{'fraction_non_alphanumeric': 0.042181130737874405, 'fraction_numerical': 0.013824068056950434, 'mean_word_length': 4.264344580936091, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum thermalization of two coupled two-level systems in eigenstate and bare-state representations

We study analytically the quantum thermalization of two coupled two-level systems (TLSs), which are connected with either two independent heat baths (IHBs) or a common heat bath (CHB). We understand the quantum thermalization in eigenstate and bare-state representations when the coupling between the two TLSs is stronger and weaker than the TLS-bath couplings, respectively. In the IHB case, we find that when the two IHBs have the same temperatures, the two coupled TLSs in eigenstate representation can be thermalized with the same temperature as those of the IHBs. However, in the case of two IHBs at different temperatures, just when the energy detuning between the two TLSs satisfies a special condition, the two coupled TLSs in eigenstate representation can be thermalized with an immediate temperature between those of the two IHBs. In bare-state representation, we find a counterintuitive phenomenon that, under some conditions, the temperature of the TLS connected with the high-temperature bath is lower than that of the other TLS, which is connected with the low-temperature bath. In the CHB case, the coupled TLSs in eigenstate representation can be thermalized with the same temperature as that of the CHB in nonresonant cases. In bare-state representation, the TLS with a larger energy separation can be thermalized to a thermal equilibrium with a lower temperature. In the resonant case, we find a phenomenon of anti-thermalization. We also study the steady-state entanglement between the two TLSs in both the IHB and CHB cases.


I. INTRODUCTION
Conventional quantum thermalization [1,2] is understood as an irreversibly dynamic process under which a quantum system immersed in a heat bath approaches a thermal equilibrium state with the same temperature as that of the bath. The thermal equilibrium state [3] of the thermalized system reads as ρ th (T ) = Z −1 S exp[−H S /(k B T )], which is merely determined by the Hamiltonian H S of the thermalized system and the temperature T of its environment, where Z S = Tr S {exp[−H S /(k B T )]} is the partition function of the thermalized system, with k B being the Bolztmann constant. During the course of a quantum thermalization, all of the initial information of the thermalized system is totally erased by its environment. Recently, much attention has been paid to quantum thermalization (e.g., Refs. [4][5][6][7][8][9][10][11][12]). Specifically, a new kind of thermalization, called canonical thermalization (e.g., Refs. [13][14][15][16][17]), has been proposed.
The conventional quantum thermalization works for the situations wherein one quantum system is connected with just one environment at thermal equilibrium. When we consider a composite quantum system, which is constructed with many subsystems and connected with many environments, the conventional quantum thermalization is no longer valid. Therefore, the density operator of the composite quantum system at thermal equilibrium can not be written as ρ th (T ). In fact, quantum thermalization of a composite quantum system is a very complex problem. On one hand, from the viewpoint of the environments, the composite quantum system could be connected with many independent environments or a common environment. At the same time, the temperatures of these environments could be the same or different. On the other hand, the coupling strengths among these subsystems can affect the physical picture to describe the quantum thermalization of the composite quantum system. When the coupling strengths among these subsystems are stronger than the system-bath couplings, the composite quantum system can be regarded as a single system, while it is regarded as many individual subsystems when the couplings among them are weaker than the system-bath couplings.
The above mentioned situations come true in recent years since quantum systems can be manufactured to be more and more complicated and small, based on the great advances in physics, chemistry, and biology [18]. Therefore, the research on thermodynamics of small systems becomes very interesting. In particular, quantum thermalization of composite quantum systems becomes an important topic since many important results in this field, such as nonequilibrium work relations and fluctuation theorem (e.g., Refs. [19][20][21][22][23]), are based on the thermal equilibrium state of the composite quantum systems. As composite quantum systems are composed of many subsystems, they could be connected either with many independent environments at different temperatures or a common environment. Therefore, it is natural to ask the following questions: (i) How do the couplings among the subsystems affect the quantum thermalization of a composite quantum system? (ii) What are the steady-state properties of a composite quantum system when it is connected with either many independent environments or a common environment?
With these questions, in this paper, we study the quantum thermalization of two coupled two-level systems (TLSs) that are immersed in either two independent heat baths (IHBs) or a common heat bath (CHB). Simple as this model is, it is illustrative. When the coupling between the two TLSs is stronger than the TLS-bath couplings, the two TLSs can be considered as an effective composite system, i.e., a four-level system, and then we understand the quantum thermalization in eigenstate representation of the composite system. In addition, when the coupling between the two TLSs is weaker than the TLS-bath couplings, we understand the quantum thermalization from the viewpoint of each individual TLS. However, due to the TLS-TLS coupling, the effective temperatures of the two TLSs should be different from those for the decoupling case.
As for the environments, there are two kinds of different situations: the IHB case and the CHB case. In the IHB case, we find that, when the two IHBs have the same temperatures, the two coupled TLSs in the eigenstate representation can be thermalized with the same temperature as those of the IHBs. However, in the case where the two IHBs have different temperatures, just when the energy detuning between the two TLSs satisfies a special condition, the two coupled TLSs in eigenstate representation can be thermalized with an immediate temperature between those of the two IHBs. In the bare-state representation, we find a counterintuitive phenomenon that, under some conditions, the temperature of the TLS connected with the high-temperature heat bath is lower than that of the other TLS, which is connected with the lowtemperature heat bath.
In the CHB case, we also study the quantum thermalization in eigenstate and bare-state representations. In the eigenstate representation, the present case reduces to the conventional quantum thermalization, [i.e., one quantum system (an effective four-level system formed by the two coupled TLSs) is thermalized by one environment (the common heat bath) in thermal equilibrium]. In addition, we also investigate the effective temperatures of the two TLSs in bare-state representation. It is found that the TLS with a larger energy separation can be thermalized with a lower temperature. In particular, in the resonant case, we find a quantum phenomenon of antithermalization when the two TLSs are connected with a common heat bath. This paper is organized as follows: In Sec. II, we present the physical models and their Hamiltonians. In Sec. III, we study the quantum thermalization of the two coupled TLSs immersed in two IHBs. In Sec. IV, we consider the case wherein the two TLSs are immersed in a CHB. We conclude this work in Sec. V. Finally, we give two appendices, A and B, for detailed derivations of quantum master equations and transition rates for the IHB and CHB cases, respectively.

II. PHYSICAL MODELS AND HAMILTONIANS
Let us start with introducing the physical models [as illustrated in Figs. 1(a) and 1(b)]: two TLSs, denoted by TLS1 and TLS2 with respective energy separations ω 1 and ω 2 , couple with each other via a dipole-dipole interaction of strength ξ. At the same time, the two TLSs couple inevitably with the  environments surrounding them. Specifically, in this paper, we consider two kinds of different cases: the IHB case and the CHB case. In the former case, the two TLSs are immersed in two IHBs, while, in the latter case, the two TLSs are immersed in a CHB. The Hamiltonian of the total system, including the two TLSs and their environments, is composed of three parts: where H TLSs is the Hamiltonian of the two coupled TLSs, H B is the Hamiltonian of the heat baths, and H I describes the interaction Hamiltonian between the two TLSs and their baths. The Hamiltonian H TLSs (with = 1) reads as The first two terms in Eq. (2) are free Hamiltonians of the two TLSs, which are described by the operators σ + l = (σ − l ) † = |e ll g| and σ z l = |e ll e| − |g ll g| (l = 1, 2), where |g l and |e l are, respectively, the ground and excited states of the lth TLS (i.e., TLSl). The last term in Eq. (2) describes a dipoledipole interaction of strength ξ between the two TLSs. We note that the Hamiltonian given in Eq. (2) has been widely studied in various physical problems, such as quantum logic gates [24], coherent excitation energy transfer [25], decoherence dynamics [26][27][28], and nonequilibrium thermal entanglement [29,30].
Aside from the dipole-dipole interaction between the two TLSs, there exist couplings between the TLSs and their environments. In general, when the couplings of a system with its environment are weak, it is universal to model the environment as a harmonic-oscillator heat bath and choose a linear coupling between the system and its environment [31]. In this paper, we consider this situation and model the environments as harmonic-oscillator heat baths. As mentioned above, we will consider two kinds of different cases: one is the IHB case and the other is the CHB case.
In the IHB case, as shown in Fig. 1(a), the two TLSs are immersed in two IHBs described by the Hamiltonian Here H (a) B and H (b) B are, respectively, the Hamiltonians of the baths for TLS1 and TLS2. The creation and annihilation operators a † j (b † k ) and a j (b k ) describe the jth (kth) harmonic oscillator with frequency ω a j (ω bk ). The interaction Hamiltonian of the two TLSs with the two IHBs reads On the other hand, in the CHB case, as shown in Fig. 1(b), the two TLSs are immersed in a CHB with the Hamiltonian where a † j and a j are, respectively, the creation and annihilation operators of the jth harmonic oscillator with frequency ω a j . The interaction Hamiltonian between the two TLSs and the CHB reads as For simplicity, we have assumed that the TLS-bath coupling strengths g 1 j , g 2 j , and g 2k are real numbers.

III. QUANTUM THERMALIZATION OF TWO COUPLED TLSS IMMERSED IN TWO IHBS
In this section, we study the quantum thermalization of the two coupled TLSs that are immersed in two IHBs. We depict the evolution of the two TLSs in terms of a quantum master equation. By solving the equations of motion of the density matrix elements to obtain steady-state solution, we study the steady-state properties of the two coupled TLSs.

A. Equations of motion and steady-state solutions
We consider the situation wherein the environments of the two TLSs are memory-less and the couplings between the TLSs and the environments are weak. Then, we may adopt the usual Born-Markov approximation in derivation of quantum master equation. At the same time, we derive the master equation in eigenstate representation of the two coupled TLSs so that we can safely make the secular approximation (equivalent to the rotating-wave approximation) to obtain a time-independent quantum master equation by neglecting the rapidly oscillating terms [1]. Therefore, our discussions are under the Born-Markov framework.
In the case of two IHBs, the evolution of the two coupled TLSs is governed by the following Born-Markov master equation in the Schrödinger picture, which will be derived in detail in Appendix A. In Eq. (8), ρ S is the reduced density operator of the two TLSs. The operators τ i j are defined by τ i j = |λ i λ j | in the eigenstate representation of the two coupled TLSs. The summation parameters (i, j) in the second line of Eq. (8) can take (i, j) = (4, 2), (3, 1), (2, 1), and (4,3). In the present model, the effective rates in Eq. (8) are defined as Γ 13 = Γ 24 = Γ 1 , Γ 31 = Γ 42 = Γ 2 , Γ 12 = Γ 34 = Γ 3 , and Γ 21 = Γ 43 = Γ 4 , with where we define The energy separations ε 1 and ε 2 are introduced as respectively, the densities of state at energy ε i of the heat baths for TLS1 and TLS2. In the following we assume γ a (ε 1 In addition, we introduce the average thermal excitation numbersn a ( for the heat baths of TLS1 and TLS2, at the respective temperatures T 1 and T 2 [1]. Based on quantum master equation (8), it is straightforward to obtain optical Bloch equations for the density matrix elements of the two TLSs in the eigenstate representation. Denoting X(t) = ( τ 11 (t) , τ 22 (t) , τ 33 (t) , τ 44 (t) ) T ("T " stands for matrix transpose), then the optical Bloch equations for the diagonal density matrix elements in the eigenstate representation can be expressed aṡ where the coefficient matrix M (IHB) is defined by From Eq. (10), we can see that the evolution of the diagonal density matrix elements decouples with off-diagonal elements.
The transient solutions of optical Bloch equation (10) can be obtained with the Laplace transform method. To study quantum thermalization, however, it is sufficient to obtain the steady-state solutions, , , , where the subscript "ss" means steady-state solutions. According to Eq. (8), we can also obtain the equations of motion for these off-diagonal density matrix elements of the two TLSs as follows: Other off-diagonal elements can be obtained via τ i j (t) * = τ ji (t) . It can be found that the steady-state solutions of the equations of motion for these off-diagonal elements are zero, Based on the steady-state solutions for these density matrix elements, we can analyze the steady-state properties of the two coupled TLSs.

B. Quantum thermalization in eigenstate representation
For the present system with two IHBs, when the coupling between the two TLSs is stronger than the TLS-bath couplings, the two coupled TLSs can be regarded as an effective four-level system connected with two IHBs. Therefore, in the eigenstate representation, the dynamic evolution process of the two-coupled TLSs approaching their steady state of thermal equilibrium can be understood as a non-equilibrium quantum thermalization: thermalization of a quantum system connected with many IHBs at different temperatures [32][33][34][35][36][37][38][39].
As the steady state of the two coupled TLSs is completely mixed in eigenstate representation, we can introduce effective temperatures to characterize the relation between any two eigenstates based on their steady-state populations. From Eq. (12) we can find the following relations Generally, it is impossible to define an effective temperature for the effective four-level system at steady state. According to Eq. (15), we can characterize the state of the effective four-level system via introducing the following two effective temperatures [40] by , In terms of Eq. (15), we get When the two IHBs have the same temperatures, namely T 1 = T 2 = T , we can show that the above two effective temperatures are equal to those of the two IHBs, that is For the nonequilibrium case of T 1 T 2 , the two effective temperatures are different for general case. We find the following relations min(T 1 , , which mean that the two effective temperatures T eff (ε 1 ) and T eff (ε 2 ) will be within the region from min(T 1 , T 2 ) to max(T 1 , T 2 ) [11]. Under some special conditions, the two effective temperatures could be equal. In this case, we consider that the effective four-level system is thermalized by the non-equilibrium environments: two IHBs at different temperatures.
In Fig. 2, we plot the two effective temperatures T eff (ε 1 ) and T eff (ε 2 ) as a function of the mixing angle θ. It shows that T eff (ε 1 ) and T eff (ε 2 ) can be equal for a special mixing angle θ 0 , which is determined by ε 1 ln(Γ 1 /Γ 2 ) = ε 2 ln(Γ 3 /Γ 4 ) . That is to say the effective four-level system formed by the two coupled TLSs can be thermalized with an effective temperature when ∆ω = ξ/ tan θ 0 . Note that this effective temperature actually is a nonequilibrium temperature [41]. We can also see from Fig. 2 that the two effective temperatures T eff (ε 1 ) and T eff (ε 2 ) are within the region from T 1 to T 2 . We point out that the mixing angle θ should be chosen to make sure that E λ 1 > E λ 2 .

C. Quantum thermalization in bare-state representation
When the coupling between the two TLSs is weaker than the TLS-bath couplings, we understand the quantum thermalization in bare-state representation. We can express the bare states with the eigenstates as Then we can obtain the relations According to Eqs. (12), (14), and (20), the steady-state solution can be obtained as In addition, the off-diagonal elements of the density matrices of the two TLSs can be expressed as Because of τ i j (t) ss = 0 (i j), we have σ + 1 ss = σ + 2 ss = 0, which implies that the steady states of the two TLSs in barestate representation are completely mixed. Based on these, it is possible to introduce two effective temperatures of the two TLSs as follows: In Fig. 3, we plot the effective temperatures T eff (ω 1 ) and T eff (ω 2 ) as a function of the mixing angle θ for various temperature distributions. From Fig. 3, we can see the following three interesting results. Firstly, when θ ≈ π/2, T eff (ω 1 ) and T eff (ω 2 ) become approximately equal. Specially, at the resonant point θ = π/2, the two TLSs have the same temperatures no matter whether the two bath temperatures are the same or not (see the cross point at θ = π/2 in figures). This result can be explained as follows: When the two TLSs are nearly in resonance with each other, the dipole-dipole interaction can induce population exchange between the two TLSs such that their temperatures are approximately equal. When θ = π/2 (ω 1 = ω 2 ), Eq. (21) reduces to then we have T eff (ω 1 ) = T eff (ω 2 ). Secondly, when θ is near to 0 and π, T eff (ω 1 ) and T eff (ω 2 ) approach approximately T 1 and T 2 , respectively. This result can be understood as follows: Near to θ ≈ 0 and π, the energy detuning |∆ω| = 2ξ/| tan θ| between the two TLSs is very large, and the population exchange between them can be neglected for a weak coupling strength ξ. At the same time, the energy separation shifts of the two TLSs are neglectable [this can be seen from the following effective Hamiltonian (25)]. Hence, the temperatures of the two TLSs should be equal to those of their IHBs, respectively.
Thirdly, in Fig. 3, there are some regions of θ where the effective temperature T eff (ω 1 ) of TLS1 could be smaller than the temperature T eff (ω 2 ) of TLS2 although the bath temperatures T 1 > T 2 . This is a counterintuitive result. Intuitively, in these cases T eff (ω 2 ) should be smaller than T eff (ω 1 ), because TLS2 is connected with a low-temperature bath while TLS1 is connected with a high-temperature bath [In Fig. 3(a), when the two bath temperatures are the same, the effective temperatures of the two TLSs should be equal]. Actually, the counterintuitive result is a net effect of the energy separations of the two TLSs, the bath temperatures, and the coupling between the two TLSs (coupling induces population exchange and energy shifts). From Figs. 3(a) to 3(d), the curves change gradually with the increase of the temperature difference T 1 −T 2 between the two baths. And the counterintuitive region decreases with the increase of the temperature difference T 1 − T 2 .
In the following, we present a microscopic explanation for Fig. 3 as a physical insight of the counterintuitive phenomenon. When the values of θ are far from the resonant point (not near to 0 and π), the two TLSs will have large detuning from each other. Under the large detuning condition ξ ≪ ∆ω, the real population exchange between the two TLSs is compressed, the dipole-dipole interaction induces energy shift to the two TLSs. In this case, we can derive an effective Hamiltonian to describe the two TLSs with the Fröhlich-Nakajima transformation approach [42,43]. Starting from the Hamiltonian Then the effective Hamiltonian reads as where the shifted energy separations are defined byω 1 = ω 1 + ξ 2 /∆ω andω 2 = ω 2 − ξ 2 /∆ω. We can see from Eq. (25) that, under the large detuning condition, there is no effective coupling between the two TLSs. The dipole-dipole interaction between the two TLSs shifts their energy separations slightly. Hence, when the two TLSs (with the shifted energy separation) are thermalized to thermal equilibrium with their baths, we have the relation exp(−ω l /T l ) = p (l) e /p (l) g (l = 1, 2) for the TLSl (with shifted energy separationω l , excited and ground state populations p (l) e and p (l) g ) in thermal equilibrium at temperature T l , and then the effective temperatures defined in Eq. (23) should be From Eq. (25), we can see that, for a positive ξ, when 0 < θ < π/2, we have ∆ω > 0, thenω 1 > ω 1 andω 2 < ω 2 . Hence the effective temperatures T eff (ω 1 ) < T 1 and T eff (ω 2 ) > T 2 .
On the other hand, when π/2 < θ < π, we have ∆ω < 0, thenω 1 < ω 1 andω 2 > ω 2 . Hence, the effective temperatures T eff (ω 1 ) > T 1 and T eff (ω 2 ) < T 2 . According to the above analysis, we can see that, when T 1 = T 2 [ Fig. 3(a)], there will exist a counterintuitive region. At the same time, when θ is near to 0 and π, the shifted energy separation ξ 2 /∆ω approaches zero, then T eff (ω l ) ≈ T l . Hence, with the increase of the bath temperature difference T 1 − T 2 , the difference between the two effective temperatures also increase, which leads to the counterintuitive region decreases. These results can be seen from Fig. 3.
In fact, the above intuitive result is based on the phenomenological master equatioṅ with The superoperator L l [ρ S ] describes the dissipation of a TLSl (l = 1, 2) immersed in a heat bath at temperature T l (n l = 1/[exp(ω l /T l ) − 1]). Therefore, Eq. (27) is not valid in the case of two coupled TLSs, especially when the coupling between the two TLSs is stronger than the TLS-bath couplings. In addition, in Eq. (27), the effects on the TLSs from the two baths are different, one is direct and the other is indirect. For example, the bath of TLS1 affects TLS1 directly, while the bath of TLS2 affects TLS1 indirectly through TLS2. On the contrary, our results given in Eq. (23) are based on quantum master equation (8), which is rigorously derived in the eigen-representation of the two coupled TLSs. Hence, the dissipation is depicted in the eigen-representation of the two coupled TLSs. In other words, the two TLSs play equivalent roles and the two baths directly affect the TLSs. The resonant case is a clear example for the equivalent role of the two TLSs. In the resonant case θ = π/2, we obtain T eff (ω 1 ) = T eff (ω 2 ) in terms of quantum master equation (8), while we get T eff (ω 1 ) > T eff (ω 2 ) from Eq. (27) when T 1 > T 2 .

D. Steady-state entanglement between the two TLSs
In the IHB case, there exists a dipole-dipole interaction between the two TLSs. Therefore, a natural question is: what the quantum entanglement is between the two TLSs after they are thermalized. As we know, during the thermalization processes (not antithermalization), all of the initial information of the two coupled TLSs is totally erased, and the steady state of the two TLSs is determined by the decay rates and bath temperatures. Hence, we need to know the steady-state entanglement in the two TLSs. We note that entanglement dynamics in similar systems has been studied [44]. In the following we apply the concurrence to quantify the steady-state entanglement in the two TLSs.
In the following, we study the steady-state concurrence of the two coupled TLSs in the case of θ = π/2 and T 1 = T 2 = T . In Fig. 4, we plot the steady-state concurrence C(ρ ss ) as a function of the bath temperature T for various values of the dipole-dipole interaction strength ξ. Figure 4 shows that a larger steady-state concurrence can be created for a larger ξ.
We also see the sudden death of the concurrence when the bath temperature increases up to a critical value. Notice that the phenomenon of threshold temperature has also been found in thermal entanglement of spin model [29,30,47,48].

IV. QUANTUM THERMALIZATION OF TWO COUPLED TLSS IMMERSED IN A CHB
In the above section, we have studied the quantum thermalization of two coupled TLSs immersed in two IHBs. However, in some cases, the two coupled TLSs can be considered to be placed in a CHB. In this section, we shall study the quantum thermalization of two coupled TLSs immersed in a CHB.
Correspondingly, the optical Bloch equation for the CHB case has the same form as Eq. (10), but the coefficient matrix is replaced by M (CHB) with the following expression We can also obtain the equations of motion for these offdiagonal density matrix elements as follows: The equations of motion for other elements can be obtained by τ i j (t) = τ * ji (t) . Obviously, the steady-state solutions of these off-diagonal density matrix elements are zero.

B. Quantum thermalization in eigenstate representation
Differently from the IHB case, for the present four-level system, we introduce three effective temperatures to characterize its state. The three effective temperatures T 12 , T 13 , and T 34 are defined according to the populations of the four levels as follows: We can show that the above introduced temperatures are the same as that of the CHB, which means the two coupled TLSs can approach a thermal equilibrium with the CHB. In other words, the two coupled TLSs in eigenstate representation can be thermalized by the CHB. According to Eqs. (18) and (42), we know that when the temperatures of the heat baths are T , irrespective of two IHBs or a CHB, the effective four-level system formed by the two coupled TLSs can be thermalized into a thermal equilibrium state with the same temperature T . In other words, based on the thermal equilibrium state at temperature T , we can not know whether the two coupled TLSs are connected with two IHBs or a CHB.

C. Quantum thermalization in bare-state representation
We also investigate the quantum thermalization of the two TLSs in the bare-state representation. In terms of Eqs. (20) and (39), we can obtain the steady-state average values of the two Pauli operators σ z 1 and σ z 2 as follows: Moreover, we have σ + 1 ss = σ + 2 ss = 0. Similar to Eq. (23) in the above section, we also introduce two effective temperatures to characterize the state of the two TLSs. In Fig. 5, we plot the two effective temperatures as a function of the mixing angle θ. We emphasize that the resonant point θ = π/2 in Fig. 5 should be taken out.
We can draw a conclusion from Fig. 5 that the TLS with a larger energy separation can be thermalized to a thermal equilibrium state with a lower temperature. This can be seen as follows: When ω 1 > ω 2 , we have 0 < θ < π/2, from Fig. 5 it is clear that the temperature of TLS1 is lower than that of TLS2. On the other hand, when ω 1 < ω 2 , we have π/2 < θ < π, Fig. 5 indicates that the temperature of TLS2 is lower than that of TLS1 in this region. The physical explanation for Fig. 5 is the same as that for Fig. 3(a).

D. Quantum anti-thermalization in the resonant case
When θ = π/2, the two TLSs are in resonance, and then the decay rates Γ 13 , Γ 31 , Γ 34 , and Γ 43 are zero under the assumption γ 1 (ε i ) = γ 2 (ε i ). Hence, the eigenstate |λ 3 decouples with other eigenstates, resulting in an anti-thermalization phenomenon. The state |λ 3 is called as a "dark state" [49]. In this case, we need to rewrite new optical Bloch equations for these density matrix elements. We obtain the equations of motion for diagonal density matrix elements τ 11 (t) = −2Γ 12 τ 11 (t) + 2Γ 21 τ 22 (t) , τ 22 (t) = 2Γ 12 τ 11 (t) − 2(Γ 21 + Γ 24 ) τ 22 (t) +2Γ 42 τ 44 (t) τ 33 and for off-diagonal density matrix elements The steady-state solutions for these density matrix elements are where τ 33 (0) = 0 is the initial population of state |λ 3 . It is obvious that the steady state of the two coupled TLSs depends on its initial state. Therefore, the two coupled TLSs exhibit a phenomenon of anti-thermalization in the sense that the heat bath can not erase totally the initial information of the two TLSs. For example, when initially the two coupled TLSs are prepared in state |λ 3 , they will stay in |λ 3 forever. However, in the subspace spanned by the three eigenstates |λ 1 , |λ 2 , and |λ 4 , the initial information of the two coupled TLSs can be totally erased, which can also be seen from Eq. (46 ) when τ 33 (0) = 0.

E. Steady-state entanglement between the two TLSs
In the CHB case, in addition to the dipole-dipole interaction between the two TLSs, the common bath can also provide a physical mechanism to entangle the two TLSs. According to Eqs. (33) and (34), the steady-state density matrix elements of the two TLSs in the CHB case have the same form as those given in Eq. (35). However, now the steady-state solutions of the eigenstate populations τ j j ss are given by Eq. (39). Correspondingly, we can obtain the concurrence between the two TLSs in terms of Eq. (36). For a given nonresonant θ, the figure in the CHB case is very similar to Fig. 4 (so it is not shown here). The phenomenon of threshold temperature also exists in the CHB case.

V. CONCLUSION AND DISCUSSIONS
In conclusion, we have studied the quantum thermalization of two coupled TLSs which are immersed in either two IHBs or a CHB. We have characterized the temperatures of the two coupled TLSs in eigenstate and bare-state representations when the coupling between the two TLSs is stronger and weaker than the TLS-bath couplings, respectively. In the IHB case, we have found that, when the two IHBs have the same temperatures, the two coupled TLSs could be thermalized in eigenstate representation with the same temperature as those of the heat baths. However, in the case where the two heat baths have different temperatures, just when the energy detuning between the two TLSs satisfies a special condition, the effective four-level system formed by the two coupled TLSs can be thermalized with an immediate temperature between those of the two heat baths. In bare-state representation, we have found a counterintuitive phenomenon that the temperature of the TLS connected with the high-temperature heat bath is lower than that of the other TLS which is connected with the low-temperature heat bath. In the CHB case, the two TLSs in eigenstate representation could be thermalized with the same temperature as that of the heat bath for nonresonant cases. In bare state representation, we have found that the TLS with a larger energy separation can be thermalized to a thermal equilibrium at a lower temperature. We have also found a phenomenon of anti-thermalization of the two TLSs in a common heat bath in the resonant case. In addition, we have studied the steady-state entanglement of the two TLSs in the IHB and CHB cases. It has been found that there exist threshold temperatures for the steady-state entanglement generation.
Finally, we present some discussions on the thermalization time over which the thermalized systems evolve from their initial states to steady states. Mathematically, the thermalization time for a system should be infinite because the long-time limit (lim t→∞ ) is needed to make sure that these density matrix elements evolve to their steady-state values. From the viewpoint of physics, we might introduce some time scales to describe a thermalization, as the half-life of an exponential decay. However, for present systems, the evolutions of these density matrix elements [i.e., the transient solution of Eq. (10)] are not purely exponential functions. At the same time, these evolutions depend on the initial conditions. Therefore, it is needed to introduce the time scales under given initial conditions, other than a universal time scale for a quantum thermalization.

Appendix A: Derivation of quantum master equation (8) for the IHB case
In this appendix, we give a detailed derivation of quantum master equation (8), which describes the evolution of the two TLSs immersed in two IHBs. In the interaction picture with respect to H 0 = H TLSs + H (IHB) B , the interacting Hamiltonian (5) becomes where τ i j = |λ i λ j | and we introduce the noise operators with A(t) = j g 1 j a j e −iω a j t and B(t) = k g 2k b k e −iω bk t . Under the Born-Markov approximation [1], the quantum master equation readsρ where Tr B stands for tracing over the degrees of freedom of the baths. We assume that the two baths are in thermal equilibrium B )] are the partition functions of the two baths, respectively. The parameters β 1 = 1/T 1 and β 2 = 1/T 2 are the inverse temperatures of the baths for TLS1 and TLS2. Through making the rotating wave approximation, we obtaiṅ and where the parameters Γ i and Λ i have been defined in Eq. (9). Based on the above one-side Fourier transforms of these correlation functions, we can obtain those for other correlation functions. By substituting them into quantum master equation (A4) and returning to the Schrödinger picture, we can obtain quantum master equation (8).
With the same method, the one-side Fourier transform for other correlation functions can also be obtained.