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

It is commonly known that the dephasing in open quantum systems is due to the establishment of bipartite correlations with ambient environments, which are typically difficult to be fully characterized. Recently, a new approach of average over disordered Hamiltonian ensemble is developed and shown to be capable of describing the nonclassicality of incoherent dynamics based on inferring the nonclassical nature of the correlations. Here we further extend the approach of Hamiltonian ensemble in the canonical form to the realm of structural disorder. Under the variable separation of the probability distribution within the Hamiltonian ensemble, the geometrical structure 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. We show that these effects are easily understood from the emergences of additional terms in the master equations, leading to rather general master equations and, consequently, going beyond the previous frameworks of pure dephasing or isotropic depolarization.

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 23 . The process nonclassicality was exemplified 24 and quantified 29 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 23 . 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. It is also worth noting that there are efforts devoting to the decomposition of pure dephasing into random unitary (RU) representation with static probability distribution 31 . However, the generator is still time-dependent, rather than the canonical form.
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 qubit dephasing dynamics in terms of HE in the canonical form, but also demonstrates the effects of symmetry breaking of the geometrical structure on the qubit 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. Finally, we stress that to understand attainable qubit dynamics is important, particulary when we try to find the right type of HE for a physical process.

Hamiltonian ensemble in the spherical coordinates
We begin with presenting the general HE in the canonical form and exploring the corresponding ensembleaveraged 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 unitarily-evolved 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 1st directional moment �n j � � = n j �(θ, φ)d� , and the 2nd 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 2nd directional moment left, due to an appropriate orthogonal transformation. This means that, for any given HE configuration admitting variable separation (3), one is possible to redefine a new set of www.nature.com/scientificreports/ 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 Eq. (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 RU. Due to the unitality, the decohering behavior of E t is of the type of depolarization, captured by the second line of Eq. (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 29,30 , 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π . Fig. 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 1st directional moments �n j � = 0 , and the 2nd directional moments �n j n k � � = δ jk /3 , satisfying the diagonal condition of Eq. (4).
Under such highly symmetric geometry, from Eq. (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 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. www.nature.com/scientificreports/ with decay rate γ (t) = −ẇ(t)/2w(t) ; namely, the three Pauli decay channels share identical decay rate. The derivation of the master equation from Eq. (5) is outlined in Ref. 23,[32][33][34] . 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 23 , 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(ω).
As a comparative study, we consider the Gaussian radial function defined as for ω ≥ 0 . Note that the standard deviation ω c controls the width of the distribution and therefore plays similar role of cutoff frequency. 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 23 , 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 Tr[ρ 2 (t)] = w 2 G (t) + 1 /2 is shown in 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): 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 Then the radial expectation www.nature.com/scientificreports/ 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 [34][35][36][37] of the type of RU 38,39 .
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 ensembleaveraged 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. Specifically, we will consider the two examples of unbalanced regimes; namely the bagel-shaped �(θ, φ) = π −2 sin θ and the dumbbell-shaped �(θ, φ) = (3/4π) cos 2 θ . Under these symmetries, it is straightforward to see that the 1st directional moments vanish again, �n j � = 0 , and the 2nd x-and 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 2nd 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 Eq. (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)/2f j (t) − k� =jḟ k (t)/2f k (t) associated to the corresponding Pauli decay channels. Moreover, due to the azimuthal symmetry, we have �n 2 . Each decay rate can be considered as a competition between ḟ j (t)/2f j (t) 's, which are determined by the 2nd directional moments. Accordingly, we will consider the two unbalanced regimes, Moreover, the profile of purity now depends on the initial state ρ 0 = ( I + � ρ 0 ·σ )/2: Bagel-shaped solid angular function. We first show the former by considering the case of a bagelshaped 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. Now 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 (17)  www.nature.com/scientificreports/ 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. 1b) 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. 1b), 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.
The numerical results of purity are shown in Fig. 3 for initial state � ρ 0 = (sin ϑ 0 , 0, cos ϑ 0 ) with ϑ 0 = 0 (solid curves), π/4 (dotted curves), and π/2 (dashed curves), respectively. For ϑ 0 = 0 , the time evolution is solely determined by f z (t) [cf. Eq. (19)], which in turn fully determines γ x (t) and γ y (t) . Therefore, the profile of purity completely reflects the behavior of γ x (t) in Fig. 2. It is obvious that positive (negative) γ x (t) results in lowering (rising) purity, respectively; meanwhile, the singularity of γ x (t) leads to a full die-out of purity followed by a revival. On the other hand, it is more involved for ϑ 0 = π/2 . Now f x (t) dominates the time evolution. This leads to a competition between γ x (t) and γ z (t) . The singular effects of γ z (t) are quenched by γ x (t) due to the regime n 2 x � > n 2 z � and therefore the purity is always finite.

Dumbbell-shaped solid angular function.
For the second regime we consider the case of a dumbbellshaped solid angular part with �(θ, φ) = (3/4π) cos 2 θ . Its visualization is shown in Fig. 4a and clearly satisfies the desired simultaneous symmetries. In this case the 2nd directional moments are �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 � . = 0 Reciprocal square Figure 3. The purities Tr[ρ 2 (t)] for initial state � ρ 0 = (sin ϑ 0 , 0, cos ϑ 0 ) with ϑ 0 = 0 (solid curves), π/4 (dotted curves), and π/2 (dashed curves), respectively. For ϑ 0 = 0 , strong relation can be observed between the purity and γ x (t) in Fig. 2. Positive (negative) γ x (t) results in lowering (rising) purity, respectively; and the singularity of γ x (t) leads to a full die-out of purity followed by a revival. While the case of ϑ 0 = π/2 is more involved. The profile of purity is a result of the competition between γ x (t) and γ z (t) . Under the regime n 2 x � > n 2 z � , the singular effects of γ z (t) are quenched by γ x (t) and therefore the purity is always finite. www.nature.com/scientificreports/ �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. 4. Figure 4b 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. 4c. 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 profile of purity shown in Fig. 5 basically obeys the same logic as that in Fig. 3. However, as here we consider the opposite regime, smaller amplitude of γ x (t) implies larger purity for ϑ 0 = 0 . The singular effects of γ z (t) are now dominant in the competition with γ x (t) and therefore lead to a full die-out of purity for ϑ 0 = π/2.
It has been pointed out that, for qubit pure dephsing, the dephasing rate is dominated by the variance of the probability distribution within HE 23 . On the other hand, provided the simultaneous symmetries considered here, the 2nd 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. This point can be understood from the decay rates γ j (t) in Eq. (20), along with f j (t) in Eq. (6). Furthermore, according to the symmetries under consideration, we have �n 2 x � = �n 2 y � � = �n 2 z � . This implies these two regimes and demonstrates striking difference in the dynamical behavior. In the following, we will further reduce the symmetry by considering the azimuthal symmetry exclusively and explore the effects of 1st directional moment.
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 � . = 0 Reciprocal square Figure 5. The purities Tr[ρ 2 (t)] for initial state � ρ 0 = (sin ϑ 0 , 0, cos ϑ 0 ) with ϑ 0 = 0 (solid curves), π/4 (dotted curves), and π/2 (dashed curves), respectively. Aforementioned relation between the purity and γ j (t) 's in Fig. 4 can also be observed. However, due to the regime n 2 x � < n 2 z � considered here, smaller amplitude of γ x (t) implies larger purity for ϑ 0 = 0 . While for the case of ϑ 0 = π/2 , the singular effects of γ z (t) are dominant and therefore lead to a full die-out of purity.

Azimuthal symmetry
In the presence of the azimuthal symmetry, it has been shown in the previous examples that the 1st x-and y-directional moments vanish, �n x � = �n y � = 0 , and the 2nd x-and y-directional moments are equal, �n 2 x � = �n 2 y � . Crucially, in the absence of reflectional symmetry about the x-y plane, the 1st z-directional moment is typically non-vanishing, �n z � � = 0 . For example, we will consider a solid angular part of three-dimensional cardioid with �(θ, φ) = (1 − cos θ)/4π . Its visualization is shown in Fig. 6a.
According to the directional moments determined by the symmetry under consideration, the action of ensemble-averaged dynamical linear map in Eq. (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)/2f z (t) . Note that the most prominent difference of Eq. (23) from the previous examples is the presence of an effective level spacing 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 threedimensional cardioid with �(θ, φ) = (1 − cos θ)/4π . Its visualization is shown in Fig. 6a. 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 (23). 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.
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 Fig. 6b,c, 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).  Figure 6. Decoherence under the HE with the geometry of 3D cardioid. (a) Visualization of the solid angular part �(θ, φ) = (1 − cos θ)/4π of the probability distribution. This geometry exhibits the azimuthal symmetry exclusively; meanwhile, this leads to the balanced regime �n 2 x � � = �n 2 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.

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 1st x-and y-directional moments, �n x � = �n y � = 0 . The three 2nd directional moments are generically different due to the lack of the azimuthal symmetry. Given the vanishing 1st x-and y-directional moments, the action of ensemble-averaged dynamical linear map in Eq. (5) is formally the same as Eq. (22)  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. 7a 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 �(θ, φ) , we have �n z � � = −1/3 , �n 2 x � � = (2 + a)/6 , �n 2 y � � = (2 − a)/6 , and �n 2 z � � = 1/3 . Along with the radial expectations cos ωt P shown in Eqs. (11), (14), and (17), as well as sin ωt P in Methods, the numerical results of γ xy (t) (25) are shown in Fig. 7.
The results for the Gaussian radial function are represented in Fig. 7b by black (solid and dashed) curves with varying lateral asymmetry a. 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 varying a. 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. 7c by blue curves with varying a. 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 The decay rate γ xy (t) for the Gaussian (black) and the exponential cutoff (red) radial functions with different lateral asymmetry a. 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 different lateral asymmetry a. 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. www.nature.com/scientificreports/ 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.

Conclusion and discussion
In summary, we have explored the incoherent dynamics of qubits raised from the ensemble average over the canonical HE of structural disorder. Under the variable separation (3) 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 23 , 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 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. Generally speaking, the more asymmetric the geometry is, the more terms emerge.
Notably, the asymmetry goes into play via the 1st and 2nd directional moments, which in turn determine the aforementioned terms in the master equation. The effective level spacing ω(t) emerges in the presence of 1st 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 2nd 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 2nd directional moments (balanced/unbalanced regimes) and the amplitudes of the decay rates. Moreover, the presence of 1st 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 23 , 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 worthwhile to note that the variable separation (3) can be further released by considering correlated radial and angular coordinates, leading to even more general joint probability distributions p(ω, θ , φ) . For example, expansion in terms of spherical harmonics Y l,m (θ, φ) has been demonstrated in Ref. 30 . Furthermore, this problem can also be investigted from the viewpoint of RU, which is highly related to HE. Since the RU representation is a useful tool in the study of open system dynamics 31,38,39 , considerable efforts have been devoted into the investigation of UR decomposition [40][41][42] . For example, master equation in the same form as Eq. (20) is derived by mixing fluctuating Gaussian noise along different angles 39 .
On the other hand, it is know that any qubit or qutrit phase damping dynamical mpas are RU 43,44 . However, neither the RU-decomposition nor the extreme points for higher dimensional cases are unclear. These problems can in general merely be numerically implemented 43 . Therefore, the canonical HE would not only provide insights into the effects of higher asymmetry beyond the framework of existing RU representation, but also establish a systematic approach to single out the attainable subset, and substantially simplify the problem of RU-decomposition. This underpins the significance of our approach.

Methods
Orthogonal transform of the 2nd 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 2nd directional moment �n j ′ n k ′ � forms a symmetric matrix, it can be diagonalized with an appropriate orthogonal matrix û = [u j ′ j ] such that 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).