Environmentally Induced Entanglement -- Anomalous Behavior in the Adiabatic Regime

Considering two non-interacting qubits in the context of open quantum systems, it is well known that their common environment may act as an entangling agent. In a perturbative regime the influence of the environment on the system dynamics can effectively be descried by a unitary and a dissipative contribution. For the two-spin Boson model with (sub-) Ohmic spectral density considered here, the particular unitary contribution (Lamb shift) easily explains the buildup of entanglement between the two qubits. Furthermore it has been argued that in the adiabatic limit, adding the so-called counterterm to the microscopic model compensates the unitary influence of the environment and, thus, inhibits the generation of entanglement. Investigating this assertion is one of the main objectives of the work presented here. Using the hierarchy of pure states (HOPS) method to numerically calculate the exact reduced dynamics, we find and explain that the degree of inhibition crucially depends on the parameter $s$ of low frequency power law behavior of the spectral density $J(\omega) \sim \omega^s e^{-\omega/\omega_c}$. Remarkably, we find that for resonant qubits, even in the adiabatic regime (arbitrarily large $\omega_c$), the entanglement dynamics is still influenced by an environmentally induced Hamiltonian interaction. We study the model in detail and present the exact entanglement dynamics for a wide range of coupling strengths, distinguish between resonant and detuned qubits, as well as Ohmic and deep sub-Ohmic environments. Notably, we find that in all cases the asymptotic entanglement does not vanish and conjecture a linear relation between the coupling strength and the asymptotic entanglement measured by means of concurrence. Further we discuss the suitability of various perturbative master equations for obtaining approximate entanglement dynamics.


I. INTRODUCTION
Entanglement -this very peculiar kind of correlation, not occurring in the classical world, is known to be a fragile property with respect to environmental influences [1][2][3][4]. On the other hand, such an environmental interaction features also the capability to induce entanglement [5][6][7][8], which can even last for arbitrary long times [9][10][11]. The investigation of these competing environmental effects [11][12][13][14][15][16][17][18][19][20][21][22] adds to the understanding of entanglement and decoherence in the context of open quantum systems which is not only of importance from a fundamental point of view but is also most relevant for fields like quantum computation, communication and metrology.
Albeit the vast amount of publications in this field, sophisticated results for the microscopic model of system + environment are rare even for the simple case of two qubits coupled to a common (sub-) Ohmic environment. Therefore we investigate the two-qubit entanglement by means of the numerically exact hierarchy of pure states (HOPS) method [23][24][25][26][27] with two major objectives in mind. First, we draw conclusions about the general applicability of various perturbative approaches. Favoring the Redfield equation (RFE) is in line with previous results [28]. As expected, our considerations confirm that in the perturbative regime the entanglement generation is well modeled by the environmentally induced Lambshift Hamiltonian mediating an effective qubit-qubit interaction. Further, it has also been argued that this environmental effect on the system may be canceled by the * richard.hartmann@tu-dresden.de so-called counterterm in the adiabatic limit [29]. Therefore, we secondly investigate the influences of the counterterm on the entanglement dynamics and find that it sensitively depends on the parameter s of the spectral density (SD). Whereas in the Ohmic case the inhibition of entanglement is well observed, in the deep sub-Ohmic regime (s < ∼ 0.3) a significantly larger cutoff frequency ω c is required to see the same effect.
The model under consideration extends the prominent spin-boson model [30] to two qubits Each qubit is modeled by the Pauli matrix σ x with tunneling frequency ω. They are coupled via σ z to a common bosonic environment where a (a † ) denotes the bosonic annihilation (creation) operator. Note, units are chosen such thath and k B become unity. For the continuous environment we choose a (sub-) Ohmic SD π λ g 2 λ δ(ω − ω λ ) = J(ω) = π 2 αω 1−s c ω s e −ω/ωc , (2) where α denotes the unitless coupling strength, ω c the cutoff frequency of the exponential cutoff and the power s allows to distinguish between the sub-Ohmic (s < 1), Ohmic (s = 1) and super-Ohmic (s > 1) regime.
For such a microscopic model time local master equations provide perturbative results in the weak coupling regime. In the case of two resonant qubits (ω A = ω B ), the widely used quantum-optical master equation (QOME) [31,32] of Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) form [33] reveals that the relaxation dynamics is accompanied by an environmentally induced Hermitian coupling between the two parties [28,32,[34][35][36]. This effective interaction easily explains the entanglement generation, even in the absence of a direct coupling [17,36,37]. However, for the more general case of two detuned qubits the QOME does not show entanglement generation [38,39]. This is consistent since the system dynamics experiences an additional timescale of the order of the detuning. As a consequence, for the QOME (the rotating wave approximation (RWA)) to be applicable, a substantially weaker coupling is required as compared to the resonant case (see [28]).
We show by means of numerical results that for a sufficiently large coupling strength two detuned qubits become entangled, too. We further investigate the question to what extent time local master equations, other than the QOME, approximate the exact entanglement dynamics of two detuned qubits sufficiently well.
In addition, we find from the exact dynamics that even in the weak coupling regime the entanglement of the asymptotic state does not vanish irrespectively of the detuning. The correct asymptotic value cannot be obtained by any of the time local master equations considered here.
To further enlighten the properties of environmentally induced entanglement we investigate the case where the environmentally mediated unitary qubit-qubit interaction is suppressed. A similar study based on GKSL type master equations has shown that if the bath induced unitary interaction is omitted, the action of the dissipator alone can result in entanglement generation as well [40].
By means of the full microscopic model it is, a priory, not clear how to implement this scenario since the exact form of the induced unitary interaction is not known. In the adiabatic limit, however, where the bath oscillators react instantaneously to the system dynamics [29], it takes the form which defines the so-called "counterterm" (Γ denotes the gamma function). Based on this reasoning, adding the counterterm H c to the microscopic Hamiltonian [Eq. (1)] is assumed to suppress the bath induced unitary interaction [19,41] and, thus, inhibit the buildup of entanglement. However, the adiabatic assumption should be questioned for the deep sub-Ohmic regime where the low frequency modes are especially important. We find and explain that for finite ω c including the counterterm does not fully suppress the entanglement generation. In addition we report the interesting phenomenon that for resonant qubits even in the limit ω c → ∞ the truly adiabatic regime is never reached. That is, a significant amount of entanglement will always be generated (and diminish afterwards) via an environmentally induced unitary interaction even in the presence of the renormalizing counter term.
The two-spin Boson model considered here with (sub-) Ohmic SD and exponential cutoff has also been investigated by means of other numerical methods. Unfortunately, the entanglement dynamics obtained using path integral Monte Carlo (PIMC) techniques [19] does not agree with our numerical and analytical results. As of missing information in Ref. [11] we were not able to connect our calculations to results obtained from the quasiadiabatic path integral (QUAPI) method. However, the time-evolving matrix product operator (TEMPO) algorithm [42], an advancement of QUAPI, gives results consistent with ours.

II. INDUCED ENTANGLEMENT DYNAMICS
As a measure of the two-qubit entanglement we use concurrence [43]. To examine its dynamics, the reduced dynamics of the two-spin boson model specified in Eq.
(1) is obtained numerically, solving the non-linear variant of the non-Markovian quantum state diffusion (MNQSD) equation [44,45] by means of the well tested HOPS approach [23][24][25]27]. To use this method for a (sub-) Ohmic environment the exact bath correlation function (BCF) with algebraic decay has to be approximated by a sum of exponentials where the accuracy of the approximation can be controlled by the number of exponential terms. Increasing the accuracy also increases the time to which the approximation follows the algebraic decay (see Fig.  1). Further details on how the HOPS method has been applied successfully on (sub-) Ohmic environments can be found in Ref. [25,27].
Throughout this work we refer to the expression c = λ 1 − λ 2 − λ 3 − λ 4 as concurrence [43,46] where, strictly speaking, only the positive values of c quantify entanglement. However, when comparing the dynamics of c it seems more convenient to also show its negative values. As usual, λ i are the decreasingly sorted eigenvalues of the matrix R = √ ρρ √ ρ withρ = σ A y σ B y ρ * σ A y σ B y . For numerical reason we calculate the eigenvalues λ i from the eigenvalues a i of the matrix ρρ via λ i = √ a i .

A. Exact Dynamics
To examine the properties of the entanglement generation we choose as initial condition the two-qubit product state |ψ 0 = |↑↑ and a zero temperature environment where |↑ (|↓ ) is the eigenvector of σ z with eigenvalue +1 (−1). Note, for the unbiased single qubit Hamiltonian ( σ z with = 0) considered here the initial state |↓↓ yields the same entanglement dynamics and is, thus, not considered explicitly. Also, the symmetry of the resonant Hamiltonian (1) (ω A = ω B ) results in a decoherence free subspace spanned by the Bell-state |Φ − ∼ |↑↓ − |↓↑ . Therefore, since the product states |↑↓ and |↓↑ contain a contribution from the decoherence free subspace, it is not too surprising that entanglement is generated [11] if they serve as initial condition. Consequently, they are not considered here but will be discussed elsewhere for detuned qubits.
As shown in Fig. 2, for the resonant case ω A = ω B a substantial amount of entanglement builds up initially (referred to the rescaled time tω A α(ω c /ω A ) 1−s ) for any weak to intermediate coupling strength as well as a sub-Ohmic and Ohmic environment. Subsequently, the entanglement diminishes. This kind of behavior can easily be explained by the interplay between the bath induced unitary interaction [28,36] and the relaxation of the system [4,47] (see also Sec. II B 1).
Our calculations show that for the slightly detuned scenario ω B /ω A = 0.95 the entanglement evolves very similar to the resonant case. As expected, differences are most prominent in the weak coupling regime. For a larger detuning ω B /ω A = 0.8 such differences become more significant.
Notably, our exact results reveal that in any case the asymptotic entanglement does not vanish (see Fig. 3). A result which cannot be obtained correctly by any of the perturbative approaches considered in the following. Further, the asymptotic value is rarely influenced by the detuning. This is particularly remarkable in the weak coupling regime where the initial dynamics is very sensitive to the detuning.
2 time tαω Ã α = 2.00 · 10 −2α = 6.32 · 10 −2α = 2.00 · 10 −1 Figure 2. The dynamics of the two-qubit entanglement quantified by concurrence is shown for a sub-Ohmic (s = 0.3) and Ohmic environment with ωc = 10ωA and various coupling strengths. By keepingα = α(ωc/ωA) 1−s constant the relaxation takes place on the same time scale for different parameters s. In addition different detuning parameters ωB/ωA are considered. The dynamics shown here was obtained by the HOPS method with an exponential approximation of the BCF which obeys an absolute error smaller than 10 −3 . The dynamics obtained using an even more accurate representation is indistinguishable with respect to the shown plots. In the same manner the convergence with respect to the hierarchy depth has bee checked. In this sense the dynamics shown here is exact.
.00 · 10 −2α = 6.32 · 10 −2α = 2.00 · 10 −1 Figure 3. The entanglement dynamics is shown with special emphasize on the asymptotic behavior (same parameters as in Fig. 2). The insets show the linear dependence of the asymptotic entanglement as a function of the coupling strength (double logarithmic plot).

B. Perturbative Master Equations
In addition to the exact calculations the properties of the dynamics obtained by perturbative master equations are of relevance, too. The motivation is given by the rather simple structure of such master equations which allows to gain insights in the relevant mechanisms. However, the suitability of these perturbative approaches has to be checked since an error estimation from within the perturbative formalism seems unfeasible. With the help of the exact dynamics obtained by the HOPS method we examine the applicability of the QOME [32,34], its variation with only a partial rotating wave approximation (PRWA) [48][49][50], the Redfield equation (RFE) [51], the very recent Markovian arithmetic geometric mean approximation (MAGMA) [52] and the coarse-graining approach [38,39,53] in the context of two independent qubits coupled to a common environment with sub-Ohmic SD.

The Rotating Wave Approximation
The great success of the QOME (Born-Markov approximation and full RWA, well known GKSL form, see Ref. [28,32,34]) to describe the dynamics of a single qubit encourages the use of the same formalism for two qubits also. Since the QOME is of GKSL form positivity of the reduced state is assured for all times and any initial condition. For the resonant case ω A = ω B it is easily seen that the unitary contribution of the QOME (Lamb shift term) [9,28,36] (4) couples the qubits such that entanglement builds up [11,17,36,37] (see dashed graphs in Fig. 2). Here S denotes the imaginary part of the half sided Fourier transform of the zero temperature BCF α bcf which is defined by means of the SD J(ω) Since the local contribution H local acts on a single qubit only it does not influence the entanglement dynamics. Importantly, as of the RWA involved in the derivation of the QOME, for detuned qubits ω A = ω B the Lamb shift term consists only of the local part and, henceforth, entanglement generation is not featured by the QOME [38]. However, as shown in Fig. 4 the exact dynamics obtained from the HOPS method reveals that for slightly detuned qubits entanglement is generated in a very similar manner as in the resonant case (see also Ref. [22,28,54,55] for more detailed discussions on  Note, the QOME and the PRWA coincide in the resonant case. Whereas the QOME is not able to predict the entanglement generation for detuned qubits, the PRWA approach agrees well with the exact dynamics (except for the superimposed oscillations). Both master equations yield a vanishing asymptotic entanglement whereas the exact value remains finite.
the severe consequences of the RWA on bipartite correlations). Assuming that the QOME provides a suitable approximation in some weak coupling regime it seems contradictory that an infinitesimal change of the system parameters results in a significant change of the dynamics. The inconsistency can be resolved by recalling that a small detuning introduces a new very slow system time scale. Only for a decay time scale (inverse of the coupling strength) much larger than the system time scale the RWA and, thus, the QOME is applicable. This means that the applicability of the QOME crucially depends on the system parameters. Only in the limit of zero coupling (in combination with a rescaled time, the so called scaling-limit) the QOME becomes exact [56]. Therefore, for a sufficiently weak coupling strength the entanglement generation will vanish for any detuning, but remain in the resonant case. This can, to some extent, be seen from the exact dynamics shown in Fig. 2.
After having pointed out the severe problems in the detuned case arising from the RWA we consider perturbative approaches which aim to circumvent the RWA. A first one, very similar to the QOME, makes use of the RWA only partially [48][49][50] which results in a master equation of GKSL kind also. The specific expression of the master equation for the model considered here can be found in Ref. [28]. The master equation using the PRWA is constructed from the same non-local Lindblad operators occurring in the QOME for resonant qubits. Henceforth, the Lamb shift term is non-local and features entanglement generation also in the detuned case.
The initial entanglement dynamics obtained from the PRWA master equation agrees quantitatively, up to some fast oscillation, with the exact dynamics obtained by the HOPS method (see Fig. 4). The fast superimposed oscillation are not captured because the PRWA still neglects some secular terms. Keeping all secular terms yields the RFE considered next. Both, the QOME and the PRWA approach yield a vanishing asymptotic entanglement which is in disagreement with the exact dynamics.

The Redfield Equation
Motivated by the detailed accuracy assessment of perturbative master equations [28] we consider not only all secular terms but also use RFE with time dependent coefficients which was shown to be the most accurate method. As shown in Fig. 5 the high degree of accuracy also holds for the sub-Ohmic environment considered here. The entanglement dynamics obtained from the RFE with time dependent coefficients matches very well the exact dynamics. Even the fast oscillations are recovered.
Small positivity violations of the reduced dynamics occur already in a regime where the accuracy is still acceptable (see Fig. 6). Consequently, using the positivity violation as an indicator for the break down of the weak  Figure 6. To justify the use of the positive matrix ρρ † in order to calculate the concurrence for the RFE the deviation of ρ and ρρ † from the exact state is shown using the Hilbert-Schmidt-norm. The deviations are nearly indistinguishable which means that ρρ † approximates the exact state as good as ρ. This is also reflected by the magnitude of the smallest negative eigenvalue of ρ which is significantly smaller than the deviation. In addition the absolute difference between the exact and the RFE concurrence is shown. The error of the concurrence is roughly of the same order of magnitude as the error of the state. The error of the concurrence as well as the density matrix is also plotted for the MAGMA which shows that ensuring positivity by the structure of the master equation results in a loss of accuracy.
coupling assumption, as proposed for a Lorentzian environment with an exponentially decaying BCF [28], is not that simple for a sub-Ohmic SD with a slow algebraic decay of the BCF. The positivity violation is of particular relevance here because the concurrence is not defined for non-positive ρρ even though ρ approximates the exact state suitable well. However, since the positivity violation is of the order of the perturbative error the matrix ρρ † is consistent with the degree of perturbation (see Fig. 6) while being positive by construction. This justifies the calculation of the concurrence using the positive matrix ρρ † obtained from the reduced dynamics of the RFE.

The Markovian Arithmetic Geometric Mean Approximation (MAGMA)
In a recent publication the failure to quantify entanglement of a non-positive state in an approximative sense has been addressed as well [52]. The proposed MAGMA modifies the RFE such that it becomes a master equation of GKSL type. Crucially, the environmentally induced unitary influence on the system is identified from the RFE before the approximation is applied. The justification of MAGMA is based on the relaxation time scale in the interaction picture, roughly given by the inverse of the coupling strength, and properties of the SD only. In contrast to the QOME the particular spectrum of the system Hamiltonian is irrelevant. It is, thus, expected that the entanglement generation is well captured for any detuning of the qubits. As shown in Fig. 5 this is indeed the case. The MAGMA master equation mimics even the fast oscillations of the entanglement dynamics, however, slightly less accurate compared to the positive state constructed from the RFE dynamics (see also Fig. 6 which shows the error). Since the MAGMA master equation is based on the RFE it is consistent that the asymptotic entanglement deviates from the exact value on the same scale as the RFE result (see Fig. 5), although the value is overestimated. This allows to conclude that the positivity issues related to the evaluation of the entanglement can be cured slightly more accurate using the Redfield formalism with ρ := ρρ † than applying the additional MAGMA.

The Coarse-Graining Master Equations
An alternative path towards a master equation for the microscopic model which assures positive dynamics is the so-called coarse-graining procedure [38,39,53]. It has been proposed precisely with the aim to overcome the limitations due to the RWA. In particular for the case of two detuned qubits entanglement generation has been shown using such a coarse-graining master equation (CGME) [38]. We confirm this result for a common sub-Ohmic environment (see Fig. 7). However, the quantitative comparison for a sub-Ohmic environment with s = 0.3 and ω c = 10ω A shows that the CGME does not reach the accuracy of the PRWA. To understand the deviation of the CGME recall that the condition for its applicability reads τ env τ τ ind. [28]. The environmentally induced timescale for the system dynamics scales with the inverse coupling strength. From the exact dynamics in Fig. 7 its follows τ ind ≈ 0.1(αω A ) −1 ≈ 60ω −1 A . Although the decay time of the BCF scales with the cutoff frequencies τ env ∼ ω −1 c the particular time at with the BCF has decayed over, for example, 2 orders of magnitude is τ env,2 ≈ 10 2/(s+1) ω −1 c . For the example considered here it follows τ env,2 ≈ 3.5ω −1 A . Obviously, a clear separation of the three timescales is not justified.

Summary
As expected, the approach with the least approximations, namely the RFE with time dependent coefficients, yields the most accurate results for the entanglement dy- τ ω A = 10 Figure 7. The entanglement dynamics of the CGME is shown for comparison (same environmental parameters as in Fig. 4).
Since the QOME is justified in the resonant case it is not surprising that for large coarse-graining parameters τ the CGME yields good results too. For the detuned case the entanglement generation is still visible, in contrast to the QOME. However, the PRWA agrees better with the exact results than the CGME (see the text for a discussion about the time scale separation required by the CGME).
namics of two qubits coupled to a common sub-Ohmic environment. Positivity violations can consistently be healed by using the positive matrix ρρ † to estimate the concurrence. The PRWA, which coincides with the QOME in the resonant case, yields the correct overall dynamics while missing faster oscillations. The CGME is less appropriate for the set of parameters considered here. Note that the recent MAGMA master equation outperforms the other GKSL-type master equations in terms of accuracy. However, as of the simple structure of the PRWA, we use this approach together with the exact results to illuminate the influence of the counterterm in the next section. The reasoning could equally well have been done based on the MAGMA master equation since due to the weak coupling both approaches yield nearly indistinguishable dynamics.

III. INFLUENCE OF THE COUNTERTERM
To recapitulate the motivation for the countertermcompensation of the environmentally induced unitary effect on the system dynamics -a particle with position q and momentum p moving in a potential V (q) and coupled to a set of harmonic oscillators with position (momentum) x λ (p λ ) is considered. The linear coupling − λ F λ (q)x λ to the environmental modes renormalizes the potential affecting the particle. It has been argued that in the adiabatic regime (instantaneous adjustment of the environmental modes to the particle position) the For a mode independent coupling operator F λ (q) = c λ F (q), as in Eq. (1), the Hamiltonian compensates the change of the potential. It is assumed that the same holds true when casting the above Hamiltonian to the more abstract form of Eq. (1). In particular the relations g λ √ 2m λ ω λ = c λ and L = −F (q) yield for the counterterm (8) For the class of (sub-) Ohmic SD, as given in Eq. (1), the counterterm evaluates to H c = αωc 2 Γ(s)L 2 . Note that for two non-interacting qubits it is primarily the environmentally induced unitary interaction which generates entanglement. Therefore, including the counterterm is expected to suppress the generation of entanglement. Its influence on the entanglement dynamics is shown in Fig. 8 (resonant qubits) and Fig. 9 (detuned qubits). Two main features are observed. First, with decreasing s (going from the Ohmic to the deep sub-Ohmic regime) the impact of the counterterm becomes less, independent of the coupling strength (see first and second row in Fig. 8). Second, for weak coupling and s = 1 the entanglement dynamics including the counterterm is very well mimicked by the action of the dissipator only of the PRWA master equation (which becomes the QOME in the resonant case). This shows that the counterterm does compensate the induced energy shift up to a certain level. However, for s = 0.3 the exact dynamics, hardly influenced by H c , does not agree with the dissipative dynamics in the perturbative regime. Consequently the above assertion about the effect of the counterterm does not hold true in the deep sub-Ohmic regime.
To understand this behavior we investigate the dynamics in the weak coupling regime. In that regime the RFE is known to be very accurate [28]. Using asymptotic coefficients it readṡ diss. only (PRWA) Figure 9. For slightly detuned qubits (ωB/ωA = 0.95) and a cutoff frequency ωc/ωA = 10 the entanglement dynamics under the influence of the counterterm is shown (same coloring as in Fig. 8). As of the detuning the PRWA master equation is used to obtain the perturbative dynamics. Similar to the case of resonant qubits, the counterterm approximately cancels the Lamb shift term in the Ohmic case. As shown in the lower left panel, the exact dynamics including the counter term (orange line) agrees with the dynamics under the action of the PRWA dissipator only (dashed black line). However, this expected behavior is not observed for a sub-Ohmic SD with small s. Including the counterterm rarely affects the entanglement dynamics (blue and orange lines, lower right panel). The influence of the cutoff frequency is shown in Fig.  11.
position of the coupling operator by means of projectors of eigenstates of the system Hamiltonian H sys | = | . It is well known that in the GKSL-formalism the real part of F (ω) = ∞ 0 dτ α bcf (τ )e iωτ = J(w) + iS(ω), which is the SD, accounts for dissipative dynamics whereas the imaginary part determines the unitary contribution. Under the assumption that the particular values S(ω i ) can be approximated by S(0) (recall ω i takes the value of all possible transition frequencies of the system Hamiltonian) the corresponding contribution in the RFE becomes where we have used that by definition L = i L ωi . Note, the proposition that S(ω i ) ≈ S(0) does not imply that J(ω i ) ≈ J(0). For zero temperature S(ω) is related to the SD by where P denotes the Cauchy principal value. We finally conclude that given S(ω i ) ≈ S(0) the unitary contribution in the perturbative treatment (Lamb shift term) is approximately canceled by the counterterm For the class of (sub-) Ohmic SD the validity of the assumption S(ω i ) ≈ S(0) depends on the parameter s in a very sensitive way. To see that, it is convenient to expand the function S(ω) by introducing x = ω/ω c . Expressing F (ω) analytically and expanding the incomplete gamma function for small x yields where e γ denotes the Euler-Mascheroni constant and Θ the Heaviside step function. The expansion shows that the lowest order behaves like |x| s in the sub-Ohmic and x log(|x|) in the Ohmic case (see Fig. 10 for examples). Thus, the convergence rate for S(x) to approach S(0) decreases with s. In addition, for positive x the term ∼ |x| s changes its sign at s = 0.5 which results in a pointed minimum at S(0) for s < 0.5.
This behavior explains that the cancellation of the Lamb shift by the counterterm fails spectacularly for small s. The dynamics of these Hermitian contributions are also shown in Fig. 8 (bottom row). As expected, the pure influence of the Lamb shift on the system Hamiltonian results in oscillatory dynamics of the concurrence (green graph). Adding the counterterm H c (purple graph) yields a significantly slower generation of entanglement for s = 1. The two contributions compensate, at least approximately, on the relaxation time scale. However, for s = 0.3 the effect of including H c is far from cancellation -entanglement is generated even slightly faster. This behavior is directly captured by the exact entanglement dynamics of the two-spin-boson model calculated using the HOPS method ( Fig. 8 and Fig. 9).
The above considerations suggest that in the limit ω c /ω A → ∞, where the approximation S(ω i ) ≈ S(0) becomes exact, the counterterm should truly cancel the Lamb shift contribution. This behavior is confirmed by the example shown in Fig. 11 where the entanglement Remarkably, for resonant qubits (right column) this expectation is not fulfilled. In that case the joint unitary dynamics Hsys + H Lamb + Hc (black solid line) builds up entanglement on the same time scale as the dissipation takes place. Consequently, the exact dynamics does not approach the dynamics of the QOME under the action of the dissipator only (dashed line). This qualitative difference originates from the time scale set by the detuning which is independent of ωc (gray vertical lines). Note, this time scale decreases on the rescaled timẽ αtωA while increasing ωc. Since the entanglement dynamics of the joint unitary part is periodic on that time scale, but the buildup takes place on the time scale set by ∆S ∝ S0 (constant on the rescaled time), the generation of entanglement is effectively suppressed. dynamics of two detuned qubits (left column) approaches the purely dissipative dynamics when increasing the cutoff frequency.
Remarkably, this does not hold true for resonant qubits (see right column in Fig. 11). To affirm that this effect remains even in the limit ω c /ω A → ∞ we argue as follows. The remaining Hermitian contribution H Lamb + H c scales in lowest order like which sets the timescale on which the unitary interaction induced by the Hamiltonian [Eq. (7)] takes place. Obviously, increasing ω c increases that timescale which could mean an effective cancellation if the damping takes place on a faster time scale. It turns out, however, that this is not the case since the damping rate γ, determined by the SD at the resonance frequency, scales with ω c in the same manner. This can be seen by expressing the coupling strength in terms of S 0 and expanding the SD in In summary we find the counter-intuitive result that for two resonant qubits and a sub-Ohmic environment the effect of the Lamb shift contribution is not compensated by the counterterm, even in the limit ω c → ∞ which is usually considered the adiabatic limit. Consequently, the entanglement generation beyond the purely dissipative contribution remains for any ω c . This, in particular, stays in contrast with the results presented in Ref. [19]. The authors claim that for sufficiently weak coupling the model including the counterterm (as in Eq. 7) does not result in any entanglement generation for the two qubits.

IV. CONCLUSIONS
In the context of an open two-qubit system we address the fundamental question: To what extent does the incorporation of the counterterm into the microscopic model counterbalance the effect of the Lamb shift Hamiltonian and, thus, inhibit the environmentally induced generation of entanglement? The answer requires some prerequisites which we provided and discussed as well. First of all we used the well tested HOPS method to obtain the exact entanglement dynamics. With the exact dynamics at hand we were able to assess the accuracy of various perturbative approaches. We found that the RFE yields the most accurate results which is not too surprising since it involves the least approximations. We argued that positivity violations can be consistently dealt with by using the positive matrix ρρ † to estimate the amount of entanglement. However, to answer our main question a master equation of GKSL kind is desirable which clearly separates the unitary from the dissipative influence of the environment. As expected, we confirm that the master equation based on the MAGMA as well as PRWA yield acceptable entanglement dynamics also for detuned qubits. Concerning the above assertion, if the Lamb shift Hamiltonian (unitary influence of the environment) is truly compensated by the counterterm the resulting dynamics should match the perturbative dynamics obtained by the master equation under the exclusive action of the dissipator in the weak coupling regime. Therefore, we compared the exact dynamics including the counter term with the dissipator only dynamics of the PRWA approach. The MAGMA approach could have been used equally well. It is clear, that for finite ω c , if at all, the compensation is only approximate. We found that for an Ohmic environment the exact dynamics matches fairly well the dissipative dynamics of the perturbative approach -the counterterm compensates the Lamb shift Hamiltonian in some sense. However, decreasing s (sub-Ohmic environment) while keeping ω c constant worsens the compensation. For small s < ∼ 0.3 the entanglement generation is rarely influenced by the counterterm. This behavior can be traced back to the power law behavior of S(ω) ∼ |ω/ω c | s , the imaginary part of the half sided Fourier transform of the BCF, and is amplified by the pointed minimum of S(ω) at ω = 0 for s < 0.5. Remarkably, for resonant qubits and a sub-Ohmic environment even in the adiabatic limit ω c → ∞ the entanglement generation cannot be modeled solely by the dissipator of the perturbative master equation. We showed that on the one hand, increasing ω c slows down the generation of entanglement, which is due to the still imperfect compensation. On the other hand the relaxation time increases in precisely the same manner. Consequently, even in the limit ω c → ∞ the entanglement dynamics is influenced by the infinitesimal non-zero difference between the Lamb shift Hamiltonian and the counterterm.
In summary, for an environment with finite ω c including the counterterm in the microscopic model does not imply that the unitary influence of the environment on the system is canceled. For the general detuned case, in order to achieve the same level of compensation while going from the Ohmic to the deep sub-Ohmic regime the cutoff frequency ω c needs to be increased. It is a special property of two resonant qubits that even tough the Lamb shift Hamiltonian and the counterterm approach each other in the large cutoff regime the entanglement dynamics is influenced significantly by the remaining infinitesimal difference.
Alongside these main findings we reported that in a perturbative regime the RFE with time dependent coefficients yields the most accurate entanglement dynamics when curing positivity issues by the replacement ρ → ρρ † . The recently proposed MAGMA which guarantees positivity yields results which qualitatively agree with the RFE while being less accurate. Applying the RWA only partially yields reasonable entanglement behavior (superimposed fast oscillations are obviously missing) also for detuned qubits. As of the rather small cutoff frequency considered here and the slow algebraic decay of the BCF the CGME is not suitable. These findings are in line with results in Ref. [28] obtained for a Lorentzian environment.
In addition we report that although the initial entanglement generation is followed by a decay, the asymptotic value of entanglement does not vanish for the initial states |↑↑ and |↓↓ . From the exact numeric solution we find a linear relation between the asymptotic entanglement and the coupling strength. None of the perturbative methods considered here capture this feature correctly.
Further, independent of the particular model and based on the RFE we argued that the Lamb shift Hamiltonian and the countertem cancel up to some error which is related to approximation S(ω i ) ≈ S(0) for all ω i where ω i denote all possible transition frequencies of the system Hamiltonian. For a (sub-) Ohmic environment this approximation becomes exact in the adiabatic limit ω c → ∞. This, however, does not necessarily imply that the dynamics is not influenced by the remaining very small unitary influence of the environment, as explained above.