Exact dynamics of interacting qubits in a thermal environment: Results beyond the weak coupling limit

We demonstrate an exact mapping of a class of models of two interacting qubits in thermal reservoirs to two separate spin-bath problems. Based on this mapping, exact numerical simulations of the qubits dynamics can be performed, beyond the weak system-bath coupling limit. Given the time evolution of the system, we study, in a numerically exact way, the dynamics of entanglement between pair of qubits immersed in boson thermal baths, showing a rich phenomenology, including an intermediate oscillatory behavior, the entanglement sudden birth, sudden death, and revival. We find that stationary entanglement develops between the qubits due to their coupling to a thermal environment, unlike the isolated qubits case in which the entanglement oscillates. We also show that the occurrence of entanglement sudden death in this model depends on the portion of the zero and double excitation states in the subsystem initial state. In the long-time limit, analytic expressions are presented at weak system-bath coupling, for a range of relevant qubit parameters.


I. INTRODUCTION
Understanding the dynamics of a dissipative quantum system is a prominent challenge in physics, as a quantum system is never perfectly isolated from a larger environment. A minimal, yet highly rich model for exploring quantum dissipation effects, is the spin-boson model, including an impurity two-level system (referred to as a spin) coupled to a thermal reservoir. This model displays a rich phase diagram in the equilibrium regime [1,2]. The non-equilibrium version of this model, referring to the case where the spin is coupled to two thermal reservoirs, has been suggested as a toy model for exploring quantum transport phenomenology through an anharmonic nanojunction [3,4]. In this case, the generic situation is one of a non-equilibrium steady-state, regardless of the initial preparation.
Interacting two-level systems are the basic element in quantum computation, thus it is paramount to extent the minimal spin-boson scenario and describe more complex modular systems, e.g., two interacting qubits immersed in a thermal environment [5]. For a schematic representation, see Fig. 1. The qubits may share their thermal environment, or may separately couple to independent baths, maintained at a nonzero temperature. The latter situation corresponds to the case where the qubits are not necessarily placed close to each other. In another relevant setup, one qubit couples indirectly to a thermal reservoir, through its interaction with the other qubit. This situation effectively corresponds to a subsystem anharmonically coupled to a harmonic bath, allowing to introduce nontrivial nonlinear effects [6]. Physical realizations include, for example, ultracold atoms in optical lattices [7], trapped ions [8], resonator-coupled superconducting qubit arrays [9,10], and electron spins in quantum dots and doped semiconductors [11]. In such systems one should consider (at least) four energy scales: the internal qubit energetics, controlling its isolated (Rabi oscillation) dynamics, qubit-bath interaction strength and the environment temperature, leading to decoherence and relaxation processes, and qubit-qubit coupling energy, admitting state transfer between qubits and a nontrivial gate-functionality.
The dissipative multi-qubit system has recently served as a simple model for resolving issues related to coherence dynamics in the time evolution of biological molecules, e.g., the Fenna-Matthews-Olson (FMO) complexes, resulting in the identification of the relevant decoherence, relaxation and disentanglement timescales [12][13][14]. It should be noted that these works have considered single-excitation states only, ignoring the contribution of the zero and the doubly excited states in the two-qubits dynamics.
In this paper, we analytically demonstrate that a class of interacting two-qubit systems immersed in separate thermal reservoirs or within a common bath, can be mapped onto two uncoupled spin-bath problems, allowing for an exact numerical solution of the qubits dynamics. For a bosonic environment and a particular system-bath interaction form, we perform those simulations using an exact numerical technique, the quasi adiabatic path-integral (QUAPI) approach [15], providing the population and coherence dynamics of the system. With this at hand, we can follow the exact dynamics of entanglement between the qubits, as quantified by Wootters' concurrence [16]. For a range of system and bath parameters, analytic results are presented, describing the system behavior in the long time limit. The model investigated here is more general than what has been typically considered before, going beyond the simple exchange interaction model [13,14].
Entanglement is associated with nonclassical correlations between two or more quantum systems [17]. Since it is a basic resource in quantum computation and information technology, it is important to understand the extent to which environmental-induced decoherence processes degrade and destroy it [18], or alternatively, generate [19] and maintain it [20]. It has been recently shown that two qubits in separate reservoirs may disentangle at finite times, as opposed to the behavior of coherences. This process is referred to as "entanglement sudden death" [18,21]. More recent theoretical and experimental studies have looked at related effects, e.g., the collapse and subsequent revival of the entanglement [22], or its delayed-sudden birth, induced by a dissipative bath [23]. Steady-state entanglement generation by dissipation has been recently observed in atomic ensembles [24]. These studies have assumed non-interacting qubits, and the system dynamics has been typically followed within quantum master equation approaches (e.g., the Redfield equation or Lindblad formalism [25]), by invoking the weak system-bath coupling approximation. The markovian limit has been further assumed in many cases, see e.g., [19,26]. Our work here departs from these studies in two substantial aspects. First, we consider a more complex model for the subsystem, introducing qubit-qubit interaction, with the motivation to examine a setup relevant for quantum computing technology. Using an exact mapping, we show how the dynamics can be followed within a simpler construction: While we take into account the zero excitation and double excitation states, under certain initial conditions their dynamics can be separated from the evolution of the single-excitation states. However, their contribution to the pair entanglement is paramount. Second, we refrain from making approximations: We study the qubits dynamics using a numerically exact method, assuming a class of initial states. The results are valid beyond the weak system-bath coupling scenario, accommodating non-markovian effects.
Our calculations display rich dynamics. Particularly, we observe the development of a stationary concurrence due to the coupling of the qubits to thermal baths. This result stands in a direct contrast to the oscillatory behavior observed in the fully coherent regime. It demonstrates that while entanglement inherently relies on the existence of quantum correlations in the system, it nevertheless requires non-vanishing decoherence and relaxation effects in order to be stabilized and become useful for quantum technologies. Other phenomena detected and explained here are entanglement delayed sudden birth, sudden death, and revival.
The paper is organized as follows. In Sec. II we describe the model of interest, and explain its mapping onto two separate spin-bath problems. In Sec. III we explain how we follow the system dynamics, and include relevant expressions for calculating the qubits concurrence. Numerical results within QUAPI are included in Sec. IV. The long-time limit is discussed in Sec. V. Sec. VI concludes.

II. MODEL
The general model to be considered here includes two interacting qubits, i = 1, 2, immersed within separate reservoirs, L and R, respectively. The formalism can be reduced to describe a single-bath scenario. The total Hamiltonian includes three terms, H S and H B = H L + H R stand for the system and reservoirs Hamiltonians, respectively. The former includes the isolated qubits with the internal energy bias ǫ i and a qubit-qubit interaction term V ss , σ p i (p = x, y, z) are the Pauli matrices for the ith spin, J is an energy parameter characterizing the exchange interaction, γ and δ set the interaction anisotropy. Our mapping holds for a dephasing-type system-bath interaction model, Here, B ν is a ν = L, R bath operator, with B L coupled to spin '1' and B R coupled to spin '2'. In our simulations below we adopt bosonic reservoirs: each thermal baths includes a collection of independent harmonic oscillators, The operators a † k and a k are bosonic creation and annihilation operators, respectively, ω k is the mode frequency. We also assume that the interaction operators constitute the reservoirs displacements from equilibrium, Here, λ k are system-bath coupling constants. The mapping described next, from the Hamiltonian (1) into two spin-bath problems, neither rely on a particular bath statistics nor on the details of the operators B ν . For example, it is valid for a model of two-qubits in fermionic or spin environments. The Hamiltonian introduced so far takes into account two independent reservoirs. We could explore a similar setup with one qubit coupled to a thermal reservoir indirectly, through its interaction with the second qubit, mimicking nonlinear effects [6], see Fig. 1. Another relevant setup includes two qubits immersed in a common thermal reservoir. The Hilbert space of the qubits is spanned by four vectors, |00 , |01 , |10 , |11 , forming the excitation basis of the Hamiltonian, where the left (right) digit indicates the state of qubit 1 (2). We now show that the Hamiltonian can be mapped onto two spin-bath type models. We begin by defining four composite system operators, In the excitation basis, these operators take the explicit form Additionally, we construct two identity-type operators, With these at hand, the Hamiltonian (1) can be written as where Here, One can easily show that the following commutators vanish [27] (m, n = x, z) The Hilbert space of two qubits can thus be factored into two direct-sum subspaces: The first, P , is spanned by |01 and |10 . The second, Q, is spanned by |00 and|11 . One can further prove that P x , [P x , P z ] and P z generate an SU P (2) group, and Q x , [Q x , Q z ] and Q z generate another SU Q (2) group. The two groups have a direct-sum structure, SU P (2)⊕ SU Q (2). We now note that P x and P z in subspace P play the role of the Pauli matrices σ x and σ z (in the space spanned by |0 and |1 ). The same principle holds for Q x and Q z in the Q subspace. Overall, in mapping Eq. (1) into Eq. (10) we replaced a model of two interacting spins coupled each to its own thermal reservoir, by a model of two separate spin-bosontype systems, where each spin in the new model is coupled to both reservoirs. The latter model is significantly simpler than the former, and we can explore its dynamics using exact simulation tools for a class of certain initial conditions. The mapping described here holds for general reservoirs and a bilinear system-bath interaction form, with an arbitrary bath operator coupled to the subsystem. The results also hold when we apply a dressing transformation [28] W , H ′ = W † HW on the two qubits Hamiltonian, or the reservoirs. For example, we may introduce a spin-orbital coupling into V ss via W = exp(i θ 2 σ z 1 ), while keeping other terms in the total Hamiltonian unchanged Another relevant case that can be simulated exactly relies on the absence of external fields, ǫ 1 = ǫ 2 = 0, taking W = exp(i π 4 (σ y 1 + σ y 2 )). This results in After this transformation, the model has turned into the anisotropic XYZ-type model with flip-flop (σ x ) coupling between the system and reservoirs. In this form, the model describes energy exchange between the qubits and the baths, unlike the original Hamiltonian [Eq. (3)] which delineates dephasing effects. When

III. DYNAMICS AND QUANTUM ENTANGLEMENT
We explain here how we time-evolve the reduced density matrix, to obtain the qubits dynamics. As an initial condition we consider a system-bath product state, where the reservoirs are maintained in a canonical-thermal state, ρ ν = e −βν Hν /Z ν , Z ν = Tr B [e −βν Hν ] is the partition function. In order to separate the dynamics into the Q and P branches, we must adopt a direct sum Q-P initial state for the qubits. Overall, the total density matrix at time t = 0 is written as For a particular example, see Eq. (23). Under this construction, the reduced density matrix follows The trace is performed over the L and R degrees of freedom. The time evolution operators, with H P and H Q given in Eq. (10). Eq. (15) establishes an important result: The dynamics of the two-qubit system proceeds in two independent branches, each equivalent to a spin-bath model. In the case of bosonic baths, the dynamics in each branch is followed next using the QUAPI technique [15], a numerically exact simulation tool that can be easily extended to include more than one thermal reservoir. The output of this calculation is the reduced density matrix of the qubits. We use this information and investigate the time evolution of entanglement between the qubits. As a side comment, we note that if the qubits were to couple to a common bath,B = 0, only the Q subspace would have become susceptible to decoherring effects, while the P subspace would be an invariant subspace, or a "decoherence free" subspace. Based on Eqs. (14)- (16), we conclude that the reduced density matrix ρ S is an X-type matrix at all times, once we organize it in the standard order of basis vectors |00 , |01 , |10 and |11 , This is an important result: The dynamics of a class of dissipative interacting qubits [Eqs.
(1)- (3)] can be reached via the solution of two spin-bath problems, and the reduced density matrix satisfies an X-form at all times. Different quantities are of interest, e.g., the timescale for maintaining coherences in the system [14]. Here, we focus on quantum correlations in the system [17], computed next in a numerically exact way beyond weak coupling. In particular, we quantify the degree of entanglement between the qubits using Wootters' concurrence [16]. For mixed states it is calculated by considering the eigenvalues of the matrix r(t) = ρ S (t)σ y 1 ⊗ σ y 2 ρ * S (t)σ y 1 ⊗ σ y 2 , given here by In terms of these eigenvalues, the concurrence is defined as It varies from C = 0 for a disentangled state to C = 1 for a maximally entangled state. In the present case it reduces to with The dynamics of concurrence for an X-state density matrix has been examined in different works. For example, in Ref. [26] it is demonstrated that the effect of entanglement sudden death should always take place in a noninteracting qubit system once coupled to a finite temperature reservoir. As we mention in the introductory section, we depart from this study and similar works in two aspects: (i) We build the reduced density matrix using an exact numerical treatment, and (ii) we consider a more general model, including quit-qubit interaction effects, with the motivation to consider a setup more relevant for quantum computation technologies.

IV. NUMERICAL RESULTS
We simulate the spin-boson dynamics in the P and Q branches (separately) using QUAPI [15], to obtain the population and coherences in each branch. With this at hand, we generate the 4 × 4 reduced density matrix ρ S (t), Eq. (17). The qubits degree of entanglement is calculated using Eq. (20). Our general description assumes two thermal reservoirs: H L , coupled to spin 1 and H R , coupled to spin 2. These reservoirs are characterized by the spectral function J ν (ω) = π k∈ν λ 2 k δ(ǫ k − ω). Specifically, we simulate Ohmic baths, J ν (ω) = πKν 2 ωe −ω/ωc ; ω c is the cutoff frequency. The dimensionless prefactor K ν is referred to as the Kondo parameter, describing the strength of the systembath interaction energy. In practice, our simulations were performed without the R reservoir, by taking K R = 0. The reason for this choice is that in the spin-boson Hamiltonian (10) the inclusion of identical reservoirs which are interacting in the same manner with the spins (same functional form for B L and B R , up to a sign), simply amounts to an additive operation, reflecting a linear scaling of the Kondo parameter when more than one reservoir is incorporated. In what follows, we thus use the short notation K ≡ K L , K R = 0 and T ≡ T L .
The energy parameters in the system are the qubit-qubit coupling, taken as J = 1, and the anisotropy parameters δ = 0.1, γ = 0.5. The qubits are assumed identical with ǫ 1 = ǫ 2 ∼ 0.1 − 0.5. For the reservoir we take as a cutoff frequency ω c = 7.5, and use temperatures at the range T = 0.1 − 1. The Kondo parameter extends from the weak coupling limit (K = 0.05) to the strong-intermediate regime, K = 0.8, where convergence of QUAPI can be achieved. The following initial condition is utilized for the qubits subsystem with 0 ≤ a ≤ 1 and ρ Q,P (0) as (maximally entangled) Bell states, The concurrence (20) can be simplified if the following conditions are simultaneously satisfied: (i) a ≤ 1/2, and (ii) ǫ 1 = ǫ 2 . The latter condition, combined with the initial state ascribing identical weight to diagonal elements in the P subspace, implies that the populations in the P (single excitation) subspace are identical at all times, (ρ S ) 01,01 = (ρ S ) 10,10 = 1−a 2 . Since (ρ S ) 00,00 (ρ S ) 11,11 ≤ a 2 at all times and the density matrix positivity condition demands that |ρ i,j | 2 ≤ ρ i,i ρ j,j , the off diagonal terms are bounded by |(ρ S ) 11,00 | ≤ a 2 . This implies that F 2 cannot be larger than zero at any instant for a ≤ 1/2. The concurrence can then be simplified to C a≤ 1 2 (t) = max(0, 2F 1 ).
This expression indicates that C is nonzero if the magnitude of the coherence in the P subspace is large, in comparison to the product of populations in the Q subspace. Thus, to understand the behavior of entanglement between the qubits one needs to follow the population and coherence dynamics in both subspaces. At the special point a = 0 the concurrence reduces to the simple form which only depends on coherences behavior, a continuous function. As a result, concurrence sudden death is eliminated, indicating that this effect is directly linked to the inclusion of zero and double excitation components in the dynamics. The qubits population behavior in time is displayed in Fig.  2, using a = 1 2 . The qubits have the same energy gap, thus in the P subspace the two states are degenerate and their population is identical at all times, independently of K. In contrast, in the Q subspace the energy difference between the states is significant, larger than the temperature, T /2ǫ < 1; the tunneling element is given by γJ = 1 2 , with γ as the anisotropy in the qubit-qubit coupling. In such a situation we expect the steady-state population of the spin-up state to be significantly smaller than the ground state population, as indeed we observe in Fig. 2(b). An interesting observation is the phenomenon of population inversion between the zero and the double excitation states before steady-state sets in. This behavior occurs roughly up to a timescale that is inversely proportional to the Kondo parameter K, independent of the temperature. For the same set of parameters Fig. 3 presents the coherence dynamics in the two subspaces. Generally, coherences are diminishing with the increase of K. Given the population and coherence dynamics, we display in Fig. 4 the concurrence, calculated using Eq. (24), manifesting a rich dynamics. The following characteristic's are of particular interest: (i) The birth-time of the concurrence, (ii) its oscillations, (iii) the occurrence of sudden death and revival, and (iv) the steady-state value. We now explain those properties.
Time-zero concurrence. The particular initial condition used here, a = 1 2 , results in C(t = 0) = 0. This is because while we are using maximally entangled states within each subspace as an initial condition, the entanglement between the two qubits themselves is zero initially, since all relevant reduced density matrix elements, necessary for evaluating Eq. (24), are identical [33].
Delayed sudden birth. When (ρ S ) 00,00 ∼ (ρ S ) 11,11 , a situation taking place at, and close to, the initial time, the concurrence should be zero, given the positivity condition that limits the value of off-diagonal elements. For small K, the time it takes the system to depart from its initial-equal population state is prolonged compared to a large-K case, thus, the concurrence birth-time is delayed with respect to the large K behavior. Interestingly, the delay in the birth time does not extend linearly with K. Rather, the delay is significant for both K = 0.05 (weak system-bath coupling) and for K = 0.5 (intermediate coupling), while it is shorter for in-between values, K ∼ 0.3; The reason is that the delay time is a nontrivial function of both the time it takes the coherences to establish, and the time it takes the population to significantly depart from the initial (equally populated) setup.
Oscillations. The oscillatory nature of C in time, best manifested for K ≤ 0.1 reflects the Rabi-type oscillations of the diagonal elements (ρ S ) 00,00 and (ρ S ) 11,11 . When these elements are similar in value, the concurrence drops, and even dies during a certain time interval, depending on the magnitude of the coherence at that time.
Steady-state value. If the two qubits are isolated from thermal effects (K=0, right panel of Fig 4), the concurrence oscillations reflect the nature of the population and coherences dynamics, depicting Rabi oscillations. The qubits behavior under a dissipative thermal bath is notably distinct: Since both population and coherences approach a constant at long time, the concurrence reaches a steady-state value as well. It predominantly reflects the magnitude of the coherence (ρ S ) 10,01 in the long time limit since the population weakly depends on K at long time, see Fig. 2(b). Interestingly, for the present a = 1 2 case the steady-state value of the concurrence is almost identical at weak system-bath coupling, K = 0.05 − 0.3. It significantly degrades around K = 0.4 − 0.5. Beyond that, it is identically zero.
Sudden death and revival. Based on Figs. 2-4, we can draw general conclusions regarding the process of entanglement sudden death. The effect is directly linked to the existence of population in the zero and double excitation states. If the dynamics within the Q space is eliminated all together (a = 0), the concurrence is only controlled by the absolute value of the coherence |(ρ S ) 01,10 |, Eq. (25). This quantity does not manifest an oscillatory behavior: Under the present initial condition it starts at a large value, touches zero at a particular time, then grows again to a certain extent. (Under a different initial condition, e.g., |(ρ S ) 01,10 (0)| = 0, the entanglement will systematically grow, up to the steady-state value). In contrast, when the double excitation state is initially populated, oscillations between the states in the Q subspace largely occur, if system-bath coupling is weak. Then, the competition between the two terms in F 2 , see Eq. (21), can result in a disentanglement over a finite time interval. Fig. 5 displays the concurrence using different initial conditions, by playing with the parameter a. This modifies the weight of the zero and double excitation states in the dynamics. When a = 0 and K = 0 the entanglement is maximal (C = 1) at all times. For finite K, keeping a = 0, it dies at the particular point at which |(ρ S ) 01,10 | = 0. Beyond that, it recovers to a value close to 1. When we include the Q states, e.g., by taking a = 0.2, we observe the effect of entanglement sudden death, over a certain time interval. The duration of this interval grows when a is further increased up to a ≤ 1/2. Beyond this point the coherence in the P subspace may dominate over the population in the Q subspace, resulting in a positive value for F 1 , eliminating entanglement sudden death. The behavior at intermediate-strong system-bath coupling, K = 0.6, is included in Fig. 6, demonstrating that temporal oscillations are washed out. The dynamics at even larger K is similar in trends, with reduced concurrence value.
The role of the temperature is displayed in Fig. 7. At high temperature the concurrence is zero. At intermediate values, T < ǫ ∼ 1 we find that its sole effect is a shift down of the qubit entanglement with increasing temperature. All other features (birth time, oscillation) stay intact. The simulation could not be performed at temperatures below T ∼ 0.1 due to convergence issues in QUAPI.
We can readily study the concurrence under different initial conditions for the P and Q subspaces, not necessarily in the form of Bell states, as long as Eq. (14) is obeyed. In particular, using a diagonal state for the time-zero reduced density matrix, similar features as those discussed above were obtained.

V. UNIVERSAL FEATURES AT LONG TIME
The long time behavior of the concurrence, representing the equilibrium limit, is displayed in Fig. 8  K, the system-bath coupling parameter, and the initial state preparation ratio a, see Eq. (22). We note that the concurrence can be significant in both the weak and strong coupling regimes, as long as the system evolves predominantly in either the P or Q subspaces. We now show that at weak coupling, K ≪ 1, and at low temperatures, T < Jγ, for a broad range of parameters (as we explain below), the following general result holds The important implication of this result is that to the lowest order in K the concurrence deviates from unity due to the occupation of the zero and doubly excited states in the system. This trend was observed before, e.g., in Ref. [29]. However, here, for the first time, it is justified analytically, based on the spin-boson model behavior [31]. We derive Eq. (26) by studying the long-time limit of F 1 , as it dictates the concurrence when a ≤ 1 2 , see Eq. (24). In the biased case, weak coupling theory (beyond the noninteracting blip approximation) provides [31] Q z = (ρ S ) 00,00 − (ρ S ) 11,11 ∼ a ǫ in the thermodynamic limit. Here ∆ 2 b = ǫ 2 + ∆ 2 ef f , ∆ ef f is a nontrivial function of K, ω c , and the bare tunneling element in the Q subspace, Jγ. In the weak coupling limit we can write ∆ ef f ∼ Jγ, thus ∆ b ∼ ǫ 2 + J 2 γ 2 . Manipulating the polarization, we obtain the relevant term, The other element in F 1 is the coherence in the P subspace. In the long time limit it satisfies [31] |(ρ S ) 01, 10 where Ω 2 = J 2 + 2J 2 Kµ; the proportionality factor obeys µ = ψ(iJ/πT ) − ln(J/T ) with ψ as the digamma function [31]. As expected, the equilibrium concurrence depends on the environmental temperature, leading to entanglement degradation at high T , as observed in Fig. 7. Considering the low temperature case, T < J, Jγ, we note that the trigonometric term in both Eq. (28) and Eq. (29) is close to unity. If we further work in the region ǫ < Jγ, the square root expression in Eq. (28) gets close to 1. Under these broad conditions, the concurrence reduces to One should note that the K dependence is more subtle than the simple linear scale attained here, since the tunneling element J should be corrected by K in a nontrivial manner [31]. The simple result (30) provides us with some basic-interesting rules for building a long-time concurrence within the range of parameters mentioned above: (i) It decays linearly with the overall population placed in the Q (zero and double excitation) subspace. (ii) The reservoir temperature does not significantly affect it. (iii) It does not depend on the qubit interaction energy. We note again that these observations are valid for a < 0.5, as long as T < J, Jγ, ǫ < Jγ and K ≪ 1.
When a > 0.5, the concurrence is determined by the competition between F 1 and F 2 , see Eq. (20). The numerics then suggests that C a> 1 2 (t → ∞) ∼ 2a − 1 holds, for similar energy parameters.
We conclude this section by emphasizing the implication of Eq. (26): One could set the steady-state entanglement in a dissipative system, by controlling the initial population in the P and Q subspaces.

VI. CONCLUSIONS
Using exact numerical tools, we simulated the time evolution of two qubits immersed in thermal environments, considering a class of initial states for the subsystem. This task was achieved by reducing the two qubits-bath model into two spinbath systems, whose dynamics could be readily followed separately. Using Wootters' formula for the concurrence [16], we quantified the degree of the qubits entanglement in time, exposing rich dynamics, including oscillations, delayed sudden birth, sudden death, and revival. Specifically, we showed that the occurrence of entanglement sudden death can be traced down to the initial population of the zero and double excitation states. The steady-state behavior was discussed in the weak coupling limit.
Our results are significant for several reasons. First, we exposed a general mapping between an interacting two-qubit system embedded in a bath, and two spin-bath models, allowing us to simulate the dynamics of the original model using a numerically exact method that was developed for studying the prominent spin-boson case. Second, based on our mapping scheme, we calculated the concurrence measure and demonstrated the essential role of the environment in generating a stationary entanglement between the (interacting) qubits. By engineering the environment and its interaction with the system one could tune the degree of disentanglement in the system [30]. Earlier studies in this field have mostly treated a simpler version of our model, ignoring qubit-qubit interaction energy, further utilizing perturbative treatments. To the best of our knowledge, our work is the first to calculate the concurrence exactly in a dissipative and interacting qubit model. Future studies will focus on the dynamics of non-classical correlations beyond the entanglement measure, evaluating quantum discord [32]. This could be done by relying on the X-form of the reduced density matrix [33,34]. With this at hand, we plan to study the dynamics of classical and quantum correlations in the qubit system, specifically, to investigate classical and quantum decoherence mechanisms and the possible transition between them [35].