Effects of symmetry breaking of the structurally-disordered Hamiltonian ensembles on the anisotropic decoherence


 It is commonly known that the dephasing in open quantum systems is due to the establishment of bipartite correlations with an environment.
Recently, a new approach of average over disordered Hamiltonian ensemble is developed and shown to capable of describing
both the incoherent dynamical behavior and the nonclassicality of dynamical processes.
Here we further extend the approach of Hamiltonian ensemble in the canonical form to the realm of structural disorder.
Under the separation of the probability distribution within the Hamiltonian ensemble,
the geometrical structural is easily visualized and can be characterized according to the degree of symmetry.
We demonstrate four degrees and investigate the effects of different types of symmetry breaking on the incoherent dynamics.
With these effects, we obtain rather general master equations, going beyond the previous frameworks of pure dephasing or isotropic depolarization.
The practicality of the Hamiltonian ensemble and the theory of process nonclassicality is significantly enhanced.


Introduction
Exposed to the inevitable interactions with the huge surrounding environments, any quantum systems generically undergo incoherent dynamical processes and gradually lose their quantumness [1][2][3][4] , constituting the primary obstacle in the developments of frontier quantum technologies [5][6][7][8][9][10] . Therefore, it is crucial to characterize, control, and eliminate the sources of decoherence. One of the main causes of the incoherent dynamical nature arises from the damage to the system-environment correlations established during their interactions. However, due to the huge environmental degree of freedom, it is infeasible to fully access these bipartite correlations. This renders a general, and precise, description of how the correlations are destroyed highly nontrivial. Consequently, their incoherent effects on the reduced system dynamics are taken into account in terms of a family of completely positive and trace-preserving (CPTP) dynamical linear maps [11][12][13][14] . Moreover, there are several alternative techniques for characterizing CPTP maps have been developed, such as the Kraus operators 15 , the process matrices 16 , and the Choi-Jamiołkowski isomorphism 17,18 .
Recently, a promising approach, referred to as Hamiltonian ensemble (HE), has been developed 19,20 . HE is initially dedicated to the investigation of disordered systems and classical noises [21][22][23][24] , which are described by an ensemble of Hermitian operators parameterized by some random variables obeying a specified probability distribution function. Irrespective of the unitarity of a single realization generated by each member Hermitian operator, the time-evolved state ρ(t) undergoes a dephasing dynamics after the ensemble-average procedure over all unitary realizations 19,20 . Furthermore, it has been pointed out that the incoherent dynamical behavior is intimately related to the properties of the probability distribution function 19 . Based on this insight, the probability distribution function encapsulated within the HE has attracted exclusive focus in the characterization of dephasing dynamics, and been promoted to be the Canonical Hamiltonian ensemble representation (CHER) of dynamical processes in the frequency domain 25 . Additionally, the CHER has been shown to be a versatile approach in the quantification of process nonclassicality 25 .
However, due to the difficulty imposed by the non-abelian algebraic structure underlying the canonical HE, most of the aforementioned works relied on significant simplifications to circumvent it. For example, the spectral disorder appealed to HEs consisting of diagonal operators, leading to the pure dephasing 19 . The process nonclassicality was exemplified 20 and quantified 25 under the same framework of Cartan subalgebra. Furthermore, the attempt going beyond the framework of pure dephasing was the unitarily invariant disorder, which was studied by incorporating the spectral disorder with the Haar measure integral, leading to the isotropic depolarization 19 . Whereas, the most drawback of this unitary invariance approach lies in the deviation from the canonical form, giving rise to the issue of double counting of each member Hermitian operator.
To cure this issue, as well as to go further beyond the above two frameworks of pure dephasing and isotropic depolarization, here we study the canonical HE of structural disorder along with different degrees of symmetry. This not only enables us to explore more general types of dephasing dynamics in terms of HE in the canonical form, but also demonstrates the effects of symmetry breaking of the geometrical structure on the incoherent dynamics. The unitarily invariant disorder is shown to be a special case of the spherical symmetry. We have reduced the continuous spherical symmetry to three lower levels, until the discrete one of simultaneous reflectional symmetries. Each symmetry breaking gives rise to an additional disturbance complicating the dynamics, including the anisotropic decay rates, the effective level spacing, and the off-diagonal decay rates. Following this line, the pattern of any further generalization can be deduced.

Hamiltonian ensemble in the spherical coordinates
We begin with presenting the general HE in the canonical form and exploring the corresponding ensemble-averaged dynamics for a qubit. Generically, any Hermitian operators acting on a qubit system are elements in the u(2) Lie algebra of the form H λ = (λ 0 I + λ xσx + λ yσy + λ zσz )/2. However, as λ 0 plays no role in describing the qubit dynamics due to the commutativity [ I,σ j ] = 0, we can restrict ourselves to the traceless operators in the su(2) Lie algebra with λ 0 = 0.
Here we consider the Hamiltonian ensemble {(p λ , H λ )} λ of canonical form parameterized by λ ∈ R 3 , where each member Hamiltonian H λ ∈ su(2) is associated with a probability p λ of occurrence. Since the su(2) Lie algebra is non-abelian, the corresponding unitary time-evolution operator U λ = exp(−i H λ t) is rather difficult to be dealt with in these parameters. Crucially, along with the change of variables to the spherical coordinates λ x = ω sin θ cos φ , λ y = ω sin θ sin φ , and λ z = ω cos θ , the HE can be recast into where n = (sin θ cos φ , sin θ sin φ , cos θ ) ∈ R 3 is the directional unit vector andσ denotes the three Pauli matrices. Accordingly, each operator U λ = cos(ωt/2) I − i sin(ωt/2) n ·σ is explicitly expressed in the spherical coordinates and leads to a unitarilyevolved single realization provided an initial state ρ 0 . In view of the member Hamiltonian encapsulated within the HE (1), as well as the single realization (2), we can observe a separation between the radial coordinate, ω, and the solid angular coordinates, (θ , φ ). Accordingly, we will further assume a separable probability distribution into two positive real functions P(ω) and Θ(θ , φ ). We will show that many interesting properties can be conveniently studied with the help of this separation; particularly, the symmetry is easily exhibited by the geometrical structure of Θ(θ , φ ). One should note that, it is p λ a legitimate probability distribution and normalized to unity p λ d 3 λ = p(ω, θ , φ )ω 2 dωdΩ = 1 with dΩ = sin θ dθ dφ ; while this is not the case for P(ω) or Θ(θ , φ ) individually. We therefore assume that ∞ 0 P(ω)ω 2 dω = 1/ξ and Θ(θ , φ )dΩ = ξ due to the separability. This guarantees the normalization condition for p(ω, θ , φ ).
Based on the separability, the ensemble-averaged dynamics under the structurally-disordered HE (1) is given by where, for convenience, we have introduced the radial expectation f (ω) P = ∞ 0 f (ω)P(ω)ω 2 dω with respect to P(ω), n j the three components of n, the 1 st directional moment n j Θ = n j Θ(θ , φ )dΩ, and the 2 nd directional moment n j n k Θ = n j n k Θ(θ , φ )dΩ. It should be noted that, in the last term of Eq. (4), there are only three square terms of the 2 nd directional moment left, due to an appropriate orthogonal transformation. This means that, for any given HE configuration, one is possible to redefine a new set of axes of the geometry by using a basis transformation to eliminate the crossing terms n j n k Θ with j = k. Consequently, without loss of generality, we can start from the dynamical linear map (4) and ignore the crossing terms in the following. This significantly simplifies the complexity of the problem. Further details on the orthogonal transformation are discussed in Methods.
On the other hand, the density matrix for a qubit system ρ = ( I + ρ ·σ )/2 is also an element in the u(2) Lie algebra parameterized by the corresponding Bloch vector ρ ∈ R 3 . Therefore, the properties of the dynamical linear map (4) can be fully understood from its action on the generators of the u(2) Lie algebra: where and ε jkl = 1 if ( j, k, l) is an even permutation of (1, 2, 3), −1 if an odd permutation, and 0 otherwise. From the first line of Eqs. (5), we can see that E t is unital. This can be understood if we note that the HE is a special case with time-independent probability distribution and member Hamiltonian of a superset of random unitary (RU). Due to the unitality, the decohering behavior of E t is of the type of depolarization, captured by the second line of Eqs. (5). From the later, we can see that the decohering behavior of E t depends highly on the structure and the symmetry of p(ω, θ , φ ) via the radial expectations and the directional moments. On the other hand, this also reflects the notion of CHER 25 , which conceives the (quasi-)distribution function as the characteristic representation of a dynamics. Accordingly, we will investigate the dynamical behavior of E t along with different degrees of symmetry of the probability distribution exhibited by the solid angular part Θ(θ , φ ).

Spherical symmetry
We first consider the case of spherically symmetric probability distribution with Θ(θ , φ ) = 1/4π. Figure 1a shows its visualization with the distance between the surface and the origin indicating the value of Θ(θ , φ ) in solid angular coordinates (θ , φ ). As Θ(θ , φ ) is a constant, the geometry forms a sphere; therefore, the HE is of spherical symmetry, which is the highest symmetry can be exhibited. It is straightforward to verify that ∞ 0 P(ω)ω 2 dω = Θ(θ , φ )dΩ = 1, the 1 st directional moments n j Θ = 0, and the 2 nd directional moments n j n k Θ = δ jk /3, satisfying the diagonal condition of Eq. (4).
Under such highly symmetric geometry, from Eqs. (5), the ensemble-averaged dynamics is given by which is a statistical mixture between the completely mixed state and the initial state ρ 0 along with a time-varying wight w(t) = (2 cos ωt P + 1)/3. Therefore, the initial state will gradually lose its coherence. Moreover, this incoherent dynamical behavior is governed by the master equation of isotropic depolarization with decay rate γ(t) = −ẇ(t)/2w(t); namely, the three Pauli decay channels share identical decay rate. The derivation of the master equation from Eqs. (5) is outlined in Ref. 19,[26][27][28] . Additionally, the decohering dynamics of this isotropic depolarization can also be well-understood by the purity It is worthwhile to note that, the spherical symmetry qualitatively reproduces the unitarily invariant disorder for qubit 19 , which is studied by means of incorporating the spectral disorder with the Haar measure distributed uniformly over whole solid angular coordinates. However, one of the drawbacks of this approach is the occurrence of double counting of each member Hamiltonian operator. For example, both of the two identical members ωσ z /2 and exp(−iπσ x /2)(−ωσ z /2) exp(iπσ x /2) contribute individually to the ensemble-average procedure. On the other hand, the HE in the canonical form rules out this circumstance by considering positive radial coordinate exclusively. Therefore, quantitative differences can be seen when we further take the radial part P(ω) into account in the following. To further clarify this situation and to exemplify the isotropic depolarization, we will carry out several types of distribution functions. As the solid angular part Θ(θ , φ ) has been specified according to the symmetry, we will characterize the distribution functions with the radial part P(ω).

Figure 1.
Decoherence under the HE with spherically symmetric probability distribution. (a) Visualization of the solid angular part Θ(θ , φ ) = 1/4π of the probability distribution. As Θ(θ , φ ) is a constant, the geometry is a sphere of radius 1/4π. (b) Time-evolution of the decay rates γ(t) for the Gaussian (black solid curve), the exponential cutoff (red dotted curve), and the reciprocal square (blue dashed cruve) radial functions. Their behaviors are different to each other, as explained in the main text. The decay rates temporarily go down to negative values in some time periods, indicating transitions between indivisibility and full divisibility of the type of random unitary. (c) The purities Tr[ρ 2 (t)] governed by the isotropic depolarization with pure initial states for the three radial functions. There is a rivival in the time period of negative decay rate. Additionally, the final value of the purities is 5/9, rather than 1/2, reflecting that the final state is not the completely mixed state.
As a comparative study, we consider the Gaussian radial function defined as for ω ≥ 0, where ω c plays the role of standard deviation and controls the width of the distribution. Note that the functional form is specified according to the normalization condition ∞ 0 P G (ω)ω 2 dω = 1 and therefore the coefficient is slightly different from the usual one. For the Gaussian radial function P G (ω), one can analytically evaluate This allows us to evaluate the mixing weight w G (t) = (2 cos ωt P G + 1)/3 and the decay rate in Eq. (8): If compared with the one obtained by the unitarily invariant disorder for qubit 19 , we can find that Eq. (12) is lesser by a factor 2, indicating that the HE in the canonical form rules out the circumstance of double counting. The numerical results are shown in Fig. 1b. It can be seen that, γ G (t) is initially increasing, then followed by a sharp descent to negative values, and finally approaching zero asymptotically from below. Additionally, the time evolution of the purity with pure initial state Fig. 1c. The purity initially decays very approaching the completely mixed state due to the increasing γ G (t); afterwards, it rises again and saturates to the value of (w 2 G (t → ∞) + 1)/2, indicating that the final state deviates from the completely mixed state. This property can also be observed in the following two examples and will be discussed latter.
Next, we consider an alternative cutoff in exponential form defined as satisfying the normalization condition ∞ 0 P EC (ω)ω 2 dω = 1. ω c is the cutoff value. With this exponential cutoff, we have the mixing weight w EC (t) = (2 cos ωt P EC + 1)/3, and the decay rate in Eq. (8):

4/13
From Fig. 1b, γ EC (t) exhibits a similar behavior to γ G (t), but even sharper oscillation. Consequently, this is also the case for the purity Tr[ρ 2 (t)] = w 2 EC (t) + 1 /2 shown in Fig. 1c. This is due to the fact that P EC (ω) possesses a long wing over high ω domain, leading to a rapid spreading of random unitary rotation and shorter coherence time.
As a third example, we consider a radial function of different type. The reciprocal square for ω ∈ [0, ω c ], is defined on a finite interval rather than infinite length. Again the functional form is specified by the normalization condition ω c 0 P RS (ω)ω 2 dω = 1. Then the radial expectation leads to the mixing weight w RS (t) = (2 cos ωt P RS + 1)/3 and the decay rate in Eq. (8): Due to the finite domain of P RS (ω), the decay rate γ RS (t) exhibits a long oscillating tail as shown in Fig. 1b. This behavior is very different from γ G (t) and γ EC (t), and renders the purity Tr[ρ 2 (t)] = w 2 RS (t) + 1 /2 oscillating as well after the initial descent.
It is evident from the above three examples that, whenever the decay rates γ(t) go down to negative values, there is a revival of the purity Tr[ρ 2 (t)] in the time period of negative decay rate. We can even clearly observe this phenomenon from the long oscillating tails in the reciprocal square example. This can be explained as a typical transition between indivisibility and full divisibility [28][29][30] of the type of random unitary.
On the other hand, although the dynamics is governed by the master equation (8) of isotropic depolarization, the final state is generically not the completely mixed state. This can be realized by observing that the purity saturates to a value of 5/9, rather than 1/2 in Fig. 1c. This does not imply the violation of unitality of the ensemble-averaged dynamics under canonical HE. Since cos ωt ω 3 /3 → 0 when t → ∞ in these examples, we have the steady mixing weight w(t → ∞) = 1/3 and, consequently, the final state ρ(∞) = 2( I/2) + ρ 0 /3, a constant mixture between the completely mixed state and the initial state ρ 0 .
We have shown that the spherically symmetric probability distribution qualitatively reproduces the unitarily invariant disorder for qubit and leads to the master equation of isotropic depolarization with identical decay rate of the three Pauli decay channels. To demonstrate the versatility of the canonical HE in characterizing the incoherent dynamics beyond isotropic depolarization, we will further reduce the degree of symmetry and explore its effects on the incoherent dynamical behaviors.

Simultaneous azimuthal and reflectional symmetries
To release the spherical symmetry, we consider the HE exhibiting simultaneously the azimuthal symmetry and the reflectional symmetry about the x-y plane. Under these symmetries, it is straightforward to see that the 1 st directional moments vanish again, n j Θ = 0, and the 2 nd xand y-directional moments are equal, n 2 x Θ = n 2 y Θ . However, in contrast to the case of spherical symmetry, they may not necessarily equal to the 2 nd z-directional moment n 2 z Θ ; meanwhile, without loss of generality, we assume that the diagonal condition of Eq. (4) is still hold, i.e., n j n k Θ = 0 for j = k.
After determining the directional moments according to the symmetries, the action of ensemble-averaged dynamical linear map in Eqs. (5) is significantly simplified as Since the spherical symmetry is no longer hold, the actions of E t on the three generators are different. Consequently, the incoherent dynamical behavior is governed by the master equation of anisotropic depolarization with decay rates γ j (t) = ḟ j (t)/2 f j (t) − ∑ k = jḟk (t)/2 f k (t) associated to the corresponding Pauli decay channels. Moreover, due to the azimuthal symmetry, we have n 2 x Θ = n 2 y Θ , f x (t) = f y (t), and, consequently, γ x (t) = γ y (t) = −ḟ z (t)/2 f z (t). Each 5/13 Figure 2. Decoherence under the HE with bagel-shaped geometry. (a) Visualization of the solid angular part Θ(θ , φ ) = π −2 sin θ of the probability distribution. The bagel-shaped geometry exhibits the azimuthal and reflectional symmetries simultaneously; meanwhile, this leads to the regime n 2 The decay rates γ x (t) (solid curves) and γ z (t) (dashed curves) for the Gaussian (black) and the exponential cutoff (red) radial functions. Each γ x (t) exhibits two singularities and finally approaches zero asymptotically. (c) The decay rates γ x (t) (solid curves) and γ z (t) (dashed curves) for the reciprocal square radial function. γ x (t) shows a similar oscillating behavior; while γ z (t) possesses more zeros. Additionally, the amplitude of γ x (t) is larger than γ z (t), reflecting the relationship n 2 x Θ > n 2 z Θ .
decay rate can be considered as a competition between ḟ j (t)/2 f j (t) 's, which are determined by the 2 nd directional moments. Accordingly we will consider the two unbalanced regimes, n 2 x Θ = n 2 y Θ > ξ /3 > n 2 z Θ and n 2 x Θ = n 2 y Θ < ξ /3 < n 2 z Θ . We first show the former by considering the case of a bagel-shaped solid angular part with Θ(θ , φ ) = π −2 sin θ . Its visualization is shown in Fig. 2a, from which it is obvious that Θ(θ , φ ) exhibits simultaneously the azimuthal symmetry and the reflectional symmetry about the x-y plane. With the specified functional form, we can explicitly compute n 2 x Θ = n 2 y Θ = 3/8, n 2 z Θ = 1/4, and ξ = Θ(θ , φ )dΩ = 1, satisfying the required relationship n 2 x Θ = n 2 y Θ > 1/3 > n 2 z Θ . In fact, this relationship can also be inferred from the visualization of Θ(θ , φ ) before explicit computations.
After specifying the solid angular part, we revisit the three radial functions, the Gaussian P G (ω) (10), the exponential cutoff P EC (ω) (13), and the reciprocal square P RS (ω) (16). We first verify that ∞ 0 P(ω)ω 2 dω = Θ(θ , φ )dΩ = 1 guarantees the normalization condition for p(ω, θ , φ ), and the radial expectation cos ωt P with respect to the three radial functions have been shown in Eqs. (11), (14), and (17), respectively. Then the incoherent dynamical behavior can be fully understood and the decay rates γ j (t) in the master equation (20) can also be computed explicitly. The analytical expressions of γ j (t) are given in Methods. We show the numerical results with respect to the radial functions in Fig. 2. Figure 2b shows γ x (t) (solid curves) and γ z (t) (dashed curves) for the Gaussian (black) and the exponential cutoff (red) radial functions. Note that, each of the sharp descents of γ x (t) for the Gaussian and the exponential cutoff under spherical symmetry [cf. Fig. (1)b] now splits into two prominent singularities under lower symmetry. This renders the behavior of γ z (t) = − ḟ x (t)/ f x (t) − γ x (t) singular as well. Finally, they again approach zero asymptotically. The results for the reciprocal square radial function are shown in Fig. 2c. In contrast to the former, both γ x (t) and γ z (t) are regular. It can be seen that γ x (t) shows similar temporal behavior to the one under spherical symmetry [cf. Fig. (1)b], whereas γ z (t) possesses more zeros. Additionally, it is interesting to note that the amplitude of γ x (t) is larger than γ z (t), reminiscent of the relationship n 2 x Θ > n 2 z Θ we are considered. Similar analogy can be observed in the latter example. In fact, this analogy provides some insights into the effects of symmetry breaking on the decay rates, and will be discussed latter.
For the second regime we consider the case of a dumbbel-shaped solid angular part with Θ(θ , φ ) = (3/4π) cos 2 θ . Its visualization is shown in Fig. 3a and clearly satisfies the desired simultaneous symmetries. In this case the 2 nd directional moments are n 2 x Θ = n 2 y Θ = 1/5, n 2 z Θ = 3/5, and ξ = Θ(θ , φ )dΩ = 1, satisfying the required relationship n 2 x Θ = n 2 y Θ < 1/3 < n 2 z Θ . Similarly, we adopt the same radial functions again. The analytic expressions are given in Methods. We show the numerical results in Fig. 3. Figure. 3b shows the results of γ x (t) (solid curves) and γ z (t) (dashed curves) for the Gaussian (black) and the exponential cutoff (red) radial functions. Under this regime γ x (t)'s are finite with amplitudes smaller than the ones under spherical symmetry [cf. Fig. 1b], whereas γ z (t)'s are still singular, reflecting the relationship n 2 x Θ < n 2 z Θ as well. On the other hand, γ x (t) for the Gaussian exhibits a similar temporal behavior to the one under spherical symmetry; while γ x (t) for exponential cutoff exhibits an additional mild rising from negative and finally approaches zero asymptotically from above. The results for the reciprocal square radial function are shown in Fig. 3c. We can see that both of the decay rates exhibits the same temporal behavior and, consequently, possess the same zeros. Meanwhile, We again see the same analogy between the amplitudes of γ j (t)'s and the relationship n 2 x Θ < n 2 z Θ .
The decay rates γ x (t) (solid curves) and γ z (t) (dashed curves) for the Gaussian (black) and the exponential cutoff (red) radial functions. γ x (t)'s are finite with smaller amplitudes, whereas γ z (t)'s are singular. Furthermore, γ x (t) for the exponential cutoff approaches zero asymptotically from above after a mild rising. (c) The decay rates γ x (t) (solid curves) and γ z (t) (dashed curves) for the reciprocal square radial function. Both of the decay rates exhibits the same temporal behavior. Finally, in these plots, we can see the analogy between the amplitudes of γ j (t)'s and the relationship n 2 x Θ < n 2 z Θ .
It has been pointed out that, for qubit pure dephsing, the dephasing rate is dominated by the variance of the probability distribution within HE 19 . On the other hand, provided the simultaneous symmetries considered here, the 2 nd directional moments play the role of variance alone the specified directions. Consequently, the asymmetry goes into play via them, as reflected by the decay rates, and gives rise to the analogy.

Azimuthal symmetry
We further reduce the symmetry by considering the azimuthal symmetry exclusively and explore the resulting master equations. In the presence of the azimuthal symmetry, it has been shown in the previous examples that the 1 st xand y-directional moments vanish, n x Θ = n y Θ = 0, and the 2 nd xand y-directional moments are equal, n 2 x Θ = n 2 y Θ . Crucially, in the absence of reflectional symmetry about the x-y plane, the 1 st z-directional moment is typically non-vanishing, n z Θ = 0.
According to the directional moments determined by the symmetry under consideration, the action of ensemble-averaged dynamical linear map in Eqs. (5) can be explicitly written as In the presence of n z Θ , we have one more additional radial expectation sin ωt P ; meanwhile, the incoherent dynamical behavior is governed by the master equation In deriving the above master equation, we have used the facts that n 2 x Θ = n 2 y Θ and f x (t) = f y (t); therefore, γ x (t) = γ y (t) = −ḟ z (t)/2 f z (t). Note that the most prominent difference of Eq.
is also altered in the presence of finite n z Θ . Both of ω(t) and the variation of γ z (t) are results of the lack of the reflectional symmetry about the x-y plane.
Instead of the unbalanced regimes considered in the two previous examples, now we demonstrate a balanced one with n 2 x Θ = n 2 y Θ = 1/3 = n 2 z Θ to simplify the complexity. We consider a solid angular part of three-dimensional cardioid with Θ(θ , φ ) = (1 − cos θ )/4π. Its visualization is shown in Fig. 4a. We can verify that the desired balanced regime is satisfied, and compute n x Θ = n y Θ = 0 and n z Θ = −1/3. With the same radial functions, we can determine the effective level spacing ω(t) and the decay rate γ j (t) in the master equation (22). In addition to cos ωt P , now we need one more radial expectation sin ωt P . The analytical expresstions with respect to the three radial functions are given in Methods. x Θ = n 2 y Θ = 1/3 = n 2 z Θ . (b) The averaged level spacing ω(t) for the Gaussian (black solid curve), the exponential cutoff (red dotted curve), and the reciprocal square (blue dashed cruve) radial functions. Its presence is a result of the lack of the reflectional symmetry about the x-y plane. (c) The decay rate γ z (t) for the Gaussian (black solid curve), the exponential cutoff (red dotted curve), and the reciprocal square (blue dashed cruve) radial functions. The line shape is similar to the one under spherical symmetry appended by an additional drop before the peaking value, which is a result of the asymmetry.
Thanks to the balanced regime under consideration, we find that γ x (t) = γ y (t) and they are exactly the same as those under spherical symmetry; while this is not the case for γ z (t) due to the breaking of the reflectional symmetry. We therefore show the numerical results of ω(t) and γ z (t) in Figs. 4b and 4c, respectively. For the Gaussian radial function (black solid curve), ω(t) begins with a negative value. After a shallow drop, a sharp peak is followed, then going down to negative again, and finally approaching zero asymptotically from below. For the exponential cutoff (red dotted), the line shape is similar to that of Gaussian but an overturned one. However, for the reciprocal square (blue dashed curve), ω(t) exhibits a long oscillating tail due to the finite domain of P RS (ω). On the other hand, it is interesting to note that, the overall line shape of γ z (t) is similar to the one under spherical symmetry, but an additional drop before reaching the peaking value. The similarity is a consequence of the balanced regime, according to the aforementioned analogy. However, the asymmetry is also influential by perturbing γ z (t).

Simultaneous reflectional symmetries
Finally, we consider a more asymmetric case by reducing the continuous azimuthal symmetry to a discrete one, i.e., the simultaneous reflectional symmetries about the x-z and y-z planes. In this case, we can only deduce the vanishment of the 1 st xand y-directional moments, n x Ω = n y Ω = 0. The three 2 nd directional moments are generically different due to the lack of the azimuthal symmetry.
Given the vanishing 1 st xand y-directional moments, the action of ensemble-averaged dynamical linear map in Eqs. (5) is formally the same as Eqs. (21) under azimuthal symmetry. Whereas the corresponding master equation possesses two more decay channels with the off-diagonal decay rates where D(t) = f x (t) f y (t) + n z 2 Θ sin ωt 2 P . It is clear that their presences are a result of the azimuthal symmetry breaking. As this is the lowest degree of symmetry considered here, Eq. (23) is the most general master equation we have demonstrated. It can be reduced to all the previous ones if the corresponding symmetries are recovered.
To exemplify this case of low symmetry, we consider a solid angular part of kneaded cardioid with Θ(θ , φ ) = (1 − cos θ )(1 + a cos 2φ )/4π, where 0 ≤ a ≤ 1 describes the degree of lateral asymmetry of the geometry. Its visualization is shown in Fig. 5a with a = 0.3. Due to the φ -dependence, the geometry appears to be subject to stress along the y-axis, and then expanding along the x-axis. Therefore the azimuthal symmetry is broken and merely the simultaneous reflectional symmetries are left. For this The decay rate γ xy (t) for the Gaussian (black) and the exponential cutoff (red) radial functions with decreasing lateral asymmetry a from 0.3 to 0.1 and from 0.7 to 0.1 (decreasing opacity), respectively. As γ xy (t) is caused by the azimuthal symmetry breaking, it is gradually vanishing when a is decreasing and, contrarily, becomes more prominent for a more asymmetric geometry. As indicated by the dashed curves, the negative peaks even split into two singularities at a = 0.3 and 0.7 for Gaussian and exponential cutoff radial functions, respectively. (c) The decay rate γ xy (t) for the reciprocal square radial function with decreasing lateral asymmetry a from 0.7 to 0.1 (decreasing opacity). Similar tendency toward vanishment with decreasing a can also be seen. However, for the reciprocal square, γ xy (t) do not exhibit singularity even if under large asymmetry.
The results for the Gaussian radial function are represented in Fig. 5b by black (solid and dashed) curves with decreasing lateral asymmetry a from 0.3 to 0.1 (decreasing opacity). As expected, γ xy (t) is gradually vanishing with decreasing a (i.e., more symmetric geometry), reflecting the fact that γ xy (t) is a result of azimuthal symmetry breaking. In contrast, γ xy (t) becomes more prominent when a is large. The negative peak even splits into two singularities at a = 0.3, as indicated by the black dashed curve. Moreover, the results for the exponential cutoff radial function are represented by red (solid and dashed) curves with decreasing a from 0.7 to 0.1 (decreasing opacity). It can be seen that the overall tendency is the same as Gaussian but an overturned one. The curves show a positive peak at an earlier time; meanwhile, the exponential-cutoff case shows singularities at a larger value of a = 0.7, as indicated by the red dashed curve.
On the other hand, the results for the reciprocal square radial function are represented in Fig. 5c by blue curves with decreasing a from 0.7 to 0.1 (decreasing opacity). The tendency toward vanishment with decreasing a can also be seen; whereas, in contrast to the other two radial functions, the deformation of the line shape with a is even, without showing singularity, even if a is large. Moreover, akin to the other decay rates with the same reciprocal square radial function, γ xy (t) also exhibits an oscillating tail, after the negative main peak. Crucially, from these results of Gaussian and reciprocal square radial functions, there is a distinct property the off-diagonal decay rate γ xy (t) can be observed as well. To guarantee the complete positivity of the dynamical linear map E t , the short time behavior of all the digonal decay rates γ j (t) associated to the Pauli decay channels is rising at beginning. Any negative values can only emerge after a period of positive values. However, γ xy (t) is possible to have a short time behavior of descent to negative at beginning. We stress that this does not imply the violation of complete positivity as the eigenvalues of the Kossakowski matrix K = [γ jk (t)] exhibit a short time behavior of rising to positive at beginning.

Conclusions
In summary, we have explored the incoherent dynamics raised from the ensemble average over the canonical HE of structural disorder. Under the separation of the probability distribution within the HE into the radial and the solid angular parts, the structural disorder can be characterized in accordance with the degree of symmetry of the solid angular geometry. The effects of asymmetry are particularly manifest via the corresponding master equation in the Lindblad form.
In this work we have considered three radial functions; namely, the Gaussian, the exponential cutoff, and the reciprocal square. We first show the most symmetric case of spherical symmetry, leading to the master equation of isotropic depolarization. It is worthwhile to note that the canonical HE of spherical symmetry not only reproduces the results of unitarily invariant disorder for qubit 19 , but also resolves the issue of double counting. Further versatility of the canonical HE in describing various types of incoherent dynamics beyond isotropic depolarization is revealed by considering asymmetric cases. In addition to the 9/13 spherical symmetry, we have also demonstrated three lower degrees of symmetry. We can observe the effects caused by the symmetry breaking of different types via the master equation in the Lindblad form, including the singularities of the decay rates, the effective level spacing, and the off-diagonal decay rates. The more asymmetric the geometry is, the more terms emerge.
Notably, the asymmetry goes into play via the 1 st and 2 nd directional moments, which in turn determine the aforementioned terms in the master equation. The effective level spacing ω(t) emerges in the presence of 1 st directional moments, which are caused by the reflectional symmetry breaking. They denote the mean values alone the specified directions. While the time dependence of ω(t) is determined by the 2 nd directional moments and the radial expectations, characterizing the overall asymmetry and the radial disorder, respectively. On the other hand, the behaviors, particularly the time dependence, of the decay rates are also intimately related to the directional moments. This can be easily understood from the analogy between the relationship of the 2 nd directional moments (balanced/unbalanced regimes) and the amplitudes of the decay rates. Moreover, the presence of 1 st directional moments will further complicate the decay rates. Finally, the highly asymmetric case gives rise to the emergence of off-diagonal decay rates. Interestingly, very similar conclusions have been drawn under the framework of qubit pure dephasing 19 , where it has been shown that ω(t) is given by the mean value and the asymmetry of the distribution, while the dephasing rate is given by the variance and the kurtosis.
Finally, it is intriguing to revisit the superset of RU. Considerable efforts have been devoted into the investigation of decomposition of unital dynamical maps into RU [31][32][33] . It is know that any qubit or qutrit phase damping dynamical mpas are RU 34,35 . However, neither the RU-decomposition nor the extreme points for higher dimensional cases are unclear. These problems can in general merely be numerically implemented 34 . Therefore, the canonical HE would provide a systematical approach to single out the attainable subset, and significantly simplify the problem of RU-decomposition.

Methods
Orthogonal transform of the 2 nd directional moments According to the last term of Eq. (2) in the main text, there should be 9 terms in the last summation of Eq. (4): Since the 9 terms of 2 nd directional moment n j n k Θ forms a symmetric matrix, it can be diagonalized with an appropriate orthogonal matrixû = [u j j ] such that n j n k Θ = 3 ∑ j,k=1 u j j δ jk n 2 j Θ u k k .
This corresponds to a basis transformationσ j = ∑ 3 j =1 u j jσ j and a rotation n j = ∑ 3 j =1 u j j n j of the directional unit vector. The later is necessary in the first summation of Eq. (4).

Decoherence under simultaneous azimuthal and reflectional symmetries
As explained in the main text, given the probability distribution p(ω, θ , φ ) = P(ω)Θ(θ , φ ) with the solid angular part Θ(θ , φ ) exhibiting simultaneously the azimuthal symmetry and the reflectional symmetry about the x-y plane, we can compute the decay rates in the master equation (20) of anisotropic depolarization, where f j (t) = cos ωt P ξ − n 2 j Θ + n 2 j Θ /ξ . In these expressions, we have used the facts that n 2 x Θ = n 2 y Θ and f x (t) = f y (t) due to the azimuthal symmetry.
To fully determine the probability distribution p(ω, θ , φ ) = P(ω)Θ(θ , φ ), we revisit the three radial functions defined in Eqs. (10), (13), and (16) in the main text. Then the radial expectation cos ωt P with respect to the three radial functions are