Frequency-renormalized multipolaron expansion for the quantum Rabi model

TL;DR

This paper introduces a frequency-renormalized multipolaron expansion method to accurately analyze the ground state of the quantum Rabi model, improving upon previous approaches with variationally optimized parameters.

Contribution

It extends the polaron and coherent state expansion techniques by incorporating frequency renormalization, enhancing the description of the quantum Rabi model's ground state.

Findings

High efficiency in describing the QRM ground state.

Potential applicability to more complex light-matter interaction models.

Improved accuracy over previous methods.

Abstract

We present a frequency-renormalized multipolaron expansion method to explore the ground state of quantum Rabi model (QRM). The main idea is to take polaron as starting point to expand the ground state of QRM. The polarons are deformed and displaced oscillator states with variationally determined frequency-renormalization and displacement parameters. This method is an extension of the previously proposed polaron concept and the coherent state expansion used in the literature, which shows high efficiency in describing the physics of the QRM. The proposed method is expected to be useful for solving other more complicated light-matter interaction models.

Equations

$$H=\frac{\Omega }{2}{\sigma }_x+\omega {\hat{a}}^{†}\hat{a}+g{\sigma }_z({\hat{a}}^{†}+\hat{a}),$$

$$H=\sum _{S_z=\pm }(h^{S_z}|S_z⟩⟨S_z|+\frac{\Omega }{2}|S_z⟩⟨\overline{S_z}|)+{\epsilon }_0,$$

$$|G⟩=\frac{1}{\sqrt{2}}({\Psi }^{+}(x)|{+}_z⟩-{\Psi }^{-}(x)|{-}_z⟩),$$

$${\Psi }^{+}(x)={\Psi }^{-}(-x).$$

$${\Psi }^{\pm }(x)=\sum _{n=1}^N{C}_n{\phi }_n^{\pm }(x),$$

$$|G⟩=\frac{1}{\sqrt{2}}\sum _{n=1}^N{C}_n({\varphi }_0({\xi }_n\omega ,x+{\zeta }_n{g}^{\prime })|{+}_z⟩$$

$$-{\varphi }_0({\xi }_n\omega ,x-{\zeta }_n{g}^{\prime })|{-}_z⟩),$$

$$\begin{split}|G\rangle=\frac{1}{\sqrt{2}}\sum_{n=1}^{N}C_{n}\big{(}\varphi^{+}_{n}|+_{z}\rangle-\varphi^{-}_{n}|-_{z}\rangle\big{)}.\end{split}$$

$$\begin{split}E_{FR-MPE}&=\langle G|H|G\rangle\\ &=\frac{1}{2}\sum_{n,m}^{N}C_{n}C_{m}\bigg{[}\big{(}\langle\varphi^{+}_{n}|h^{+}|\varphi_{m}^{+}\rangle+\langle\varphi^{-}_{n}|h^{-}|\varphi^{-}_{m}\rangle\big{)}\\ &-\frac{\Omega}{2}\big{(}\langle{\varphi^{+}_{n}}|{\varphi^{-}_{m}}\rangle+\langle{\varphi^{-}_{n}}|{\varphi_{m}^{+}}\rangle\big{)}\bigg{]}-\frac{1}{2}\omega({g^{{}^{\prime}}}^{2}+1)\\ &=\sum_{n,m}^{N}C_{n}C_{m}\langle\varphi^{+}_{n}|h^{+}|\varphi_{m}^{+}\rangle-\frac{\Omega}{2}\sum_{n,m}^{N}\langle\varphi^{+}_{n}|\varphi^{-}_{m}\rangle\\ &-\frac{1}{2}\omega({g^{{}^{\prime}}}^{2}+1).\\ \end{split}$$

$$\begin{split}h_{nm}^{+}&=\langle\varphi^{+}_{n}|h^{+}|\varphi^{+}_{m}\rangle\\ &=\frac{1}{2}\omega\langle\varphi_{n}^{+}|\big{(}\hat{p}^{2}+(\hat{x}+g^{{}^{\prime}})^{2}\big{)}|\varphi_{m}^{+}\rangle\\ &=\frac{1}{2}\omega\langle\varphi_{n}^{+}|(-\frac{\partial^{2}{}}{\partial{x^{2}}}+{x}^{2}+2{x}{g^{{}^{\prime}}}+{g^{{}^{\prime}}}^{2})|\varphi_{m}^{+}\rangle.\\ \end{split}$$

$$\begin{array}{rl}⟨{\phi }_n^{+}|(-\frac{{\partial }^2}{\partial x^2})|{\phi }_m^{+}⟩=& -⟨{\phi }_n^{+}|(4D^2{x}^2+4DEx+E^2\\ & +2D)|{\phi }_m^{+}⟩,\end{array}$$

$$D$$

$$h_{nm}^{+}$$

$$(2g^{{}^{\prime }}-4DE)x-E^2-2D+{g^{{}^{\prime }}}^2)|{\phi }_m^{+}⟩,$$

$$S_{nm}$$

$$\displaystyle=\sqrt{2}\big{[}\frac{\xi_{n}\xi_{m}}{(\xi_{n}+\xi_{m})^{2}}\big{]}^{\frac{1}{4}}e^{\big{(}-\frac{(\zeta_{n}-\zeta_{m})^{2}{g^{{}^{\prime}}}^{2}\xi_{n}\xi_{m}}{2(\xi_{n}+\xi_{m})}\big{)}},$$

$${⟨\hat{x}⟩}_{nm}$$

$$=S_{nm}\frac{-{\xi }_m{\zeta }_m-{\xi }_n{\zeta }_n}{{\xi }_n+{\xi }_m}{g}^{{}^{\prime }},$$

$${⟨{\hat{x}}^2⟩}_{nm}$$

$$\displaystyle=S_{nm}[\frac{1}{\xi_{n}+\xi_{m}}+\big{(}\frac{{(\xi_{m}\zeta_{m}+\xi_{n}\zeta_{n}){g^{{}^{\prime}}}}}{\xi_{n}+\xi_{m}}\big{)}^{2}].$$

$$S_{n\overline{m}}$$

$$\displaystyle=\sqrt{2}\big{[}\frac{\xi_{n}\xi_{m}}{(\xi_{n}+\xi_{m})^{2}}\big{]}^{\frac{1}{4}}e^{\big{(}-\frac{(\zeta_{n}+\zeta_{m})^{2}{g^{{}^{\prime}}}^{2}\xi_{n}\xi_{m}}{2(\xi_{n}+\xi_{m})}\big{)}},$$

$$⟨G|H|G⟩$$

$$\displaystyle=\sum_{n,m}^{N}C_{n}C_{m}\big{(}\langle{\psi^{+}_{n}}|{\psi^{+}_{m}}\rangle\big{)}$$

$$=1.$$

Full text

I Introduction

Quantum Rabi model (QRM) describes a two-level system interacting with a single-mode bosonic field Rabi (1936, 1937). It plays a fundamental role in many fields of physics, such as superconducting circuit quantum electrodynamics (QED) Niemczyk et al. (2010); Wallraff et al. (2004, 2005), quantum optics Raimond et al. (2001); Scully and Zubairy (1997), quantum information Ashhab and Nori (2010); Zhou et al. (2016); Zhao et al. (2016), and quantum computation Romero et al. (2012), condensed matter physics Holstein (1959).

Experimentally, the model has been first realized in the cavity QED systems Xiang et al. (2013), in which the coupling strength is quite weak, corresponding to the so-called weak coupling regime. In this regime, the rotating-wave approximation (RWA) has been widely employed, which leads to an analytically solvable model, namely, the Jaynes-Cummings (JC) model Jaynes and Cummings (1963). The JC model is a basic model in quantum optics which is successful in the understanding of a range of experimental phenomena, such as the well known quantum Rabi oscillation Brune et al. (1996) and vacuum Rabi mode splitting Thompson et al. (1992).

Recently, with the advancement of quantum technology Leibfried et al. (2003); Englund et al. (2007), the so-called strong coupling Wallraff et al. (2004), ultra-strong coupling Peropadre et al. (2010); Niemczyk et al. (2010); Forn-Díaz et al. (2010); Andersen and Blais (2017) and even the deep strong coupling Casanova et al. (2010) regimes have been experimentally realized in many devices. As a result, the RWA widely used in the literature is no longer valid in these strongly coupling regimes, and thus a full QRM has to be reconsidered in order to describe well the physics observed in these strongly coupling regimes.

It turns out that, despite its simple form, it is not an easy task to fully solve and understand the QRM. Therefore, many approximate methods including adiabatic approximation Irish et al. (2005), general rotating-wave approximation (GRWA) Irish (2007) and its extensions Zhang et al. (2011); Zhang (2016); Yu et al. (2012), unitary transformationYan et al. (2015), the variational technique Liu et al. (2015), to name just a few, have been proposed. However, it has been shown that each of these approximate methods may be valid in certain limited regime, and an approximate method which is valid in whole parameter regime of the model is still favorable and deserves efforts to develop and improve. In 2011, a remarkable mathematical progress on the integrability of the QRM has been obtained by Braak Braak (2011) and thus the exact spectra of the QRM have been determined in an analytical way. The exact spectra of the QRM has also been formulated by using Bogoliubov operator technique Chen et al. (2012). However, in order to explore the full physics of the QRM, it is obviously not enough to only know the spectra of the model, one still needs to know exactly the wavefunction of the model. Therefore, it is still important to explore a simple and straightforward method to study the QRM. On the one hand, this method should be valid in whole parameter regime of the QRM, on the other hand, it should also be convenient to formulate the wavefunction of the QRM. This is our motivation of the present work.

Very recently, we have introduced a trial wavefunction based on the concept of polaron and anti-polaron picture in order to explore the phase diagram of the QRM Ying et al. (2015). An important feature of this trial wavefunction is that it provides a unified framework to accurately describe the physics of the QRM both in the weak- and strong coupling regimes. However, in the crossover regime, some errors relatively are still not negligible, particularly, at a low oscillator frequency, as shown in Fig. 1. Therefore, some improvements are still desirable in order to capture more accurate physics in the crossover regime.

How to further improve the calculation accuracy in the crossover regime, particularly, in the low oscillator frequency case? One notes that a variational coherent-state expansion method has been proposed by Bera et al. Bera et al. (2014a, b) in the study of the spin-boson model. It was shown that the more the polarons are used, the more accurate the result becomes. Following the idea of the multipolaron expansion, here we propose a variational frequency-renormalized multipolaron expansion method to improve the performance of the trial wavefunction based on the polaron and anti-polaron picture. The key difference is that in contrast to Bera et al.’s multipolaron expansion, we introduce the frequency renormalization feature. As shown later, the frequency renormalization introduced shows a high efficiency in calculating the energy and wavefunction of the QRM. Therefore, it is expected that our frequency-renormalized multipolaron expansion method is useful in solving more complicated models related to light-matter interaction.

The paper is organized as follows. In Sec. II, the QRM is introduced. In Sec. III, we construct our variational method based on the frequency-renormalized multipolaron expansion as the trial ground state wavefunction. In Sec. IV, we present the results based on the proposed method, and compare them with those obtained by the multipolaron expansion without introducing frequency-renormalized feature. Sec. V is devoted to a brief conclusion.

II The QRM model

Following the notation in Refs. Albert (2012); Irish and Gea-Banacloche (2014); Ashhab and Nori (2010); Zhang et al. (2011, 2013); Liu et al. (2013), the Hamiltonian of the QRM model reads ($\hbar=1$):

[TABLE]

where $\Omega$ is the qubit energy level splitting, $\sigma_{x,z}$ is the Pauli matrix to describe the qubit, $\hat{a}^{\dagger}$ and $\hat{a}$ are the bosonic creation and annihilation operators, respectively, of the bosonic mode with frequency $\omega$, and $g$ denotes the coupling strength between the qubit and the bosonic mode.

In terms of the quantum harmonic oscillator with dimensionless formalism $\hat{a}^{\dagger}=(\hat{x}-i\hat{p})/\sqrt{2}$, $\hat{a}=(\hat{x}+i\hat{p})/\sqrt{2}$, where $\hat{x}=x$ and $\hat{p}=-i\frac{\partial}{\partial x}$ are the position and momentum operators respectively, the model can be rewritten as Irish and Gea-Banacloche (2014)

[TABLE]

where, $S_{z}=\pm_{z}$ represents the eigenvalue of $|\pm_{z}\rangle$ in z-direction, $\bar{S_{z}}=-{S_{z}}$, $h^{\pm}={\omega}(\hat{p}^{2}+(\hat{x}\pm g^{\prime})^{2})/2$, $\varepsilon_{0}=-\omega({g^{\prime}}^{2}+1)/2$ is a constant, and $g^{\prime}=\sqrt{2}g/\omega$. For simplicity, we take $\Omega=1$ as the units of energy here and hereafter.

III Frequency-renormalized multipolaron expansion method

Consider the parity operator $\Pi=\sigma_{x}e^{a^{\dagger}a}$, one has $[H,\Pi]=0$. Such a $\mathbb{Z}_{2}$ symmetry leads to a decomposition of the state space into just two subspaces with odd and even parity, respectively. Since the ground state has an odd parity, one takes the trial wavefunction of the QRM in the position representation as

[TABLE]

where $\Psi^{\pm}$ is the wavefunction associated with the spin state $|\pm_{z}\rangle$, and we have

[TABLE]

Based on the polaron picture Ying et al. (2015); Bera et al. (2014a, b), they can be expanded as

[TABLE]

where $C_{n}$ is the coefficient and $N$ is the number of polaron used. Here $\varphi_{n}^{\pm}(x)=\phi_{0}(\xi_{n}\omega,x{\pm}\zeta_{n}g^{\prime})$ denotes the $n$th polaron which is given by deformed oscillator ground state wavefunction $\phi_{0}(\omega,x)$ with the frequency renormalization parameter $\xi_{n}$ and shifted position parameter $\zeta_{n}$.

Thus, Eq. (4) can be rewritten as

[TABLE]

which is the starting point of the present work. Due to the deformed polaron introduced, we call Eq. (III) as frequency-renormalized multi-polaron expansion (FR-MPE). Actually, if one takes $N=2$, Eq. (III) recovers the previous polaron and anti-polaron wavefunction which has been used to explore the ground state phase diagram of the Rabi model Ying et al. (2015). On the other hand, if one takes $\xi_{n}=1$, i.e., the frequency-renormalization factor is not considered, Eq.(III) is the single mode version of the coherent state expansion used in the study of the spin-boson model Bera et al. (2014a, b), since a coherent state is a displaced oscillator state in the $\hat{x}$ representation.

The variational parameters introduced in Eq.(III), i.e., $C_{n}$, $\xi_{n}$ and $\zeta_{n}$ ($n=1,\cdots,N$), can be determined by minimizing the ground state energy $E_{G}=\langle G|H|G\rangle$ (see Appendix A for a detailed derivation), subject to the constraint of wavefunction normalization $\langle G|G\rangle=1$, under which the number of variational parameters will be $3\times N-1$. In order to determine the variational parameters, we first adopt simulated annealing algorithm Kirkpatrick (1984); Hwang (1988) to search the rough values of the variational parameters. Then we use pattern search algorithm Hooke and Jeeves (1961); Davidon (1991) to refine these variational parameter values in order to further minimize the ground state energy. The combination of these two algorithms is found to be sufficient to determine the variational parameters with high efficiency and high precision.

IV numerical results and discussion

IV.1 The result with $N=4$

First of all, we present some results to show the high precision of our method. As shown in Fig. 1, the lower the oscillator frequency is, the larger the error near $g_{c}$ is. Moreover, the main error is located around $g_{c}$. Therefore, here we consider the case of $\omega=0.01$ and take $N=4$. Meanwhile, we compare the obtained results with those with numerical exact diagonalization. In Fig. 2 we show various physical quantities including (a) the ground state energy, (b) the spin polarization $\langle\sigma_{x}\rangle$, (c) the correlation $\langle\sigma_{z}(a^{\dagger}+a)\rangle$ and (d) the mean photon number $\langle a^{\dagger}a\rangle$ as a function of the coupling strength. Due to the afore-mentioned reason, hereafter we limit ourself to the region around $g_{c}$. It is found that the agreement is quite good, which will be further discussed in the next subsection. This result confirms the high precision of our method, even in the low oscillator frequency regime. In addition, our method also exhibits high efficiency since only two pairs of polaron and anti-polaron ($N=4$) have been considered here, which indicates that the deformed polaron picture Ying et al. (2015) is a good starting point to capture the physics of the QRM.

IV.2 The high efficiency of the frequency-renormalized multipolaron expansion

In order to confirm the high efficiency of the frequency-renormalized multipolaron expansion method, we compare the result with $N=4$ with that of $N=2$ with $\omega=0.01$ in Fig. 3. As also mentioned above, the case of $N=2$ recovers the previous polaron and anti-polaron picture Ying et al. (2015). From Fig. 3 it is noted that the results for $N=4$ have a vanishing small error in comparison to that of $N=2$. Technically, we just increase an additional pair of polaron and anti-polaron beyond that of $N=2$, leading to a dramatic improvement, which shows the high efficiency of the polaron and anti-polaron basis. Physically, it is not difficult to understand why the polaron and anti-polaron basis is so high efficient. This is because that the polaron and anti-polaron as a basic ingredient in describing the ground state wavefunction of the QRM is able to capture the essential physics of the model. The anti-polaron component originates naturally from the tunneling feature Irish and Gea-Banacloche (2014); Ying et al. (2015) of the QRM. By the same reason, the additional pair of polaron and anti-polaron also originate from the high-order effect of the tunneling feature, which reflects the many-body effect in the QRM. This point is more clear in the discussion of the ground state wavefunction below.

Figure 4 shows the ground state wavefunction of the QRM for $N=4$ and its polaron components. As a comparison, we also present the result for $N=2$. For simplicity, we only provide the $\Psi^{+}$-component, and the $\Psi^{-}$ has a similar behavior due to the odd parity. From Fig. 4(a), although the case of $N=2$ captures the nature of the wavepacket separation for $g/g_{c}=1.05$, the error from the numerical exact result is still obvious, as shown in Fig. 4(c) as yellow dashed line. In particular, in the region around $x=0$, the error is more obvious. This is because that after the wavepacket becomes separated, a pair of polaron and anti-polaron is not sufficient to describe the region away from the positions of the polaron and antipolaron, for example, the region of $x=0$. Therefore, according to the idea of the frequency-renormalized multipolaron expansion, one increases an additional pair of polaron and anti-polaron, not only the region away from the positions of the polaron and anti-polaron, the whole wavefunction has a vanishing small error in comparison with the exact one, as shown in Fig. 4(c) as blue dash-dotted line. This result indicates that the frequency-renormalized multipolaron expansion is highly efficient in calculating the wavefunction of the QRM, as a result, is able to calculate accurately the physical observables, as shown in Fig. 3.

IV.3 The importance of the frequency renormalization

Based on the scheme of multipolaron expansion introduced by Bera et al. in their coherent state expansion (CSE) method Bera et al. (2014a, b), one of the important features of our method is to introduce the frequency renormalization feature. The physical background of the introduction of the frequency renormalization is the change of the effective potentials induced by the tunneling between these two energy levels, which is the origin of the anti-polaron. The existence of the effective potentials would, of course, modify the frequency of the oscillator states in the QRM. In order to show the importance of the frequency renormalization, in the following we compare our method with the frequency renormalization with the simple multipolaron expansion without the frequency renormalization. The latter is called the CSE method below, which has been used to study the spin-boson model. Here we employ it to the QRM.

Firstly, we consider the ground state wavefunction for $\Psi^{+}$-component obtained by these two methods, as shown in Fig. 5, and also present the numerical exact wavefunction for comparison. It is quite obvious that the result with the frequency renormalization agrees much better with the exact one than that without the frequency renormalization. The comparison indicates that the frequency renormalization is indeed an important feature in describing the ground state physics of the QRM. This feature is even more important than simply increasing the number of the polaron and anti-polaron. This can be clearly seen from Fig. 6, showing the comparison of the ground state energy and the other physical observables between the results with frequency renormalization for $N_{FR-MPE}=2,4$ and without frequency renormalization for $N_{CSE}=4,6$, respectively. Obviously, our method has much higher efficiency than that without the frequency renormalization. Moreover, the corresponding variational parameters in our scheme are less than or equal to those of the CSE.

V Conclusions

Based on the polaron and anti-polaron picture, we proposed a frequency renormalization multipolaron expansion method to improve significantly the ground-state wavefunction of QRM, as a result, the ground-state energy and other physical quantities have also been significantly improved, in particular, near the crossover region. In comparison to the coherent state expansion given by Bera et al. Bera et al. (2014a, b), our method shows higher efficiency, in which the frequency renormalization plays an important role. Physical origin of the frequency renormalization is due to the deformed potential induced by the tunneling between two energy levels. This method can be applied to other more complicated quantum models like the spin-boson model Bera et al. (2014a, b), multi-qubit QRM Benatti et al. (2010); Mao et al. (2016), anisotropic QRM Xie et al. (2014), and biased or asymmetric Rabi model Liu et al. (2017) and so on. The further works are in progress.

Acknowledgments

We acknowledges useful discussion with Gao-Yang Li, Fu-Zhou Chen, Chong Chen and Chen Cheng. This work was supported by NSFC (Grants No. 11325417, No. 11674139) and PCSIRT (Grant No. IRT-16R35). ZJY also acknowledges the financial support of the Future and Emerging Technologies (FET) Programme within the Seventh Framework Programme for Research of the European Commission, under FET-Open Grant No.: 618083 (CNTQC).

Appendix A The ground state energy of FR-MPE

In this appendix, we present the main steps to get the ground state energy. The ground state $|G\rangle$ of FR-MPE is

[TABLE]

The ground state energy of FR-MPE is given by

[TABLE]

A.0.1 Calculation of each term

For the first term in Eq. (8)

[TABLE]

For simplicity we have assumed the unit $\hbar=m=1$.

We first give the first term in Eq. (9).

[TABLE]

where, the coefficients $D$ and $E$ are introduced for the simplicity of the formulation, they are defined as

[TABLE]

So, the expression of Eq. (9) is

[TABLE]

in which,

[TABLE]

[TABLE]

and

[TABLE]

Until now, we can get the first term in Eq. (8) completely.

For the second term in Eq. (8)

[TABLE]

So we finally can get the ground state energy.

A.0.2 Normalization condition

Besides of the above formulation, we still have the normalization condition which describes the relationships between the parameters.

[TABLE]

The derivation of the ground state energy of CSE and FR-MPE are nearly the same. Since in the $\hat{x}$ representation, a coherent state is a displaced ground state of the oscillator. We can get the CSE result by simply setting frequency renormalization factor $\xi=1$ in FR-MPE. So here we present the derivation of FR-MPE result only.