Effective dynamics of qubit networks via phase-covariant quantum ensembles
Sean Prudhoe Unnati Akhouri Tommy Chin Sarah Shandera






















1-4mu l1-4mu l
1-4.5mu l1-5mu l
⟨
⟩
equation






Institute of Gravitation and the Cosmos, The Pennsylvania State University,

Whitmore 321, University Park, PA, USA
Department of Physics, The Pennsylvania State University,

Davey 104, University Park, PA, USA

Effective dynamics of qubit networks
May 1, 2024
Sean Prudhoe, Unnati Akhouri, Tommy Chin, and Sarah ShanderaEffective dynamics of qubit networks via phase-covariant quantum ensembles



We study ensembles of phase-covariant channels. We show that such ensembles arise naturally from familiar spin-chain models (e.g., XXZ) with a special class of initial states, and that the disorder-averaged map of disordered spin chains is phase-covariant under a weak symmetry constraint on the distribution. We use those examples to motivate a broader class of phase-covariant ensembles, which include both unital and non-unital channels. We demonstrate the physical properties captured by the late-time limit of the average map over the ensemble.










§ INTRODUCTION


Most of the quantum systems used for measurement and computation are open, rather than isolated, systems. We are frequently interested in understanding the dynamics of an ensemble of such systems, such as the set of evolutions of some quantum state under many trials of a noisy circuit <cit.>. An ensemble of noisy circuits is an approximate description of some more complex, closed system, and more generally ensembles of open systems provide an alternative description of closed systems. In that case, the ensemble is generated by a partition of the closed system into a set of subsystems, with the evolution of each subsystem described by an open-systems equation. The set of sub-system evolutions is constrained by conservation laws and symmetries of the closed, unitary evolution. One might in practice be primarily interested only in the evolution of some special subsystems, for example the computational qubits or the probe systems for a measurement task. For other questions, such as studying thermalization, it may be useful to define a typical or average open evolution together with the deviations from the average.  

In this work, we derive and characterize ensembles of open-system qubit evolutions. We consider cases where either all members of the ensemble are phase-covariant dynamical maps <cit.>, or ensembles whose average map is phase covariant. Prior work on random ensembles of open-system dynamics has largely been restricted to evolution that is time-homogeneous, Markovian, and unital <cit.>. Phase-covariant maps provide a well-characterized way to consider more complex open-system evolution with none of those restrictions <cit.>, making them appealing for phenomenological studies of thermalization, quantum homogenization, dephasing process, quantum metrology and quantum optics <cit.>. For computational questions, they carry the advantages of non-unitality <cit.>. 

The most general phase-covariant qubit channel can be parameterized as 

    Λ_ PC=[        1        0        0        0;        0  λ_1cosθ -λ_1sinθ        0;        0  λ_1sinθ  λ_1cosθ        0;      τ_3        0        0      λ_3 ] ,

where λ_1, λ_3, θ, and τ_3 are all real numbers and complete positivity imposes some restrictions on the parameters, discussed below. The action of this class of channels includes translations and xy-symmetric deformations of the Bloch sphere that commute with rotations about the z-axis. Phase covariant maps provide agnostic noise models, since they describe a combination of pure dephasing with energy absorption and emission. An additional useful property of phase-covariant channels is that they need not have the maximally mixed state as a fixed point. That is, they need not be unital. The state that is invariant under the action of the phase-covariant map, Eq.(<ref>), is ρ_* =  1/2(1+ (τ_3/1-λ_3) Z),
where Z is the Pauli matrix. This state is different from the maximally mixed state (the Gibbs state at infinite temperature) whenever τ_3≠ 0, but it is still diagonal and can be understood as a classical ensemble or Gibb's state with β_* = log[ 2/1-a_*]>0, where a_* = τ_3/1-λ_3. This simple steady state is useful for studying thermalization beyond the Markovian regime <cit.>. 

Evolution under a phase-covariant dynamical map (a smooth, time parameterized sequence of phase-covariant channels) can be generated for a system with free Hamiltonian H_ S by coupling it via an interaction term H_ SE to an environment with free Hamiltonian H_ E in initial state ρ_ S⊗ρ_ E, with the restrictions 

    [H_ S+H_ E,H_ SE]=[H_ E,ρ_ E]=0 .

These relations imply that (i) no energy builds up within interactions and (ii) the initial environment state is a Gibbs state with respect to H_ E, and plays the role of a classical reservoir in the case that the reduced dynamics is Markovian. But, these maps are generically non-Markovian, allowing much richer dynamics <cit.>.

By choosing Hamiltonians for N coupled spins with initial states that obey the conditions above, we can construct ensembles of open-system evolution where each member of the ensemble is a phase-covariant map for a single spin degree of freedom. Except in special cases, the full set of maps is time-consuming to derive for large N. However, we can take a useful limit that effectively averages over the oscillatory behavior associated with the finite system size. We will find that in useful cases this average map provides a stationary channel that depends only on the initial state and the symmetries of the interaction terms in the Hamiltonian. 

To construct these channels from example systems, we first average over the set of individual dynamical maps at fixed time t associated with the evolution of the ith subsystem, Λ_i(t'),  and then take the limit of a long time-average:

    Λ^∞_N≡lim_t→∞1/t∫_0^t[1/N∑_i=1^N
        Λ_i(t')]dt' .

In useful cases, the limit of the long-time and ensemble averaged map will be a time-independent channel. Furthermore, for disordered Hamiltonians, Λ_N(t) (and the map averaged over the N sub-systems at any time) may be phase-covariant even when individual maps are not. This is most directly seen by working with the distribution of map parameters, inherited from the Hamiltonian disorder, over the N sub-systems. Then Λ^∞_N may again be a time-independent, phase-covariant channel. We will demonstrate these points with specific examples below.

In the rest of the paper, we establish a procedure for generating ensembles of phase-covariant maps with a long-time average that is time-independent and with fluctuations about the average that mimic those found in finite-size spin networks. We begin in Section 2 by constructing ensembles of maps for spin networks where every division into sub-system and environment satisfies Eq.(<ref>). That is, every map is phase-covariant and differences in the map come from the inhomogeneity of the initial state. We derive the Λ^∞_N for those ensembles. In Section 3 we provide a perturbative calculation justifying the use of Λ^∞_N as defined in Eq.(<ref>) as a large-N limit for at least some kinds of systems. Next, in Section 4, we consider a case where different members of the ensemble may also have different dynamics, for example, as generated by different values for couplings within the full, closed Hamiltonian. There we demonstrate how moving from a distribution over Hamiltonian parameters to a distribution over map parameters facilitates the construction of additional scenarios whose average dynamics are phase covariant.  We put these pieces together in Section 5, and provide a prescription for generating ensembles of open-system evolution that we conjecture to be useful for studying the ensemble of subsystems of large-N, closed systems. We conclude with discussion in Section 6.


 


§ PHASE-COVARIANT ENSEMBLES FROM XXZ NETWORKS
 



In this section we consider several examples of XXZ spin-networks in external magnetic fields, which we use as benchmark examples of the phase-covariant ensembles. We are interested in extracting features that can be usefully generalized to a construction that can be carried out using only the parameters of the maps in Eq.(<ref>), without reference to an explicit Hamiltonian, although we will see how properties of the Hamiltonian are recognizable in the maps. The phase-covariant restrictions on XXZ models have a simple late-time averaged behavior, described by a single phase-covariant channel (since these channels are closed under convex linear combinations <cit.>). Besides characterizing this channel, we use these models to extract a typical scaling with N of how maps vary around the average, which can be used to inform more general probability distributions over map parameters.  

To begin, the dynamical map associated to the i^ th qubit in the network is defined as

    Λ_i(t)[ρ_i(0)]   = tr_j≠ i[e^-iHtρ_i(0)⊗ρ_ E_i(0)e^iHt]

 where H is the full Hamiltonian describing the network dynamics. Λ_i(t) depends on time, the Hamiltonian parameters, and the initial state of the remaining qubits. To obtain phase-covariant dynamical maps for a given qubit network, we must ensure the thermodynamic constraints from Eq.(<ref>) are satisfied. Therefore, we choose initial global states that are totally uncorrelated and diagonal in the computational basis (as appropriate for XXZ type Hamiltonians)

    ρ(0)=1/2^N⊗_i=1^N(1_i+z_iZ_i)  .

The initial state must be totally uncorrelated so that for all possible choices of focal qubit the reduced dynamics is completely positive. Then, all maps are of the form Eq.(<ref>) and satisfy the following inequalities <cit.>
 
    |λ_3,i|+|τ_3,i|   ≤1
    
    4λ^2_1,i+τ^2_3,i   ≤(1+λ_3,i)^2 .






 §.§ XXZ-networks



In this section, we consider homogeneous XXZ-networks with two different topologies. The first has nearest-neighbor interactions and periodic boundary conditions (a ring), with Hamiltonian 

    H^ ring_ XXZ=h∑_i=1^NZ_i+J_⊥/2∑_i=1^N(X_iX_i+1+Y_iY_i+1)+J_∥/2∑_i=1^NZ_iZ_i+1 .

Here X,Y,Z are the usual Pauli matrices and the qubit indices take integer values mod N. The coupling constants J_∥ and J_⊥ determine the class of the model. For J_⊥=0, H_ XXZ reduces to the (anti)-ferromagnetic Ising model with J_∥<0 (J_∥>0). When J_⊥=J_∥ the model reduces to the XXX model. These models belong to  different universality classes as they flow to distinct fixed points under renormalization. 


Second, we consider completely connected networks, with all possible pairwise interactions:

    H^cc_ XXZ=h∑_i=1^NZ_i+J_⊥/2∑_i,j(X_iX_j+Y_iY_j)+J_∥/2∑_i,jZ_iZ_j .


Technically, the computation of the dynamical maps via Eq.(<ref>) involves exponentiating the Hamiltonian and tracing over most of the system. This is easiest to do if the Hamiltonian is first diagonalized. Both topologies have a number of symmetries that can be exploited to perform this diagonalization. A detailed description can be found in the Appendix, but we summarize the key points here as they provide a helpful organization to understand some of the structure in the maps.  

First, H_ XXZ for either connectivity has a U(1) symmetry generated by the charge operator Q_N=∑_iZ_i. Therefore H_ XXZ is block-diagonal in an excitation basis,

    𝒫_NH_ XXZ𝒫_N^†=H^0_ XXZ⊕ H^1_ XXZ⊕...⊕ H^N_ XXZ ,

with dim H^q_ XXZ= N q and 𝒫_N is a 2^N×2^N permutation matrix taking the computational basis to a chosen excitation basis. There are many choices of excitation basis, as one can always reorder computational basis states within a given q-block to obtain a new basis. 



Next, one can leverage the translation symmetry of the models. For homogeneous couplings and the two topologies considered here, the Hamiltonian is invariant under relabelings like i→ i-1. That is, the translation operator over N spins, T_N, satisfies [T_N,H_ XXZ]=0. Thus the next stage in block-diagonalization for both models is achieved by using a Fourier basis, with blocks labelled by a:

    Γ_N𝒫_N(H_ XXZ)𝒫_N^†Γ_N^†=⊕_q=0^N⊕_a=0^N-1H^q,a_ XXZ .

See the Appendix for additional details and examples. 

However, for N>3, even these symmetries are not enough to achieve exact diagonalization (see, for example, Eq.(<ref>) in the Appendix). The additional degeneracy in the eigenstates of T_N increases with N, and so eigenstates of H_ XXZ are superpositions of these Fourier modes. That is, if eigenstates of the Hamiltonian are labelled by a complete set of quantum numbers {q,a,l}, and the index k runs over the Fourier modes in each block a, then

    |E^a_q;l⟩   =∑_k=k_ min^k_ max C_q^kl|ℱ^a_q;k⟩ .

Whenever k_ max=k_ min, the stationary states are just Fourier modes i.e. 
|E^a_q;k_ min⟩=|ℱ^a_q;k_ min⟩. For small (k_ max-k_ min), the time evolution may be computed analytically. That is the procedure we will follow in the examples below, but the maximum degree of degeneracy amongst Fourier modes grows linearly with N, and it becomes more difficult to perform the diagonalization.  

To make this a bit more explicit, consider the q=1 block for an XXZ-network with N qubits. The block has dimension N and T_N is non-degenerate within it, so diagonalization of H^q=1_ XXZ is achieved through an N-dimensional discrete Fourier transformation. Explicitly the eigenstates are given as discrete plane waves

    |ℱ^a_1;0⟩≡1/√(N)∑_k=0^N-1e^i kθ_aσ^+_k|1...1⟩ ,

where θ_a=2 π a/N and σ^+_k flips the kth spin from 1 to 0 and acts as the identity every where else. The q=N-1 block has essentially identical eigenstates, expect with the replacement

    σ^+_k|1...1⟩→σ^-_k|0...0⟩ .
 
For other values of q, the translation operator becomes more degenerate. In these cases an appropriately sized discrete Fourier transformation must be performed to reach the next step of block-diagonalization. 

In the remainder of this subsection, we go through the computation of the dynamical maps for three qubits in detail. Then, in the rest of the section, we report map parameters and average maps for N>3, leaving details of the process to the Appendix. 





  §.§.§ 3-qubit XXZ-network



We shall now demonstrate how we compute the ensemble of phase-covariant maps and the average map, and characterize the fluctuations, for the example of a 3-qubit network. In this case, the ring and completely connected topologies are equivalent. We begin by diagonalizing H_ XXZ, going first to the canonically ordered excitation basis

    𝒫_3H_ XXZ𝒫_3^†= H^0_ XXZ⊕ H^1_ XXZ⊕ H^2_ XXZ⊕ H^3_ XXZ .
 
This example is especially simple as the only interacting blocks, H^1_ XXZ and H^2_ XXZ, are both diagonalized through discrete Fourier transformation. The other blocks are non-interacting, and consist of the computational basis states |1... 1⟩ (q=0) and |0...0⟩ (q=N).
Using the exact diagonalization of the Hamiltonian, we determine explicit expression for the parameters appearing in the phase-covariant dynamical maps when the initial state is of the form Eq.(<ref>). We find  

    λ_1,i(t,3)=√(α_i^2(t,3)+β_i^2(t,3))
       θ_i(t,3)=2ht+ϕ_i(t,3) 
       λ_3,i(t,3) = 1/9(5+4cos3 J_⊥t)
       τ_3,i(t,3) = 2/9(∑_k ≠ iz_k)(1-cos3 J_⊥t)

where we have made the definitions

    α_i(t,3)   =1/18((1-∏_k≠ iz_k)(7+2cos3J_⊥t)
       +(1+∏_k≠ iz_k)(3cos(2J_∥-2J_⊥)t+6cos(J_∥+2J_⊥)t)) 
    β_i(t,3)   =1/18(∑_k ≠ iz_k)(3sin(2J_∥-2J_⊥)t+6sin(J_∥+2J_⊥)t) 
    ϕ_i(t,3)   =  tan^-1(β_i(t,3)
            /α_i(t,3)) .

Notice that the strength of the external magnetic field, h only enters into the rotation angle θ. 

There are two relevant time scales in these maps, determined by J_∥ and J_⊥. For simplicity, we will illustrate results in Figure 1 for J_∥= J_⊥, and define a single time-scale 

    t_J=2π/J_⊥ .

While a simplifying assumption, we do not expect this reduction to drastically change the primary aspects of our results for the average phase-covariant maps computed below. For example, the late-time average values of both λ_3 and τ_3 are independent of the ratio J_⊥/J_∥. Further, if we avoid certain commensurate values of Hamiltonian parameters, the steady-state values of λ_1 are also independent of the coupling ratio. 


We now turn to the computation of the steady-channel, Λ_3^∞ and the normalized fluctuations. For generic values of h (i.e. not commensurate with J_⊥ or J_∥), Λ_xx(t), Λ_xy(t), Λ_yx(t), and Λ_yy(t), vanish under long-time averaging. Therefore we assume that h is generic, precluding this situation. Choosing to define the time average of λ_1(t) as

    λ_1(t)=√(Λ^2_xx(t)+Λ^2_xy(t))

we then find

    λ_1(t=∞,N=3)=0  ,
       λ_3(t=∞,N=3)=5/9 , 
       τ_3(t=∞,N=3)=4/9𝔢_1(z_i)/3 .

Anticipating results for larger N, we express τ_3^∞ in terms of elementary symmetric polynomials in the initial state 𝔢_k(z_i). These polynomials are defined for any N by the equation

    ∏_i=1^N(z+z_i)=𝔢_0(z_i)z^N+𝔢
        _1(z_i)z^N-1+...+𝔢_N(z_i) .

Examples and related applications of the polynomials can be found in <cit.>. 

Before we move on to the next topic, it is useful to point out a few features of the dynamical map ensembles for the various N we consider. 

Figure 2 emphasizes the non-time-local nature or non-invertibility of our ensembles, following from the existence of zeros in the λ_3(t) plots. The top left panel shows that while not every individual map is time-local for the N=4 hierarchy state, the average dynamics will be. While not plotted here, the Néel serves as a maximal example of locality breaking, as even the average dynamical map is time-non-local.

We may also use these plots to draw a point of distinction between the two topologies of networks. In the completely connected case λ_3(t) and τ_3(t) are always periodic over some multiple of t_J, while in the case of the XXX-ring the dynamics need not be periodic, as can be seen from the lower panel of Figure 2. Of course if N is a perfect square then periodicity occurs, as can be seen for the case of N=4. Equations (<ref>) and (<ref>) in the Appendix makes this point clear, showing that for the XXX-ring √(N)J_⊥ appears in the various frequencies comprising the dynamical map components.  




 §.§ Long-time averaged channels



The individual maps and steady-channels for N>3 require a similar, but longer computation. the two topologies of XXZ-clusters for. We also discuss features of the reduced dynamics found for the computed ensembles, i.e. whether they are time-local for certain initial states and other average properties. Details of the map components as a function of time, and some additional plots of that time-dependence, can be found in the Appendix.




  §.§.§ Complete connectivity
 

Turning our attention to the completely-connected XXZ-network we compute the late-time network averaged quantum channels. For N=4 we find 

    λ^∞_1(N=4)   =0 
    λ^∞_3(N=4)   =7/16𝔢_0(z_i)/1+3/16𝔢_2(z_i)/6
    τ^∞_3(N=4)   =9/16𝔢_1(z_i)/4-3/16𝔢_3(z_i)/4  .

 For N=5 we find

    λ^∞_1(N=5)   =0 
    λ^∞_3(N=5)   =7/15𝔢_0(z_i)/1+16/75𝔢_2(z_i)/10
    τ^∞_3(N=5)   =8/15𝔢_1(z_i)/5-16/75𝔢_3(z_i)/10 .

and for N=6 we find
 
    λ^∞_1(N=6)=0 
       λ^∞_3(N=6)=59/144𝔢_0(z_i)/1+5/12𝔢_2(z_i)/15-5/48𝔢_4(z_i)/15
       τ^∞_3(N=6)=85/144𝔢_1(z_i)/6-5/12𝔢_3(z_i)/20+5/48𝔢_5(z_i)/6 .

The relationship appearing between λ^∞_3 and τ^∞_3 are in part due to constraints that must be satisfied. Namely if all z_i=± 1, the total global state is a fixed point of the unitary evolution therefore

    τ^∞_3|_z_i=± 1±λ^∞_3|_z_i=± 1=± 1  .

 Our results taken together with these constraint leads us to make the following conjectures for the values of any N where for odd values

    λ^∞_1( N = odd)=0 
       λ^∞_3(N = odd)=a_0𝔢_0(z_i)/1+a_2𝔢_2(z_i)/N 2+...+a_N-3𝔢_N-3(z_i)/N N-3
       τ^∞_3(N=  odd)=(1-a_0)𝔢_1(z_i)/N 1-a_2𝔢_3(z_i)/N 3-...-a_N-3𝔢_N-2(z_i)/N N-2

and for even values

    λ^∞_1(N = even)=0 
       λ^∞_3(N = even)=a_0𝔢_0(z_i)/1+a_2𝔢_2(z_i)/N 2+...+a_N-2𝔢_N-2(z_i)/N N-2
       τ^∞_3(N=  even)=(1-a_0)𝔢_1(z_i)/N 1-a_2𝔢_3(z_i)/N 3-...-a_N-2𝔢_N-1(z_i)/N N-1 .

Another property of these long time averages, and in fact of the dynamical map components in general, is that λ_3 is an even function of {z_i} while τ_3 is an odd function. This is a consequence of the symmetry generated by  

    𝒫=∏_i=1^NZ_i.

which is a symmetry for both models considered, therefore λ_3(t) and τ_3(t) have definite values of 𝒫 in the ring topology as well. 


 To conclude the discussion, notice that the average map depends only on the initial state and symmetry of the Hamiltonian. We note that in this case the long time averages are invariant under any permutation of the initial z_i. Comparing the fully connected case to the ring, shown next, illustrates this dependence. In addition, notice that the invariant state of the average map is not the average of the initial states, but any explicit dependence on coupling strength is gone. As expected from the individual maps at any time, there is also no dependence on magnetization. In the case of non-phase covariant maps, there will be terms depending on h, and, just as in the case of λ_1, special values of the coupling lead to different long time average values.   




  §.§.§ Ring connectivity


 Computing the late-time and ring-averaged quantities we find for N=4

    λ^∞_1(N=4)   =0 
    λ^∞_3(N=4)   =7/16+3/16𝔢_2(z_i)-4z_1z_3-4z_2z_4/6
    τ^∞_3(N=4)   =9/16𝔢_1(z_i)/4+1/16𝔢_3(z_i)/4 ,

and for N=5

    λ^∞_1(N=5)   =0 
    λ^∞_3(N=5)   =71/225+2/45𝔢_2(z_i)/10
    τ^∞_3(N=5)   =154/225𝔢_1(z_i)/5-2/45𝔢_3(z_i)/10 .

Comparing to the complete connectivity, Eq.(<ref>) and Eq.(<ref>), demonstrates how the interaction structure appears in the maps. For example taking N=4, the long-time average is only invariant under cyclic permutations of the initial state while in the fully connected case the invariance was with respect to the full permutation group. However, for N=5 the full permutation symmetry is recovered in the infinity time limit. But we do not expect this to be a generic feature for larger N, and in principle only cyclic polynomials of the initial state should appear. As before, the invariant state of the average map does not correspond to the average initial state.



 §.§ Fluctuations about the average channel

 We now consider physical data to constrain the variance we expect in the distributions at late-times. We do this by considering the maximum distance from the finite-time average to the steady value. We define the fluctuations at finite-time about these average quantities, estimated by the normalized fluctuations

    Δλ_i   ≡λ_i(t)-λ_i^∞/λ_i^∞
    Δτ_3   ≡τ_3(t)-τ_3^∞/τ_3^∞  .

These quantities are used to calibrate the second moments of the phase-covariant measures developed in Section 5. The normalized fluctuations are generally initial-state dependent, but to determine an approximate size of the fluctuations, we consider a class of initial states with essentially initial state independent fluctuations. To see that such a class of states exists, consider the completely connected network and assume that each qubit satisfies 0≤ z_i≤𝒵<1. Since

    𝔢_k(z_i)/N k≤𝒵^k ,

then for small enough 𝒵 only the first terms in λ_3(t) and τ_3(t) are relevant. In such a class the fluctuations about the long time average will be initial state independent. 

Figures 3 and 4 show the fluctuations about the long time averages for the completely and ring connected networks, using only the lowest order terms in the initial state from λ_3(t) and τ_3(t). We confirm in Section 5 that these approximations provide a reasonable estimate of the fluctuation size, even taking a generic initial state. 

These plots show that the late-time fluctuations go as 

    Δζ≈C_ζt_J/Nt ,

where t_J=2π/J_⊥ is the dynamical timescale defined in Eq.(<ref>) and C_ζ is an order one constant. The completely connected network, for a fixed N, approaches the steady channel more quickly than in the case of the ring. Comparing the N=6 completely connected network and the N=5 ring network, we see that the fluctuations persist at least four times longer in the case of the ring. Note that the scale of the vertical axis of all plots are identical. 








§ PHYSICAL RELEVANCE OF THE LONG-TIME AVERAGED CHANNEL





 Consider a qubit network that consists of N_ CL clusters, each containing N qubits, where we assume  N_ CL≫ N. For simplicity, each cluster is assumed to have the same connectivity. The spin-averaged dynamical map for the Ith cluster is 
 
    ⟨Λ_ I(t)⟩   =1/N∑_i=1^N
        Λ_ I_i(t)  .

 Suppose at some very early time there are no interactions and the state of the network is entirely uncorrelated. Interactions within each individual cluster are turned on via a quenching mechanism, but that not all clusters are quenched at the same time. Let the time-scale for over which the cluster quenching is (uniformly) staggered be longer than that set by the interaction time-scale within each cluster. The time-dependent Hamiltonians considered (and initial states) are such that the conditions for phase covariance are satisfied for each cluster. For example, this scenario can be modeled by a Hamiltonian of the form
 
    H(t)= ∑_ I=1^N_ CL( h_ I∑_i=1^N Z_ I_i+Jθ(t-t_ I)∑_i,j(X_ I_iX_ I_j+Y_ I_iY_ I_j+Z_ I_iZ_ I_j))

 where θ(t) is the Heaviside step function, and t_ I are the various quench times associated with the clusters. If the timing of the quench varies over a time-scale longer than that set by J, and the large number of clusters ensures that all oscillation scales of the single-cluster Hamiltonian are smoothly sampled, then the clusters will be at varying points in their evolution and the average map at late times will be given by
 
    lim_N_ CL→∞1/N_ CL∑_ I⟨Λ_ I(t)⟩  =1/t∫_0^tdt'⟨Λ_ K(t')⟩=Λ_N(t) .

In fact, as long as there are many clusters with any given initial state, which experience the quench at offset times, one can imagine a further average over steady channels with different initial states to match the full complexity of the initial system of size ≫ N. That is, one could average over the initial states appearing in the steady channels.

Under such a quenching scenario, even for relatively small clusters, we obtain average dynamical maps that approach a steady-channel, denoted as Λ_N^∞, in the limit t →∞. We are then able to construct distributions of channels that have the steady channel as their average. By studying the approach to the steady channel, we able to constrain the second moment of the distributions. 






§ PHASE-COVARIANCE ON AVERAGE FROM DISORDERED HAMILTONIANS
 




While in the staggered quench example above we imposed that each cluster had the same Hamiltonian and an initial state from the phase-covariant set, the regime of applicability of the phase-covariant steady channel is actually larger. In this section we show that constrained Hamiltonian disorder, or noise, can be rather naturally consistent with phase-covariant, average reduced dynamics. We show this by considering ensembles of qubit pairs, which serve as a simple model of the paradigm set forth at the end of the previous section.





 §.§ An ensemble with non-phase-covariant elements



 To begin we start with a physical setting where distributions of maps that are not necessarily phase-covariant appear. We consider a Hamiltonian that governs the time evolution of N_ CL pairs of spins, where each pair evolves independently, and conserves the total excitation number Q=∑ Z_i. This can be compactly written as

    H   =∑_i=1^2N_ CL h_i Z_i+∑_i=1^2N_ CLJ_i(i+1)/2[1+(-1)^i+1/2][X_iX_i+1+Y_iY_i+1] 
       =∑_i=1^2N_ CL1+(-1)^i+1/2H^i,i+1_ XX=H^1,2_ XX+H^3,4_ XX+...+H^2N-1,2N_ CL_ XX .

The 2N_ CL qubits are taken to be in a completely uncorrelated initial state. For each pair, we allow this example to break both phase-covariant conditions from Eq.(<ref>). First, allowing h_i≠ h_j means that the interaction term will not commmute with the free Hamiltonian. In addition, we do not impose any restriction on the initial state. 

As the Hamiltonian is block diagonal in the qubit pairs, we can determine both the 1-spin and 2-spin dynamical maps obtained for a single pair. The 2-qubit maps describe unitary evolution, and this example may be viewed as an generalization of the noisy qubit models proposed in <cit.> to a larger system size. However, at the single-spin level the evolution is non-unitary and richer than what can be found with noisy Hamiltonians alone.

To write the single-spin dynamical map, it is convenient to decompose the unitary evolution in terms of 

    𝒰^ij_lk≡1/4 tr[(σ^i_ 1⊗σ^j_ 2)e^-i(H^1,2_ XX)t(σ^l_ 1⊗σ^k_ 2)e^i(H^1,2_ XX)t]  .

The Hamiltonian parameters enter in the unitary evolution in the following combinations:

    h_12=h_1+h_2/2 ,     ω_12= sgn(Δ_12)√(Δ_12^2+J_12^2)
       Δ_12=h_1-h_2 ,     tanϕ_12=J_12/Δ_12 .

Further details can be found in the Appendix. Then the evolution of the first spin in the first pair, for example, is given by

    ρ_1(t) = tr_2[ e^-i (H^1,2_ XX) tρ_1(0)⊗ρ_2(0)  e^i (H^1,2_ XX) t] ≡Λ_1(t,0) ρ_ 1(0)  .

Using the Bloch parameterization for the second spin (which provides the initial environment state for the first spin) as ρ_2=1/2(1+x_2X_2+y_2Y_2+z_2Z_2), and defining partial components

    Λ_1^ijk=𝒰^i0_jk ,
 
the dynamical map in the Pauli basis is

    Λ_1(t,0)=
       [                        1                        0                        0                        0;                        0                  Λ_1^xx0                  Λ_1^xy0 (x_2Λ_1^xzx +y_2Λ_1^xzy);                        0                  Λ_1^yx0                  Λ_1^yy0 (x_2Λ_1^yzx +y_2Λ_1^yzy);               z_2Λ_1^z0z (x_2Λ_1^zxx +y_2Λ_1^zxy) (x_2Λ_1^zyx +y_2Λ_1^zyy)                  Λ_1^zz0 ]

and is non-unital so long as z_ 2(0)≠0. The dynamical map Λ_2(t,0) for the second spin can be obtained directly from Eq.(<ref>). From the Hamiltonian it is easy to see that we need only exchange the values of h_i, so Λ_2(t,0) is directly obtained by exchanging the initial states and sending Δ_12→ -Δ_12. In terms of map components, this is accomplished by

    Λ_2^ijk(t,Δ_12)=𝒰^0i_kj(t,Δ_12)=𝒰^i0_jk(t, -Δ_12)=Λ_1^ijk(t, -Δ_12) .

The dynamical maps associated to the second qubit pair are obtained through the replacement 12 → 34 in the initial states, h_i and J_i,i+1. In this way we obtain an ensemble of non-unital 1-qubit dynamical maps and unitary 2-qubit dynamical maps. 

Given the collection of disjoint qubit pairs, and the associated 2-qubit unitary channels, it is possible to construct their ensemble averages. Both averages will be non-unitary, since convex combinations of unitary channels may be non-unitary. This provides a useful means to efficiently simulate non-unitary evolution <cit.>. If each pair is assumed to be in the same uncorrelated 2-qubit state, the following non-unitary ensemble averaged 2-qubit dynamical map may be defined 

    𝒰(t)[ρ_1(0)⊗ρ_2(0)]=1/N_ CL∑_ clusters(e^-i (H^i,i+1_ XX) tρ_1(0)⊗ρ_2(0) e^i (H^i,i+1_ XX) t)  .


At this point we have a physical representation of the sort of noisy Hamiltonians that have appeared in numerous works. Now conceptually we may view the the previous ensemble as an instantiation of some noise distribution over Hamiltonian parameters. Therefore in the limit of large N_ CL the network average is equivalent to 

    𝒰(t)[ρ_1(0)⊗ρ_2(0)]=⟨ e^-i H_ XX tρ_1(0)⊗ρ_2(0) e^i H_ XXt⟩_ disorder

where we have dropped the superscripts on H_ XX to emphasis the conceptual difference with the network of clusters. 

 Now instead of considering a disorder distribution directly over the parameters in the Hamiltonian, we consider the distributions directly over its eigenvalues and eigenvectors. That is, we take the fundamental parameters to be those on the left hand side of Eq.(<ref>), e.g. h (from h_12), ω, and ϕ, rather than (h_1,h_2,J_⊥). This approach allows for a closer connection to the dynamical map components. For example, the following component

    Λ^xzy_1∝sinϕ ,

which is compactly expressed in terms of the mixing angle ϕ from H_ XX eigendecomposition. Notice that such a component vanishes if the noise distribution is an even function of ϕ, where tanϕ=J/Δ. The other dynamical map components that break phase-covariance also vanish under the same assumption. We understand this as the phase-covariant constraints in Eq.(<ref>) being satisfied on average. 


We can apply a similar procedure to other non-phase-covariant maps. In particular, for more general non-unital quantum channels with non-zero shift in the x and y direction (τ_1≠0 and τ_2 ≠0), the invariant state can have coherences and therefore cannot be represented as a Gibbs state in the computational basis. However, ensembles of such general non-unital maps can be generated such that the shift in the x and y directions have mean zero, resulting in a Gibbs-like state for the ensemble-averaged map.



 §.§ Phase-covariance via disorder averaging

We begin by considering independent Gaussian distributions for h and ω,

    p_h(h)= 1/√(2πσ_h^2)e^-(h-B)^2/2σ_h^2   ,      p_ω(ω)= 1/√(2πσ_ω^2)e^-(ω-Ω)^2/2σ_ω^2 ,

and two different examples of probability distributions over ϕ,

    p_ G(ϕ)= 1/√(2πσ_ϕ^2)e^-ϕ^2/2σ_ϕ^2   ,      
       p_ T(ϕ)= tanh[a_ϕϕ]^2/(π-2/a_ϕtanh[a_ϕπ/2]) ,

where for p_ T(ϕ) the range is restricted to be -π/2≤ϕ≤π/2. Both distributions have ⟨ϕ⟩=0 to ensure that the 1-qubit dynamical maps are phase-covariant. We consider multiple ϕ-distributions as the average values of τ_3 and λ_3 are strongly tied to this choice. Certain ϕ-distributions restrict the range of allowed average values of map parameters, and we see from 
 Figure 5 that the Gaussian distribution cannot generate large enough τ_3 to be compatible with all examples from Section 2. 

Using these distributions, the partial components of the effective focal dynamical map are 

    Λ_1^xx0=Λ_1^yy0=(cosΩ tcos2Bt-e^-π^2/2φ^2sinΩ tsin2Bt)e^-(σ_ω^2+4σ_h^2)t^2/2
       Λ_1^yx0=-Λ_1^xy0= (cosΩ tsin2Bt+e^-π^2/2φ^2sinΩ tcos2Bt)e^-(σ_ω^2+4σ_h^2)t^2/2
       Λ_1^z0z=1-Λ_1^zz0=1/4(1-e^-2π^2/φ^2)(1-cos(2Ω t)e^-2σ_ω^2t^2) .

Figure 5 shows the non-unitary phase-covariant parameters of the map Λ_1(t,0). For plotting purposes we have defined the dimensionless quantity χ=Ω/2πt. 


We have introduced the truncated distribution, p_ T(ϕ), to demonstrate that disorder distributions can capture the larger values of τ_3 found in Section 2. The truncated hyperbolic tangent distribution has a maximum allowed value of τ_3 given by  

    lim_a_ϕ→ 01/π-2/a_ϕtanh[a_ϕπ/2]∫_-π/2^π/2sin^2[ϕ]tanh^2[a_ϕϕ] dϕ = 6+π^2/2π^2≈ 0.8  .







§ PHASE-COVARIANT DYNAMICS FROM CHANNEL ENSEMBLES





In the final section of this work we construct measures over phase-covariant channels that capture key features of the long-time averaged channels derived from XXZ networks. However, the time-evolving ensemble can be considered on its own, characterized by its thermodynamic properties without reference to a particular Hamiltonian. From the Gaussian examples considered in Section 4, we are motivated to considered a phase-covariant measure consisting of truncated Gaussian distributions over the non-unitary parameters. We use the physical constraints from Section 2 to fix the mean (μ) and scale (σ) parameters. 




 §.§ Phase-covariant measures



The key point in defining a measure over channels is determining the domain of integration. This depends on the constraints imposed by complete positivity, which are especially simple in the case of phase-covariant qubit channels. Then, a phase-covariant measure may be defined by supplying a joint probability distribution p_N(λ_1,τ_3,λ_3;t) which satisfies the normalization condition

    1=   ∫_Γ_ PC p_N(λ_1,τ_3,λ_3;t) dλ_1 dτ_3 dλ_3
     =   ∫_0^1∫_-(1-λ_3)^1-λ_3∫_-1/2√((1+λ_3)^2-τ_3^2)^1/2√((1+λ_3)^2-τ_3^2)p_N(λ_1,τ_3,λ_3;t)dλ_1 dτ_3 dλ_3
    +   ∫_-1^0∫_-(1+λ_3)^1+λ_3∫_-1/2√((1+λ_3)^2-τ_3^2)^1/2√((1+λ_3)^2-τ_3^2)p_N(λ_1,τ_3,λ_3;t)dλ_1 dτ_3 dλ_3 ,

where the order of integration follows from the form of the constraints in Eq.(<ref>). For more general quantum systems it may not be possible to determine the boundary of integration, but restricting the reduced dynamics to be G-covariant potentially makes the problem tractable for larger quantum channel. 

The above measure allows us to define the average dynamical map,    

    ⟨Λ_N⟩(t)=∫_Γ_ PC p_N(λ_1,τ_3,λ_3;t) Λ(λ_1,τ_3,λ_3)dλ_1 dτ_3 dλ_3 .

 We choose to introduce time-dependence into the probability distribution, instead of directly into the dynamical map, as the time-dependence of the dynamical map becomes quite complicated for large clusters. 

This form of the map ensemble can be usefully compared to prior explorations of random open system dynamics, and to compare the phase-covariant class to the larger class of channels. 
  
  We compare the distribution of eigenvalues for a non-unital phase-covariant ensemble to those of the more general class of dynamical maps, treated for example in <cit.>. To do so, it is helpful to diagonalize the phase-covariant channel and define
 
    D_Λ_ PC=
         [          1          0          0          0;          0  λ_1 e^i θ          0          0;          0          0 λ_1 e^-i θ          0;          0          0          0        λ_3;            ] ,

where μ_±=λ_1 e^± i θ. Note that τ_3 does not appear in the eigenvalues but only in the eigenvectors. The dynamical map parameters (restricted by complete positivity) are uniformly sampled and plotted on the left panels of Figures 6 and 7.

To make a comparison with the more general class of dynamical maps in <cit.>, we consider a minimal breaking of the phase-covariance through the introduction of a new parameter λ_2,


    Λ_ PC[        1        0        0        0;        0  λ_1cosθ -λ_2sinθ        0;        0  λ_1sinθ  λ_2cosθ        0;      τ_3        0        0      λ_3;          ] .

The eigenvalues of the non-phase-covariant map are

    μ_± = (λ_1+λ_2)cosθ±√((λ_1+λ_2)^2cos^2θ-4λ_1λ_2)/2.

The dynamical map parameters are uniformly sampled (again restricting to complete positivity) and the eigenvalues of the non-phase-covariant ensemble are plotted on the right panels of Figures 6 and 7. This simple breaking of phase covariance introduces additional structure seen in Figure 6 (right), several of which match those in Figure 1 of the slightly more general treatment of <cit.>. For example, there is agreement in the empty parts of the plane, and the apparent `repulsion' of the points away from the real axis. Note also that a significant number of μ_± eigenvalues become real when phase covariance is broken.
 

 
Finally, we end by defining the uniform distribution over phase-covariant channels

    vol(Γ_ PC)= vol(Γ^+_ PC)+vol(Γ^-_ PC)=(16/9-π/6)+π/6=16/9

where we have separated the integral into two regions based on the sign of λ_3. The positive region has a larger volume than that of the negative region, which is confirmed in Figure 7 (left). Therefore when sampling the phase-covariant channels using a uniform distribution, we expect there to be a bias towards positive λ_3. Our results are consistent with, although scaled differently from, those found in <cit.>.



 §.§ Example distribution: Truncated-Gaussians

 There are of course many other non-uniform distributions to use to define a phase-covariant measure, and making an appropriate choice allows one to more easily capture the behavior of the long time average channels determined in Section 2. Although a time-dependent uniform distribution could be used, this is a more difficult to define from the data gathered from the average and late-time fluctuations. Ostensibly, the results of Section 4 indicate that using truncated Gaussian distributions could be a fruitful approach. We consider phase-covariant measures defined in this way for the remainder of this section. 


A truncated Gaussian distribution over a single random variable (X) is given by 

    p_ TG(X;η⃗)=e^-(X-μ)^2/2σ^2/𝒩(η⃗)

where η⃗ stands in for the parameters of the distribution and 𝒩 is the normalization factor

    𝒩(η⃗)=∫_a^be^-(X-μ)^2/2σ^2dX  .

There are four total parameters associated to a two-sided truncated Gaussian, the mean parameter (μ) , the scale parameter (σ), and the end points of the interval that define the truncation. For our purposes the end points of each distribution are set by the complete-positivity constraints. For example the truncated Gaussian associated with λ_3 may take values between -1 and 1, while the other distributions have intervals depending on the values of λ_3. 

The remaining parameters we set through the use of the physical constraints found in the previous section. For example the scale parameter (assumed real) takes the form

    σ=C/N t_ ref/t ,

where the constant C is O(1). We define the mean of each parameter by

    μ_1(N)=λ^∞_1(N)=0
       μ_λ(N)=λ^∞_3(N) 
       μ_τ(N)=τ^∞_3(N) ,

where we have seen that μ_λ≈ a_0 and μ_τ≈ (1-a_0)z are approximations that are valid for a non-empty set of initial states. Notice that these parameters do not correspond to the mean and variance of the truncated Gaussians, and in general they will be larger than these values. 

In order to construct instances of these distributions and how they change over time we use the following procedure. First, we choose a sampling time, set by the time-scale t_ ref=t_J=2π/J_⊥. We then sample the distribution p_ TG(λ_3,η⃗_3) at each time step, where the scale parameter is changed according to Eq.(<ref>). In this way we obtain an ordered list of values that λ_3 during each time-step. From these values of λ_3, we may construct the distribution p_ TG(τ_3,η⃗_τ), as recall the end-point parameters for a given time-step depend explicitly on obtained values of λ_3. The procedure then carries on in the same manner as for λ_3, a successive list of values for τ_3 are generated for each time-step. Finally, using the outcomes of the previous steps a list of values for λ_1 may be generated from the distribution p_ TG (λ_1,η⃗_1). Figure 8 illlustrates the results of this procedure. One could define a smaller sampling time, as it is clear that the sampling procedure does not capture the smaller scale fluctuations. However, these fluctuations do not carry additional physical information about the convergence toward the steady channel, for single-qubit dynamics, at least. 

Care must be taken into the above procedure, as the endpoints for the distribution over τ_3 depends explicitly on the sign of λ_3. One way around the potential issues is to simply take the τ_3-interval to be set by the sign of the late-time limit of λ_3. This guarantees that at late-enough times, no positivity breaking occurs. Of course one may just consider the sign of λ_3 in the above procedure and take care to use the correct interval for τ_3 at a given time-step.







§ CONCLUSIONS


In this work we have defined phase-covariant ensembles in two different ways. The first approach was to construct ensembles directly from a time-independent Hamiltonian, and consider the individual reduced dynamics for all possible subsystems of the smallest size. We constructed examples where each dynamical map in the ensemble was phase-covariant, and examples where only the network averaged reduced dynamics was phase-covariant. From these examples, we are able to conjecture the form of the map parameters in terms of the initial state and connectivity, as shown in Eq.(<ref>) and Eq.(<ref>) for completely connected maps and Eq.(<ref>) and Eq.(<ref>) for one-dimensional rings with only nearest neighbor interactions. While we are unable to conjecture the general form of the long-time average for the ring connectivity, we do see that in general only the translation symmetry should be imposed. 

The resulting time-independent, phase-covariant channels generate a notion of thermalization associated with each Hamiltonian and initial state, since under many repeated applications of the steady channel, every single-spin state will be driven toward the Gibbs state associated with β_* = log[ 2/1-a^∞_*], where a^∞_* = τ^∞_3/1-λ^∞_3. The repeated application of such maps has also found a role in the simulation of open quantum systems <cit.>, although so far restricted to systems described by stationary Lindblads. Our results may open up a similar route to efficient simulation of a broader class of open systems. This is especially interesting as quantum simulators may have an exponential advantage for simulating open-system evolution <cit.>. The procedure in Section 5 could be used in this direction.

The phase-covariant dynamics of 1-qubit subsystems do not seem to differ based on the integrable vs. non-integrable nature of the spin-chain Hamiltonian considered. In principle, such a characterization should be crucially important in determining the nature of a random distribution over channels, similar to how it constrains the form of random unitary dynamics in various large systems. It will be interesting to expand this study to both larger subsystems and more explicit breaking of phase-covariance.








§ ACKNOWLEDGEMENTS

This work was supported by NSF PHY-2310662.
















§ APPENDIX A


A.equation

In this appendix we lay out the symmetries that are useful in computing exact dynamical maps. The first of these symmetries is the charge, or equivalently the excitation number, operator

    Q_N=∑_iZ_i .
 
Throughout the article the excitation number is labeled as q, but note this is not the eigenvalue of Q_N. Each block with a given value of q has dimension m=N q, simply the number of ways to pick q qubits to be excited. 
  
  
  Hamiltonians satisfying [H,Q_N]=0 will be block diagonal in any excitation eigenbasis. Since Q_N is already diagonal in the computational basis, block diagonalizing the Hamiltonian at most requires a permutation matrix that reorders computational basis states by excitation number. For this article we chose 𝒫_N so that the computational states within the same q-blocks are ordered based on their binary representation. For example when N=3, in the q=1 block the computational state |001⟩ comes before |100⟩. The next symmetry available to the models we consider is the translation operator (T_N). The translation action on computational basis states is

    T_N|a_1a_2...a_N⟩=|a_Na_1...a_N-1⟩

and the adjoint action on Pauli matrices is

    T_Nσ^α_iT_N^†=σ^α_i+1 .

where a tensor product of the Pauli with the identity operator on all other spins is implied. Since [Q_N,T_N]=0 the translation operator is block diagonal in any excitation-ordered computational basis. Further noting that T_N^N=1, it is not hard to see that its eigen-spectrum consists entirely of N^ th roots of unity.  Combining these results we may label the eigenstates of  T_N as

    T_N|ℱ^a_q;k⟩=e^i2π a/N|ℱ^a_q;k⟩
       Q_N|ℱ^a_q;k⟩=(2q-N)|ℱ^a_q;k⟩

where q is the excitation number, and k is a label that removes any remaining degeneracy of the eigenstates. In a given q-block, Fourier modes corresponding to subgroups of the cyclic group ℤ_N may appear, which happens when N is a composite integer. In the remainder of this appendix we give constructions of the Fourier modes for N=3 and N=4. 


 §.§ 3-qubits


    |ℱ^0_0;0⟩=|111⟩   |ℱ^0_3;0⟩=|000⟩


    |ℱ^0_1;0⟩=1/√(3)(|011⟩+|101⟩+|110⟩)
       |ℱ^1_1;0⟩=1/√(3)(|011⟩+e^2π i/3|101⟩+e^4π i/3|110⟩) 
       |ℱ^2_1;0⟩=1/√(3)(|011⟩+e^4π i/3|101⟩+ e^2π i/3|110⟩)


    |ℱ^0_2;0⟩=1/√(3)(|001⟩+|010⟩+|100⟩)
       |ℱ^1_2;0⟩=1/√(3)(|001⟩+e^2π i/3|010⟩+e^4π i/3|100⟩) 
       |ℱ^2_2;0⟩=1/√(3)(|001⟩+e^4π i/3|010⟩+ e^2π i/3|100⟩)



 §.§ 4-qubits


    |ℱ^0_0;0⟩=|1111⟩   |ℱ^0_4;0⟩ =|0000⟩


    |ℱ^0_1;0⟩=1/√(4)(|0111⟩+|1011⟩+|1101⟩+|1110⟩)
       |ℱ^1_1;0⟩=1/√(4)(|0111⟩+i|1011⟩-|1101⟩-i|1110⟩) 
       |ℱ^2_1;0⟩=1/√(4)(|0111⟩-|1011⟩+ |1101⟩-|1110⟩)
       |ℱ^3_1;0⟩=1/√(4)(|0111⟩-i|1011⟩-|1101⟩+i|1110⟩)


    |ℱ^0_3;0⟩=1/√(4)(|0001⟩+|0010⟩+|0100⟩+|1000⟩)
       |ℱ^1_3;0⟩=1/√(4)(|0001⟩+i|0010⟩-|0100⟩-i|1000⟩) 
       |ℱ^2_3;0⟩=1/√(4)(|0001⟩-|0010⟩+ |0100⟩-|1000⟩)
       |ℱ^3_3;0⟩=1/√(4)(|0001⟩-i|0010⟩-|0100⟩+i|1000⟩)


    |ℱ^0_2;0⟩ =1/√(4)(|0011⟩+|0110⟩+|1001⟩+|1100⟩) 
       |ℱ^1_2;0⟩= 1/√(4)(|0011⟩+i|0110⟩-|1001⟩-i|1100⟩)
       |ℱ^2_2;0⟩=1/√(4)(|0011⟩-|0110⟩+|1001⟩-|1100⟩) 
       |ℱ^3_2;0⟩= 1/√(4)(|0011⟩-i|0110⟩-|1001⟩+i|1100⟩) 
       |ℱ^0_2;1⟩ =1/√(2)(|0101⟩+|1010⟩)
       |ℱ^2_2;1⟩=1/√(2)(|0101⟩-|1010⟩)



§ APPENDIX B

B.equation


In this Appendix we present the exact formulae for the dynamical maps considered in this work. We have the following form of phase-covariant maps

    Λ(t)=[       1       0       0       0;       0 Λ_xx(t) Λ_xy(t)       0;       0 Λ_yx(t) Λ_yy(t)       0; Λ_z0(t)       0       0 Λ_zz(t) ]

where for the Hamiltonians considered, the various components all may be decomposed as 

    Λ_xx(t)=α(t) cos2ht-β(t) sin2ht=λ_1(t)cos(2ht+ϕ(t))
       Λ_yx(t)=β(t)cos2ht+α(t) sin2ht=λ_1(t)sin(2ht+ϕ(t))
       Λ_zz(t)=λ_3(t) 
       Λ_z0(t)=τ_3(t) ,

with the remaining components determined by the relations Λ_yy=Λ_xx and Λ_xy=-Λ_yx. 
To perform the minimum amount of computations, we only provide explicit formulas for τ_3, λ_3, α, and β for each model. 



 §.§ Completely-connected XXZ-network

We present the exact dynamical map components for the completely-connected and homogeneous XXZ-network. We present all the non-unitary parameters defined in Eq.(<ref>). 


  §.§.§ 4-qubits

For 4-qubits the Hamiltonian has the block diagonal structure

    𝒫_4H_ XXZ𝒫_4^†= H^0_ XXZ⊕ H^1_ XXZ⊕ H^2_ XXZ⊕ H^3_ XXZ⊕ H^4_ XXZ

The only non-trivial diagonalization occurs in the q=2 block. Restricting to this sector, we now invoke the translation symmetry of the model

    Γ_4 H^2_ XXZΓ_4^†=H^2,0_ XXZ⊕ H^2,1_ XXZ⊕ H^2,2_ XXZ⊕ H^2,3_ XXZ

where Γ_4 is a unitary transformation that diagonalizes T_4. The a=1,3 blocks are singlets, while a=0,2 are doublets (see Eq.(<ref>)). Only for a=0 do the Fourier modes mix, although the mixing angle is parameter independent. But obviously the eigenvalues depend on the Hamiltonian paramters. We find the non-unitary parameters of the dynamical map to be 

    λ_3,i(t,4) =   (7/16 +1/4cos2J_⊥t +5/16cos4J_⊥t) 
     +   (1/16-1/12cos2J_⊥t+1/48cos4J_⊥t) ∑_k<l 
     k,l≠ iz_kz_l


    τ_3,i(t,4) =   (3/16-1/12cos2J_⊥t-5/48cos4J_⊥t)∑_k≠ iz_k
    -   (3/16-1/4cos2J_⊥t+1/16cos4J_⊥t)∏_k≠ iz_k



    α_i(t,4)=   (3/16cos(3J_⊥ +J_∥)t +1/16cos3(J_⊥ -J_∥)t
     +   3/32cos(J_⊥ +3J_∥)t +1/4cos(J_⊥+J_∥)t
    +   3/8cos(J_⊥ -J_∥)t +1/32cos(5J_⊥ -J_∥)t)
     +   (3/16cos(3J_⊥ +J_∥)t+1/16cos3(J_⊥ -J_∥)t
     -   1/32cos(J_⊥ +3J_∥)t -1/12cos(J_⊥ +J_∥)t
     -   1/8cos(J_⊥ -J_∥)t -1/96cos(5J_⊥ -J_∥)t)∑_k<l 
     k,l≠ iz_kz_l



    β_i(t,4)=   (3/16sin(3J_⊥ +J_∥)t -1/16sin3(J_⊥ -J_∥)t
     +   1/32sin(J_⊥ +3J_∥)t +1/12sin(J_⊥ +J_∥)t
    -   1/8sin(J_⊥ -J_∥)t -1/96sin(5J_⊥ -J_∥)t)∑_k≠ iz_k
    +   (3/16sin(3J_⊥ +J_∥)t -1/16sin3(J_⊥ -J_∥)t
     -   3/32sin(J_⊥ +3J_∥)t -1/4sin(J_⊥ +J_∥)t
    +   3/8sin(J_⊥ -J_∥)t +1/32sin(5J_⊥ -J_∥)t)∏_k≠ iz_k



  §.§.§ 5-qubits

As the exact diagonalization is more complicated for the 5-qubit network, we shall go through the steps more explicitly than the 4-qubit case. Again we may change to the canonical excitation basis

    𝒫_5H_ XXZ𝒫_5^†= H^0_ XXZ⊕ H^1_ XXZ⊕ H^2_ XXZ⊕ H^3_ XXZ⊕ H^4_ XXZ⊕ H^5_ XXZ .

As with the ring topology, we further go to the basis of translation eigenstates. Again the only excitation number blocks that are not diagonalized trivially at this stage are

    Γ_5H^2_ XXZ⊕ H^3_ XXZΓ^†_5=⊕ _a=0^4 H^2,a_ XXZ⊕ H^3,a_ XXZ .

And as for the ring-topology, we obtain a set of ten 2×2 matrices that require further computation to diagonalize. The q=2 blocks are found to be 

    H^2,0_ XXZ   =
    [ -h-J_∥+2J_⊥        4J_⊥;            ;        4J_⊥ -h-J_∥+2J_⊥ ]
    
    H^2,1_ XXZ=   [ -h-J_∥+(√(5)-1/2)J_⊥       -J_⊥e^-6π i /5;                     ;        -J_⊥e^6π i /5 -h-J_∥-(√(5)+1/2)J_⊥ ]
    
          H^2,2_ XXZ=   [ -h-J_∥-(√(5)+1/2)J_⊥        -J_⊥e^-2π i/5;                     ;         -J_⊥e^2π i/5 -h-J_∥+(√(5)-1/2)J_⊥ ]

With the remaining blocks determined by the relations H^2,2_ XXZ=(H^2,3_ XXZ)^T and H^2,1_ XXZ=(H^2,4_ XXZ)^T. 
The matrix representations for the q=3 block are obtained through the replacement h→ -h, affecting only the eigenvalues and not the form of the associated eigenvectors. Again to due to large degree of symmetry of the model, the eigenvectors are parameter independent. The resulting map parameter functions are:

    λ_3,i(t,5) =   (7/15 +1/3cos3J_⊥t +1/5cos5J_⊥t) 
     +   (8/225-1/18cos3J_⊥t+1/50cos5J_⊥t)∑_k<l 
     k,l≠ iz_kz_l


    τ_3,i(t,5) =   (2/15-1/12cos3J_⊥t-1/20cos5J_⊥t)∑_k≠ iz_k
    -   (4/75-1/12cos3J_⊥t+3/100cos5J_⊥t)∑_k<l<j 
     k,l,j≠ iz_kz_lz_j


    α_i(t,5)=   (157/600+1/12cos3J_⊥t+3/100cos5J_⊥t+3/50cos(3J_⊥+2J_∥)t
    +   1/40cos4(J_⊥-J_∥)t+1/100cos(7J_⊥-2J_∥)t
    +   9/50cos2(J_⊥-J_∥)t+1/4cos(J_⊥+2J_∥)t
     +   1/10cos(J_⊥+4J_∥)t) 
     -   (157/1800+1/36cos3J_⊥t+1/100cos5J_⊥t-1/40cos4(J_⊥-J_∥)t
    -   1/10cos(J_⊥ +4J_∥)t)∑_k<l 
     k,l≠ iz_kz_l
     +   (157/600+1/12cos3J_⊥t+3/100cos5J_⊥t+1/10cos(J_⊥+4J_∥)t
    -   1/4cos(J_⊥+2J_∥)t-3/50cos(3J_⊥+2J_∥)t
    -   1/100cos(7J_⊥-2J_∥)t-9/50cos2(J_⊥-J_∥)t
    +   1/40cos4(J_⊥-J_∥)t)   ∑_j<k<l<m 
     j,k,l,m≠ iz_jz_kz_lz_m


    β_i(t,5)=   (3/100sin(3J_⊥+2J_∥)t-1/40sin4(J_⊥-J_∥)t
    -   1/200sin(7J_⊥-2J_∥)t-9/100sin2(J_⊥-J_∥)t
    +   1/8sin(J_⊥+2J_∥)t+1/10sin(J_⊥+4J_∥)t) ∑_k≠ iz_k
     -   (3/100sin(3J_⊥+2J_∥)t+1/40sin4(J_⊥-J_∥)t
    -   1/200sin(7J_⊥-2J_∥)t-9/100sin2(J_⊥-J_∥)t
    +   1/8sin(J_⊥+2J_∥)t-1/10sin(J_⊥+4J_∥)t)∑_k<l<j 
     j,k,l≠ iz_jz_kz_l




  §.§.§ 6-qubits

Repeating the procedures of the previous sections we may block diagonalize the 6-qubit XXZ-network. Going to the Fourier basis we have the remaining blocks to diagonalize

    Γ_6(H^2_ XXZ⊕ H^3_ XXZ⊕ H^4_ XXZ)Γ^†_6=⊕_a=0^5H^2,a_ XXZ⊕ H^3,a_ XXZ⊕ H^4,a_ XXZ

which provides more of a challenge than the previous examples as now 4× 4 blocks and 3× 3 blocks appear. But the example still retains a degree of simplicity, as the eigenvectors do not contain any parameter dependence. The computation may be further reduced by noting that once H^2,a_ XXZ is diagonalized, a similar transformation may be used to diagonalize H^4,a_ XXZ. As for the size of the various blocks we have for q=2

    dimH^2,0_ XXZ= dimH^2,2_ XXZ= dimH^2,4_ XXZ=3 
        dimH^2,1_ XXZ= dimH^2,3_ XXZ= dimH^2,5_ XXZ=2

where the same holds for q=4. All of these blocks create superpositions amongst the Fourier modes except H^2,3_ XXZ, which is comprised of two singlet states |ℱ^3_2,0⟩ and |ℱ^3_2,1⟩. The same conclusion applies for q=4. For q=3 we have

    dimH^3,0_ XXZ= dimH^3,3_ XXZ=4 
        dimH^3,1_ XXZ= dimH^3,2_ XXZ= dimH^3,4_ XXZ= dimH^3,5_ XXZ=3

and there are mixtures created among all the Fourier modes except |ℱ^2_3,0⟩ and |ℱ^4_3,0⟩.

    λ_3,i(t,6)    =(59/144 +5/32cos2J_⊥t +5/16cos4J_⊥t+35/288cos6J_⊥t) 
       +(1/24 -1/32cos2J_⊥t -1/40cos4J_⊥t+7/480cos6J_⊥t)∑_k<l 
     k,l≠ iz_kz_l
       -(1/48 -1/32cos2J_⊥t +1/80cos4J_⊥t-1/480cos6J_⊥t)   ∑_j<k<l,m 
    j, k,l,m≠ iz_jz_kz_lz_m


    τ_3,i(t,6) =   (17/144 -1/32cos2J_⊥t -1/16cos4J_⊥t-7/288cos6J_⊥t)∑_k≠ iz_k
    -   (1/24 -1/32cos2J_⊥t -1/40cos4J_⊥t+7/480cos6J_⊥t)∑_k<l<j 
     k,l,j≠ iz_kz_lz_j
    +   (5/48 -5/32cos2J_⊥t +1/16cos4J_⊥t-1/96cos6J_⊥t)∏_j≠ iz_j

The expressions for α(t) and β(t) are lengthy, therefore we shall express our results using the partial components. The remaining non-zero partial components are obtained using the permutation symmetry, for example α^zz000_i=α^z0z00_i=α^00zz0_i, etc. We find 

    α^00000_i(t,6)=   5/96cos(J_⊥+5J_∥)t+1/96cos5(J_⊥-J_∥)t
    +   3/16cos(J_⊥+3J_∥)t+25/288cos3(J_⊥-J_∥)t
    +   1/288cos(9J_⊥-3J_∥)t+1/48cos(5J_⊥+J_∥)t
    +   3/32cos(3J_⊥+J_∥)t+5/32cos(J_⊥+J_∥)t
    +   5/16cos(J_⊥-J_∥)t+1/32cos(5J_⊥-J_∥)t
    +   1/96cos(7J_⊥-J_∥)t+5/144cos3(J_⊥+J_∥)t


    α^zz000_i(t,6)=   5/96cos(J_⊥+5J_∥)t+1/96cos5(J_⊥-J_∥)t
    +   3/80cos(J_⊥+3J_∥)t-1/240cos(5J_⊥+J_∥)t
    -   3/160cos(3J_⊥+J_∥)t-1/32cos(J_⊥+J_∥)t
    +   5/288cos3(J_⊥-J_∥)t+1/1440cos(9J_⊥-3J_∥)t
    -   1/16cos(J_⊥-J_∥)t-1/160cos(5J_⊥-J_∥)t
    -   1/480cos(7J_⊥-J_∥)t+1/144cos3(J_⊥+J_∥)t


    α^zzzz0_i(t,6)=   5/96cos(J_⊥+5J_∥)t+1/96cos5(J_⊥-J_∥)t
    -   9/80cos(J_⊥+3J_∥)t-5/96cos3(J_⊥-J_∥)t
    -   1/480cos(9J_⊥-3J_∥)t+1/240cos(5J_⊥+J_∥)t
    +   3/160cos(3J_⊥+J_∥)t+1/32cos(J_⊥+J_∥)t
    +   1/16cos(J_⊥-J_∥)t+1/160cos(5J_⊥-J_∥)t
    +   1/480cos(7J_⊥-J_∥)t-1/48cos3(J_⊥+J_∥)t


    β^z0000_i(t,6)=   5/96sin(J_⊥+5J_∥)t-1/96sin5(J_⊥-J_∥)t
    +   9/80sin(J_⊥+3J_∥)t-5/96sin3(J_⊥-J_∥)t
    -   1/480sin(9J_⊥-3J_∥)t+1/240sin(5J_⊥+J_∥)t
    +   3/160sin(3J_⊥+J_∥)t+1/32sin(J_⊥+J_∥)t
    -   1/16sin(J_⊥-J_∥)t-1/160sin(5J_⊥-J_∥)t
    -   1/480sin(7J_⊥-J_∥)t+1/48sin3(J_⊥+J_∥)t


    β^zzz00_i(t,6)=   5/96sin(J_⊥+5J_∥)t-1/96sin5(J_⊥-J_∥)t
    -   3/80sin(J_⊥+3J_∥)t+5/288sin3(J_⊥-J_∥)t
    +   1/1440sin(9J_⊥-3J_∥)t-1/240sin(5J_⊥+J_∥)t
    -   3/160sin(3J_⊥+J_∥)t-1/32sin(J_⊥+J_∥)t
    +   1/16sin(J_⊥-J_∥)t+1/160sin(5J_⊥-J_∥)t
    +   1/480sin(7J_⊥-J_∥)t-1/144sin3(J_⊥+J_∥)t


    β^zzzzz_i(t,6)=   5/96sin(J_⊥+5J_∥)t-1/96sin5(J_⊥-J_∥)t
    -   3/16sin(J_⊥+3J_∥)t+25/288sin3(J_⊥-J_∥)t
    +   1/288sin(9J_⊥-3J_∥)t+1/48sin(5J_⊥+J_∥)t
    +   3/32sin(3J_⊥+J_∥)t+5/32sin(J_⊥+J_∥)t
    -   5/16sin(J_⊥-J_∥)t-1/32sin(5J_⊥-J_∥)t
    -   1/96sin(7J_⊥-J_∥)t-5/144sin3(J_⊥+J_∥)t



 §.§ Ring-connected XXZ-network

Turning to the more difficult model, the homogeneous XXZ-ring. Due to the more complicated structure of this model we only investigate the XXZ model for N=4, and use the XXX model for N=5, where we only present the computations for λ_3(t) and τ_3(t). 


  §.§.§ 4-qubits

For 4-qubits the Hamiltonian has the block diagonal structure

    𝒫_4H_ XXZ𝒫_4^†= H^0_ XXZ⊕ H^1_ XXZ⊕ H^2_ XXZ⊕ H^3_ XXZ⊕ H^4_ XXZ

 The only non-trivial diagonalization occurs in the q=2 block. Restricting to this sector, we may reduce it to smaller blocks by invoking the translation symmetry

    Γ_4 H^2_ XXZΓ_4^†=H^2,0_ XXZ⊕ H^2,1_ XXZ⊕ H^2,2_ XXZ⊕ H^2,3_ XXZ

where Γ_4 is a unitary transformation that diagonalizes T_4. The a=1,3 blocks are singlets, while a=0,2 are both doublets. Mixing only occurs within the a=0 block, which in this case does depend on the parameters of the Hamiltonian.  
Theses stationary states of H^2,2_ XXZ are determined to be 

    |E^0_2;0⟩=cosϕ_×/2|ℱ^0_2;0⟩+sinϕ_×/2|ℱ^0_2;1⟩
       |E^0_2;1⟩=sinϕ_×/2|ℱ^0_2;0⟩-cosϕ_×/2|ℱ^0_2;1⟩
       ω_±=J_∥± sgn(J_∥) J_×


where the Hamiltonian dependent parameters we have introduced are

    J_×=√(J_∥^2+8J_⊥^2)
       ϕ_×=- tan^-12^3/2J_⊥/J_∥ .





 
    λ_3,i(t,4)    =(7/16+1/4cos2J_⊥t+1/16cos4J_⊥t+1/4cos J_∥tcos J_×t
       +1/4cosϕ_×sin J_∥tsin J_×t)
       +(1/16+cos4J_⊥t-1/8cos J_∥tcos J_×t
       -1/8cosϕ_×sin J_∥tsin J_×t )(z_i-1z_i+2+z_i+1z_i+2)
       -(3/16-1/4cos2J_⊥t+1/16cos4J_⊥t)(z_i-1z_i+1)


    τ_3,i(t,4)    =(1/16-1/4cos2J_⊥t-1/16cos4J_⊥t+1/4cos J_∥tcos J_×t
       +1/4cosϕ_×sin J_∥tsin J_×t)∏_k ≠ iz_k
       +(3/16-1/16cos4J_⊥t-1/8cosϕ_×sin J_∥tsin J_×t
       -1/8cos J_∥tcos J_×t) (z_i-1+z_i+1)
       +(3/16-1/4cos2J_⊥t+1/16cos4J_⊥t)z_i+2

As seen in the previous section, the expression for α(t) and β(t) are more cumbersome than the other non-unitary parameters. Therefore we again split them into their partial components. For the ring topology we only need to consider the reflection symmetry to generate additional partial components as will be clear in the following equations. We find

    α^000_i(t,4)   =1/4+1+3cos2J_∥t/16cos^2J_⊥t +1+3cos2J_⊥t/16cos J_∥tcos J_×t 
       -√(2)/8sinϕ_×sin2J_⊥tcos J_∥tsin J_×t
       -(1+cosϕ_×sin J_⊥tsin J_×t)/8sin^2J_∥t
 

    α^0zz_i(t,4)   =cos^2J_⊥t/8(cos2J_∥t-cos J_∥tcos J_×t +cosϕ_×sin J_∥tsin J_×t) 
       +√(2)/16sinϕ_×cos J_∥tsin2J_⊥tsin J_×t= α^zz0_i(t,4)



    α^z0z_i(t,4)=   -1/8+cos^2J_⊥t/4(1+3/2cos2J_∥t)-sin^2J_⊥t/8cos J_∥tcos J_×t
       +1+3cos2J_⊥t/16cosϕ_×sin J_∥tsin J_×t


    β^z00_i(t,4)=   -cos^2J_⊥t/8(3sin2J_∥t+sin J_∥tcos J_×t+cosϕ_×cos J_∥tsin J_×t)
       +√(2)/16sinϕ_×sin J_∥tsin2J_⊥tsin J_×t=β^00z_i(t,4)


    β^0z0_i(t,4)   =1+3cos2J_⊥t/16cosϕ_×cos J_∥tsin J_×t-cos^2J_⊥t/8sin2J_∥t
       +sin^2J_⊥t/8sin J_∥tcos J_×t


    β^zzz_i(t,4)=   sin^2J_⊥t/8cosϕ_×tcos J_∥tsin J_×t-cos^2J_⊥t/8sin2J_∥t
    +   1+3cos2J_⊥t/16sin J_∥tcos J_×t
    -   √(2)/16sinϕ_×sin J_∥tsin2J_⊥tsin J_×t




  §.§.§ 5-qubits

As with the previous examples we change from the computational basis to the canonically  ordered excitation basis. We achieve this by using a permutation matrix

    𝒫_5H_ XXZ𝒫_5^†= H^0_ XXZ⊕ H^1_ XXZ⊕ H^2_ XXZ⊕ H^3_ XXZ⊕ H^4_ XXZ⊕ H^5_ XXZ

and we are interested in finding the stationary states of the middle two blocks, as the others may be found trivially. Now we can go further by changing to a translation eigenbasis

    Γ_5H^2_ XXZ⊕ H^3_ XXZΓ^†_5=⊕ _a=0^4 H^2,a_ XXZ⊕ H^3,a_ XXZ

where the q=2 blocks take the form

    H^2,a_ XXZ=[                   -h+1/2J_∥ (2cos6π a/5J_⊥) e^-6π i a/5;                            ;   (2cos6π a/5J_⊥)e^6π i a/5     -h-3/2J_∥+2cos6π a/5J_⊥ ]

and the q=3 blocks are obtained by sending h→ -h. We thus obtain the following parameterization of the non-trivial stationary states

    tanψ_a=6π a/5
       tanϕ_a=2cos6π a/5J_⊥/(J_∥-cos6π a/5J_⊥)
        J_a=√((J_∥-cos6π a/5J_⊥)^2+4cos^26π a/5J^2_⊥)

Due to the length of the expression obtained, we must report the non-unitary parameters by the partial components defined in the main text. Setting J_⊥=J_∥ we find
 
    λ^0000_3,i(t;5)   =71/225+13/90cos(3-√(5)/2)J_⊥t+1/10cos(5-√(5)/2)J_⊥t
       +1/45cos3(1-√(5)/2)J_⊥t+32/225cos√(5)J_⊥t
       +13/90cos(3+√(5)/2)J_⊥t+1/10cos(5+√(5)/2)J_⊥t
       +2/225cos2√(5)J_⊥t+1/45cos3(1+√(5)/2)J_⊥t

    
    λ^zz00_3,i(t,5)=   -6√(5)+25/900cos(3+√(5)/2)J_⊥t+6√(5)-25/900cos(3-√(5)/2)J_⊥t
    +   1/100cos(5+√(5)/2)J_⊥t+1/100cos(5-√(5)/2)J_⊥t
    -   √(5)/450cos3(1-√(5)/2)J_⊥t+√(5)/450cos3(1+√(5)/2)J_⊥t
    -   2/225cos2√(5)J_⊥t+2/45cos√(5)J_⊥t= λ^00zz_3,i(t,5)
    
    
    λ^z0z0_3,i(t,5)=   1/100cos(5-√(5)/2)J_⊥t-6√(5)+25/900cos(3-√(5)/2)J_⊥t
       +√(5)/450cos3(1-√(5)/2)J_⊥t-√(5)/450cos3(1+√(5)/2)J_⊥t
       +6√(5)-25/900cos(3+√(5)/2)J_⊥t+1/100cos(5+√(5)/2)J_⊥t
       -2/225cos2√(5)J_⊥t+2/45cos√(5)J_⊥t= λ^0z0z_3,i(t,5)
     
     
    λ^0zz0_3,i(t,5)   =1/45+3√(5)-5/300cos(3-√(5)/2)J_⊥t+1/75cos2√(5)J_⊥t
       +1-√(5)/100cos(5-√(5)/2)J_⊥t-5-√(5)/450cos3(1-√(5)/2)J_⊥t
       -3√(5)+5/300cos(3+√(5)/2)J_⊥t+1+√(5)/100cos(5+√(5)/2)J_⊥t
       -5+√(5)/450cos3(1+√(5)/2)J_⊥t

    
    λ^z00z_3,i(t,5)   =1/45-3√(5)+5/300cos(3-√(5)/2)J_⊥t+1+√(5)/100cos(5-√(5)/2)J_⊥t
       -5+√(5)/450cos3(1-√(5)/2)J_⊥t+3√(5)-5/300cos(3+√(5)/2)J_⊥t
       +1-√(5)/100cos(5+√(5)/2)J_⊥t+1/75cos2√(5)J_⊥t
       +5-√(5)/450cos3(1+√(5)/2)J_⊥t

 
    τ^z000_3,i(t;5)   =77/450-1+√(5)/40cos(5+√(5)/2)J_⊥t-1-√(5)/40cos(5-√(5)/2)J_⊥t
       -1-√(5)/180cos3(1-√(5)/2)J_⊥t-1+√(5)/180cos3(1+√(5)/2)J_⊥t
       -13+3√(5)/360cos(3+√(5)/2)J_⊥t-13-3√(5)/360cos(3-√(5)/2)J_⊥t
       -1/450cos2√(5)J_⊥t-8/225cos√(5)J_⊥t=τ^000z_3,i(t;5)

 
    τ^0z00_3,i(t;5)   =77/450-1-√(5)/40cos(5+√(5)/2)J_⊥t-1+√(5)/40cos(5-√(5)/2)J_⊥t
       -1+√(5)/180cos3(1-√(5)/2)J_⊥t-1-√(5)/180cos3(1+√(5)/2)J_⊥t
       -13-3√(5)/360cos(3+√(5)/2)J_⊥t-13+3√(5)/360cos(3-√(5)/2)J_⊥t
       -1/450cos2√(5)J_⊥t-8/225cos√(5)J_⊥t=τ^00z0_3,i(t;5)

 
    τ^zzz0_3,i(t;5)=   -1/90-3-√(5)/200cos(5+√(5)/2)J_⊥t-3+√(5)/200cos(5-√(5)/2)J_⊥t
       +5+3√(5)/900cos3(1-√(5)/2)J_⊥t+5-3√(5)/900cos3(1+√(5)/2)J_⊥t
       +65+3√(5)/1800cos(3+√(5)/2)J_⊥t+65-3√(5)/1800cos(3-√(5)/2)J_⊥t
       +1/450cos2√(5)J_⊥t-2/45cos(√(5)J_⊥t)=τ^0zzz0_3,i(t;5)

 
    τ^z0zz_3,i(t;5)=   -1/90-3+√(5)/200cos(5+√(5)/2)J_⊥t-3-√(5)/200cos(5-√(5)/2)J_⊥t
       +5-3√(5)/900cos3(1-√(5)/2)J_⊥t+5+3√(5)/900cos3(1+√(5)/2)J_⊥t
       +65-3√(5)/1800cos(3+√(5)/2)J_⊥t+65+3√(5)/1800cos(3-√(5)/2)J_⊥t
       +1/450cos2√(5)J_⊥t+2/45cos(√(5)J_⊥t)=τ^zz0z_3,i(t;5)



 §.§ Disconnected XX-network

 
 The last dynamical map components we determine are those associated to the disconnected XX network presented in Section 4. The eigen-parameterization used for H^1,2_ XX is 

    h_12=h_1+h_2/2  ω_12= sgn(Δ_12)√(Δ_12^2+J_12^2)
       Δ_12=h_1-h_2  tanϕ_12=J_12/Δ_12

Using these parameters, the Pauli components of the unitary channel may be expressed as

    𝒰^x0_zy=𝒰^zx_y0=𝒰^0x_yz=𝒰^xz_0y=-𝒰^y0_zx=-𝒰^zy_x0=-𝒰^0y_xz=-𝒰^yz_0x=sinϕ_12sinω_12 tcos2 h_12 t
       𝒰^xx_0z=𝒰^yy_z0=𝒰^z0_xx=𝒰^z0_yy=-𝒰^xx_0z=-𝒰^yy_0z=-𝒰^0z_xx=-𝒰^0z_yy=cosϕ_12sinϕ_12sin^2ω_12 t
       𝒰^xy_z0=𝒰^yx_z0=𝒰^z0_yx=𝒰^0z_xy=-𝒰^xy_0z=-𝒰^yx_z0=-𝒰^z0_xy=-𝒰^0z_yx=1/2sinϕ_12sin2ω_12 t
       𝒰^0x_xz=𝒰^0y_yz=𝒰^xz_0x=𝒰^yz_0y=𝒰^x0_zx=𝒰^y0_zy=𝒰^zx_x0=𝒰^zy_y0=sinϕ_12sinω_12 tsin2 h_12 t
       𝒰^y0_x0=-𝒰^x0_y0=𝒰^yz_xz=-𝒰^xz_yz=cosω_12 tsin2h_12t+cosϕ_12sinω_12 tcos2h_12t
       𝒰^0y_0x=-𝒰^0x_0y=𝒰^zy_zx=-𝒰^zx_zy=cosω_12 tsin2h_12t-cosϕ_12sinω_12 tcos2h_12t
       𝒰^x0_x0=𝒰^y0_y0=𝒰^xz_xz=𝒰^yz_yz=cosω_12 tcos2h_12t-cosϕ_12sinω_12 tsin2h_12t

with the remaining components of the unitary channel given by

    𝒰^0x_0x=𝒰^0y_0y=𝒰^zx_zx=𝒰^zy_zy=cosω_12 tcos2h_12t+cosϕ_12sinω_12 tsin2h_12t
       𝒰^yx_xx=𝒰^yy_xy=-𝒰^xx_yx=-𝒰^xy_yy=1/2(sin4h_12t+cosϕ_12sin2ω_12 t)
       𝒰^xy_xx=𝒰^yy_yx=-𝒰^xx_xy=-𝒰^yx_yy=1/2(sin4h_12t-cosϕ_12sin2ω_12 t) 
       𝒰^yy_xx=𝒰^xx_yy=1/2(-cos4h_12t+cos^2ϕ_12cos2ω_12 t+sin^2ϕ_12) 
       𝒰^xx_xx=𝒰^yy_yy=1/2(cos4h_12t+cos^2ϕ_12cos2ω_12 t+sin^2ϕ_12) 
       𝒰^xy_xy=𝒰^yx_yx=1/2(cos4h_12t+cos2ω_12 t) 
       𝒰^xy_yx=𝒰^yx_xy=1/2(cos4h_12t-cos2ω_12 t)  
       𝒰^z0_z0=𝒰^0z_0z=1-sin^2ϕ_12sin^2ω_12 t
       𝒰^z0_0z=𝒰^0z_z0=sin^2ϕ_12sin^2ω_12 t

Constructing the 1-qubit dynamical map, a non-unital dynamical map is obtained as given in Section 2. But enforcing a noise-distribution over ϕ_12 that is an even function, the noise-averaged dynamical map takes the form Eq.(<ref>).

unsrt
 
