Going beyond Local and Global approaches for localized thermal dissipation
Donato Farina Giulio De Filippis Vittorio Cataudella Marco Polini Vittorio Giovannetti
quant-ph cond-mat.mes-hall





^1NEST, Scuola Normale Superiore, I-56126 Pisa, Italy.

^2Istituto Italiano di Tecnologia, Graphene Labs, Via Morego 30, I-16163 Genova, Italy.

^3SPIN-CNR and Dip. di Fisica - Università di Napoli Federico II, I-80126 Napoli, Italy.

^4INFN, Sezione di Napoli - Complesso Universitario di Monte S. Angelo, I-80126 Napoli, Italy.

^5NEST, Scuola Normale Superiore and Istituto Nanoscienze-CNR, I-56127 Pisa, Italy.





Identifying which master equation is preferable for the description of a multipartite open quantum system is not trivial and has led in the recent years to the local vs. global debate in the context of Markovian dissipation. 
We treat here a paradigmatic scenario in which the system is composed of two interacting harmonic oscillators A and B, with only A interacting with a thermal bath - collection of other harmonic oscillators - and we study the equilibration process of the system initially in the ground state with the bath finite  temperature.

We show that the completely positive version of the Redfield equation obtained using coarse-grain and an appropriate time-dependent convex mixture of the local and global solutions give rise to the most accurate semigroup approximations of the whole exact system dynamics, i.e. both at short and at long time scales, outperforming the  local and global approaches.



Going beyond Local and Global approaches for localized thermal dissipation
    D. Farina^1,2, G. De Filippis^3,4, V. Cataudella^3,4, M. Polini^2 and V. Giovannetti^5
    December 30, 2023
==========================================================================================










§ INTRODUCTION


In the rising field of quantum technology <cit.>, considering a quantum system isolated from its surroundings is a non-realistic idealization.
In the majority of the implementations of quantum information algorithms <cit.> and quantum computation <cit.>,

the interaction with the environment is detrimental for quantum resources, becoming a crucial ingredient to monitor, with the scope of reducing its effects or with the aim of accounting for it by applying quantum error correction methods.
Interestingly, in more rare cases the environment itself acts as a mediator for the production of quantum correlations into the system <cit.>.


Unfortunately
 our ability in accounting for environmental effects is severely limited by  the difficulty of 
keeping track  of  the exact dynamics of the entire system-environment compound: a  problem which is made computationally hard  by the large number of degrees of freedom involved in
the process.
For this reason, effective models for the way the environment acts on the reduced system density matrix have been developed, leading to the master equation (ME) formalism <cit.>. 
The lowest level of approximation contemplates the assumption of weak system-environment coupling (Born approximation)
and time-divisibility for the system dynamics (Markov approximation).
This  leads to the Redfield equation <cit.> which regrettably, while being able to capture some important features of the model <cit.>, does not ensure positive (and hence completely positive) evolution <cit.>.
In quantum mechanics, the positivity
 of density matrices – i.e. the fact that all their eigenvalues are non-negative –
 is an essential  property  imposed by 
  the probabilistic interpretation of the theory <cit.>. 
Allowing for mathematical  structures that do not comply with such 
requirement paves the way to a series of inconsistencies that include
negative probabilities of measurements outcomes, violation of the uncertainty relation, 
and ultimately the non-contractive character of the underlying dynamics.
Ways to correct or to circumvent the pathology exhibited by the Redfield equation 
typically relay on 
 the full <cit.> or the partial <cit.> implementation of the secular approximation: a coarse-grain temporal  average of the system dynamics which, performed in conjunction with the above mentioned
Born and Markov approximations,  leads to a more reliable differential equation for the system
density matrix known as the 
 Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) master equation <cit.>.

 
The situation becomes more complicated when the system is composed of two or more interacting subsystems that are locally coupled to possibly independent reservoirs <cit.>. In this case, a brute force application of a full secular approximation
leads to the so called global ME, a GKSL equation obtained
under the implicit assumption  that  the environment will perceive the composite system 
as a unique body irrespectively from the local structure of their mutual  interactions.
While formally correct in terms of the positivity and complete positivity requirements and
predicting long term  behaviours  which are thermodynamically consistent, 
the resulting ME  is prone to introduce errors in the short term description of the dynamical 
process.
A suitable alternative  is provided by 
the so called local ME approach where,
contrarily to the Global ME, 
 each subsystem is assumed to independently interacts with its own environment, keeping track of the local nature of the microscopic interaction.
 Despite in certain situations it can imply the breaking of the second law of thermodynamics <cit.>, it allows for a more precise description of the short term dynamics of the composite system. A local approach is usually allowed when the subsystems interact weakly between each other <cit.>. As well as the Global ME, the local ME can be microscopically derived <cit.> and is in GKSL form.
Notably, such master equation has recently acquired full dignity showing that it 
exactly  
describes the dynamics induced by an engineered bath schematized by a collisional model
<cit.>. 
Furthermore, even under a more conventional description of the environment, 
thermodynamics inconsistencies only 
occur  at the order of approximation where the local approach is not guaranteed to be valid and, eventually, it is possible to completely cure such inconsistencies by implementing a perturbative treatment   
around the local approximation <cit.>.



The scope of the present work is to test the effectiveness of different classes of MEs  to describe  the system dynamics, particularly focusing on alternative approaches beyond those adopted in deriving the local and global  MEs  and using as benchmark a model that we are able to solve exactly.
Differently from previous studies <cit.>, where the focus was on the steady state properties of a bipartite system with each subsystem coupled to a different thermal reservoir, 
we deal with a bipartite system  asymmetrically coupled to a single thermal bath and analyze its whole dynamics including both the transient and asymptotic regime.
More specifically, in our case the system of interest  is composed of two interacting harmonic oscillators A and B, with only A  microscopically coupled with an external bosonic thermal bath
described as a collection of extra harmonic oscillators. 
About the exact dynamics benchmark, the unitary evolution of the joint system+environment compound has been calculated by restricting ourself to exchange interactions and gaussian states <cit.>. 
Our analysis leads to the conclusion that
the completely positive version of the Redfield equation  obtained as described in <cit.> by applying the secular approximation via coarse-grain averaging in a partial and tight way, provides a  semigroup description of the system dynamics that 
outperform both the local and Global ME approaches. 
We also observe that analogous advantages can be obtained by adopting a phenomenological 
description of the system dynamics,  constructed in terms of 
an appropriate time-dependent convex mixture of the local and Global ME solutions. 


Despite the selected model has been chosen primarily for its minimal character, possible implementations of the set up we deal with can be found in cavity (or in circuit) quantum electrodynamics. 
An example is the open Dicke model <cit.> 
for large enough number of two-level atoms inside the cavity <cit.> and 
assuming that the interaction of the cavity mode with the radiation field is more relevant than the direct coupling of the radiation field with the atoms.

Alternatively, our bipartite system may directly describe coupled cavities in an array  <cit.> in the instance of two cavities.

About the kind of dynamics we chose, 
it may be of interest 
for ground state storage in quantum computation
<cit.> or, conversely, for thermal charging tasks <cit.>.



The paper is organized as follows. 

In Sec. <ref> we introduce the  model.
The different approximations are described in Sec <ref>.
In Sec. <ref>  we integrate the dynamical evolution under the various approximations and 
present a comparison between the various results. 
In Sec <ref> we draw the conclusions and we discuss possible future developments.
Details on the approximation methods and on the evaluation of the exact dynamics are reported in the Appendix.  














§ THE MODEL


The model we consider is schematically described in Fig. <ref>. 
It consists into a bipartite system 𝒮 composed of two resonant bosonic modes A and B of frequency ω_0 and described by the ladder operators a,a^† and b, b^†,
 that interact through an excitation preserving coupling characterized by 
an intensity parameter g≥ 0. 
Accordingly, setting ħ=1, the free Hamiltonian of  𝒮 reads 


    H_ S:=H_ S, 0+H_ S, g ,   
       H_ S, 0:=ω_ A a^† a + ω_ B b^† b  ,
         with    ω_ A=ω_ B:=ω_0 ,    
       H_ S, g:= g(a^† b +  h.c.) ,

which can also be conveniently expressed as

    H_ S   =   ω_+ γ_+^†γ_+ + ω_- γ_-^†γ_- ,

with  

    ω_±   :=   ω_0± g  ,     γ_±:=1/√(2) (a± b) ,

being, respectively, the associated eigenmode frequencies and operators
 <cit.>, the last obeying the commutation rules


    [γ_-,γ_+]_-=[γ_-,γ^†_+]_-=0 , [γ_±,γ^†_±]_-=1 .


 Through the exclusive mediation of  subsystem A, we then assume
S   to be connected with an
 external environment E formed
of a collection of a large number M of independent bosonic modes,   
no direct coupling being instead allowed between  B and E. 
Indicating with  c_k, c_k^† the ladder operators of the k-th mode of 
 E,
we hence express the full Hamiltonian of the joint system 𝒮+ℰ as 

    H:=H_ S+H_ E+H_1 ,

 with

    H_ E:=∑_k=1^M ω_k c_k^† c_k ,      
    H_1:=∑_k=1^M γ_k ( a^† c_k +   h.c.) ,

being respectively the free Hamiltonian of the environment 
and 
the exchange coupling between A and E.
More in details, in our analysis
 we shall assume the frequencies ω_k of the environmental modes to be equally spaced  with a cut-off value ω_c> ω_0, i.e. 

    ω_k := k/Mω_ c ,  k ∈{1,..., M} ,

and take  the  system-environment coupling constants γ_k 
to have the form 

    γ_k:= √(κ(ω_0) (ω_k/ω_0)^αω_c/2π M) ,

with  κ(ω_0) controlling the effective strength of the interaction between
A and E. The parameter α≥ 0 appearing in Eq. (<ref>) gauges the 
 bath's dispersion relation by imposing the following form for the (rescaled) spectral density of the reservoirs modes <cit.> 

    κ(ω):=  2 π∑_k=1^M γ_k^2 δ(ω-ω_k)=
    κ(ω_0) (ω/ω_0)^α Θ (ω_ c-ω) ,

with Θ(x) being the Heaviside step function (α=1, α>1 and α<1 being associated to the Ohmic, super-Ohmic, and sub-Ohmic scenarios respectively <cit.>). 
Finally we shall assume the joint S +  E system to be initialized into a factorized
state 

    ρ_ SE(0)=ρ_ S(0) ⊗ρ_ E(0) ,
  where the
bath is in a thermal state of temperature
 1/β>0:


    ρ_ E(0)   :=   e^-β H_ E/ tr[e^-β H_ E] = ρ_1(β)⊗…⊗ρ_M(β) ,
    ρ_k(β)   :=   e^-βω_k c^†_k c_k/ tr[e^-βω_k c^†_k c_k] .
 
























§ APPROXIMATED EQUATIONS FOR S



In this section we review the different ME approaches one can use to effectively describe 
the evolution of the system S by integrating away  the degrees of freedom 
of  the environment E. We shall start our presentation by  
introducing the coarse-grained  regularized version of  the Redfield equation <cit.>, 
which includes  the Global ME as a special case. We then introduce the local ME approach and finally 
discuss the phenomenological approach which employs  convex combinations of local and global 
ME solutions.  Since most of the derivations of the above expressions are 
 discussed in details elsewhere (see e.g. <cit.>) here we just give an overview of the methods involved and refer the interested reader to the Appendix <ref> for further details.



 §.§ From CP-Redfield ME to Global ME

The starting point of this section is the Redfield equation 
 which one obtains by 
expressing the dynamical evolution of the joint system in the interaction picture, and 
enforcing 
the Born and, then, the Markov approximations <cit.>.  The Born approximation assumes that the 𝒮+ℰ coupling is weak in such a way that the state of ℰ is negligibly influenced by the presence of 𝒮, while the Markov approximations assumes invariance of the interaction-picture system state over time-scales of order τ_ E, the last being the time over which ℰ loses the information coming from 𝒮 and can be estimated from width of the bath correlation functions (see Appendix <ref>). 

As anticipated in the introductory section, the Redfield equation does not ensure completely positive evolutions and in certain cases neither positive evolution, hence preventing one from framing the obtained
results with the  probabilistic interpretation of quantum mechanics. 
To cure this issue we refer to the version of the partial secular approximation described in Ref. <cit.>.

Performing a coarse-grain averaging on the Redfield equation in interaction picture over a time interval Δ t that is much larger than the typical time scale of the system state in interaction picture, is a way to appropriately smooth the non-secular terms responsible of the non-positive  character, even in a tight way.
As schematically pictured in Fig. <ref>, 
by moving the parameter Δ t
 along the interval [0,∞[
the reported 
 technique is also capable to formally connect  the original Redfield equation (Δ t=0) and the full secular approximation (Δ t=∞) in a continuous
 way. 
Expressed in Schrödinger picture, 
 the coarse-grained Redfield equation for the evolution of ρ_ S for
 fixed coarse-graining time Δ t, 

reads


    ρ̇_̇ ̇Ṡ(t)=
    -i [H_ S+H_ LS^(Δ t), ρ_ S(t)]_-
    
    
              
    + ∑_σ , σ'=±
    S_σσ'^(Δ t){γ^(1)_σσ'(γ_σ^†ρ_ S(t) γ_σ'
    -1/2[ γ_σ'γ_σ^†,ρ_ S(t)]_+)
    
    
                   +
    
    γ^(2)_σ' σ( γ_σ'ρ_ S(t) γ_σ^†
    -1/2[ γ_σ^†γ_σ' ,ρ_ S(t)]_+)},
 
where hereafter we shall use the symbols [⋯, ⋯]_∓ to represents commutator and  anti-commutator relations, where γ_± are the eigenmode operators of H_ S introduced in Eq. (<ref>), 
and where 

    H_ LS^(Δ t):=
    ∑_σ , σ'=± S_σσ'^(Δ t)
    (η_σσ'^(1)  +
    
    η_σ' σ^(2))
    γ^†_σγ_σ' ,

being the so called Lamb-shift Hamiltonian correction term.

As indicated by the notation,
 the dependence of Eq. (<ref>) 
 upon the coarse-graining time interval Δ t
 is carried out by the tensor S_σσ'^(Δ t) 
of components 
 
    S_σσ'^(Δ t)   :=    sinc((σ-σ')gΔ t2)
    
       =   δ_σσ' + (1-δ_σσ')   sinc(g Δ t) ,


with
sinc(x):= sin(x)/x being the cardinal sinus.
The functional dependence of the right-hand-side of
(<ref>) 
upon the bath temperature
  is instead  carried on by
the tensors
γ^(i)_σσ' and η^(i)_σσ'. 
 Specifically, for 
 σ, σ' ∈{+, -}  and i∈{ 1,2}, 
these elements fulfill the constraints

    γ^(i)_σσ':=γ^(i)_σσ+γ^(i)_σ' σ'/2+
    i (η^(i)_σσ-η^(i)_σ' σ') ,
    
    
    η^(i)_σσ':=-i γ^(i)_σσ-γ^(i)_σ' σ'/4+
    η^(i)_σσ+η^(i)_σ' σ'/2 ,


 which allows one to express all of them in terms of their diagonal (σ= σ')
 components 





    γ_σσ^(1)   :=   1/2 κ(ω_σ)
     𝒩(ω_σ) ,
    
    γ_σσ^(2)   :=   1/2 κ(ω_σ)
     [1+𝒩(ω_σ)] ,
    
    η^(1)_σσ   :=   1/2_0^∞ d ϵ 1/2 πκ(ϵ) 𝒩(ϵ)/ϵ-ω_σ ,
    η^(2)_σσ   :=   
    -1/2_0^∞   d ϵ 1/2 πκ(ϵ) [1+𝒩(ϵ)]/ϵ-ω_σ ,

with  κ(ω) the spectral density of the reservoirs 
defined in Eq. (<ref>) and with 

    𝒩(ω_k):=[  c_k^† c_k ρ_k(β)] =1/e^βω_k-1 ,

being the Bose-Einstein factor of the mode k of the thermal bath.



For  gΔ t→ 0, S_σσ'^(Δ t)  assumes constant  value 1 for all σ and σ': this corresponds to the pathological case of the (uncorrected) Redfield equation in which both the diagonal (secular) and the off-diagonal 
(non-secular) σ, σ'
terms  of the right-hand-side of 
Eq. (<ref>) 
 contribute   at the same level to the dynamical evolution of ρ_ S(t)
 paving the way to unwanted non-positive effects. 
As gΔ t increases 
the off-diagonal component  S_+-^(Δ t)=  sinc(g Δ t)  acts as the smoothing factor for the non-secular (σ≠σ') part of the ME,
which 
gets progressively depressed as 
 the coarse-grain time interval Δ t gets comparable or even larger 
 than the inverse of the  energy scale g of the system.
Following Ref. <cit.>  one can then 
show that the model 
admits a (finite) threshold value for Δ t above which Eq. (<ref>) acquires  
 the explicit  GKSL form that is necessary and sufficient to ensure
complete positivity of the resulting evolution. 
Specifically, as discussed in details  in Appendix <ref>, such threshold is triggered by the inequality

    |S_+ -^(Δ t)|≤min_i ∈{1,2}√(γ_++^(i)γ_–^(i)|γ_+, -^(i)|^2
    ) .

In the following, the  equation (<ref>) at positivity threshold, i.e. with the choice of S_+ -^(Δ t) tightly saturating the bound of Eq. (<ref>), will be called CP-Redfield.

A numerical study of the condition (<ref>)
for some selected values of the system parameters is presented
in  
Fig. <ref>.

This plot makes it clear that the low temperature regime 
(𝒩(ω_0)≪ 1) 
constraints one to take very small values of |S_+ -^(Δ t)| to guarantee the completely positive character of the evolution <cit.>, while just a tiny correction is needed at high temperatures. These  facts are 
in full agreement with the observation  <cit.> that the non-positivity character of the Redfield equation is  enhanced at low temperature as a signature of the deviations from 
the Born-Markov assumptions underlying it <cit.>. We stress that, in this context, non-positivity is originated by the multipartite nature of 𝒮: indeed, as g→ 0, the right-hand side of Eq. (<ref>) tends to 1 and consequently the non-positivity of the Redfield ME disappears in this limit.
Notice finally that irrespectively from the value of g, Eq. (<ref>) is 
trivially fulfilled 
in the asymptotic 
gΔ t
→∞  limit where |S_+ -^(Δ t)| approaches the value zero
leading to 

    S_σσ'^(∞)=δ_σσ' .
 
 This condition identifies 
the full secular approximation
of Eq. (<ref>) 
 that transforms such equation into   the Global ME of the model which, 
 for the sake of completeness, we report here in its explicit form 

    ρ̇_ S(t)=
    -i [H_ S+
    H_ LS^ (glob)
    , ρ_ S(t)]_-
    
    
            + 
    
    ∑_σ=±{1/2κ(ω_σ)𝒩(ω_σ)
    (γ_σ^†ρ_ S(t) γ_σ
    -1/2[ γ_σγ_σ^†, ρ_ S(t)]_+)
    
    
            
    +
    
    1/2κ(ω_σ)[1+𝒩(ω_σ)]
    ( γ_σρ_ S(t) γ_σ^†
    -1/2[ γ_σ^†γ_σ ,ρ_ S(t)]_+)} ,

with

    H_ LS^ (glob)   :=    H_ LS^(∞)=∑_σ= ±δω_σγ^†_σγ_σ ,
    δω_σ   :=   η_σσ^(1)  
    +
    η_σσ^(2)
    =
    1/4 π_0^∞ dϵ
    κ(ϵ)  /ω_σ-ϵ,

being  the secular component of the 
Lamb-shift term <cit.>.



 §.§ Local ME
 


The Local ME for S is a GKSL equation characterized by 
Lindblad operators which act locally on the mode A.
Explicitly it is given by

    ρ̇_̇ ̇Ṡ(t)=
    -i[H_ S+H_ LS^ (loc), ρ_ S(t)]_-
    
               + κ(ω_0)  𝒩(ω_0) 
    (a^†ρ_ S(t) a-1/2[ a a^†, ρ_ S(t)]_+
    ) 
                + κ(ω_0) (1+𝒩(ω_0)) 
     (a ρ_ S(t) a^†-1/2[a^† a, ρ_ S(t)]_+ 
    ) ,

with κ(ω_0) and 𝒩(ω_0)
defined as in the previous section and where now the Lamb-shift term is 
expressed as a modification of the local  Hamiltonian of the A mode only, i.e. 

    H_ LS^ (loc)   :=   δω_ A  a^† a ,
    δω_ A   :=   1/2 π_0^∞ d ωκ(ω)/ω_0-ω.

 Effectively Eq. (<ref>) can be obtained 
 starting  from a  Hamiltonian model for the global system S+ E
 where one initially completely
neglects the presence of the B mode, enforces the  same approximations
 that leads one to 
  (<ref>) (i.e.
the Born, Markov, and full secular approximation),
and finally introduces B and its coupling with A 
as an additive Hamiltonian contribution in the resulting expression. 
More formally as shown e.g. in Ref. <cit.>, Eq. (<ref>)
can be derived in the  weak internal coupling limit g τ_ E≪ 1 (τ_ E being the bath memory time scale, see Appendix <ref> for details) which allows one to treat the interaction between A and B as a perturbative correction with respect to the direct A-E coupling – see Appendix <ref> for
more on this.






 §.§ Convex mixing  of Local and Global solutions
 
As we shall explicitly see in the next section (see Eq. (<ref>)), the main advantage offered by the Global ME (<ref>) is that it provides
an accurate description of the steady state of S  at least in the infinitesimally small S+ E 
 coupling regime where  on pure thermodynamic considerations one expects
 independent thermalization of the eigenmodes γ_± of the system.
On the contrary  the steady state predicted by the Local ME (<ref>) is wrong (even if increasingly accurate as g/ω_0 → 0) because it implies the thermalization of the subsystems A and B regardless of the presence of the internal coupling H_ S, g.
Conversely, the Local ME has the quality to predict Rabi oscillations between A and B at shorter time scales, that are completely neglected when adopting the Global ME.

In view of these observations 
a reasonable way of keeping local effects during the transient still maintaining an accurate steady state solution is to adopt an appropriate phenomenological ansatz
describing the evolution of S in terms of quantum trajectories that interpolate between the solutions ρ^( glob)_ S(t) and ρ^( loc)_ S(t)  of the Global and Local ME, see
Fig. <ref>.
The simplest of these construction is provided by the following 
time-dependent mixture


    ρ^( mix)_ S(t):= 
    e^- 𝒢 tρ^( loc)_ S(t) + (1-e^- 𝒢 t) 
    ρ^( glob)_ S(t) .

In this expression 𝒢>0 is an effective rate, whose inverse fix the time scale of the problem that determines when global thermalization effect start dominating the
system dynamics. 
Accordingly Eq. (<ref>) 
allows us to keep local effects for short time scales t≲𝒢^-1 and the correct thermalization of the eigenmodes of the system at longer time scales t≫𝒢^-1.
The above formula can be interpreted as follows: the environment needs a finite amount of time to become aware of the presence of the part B because of its short time correlations (Markovian hypothesis). 
The specific value of  𝒢 is a free variable in this model and works as a fitting parameter: its value can be even estimated  quite roughly because of the relatively large time interval at intermediate time scales where the Global and Local approximations look alike (more on this later). 
It is finally worth observing that from the complete positivity properties of both the solutions
of the Global and Local ME, it follows that (<ref>) 
also fulfills such requirement (indeed convex combinations of completely positive
transformations are also completely positive). On the contrary 
at variance with the original expressions (<ref>) and (<ref>), as well as the CP-Redfield 
expression (<ref>), Eq. (<ref>) will typically exhibit a non Markovian character and will not  be possible to present it 
in the form of a GKSL differential equation.
This property is a direct consequence of the fact 
that the set of Markovian evolutions is not closed under  convex combinations
 <cit.>. 




§ DYNAMICS
 


In the study of the approximated equations introduced in the previous section, as well as for their comparison with the 
exact  solution of the S+  E dynamics,
an important simplification  arises from  the choice we made in fixing the
initial condition of E.
Indeed thanks to Eqs. (<ref>), (<ref>) 
the  resulting CP-Redfield, Global, and Local MEs, happen to be  Gaussian processes <cit.> which
admit complete characterization only in terms of the first and second moments of the field operators γ_±
 (notice that while the mixture (<ref>) does not
fit into the set of Gaussian processes – formally speaking it belongs to the convex-hull of such set –  we can still resort to the above simplification by
 exploiting the fact that ρ^( mix)_ S(t) is explicitly given by the sum of the Global
and Local ME solutions).
Accordingly in studying the dynamics of our approximated schemes we can just focus on 
the functions
⟨γ_σ⟩(t):=
 [ γ_σρ_ S(t)],
⟨γ_σγ_σ'⟩(t):=
 [ γ_σγ_σ'ρ_ S(t)], and ⟨γ_σ^†γ_σ'⟩(t):=
 [ γ_σ^†γ_σ'ρ_ S(t)] whose temporal dependence can be determined
by solving a restricted set of coupled linear differential equations.
 We also observe that
since the full Hamiltonian (<ref>) conserves the total number of excitations in the S +  E model,  coupling between excitations  conserving and  non-conserving  moments are prevented  <cit.> yielding further
 simplification in the analysis. 





 
 
Having clarified these points, in what follows we shall focus on the special case where 
the input state of S is fixed assuming that 
 both A and B are initialized  in the  ground states of their local Hamiltonians,
i.e. 

    ρ_ S(0)   =   |0⟩_ A⟨0|⊗|0⟩_ B⟨0| ,
 
with |0⟩ representing the zero Fock state of the corresponding mode. 
Under these conditions the input state is Gaussian <cit.> 
and,
evolved under CP-Redfield, Global, Local and
the exact dynamics, will remain Gaussian at all time. Furthermore 
all the first order moments and all the non-excitation-conserving second order
terms exactly nullify, i.e.

    ⟨γ_σ⟩(t) =0  ,     ⟨γ_σγ_σ'⟩(t) =0  ,
 
 leaving only a restricted set of equations to be explicitly integrated. 
In particular, for the case of the coarse-grained Redfield equation
(<ref>) we get 
 

    d/dt⟨γ_+^†γ_+⟩(t)   =   
    -1/2κ(ω_+)[⟨γ_+^†γ_+⟩(t)-𝒩(ω_+)]
    
          +  S_+-^(Δ t)×{2  Im ((η^(1)_+-+η^(2)_-+) ⟨γ_-γ_+^†⟩(t) )
    +
    
    
      Re [(γ^(1)_+--γ^(2)_-+) ⟨γ_-γ_+^†⟩(t)]}
      ,
    
    
    
    
    d/dt⟨γ_-^†γ_-⟩(t)   =   
    -1/2κ(ω_-)[⟨γ_-^†γ_-⟩(t)-𝒩(ω_-)]
    
          +   S_+-^(Δ t)×{-2  Im((η^(1)_+-+η^(2)_-+) ⟨γ_-γ_+^†⟩(t))
    +
    
    
       Re [(γ^(1)_+--γ^(2)_-+) ⟨γ_-γ_+^†⟩(t)]
    } ,
    
    
    
    
    d/dt⟨γ_-γ_+^†⟩(t)   =   {i(ω_+ + δω_+-ω_-  -  δω_-)-1/4[κ(ω_+)+κ(ω_-)]}⟨γ_-γ_+^†⟩(t)
          +   S_+-^(Δ t)×{
    i (η^(1)_-++η^(2)_+-)   
    [⟨γ_-^†γ_-⟩(t)-
    ⟨γ_+^†γ_+⟩(t) ]+
    
    
    γ^(1)_-++1/2 (γ^(1)_-+-γ^(2)_+-) 
    [⟨γ_-^†γ_-⟩(t)+ ⟨γ_+^†γ_+⟩(t)] 
     } ,


with
initial values 


    ⟨γ_+^†γ_+ ⟩(0)   =   ⟨γ_-^†γ_- ⟩(0)=
    ⟨γ_- γ_+^†⟩(0)=0 ,

imposed by  (<ref>).
In particular in the case of full secular approximation (S_+-^(Δ t)=0) the above set of equations become

    d/dt⟨γ_+^†γ_+⟩(t)   =   
    -1/2κ(ω_+)[⟨γ_+^†γ_+⟩(t)-𝒩(ω_+)] , 
    
    
    d/dt⟨γ_-^†γ_-⟩(t)   =   
    -1/2κ(ω_-)[⟨γ_-^†γ_-⟩(t)-𝒩(ω_-)] , 
    
    
    
    d/dt⟨γ_-γ_+^†⟩(t)   =   {i(ω_+ + δω_+-ω_-  -  δω_-)
           -1/4[κ(ω_+)+κ(ω_-)]}⟨γ_-γ_+^†⟩(t) ,

which yield the evolution of the moments for the  Global ME (<ref>). 
Similar considerations hold true for  the Local ME (<ref>). In this case 
following Refs. <cit.> we get 
 

    d/dt⟨γ_+^†γ_+ ⟩ (t)    =    -1/2κ(ω_0) [⟨γ_+^†γ_+ ⟩ (t) - 𝒩(ω_0)+   Re⟨γ_- γ_+^†⟩(t)]  
    + δω_ A  Im⟨γ_- γ_+^†⟩(t)       , 
    
    d/dt⟨γ_-^†γ_- ⟩ (t)    =    
    -1/2κ(ω_0) [⟨γ_-^†γ_- ⟩(t) - 𝒩(ω_0)
    +   Re⟨γ_- γ_+^†⟩(t)] 
    - δω_ A  Im⟨γ_- γ_+^†⟩(t)  
          , 
    
    d/dt⟨γ_- γ_+^†⟩ (t)    =    [i 2 g - 1/2κ(ω_0)] ⟨γ_- γ_+^†⟩(t) + 
    κ(ω_0) /2{𝒩(ω_0)
    
    -1/2   [⟨γ_+^†γ_+ ⟩ (t) + ⟨γ_-^†γ_- ⟩(t)]}
    +i  δω_ A/2
       [ ⟨γ_-^†γ_- ⟩(t)- ⟨γ_+^†γ_+ ⟩ (t)]    
         ,


which, for a direct 
 comparison with Eq. (<ref>),
  we 
 express here in terms of the eigenmodes γ_±.

 
 
 






 §.§ Evolution of the second moments
   

A closer look at Eq. (<ref>) reveals that in this case 
one has that for large enough t we get 

    ⟨γ^†_±γ_±⟩|_   ( glob) (∞)= 𝒩(ω_±) ,     ⟨γ_ -γ^†_ +⟩|_( glob) (∞)=0 .

This enlightens the fact that,
as anticipated at the beginning of Sec. <ref>, the Global ME (<ref>) 
imposes S to asymptotically converge  toward 
 the Gibbs thermal state 

    ρ^( glob)_ S(∞):=e^-β H_ S/ tr[ e^-β H_ S] ,

in agreement with what one would expect  from purely thermodynamics considerations under 
weak-coupling conditions for the system-environment 
interactions.
On the contrary the steady state predicted by the Local ME is wrong (even if increasingly accurate as g/ω_0 → 0) because it implies the thermalization of the subsystems A and B regardless of the presence of the internal coupling H_ S, g. Indeed 
from Eq. (<ref>) we get 

    ⟨γ_±^†γ_±⟩|_( loc)(∞)=𝒩(ω_0) ,    ⟨γ_- γ_+^†⟩|_( loc)(∞)=0
 
or equivalently 

    ⟨ a^† a ⟩|_( loc)(∞)   =   ⟨ b^† b ⟩|_( loc)(∞)=𝒩(ω_0) ,
    ⟨ a b^†⟩|_( loc)(∞)   =   0 ,

 which identifies 


    ρ^( loc)_ S(∞):=e^-β H_ S,0/ tr[ e^-β H_ S,0] ,

 as the new fixed point for the dynamical evolution


  (see also Appendix <ref>).
The discrepancy between the above expressions and
 Eqs. (<ref>), (<ref>) is even accentuated in the low temperature regime   
 βω_0 ≫ 1,
where in particular the ratio 
𝒩(ω_-)/𝒩(ω_0) ≃ e^β g explodes exponentially. 


 
The situation gets reversed  at shorter time scales. Here 
the Local ME correctly presents coherent energy exchanges between A and B 
which instead the Global approach completely neglects. 
Indeed from Eq. (<ref>) 
 it follows that the Global ME predicts Im [⟨ ab^†⟩(t)]=0, the term  
 being responsible of the Rabi oscillations between A and B.
The Local ME on the contrary 
– when the Lamb-shift correction can be neglected – 
gives Re[ ⟨ ab^†⟩(t)]=0, the latter  being proportional to the average internal interaction energy ⟨ H_ S, g⟩.






The above observations are confirmed by the numerical study we present  in the remaining of the section
(see however also the material presented in Appendix <ref>).
In particular, in panels (a) and (b) of Fig. <ref> the temporal evolution of the second order moments obtained
by solving Eq. (<ref>) and (<ref>) are compared with 
 the exact values of the corresponding quantities obtained by numerical
integration of the exact S + E Hamiltonian model along the lines detailed
in Appendix <ref>.
In panel (c) of such figure we also  present 
 the results obtained by using the effective model of Sec. <ref>, where according to Eq. (<ref>)
 the expectation values of the relevant quantities are computed as 

    ⟨γ_σ^†γ_σ'⟩|_( mix)(t)
       =    e^- 𝒢 t⟨γ_σ^†γ_σ'⟩|_( loc)(t)
            + (1-e^- 𝒢 t) 
    ⟨γ_σ^†γ_σ'⟩|_( glob)(t) ,
 
with ⟨γ_σ^†γ_σ'⟩|_( loc)(t) and 
⟨γ_σ^†γ_σ'⟩|_( glob)(t) representing  the solutions of
Eq. (<ref>) and 
Eq. (<ref>) respectively. 
In our analysis the system  parameters  have been set in order to enforce 
 S +  E weak-coupling conditions
(ω_0,ω_±≫κ(ω_0))  to make sure that 
that 
the long term prediction 
(<ref>) of 
the Global ME provides a proper description of the system dynamics.
By the same token, the temperature of the bath  has been fixed to be relatively high,
i.e. 1/β≈10.5 ω_0, to avoid 
 to enhance 
 correlation effects between the bath and the system which are not included in 
 the Born and Markov approximations needed to derive both the Global and the Local ME <cit.> (a study of the impact of
 low temperature effects on the S +  E correlations is presented in  Appendix <ref>). Finally regarding the value of the phenomenological parameter G entering in (<ref>) we set it be equal to 0.4 κ(ω_0)
finding a relatively good agreement with the exact data at all times. 






 
The convex combination (<ref>) is not the only way of keeping the best from both the local and the global approximations.
Indeed, by making a step back, one can consider the coarse-grained Redfield equations (<ref>) once that the pathology related to their non-positivity has been cured.
A detailed study of the performances of this approach is presented in Fig. <ref>. Here, for the same values of the parameters used in Fig. <ref>, in panel (a) we exhibit the plots associated with the CP-Redfield equation obtained  by fixing S_+ -^(Δ t)
in such a way to saturate the positivity bound (<ref>), i.e.
S^(Δ t)_+-=0.989.
As in the case of  panel (c)  of Fig. <ref>, we notice that
CP-Redfield is in a good agreement with the exact data both at long and short
time scales. As a check in panel (b) of Fig. <ref>
we also present the
(uncorrected)
Redfield equation obtained by setting in Eq. (<ref>) Δ t=0, corresponding to have S_+ -^(Δ t)=1 which for the system parameters we choose
 gives a clear violation of the positivity bound (<ref>).
Interestingly enough, despite the fact that
the resulting equation does not guarantee complete positivity of the associated
evolution, we notice that also in this case one has an apparent good agreement with the exact results for all times. 

In particular both CP-Redfield and Redfield equations appear to be able to capture a non-weak coupling correction to the asymptotic value of 2  Re⟨γ_- γ_+^†⟩ (t)=
⟨ a^† a⟩(t)-⟨ b^† b⟩(t) ,
an effect that is present in the exact model due to the fact that the
subsystem A remains slightly correlated with the bath degrees of freedom,  but which
is not present when adopting neither Global,  Local, or Mixed approximations (see Fig. <ref>). 

An evidence of this can be
obtained by observing that  from Eq. (<ref>) we have


    2  Re⟨γ_- γ_+^†⟩ (∞)=
    S_+ -^(Δ t)ω_+ -ω_-
    
    
    
    
    
    
    
    
    
          ×_0^∞ dϵ κ(ϵ)/2 π (𝒩(ϵ)-𝒩(ω_+) ϵ- ω_+ - 𝒩(ϵ)-𝒩(ω_-) ϵ- ω_-)+O[κ(ω_0)^2]


which is exactly null for the Global ME (S_+ -^(Δ t)=0), but which 
 is different from zero (and in  good agreement with the exact result) both
 for the uncorrected Redfield equation (S_+ -^(Δ t)=1) and 
 CP-Redfield (S^(Δ t)_+-=0.989).








Despite the  apparent success of the uncorrected Redfield equation reported above, a clear signature of its non-positivity can still be spotted by looking 
at a special functional  of the second order moments of the model, i.e. the quantity 
    λ_c (t) :=   1/2 min{ eigenvalues [ Γ_ S(t)+i Ξ_ S ] } .


In the above  definition 
 Γ_ S(t) and Ξ_ S are 
 respectively  the covariance matrix and the symplectic form of the two-mode system S.
 Expressed in terms of the eigenoperators γ_± their elements are given by 
 
    [Γ_ S(t)]_ij:= ⟨[ Γ_i -⟨Γ_i⟩(t),Γ_j^† -⟨Γ^†_j⟩(t)]_+⟩(t),

and 
 
    [Ξ_ S]_ij :=  -i ⟨[ Γ_i, Γ_j^†]_- ⟩(t) =
    -i 
    [  1  0  0  0;  0 -1  0  0;  0  0  1  0;  0  0  0 -1 ],

 with Γ_i being the i-th component of the operator vector Γ:=(γ_+,γ_+^†,γ_-,γ_-^†)^T. 
 In particular due to the choice of the input state we made in Eq. (<ref>),
 we get 
  
 
    Γ_ S(t) =
    [ 2 ⟨γ_+^†γ_+ ⟩ (t)+1                   0  2 ⟨γ_- γ_+^†⟩(t)^*                   0;                   0 2 ⟨γ_+^†γ_+ ⟩ (t)+1                   0    2 ⟨γ_- γ_+^†⟩(t);    2 ⟨γ_- γ_+^†⟩(t)                   0 2 ⟨γ_-^†γ_- ⟩ (t)+1                   0;                   0  2 ⟨γ_- γ_+^†⟩(t)^*                   0 2 ⟨γ_-^†γ_- ⟩ (t)+1 ] ,

and hence
 
    Γ_ S(t) + i Ξ_ S/2
     =
    [ ⟨γ_+^†γ_+ ⟩ (t)+1                 0  ⟨γ_- γ_+^†⟩(t)^*                 0;                 0   ⟨γ_+^†γ_+ ⟩ (t)                 0    ⟨γ_- γ_+^†⟩(t);    ⟨γ_- γ_+^†⟩(t)                 0 ⟨γ_-^†γ_- ⟩ (t)+1                 0;                 0  ⟨γ_- γ_+^†⟩(t)^*                 0   ⟨γ_-^†γ_- ⟩ (t) ] ,
    λ_c(t)=1/2{⟨γ_+^†γ_+ ⟩(t)+
    ⟨γ_-^†γ_- ⟩(t)-
    √([⟨γ_+^†γ_+ ⟩(t)-
    ⟨γ_-^†γ_- ⟩(t)]^2+4 |⟨γ_- γ_+^†⟩(t)|^2)} .

  
  When evaluated on a proper state of the system, the  
 Robertson-Schrödinger uncertainty relation inequality <cit.> 
 forces  the spectrum of the above matrix to be non-negative –
 see Appendix <ref> for details.
Accordingly 
 when  ρ_ S(t) is positive semi-definite (i.e. it is a physical state) one must have λ_c (t)≥ 0.
The temporal evolution  of λ_c (t) is reported in 
 Fig. <ref> for the various approximation methods and for the exact dynamics: one notice that while Global, Local, and CP-Redfield always complies with 
 the positivity requirement, the uncorrected 
 Redfield equation exhibit negative values of λ_c(t) at short time scales.

Analytically, this can  
 be seen from the short time scale trend of λ_c (t) , which from 
 Eq. (<ref>) can be determined as


    λ_c(δ t)   ≃   (γ_–^(1)+γ_++^(1)2)
    [1-
          √(1+
    4
    (S_+ -^(Δ t)^2-
    γ_++^(1)γ_–^(1)|γ_+-^(1)|^2)
    |γ_+-^(1)|^2(γ_–^(1)+γ_++^(1))^2) ] δ t ,


which tightly gives λ_c(δ t)≥  0 if and only if the complete positivity constraint (<ref>) is fulfilled.

 Notice also that while  none of the approximated methods
 are able to follow the whole exact behaviour λ_c(t), CP-Redfield and Global provide good agreement in the long time limit, while CP-Redfield and Local correctly 
predict λ̇_̇ċ(0)=0. 



























 












 §.§ Fidelity Comparison
 

In this section we further discuss the difference between the various
approximation methods, as well as their relation with the exact solution,
evaluating the temporal evolution of the Uhlmann fidelity <cit.> between the associated
density matrices of S. 
We remind that given ρ^(1)_ S and ρ^(2)_ S two quantum states of the system their fidelity
is defined as the positive functional  

    ℱ(ρ^(1)_ S, ρ^(2)_ S):=√(ρ^(1)_ S)√(ρ^(2)_ S)_1 ,

with Θ_1 := [ √(Θ^†Θ)] being the trace norm
of the operator Θ. This quantity provides a bona-fide estimation of how close
the two density matrices are, getting its maximum value 1 when ρ^(1)_ S=ρ^(2)_ S,
and achieving zero value instead  when the support of
ρ^(1)_ S and ρ^(2)_ S  are orthogonal, i.e. when they are perfectly distinguishable. 
In the case of two-mode Gaussian states <cit.> with null first order moments, a relatively simple closed expression for ℱ(ρ^(1)_ S, ρ^(2)_ S) is known in terms of the covariance matrices of the two density matrices  <cit.>.
Specifically, in the eigenmode representation, one has 

    ℱ^2(ρ^(1)_ S, ρ^(2)_ S)=1/√( b)+√( c)-√((√( b)+√( c))^2- a) ,


with 

    a   :=   2^-4 [Γ^(1)_ S+Γ^(2)_ S] ,
    
    
     b   :=   2^-4[  Ξ_ S Γ^(1)_ S Ξ_ S Γ^(2)_ S -   1_4 ]
     ,
    
     c   :=    2^-4[Γ^(1)_ S+i Ξ_ S][Γ^(2)_ S+i Ξ_ S]
     ,

where  Γ^(1)_ S, Γ^(2)_ S being the covariance matrices of ρ^(1)_ S and ρ^(2)_ S defined in 
 (<ref>), and with 
Ξ_ S the
symplectic form given in  Eq. (<ref>) – see final part of Appendix <ref> for details.
In what follows we shall make extensive use of the
identity (<ref>)
 thanks to the fact that 
for the input
state (<ref>) we are considering in our analysis, 
the density matrix  of S remains Gaussian at all times when evolved 
under 
Global, Local, CP-Redfield ME, as well as under the exact integration of the full
S + E Hamiltonian model. The same property unfortunately does not hold for the convex mixture (<ref>) which is explicitly non-Gaussian (indeed it is a convex combination of Gaussian states). In this case
hence the result of  <cit.> can not be directly applied to
compute ℱ(ρ^( mix)_ S(t) ,ρ^( exact)_ S(t)).
Still  the concavity property <cit.> of the ℱ 
can be invoked to compute the following  lower bound  

    ℱ(ρ^( mix)_ S(t) ,ρ^( exact)_ S(t))
    ≥
    e^- 𝒢 tℱ(ρ^( loc)_ S(t) ,ρ^( exact)_ S(t))
    
                       +
    (1-e^- 𝒢 t) ℱ(ρ^( glob)_ S(t) ,ρ^( exact)_ S(t)),
  
with  the  right-hand-side  being provided by Gaussian terms.
Finally 
the non-positivity of the (uncorrected) Redfield equation also gives rise to problems in the evaluation of the associated fidelity (as a matter of fact, in this case 
the quantity ℱ(ρ^( red)_ S(t) ,ρ^( exact)_ S(t)) is simply ill defined). 
Aware of this fundamental limitation, but also of the  fact that the departure from the
positivity condition  of the solution ρ^( red)_ S(t) of the Redfield equation 
is small, in our analysis we decided to present the real part of ℱ^2(ρ^( red)_ S(t) ,ρ^( exact)_ S(t)). 


 
To begin, in Fig. <ref> we present the value of
 ℱ^2(ρ^( loc)_ S(t), ρ^( glob)_ S(t)): as clear from the plot, this quantity 
is sensibly different from 1  at short and at long time scales (confirming the observation
of the previous section) while it is ∼ 1 at intermediate time scales.
In Fig. <ref> instead we proceed with the comparison of the
approximate solutions with the exact one. The reported plots  
confirm that the convex combination of the local and global solutions (<ref>) is an effective ansatz to approximate the system evolution, giving a (lower) bound for the fidelity computed as in Eq. (<ref>)
that is close to 1 both at short and at long time scales.
On the same footing we find the CP-Redfield equation which still remaining positive brings all the main qualities of the (full) Redfield ME.
For completeness,  in Fig. <ref> we report two situations in which the Global ME  and the local ME work extremely bad respectively. 
In Panel (a) we consider weaker internal coupling g such that the local ME gives a satisfying result for the whole dynamics while the inadequacy of the Global ME during the transient is accentuated;

In Panel (b) we decrease instead the temperature accentuating the inadequacy of the local ME in the steady prediction.

In both the Panels we report the curve corresponding to the CP-Redfield approximation. The last 
follows either the local or the global curve depending on which one performs better in the two instances.









 



 

















§ CONCLUSIONS





In the study of multipartite Markovian open quatum systems it has been widely discussed in literature whether the local dissipator or the global dissipator is more adapt to effectively  reproduce the system dynamics <cit.>.

Here we have treated a case where the system is composed of two interacting harmonic oscillators A and B, with only A interacting with a thermal bath - collection of other harmonic oscillators - and we have analyzed the equilibration process of the system initially in the ground state with the finite bath temperature.

We have shown that the  “completely positive Redfield” equation — i.e. the cured version of the Redfield equation by means of coarse-grain averaging <cit.> — and an appropriate time-dependent convex mixture of the local and global solutions (<ref>) give rise to the most accurate semigroup approximations of the exact system dynamics, both during the time transient and for the steady state properties, going beyond the pure local and global approximations.

The convex mixture of the local and global channels has been introduced phenomenologically for allowing at the same time coherent local energy exchange at short time scales between A and B and the steady state  expected from the thermodynamics at long time scales, i.e. the global Gibbsian state. 

Future developments on this route may concern the search of a microscopic derivation of this (non-Markovian) quantum channel.

D. F and V. G. acknowledge support by MIUR via PRIN 2017 (Progetto di Ricerca di Interesse Nazionale): project QUSHIP (2017SRNBRK). 


















§ DERIVATION OF THE COARSE-GRAINED REDFIELD AND LOCAL ME



In this section we provide details about the derivations of the coarse-grained Redfield (<ref>) and Local (<ref>) MEs.
For (<ref>) we make use of Refs. <cit.> and of the method to correct the non-positivity of the Redfield equation given in Ref. <cit.>, while for (<ref>) we follow the approach of Ref. <cit.>. 



Expressed in interaction picture the evolution of the joint state of S +  E induced by the
 Hamiltonian (<ref>)  is given by the 
 
 Liouville-von Neumann equation
 


    ρ̇̃̇_ SE(t)   =   -i [  H̃_1(t), ρ̃_ SE(t)]_- ,
 
where given U_ 0(t):= e^i (H_ S+H_ E) t we have 

    H̃_1(t)   :=   U_ 0 (t) H_1  U_ 0^† (t) =a^†(t) C(t)+ h.c.









with a^†(t):= e^i H_ S,gt a^† e^-i H_ S,gt and 
C(t):= ∑_k γ_k c_k e^-i(ω_k - ω_0)t.



Tracing out the environment degrees of freedom, Eq. (<ref>) can be written as



    ρ̇̃̇_̇̃̇ ̇̃̇Ṡ̃̇(t)   =   -i  Tr_ E[ H̃_ 1(t), ρ̃_ SE(0)]_-
    
          -∫_0^t  Tr_ E[ H̃_ 1(t),[ H̃_ 1(t'), ρ̃_ SE(t')]_-]_-dt' 
     .
 
We assume now weak system-environment coupling such that the environment stays in its own Gibbs state (<ref>)  (invariant in interaction picture) for all the system dynamics and the SE state can be approximated by the tensor product 


    ρ̃_ SE(t)≃ρ̃_ S(t) ⊗ρ_ E(0).











Equation (<ref>) means that the environment, being a macroscopic object, can be considered insensitive to the interaction with the system (Born approximation <cit.>). On the contrary the system state is affected by the coupling with the environment. Being the first moments null over a thermal state, the first commutator in (<ref>) is zero and by inserting the tensor product (<ref>) such equation becomes


    ρ̇̃̇_ S(t)≃
          ∫_0^t dt'  c^(1)(t-t')
    
    (a^† (t') ρ̃_ S(t') a(t)-
    a(t) a^†(t') ρ̃_ S(t'))+
    
          
    c^(2)(t-t') (a(t') ρ̃_ S(t') a^†(t) -
    a^†(t) a(t') ρ̃_ S(t')) + h.c. ,


where c^(1)(τ) and c^(2)(τ) are bath correlation functions defined as


    c^(1)(τ)   :=   <C^†(τ) C >=∑_k γ_k^2 𝒩(ω_k) e^i (ω_k-ω_0)τ , 
    
    
    
    c^(2)(τ)   :=   <C(τ) C^†>=∑_k γ_k^2 [1+𝒩(ω_k)] e^-i (ω_k-ω_0)τ .


Next step is the Markovian assumption τ_ E≪δ t,  where δ t is the typical time scale of the state in interaction picture and τ_ E is the bath memory time scale, i.e. the characteristic width of the bath correlation functions (<ref>).

Such time scale separation allows to replace in Eq. (<ref>) the upper integration bound with +∞ and to neglect the τ:= t-t' dependence of the state ρ̃_ S, 

leading to the Redfield equation (interaction picture):


    ρ̇̃̇_ S(t)≃∫_0^∞ dτ[ c^(1)(τ)
    
              ×(a^† (t-τ) ρ̃_ S(t) a(t)-
    a(t) a^†(t-τ) ρ̃_ S(t))
    
    
          +
    
    c^(2)(τ) (a(t-τ) ρ̃_ S(t) a^†(t) -
    a^†(t) a(t-τ) ρ̃_ S(t)) ]+ h.c.
 

As described in Ref. <cit.>,
if the bath correlation functions are narrow enough with respect to the internal coupling time scale, i.e. g τ_ E≪ 1,
in Eq. (<ref>)
one can approximate a(t-τ)≈ a(t) 

obtaining the interaction picture version of the Local ME (<ref>), which is in Lindblad form without the need of any secular approximation. Alternatively, passing to the eigenmode basis of Eq. (<ref>), 
Eq. (<ref>) can be equivalently written as  


    ρ̇̃̇_ S(t)   =   
    1/2∑_σ, σ'[ 
    Ω^(1)_σ 
    e^i (σ-σ')g t(γ_σ^†ρ̃_ S(t) γ_σ'-
    γ_σ'γ_σ^†ρ̃_ S(t)) 
    
          +
    Ω^(2)_σ'
    e^i (σ-σ')g t(
     γ_σ'ρ̃_ S(t) γ_σ^† -
    γ_σ^†γ_σ'ρ̃_ S(t)) ]
    +
     h.c.
 


where


    Ω^(1)_σ   :=   ∫_0^∞ dτ  c^(1)(τ) e^-i σ g τ , 
    Ω^(2)_σ'   :=   
    ∫_0^∞ dτ c^(2)(τ) e^i σ' g τ .
 
Next step is to

perform a coarse-grain average on Eq. (<ref>) over a time interval Δ t≪δ t, 

which amounts in applying the following substitution 


    e^i (σ-σ')g t   ⟶   1/Δ t∫_t-Δ t/2^t+Δ t/2 ds  e^i (σ-σ')g s
       =   e^i (σ-σ')g t sinc((σ-σ')g Δ t/2),


without affecting the system state in interaction picture.
Equation (<ref>) is eventually obtained by passing to the Schrödinger picture. 

Indeed the Lamb-shift and the dissipator coefficients of Eqs. (<ref>) and 
(<ref>) are related to the quantities Ω^(i)_σ as 


    γ^(i)_σσ'=1/2 
    (Ω^(i)_σ+Ω^(i)_σ'^*) ,




    η^(i)_σσ'=1/4i 
    (Ω^(i)_σ-Ω^(i)_σ'^*) .
















§ COMPLETELY POSITIVE MAP REQUIREMENT FOR THE COARSE-GRAINED REDFIELD EQUATION
 
To discuss the complete positivity condition for the coarse-grained Redfield equation let us observe that 
its  dissipator is given by  
the last two
lines in the right-hand-side of  Eq. (<ref>). 
Following  Ref. <cit.> we 
write them as

    ∑_i, σ, i', σ'γ_i' σ' , i σ(𝒜_i', σ'^†ρ_ S(t) 𝒜_i, σ
    - 1/2[ 𝒜_i, σ𝒜_i', σ'^† , ρ_ S(t) ]_+ ),


with 𝒜_1, σ=γ_σ, 𝒜_2, σ=γ_σ^†,
and  γ_i' σ' , i σ being the  elements of the  
4 × 4 hermitian matrix 


    γ_I , J  
    =
    
     [            γ^(1)_+ + γ^(1)_+ -S_+ -^(Δ t)                    0                    0; γ^(1)_- +S_+ -^(Δ t)            γ^(1)_- -                    0                    0;                    0                    0            γ^(2)_+ + γ^(2)_+ -S_+ -^(Δ t);                    0                    0 γ^(2)_- +S_+ -^(Δ t)            γ^(2)_- -;                      ].

Complete positivity of  the evolution described by Eq. (<ref>) can now be guaranteed by the imposing positiveness of the spectrum 
of  (<ref>), a condition which  by explicit diagonalization  leads to 
Eq. (<ref>).



§ COVARIANCE MATRICES
 

Expressed in terms of the  system canonical coordinates 
 
    x_ A:=(a+a^†)/√(2) ,            p_ A:=(a-a^†)/(√(2) i) , 
     
    x_ B:=(b+b^†)/√(2) ,            p_ B:=(b-b^†)/(√(2) i) ,
 
the covariance matrix Σ_ S  associated with the quantum state ρ_ S of the two-mode system S is defined as the 4× 4 real hermitian  matrix
 
    [Σ_ S]_αβ:= ⟨[ r_ S,α - ⟨r_ S,α⟩ , 
    r_ S,β  - ⟨r_ S,β⟩]_+
    ⟩ ,

 where  as usual we adopt the shorthand notation ⟨⋯⟩ : =  [⋯ρ_ S],  and where r_ S,α is the α-th component of the operator vector 
 r_ S:=(x_ A, p_ A, x_ B, p_ B)^T.
 In this notation 
 the symplectic form of the system is defined by the matrix   Ω_ S of elements

    Ω_ S:=
    [  0  1  0  0; -1  0  0  0;  0  0  0  1;  0  0 -1  0 ] ,

which embodies the 
canonical  commutation rules of the model via the identity 

⟨ [r_ S,α,r_ S,β]_- ⟩=i [Ω_ S]_αβ.
From the Robertson-Schrödinger uncertainty relations <cit.> it hence follows, 
that for all choice of  ρ_ S we must have that the matrix 
Σ_ S+i Ω_ S is non-negative or equivalently that the following inequality must hold

    min{ eigenvalues [ Σ_ S+i Ω_ S ] }≥ 0 .

Equation (<ref>) is at the origin 
 of the study we presented in Fig. <ref>.
We notice indeed that introducing
 the unitary matrix 
 
    𝒱 := 1/2[   1   1   1   1;  -i   i  -i   i;   1   1 - 1  -1;  -i   i   i  -i ] ,

from  Eq. (<ref>)  the following identity holds, 
 
    r_ S   =    VΓ_ S ,
 
 with Γ_ S the operator vector introduced in Eq. (<ref>),
 which in turns implies
 
    Σ_ S =  VΓ_ S V^† ,     Ω_ S =  VΞ_ S V^† ,
 
    with Ξ_ S as in Eq. (<ref>). 
Accordingly, we get 

    Σ_ S+i Ω_ S =  V( Γ_ S+i Ξ_ S)  V^† ,
 
which finally allows us to translate Eq. (<ref>) into the positivity condition for the quantity
 λ_c(t) introduced in Eq (<ref>). 


Notice finally that the unitary relations (<ref>) are also at the origin of 
Eqs. (<ref>) and (<ref>) which we derived from <cit.>
via the
 identities

    [Γ^(1)_ S+Γ^(2)_ S]    =    [Σ^(1)_ S+Σ^(2)_ S]  , 
    [  Ξ_ S Γ^(1)_ S Ξ_ S Γ^(2)_ S -   1_4 ]   =   [  Ω_ S Σ^(1)_ S Ω_ S Σ^(2)_ S -   1_4 ] , 
    [Γ^(j)_ S+i Ξ_ S]   =   [Σ^(j)_ S+i Ω_ S] ,
 
where for j=1,2,  Γ^(j)_ S and Σ^(j)_ S represent the covariance matrices
(<ref>) and (<ref>) of the matrices ρ_ S^(j).



§ THE EXACT MODEL



In this section, following a procedure similar to <cit.>, we discuss  how to explicitly solve the exact dynamics of the Hamiltonian model for the joint system S +  E.


Passing to the canonical variables  of the full model, i.e. introducing the operators
x_ A=(a+a^†)/√(2), p_ A=(a-a^†)/(√(2) i),
x_ B=(b+b^†)/√(2), p_ B=(b-b^†)/(√(2) i) as in Eq. (<ref>) and
x_k=(c_k+c_k^†)/√(2), p_k=(c_k-c_k^†)/(√(2) i),
the Hamiltonian (<ref>) of 𝒮+ E can be written as 


    H=1/2r^T ℋr +  const .


The vector operator r is the generalization of r_ S introduced in Sec. <ref>
that now contains the canonical coordinates of all the S +  E modes, i.e. 


    r   =   (x_ A, p_ A, x_ B, p_ B, x_1, p_1, ... , x_M, p_M)^T ,


and ℋ is a real symmetric (2 M + 4 )× (2 M + 4 )  matrix, having non null elements only on the diagonal and on the first two  rows and on the first two  columns. This is because only the sub-system A is microscopically attached to the thermal bath:


    ℋ=
    [ ω_ A    0    g    0  γ_1    0       γ_M    0;    0 ω_ A    0    g    0  γ_1         0  γ_M;    g    0 ω_ B    0    0    0         0    0;    0    g    0 ω_ B    0    0         0    0;  γ_1    0    0    0  ω_1    0         0    0;    0  γ_1    0    0    0  ω_1         0    0;    ⋮    ⋮    ⋮    ⋮    ⋮    ⋮    ⋱    ⋮    ⋮;  γ_M    0    0    0    0    0       ω_M    0;    0  γ_M    0    0    0    0         0  ω_M;      ] .


Exploiting the above construction the expectation value of r can now be shown to  evolve in time as <cit.>


    ⟨r(t) ⟩ : = [ rρ_ SE(t)] =e^Ωℋt⟨ r(0) ⟩ ,

where Ω is the symplectic form of the entire model, i.e. the (2 M + 4 )× (2 M + 4 ) matrix 

    Ω:= ⊕_i=1^M+2[  0  1; -1  0 ] ,
 
whose elements embody  the canonical  commutation rules of entire S +  E system via the identity 

⟨ [r_α,r_β]_- ⟩ =i Ω_αβ.
Similarly  the covariance matrix of elements 


    Σ_αβ(t)   :=   [ [ r_α - ⟨r_α(t) ⟩ , 
    r_β - ⟨r_β(t) ⟩]_+ ρ_ SE(t)]  
      
       =   ⟨[ r_α(t) - ⟨r_α(t)⟩ , 
    r_β(t) - ⟨r_β(t)⟩]_+⟩,


can be shown to evolve as 


    Σ(t)= e^ΩℋtΣ(0) e^ℋΩ^Tt .

For future reference it is worth stressing  that the 4 × 4 principal minor of the  matrix Σ(t) 
(i.e. the sub-matrix obtained from the latter by taking the upper left 4 × 4 part) 
corresponds to the covariance matrix Σ_ S(t)  of the S  system alone, 
whose elements can be formally expressed as in Eq. (<ref>). 

 

In the evaluation of    Eqs. (<ref>), (<ref>) one can  resort to the exact diagonalization of the Hermitian matrix ℳ
defined as


    ℳ:= i Ωℋ .


Calling (g_1, , g_2M + 4) the eigenvalues of ℳ and


    V_αβ:=[𝗀^(β)]_α


the unitary matrix whose columns are the normalized eigenvectors 𝗀^(α)  corresponding to the eigenvalues 
g_α, 


the diagonal form of the matrix ℳ is obtained as:


    diag(g_1, , g_2M + 4)=
    V^†ℳ V .
 

Accordingly, we can now  rewrite Eqs. (<ref>), (<ref>) in the form


    ⟨r(t)⟩   =    V E_-(t) V^†⟨r(0)⟩
    Σ(t)    =    V E_-(t) V^†Σ(0) V E_+(t) V^† ,


with 


    E_∓(t)=diag(e^∓ i g_1 t, , e^∓ i g_2M + 4 t) .


In summary, the exact dynamics is obtained thanks to the numerical diagonalization of the matrix
ℳ of Eq. (<ref>) and by performing the matrix multiplications in Eq. (<ref>). 


Regarding the initial conditions, we observe that in the case of the input state we have selected in Eqs. (<ref>), (<ref>) and 
 (<ref>), the initial covariance matrix reads as the direct sum



    Σ(0)=
    [               1_2                 0                 0                                   0;                 0               1_2                 0                                   0;                 0                 0   [2 𝒩(ω_1)+1]1_2                                   0;                 ⋮                 ⋮                 ⋮                 ⋱                 ⋮;                 0                 0                 0                   [2 𝒩(ω_M) + 1]1_2 ]



with 1_2 being the 2× 2 identity matrix and 𝒩(ω_k) being the Bose-Einstein mean occupation numbers
introduced in Eq. (<ref>).

Regarding the first order moments instead, since  ⟨r(0) ⟩=0 the evolution law of  Eq. (<ref>) leads to ⟨r(t) ⟩=0 for all t≥ 0. 


 







 §.§ Memory and recurrence time scales


 







When resorting to numerical methods in solving the exact Hamiltonian model one should be aware of the fact that since it
 involves a finite number of parties (i.e. the system modes A and B and
the M environmental modes), it will be characterized by a recurrence time scale T_ rec that, due to the various approximation involved in their
derivation, leave no trace in the corresponding ME expressions. 
An estimation of such quantity   can be retrieved directly from the periodicity of the correlation functions  of Eq. (<ref>) which leads
us to (see Fig. <ref>a):


    T_ rec=2 π M /ω_c .


The choice of the parameters ω_c=3 ω_0 and M≈ 400 <cit.> ensures that the discretization  does not play any role in the time window we have considered for all the plots.


 
The width of the correlation functions (<ref>)  also  plays an important role in the model: 
it yields the time  τ_ E which takes for the information that emerges from the system 
to get lost into the environment and never coming back <cit.>.

Such time scale can't be resolved by any approximation we have discussed so far, because of the Markovian assumption which is present in all of them.
The estimation of this time scale is given by the half width at half maximum (see Fig. <ref>b) of  |c^(1)(τ)| and |c^(2)(τ)|. 
For  𝒩(ω_0)=10 we get 


    τ_ E≈ 3.8 / ω_c .







 §.§ Low temperature effects
  

It is well known that  in the low temperature regime 
correlation effects between the bath and the system tent to arise,  
challenging the 
Born approximation used in the derivation of the Markovian MEs
 <cit.>.
An evidence of this fact is presented in  Fig. <ref>
where the time evolution of 
 the average components of the Hamiltonian (<ref>) 
 are presented for two different choices of the parameter 1/β. 








§ ON THE THERMALIZATION OF THE SYSTEM EIGENMODES


 

We show here the dual counterparts of the moments reported in Figs. <ref> and <ref> in the basis of the eigenmodes (<ref>),

making clearer when these eigenmodes reach the correct thermalization or not depending on the implemented approximation.

The second order moments in the  a, b basis and the ones in the  γ_+, γ_- basis are related each other as 


    1/2( ⟨ a^† a⟩-⟨ b^† b⟩)
       =    Re⟨γ_- γ_+^†⟩ ,
    
     Im ⟨ a b^†⟩   =    Im ⟨γ_- γ_+^†⟩ ,
    
     Re ⟨ a b^†⟩   =   1/2( ⟨γ_+^†γ_+⟩-⟨γ_-^†γ_- ⟩) ,
    
    
    ⟨ a^† a⟩+⟨ b^† b⟩   =   ⟨γ_+^†γ_+⟩+⟨γ_-^†γ_- ⟩ .


The steady state (<ref>) is what one expects from thermodynamics. It implies

⟨γ^†_±γ_±⟩ (∞)= 𝒩(ω_±) , 
⟨γ_ -γ^†_ +⟩ (∞) =0 .

This result is captured by applying the global approximation (see Eqs. (<ref>)), which under the initial conditions (<ref>) gives 


    ⟨γ_-γ_+^†⟩|_( glob)(t)   =   0,
    ⟨γ_±^†γ_±⟩|_( glob)(t)   =   𝒩(ω_±)(1-e^-1/2κ (ω_±) t) .
 

On the other hand, the local approximation fails just about the steady state properties.

As discussed in <cit.>, under the same initial conditions and when the Lamb-shift correction δω_ A can be neglected,  the local ME (see Eqs. (<ref>)) leads to


    Re⟨γ_- γ_+^†⟩|_( loc)
    (t)=  𝒩(ω_0) e^-  κ(ω_0) t/2/ϵκ(ω_0)   sin(ϵ t/2)  
         ,
    
           Im⟨γ_- γ_+^†⟩|_( loc)(t)=
    4 𝒩(ω_0) κ(ω_0) g
    e^-κ(ω_0) t/2/ϵ^2 
    [1-cos(ϵ t / 2)] , 
    
          ⟨γ_±^†γ_±⟩|_( loc)(t)
    = 𝒩(ω_0) 
              ×{ 1-e^-  κ(ω_0) t/2/ϵ^2[16 g^2-κ(ω_0)^2 cos(ϵ t/2)]} ,

with 
ϵ:= √((4 g)^2-κ(ω_0)^2).

Using the relations (<ref>-<ref>), the above equations imply in the a, b basis:


    ⟨ a^† a ⟩|_( loc)
    (t)=  𝒩(ω_0) { 1-e^-  κ(ω_0) t/2/ϵ^2
                ×[16 g^2-κ(ω_0) ϵsin(ϵ t/2)-κ(ω_0)^2 cos(ϵ t/2)]} ,
          ⟨ b^† b ⟩|_( loc)(t)= 
          𝒩(ω_0) {1-e^-  κ(ω_0) t/2/ϵ^2
                ×[16 g^2+κ(ω_0) ϵsin(ϵ t/2)-κ(ω_0)^2 cos(ϵ t/2)]} ,
           Im⟨  a b^†⟩|_( loc)(t)=
    4 𝒩(ω_0) κ(ω_0) g
    e^-κ(ω_0) t/2/ϵ^2 
    [1-cos(ϵ t / 2)] ,
 

    Re⟨  a b^†⟩|_( loc)(t)=0 ,





i.e.


    ⟨ H_ S,g⟩|_( loc)(t)=0  .


In Fig. <ref>
we plot the moments in the eigenmodes basis, comparing the results obtained by the Global, Local, convex mixture, Redfield, CP-Redfield approximations with the ones predicted by the exact dynamics, by including this time also the Lamb-shift contributions. In the local case for instance the Lambshift implies a tiny splitting between ⟨γ_+^†γ_+⟩|_( loc)(t) and ⟨γ_-^†γ_-⟩|_( loc)(t) at short time scales (connected to a small but non-vanishing Re⟨ a b^†⟩|_( loc)(t),  see Eq. (<ref>)).
Again the convex mixture of the local and global approximations of Eq. (<ref>) and the CP-Redfield equation yield a very good approximation either of the transient than of the steady state properties.












   





























10

riedel2017european
M. F. Riedel, D. Binosi, R. Thew, T. Calarco, https://doi.org/10.1088
  2, 030501 (2017).

nielsen-chuang-book
M. A. Nielsen, I. L. Chuang, Quantum Computation and Quantum Information:
  10th Anniversary Edition  (Cambridge University Press, New York, NY, USA,
  2011), 10th edn.

arute2019supremacy
F. Arute, et al. , https://doi.org/10.1038/s41586-019-1666-5Nature  574, 505 (2019).

benatti2003environment
F. Benatti, R. Floreanini, M. Piani, https://doi.org/10.1103/PhysRevLett.91.070402Phys. Rev. Lett.  91, 070402
  (2003).

lindblad1976generators
G. Lindblad, https://doi.org/10.1007/BF01608499Commun. Math. Phys.  48, 119 (1976).

gorini1976completely
V. Gorini, A. Kossakowski, E. C. G. Sudarshan, https://doi.org/10.1063/1.522979J. Math. Phys.  17,
  821 (1976).

redfield1957theory
A. G. Redfield, https://doi.org/10.1147/rd.11.0019IBM J. Res. Dev.  1, 19 (1957).

breuer2002theory
H.-P. Breuer, F. Petruccione, The theory of open quantum systems 
  (Oxford University Press on Demand, 2002).

jeske2013derivation
J. Jeske, J. H. Cole, https://doi.org/10.1103/PhysRevA.87.052138Phys. Rev. A  87, 052138 (2013).

lim2017signatures
J. Lim, et al. , https://doi.org/10.1063/1.4973975J. Chem. Phys.  146, 024109 (2017).

Purkayastha2016
A. Purkayastha, A. Dhar, M. Kulkarni, https://link.aps.org/doi/10.1103/PhysRevA.93.062114Phys. Rev. A  93, 062114 (2016).

gaspard1999slippage
P. Gaspard, M. Nagaoka, https://doi.org/10.1063/1.479867J. Chem. Phys.  111, 5668 (1999).

argentieri2014violations
G. Argentieri, F. Benatti, R. Floreanini, M. Pezzutto, https://doi.org/10.1209
  Lett.  107, 50007 (2014).

ishizaki2009adequacy
A. Ishizaki, G. R. Fleming, https://doi.org/10.1063/1.3155214J. Chem. Phys.  130, 234110 (2009).

benatti2003nonpositive
F. Benatti, R. Floreanini, M. Piani, https://doi.org/10.1103/PhysRevA.67.042110Phys. Rev. A  67, 042110
  (2003).

wilkie2001dissipation
J. Wilkie, https://doi.org/10.1063/1.1365955J. Chem. Phys.  114, 7736 (2001).

suarez1992memory
A. Suárez, R. Silbey, I. Oppenheim, https://doi.org/10.1063/1.463831J. Chem. Phys.  97, 5101
  (1992).

dumcke1979proper
R. Dümcke, H. Spohn, https://doi.org/10.1007/BF01325208Z. Phys. B Con. Mat.  34, 419 (1979).

benatti2005open
F. Benatti, R. Floreanini, https://doi.org/10.1142/S0217979205032097Int. J. Mod. Phys. B  19, 3063 (2005).

schaller2008preservation
G. Schaller, T. Brandes, https://doi.org/10.1103/PhysRevA.78.022106Phys. Rev. A  78, 022106 (2008).

cresser2017coarse
J. D. Cresser, C. Facer, https://arxiv.org/abs/1710.09939arXiv preprint arXiv:1710.09939   (2017).

seah2018refrigeration
S. Seah, S. Nimmrichter, V. Scarani, https://doi.org/10.1103/PhysRevE.98.012131Phys. Rev. E  98, 012131
  (2018).

jeske2015bloch
J. Jeske, D. J. Ing, M. B. Plenio, S. F. Huelga, J. H. Cole, https://doi.org/10.1063/1.4907370J. Chem.
  Phys.  142, 064104 (2015).

rivas2017refined
Á. Rivas, https://doi.org/10.1103/PhysRevA.95.042104Phys. Rev. A  95, 042104 (2017).

farina2019psa
D. Farina, V. Giovannetti, https://link.aps.org/doi/10.1103/PhysRevA.100.012107Phys. Rev. A  100, 012107 (2019).

cattaneo2019psa
M. Cattaneo, G. L. Giorgi, S. Maniscalco, R. Zambrini, https://doi.org/10.1088
  21, 113045 (2019).

levy2014local
A. Levy, R. Kosloff, https://doi.org/10.1209

hofer2017markovian
P. P. Hofer, et al. , https://doi.org/10.1088

adesso2017loc-vs-glob
J. O. González, et al. , https://doi.org/10.1142/S1230161217400108Open Syst. Inf. Dyn.  24,
  1740010 (2017).

rivas2010markovian
A. Rivas, A. D. K. Plato, S. F. Huelga, M. B. Plenio, https://doi.org/10.1088
  12, 113032 (2010).

de2018reconciliation
G. De Chiara, et al. , https://doi.org/10.1088

trushechkin2016perturbative
A. S. Trushechkin, I. V. Volovich, https://doi.org/10.1209
  (2016).

serafini2017quantum
A. Serafini, Quantum Continuous Variables: A Primer of Theoretical
  Methods  (CRC Press, 2017).

dicke1954coherence
R. H. Dicke, https://doi.org/10.1103/PhysRev.93.99Phys. Rev.  93, 99 (1954).

emary2003chaos
C. Emary, T. Brandes, https://doi.org/10.1103/PhysRevE.67.066203Phys. Rev. E  67, 066203 (2003).

hartmann2008quantum
M. J. Hartmann, F. G. S. L. Brandao, M. B. Plenio, https://doi.org/10.1002/lpor.200810046Laser Photonics Rev.  2, 527 (2008).

farina2019charger
D. Farina, G. M. Andolina, A. Mari, M. Polini, V. Giovannetti, https://doi.org/10.1103/PhysRevB.99.035421Phys. Rev.
  B  99, 035421 (2019).

hovhannisyan2020charging
K. V. Hovhannisyan, F. Barra, A. Imparato, https://arxiv.org/abs/2001.07696arXiv preprint
  arXiv:2001.07696   (2020).

OHM
A. J. Leggett, et al. , https://link.aps.org/doi/10.1103/RevModPhys.59.1Rev. Mod. Phys.  59, 1 (1987).

cheng2005markovian
Y. C. Cheng, R. J. Silbey, https://doi.org/10.1021/jp051303oJ. Phys. Chem. B  109, 21399 (2005).

hartmann2020accuracy
R. Hartmann, W. T. Strunz, https://doi.org/10.1103/PhysRevA.101.012103Phys. Rev. A  101, 012103 (2020).

CONVEX
M. M. Wolf, J. Eisert, T. S. Cubitt, J. I. Cirac, https://link.aps.org/doi/10.1103/PhysRevLett.101.150402Phys. Rev. Lett.  101, 150402 (2008).

cattaneo2019simmetry
M. Cattaneo, G. L. Giorgi, S. Maniscalco, R. Zambrini, https://arxiv.org/abs/1911.01836arXiv preprint
  arXiv:1911.01836   (2019).

PhysRevA.86.022340
P. Marian, T. A. Marian, https://link.aps.org/doi/10.1103/PhysRevA.86.022340Phys. Rev. A  86, 022340 (2012).

  













