Local versus global master equation with common and separate baths: superiority of the global approach in partial secular approximation

Open systems of coupled qubits are ubiquitous in quantum physics. Finding a suitable master equation to describe their dynamics is therefore a crucial task that must be addressed with utmost attention. In the recent past, many efforts have been made toward the possibility of employing local master equations, which compute the interaction with the environment neglecting the direct coupling between the qubits, and for this reason may be easier to solve. Here, we provide a detailed derivation of the Markovian master equation for two coupled qubits interacting with common and separate baths, considering pure dephasing as well as dissipation. Then, we explore the differences between the local and global master equation, showing that they intrinsically depend on the way we apply the secular approximation. Our results prove that the global approach with partial secular approximation always provides the most accurate choice for the master equation when Born–Markov approximations hold, even for small inter-system coupling constants. Using different master equations we compute the stationary heat current between two separate baths, the entanglement dynamics generated by a common bath, and the emergence of spontaneous synchronization, showing the importance of the accurate choice of approach.


Introduction
Open quantum systems of two coupled qubits are of fundamental importance in many disparate fields, being for instance at the basis of the realization of multi-qubit gates for quantum computation [1][2][3], distributed quantum sensing and metrology [4,5], and entanglement generation [6][7][8]. Such systems have been experimentally simulated in a variety of platforms, including trapped ions [9,10], superconducting qubits [11], or cavity QED arrays [12]. They are also useful in the context of quantum thermodynamics as they possess the minimum ingredients to realize thermal machines [13][14][15]. Furthermore, in spite of their simplicity, they allow for the observation of fundamental effects such as Dicke superradiance [16] or spontaneous quantum synchronization [17]. The derivation of the master equation describing the evolution of the qubits, and the subsequent search for an easy path to solve it, is therefore of the greatest importance.
While partial results investigating specific cases are available in the literature [18][19][20][21][22][23][24], a general description of the problem is still missing. In this paper we provide a comprehensive analysis based on a miscroscopic derivation in the case of two qubits, addressing: the presence of dissipative as well as dephasing baths, which can be common and/or separate, and considering a sufficiently general interaction between the qubits not limited to a Hamiltonian in rotating wave approximation (RWA), and also allowing for frequency detuning. As this is often

Deriving the master equation
The aim of the present work is to address a general Markovian master equation for two qubits that can be detuned, exchange energy and are coupled to thermal baths: we consider both dephasing and dissipative interactions, and both separate and common baths, as in the pictorial representation in figure 1(a).

Full Hamiltonian
Let us start by writing the free Hamiltonian of the system, in which we have set  = 1: where w ¢ 1 and w ¢ 2 are the frequencies of respectively the first and second qubit, and l¢ is the qubit-qubit coupling constant. We note that for the sake of generality we do not approximate the interaction by σ 1 + σ 2 -+h.c., as in RWA. The generality of equation ( 1 and l l w = ¢ ¢ 1 . Through all the work, we assume that the renormalized qubit frequencies are of the same order, ω 2 =O(1).
We now write the most general microscopic Hamiltonian of two coupled qubits interacting with common and separate thermal baths (consistently renormalized by the frequency of the first qubit): where following the convention of quantum optics the summation over k in the limit of infinite size bath represents as usual an integral over all the dense frequencies, and α=l 1 , l 2 , c indicates the specific bath. and the dissipative and dephasing couplings are mediated by the coefficients g x and g z . For instance, ( ) g x l 1 is the dimensionless coupling constant describing the strength of the dissipative interaction between the first qubit and the respective local bath, and so on. For simplicity, we take the coupling constants real. Notice that we are using the standard denomination for 'dissipation' and 'dephasing', where the former refers to a coupling through σ x , inducing both loss of energy and decoherence, while the latter denotes a coupling through σ z , causing, at least in the uncoupled case, pure decoherence but no energy leak. Markovian master equations can be derived in the weak coupling limit of the qubit-bath interaction. Therefore, we introduce a constant μ such that each coupling strength appearing in the interaction Hamiltonian is at maximum of the order of μ, i.e. m = ( ) ( ) g O j a " = a l l c c , , , 1 2 1 2 and j=x, z, and we assume μ=1 (consistently in units of ω 1 ). The coupling coefficients f k,α define the spectral density J α (ω) of each bath through: and we notice that the distinct dephasing and dissipative (and 'small' O(μ)) coupling ( ) g j a are not included in equation (7). One may wonder why, aiming at a complete description of any possible two-qubit system, we have considered the same bath inducing both dissipation and dephasing (in fact we could consider 6 instead of 3 baths). Assuming that different effects are due to different phenomena, a description employing a distinct bath for each of them should be necessary. Moreover, many uncorrelated environments may interact locally on each qubit, as it happens for example with a transmon qubit [30], so why shall we describe them through a single bath, as done in equation (3)? We anticipate that this assumption simplifies the notation and actually does not limit the following analysis, as we will be discussing in section 2.3.

Bloch-Redfield master equation in the secular approximation
In this section we will illustrate how to obtain a Markovian master equation starting from the microscopic Hamiltonian, stressing the validity of each approximation in order to get to a global Bloch-Redfield master  (3), where the environments are characterized by the inverse temperatures β ( c) , b ( ) l1 and b ( ) l2 , respectively for the common bath, the local bath on the first qubit and the local bath on the second qubit. An additional direct coupling between the qubits is mediated by the coupling constant λ. (b): Diagram of the states of the system Hamiltonian equation (2), setting λ=0 and with all the possible emission frequencies.
equation in the (partial) secular approximation. The possibility for a local master equation will be discussed in section 3. Let where ρ(t) and H I (t) denote the overall density matrix and the interaction Hamiltonian in the interaction picture representation (see [31] for details). By integrating equation (8), inserting it once again in equation (8) and taking the partial trace as usual, we obtain an integro-differential equation for the reduced density matrix of the system , if the environment is in a thermal state. We now set the validity of essential approximations in order to get to a Markovian master equation: Born approximation-The interaction between system and environment is so weak that the state of the latter is almost not perturbed by the coupling with the system. If the initial state of the overall system is the product state ρ(0)=ρ S (0)⊗ρ B , the evolved state at a certain time t is assumed product as well: The approximation (10) can be considered as a heuristic and intuitive way to obtain an important result, mathematically proven through the method developed by Nakajima [32] and Zwanzig [33]. Indeed, it can be shown [34] that, by inserting equation (10) in equation (9), we are neglecting terms of the order of O(μ 3 ), where μ is the coupling constant defined in the previous section. Therefore We point out that, while the full state ρ(t) is not expected to remain factorized (as in equation (10)) for long times [34], equation (11) is an exact result holding for any time t.
Let us now decompose the interaction Hamiltonian in the interaction picture in the following way: where A β (t) are system operators, while B β (t) are bath operators 4 . If we make the change of variable t = -¢ t t and insert equation (12) in (11), after some algebra we obtain: Tr 0 B B , with the assumption that the bath is stationary, i.e. [ρ B , H B ]=0. We are now ready to perform the next fundamental approximation.
Markov approximation-We assume that the bath operators have a very short correlation time, and the correlation functions decay as tb Remembering the weak coupling limit we then set τ B =τ R , i.e. the system will relax slowly with respect to the bath correlation functions, being τ R the timescale over which the state in the interaction picture changes appreciably. Considering the highest order appearing in equation (11), it is usually heuristically set where we remind that μ is the qubit-bath coupling constant renormalized by the frequency of the first qubit. The validity of the assumption needs often to be checked. For instance, in the limit of very high temperatures this might not be fulfilled, since a huge number of excitations would be available to interact with the system, making the decay rates very high as well. Nonetheless, there may exist the case in which, in the limit for the temperature  ¥ T , the autocorrelation functions of the bath decay faster than the relaxation time τ R , and therefore the Markov approximation is still valid. In this scenario, for  ¥ T the autocorrelation functions of the bath are proportional to a Dirac delta,  t d t µ bb¢ ( ) ( ), and we recover the so-called singular-coupling limit [31]. If now we calculate the integral in equation (13) for a sufficiently large time t * ?τ B , such that t * is still way smaller than the time τ R at which the state of the system in interaction picture changes appreciably, then we can safely replace ρ S (t−τ) with ρ S (t) in the same equation, since the dynamics of ρ S (t) is way slower than the decay 4 The bath operators B β in equation (12) should not be confused with B ( α) defined in equations (4) and (6) of  t bb¢ ( ). For the same reason, we can extent the integral till infinity, since the added part will give a negligible contribution. This is the Markov approximation, which sets a resolution on the timescale of the dynamics for t * , such that This corresponds to defining a certain coarse-grained timescale of the evolution; indeed, the Markovian master equation can alternatively be derived by making averages on these coarse-grained time intervals, as recently discussed in [35,36]. Finally we write: Unfortunately, to the best of our knowledge a precise order for the remainder neglected in the Markov approximation equation (16) has not been reported in general. An interesting bound is however provided in a recent paper [37], where instead of equation (15) the authors consider the tighter condition t m - B 1 . In general, assuming the condition in equation (15), the approximated master equation neglects a remainder of order higher than O(μ 2 ) (that from now on, we will drop); a special care in checking the validity of the Markov approximation in each specific case is anyway indispensable.
We will now further decompose the interaction Hamiltonian equation (12) by introducing the jump operators associated to each system operator A β : where   ñ {| } is the basis of the eigenvectors of the system Hamiltonian H S . The following properties hold: By writing equation (16) with the time-evolved jump operators, we get to the Bloch-Redfield master equation where we have introduced the one-side Fourier transform of the bath correlation functions . 20 t 0 i Secular approximation-The evolution of the state of the system ρ S (t) has, in the interaction picture, a typical relaxation timescale of the order of the square of the inverse of the coupling strength μ, as stated in equation (14). If there exist values of w¢ and ω in equation (19) being coarse-grained in time as from equation (15), i.e.
then the terms in equation (19) oscillating with frequency w w ¢will not give any significant contribution to the system evolution, since by integrating equation (19) for a time t * such that the fastoscillating quantities vanish. Equation (21) corresponds to a refinement of the coarse-grain condition written in equation (15). Indeed, a slightly different approach to the derivation of the master equation makes use of a unique coarse-grained average, including both the Markov and the secular approximation (see for instance [35,36]). Notice that the interaction picture is particularly suited to distinguish the terms bringing a negligible contribution to the evolution of the system.
Neglecting the fast oscillating terms in the interaction picture is usually referred to as secular approximation. Unfortunately, it is easy to run into a nomenclature issue in the literature: in many works we can find the name 'secular approximation' for the removal of all the terms in equation (19) for which w w ¢ ¹ , without questioning the validity of equation (21). This is of course feasible for all the systems in which the relevant frequencies are well-spaced, i.e. w w t m ¢ -» - | | R 1 2 for any w w ¢, , but it might lead to confusion in other cases, as we will discuss in section 3.
To avoid confusion, we will call full secular the approximation for which we neglect all the terms in equation (19) with w w ¢ ¹ , while we will employ the name partial secular approximation for the cases in which we keep some slowly rotating terms with w w ¢ ¹ , for which the relation in equation (21) would actually fail. For the sake of clarity, when discussing concrete examples throughout the paper we will name the master equation with cross terms often retained in the secular approximation (21) as master equation in partial secular approximation (also in regimes where these terms could be neglected).
After some algebra, the Bloch-Redfield master equation equation (19) may be rewritten in the Schrödinger picture as where we have introduced the Lamb-Shift Hamiltonian: and the dissipator of the master equation, responsible for the energy losses of the system: , .
Notice that, prior to the secular approximation, the Lamb-Shift Hamiltonian is not Hermitian and contains imaginary terms as well, and we do not have a 'purely dissipative' dissipator. By employing the full secular approximation and coming back to the Schrödinger picture, the Lamb-Shift Hamiltonian and dissipator read: The master equation (22) with Lamb-Shift Hamiltonian and dissipator given by equation (26) is written in the GKLS form [38][39][40], and it therefore generates a dynamical semigroup, i.e. a perfectly Markovian evolution. On the other hand, this is a strong condition which is not necessary to get a GKLS form of the master equation. In fact, some very recent papers [35,36] have shown that performing an accurate partial secular approximation leads to a GKLS master equation as well (see also [41,42]). An interesting observation is that the partial secular approximation condition (21), is sufficient to remove fast oscillating terms leading to a dissipator (24) where , with ω and ω′ with the same sign (which would produce a squeezing-like effect), are prevented. In fact the fastest terms, more susceptible to fulfill the condition (21) will oscillate at frequency w w --¢ | ( )| (again with ω and ω′ with the same sign): if this is the case, then all terms will be consistently neglected. Finally, let us term the master equation (19) in partial secular approximation, which may be rewritten in the form of equation (22), 'global master equation with partial secular approximation'. We stress the fact that this master equation is derived within the Born-Markov approximations. The latter can be more delicate to assess and unphysical effects can arise as signatures of an inaccurate Markovian description of an intrinsically non-Markovian evolution [43]. On the other hand, in general it is immediate to establish the validity of the condition (21) to get an equation in the partial secular approximation. Also the inaccurate secular approximation, i.e. out of the validity region (21), can lead to unphysical effects, as can be displayed by a full secular master equation with respect to a partial secular one.

Diagonalizing the system Hamiltonian and finding the jump operators
As discussed in the previous section, diagonalizing the system Hamiltonian H S is a necessary step to derive the Markovian master equation, since it allows us to write in the correct form the jump operators defined in equation (17) which describe the effects of the interaction with the baths.
In the interaction Hamiltonian we can find the system operators σ j x and s j z coupled to the bath operators, with j=1, 2. Their decomposition in terms of jump operators is readily written according to equation (17): (28) describes the possible emission and absorption processes of the system, depicted in figure 1 A self-consistent secular approximation depends on the detuning between the qubits. In the case in which there is a small detuning, such that ω 1 −ω 2 is not way greater than μ 2 , we cannot employ the full secular approximation, but we need to rely on a partial secular approximation in which we keep in equation The validity of the full secular approximation for big detuning and its breakdown in the opposite scenario are respectively shown in figures 2(a) and (b). The most evident difference between partial and full secular descriptions in the regime in which the latter fails (figure 2(b)) is the presence of so-called quantum beats, i.e. of oscillations of the population of the excited state of the qubit. The quantum beats are a well-known phenomenon occurring during a superradiant emission [44] and predicted also for two nonidentical atoms [45] as in our case. The master equation with partial secular approximation correctly describes them, while the full secular one is too crude and not able to reproduce the beats, leading to a completely smooth evolution. Further cases in which the full secular approximation is not suitable will be discussed in section 4.
Starting from equation (19) and employing the notation of equations (23) and (24)  By looking at the jump operators in equation (28) and at the partial secular approximation performed on equation (19) we can now address the claim in section 2.1 about the simplifying choice of considering one single bath inducing both dephasing and dissipation. The point is that considering two distinct baths rather than a single one is in general not needed (unless dephasing and dissipation need to be considered with different spectral density or baths temperatures), and would lengthen all expressions. More in detail, considering multiple equivalent and independent baths could at most affect the coefficients w G bb¢ ( ) (see appendix B for their specific form): then we would have w G = bb¢ ( ) 0 for any system operators A β and b¢ A coupled to distinct baths through  figure 2(a), the full secular approximation (dashed red) provides a correct way to describe the evolution, although the tiny oscillations given by the partial secular approximation in which we keep the cross terms (solid blue) can be observed by zooming to a proper small time interval (inset), which anyway cannot be resolved in the timescale defined by the coarse-graining. In figure 2(b), due to the small detuning, the full secular approximation (dashed red) fails and it leads to a completely different evolution with respect to the partial secular (solid blue).
B β and b¢ B , since there are no correlations between the two baths. On the contrary, if B β and b¢ B are operators of the same bath, then w G bb¢ ( ) does not vanish a priori and there could be a case in which Let us for instance consider the coupling with the local bath on the first qubits: What if the dissipation would be induced by a bath different than the dephasing one? Looking at equation (28) we can see that the operators s x 1 and s z 1 may in theory couple in equation , but their corresponding terms would vanish because of the partial secular approximation, since w It is easy to recognize that this argument holds for any case where dephasing and dissipation could arise from different baths. Therefore, for equivalent but independent baths, the simplified Hamiltonian equation (3) can be assumed. Otherwise, considering 6 baths (instead of 3) would lead to different values of the bath correlation functions in equation (B.1), but not change the structure of the master equation 5 . We will reach the same conclusion in presence of qubits coupling, apart from a singular case (corresponding to a very specific parameter choice, when the condition in equation (C.4) holds).
Following the same path, we can readily see that if different sources, associated to different baths, would induce, let us say, independent dissipations on the same qubit, by assuming a single dissipative local bath we are not losing generality, since the effects of the multiple baths would not change the form of the master equation, but at most the value of the coefficients: the effects of independent baths would just sum, i.e. the final decay rate would be the sum of the decay rates given by each single bath. The argument we have just discussed is reflected in the values of the coefficients in equation (B.1).

Direct coupling
The case in which we have a direct qubit-qubit coupling should not in principle be more complex, since all we need to do is to diagonalize a 4×4 matrix, find the corresponding eigenvalues and eigenvectors and work in the new basis. The system Hamiltonian H S now reads: and the corresponding matrix in the canonical basis ñ ñ ñ ñ {| | | | } 11 , 10 , 01 , 00 is written as We can easily diagonalize equation (32) by finding the eigenvalues: with associated eigenvectors [19] q where the parameters θ and f are given by Once we know the spectral decomposition of the Hamiltonian, we can proceed to calculate the jump operators associated with each system operator appearing in the interaction Hamiltonian, i.e. s j x and s j z with = j 1, 2. The explicit form of each jump operator is given in appendix C. With the aim at a complete description of the problem, we also consider the possibility that two of the eigenstates of the system are almost degenerate, which would happen if w - 1 and λ=1, as it can be seen from equation (33). In this case, some additional terms beyond the full secular approximation need to be consistently kept, as fully observed in appendix C.

Local versus global: an in-depth discussion
A debate about the validity of the local rather than the global description of an open quantum system has arisen since the early era of the field: to the best of our knowledge, the first discussions about how to derive a global master equation accounting for the inter-system interactions date back to the early seventies [46][47][48]. Twenty years later, Cresser observed the failure of the local approach to describe a lossy Jaynes-Cummings model [49], terming 'phenomenological master equation' what it is nowadays usually called 'local master equation'. Quite the same issue has been addressed in some more recent papers [50,51], extending the analysis to three coupled Josephson junctions [52] or coupled harmonic oscillators [34], while the 3-level atom has been investigated in [53]. Some comments about the validity of the local approach to describe energy transport in chains of harmonic oscillators or spins appear in [54][55][56].
In the past few years, a renewed interest in the topic has grown, also because of a paper in 2014 suggesting that the local approach was breaking the second law of thermodynamics in a thermal machine composed of two quantum nodes [25]; this violation was later shown to be beyond the order of the employed approximation, thus only apparent [26]. Related discussions date back to 2002 [57,58]. Furthermore, a recent paper has shown that the local master equation reconciles with the laws of thermodynamics when analyzing a suitable associated collisional model [59]. Connections between the local evolution of an open system and its thermodynamic microscopic model were previously addressed in [60,61]. The investigation of the local versus global problem in different scenarios is nowadays quite active [24, 27-29, 59, 62-71]. For instance the failure of the local approach when studying two coupled qubits is claimed in [24,63,66]. On the contrary, two distinct works have tested the validity of the local description applied to the calculation of thermodynamics quantities in quantum heat engines [27,28], showing its goodness in a quite large range of parameters of the coupling constant, and claiming that the global approach fails when the two subsystems are weakly coupled. More precisely, this assertion is due to a restrictive consideration of the global master equation as limited by a full secular approximation, which also the authors recognize as responsible for the breakdown of the master equation. In [27] the possibility for a partial secular approximation is also suggested in order to cure such deficiency, and many other papers have pointed out why the full secular approximation may not be valid in some parameters ranges of different scenarios [25,34,56,64,69,70]. In the following, we therefore analyze in detail both local and global approach and show that the partial secular approximation allows to derive a global master equation that never leads to unphysical results, given that in the limit l  0 it coincides with the local master equation. The discussion in sections 3.1 and 3.2 are generally valid also beyond the 2-qubit system, while section 3.3 addresses the validity of the local approach and full secular approximation in the specific case of two coupled spins.
3.1. Setting the nomenclature Let us start by setting a common nomenclature for local and global approach. We discuss the case of two subsystems, but generalizations to multipartite systems are straightforward. We can thus consider H S =H 1 +H 2 +H 12 with no need of specifying the nature of the subsystems and their interaction (for two spins we are considering H j =ω j /2σ j z , and ls s = H x x 12 1 2 ). We now recall the local and global master equations for an open quantum systems. Local master equation.The local approach is an approximation that consists in calculating the jump operators in equation (17)  2 , i.e. neglecting the interaction between the subsystems when computing the effects of the environment. This clearly leads to two separate sets of local jump operators which (non-trivially) act only on the first or second subsystem. If a full secular approximation is applied, the direct coupling between the subsystems only appears in the commutator [H S , ρ S ] of the Bloch-Redfield master equation (22), thus it only influences the unitary part of the evolution. Intuitively, the local approach is expected to provide us with a valid approximated master equation only when the coupling constant between the subsystems is sufficiently small [26-28, 56, 59, 72].
Global master equation.The global approach consists in considering the full (exact) system Hamiltonian, interacting term included, when calculating the jump operators. Hence, the global master equation is the Bloch-Redfield one (19) without further approximations. The jump operators appearing on the right term of the equation are not local anymore, i.e. since they are obtained after the diagonalization of H S , they can act on both the first and the second subsystem. The global master equation is in general more precise than the local one, as the latter relies on a further approximation. It might however be too involved to be solved, due to the nonlocality of the jump operators [27]. To simplify its form, one may rely on the standard secular approximation, provided that the condition in equation (21) is fulfilled.
In the recent literature it is sometimes used the term 'global master equation' to indicate the result of the global approach described above and after having performed an indiscriminate full secular approximation. We believe that this nomenclature may lead to confusion: per se, the fact of being 'global', i.e. to lead to jump operators which act jointly on both the subsystems, is not related to the secular approximation. The reported appearance of unphysical currents is not due to the non-locality of the jump operators (global approach), being instead the result of the indiscriminate application of the full secular approximation when only the partial one was justified [27,28,56]. The attempt to apply the full secular approximation for any different frequencies w w ¢ ¹ , even when the condition in equation (21) is not fulfilled, leads to inconsistencies: all the approximations listed in section 2.2 are indeed valid in well-defined parameter regimes. From the formal derivation in the previous section, in the framework of Born-Markov approximations, a global master equation with a justified partial secular approximation is in general more accurate (or less approximated) than any local one. In reference [27] this was also suggested and named partial Markovian Redfield master equation. An important exception is a master equation derived in the singular-coupling limit [28,31], that would lead to a local master equation. In this particular case, the global master equation with partial secular approximation, even if accurate, would be unnecessarily more complicate than the local one. For instance, this is the case when addressing a Markovian scenario with very high temperature in which the autocorrelation functions of the bath decay faster than the system itself [29].
One may argue that the full secular approximation is anyway preferable to the partial one, since it is generally introduced to obtain a GKLS master equation [28] such as the one in equation (26), free from any unphysical behavior. It is indeed known that the Bloch-Redfield equation (19) may in some cases violate the positivity of the dynamical map [73]. However, if the full secular approximation is not well justified from a microscopic model (because equation (21) does not apply) then a global full secular master equation needs to be considered as a phenomenological one, as the correspondence with the microscopic model is lost. On the other hand, as the approximations we have performed to obtain equation (19) are correct up to the order of the remainders, unphysical departures are expected to be consistently small [43]. For an in-depth discussion about Bloch-Redfield equation and complete positivity we refer the reader to the broad literature on the topic [64,[73][74][75][76][77][78][79][80][81].
For our purpose, we conclude stressing that the GKLS form of the master equation is also guaranteed by the correct application of the partial secular approximation, as recently discussed in [35,36]. This is therefore a preferable approach, being well related to a microscopic model instead of being phenomenological.

Accuracy of the local master equation
In order to assess the accuracy of the local master equation, let us write the interaction Hamiltonian as H 12 =λ V, where λ is a 'small' parameter which we consider as a perturbation order, i.e. λ=1 . Following [26], we apply standard perturbation theory to find the zeroth order eigenvectors and eigenvalues 6 . Within the assumption of not-degenerate system Hamiltonian, we write the eigenvalues as the infinite perturbation expansion [82]: . . where ω (0) is the frequency given by differences of the unperturbed eigenvalues - If we employ the local master equation to compute physical quantities, clearly we can resolve them only up to the order of O(μ 2 ). Any quantity of the order of O(μ 2 λ) or smaller, then is null in the framework of the local approach. This is the reason why the violation of the second law of thermodynamics [25] is only an apparent one [26], given that it is of the order of O(μ 2 λ 2 ).

Local versus global approach for two coupled qubits
We will now address the local versus global comparison focusing on the case of two spins as in equation (2) and interaction ls s = H x x 12 1 2 . As discussed in the previous sections, the global master equation with partial secular approximation is always valid up to the errors induced by the Born-Markov approximations. On the contrary, the local master equation and the full secular approximation are accurate only for some parameter regimes, which we will investigate starting from the derivation of the master equation. We recall the fact that the local approach is always valid in the singular-coupling limit, which is more restrictive.
We will first present the local master equation for two coupled qubits, studying the ranges of parameters in which each approximation works, in the presence of common or separate baths, and then show how the difference between master equations emerges through some illustrative examples.

Local master equation
The local approach is valid in the case in which λ=1. In order to derive the master equation we must find the jump operators, thus the eigenvalues and eigenvectors of the Hamiltonian. Let us start from the situation where the local Hamiltonian is not degenerate, i.e. λ=ω − . In this case we have the following eigenvalues: Equations (41) and (42) are the zeroth order expressions for the infinite perturbative expansion in equation (38). By inserting them in equation (40), we see that the local master equation is exactly the master equation (29) found for decoupled qubits, but with free system Hamiltonian H S including a coupling term, i.e. H S given by equation (31) instead of equation (27). So, the difference between the local and the global master equation is of the order of O(μ 2 λ).
In the degenerate regime ω − =0, there is an apparent freedom in the choice of the basis with respect to which the perturbative expansion must be performed, as any linear combination of ñ |e 1 and ñ |e 2 could in principle be selected. This apparent freedom is actually removed by the interaction Hamiltonian. For instance, in the case of separate baths, deriving the master equation starting from any possible choice of the basis would in any case lead to local jump operators, as in the absence of degeneracy. Working near degeneracy, that is, , but the aforementioned selection rule would still apply and the final master equation would not change.

Comparison between master equations
In this section we will investigate the limits of validity of the full secular approximation and of the local master equation depending on the parameters of the system and show some examples of the different system dynamics generated by them. In section 4 we will provide further examples focused on physical quantities computed through distinct master equations. The results are summarized in table 1 presented in the concluding remarks. The parameters we vary in order to study each master equation are the qubit-qubit coupling constant λ and the detuning between the qubits ω − , as can be seen in table 1. All the remaining parameters will be fixed as follows, paying attention to the conditions for the Born-Markov approximations.
• We choose as qubit-bath coupling constant μ=10 −2 . This means that the timescale of the evolution of the system will be τ R =O(μ −2 )=10 4 . This quantity is important to check the validity of each approximation, as shown in table 1. The remainder given by the Born approximation will be, according to equation (11), of the order of O(μ 3 )=10 −6 .
• Both the common and the separate baths will have an Ohmic spectral density, i.e.
where J(ω) is defined in equation (7) and Ω is a cutoff frequency which we have set as large as Ω=20.
Regarding the inverse of the temperature of each bath, we have chosen . An unbalance between the local baths is important in quantum thermodynamics, in order to study the heat transport between them. We finally have to check that these baths with the chosen temperatures satisfy the condition for the Markov approximation equation (15). The timescale of the decay of the bath autocorrelation functions for an Ohmic spectral density reads [31] , and τ R =O(μ 4 ), therefore τ B =τ R and the Markov approximation is valid for our choice of parameters.
• We choose as initial state of the system the factorized state r r r = Ä OV OV 0 , where ρ OV is the completely overlapped state r = ñá + ñá + ñá + ñá We will evaluate the dynamics using four different master equations, namely the global master equation in partial secular approximation (GP), which we will consider as the most correct one according to the discussion in section 2.2, the global master equation in full secular approximation (GF), the local master equation in partial secular approximation (LP), and the local master equation in full secular approximation (LF). We remind that here we are using 'partial secular approximation' to refer to a master equation in which we keep the cross terms which, in some scenario, may be slow-rotating and not negligible, even in regimes in which such terms may actually be eliminated. Each master equation leads to a (a priori) different evolution. For simplicity, we focus on three different figures of merit. The first one is the mean value of s z 1 as a function on time, i.e. s á ñ ( ) t z 1 . The second is the fidelity [83] between the state obtained through the global master equation with partial secular (most accurate one) and a state computed with another master equation, i.e.  r ( ( ) ·) t , GP , as a function on time. The third and last figure of merit is the steady state of the system, i.e. the state obtained for  ¥ t . Notice that, while the fidelity is quite a general and reliable indicator, both the population of the first qubit and the steady state may not display differences between two master equations, even if the latter are substantially different.
For convenience, we first study the scenario with separate baths only, and then the one in the presence of a common bath, addressing in both cases the local and global master equation separately, and focusing on dissipative couplings with the environments, as we do not expect any qualitative difference if we also add dephasing baths.

Separate baths
• In the case of separate baths, the local master equation presents local dissipators on each qubit, and the interaction between the subsystems only comes into play in the unitary part of the evolution. This means that the full secular approximation always coincides with the partial one and is valid for any value of the detuning ω − , provided that λ=1, which fixes the validity of the local master equation.
• The full secular approximation with global approach may break down for some range of values. Indeed, the global approach makes use of the basis of eigenmodes to build the jump operators, which is composed of entangled states, therefore a single local bath coupled to s x 1 induces dissipation on the second qubit as well, and the secular approximation comes into play. Let us look at the jump frequencies presented in equation (C.1). We have to identify the frequency differences which may be comparable with If the qubits have small detuning, i.e. w - | | 1, and the qubit-qubit coupling constant is small as well, λ=1, the secular approximation on these frequencies equation (21) does not apply if for instance Still, if the basis of eigenvectors of the system Hamiltonian is quasi-local, in the sense that it is well described by the states in equation (42), the cross terms between the qubits arising with separate baths are very small, of the order of O(μ 2 λ). Thus we can neglect them, and the full secular approximation is still valid even in the global case. This final condition of validity reads λ=ω − , highlighting the importance of the ratio between detuning and qubit-qubit coupling constant.
We show an example in figure 3, considering s á ñ ( ) t z 1 and the fidelity between evolved states as a function of time, and in figure 4, focusing on the steady state. Being the baths separate, the local master equation with full secular approximation coincides with the partial secular one. In figure 3, since λ is very small, the local approach provides a reliable description of the dynamics independently of the value of the detuning. On the contrary, for identical qubits (figures 3(a) and (b), with ω − =0), the global master equation with full secular approximation fails, while this approximation in the global approach is justified for λ=ω − (figures 3(c) and (d)), despite the detuning being small. Looking at the stationary states and also allowing for baths at different temperatures, we show in figure 4 the predicted parameters regimes of failure of the local master equation and of the full secular approximation in the global one. While λ=ω − , both the approaches are reliable, but as soon as λ gets close to ω − the GF fails; this clearly starts from smaller values of λ on the figure 4 left than on the right, since in the former case the detuning is smaller. As far as λ increases toward 1, the global approach with full secular recovers validity, since it fulfills the condition for the full secular approximation. As λ becomes of the order of the qubit frequency O(1), the local m.e. loses reliability.

Common bath
• If the bath is common, in general the local master equation does not display local dissipators. Indeed, if the detuning ω − is small, i.e. not way larger than t m = -( ) O R 1 2 , we obtain cross terms in the master equation which have the effect of exchanging excitations between the qubits (see equation (29)). Therefore, the full secular approximation for the local master equation is valid only if w m - 2 . This means that the claim about the goodness of an indiscriminate full secular approximation when following the local approach is not general, but limited to the separate baths scenario only. Of course, the local master equation is accurate only if λ=1.
• For the global master equation in the presence of a common bath, the same discussion about the case with separate baths hold, with the only difference that a local basis (λ?ω − ) does not allow us to perform the full secular approximation anymore. Therefore, the condition for the global master equation with full secular approximation reads ω − ?O(μ 2 ) or λ?O(μ 2 ).  first row), the local approach with partial secular approximation provides a reliable description of the dynamics. On the contrary, both the local and global approach with full secular approximation fail, since the detuning is very small as well. In the scenario of λ being of the order of the qubit frequency (second row), the local approach always fails, while the global approach with full secular approximation is reliable in spite of the small detuning, since λ?μ 2 .

Common bath: entanglement, quantum beats and synchronization
We will now show how, in the critical regime of small qubit-qubit coupling constant λ, small detuning and a common bath, the full secular approximation leads to unphysical results, while the partial secular correctly describes several physical phenomena, both in the local and in the global approach. Since we are considering the limit λ=1, the same considerations hold for the case of uncoupled qubits discussed in section 2.3.1. In particular, we consider the case addressed in figures 5(a) and (b) for the coupled case, i.e. ω − =10 −2 , λ=10 −4 and the rest of parameters as discussed in section 3.3.2; to include a discussion about the uncoupled case, we consider the scenario of figure 2(b) with ω − =10 −2 . These two situations are almost equivalent due to the very small λ, as can be seen by comparing figures 2(b) and 5(a).
We first focus on the dynamics of entanglement obtained through different master equations. It is indeed wellknown that a common bath may generate entanglement between non-interacting qubits immersed in it [84], even when the reservoir is Markovian [7]. This phenomenon has been predicted using diverse methods of obtaining a master equation, such as a coarse-graining procedure [85] or employing master equations originally derived for quantum optics [86]. As entanglement measure we choose the negativity  [84]; in the consider case of two qubits, a non-zero value of the negativity is a necessary and sufficient condition to have entanglement. We plot in figure 6 the negativity as a function of time for both uncoupled and coupled case, which do not show a visible difference consistently with the very small qubit-qubit coupling constant. Both in the local and global approach, the master equation with partial secular approximation correctly displays entanglement creation, sudden death and subsequent sudden birth [87], that resemble dynamics already observed in similar scenarios [86]. On the contrary, the full secular approximation completely misses the detection of entanglement, since it does not include a qubit-qubit coupling mediated by the common bath.
As already discussed in section 2.3.1, another physical phenomenon that the full secular approximation is not able to reproduce are the quantum beats, i.e. the oscillations in the dynamics of the qubit populations that can be observed in figures 5(a) and 2(b). The quantum beats are known since the studies on superradiance in the eighties [44], and appear during the evolution of two slightly-detuned atoms because of the tiny difference between their frequencies in the phase of the emission power [45]. The full secular approximation does not make the two qubits 'communicate', and therefore it does not detect the detuning and leads to an incorrect smooth decay of the population of the excited state. Finally, we mention another important physical effect missed by the full secular approximation. Quantum synchronization is a paradigmatic phenomenon investigated in disparate fields in the recent years (for a review see [88]). In particular, spontaneous synchronization of the spinning frequencies of two uncoupled qubits in a common bath has been predicted using a Bloch-Redfield master equation without any further secular approximation [17]. Using the master equation (29) in partial secular approximation we have observed quantum synchronization starting from a time t≈6000, while the same master equation in full secular approximation never displays synchronization of the qubit frequencies, and it is therefore not suitable to analyze such a phenomenon.

Steady state heat current incoming from separate baths
We now consider the case of two coupled qubits and separate baths addressed in figure 4: for small values of the detuning and varying the coupling constant λ, we compute the steady state heat currents coming from the hot and cold reservoir. As in section 3.3.2, we assume that the inverse temperatures of the baths are b = where H S is the system Hamiltonian equation (31) and  2 is the dissipator generated by the hotter bath coupled to the second qubit. Analogously we can define the heat current incoming from the colder reservoir as , and since r ¥ is the steady state we have J 1 +J 2 =0. Note that equation (44) is widely used in the literature [27,28,89], but cannot be associated to a heat current observable. A definition based on a current observable can be found in [90], which also discusses some issues regarding its consistency with different master equations. Figure 7 depicts the stationary heat current incoming from the hot bath as a function of the coupling constant λ, with detuning ω − =10 −3 (a) and ω − =10 −5 (b), for distinct master equations. We see that there is a region of λ in which the global master equation with full secular approximation fails, highly overestimating the value of the current: this is the effect that was extensively observed and discussed in the recent works on the topic [27,28]. The full secular approximation breaks down because ω − =1 and λ=1, making small the energy difference in equation (43). However, we also observe that such region depends on the ratio ω − /λ, as discussed in section 3.3.2 and displayed in table 1: if ω − ?λ, the eigenmodes basis of the system is almost local and the global approach with full secular is accurate as well. For this reason, the range in which it fails is more narrow for ω − =10 −3 than for ω − =10 −5 . If the detuning were null, the global approach with full secular approximation would provide a non-zero heat current even for l  0 [28], which is clearly non-physical. On the contrary, the global approach with partial secular approximation and the local approach (for which partial and full secular coincide) provide a correct description of the incoming heat current for small values of λ. Increasing the coupling constant, we observe that the current increases as well till reaching the value given by the full secular approximation, which starting from here recovers its validity. For big values of λ the global master equations describe a current decreasing toward 0. This is correct, since if λ?1 the only relevant part of the system Hamiltonian H S is ls s x x 1 2 , and therefore it commutes with the interaction Hamiltonian: Hence, the dissipator associated to each bath does not induce an energy exchange in any stationary state (which now depends on the initial conditions), and no heat current is produced. On the contrary, the local master equation is written in a basis which is different from the diagonal basis of H S , and thus . We confront the value obtained through the global master equation (m.e.) with partial secular approximation (solid blue), the global m.e. with full secular approximation (dashed red), the local m.e. with partial secular approximation (dotted-dashed orange), and the local m.e. with full secular approximation (dotted green). For both local and global approach, the full secular approximation incorrectly shows no entanglement during the evolution of the qubits, while both the global and local m.e. with partial secular approximation provide the phenomena of entanglement creation, sudden death and sudden birth. Since λ is very small, we observe no remarkable difference between the coupled and uncoupled case.
indicates a fictitious non-zero heat current even for λ?1. Note that the global approach with partial secular approximation is appropriate in all considered parameter regimes.
The reader can verify that the failure of the local master equation or of the full secular approximation to describe the incoming heat current in different regimes correctly reproduces the parameters ranges summarized in table 1.

Concluding remarks
In the present work we have extensively addressed the derivation of the Markovian master equation for two qubits interacting with thermal baths, in order to assess the validity of local and global master equations. A comprehensive description is achieved considering all the possible scenarios: presence of separate as well as a common bath, including both dissipative and dephasing interaction, and taking into account the possibility of a direct coupling between the qubits and of detuning. The Markovian master equation, in the form of a Bloch-Redfield master equation, has been derived in section 2 reviewing all the necessary approximations, carefully stating the condition for the validity of the Born and the Markov approximations. We have obtained two general master equations: equation (29) in the case without a direct qubit-qubit coupling, and equation (36) when a direct coupling between the qubits must be taken into account. . We confront the value obtained through the global master equation (m.e.) with partial secular approximation GP (solid blue), the global m.e. with full secular approximation GF (dashed red), and the local m.e. with full secular approximation LF (dotted-dashed orange), which coincides with the local m.e. with partial secular approximation. When λ=1 both the global m.e. with partial secular approximation and the local master equation correctly reproduce the value of the heat current, which decreases toward zero for l  0. On the contrary, the global m.e. with full secular approximation highly overestimates the value of the current when the detuning is not way higher than the coupling constant. Starting from λ=O(1), the local approach incorrectly produces a stationary non-zero value of the heat current, while the global approach accurately describes the decay of the current toward zero. Table 1. Conditions of validity of each master equation (m.e.), depending on the values of the detuning between the qubits ω − , of the qubit-qubit coupling constant λ and of the qubit-bath coupling constant μ. Each possible scenario is taken into account: local or global master equation, partial or full secular approximation, presence of a common as well as separate baths. We recall that all the constants are dimensionless quantities, according to the renormalization discussed in equations (1)- (3), and that the general condition for the validity of the secular approximation is presented in equation (21). The table only deals with scenarios in which the Born-Markov approximations hold. Furthermore, we recall that the local approach is always valid in the singular-coupling scenario, independently of the system parameters.
Validity of the master equation These preliminary steps put ourselves in the conditions of determining the validity of the final possible approximation, namely the secular one. We have established the requirements under which one can apply a full secular or only a partial secular approximation. This assessment is especially relevant in the context of the feasibility of a local approach to the master equation with respect to the global one, which may simplify the computation of the solution in many cases. This renewed problem is addressed in section 3.
Many works have already proven that the global master equation may fail in some scenarios because of the breaking of the (full) secular approximation. In this paper, we have shown how to overcome such problem by applying an accurate partial secular approximation which does not remove slowly-rotating terms. In particular, in section 2 we have provided an extensive mathematical derivation of such equation, termed global master equation with partial secular approximation, proving that it is always the most correct choice in any parameters scenario in which the Born and Markov approximations are valid. Then, in section 3 we have shown how to derive the local master equation, and focused on the comparison between it and the global approach for the case of two coupled qubits, with partial or full secular approximation. The local master equation is always accurate in the singularcoupling limit, as already extensively proven [28,31]. If this limit is not assumed, the local approach is valid only for small values of the qubit-qubit coupling constant λ=1. The feasibility of the full secular approximation must be checked: if a common bath is present, the detuning between the qubits plays a fundamental role, since a small value of it would make the local master equation with full secular approximation fail; in the case of the global master equation, the value of λ comes into play as well. If the qubits interact with separate baths only, some subtleties emerge: for the local master equation the full secular approximation is always valid, while in the global approach care must be taken, since if both λ and the detuning are very small, the condition for the approximation may break down. We have shown that the value of the qubit-qubit coupling constant λ is not the only important actor here: the ratio between λ and the detuning ω − must be considered as well, since if ω − ?λ the full secular approximation in the global approach recovers its validity. These results are summarized in table 1, where we have highlighted the scenarios with local or global approach, partial or full secular approximation, and common or separate baths. At the end of section 3, we have discussed the consequences of the local versus global issue in the case of two qubits, and we have compared the results of the system evolution obtained through different approaches by depicting them in figures 3-5.
In section 4 we have shown how several physical quantities change when been computed through different master equations. In particular, in the case of weak qubit-qubit coupling constant, small detuning and a common bath, the full secular approximation is not able to detect important phenomena such as quantum beats and quantum synchronization, and does not produce entanglement during the evolution as depicted in figure 6. In the scenario with two separate baths with unbalanced temperatures, the global master equation with full secular approximation leads to an unphysical stationary heat current in the region in which the detuning ω − is not way larger than λ. On the contrary, after a transient in which the stationary heat current is correctly detected, it is the local master equation which predicts a non-zero fictitious current when λ?1. The global master equation with partial secular approximation reproduces the correct physical results in any scenario and any range of parameters.
Our discussion remains valid while considering general interacting bipartite systems and can also be immediately extended to multipartitite scenarios, while more challenging will be to explore local versus global approaches in non-Markovian situations. Also, this analysis is relevant for extended systems experiencing dissipation only in some of their parts, so that global master equations need in principle to be considered [91]. Our general conclusions are especially timely, because of the renewed interest on the topic that has lead to several results, sometimes contradictory or only partial. Beyond being a fundamental instrument for the appropriate description of coupled qubits in contact with environments, further implications may be foreseen in the context of thermodynamics, computations and information, considering the differences arising between phenomenological approaches and the microscopically derived global master equation with partial secular approximation. The jump operators with negative frequencies are obtained by employing the property in equation A . Notice that, once again, the jump operators associated to s j x and s j z have different frequencies, whose difference is not 'small' in the sense of the condition for the secular approximation equation (21). Actually, there may be a singular case in which two frequencies of different bath operators assume the same value, namely ω II and ω IV : Figure C1. Diagram of the states of the system Hamiltonian equation (31), with all the possible emission frequencies.