Decoherence of two qubit systems: A random matrix description

We study decoherence of two non-interacting qubits. The environment and its interaction with the qubits are modelled by random matrices. Decoherence, measured in terms of purity, is calculated in linear response approximation. Monte Carlo simulations illustrate the validity of this approximation and of its extension by exponentiation. The results up to this point are also used to study one qubit decoherence. Purity decay of entangled and product states are qualitatively similar though for the latter case it is slower. Numerical studies for a Bell pair as initial state reveal a one to one correspondence between its decoherence and its internal entanglement decay. For strong and intermediate coupling to the environment this correspondence agrees with the one for Werner states. In the limit of a large environment the evolution induces a unital channel in the two qubits, providing a partial explanation for the relation above.


Introduction
In recent experiments [1,2,3,4] it has been demonstrated that it is possible to protect ever larger entangled quantum systems, often arrays of qubits, ever more efficiently from decoherence. This improved protection of states makes a more detailed analysis of decoherence desirable. In the context of fidelity decay, a random matrix description [5,6] is accessible and very effective in describing experiments [7,8,9]. A close connection between the dynamics of fidelity decay and decoherence has been shown in some instances [10,11], which suggests to apply methods successful in one field to the other. Based on some previous work we shall therefore analyse in detail the consequences of using random matrix theory (RMT) to model decoherence [6,12,13] for the case of two non-interacting qubits or a single qubit as the central system.
Two perspectives make such a random matrix treatment particularly attractive. First, reduction of decoherence may, in some instances, be achieved by isolating some "far" environment (including spontaneous decay) to a degree that it can, to first approximation, be neglected. Then it can happen that the Heisenberg time of the relevant "near" environment is finite on the time scale of decoherence. In such a case Decoherence of two qubit systems: A random matrix description 2 it becomes relevant that RMT shows, in linear response approximation, a transition from linear to quadratic decay at times of the order of the Heisenberg time. This behaviour is seen with spin chain environments [14], and is essential for the success of the theory in describing the above mentioned experiments of fidelity decay. Note also that the concept of a two stage environment has been used for basic considerations [15]. Second, the long term goal must be to describe in one theory the decay of fidelity that includes undesirable deviations of the internal Hamiltonian of the central system, together with decoherence.
This study is focussed on weak coupling of one or two qubits to an environment, and thus we can use the correlation function approach proposed for purity decay in echo-dynamics [16], treating the coupling as the perturbation [6]. The linear response approximation will be sufficient and in this approximation the ensemble averages, which we have to take in any RMT model, are feasible though somewhat tedious. Exact solutions, which exist in some instances for the decay of the fidelity amplitude [17,18], seem to be out of reach at present, because they would require the evaluation of four-point functions.
Assume that the qubit(s) are initially in a pure state, and evolve under their own local Hamiltonians. The qubit(s) are coupled independently via random matrices to a large environment in turn described by another random matrix. The coupling to the environment gives rise to decoherence. Averaging both the coupling and the environment Hamiltonian over the RMT ensembles yields the generic behaviour. We then focus on the correspondence between decoherence and entanglement decay for a Bell state. We find a relation between this two quantities, that can be partially understood assuming unitality. This relation in most cases coincides with the corresponding one obtained for Werner states. For that part, it is essential that the two qubits do not interact. Otherwise, the coupling between the qubits would act as an additional sink (or source) for internal entanglement -a complication we wish to avoid. However, a residual coupling between the qubits should be taken into account at some point [19]. A preliminary study of some aspects of this system has been presented in [13], but here we shall present a more general picture: First, we allow for local onequbit Hamiltonians. They have a considerable effect on decoherence. Second, we cover exhaustively all possible initial states. Third, we study the consequences of the whole system being time-reversal invariant. Indeed, we present the mathematical procedure in detail using both the Gaussian unitary (GUE) and the Gaussian orthogonal (GOE) ensembles [20,21] for the description of the environment and the coupling. The two ensembles correspond to time reversal invariance (TRI) breaking and conserving dynamics respectively. All these cases are treated on equal footing.
In section 2 we shall state the model, recall the linear response formalism for echo dynamics, and show how it can be adapted to forward evolution. In section 3 we shall discuss decoherence of a single qubit, which already shows some interesting features by itself and serves as a warm up for the more complicated case of two qubits. The details of the calculations are given in the appendices. This shall be treated in section 4, starting from the situation where one qubit is a spectator, exerting influence only through entanglement. More general cases are analysed and reduced to the spectator situation. Most analytic results will be accompanied by Monte Carlo simulations. In section 5 we investigate the correspondence between decoherence (purity) and entanglement (concurrence) as both evolve in time. We end with conclusions in section 6.
Decoherence of two qubit systems: A random matrix description 3

General considerations
In this section, we describe the general framework for our study of decoherence. We consider one-and two-qubit systems, coupled to an environment which is modelled using random matrix theory. The one-qubit case is introduced in section 2.1, the twoqubit case in section 2.2. The quantities studied in this paper are defined in section 2.3 and the tools used to solve the problem are introduced in section 2.4.

The model for one qubit
We describe decoherence by considering explicitly the additional degrees of freedom (henceforth called "environment") which are interacting with the qubit. We follow the unitary time evolution of a pure, initially separable, state in the product space H = H 1 ⊗ H e , where H 1 (of dimension two) and H e (of dimension N e ) denote the Hilbert spaces of the qubit and the environment, respectively. As time evolves, the qubit and the environment become more and more entangled, which means that after tracing out the environmental degrees of freedom, the state of the qubit becomes more and more mixed. For such a setup, the Hamiltonian is of the following form H λ = H 1 ⊗ 1 1 e + 1 1 1 ⊗ H e + λV 1,e ≡ H 1 + H e + λV 1,e . (1) Here, H 1 represents the Hamiltonian acting on the qubit, H e the Hamiltonian of the environment, and V 1,e the coupling between the qubit and the environment. The real parameter λ controls the strength of the coupling. For notational ease, and if there is no danger of confusion, we omit extensions with the identity, as it is done in (1). We describe both the coupling and the dynamics in the environment within random matrix theory. To this end, H e and V e,1 are chosen from either the GUE or the GOE, depending on whether we wish to describe a TRI breaking or TRI conserving situation. The Hamiltonian H 1 implies another free parameter of the model, namely the level splitting ∆ of the two level system representing the qubit.
We shall study the time evolution of an initially pure and separable state where |ψ 1 ∈ H 1 and |ψ e ∈ H e . At any time t, the state of the whole system is thus |ψ(t) = exp(−itH λ )|ψ(0) , and the state of the single qubit is tr e |ψ(t) ψ(t)|. While the state of the qubit |ψ 1 is yet another free parameter in our model, we assume the state of the environment |ψ e to be random. This means that the state is chosen from an ensemble which is invariant under unitary transformations. In practice, this means that the coefficients are chosen as complex random Gaussian variables, and subsequently the state is normalized.

The model for two qubits
For the two qubit case, we need to introduce the additional Hilbert spaces H 2 of dimension two and H e ′ of dimension N e ′ , which correspond to the second qubit and its environment respectively. As in the one-qubit case, the dynamics in the space of the two qubits H 1 ⊗H 2 (henceforth called "central system") is obtained by considering the unitary evolution in whole Hilbert space H 1 ⊗ H 2 ⊗ H e ⊗ H e ′ and subsequently tracing-out the environmental degrees of freedom (henceforth called "environment").
Typically, we assume that at t = 0 the two qubits (and hence the state |ψ 12 ) to be entangled. This is the main new ingredient, as compared to the one-qubit case. In this work we shall consider three different dynamical scenarios, all explicitly excluding any interaction between the two qubits: where λ, H 1 , H e and V 1,e are defined as in (1). This situation is shown schematically in figure 1(a). If we choose an initial state where the two qubits are already entangled, this provides the simplest situation which allows to study entanglement decay. If the two qubits are initially not entangled, the process reduces effectively to the one-qubit case, described previously section 2.1. The special case H 1 = H 2 = 0 has been considered in Ref. [13].
(b) The separate environment Hamiltonian: This case differs from the spectator model only inasmuch as the second qubit is also coupled to an environment. Both environments are assumed to be non-interacting. The total Hamiltonian then reads where V 2,e ′ and H e ′ describe the coupling to -and the dynamics in the additional environment. Both quantities are chosen independently from the respective random matrix ensembles, in perfect analogy with V 1,e and H e . The real parameters λ 1 and λ 2 fix the coupling strengths to either environment. This model, see figure 1(b), may describe two qubits that are ready to perform a distant teleportation, where each of them is interacting only with its immediate surroundings. It can also represent a pair of qubits that, although close together, interact with different and independent degrees of freedom.
Decoherence of two qubit systems: A random matrix description 5 (c) The joint environment Hamiltonian: The third case, shown in figure 1(c), describes a situation in which both qubits are coupled to the same environment, even though the coupling matrices are still independent. In that case, the total Hamiltonian reads where V 2,e describes the coupling of the second qubit to the environment. It is chosen independently from the same random matrix ensemble as V 1,e .

The measures
As a measure of decoherence we use purity. This measures the degree of mixedness of a density matrix ρ representing a physical system. Purity is defined as It reaches a maximum value of one when the state is pure, i.e. when ρ = |ψ ψ| for some |ψ . Otherwise it is less than one and reaches a minimum when the state is completely mixed. We use purity instead of e.g. the von Neumann entropy because it is simpler to handle from an algebraic point of view. Note that for one qubit both quantities are equivalent, and for two qubits they contain similar information [22]. The other measure considered here is concurrence. It quantifies the degree of entanglement of two qubits in a pure or mixed state. It is a measure used extensively in the literature [23] that is straightforward to calculate. It is closely related to the entanglement of formation, which measures the minimum number of Bell pairs needed to create an ensemble of pure states representing the state to be studied [24]. Given a density matrix ρ representing the state of two qubits, concurrence is defined as where Λ i are the eigenvalues of the matrix ρ(σ y ⊗ σ y )ρ * (σ y ⊗ σ y ) in non-increasing order. The superscript * denotes complex conjugation in the computational basis and σ y is a Pauli matrix [25]. The concurrence has a maximum value of one for Bell states, and a minimum value of zero for separable states. Furthermore, it is invariant under bilocal unitary operations. Both, purity and concurrence, provide a measure of entanglement, but in the present context they quantify completely different aspects of the problem. Purity will be used to quantify the entanglement between the pair of qubits and the environment i.e. decoherence. Concurrence will be used to measure the entanglement within the pair. Although they are not fully independent, information about one provides little knowledge about the other. E.g. a state with purity one can have full entanglement (a Bell pair) or null entanglement (a separable pure state).

Echo dynamics and linear response theory
We shall calculate the value of purity as a function of time analytically, in a perturbative approximation. To perform this task it is useful to consider each of the above Hamiltonians [(1), (4), (5), and (6)], as composed by an unperturbed part H 0 and a perturbation λV with λ a parameter. The unperturbed part corresponds to the operators that act on each individual subspace alone whereas the perturbation corresponds to the coupling among the different subspaces; for example in the one qubit case H 0 = H e + H 1 and V = V e,1 . For the case of two qubits, λ can be taken as max{λ 1 , λ 2 }. We shall use the tools developed in Appendix A for a linear response formalism in echo dynamics. In order to do so we must formally write the problem in this language.
We write the Hamiltonian as and introduce the evolution operator and the echo operator defined by respectively ( = 1). For the calculation of purity at a given time t, we replace the forward evolution operator U λ by the corresponding echo operator M λ . Even though the resulting states are different, i.e. ρ(t) = tr e,e ′ U λ (t)ρU † λ (t) = ρ M (t) = tr e,e ′ M λ (t)ρM † λ (t), they are still related by the local (in the two qubits and the environment) unitary transformation U 0 (t). Since such transformations do not change the entanglement properties, it holds that P (t) = P [ρ(t)] = P [ρ M (t)]. This step is crucial, since the echo operator admits a series expansion with much larger range of validity (both, in time and perturbation strength). The numerical simulations are all done with forward evolution alone as they require less computational effort.
The Born expansion of the echo operator up to second order reads with andṼ (t) = U 0 (t) † V U 0 (t) being the coupling in the interaction picture. Using this expansion we calculate the purity of the central system, averaged over the coupling and the Hamiltonian of the environment.

One qubit decoherence
We first study the GUE case with and without an internal Hamiltonian governing the qubit. The next step is to work out the GOE case. There we concentrate on the case with no internal Hamiltonian governing in the qubit since we want to keep the discussion as simple as possible to focus on the consequences of the weaker invariance properties of the ensemble.
Here, we study the evolution under the Hamiltonian where the initial state is the product of a fixed pure state |ψ 1 ∈ H 1 and a random pure state |ψ e ∈ H e ; see (2). At any later time, the state of the qubit and its purity are given by In Appendix A, we compute the average purity P (t) as a function of time in the linear response approximation (11), following the steps outlined in section 2.4. The Decoherence of two qubit systems: A random matrix description 7 average is taken with respect to the coupling V 1,e , the random initial state |ψ e and the spectrum of H e . In the limit of N e → ∞, we obtain [(A.7), (A.10)] where χ GOE = 1 for the TRI case, and χ GOE = 0 for the non-TRI case. The correlation functions C 1 (τ ), S 1 (τ ), S ′ 1 (τ ), andC(τ ) are defined in Appendix B.C(τ ) deserves special attention, since all the dependence of the environment is via the function where the E j 's are the eigenenergies of H e and τ H is the corresponding Heisenberg time. The two-point form factor, b (β) 2 (τ ), is known analytically for the GUE (β = 2) and the GOE (β = 1), which are the two cases treated here [21].

The GUE case
We are now in the position to give an explicit formula for P (t) in the GUE case. This formula will generally depend on some properties of the initial condition |ψ 1 . We wish to write |ψ 1 in the most general way, but grouping together cases, equivalent due to the invariance properties of the problem.
Recall that H = H e + H 1 + λV represents an ensemble of Hamiltonians in which H e ∈ GUE and V ∈ GUE, whereas H 1 (together with the initial condition |ψ 1 ) remains fixed throughout the calculation. The operations under which the ensemble is invariant are local (with respect to the partitioning of the Hilbert space into H 1 and H e ), unitary (due to the invariance properties of the GUE), and leave H 1 invariant. Hence the transformation matrices must be of the form U ⊗ exp(iαH 1 ) with α a real number and U a unitary operator acting on H e . This freedom allows to choose a convenient basis to solve the problem. On one hand, it allows to write H 0 in diagonal form (as done in Appendix A), and on the other hand, we can use it to represent the initial state of the qubit in such a way that there is no phase shift between the two components of the qubit. This can be achieved by appropriately choosing α (see also the discussion in section 4.1). We thus write, without loosing generality where |0 and |1 denote the eigenstates of H 1 . Notice that if φ = 0 and φ = π/2, |ψ 1 is an eigenstate of H 1 . Finally, we choose the origin of the energy scale in such a way that the Hamiltonian of the qubit can be written as Hence, φ is the only relevant parameter describing the initial state. We obtain the average purity from the general expression in (16) and (17). For a pure initial state ρ 1 = |ψ 1 ψ 1 | the relevant correlation functions Re C 1 (τ ), S 1 (τ ) and C(τ ) are given in (B.10), (B.12), and (B.5), respectively. Using the symmetry of the resulting integrand with respect to the exchange of τ and τ ′ , we find with Decoherence of two qubit systems: A random matrix description 8 quantifying the "distance" between φ and the eigenbasis of H 1 . Let us consider following two limits for H 1 . The "degenerate limit", where the level splitting ∆ is much smaller than the mean level spacing d e = 2π/τ H of the environmental Hamiltonian, and the "fast limit", where the level splitting is much larger. In the latter case, the internal evolution of the qubit is fast compared with the evolution in the environment.
The degenerate limit leads to the known formula [13] with The result does not depend on the initial state of the qubit. Due to the degeneracy all states are eigenstates of H 1 and thus equivalent. The leading term of the purity decay is linear before the Heisenberg time and quadratic after the Heisenberg time. Similar features were already observed in fidelity decay and purity decay in other contexts [6].
In the fast limit (∆ ≫ d e ), purity is obtained from (20) by replacing cos ∆τ ′ with one when it is multiplied with the δ function [see (18)], and with zero everywhere else. For finite N e care must be taken, since we are assuming Zeno time (which is given by the "width of the δ-function") to be much smaller than all other time scales, such that ∆ ≪ N e d e . The resulting expression is Typically (depending on the initial state), this formula again displays a dominantly linear decay below the Heisenberg time, and a dominantly quadratic decay above, similar to (22). In figure 2 we compare numerical simulations of the average purity P (t) (symbols) with the corresponding linear response result (dashed lines) based on (20). The numerical results are obtained from Monte Carlo simulations with 15 different Hamiltonians and 15 different initial conditions for each Hamiltonians. We wish to underline two aspects. First, the energy splitting in general leads to an attenuation of purity decay. Even though a strict inequality only holds for the limiting cases, P F (t) > P D (t) (for t = 0), we may still say that increasing ∆ tends to slow down purity decay. This result is in agreement with earlier findings on the stability of quantum dynamics [26]. Second, for the fast limit and an eigenstate of H 1 (g φ = 1) we find linear decay even beyond the Heisenberg time. A similar behaviour has been obtained in [6], but there an eigenstate of the whole Hamiltonian was required.
In [13] it was shown that exponentiation of the linear response result leads to very good agreement beyond the validity of the original approximation. We use the formula (C.1) where P LR (t) is truncation to second order in λ of the expansion (20), and P ∞ = 1/2 the estimated asymptotic value of purity for t → ∞, see Appendix C. From figure 2 we see that the exponentiation indeed increases the accuracy of the bare linear response approximation.
Decoherence of two qubit systems: A random matrix description   (20) and (25), where P∞ is given in Appendix C. Note that the level splitting ∆ tends to slow down decoherence.

The GOE case
We now drop H 1 leaving H 0 = H e , resulting in where H e ∈ GOE acts on H e and V ∈ GOE on H e ⊗ H 1 . The resulting ensemble of Hamiltonians is invariant under local orthogonal transformations. In the environment, this allows again to diagonalize H e . In the qubit it allows rotations of the kind exp(iασ y ) ∈ O(2). If such transformations are represented on the Bloch sphere, they become rotations around the y axis. Hence, they can take any point on the Bloch sphere onto the xy-plane. Supposing this point represents the initial state, it shows that we may assume |ψ 1 to be of the form In this expression, γ ∈ [−π/2, π/2] denotes the angle of the point representing the initial state with the xz-plane (see figure 3). In order to obtain the linear response expression for P (t) we again make use of (16) and (17). However, apart from the correlation functions used in the GUE case, we have now to consider in addition S ′ 1 (τ ), as given in (B.17). The special case H 1 = 0 can simply be obtained by setting ∆ = 0. This yields After evaluating the double integral in (16), we obtain where B (1) is the double integral of the form factor. It can be computed analytically, but the resulting expression is very involved [5]. For our purpose it is sufficient to note that for t ≪ τ H , B 2 (t) grows only logarithmically.  In figure 4 we show P (t) for γ = 0 (green squares) and for γ = π/2 (blue circles). In contrast to the GUE case (degenerate limit), the average purity depends on the initial state via the angle γ. The fastest decay of purity is observed for γ = π/2, where the image under the time reversal operation becomes orthogonal to the initial state. The slowest decay is observed for γ = 0, which characterizes states which remain unchanged under the time reversal symmetry operation. These statements can be directly translated to the von Neumann entropy S. For one qubit it has a one to one Decoherence of two qubit systems: A random matrix description 11 relation with purity

PSfrag replacements
Observe the entropy scale on figure 4. Due to the dependence of (29) on γ, different initial conditions in the Bloch sphere will yield different behaviours of purity. However for a fixed value of γ we numerically check that there is self averaging (not shown), meaning that different members of the ensemble behave the same as the average, in the large N e limit. As for one qubit the von Neumann entropy is a function of purity, this behaviour will also be observed in this quantity. The consequences of the weaker invariance properties of the GOE will be analysed in detail in a more general framework in a latter paper.

Two qubit decoherence
In this section, we address the question whether entanglement within a given system affects its decoherence, if the system is coupled to an environment. We discuss the three different situations described in section 2.2.

The spectator Hamiltonian
The spectator Hamiltonian discussed in section 2.2, reads In order to study the decoherence of the initial state we use the fact that the dynamics in H 2 also decouples from that in H 1 and H e , i.e. that H 2 commutes with all other terms in the Hamiltonian. The quantum echo of ̺ 0 after time t is Since ̺ M (t) remains a pure state in H 1 ⊗ H 2 ⊗ H e for all times, As the echo operator acts as the identity on the second qubit, where ρ 1 = tr 2 |ψ 12 ψ 12 |. We may therefore compute the purity of the spectator model, without ever referring explicitly to the second qubit. Any dependence of the decay of the purity on the entanglement between the two qubits is encoded into the initial density matrix ρ 1 . This also implies that we can use the results obtained in Appendix A, and hence (16) and (17) remain valid. The only difference is that for the correlation functions C 1 (τ ), S 1 (τ ), and S ′ 1 (τ ), we now have to insert the respective expressions which apply for mixed initial states of the first qubit. These expressions are given in Appendix B.

PSfrag replacements
. We see a plot to help visualize the way the initial condition is parametrized. On the left, qubit one, has an internal Hamiltonian. Its eigenvectors (|0 and |1 ) are represented in blue. The z axis is chosen parallel to the vector |0 . The x axis is chosen in such a way as to make both |0 1 and |1 1 have real coefficients i.e. such that the xz plane contains |0 1 and |1 1 . On the right we represent the second qubit where we have absolute freedom to choose the basis (even if an internal Hamiltonian is present), and thus we choose it according to the natural Schmidt decomposition.

The GUE case:
We again wish to write the initial condition in its simplest form. We must respect the structure of (4), but take advantage of all its invariance properties. Given a fixed H 1 , the ensemble of Hamiltonians is invariant under local operations of the form U Ne ⊗exp iαH 1 ⊗U 2 where U Ne ∈ U(N e ) is any unitary operator acting on the environment, α a real number, and U 2 ∈ U(2) is any unitary operator acting on the second qubit. The freedom in the environment allows to write H e in diagonal from, whereas the freedom within the qubits allows to choose a basis {|0 , |1 } ⊗ {|0 , |1 } in which the state can be written as |ψ 12 = cos θ(cos φ|0 + sin φ|1 )|0 + sin θ(sin φ|0 − cos φ|1 )|1 , (37) and still, To find this basis we start using the Schmidt decomposition to write with {|0 i , |1 i } being an orthonormal basis of particle i. For the first qubit, we fix the z axis of the Bloch sphere (containing both |0 and |1 ) parallel to the eigenvectors of H 1 , and the y axis perpendicular (in the Bloch sphere) to both the z axis and |0 1 . The states contained in the xz plane are then real superpositions of |0 and |1 , which implies that |0 1 = cos φ|0 + sin φ|1 and |1 1 = sin φ|0 − cos φ|1 for some φ. In the second qubit it is enough to set |0 = |0 2 and |1 = |1 2 . This freedom is also related to the fact that purity only depends on tr 2 |ψ 12 ψ 12 |. A visualization of this procedure is found in figure 5. The angle θ ∈ [0, π/4] measures the entanglement (C(|ψ 12 ψ 12 |) = sin 2θ) whereas the angle φ ∈ [0, π/2] is related to an initial magnetization. The general solution for purity using this parametrization is where the geometric factors g (1) θ,φ ∈ [0, 1/2], and g (2) θ,φ ∈ [1/2, 1] are expressed as in terms of the functions g φ and g θ , defined in (21). Both geometric factors are shown in figure 6. Equation (39) is obtained from (16) and (17) by insertion of the (B.10) and (B.16) for Re C 1 (τ ) and S 1 (τ ), respectively.

PSfrag replacements
θ,φ , and g (2) θ,φ from left to right respectively. For g (1) θ,φ we see that for pure eigenstates of H 1 its value is zero. This leads to a higher qualitative stability of this kind of states.
We consider again two limits for ∆. In the degenerate limit (∆ ≪ 1/τ H ) purity decay is given by where f τH (t) is defined in (23). The result is independent of φ since a degenerate Hamiltonian is, in this context, equivalent to no Hamiltonian at all. The θ-dependence in this formula shows that an entangled qubit pair is more susceptible to decoherence than a separable one.
For initial states chosen as eigenstates of H 1 we find linear decay of purity both below and above Heisenberg time. In order for ρ 1 to be an eigenstate of H 1 it must, first of all, be a pure state (in H 1 ). Therefore this behaviour can only occur if θ = 0 or θ = π/2. Apart from that particular case, we observe in both limits, the fast as well as the degenerate limit, the characteristic linear/quadratic behaviour before/after the Heisenberg time similar to the one qubit case. In figure 7 we show numerical simulations for P (t) . We average over 30 different Hamiltonians each probed with 45 different initial conditions. We contrast Bell states (φ = π/4, θ = π/4) with separable states (φ = π/4, θ = 0), and also systems with a large level splitting (∆ = 8) in the first qubit with systems having a degenerate Hamiltonian (∆ = 0). The results presented in this figure show that entanglement generally enhances decoherence. This can be anticipated from figure 6, since for fixed φ, increasing the value of θ (and hence entanglement) increases both g (1) θ,φ and g (2) θ,φ . At the same time we find again that increasing ∆ tends to reduce Decoherence of two qubit systems: A random matrix description  the rate of decoherence, while a strict inequality only holds among the two limiting cases (just as in the one qubit case). From g (2) θ,φ = 2 − g (1) θ,φ − g θ , it follows that (P F − P D )/λ 2 = g (2) θ,φ [f τH (t) − 2tτ H ] ≥ 0. Therefore, for fixed initial conditions and t greater than 0, P F (t) > P D (t). This is the second aspect illustrated in figure 7.
In order to extend the formulae to longer times/smaller purities we exponentiate them using the results of Appendix C. The numerical simulations (see figure 7) agree very well with that heuristic exponentiated linear response formula. In one case (blue rhombus) where the agreement is not so perfect, we found that it is the inaccurate estimate P ∞ = 1/4 of the asymptotic value of purity, which leads to the deviations. 4.1.2. The GOE case: Let us consider the GOE average of both H e and V 1,e . We have to take into consideration the weaker symmetry properties of the GOE ensemble as compared with the GUE. When averaging H e and V 1,e over the GOE, we are again confronted with the fact that the invariance group is considerably smaller than in the GUE case. In this context the initial entanglement between the two qubits has a crucial importance since it "transports" the invariance properties from the spectator to the coupled qubit.
For the sake of simplicity, we focus on the degenerate limit setting H 1 = 0. Note that on the basis of the results in Appendix A, the general case can be treated similarly and the corresponding result will be presented at the end of this subsection.
We first specify the operations under which the spectator Hamiltonian (32), considered as a random matrix ensemble, is invariant. As both, the internal Hamiltonian of the environment and the coupling, are selected from the GOE the invariance operations form the group and have the structure O Ne ⊗ exp(iασ y ) ⊗ U 2 , with O Ne being an orthogonal matrix (acting in H e ), α a real number, and U 2 a unitary operator acting on the spectator qubit.
The direct product structure of the invariance group oblige us to respect the identity of each particle, but allows to analyse each qubit separately. For instance, if we would replace the random coupling matrix V 1,e with one, which involves both qubits, the invariance group would be O(N e ) × O(4). As a consequence, purity decay would become independent of the entanglement within the qubit pair: for any entangled state one can find a orthogonal matrix which maps the state onto a separable one.
We write the initial condition |ψ 12 as in (38). For the coupled particle follows the same analysis made in section 3.2. We can thus write |0 1 = 2 −1/2 [|0 +exp(iγ)|1 ] and, in order to respect orthogonality, For the second qubit we have the same complete freedom as in (4.1.1). We thus select |0 2 = |0 and |1 2 = exp(−iζ)|1 to erase the relative phase in the first qubit and finally write the initial state as The average purity is still given by the double integral expression in (16). However, in the present case the mixed initial state ρ 1 = tr 2 |ψ 12 ψ 12 | must be used. For ∆ = 0, the resulting integrand reads where C 1 (τ ), S 1 (τ ), and S ′ 1 (τ ) are given in (B.10), (B.16), and (B.20), respectively. Evaluating the double integral, we obtain where B (1) 2 (t) is given in (30). As in the GUE case, this result depends on the entanglement of the initial state, and, as in the one-qubit GOE case, it also depends on γ. Again it turns out that Bell states are more susceptible to decoherence than separable states. Note however, that purity as a function of θ is not monotonous. Hence, a finite increase of entanglement does not guarantee that the purity decreases everywhere in time. For separable states, g θ = 1, we retrieve formula (29). However, for completely entangled states, θ = π/4, and the dependence on γ is lost. This is understood from a physical point of view, noticing that any bilocal unital channel (in particular unitary operations) on a Bell state can be reduced to a single local unitary operation acting on a single qubit i.e. for any local unitary operation U 1,2 , there exists (|Bell is any 2 qubit pure state with C = 1, e.g. |00 + |11 ). We can then say that the invariance properties in the second qubit are inherited from the first qubit via entanglement.
Let us obtain the standard deviation for the different possible initial conditions in the qubits. We want to analyse the situation separately for a fixed value of concurrence. Then, as the invariant measure of the ensemble of initial conditions, and fixing the amount of entanglement, we shall use the tensor product of the invariant measures in each of the qubits. Since there is no dependence of (47) on the second qubit, the appropriate invariant measure is trivially inherited from the invariant measure for a single qubit. The resulting value for the standard deviation is Based on Appendix A and Appendix B we can also obtain the average purity for ∆ = 0. The parametrization of the initial states is more complicated since two preferred directions arise, one from the eigenvectors of the internal Hamiltonian and the other from the invariance group. The result can be expressed in the form given in (16), with The angle η is related to a phase shift between the components of any of the eigenvectors of the initial density matrix ρ 1 .

The separate environment Hamiltonian
We proceed to study purity decay with other configurations of the environment. Consider the separate environment configuration, pictured in figure 1(b). The corresponding uncoupled Hamiltonian is and the coupling is From now on we assume that the internal Hamiltonians of the environment and the couplings are chosen from the GUE. Generalization for the GOE can be obtained along the same lines using the corresponding results of section 4.1.2. The initial condition has a separable structure with respect to both environments, see (3). Next we calculate the coupling in the interaction picture. It separates into two parts acting on different subspaces λṼ = λ 1Ṽ (1) + λ 2Ṽ (2) , wherẽ Notice thatṼ (1) (Ṽ (2) ) does not depend on H e ′ (H e ). Since V (1) and V (2) are uncorrelated quadratic averages separate as This leads to a natural separation of each of the contributions to purity where P (i) spec (t) denotes the average purity with particle i being a spectator, as given in section 4.1. In this way, the problem reduces to that of the spectator model. The respective expressions in section 4.1 may be used. For instance, if we assume broken TRI, we obtain from (39) whereC 1 andC 2 are the correlation functions of the corresponding environments defined in exact correspondence with (18), for H e and H e ′ respectively. If in one or both of the qubits, the level splitting in the internal Hamiltonians is very large/small compared to the Heisenberg time in the corresponding environment (denoted by τ e and τ e ′ for H e and H e ′ , respectively) the degenerate and/or fast approximations may be used. As an example, if ∆ 1 ≪ 1/τ e and ∆ 2 ≪ 1/τ e ′ we find whereas if ∆ 1 ≫ 1/τ e and ∆ 2 ≫ 1/τ e ′ P F (t) = 1 − λ 2 1 g θ,φ1 f τe (t) + 2τ e g (2) θ,φ1 t − λ 2 2 g (1) θ,φ2 f τ e ′ (t) + 2τ e ′ g (2) θ,φ2 t .
It is interesting to note that if we have two separate but equivalent environments (i.e. both Heisenberg times are equal), we get exactly the same result as for a single environment. Also notice that the Hamiltonian of the entire system separates and thus the total entanglement of the two subsystems (H 1 ⊗ H e and H 2 ⊗ H e ′ ) becomes time independent.

The joint environment configuration
The last configuration we shall consider is the one of joint environment; see figure 1(c).
Note the slight difference with (53). However still V (1) and V (2) are uncorrelated, enabling us to write again (54). From now on, the calculation is formally the same as in the separate environment case. Hence we can inherit the result (56) directly, taking into account that since they come from the same environmental Hamiltonian, the two correlation functions are the same. In case any of the Hamiltonians fulfils the fast or degenerate limit conditions, the corresponding expressions to (57) and (58) can be written. As an example, if the first qubit has no internal Hamiltonian, and the second one has a big energy difference, the resulting expression for purity decay is where τ H is the Heisenberg time of the joint environment. Monte Carlo simulations showing the validity of the result were done with satisfactory results, comparable to those obtained in figure 7. The parameter range checked was similar to that in the figure.

A relation between concurrence and purity
In the previous section we have found that purity decay is very sensitive to the initial degree of entanglement between the two qubits. In other words, we related the behaviour of purity in time to the degree of entanglement (measured by concurrence) at t = 0. In this section, we study the relation between purity and entanglement (again measured by concurrence) as both are evolving in time. We focus on initial states which are Bell states such that purity and concurrence are both maximal and equal to one at t = 0. Since concurrence is defined in terms of the eigenvalues of a hermitian 4 × 4-matrix, an analytical treatment even in linear response approximation is much more involved than in the case of purity. Instead, we work with a phenomenological relation between purity decay and concurrence discovered in [27] and further studied in [13]. It proves to be valid in a wide parameter range, and it allows to obtain an analytic prediction for concurrence decay.
We study the relation between concurrence and purity on the CP -plane, where we plot concurrence against purity with time as a parameter. Since the initial state is a Bell state, the starting point is always at (C, P ) t=0 = (1, 1). In this plane not all points are allowed for physical states; two constrains appear. The first one follows from the ranges of C ∈ [0, 1] and P ∈ [1/4, 1]. The second restriction follows from the "monogamy" of entanglement: if a qubit is very entangled with another qubit, it cannot be very entangled with the rest of the universe (i.e. the environment), implying some degree of purity. Conversely if the qubit is very entangled with the universe, it cannot be entangled with the other qubit. This statement can be made quantitative, which leads to the concept of maximally entangled mixed states. They define the upper bound of the set of admissible states (the corresponding area on the CP -plane is shown in figure 8(a) in grey). Another region of interest corresponds to those states, which form the image of a Bell state under the set of bi-local unital operations [28] (red area in the same figure). Finally, we have the Werner states ρ W = α 1 1 4 + (1 − α)|Bell Bell|, 0 ≤ α ≤ 1, which define a smooth curve on the CP -plane (black solid line). The analytic form of this curve is [13] C W (P ) = max 0, and will be referred to as the Werner curve. Note that states belonging to the Werner curve are not necessarily Werner states. The spectator and separate environment configuration are strictly bi-local, whereas the joint environment configuration becomes bi-local in the large N e limit. None of the three types of environmental couplings are strictly unital. In this limit we find from our numerical simulations that all CP -curves fall into the region defined by the bi-local unital operations.
For figure 8(b), we perform numerical simulations of the spectator Hamiltonian assuming broken time reversal symmetry (GUE case). We compute the average purity P (t) as well as the average concurrence C(t) as a function of time. This average is over 15 different Hamiltonians and 15 different initial states for each Hamiltonian. We fix the level splitting in the coupled qubit to ∆ = 1 and consider two different values λ = 0.02 and 0.14 for the coupling to the environment. figure 8(b) shows the resulting CP -curves for different dimensions of H e . Observe that for both values of λ, the curves converge to a certain limiting curve inside the red area (defined by the set of bi-local unital channels) as dim(H e ) tends to infinity. While for λ = 0.02, this curve is at a finite distance of the Werner curve, for λ = 0.14 it practically coincides with C W (P ). To check this statement in more detail, consider the numerical CP -curve C num (P ) obtained from our simulations and define its distance E to the Werner curve as The behaviour of E is shown in figure 9(a). For λ = 0.14 (black dots), the error goes to zero, which means that the corresponding CP -curves indeed converge to the Werner curve. In fact from a comparison with the black solid line we may conclude that the deviation E is inversely proportional to the dimension of H e . By contrast, for λ = 0.02 (red dots), the error tends to an approximately constant value, in line with the assertion that the numerical results converge to a different curve. In figure 9(b) we plot the error E as a function of λ, for different dimensions of H e as indicated in the figure legend. The results suggest that the convergence to the Werner curve occurs as long as λ > λ c which is of the order of 0.1 for ∆ = 1. In fact for all 0.001 ≤ ∆ ≤ 100 studied, we found a critical value λ c such that the CP -curves converge to the Werner curve as long as λ λ c . For ∆ = 0 all studied couplings led to convergence to the Werner curve.
In the presence of TRI the CP -curves again converge to the Werner curve, for values of the coupling greater than a critical value λ c . As in the GUE case, this critical value depends on ∆. However, in contrast to the GUE case, we find that λ c remains finite for ∆ = 0. In the other configurations considered (the joint and the separate environment), the behaviour is similar. In those cases it is the largest (of the two) coupling strengths which determines whether the CP -curves converge to the Werner curve, or not.
A partial explanation for this particular behaviour in the CP plane is provided in [28]. The authors find that for unital channels, the possible points in the CP plane that can be accessed are limited to a small 2-dimensional region, see figure 8(a).
We test the unitality of the channels induced by the non-unitary evolution in the  Figure 9. In (a) we show E defined in (64) which measures the distance between the numerical CP -curves and the Werner curve as a function of dim Ne. The same data is used as in figure 8. The red dots refer to the case λ = 0.02 where the error E apparently remains finite. The black dots refer to the case λ = 0.14, where E seems to tend to zero. The black solid line, which is proportional to N −1 e , is meant as a guide to the eye. In (b) we plot E as a function of the coupling λ, for various values of Ne. We may say that the line λ ≈ 0.1 separates two regimes. For λ > 0.1 we observe the convergence to the Werner curve, whereas for λ < 0.1 the limiting curve is a different one. PSfrag replacements Figure 10. We evaluate the unitality condition for the spectator configuration (i.e. on one of the qubits). We vary the size of the environment and test different times, that would lead to different values of concurrence as shown in the inset, for the spectator case with the corresponding coupling. A line with slope −1/2 is also included for comparison.
following way. Since we want to test E(1 1) = 1 1, we prepare an initial pure condition that leads to a completely mixed state in the qubit (let us focus on the spectator case), i.e. , with ψ (i) e |ψ (j) e = δ ij . We let the state evolve with a particular member of the ensemble of Hamiltonians defined in (1). Afterwards we evaluate the Euclidean distance of the resulting mixed state in the qubit from the fully mixed state, within the Bloch sphere. The average distance is plotted as a function of the size of the environment in figure 10. We conclude that the unitality condition is approached algebraically fast as the size of the environment increases. This by no means explains the generic accumulation to the Werner curve, but certainly restricts the possibilities; furthermore, for high purities/concurrences the area converges to the curve C = P , hence providing a satisfactory explanation in this regime. Sufficiently close to P = 1, the above arguments imply a one to one correspondence between purity and concurrence, which simply reads C ≈ P . This allows to write an approximate expression for the behaviour of concurrence as a function of time using the appropriate linear response result for the purity decay. The corresponding expressions have been discussed in detail in the previous section. (66) has similar limits of validity as the linear response result for the purity. Thus, in a sense, we may call it a linear response expression for concurrence decay. In those cases, where the coupling is beyond the critical coupling λ c and where the exponentiated linear response expression (25) holds for the average purity, we can write down a phenomenological formula for concurrence decay, which is valid over the whole range of the decay In figure 11 we show random matrix simulations for concurrence decay in the joint environment configuration. We consider small couplings λ 1 = λ 2 = 0.01 which lead to the Gaussian regime for purity decay, as well as strong couplings λ 1 = λ 2 = 0.1 which Decoherence of two qubit systems: A random matrix description 22 lead to the Fermi golden rule regime. We find good agreement with the prediction (67), except for the Gaussian regime when we switch-on a fast internal dynamics in both qubits (∆ 1 = ∆ 2 = 10). These level splitting are so large, that the coupling strengths are much smaller than λ c , which leads to the deviations observed. Note that, in the separate environment case, the entanglement of the two subsystems defined on the spaces H 1 ⊗ H e and H 2 ⊗ H e ′ is constant in time. For the initial conditions we use, the entanglement stems entirely from the two qubits. The concurrence of these two qubits decays because the constant entanglement spreads over all constituents of the two large subspaces.

Conclusions and outlook
We have derived general linear response formulae for the decay of purity for one and two qubits. The initial states are grouped in relevant families and the dependence on these families is analysed in detail. In particular we see that time reversal invariance will generally increase the number of families we have to consider, though entanglement will undo some of this in the two-qubit case. We do not discuss the well-known quadratic decay at very short (Zeno) times, though it obviously occurs in our model [6]. We find the linear decay resulting from the usual master equation treatment at short times, but we see the crossover to quadratic decay as we approach the Heisenberg time. This crossover is important if decoherence is dominated by a few degrees of freedom coupled to the central system, leading to a finite Heisenberg time in the relevant environment. We show that the quadratic term is absent if we choose an eigenstate of the intrinsic Hamiltonian of the central system.
In general we can say that entanglement between the qubits accelerates decoherence mainly by affecting the prefactors in the linear response expression. Internal dynamics introduces an additional term but also affects the prefactors. It tends to stabilize the system though the interplay with entanglement does not make this statement strict for small changes.
The analysis of concurrence decay and its relation to purity decay largely confirms previous findings [13,27] in this more general setting. In particular the time-evolution follows the purity -concurrence relation of Werner states. We use this relation to give an analytic, albeit heuristic, formula for concurrence decay. Actually the relation shows small deviations for very weak couplings, and small concurrences. We obtain some understanding of these findings by showing that the process becomes unital in the limit of a large environment.
We thus have given a rather complete description of the time-evolution of decoherence and entanglement of two non-interacting qubits. From a technical point of view we find a remarkable sequence of argumentation. In linear response approximation, all two-qubit configurations can be reduced to the spectator case. This case is equivalent to a single qubit, initially in a mixed state yet not entangled with the environment. The evolution of the mixed central system is obtained in the appendix in the linear response approximation. This result in turn allows to derive the general formulae for all cases mentioned above. While this allowed a simple unified representation of the two-qubit problem in the framework of an RMT model, actually one readily suspects, that this result is more general, and this is indeed the case as we shall see in a forthcoming paper [29]. Another line, that could be easily followed, is to consider a random interaction between the two qubits and a joint coupling to the environment. While such a study is easy and has some interest, the real problem to tackle is to introduce time dependent operations such as gates, while the decoherence modelled by RMT acts. Another open problem that can be tackled in this framework is a situation of a near and far environment the first coupled weakly and the second extremely weakly. The first should probably have a finite Heisenberg time on a time scale relevant to the system, while the second more likely would have very dense spectrum and thus an infinite Heisenberg time. Finally, recall that we used formally echo dynamics for purely technical reasons. Yet this could be used to great advantage if we simultaneously want to study internal errors in the implementation of a quantum algorithm on a physical device as well as decoherence effects. As RMT is extremely successful in describing fidelity decay it seems to be the optimal approach. linear operators. For arbitrary operators A, B acting both on H ′ , the purity form has the following property: The purity defined in (7) can be expressed in terms of the purity form as Note that on the RHS of this equation, we first take the trace over the central system, i.e. we consider the purity of the state of the environment after tracing out the central degrees of freedom. With this little twist we can describe the one-qubit decoherence in section 3 and the spectator model in section 4.1 at the same time.
Purity decay in the linear response approximation To average purity, we use (11) to compute ̺ M (t) ⊗ ̺ M (t) in linear response approximation Keeping terms only up to second order in λ: In the next step, we perform the ensemble average over the coupling. As a result, the terms linear in λ vanish. For the remaining terms, we use the properties in (A.2) and the fact that I(t) is hermitian, to obtain In order to pull out the same double integral of A J and A I , we use the time ordering symbol T . It allows to write for the average purity as a function of time whereṼ andṼ ′ are short forms for the coupling matricesṼ (τ ) andṼ (τ ′ ), respectively. The arguments τ and τ ′ of the coupling matrices are interchanged in the second and fourth term of A JI . This does not change the value of the integral of course, but it facilitates to handle common terms in the following considerations.
The general result If the coupling matrix is taken either from the Gaussian unitary (GUE) or orthogonal (GOE) ensemble, we find . (A.9) These expressions are derived in the following two subsections, for the GUE case in Appendix A.1, for the GOE case in Appendix A.2. The correlation functions C x , S x , C ′ x and S ′ x with x ∈ {1, e} are defined and discussed in Appendix B.
The limit dim(H e ) → ∞: If the dimension of the Hilbert space of the environment goes to infinity, the corresponding correlation functions simplify, as discussed in Appendix B. We are then left with The integrand A JI in (A.8) consists of four terms. These terms will be considered, one after the other. We will always first average the argument of the purity form.
Only then we perform the partial traces over subsystem H 1 , and therefore the final trace over the environment. The averaging is only over the coupling, and it is done by applying two simple rules, as described below.
For the coupling V 1,e in the product eigenbasis of H 0 , we use either random Gaussian unitary (GUE) or orthogonal (GOE) matrices. Their statistical properties are completely characterized by the following second moments (for notational ease we ignore the subscript 1,e for a moment) where χ GOE = 1 if V is taken from the GOE, while χ GOE = 0 if V is taken from the GUE. In the interaction picture we then find: We first compute the average where we have used the averaging rule, (A.12), as well as the fact that the timeordering operator T requests to exchange τ with τ ′ whenever τ < τ ′ . The indices i, k and m denote basis states in H 1 , while the indices j, l and n denote basis states in H e . We can rewrite that expression in a more compact form by employing the diagonal matrices C x (τ ), defined in (B.1): T ṼṼ ′ ̺ 0 ⊗̺ 0 = C 1 (|τ −τ ′ |)⊗C e (|τ −τ ′ |) ρ 1 ⊗|ψ e ψ e | ⊗ ρ 1 ⊗|ψ e ψ e | .(A.14) Now we apply the partial traces over subsystem H 1 The final trace over the environment yields We first compute the average where we have used (A.12). We apply the partial traces over subsystem H 1 The final trace over the environment yields We first compute the average Then we apply the partial traces over subsystem H 1 (A.21) The final trace over the environment yields We first compute the average We apply the partial traces over subsystem H 1 (A.24) The final trace over the environment yields In the GOE case, the average over the coupling yields, besides the previously considered GUE-term an additional one; see (A.12). In the following we will redo the calculation of the previous subsection, but consider only that additional term. As a reminder, we add the subscript 2 to the brackets which denote the ensemble average.
where we have used the averaging rule, (A.12). Here the time-ordering operator has no effect. Now we apply the partial traces over subsystem H 1 The final trace over the environment yields We first compute the average where we have used (A.12). We apply the partial traces over subsystem H 1 The final trace over the environment yields l|ψ e e iE l (τ +τ ′ ) l|ψ e ψ e |j e −iEj (τ +τ ′ ) ψ e |j We first compute the average Then we apply the partial traces over subsystem H 1 (A.33) The final trace over the environment yields We first compute the average We apply the partial traces over subsystem H 1 (A.36) The final trace over the environment yields 2 (t) is the two-point spectral form factor with time measured in units of the Heisenberg time τ H for the corresponding spectral ensemble. In the limit of large dimension N e = dim(H e ) we find for the other correlation functions: S e (τ ) = tr e −iHeτ |ψ e ψ e |e iHeτ |ψ e ψ e | (B.6) S ′ e (τ ) = tr e −iHeτ |ψ e ψ e |e iHeτ |ψ * e ψ * e | . (B.7) For random pure states, as considered in the present work, both correlation functions are at most of order one at τ = 0. For τ > 0 they drop very quickly and soon become of order N −1 e . This happens on the same time scale, where C e (τ ) drops from values of the order N e to values of the order one. In that sense we consider these correlation functions to contribute only O(N −1 e ) corrections to the result given in (A.10).

Appendix C. The exponentiation
The linear response formulae derived throughout the article are expected to work for high purities or, equivalently, when the first two terms in the series (11) approximate the echo operator well; it is of obvious interest to extend them to cover a larger range. The extension of fidelity linear response formulae has, in some cases, been done with some effort using super-symmetry techniques. This has been possible partly due to the simple structure of fidelity, but trying to use this approach for a more complicated object such as purity seems to be out of reach for the time being. Exponentiating the formulae obtained from the linear response formalism has proven to be in good agreement with the exact super-symmetric and/or numerical results for fidelity, if the perturbation is not too big. The exponentiation of the linear response formulae for purity can be compared with Monte Carlo simulations in order to prove its validity. We wish to explain the details required to implement this procedure in this appendix.
Given a linear response formula P LR (t) (for which P (0) = 1), and an expected asymptotic value for infinite time P ∞ , the exponentiation reads This particular form guaranties that P ELR equals P LR for short times, and that lim t→∞ P ELR (t) = P ∞ The particular value of P ∞ will depend on the physical situation; in our case it will depend on the configuration and on the initial conditions. Let us write the initial condition as |ψ 12 (θ) = cos θ|0 102 + sin θ|1 112 (C.2) for some rotated qubits |0 i , |1 i . In each of the qubits interacting with the environment we assumed that for long enough times, the totally depolarizing channel E d will act (recall that the totally depolarizing channel is defined as E d [ρ] = 1 1/ tr 1 1 for any density matrix ρ). Hence, for the spectator configuration, the final value of purity is whereas for both the joint and the separate environment configuration we use the value Good agreement is found with Monte Carlo simulations for moderate and strong couplings.