{'fraction_non_alphanumeric': 0.04127888572333017, 'fraction_numerical': 0.010794555238999684, 'mean_word_length': 4.277480788506515, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum thermalization and thermal entanglement in the open quantum Rabi model

We study quantum thermalization and thermal entanglement in the open quantum Rabi model (QRM), in which a two-level system and a single-mode bosonic field are coupled to either two individual heat baths or a common heat bath. By treating the QRM as an effective multilevel system and deriving global quantum master equations in the eigenstate representation of the QRM, we study the physical conditions for quantum thermalization of the open QRM. It is found that, in the individual heat-bath case, the QRM can only be thermalized when either the two heat baths have the same temperature or the QRM is only coupled to one of the two baths. In the common heat-bath case, differently, the QRM can always be thermalized. We also study thermal entanglement of the QRM in both the resonant- and non-resonant coupling cases. The logarithmic negativity for the thermal state of the QRM is obtained in a wide parameter space, ranging from the low- to high-temperature limits, and from the weak- to deep-strong-coupling regimes. This work paves the way towards the study of quantum effects in nonequilibrium ultrastrongly-coupled light-matter systems.


I. INTRODUCTION
Quantum thermalization [1,2], as one of the most important topics in quantum thermodynamics, is understood as an irreversible dynamic process via which a quantum system immersed in its environment approaches a thermal equilibrium state at the same temperature as the environment. Until now, considerable studies have been done on both quantum thermalization of coupled quantum systems and thermal entanglement [3] between the subsystems in coupled systems. For instance, quantum thermalization and thermal entanglement of two coupled two-level atoms have been studied in the dressedstate representation [4]. Recently, quantum thermalization and thermal entanglement in the open Jaynes-Cummings model (JCM) have also been studied [5]. In addition, thermal entanglement has been studied in various physical systems, including coupled spins [6][7][8][9][10] and coupled oscillators [11].
In general, the thermal state of a coupled quantum system should be entangled because the thermal state takes the form as ρ th = Z −1 sys exp(−βH sys ) with H sys = H A + H B + H I , where H A and H B are, respectively, the free Hamiltonians of the subsystems A and B, H I is the interaction Hamiltonian between A and B, and Z sys = Tr[exp(−βH sys )] is the partition function with β = 1/(k B T ) being the inverse temperature (k B is the Boltzmann constant). However, a counterintuitive phenomenon of vanishing thermal entanglement in the open JCM [5] has been found and proved. Usually, the JCM [12] is obtained by discarding the counter-rotating terms in the quantum Rabi model (QRM) [13] with the rotating-wave approximation (RWA), when the interaction between a two-level system (TLS) and a single-mode bosonic field does not enter the ultrastrong-coupling regime [14]. Therefore, a natural and important question is what happens with the thermalization and thermal entanglement of the open QRM in the ultrastrong-coupling regime. In particular, great advances have been made in the enhancement of the coupling strength of the QRM in the last decade, and the ultrastrong couplings [15][16][17] even deep-strong couplings [18] in the QRM have been realized in various physical systems, such as superconducting quantum circuits [17,[19][20][21][22] and semiconductor quantum wells [23][24][25]. Therefore, the study of quantum thermalization and thermal entanglement in the open QRM becomes an urgent topic to be addressed.
In this paper, we study quantum thermalization of the open QRM, which is coupled to either two individual heat baths (IHBs) or one common heat bath (CHB). Concretely, we derive global quantum master equations [26][27][28][29][30][31] to describe the evolution of the QRM. The global quantum master equations are valid even when the interaction between the TLS and the bosonic mode is much stronger than the system-bath couplings. Note that quantum statistics based on the global quantum master equations have been studied in coupled atom systems [32][33][34][35][36][37][38][39][40][41][42][43][44][45] and coupled harmonic-oscillator systems [46][47][48]. To characterize the thermalization of the QRM, we introduce the effective temperatures associated with any two eigenstates, and evaluate the thermalization of the QRM by inspecting whether its steady-state density matrix can be expressed as a thermal equilibrium state. Furthermore, we study thermal entanglement [6] between the two subsystems of the QRM by calculating the logarithmic negativity of the thermal state. We obtain the thermal entanglement of the QRM in both the resonant-and nonresonant-coupling cases. Concretely, we find the dependence of the logarithmic negativity on the coupling strength and the temperature. We also explain the thermal entanglement in the low-temperature limit based on the eigenstate entanglement. The rest of this paper is organized as follows. In Sec. II, we introduce the open QRM in both the resonant-and nonresonant-coupling cases, and present the Hamiltonians of the system. In Secs. III and IV, we study quantum thermalization of the QRM in both the IHB and CHB cases, respectively. In Sec. V, the thermal entanglement of the QRM is analyzed in both the resonantand nonresonant-coupling cases. A conclusion is given in Sec. VI.

II. MODEL AND HAMILTONIAN
We consider a quantum Rabi model which describes the interaction of a TLS with a single-mode bosonic field through the dipole coupling (see Fig. 1). The two levels of the TLS are denoted as the ground state |g and the excited state |e , and the energy separation between these two levels is ω 0 . We introduce the Pauli operators σ x = |e g| + |g e|, σ y = i(|g e| − |e g|), and σ z = |e e| − |g g| to describe the TLS. For the singlemode bosonic field, we assume its resonance frequency as ω c and denote its annihilation (creation) operator as a (a † ). The Hamiltonian of the QRM takes the form as Here, the first and second terms describe the free Hamiltonian of the TLS and the single-mode bosonic field, respectively. The last term describes the interaction between the TLS and the field, with g being the coupling strength. When g/ω 0 , ω c > 0.1, the QRM enters the ultrastrong-coupling regime. The QRM is in the deepstrong-coupling regime when g/ω 0 , ω c > 1. For the QRM, the total excitation operator N = a † a + |e e| is no longer a conserved quantity. However, there exists a parity operator P = −σ z (−1) a † a , which commutes with the Rabi Hamiltonian H QRM , and hence we can divide the whole Hilbert space of the system into odd and even parity subspaces [49,50]. Generally, it is hard to analytically solve the eigensystem of the QRM with the element functions. However, the eigenenergy spectrum of the QRM has been obtained with several methods under the assistant of numerical calculations [13,[51][52][53]. In this work, we will numerically solve the eigensystem of the QRM and study the steady-state populations of these eigenstates in a sufficiently large Hilbert space. For expressional convenience, we denote the eigensystem of the Rabi Hamiltonian as H QRM |ε m = ε m |ε m (m = 1, 2, 3, ..., ∞), where |ε m are eigenstates with the corresponding eigenvalues ε m .
In Fig. 2, we show the lowest eight eigenstate energies of the QRM as a function of the coupling strength in both resonant and nonresonant-coupling cases. We can see several features in the energy spectrum. (i) The entire pattern of these levels decreases with the increase of the coupling strength g/ω 0 . For a larger value of the ratio ω c /ω 0 , the slope is smaller. (ii) These energy levels can be paired off from lower states to upper states. Except for the first pair of states, the two states in other pairs intersect with each other, and the crossing levels have different parities. With the increase of the ratio ω c /ω 0 , the horizontal distance between two neighboring crossing points increases. (iii) In the absence of the coupling, i.e., g = 0, these states are reduced to the bare states, and hence the states with the same excitation are degenerate in the resonant case [i.e., panel (b)]. (iv) When the system enters the deep-strong coupling regime, the states in the same pair become near degenerate. Note that the information of the energy spectrum will be useful for analyzing the thermal entanglement.
To characterize the interactions of both the TLS and the bosonic mode with the environments, we consider two different cases: (i) The TLS and the bosonic mode are coupled to individual heat bath alone [ Fig. 1(a)]. (ii) Both the TLS and the bosonic mode are coupled to a common heat bath [ Fig. 1(b)]. In the IHB case, the free Hamiltonians of the two heat baths are expressed as H where the creation and annihilation operators C † q (A † k ) and C q (A k ) describe the qth (kth) mode with frequency ω q (ω k ) in the heat bath of the TLS (bosonic mode). The interaction Hamiltonian between the QRM and the two heat baths can be expressed as (2) where λ q (η k ) is the coupling strength between the TLS (bosonic mode) and the qth (kth) mode of its heat bath.
In the CHB case, both the TLS and the bosonic mode are connected with a common heat bath described by the Hamiltonian H where B † p and B p are, respectively, the creation and annihilation operators of the pth mode with resonance frequency ω p in the common bath. The interaction Hamiltonian between the QRM and the CHB reads where λ p (η p ) is the coupling strength between the TLS (bosonic mode) and the pth mode of the CHB.
In both the IHB and CHB cases, the Hamiltonian of the whole system including the QRM and the heat baths can be written as H = H QRM +H are, respectively, the bath Hamiltonians and the interaction Hamiltonians between the QRM and its baths, with s = IHB and CHB corresponding to the individual heat-bath and common heat-bath cases, respectively.

III. QUANTUM THERMALIZATION OF THE QRM IN THE IHB CASE
In this section, we study the quantum thermalization of the QRM in the IHB case. We derive a global quantum master equation to govern the evolution of the open QRM. We also investigate the effective temperature associated with any two eigenstates to evaluate the quantum thermalization of the QRM.

A. Global quantum master equation
To include the dissipation in this system, we derive the quantum master equation within the Born-Markov framework, which is valid under the assumption of weak system-bath coupling and short bath correlation time.
In particular, the quantum master equation is derived in the eigenstate representation of the Rabi Hamiltonian. In the interaction picture with respect to , the evolution of the QRM is governed by the global quantum master equation where the Lindblad superoperators D[|ε n ε m |]ρ S (t) and D[|ε m ε n |]ρ S (t) are defined by for o = |ε n ε m | or o † = |ε m ε n |. In Eq. (4), the decay rates related to the dissipation channels of the TLS and the bosonic mode are, respectively, defined by where ̺ σ (ω q ) and ̺ a (ω k ) are, respectively, the spectral density functions of the heat baths associated with the TLS and the bosonic mode, and ω m,n = ε m −ε n denotes the energy separation between the two eigenstates |ε m and |ε n of the QRM. In our calculations, we suppose that the decay rates γ σ (ω m,n ) = γ σ and γ a (ω m,n ) = γ a , which means that the decay rates related to all the transitions caused by the same subsystem are identical. The transition coefficients in Eq. (4) are defined by χ σ,m,n = ε m |σ x |ε n and χ a,m,n = ε m |(a + a † )|ε n . In addition, the average thermal-excitation numbers in Eq. (4) are defined by n l=σ,a (ω m, where ω m,n denotes the energy separation between the involved two eigenstates |ε m and |ε n . The parameters T σ and T a are, respectively, the temperatures of the heat baths connected to the TLS and the bosonic mode. In terms of the transformation ρ S (t) = e −iHQRMt/ ρ S (t)e iHQRMt/ and quantum master Eq. (4), we can obtain the quantum master equation in the Schrödinger picture as where ρ S (t) is the reduced density matrix of the QRM in the Schrödinger picture, and the dissipator L IHB [ρ S (t)] is obtained by substitutingρ S (t) with ρ S (t) in Eq. (4). We should point out that the global quantum master equation (7) does not work at the degenerate points in the energy spectrum. This is because we have used the secular approximation in the derivation of the Lindblad dissipators by discarding the related crossing terms. Theoretically, these crossing terms should be kept at the degenerate points. For keeping the uniformity of the systematic description of the open Rabi system, we derive the global quantum master equation in the full parameter space, and add this notice to avoid these degenerate points. Note that this notice also works for the open Rabi model in the CHB case.

B. Quantum thermalization
In the non-equilibrium open-system case, the QRM is connected with two heat baths, which could be at different temperatures. To evaluate the thermalization, we check whether the steady-state density matrix of the QRM can be expressed as a thermal state } is the partition function of the QRM, with T being the temperature of the thermalized system. In the quantum thermalization process, the environment erases all the initial-state information of the thermalized system. For the QRM, it is valid when the coupling enters the ultrastrong even deep-strong coupling regimes. Therefore, we should treat the QRM as an effective multiple-level system, i.e., working in the eigenstate representation of the QRM. As a result, the thermalization of the QRM in the IHB case can be understood as the thermalization of a multiple-level system connected with two heat baths, which would be at either different temperatures or the same temperature. By solving the global quantum master equation (7) in the eigenstate representation of the QRM, we find that the steady state of the QRM is a completely mixed state in this representation. Motivated by this feature, we introduce effective temperatures associated with any two eigenstates based on their populations. If all the effective temperatures are the same, then the steady-state density matrix of the QRM can be written as a thermal state. In this case, we regard it as the quantum thermalization of the QRM.
We denote the populations of the states |ε m (|ε n ) as p m (p n ). Then, we define the effective temperature as According to these effective temperatures, we can characterize the thermalization of the QRM. Based on Eq. (7), we can obtain the equations of motion for the density matrix elements ε n |ρ ss |ε m . For investigating the thermalization, we focus on the steady state of the system, which can be solved by setting d dt ρ S (t) → 0, then the steady-state density matrix elements ε n |ρ ss |ε m can be obtained. In terms of the population p n = ε n |ρ ss |ε n of the eigenstate |ε n , then the effective temperatures T eff (ω m,n ) can be calculated accordingly. In our numerical simulations, we need to truncate the dimension up to n d of the Hilbert space of the bosonic field so that n d n=1 ε n |ρ ss |ε n ≈ 1. In Fig. 3, we plot the effective temperatures k B T eff (ω m,n )/ ω 0 as functions of the energy-level indexes m and n when the scaled bath temperatures k B T σ / ω 0 and k B T a / ω 0 take various values. Figures 3(a) and 3 Figure 3(a) shows that the QRM can be thermalized under the circumstance of k B T σ / ω 0 = k B T a / ω 0 . In this case, the temperature of the thermalized QRM equals to that of the two baths, and the density matrix of the QRM can be written as the thermal state. In particular, the thermalization of the QRM is independent of the nonzero values of the decay rates. In Figs. 3(b)-3(d), the temperatures of these adjacent energy levels of the QRM are not the same, then the steady state of the QRM cannot be expressed as a thermal state, and hence the QRM cannot be thermalized in this case. It should be emphasized that the system-bath coupling configuration is of great importance for the thermalization of the QRM. In the case of unequal temperatures of the two baths, the QRM cannot be thermalized. When one of the system-bath couplings is turned off, the thermalization of the QRM will have different results. In Fig. 4, we show the effective temperatures k B T eff (ω m,n )/ ω 0 as functions of the energy-level indexes m and n. Figures 4(a) and 4(b) correspond to the cases where the bath of the bosonic mode is decoupled [ Fig. 4(a)] and the bath of the bosonic mode is at zero temperature [ Fig. 4(b)], respectively. Similar results can also be found when either the TLS decouples from its bath or the bath of the TLS does not excite the system [Figs. 4(c) and 4(d)]. These results indicate that the QRM can be thermalized when one of the two baths is decoupled from the QRM. When the couplings of the QRM with the two baths exist simultaneously, the QRM cannot be thermalized even when one of the baths is at zero temperature (the temperature of the other bath is not zero). Consequently, it can be concluded that the QRM can be thermalized when only one heat bath coupled to the QRM and cannot be thermalized when the two heat baths are at two different temperatures.

IV. QUANTUM THERMALIZATION OF THE QRM IN THE CHB CASE
In this section, we investigate whether the QRM can be thermalized in the common heat-bath case. Concretely, we derive a global quantum master equation to govern the evolution of the QRM as with the dissipator L CHB [ρ S (t)] given by where the Lindblad superoperator D[o] is defined in Sec. III. We define the decay rates in Eq. (10) as which correspond to the dissipations through the TLS, the bosonic mode, and the cross effect between the two subsystems, respectively. Similar to the IHB case, we obtain the steady state of quantum master equation (9), and find that the density matrix of the QRM in the eigenstate representation is diagonal. Based on the steady-state populations of these eigenstates |ε n , we calculate the effective temperatures assosiated with any two eigenstates |ε m and |ε n . We also checked the thermalization of the QRM in the CHB case when the two decay rates take different values. It has been found that the thermalization in the CHB case is independent of the values of the two decay rates (the two decay rates cannot be zero at the same time). Here, we do not present the figure because the results are the same as Fig. 3(a).

V. THERMAL ENTANGLEMENT IN THE QRM
In this section, we study thermal entanglement in the QRM. When the QRM is thermalized at temperature T , its density matrix can be written as ρ th (T ) = Z −1 QRM exp(−βH QRM ) with β = 1/(k B T ). We study quantum entanglement between the TLS and the singlemode bosonic field by calculating the logarithmic negativity of the thermal state. Note that here we use the logarithmic negativity to quantify the quantum entanglement between the TLS and the bosonic mode because the thermal state is a mixed state. We should point out that, in the zero-temperature case, i.e., T = 0, the thermal state is reduced to the ground state of the QRM. In some previous papers [54][55][56], quantum entanglement in the ground state of the QRM has been studied using the von-Neumann entropy. It has been found that the entropy of entanglement monotonically increases and exhibits a saturation effect with the increase of the coupling strength g/ω 0 at the resonant case ω c = ω 0 . We remind that the logarithmic negativity does not reduce to the entropy of entanglement, and that the logarithmic negativity could be zero even the state is entangled. However, for quantifying quantum entanglement of the QRM in thermal state, the logarithmic negativity is a computable quantity.
For the QRM in the thermal state ρ th (T ), its logarithmic negativity can be calculated by where "T a " denotes partial transpose with respect to the bosonic mode. The trace norm is defined by ||ρ Ta th (T )|| 1 = Tr (ρ Ta th (T )) † ρ Ta th (T ) . According to the relation H QRM |ε k = ε k |ε k , the thermal state of the QRM can be written as ρ th (T ) = Z −1 QRM ∞ k=1 exp(−βε k )|ε k ε k |. For calculation of the logarithmic negativity, we need to express the density matrix ρ th (T ) with the bare states |e, m and |g, n . Therefore, we expand the eigenstates of the QRM with the bare basis vectors |e, m and |g, n as |ε k = ( ∞ m=0 |e, m e, m| + ∞ n=0 |g, n g, n|)|ε k . Before discussing the thermal entanglement, we first analyze the parameter space of the thermal state, such that the thermal entanglement can be discussed clearly in the full parameter space. For the QRM in the thermal state, its density matrix can be expressed as It can be seen from Eq. (13) that the density matrix ρ th (T ) is determined by three ratios: (i) The scaled inversed temperature ω 0 β, which can be used to characterize the temperature of the thermalized system. When both the coupling strength g and the resonance frequency ω c are either smaller than or of the same order of ω 0 , then the ω 0 can be used to characterize the energy scale of the system, and the ratio ω 0 β can be used to describe the temperature scale of the system. The relations ω 0 β ≫ 1 and ω 0 β ≪ 1 stand for the low-and hightemperature limits, respectively. (ii) The frequency ratio ω c /ω 0 , which describes the frequency mismatch between the TLS and the bosonic mode. Here ω c /ω 0 = 1 and ω c /ω 0 = 1 correspond to the resonant-and nonresonantcouplings between the TLS and the bosonic mode, respectively. (iii) The scaled coupling strength g/ω 0 , which is the ratio of the coupling strength over the resonance frequency of the TLS. Based on the appoint in ultrastrong coupling, g/ω 0 > 0.1 and g/ω 0 > 1 are considered as the conditions for characterizing the ultrastrongcoupling and deep-strong-coupling regimes of the QRM, respectively [15,16].
In the resonant-coupling case ω c /ω 0 = 1, the density matrix of the thermal state just depends on the two ratios ω 0 β and g/ω 0 . Therefore, we can clearly investigate the dependence of the logarithmic negativity N on the ratios ω 0 β and g/ω 0 . In this way, the dependence of the thermal entanglement in the QRM on the system parameters can be analyzed clearly. In the nonresonantcoupling case, we can investigate the dependence of the thermal entanglement on the two ratios ω c /ω 0 and g/ω 0 when the ratio ω 0 β takes different values.
In Fig. 5, we plot the logarithmic negativity N of the QRM in the thermal state as a function of ω 0 β and g/ω 0 in the resonant case ω 0 = ω c . Here the thermal state ρ th (T ) just depends on the two ratios ω 0 β and g/ω 0 . Therefore, Fig. 5(a) shows the contour map of the thermal entanglement in the overall view. In a macroscopic view, thermal entanglement shows a wedge-ridge pattern. The peak value of the entanglement is approximately located around the coupling g/ω 0 ≈ 1. In addition, we can see that the thermal entanglement disappears in the hightemperature limit ω 0 β ≪ 1. In the low-temperature limit ω 0 β ≫ 1, the thermal entanglement approaches a stable value, which is the ground state entanglement of the QRM. This point can be seen more clearly in panel (b). Here, we see that for a given g/ω 0 , the thermal entanglement increases with the increase of ω 0 β, and then the thermal entanglement tends to be a stable value in strong, ultrastrong, and deep-strong coupling regimes. We note that this feature can be understood by analyzing the dependence of the eigenstate population p |εm on the ω 0 β. In addition, we can see that the stable value of the logarithmic negativity N is in a monotonous order with respect to the coupling strength g/ω 0 . In particular, for a larger coupling strength g/ω 0 , the value of ω 0 β corresponding to the turning point is larger. To see the dependence of the logarithmic negativity N on the coupling strength, in Fig. 5(c) we plot the N as a function of g/ω 0 at different values of ω 0 β. Here we can see that, for ω 0 β = 0.01, the system is in the hightemperature limit, and hence the thermal engtanglement disappears. For ω 0 β = 600, the N increases with the increase of g/ω 0 . When ω 0 β = 40 and 80, however, we find that the N experiences a nonmonochromatic change, it increases firstly and then decreases. These features can be explained by analyzing the dependence of the populations of these eigenstates on the value of ω 0 β.
In Fig. 6, we show the populations p |εm as functions of ω 0 β at different values of g/ω 0 . Here, we can see that, for a given value of g/ω 0 , the ground state probability p |ε1 increases with the increase of the inverse temperature ω 0 β. Accordingly, the populations p |εm for m = 2, 3, 4 decrease with the increase of ω 0 β. For a larger value of g/ω 0 , a larger value of ω 0 β is needed to make sure the system reaching its ground state. For example, ω 0 β ≈ 600 is needed to make sure the system reaching the ground state for g/ω 0 = 1.5 as shown in Fig. 6(d). For the cases of ω 0 β = 40 and 80 in Fig. 5(c), when the coupling strength g/ω 0 is larger than a value around 1, the system will be in a mixed state involving the ground state and the first excited state. Here, the ground state |ε 1 and the first excited state |ε 2 are near degenerate, as shown in Fig. 2(b), and hence the entanglement is not a monochromatic function of the coupling strength.
For the case of ω 0 β = 600, when the coupling strength changes from g/ω 0 = 0.01 to 1.5, the system is always in the ground state, then the engtanglement in the ground state increases with the increase of the coupling strength, this feature matches the result when the entanglement is measured by the von Neumann entropy [54][55][56].
In Fig. 7, we plot the logarithmic negativity N of the QRM in the thermal state as a function of the ratios ω c /ω 0 and g/ω 0 , when the bath temperature takes different values. From Figs. 7(a) to 7(d), we take ω 0 β = 600, 30, 2, and 0.01, corresponding to a gradual change of the thermalized temperature from the low-temperature limit to the high-temperature limit. From Fig. 7, we can see three features of the thermal entanglement. (i) Thermal entanglement exists in the low-temperature and finittemperature cases, and there is no thermal entanglement in the high-temperature limit. This result is understandable because the thermal noise is harmful to quantum effect, and hence the quantum entanglement disappears in the high-temperature limit. In the zero-temperature limit, the thermal equilibrium state is reduced to the ground state, and the ground-state entanglement in the QRM has been studied in some previous works [54][55][56][57]. In the high-temperature limit, the thermal state becomes a completly mixed state, i.e., the density operator in the eigenstate representation becomes an identical matrix. Therefore, the thermal entanglement disappears in the high-temperature limit. (ii) In the finit-temperature case, the thermal entanglement exhibits a wedge-ridge pattern. For a given value of ω c /ω 0 , the logarithmic negativity N increases with increase of the coupling strength g/ω 0 . For a given coupling g/ω 0 , the logarithmic negativity N first increases and then decreases when the ratio ω c /ω 0 changes from 1 to a much larger value. (iii) At a relatively low temperature, the peak value of the logarithmic negativity N appears approximately around the resonance-coupling case, i.e., ω c /ω 0 ≈ 1. In particular, the slope of the ridge line corresponding to the peak value of entanglement increases slightly with the increase of the coupling strength g/ω 0 . For a higher thermalized temperature, the ridge line moves slightly towards the increasing direction of the value ω c /ω 0 .

VI. DISCUSSIONS ON THE EXPERIMENTAL IMPLEMENTATION
In this section, we present some discussions on the experimental implementation of this scheme. Concretely, we analyze the experimental feasibility based on four aspects. (i) The realization of the QRM with practical systems. The physical mode considered in this paper is general, and it can be realized in realistic physical setups with which the QRM can be implemented. Currently, the QRM can be realized in many systems, including various cavity-QED systems [14,58], circuit-QED systems [59], and other hybrid systems [60,61] described by the interaction between a two-level system (artificial atoms, spins, charge states etc.) and a single-mode bosonic field (cavity field, mechanical resonator, LC resonator etc.). (ii) The parameter condition of the physical model. In this work, the ultrastrong-and deep-strong-coupling regimes are considered. Therefore, the candidate physical system should enter these two regimes. Currently, both the ultrastrong and deep-strong couplings have been realized in many physical platforms [15,16]. (iii) The control of the effective bath temperatures. To study the thermalization, the candidate systems should be chosen such that the effective bath temperature can be tuned to control the equilibrium and nonequilibrium environments. Experimentally, the control of the effective bath temperature can be realized in circuit-QED systems [62]. (iv) The measurement of the eigenstate populations and eigenstate entanglement. In realistical physical systems, the measurement of the two-level system and the bosonic mode can be realized. However, it remains unexplored for measuring the population and entanglement of the eigenstates for the QRM. Therefore, to experimentally implement the present scheme, the experimental measurement of the eigenstate population and entanglement need to be developed. Based on the above analyses, we can conclude that the present scheme should be experimentally accessible once the experimental measurement of the eigenstate population and entanglement can be realized. Thanks to the recent great advances in ultra-strong couplings, much effort is devoting to the study of the measurement of the ultrastrongly coupled systems in the eigenstate representation. This will provide a hopeful prospect for the experimental realization of this scheme.

VII. CONCLUSION
In conclusion, we have studied the quantum thermalization of the open QRM, which is connected with either two IHBs or a CHB. We have derived the global quantum master equations in the eigenstate representation of the QRM to govern the evolution of the QRM in both the IHB and CHB cases. Based on the steady-state populations of the eigenstates, we have studied the quantum thermalization of the QRM by checking whether all the effective temperatures between any two eigenstates are the same. In the IHB case, we have found that, when the two IHBs have the same temperature (different temperatures), the QRM can (cannot) be thermalized. When one of the two IHBs is decoupled from the QRM, then the QRM can be thermalized with the coupled bath. In the CHB case, the system can always be thermalized. Additionally, we have studied quantum entanglement between the TLS and the bosonic field when the QRM is in a thermal state. The logarithmic negativity for the thermal state of the QRM has been calculated in both the resonant-and nonresonant-coupling cases. The dependence of the thermal entanglement on the system parameters and the thermalized temperature has been found and discussed in detail. This work will give an insight to quantum information processing at finite temperature.