A continued fraction based approach for the Two-photon Quantum Rabi Model

We study the Two Photon Quantum Rabi Model by way of its spectral functions and survival probabilities. This approach allows numerical precision with large truncation numbers, and thus exploration of the spectral collapse. We provide independent checks and calibration of the numerical results by studying an exactly solvable case and comparing the essential qualitative structure of the spectral functions. We stress that the large time limit of the survival probability provides us with an indicator of spectral collapse, and propose a technique for the detection of this signal in the current and upcoming quantum simulations of the model.

A continued fraction based approach for the two-photon Quantum Rabi Model elena Lupo 1,2 , Anna Napoli 1,3 , Antonino Messina 3,4 , Enrique solano 5,6,7 & Íñigo L. Egusquiza 8 We study the two photon Quantum Rabi Model by way of its spectral functions and survival probabilities. This approach allows numerical precision with large truncation numbers, and thus exploration of the spectral collapse. We provide independent checks and calibration of the numerical results by studying an exactly solvable case and comparing the essential qualitative structure of the spectral functions. We stress that the large time limit of the survival probability provides us with an indicator of spectral collapse, and propose a technique for the detection of this signal in the current and upcoming quantum simulations of the model.
The Quantum Rabi Model (QRM) and the Two-Photon Quantum Rabi Model (2γQRM) represent two basic models for the description of the interaction of light and matter. The first one describes a two-level system bilinearly coupled to a quantized bosonic field mode; 2γQRM is one of its simplest generalizations, in which the interaction term is now quadratic in the annihilation and creation bosonic operators. The bilinear QRM for light-matter interaction appeared more than 80 years ago [1][2][3] . Yet interest in this model has never waned and, rather, it has even grown recently. This growth mainly stems from its potential application to platforms used for quantum technologies 4 . The QRM depends on two independent parameters, and the dynamical properties of the atom mode are qualitatively very different in different regions of the parameter space. Most of the experimental Cavity Quantum Electrodynamics (CQED) setups are characterized by physical conditions inside the weak-coupling regime, in which the Quantum Rabi model can be effectively simplified to the exactly treatable Jaynes-Cummings model. So as to best describe new, more advanced quantum devices, such as superconducting circuits or trapped ions systems, for instance, the description of the QRM must be extended to the appropriate regions of the parameter space, for which the Jaynes-Cummings approximation fails.
Alternatively, one can view these newer platforms as 'Quantum Simulators' [5][6][7] , in which one can realize models that had been previously discarded as 'unphysical' . In fact, coupling constant values much higher than the ones typical of CQED setups have been measured in the last years, even reaching the so-called Ultrastrong Coupling (USC, ω ω .   g 0 1 ) and the Deep Strong Coupling (DSC, ω  g ) regimes in the context of circuit Quantum Electrodynamics cQED [8][9][10] .
In this vein of Quantum Simulation, other possibilities have appeared. For instance, the interaction Hamiltonian for trapped ions is non-linear, thus allowing this system to be exploited in order to investigate the dynamics of various QRM generalizations 7,11 . For these reasons, the interest in the QRM and its variants has been rekindled, and a strong effort to construct their solutions and to clarify the relative dynamical properties is under way [12][13][14][15][16][17][18] . A major role in these new developments has been played by the analytic solutions of the QRM, found first in 2011 15 and based on its representation in the Bargmann space of the holomorphic functions 19 . Other approaches, exploiting a suitable Bogoliubov transformation 17 , or an expansion in the basis of Heun functions, have also been proposed 20,21 .
Among the many generalizations of the QRM, the Two-photon Quantum Rabi Model (2γQRM) is of particular interest. It was introduced as an effective model for a three-level system interacting with a bosonic mode in which the intermediate level can be adiabatically eliminated [22][23][24][25] . Even if the original phenomenological model was treated in the Rotating Wave Approximation 22 , some work on the influence of the counter-rotating terms has been carried out in the past [26][27][28][29] . The more recent possibility of realizing the 2γQRM in quantum simulators has sparked a new flurry of studies. In particular one should notice the recent proposal for its implementation in trapped ions systems and superconducting circuits 11,[30][31][32] . Moreover, after Braak's solution for the QRM, the same approach was applied to the determination of the 2γQRM spectrum using G-functions 33,34 . Alternatively the search of the analytic solution has also been presented as an expansion in the generalized squeezed number states 17,35,36 . An important feature of the model, namely the collapse of the discrete spectrum into a continuum at a value of the coupling constant g = ω/2 37 , has also been made evident both with squeezed states 38 and with the Bargmann space 19 approach.
Nonetheless, useful as these analytical approaches are for the spectrum as a function of the coupling strength between the fermionic and bosonic subsystems, an analytical form of the eigenstates, and thus of all quantities of interest, is still to be obtained 18 . One such quantity, of particular relevance from an experimental point of view, is the spectral function, defined as ρ ) , in terms of a generic state |Ψ〉 of the system. In fact ρ(E, |Ψ〉) contains all the information useful for generating the time evolution of |Ψ〉, namely those eigenvalues of the Hamiltonian whose eigenfunctions have an overlap with |Ψ〉 and the relative transition probabilities.
In this paper we put forward the spectral analysis of a factorized state σ σ | 〉 ≡ | 〉| 〉 n n , of the 2γQRM, n being the eigenvalue of the number operator † a a and σ being the eigenvalue of the spin operator σ z . This approach, valid in each point of the parameter space, is an alternative to the Bargman solution of the model. To achieve this goal the relevant matrix element of the resolvent is presented in continued fraction form. We have thus direct access to two complementary quantities of interest: the spectral density and the survival probability. The structure of our approach allows us clean access to the dynamics of the system near the collapse point of the 2γQRM corresponding to the value of g = 0.5ω.
The paper is organized as follows: we present the model and apply a unitary transformation such that the eigenstates factorize in a bosonic and a spin part 39 ; then we exploit the connection between the resolvent of a tridiagonal matrix and continued fractions to obtain a numerical determination for the spectral function of factorized states σ | 〉 n, ; finally we use the previous results to study the survival probability of the vacuum state of the system.

the two-photon model
The Two-photon Quantum Rabi Model (2γQRM) presents an interaction which is non linear in the bosonic operators. The Hamiltonian can be expressed as: where ω is the frequency of the bosonic mode, ω 0 the atomic frequency and g the coupling constant between the two subsystems. Here and subsequently we set ℏ to 1. The spectrum of this model has been numerically calculated in many works 11,17,27,[33][34][35][36][37] , and a link with the squeezed number states has been pointed out 17,35,36 . This link provides us with a better understanding of the spectrum collapse at g = ω/2 11,37 , as follows. Define, as usual, the squeezing , and consider the 2γQRM for ω 0 = 0, written in the basis for which σ x is diagonal. Under squeezing transformations with squeezing parameters β = ± ω ± − ( ) tanh g 1 2 1 2 one of the diagonal elements of the Hamiltonian becomes a harmonic oscillator. Clearly. the limit |g|→ω/2 is the limit of infinite squeezing. This entails, in what regards the spectrum, the collapse of eigenvalues into a continuum in the limit g → 0.5ω, and the (generalized) eigenstates are no longer normalizable 11,37 . As in 11 , one can rewrite the Hamiltonian (1) in terms of the position and momentum operators of the oscillator, = + ω † x a a ( ) , with unit mass: For g < ω/2 the effective potential makes the system stable, while at the point g = ω/2 one of the two quantities x 2 or p 2 disappears and the spectrum collapses into a continuum. This is immediately obvious if ω 0 = 0. Were this parameter different from zero, isolated eigenstates would appear. In the context of an analysis of the asymptotic behaviour of solutions in Bargmann space, the collapse point coincides with the limit situation, for which the eigenfunction is no longer normalizable 11,16,37 .
As is well known, the QRM Hamiltonian commutes with a parity operator, and its eigenvalues can be arranged in parity subspaces. In the case of the 2γQRM the symmetry is  4 , since the Hamiltonian commutes with Π σ = − π † e i a a z 4 2 , whose eigenvalues are the quartic roots of unity {±1, ±i}. It follows that the full Hilbert space is organized in four infinite-dimensional chains: www.nature.com/scientificreports www.nature.com/scientificreports/ We denote the corresponding four infinite-dimensional subspaces S w , with w ∈ {±1, ±i}. For instance, the vacuum state | − 〉 ≡ | 〉| − 〉 0, 0 1 belongs to the subspace S +1 . Explicitly, We shall now apply a transformation which factorizes the state |Ψ w 〉 into a bosonic and an atomic part, following the procedure of 39 for the QRM. In other words 12 , we use the parity basis. This factorization is indeed achieved with the rotation

This rotation transfoms the Hamiltonian into
The coupling term is now expressed as diagonal in the bosonic number operator. Under this rotation the subspaces of constant 4-parity become: i i 1 1 and the Hamiltonian projected into each subspace is a quadratic of the bosonic creation and annihilation operators,

the spectral Function and the Resolvent
In this section we derive an expression of the spectral function of a factorised state |n, σ〉 in a continued fraction form. The spectral function ρ |Ψ〉 |Ψ′〉 E ( , , ) is the matrix element ρ(E) Ψ,Ψ′ of the microcanonical density operator defined in the following way: where H is the Hamiltonian of the model and |ε λ 〉 is the eigenstate of H related to the eigenvalue Other than in quantum statistical mechanics, it appears in relation with the resolvent (E − H) −1 , whose spectral representation is: H 0 From the definition (9) one can see that the diagonal element ρ(E, |Ψ〉) contains all the spectral information useful in the study of the state |Ψ〉 of the system. It can be in fact interpreted as the probability distribution of the state |Ψ〉 to be in a particular eigenstate of the Hamiltonian: 2 and it is instrumental in studying the time evolution of the state. Its numerical computation can be rather involved if attacked in terms of Bargmann functions, though. Here we address this issue by making use of the connection of the spectral function to the resolvent of the system. The , with ε > 0 and P principal part, determine www.nature.com/scientificreports www.nature.com/scientificreports/ In the factorized states basis σ | 〉 n, , the resolvent of the QRM and 2γQRM can readily be expressed in continued fraction form (see 40,41 or Appendix A), which makes a numerical calculation of the spectral function ρ σ | 〉 E n ( , , ) accessible. Notice that the use of continued fractions has been a staple in the treatment of the QRM, in different guises and forms 42 . Taking the rotated Hamiltonian (6), the element of the resolvent related to the state σ | 〉 n, is in the form: where we set z = E − iε for the numerical calculation of (12) and the coefficients A j and R j depend on the subspace which σ | 〉 n, belongs to: if the state has the form ± − n 2 , ( 1) n (i.e. it belongs to the subspace  S 1 ) the coefficients are A j = 2jω ± (−1) j ω 0 /2 and = − R g j j 2 (2 1) j ; if the state is in the form | + ± − 〉 n 2 1, ( 1) n (i.e. it belongs to the subspace S ±i ) the coefficients are A j = (2j + 1)ω  (−1) j ω 0 /2 and = + R g j j 2 (2 1) j . Equation (13) allows us to calculate the spectral function of any factorized state σ | 〉 n, of the 2γQRM. In this work we show the results related to the positive parity subspace S +1 . Since we are working in the rotated basis, from now on we use |2n, −〉 as notation for the state belonging to +  S 1 . The convergence of the continued fraction has been determined through Pringsheim's Theorem, under the condition that g < ω/2 (see Appendix B). The actual computation of the continued fraction expansion involves a truncation in Fock space for each truncation of the continued fraction. In Fig. 1 we report the numerical determination of the the spectral density for the vacuum state of the 2γQRM at different values of g/ω. Notice that the parameter ε has to be fixed for the numerical evaluation. Its value is chosen in such a way it doesn't affect the ratio between the peaks, and a smaller value would result only in a common scaling factor that does not bring improvement in the determination of the spectral function.
The method at hand, namely the numerical computation by continued fractions of spectral functions, allows us to insert much higher truncation numbers than with a direct simulation with truncation in Fock space, even very close to the collapse point g/ω = 0.5, where the spectrum will no longer be purely discrete. Fig. 2 shows the spectral density as we approach the special value g/ω = 0.5, making apparent this change of the spectrum into an isolated discrete value and a continuum.
We now apply our technique to the collapse point g/ω = 0.5, even though Pringsheim's theorem only guarantees convergence in the discrete case g/ω < 0.5. In fact the continued fraction approach allows only a discrete approximation of a continuum spectrum, but this is done at very high truncation numbers. In Fig. 3 the spectral function of the vacuum state for g/ω = 0.5 is calculated at different values of ω 0 . A first point of note is that the presence of an isolated ground state is linked to the atomic frequency ω 0 being different from zero. Secondly, observe that the energy difference between the ground state and the continuum (Fig. 4) is not linear in ω 0 , as observed also for g/ω < 0.5 in previous papers 33,34,37 .
In the case ω 0 = 0 the spectral function ρ | − 〉 E n ( , 2 , ) 0 , with | − 〉 ∈ +  n S 2 , 1 , can be calculated analytically. Consider the Hamiltonian of the 2γQRM projected in +  S 1 for a coupling value g = ω/2:     www.nature.com/scientificreports www.nature.com/scientificreports/ x p 2 2 0 ( 1 )/4 2 2 In the case ω 0 = 0 the Schrödinger equation takes the form ω ω − |Ψ 〉 = |Ψ 〉 x x E x ( / 2) ( ) ( ) 2 2 and the eigenstates coincide with the position operator eigenstates |x〉. Therefore, the spectral function related to the state |2n, −〉 can be expressed in terms of the Hermite polynomials H m (ξ): We can now contrast and calibrate the numerical results at ω 0 ≠ 0 for the first six states of the subspace + S 1 with the corresponding analytical expression (16), in Fig. 5. Clearly the qualitative structure is well tracked by our numerical procedure, setting aside the divergence of ρ 0 at E = −ω/2. In particular, notice the number of nodes in the corresponding spectral functions. Moreover, in Fig. 6 we plot the ratio between the two quantities. Even if the polynomial trend of the truncated continued fraction can not track the exponential trend of (16), in a range of high energies for which the Hermite trend contributes mostly, we can notice a constant value which is due to the atomic term in the Hamiltonian (8) becoming progressively less relevant.

the survival probability of the vacuum state
The results of the previous section can be exploited for the determination of an important dynamical quantity: the survival probability, that is, the probability of finding the system in its initial state after a time evolution of interval t.
The connection between the spectral function ρ(E, |Ψ〉) and the survival probability is given through the survival amplitude by Fourier transform, Figure 5. Comparison between the the spectral functions related to the first six states belonging to +  S 1 , between the cases ω 0 = 0.8ω and ω 0 = 0 (whose analytic form is known) at the collapse point g/ω = 0.5. In all cases ε = 0.0005, while the truncation number exploited for the continued fraction is N = 24000. Figure 6. Plot of the ratio between the two cases compared in Fig. 5. It can be seen that for high energies the ratio between the spectral function in the case ω 0 ≠ 0 and the exact case ω 0 = 0 is constant. That is, with the integration on the domain defined by ρ(E, |Ψ〉).
Let us now focus on the vacuum state of the 2γQRM. It is of interest since it can be prepared as the ground state in the decoupled or strong coupling regime ( ω .  g/ 01), and then adiabatically moved to larger couplings. In terms of the eigenenergies E λ and the transition probabilities ε 〈 | − 〉 λ 0, 2 , which can be derived from its spectral function, the survival probability of the vacuum state |0, −〉 can be written as: This connection provides us with a numerical technique to compute the survival probability, through numerical computation of the spectral function. The fact that we do not use matrix inversion, diagonalization, nor exponentiation in the process means that the point of the truncation can be much higher than what could be reasonably achieved with Fock space expansions for the survival probability. This numerical advantage allows us, in particular, an analysis of the survival probability for the 2γQRM close to the collapse point g = ω/2.
In Fig. 7 we report the numerical determination of the survival probability at different values of g/ω. Near the collapse point g = ω/2 interference effects become predominant. This was to be expected from the spectral density depicted in Fig. 2, since the density of eigenstates means that small frequencies (small energy differences) will play a major role in the survival probability. Indeed the long time behaviour of the survival probability becomes flatter, as seen in the last graph of Fig. 7.
We also compute the survival probability for g/ω = 0.5 at different values of ω 0 , as portrayed in Fig. 8. We again see that the survival probability for |0, −〉 presents a dominant constant value, dependent on the atomic parameter, after a short transient. This can be understood by looking at the form of the Survival Probability − P t ( ) 0, in terms of the spectral function, iEt 0, 0 2 /2 2 and application of the Riemann-Lebesgue lemma. Indeed, we know that ρ | − 〉 E ( , 0, ) is integrable -in fact, as pointed out above, it is normalized to 1. Therefore the Fourier transform above tends to zero at infinity. To be more precise, only the discrete part of the spectrum contributes to the long time behaviour,  Notice the asymptotic 1/t behaviour, that is due to the 1/(E + ω/2) 1/2 divergence in the integrand.
As ω 0 grows, a discrete point will appear in the spectrum, and thus a constant term in the long time behaviour of the survival probability. The subleading term will be generically of form 1/t, since the leading behaviour of the Fourier transform of the continuum part will be 1/t or faster decay.

Conclusions and perspectives
In this work we have studied numerically spectral functions for the Two Photon Quantum Rabi Model (2γQRM) and the corresponding survival probabilities. These two quantities are more readily amenable to numerical treatment than direct diagonalization of the Hamiltonian, as is shown by the much higher truncation numbers we can achieve in this approach.
Since there are indeed several proposals for quantum simulation implementation of the 2γQRM 11,30-32 , our improved numerical approach will prove beneficial for their analysis.
This improvement of numerics has allowed us to investigate further the collapse point, at which the spectrum becomes continuous. This is indeed the result recovered both from spectral functions and from survival probabilities.
As all numerics are suspect in the environment of a drastic structural change, such as the spectral collapse at hand, we have proposed an independent check by comparing spectral functions at the collapse point for the exactly solvable case with ω 0 = 0, expressed in terms of Hermite polynomials, with those corresponding to ω 0 ≠ 0. The qualitative structure, in particular the number of modes and the large energy/short time behaviours, is maintained as expected, thus providing us with a calibration tool.
In particular we note that a signature of the collapse of the spectrum into a purely continuous one would be that all survival probabilities necessarily tend to zero. In the case at hand there is a remaining relevant discrete point in the spectrum, and the long time limit of the survival probability is a constant, determined by the projection of the initial state onto the corresponding proper eigenstate.
The direct measurement of such a phenomenon in the survival probability might not be immediately possible in the different platforms in which the 2γQRM is a good description of the dynamics for some range of the parameters. However, there are alternatives to detect the spectral collapse, one of which we now put forward. Consider thus that there is another eigenstate of the full system which can be coupled to the discrete element of the +  S 1 subspace. In such a situation, the long term behaviour of the survival probability for any state in the +  S 1 subspace will be given by coherent Rabi oscillations, providing us with a target for detection.
Notice a very recent alternative proposal to investigate the spectral collapse, in this case for the 2γQRM with full quadratic coupling, studying a two-time correlation for the output field in a driven system 43 .
In summary, we have investigated further the rich phenomenology of the 2γQRM, with emphasis on the numerically computable spectral functions and survival probability, and we suggest new avenues for the exploration of the spectral collapse.

Appendix A: the continued fraction form of the resolvent
In each subspace of defined four-parity Π 4 the rotated Hamiltonian ∼ H is tridiagonal in the basis of Fock states. For instance, the 2γQRM Hamiltonian projected in the subspace of positive parity +  S 1 is (see Eq. (8)): is the same as (25), the derivation of the diagonal element of the resolvent does not change from the one in the subspace From the theory of linear algebra the inverse of a square matrix Q is the matrix of elements Q ij = Δ ij /det(Q), where Δ ij = (−1) i+j det(M ij ) is the (i, j)-cofactor and M ij is the first minor, obtained by eliminating the i-th row and the j-th column. We use the following notation: D 0 as the determinant of − ∼ + z H ( ) 1 , D k as the corresponding determinant of the matrix resulting from eliminating the first k rows and k columns, and ∼ D k as the determinant of the matrix given by the restriction to the first k + 1 rows and columns of the same matrix. So we have: Since the sub-matrix obtained by removing the n-th row and n-th column of − ∼ + z H ( ) 1 has a block-tridiagonal form, applying the formula for the inverse matrix we have: Applying iteratively the two formulas (30) the two quantities D n /D n−1 and ∼ ∼ + D D / n n 1 can be expressed in a continued fraction and a finite continued fraction form respectively: Substituting them in equation (30) we obtain the continued fraction form of a diagonal element of the resolvent (13). whence ω δ < + g 2(2 ) , 2 2 thus inside the normal region before the collapse. The asymptotic convergence does not guarantee convergence of the resolvent if the energy is one of the real eigenvalues; it is enough for our purposes, though, since in order to determine the relevant spectral function we have to compute the limit of the imaginary part of the resolvent.