Numerical approach for the evolution of spin-boson systems and its application to the Buck【-逻*辑*与-】n
本站小编 Free考研考试/2022-01-02
Xue-Ying Liu1,2, Xue-Zao Ren2, Chen Wang1, Xian-Long Gao,1,4, Ke-Lin Wang31Department of Physics, Zhejiang Normal University, Jinhua 321004, China 2School of Science, Southwest University of Science and Technology, Mianyang 621010, China 3Department of Modern Physics, University of Science and Technology of China, Hefei 230026, China
First author contact:4Author to whom any correspondence should be addressed. Received:2020-02-19Revised:2020-03-15Accepted:2020-03-20Online:2020-05-18
Abstract We study the evolution properties of spin-boson systems by a systematic numerical iteration approach, which performs well in the whole coupling regime. This approach evaluates a set of coefficients in the formal expression of the time-dependent Schrödinger equation by expanding the initial state in Fock space. This set of coefficients is unique for the spin-boson Hamiltonian studied, allowing one to calculate the time evolution from different initial states. To complement our numerical calculations, we apply the method to the Buck–Sukumar model. We find that when the ground-state energy of the model is unbounded and no ground state exists in a certain parameter space, the time evolution of the physical quantities is naturally unstable. Keywords：Buck–Sukumar model;spin-boson systems;light-atom coupling;unbounded ground state;time-dependent Schrödinger equation
PDF (1003KB)MetadataMetricsRelated articlesExportEndNote|Ris|BibtexFavorite Cite this article Xue-Ying Liu, Xue-Zao Ren, Chen Wang, Xian-Long Gao, Ke-Lin Wang. Numerical approach for the evolution of spin-boson systems and its application to the Buck–Sukumar model. Communications in Theoretical Physics, 2020, 72(6): 065502- doi:10.1088/1572-9494/ab8a0f
1. Introduction
The interaction between light and matter is a central and fundamental problem in quantum optics. The simplest light–matter interaction model is the Jaynes–Cummings (JC) model, proposed by Jaynes and Cummings within the rotation-wave approximation (RWA). The JC model describes the interaction between a two-level atom and a single-mode field [1], which is an exactly solvable model with applications in many fields [2–6]. The JC model is experimentally realized by the squeezed radiation field in a 85Rb-atom micromaser [7] and by the collapse and revival of atomic inversion in a 138Ba-atom microlaser [8].
The JC model has been generalized in different ways. For example, the Dicke model (also known as the Tavis–Cummings model) couples the single mode with multiple two-level atoms. The coupling of multi-modes and multi-atoms can be constructed by Yang–Baxter algebra and is integrable through the Bethe ansatz [9]. An extra counter-rotating term in the JC model gives the famous Rabi model, which has been only analytically solved by Braak [10] till date. The JC model can also be generalized by including nonlinear couplings between a single atom and the radiation field [11–14]. An example is the intensity-dependent JC model, which exhibits significantly different dynamics in the absence and presence of RWA, owing to the dramatically enhanced field-squeezing effect in the latter case [15]. The Buck–Sukumar (BS) model [16] is a more general JC model considering the interaction of a two-level system with a field under intensity-dependent coupling, which is still analytically solvable. In contrast, the BS model with the counter-rotating term (sometimes called the nonlinear Rabi model) is non-analytic. Although the nonlinear term in the intensity-dependent JC or the BS model involves multiphoton interactions that are not physically realizable in current cavity- or circuit-QED setups, they are still of interest to the quantum optics and cold-atom fields [17, 18]. For instance, the cold atoms in optical lattices may provide a means of realizing the nonlinear light-atom coupling appearing in BS-type Hamiltonians, via the engineering of adequate optical lattices [19].
The dynamical nature of the quantum Rabi-like model has been recently discussed. Through time-dependent correlation functions, the dynamical correlation functions of the quantum Rabi model involve the JC scheme in the ultrastrong-coupling and deep strong-coupling regimes [20, 21]. The nonlinear dynamics of trapped-ion models have been proposed for blocking the propagation of quantum information along the Hilbert space of the JC and quantum Rabi models [22]. Exploiting the parity symmetry of the Rabi model, Hu et al [23] recursively evolved the corresponding positive and negative parities using the expansion coefficients in the dynamical equations.
In this paper, quantum spin-boson systems expressed by time-independent Hamiltonians $\hat{H}$ are evolved by a non-perturbation numerical approach. When the stationary eigenvalues of the Hamiltonians are known, the non-perturbation numerical approach recovers the existing results [4, 24] obtained by evolving the known eigenvalues and eigenvectors of $\hat{H}$ [10, 25] from the initial state. However, when the eigenvalues are unknown, our method can evolve the properties from different initial states. Herein, we apply this method to the BS model and discuss the evolution of the mean photon number, the atomic inversion, and the ground-state properties in the absence and presence of RWA. Ng et al proved that the BS Hamiltonian is non-real positive-definite when the strengths of the rotating- and counter-rotating-wave terms are equal [15]. Such a system has no physical meaning and no definite time evolution. However, Rodríguez-Lara et al [26] temporally evolved the equivalent nonlinear JC photonic lattice using quantum optics-based methods via the analogy between the transport of single-photon states and propagation of a classical field. We resolve this obvious contradiction through numerical calculations.
The remainder of this paper is organized as follows. Section 2 introduces the BS model, and section 3 presents our non-perturbation numerical method. Sections 4 and 5 present the BS results with and without RWA, respectively. Conclusions are presented in section 6.
2. The model
The BS model generalizes the JC model to allow nonlinear couplings [15]. The Hamiltonian (ℏ=1) of the BS model is given by$\begin{eqnarray}\begin{array}{rcl}\hat{H} & = & {\omega }_{f}{\hat{a}}^{\dagger }\hat{a}+\displaystyle \frac{{\omega }_{0}}{2}{\hat{\sigma }}_{z}+{g}_{-}(\hat{a}\sqrt{{\hat{a}}^{\dagger }\hat{a}}{\hat{\sigma }}_{+}+\sqrt{{\hat{a}}^{\dagger }\hat{a}}{\hat{a}}^{\dagger }{\hat{\sigma }}_{-})\\ & & +{g}_{+}(\hat{a}\sqrt{{\hat{a}}^{\dagger }\hat{a}}{\hat{\sigma }}_{-}+\sqrt{{\hat{a}}^{\dagger }\hat{a}}{\hat{a}}^{\dagger }{\hat{\sigma }}_{+}),\end{array}\end{eqnarray}$where ωf and ω0 are the frequencies of the field and the two-level transition, respectively; ${\hat{a}}^{\dagger }(\hat{a})$ is the creation (annihilation) operator of the field; and ${\hat{\sigma }}_{i}$ ($i=z,+,-$) are Pauli matrices. The nonlinear coupling is split into a rotating-wave term described by the coupling strength g− and a counter-rotating term described by the coupling strength g+.
When g+=0, we restore the original BS model without the counter-rotating term. In this case, the energy is lower-bound only if ${g}_{-}\leqslant 1/(2S)$, where S denotes the eigenvalues of the constant of motion of the Casimir operator ${\hat{S}}^{2}$. In terms of the Pauli operators, ${\hat{S}}^{2}={\hat{S}}_{\mp }{\hat{S}}_{\pm }+{\hat{S}}_{z}^{2}\pm {\hat{S}}_{z}$. When ${g}_{-}\gt 1/(2S)$, the eigenenergies can be arbitrarily low. As there is no stable ground state, the Hamiltonian is unbounded and no ground state exists. In this sense, the model appears to be incomplete for ${g}_{-}\gt 1/(2S)$. The ground-state energy decreases indefinitely when the eigenvalues of $\hat{C}\,={\hat{a}}^{\dagger }\hat{a}+2{\hat{S}}_{z}$, denoted as c, increase above a critical value. This effect illustrates a dramatic phase transition between the normal phase and the super-radiant phase [27].
When the counter-rotation term is included (${g}_{+}\ne 0$), the Hamiltonian cannot be analytically diagonalized, but the eigenstates and eigenenergies can be computed by a numerical approach. After expanding the initial state in terms of the numerically determined eigenstates and eigenenergies, one can temporally evolve the physical quantities. The following section proposes another formal and systematic numerical method, which does not require precalculation of the stationary eigenstates and eigenenergies.
3. Non-perturbative time evolution
The dynamics of a Hamiltonian system are given by the time-dependent Schrödinger equation$\begin{eqnarray}{\rm{i}}{\hslash }\displaystyle \frac{\partial }{\partial t}| t\rangle =\hat{H}| t\rangle ,\end{eqnarray}$which is formally solved as$\begin{eqnarray}| t\rangle ={{\rm{e}}}^{-{\rm{i}}\hat{H}t}| 0\rangle =\sum _{n=0}^{\infty }\displaystyle \frac{1}{n!}{\left(-{\rm{i}}\hat{H}t\right)}^{n}| 0\rangle .\end{eqnarray}$Here, $| 0\rangle $ denotes the initial state $| t=0\rangle $. Defining$\begin{eqnarray}| {B}_{n}\rangle =\displaystyle \frac{1}{n!}{\left(-{\rm{i}}\hat{H}t\right)}^{n}| 0\rangle ,\end{eqnarray}$we notice that $| {B}_{n}\rangle $ can be concisely iterated as$\begin{eqnarray}| {B}_{n+1}\rangle =-\displaystyle \frac{{\rm{i}}t}{n+1}\hat{H}| {B}_{n}\rangle .\end{eqnarray}$$\hat{H}$ can be easily operated on the initial state $| 0\rangle $ when $| 0\rangle $ is expanded in the Fock state space $\{| p\rangle \}$, p=0, 1, 2, … or the coherent space $\{| \alpha \rangle \}$.
In the following analysis, we illustrate the proposed numerical recipe on a BS model. The time-dependent wavefunction of the two-level system is expanded as an excited state $| e\rangle $ and the ground state $| g\rangle $, namely, as $| t\rangle =| t{\rangle }_{1}| e\rangle +| t{\rangle }_{2}| g\rangle $, where $| t{\rangle }_{1}$ and $| t{\rangle }_{2}$ are the time-dependent wavefunctions in the $| e\rangle $ and $| g\rangle $ states, respectively. Furthermore, $| t\rangle $ is expanded as the following Taylor series:$\begin{eqnarray}\begin{array}{rcl}| t\rangle & = & | t{\rangle }_{1}| e\rangle +| t{\rangle }_{2}| g\rangle \\ & = & \displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{1}{n!}{\left(-{\rm{i}}\hat{H}t\right)}^{n}(| 0{\rangle }_{1}| e\rangle +| 0{\rangle }_{2}| g\rangle ).\end{array}\end{eqnarray}$Defining$\begin{eqnarray}| {B}_{n}\rangle =\displaystyle \frac{1}{n!}{\left(-{\rm{i}}\hat{H}t\right)}^{n}(| 0{\rangle }_{1}| e\rangle +| 0{\rangle }_{2}| g\rangle ),\end{eqnarray}$we recover the iteration relation equation (5). Expanding $| {B}_{0}\rangle $ as $| 0{\rangle }_{1}| e\rangle +| 0{\rangle }_{2}| g\rangle $, equation (7) can be rewritten as the following operator iteration relation:$\begin{eqnarray}| {B}_{n}\rangle =\displaystyle \frac{1}{n!}{\left(-{\rm{i}}\hat{H}t\right)}^{n}| {B}_{0}\rangle .\end{eqnarray}$
Once $| {B}_{n}\rangle $ is known, the time evolution of the state is clearly given by$\begin{eqnarray}| t\rangle =\sum _{n}| {B}_{n}\rangle .\end{eqnarray}$The initial states can then be expanded in the orthogonal complete Fock basis $\{| p\rangle \}$ as$\begin{eqnarray}| {B}_{n}\rangle =\sum _{p=0}^{\infty }({f}_{p}^{1,n}(t)| e\rangle +{f}_{p}^{2,n}(t)| g\rangle )\otimes | p\rangle .\end{eqnarray}$Substituting $| {B}_{0}\rangle ={\sum }_{p=0}^{\infty }({f}_{p}^{1,0}| e\rangle +{f}_{p}^{2,0}| g\rangle )\otimes | p\rangle ,$ into equation (8), applying $\hat{a}| p\rangle =\sqrt{p}| p-1\rangle $ and ${\hat{a}}^{\dagger }| p\rangle =\sqrt{p+1}| p+1\rangle $ (which appear in ${\hat{H}}^{n}$ over the Fock basis $| p\rangle $), and comparing the result with equation (10), we obtain the following iterative relationship between the corresponding expansion coefficients ${f}_{p}^{i,n}(i=1,2;p=0,1,2,\ldots )$,$\begin{eqnarray}\left[\begin{array}{c}{f}_{0}^{1,n}(t)\\ {f}_{1}^{1,n}(t)\\ .\\ .\\ .\\ {f}_{0}^{2,n}(t)\\ {f}_{1}^{2,n}(t)\\ .\\ .\\ .\end{array}\right]\,=\,\displaystyle \frac{{\left(-{\rm{i}}t\right)}^{n}}{n!}{Q}^{n}\left[\begin{array}{c}{f}_{0}^{1,0}(0)\\ {f}_{1}^{1,0}(0)\\ .\\ .\\ .\\ {f}_{0}^{2,0}(0)\\ {f}_{1}^{2,0}(0)\\ .\\ .\\ .\end{array}\right],\end{eqnarray}$where Q is the transfer matrix, given by$\begin{eqnarray}Q=\left[\begin{array}{cc}{Q}_{11} & {Q}_{12}\\ {Q}_{21} & {Q}_{22}\end{array}\right],\quad \end{eqnarray}$and$\begin{eqnarray}\begin{array}{rcl}{Q}_{11} & = & ({\omega }_{0}/2+{\omega }_{f}p){\delta }_{p,{p}^{{\prime} }},\\ {Q}_{22} & = & (-{\omega }_{0}/2+{\omega }_{f}p){\delta }_{p,{p}^{{\prime} }},\\ {Q}_{12} & = & {g}_{+}(p+1){\delta }_{p+1,{p}^{{\prime} }}+{g}_{-}p{\delta }_{p-1,{p}^{{\prime} }},\\ {Q}_{21} & = & {g}_{-}(p+1){\delta }_{p+1,{p}^{{\prime} }}+{g}_{+}p{\delta }_{p-1,{p}^{{\prime} }}.\end{array}\end{eqnarray}$This formulation transforms the operator iteration relations between the $| {B}_{n}\rangle $ into matrix relations between the ${f}_{p}^{i,n}(i=1,2)$. As the Hamiltonian multiplication ${\hat{H}}^{n}$ is absorbed in the matrix multiplication Q as Qn, the number of numerical calculations is greatly reduced, gaining an additional advantage. In this formulation, the matrix multiplication Qn is uniquely determined by the Hamiltonian and is independent of the ${f}_{p}^{i,0}(i=1,2;p=0,1,2,\ldots )$ determined by the initial states. The time-dependent properties of the system under different initial states depend on the same Qn. Therefore, the matrix multiplication Qn can be calculated and stored prior to defining the arbitrary initial states.
To avoid error accumulation during the time evolution, we discrete the time t into K steps, where each unit period Δt satisfies $t=K{\rm{\Delta }}t$. The term$\begin{eqnarray}M({\rm{\Delta }}t)={\left[\displaystyle \frac{{\left(-{\rm{i}}{\rm{\Delta }}t\right)}^{n}}{n!}{Q}^{n}\right]}^{K},K=0,1,2,\ldots \end{eqnarray}$is saved after each calculation. When $M({\rm{\Delta }}t)$ acts on ${f}_{p}^{i,0}(i=1,2;p=0,1,2,\ldots )$ for different K, ${f}_{p}^{i,n}(t)(i=1,2;p=0,1,2,\ldots )$ at each time is obtained as$\begin{eqnarray}\begin{array}{c}\left[\begin{array}{c}{f}_{0}^{1,n}(t)\\ {f}_{1}^{1,n}(t)\\ .\\ .\\ .\\ {f}_{0}^{2,n}(t)\\ {f}_{1}^{2,n}(t)\\ .\\ .\\ .\end{array}\right]\,=\,M({\rm{\Delta }}t)\left[\begin{array}{c}{f}_{0}^{1,0}(0)\\ {f}_{1}^{1,0}(0)\\ .\\ .\\ .\\ {f}_{0}^{2,0}(0)\\ {f}_{1}^{2,0}(0)\\ .\\ .\\ .\end{array}\right].\end{array}\end{eqnarray}$
In summary, once the initial state is given, ${f}_{p}^{i,0}$, $i=1,2;$ $p=0,1,2,\ldots $ can be determined, and the subsequent ${f}_{p}^{i,n}(t),i=1,2;p=0,1,2,\ldots $ can be calculated by equations (11)–(15). Numerically, we must truncate n and p. Different combinations of the truncated forms, denoted as N and P, respectively, were tested in the present study. The results with truncation errors below 10−9 are presented below.
Now, assume a coherent initial state of the system:$\begin{eqnarray}| 0\rangle ={{\rm{e}}}^{-{\alpha }^{2}/2}\sum _{p}\displaystyle \frac{{\alpha }^{p}}{\sqrt{p!}}| p\rangle (\sin \theta | e\rangle +\cos \theta | g\rangle ),\end{eqnarray}$where the amplitude α of the coherent state and the mean photon number of the initial state are related: ${\alpha }^{2}=\langle p\rangle $. Given $| {B}_{0}\rangle =| 0{\rangle }_{1}| e\rangle +| 0{\rangle }_{2}| g\rangle $ and equation (10), ${f}_{p}^{i,0}(0)$ for $i=1,2$ are obtained as$\begin{eqnarray}\begin{array}{rcl}{f}_{p}^{1,0}(0) & = & {{\rm{e}}}^{-{\alpha }^{2}/2}\displaystyle \frac{{\alpha }^{p}}{\sqrt{p!}}\sin \theta ,\\ {f}_{p}^{2,0}(0) & = & {{\rm{e}}}^{-{\alpha }^{2}/2}\displaystyle \frac{{\alpha }^{p}}{\sqrt{p!}}\cos \theta .\end{array}\end{eqnarray}$
The resulting wavefunction $| t\rangle $, given by$\begin{eqnarray}| t\rangle =\sum _{n}\sum _{p}({f}_{p}^{1,n}(t)| p\rangle | e\rangle +{f}_{p}^{2,n}(t)| p\rangle | g\rangle ),\end{eqnarray}$can be numerically calculated by the above procedure. In the following analysis, we investigate two typical measurements of interest in the quantum optics community. The first is the mean photon number $\langle \hat{n}\rangle $, expressed as$\begin{eqnarray}\begin{array}{rcl}\langle \hat{n}\rangle & = & \langle t| {\hat{a}}^{\dagger }\hat{a}| t\rangle \\ & = & \displaystyle \sum _{m}\displaystyle \sum _{q}\left[{f}_{q}^{1,m* }(t)\langle q| \langle e| +{f}_{q}^{2,m* }(t)\langle q| \langle g| \right]{\hat{a}}^{\dagger }\hat{a}\\ & & \times \displaystyle \sum _{n}\displaystyle \sum _{p}[{f}_{p}^{1,n}(t)| p\rangle | e\rangle +{f}_{p}^{2,n}(t)| p\rangle | g\rangle ]\\ & = & \displaystyle \sum _{m,n}\displaystyle \sum _{p}p[{f}_{p}^{1,m* }(t){f}_{p}^{1,n}(t)+{f}_{p}^{2,m* }(t){f}_{p}^{2,n}(t)].\end{array}\end{eqnarray}$
The second is the atomic inversion $\langle {\hat{\sigma }}_{z}\rangle $:$\begin{eqnarray}\begin{array}{rcl}\langle {\hat{\sigma }}_{z}\rangle & = & \langle t| {\hat{\sigma }}_{z}| t\rangle \\ & = & \displaystyle \sum _{m}\displaystyle \sum _{q}[{f}_{q}^{1,m* }(t)\langle q| \langle e| +{f}_{q}^{2,m* }(t)\langle q| \langle g| ]\\ & & \times (| e\rangle \langle e| -| g\rangle \langle g| )\\ & & \times \displaystyle \sum _{n}\displaystyle \sum _{p}[{f}_{p}^{1,n}(t)| p\rangle | e\rangle +{f}_{p}^{2,n}(t)| p\rangle | g\rangle ]\\ & = & \displaystyle \sum _{m,n}\displaystyle \sum _{p}[{f}_{p}^{1,m* }(t){f}_{p}^{1,n}(t)-{f}_{p}^{2,m* }(t){f}_{p}^{2,n}(t)].\end{array}\end{eqnarray}$Physically, the mean photon number and the atomic inversion may periodically collapse and revive, as observed in the JC model.
For comparison, we briefly describe the time evolution from the known eigenvalues and eigenvectors of $\hat{H}$ (hereafter, this method is abbreviated as the TEEE method). The stationary state of equation (1) is assumed as $| \psi \rangle \,={\sum }_{p=0}^{\infty }[{a}_{p}| p\rangle | e\rangle +{b}_{p}| p\rangle | g\rangle $, where $| p\rangle $ is the Fock state in the Fock representation. Solving $\hat{H}| \psi \rangle =E| \psi \rangle $, we get$\begin{eqnarray}\begin{array}{rcl} & & \left({\omega }_{f}p+\displaystyle \frac{{\omega }_{0}}{2}\right){a}_{p}+{g}_{-}(p+1){b}_{p+1}+{g}_{+}{{pb}}_{p-1}={{Ea}}_{p},\\ & & \left({\omega }_{f}p-\displaystyle \frac{{\omega }_{0}}{2}\right){b}_{p}+{g}_{+}(p+1){a}_{p+1}+{g}_{-}{{pa}}_{p-1}={{Eb}}_{p}.\end{array}\end{eqnarray}$After numerically solving the above equations, we obtain the eigenenergy spectrum $\{{E}_{j}\}$ and the corresponding stationary-state set $\{{a}_{p}^{j},{b}_{p}^{j}\}$, where $j=0,1,2,\ldots ,2P+2$ with P being the truncation of the integer number p. Then, expanding the given initial state $| 0\rangle $ in terms of the stationary-state set $| 0\rangle ={\sum }_{j=0}^{2P+2}{F}_{j}| {\psi }_{j}\rangle $, where Fj is the expansion coefficient, the time evolution of the wavefunction is obtained as$\begin{eqnarray}\begin{array}{rcl}| t\rangle & = & {{\rm{e}}}^{-{\rm{i}}\hat{H}t}| 0\rangle \\ & = & \displaystyle \sum _{j=0}^{\infty }{F}_{j}{{\rm{e}}}^{-{\rm{i}}{E}_{j}t}\displaystyle \sum _{p=0}^{\infty }\left[{a}_{p}| p\rangle | e\rangle +{b}_{p}| p\rangle | g\rangle \right]\\ & \approx & \displaystyle \sum _{j\,=\,0}^{2P+2}{F}_{j}{{\rm{e}}}^{-{\rm{i}}{E}_{j}t}\displaystyle \sum _{p\,=\,0}^{P}\left[{a}_{p}| p\rangle | e\rangle +{b}_{p}| p\rangle | g\rangle \right].\end{array}\end{eqnarray}$As in equations (19) and (20), we can calculate the time evolutions of the mean photon number $\langle \hat{n}\rangle =\langle t| {\hat{a}}^{\dagger }\hat{a}| t\rangle $ and the atomic inversion $\langle {\hat{\sigma }}_{z}\rangle =\langle t| {\hat{\sigma }}_{z}| t\rangle $.
Here, we emphasize that both the Taylor expansion in equation (3) and the basis expansion in equation (10) include an infinite sum of terms Therefore, the accuracy of the time evolution can be systematically refined by including more Taylor series and Fock states, respectively. The expansion converges quickly due to the factorial of n in the denominator. This approach, which we call the non-perturbation time-evolution method, offers several advantages over expanding the eigensolutions (i.e. the TEEE method). First, this method is directly applicable when the eigenenergy and eigenstate are unknown. Secondly, it systematically obtains the long-time evolution. Finally, the expansion coefficients M(t) in equation (15) need not be recalculated in the time evolutions of different initial states.
We test our numerical method for the dynamical properties of the Rabi model and BS model compared with the TEEE one. We found both of the methods have almost the same achievements. Furthermore, we find that our method may have advantage over the spin-boson model of multi-modes, which is the case for the Rabi-dimer model [28] discussed in appendix A. The performance of the numerical method on a Hamiltonian with anti-Hermitian terms (which models open quantum systems) was also evaluated in appendix B. In the following, we concentrate on the method to the application of the BS model.
4. BS model with RWA $({g}_{+}=0)$
We now apply the above numerical method to the BS model and compare the results with those of the TEEE method and Rodríguez-Lara et al [26]. We first study the BS model with RWA (i.e. g+=0) at resonance, setting ω0=ωf . In this case, the BS model is exactly solvable. Its stationary state can be expressed as $| \psi \rangle ={a}_{p}| p\rangle | e\rangle +{b}_{p}| p+1\rangle | g\rangle $ due to the conservation of the excitation number ${\hat{P}}_{\mathrm{exc}}={\hat{a}}^{\dagger }\hat{a}+\tfrac{{\hat{\sigma }}_{z}}{2}+\tfrac{1}{2}$. By solving the time-independent Schrödinger equation $\hat{H}| \psi \rangle =E| \psi \rangle $, we can obtain the eigenenergy ${E}_{p}^{\pm }\,=(p+\tfrac{1}{2}){\omega }_{f}\pm {g}_{-}(p+1)$. Because that the state $| 0\rangle | g\rangle $, called the dark state, is also its eigenstate, the corresponding eigenenergy is ${E}_{d}=-\tfrac{{\omega }_{f}}{2}$. Therefore, the ground state energy is Ed for ${g}_{-}\leqslant {\omega }_{f}$ and ${E}_{p}^{-}=({\omega }_{f}-{g}_{-})p+\tfrac{{\omega }_{f}}{2}-{g}_{-}$ otherwise. As a result, the energy is unbounded in the region ${g}_{-}\gt {\omega }_{f}$, due to the fact that the integer number p, is theoretically infinite, indicating a quantum phase transition [29] at ${g}_{-}={\omega }_{f}$.
In the following, we show the time evolutions of the mean photon number $\langle \hat{n}\rangle $ and the atomic inversion $\langle {\hat{\sigma }}_{z}\rangle $ in figure 1, starting from a coherent initial state with α=5 . We choose the following initial state$\begin{eqnarray}\begin{array}{rcl}| t=0\rangle & = & | {\alpha }_{+},g\rangle +| {\alpha }_{-},e\rangle \\ & = & | {\alpha }_{+}\rangle | g\rangle +| {\alpha }_{-}\rangle | e\rangle ,\end{array}\end{eqnarray}$where$\begin{eqnarray}\begin{array}{rcl}| {\alpha }_{\pm }\rangle & = & | \alpha \rangle \pm | -\alpha \rangle \\ & = & {{\rm{e}}}^{-| \alpha | /2}\left[\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{{\alpha }^{n}}{\sqrt{n!}}| n\rangle \pm \displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{{\left(-1\right)}^{n}{\alpha }^{n}}{\sqrt{n!}}| n\rangle \right].\end{array}\end{eqnarray}$Physically, the atomic inversion $\langle {\hat{\sigma }}_{z}\rangle $ is equivalent to the atomic intensity difference between the two levels $\langle e\rangle $ and $\langle g\rangle $. The atomic inversion collapses and revives with the variable photon number. Firstly, we study the dynamics for ${g}_{-}\lt {\omega }_{f}$, where the ground state exists and the dynamics is believed to be stable. Under the truncations N=20 and P=50, our results accorded with those of Rodríguez-Lara et al [26] and the TEEE method (truncated with P=60). Similar agreement was observed in other parameter spaces (data not shown). Note that our method can be extended to very long-time evolution. Our tests confirmed highly accurate time-evolution behavior up to t=600 (in units of ωf), and the evolution time is extendable if necessary.
Figure 1.
New window|Download| PPT slide Figure 1.Evolutions of (a) mean photon number $\langle \hat{n}\rangle $ and (b) atomic inversion $\langle {\hat{\sigma }}_{z}\rangle $ in the BS model with RWA. The initial state is the coherent state $| t=0\rangle =| {\alpha }_{+},g\rangle +| {\alpha }_{-},e\rangle $ with α=5 in the BS model when g+=0 at resonance ωf=ω0=1. The coupling strength is g−=0.1ωf . For comparison, the results of equations (11)–(15) (black line) and the TEEE method (blue points) are also plotted.
In the strong coupling case, the counter-rotating wave term cannot be ignored. So in the next section, we will analyze the RWA-void BS model.
5. BS model without RWA $({g}_{+}\ne 0)$
In this section, we study the effects of the counter-rotating term in the BS model. In this case, the quantum phase transition happens at the critical points ${g}_{-}={g}_{+}=\tfrac{{\omega }_{f}}{2}$. So the BS model has a different phase transition points for the case with and without RWA and the transition positions. The region of the unbounded ground state shrinks due to the RWA.
We now confirm that our numerical method is also applicable to the evolution of the RWA-void BS model, i.e. when ${g}_{+}\ne 0$. In this case, the stationary state of the BS model cannot be solved analytically. As before, the non-perturbation method is implemented by equations (11)–(15). While our method must only recalculate M(t) for the different Qs in equation (13), the TEEE method must numerically recalculate the eigenvalues and eigenvectors of the Hamiltonian. Figure 2 shows the evolutions of the atomic inversion $\langle {\hat{\sigma }}_{z}\rangle $ and mean photon number $\langle \hat{n}\rangle $ in RWA-void BS model, starting from the vacuum initial state $| t=0\rangle =| 0,e\rangle $ with parameters ωf=1, ${\omega }_{0}=3/4{\omega }_{f}$ and ${g}_{-}={g}_{+}=0.4{\omega }_{f}$ (corresponding to $\alpha =0$ and θ=π/2 in equation (17)).
Figure 2.
New window|Download| PPT slide Figure 2.Evolution of (a) mean photon number $\langle \hat{n}\rangle $ and (b) atomic inversion $\langle {\hat{\sigma }}_{z}\rangle $ in the RWA-void BS model. The initial state is a vacuum $| t=0\rangle =| 0,e\rangle $ corresponding to $\alpha =0$ and θ=π/2 in equation (17) (meaning that no photons are excited at t=0). Other parameters are ωf=1 , ω0=3/4ωf and ${g}_{-}={g}_{+}=0.4{\omega }_{f}$. The results of equations (11)–(15) (black line) and the TEEE method (blue points) are plotted for comparison.
Our results (truncated by N=30 and P=50) agree well with those of Rodríguez-Lara et al [30] and the TEEE method (truncated by P=60), demonstrating the applicability of our method to the RWA-void BS model. Note that the collapse and revival phenomena are absent in figure 2, as they were removed by the counter-rotating-wave term.
We subsequently consider the RWA-void BS model under the resonance condition ωf=ω0=1 with parameter values ${g}_{-}={g}_{+}=2{\omega }_{f}$. In this case, the Hamiltonian has no ground state and therefore has no physical meaning [15]. The aim is to confirm the behaviors of the mean atomic inversion $\langle {\hat{\sigma }}_{z}\rangle $ and the mean photon number $\langle \hat{n}\rangle $, and the applicability of the numerical method. Rodríguez-Lara et al [26] calculated the physics of this system via quantum optics methods. They assumed an analogy between the transport of single-photon states and the propagation of a classical field. The numerical propagation was stable on a 2000×2000 photonic lattice. However, the ground-state energy of the Hamiltonian is unbounded in this parameter space. In the naive sense, this result implies an unstable mean atomic inversion and mean photon number. In the following, we first reproduce Rodríguez-Lara et al's results by the TEEE method, then reveal the contradiction between the two arguments.
In practice, equation (22) is truncated by the limit of P instead of being implemented in infinite dimensional Fock space. The numerical procedures of the TEEE method should therefore converge as the truncation number increases. In this case, the final results can be reliably obtained while satisfying the given precision requirements. As shown in panels (a) and (b) of figure 3, the mean photon number and mean atomic inversion results that were truncated by P=2000 matched those obtained by the TEEE method of Rodríguez-Lara et al [26] with the same truncation factor. However, when the truncation limit was increased to 3000 (see figures 3(c) and (d)), the results of our method and TEEE were very different (In order to see clearly the difference, we showed only the results for $t\in [0,40]$ comparing to figures 3(a) and (b)). Hence, whether the truncation requires to be increased or whether the discrepancy naturally arises from the non-physicality of the unbounded ground-state energy due to the incompleteness of the model must be elucidated in further analysis.
Figure 3.
New window|Download| PPT slide Figure 3.Evolutions of (a) mean photon number $\langle \hat{n}\rangle $ and (b) atomic inversion $\langle {\hat{\sigma }}_{z}\rangle $ for P=2000. (c) and (d) represent the same data for (a) and (b), respectively, but for P=2000 and 3000 simultaneously. The BS model begins from a vacuum state $| t=0\rangle =| 0,e\rangle $ with ωf=ω0=1 and g−=g+=2ωf.
To resolve the above doubt, figure 4(a) shows the energy-level differences (${\rm{\Delta }}E={E}_{j}-{E}_{0}$) between the jth ($j=1,2,3,\ldots $) excited state and the ground state in the RWA-void BS model. Figure 4(b) compares the (${\rm{\Delta }}E={E}_{j}-{E}_{0}$ in the models truncated by P=2000 and P=3000. Comparison of the two results clearly shows the invalidity of simply increasing the truncation limit, indicating that the non-convergence of the physical quantities originates from the incomplete Hamiltonian and the unbounded ground-state energy.
Figure 4.
New window|Download| PPT slide Figure 4.Top panel: difference between the jth ($j=1,2,3,\ldots $) excited-state energy and the ground-state energy (${\rm{\Delta }}E={E}_{j}-{E}_{0}$) in the RWA-void BS model with parameters ${\omega }_{f}={\omega }_{0}=1$ and ${g}_{-}={g}_{+}=2{\omega }_{f}$, calculated by the TEEE method truncated by P=2000 (red line) and P=3000 (blue line). Bottom panel: energy-level differences for P=2000 (red line in bottom left) and P=3000 (blue line in bottom right panel), respectively.
To further confirm this judgment, we studied the ground-state energy as a function of truncation. Figure 5(a) shows the results for ${g}_{-}={g}_{+}=0.4{\omega }_{f}$. After an initial rapid decrease, the ground-state energy stabilized at P>6. Moreover, as shown in figure 5(b), the ground-state energy decreased linearly with increasing truncation limit for ${g}_{-}={g}_{+}=2{\omega }_{f}$. In other words, the investigated system lacked a ground state and was physically infeasible.
Figure 5.
New window|Download| PPT slide Figure 5.Ground-state energy versus truncation factor in the BS model with parameters (a) ωf=1 and ${\omega }_{0}=3/4{\omega }_{f}$, ${g}_{-}={g}_{+}=0.4{\omega }_{f}$, and (b) ${\omega }_{f}={\omega }_{0}=1$ and ${g}_{-}=g+=2{\omega }_{f}$.
6. Conclusion
We proposed a systematic numerical method that calculates the time evolution of the physical quantities in the BS model of the interactions between a two-level atomic system and a field. Applying this method, we successfully transformed the Hamiltonian multiplication of ${\hat{H}}^{n}$ into a matrix multiplication of Qn. This transformation greatly reduced the numerical calculation involved. Moreover, the matrix multiplication Qn is independent of the initial states, meaning that the time-dependent properties of the system can be computed for arbitrary initial states without recalculating the matrix multiplication Qn.
The performance of the method was tested in numerical comparisons with the TEEE method. The method performed well in the whole coupling regime, and captured the collapses and revivals of the atomic inversion along with the photon-number variations. In the parameter space of a non-real positive-definite Hamiltonian, we also resolved the uncertainty regarding whether the evolution was stabilized by increasing the truncation limit over the dimension of the Fock space or whether it was a natural consequence of the unbounded ground-state energy (a non-physical situation caused by the incompleteness of the Hamiltonian). We found that when the ground-state energy of the Hamiltonian is unbounded in a certain parameter space, the time evolution of the physical quantities is unstable for the case without RWA.
Acknowledgments
This work was supported by the National Natural Science Foundation of China under Grant Nos. 11835011 and 11774316.
JaynesE TCummingsF W1963 Comparison of quantum and semiclassical radiation theories with application to the beam maser Proc. IEEE 51 89 DOI:10.1109/PROC.1963.1664 [Cited within: 1]
Gea-BanaclocheJ1991 Atom- and field-state evolution in the Jaynes–Cummings model for large initial fields Phys. Rev. A 44 5913 DOI:10.1103/PhysRevA.44.5913 [Cited within: 1]
ShoreB W2007 Sir Peter Knight and the Jaynes–Cummings model J. Mod. Opt. 54 2009 DOI:10.1080/09500340701303501
HeSWangCChenQ HRenX ZLiuTWangK L2012 First-order corrections to the rotating-wave approximation in the Jaynes–Cummings model Phys. Rev. A 86033837 DOI:10.1103/PhysRevA.86.033837 [Cited within: 1]
MirzaeeMBatavaniM2015 Atom-field entanglement in the Jaynes–Cummings model without rotating wave approximation Chin. Phys. B 24040306 DOI:10.1088/1674-1056/24/4/040306
KunduA2004 Quantum integrability and Bethe ansatz solution for interacting matter-radiation systems J. Phys. A: Math. Gen. 37 L281 DOI:10.1088/0305-4470/37/25/L03 [Cited within: 1]
NgK MLoC FLiuK L2000 Exact eigenstates of the intensity-dependent Jaynes–Cummings model with the counter-rotating term Physica A 275 463 DOI:10.1016/S0378-4371(99)00401-X [Cited within: 4]
YuL BTongN HXueZ YWangZ DZhuS L2012 Simulation of the spin-boson model with superconducting phase qubit coupled to a transmission line Sci. China A 55 1557 DOI:10.1007/s11433-012-4863-x [Cited within: 1]
FengY LZhangKFanJ TMeiFChenGJiaS T2018 Superfluid-superradiant mixed phase of the interacting degenerate Fermi gas in an optical cavity Sci. China A 61123011 DOI:10.1007/s11433-018-9271-5 [Cited within: 1]
WolfF AValloneFRomeroGKollarMSolanoEBraakD2013 Dynamical correlation functions and the quantum Rabi model Phys. Rev. A 87023835 DOI:10.1103/PhysRevA.87.023835 [Cited within: 1]
LiuTWangK LFengM2009 The generalized analytical approximation to the solution of the single-mode spin-boson model without rotating-wave approximation Europhys. Lett. 86 54003 DOI:10.1209/0295-5075/86/54003 [Cited within: 1]
Rodríguez-LaraB MSoto-EguibarFCárdenasA ZMoya-CessaH M2013 A classical simulation of nonlinear Jaynes–Cummings and Rabi models in photonic lattices Opt. Express 21 12888 DOI:10.1364/OE.21.012888 [Cited within: 5]
CordeiroFProvidênciaCda ProvidênciaJNishiyamaS2008 The two-photon Jaynes–Cummings model: a coherent state approach Adv. Stud. Theor. Phys. 2 181 DOI:10.1088/1751-8113/40/40/010 [Cited within: 1]
CordeiroFProvidênciaCda ProvidênciaJNishiyamaS2007 The Buck–Sukumar model described in terms of ${su}(2)\bigotimes {su}(1,1)$ coherent states J. Phys. A: Math. Theor. 40 12153 DOI:10.1088/1751-8113/40/40/010 [Cited within: 1]
ChenQ HLiuTZhangY YWangK L2010 Quantum phase transition in coupled two-level atoms in a single-mode cavity Phys. Rev. A 82053841 DOI:10.1103/PhysRevA.82.053841 [Cited within: 1]
Rodríguez-LaraB M2014 Intensity-dependent quantum Rabi model: spectrum, supersymmetric partner, and optical simulation J. Opt. Soc. Am. B 31 1719 DOI:10.1364/JOSAB.31.001719 [Cited within: 1]
SergiAZloshchastievK G2013 Non-Hermitian quantum dynamics of a two-level system and models of dissipative environments Int. J. Mod. Phys. B 271350163 DOI:10.1142/S0217979213501634 [Cited within: 1]
MoiseyevN2011Non-Hermitian Quantum MechanicsCambridgeCambridge University Press [Cited within: 1]
Echeverri-ArteagaaSVinck-PosadaaHGómezE A2019 A comparative study on different non-Hermitian approaches for modeling open quantum systems Optik 180 505 DOI:10.1016/j.ijleo.2018.11.133 [Cited within: 1]
Echeverri-ArteagaSVinck-PosadaHGómezE A2018 A comparative study on the reliability of non-Hermitian effective Hamiltonian approach for modeling open quantum systems Optik 171 413 DOI:10.1016/j.ijleo.2018.11.133 [Cited within: 1]
HwangM JKimM SChoiM S2016 Recurrent delocalization and quasi-equilibration of photons in coupled circuit QED systems Phys. Rev. Lett. 116153601 DOI:10.1103/PhysRevLett.116.153601 [Cited within: 1]
,1,4, Ke-Lin Wang31Department of Physics, 
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
