Quantum Zeno effect in the multimode quantum Rabi model

TL;DR

This paper investigates the quantum Zeno and anti-Zeno effects in the multimode quantum Rabi model, deriving analytical decay rates and exploring how cavity mode squeezing influences these quantum phenomena.

Contribution

It provides an analytical expression for decay rates in the MQRM and explores the effects of squeezing, temperature, and high-frequency modes on the QZE and QAZE.

Findings

Crossover from QZE to QAZE due to energy backflow.

Decay rate depends on squeezing angle, strength, and temperature.

Numerically exact methods confirm the analytical results.

Abstract

We study the quantum Zeno effect (QZE) and quantum anti-Zeno effect (QAZE) of the multimode quantum Rabi model(MQRM). We derive an analytic expression for the decay rate of the survival probability where cavity modes are initially prepared as thermal equilibrium states. A crossover from QZE to QAZE is observed due to the energy backflow induced by high frequency cavity modes. In addition, we apply a numerically exact method based on the thermofield dynamics(TFD) theory and the matrix product states(MPS) to study the effect of squeezing of the cavity modes on the QZE of the MQRM. The influence of the squeezing angle, squeezing strength and temperature on the decay rate of the survival probability are discussed.

Equations

$$H=H_0+H_B+H_{\text{int}}$$

$$H_0=\frac{\Delta }{2}{\sigma }_z,H_B=\sum _{m=0}^{M-1}{\omega }_m{a}_m^{†}{a}_m$$

$$H_{\text{int}}=\sum _{m=0}^{M-1}{g}_m{\sigma }_x(a_m^{†}+a_m)$$

$$\rho (0)=|\uparrow ⟩⟨\uparrow |\otimes {\rho }_{\text{B}},{\rho }_{\text{B}}=S(\xi )\left(Z(\beta {)}^{-1}{e}^{-\beta H_B}\right)S^{†}(\xi )$$

$$S(\xi )=\prod _m{S}_m({\xi }_m)=\prod _m\exp \left(-\frac{1}{2}{\xi }_m^{\ast }{a}_m^2+\frac{1}{2}{\xi }_m(a_m^{†}{)}^2\right)$$

$$H^{\prime }=S^{†}(\xi )HS(\xi )=H_0+H_B^{\prime }+H_{\text{int}}^{\prime }$$

$$H_B^{\prime }=\sum _{m=0}^{M-1}{\omega }_m\left(A{a}_m^{†}{a}_m+B{a}_m^{†2}+B^{\ast }{a}_m^2\right)$$

$$H_{\text{int}}^{\prime }=\sum _{m=0}^{M-1}{g}_m\left(K{a}_m^{†}+K^{\ast }{a}_m\right){\sigma }_x$$

$$\hat{H}=H_0+{\hat{H}}_B+{\hat{H}}_{\text{int}}$$

$$\displaystyle\hat{H}_{B}=\sum_{m=0}^{M-1}\omega_{m}\big{[}A\left(a^{\dagger}_{m}a_{m}-b_{m}^{\dagger}b_{m}\right)+\left(Ba_{m}^{{\dagger}2}-Bb_{m}^{2}+\text{H.c}\right)\big{]}$$

$${\hat{H}}_{\text{int}}=\sum _{m=0}^{M-1}{g}_m(K{a}_m^{†}\cosh {\theta }_m+Kb\sinh {\theta }_m){\sigma }_x+\text{H.c}$$

$$|\hat{\psi }(0)⟩=|\uparrow ⟩\otimes \prod _{m=0}^{M-1}|0{⟩}_{a_m}|0{⟩}_{b_m}$$

$$⟨O(t)⟩=⟨\hat{\psi }(t)|O|\hat{\psi }(t)⟩$$

$$P_{\text{sur}}(n\tau )=\text{Tr}{\left[P{e}^{-i\hat{H}\tau }\rho (0)e^{i\hat{H}\tau }P\right]}^n=P_{\text{sur}}^n(\tau )$$

$$P_{\text{sur}}(n\tau )=\exp (-\gamma (\tau )t)$$

$$\gamma (\tau )=-\frac{1}{\tau }\ln \left[P_{\text{sur}}(\tau )\right]$$

$$P_{\text{sur}}(\tau )$$

$$=\text{Tr}\left[\sum _{n=0}^{\infty }\frac{(-i\tau {)}^2}{n!}[\hat{H},\rho (0){]}_n|e⟩⟨e|\right]$$

$$P_e(\tau )\approx 1-{\left(\frac{\tau }{{\tau }_z}\right)}^2$$

$$\displaystyle P(t=N\tau)=\bigg{[}1-\left(\frac{\tau}{\tau_{z}}\right)^{2}\bigg{]}^{N}\approx e^{-\gamma(\tau)t}$$

$${\gamma }_{\text{2nd}}(\tau )=\tau \sum _m{g}_m^2\left({\cosh }^2({\theta }_m)+{\sinh }^2({\theta }_m)\right)$$

$${\gamma }_{\text{th}}(\tau )=$$

$$+$$

$$P_{\text{sur}}(\tau )=(⟨{\sigma }_z(\tau )⟩+1)/2=E_{\text{TLS}}(\tau )/\Delta +1/2$$

$$\frac{d\gamma (\tau )}{d\tau }<0\iff -\frac{1}{\tau }\left(\frac{d\ln (P_{\text{sur}})}{d\tau }-\frac{\ln (P_{\text{sur}})}{\tau }\right)<0$$

$$E_{\text{TLS}}(\tau )/\Delta =1/2-a_1(g\tau )-a_2(g\tau {)}^2+O(g\tau {)}^2$$

$$\displaystyle\frac{d\gamma(\tau)}{d\tau}\approx\frac{d\gamma}{d\tau}\bigg{|}_{\tau=0}+\frac{4}{3}\tau a_{1}g\left(\frac{d\gamma}{d\tau}\bigg{|}_{\tau=0}+\frac{1}{2}a_{2}g^{2}\right)+O(\tau)<0$$

$${\partial }_t\rho (t)=-i[H,\rho (t)]$$

$$\rho (0)=\frac{1}{Z(\beta )}\exp \left(-\beta H\right)$$

$$i\frac{\partial }{\partial t}|\psi (t)⟩=\hat{H}|\psi (t)⟩$$

$$|I⟩=\sum _k|k⟩|\tilde{k}⟩$$

$$|\psi (t)⟩=\rho (t{)}^{1/2}|I⟩$$

$$|\psi (0)⟩=\rho (0{)}^{1/2}|I⟩=Z(\beta {)}^{-1/2}\exp \left(-\frac{1}{2}\beta H\right)|I⟩$$

$$⟨A(t)⟩=⟨\psi (t)|A|\psi (t)⟩=\text{Tr}\left\{\rho (t)A\right\}$$

Full text

I Introduction

Frequent measurements on a quantum system may cause its dynamical evolution slow down or accelerated. This phenomenon, known as the quantum Zeno effect(QZE) or quantum anti-Zeno effect (QAZE)misra1977zeno ; facchi2008quantum ; francica2010quantum ; schafer2014experimental , has been observed in experiments of trapped ionitano1990quantum , cavity QEDbernu2008freezing and nuclear spin ensemblesalvarez2010zeno . QZE has also been considered as a powerful strategy to implement the quantum control, including quantum communicationpouyandeh2014measurement ; bayat2015measurement , quantum information protectionbarenco1997stabilization , decoherence suppressionbeige2000quantum , purification and coolingerez2008thermodynamical .

Recently, the QZE has been observed in the circuit-QED system where a superconducting flux qubit is coupled to a transmission linekakuyanagi2015observation . Compared to the traditional optical experimental domain, the time scales involved in the circuit-QED system are much larger and can be resolved with ordinary electronics. Moreover, the ultrastrong-coupling strength between the the qubit and the resonator has already been realized in circuit-QED systemsforn2017ultrastrong , making it possible to observe QZE and QAZE in the strong coupling regime. As one of consequences in this regime, additional modes of the electromagnetic resonator become increasingly relevantbosman2017multi . To capture such multiply modes effect, the widely used quantum Rabi modelrabi1936process , describing an qubit interacting with a single electromagnetic mode, has to be generalized to its multimode version, known as multimode quantum Rabi model(MQRM)sundaresan2015beyond ; gely2017convergence . Recent studies have shown that such higher-lying electromagnetic modes of the resonator has a profound impact on various quantum optical phenomenons in the strong coupling regimesanchez2018resolution ; de2014light ; garcia2015light ; sundaresan2015beyond ; gely2017convergence . Thus, it is an interesting question that what role do multiple modes play in the quantum Zeno dynamics of the MQRM. In addition, previous studies have shown that the squeezing of the resonator mode has a significant influence in the quantum Zeno dynamics of a qubit in the quantum Rabi modellizuain2010zeno . It is natural to raise a question that how does the squeezing of multiple resonator modes affect the QZE and QAZE of the MQRM.

To address the above problems, we exploit a numerical exact method based on the matrix product states and time-dependent variational principlehaegeman2016unifying to study the QZE and QAZE of the MQRM. By restricting the evolution in the single excitation subspace under the framework of TFD, we derive an analytical expression for the decay rate of survival probability when multiply modes are initially prepared as a thermal equilibrium states. We observe a crossover between the QZE to the QAZE under repeated projective measurements. By numerically calculating the energy transport between the TLS and multiply modes, we show that this crossover is attributed to the energy back flow from the high frequency modes to the qubit. Moreover, we generalize the initial state of multiple modes to a squeezed thermal state and study effects of squeezing phase angle and amplitude to the QZE of the MQRM. We find that the decay of survival probability is accelerated by non-vanishing squeezing strengthes. Squeezing angles also significantly affect the decay rate. Particularly, high frequency modes in the MQRM cause a positive shift on the critical squeezing angles where the decay rate reaches its extremum.

This paper is organized as follows. In Sec II, we briefly introduce the MQRM and the numerical exact method we used to obtain the evolution of the MQRM at finite temperature. In Sec III, we study QZE the MQRM at finite temperature. We discuss the crossover of the survival probability decay rate from QZE to QAZE and its relation to the energy transport between the qubit and multiple cavity modes. Then, we extend the initial state of multiple cavity modes to squeezed thermal states and discuss the effect of the squeezing on the QZE in Sec IV. We close this paper with a short summary in Sec V.

II Model and Methodology

II.1 The multimode quantum Rabi model

In this paper, we focus on the MQRM which describes a qubit or a two-level system(TLS) interacting with multiple quantized photonic modes(e.g., cavity or resonator). The Hamiltonian of MQRM can be written asgely2017convergence ; malekakhlagh2017cutoff ($\hbar=1$):

[TABLE]

where $\sigma_{x},\sigma_{z}$ are standard pauli operators and $a_{m}^{\dagger}$,$a_{m}$ are the creation and annihilation boson operators for the $m$th mode with the frequency $\omega_{m}=(m+1)\omega_{0}$ and coupling strength $g_{m}=\sqrt{m+1}g$ to the TLS. $M$ is the total number of modes and depends on the specific physical realization. If not mentioned, we set $M=15$ throughout this study which is already feasible in the current circuit-QED experimentsundaresan2015beyond . Since we are interested in the resonant situation, $\Delta=\omega_{0}$ is assume in this study.

To study the QZE of the MQRM, we focus on the survival probability of the TLS under successive ideal projective measurements with operator $\sigma_{z}$. The initial state is assumed as the product state: the TLS is at the exited state $|\uparrow\rangle$ and cavity modes are squeezed thermal states:

[TABLE]

here $Z(\beta)=\text{Tr}e^{-\beta H_{B}}$ is the partition function. $\beta=1/(k_{B}T)$ where $T$ is the temperature and $k_{B}$ is the Boltzmann constant. $S(\xi)$ is the squeezing operator:

[TABLE]

where $\xi_{m}=r_{m}{e^{i\phi_{m}}}$ is the squeezing parameter for the $m$th mode. For simplicity, we assume a uniform squeezing for all $M$ modes, i.e. $\xi_{m}=\xi=re^{i\phi}$.

II.2 MPS-based numerical method

In the following, we briefly introduce the numerically exact method based on TFD and MPS to simulate the dynamics. We first remove the squeezing operator in the initial state (2) by applying the unitary transform $S(\xi)$ to the Hamiltonian (1):

[TABLE]

where $A=\cosh^{2}r+\sinh^{2}r,B=e^{i\phi}\cosh r\sinh r,K=\cosh r+e^{i\phi}\sinh r$.

According to the TFD methodsuzuki1991density ; gelin2017thermal , the evolution of this transformed Hamiltonian (4) from the initial state of thermalized cavity modes is equivalent to a Schrödinger equation of a modified Hamiltonian $\hat{H}$ defined in an enlarged Hilbert space(see Appendix):

[TABLE]

The the vacuum initial state is

[TABLE]

where $b_{m}^{\dagger},b_{m}$ are boson operators of the fictitious modes and $\theta_{m}=\text{arctanh}(e^{-\beta\omega_{m}/2})$ . The evolution of the expectation value of an arbitrary operator $O$ that affects in original Hilbert space can be straightforwardly calculated :

[TABLE]

Finally, the TFD Schrödinger equation is simulated by MPS-based numerically exact method schollwock2011density ; haegeman2016unifying ; wall2016simulating which has been widely used and proved to be highly efficient in solving quantum many body dynamicsschroder2016simulating ; kloss2018time .

II.3 Survival probability and effective decay rate

Now we turn to the QZE and QAZE of the MQRM. The QZE misra1977zeno can be described by the survival probability $P_{\text{sur}}(n\tau)$ which is defined as the probability of finding the initial state after $n$ successive measurements with equal time interval $\tau$. The measurement considered in this paper is assumed to be an ideal projecting measurements of operator $\sigma_{z}$ followed by a post selection regarding to the positive measurement result.The survival probability can be written as:

[TABLE]

where $P=|\uparrow\rangle\langle\uparrow|$ is measurement projecting operator. In the short interval time limit $\tau\rightarrow 0$, one can further write $P_{\text{sur}}(t=n\tau)$ in an exponentially decay formfacchi2001quantum ; maniscalco2006zeno :

[TABLE]

where $\gamma(\tau)$ is an effective decay rate:

[TABLE]

Eq.(10) is valid when the measurement interval $\tau$ is relatively short where the measurement disturbance to the environment(here is cavity modes) can be neglected thus the state after each projecting measurement collapse to the identical state as Eq. (6)he2017zeno . Therefore, we restrict $g\tau<1$ in the following discussion.

The effective decay rate $\gamma(\tau)$ is a crucial quantity to characterize the QZE and the QAZEzheng2008quantum ; cao2010dynamics : ${\partial\gamma(\tau)}/{\partial\tau}>0$ means that the system is more severely slowed-down by faster repeated measurements, indicating the occurrence of QZE; on the contrary, ${\partial\gamma(\tau)}/{\partial\tau}<0$ can be regarded as the characteristic of QAZE since the decay is accelerated by frequent measurements. Compared to the original criterion to classify QZE by using $\gamma(\tau)/\gamma_{0}<0$ where $\gamma_{0}=\gamma(\tau\rightarrow\infty)$ is the natural decay ratezheng2008quantum , new definitions through ${\partial\gamma}/{\partial\tau}$ retains the core physical picture of QZE and QAZE without calculating $\gamma_{0}$ that may not exist in some modelschaudhry2014zeno ; zhang2015zeno ; wu2017quantum ; he2017zeno . Throughout this paper, we use this new criterion to classify QZE and QAZE, the potential crossover point is denoted as $\tau_{c}$, i.e. {\partial\gamma(\tau)}/{\partial\tau}\big{|}_{\tau_{c}}=0 .

III Results and discussion

III.1 Thermal equilibrium initial state

We first consider the QZE when cavity modes are prepared as the thermal equilibrium state. The survival probability of Hamiltonian (5) with the product initial state (6) can be written aslizuain2010zeno :

[TABLE]

where $\rho_{0}=|\hat{\psi}(0)\rangle\langle\hat{\psi}(0)|$ and $[\hat{H},\rho(0)]_{n}=[\hat{H},[\hat{H},\rho(0)]_{n-1}]$. A commonly used approximation to obtain the decay rate of the survival probability is to keep terms in (11) up to $\tau^{2}$ , this leads to the decay of $P_{\text{sur}}(\tau)$ in a quadratical form:

[TABLE]

where $\tau_{z}$ is known as ”Zeno time”schulman1994characteristic . Thus for $N$ repeated projective measurement with equal interval $\tau$, the survival probability can be approximated as an exponential decay:

[TABLE]

with the effective decay rate $\gamma_{\text{2nd}}(\tau)$:

[TABLE]

Since Eq.(14) is independent of $\Delta$ and $\omega_{m}$ , its validity is limited to the case where the measurement interval $\tau$ is much shorter than the typical time scale of each mode , that is $\tau<1/\omega_{m}$. When high frequency modes are included, it severely precludes the availability of $\gamma_{\text{2nd}}(\tau)$ compared to the Rabi model where only single resonant mode is considered.

Instead, we employ an alternative method based on the thermofield dynamics(TFD)umezawa1982thermo ; suzuki1991density . This method transforms the evolution of a mixed thermalized initial state into another evolution with a pure initial state in an enlarged Hilbert space whose total number of degrees of freedom is double compared to the original one. Recently, this method has received much attention and has been widely used to study finite temperature dynamics of quantum electron-vibrational systemsborrelli2016quantum ; borrelli2017simulation ; chen2017finite and open quantum systemswang2017finite ; he2018zeno . Details of TFD can be found in Refborrelli2016quantum ; gelin2017thermal .

By restricting the evolution of the TFD Schrödinger equation to the single expiation subspace, we can obtain an analytic expression for the decay rate of the survival probability(see Appendix):

[TABLE]

where $\text{sinc(x)}\equiv\sin(x)/x$. Different from the $\gamma_{\text{2nd}}$ in Eq(14), $\gamma_{\text{th}}$ shows apparent dependence of frequencies $\omega_{k}$ and $\Delta$ and has a leading scaling of $g_{k}^{2}/(\Delta-\omega_{k})^{2}$ which is similar to the decay rate of the spontaneous emission in the multimode Percell effect due to non-resonant modeshouck2008controlling .

In Figure(1), the analytical result $\gamma_{\text{th}}(\tau)$ agree well with numerical exact ones in a large range of parameters including strong coupling and high temperature regimes. Specifically, for zero temperature case( Figure(1.a)) the decay rate $\gamma(\tau)$ for a single mode Rabi model shows a monotonic increasing with the increase of measurement interval $\tau$ for both weak($g/\omega_{0}=0.01$) and strong($g/\omega_{0}=0.2$) coupling strengthes, clearly demonstrates a pure Zeno effect. However, by including the high-lying photonic modes($M=15$), the decay rate of survival probability presents a nonmonotonic behavior, which increases with $\tau$ before $\omega_{0}\tau<0.2$ and then decreases, indicating a crossover from QZE to QAZE. Similar transitions can be observed for finite temperature situations, as shown in (Figure(1.b)). Thus we conclude that these high-frequency cavity modes in the MQRM may induce the QAZE. In the following, we show that such QAZE is attributed to the energy backflow from high-frequency modes to the TLS.

According to Eq.(10), the decay rate $\gamma(\tau)$ can be related to the population of the TLS in the excited state $|\uparrow\rangle$. Since the energy of the TLS is $E_{\text{TLS}}(t)=\Delta\langle\sigma_{z}(t)\rangle/2$ which is proportional to the energy of the TLS:

[TABLE]

thus we may understand this crossover from the view of energy transport between multiple mods and the TLS. Since the QAZE is classified by $d\gamma(\tau)d\tau<0$, this requires:

[TABLE]

By considering the assumption that the measurement interval is short ($g\tau<1$), $E_{\text{TLS}}(\tau)$ can be approximately expressed as:

[TABLE]

where $a_{1}>0$ since the TLS is initially prepared in exited state and $E_{\text{TLS}}$ is maximized at $\tau=0$. Inserting Eq.(18) into Eq.(17), we obtain:

[TABLE]

For arbitrary system, ${d\gamma}/{d\tau}|_{\tau=0}>0$ which implicates the pure Zeno effet when $\tau\rightarrow 0$, the occurrence of the QAZE requires $a_{2}<0$, i.e. ${d^{2}E_{\text{TLS}}}/d^{2}\tau>0$. As depicted in Figure(2).a, the crossover point $\tau_{c}$ is close to the pole of the energy transport rate as $d^{2}E_{\text{TLS}}/d^{2}\tau=0$. As shown in $(\ref{Fig2})$, by further calculating the evolution of the energy operator for each mode $E_{m}(\tau)=\omega_{m}\langle\psi(\tau)|a_{m}^{\dagger}a_{m}|\psi(\tau)\rangle$, we observe an energy backflow from the higher frequency parts of cavity modes to the TLS during the crossover point $\tau_{c}$. It is this energy backflow induced by high frequency modes that leads to $d^{2}E_{\text{TLS}}/d^{2}\tau>0$ and makes it possible to observe QAZE. On the contrary, the low frequency part of modes always absorbs energy from the TLS during $g\tau<1$, leading to the pure QAZE in the single mode Rabi model as shown in Figure(1).a.

III.2 Squeezed thermal initial states

In this section, we generalize the initial state of the cavity modes to the squeezed thermal states.

From the view of the Hamiltonian (5), effects of the squeezing come from two aspects: 1).The renormalization of mode frequencies and effective coupling strengthes by factors $A$ and $K$ respectively. 2).The nonharmonic terms in $\hat{H}_{B}$ with the factor $B$.

The effect of 1) can be included in the analytic expression (15) by replacing $g_{m}^{2}\rightarrow|K|^{2}g_{m}^{2}$ and $\omega_{m}\rightarrow A\omega_{m}$. By increasing $r$, the effective coupling strength is increased, leading to a larger decay rate. Moreover, due to the factor $|K|^{2}$, $\gamma(\tau)$ depends on the squeezing angle $\phi$. Since $|K|^{2}=\cosh(2r)+\cos\phi\sinh(2r)$, the critical angle where $\gamma(\tau)$ takes the extreme value may be $\phi_{\text{max}}=0$ and $\phi_{\text{min}}=\pi$ respectively. Same critical angles have been observed on the single-mode rabi model with initial squeezed fieldlizuain2010zeno and the TLS in a squeezed bath under RWAmundarain2006zeno . However, this is not the case in the current MQRM where effects of both multimodes and non-rotating terms are taken into consideration. Specifically, by increasing the number of modes from $M=1$ to $M=15$, critical angles $\theta_{\text{min}}$ and $\theta_{\text{max}}$ gradually increase with a fixed angle difference $\theta_{\text{min}}-\theta_{\text{max}}=\pi$ as shown in Fig(3.a). Interestingly, as depicted Fig(3.b), such critical angles are independent of squeezing strength $r$ and coupling strength $g$. This independence agrees with results of the TLS in squeezed bathmundarain2006zeno . Surprisingly, the scaling of $\gamma(\tau)\sim g^{2}$ for thermal equilibrium initial states shown in Eq(15) may still holds for squeezed thermal bath by noting that $\gamma(\tau)/g^{2}$ as a function of $\phi$ for different coupling strengthes $g$ in Fig(3.b) collapses to the same curve.

The measurement interval $\tau$ also affects such critical angles, as shown in Fig(3.c). Note when $\tau\rightarrow 0$, the shift of critical angles approaches zero. The crossover from QZE to QAZE can still be observed for the squeezed thermal initial state. In addition, a specific squeezing angle can either highlights or downplay this transition: as shown in Fig(3.c), during transition point $\tau_{c}$, $\gamma(\tau)$ are significantly increased for $\phi=\phi_{\text{max}}$ and suppressed for $\phi=\phi_{\text{min}}$, while with the increase of $\tau$, the difference between $\gamma(\tau)$ at two critical angles is narrowed.

Finally, the effect of temperature on the decay rate for the squeezed thermal initial state is described by $\theta_{m}$ in $\hat{H}_{\text{int}}$ where higher temperature leads to larger effective coupling strength and consequently larger decay rate, as shown Fig(4.b). Moreover, increasing temperature also causes reduction of the shift of critical angles. This is not a surprising result. As explained above, the shift of critical angles are related to the $Ba_{m}^{{\dagger},2}+B^{*}a_{m}^{2}$ in $H_{B}$. However, with the increase of temperature, $\hat{H_{\text{int}}}$ becomes dominant compared to the nonharmonic term, thus can be neglected in the high-temperature limit and the relation between $\gamma(\tau)$ and $\phi$ are only determined through $|K|^{2}$.

IV Conclusion

In summary, we study the QZE and QAZE of the MQRM where cavity modes are prepared in thermal equilibrium state. We derive an analytic expression of the decay rate of the survival probability which shows a transition from QZE to QAZE. By calculating the evolution of energy flow between the TLS and cavity modes, we show that such QZE-QAZE transition is related to the energy backflow from high frequency modes to the TLS. Furthermore, we generalize the initial state of cavity modes to squeezed thermal states and study effects of squeezing angle and strength by applying a numerically exact method. We find that the decay of survival probability is accelerated by non-vanishing squeezing strengthes. Moreover, squeezing angles also significantly affect the decay rate: compared to $\theta_{\text{max}}=0$ and $\theta_{\text{min}}=\pi$ in which decay rates take maximal and minimal values in the single-mode Rabi model, high frequency modes in the MQRM cause a positive shift on these critical squeezing angles. Finally, with the increase of the temperature, the decay of the survival probability is accelerated and the shift of critical squeezing angles are reduced.

V Acknowledgement

Shu He is supported by the National Natural Science Foundation of China under Grant No. 11804240. Chen Wang is supported by the National Natural Science Foundation of China under Grant Nos. 11704093 and 11547124. Liwei Duan and Qing-Hu Chen are supported by the National Natural Science Foundation of China under Grant Nos.11674285 and 11834005.

Appendix A A Brief introduction to TFD

We briefly introduce the theory of thermofield dynamics in this Appendix.

The time evolution of an arbitrary Hamiltonian $H$ at finite temperature can be described by the Liouville-von Neumann equation($\hbar=1$):

[TABLE]

where the initial state is assumed to be at thermal equilibrium:

[TABLE]

here $Z(\beta)=\text{Tr}e^{-\beta H}$ is the partition function and $\beta=\left(k_{B}T\right)^{-1}$. $k_{B}$ is the Boltzmann constant and $T$ is the temperature.

According to Refsuzuki1986thermo , the time evolution of such density matrix (20) can be equivalently described by a modified Hamiltonian $\hat{H}=H-\tilde{H}$ and corresponding Schrödinger equation defined in an enlarged Hilbert space($\hbar=1$):

[TABLE]

$\tilde{H}$ represents the fictitious Hamiltonian which can be derived from the original Hamiltonian $H$suzuki1991density .

By defining the vector $|I\rangle$:

[TABLE]

where $|k\rangle$($|\tilde{k}\rangle$) are arbitrary basis vectors of the original(fictitious) space, the density matrix $\rho(t)$ and wave function $|\psi(t)\rangle$ have the following relationsuzuki1985thermo :

[TABLE]

Particularly, the initial state at thermal equilibrium is:

[TABLE]

The expectation value for an arbitrary operator $A$ defined in the original Hilbert space $\{|k\rangle\}$ can be calculated by:

[TABLE]

Thus the evaluation of $\langle A(t)\rangle$ by using the TFD wave function $|\psi(t)\rangle$ is equivalent to the that by using the density matrix $\rho(t)$.

Now let us turn to the Hamiltonian of the MQRM (1). To describe its evolution from the initial state $(\ref{IniOrigSq})$ under the framework of the TFD, we first remove the squeezing operator in the initial state by applying a unitary transform $S(\xi)$ defined in (3) to the Hamiltonian (1):

[TABLE]

where $A=\cosh^{2}r+\sinh^{2}r,B=e^{i\phi}\cosh r\sinh r,K=\cosh r+e^{i\phi}\sinh r$ , this is the Hamiltonian (4) in the main text.

We can choose the fictitious Hamiltonian $\tilde{H}$ as :

[TABLE]

where $b_{m}^{\dagger}$($b_{m}$) is the creation (annihilation) boson operator defined in the fictitious Hilbert space corresponding to the the mode $a_{m}$. This leads to the total Hamiltonian for the TFD schrödinger equation:

[TABLE]

The corresponding initial TFD wave function can be written as:

[TABLE]

where $|\textbf{0}(\beta)\rangle$ is often referred to as the thermal vacuum state:

[TABLE]

This thermal vacuum state can be further transformed to the ground state(vacuum state) of the modes $\{a_{m}\}$ and $\{b_{m}\}$ by applying the so-called Bogoliubov thermal transformationtakahashi1996thermo :

[TABLE]

where $|\textbf{0}\rangle=\prod_{m=0}^{M-1}|0\rangle_{a_{m}}|0\rangle_{b_{m}}$ and the transformation operator $G$ is defined as:

[TABLE]

where

[TABLE]

Instead of solving the TFD schrödinger equation (29) with the initial state (31) , one can apply the inverse Bogoliubov thermal transformation to the Hamiltonian $\hat{H}$. By considering the following new relations:

[TABLE]

We have:

[TABLE]

$\hat{H}_{\theta}=e^{iG}\hat{H}e^{-iG}$ is the final Hamiltonian (5) in the main text:

[TABLE]

and corresponding initial state (6):

[TABLE]

Appendix B Derivation of Eq.(15)

In this Appendix, we show the detailed derivation to the analytic expression of the decay rate $\gamma_{\text{th}}(\tau)$.

The evolution of the MQRM from a product initial state where resonator modes are at thermal equilibrium is govern by the TFD Schrödinger equation of the Hamiltonian (37) with $r=0$(without squeezing):

[TABLE]

To obtain an analytic expression, we further apply an approximation by restricting the evolution in the single excitation(SE) Hilbert space. This can be achieved by writing the wave function at arbitrary time $t$ as:

[TABLE]

Thus the Hamiltonian (39) can be simplified as:

[TABLE]

where $\sigma_{\pm}=\frac{1}{2}\left(\sigma_{x}\pm i\sigma_{y}\right)$.

Inserting the wave function (40) into the TFD Schrödinger equation, we have:

[TABLE]

applying the following transformation:

[TABLE]

we have:

[TABLE]

Integrating last two equations in (44) and substituting into the first one, we have:

[TABLE]

The solution of the integro-differential equation above can be obtained iteratively. Since we are limited to the frequent measurement limit where $g\tau<1$, only the short time behavior of $\tilde{\chi}(t)$ is concerned. This can be obtained by taking first iteration and inserting the initial condition $\chi(0)=1$:

[TABLE]

where

[TABLE]

Substituting (46) into the definition of the decay rate of the survival probability (10) , we finally obtain:

[TABLE]

where $\text{sinc(x)}\equiv\sin(x)/x$.

Appendix C MPS-based numerically exact method

As derived in the previous Appendix, the evolution of the MQRM from a squeezed thermal states can be described by the TFD schrödinger equation of the Hamiltonian (5) with the initial state (6). To implement MPS-based numerical method, the Hamiltonian is rewritten in a compact form of the following matrix product operator(MPO)wall2016simulating :

[TABLE]

where $I$ is the identity operator of the TLS and the boson space

[TABLE]

Note for the open boundary condition, only first row of $\mathcal{W}_{b}^{[1]}$ and first column of $\mathcal{W}_{a}^{[M]}$ are used. With this MPOs, we employ the recently proposed time-dependent variational principle(TDVP) methodhaegeman2016unifying to simulate the evolution. For numerical results throughout this study, the convergence can be achieved by setting the maximal bound dimensions of the MPS $D_{\text{max}}\leq 15$ and the cutoff boson number for each cavity mode $N_{\text{max}}\leq 80$. The time step of simulations throughout this paper is set to $g\delta t=10^{-3}$.