Quantum phase transitions in the spin-boson model without the counterrotating terms

TL;DR

This paper investigates quantum phase transitions in the spin-boson model without counterrotating terms, revealing the presence of both first- and second-order transitions and challenging previous understandings under the rotating-wave approximation.

Contribution

It provides a numerically exact analysis of the spin-boson model without counterrotating terms, discovering new types of quantum phase transitions and their dependence on bath exponents.

Findings

Second-order QPT observed in sub-Ohmic bath

First-order QPTs appear before critical points

Multiple first-order QPTs in Ohmic bath at strong coupling

Abstract

We study the spin-boson model without the counterrotating terms by a numerically exact method based on variational matrix product states. Surprisingly, the second-order quantum phase transition (QPT) is observed for the sub-Ohmic bath in the rotating-wave approximations. Moreover, first-order QPTs can also appear before the critical points. With the decrease of the bath exponents, these first-order QPTs disappear successively, while the second-order QPT remains robust. The second-order QPT is further confirmed by multi-coherent-states variational studies, while the first-order QPT is corroborated with the exact diagonalization in the truncated Hilbert space. Extension to the Ohmic bath is also performed, and many first-order QPTs appear successively in a wide coupling regime, in contrast to previous findings. The previous pictures for many physical phenomena for the spin-boson model in the rotating-wave approximation have to be modified at least at the strong coupling.

Equations

$$H_{SB}=\frac{\Delta }{2}{\sigma }_z+\sum _k{\omega }_k{a}_k^{†}{a}_k+\sum _k{g}_k\left(a_k^{†}+a_k\right){\sigma }_x,$$

$$J(\omega )=\pi \sum _k{g}_k^2\delta (\omega -{\omega }_k)=2\pi \alpha {\omega }_c^{1-s}{\omega }^s\Theta ({\omega }_c-\omega ),$$

$$H_{SB}^{RWA}=\frac{\Delta }{2}{\sigma }_z+\sum _k{\omega }_k{a}_k^{†}{a}_k+\sum _k{g}_k\left(a_k^{†}{\sigma }_{-}+a_k{\sigma }_{+}\right).$$

$$\widehat{H_{SB}}=\frac{\Delta }{2}{\sigma }_z+\frac{ϵ}{2}{\sigma }_x+\sum _k{\omega }_k{a}_k^{†}{a}_k+\frac{1+\lambda }{2}\sum _k{g}_k\left(a_k^{†}{\sigma }_{-}+a_k{\sigma }_{+}\right)+\frac{1-\lambda }{2}\sum _k{g}_k\left(a_k{\sigma }_{-}+a_k^{†}{\sigma }_{+}\right),$$

$$\hat{H}$$

$$H_{\text{chain}}$$

$$c_0$$

$${ϵ}_n$$

$$t_n$$

$${\xi }_s$$

$$A_n$$

$$C_n$$

$$N_n^2$$

$${\left(a:q\right)}_n=\left(1-a\right)\left(1-aq\right)...\left(1-a{q}^{n-1}\right)$$

$$|\psi ⟩=\sum _{\{N_n\}=1}^{d_n}M\left[N_1\right]...M\left[N_L\right]|N_1,...,N_L⟩,$$

$$|{\psi }_0⟩$$

$$|{\psi }_1⟩$$

$$|{\psi }_2⟩$$

$$|{\psi }_3⟩$$

$${|\psi ⟩}^{\le N}=\sum _{l=0}^N|{\psi }_l⟩.$$

$$H^T$$

$$|\psi^{T}\rangle=\left(\begin{array}[]{c}\sum_{n=1}^{N_{c}}A_{n}\exp\left[\sum_{k=1}^{L}f_{n,k}\left(a_{k}^{{\dagger}}-a_{k}\right)\right]|0\rangle\\ \sum_{n=1}^{N_{c}}B_{n}\exp\left[\sum_{k=1}^{L}h_{n,k}\left(a_{k}^{{\dagger}}-a_{k}\right)\right]|0\rangle\end{array}\right),$$

$$E=\frac{⟨{\psi }^T|H^T|{\psi }^T⟩}{⟨{\psi }^T|{\psi }^T⟩},$$

$$⟨{\psi }^T|H^T|{\psi }^T⟩$$

$$⟨{\psi }^T|{\psi }^T⟩$$

$$F_{m,n}$$

$$G_{m,n}$$

$${\Gamma }_{m,n}$$

$${\alpha }_{m,n}$$

$${\beta }_{m,n}$$

$${\gamma }_{m,n}$$

$$\frac{\partial E}{\partial A_n}=\frac{\partial E}{\partial B_n}=\frac{\partial E}{\partial f_{ij}}=\frac{\partial E}{\partial h_{ij}}=0$$

$$\sum _n\left(2A_n{F}_{i,n}\left({\alpha }_{i,n}-E\right)-{\Gamma }_{i,n}{B}_n{\gamma }_{i,n}\right)=0,$$

$$\sum _n\left(2B_n{G}_{i,n}\left({\beta }_{i,n}-E\right)-{\Gamma }_{n,i}{A}_n{\gamma }_{n,i}\right)=0,$$

$$\sum _n\{-{\Gamma }_{i,n}{B}_n\left(h_{n,j}{\gamma }_{i,n}+\lambda g_j\right)$$

$$+A_n{F}_{i,n}\left[2\left({\alpha }_{i,n}+{\omega }_j-E\right)f_{n,j}+g_j\right]\}=0,$$

$$\sum _n\{-{\Gamma }_{n,i}{A}_n\left(f_{n,j}{\gamma }_{n,i}-\lambda g_j\right)$$

$$+B_n{G}_{i,n}\left[2\left({\beta }_{i,n}+{\omega }_j-E\right)h_{n,j}-g_j\right]\}=0.$$

$$E=-P_{+}{\log }_2{P}_{+}-P_{-}{\log }_2{P}_{-}$$

Full text

I Introduction

The spin-boson model describes a qubit (two-level system) coupled to a quantum environment represented by a continuous bath of bosonic modes. It is a paradigmatic model in many fields, ranging from quantum optics Scully to condensed matter physics Leggett to open quantum systems Breuer ; weiss1 . The Hamiltonian is given by

[TABLE]

where $\sigma_{i}$ ($i=x,y,z$) is the Pauli matrices, $\Delta$ is the qubit frequency, $a_{k}$ ($a_{k}^{{\dagger}}$) is the bosonic annihilation (creation) operator which can annihilate (create) a boson with frequency $\omega_{k}$, and $g_{k}$ denotes the coupling strength between the qubit and the bosonic bath, which is usually characterized by the power-law spectral density $J(\omega)$,

[TABLE]

where $\alpha$ are a dimensionless coupling constant, $\omega_{c}$ is the cutoff frequency, and $\Theta(\omega_{c}-\omega)$ is the Heaviside step function. The bath exponent $s$ classifies the reservoir into super-Ohmic $\left(s>1\right)$, Ohmic $\left(s=1\right)$, and sub-Ohmic $\left(s<1\right)$ types. In many theoretical studies, due to the weak coupling strength in the real quantum optical and quantum dissipative systems, the counterrotating terms involving higher excited states, $a_{k}^{{\dagger}}\sigma_{+}$ and $a_{k}\sigma_{-}$, can be neglected, a condition which is the so-called rotating-wave approximation (RWA); thus the full Hamiltonian (1) can be reduced to the following RWA form

[TABLE]

It is generally believed that the RWA is a reasonably good approximation because the counterrotating terms violate energy conservation, leading to virtual processes, and thus are suppressed.

The total excitation number of the spin-boson model is $\hat{N}=\sum_{k}a_{k}^{{\dagger}}a_{k}+\sigma_{+}\sigma_{-}$, where $\sigma_{\pm}=\left(\sigma_{x}\pm i\sigma_{y}\right)/2$. With (without) the RWA, the system possesses a $U(1)$ ($Z_{2}$) symmetry, which is characterized by the action of the operator $\Re\left(\theta\right)=\exp\left(i\theta\hat{N}\right)$, where the arbitrary angle $\theta$ corresponds to $U(1)$ symmetry and the special value of \theta=\pi\  corresponds to $Z_{2}$ (parity) symmetry. The parity operator $\hat{\Pi}=\Re\left(\pi\right)$ has two eigenvalues  $\pm 1$.

With the advance of modern technology, various qubit and oscillator coupling systems can be engineered in many solid-state devices, such as superconducting circuits Niemczyk ; Yoshihara , cold atoms Dimer , and trapped ions Cirac . Recently, the spin-boson model has been realized by the ultrastrong coupling of a superconducting flux qubit to an open one-dimensional (1D) transmission line  Forn2 . The counterrotating terms can be strongly suppressed in some proposed schemes AnistropicDicke ; Keeling ; Fanheng . In some systems, the anisotropy appears quite naturally, because they are controlled by different input parameters  Grimsmo1 . On the theoretical side, the anisotropic models where the rotating and counterrotating terms are different have attracted considerable attentions. Rich quantum phase transitions (QPTs) in the anisotropic Dicke models including the Tavis-Cummings model TC without counterrotating terms have been reported AnistropicDicke in the thermodynamic limit, i.e., infinite qubits. It was found in JC_Plenio that the Jaynes-Cummings (JC) model  JC , the quantum Rabi model in the RWA, could also undergo the second-order QPT in the extreme model parameter limit, $\Delta/\omega\rightarrow\infty$ where $\Delta$ and $\omega$ are the frequencies of qubit and cavity. It was demonstrated that the ratio of frequencies $\Delta/\omega$ plays the same role as the qubits number, and the second-order QPT can be identified by finite frequency scaling analysis.

Since both the Dicke and Rabi models in the RWA undergo the second-order QPT, can the spin-boson model in the RWA display the second-order QPT? Are some extreme conditions for model parameters required to realize the second-order QPT? As we know from the literature, the spin-boson Hamiltonian with the RWA is usually treated in sub-space with fixed total excitations $N$. In this way the second-order QPT is absolutely excluded. The constrained condition for the given sub-space should be relaxed to detect more physical phenomena.

In the sub-Ohmic spin-boson model without the RWA, the second-order QPT from the delocalized phase, where spin has the equal probability in the two states, to localized phase, in which spin prefers to stay in one of the two states, has been studied extensively Bulla ; QMC ; ED ; Zhang ; guo2012critical ; Frenzel ; CRduan . Many advanced numerical approaches in quantum many-particle physics have been applied and extended to this model, such as the numerical renormalization group Bulla , quantum Monte Carlo simulations QMC , sparse polynomial space approach ED , exact diagonalization in terms of shift bosons Zhang , and various matrix product state approaches guo2012critical ; Frenzel . Some analytical approaches based on the polaronic unitary transformation, also known as the Silbey-Harris ansatz  Silbey , have been also developed for this model zheng ; Chin ; mD1 ; multi_CS2 ; yzhao ; Blunden ; ZH ; cao ; HeShu1 . The single-coherent-states ansatz Silbey was improved by simply adding other coherent states on an equal footing mD1 and by superpositions of two degenerate single coherent states ZH , which are generally termed as the multi-coherent-states (MCS) ansatz. Actually, the MCS in the single-mode model was proposed much earlier by Ren and Chen tcs .

In this paper, we will extend the variational matrix product state (VMPS) approach  guo2012critical to study the spin-boson model in the RWA for all values of the bath exponents. The MCS variational approach and exact diagonalization within truncated Hilbert space are also employed to provide independent checks in different regimes. The paper is organized as follows. In Sec. II, we briefly introduce the generalized spin-boson model. Some methodologies including the VMPS, the MCS variational approaches, and exact diagonalization in truncated Hilbert space are described. The rich phase transitions revealed by the VMPS method are presented in Sec. III, where the MCS variational approaches and the exact diagonalization are also applied to provide further evidence. The quantum criticality based on VMPS studies is also analyzed. Finally, conclusions are drawn in Sec. IV.

II Generalized Model Hamiltonian and methodologies

Based on Hamiltonian (1), the generalized spin-boson Hamiltonian can be written as ($\hbar=1$)

[TABLE]

where $\epsilon$ is the energy bias applied in a two-level system ($\epsilon=0$ except special statements) and \frac{1+\lambda}{2}\and \frac{1-\lambda}{2}\ are the weights of the rotating-wave and counterrotating terms, respectively. In this sense $\lambda$ is the anisotropy constant of this model. Obviously, $\lambda=1$ ($\lambda=0$) corresponds to the spin-boson model with (without) RWA.

For later use, Hamiltonian (4) can be rewritten as

[TABLE]

In the following, three methods are introduced to study this generalized model.

VMPS approach.-. As is well known, the VMPS approach works efficiently in one-dimensional chain models  MPS_intro1 ; MPS_intro2 . To apply VMPS in the spin-boson model, we therefore transform the model into a 1D chain model. We first perform the logarithmic discretization of the spectral density of the continuum bath Bulla with discretization parameter $\Lambda>1$; then, by using orthogonal polynomials $b_{n}^{{\dagger}}=\sum_{k}U_{nk}a_{k}^{{\dagger}}$ ($b_{n}=\sum_{k}U_{nk}a_{k}$) as described in Ref. Plenio_Chin , the spin-boson models can be mapped into the representation of a one-dimensional semi-infinite chain with nearest-neighbor interaction  Friend . Thus, Hamiltonian (5) can be written as:

[TABLE]

where $b_{n}^{{\dagger}}$($b_{n}$) is the creation (annihilation) operator for a new set of boson modes in a transformed representation with $\epsilon_{n}$ describing frequency on chain site $n$, $t_{n}$ describing the nearest-neighbor hopping parameter, and $c_{0}$ describing the effective coupling strength between the spin and the new effective bath. All of the parameters mentioned above, such as $t_{n},\epsilon_{n}$ and $c_{0}$, are determined by the logarithmic discretization parameter $\Lambda$, the cutoff frequency $\omega_{c}$, and a specific form of the spectral density, which are expressed below

[TABLE]

where

[TABLE]

with

[TABLE]

For details, one may refer to Ref. Plenio_Chin .

We now briefly introduce the VMPS approach VMPS1 ; VMPS2 ; MPS_intro1 . For the transformed spin-boson model of a 1D chain with $L$ sites, the ground-state wave function of Hamiltonian (6) can be depicted as

[TABLE]

where $N_{n}$ is the physical dimension of each site $n$ with truncation $d_{n}$, and $D_{n}$ is the bond dimension for matrix $M$ with the open boundary condition, bounding the maximal entanglement in each subspace. $M$ on each site is optimized through sweeping the 1D chain iteratively, where accuracy of numerical results is determined by values of $d_{n}$ and $D_{n}$. In order to deal with the spin-boson model near the quantum critical region effectively, we apply an optimized boson basis through an additional isometric map with $d_{opt}\ll d_{n}$ like in Refs. guo2012critical ; Friend . In this way, we can effectively reduce the local boson basis and improve the maximum boson number in the quantum critical region.

We have recovered all the results in Ref. guo2012critical and confirmed the classical mean-field behavior for $s<1/2$ of the sub-Ohmic spin-boson model. In the present paper, we will extend it to study the RWA spin-boson model, mainly focusing on phase transitions in the sub-Ohmic baths, namely $\lambda=1$, $0<s<1$. For the data presented below, we typically choose the model parameters as $\Delta=0.1$, $\omega_{c}=1$, and $\epsilon=0$, if there are no special statements. Some parameters of orthogonal polynomials transformation and VMPS are the same as those in Ref. guo2012critical , e.g.  the logarithmic discretization parameter $\Lambda=2$, the length of the semi-infinite chain $L=50$, and optimized truncation numbers $d_{opt}=12$. In addition, we adjust the bond dimension to achieve better convergence of the results. In this paper, we choose $D_{max}=20,40$ for $s=0.3,0.7$, respectively, which is sufficient to obtain the converged results, as demonstrated in Appendix A in detail.

Exact diagonalization in truncated Hilbert space.-. It is well known that the spin-boson model with the RWA possesses $U(1)$ symmetry, and the total excitation number $\hat{N}$ is conserved because $\left[\hat{N},H\right]=0$. It has been reported in Ref. Tong that the total excitation number $N$ in the ground state of the RWA spin-boson model jumps from [math] to $1$ at a critical coupling strength. This instability can be also called the first-order phase transition, because the first derivative of the ground-state energy is discontinuous. Due to the possible sequence of instabilities with the coupling strength, we can truncate the Hilbert space up to a finite $N$ excitation number. To this end, we first separate the Hilbert space into several subspaces with different excitation numbers $l=0,1,2...N$. The wave function in $l$-subspace $\left|\psi_{l}\right\rangle$ can be written explicitly with $l$ excitations. E. g. the wavefunctions for $l=0,1,2$ and $3$-subspace are listed in the following

[TABLE]

Then the wavefunction in the truncated Hilbert space up to $N$ excitations can be expressed as

[TABLE]

For $N=0$, $\left|\psi\right\rangle^{\leqslant 0}=\left|0\right\rangle\left|\downarrow\right\rangle$, the ground-state energy is $E_{0}=-\frac{\Delta}{2}$. However, it is very difficult to obtain the analytical solution for $\left|\psi\right\rangle^{\leqslant N}$ up to $N>0$ excitations, and numerically exact diagonalizations are then required to obtain the converged lowest energy. This approach is called $N$ED below.

MCS ansatz.-. We also apply the MCS ansatz tcs ; mD1 ; multi_CS2 to the spin-boson model in the RWA. To facilitate the variational study and visualize the symmetry breaking explicitly, we rotate the Hamiltonian (5) around the y axis by an angle $\pi/2$ with $\epsilon=0$, which gives

[TABLE]

The trial state $|\psi^{T}\rangle$ is written in the basis of the spin-up state $|\uparrow\rangle$ and spin-down state $|\downarrow\rangle$

[TABLE]

where $A_{n}$ ($B_{n}$) are related to the occupation probabilities of the spin-up (spin-down) state in the $n$th coherent state; $N_{c}$ and $L$ are numbers of coherent states and total bosonic modes, respectively; and $f_{n,k}$ ($h_{n,k}$) represents bosonic displacement of the $n$th coherent state and $k$ th bosonic mode. The symmetric MCS ansatz ($A_{n}=B_{n}$ and $f_{n,k}=-g_{n,k}$) can only be applied to the delocalized phase, so one can easily detect the symmetry breaking.

The energy expectation value can be calculated as follows

[TABLE]

where

[TABLE]

with

[TABLE]

Minimizing the energy expectation value with respect to variational parameters gives the following self-consistent equations

[TABLE]

which leads to

[TABLE]

In practice, these parameters can be obtained by solving the coupled equations self-consistently, which in turn give the ground-state energy and wave function. The ground state with zero excitation $\left|0\right\rangle\left|\downarrow\right\rangle$ to Hamiltonian (5) can be contained in the MCS wave function to Hamiltonian (9) by setting the constrained coefficients: $A_{n}=B_{n}$ and$\ f_{n,k}=h_{n,k}=0$ for all $n$. However the state with nonzero total excitations cannot be included in the  MCS ansatz due to the finite number of coherent states. It has  been demonstrated that this wave function can describe the localized phase of the spin-boson model Blunden . The number of the coherent states in the practical calculations in this paper is $N_{c}=6$, which is sufficient to judge the existence of the second-order QPT.

For the three approaches described above, discretization of the energy spectrum of the continuum bath should be performed at the very beginning in the practical calculations. The same logarithmic discretization is taken for different approaches if comparison is made in the data presented below.

III Results and discussions

We first describe our main results by VMPS method to the sub-Ohmic $\left(0<s<1\right)$ spin-boson model with the RWA described by Eq. (5) for $\lambda=1$. We observe a second-order QPT in the spin-boson model under RWA, which was unnoticed in literature, to the best of our knowledge. This surprising observation is further confirmed by the MCS variational approach. The symmetry breaking was unambiguously found above the critical point in this wave function based approach. The $N$ED results are also given, which should be exact at weak coupling and can be regarded as a benchmark in this regime. Besides the second-order QPT, we also find a few first-order QPTs before the critical point of the second-order QPT for large bath exponent $s$. As $s$ decreases, the first-order QPTs disappear successively, but the second-order phase transition remains robust no matter how small $s$ is. We will discuss those phenomena based on various numerical calculations in the following subsections.

III.1 Magnetization $\left|\left\langle\sigma_{x}\right\rangle\right|$

In terms of Hamiltonian (5), magnetization $\left|\left\langle\sigma_{x}\right\rangle\right|$ can be regarded as the order parameter in this model. In the second-order phase transition, due to the symmetry breaking, the order parameter changes from zero to nonzero at the critical points. In Fig. 1 we present the VMPS results for $\left|\left\langle\sigma_{x}\right\rangle\right|$ as a function of the coupling strength $\alpha$ for the spin-boson model in the RWA (i.e.$\ \lambda=1$) for two typical values of $s=0.3$ and $0.7$. Surprisingly, magnetization changes abruptly from zero to nonzero for both cases. The critical points of the second-order QPTs are $\alpha_{c}=0.016,0.685$ for s=0.3,0.7\ respectively.

However, the $N$ED within excitation numbers up to $N$ shows that the order parameter $\left|\left\langle\sigma_{x}\right\rangle\right|$ remains zero for all coupling strengths. Note that the VMPS has provided convincing results for the full spin-boson model guo2012critical . $N$ED cannot yield consistent results with the VMPS ones, because the total excitation number may not conserve with $\alpha$ due to the unexpected symmetry breaking.

In the finite-size Dicke model with the RWA Buzek ; Zhangyy , the system usually undergoes the first-order phase transition, i.e., sequence of instabilities, among the phases within different conserved excitation numbers as the coupling strength increases. In the limit $\Delta/\omega\rightarrow\infty$, the second-order quantum phase transitions have been observed in the quantum Rabi model with the RWA JC_Plenio . It has been also reported that the Dicke model under the RWA displays a second-order QPT in the thermodynamic limit AnistropicDicke . Contrary to the quantum Rabi (Dicke) model in the RWA, the spin-boson model under the RWA could undergo the second-order QPT even for a finite value of $\Delta/\omega$ ( one qubit).

III.2 Magnetization $\left|\left\langle\sigma_{z}\right\rangle\right|$ and the ground-state energy

In the original spin-boson Hamiltonian, $\sigma_{z}$ describes a tunneling two-level system  zheng ; Chin . The magnetization along the $z-$direction $\left|\left\langle\sigma_{z}\right\rangle\right|$, simply the order parameter along z direction, is the renormalized factor of the tunneling amplitude $\Delta$. In this subsection, we examine the magnetization along the z direction as well as the ground-state energy by both VMPS and $3$ED, which are exhibited in Fig. 2.

Obviously, two observables begin to deviate only after the critical point $\alpha_{c}$. Especially we find that the ground-state energies by the VMPS become lower than those by $3$ED after the critical points, indicating again the invalidity of ED method at strong coupling. Practically, we cannot perform ED with a very large total excitation number due to the huge Hilbert space. In this paper, our exact diagonalization is only performed up to $\ N=3$. Of course, if one can really perform true exact diagonalization without the limitation of total excitation numbers, the true ground state could be also correctly described in the ED method.

It is interesting to note from Figs. 1 and 2 that $\left|{\left\langle\sigma_{x}\right\rangle}\right|$ and $\left|{\left\langle\sigma_{z}\right\rangle}\right|$ exhibit different behaviors with the increase of the coupling strength. We will discuss this issue based on parity symmetry breaking and total excitation numbers in the next subsection.

III.3 Parity symmetry breaking

As stated before, the spin-boson model with (without) the RWA possesses a $U(1)$ ($Z_{2}$) symmetry. $U(1)$ is a higher symmetry than $Z_{2}$, so in the RWA spin-boson model, the system also has $Z_{2}$ symmetry, i.e., parity symmetry, like in the full spin-boson model. In this section, we study the behavior of the expectation value of the parity $\hat{\Pi}=\exp\left(i\pi\hat{N}\right)$.

The finite order parameter above the critical points obtained by VMPS in the previous subsection displays spontaneous symmetry breaking, while parity is generally just the criterion to determine whether symmetry is broken. The upper panel in Fig. 3 gives the expectation value of parity \left\langle\hat{\Pi}\right\rangle\obtained by both methods for $s=0.3,0.7$. Before the critical points $\alpha_{c}$, both methods yield the same results for the parity. It is interesting to find that in this regime, the value of parity for all values of $s<1$ is either $1$ or $-1$, which corresponds to even or odd parity respectively. It follows that the symmetry is not broken in this regime. Nevertheless, the average parity $\left\langle\hat{\Pi}\right\rangle$ becomes zero due to quantum fluctuations above the critical points. It is not the eigenstate of the parity operator, indicating spontaneous parity symmetry breaking. Note also that the parity jumps for a few times before the critical coupling for $s=0.7$, and remains unchanged for $s=0.3$.

We also present the average value of the total excitation number $\left\langle\hat{N}\right\rangle$ in the ground-state as a function of the coupling strength$\ \alpha$ in the lower panel of Fig. 3. Below the critical points, both approaches give the same $\left\langle\hat{N}\right\rangle$, which jumps between different plateaus with different integers, just like in the JC model JC_Plenio . For $s=0.3$, $\left\langle\hat{N}\right\rangle=0$ remains until the second-order critical point $\alpha_{c}$, while, for $s=0.7$, $\left\langle\hat{N}\right\rangle$ increases from [math] to $2$ one by one before $\alpha_{c}$. It follows that the lowest energies in different coupling regime belong to the energy levels in different $\left\langle\hat{N}\right\rangle$ subspaces, leading to level crossing at some coupling strengths. Thus the first derivative of the ground-state energy with respect to the coupling strength must be discontinuous at these coupling strengths, so instability of $\left\langle\hat{N}\right\rangle$ before $\alpha_{c}$ just corresponds to the first-order phase transitions. In addition, the jump of $\left\langle\hat{N}\right\rangle$ can also account for the back and forth of the parity between $1$ and $-1$.

One can also note that the total excitation number$\ \left\langle\hat{N}\right\rangle$ by VMPS increases abruptly at the critical coupling $\alpha_{c}$. The $N$ED approach can only describe the phase with excitation number less than or equal to $N$. The total excitation number is not limited in VMPS, so it in principle can describe all phases. The total excitation numbers are not conserved in the ground state above $\alpha_{c}$, because of the symmetry breaking. At the critical points, $\left\langle\hat{N}\right\rangle$ does not jump to a plateau with the finite larger integer, different from the first-order phase transitions.

The instability of total excitation number $\left\langle\hat{N}\right\rangle$ at the phase transitions can account for the rich behavior of $\left|\langle\sigma_{z}\rangle\right|$ here. For any value of $s$, in the weak-coupling regime, the corresponding ground state is the spin-down state with photonic vacuum for $\left\langle\hat{N}\right\rangle=0$, i.e., $|\downarrow\rangle|0\rangle$, thus $\langle\sigma_{z}\rangle=-1$, as just demonstrated in the upper panel of Fig. 2 for $s=0.3,0.7$. Once $\left\langle\hat{N}\right\rangle\neq 0$ at the first phase transition no matter whether it is of the first or the second order, the spin state in the ground state consists of both spin-up and spin-down states, which is drastically different from the spin state only including the spin-down state at $\left\langle\hat{N}\right\rangle=0$, leading to the jump of $\left|\langle\sigma_{z}\rangle\right|$ and $\left|\langle\sigma_{z}\rangle\right|\neq 1$, as exhibited in Fig. 2. At the later phase transitions between two finite $\left\langle\hat{N}\right\rangle\neq 0$, the ground states do not change drastically because both have spin-up and spin-down states, therefore $\left|\langle\sigma_{z}\rangle\right|$ will exhibit kinks, as shown in Fig. 2 for $s=0.7$. The drop of $\left|\langle\sigma_{z}\rangle\right|$ to zero is an artefact due to the plot using magnetization, which is the absolute  value of $\langle\sigma_{z}\rangle$ in our paper. In the weak-coupling limit, $\langle\sigma_{z}\rangle$ is actually $-1$ as stated above. In the strong-coupling limit, we always  find $\langle\sigma_{z}\rangle>0$ in the ground-state, so,  in between $\langle\sigma_{z}\rangle$ must cross zero, resulting in the drop of $\left|\langle\sigma_{z}\rangle\right|$ to zero at some coupling that is not at any phase transition point. $\left|\langle\sigma_{x}\rangle\right|$ is the order parameter of the second-order QPT, so it only becomes non-zero for couplings larger than the critical one.

In short, we find rich phase transitions in the sub-Ohmic spin-boson model with the RWA. First,  the second-order QPT occurs for any finite model parameters, similar to its counterpart without the RWA. In both the Rabi model and the Dicke model with the RWA, the second-order QPTs cannot happen for finite ratio $\Delta/\omega$  or finite qubit number. Second, for larger bath exponents, e.g., $s=$ $0.7$, both the first- and second-order QPTs occur subsequently with the coupling strength. Several first-order phase transitions before the critical points indicate a sequence of instabilities. The first-order phase transition is absent for small bath exponent $s$, such as for $s=0.3$.

III.4 Evidence for the second-order QPT by MCS variational studies

Since this model at $s=0.3$ exhibits only one second-order QPT from the zero excitation ground state, we can employ the MCS variational approach to provide additional evidence, because the zero excitation can be realized and the localized phase can be also described in the trial wave function (10). In Fig. 4, we list results by MCS, VMPS, and $N=3$ED approaches for $s=0.3$. The logarithmic discretization parameter $\Lambda=2$ and $L=20$ bosonic modes are taken for all three approaches here. Note that the number of bosonic modes here is smaller than those in other figures due to the computational difficulties in the MCS approach, but it does not influence the essential results at all. The MCS approach is used here to account for the existence of the second-order QPT qualitatively, not for the precise location of the critical points.

Before the critical points, the results by all three methods are the same. Nonzero order parameter $|\left\langle\sigma_{x}\right\rangle|$ by MCS approach appears above the critical points, providing further convincing evidence of the spontaneous symmetry breaking. The MCS ground-state energies are lower than those by $N$ED method after the critical point, indicating that the state with broken symmetry is more stable. The total excitation number by MCS method increases suddenly above the critical point, because of no limitation of total excitation number in the coherent state. All these findings in the MCS variational study provide strong evidence of the second-order QPT in the spin-boson model with RWA. As found recently by Blunden-Codd et al. Blunden that a very accurate wave function can be only obtained by at least $100$ coherent states, by $N_{c}=6$ coherent states the MCS results for the order parameter and energy still slightly deviate from those by VMPS above the critical points.

III.5 The critical exponent for magnetization

Now we will study the nature of the second-order QPT in the sub-Ohmic spin-boson model under RWA. The most important question is whether the RWA changes the universality of the second-order QPT. The field relevant order-parameter critical exponent $\delta$ can be determined through the displayed power-law behavior $\left\langle\sigma_{x}\right\rangle\propto\epsilon^{1/\delta}$ at the critical coupling strength $\alpha=\alpha_{c}$. Previously, various critical exponents in the full spin-boson model have been calculated with different numerical approaches QMC ; ED ; Zhang ; Chin ; guo2012critical . It is generally accepted that the exponent $\delta$ takes the mean-field value $1/\delta=1/3$ for $s<1/2$, and the non-classical one $1/\delta=\left(1-s\right)/\left(1+s\right)$ by the exact hyperscaling for $s>1/2$ .

We present the magnetization by the VMPS method as a function of bias $\epsilon$ in a log-log plot for both $\lambda=1$ (RWA) and $\lambda=0$ (non-RWA) in Fig. 5. A very nice power-law behavior over three decades is demonstrated in all cases. The results for the full model are nearly the same as those in Ref. guo2012critical . It is interesting to find that $1/\delta$ is around $1/3$ for $s<1/2$ for both RWA and non-RWA cases. Very surprisingly, even for $s=0.7$ in the RWA model, $1/\delta$ is still close to that in the full model within the statistical errors, also indicating hyperscaling. Our results suggest that counterrotating terms would have no effect on critical exponent $\delta$ in the spin-boson model.

III.6 Extensions to the Ohmic bath.

Now we turn to the Ohmic spin-boson model under the RWA. By both VMPS and $3$ED methods, the above observables are also calculated. In Fig. 6, we collect magnetization $\left\langle\sigma_{z}\right\rangle$, parity $\left\langle\hat{\Pi}\right\rangle$, and total excitations $\left\langle N\right\rangle$ as a function of the coupling strength $\alpha$ in the range $(0,2)$. During this wide regime, one can find that the model undergoes a few first-order QPTs with the increment $1$ of the total excitation number. The two approaches almost give the same results. It follows that in the Ohmic bath, the $3$ED approach is suited to the wide coupling regime where no second-order QPT occurs. This result is obviously different from that in Ref. Tong where only $N=1$, i.e., single excitation, is considered. We have also studied the super-Ohmic bath $s=3/2$, and find that there is only one first-order QPT in the regime $\alpha=[0,2]$ (not shown here).

IV Conclusion

In this paper, we study the spin-boson model in the RWA by the VMPS method, MCS variational ansatz, and exact diagonalizations within the truncated Hilbert space. Surprisingly, we find the second-order QPT in the RWA model for any bath exponent $s<1$ for the first time, to the best of our knowledge. A rich picture for the quantum phase transitions is observed. Besides the second-order phase transition, the first-order phase transition also appears in the same model, which could only vanish for small bath exponents. The coexistence of both first- and second-order QPTs in the same model has never been observed in other spin-boson-like models, such as the quantum Rabi and Dicke models in the RWA.

Within the statistical error, for all values of the bath exponents $s<1$, the critical exponent $\delta$ is found to be nearly the same as those in the full model. It is then suggested that the counterrotating terms would have almost no effect on critical exponent $\delta$ in the critical regime. The analytical argument about the quantum-to-classical mapping  map ; berry ; Kirchner in the spin-boson model in the RWA would be helpful to account for this robust nature uniquely from the rotating-wave terms.

ACKNOWLEDGEMENTS We acknowledge useful discussions with Stefan Kirchner. This work is supported by the National Science Foundation of China (Nos. 11834005, 11674285), the National Key Research and Development Program of China (No. 2017YFA0303002),

∗ Email:qhchen@zju.edu.cn

Appendix A Convergence of results by the variational matrix product state

approach

We provide evidence for the full convergence of our VMPS results here.

The most important point of our paper is whether the second-order QPT occurs in the spin-boson model in the RWA. To demonstrate this point convincingly, we have checked the magnetization as a function of the coupling strength with increasing bond dimension (D). The results are shown in Fig. A1. One can see that the magnetization is converged even for the smallest one $D=12$, indicating that the second-order QPT in our model is reliable.

The field relevant order-parameter critical exponent $\delta$ should be determined precisely to study the universality class of the second-order QPT in this model. In doing so, we should check its convergence for the VMPS parameters carefully. Figures A2 and  A3 exhibit the nice power-law curves with the same fitting exponent for the selected bond dimension ($D$), optimized physical dimension ($d$), and discrete logarithm ($\Lambda$) for the same values of the bath exponent $s$. It is clearly demonstrated that we have obtained a convincing critical exponent of the second-order QPT in this paper.

Finally we turn to the convergence of the entanglement as a function of the inverse of the bond dimension 1/D. In the spin-boson model, the entanglement entropy can be rewritten as Chin

[TABLE]

where $P_{\pm}=\left(1\pm\sqrt{\left\langle\sigma_{x}\right\rangle^{2}+\left\langle\sigma_{y}\right\rangle^{2}+\left\langle\sigma_{z}\right\rangle^{2}}\right)/2$, which describes the correlation between the spin and the bosonic bath.

In Fig. A4, we plot the entanglement entropy as a function of 1/D for s=0.3 (upper) and s=0.7 (lower) at fixed coupling strength. It is clearly shown that the entropy saturates even before $D=20$ for both cases, providing further strong evidence for the convergence of our results for $D\geq 20$, i.e., the value selected in our calculation. Thus we believe that our algorithm has converged to the global minimum, the ground state is the correct one, and the results given here should be reliable.