Creation and protection of entanglement in systems out of thermal equilibrium

We investigate the creation of entanglement between two quantum emitters interacting with a realistic common stationary electromagnetic field out of thermal equilibrium. In the case of two qubits we show that the absence of equilibrium allows the generation of steady entangled states, which is inaccessible at thermal equilibrium and is realized without any further external action on the two qubits. We first give a simple physical interpretation of the phenomenon in a specific case and then we report a detailed investigation on the dependence of the entanglement dynamics on the various physical parameters involved. Sub- and super-radiant effects are discussed, and qualitative differences in the dynamics concerning both creation and protection of entanglement according to the initial two-qubit state are pointed out.


Introduction
Quantum systems may present correlations of both quantum and classical nature. Entanglement captures quantum correlations due to the non-separability of the system state [1][2][3]. The presence of these correlations is connected to the rise of non-local effects in quantum theory [4] and has been recognized as a key resource in several fields of quantum technology, including quantum computing [5], quantum cryptography [6], quantum teleportation [7] and quantum metrology [8]. A main obstacle to the concrete exploitation of quantum features in the above applications is the detrimental effect of environmental noise [9]. The unavoidable coupling with degrees of freedom of the surrounding environment generally leads to a decay of quantum coherence properties [10], preventing the possible exploitation of quantum correlations present in the system. Considerable effort has been made to understand the effects of environmental noise on the dynamics of correlations present in an open quantum system [11][12][13][14], and to contrast the natural fragility of quantum coherence properties [15][16][17]. Reservoir engineering methods have pointed out the possibility of changing the perspective from reducing the coupling with the environment to modifying the environmental properties in order to manipulate the system of interest, thanks to its proper dissipative dynamics [18][19][20]. Other approaches exploit the effect of measurements and feedback to drive the systems toward a target state [21,22].
A possible way to create quantum correlations between two systems is to make them interact with a common environment [23], which can also cause a revival of entanglement [24].
In the case of two emitters in a common vacuum or thermal electromagnetic field, in the absence of matter close to them, the mediated interaction plays a role over distances of the order of the common transition wavelength [25,26]. It has been found that the presence of plasmonic waveguides near the emitters can allow a mediated interaction over larger distances [27] whose effect on the entanglement dynamics has been discussed [28]. However, at thermal equilibrium, the dynamic creation of entanglement eventually ceases at some time and the system thermalizes toward a thermal state which is a classical mixture. Steady entanglement can be instead generated by adding the action of an external driving laser [29].
The influence of several independent reservoirs at different temperatures, whose emission does not depend on their internal structure (material or geometry), has been considered in several contexts, including generation of entanglement in non-equilibrium steady states, both in the case of a few spins [30][31][32] and of a chain of spins [33][34][35] and in the context of quantum thermal machines [36,37].
However, in a realistic configuration, the actual reflection and transmission properties of the bodies surrounding the quantum emitters should be taken into account, and may become particularly relevant if the emitters are placed close to the bodies (near-field effects). New possibilities emerging in such realistic systems out of thermal equilibrium have been recently pointed out in different contexts ranging from heat transfer [38,39], to Casimir-Lifshitz forces [40][41][42][43][44][45]. There, radiation fields out of thermal equilibrium in configurations of a quite general nature have been characterized in terms of the correlators of the total field depending on the scattering matrices of the bodies composing the total system [46,47]. In the case of single emitters in such environments, new tools exploiting the absence of thermal equilibrium to manipulate the atomic dynamics realizing inversion of population and cooling of internal atomic temperature have been pointed out [48,49]. Recently, the case of two quantum emitters has also been analyzed, pointing out a new remarkable mechanism to generate and protect entanglement in a steady way in systems out of thermal equilibrium [50].
In this paper, we report a detailed investigation of this phenomenon by studying the internal dynamics of a system composed of two quantum emitters (real atoms or artificial ones such as quantum dots or superconducting qudits) placed in front of an arbitrary body embedded in a thermal radiation whose temperature is different from that of the body. The paper is organized as follows. In section 2, we describe the physical model under investigation and we derive a master equation for the general case of two N -level emitters. In section 3, we derive closedform expressions for the functions governing the dynamics, in terms of the scattering matrices of the body and valid for arbitrary geometrical and material properties. In section 4, we develop these expressions in the case when the body is a slab of finite thickness. From section 5 onwards we specialize our analysis to the case of a two-qubit system, comparing cases in and out of thermal equilibrium. We point out the occurrence of peculiar phenomena emerging out of thermal equilibrium such as the generation of steady entanglement and a simple interpretation for this phenomenon is presented for a particularly interesting case. The general case of arbitrary values of the parameters is then discussed in section 6. In section 7, we draw our conclusions.

The model
We consider a system made of two quantum emitters q = 1, 2 interacting with an environment consisting of an electromagnetic field which is stationary and out of thermal equilibrium. This is generated by the field emitted by a body (M) at temperature T M of arbitrary geometry and Physical configuration: two quantum emitters close to an arbitrary body whose temperature T M is kept fixed and different from that of the surrounding walls T W . The two emitters are placed at R 1 = (r 1 , z 1 ) and R 2 = (r 2 , z 2 ), where r 1 and r 2 are vectors in the x y plane. In the figure, we choose the x-axis along the direction r 1 − r 2 and x 2 = 0, naming r 12 = |r 1 − r 2 | = x 1 (this choice of the reference system is used in section 4 in the specific case when the body is a slab). dielectric response and by the field emitted by the far surrounding walls (W) at temperature T W , which is eventually transmitted and reflected by the body itself (see figure 1).
T M and T W are kept fixed in time, realizing a stationary configuration for the electromagnetic field. The surrounding walls have an irregular shape and are distant enough from the body and the emitters so that their field can be treated at the emitters' locations, in the absence of the body, as blackbody radiation independent of their composition. This is not true for the field emitted by the body M which cannot be treated as a blackbody since its radiation depends on its actual properties such as its geometry and its dielectric function. The total Hamiltonian has the form where H S and H E are the free Hamiltonians of the two emitters and of the environment. The interaction between the emitters and the field, in the multipolar coupling and in dipole approximation, is [51] where D q is the electric-dipole operator of emitter q and E(R q ) is the electric field at its position R q . We first consider the general case in which each emitter has N q internal levels |n q where n ∈ {1, . . . , N q } of frequency ω q n (ordered by increasing energy). Given two arbitrary levels n and m, their frequency difference is indicated by ω q nm = ω q n − ω q m and the transition matrix element of the dipole operator by d q mn = q m|D q |n q . The free Hamiltonian of the two emitters is where ( q n ) = |n qq n| are the projectors associated with each eigenvalue q n =hω q n (possibly degenerate) of H q . The dipole operator of emitter q in the interaction picture, where σ q mn = |m qq n| and ω q nm 0. By moving to the interaction picture, we obtain for H I where the time-dependent electric field is given by In the following each mode of the field is identified by the frequency ω, the transverse wave vector k = (k x , k y ), the polarization index p (taking the values p = 1, 2 corresponding to transverse-electric (TE) and transverse-magnetic (TM) polarizations respectively) and the direction of propagation φ = ±1 (shorthand notation φ = ±) along the z-axis (see figure 1). In this approach, the total wavevector takes the form K φ = (k, φk z ), where the z component of the wavevector k z is a dependent variable given by k z = ω 2 c 2 − k 2 , where k = |k|. The explicit expression of the field at an arbitrary point R is where a single-frequency component reads where E φ p (k, ω) is the field amplitude operator associated with the mode (ω, k, p, φ). For the TE and TM polarization vectors appearing in (7), we adopt the following standard definitions: wherex,ŷ andẑ are the unit vectors along the three axes andk = k/k.

Master equation
The starting point to study the dynamics of the two emitters is, in the interaction picture, the von Neumann equation for the total density matrix ρ tot (t): The reduced density matrix of the two emitters is given by ρ = Tr E [ρ tot ], where Tr E denotes the trace over the degrees of freedom of the environment. To derive a master equation for ρ we follow the procedure described in [9] for the case of one emitter by extending it to our system made of two emitters. We name ω q an arbitrary transition frequency of emitter q (positive and negative). In general, several transitions can be characterized by the same frequency ω q both because of degeneracy and/or the occurrence of equidistant levels. We rewrite each Cartesian component of the dipole operator, [D q ] i (i = {x, y, z}), as where A q i (ω q ) and A q † i (ω q ) turn out to be the eigenoperators of H q with frequencies −ω q and +ω q , respectively, i.e.
In the central term of (10), the first sum is over all of the frequencies ω q while the second is over all of the couples of energy eigenvalues q n and q m of H q such that q n − q m =hω q . Following [9], it is useful to rewrite H I (t) of (5) in terms of the eigenoperators A q i (ω q ) as From (10) it follows that the vector where the sum is over all of the couples n and m such that ω q nm = ω q . By applying to the case of two emitters the standard procedure for the microscopic derivation of a master equation reported in [9], under Born, Markovian and rotating-wave approximations 4 , one can obtain (by using also the condition E i (R, t) = 0) in the Schrödinger representation d dt where ω 0, being terms with positive or negative ω associated, respectively, with downward and upward transitions. In the above equation, for q = q the sum q,q ,ω is over all of the common frequencies ω q = ω q = ω (this condition derives from the rotating wave approximation) while for q = q it is over all of the transition frequencies of each emitter, and γ qq ii (ω), and s qq ii (ω) are defined by where the field correlation functions enter in the function qq ii (ω). It follows that The initial state of the total system in (13) is assumed to be factorized, ρ tot (0) = ρ(0)ρ E . In the case ρ E is a stationary state of the environment ([H E , ρ E ] = 0) the correlation functions are homogeneous in time, that is does not depend on time. The functions defined in (14) appearing in the master equation (13) depend thus only on the field correlation functions E i (R q , s)E i (R q , 0) , whose computation out of thermal equilibrium will be the subject of sections 3 and 4. We now explicitly write the master equation (13) in the case of absence of degenerate and equidistant levels in each emitter, when the definition of the eigenoperators A q (ω q ) (12) reduces to A q (ω q ) = d q mn σ q mn (each ω q corresponds to only one couple of energy eigenvalues { q m , q n }). For this purpose, we develop the sum over ω in (13), which for each |ω| runs over ω (downward transitions) and −ω (upward transitions), as ω f (ω) = ω>0 f (ω) + ω>0 f (−ω). From now on ω always indicates a positive frequency and we drop '> 0' in the sums over ω. By introducing this new convention and using the explicit form for A q i (ω) (12), we can recast (13) as where the sum q,q ,ω in the second line is relative to all of the transition frequencies of each emitter for q = q and only to the common transition frequencies for q = q , (m, n) and (m , n ) individuate, respectively, the transition of each emitter corresponding to the frequency ω, and we have defined the functions We remark that it holds [ qq (ω)] * = q q (ω) and [ qq (ω)] * = q q (ω). In (16), function qq (ω) represents a coherent (dipole-dipole) interaction between the emitters mediated by the field while dissipative effects enter through the functions. In particular, qq (±ω) are individual (q = q ) and common field-mediated collective (q = q ) emitter transition rates, related to both quantum and thermal fluctuations of the electromagnetic field at the emitters' position.

Emitters close to an arbitrary body
Here, we derive the field correlation functions needed to compute the functions in (17) for non-equilibrium configurations in the case of an arbitrary body and multilevel emitters. These functions will depend on the two temperatures T M and T W and on the material and geometrical properties of the body as well. We follow the derivation discussed in [47] in the more general case of two bodies and three temperatures and the derivation relative to a single quantum emitter in the presence of a single body and two temperatures [49]. Here, we extend the latter derivation to the case of two quantum emitters. Some of the computations involved are reported in appendix A.
The starting point is to decompose on the right side of the body where the emitters are located, the amplitude operators of total field modes propagating in the two directions z > 0 and z < 0 in terms of the fields emitted by the surrounding walls (W) and by the body (M). For a given set (ω, k, p), we have for the two directions where we made the dependence on ω, k and p implicit. The total field E − propagating toward the body (i.e. toward the left) is equal to the field emitted by the walls E (W)− coming from the left, while the total field E + propagating toward the right results from the field E (M)+ directly produced by the body, the transmission through the body of the field E (W)+ emitted by the walls coming from the left and the reflection by the body of the field E (W)− coming from the right (see figure 2). The operators R and T are the reflection and transmission scattering operators associated with the right side of the body, whose explicit definition can be found for example in [47]. They connect any outgoing (reflected or transmitted) mode of the field to the entire set of incoming modes. By using (18) one can write the total field correlators in terms of the correlators of the fields emitted by each source.
The source fields have been characterized as in [47] by assuming that for the body M and the walls W a local temperature that remains constant in time can be defined and that the emission process of the body is essentially not influenced by the presence of the external radiation impinging on the body itself. This assumption leads to the hypothesis that the part of the total field emitted by the body is the same as it would be if the body were at thermal equilibrium with the environment at its own temperature so that the correlators of the field emitted by each body can still be deduced by using the fluctuation-dissipation theorem at its local temperature.
Under this assumption, the following symmetrized correlation functions [ AB sym = ( AB + B A )/2] have been derived where we have introduced In the above equation (pw) and (ew) are the projectors on the propagative (c k < ω, corresponding to a real k z ) and evanescent (c k > ω, corresponding to a purely imaginary k z ) sectors respectively. By combining (18) and (19), in appendix A a general expression for the total correlation functions in frequency space has been derived in (A.7). This expression can be used to compute the functions γ qq ii (ω), γ qq ii (−ω) and s qq ii (ω) entering in (17), by exploiting their connection with the correlation functions between frequency components of the total field given in (A.8).
To move to the final expression of the functions in (17), we first rewrite the antinormally ordered correlation functions (A.7) as from which the normally ordered correlation functions are obtained by replacing [1 + n(ω, T i )] with n(ω, T i ) and by taking the complex conjugate (this procedure derives from Kubo's prescription as explained in appendix A) and where we have introduced two α functions which do not depend on the temperatures and on the dipoles, and do depend on the geometrical and material properties of the body through the operators R and T : ii is real as expected, being proportional to the imaginary part of the Green function (see (C.3)). Now we can compute the transition rates in (17), by using (A.8), (21) and (22), where q 0 (ω) = |d q mn | 2 ω 3 3hπ 0 c 3 is the vacuum spontaneous-emission rate of transition |n q → |m q of emitter q and we have introduced the new functions depend on the choice of emitters' dipoles. In the case of two qubits, which will be treated in sections 5 and 6, there is only one transition for each emitter and the above equations (24) and (25) hold with the notation d q mn = d q . With regards to the function qq (ω) we obtain, by using (A.8), (21) and (22), where we used the properties [α qq ii . It follows that qq (ω) does not depend on the presence or the absence of thermal equilibrium, being independent of the temperatures. By using the relation between α functions of (23) and the Green's function of the system in (C.3) derived in appendix C, the integration over frequencies in (26) can be achieved by using the Kramers-Kronig relations connecting real and imaginary parts of Green's function:

Emitters close to a slab
We now specialize the derivation of the previous section to the case when the body is a slab of finite thickness δ, defined by the two interfaces z = 0 and −δ (see figure 3). In this simple case, explicit expressions for the transmission and reflection operators can be exploited [46,47]. Owing to the translational invariance of a planar slab with respect to the x y plane, the slab reflection and transmission operators, R and T , are diagonal and equal to where the Fresnel reflection and transmission coefficients modified by the finite thickness δ are given by (we recall that p = 1, 2 corresponding to the TE and TM polarizations) In the previous equations we have introduced the z component of the K vector inside the medium ε(ω) being the dielectric permittivity of the slab, the ordinary vacuum-medium Fresnel reflection coefficients as well as both the vacuum-medium (denoted with t) and medium-vacuum (denoted witht) transmission coefficients After replacing the matrix elements (28) in (23) we obtain for the α functions (we choose the x-axis along the vector r 1 − r 2 whose coordinates in the plane x y are then (r 12 , 0), being where, using the fact that ρ p (k, ω) and τ p (k, ω) are independent of θ (the angle formed by k and the x-axis in the plane x y), we have introduced the angular integrals where r 21 = −r 12 . The matrix elements different from zero are, for where J n (x) is the nth-order Bessel function of the first kind. For ii become diagonal and reduce to the vectors defined in (55) of [49] in the case of a single emitter.
To simplify the functions [α ii . By using the angular integrals (34), equation (33) can thus be rewritten as where we have introduced the integral matrices For q = q , α qq W (ω) and α qq M (ω) coincide with the functions defined in (56) of [49] in the case of a single emitter. For q = q , in the limit R q → R q , [α qq W (ω)] ii and [α qq M (ω)] ii tend to their values in the case of a single emitter placed in R q . In the limit, the distance between the two emitters goes to infinity, both α qq W (ω) and α qq M (ω) go to zero. In the limit of |r qq | → ∞, this is due to the fact that the functions [N qq p (k, ω)] φφ go to zero (for |x| → ∞, it is J 0 (x) → 0, J 1 (x) → 0 and J 2 (x) → 0). In the limit |z q − z q | → ∞, this is due to the presence of an oscillating function whose frequency goes to infinity in the integrals A(ω), B(ω) and C(ω) (this can be seen explicitly by integrating by parts) and to the presence of the exponential function going to zero in the integral D(ω).
Concerning the function qq (ω), in order to develop its expression in (27) one has to compute the real part of the Green function in terms of the scattering operators R and T . This is done in appendix C where a free term qq 0 (ω) (see (C.11)) remaining in the absence of matter has been isolated from a reflected part qq R (ω) (see (C.12)), qq (ω) = qq 0 (ω) + qq R (ω). By using the expressions for R and T in the case of a slab and the angular integrals in (34), one can derive, starting from (C.12): Equation (27) can be thus cast in the form where we have introduced the integral matrices We observe that the limit case when the body is absent is discussed in appendix B, where known expressions for qq (ω) and qq (ω) are retrieved.

Two-qubit system
From now on we specialize our investigation to the case of two emitters (qubits) characterized by two internal levels |1 ≡ |g and |2 ≡ |e with the same transition frequency ω = ω 1 e − ω 1 g = ω 2 e − ω 2 g . In this case, master equation (16) reduces to where we used [σ q gg , ρ] = −[σ q ee , ρ], so that and where the functions S qq (±ω), qq (ω) and qq (±ω) are defined in (17) for the specific case {m, n} = {1, 2} (in the following we use the notation d q 12 = d q ). We observe that master equation (41) can also describe the case in which the emitters' frequencies are close enough (but not identical) so that rotating wave approximation used in the derivation of (13) still holds. This typically occurs when the frequency difference is much smaller than the average frequency [26]. The operator δ S (42) represents a shift of energy levels, the renormalized transition frequencies being equal toω 1 e − ω 1 g = ω + S 11 (ω) − S 11 (−ω) and ω 2 e − ω 2 g = ω + S 22 (ω) − S 22 (−ω). In the following, these shifts will not play any role. To discuss the properties of (41), we will use two different bases, the decoupled basis {|1 ≡ |gg , |2 ≡ |eg , |3 ≡ |ge , |4 ≡ |ee } and the coupled basis where we have introduced the collective antisymmetric |A and symmetric states |S . The coupled basis is the one diagonalizing the effective Hamiltonian, H S + δ S + q =q h qq (ω)σ q † ge σ q ge , appearing in the first line of (41). In particular, the sign of 12

X states
In the decoupled basis we can distinguish elements along the two main diagonals of the twoqubit density matrix from the remaining ones because they are not connected through master equation (41). We thus focus our attention on the class of X states, having non-zero elements only along the main diagonal and anti-diagonal of the density matrix (we use the notation Bell, Werner and Bell diagonal states belong to this class of states [55]. X-structure density matrices are found in a wide variety of physical situations and are also experimentally achievable [56]. For example, X states are encountered as eigenstates in all of the systems with odd-even symmetry such as in the Ising and the XY models [57]. Moreover, in many physical evolutions of open quantum systems an initial X structure is maintained in time [58], as it is in our case. Terms outside the two main diagonals initially populated would be eventually washed off asymptotically. In the following, the two-qubit state will always have an X structure.

Concurrence
We shall quantify the entanglement in the two-qubit dynamics by evaluating the concurrence, C(t) (C = 0 for separable states, C = 1 for maximally entangled states) [59]. For the X states it takes the form [58] C(t) = 2 max{0, K 1 (t), K 2 (t)}, The master equation (41) always induces an exponential decay for ρ 14 (t), so that in the steady state only K 1 (t) could be responsible for having C(∞) > 0.
To discuss new phenomena emerging out of thermal equilibrium, it will be instructive to rewrite K 1 (t) in terms of the populations in the coupled basis (we use the notation ρ IJ = I |ρ|J and ρ I = I |ρ|I ): We will see that out of thermal equilibrium, it is always ρ AS (∞) = 0, but ρ S (∞) and ρ A (∞) can differ, so that K 1 (∞) could be positive.

Thermal equilibrium
When T W = T M ≡ T , master equation (41) describes the thermalization toward the thermal equilibrium state, which is diagonal with the four steady populations given by where Z eq = [1 + 2n(ω, T )] 2 . By moving to the coupled basis, the thermal state remains As a mathematical remark, we note that the thermal state is always reached asymptotically except if the identities 11 (±ω) = 22 (±ω) = 12 (±ω) = 21 (±ω) ≡ (±ω) are strictly verified. In this peculiar case, both in and out of thermal equilibrium, the steady state depends upon the initial state and may be entangled. In particular, it is diagonal in the coupled basis with populations equal to where Z = (−ω) 2 + (ω) (−ω) + (ω) 2 . Apart from this case, at thermal equilibrium the steady state is always a thermal state, thus not entangled. We can see it by looking at the concurrence (44) which is zero being ρ 23 (∞) = 0. This can also be seen in the coupled basis, where ρ AS (∞) = 0 and ρ S (∞) = ρ A (∞), so that K 1 (∞) (45) is negative.

Out of thermal equilibrium: an instructive case
When T W = T M , qualitative differences emerge in the dynamics and in the steady states. To highlight these new features, we first consider a simple case where a clear physical interpretation in terms of |S and |A is available. This is the case when 11 (±ω) = 22 (±ω) ≡ (±ω) and 12(21) (±ω) are real. These conditions are verified, for example, in the case of identical qubits, with d 1 = d 2 ≡ d, placed in equivalent positions with respect to the body (in the case of a slab, z 1 = z 2 ) and with d real and having components different from zero either only along the z-axis or only along the plane x y. In this case, master equation (41) gives in the coupled basis a set of rate equations for the populations, which are decoupled from the other density matrix elements: Here, the coefficient 0 (ω) has been absorbed by the time variable in the derivative, which is now dimensionless, and we have used the relations where α W(M) (ω) ≡ α 11 W(M) (ω) = α 22 W(M) (ω). We remark that qq (ω) does not enter in the rate equations (48), which are schematically represented in figure 4. We observe that with each decay channel from |E to |G we can associate distinct effective temperatures T S and T A confined between T W and T M in correspondence to the effective number of photons n S and n A , which have the property of being confined between n(ω, T W ) and n(ω, T M ) [49].
Concerning the coherences in the second diagonal: which give for each coherence an exponential decay modulating oscillations due to 12 (ω). The stationary solution of (48) is where Z neq is the sum of the elements of the vector in the second line of the above equation. Out of equilibrium ρ 23 (∞) neq is different from zero and is given by where we see easily how it tends to zero at thermal equilibrium when n S = n A . By using (52) in (44) and (45), we obtain for the steady concurrence: Simplifying S , C(∞) neq becomes the function of the three dimensionless quantities A / S , n S and n A . This dependence is discussed in figure 5, where C = C(∞) neq is depicted as a function of n A for n S = 0.001 and for several values of A / S , as indicated in the legend. We observe that by decreasing the value of A / S , higher values of C are reachable at higher values of n A . The maximum value of C is 1/3, which can be obtained in the limits A / S → 0, n S → 0 and n A → ∞. In particular, the corresponding maximally entangled state, which is a statistical mixture of the ground and of the antisymmetric state with weights, respectively, equal to 2/3 and 1/3, has also been found in [35]. For smaller values of n S the behavior remains almost identical, while by increasing its value, the values of C decrease progressively. We remark that identical behavior is found in the opposite case, i.e. when S / A → 0, the case in which the roles of states |S and |A are inverted. This can be achieved by looking for values of the various parameters such that α 12 M(W) (ω) is negative and very close to α M(W) (ω) in order to make the ratio S / A very small. Figure 5 describes the generation of steady entangled states emerging only in the absence of thermal equilibrium. Two main conditions must be fulfilled, the first being to have small values for A / S and the second to realize (quite) different effective temperatures for the two decay channels, which can be achieved only in the absence of thermal equilibrium.

Numerical investigation
Here, we report the numerical investigation concerning the case treated in section 4 when the body close to the emitters is a slab of finite thickness δ. According to (13), a relevant parameter involved in our investigation concerning the role of the body is the value of the dielectric permittivity at the common transition frequency of the two qubits. As a material we choose silicon carbide (SiC) whose dielectric permittivity ε(ω) is described by using a Drude-Lorentz model [53] characterized by a resonance at ω r = 1.495 × 10 14 rad s −1 and where ε ∞ = 6.7, ω l = 1.827 × 10 14 rad s −1 and = 0.009 × 10 14 rad s −1 . This model implies a surface phonon-polariton resonance at ω p = 1.787 × 10 14 rad s −1 . A relevant length scale in this case is c/ω r 2 µm while a reference temperature ishω r /k B 1140 K. We will assume that ε(ω) does not vary much in the interval of temperatures considered. In the following study, we explore a region of parameters much wider than that allowing the analytical description of section 5.4.

Steady configurations
We first focus on the properties of steady states. In particular, we are interested in the amount of entanglement present asymptotically, which is quantified by the concurrence (44). This analysis is supported by an analytical solution of the steady state of (41), which is not reported here, since it is particularly cumbersome. In figure 6(a), we plot the maximum of steady concurrence obtained for an interval of transition frequencies ranging from 0.3 ω r to 1.7 ω r , in the case of δ = 0.01 µm. In our numerical sample, z 1 and z 2 may vary between 0.05 and 50 µm and r 12 between 0 and 15 µm. The two temperatures range in an interval such that the associated number of photons is between 0 and 3 (we checked that larger values are not needed). The red curve is relative to the case of dipoles oriented along the z-axis, while the green curve to the case of dipoles oriented along the x-axis (in figure 9 we will show that this is the best choice if we limit ourselves to directions lying in the x y plane). Higher values of concurrence are obtained immediately before/after the resonance frequency ω r . The red configuration always gives better results except around the surface phonon-polariton frequency ω p 1.2ω r where choosing the dipole directions along the x-axis is the best choice. The values of the parameters corresponding to each maximum vary with frequency. The best configuration is always characterized by values of n(ω, T W ) close to zero and n(ω, T M ) between 1 and 3. Smaller values of n(ω, T M ) are needed in the green curve. The zone where to place the qubits is around 1 µm from the slab at 0.3ω r , gradually decreasing (specially after ω r ) down to 0.25 µm at 1.7ω r . For the red curve the best choice is always z 1 close to z 2 (in our numerical sample, we limit the minimal distance at the order of 0.1 µm) and r 12 = 0, while for the green curve it is z 1 = z 2 and r 12 small (of the order of 0.01 µm). This means that the best configuration is when the interatomic axis is aligned with the dipoles direction. For ω around ω p we point out the occurrence of larger values of C for different values of δ, points indicated with a cross above the green curve. In the absence of large values of C in correspondence to the canonical choice of the parameters described above, small values of C become evident for a different set of parameters. This corresponds to larger values of δ (of the order of 1 µm or more), z 1 2 µm, z 2 4 µm and r 12 0.5 µm. In part (b), we plot the dependence of C on δ for several values of ω as indicated in the figure. The maximum of C is always obtained close to δ = 0.01 µm, which is the value chosen in part (a), except around ω r where much smaller values of δ are required. This explains why in part (a) concurrence decreases around ω r for δ = 0.01 µm.
In figure 7, we plot the steady concurrence as a function of the position of the second qubit z 2 and of the slab temperature T M for four different values of r 12 . From (a) to (d) the two-qubit distance [r 2 12 + (z 1 − z 2 ) 2 ] 1/2 increases leading to a progressive decrease of the values of concurrence generated. A maximum of C 0.24 is obtained in part (a) for z 2 1.3 µm and T M 1100 K. The white lines correspond to the case z 2 = z 1 for which equation (54) holds for concurrence. In part (a), the maximum along the white curve is C ∼ 0.222 in correspondence to A / S ∼ 4.6 × 10 −7 , n S ∼ 0.02 and n A ∼ 1.56 which correspond to effective temperatures for the two decay channels T S ∼ 90 K and T A ∼ 690 K. We observe that very high temperatures are considered in this plot only to highlight the entire region where steady entanglement is present. At unphysical temperatures (e.g. above the melting temperature), the plot is only indicative of what would occur if a different material was chosen such that similar values of ε(ω) (for SiC it is ε(0.3 ω r ) ∼ 10.3 + 0.00721i) were encountered at lower frequencies. In this case, similar behavior for steady concurrence at lower temperatures is expected. However, we remark that in our case values of C higher than 0.14 are already present at T M 500 K in part (a).
In figure 8, we discuss the behavior of n S , n A and A / S , appearing in (54), as a function of z. We plot in part (a) n S and n A as a function of z = z 1 = z 2 and compare them with the value of n(ω, T ) computed at the minimal (here T W = 30 K) and maximal (here T M = 1200 K) temperature considered. The temperatures and the other parameters are equal to the ones giving the maximum of concurrence in figure 7(b) along the white line. In the inset (part (b)) we plot  In figure 9, we analyze the dependence of steady concurrence on dipole orientations. In general, higher results are obtained when the two dipoles are parallel. We use again the set of parameters corresponding to the maximum in figure 7(b) along the white line, which is obtained for dipoles along the z-axis. We show how concurrence decreases by changing the dipole directions toward the x-axis (part (a)) always lying on the x z plane and then toward the y-axis (part (b)) always lying on the x y plane (see insets in figures 9(a) and (b)). From part (b) it emerges that aligning the dipoles direction with the interatomic axis (which here is the x-axis) is the optimal choice in the x y plane inducing a lack of symmetry in this plane between the x and y directions.

Dynamics
Here, we discuss the dynamic behavior of two-qubit density matrix elements and concurrence out of thermal equilibrium, also making comparisons with the thermal equilibrium case. This analysis is performed by solving numerically the evolution governed by (41).
In figure 10, we plot several density matrix elements and concurrence as a function of dimensionless time 0 (ω)t. Parts (a) and (b) concern the case of maximally entangled initial states, respectively, the antisymmetric state |A in (a) and the symmetric state |S in (b) (see the values of the parameters in the caption of the figure). Quite different dynamic behavior is pointed out. While starting from |A entanglement is just preserved at a high value, starting from |S concurrence first decreases (going to zero) mainly because of the decrease of ρ S and then revives because of the increase of ρ A . Dynamic creation of entanglement is yet more evident in part (c) where the initial state is the factorized state |2 . In this case, concurrence is initially zero and increases because of the mediated interaction between qubits. Oscillations of C and ρ 23 are linked to the behavior of ρ AS(SA) which rapidly oscillate (see (51)) because of the large value of 12 (ω) which here is equal to 12 (ω)/ 0 (ω) −2.3 × 10 4 . The oscillations are present because ρ AS(SA) are initially populated (ρ AS(SA) = 1/2) while they are not in the cases plotted in parts (a) and (b). States |S and |A are initially equally populated and become different asymptotically. In figure 11, we compare the evolution of concurrence out of thermal equilibrium with the evolutions at equilibrium at the minimal temperature T min = T W = T M = 30 K and at the maximal temperature T max = T W = T M = 1300 K. Two initial maximally entangled configurations are compared, the antisymmetric state |A in part (a) and the symmetric state |S in part (b). At thermal equilibrium, concurrence vanishes on shorter times by increasing the temperature, while out of equilibrium steady entanglement is present. At equilibrium, a larger decay time is observed by starting from the antisymmetric state (see also figure 12 on this subject). In part (b), out of equilibrium, concurrence decays on the same equilibrium timescale, the two-qubit state becoming separable, but it reemerges successively. Both in (a) and (b) a large amount of the initial entanglement is thus asymptotically preserved. In part (c) the initial state is the factorized state |2 . The main difference here is that concurrence presents strong oscillations (see comment on part (c) of figure 10). At thermal equilibrium, entanglement eventually vanishes on a timescale similar to the one of part (a), while out of equilibrium it is maintained after its creation.
In figure 12, we discuss the dependence of super-and sub-radiant effects from the presence/absence of thermal equilibrium. Here, super-and sub-radiance are connected to the occurrence of a decay rate larger or smaller than the one observed in the case of independent qubits, phenomenon due to the interaction of the qubits with a common environment and which depend on the nature of their initial state [26]. In particular, we compare the evolution of the ground state population ρ G starting from |S and |A for two different values of r 12 (0.25 and 15 µm) at thermal equilibrium at T W = T M = 100 K (part (a)) and at T W = T M = 800 K (part (b)) and out of thermal equilibrium for T W = 100 K, T M = 800 K (part (c)). The figure evidences super-radiant behavior when the initial state is |S and sub-radiant when it is |A . The faster or slower increase of ρ G is due to the role of 12 (±ω) which in the two channel decay rates of (49)  figure 4, tend to the same value (+(−)ω)/ 0 (ω), which is the decay rate in the case of single emitters. At thermal equilibrium, the asymptotic state is independent of the values of r 12 while this is not the case out of thermal equilibrium, as pointed out in part (c).
We finally remark that relevant differences are expected when the Markovian and the rotating wave approximation, here adopted, are not valid. In the non-Markovian regime another source of oscillations in the dynamics of concurrence typically emerges [13], while the effect of counter-rotating terms is known to modify the creation of entanglement between the two emitters [60].

Conclusions
In this paper, we have investigated a system made of two quantum emitters interacting with a common stationary electromagnetic field out of thermal equilibrium generated by an arbitrary body and by the surrounding walls held at fixed different temperatures. The environmental field is characterized by means of its correlation functions out of equilibrium which also depend on the scattering properties of the body. We have derived the expressions in the absence of thermal equilibrium of the various functions governing the dissipative dynamics of the two emitters and compared them with the ones holding at thermal equilibrium. This has been done in the case of emitters characterized by an arbitrary number of levels. We have then specialized our investigation to the case of two qubits discussing the new features emerging out of thermal equilibrium.
For a restricted parameter region we have analytically shown that the absence of equilibrium may lead to the generation of steady entangled states. This phenomenon has been interpreted in terms of different effective temperatures associated with two decay channels connecting the total excited and ground states via the symmetric and antisymmetric states, respectively. The two-qubit dynamics can be directed toward mixed states where the antisymmetric contribution is larger than the symmetric one (or vice versa), resulting in the presence of steady entanglement. In this specific case a value of 1/3 has been found as maximum for the concurrence, quantifying the steady entanglement.
We have then numerically investigated the general dependence of steady states and dynamics on the various parameters, without any restriction on the decay rates, in the case the body placed in proximity of the two qubits is a slab made of SiC. The dependence of steady entanglement on the two-qubit distance, their common transition frequency with respect to the slab resonances, the slab thickness, the dipoles orientations and the two involved temperatures have been discussed. Values of concurrence up to 0.24 have been found. Protection and/or generation of entanglement according to the nature of the two-qubit initial state, entangled or not, has been pointed out, also comparing entanglement dynamics in the presence or absence of thermal equilibrium. Higher values of steady concurrence are found for transition frequencies far from the slab resonances (ω/ω r = 0.3) and small thickness (δ 0.01 µm). Remarkably, steady entanglement can be obtained by starting from configurations at thermal equilibrium and by increasing one of the two temperatures involved in the environment of the two qubits.
The possibility to observe the effects we discussed could be explored, for example, for emitters made by trapped atoms [42] or by artificial atoms such as quantum dots or superconducting qubits, placed in proximity to a substrate held at a temperature different from that of the cell surrounding the emitters and the substrate.  which does not depend anymore on the reference system chosen to derive (B.2). In order to compare the previous result with the known expressions, let us consider the case of two qubits in vacuum with dipoles parallel between them (directiond) with different modulus |d q | = |d q |.
In this case, the previous equation reduces to the form (see for instance [25,26]) (B.4) operators, and a reflected one, G (R) ii (R, R , ω), proportional to R. With regards to the imaginary part of G ii (R q , R q , ω) it is possible to check starting from (C.4) that equation (C.3) is verified.