Markovian Master Equations: A Critical Study

We derive Markovian master equations of single and interacting harmonic systems in different scenarios, including strong internal coupling. By comparing the dynamics resulting from the corresponding Markovian master equations with exact numerical simulations of the evolution of the global system, we precisely delimit their validity regimes and assess the robustness of the assumptions usually made in the process of deriving the reduced dynamics. The proposed method is sufficiently general to suggest that the conclusions made here are widely applicable to a large class of settings involving interacting chains subject to a weak interaction with an environment.


Introduction
It is widely assumed that one of the crucial tasks currently facing quantum theorists is to understand and characterize the behaviour of realistic quantum systems. In any experiment, a quantum system is subject to noise and decoherence due to the unavoidable interaction with its surroundings. The theory of open quantum systems aims at developing a general framework to analyze the dynamical behaviour of systems that, as a result of their coupling with environmental degrees of freedom, will no longer evolve unitarily. If no assumptions are made concerning the strength of the system-environment interaction and the time-correlation properties of the environment, the dynamical problem may become intractable, despite that the functional forms of very general evolutions can be derived [1]. However, there exists a broad range of systems of practical interest, mostly in quantum optics and in the solid state physics, where it is possible to account for the observed dynamics by means of a differential equation for the open system's density matrix derived in the context of Markovian processes. Such a differential equation, the so-called Markovian (or Kossakowski-Lindblad) master equation, is required to fulfill several consistency properties such as being trace preserving and satisfying complete positivity [2,3,4,5,6,8,9,10].
However, from the theoretical point of view, the conditions under which these type of equations are derived are not always entirely clear, as they generally involve informal approximations motivated by a variety of microscopic models. This leaves open the range of validity of these equations, and which in some circumstances can lead to non physical evolutions. The situation becomes even worst as the complexity of the open system increases. In particular, it is not an easy question to decide whether the dynamics of a composite, possibly driven, quantum system can be described via a Markovian master equation, and if so, in what parameter regime. Actually, several groups have recently put forward operational criteria to check for deviations from Markovianity of real quantum evolutions [11,12,13,14].
The main propose of this work is to study such interacting open quantum systems, and show that there are Markovian master equations close to the real dynamics, characterizing the range of validity of each one. To this aim we have chosen a system consisting of quantum harmonic oscillators, as one can easily follow the exact dynamics using numerical simulations of a particular, but wide class of simple states, the so-called Gaussian states. Moreover, the proposed method is general enough to be applicable to non-harmonic systems and, in particular, when the coupling between oscillators is sufficiently weak so that their local dynamics is effectively two-dimensional, we expect the conditions obtained for strict Markovianity to be directly applicable to systems of interacting qubits.
The damped harmonic oscillator is the canonical example used in most references to discuss both Markovian and non-Markovian open system dynamics (see for instance [3,15,16,17,18,19,20] and references therein) and exact solutions in the presence of a general environment are known [21]. The dynamics of coupled damped oscillators, including those interacting with a semiclassical field, are significantly less studied, with most analysis focusing on evaluating the decoherence of initially entangled states provided that certain dynamical evolution, Markovian or not, is valid [22]. Recently, an exact master equation for two interacting harmonic oscillators subject to a global general environment was derived [23]. Here we will focus on the derivation of Markovian master equations for interacting systems. We will focus on a scenario where two harmonic systems are subject to independent reservoirs and present a detailed study based on the numerical simulation of the exact dynamics. The advantage of this approach is that it allows us to compute not only quantities for the damped system but also for the environment. This enables us to check the rigour of some of the assumptions usually made in obtaining a Markovian master equation and assess their domain of validity.
We have extensively studied three damped systems. For completeness, we start our analysis by considering a single harmonic oscillator (section 2) and subsequently move to the core of our study by analyzing the dynamics of two interacting harmonic oscillators (section 3), finding Markovian master equations for both weak and strong internal coupling. We finally address the dynamics of an harmonic oscillator driven by a semiclassical field (section 4), where different Markovian master equations have been obtained and studied depending on the values of the external Rabi frequency and the detuning from the oscillator's natural frequency. To make the reading more fluent, details of the simulations and the derivation procedure are left for the appendices.
In the following two introductory sub-sections, and with the aim of setting up the notation and making the presentation as self contained as possible, we present a brief discussion of how Markovian master equations are obtained in the weak coupling limit (section 1.1), and present a short review of the properties of the harmonic oscillator Gaussian states, which will be used in subsequent sections (section 1.2).

Markovian Master Equations
To derive Markovian master equations we follow the approach of projection operators initiated by Nakajima [24] and Zwanzig [25], see also [3,15,16] for instance. In this method we define in the Hilbert space of the combined system and environment H = H S ⊗ H E two orthogonal projection operators, Pρ = Tr E (ρ) ⊗ ρ E and Q = 1 − P.
Here ρ ∈ B(H) is the combined state and ρ E ∈ B(H E ) is a fixed state of the environment, which we choose to be the real initial (thermal, k B = 1) state, Note that Pρ gives all the necessary information about the reduced system state ρ S , so to know the dynamics of Pρ implies that one knows the time evolution of the reduced system.
We then assume that the dynamics of the whole system is given by the Hamiltonian H = H S +H E +αV , where H S and H E are the individual Hamiltonians of the system and environment respectively and V describes the interaction between them with coupling strength α. Working in the interaction picture ( = 1), and analogously forṼ (t), we obtain the evolution equation For the class of interactions that we are interested in as can be easily checked by applying it over an arbitrary state ρ ∈ B(H). It is not difficult to redefine the interaction Hamiltonian such that this always holds, see for example [10,20]. Our aim is to obtain a time-evolution equation for Pρ under some approximation, in such a way that it describes a quantum Markovian process. To this end, we apply the projection operators on equation (1), introducing the identity 1 = P + Q between The solution of the second equation can be written formally as This is nothing but the operational version of the variation of parameters formula for ordinary differential equations (see for example [26,27]), where the solution to the homogeneous equation d dt is given by the propagator where T is the time-ordering operator. Inserting the formal solution for We now assume that the initial state of the system and bath are uncorrelated, so that the total density operator is factorised into ρ(t 0 ) = ρ S (t 0 ) ⊗ ρ th . From this we find Qρ(t 0 ) = 0, which was guaranteed by our choice of P as projecting onto the initial state, and then by using (2) we finally arrive at with kernel Equation (5) is still exact. We now consider the weak coupling limit, by taking the kernel at lowest order in α, so that by again using condition (2) we get a Born approximation for (5): Note that we are not asserting here that the state of the bath is always ρ th , the term ρ S (s) ⊗ ρ th appears just as a result of the application of the projection operator (see discussion in section 2.3.3). Now we take the initial time t 0 = 0 and an elementary change of variable s by t − s in the integral yields We expect this equation to be valid in the limit α → 0, but in such a limit the change inρ S becomes smaller and smaller and so if we want to see dynamics we need to rescale the time by a factor α 2 [2,4,5] otherwise the right side of the above equation goes to zero. Thus in the limit α → 0 the integration is extended to infinity. However in order to get a finite value for the integral, the functions Tr E [Ṽ (t), [Ṽ (t − s), ρ B ]] must decrease appropriately. In particular this implies that they should not be periodic, which requires that the number of degrees of freedom in the environment must be infinite, as otherwise there will be a finite recurrence time. Moreover, asρ S changes very slowly in the limit α → 0, we can take it as a constant inside width τ B around s = 0 where is not zero, and so finally we obtain These informal arguments contain the basic ideas behind the rigorous results obtained by Davies [4,5].
Since we have started from a product state ρ(t 0 ) = ρ S (t 0 ) ⊗ ρ th , we require, for consistency, that our evolution equation generates completely positive dynamics. The last equation does not yet warrant complete positivity in the evolution [8], and so we need to perform one final approximation. To this end, note that the interaction Hamiltonian may be written as: where each A k can be decomposed as a sum of eigenoperators of the superoperator where This kind of decomposition can always be made [3,10]. On the other hand, by taking the Hermitian conjugate, Now, substituting the decomposition in terms of A k (ν) forṼ (t − s) and A † k (ν) forṼ (t) into equation (8) gives, after expanding the double commutator, where we have introduced the quantities with the last step being justified because ρ th commutes with exp(iH E t).
In equation (12) the terms with different frequencies will oscillate rapidly around zero as long as |ν ′ − ν| ≫ α 2 , so in the weak coupling limit these terms vanish to obtain Now we decompose the matrices Γ k,ℓ (ν) as a of sum Hermitian and anti-Hermitian parts where the coefficients form Hermitian matrices. In terms of these quantities (14) becomes where is a Hermitian operator which commutes with H S , as a consequence of (11). This is usually called the Shift Hamiltonian, since it produces a renormalization of the free energy levels of the system induced by the interaction with the environment. The dissipator is given by Returning to Schrödinger picture, the time-evolution equation is then just Note that the matrices γ k,ℓ (ν) are positive semidefinite for every ν, this is a consequence of the Bochner's theorem [28], that is, it is easy to check that the correlation functions Tr B † k (s)B ℓ ρ th are functions of positive type, and γ k,ℓ (ν) are just the Fourier transform of them. With this final remark we conclude that the equation (15) generates a completely positive semigroup [6] and so defines a proper Markovian master equation, i.e. a completely positive semigroup.

Gaussian States
We saw in the last section that to avoid a finite recurrence time, the number of environment degrees of freedom should strictly tend to infinity. However, in practice, the recurrence time grows very rapidly with the size of the environment and so one can still test the validity of such equations with only a finite, yet still large environment model, as long as the domain of interest is restricted to early times. The prototypical example of which is afforded by a collection on n harmonic oscillators. In fact, such models are often explicitly included in master equation derivations both due to their easy handling and due to realistic physical justification. Phenomenologically speaking, they correctly describe both quantum Brownian motion and the derivation of Langevin style equations from first principles [16]. However, they also provide a convenient numerical testing ground as the number of variables needed to model such systems scales polynomially in the number of degrees of freedom. This is because the harmonic oscillator falls into a class of quantum states known as Gaussian states, which are entirely characterised by their first and second moments. We now review some of their basic properties [7].
For any system of n canonical degrees of freedom, such as n harmonic oscillators, or n modes of a field, we can combine the 2n conjugate operators corresponding to position and momentum into a convenient row vector, The usual canonical commutation relations (CCR) then take the form where the skew-symmetric real 2n × 2n matrix σ is called the symplectic matrix. For the choice of R above, σ is given by, One may also choose a mode-wise ordering of the operators, R = (x 1 , p 1 , ..., x n , p n ) T , in which case the symplectic matrix takes on the form, Canonical transformations of the vectors S : R → R ′ are then the real 2n−dimensional matrices S which preserve the kinematic relations specified by the CCR. That is, the elements transform as R ′ a = S ab R b , under the restriction, This condition defines the real 2n-dimensional symplectic group Sp(2n, R). For any element S ∈ Sp(2n, R), the transformations −S, S T and S −1 are also symplectic matrices, and the inverse can be found from S −1 = σS T σ −1 . The phase space then adopts the structure of a symplectic vector space, where (19) expresses the associated symplectic form. Rather than considering unitary operators acting on density matrices in a Hilbert space, we can instead think of all the quantum dynamics taking place on the symplectic vector space. Quantum states are then represented by functions defined on phase space, the choice of which is not unique, and common examples include the Wigner function, Q-function and the P -function. Often one has a particular benefit for a given physical problem, however for our purposes we shall consider the (Wigner) characteristic function χ ρ (ξ), which we define through the Weyl operator as Each characteristic function uniquely determines a quantum state. These are related through a Fourier-Weyl transform, and so the state ρ can be obtained as We then define the set of Gaussian states as those with Gaussian characteristic functions. Equivalent definitions based on other phase space functions also exist, but for our choice we consider characteristic functions of the form, where C is a 2n×2n real matrix and D ∈ R 2n is a vector. Thus, a Gaussian characteristic function, and therefore any Gaussian state, can be completely specified by 2n 2 + 3n real parameters. The first moments give the expectation values of the canonical coordinates d j = Tr[R j ρ] and are related to D by d = σ −1 D, while the second moments make up the covariance matrix defined by These are related to C by the relation C = σ T Cσ. It is often the case that only the entanglement properties of a given state are of interest. As the vector d can be made zero by local translations in phase space, one can specify the state entirely using the simpler relation, However, in this work we shall predominantly use the relation (25). Using this convention, we mention two states of particular interest; the vacuum state, and the n-mode thermal state. Both take on a convenient diagonal form. In case of the vacuum this is simply the identity C = 1 2n , while for the thermal state the elements are given by where ω j is the frequency of the j th mode, and the equilibrium temperature is given by T .

Operations on Gaussian States
We now consider Gaussian transformations.
As the R j are hermitian and irreducible, given any real symplectic transform S, the Stone-Von Neumann theorem tells us there exists a unique unitary transformation U S acting on H such that U S W ξ U † S = W Sξ . Of particular interest are those operators, U G , which transform Gaussian states to Gaussian states. To this end, we consider the infinitesimal generators G, of Gaussian unitaries U G = e −iǫG = 1 − iǫG + O(ǫ 2 ). Then to preserve the (Weyl) canonical commutation relations, the generators G must have the form G = 2n j,k=1 g jk (R j R k − R k R j )/2 [7]. It follows that Hamiltonians quadratic in the canonical position and momentum operators (and correspondingly the creation and annihilation operators) will be Gaussian preserving, in particular, the Hamiltonian for n simple harmonic oscillators, H = n j=1 ω j a † j a j . It is for this reason that harmonic oscillators provide such a useful testing ground for many body systems.
An additional, though simple, property worth highlighting is the action of the partial trace. Using the expression for the density matrix (23), it is straightforward to see the effect of the partial trace operation on the characteristic function. If we take a mode-wise ordering of the vector R = (R 1 , R 2 ), where R 1 and R 2 split two subspaces of n 1 and n 2 conjugate variables corresponding to partitions of the state space of ρ into H = H 1 ⊗ H 2 , then the partial trace over H 2 is given by That is, we need only consider the characteristic function χ(ξ 1 ) associated to the vector R 1 . At the level of covariance matrices, we simply discard elements corresponding to variances including any operators in R 2 , and so the partial trace of a Gaussian state will itself remain Gaussian. Finally, we make some remarks regarding closeness of two Gaussian states. Given ρ 1 and ρ 2 the fidelity between them is defined as F (ρ 1 , ρ 2 ) = Tr √ ρ 1 ρ 2 √ ρ 1 2 , and is a measure of how close both quantum system are each other. Actually a distance measure can be defined as . This distance will be very useful for quantifying how well the dynamics generated by a Markovian master equation approximate the real one.
In general the fidelity is quite difficult to compute, however in the case of Gaussian states Scutaru has given closed formulas in terms of the covariance matrix [30]. For example, in case of one mode Gaussian states ρ G1 and ρ G2 , with covariance matrices C (1) and C (2) and displacement vectors d (1) and d (2) respectively, their fidelity is given by the formula

Damped Harmonic Oscillator
We will first consider a single harmonic oscillator damped by an environment consisting of M oscillators (see figure 1). We want to know under which conditions the Markovian master equation that we derived in the previous section for the evolution of the damped oscillator is valid. To this aim we will approach the exact dynamical equations of the whole system when M is large; these will be solved via computer simulation, and we can then compare this solution with the one obtained using a master equation. The Hamiltonian for the whole system will be given by ( = 1) Note that the coupling to the bath has been considered in the rotating wave approximation (RWA), which is a good description of the real dynamics for small damping Ω ≫ max{g j , j = 1, . . . , M} (e.g. in the weak coupling limit) [31]. For definiteness, in this paper we have chosen to distribute the environmental oscillators according to an Ohmic spectral density with exponential cut-off. In the continuous limit, this has the form [19] where α is a constant which modifies the strength of the interaction and ω c is the so-called cutoff frequency. Clearly J(ω) increases linearly for small values of ω, decays exponentially for large ones, and has its maximum at ω = ω c . Of course any other choice of spectral density could have been taken, but this in turn would require a re-analysis of the master equations' range of validity.

Exact Solution
The exact solution of this system can be given in terms of the time-evolution of the collection {a, a j } in the Heisenberg picture [17]. From (30) we have and so by writing A = (a, a 1 , a 2 , . . . , a M ) T , the system of differential equations may be expressed as where W is the matrix and the solution of the system will be given by Analogously, the evolution of the creation operator will be We can also compute the evolution of position and momentum operators X = and similarly in these expressions, T R and T I are the self-adjoint matrices defined by So, the time-evolution of the vector R = (x, x 1 , . . . , x M , p, p 1 , . . . , p M ) T will be given by note that the size of M is 2(M + 1) × 2(M + 1). Due to the linearity in the couplings in H, an initial (global) Gaussian state ρ G will remain Gaussian at all times t, and so we can restrict our attention to the evolution of its covariance matrix Particularly, since we are interested in just the first oscillator, we only need the evolution of the 2 × 2 submatrix {C ij ; i, j = 1, M + 2}. The evolution of pairs of position and momentum operators is and similarly for products of expectation values R i (t) R j (t) . So the elements of the covariance matrix at time t will be and for the first oscillator we have here (·, ·) denotes the scalar product, and the vectors M 1 and M M +2 are given by More details of how this exact solution is simulated in order to approach the Markovian master equation description are given in Appendix A.

Markovian Master Equation
The damped harmonic oscillator is a standard example for the derivation of master equations (see for example [3,17,18,20]). The Markovian master equation (15) is given by d dt whereΩ is a renormalized oscillator energy arising for the coupling to the environment (here P.V. denotes the Cauchy principal value of the integral),n is the mean number of bath quanta with frequency Ω, given by the Bose-Einstein distribution and γ is the decay rate, which is related to the spectral density of the bath J(ω) = Note that the shift ∆ is independent of the temperature, and although its effect is typical small (e.g. [3,18]) we will not neglect it in our study. For an ohmic spectral density the frequency shift is where Ei is the exponential integral function defined as In addition, note that the equation (49) is Gaussian preseving [32], as it is the limit of a linear interaction with an environment and so the total system remains Gaussian while the partial trace also preserves Gaussianity.

Study of the Approximations
As a first step, we have plotted the variance of the x coordinate for two different initial states of the system, these are a thermal and a squeezed state, see figure 2. The last plot clearly illustrates the closeness of the results for the Markovian master equation, when compared to the effect of the Lamb shift. To explore this further, we now study several effects which pertain to the validity of this equation, by calculating the distance (in terms of the fidelity) between the simulated state ρ

Discreteness of the bath
Due to the finite number of oscillators in the bath, we can only simulate inside a bounded time scale free of the back-action of the bath. This produces revivals in the visualized dynamical quantities for times t < τ R , where τ R is the recurrence time of the bath. Of course, the time after which these revivals arise increases with the number of oscillators in the bath, and roughly speaking it scales as τ R ∝ M. This behaviour is shown in figure 3, where the distance between the simulation and the Markovian master equation for a system initially in a thermal state with temperature T S = 30 is plotted as a function of the time and the number of oscillators.

Temperature
It is sometimes claimed that for ohmic spectral densities the Markovian master equation (49) is not valid at low temperatures [18,19]. Of course, one must make clear the context in which this claim is made, and so for definiteness, let us focus on the validity with respect to the bath temperature. A detailed discussion of this situation can be found in the book by Carmichael [18]. There the argument is based on the width of the correlation function C 12 g j a j , which increases for an Ohmic spectral density as the bath temperature decreases. More specifically, in the derivation of the Markovian master equation two kinds of correlation functions appear,  We may call C 12 (s) ≡ C(−s, T ) and C 21 (s) ≡ C(s, T ) + C 0 (s), and so in the continuous limit is the so-called Hurwitz Zeta function, which is a generalization of the Riemann zeta function ζ(z) = ζ(z, 1) [33].  In the left plot of figure 4, the absolute value of C(s, T ) is plotted for different temperatures. Note that the spreading of the correlation function is mainly caused by its "height" decrease, that is, in the limit T → 0, C(s, T ) → 0. So one may also expect that the contribution of these correlations to the motion becomes less important as T → 0, in such a way that the problem of the infinite width can be counteracted, and this is indeed what seems to happen. To visualize this more carefully we have plotted in the right of figure 4 the full weight at half height (FWHH) for both C 0 (s) and C(s, T ). In order to make valid the Markovian approximation, the typical time scale for the evolution of the system due to its interaction with the bath τ S must be large in comparison with the decay time τ B of the correlation functions. Loosely speaking, this can be characterized by the FWHH.
From figure 4 one sees that for small temperatures τ B (i.e. FWHH) is quite large, so it is expected that the Markovian approximation breaks down for values of T such that τ S τ B . However if α is small enough this will happen for values where the contribution of C(s, T ) to the convolution integrals is negligible in comparison with the contribution of C 0 (s), whose FWHH will remain constant and small with respect to τ S . As a rough estimation, using the parameters in figure 2, we find that to get a value of the FWHH comparable with τ S ∼ 1/ √ α ∼ 22.4, we need a temperature of at least T ∼ 0.05. Both contributions enter in the Markovian master equation derivation via some convolution with the quantum state and one oscillating factor. We may get a very informal idea of how both contributions matter by looking at their maximum values at s = 0, for example C(s = 0, T = 0.05) = 3.27391 × 10 −7 and C 0 (s = 0) = 0.018, and so it is clear that C(s, T = 0.05) will not have a large effect on the dynamics. For large temperatures the FWHH of C(s, T ) remains small though now larger than C 0 (s), so it is expected that in the limit of high temperatures the accuracy of the Markovian master equation stabilizes to a value only a little worse than for T = 0.
All of these conclusions are illustrated in figure 5, where the fidelity between the state from the simulation and that from the Markovian master equation is plotted. The behaviour at very early times is mainly related to the choice of the initial state of the system, and reflects how it adjusts to the state of the bath under the Markovian evolution [34], different tendencies have been founded depending on the choice of initial state. However the behaviour with temperature is visible at longer times (since τ B ∼ FHWW increases with T ) which is in agreement with the conclusions drawn from the correlation functions (see small subplot). At zero temperature (blue line) the results are in closest agreement, however, as the temperature is increased to T = 0.1 the correlation function broadens, which leads to a degradation (albeit small) in the modelling precision. As the temperature increases further, the influence of this correlation function becomes more important and the FWHH decreases to a limiting value (see the plot on the right of figure 4), this convergence is reflected by the red, cyan and purple lines which show that the accuracy at large temperatures stabilizes to only a little worse than that at T = 0, as was expected from figure 4.
In summary, the Markovian master equation (49) does not properly describe the stimulated emission/absorption processes (the ones which depend on C(s, T )) for low temperatures, however the temperatures when this discrepancy is apparent are so small that the contribution from stimulated process are negligible in comparison with spontaneous emission, and so the discrepancy with the Markovian master equation is never large.

Assumption of factorized dynamics
In the derivation of the Markovian master equation, one can arrive at equation (7) by iterating the Von-Neumann equation (1) twice and assuming that the whole state factorizes as ρ(t) ≈ ρ S (t) ⊗ ρ th at any time ( [3,10,18,20]). This assumption has to be understood as an effective model for arriving at equation (7) without the use of projection operator techniques, however it does not make sense to assume that the physical state of the system is really a factorization for all time. Taking advantage of the ability to simulate the entire system we have plotted the distance between the simulated whole state ρ(t) and the ansatz ρ S (t) ⊗ ρ th as a function of time, see figure 6. On the left we have plotted the distance for M = 350 oscillators in the bath, actually we have checked from several simulations that the results turn out to be independent of the number of oscillators as long as the maximum time is less than the recurrence time of the system. From figure  3 we see that t = 50 is less than the recurrence time for M = 175, and so we have used this value and plotted the distance for different coupling strengths on the right. It is clear that this distance is monotonically increasing in time (strictly, in the limit of an environment with infinite degrees of freedom), and the slope decreases with coupling strength. In section 1.1 we pointed out that the weak coupling approach make sense if the coupling is small and the environment has infinite degrees of freedom. This fits with the usual argument to take ρ ≈ ρ S (t) ⊗ ρ th in more informal derivation of Markovian master equations, that is "the state of the environment is not so affected by the system", but we stress again that this is an effective approach, without any physical meaning on the real state ρ.

Two Coupled Damped Harmonic Oscillators
We now consider two coupled harmonic oscillators, which for simplicity we take to have the same frequency Ω 1 = Ω 2 = Ω, and each locally damped by their own reservoir (see figure 7), the Hamiltonian of the whole system is where the free Hamiltonians are given by with the couplings to the baths, and the coupling between oscillators, V 12 = β(a † 1 a 2 + a 1 a † 2 ). Again we have employed the rotating wave approximation, and so we assume Ω ≫ β. For the case of Ω ∼ β we must keep the antirotating terms a 1 a 2 and a † 1 a † 2 . However note that the eigenfrequencies of the normal modes become imaginary if ω < 2β (see for example [35]) and the system then becomes unstable, so even when keeping the antirotating terms, we must limit β if we wish to keep the oscillatory behaviour.

Exact Solution
For the exact solution, the extension to two oscillators follows closely that of a single damped harmonic oscillator. Again, we work in the Heisenberg picture, and wish to solve for the vector A = (a 1 , a 11 , . . . , a 1M , a 2 , a 21 , . . . , a 2M ) T , given the differential equation, where W is now given by the matrix The simulation process is then analogous to that of section 2.1.

Markovian Master Equations
Unfortunately, the derivation of a Markovian master equation for coupled systems introduces a number of additional complications. If the oscillators are uncoupled β = 0, it is obvious that the Markovian master equation for their joint density matrix will be a sum of expressions like (49), where here each frequency shift, decay rate and number of quanta are individually computed via equations (51), (52) and (50) for each bath j. However for finite intercoupling we split the analysis in two subsections.

Small intercoupling β
If β is sufficiently small to not affect the shift and decay rates, one can expect a Markovian master equation of the form an example of which for coupled subsystems can be found in [36], and we have given the details of a derivation based on projection operators in Appendix B.1. In addition, this kind of approximation is often made in other contexts such as with damped systems driven by a classical field [18]. Such a case will be analyzed in detail in section 4.

Large intercoupling β
To go further we must work in the interaction picture generated by the Hamiltonian H 0 = H free + V 12 and apply the procedure described in section 1.1. The details of the derivation are left for Appendix B.2, what is important however, is that the non-secular terms oscillate with a phase e ±2iβt so in order to neglect them we must impose β ≫ α, therefore the resultant equation is, in some sense, complementary to (58) valid if α β. The final Markovian master equation in this regime takes the form where γ j , ∆ j andn j are evaluated according to the spectral density and temperature of the bath j and Ω ± = Ω ± β.

Study of the Approximations
By virtue of the derivation, equations (58) and (59) preserve both complete positivity and Gaussianity (because they arise from a linear interaction with the environment).
Thus we can test their regimes of validity using simulations of Gaussian states, and the appropriate fidelity formulas. In figure 8 we have plotted the fidelity between both states for the Markovian master equation (58) (left side) and for (59) (right side).  From these results one concludes that when modeling a system with multiple baths at different temperatures equations (58) and (59) are each accurate in their theoretically applicable regimes. However, for baths at the same temperature, it seems both equations give good results. A natural, and important, question is to ask is whether an intermediate range of couplings exist, such that neither (58) or (59) give useful results. In figure 9 the fidelity between the simulation and the Markovian master equation  We see that for the parameters shown on the plot, there is a small range between β ∼ 0.01 − 0.02 where neither Markovian master equation obtains a high precision. However, note that this range becomes smaller as the coupling with the bath decreases, and so generally both master equations cover a good range of values of β.

Baths with the same temperature
We now examine the role of the bath temperatures in more detail. Since the simulations seem to produce good results for both Markovian master equations when the temperature of the local baths are the same, regardless of the strength of the intercoupling, it is worth looking at why this happens. In the case of equation (59) it is reasonable to expect that this will remain valid for small β, because when β → 0 this equation approaches (58) if the bath temperatures and spectral densities are the same. That is, the off-diagonal terms of the matrices K (E) and K (A) do not contribute much,β ∼ β and the rest of coefficients become approximately equal to those in (58.) Note this only happens under these conditions.
Essentially the same argument applies to equation (58) in the large β limit. On the one hand, for a relatively small value of β (= 0.1) in comparison to ω, the off-diagonal elements of the matrices K (E) and K (A) in the master equation (59) are unimportant in comparison with the diagonals. On the other hand, the diagonal terms are also alike for the same reason, and so both master equations will be quite similar. However note that at later times the behaviour of both equations start to differ, and the steady states are not the same. By construction, the steady state of equation (59) is the thermal state of the composed system [3,4], whereas that of master equation (58) is not (although it tends to the thermal state as β → 0 of course). Surprisingly the divergences between both equations, even for large times, are actually very small, see figure 10. In some cases, while the steady state of (58) is not strictly thermal, the fidelity with that of (59) is more than 99.999%.

Driven Damped Harmonic Oscillator
One situation which is also interesting to analyze is that of adding a driving term in the Hamiltonian of the damped oscillator. At this stage we consider again one single oscillator, damped by a thermal bath and driven by a coherent field (figure 11). This is described by a semiclassical Hamiltonian in the rotating wave approximation: here ω L is the frequency of the incident field and r the Rabi frequency.

Exact Solution
To obtain the exact solution of this system let us consider for a moment the Schrödinger picture,

We solve this equation by means of the unitary transformation
where H 0 = (Ω − ω L )a † a + r(a + a † ) + M j=1 (ω j − ω L )a † j a j + M j=1 g j (a † a j + aa † j ) is time-independent. Returning to the Schrödinger picture, the evolution of the states is then, In order to avoid differential equations with time-dependent coefficients, we can study the evolution in a X-P time rotating frame; in that frame the annihilation (and creation) operatorsã = e −iHrott ae iHrott will evolve according tõ which is quite similar to (32) but with the additional time-independent term r. Following the notation of section 2.1 we can write here b = (r, 0, . . . , 0) T and W 0 is found from (35) as W − ω L 1. The solution of this system of differential equations is If W 0 is invertible this equation can be written as Analogously to (38) and (39) we find where T 0 R and T 0 I are as in (40) for W 0 . Thus, by writing we find that the position and momentum expectation values evolve as Note that in this case the first moments of the state change, despite R(0) = 0. To calculate the evolution of the covariance matrix, we proceed in the same way as before, and analogously for the solutions for R j (t)R i (t) and R i (t) R j (t) . Combining these terms, we find the B cancel and so, in a similar fashion to (44), (45) and (46), where, of course, M 0 1 and M 0 2 are as in (48) for M 0 .

Markovian Master Equations
In order to derive a Markovian master equation for this system we must take account of two important details. First, since the Hamiltonian is time-dependent the generator of the master equation must also be time-dependent, whose solution defines a family of propagators E(t 2 , t 1 ) such that ρ S (t 2 ) = E(t 2 , t 1 )ρ S (t 1 ), These can be written formally as a time-ordered series where T is the well-known time-ordering operator. Similarly to the case of timeindependent equations it can be shown that the family E(t 2 , t 1 ) is completely positive for all (t 2 ≥ t 1 ) if and only if L t has the Kossakowski-Lindblad form for any time t [13]. The second problem is that there is an absence of rigorous methods to arrive at a completely positive master equation in the Markovian limit when the Hamiltonian is time-dependent, with the exception of adiabatic regimes of external perturbations [37]. Fortunately in this case, due to the simple periodic time-dependence of the Hamiltonian, we will be able to obtain Markovian master equations valid for large (to some degree) Rabi frequencies, even though the complexity of the problem has increased. In our derivation, we will distinguish between three cases: these will be when the Rabi frequency is very small; when the driving is far off resonance (|ω L − Ω| ≫ 0) and finally the identical case without the secular approximation.
The details of the derivation are left for the Appendix B.3, but in these three cases we find a Markovian master equation with the structure d dt where D is given by (57),Ω = Ω + ∆ is the same as for a single damped oscillator, andr is a renormalized Rabi frequency due to the effect of the bath. Note that as the incident field alters the position operator of the oscillator, which in turn couples to the bath, one should expect that the field is itself also effected by the environment. For small Rabi frequencies an argument similar to section 3.2.1 gives simplȳ whereas, when the driving field is far from resonance, |ω L − Ω| ≫ 0, we obtain Finally, if we neglect the secular approximation, this regime yields Without entering into the details of the derivation, one sees that equations (74) and (75) are problematic on resonance |Ω − ω L | ∼ 0. This is due to two approximations, one is the secular approximation in (74), and the other is the second order in the perturbative series. In the derivation in Appendix B.3 it is clear why in this case the series diverges for |Ω − ω L | ∼ 0.

Study of the Approximations
Note that in this case the range of validity of each equation is now more ambiguous than in previous sections where we have dealt with undriven systems. Which one is more appropriate is going to be discovered by simulation, although one could suppose that the more elaborate equations (74) and (75) would provide the better approximation.
However, there is still the question of how effective they are, and whether the additional effort required to obtain them is worthwhile in comparison to the simpler equation (73).
In addition note that in every case the covariance matrix is unaffected by the driving term, which only produce a change in the first moments. Furthermore, as the fidelity is invariant under unitary operations, we are always free to work in the frame rotating with the field. Therefore, all calculations can be performed with the rotating observables. In figure 12 the fidelities are plotted for close to and far from resonance. Compare the amount of disagreement with the fidelity of a single damped oscillator in figure 5. For global features, the more elaborate equation (75) works better in both cases, although the difference with (73) is very small. As expected, the choice of (74) is preferable to the choice of (73) when out of resonance, but gives quite poor results when close to resonance. However, when off resonance the difference among the three choices is essentially small.
Given these results, it is worthwhile to look at how the fidelities at one fixed time vary as a function of the detunning, this is done in figure 13 (note we choose a large value for the time, so we avoid the potentially confusing effect due to the oscillatory behaviour depicted in figure 12).
Here we see that both (74) and (75) fail close to resonance, as was expected from the perturbative approach. Equation (73) gives good results due to the small Rabi frequency, however note in comparison to (75) the accuracy quickly drops off as we move away from ω L − Ω = 0. A similar effect can be seen when compared to (74) for larger detunnings.
Finally, in figure 14 we test the dependency of the fidelities on the strength of the In summary, for the case of a driven damped harmonic oscillator the difference in accuracy among Markovian master equations is generally small. Equations (74) and (75) work better except in the case of resonance, where (73) gives more accurate results, as long as the Rabi frequency is small. The justification to use one equation over another will depend on the context and the accuracy which one wants to obtain, but given that the differences are so small the simplest choice (73) seems to be the more "economical" way to describe the dynamics.

Conclusions
We have obtained and studied the range of validity of different Markovian master equations for harmonic oscillators by means of exactly simulating the dynamics, and comparing the predictions with those obtained from evolving the system using the master equations. In particular, • We have shown that the system-environment state factorization assumption for all times has to be understood in general as an effective model by deriving the same equation using the projection operator technique.
• We analysed two strategies for finding completely positive Markovian master equations for two harmonic oscillators coupled together under the effect of local baths, indicating that both are complementary in their range of validity. Moreover, when the temperature of the local baths is the same the difference between them is quite small.
• In the same spirit, we derived time inhomogeneous completely positive Markovian master equations for a damped oscillator which is driven by an external semiclassical field. We studied the validity of each one and pointed out that completely positive dynamics can be obtained even without secular approximation (for these kinds of inhomogeneous equations).
Despite the fact that we have focused on harmonic oscillator systems, the proposed method is general and we expect that non-harmonic systems should behave in a similar manner with respect to the validity of the equations. This suggest that the general conclusions made here are widely applicable to any other settings involving a weak interaction with an environment. In this regard, we hope that the present study may help in providing a better understanding and a transparent description of noise in interacting systems, including those situations where the strength of the internal system interaction is large. There are currently many quantum scenarios open to the use of these techniques, including realizations of harmonic and spin chains in systems of trapped ions [38], superconducting qubits [39] and nitrogen-vacancy (NV) defects in diamond [40].
Moreover, interacting systems subject to local reservoirs have been recently treated under the assumption of weak internal system interaction in theoretical studies ranging from the excitation transport properties of biomolecules [41] to the stability of topological codes for quantum information [42].

Acknowledgments
A.R. acknowledges Alex Chin for fruitful discussions. This work was supported by the STREP projects CORNER and HIP, the Integrated projects on QAP and Q-ESSENCE, the EPSRC QIP-IRC GR/S82176/0 and an Alexander von Humboldt Professorship.

Appendix A. Details of the simulation
In order to make an appropriate comparison between the exact evolutions, such as those in sections 2.1, 3.1 and 4.1, and the corresponding master equations, we must make a careful choice of a number numerical parameters. In practice, however, this is not a difficult issue. The essential ingredient is to choose the couplings to the bath according to the desired spectral density. Throughout this paper, we have made the choice (31), The first step in picking g j is to remove the Dirac delta functions by integrating over a frequency range bounded by a frequency cut-off ω max , which means g 2 j ′ ≈ αω j ′ e −ω j ′ /2ωc ∆ω j ′ due to the decomposition of the integral in terms of Riemann sums. We should also take care to set the range of oscillators, ω max , large enough to cover (31) significantly. For example, if we take ω 1 = c, with c small, then one possible convention is to take ω max such that J(ω max ) = J(c), and so we neglect all possible oscillators with coupling constant less than J(c)∆ω 1 . Another polisher convention is to take ω 1 and ω max such that However, in practice this choice is not really a crucial point.
whereṼ SB (t) =Ṽ 1B1 (t) +Ṽ 2B2 (t) and for simplicity we have assumed that the strength of the coupling to each bath is identical (the reader will note afterwards that this is not a crucial assumption). We now define the projector Pρ(t) = Tr (B1,B2) [ρ(t)] ⊗ ρ th1 ⊗ ρ th2 , along with Q = 1 − P. The application of the projection operators on (B.1) yields and so (c.f. section 1.1) we find a formal solution to the second equation as where here the kernels are K 1 (t, s) = αβPV SB (t)G(t, s)QV 12 (s)P = 0, The first vanishing because V 12 (s) commutes with P and QP = 0. If we consider the second kernel, weak coupling implies α β, and so to second order in α and β this becomes which has exactly the same form as (6) and therefore the equation of motion becomes Finally we note that note that RWA approximation implies Ω ≫ β so both normal mode frequencies are positive. We can reexpress the interactions with the baths in terms of these new operators, the benefit of this is that it allows us to easily deal with the interaction picture with respect to H 0 = H 12 + H B1 + H B2 . By following the method of section 1.1 we obtain the analog of (8), where we have noted PV 1B1 V 2B2 P = PV 2B2 V 1B1 P = 0. Each of the above terms correspond, essentially, to one of a pair of two free harmonic oscillators with frequencies Ω + and Ω − , coupled to a common bath. Consequently, we can deal with them separately. Starting with the first term we decompose the interaction in to eigenoperators of [H 12 , ·] (see (11)) Notice the A 1 operator can be written as and so we can write (B.8) as which in interaction picture becomes Now, for the first element of (13) we have where the mean number of quanta in the first bathn 1 (ω 1j ) with frequency ω 1j , is given by the Bose-Einstein distribution (51). Going to the continuous limit we take M → ∞ and introduce the spectral density of the first bath J 1 (ω) = j g 2 1j δ(ω − ω 1j ), Next we perform the secular approximation; the cross terms ν ′ = ν in the above expression, which go as e ±2βti , can be neglected provided that 2β is large in comparison with the inverse of the relaxation rate (β ≫ α) and so we obtain Taking t 0 = 0 without lost of generality, the time-evolution equation forρ(t) = U † (t, 0)ρ(t)U(t, 0) iṡρ (t) = −i[Ṽ (t),ρ(t)], (B.18) so by following the analogous procedure for time-independent generators, one immediately deals with the problem that is not clear whether there exists a similar eigenoperator decomposition forṼ (t) = U † (t, 0)V U(t, 0) (V = M j=1 g j (a † a j + aa † j )) as in (10) and (11). Note however that the operatorÃ 1 (t) =ã(t) satisfies a differential equation with periodic terms iȧ(t) = [ã(t), H 0 (t)] = Ωã(t) + re −iω L t . (B.19) This kind of equation can be studied with the well-established Floquet theory (see for example [26,27]), particularly it is possible to predict if its solution is a periodic function.
In such a case, the operator in the new picture would have a formal decomposition similar to that in (10) and (11), such thatÃ k (t) = ν A k (ν)e iνt . This would then allow us to follow a similar procedure to that for time-independent Hamiltonians. Note that the importance of such a decomposition is that the operators A k (ν) are themselves time-independent. Such ideas have already been used before in, for instance, [45,46].
The solution to equation (B.19), with the initial conditionã(0) = a and for Ω = ω L is given byã Before continuing note that in the perturbative series of (B.18), the "strength" of the interactionṼ (t) is now not solely dependent on the coupling with the bath. This is because the operators A(ν) depend linearly on r w L −Ω , so when this ratio becomes large we expect that the approximation breaks down, i.e. for r ≫ 1 or very close to resonance |w L − Ω| ≈ 0.
Next we assume that the detunning is large enough |ω L − Ω| ≫ α, |ω L − Ω| 2 ≫ αr in order to make the secular approximation and after some tedious, but straightforward, algebra we find the master equation in the interaction picture to be d dtρ where D(·) has again the form of (57). Finally, on returning to the Schrödinger picture we have, and so one can see thatã(t) is not a periodic function, so the desired decomposition as a sum of exponentials with time-independent coefficients does not exist. On the other hand, the decomposition (B.20) tends to (B.24) in the limit ω L → Ω, so we may attempt to work with this decomposition and wonder whether on resonance the new master equation holds in this limit as well (in fact, we have shown that this is not true in section 4.3). The only problem to deal with is the possible lack of positivity due to the absence of the secular approximation. However, note that in this particular case only a commutator term arises from the cross terms in the analog of equation (B.15), so positivity is not lost. In fact, we obtain an equation similar to (B.22) except for an additional correction to the Rabi frequency: Note that to first order in r and α we again obtain the equation (56).