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

Let us start by writing the free Hamiltonian of the system, in which we have set ${\hslash }=1$:

$$\begin{eqnarray}&&{H}_{S}^{{\prime} }=\displaystyle \frac{{\omega }_{1}^{{\prime} }}{2}{\sigma }_{1}^{z}+\displaystyle \frac{{\omega }_{2}^{{\prime} }}{2}{\sigma }_{2}^{z}+\lambda ^{\prime} {\sigma }_{1}^{x}{\sigma }_{2}^{x},\end{eqnarray} \tag{ 1 }$$

where ${\omega }_{1}^{{\prime} }$ and ${\omega }_{2}^{{\prime} }$ are the frequencies of respectively the first and second qubit, and $\lambda ^{\prime}$ is the qubit–qubit coupling constant. We note that for the sake of generality we do not approximate the interaction by σ₁⁺ σ₂^- + h.c., as in RWA. The generality of equation (1) is further discussed in appendix A.

In order to work with dimensionless units, we renormalize the above Hamiltonian by ${\omega }_{1}^{{\prime} }$, i.e. by the frequency of the first qubit:

$$\begin{eqnarray}&&{H}_{S}=\displaystyle \frac{{\omega }_{1}}{2}{\sigma }_{1}^{z}+\displaystyle \frac{{\omega }_{2}}{2}{\sigma }_{2}^{z}+\lambda {\sigma }_{1}^{x}{\sigma }_{2}^{x},\ \ \mathrm{with}\ {\omega }_{1}\,=\,1,\end{eqnarray} \tag{ 2 }$$

${\omega }_{2}={\omega }_{2}^{{\prime} }/{\omega }_{1}^{{\prime} }$ and $\lambda =\lambda ^{\prime} /{\omega }_{1}^{{\prime} }$. Through all the work, we assume that the renormalized qubit frequencies are of the same order, ω₂ = 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):

$$\begin{eqnarray}&&H={H}_{S}+{H}_{B}^{({l}_{1})}+{H}_{B}^{({l}_{2})}+{H}_{B}^{(c)}+{H}_{I},\end{eqnarray} \tag{ 3 }$$

with

$$\begin{eqnarray}&&{H}_{I}=\left({g}_{x}^{({l}_{1})}{\sigma }_{1}^{x}+{g}_{z}^{({l}_{1})}{\sigma }_{1}^{z}\right){B}^{({l}_{1})}+\left({g}_{x}^{({l}_{2})}{\sigma }_{2}^{x}+{g}_{z}^{({l}_{2})}{\sigma }_{2}^{z}\right){B}^{({l}_{2})}+\ \left({g}_{x}^{({c}_{1})}{\sigma }_{1}^{x}+{g}_{x}^{({c}_{2})}{\sigma }_{2}^{x}+{g}_{z}^{({c}_{1})}{\sigma }_{1}^{z}+{g}_{z}^{({c}_{2})}{\sigma }_{2}^{z}\right){B}^{(c)}.\end{eqnarray} \tag{ 4 }$$

${H}_{B}^{({l}_{1})}$ and ${H}_{B}^{({l}_{2})}$ refer to the free Hamiltonian of the local bath respectively on the first and on the second qubit, while ${H}_{B}^{(c)}$ is the free Hamiltonian of the common bath. They read:

$$\begin{eqnarray}&&{H}_{B}^{(\alpha )}=\displaystyle \sum _{k}{{\rm{\Omega }}}_{k,\alpha }{a}_{k,\alpha }^{\dagger }{a}_{k,\alpha },\end{eqnarray} \tag{ 5 }$$

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₁, l₂, c indicates the specific bath.

The bath operators appearing in the interaction Hamiltonian H_I are given by

$$\begin{eqnarray}&&{B}^{(\alpha )}=\displaystyle \sum _{k}{f}_{k,\alpha }\left({a}_{k,\alpha }^{\dagger }+{a}_{k,\alpha }\right),\end{eqnarray} \tag{ 6 }$$

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. ${g}_{j}^{(a)}=O(\mu )$ $\forall a={l}_{1},{l}_{2},{c}_{1},{c}_{2}$ and j = x, z, and we assume μ ≪ 1 (consistently in units of ω₁).

The coupling coefficients f_k,α define the spectral density J_α(ω) of each bath through:

$$\begin{eqnarray}&&{J}_{\alpha }(\omega )=\displaystyle \sum _{k}{f}_{k,\alpha }^{2}\delta (\omega -{{\rm{\Omega }}}_{k,\alpha }),\end{eqnarray} \tag{ 7 }$$

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.

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 equation in the (partial) secular approximation. The possibility for a local master equation will be discussed in section 3.

Let us work in the interaction picture according to the free Hamiltonian ${H}_{0}={H}_{S}+{H}_{B}$, where the full bath Hamiltonian is ${H}_{B}={H}_{B}^{({l}_{1})}+{H}_{B}^{({l}_{2})}+{H}_{B}^{(c)}$. The Von-Neumann equation thus reads

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\rho (t)=-{\rm{i}}[{H}_{I}(t),\rho (t)],\end{eqnarray} \tag{ 8 }$$

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 ${\rho }_{S}(t)={\mathrm{Tr}}_{B}[\rho (t)]$:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\rho }_{S}(t)=-{\int }_{0}^{t}{\rm{d}}{t}^{\prime} \,{\mathrm{Tr}}_{B}[{H}_{I}(t),[{H}_{I}(t^{\prime} ),\rho (t^{\prime} )]],\end{eqnarray} \tag{ 9 }$$

where [H_I(t), ρ(0)] = 0, 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:

$$\begin{eqnarray}&&\rho (t)\approx {\rho }_{S}(t)\otimes {\rho }_{B}.\end{eqnarray} \tag{ 10 }$$

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(μ³), where μ is the coupling constant defined in the previous section. Therefore

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\rho }_{S}(t)=-{\int }_{0}^{t}{\rm{d}}{t}^{\prime} {\mathrm{Tr}}_{B}[{H}_{I}(t),[{H}_{I}(t^{\prime} ),{\rho }_{S}(t^{\prime} )\otimes {\rho }_{B}]]+O({\mu }^{3}).\end{eqnarray} \tag{ 11 }$$

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:

$$\begin{eqnarray}&&{H}_{I}(t)=\displaystyle \sum _{\beta }{A}_{\beta }(t)\otimes {B}_{\beta }(t),\end{eqnarray} \tag{ 12 }$$

where A_β(t) are system operators, while B_β(t) are bath operators⁴ . If we make the change of variable $\tau =t-t^{\prime}$ and insert equation (12) in (11), after some algebra we obtain:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\rho }_{S}(t)=-\displaystyle \sum _{\beta ,\beta ^{\prime} }{\displaystyle \int }_{0}^{t}{\rm{d}}\tau \left({{ \mathcal B }}_{\beta \beta ^{\prime} }(\tau )[{A}_{\beta }(t),{A}_{\beta ^{\prime} }(t-\tau ){\rho }_{S}(t-\tau )]+{\rm{h}}.{\rm{c}}.\right)+\ O({\mu }^{3}),\end{eqnarray} \tag{ 13 }$$

having introduced the bath correlation function ${{ \mathcal B }}_{\beta \beta ^{\prime} }(\tau )=\langle {B}_{\beta }(\tau ){B}_{\beta ^{\prime} }(0){\rangle }_{B}=\mathrm{Tr}[{B}_{\beta }(\tau ){B}_{\beta ^{\prime} }(0){\rho }_{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 $\left|{B}_{\beta \beta ^{\prime} }(\tau )\right|\sim {{\rm{e}}}^{-\tau /{\tau }_{B}}$. 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

$$\begin{eqnarray}&&{\tau }_{R}=O({\mu }^{-2}),\end{eqnarray} \tag{ 14 }$$

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\to \infty$, 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\to \infty$ the autocorrelation functions of the bath are proportional to a Dirac delta, ${{ \mathcal B }}_{\beta \beta ^{\prime} }(\tau )\propto \delta (\tau )$, 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 of ${{ \mathcal B }}_{\beta \beta ^{\prime} }(\tau )$. 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

$$\begin{eqnarray}&&{\tau }_{B}\ll {t}^{* }\ll {\tau }_{R}=O({\mu }^{-2}).\end{eqnarray} \tag{ 15 }$$

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:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\rho }_{S}(t)=-\displaystyle \sum _{\beta ,\beta ^{\prime} }{\displaystyle \int }_{0}^{\infty }{\rm{d}}\tau \left({{ \mathcal B }}_{\beta \beta ^{\prime} }(\tau )[{A}_{\beta }(t),{A}_{\beta ^{\prime} }(t-\tau ){\rho }_{S}(t)]+{\rm{h}}.{\rm{c}}.\right)+\ o({\mu }^{2}).\end{eqnarray} \tag{ 16 }$$

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 ${\tau }_{B}\ll {\mu }^{-1}$. In general, assuming the condition in equation (15), the approximated master equation neglects a remainder of order higher than O(μ²) (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_β:

$$\begin{eqnarray}&&{A}_{\beta }(\omega )=\displaystyle \sum _{\epsilon ^{\prime} -\epsilon =\omega }| \epsilon \rangle \langle \epsilon | {A}_{\beta }| \epsilon ^{\prime} \rangle \langle \epsilon ^{\prime} | ,\end{eqnarray} \tag{ 17 }$$

where $\{| \epsilon \rangle \}{}_{\epsilon }$ is the basis of the eigenvectors of the system Hamiltonian H_S. The following properties hold:

$$\begin{eqnarray}&&{A}_{\beta }^{\dagger }(\omega )={A}_{\beta }(-\omega ),\ \displaystyle \sum _{\omega }{A}_{\beta }(\omega )=\displaystyle \sum _{\omega }{A}_{\beta }^{\dagger }(\omega )={A}_{\beta }.\end{eqnarray} \tag{ 18 }$$

By writing equation (16) with the time-evolved jump operators, we get to the Bloch–Redfield master equation

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\rho }_{S}(t)=\displaystyle \sum _{\omega ,\omega ^{\prime} }\displaystyle \sum _{\beta ,\beta ^{\prime} }{{\rm{e}}}^{{\rm{i}}(\omega ^{\prime} -\omega )t}{{\rm{\Gamma }}}_{\beta \beta ^{\prime} }(\omega )\left({A}_{\beta ^{\prime} }(\omega ){\rho }_{S}(t){A}_{\beta }^{\dagger }(\omega ^{\prime} )-{A}_{\beta }^{\dagger }(\omega ^{\prime} ){A}_{\beta ^{\prime} }(\omega ){\rho }_{S}(t)\right)+{\rm{h}}.{\rm{c}}.,\end{eqnarray} \tag{ 19 }$$

where we have introduced the one-side Fourier transform of the bath correlation functions

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\beta \beta ^{\prime} }(\omega )={\int }_{0}^{\infty }{\rm{d}}{t}^{\prime} {{\rm{e}}}^{{\rm{i}}\omega t^{\prime} }{{ \mathcal B }}_{\beta \beta ^{\prime} }(t^{\prime} ).\end{eqnarray} \tag{ 20 }$$

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 $\omega ^{\prime}$ and ω in equation (19) being coarse-grained in time as from equation (15), i.e.

$$\begin{eqnarray}&&\exists \,{t}^{* }{\rm{such}}\,{\rm{that}}{\left|\omega ^{\prime} -\omega \right|}^{-1}\ll {t}^{* }\ll {\tau }_{R}=O({\mu }^{-2}),\end{eqnarray} \tag{ 21 }$$

then the terms in equation (19) oscillating with frequency $\omega ^{\prime} -\omega$ will not give any significant contribution to the system evolution, since by integrating equation (19) for a time t* such that ${\left|\omega ^{\prime} -\omega \right|}^{-1}\ll {t}^{* }\ll {\tau }_{R}$ the fast-oscillating 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 $\omega ^{\prime} \ne \omega$, 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. $\left|\omega ^{\prime} -\omega \right|\gg {\tau }_{R}^{-1}\approx {\mu }^{2}$ for any $\omega ^{\prime} ,\omega$, 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 $\omega ^{\prime} \ne \omega$, while we will employ the name partial secular approximation for the cases in which we keep some slowly rotating terms with $\omega ^{\prime} \ne \omega$, 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

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\rho }_{S}(t)=-{\rm{i}}[{H}_{S}+{H}_{{LS}},{\rho }_{S}(t)]+D[{\rho }_{S}(t)],\end{eqnarray} \tag{ 22 }$$

where we have introduced the Lamb-Shift Hamiltonian:

$$\begin{eqnarray}&&{H}_{{LS}}=\displaystyle \sum _{\omega ,\omega ^{\prime} }\displaystyle \sum _{\beta ,\beta ^{\prime} }{S}_{\beta \beta ^{\prime} }(\omega ,\omega ^{\prime} ){A}_{\beta }^{\dagger }(\omega ^{\prime} ){A}_{\beta ^{\prime} }(\omega ),\end{eqnarray} \tag{ 23 }$$

and the dissipator of the master equation, responsible for the energy losses of the system:

$$\begin{eqnarray}&&{ \mathcal D }({\rho }_{S})=\displaystyle \sum _{\omega ,\omega ^{\prime} }\displaystyle \sum _{\beta ,\beta ^{\prime} }{\gamma }_{\beta \beta ^{\prime} }(\omega ,\omega ^{\prime} )\left({A}_{\beta ^{\prime} }(\omega ){\rho }_{S}{A}_{\beta }^{\dagger }(\omega ^{\prime} )-\displaystyle \frac{1}{2}\{{A}_{\beta }^{\dagger }(\omega ^{\prime} ){A}_{\beta ^{\prime} }(\omega ),{\rho }_{S}\}\right),\end{eqnarray} \tag{ 24 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}{S}_{\beta \beta ^{\prime} }(\omega ,\omega ^{\prime} ) & = & \displaystyle \frac{{{\rm{\Gamma }}}_{\beta \beta ^{\prime} }(\omega )-{{\rm{\Gamma }}}_{\beta ^{\prime} \beta }^{* }(\omega ^{\prime} )}{2i},\\ {\gamma }_{\beta \beta ^{\prime} }(\omega ,\omega ^{\prime} ) & = & {{\rm{\Gamma }}}_{\beta \beta ^{\prime} }(\omega )+{{\rm{\Gamma }}}_{\beta ^{\prime} \beta }^{* }(\omega ^{\prime} ).\end{array}\end{eqnarray} \tag{ 25 }$$

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:

$$\begin{eqnarray}\begin{array}{rcl}{H}_{{LS}} & = & \displaystyle \sum _{\omega }\displaystyle \sum _{\beta ,\beta ^{\prime} }{S}_{\beta \beta ^{\prime} }(\omega ,\omega ){A}_{\beta }^{\dagger }(\omega ){A}_{\beta ^{\prime} }(\omega ),\\ { \mathcal D }({\rho }_{S}) & = & \displaystyle \sum _{\omega }\displaystyle \sum _{\beta ,\beta ^{\prime} }{\gamma }_{\beta \beta ^{\prime} }(\omega ,\omega )\left({A}_{\beta ^{\prime} }(\omega ){\rho }_{S}{A}_{\beta }^{\dagger }(\omega )-\displaystyle \frac{1}{2}\{{A}_{\beta }^{\dagger }(\omega ){A}_{\beta ^{\prime} }(\omega ),{\rho }_{S}\}\right).\end{array}\end{eqnarray} \tag{ 26 }$$

The master equation (22) with Lamb-Shift Hamiltonian and dissipator given by equation (26) is written in the GKLS form [38–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 terms $\left({A}_{\beta ^{\prime} }(\omega ){\rho }_{S}{A}_{\beta }(\omega ^{\prime} )-\tfrac{1}{2}\{{A}_{\beta }(\omega ^{\prime} ){A}_{\beta ^{\prime} }(\omega ),{\rho }_{S}\}\right)$, 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 $| \omega -(-{\omega }^{{\prime} })|$ (again with ω and ω' with the same sign): if this is the case, then all terms ${A}_{\beta ^{\prime} }(\omega ){\rho }_{S}{A}_{\beta }^{\dagger }(-\omega ^{\prime} )\equiv {A}_{\beta ^{\prime} }(\omega ){\rho }_{S}{A}_{\beta }(\omega ^{\prime} )$ 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.

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.

2.3.1. No direct coupling

Let us start with the simplest case, i.e. in the absence of a direct coupling between the qubits: setting λ = 0 in equation (2), the system Hamiltonian reads:

$$\begin{eqnarray}&&{H}_{S}=\displaystyle \frac{{\omega }_{1}}{2}{\sigma }_{1}^{z}+\displaystyle \frac{{\omega }_{2}}{2}{\sigma }_{2}^{z},\end{eqnarray} \tag{ 27 }$$

which is already diagonal in the 'canonical basis' $\{| 00\rangle ,| 01\rangle ,| 10\rangle ,| 11\rangle \}$, with eigenvalues respectively E₀ = −ω₊ /2, ${E}_{1}=-{\omega }_{-}/2$, ${E}_{2}={\omega }_{-}/2$, ${E}_{3}={\omega }_{+}/2$, where ${\omega }_{\pm }={\omega }_{1}\pm {\omega }_{2}$ and (without losing generality) we set ω₁ > ω₂.

In the interaction Hamiltonian we can find the system operators σ_j^x and ${\sigma }_{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):

$$\begin{eqnarray}&&{\sigma }_{j}^{x}({\omega }_{j})={\sigma }_{j}^{-},{\sigma }_{j}^{x}(-{\omega }_{j})={\sigma }_{j}^{+}\ \Rightarrow \ {\sigma }_{j}^{x}={\sigma }_{j}^{-}+{\sigma }_{j}^{+},{\sigma }_{j}^{z}(0)={\sigma }_{j}^{z}\ \Rightarrow \ {\sigma }_{j}^{z}={\sigma }_{j}^{z},\end{eqnarray} \tag{ 28 }$$

with j = 1, 2. Equation (28) describes the possible emission and absorption processes of the system, depicted in figure 1(b).

A self-consistent secular approximation depends on the detuning between the qubits. In the case in which there is a small detuning, such that ω₁ − ω₂ is not way greater than μ², we cannot employ the full secular approximation, but we need to rely on a partial secular approximation in which we keep in equation (19) slow terms with $\omega ^{\prime} -\omega =\pm ({\omega }_{1}-{\omega }_{2})$. 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 non-identical 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.

Zoom In Zoom Out Reset image size

Starting from equation (19) and employing the notation of equations (23) and (24), the master equation finally reads:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\rho }_{S}(t) & = & -{\rm{i}}[{H}_{S}+{H}_{{LS}},{\rho }_{S}(t)]+\mathop{\underbrace{\displaystyle \sum _{j,k=1,2}{\gamma }_{{jk}}\left({\sigma }_{j}^{-}{\rho }_{S}(t){\sigma }_{k}^{+}-\displaystyle \frac{1}{2}\{{\sigma }_{k}^{+}{\sigma }_{j}^{-},{\rho }_{S}(t)\}\right)}}\limits_{{\rm{emission}}}\\ & & +\ \mathop{\underbrace{\displaystyle \sum _{j,k=1,2}{\tilde{\gamma }}_{{jk}}\left({\sigma }_{j}^{+}{\rho }_{S}(t){\sigma }_{k}^{-}-\displaystyle \frac{1}{2}\{{\sigma }_{k}^{-}{\sigma }_{j}^{+},{\rho }_{S}(t)\}\right)}}\limits_{{\rm{absorption}}}+\mathop{\underbrace{\displaystyle \sum _{j,k=1,2}{\eta }_{{jk}}\left({\sigma }_{j}^{z}{\rho }_{S}(t){\sigma }_{k}^{z}-\displaystyle \frac{1}{2}\{{\sigma }_{k}^{z}{\sigma }_{j}^{z},{\rho }_{S}(t)\}\right)}}\limits_{{\rm{dephasing}}},\end{array}\end{eqnarray} \tag{ 29 }$$

where the Lamb-shift Hamiltonian is given by:

$$\begin{eqnarray}&&{H}_{{LS}}=\displaystyle \sum _{{jk}}\left({s}_{{jk}}{\sigma }_{k}^{+}{\sigma }_{j}^{-}+{\tilde{s}}_{{jk}}{\sigma }_{k}^{-}{\sigma }_{j}^{+}\right)+2{s}_{0}{\sigma }_{1}^{z}{\sigma }_{2}^{z}=\displaystyle \frac{{s}_{1}}{2}{\sigma }_{1}^{z}+\displaystyle \frac{{s}_{2}}{2}{\sigma }_{2}^{z}+{s}_{+}{\sigma }_{1}^{+}{\sigma }_{2}^{-}+{s}_{-}{\sigma }_{1}^{-}{\sigma }_{2}^{+}+2{s}_{0}{\sigma }_{1}^{z}{\sigma }_{2}^{z},\end{eqnarray} \tag{ 30 }$$

and the coefficients of the master equation are presented in appendix B.

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 ${{\rm{\Gamma }}}_{\beta \beta ^{\prime} }(\omega )$ (see appendix B for their specific form): then we would have ${{\rm{\Gamma }}}_{\beta \beta ^{\prime} }(\omega )=0$ for any system operators A_β and ${A}_{\beta ^{\prime} }$ coupled to distinct baths through B_β and ${B}_{\beta ^{\prime} }$, since there are no correlations between the two baths. On the contrary, if B_β and ${B}_{\beta ^{\prime} }$ are operators of the same bath, then ${{\rm{\Gamma }}}_{\beta \beta ^{\prime} }(\omega )$ does not vanish a priori and there could be a case in which ${A}_{\beta }\ne {A}_{\beta ^{\prime} }$ but ${{\rm{\Gamma }}}_{\beta \beta ^{\prime} }(\omega )\ne 0$. Let us for instance consider the coupling with the local bath on the first qubits: ${H}_{I}^{({l}_{1})}=\left({g}_{x}^{({l}_{1})}{\sigma }_{1}^{x}+{g}_{z}^{({l}_{1})}{\sigma }_{1}^{z}\right){B}^{({l}_{1})}$.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 ${\sigma }_{1}^{x}$ and ${\sigma }_{1}^{z}$ may in theory couple in equation (19) with a non-zero coefficient ${{\rm{\Gamma }}}_{{xz}}^{({l}_{1})}({\omega }_{1})$ or ${{\rm{\Gamma }}}_{{xz}}^{({l}_{1})}(0)$, but their corresponding terms would vanish because of the partial secular approximation, since $\left|\omega ^{\prime} -\omega \right|=\left|\pm {\omega }_{1}\right|\gg 1/{\tau }_{R}$. 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⁵ . 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).

2.3.2. 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:

$$\begin{eqnarray}&&{H}_{S}=\displaystyle \frac{{\omega }_{1}}{2}{\sigma }_{1}^{z}+\displaystyle \frac{{\omega }_{2}}{2}{\sigma }_{2}^{z}+\lambda {\sigma }_{1}^{x}{\sigma }_{2}^{x},\end{eqnarray} \tag{ 31 }$$

and the corresponding matrix in the canonical basis $\{| 11\rangle ,| 10\rangle ,| 01\rangle ,| 00\rangle \}$ is written as

$$\begin{eqnarray}{H}_{S}=\left(\begin{array}{cccc}{\omega }_{+}/2 & 0 & 0 & \lambda \\ 0 & {\omega }_{-}/2 & \lambda & 0\\ 0 & \lambda & -{\omega }_{-}/2 & 0\\ \lambda & 0 & 0 & -{\omega }_{+}/2\end{array}\right),\end{eqnarray} \tag{ 32 }$$

with ω_± = ω₁ ± ω₂.

We can easily diagonalize equation (32) by finding the eigenvalues:

$$\begin{eqnarray}&&{E}_{0}=-\sqrt{{\lambda }^{2}+{\omega }_{+}^{2}/4},\ {E}_{1}=-\sqrt{{\lambda }^{2}+{\omega }_{-}^{2}/4},{E}_{2}=+\sqrt{{\lambda }^{2}+{\omega }_{-}^{2}/4},\ {E}_{3}=+\sqrt{{\lambda }^{2}+{\omega }_{+}^{2}/4},\end{eqnarray} \tag{ 33 }$$

with associated eigenvectors [19]

$$\begin{eqnarray}\begin{array}{rcl}| {e}_{0}\rangle & = & -\sin \theta | 11\rangle +\cos \theta | 00\rangle ,\ | {e}_{1}\rangle =-\sin \phi | 10\rangle +\cos \phi | 01\rangle ,\\ | {e}_{2}\rangle & = & +\cos \phi | 10\rangle +\sin \phi | 01\rangle ,\ | {e}_{3}\rangle =+\cos \theta | 11\rangle +\sin \theta | 00\rangle ,\end{array}\end{eqnarray} \tag{ 34 }$$

where the parameters θ and ϕ are given by

$$\begin{eqnarray}&&\sin 2\theta =\displaystyle \frac{\lambda }{{E}_{3}},\cos 2\theta =\displaystyle \frac{{\omega }_{+}/2}{{E}_{3}},\sin 2\phi =\displaystyle \frac{\lambda }{{E}_{2}},\ \cos 2\phi =\displaystyle \frac{{\omega }_{-}/2}{{E}_{2}}.\end{eqnarray} \tag{ 35 }$$

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. ${\sigma }_{j}^{x}$ and ${\sigma }_{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 ${\omega }_{-}\ll 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.

With these prescriptions, the master equation reads:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\rho }_{S}(t) & = & -{\rm{i}}[{H}_{S}+{H}_{{LS}},{\rho }_{S}(t)]+\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{j,k=I,{II}}{m,n=1,2}}{\gamma }_{{jk}}^{{mn}}\left({\sigma }_{m}^{x}({\omega }_{j}){\rho }_{S}(t){\sigma }_{n}^{x}(-{\omega }_{k})-\displaystyle \frac{1}{2}\{{\sigma }_{n}^{x}(-{\omega }_{k}){\sigma }_{m}^{x}({\omega }_{j}),{\rho }_{S}(t)\}\right)\\ & & +\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{j,k=I,{II}}{m,n=1,2}}{\tilde{\gamma }}_{{jk}}^{{mn}}\left({\sigma }_{m}^{x}(-{\omega }_{j}){\rho }_{S}(t){\sigma }_{n}^{x}({\omega }_{k})-\displaystyle \frac{1}{2}\{{\sigma }_{n}^{x}({\omega }_{k}){\sigma }_{m}^{x}(-{\omega }_{j}),{\rho }_{S}(t)\}\right)\\ & & +\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{j,k=0,\pm {IV}}{m,n=1,2}}{\eta }_{{jk}}^{{mn}}\left({\sigma }_{m}^{z}({\omega }_{j}){\rho }_{S}(t){\sigma }_{n}^{z}(-{\omega }_{k})-\displaystyle \frac{1}{2}\{{\sigma }_{n}^{z}(-{\omega }_{k}){\sigma }_{m}^{z}({\omega }_{j}),{\rho }_{S}(t)\}\right)\\ & & +\ \displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{j=\pm {III}}{m,n=1,2}}{\zeta }_{j}^{{mn}}\left({\sigma }_{m}^{z}({\omega }_{j}){\rho }_{S}(t){\sigma }_{n}^{z}(-{\omega }_{j})-\displaystyle \frac{1}{2}\{{\sigma }_{n}^{z}(-{\omega }_{j}){\sigma }_{m}^{z}({\omega }_{j}),{\rho }_{S}(t)\}\right),\end{array}\end{eqnarray} \tag{ 36 }$$

where the jump operators and relative frequencies ${\sigma }_{m}^{x}({\omega }_{j})$ are defined in equations (C.1) and (C.2), and we are using the short notation ${\omega }_{-{IV}}=-{\omega }_{{IV}}$, ${\omega }_{-{III}}=-{\omega }_{{III}}$ and ω₀ = 0. The Lamb-Shift Hamiltonian is given by:

$$\begin{eqnarray}\begin{array}{rcl}{H}_{{LS}} & = & \displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{j,k=I,{II}}{m,n=1,2}}\left({s}_{{jk}}^{{mn}}{\sigma }_{n}^{x}(-{\omega }_{k}){\sigma }_{m}^{x}({\omega }_{j})+{\tilde{s}}_{{jk}}^{{mn}}{\sigma }_{n}^{x}({\omega }_{k}){\sigma }_{m}^{x}(-{\omega }_{j}\right)+\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{j,k=0,\pm {IV}}{m,n=1,2}}{r}_{{jk}}^{{mn}}{\sigma }_{n}^{z}(-{\omega }_{k}){\sigma }_{m}^{x}({\omega }_{j})\\ & & +\ \displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{j=\pm {III}}{m,n=1,2}}{u}_{j}^{{mn}}{\sigma }_{n}^{z}(-{\omega }_{j}){\sigma }_{m}^{z}({\omega }_{j}).\end{array}\end{eqnarray} \tag{ 37 }$$

The coefficients of the master equation are listed in equation (B.2) in appendix B.

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₁ + H₂ + H₁₂ 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 ${H}_{12}=\lambda {\sigma }_{1}^{x}{\sigma }_{2}^{x}$). 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) using as free system Hamiltonian ${H}_{S}^{{\rm{local}}}={H}_{1}+{H}_{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 non-locality 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 $\omega ^{\prime} \ne \omega$, 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–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.

In order to assess the accuracy of the local master equation, let us write the interaction Hamiltonian as H₁₂ = λ V, where λ is a 'small' parameter which we consider as a perturbation order, i.e. λ ≪ 1. The system Hamiltonian reads ${H}_{S}={H}_{1}+{H}_{2}+{H}_{12}={H}_{S}^{{\rm{local}}}+\lambda V$. Following [26], we apply standard perturbation theory to find the zeroth order eigenvectors and eigenvalues⁶ . Within the assumption of not-degenerate system Hamiltonian, we write the eigenvalues as the infinite perturbation expansion [82]:

$$\begin{eqnarray}\begin{array}{rcl}{E}_{n} & = & {E}_{n}^{(0)}+\lambda {E}_{n}^{(1)}+{\lambda }^{2}{E}_{n}^{(2)}+...\\ | {e}_{n}\rangle & = & | {e}_{n}^{(0)}\rangle +\lambda | {e}_{n}^{(1)}\rangle +{\lambda }^{2}| {e}_{n}^{(2)}\rangle +...,\end{array}\end{eqnarray} \tag{ 38 }$$

where ${E}_{n}^{(0)}$ and $| {e}_{n}^{(0)}\rangle$ are respectively the eigenvalues and eigenvectors of the unperturbed Hamiltonian, in our case of the local Hamiltonian ${H}_{S}^{{\rm{local}}}$. A jump operator equation (17) using the expansions in equation (38) will read

$$\begin{eqnarray}&&{A}_{\beta }(\omega )=\displaystyle \sum _{{E}_{n}^{{\prime} }-{E}_{n}=\omega }| {e}_{n}\rangle \langle {e}_{n}| {A}_{\beta }| {e}_{n}^{{\prime} }\rangle \langle {e}_{n}^{{\prime} }| =\ \mathop{\underbrace{\displaystyle \sum _{{E}_{n}^{{\prime} }-{E}_{n}=\omega }| {e}_{n}^{\left(0\right)}\rangle \langle {e}_{n}^{\left(0\right)}| {A}_{\beta }| e{{\prime} }_{n}^{\left(0\right)}\rangle \langle e{{\prime} }_{n}^{\left(0\right)}| }}\limits_{{A}_{\beta }^{\left(0\right)}({\omega }^{(0)})}+O(\lambda ),\end{eqnarray} \tag{ 39 }$$

where ${A}_{\beta }^{\left(0\right)}({\omega }^{(0)})$ are the local jump operators appearing in the local master equation.

While [26] considers the GKLS master equation in full secular approximation, here we derive the Bloch–Redfield local master equation. Inserting equation (39) in equation (19) and coming back to the Schrödinger picture we obtain:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\rho }_{S}(t) & = & -{\rm{i}}[{H}_{S}^{{\rm{local}}}+\lambda V,{\rho }_{S}(t)]+\ \displaystyle \sum _{{\omega }^{(0)},\omega {{\prime} }^{(0)}}\displaystyle \sum _{\beta ,\beta ^{\prime} }{{\rm{\Gamma }}}_{\beta \beta ^{\prime} }({\omega }^{(0)})\left({A}_{\beta ^{\prime} }^{(0)}({\omega }^{(0)}){\rho }_{S}(t){\left({A}_{\beta }^{\left(0\right)}\right)}^{\dagger }(\omega {{\prime} }^{(0)})\right.\\ & & -\ \left.{\left({A}_{\beta }^{\left(0\right)}\right)}^{\dagger }(\omega {{\prime} }^{(0)}){A}_{\beta ^{\prime} }^{\left(0\right)}({\omega }^{(0)}){\rho }_{S}(t)\right)+\ h.c.+O({\mu }^{2}\lambda )+o({\mu }^{2}),\end{array}\end{eqnarray} \tag{ 40 }$$

where ω⁽⁰⁾ is the frequency given by differences of the unperturbed eigenvalues ${E}_{n}^{\left(0\right)}-{E}_{m}^{\left(0\right)}$. The error we are making by employing the local master equation is of the order of O(μ²λ), since the leading order of the master equation after the Born approximation is O(μ²). Moreover, we explicitly write the order of the remainder after the Born–Markov approximations o(μ²), to stress that there are already some neglected terms which may be larger than the error given by the local master equation O(μ²λ).

If we employ the local master equation to compute physical quantities, clearly we can resolve them only up to the order of O(μ²). Any quantity of the order of O(μ²λ) 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(μ²λ²).

We will now address the local versus global comparison focusing on the case of two spins as in equation (2), with 'local' Hamiltonian ${H}_{S}^{{\rm{local}}}=\tfrac{{\omega }_{1}}{2}{\sigma }_{1}^{z}+\tfrac{{\omega }_{2}}{2}{\sigma }_{2}^{z}$ and interaction ${H}_{12}=\lambda {\sigma }_{1}^{x}{\sigma }_{2}^{x}$. 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.

3.3.1. 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:

$$\begin{eqnarray}&&{E}_{0}=-\displaystyle \frac{{\omega }_{+}}{2}+O(\lambda ),\ {E}_{1}=-\displaystyle \frac{{\omega }_{-}}{2}+O(\lambda ),{E}_{2}=+\displaystyle \frac{{\omega }_{-}}{2}+O(\lambda ),\ {E}_{3}=+\displaystyle \frac{{\omega }_{+}}{2}+O(\lambda ),\end{eqnarray} \tag{ 41 }$$

with corresponding eigenvectors

$$\begin{eqnarray}&&| {e}_{0}\rangle =| 00\rangle +O(\lambda ),\ | {e}_{1}\rangle =| 01\rangle +O(\lambda ),| {e}_{2}\rangle =| 01\rangle +O(\lambda ),\ | {e}_{3}\rangle =| 11\rangle +O(\lambda ).\end{eqnarray} \tag{ 42 }$$

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(μ²λ).

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}\rangle$ and $| {e}_{2}\rangle$ 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, assuming $\lambda \gg {\omega }_{-}\ne 0$, we would have for instance $| {e}_{2}\rangle =1/\sqrt{2}(| 01\rangle +| 10\rangle )+O(\lambda )$, but the aforementioned selection rule would still apply and the final master equation would not change.

3.3.2. Comparison between master equations

3.3.3. 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 ${\sigma }_{1}^{x}$ 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 ${\tau }_{R}^{-1}=O({\mu }^{2})$, always avoiding the singular parameter choice expressed in equation (C.4). Critical cases are:
  $$\begin{eqnarray}&&{\omega }_{{II}}-{\omega }_{I}={\omega }_{{IV}}-(-{\omega }_{{IV}})=2({\omega }_{{IV}}-{\omega }_{0})=2\sqrt{{\lambda }^{2}+{\omega }_{-}^{2}/4}.\end{eqnarray} \tag{ 43 }$$
  If the qubits have small detuning, i.e. $\left|{\omega }_{-}\right|\ll 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 $\left|{\omega }_{{II}}-{\omega }_{I}\right|\lessapprox {\mu }^{2}$. 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(μ²λ). 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 $\langle {\sigma }_{1}^{z}(t)\rangle$ 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.3.4. 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 ${\tau }_{R}^{-1}=O({\mu }^{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 ${\omega }_{-}\gg {\mu }^{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(μ²) or λ ≫ O(μ²).

Figure 5 shows some relevant examples through $\langle {\sigma }_{1}^{z}\rangle$ and the fidelity of the states obtained with different master equations compared with the GP: if λ is very small (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 λ ≫ μ².

Zoom In Zoom Out Reset image size

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⁻², λ = 10⁻⁴ 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⁻². 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 well-known 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 ${ \mathcal N }$ [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.

Zoom In Zoom Out Reset image size

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.

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 ${\beta }^{({l}_{1})}=1$ (colder), ${\beta }^{({l}_{2})}=0.1$ (hotter). If ${\rho }_{\infty }$ is the steady state of the open system, the heat current from the hotter reservoir is defined as [89]:

$$\begin{eqnarray}&&{J}_{2}=\mathrm{Tr}\left({H}_{S}{{ \mathcal D }}_{2}[{\rho }_{\infty }]\right),\end{eqnarray} \tag{ 44 }$$

where H_S is the system Hamiltonian equation (31) and ${{ \mathcal D }}_{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 ${J}_{1}=\mathrm{Tr}\left({H}_{S}{{ \mathcal D }}_{1}[{\rho }_{\infty }]\right)$, and since ${\rho }_{\infty }$ is the steady state we have J₁ + J₂ = 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⁻³ (a) and ω₋ = 10⁻⁵ (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⁻³ than for ω₋ = 10⁻⁵. If the detuning were null, the global approach with full secular approximation would provide a non-zero heat current even for $\lambda \to 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 $\lambda {\sigma }_{1}^{x}{\sigma }_{2}^{x}$, and therefore it commutes with the interaction Hamiltonian: ${\mathrm{lim}}_{\lambda \to \infty }[{H}_{S},{H}_{I}]=0$. 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 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.

Zoom In Zoom Out Reset image size

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.

Following the standard convention in quantum information theory, we have chosen to write throughout all the paper the free Hamiltonian of a single qubit as H₁ = ω₁/2σ₁^z. Historically, this has not been always the preferred choice, since for instance in the seminal papers by Caldeira and Leggett [92] a different notation was employed, in particular:

$$\begin{eqnarray}&&{H}_{1}^{{\prime} }=\displaystyle \frac{{\epsilon }_{1}}{2}{\sigma }_{1}^{z}+\displaystyle \frac{{{\rm{\Delta }}}_{1}}{2}{\sigma }_{1}^{x},\end{eqnarray} \tag{ A.1 }$$

which regards the two level system as two levels with detuning separated by a potential barrier, with the possibility of hopping through it via σ₁^x.

H₁ and ${H}_{1}^{{\prime} }$ are connected by a unitary transformation that, in absence of interactions between the qubits, will change the dephasing and dissipative character of each bath, and in this work both possibilities are taken into account. On the other hand, in presence of the coupling between qubits, the local unitary transformations will in general transform the coupling ${H}_{12}=\lambda {\sigma }_{1}^{x}{\sigma }_{2}^{x}$ into a more complex interaction term.

In this section we address further possible forms of the qubit–qubit interaction. We remind the system Hamiltonian:

$$\begin{eqnarray}&&{H}_{S}=\displaystyle \frac{{\omega }_{1}}{2}{\sigma }_{1}^{z}+\displaystyle \frac{{\omega }_{2}}{2}{\sigma }_{2}^{z}+{H}_{12},\end{eqnarray} \tag{ A.2 }$$

with ${H}_{12}=\lambda {\sigma }_{1}^{x}{\sigma }_{2}^{x}$. This choice is justified by the fact that such interaction is the standard one employed in the framework of quantum information, to which this paper is mostly devoted; we indeed find the ${\sigma }_{1}^{x}{\sigma }_{2}^{x}$ interaction (or fully equivalently ${\sigma }_{1}^{y}{\sigma }_{2}^{y}$) in many experimental platforms, such as superconducting qubits [93–95] or coupled atomic dipoles [96]. Anyway, let us examine possible alternatives and the connection with this Ising-like coupling.

A.2.1. Heisenberg-type interaction

Let us here consider two qubits coupled through an Heisenberg-type interaction, i.e. ${H}_{12}={\sum }_{k=x,y,z}{\lambda }_{k}{\sigma }_{1}^{k}{\sigma }_{2}^{k}$, quite common in many physical systems. This interaction conserves the parity of the number of excitations and, for this reason, the system Hamiltonian in the canonical basis has the same structure of equation (31):

$$\begin{eqnarray}{H}_{S}=\left(\begin{array}{cccc}{\omega }_{+}/2+{\lambda }_{z} & 0 & 0 & {\lambda }_{x}-{\lambda }_{y}\\ 0 & {\omega }_{-}/2-{\lambda }_{z} & {\lambda }_{x}+{\lambda }_{y} & 0\\ 0 & {\lambda }_{x}+{\lambda }_{y} & -{\omega }_{-}/2-{\lambda }_{z} & 0\\ {\lambda }_{x}-{\lambda }_{y} & 0 & 0 & -{\omega }_{+}/2+{\lambda }_{z}\end{array}\right),\end{eqnarray} \tag{ A.3 }$$

with ω_± = ω₁ ± ω₂.

The parity symmetry makes the eigenvectors remain of the same form of the eigenvectors in equation (34). The only thing that may change is the value of the angles θ and ϕ in equation (35), according to the different values of λ_k appearing in equation (A.3). This means that the general master equation preserves the structure in equation (36): the value of the coefficients marks the only difference between Heisenberg- and Ising-type interaction.

A.2.2. Rotating wave approximation

Another very common case in the literature is to consider the qubit–qubit coupling in RWA, i.e. ${H}_{12}=\lambda ({\sigma }_{1}^{-}{\sigma }_{2}^{+}+{\sigma }_{1}^{+}{\sigma }_{2}^{-})$. Even if in some cases this can be justified, we notice that without counter-rotating terms the Hamiltonian not only conserves the parity of the number of excitations, but also conserves the number of excitations itself. A standard rotation is sufficient to diagonalize it, leading to the eigenvalues:

$$\begin{eqnarray}&&{E}_{3}=+{\omega }_{+}/2,{E}_{2}=+\sqrt{{\lambda }^{2}+{\omega }_{-}^{2}/4},{E}_{1}=-\sqrt{{\lambda }^{2}+{\omega }_{-}^{2}/4},{E}_{0}=-{\omega }_{+}/2,\end{eqnarray} \tag{ A.4 }$$

with associated eigenvectors

$$\begin{eqnarray}&&| {e}_{3}\rangle =| 11\rangle ,| {e}_{2}\rangle =\cos \phi | 10\rangle +\sin \phi | 01\rangle ,| {e}_{1}\rangle =\sin \phi | 10\rangle -\cos \phi | 01\rangle ,| {e}_{0}\rangle =| 00\rangle ,\end{eqnarray} \tag{ A.5 }$$

where the angle ϕ is the same as in equation (35):

$$\begin{eqnarray}&&\sin 2\phi =\displaystyle \frac{\lambda }{{E}_{2}},\ \cos 2\phi =\displaystyle \frac{{\omega }_{-}/2}{{E}_{2}}.\end{eqnarray} \tag{ A.6 }$$

We see that, unless we can neglect the coupling (λ ≪ ω₊), the geometry of the Hamiltonian in the RWA is non-trivially different from the one used in the work. In the master equation, the absence of counter-rotating terms eliminates all the 'double emission' or 'double absorption' jump operators, namely ${\sigma }_{j}^{z}(\pm {\omega }_{{III}})$ in equation (C.2). Without these terms important effects do not arise, such as stationary entanglement [23], or are not properly described, as the refrigerator performance analyzed in [70].

The coefficients of the master equation in equation (29) read:

$$\begin{eqnarray}\begin{array}{rcl}{\gamma }_{{jk}} & = & \left\{\begin{array}{ll}({g}_{x}^{({c}_{1})}{g}_{x}^{({c}_{2})})\left[{{\rm{\Gamma }}}^{(c)}({\omega }_{j})+{\left({{\rm{\Gamma }}}^{(c)}({\omega }_{k})\right)}^{* }\right] & {\rm{if}}\ j\ne k,\\ \displaystyle \sum _{\alpha ={l}_{j},{c}_{j}}{\left({g}_{x}^{(\alpha )}\right)}^{2}\left[{{\rm{\Gamma }}}^{(\alpha )}({\omega }_{j})+{\left({{\rm{\Gamma }}}^{(\alpha )}({\omega }_{j})\right)}^{* }\right] & {\rm{if}}\ j=k,\end{array}\right.\\ {\tilde{\gamma }}_{{jk}} & = & \left\{\begin{array}{ll}({g}_{x}^{({c}_{1})}{g}_{x}^{({c}_{2})})\left[{{\rm{\Gamma }}}^{(c)}(-{\omega }_{j})+{\left({{\rm{\Gamma }}}^{(c)}(-{\omega }_{k}\right)}^{* }\right] & {\rm{if}}\ j\ne k,\\ \displaystyle \sum _{\alpha ={l}_{j},{c}_{j}}{\left({g}_{x}^{(\alpha )}\right)}^{2}\left[{{\rm{\Gamma }}}^{(\alpha )}(-{\omega }_{j})+{\left({{\rm{\Gamma }}}^{(\alpha )}(-{\omega }_{j}\right)}^{* }\right] & {\rm{if}}\ j=k,\end{array}\right.\\ {\eta }_{{jk}} & = & \left\{\begin{array}{ll}({g}_{z}^{({c}_{1})}{g}_{z}^{({c}_{2})})\left[{{\rm{\Gamma }}}^{(c)}(0)+{\left({{\rm{\Gamma }}}^{(c)}(0)\right)}^{* }\right] & {\rm{if}}j\ne k,\\ \displaystyle \sum _{\alpha ={l}_{j},{c}_{j}}{\left({g}_{z}^{(\alpha )}\right)}^{2}\left[{{\rm{\Gamma }}}^{(\alpha )}(0)+{\left({{\rm{\Gamma }}}^{(\alpha )}(0)\right)}^{* }\right] & {\rm{if}}\ j=k,\end{array}\right.\\ {s}_{{jk}} & = & \left\{\begin{array}{ll}({g}_{x}^{({c}_{1})}{g}_{x}^{({c}_{2})})\left[{{\rm{\Gamma }}}^{(c)}({\omega }_{j})-{\left({{\rm{\Gamma }}}^{(c)}({\omega }_{k})\right)}^{* }\right]/2{\rm{i}} & {\rm{if}}\ j\ne k,\\ \displaystyle \sum _{\alpha ={l}_{j},{c}_{j}}{\left({g}_{x}^{(\alpha )}\right)}^{2}\left[{{\rm{\Gamma }}}^{(\alpha )}({\omega }_{j})-{\left({{\rm{\Gamma }}}^{(\alpha )}({\omega }_{j})\right)}^{* }\right]/2{\rm{i}} & {\rm{if}}\ j=k,\end{array}\right.\\ {\tilde{s}}_{{jk}} & = & \left\{\begin{array}{ll}({g}_{x}^{({c}_{1})}{g}_{x}^{({c}_{2})})\left[{{\rm{\Gamma }}}^{(c)}(-{\omega }_{j})-{\left({{\rm{\Gamma }}}^{(c)}(-{\omega }_{k}\right)}^{* }\right]/2{\rm{i}} & {\rm{if}}\ j\ne k,\\ \displaystyle \sum _{\alpha ={l}_{j},{c}_{j}}{\left({g}_{x}^{(\alpha )}\right)}^{2}\left[{{\rm{\Gamma }}}^{(\alpha )}(-{\omega }_{j})-{\left({{\rm{\Gamma }}}^{(\alpha )}(-{\omega }_{j}\right)}^{* }\right]/2{\rm{i}} & {\rm{if}}\ j=k,\end{array}\right.\\ {s}_{0} & = & ({g}_{z}^{({c}_{1})}{g}_{z}^{({c}_{2})})\displaystyle \frac{{{\rm{\Gamma }}}^{(c)}(0)-{\left({{\rm{\Gamma }}}^{(c)}(0)\right)}^{* }}{2{\rm{i}}},\ {s}_{1}={s}_{11}-{\tilde{s}}_{11},\\ {s}_{2} & = & {s}_{22}-{\tilde{s}}_{22},{s}_{+}={s}_{12}+{\tilde{s}}_{21},{s}_{-}={s}_{21}+{\tilde{s}}_{12}.\end{array}\end{eqnarray} \tag{ B.1 }$$

Γ^(α)(ω) are the one-side Fourier transforms of the baths correlation functions defined in equation (20), where with abuse of notation α = c₁ = c₂ = c. The mean value is performed on the thermal state ${\rho }_{B}^{(\alpha )}$ of each bath, at a given inverse temperature β^(α) [31]. g_j^(α) are the coupling constants expressing the strength of the qubit-bath interactions. In the Lamb-Shift Hamiltonian equation (30), s₁/2σ₁^z and ${s}_{2}/2{\sigma }_{2}^{z}$ lead to a renormalization of the qubit frequencies. Moreover, note that we have neglected any multiple of the identity appearing in the Lamb-Shift Hamiltonian.

Analogously to equation (B.1), we list the coefficients appearing in the master equation in equation (36):

$$\begin{eqnarray}\begin{array}{rcl}{\gamma }_{{jk}}^{{mn}} & = & \left\{\begin{array}{ll}({g}_{x}^{({c}_{1})}{g}_{x}^{({c}_{2})})\left[{{\rm{\Gamma }}}^{(c)}({\omega }_{j})+{\left({{\rm{\Gamma }}}^{(c)}({\omega }_{k})\right)}^{* }\right] & {\rm{if}}\ m\ne n,\\ \displaystyle \sum _{\alpha ={l}_{m},{c}_{m}}{\left({g}_{x}^{(\alpha )}\right)}^{2}\left[{{\rm{\Gamma }}}^{(\alpha )}({\omega }_{j})+{\left({{\rm{\Gamma }}}^{(\alpha )}({\omega }_{k})\right)}^{* }\right] & {\rm{if}}\ m=n,\end{array}\right.\\ {\tilde{\gamma }}_{{jk}}^{{mn}} & = & \left\{\begin{array}{ll}({g}_{x}^{({c}_{1})}{g}_{x}^{({c}_{2})})\left[{{\rm{\Gamma }}}^{(c)}(-{\omega }_{j})+{\left({{\rm{\Gamma }}}^{(c)}(-{\omega }_{k}\right)}^{* }\right] & {\rm{if}}\ m\ne n,\\ \displaystyle \sum _{\alpha ={l}_{m},{c}_{m}}{\left({g}_{x}^{(\alpha )}\right)}^{2}\left[{{\rm{\Gamma }}}^{(\alpha )}(-{\omega }_{j})+{\left({{\rm{\Gamma }}}^{(\alpha )}(-{\omega }_{k}\right)}^{* }\right] & {\rm{if}}\ m=n,\end{array}\right.\\ {\eta }_{{jk}}^{{mn}} & = & \left\{\begin{array}{ll}({g}_{z}^{({c}_{1})}{g}_{z}^{({c}_{2})})\left[{{\rm{\Gamma }}}^{(c)}({\omega }_{j})+{\left({{\rm{\Gamma }}}^{(c)}({\omega }_{k})\right)}^{* }\right] & {\rm{if}}\ m\ne n,\\ \displaystyle \sum _{\alpha ={l}_{m},{c}_{m}}{\left({g}_{z}^{(\alpha )}\right)}^{2}\left[{{\rm{\Gamma }}}^{(\alpha )}({\omega }_{j})+{\left({{\rm{\Gamma }}}^{(\alpha )}({\omega }_{k})\right)}^{* }\right] & {\rm{if}}\ m=n,\end{array}\right.\\ {\zeta }_{j}^{{mn}} & = & \left\{\begin{array}{ll}({g}_{z}^{({c}_{1})}{g}_{z}^{({c}_{2})})\left[{{\rm{\Gamma }}}^{(c)}({\omega }_{j})+{\left({{\rm{\Gamma }}}^{(c)}({\omega }_{j})\right)}^{* }\right] & {\rm{if}}\ m\ne n,\\ \displaystyle \sum _{\alpha ={l}_{m},{c}_{m}}{\left({g}_{z}^{(\alpha )}\right)}^{2}\left[{{\rm{\Gamma }}}^{(\alpha )}({\omega }_{j})+{\left({{\rm{\Gamma }}}^{(\alpha )}({\omega }_{j})\right)}^{* }\right] & {\rm{if}}\ m=n,\end{array}\right.\\ {s}_{{jk}}^{{mn}} & = & \left\{\begin{array}{ll}({g}_{x}^{({c}_{1})}{g}_{x}^{({c}_{2})})\left[{{\rm{\Gamma }}}^{(c)}({\omega }_{j})-{\left({{\rm{\Gamma }}}^{(c)}({\omega }_{k})\right)}^{* }\right]/2{\rm{i}} & {\rm{if}}\ m\ne n,\\ \displaystyle \sum _{\alpha ={l}_{m},{c}_{m}}{\left({g}_{x}^{(\alpha )}\right)}^{2}\left[{{\rm{\Gamma }}}^{(\alpha )}({\omega }_{j})-{\left({{\rm{\Gamma }}}^{(\alpha )}({\omega }_{k})\right)}^{* }\right]/2{\rm{i}} & {\rm{if}}\ m=n,\end{array}\right.\\ {\tilde{s}}_{{jk}}^{{mn}} & = & \left\{\begin{array}{ll}({g}_{x}^{({c}_{1})}{g}_{x}^{({c}_{2})})\left[{{\rm{\Gamma }}}^{(c)}(-{\omega }_{j})-{\left({{\rm{\Gamma }}}^{(c)}(-{\omega }_{k}\right)}^{* }\right]/2{\rm{i}} & {\rm{if}}\ m\ne n,\\ \displaystyle \sum _{\alpha ={l}_{m},{c}_{m}}{\left({g}_{x}^{(\alpha )}\right)}^{2}\left[{{\rm{\Gamma }}}^{(\alpha )}(-{\omega }_{j})-{\left({{\rm{\Gamma }}}^{(\alpha )}(-{\omega }_{k}\right)}^{* }\right]/2{\rm{i}} & {\rm{if}}\ m=n,\end{array}\right.\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{rcl}{r}_{{jk}}^{{mn}} & = & \left\{\begin{array}{ll}({g}_{z}^{({c}_{1})}{g}_{z}^{({c}_{2})})\left[{{\rm{\Gamma }}}^{(c)}({\omega }_{j})-{\left({{\rm{\Gamma }}}^{(c)}({\omega }_{k})\right)}^{* }\right]/2{\rm{i}} & {\rm{if}}\ m\ne n,\\ \displaystyle \sum _{\alpha ={l}_{m},{c}_{m}}{\left({g}_{z}^{(\alpha )}\right)}^{2}\left[{{\rm{\Gamma }}}^{(\alpha )}({\omega }_{j})-{\left({{\rm{\Gamma }}}^{(\alpha )}({\omega }_{k})\right)}^{* }\right]/2{\rm{i}} & {\rm{if}}\ m=n,\end{array}\right.\\ {u}_{j}^{{mn}} & = & \left\{\begin{array}{ll}({g}_{z}^{({c}_{1})}{g}_{z}^{({c}_{2})})\left[{{\rm{\Gamma }}}^{(c)}({\omega }_{j})-{\left({{\rm{\Gamma }}}^{(c)}({\omega }_{j})\right)}^{* }\right]/2{\rm{i}} & {\rm{if}}\ m\ne n,\\ \displaystyle \sum _{\alpha ={l}_{m},{c}_{m}}{\left({g}_{z}^{(\alpha )}\right)}^{2}\left[{{\rm{\Gamma }}}^{(\alpha )}({\omega }_{j})-{\left({{\rm{\Gamma }}}^{(\alpha )}({\omega }_{j})\right)}^{* }\right]/2{\rm{i}} & {\rm{if}}\ m=n.\end{array}\right.\end{array}\end{eqnarray} \tag{ B.2 }$$