# Mean-field evolution of open quantum systems: an exactly solvable model

M. Merkli, G. P. Berman

## Abstract

We consider quantum particles coupled to local and collective thermal quantum environments. The coupling is energy conserving, and the collective coupling is scaled in the mean-field way. There is no direct interaction between the particles. We show that an initially factorized state of the particles remains factorized at all times, in the limit of large particle number. Each single-particle factor evolves according to an explicit, nonlinear, dissipative and time-dependent Hartree–Lindblad equation. The model is exactly solvable; we do not make any weak coupling or any Markovian approximations, and our results are mathematically rigorous.

## 1. Introduction and results

Studying the effects of noise on quantum systems is of central importance in modern quantum theory. Quantum information processing is based on manipulation of superposition and entanglement of basic quantum bits forming a quantum processor. As any such system is subject to noise, e.g. due to contact with thermal environments, it is crucial to understand the mechanisms of noise effects, and to quantify them mathematically. It is known that interactions with environments generically destroy phase coherence (decoherence), as well as quantum correlations (entanglement). The literature on this subject is huge; we refer only to Yu & Eberly (2004, 2009), Merkli (2011) and Merkli et al. (2011a) and references therein. Another (more positive) aspect is that, starting with disentangled states, one can create and control entanglement by coupling quantum systems to a common thermal noise (Beige et al. 2000; Basharov 2002; Braun 2002; Jakóbczyk 2002; Benatti et al. 2003; Ficek & Tanaś 2003, 2008; Benatti & Floreanini 2005; Paternostro et al. 2006; Zhang & Yu 2007; Cole 2010; Merkli et al. 2012). This opens up the possibility of manipulating quantum bits by controlling their surroundings. In the works cited earlier, the focus is on the entanglement of two subsystems interacting indirectly via a joint quantum thermal noise reservoir. (See Merkli et al. (2012), for a numerical investigation of arbitrary numbers of subsystems, where it is shown that creation of two-spin entanglement is suppressed if N is large.)

In this work, we consider entanglement creation in complex open quantum systems, i.e. for a large number N of particles interacting indirectly via a common thermal environment. As the number of particles increases, the total energy of the system increases (linearly in N). It is plausible that the interaction to the common heat bath creates an effective interaction between particles. The size of the interaction energy—for instance, an induced two-body interaction—will be proportional to the number of pairs of particles, hence of order N2. Therefore, in order to have a balanced competition between individual energy and interaction energy, one scales the interaction with a negative power of N. As we show, there is exactly one ‘correct’ scaling, having the property that the dynamics is non-trivial in the limit of infinitely many particles. In many-body quantum physics, an N-dependent scaling of particle interactions is called a mean-field scaling. So far, mean-field models of directly interacting particles (coupled via two-body potentials) have been considered. For these systems, it is known that in the mean-field picture, each particle feels the same averaged effect of all other particles (Spohn 1980). Each particle thus evolves independently, although according to a more complicated dynamics containing interaction effects. The equation governing the evolution of the state of a single particle turns out to be nonlinear in the state and is called the nonlinear Hartree equation. The emerging picture is that high complexity (large numbers N of particles) favours particle independence, but complicates the individual dynamics.

In this work, we use ideas from mean-field theory in the context of indirect interactions in complex open quantum systems. Here, a large number N of particles do not interact directly, but via a common thermal quantum noise. One motivation to consider such interactions is the question of whether a common noise can create entanglement in a particle system, the answer to which is shown here to be negative if N is large. We are not aware of any previous work considering this type of interaction. We show that for disentangled initial states, and for all times, the state of the particles is very close to a product state, the difference being of the order 1/N. We identify the dynamical equation of a single particle. It is a time-dependent effective evolution equation containing nonlinear Hartree terms as well as dissipative Lindblad contributions. The dynamics is not Markovian, as is visible in our dynamical equations via the dependence on time, of both the nonlinear and the dissipative contributions. Our model has purely energy-conserving (purely dephasing, non-demolition) particle–reservoir interactions. As a consequence, we can solve the dynamical equations explicitly and obtain exact results, for all times and for all strengths of interaction between particles and reservoirs. This kind of interactions are physically justified for systems (or time spans) in which decoherence effects dominate relaxation effects (see also references cited earlier). A rigorous treatment of systems including energy-exchange interactions between particles and reservoirs is technically much more involved. We are going to address this problem elsewhere.

### (a) Explanation of main result

We consider N quantum particles interacting with local and collective thermal environments. Each particle is in contact with its own, independent environment and coupled to another one, common to all particles. There is no direct interaction, but the collective environment induces interaction among all particles. Our main result is the following: starting in a factorized (product) initial state for N particles, ρ0⊗⋯⊗ρ0, the reduced density matrix for any n particles (nN fixed), at any given moment in time t, approaches a product state as . We summarize this as where the single-particle density matrix ρt is the solution of the equation See theorem 1.1 for a precise statement. The first commutator term represents the free dynamics, A being the Hamiltonian of a single particle. The trace term (the trace is taken on the second factor of the space of two particles) is due to the indirect interaction of particles via the common environment. It is called the Hartree term and is quadratic in the density matrix. The operator Weff(t) is an effective two-body interaction. The interaction with the collective thermal quantum noise creates k-body interactions in the particle system, for all k≥2. The lowest order is a pair interaction, and higher-order corrections to the Hartree–Lindblad equation, which are proportional to powers in 1/N, involve many-body effective interactions (we are going to explain the details of this in a subsequent work). The last term in the evolution equation is a time-dependent Lindblad term, describing the effect of the local environment. Both the effective two-body potential and the Lindblad operator are given by explicit time-dependent operators.

The convergence to a product state implies that complexity (large N) prevents the particles from ever becoming entangled, even though they interact via a common environment. We give an estimate on the difference between the reduced n-body density matrix and the limit product state in theorem 1.5, provided N is large enough, and with an explicit Cn(t). An upper bound is Cn(t)≤{2(1+t+{2t)en[1+n{2t(1+{2t)], where is the particle–reservoir interaction strength (see, however, theorem 1.5 for a better bound).

### (b) Relation to previous work

There are previous works on mean-field open quantum systems, see Breuer & Petruccione (2003) and references therein. However, in all of them, the dynamics of the N particles is considered to be given by a master equation. By using a master equation, one implicitly assumes that a Markovian (weak coupling, van Hove) limit of the dynamics has been taken. The mean-field limit dynamics, originating from the master equation in the limit is thus not guaranteed to approximate the true dynamics (which is not Markovian). The technical issue is that of an interchange of two limits: that of a large number of particles and that of small coupling. In this work, we do not encounter this problem, as we are able to solve the dynamical equations exactly. As a consequence, the single-particle dynamics we derive exhibits the correct non-Markovian behaviour, expressed through a time dependence of the nonlinear and dissipative parts in the equation. This feature of the dynamics cannot be seen in an analysis based on master equations, as the mean-field limit of a Markovian dynamics will still be Markovian.

Another difference between our and previous works is the model itself. In previous works, particles are taken to interact directly (via pair-potentials). In contrast, we consider here particles that do not interact directly, but only indirectly, via a coupling to a common noise (collective quantum heat bath reservoir). One of our motivations for considering such an interaction is the question of whether a common noise can produce entanglement in a particle system. We show here that entanglement creation is disabled by complexity (in the mean-field limit).

### (c) Mathematical description

The Hilbert space of pure states of each particle is ℋ, with The Hamiltonian of particle j is denoted by Aj. To each particle is associated a local environment, and there is another, collective environment. Each of them is given by a spatially infinitely extended free Bose gas initially in thermal equilibrium (a heat bath). Without further effort, one could consider individual temperatures for each reservoir, but we assume they are all at a common temperature T=1/β>0. The free-field Hamiltonian of a reservoir is given by acting on the (symmetric) Fock space over the one-Boson space , (momentum representation, where |k| is the energy of a Boson). In some physics literature, one would write and understand that an infinite-volume (continuous momentum) limit is performed. The creation and annihilation operators satisfy the usual canonical commutation relations [a(k),a*(l)]=δ(kl). The equilibrium state of a reservoir is determined by 〈a*(k)a(l)〉β=δ(kl)/(eβ|k|−1) (plus Wick's rule). We refer to Attal et al. (2006) for more details on the mathematical description of reservoirs. Let be a ‘form factor’. We define the field operator where and a(f) is the adjoint of a*(f). Particle j interacts with its local reservoir through an interaction {jV jφj(fj), and to a collective reservoir via {Wjφ(f). The operator φj(fj) acts on the jth reservoir as φ(fj). We have introduced coupling strengths { and {j that are arbitrary real numbers (not necessarily small). The Hilbert space of the total system is with N copies of the single-particle space ℋ and N+1 copies (N local plus one collective) of the single reservoir space ℱ (Fock space). The total Hamiltonian is 1.1 1.2 1.3 Each of the Kj is the same operator K, but acting on the space of the jth reservoir. The collective interaction is homogeneous, in that Wj is a fixed operator W, acting on the space of the jth particle. It is assumed that the interactions are purely dephasing (energy conserving, of non-demolition type), in the sense that they are diagonalized jointly with the free Hamiltonians Aj. This kind of interaction produces phase processes, such as decoherence and entanglement (Palma et al. 1996). It does not describe relaxation processes (populations, or the diagonal of the density matrix in the energy basis, are constant).

The time-dependent reduced n-body density matrix is 1.4 Here, ρ0 is the initial single-particle density matrix and ρR is the product state of all N local reservoirs and the collective one, each in its own equilibrium (at a fixed temperature 1/β). The symbol Tr[n+1,N] means that we take the trace over all degrees of freedom of particles n+1,n+2,…,N, and over all reservoirs.

The following two functions are key dynamical quantities: 1.5 and 1.6 We denote by Sj(t) the quantity (1.5) with f replaced by fj, and similarly we define Γj(t).

### Theorem 1.1 (Convergence to Hartree–Lindblad dynamics)

For any and n≥1, we have The single-particle density matrix ρj,t satisfies the time-dependent Hartree–Lindblad equation 1.7 where the effective two-particle interaction and the Lindblad operator are

### Remark 1.2

The coupling constants ϰ and ϰj can take arbitrary values in (they do not have to be small). In the outline of our results at the beginning of the introduction, we have taken them to be equal to one, for simplicity of the explanations.

### Remark 1.3

Equation (1.7) has three terms on the right-hand side. The first one originates from the free individual dynamics of spin j. The second one describes the effect of all spins on spin j, through an effective interaction operator Weff, and is quadratic in the density matrix. This effective interaction is due entirely to indirect coupling of the spins, mediated by a common heat bath. The third term is due to the local reservoirs. Note that the noise effects are independent, in the sense that each term on the right-hand side of (1.7) can be switched on and off individually by its own coupling constant.

### Remark 1.4

Equation (1.7) is time dependent. This means that the dynamics is not Markovian, a well-known property of open quantum systems.

The convergence of the density matrix as is given in the ‘trace norm’, that is, the trace of the absolute value of an operator. The absolute value of an operator is defined by . The choice of the trace norm is motivated as follows. Suppose that A is an operator acting on n particles, where nN is a fixed number. By theorem 1.1, we have as . We use here that |TrAB|≤∥A∥ Tr|B|, where ∥⋅∥ is the operator norm, (the supremum is taken over all normalized vectors in the Hilbert space on which the operator A acts). This means that for large N, the expectation of the observable A is obtained as an average in a product state. In particular, (any measure of) entanglement between arbitrary n particles vanishes for all times, in the limit of large complexity. The following result estimates how quickly this happens as complexity grows.

### Theorem 1.5 (Convergence speed)

Let η=4ϰ2n∥W∥2 and recall that d is the dimension of the single-particle Hilbert space. For and 1≤n≤N, such that N> η|S(t)|, we have where Cn(t)=ηd2n((2n+2η|S(t)|+1)eη|S(t)|(n+2η|S(t)|+1)+n|Γ(t)|).

For large times, we have and , where is the (collective) reservoir spectral density, see Merkli et al. (2008a) and Merkli et al. (2011b).

## 2. Illustration: N spins (qubits)

We consider (1.1)–(1.3) for N spins, each one having pure state space and each Hamiltonian Aj being equal to ω>0 being the frequency of a spin (in units of ). The noise interaction operators are the same for each spin, i.e. V j=V and Wj=W, with The coupling constants are taken to be homogeneous, ϰjl (local) for all j and the collective coupling constant is denoted ϰ=ϰc. We assume the local form factors in (1.2) to be all equal, fj=fl and denote the collective one in (1.3) by f=fc. According to theorem 1.1, the single-spin density matrix evolves according to the Hartree–Lindblad equation (1.7), 2.1

This equation is the same for each spin j; hence, we do not write this index. Here, Sc(t) and Γl(t) are the quantities (1.5) and (1.6) with f the collective and local form factors fc and fl, respectively. Denote by [ρt]m,n the density matrix of ρt in the ordered energy basis ϰφ1,φ2}, where Szφ1=(1/2)φ1, Szφ2=−(1/2)φ2. Let ρ0 be the initial condition ρt|t=0 and denote its populations by [ρ0]1,1=p∈[0,1], [ρ0]2,2=1−p. The solution to (2.1) is given by (see also (3.11)) 2.2 and 2.3 The populations are constant because all interactions commute with the single-particle Hamiltonian A. The single-particle off-diagonal evolution (2.3) has two oscillatory factors and a decaying one. The first oscillatory factor is due to the free evolution of a single spin, the other one is induced by the collective interaction and depends on the initial state. The relaxation or decay term originates purely from the local coupling. One can view the collection of all spins as some environment for each fixed single spin. Equation (2.3) shows that the effect of such an environment is entirely different from the usual one, induced by thermal noises (an infinite Bose gas initially in equilibrium). Indeed, the ‘other spin environment’ only contributes with a phase factor to the dynamics of the off-diagonal single-particle density matrix elements, whereas the thermal reservoir generates a decaying factor . We phrase this as follows.

In the mean-field limit, the cumulative effect of all spins on a fixed single one gives a modulation of the dephasing process, it does not accelerate decoherence. Decoherence is driven entirely by the coupling to local reservoirs.

One can view the spins as magnetic entities. Denote by Sx,Sy and Sz the three spatial components of a single spin (spin operators, proportional to the Pauli matrices). The Hamiltonian A introduced at the beginning of this section is the energy of a spin placed in a magnetic field in the z-direction, of magnitude Bz=ω/γ, where γ is the gyromagnetic ratio. Define the transversal spin operator by The average of S at time t is precisely Tr(ρtS)=[ρt]2,1, the complex conjugate of (2.3). It follows that the effect of the collective coupling, in the mean-field way, is directly visible in the angular speed of precession of the spins around the z-axis. When the spins are coupled merely to their local reservoirs, this speed is simply ω. However, the effect of the presence of a large number of other spins, coupled indirectly via the collective reservoir, and in the mean-field scaling limit, becomes time dependent and is given by Note that when Sc(t) approaches its linear limit at, the precession speed becomes .

## 3. Proof of theorems 1.1 and 1.5

### (a) Extraction of the main term

There are rank-one projections , acting on the jth particle space, such that 3.1 3.2 and 3.3 where , and w(m) are the (possibly repeated) eigenvalues of the operators Aj, V j and W, respectively.

Let and let pk be the population of level k in ρ0, i.e. . We have 3.4 where We simplify the trace in (3.4). By independence of the reservoirs, we have and similarly for the second exponential in the trace in (3.4). That trace thus becomes the product 3.5 where 〈Xβ is the average of X in the Bosonic equilibrium state at temperature 1/β. One has explicitly (Merkli et al. 2008b), for , 3.6 where S(t) and Γ(t) are given in (1.5) and (1.6). Using (3.5) and (3.6), we obtain 3.7 We need to insert this expression into (3.4). The result is 3.8 We separate (3.8) into a main term and a remainder, as , as follows. Define 3.9 Note that this term is simply (3.8) with ‘N replaced by ’. In this limit, of course , and, by using that and the expansion , see proposition 3.1 for details, we obtain that in the limit , the term [⋯ ]Nn becomes the second last exponential in (3.9). We thus have, for n and t fixed, 3.10 where R:=ρn,N(t)−ρn,N(t) is the reminder term, an operator depending on N,n,t. We estimate the size (trace norm) of R in §3b. Here, we continue the analysis of the main term .

We write the second last exponential in (3.9) as the product and noting that we see that the right-hand side of (3.9) is of the product form ρ1,t⊗⋯⊗ρn,t, where 3.11 We take the derivative with respect to time of the last equation to obtain 3.12 The sum over m is where the last equality holds because (d/dt)Tr(ρj,tW)=0 (use, for example, (3.12) to see this). It follows that 3.13 The trace is taken over the second space. Combining (3.12) and (3.13) yields the evolution equation in theorem 1.1.

### (b) Control of convergence speed

We now investigate the speed of convergence of R, as . From the definition of R, (3.10), it follows that 3.14 where and, setting 3.15 3.16 and 3.17 Note that |e|≤1. R is an operator on the n-fold tensor product ℋn:=ℋ⊗⋯⊗ℋ, where ℋ is the Hilbert space of a single particle. Because the space of bounded linear operators on ℋn, denoted by ℬ(ℋn) (and having the usual operator norm ∥⋅∥), is the dual space of the Banach space L1(ℋn) of trace-class operators on ℋn (with norm ∥x1=Tr|x|), we have 3.18 Let B be a bounded operator on ℋn. By cyclicity of the trace, we have 3.19 Since |TrXY |≤∥Y ∥Tr|X|, and because a density matrix has trace one, we obtain from (3.19) the bound 3.20 where .

We estimate the supremum of |f2|. For x real, we have . Similarly, for x≥0, we have . Furthermore, the last term [⋯ ]Nn in the expression (3.17) for f2 has modulus less than or equal to 1 (because ). We obtain 3.21 In order to estimate |f1|, we establish the following result.

### Proposition 3.1

Let xm be as in (3.15) and set ξ=4ϰ2nW2|S(t)|. Suppose that N>2ξ. Then, we have where satisfies |R′|≤(ξ/N)(n+2ξ+1) e(ξ/N)(n+2ξ+1).

We give a proof of proposition 3.1 below. Combining the result of the proposition with (3.21) gives where η=4ϰ2nW2. Using this bound in (3.20), we arrive at 3.22 This is the bound given in the theorem.

### Proof of proposition 3.1.

The left-hand side equals . We have 3.23 where Because xm is real, we have (as mentioned earlier) |eixm/N−1|≤|xm|/N, so that (by the definition of ξ and the condition on N in the proposition). We can thus use the Taylor expansion , valid for |a|<1 in (3.23), . The modulus of the last sum is bounded by the value of the geometric series, (1−|a|)−1<2. It follows that , with |R1|<2. Therefore, 3.24 with |R2|<2(ξ2/N). Next, we have 3.25 where |R3|≤(n+1)(ξ/N). We combine (3.24) and (3.25), 3.26 with |R4|≤(ξ/N)(n+2ξ+1). If follows from (3.26) that . Finally, we use the bound |eR4−1|≤|R4| e|R4| (note that R4 is complex, not real) to conclude that 3.27 where |R5|≤(ξ/N)(n+2ξ+1) e(ξ/N)(n+2ξ+1). This completes the proof of proposition 3.1 and hence that of theorem 1.1. ■

## Acknowledgements

M.M. acknowledges the support of NSERC, the Natural Sciences and Engineering Research Council of Canada, under Discovery grant no. 205247. The work by G.P.B. was carried out under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy at Los Alamos National Laboratory under contract no. DE-AC52-06NA25396. Both authors thank the Institut Henri Poincaré for support during the final stages of this work (programme ‘Research in Paris’).