Spectral determinant of the two-photon quantum Rabi model

The various generalized spectral determinants (G-functions) of the two-photon quantum Rabi model are analyzed with emphasis on the qualitative aspects of the regular spectrum. Whereas all of them yield at least a subset of the exact regular eigenvalues, only the G-function proposed by Chen et al. in 2012 exhibits an explicitly known pole structure which dictates the approach to the collapse point. We derive this function rigorously employing the $\mathbb{Z}_4$-symmetry of the model and show that its zeros correspond to the complete regular spectrum.


Introduction
The quantum Rabi model (QRM) [1] is a well-known minimalistic way to describe the interaction of matter (fermions, discrete degrees of freedom) with light (bosons, continuous degrees of freedom).It consists of a spin 1/2 coupled linearly to a single mode of the radiation field, H R = ωa † a + ∆σ z + g(a + a † )σ x . (1) The Pauli matrices σ x,z describe the spin 1/2 and a (a † ) are the annihilation (creation) operators of the bosonic mode.The Hilbert space is thus H = L 2 (R)⊗C 2 .This model is integrable due to its manifest Z 2 -symmetry for all values of frequency ω and coupling g [2].The small and large coupling regions are also accessible via perturbation theory.
For small g the rotating wave approximation is feasible [1], for large g the adiabatic approximation [3].The validity of the adiabatic approximation for large coupling is due to the boundedness of 1 1 ⊗ σ z and the fact that the operator a † a + g(a + a † ) + g 2 is unitarily equivalent to a † a, the Hamiltonian of the uncoupled field mode.Therefore, the qualitative features of the spectral graph as function of g do not change much in going from small to large coupling [4], as long as neither ω nor ∆ become singular [5,6].Especially the average distance between adjacent levels with the same parity (the eigenvalue of the operator PR = − exp(iπa † a) ⊗ σ z generating the Z 2 -symmetry) is always ω, independent of g and ∆ [7].
There are many generalizations of the QRM which describe various implementations within cavity and circuit QED as well as quantum simulation platforms, e.g. the anisotropic [8] and the asymmetric QRM [9], the Rabi-Stark model [10,11] and the QRM with non-linear coupling [12], the so-called two-photon quantum Rabi model (2pQRM) which is the subject of the present investigation.The 2pQRM has received lot of attention recently due to the possibilities to realize it in superconducting circuits [13,14] or via quantum simulation [15].After a trivial unitary transformation exchanging σ x and σ z , the Hamiltonian of the 2pQRM reads H 2p = ωa † a + (a 2 + a †2 )σ z + ∆σ x . ( For convenience, we measure ω and ∆ in units of g which is set to 1.In contrast to (1), the coupling term in (2) is not relatively bounded with respect to a † a [16] and the spectrum undergoes a dramatic change when ω approaches the critical value ω c = 2 from above.For ω less than 2, H 2p is no longer bounded from below [12].Exactly at ω = 2, the spectrum contains a discrete and a continuous part while it is pure point for ω > 2. The change from discrete to continuous spectrum at ω c has been termed "spectral collapse" [15], although this is a misnomer because the discrete spectrum does not collapse to a single point at ω c .However, exactly this behavior is seen in all numerical evaluations of the spectrum, simply because the Hilbert space used in these calculations has finite dimension.The continuous part of the spectrum at ω c covers the interval [E 0 , ∞[ with E 0 = −ω c /2 = −1 which follows from a direct treatment of this case [17,18].The critical point ω c = 2 allows for a mapping of the eigenvalue problem to a 1D Schrödinger equation in an energy-dependent potential.Similar simplifications at criticality appear also in the Rabi-Stark model [10,19].On the other hand, it is not clear from this calculation whether the Hamiltonian exactly at ω c is the limit of H 2p (ω) for ω → ω + c or a singular special case.To answer this question one has to compute the full spectrum for all ω > ω c and deduce its qualitative behavior in approaching the "collapse" point.This can be done without much effort for ∆ = 0, where the collapse phenomenon already happens.In this case H 2p is equivalent to a linear combination of generators of su (1,1) and can be diagonalized by standard methods [12,20].The general situation is much more complicated.One reason is the fact that the spectral problem for the 2pQRM and the linear QRM cannot be solved by polynomial ansatz functions which yield f.e. the discrete spectrum of the harmonic oscillator and the hydrogen atom.The simple harmonic oscillator with Hamiltonian has a spectrum of pure point type.The eigenfunctions can be obtained with the ansatz ψ n (x) = P n (x)e −ωx 2 /2 , where P n (x) is a polynomial of degree n, n ∈ IN 0 .Functions of this type are clearly normalizable in L 2 (R), so each ψ n (x) is an eigenfunction of H osc with eigenvalue E n = nω.But it is not easy to prove that the discrete spectrum of H osc consists only of the E n , meaning that the set {ψ n } is complete in L 2 (R) (supposing in addition the absence of a continuous spectrum).Moreover, the fact that a polynomial ansatz for the eigenfunctions indeed works is due to the simplicity of the differential equation (H osc − E)ψ = 0 which belongs to the hypergeometric class [21].
In most cases, the eigenfunctions cannot be obtained by simple ansatz functions and the spectral problem becomes difficult because the normalizability condition in for the eigenunction ψ(x) requires knowledge of ψ(x) for all x, it is non-local.
Fortunately, there exists a Hilbert space B, isomorphic to L 2 (R), which consists of analytic functions f (z) of the complex variable z [22].In B, the multiplication operator z is adjoint to the derivative d/dz.All elements f (z) in B satisfy two conditions: The first Bargmann condition means they are holomorphic in all of C, which is a local property of f (z).The second Bargmann condition requires finiteness of f , the Bargmann norm of f (z), obtained form the scalar product, This condition is again non-local.But in the case of the harmonic oscillator (and the QRM) the first, local Bargmann condition is sufficient to determine the whole discrete spectrum and prove its completeness.Writing H osc in terms of a and a † , we have in B, because a = d dz and a † = z.All solutions of (H osc − E)ψ(z) = 0 read ψ E (z) = z E/ω .They have finite Bargmann norm for all real E ≥ 0. So, the second Bargmann condition does not determine the discrete spectrum of H osc .But the first Bargmann condition requires E/ω to be a non-negative integer, yielding the spectrum of the harmonic oscillator.Moreover, the set {z n } is clearly a basis for the entire functions in C, so it is complete in B. The absence of a continuous spectrum follows at once.Note that the ψ En (z) are not obtained by a special ansatz because the Schrödinger equation has no other solutions than z E/ω .We show in the following that also the second Bargmann condition, despite its non-locality, can be used to determine the exact spectrum of Hamiltonians with a single continuous degree of freedom, in this case H 2p , by analyzing the singularity structure of the corresponding ordinary differential equation in the complex domain.This is done for all coupling regimes in section 2. In section 3.1, we discuss generalized spectral determinants proposed in [23] and [24] for the case ω > 2. In section 3.2 we derive the G-function proposed in [25,26].Conclusions are given in section 4.

Singularities of the eigenvalue equation
We write an eigenfunction of H 2p , |ψ ∈ B ⊗ C 2 as the column vector |ψ = (φ 1 (z), φ 2 (z)) T .The eigenvalue equation H 2p |ψ = E|ψ is the following coupled system 2 − ωzφ where φ (n) denotes the n-th derivative of φ(z) with respect to z.The system (7,8) has no singular points for |z| < ∞ but a strong irregular singularity at z = ∞.Upon elimination of φ 2 , the equation for φ 1 is of fourth order, The analysis of the irregular singularity at z = ∞ according to [27] shows that the singularity is unramified with rank two, or, following the notation of [21], "s-rank" three.The class of the singularity is four and therefore maximal [27].This means that for z → ∞ all solutions are entire and have the asymptotic form where c 0 = 0 and γ is a solution of the equation The four solutions of (11), called characteristic exponents of second kind with order two, determine the growth type of the entire function φ 1 (z) which also has order two [28].The coefficient β of the linear term in the exponential factor of ( 10) is an exponent of second kind with order one [21].In general, the formal expansion (10) of a solution φ 1 (z) with a fixed γ j is only valid in a certain sector S l ("Stokes sector") of the complex plane, while in the complementary sectors the asymptotic expansion requires γ k = γ j [21].The growth type of φ 1 (z) is then given by the sector and associated γ with maximal Re(γz 2 ).Now the spectral problem consists in the task to find those solutions of (9) which have finite Bargmann norm because the first Bargmann condition [22] is satisfied by all of them.The opposite case is realized in the linear QRM: The generic solutions are not entire in C but their Bargmann norm is always finite as for the harmonic oscillator [2].
If |φ 1 | < ∞, φ 1 (z) is an element of B, which is a necessary condition for |ψ to be an eigenvector of H 2p .To decide whether |φ 1 | is finite, it is sufficient to know the asymptotic behavior of φ 1 (z) in each Stokes sector S l , where θ(x) denotes the Heaviside step function.Inequality (13) holds for arbitrary R < ∞ because φ 1 (z) is entire.Now, for z ∈ S l and |z| > R, for sufficient large but finite R which means that the inequality on the right of (13) will be satisfied if If ( 15) is satisfied for all Stokes sectors S l , φ 1 will be normalizable.Because the γ j are non-degenerate, we have β j = 0 for all j, but this has no bearing on the normalizability properties of φ 1 (z).
The fact that all asymptotics have |γ| = 1 entails that in this region the spectrum has a "pure point" characteristic [29] and no continuous part if the continuous spectrum requires asymptotics with |γ| = 1.That this is likely correct follows from consideration of the known "generalized eigenstates" which are associated with the continuous spectrum [22], with the real phase angle 0 ≤ ϑ < 2π and x ∈ R.These states satisfy the orthogonality relation for fixed ϑ, needed to constitute a spectrum continuous in the parameter x.In fact, the image of , a generalized eigenstate of the position operator x, while the image of φ π/2 p reads ψ p (x) = (2π) −1/2 exp(ipx), a plane wave associated with the eigenbasis of the momentum operator p. Fourier transformation in L 2 (R) corresponds to the unitary transform z → iz in B, embedded in the continuous set of isometries of B, z → e iϑ z, each corresponding to a family of generalized eigenfunctions parameterized by the angle ϑ.Each of them is characterized by its exponent of the second kind γ = −e 2iϑ of order two, the generalized eigenvalue is proportional to the exponent of order one.The position basis (ϑ = 0) has γ = −1 and the momentum basis (ϑ = π/2) γ = 1.If the conjecture above is correct, the spectrum of H 2p in the region ω > 2 has pure point characteristic.

ω < 2
All γ j are located on the unit circle, |γ j | = 1 and are non-degenerate.
Because the γ j are all different, the β j vanish as in case I.As the arguments of the γ j differ, they do not have common Stokes sectors.The convergence of the integral (15) in a certain sector S l of the complex plane where ( 10) is asymptotically valid depends on whether S l contains the line {z|arg(±z) = −arg(γ)/2}.If not, the integral converges.If the line is contained in S l , the convergence depends on the value of the exponent of the first kind, ρ(S l ).As example, lets consider the sector The integral in (15) reads then (γ = e iϕ ) where the arc from x − ix to x + ix has been replaced by a straight line.We have with µ = max(0, −Im(ρ)π/2).It follows that (15) will be satisfied if Re(ρ) < −1/2.
In the present case, ρ is a function of γ, ω and E, With γ = e iϕ we find, This means Re(ρ) = −1/2 for γ 1,4 and Re(ρ) = −5/2 for γ 2,3 .Normalizable states exist if the exponents γ 1 , γ 4 are not present in the Stokes sectors containing their critical lines.The other possibility, ρ = −1/2 leads to a continuous spectrum, but the generalized eigenstates are not of the form given in (17), because the exponent β vanishes.States with these asymptotics appear in the harmonic oscillator with inverted potential [30] which is self-adjoint despite the continuous spectrum being unbounded from below, spanning the whole real axis [31,32,33].As H 2p can be written in terms of the inverted oscillator and its Fourier transform for ∆ = 0, it follows from the Kato-Rellich theorem [16] that H 2p is self-adjoint for ω < 2, contrary to a conjecture by Ng et al. [12].From the two possible values of Re(ρ), a discrete spectrum cannot be ruled out for ω < 2. However, its determination is not subject of the present paper.

ω = 2
Here, γ takes the values +1 and −1, both of them doubly degenerate.It follows from (9) together with (10) that now β is determined in terms of γ, For γ = 1, that means and for γ = −1, In both cases, the two-fold degeneracy of γ is lifted by the two possible values for β(γ), making all four asymptotic solutions linearly independent.There are now two different parts of the spectrum, E > −1 and E < −1.The threshold energy is E 0 = −1, in accord with [17,18].The four Stokes sectors are given in (16).The "critical" line with arg(±z) = −arg(γ)/2 is the real axis for γ = 1 and the imaginary axis for γ = −1.Let's first look at γ = 1 and E < −1.The integral in (15) reads in sector S 0 , as β is real and will only converge if β < 0, independent of ρ.In sector S 1 , we have instead which converges unconditionally.Analogously, the integral converges in sector S 2 if β > 0 and independent of β in S 3 .For γ = −1 and E < −1, the same description obtains with the real and imaginary axes interchanged (β ∈ iR).It follows that the spectral problem is equivalent to a lateral connection problem with non-trivial Stokes multipliers [21], because asymptotic solutions in sectors S 0 and S 2 (S 1 and S 3 ) must have different β in case γ = 1 (γ = −1) for a state with E < −1 to be normalizable.
There is seemingly also a second possibility, namely that γ = 1 in sectors S 1,3 and γ = −1 in sectors S 0,2 which also corresponds to off-diagonal Stokes matrices but involving the exponents of second kind with order two.We will see later that this case cannot occur, but the first one leads to a normalizable eigenstate if the parameter E is determined such that the lateral connection problem is solved by φ 1 (z).Indeed we know from the direct calculation [17,18] (see also [26]) that there are normalizable states below E < −1, contrary to the statement made in [24], where it is claimed that such states cannot exist for any E.
For E > −1 and γ = 1, β is imaginary and we have in sector S 0 , An estimate similar to (20) yields now that I 0 (R) will be finite if Re(ρ) < −1/2.The same condition applies to S 2 , while the integral converges in sectors S 1,3 .For g = −1 (real β) normalizability of a state with E > −1 imposes the same condition on ρ and S 0,2 replaced by S 1,3 .We find for γ = 1, ρ(1) = −8/5 and for γ = −1, ρ(−1) = 0.According to the argument for the case ω < 2, one would conclude that there are normalizable solutions with E > −1, as the integral in the critical sector converges, provided γ = 1 determines the asymptotics of φ 1 (z).This does not happen because no eigenstate of H 2p at ω = 2 can have φ 1 -asymptotics with γ = 1.Only γ = −1 is a valid exponent of second kind with order two for φ 1 (z) satisfying the system (7), (8).To see this, let's use the simplification of the system at ω = 2 mentioned in the introduction, namely The system ( 7),( 8) becomes It can be easily analyzed in the space L 2 (R) [17,18], but the asymptotics of the eigenstates can be inferred already from the special case ∆ = 0. We find that φ 1 (z) is an eigenstate of the position operator x with eigenvalue x 0 = √ E + 1/ √ 2 which agrees with the value of β in (26), the value ρ = 0 found above for γ = −1 and the form of the corresponding generalized eigenstate of x in (17).Therefore, only γ = −1 may occur in the asymptotic form of φ 1 .This is not changed for ∆ = 0 because φ 1 satisfies then a Schrödinger equation on the real line with a potential approaching zero asymptotically [17,18], therefore the asymptotic form of the scattering states are the same as for ∆ = 0. Likewise, for φ 2 , which is the image of φ 1 under Fourier transformation (or z → iz), its asymptotic form must have γ = 1.Because the formal solutions φ 1 (z) with γ = −1, real β and ρ = 0 cannot be made normalizable by adjusting E but are asymptotically generalized eigenstates of x for all E > −1, they must correspond to the continuous spectrum of H 2p in this range of energy.Moreover, as γ is fixed to −1 for φ 1 (z) in each Stokes sector, the option to have normalizable states for E > −1, namely γ = −1 in sectors S 0,2 and γ = 1 in sectors S 1,3 , is ruled out.It follows that no bound states are embedded into the continuum as it appears in the Rabi-Stark model [19].The same argument applies for E < −1, however, here we must assume ∆ = 0 to provide a potential allowing for bound states.The asymptotic form φ2 (x) of a solution to (33) ) which corresponds to (25) with two different β in sectors S 0 and S 2 , the first case mentioned above.The second case, different γ in sectors S 0,2 and S 1.3 is not realized.Therefore, the exponent of second kind with order two is unique throughout the complex plane and the Stokes phenomenon manifests only in the subleading exponent of order one.

Pure point spectrum for ω > 2
For ω > 2, H 2p has only discrete eigenvalues and normalizable eigenfunctions, because the exponents of second kind with order two are non-degenerate solutions of (11).The spectral condition amounts to the requirement that the formal solution of ( 7),(8) has only normalizable asymptotics in all four Stokes sectors (16).A sufficient condition for this to happen is clearly that only exponents γ 1,2 appear in the asymptotic expansion of φ 1 (z) which occurs only for a discrete set E n , n = 0, 1, 2 . .., of the parameter E.
We call any function G H (E) a "G-function" for H, if it has only zeros on the real axis, coinciding with all (respectively the regular subset of) the discrete eigenvalues where h(E) may be any function bounded from below such that the right hand side of ( 34) converges for all finite E. The spectral problem can be considered "solved" if one finds a G-function for H, which is, of course, not unique.Note that h(E) does not need to be bounded from above -which will turn out to be rather useful in section 3.2.
It is often convenient to define G-functions for different sectors of the spectrum related to invariant subspaces under the symmetries of H.In the case of H 2p , we have a Z 4 -symmetry generated by P = exp(i(π/2)a † a)σ x .The four eigenvalues of P are ±1, ±i and label four invariant subspaces of H.The Z 2 -subgroup generated by P 2 acts only in the bosonic part of H with eigenvalues ±1 corresponding to even and odd functions of z.We denote each invariant subspace as H ab with a, b ∈ {+, −} corresponding to the eigenvalues of P 2 and P : a = + [−] for φ 1,2 (z) even [odd] and b = ± if φ 1 (iz) = ±φ 2 (z) for even φ 1,2 (z), resp.φ 1 (iz) = ±iφ 2 (z) for odd φ 1,2 (z).These relations between φ 1 (z) and φ 2 (z) apply to their asymptotic forms (14) and therefore relate the asymptotics of φ 1 (z) in Stokes sectors S 0 and S 2 to the asymptotics of φ 2 (z) in sectors S 1 and S 3 and vice versa.The asymptotic relations among Z 4 -eigenstates will be of importance in the following.
The eigenstates of H 2p belonging to different subspaces may have degenerate eigenvalues at certain points in the parameter space and manifest as level crossings in the spectral graph.In the 2pQRM the level crossings are of two different types: The first, the so-called Juddian case, describes level crossings between H ++ and H +− , respectively between H −+ and H −− .It has been investigated in [34,35].The two-fold degeneracy of the Juddian solutions corresponds to a two-dimensional representation of Z 4 /Z 2 ∼ Z 2 .The second type describes also two-fold degeneracies but between subspaces with different eigenvalues of P 2 .They have been studied in [36].In both cases the eigenfunctions can be given in terms of either elementary functions or functions belonging to the hypergeometric class.The degenerate eigenvalues belong to the exceptional spectrum [2] and are not the subject of the present paper which concerns the regular spectrum comprising the non-degenerate eigenvalues.It has to be noted that the non-degenerate spectrum may contain certain eigenvalues at isolated points in parameter space (the so-called non-degenerate exceptional spectrum) which have to be dealt with separately, see [7] and section 3.2.We confine the following discussion to the subspace H ++ , containing only even functions of z.

The G-functions based on recurrence relations
The eigenvalue +1 of Z 4 /Z 2 means that a solution of ( 7), (8) belonging to H ++ satisfies φ 1 (iz) = φ 2 (z).If φ 1,2 (z) are even, their expansion in powers of z around z = 0 reads with arbitrary a 0 = 0. From (7) we deduce the following three-term recurrence relation for the a n , and a n = 0 for n < 0. This recurrence relation is the starting point for the G-functions proposed in [23] and [24].A Poincaré analysis of (36) with the asymptotic ansatz leads to the equation for x, with solutions As |x − | > |x + |, the Perron-Kreuser theorem [37] states that there exist a minimal solution to the recurrence (36) with and (37) suggests that the asymptotic form of φ 1 (z) for |z| → ∞ reads in this case where r(z) grows at most linear in z but is undetermined by the asymptotics of a n .A function with these leading asymptotics in all Stokes sectors corresponds to a normalizable solution of (9).Because Pincherle's theorem states that the minimal solution of a three-term recurrence relation, if it exists, is determined by a convergent continued fraction [37] (see also [38]), one may write down a G-function for the spectral problem in H ++ as follows [23], with and (E − ∆)/2 = a 1 /a 0 , as determined from the initial condition for (36).The vanishing of G Z (E) at some point E s means that the minimal solution of (36) whose first terms are given by the continued fraction V 1 (E s ) = a min 1 /a min 0 matches the initial condition a n = 0 for negative n, which is dictated by the analyticity of φ 1 (z) at z = 0.The state is therefore normalizable because the asymptotic behavior of φ 1 (z) apparently contains only the exponent γ 1 , as guaranteed by the minimality of the solution to (36).It follows that E s belongs to the pure point spectrum of Zhang's approach is certainly the simplest and most elegant way to obtain a G-function for H 2p in the regime ω > 2. Following Okubo [39], Maciejewski and Stachowiak [24] have proposed a method which is based as well on the determination of those solutions of a three-term (matrix) recurrence relation which correspond to normalizable asymptotics for φ 1,2 because they contain only exponents γ 1,2 .After "gluing" them to the initial conditions as above, they obtain another G-function ((43) in [24]) which is given in terms of a contour integral over the solution of an auxiliary differential equation.
Coming back to (38), we have |x + | = |x − | for ω ≤ 2, a minimal solution does not exist and the continued fraction does not converge.One would thus conclude with the authors of [24] that in this regime there are no normalizable solutions of ( 7), (8) satisfying in addition φ 1 (iz) = φ 2 (z).However, as we have seen in section 2, the pure point spectrum is not empty for ω = 2 and normalizable states could exist also for ω < 2. The discrepancy is resolved by noting that the Poincaré analysis of the asymptotic behavior of the recurrence relation only yields the leading, most divergent term in the asymptotic expansion and says nothing about the subleading terms which determine normalizability in the case |γ| = 1.The formal solution (35) may or may not be normalizable even if no minimal solution of (36) exists.Another consequence is the missing Stokes phenomenon: The minimal solution of the recurrence relation for φ 1 (z) fixes the allowed exponent to γ 1 , independent of the argument of z. φ 1 (z) should approach zero in sectors S 0,2 whereas φ 2 (z) diverges there, obviously contradicting (7) if ∆ = 0.In fact, the form (41) of the asymptotic behavior of φ 1 (z) is only valid in sectors S 1,3 where the term exp(γ 1 z 2 /2) diverges for |z| → ∞.In the sectors S 0,2 , where it is recessive, subdominant contributions in the recurrence (36) determine the asymptotic behavior of φ 1 (z), given in these sectors as φ 1 (z) ∼ exp(γ 2 z 2 /2).Because the asymptotics of the recurrence relation does not provide the correct behavior of φ 1 (z) in all Stokes sectors, one may ask whether indeed all normalizable states satisfying ( 7), ( 8) are related to the minimal solution of (36).For certain values of E, there could be an entire function φ(z) behaving asymptotically like exp(γ 4 z 2 /2) in sectors S 0,2 and as exp(γ 3 z 2 /2) in S 1,3 .This function would be normalizable although it belongs to a dominant solution of (36).If this happens, there would be elements of the spectrum not given by zeros of G Z (E).The same argument applies to the G-function of [24].The phenomenon occurs (for different reasons) in the QRM with linear coupling, where the non-degenerate exceptional spectrum cannot be obtained by Schweber's G-function based on continued fractions [40].
In the present case of the 2pQRM, such an exceptional eigenstate cannot occur, at least not in the non-degenerate part of the spectrum, due to the relation φ 1 (z) = αφ 2 (iz) with α ∈ {±1, ±i}.It entails that if φ 1 (z) has strict recessive behavior in a given sector with exponent γ, φ 2 (z) will have exponent −γ the same sector and would not be normalizable if γ is either γ 3 or γ 4 .Indeed, the Z 4 -symmetry of the solution shows that exponents γ 1 and γ 2 appear in the asymptotics of both φ 1 and φ 2 in all sectors but not γ 3 or γ 4 if the state is normalizable.Thus we may conclude that G Z (E) yields the complete spectrum for ω > 2. The argument is less clear for the G-function based on Okubo's method [24].Because the factorial power series in x + is actually a double series (see (34) in [24], x + is called u 0 ), it is not obvious that the asymptotics have the form φ 1 (z) ∼ exp(±x + z 2 ) in all sectors.
Zhang's G-function is shown in Figure 1.It has poles and zeros which move closer together for higher values of E, as typical for a continued fraction.The resolution of higher zeros poses thus known numerical problems but the major drawback of G Z (E) is the lack of qualitative knowledge about the distribution of its zeros, constituting the regular spectrum of H 2p in H ++ .The poles are all of first order and G Z (E) is monotonous, but the location of the poles is just as unknown as the zeros -they do not allow to restrict the distribution of the zeros as the G-functions for the QRM which possess poles at known positions, namely E = nω, with n a non-negative integer [2,41,42].The G-function proposed in [24] has the same problem: The exact eigenvalues of H 2p can be extracted by a numerical procedure but this has no advantage compared with direct diagonalization in a truncated Hilbert space because the graph of the G-function has no properties, say the position of maxima or minima, which are known in closed form (see Figures 1 and 2 in [24]).On the other hand, exact diagonalization in a state space of finite dimension suggests the collapse of the complete spectrum to a single point at ω c , a numerical artifact of the truncation.But the G-functions based only on the recurrence relation (36) offer no option to study the pure point spectrum arbitrarily close to ω c in a qualitative way.

The G-functions based on the Z 4 -symmetry
The recurrence (36) implements the Z 4 -symmetry directly via the relation (35) between the series expansions of φ 1 (z) and φ 2 (z).It is therefore automatically satisfied by all solutions of (36), independent of their normalizability.The question arises whether the Z 4 -symmetry can be used to discern normalizable from non-normalizable solutions as in the linear QRM.In the latter case, the two regular singular points of the eigenvalue equation are located at z = ±g and mapped onto each other by the Z 2 -symmetry z → −z of that model.This can be used to construct a G-function which has only zeros at those energies corresponding to solutions φ 1 (z),φ 2 (z) which are analytic in the whole complex plane by using Frobenius expansions which are analytic in the vicinity of only one of the singular points -the symmetry requirement φ 2 (z) = φ 1 (−z) is then equivalent to analyticity at the other singular point [2].
The situation is different for the 2pQRM: Here there is only a single singular point at z = ∞ but with higher rank, preventing all formal solutions to be normalizable, although all are entire.The invariant points of the map z → iz are z = 0, ∞ as in the QRM and z = 0 is an ordinary point of the system (7), (8).The symmetry generator acts trivially on the point z = ∞ and maps the asymptotics of φ 1,2 (z) onto each other: γ 1 ↔ γ 2 and γ 3 ↔ γ 4 .Clearly, because z → iz is an isometry of B, the normalizable forms are mapped onto themselves, as well as the non-normalizable ones.
Travěnec [43] has attempted to use the Z 4 -symmetry to derive a G-function for H 2p in analogy with the linear QRM as follows.First, separate the asymptotic factor exp(γ 1 z 2 /2) from the expansions of φ 1,2 , The coefficients a j n satisfy the coupled recurrence relations with initial condition a 1 0 = a 2 0 = 1.The symmetry φ 1 (iz) = φ 2 (z) is not assumed beforehand in the system (45), (46).The G-function is then defined as with z 0 = 0 an arbitrary point in the complex plane.The zeros of G T (E) should not depend on the choice of z 0 .The problem with this approach is that the implementation of the symmetry condition in (47) cannot discern between normalizable and nonnormalizable solutions.(Non-)normalizability is an invariant property of any function under z → iz, as we see from the behavior of the possible asymptotic forms.Moreover, the initial condition φ 1 (0) = φ 2 (0) enforces the symmetry for any solution of ( 45), (46).The function G T (E) vanishes therefore identically [44].Nevertheless, a straightforward numerical implementation of (47) yields a function which does not vanish identically and possesses zeros at the exact eigenenergies of H 2p in H ++ .
Figure 2 shows a plot of G T (E) for the same parameters as in Figure 1.The zeros of G Z (E) and G T (E) coincide, but the non-zero values of G T (E) are extremely large.Both observations, the numerical validity of G T (E) and its large values, have the same origin, namely the fact that a three-term recurrence relation has a unique minimal solution.The numerical evaluation of ( 45), (46) will give valid results only if the initial conditions including the parameter E are fine-tuned to yield the minimal solution [37,45].The evaluation of dominant solutions is always numerically unstable, leading to the exponential accumulation of rounding errors which produce the huge non-zero values of G T (E) for any value of E not belonging to the spectrum.Only if E coincides with some E n , the recurrence yields the (normalizable) minimal solution and the numerics will reproduce G T (E) = 0, the value of G T (E) for any E if evaluated exactly.
The same mechanism appears also in the G-function for the oscillator with quartic anharmonicity proposed by Lay [46].Travěnec's G-function is thus numerically exact in the same sense as Lay's method to compute the spectrum of the quartic oscillator.But as it is only non-zero due to numerical instabilities, G T (E) cannot be used to draw qualitative inferences about the spectrum in an analytical way because it has no exactly known properties.In this respect, there is no difference between Travěnec's [43], Zhang's [23] or Maciejewski's G-functions [24].The root of the problem is the invariance of the singularity of ( 7),( 8) and its asymptotic solutions under the symmetry operation z → iz.Therefore, this system of differential equations has to be transformed in such a way that the symmetry acts non-trivially on its singularity structure.
The bosonic Bogoliubov transformation is well established in quantum optics to describe squeezed light [47] and was used by Chen et al. in [25] to derive a G-function for the 2pQRM.The derivation was heuristic in the sense that it did not establish a = ch(θ)a † + sh(θ)a.Applying I θ in the Bargmann space B to the system ( 7),( 8) leads to 2 + ω 2 zϕ (1) The transformed system (51),(52) has a singularity structure different from ( 7), (8) but it is still invariant under z → −z, so the solutions are even or odd functions of z.Upon elimination of ϕ 2 (z), we obtain an equation of third order for ϕ 1 (z), with abbreviations The equation (54) has a regular singular point at z = 0 and an unramified irregular singular point at z = ∞ of rank two and class two.The indicial equation at z = 0 has three solutions corresponding to exponents ρ = 0, 1, E 1 /ω 1 in the Frobenius expansion around z = 0.The first two correspond to even and odd solutions analytic at z = 0.The third solution with ρ = E 1 /ω 1 is only entire if E 1 /ω 1 is a non-negative integer.This is in close analogy to the linear QRM, where the exceptional spectrum is characterized by E/ω ∈ IN 0 [7].For now we assume E 1 /ω 1 / ∈ Z.The other cases will be discussed later.The irregular singular point is z = ∞, as in (9).Rank two entails that the asymptotics of ϕ 1 (z) have again the form (10).The equation for γ reads now with solutions 0, ω/2 and 2/ω.The first solution γ = 0 corresponds to the fact that the class of the singularity is not maximal [27].The indicial equation has degree one, meaning that there is a solution ϕ 1 (z) where z = ∞ could be a regular singular point.
The asymptotic expansion has then the form Plugging this expansion into (54), we find that all odd coefficients vanish, c 2m+1 = 0, for m ∈ IN 0 , and the exponent of first kind ρ has the unique value ρ = E 1 /ω 1 .The even coefficients c 2m are determined recursively for n ≥ 2 and c −2 = 0.The coefficients in the three-term recurrence relation (59) are diverging for n → ∞, therefore the convergence radius of (58) is zero.There is, as expected, no solution of (54) regular at z = ∞ if ρ / ∈ Z.Nevertheless, there is a certain Stokes sector S * where (58) is asymptotically valid for z ∈ S * and |z| → ∞.This means that in S * the branching behavior of ϕ 1 (z) is correctly described by the exponent ρ / ∈ Z which is the same as the non-integer exponent ρ = E 1 /ω 1 of the third Frobenius solution in the vicinity of z = 0.The global behavior of this solution can be described by a branch-cut running from z = 0 to z = ∞ as shown in Figure 3.The two other solutions have non-zero exponents of second kind γ and β = 0.One of them (γ = ω/2) has not normalizable asymptotics and the other (γ = 2/ω) is normalizable if it is asymptotically valid in S 0,2 .For generic values of E, both solutions analytic at z = 0 will exhibit the exponent ω/2 in sectors S 0,2 and are thus not elements of B. If E belongs to the spectrum of H 2p in H ++ , the even solution of (54) will be entire and normalizable with asymptotics described by γ = 2/ω in S 0,2 .The asymptotic exponent ρ for γ = 0 depends on E and γ but has no influence on the normalizability of ϕ 1 .Because γ = 0, its behavior close to z = ∞ is determined by the essential singularity of the factor exp((γ/2)z 2 ) which is isolated if ϕ 1 is entire.Up to now, we have not employed the Z 4 -symmetry which is implemented in B by the operator T as follows T is therefore an element of SL 2 (C) like the Bogoliubov transformation I θ .To compute and to find the normal ordered form of any group element g in SL 2 (C), it is convenient to use the two-dimensional representation with We find then and To compute the image of ϕ 1 (z) under T θ , we write such that f (z) → c 0 = 0 and g(z)e −(γ/2)z 2 → 0 for |z| → ∞ in S 0,2 .The asymptotic behavior of the image T θ [ϕ 1 ](z) is determined by the diverging factor in (65), which can be expressed through an element of SL 2 (C) acting on the constant function, Using now the two-dimensional representation of T θ and X, we obtain with It follows Solutions of (54) are mapped by T θ onto solutions of which is obtained from (51),(52) after elimination of ϕ 1 .It has the same singular points as (54), the exponents of first kind at z = 0 are 0, 1 and E 1 /ω 1 + 2. As usual, the non-integer exponent differs from the one for ϕ 1 by an integer, so that these solutions for ϕ 1 and ϕ 2 fulfill together the system (51),(52).Clearly, the exponents of second kind at z = ∞ are the same as in (54), with ρ = E 1 /ω 1 − 2. Therefore, the branching structure of the third solution of (71) with γ = 0 has the form shown in Fig. 3 for the associated solution of (54).
Lets assume now that the exponent of φ 1 (z) in ( 65) is γ = 2/ω.The exponent of its image under T θ given by ( 69) is zero which means that the dominant part of φ 1 (z) is mapped onto the recessive part g(z) of ϕ 2 (z) in the decomposition (65).Because T θ is an isometry, normalizable solutions are mapped onto normalizable ones, so the part diverging as exp(z 2 /ω) is mapped onto a less diverging part of a normalizable solution of (71) which must be analytic at z = 0. On the other hand, if the exponent of ϕ 1 (z) is γ = ω/2, the exponent of the image T θ [ϕ 1 ](z) would be infinite according to (69).This indicates that the corresponding solution of (71) cannot be entire, otherwise it would be described by either γ = ω/2 or γ = 2/ω in sectors S 0,2 .The only possibility is that non-normalizable solutions of (54) which are analytic at z = 0 are mapped by T θ onto solutions of (71) which are not analytic at z = 0, with non-integer exponent ρ = E 1 /ω 1 + 2. Recall that those solutions possess an essential but non-isolated singularity at z = ∞, because the asymptotic expansion (58) is valid in the Stokes sector S * .Their behavior at infinity cannot be captured by an entire function of order two, indicated by the diverging η x (γ).They are not elements of B already due to their non-analyticity at z = 0.
It follows that T θ maps the set of normalizable solutions of (54) onto the corresponding set of solutions of (71).These solutions are all analytic at z = 0 and have exponent γ = 2/ω in sectors S 0,2 .The non-normalizable solutions which are analytic at z = 0 but have exponent γ = ω/2 are mapped onto solutions which are not analytic at z = 0 but have there a branch-cut with exponent E 1 /ω 1 + 2. We can now follow the same strategy as for the linear QRM and construct a pair of even functions ϕ 1 (z), ϕ 2 (z) satisfying the system (51),(52) and both analytic at z = 0.If ϕ 2 (z) and ϕ 1 (z) have the expansions it follows from (51) that the coefficients are related by ān if E 1 /ω 1 is not an integer.The coefficients of ϕ 1 are uniquely fixed in terms of the coefficients of ϕ 2 if both are analytic at z = 0 without invoking the symmetry.From (52) one deduces the three-term recurrence relation for the a n , The initial conditions can be chosen as a −2 = 0, a 0 = 1.A Poincaré analysis of (74) shows that the convergence radius of the expansions in (72) is infinite, as it should be.Both ϕ 1 and ϕ 2 are entire and their asymptotics in S 0,2 have either the exponent γ = 2/ω if they are normalizable or the exponent ω/2 if not.The symmetry imposes now the additional relation between ϕ 1 and ϕ 2 , for all z ∈ C. If φ 1 (z) has asymptotics with exponent γ = ω/2 in S 0,2 , rendering it non-normalizable, it will be mapped onto a function which is not analytic at z = 0 and (75) cannot be satisfied, because ϕ 2 (z) is analytic by construction.It follows that the solutions of ( 73) and (74) cannot satisfy (75) for generic values of E because they both have asymptotics with γ = ω/2.If E is element of the regular spectrum of H 2p and corresponds thus to a certain eigenvalue of the Z 4 -symmetry, condition (75) will be satisfied because the image of ϕ 1 under T θ is a normalizable solution of (71) which must be analytic at z = 0. Therefore, E is an element of the spectrum with eigenstate in H ++ if and only if ϕ 1,2 as determined by ( 73),( 74) satisfy (75).This argument is independent from the presence of the Stokes phenomenon because only sectors S 0,2 determine the normalizability as ω/2 > 0 and 2/ω > 0. The asymptotics in sectors S 1,3 play no role, in contrast to the original situation described by (9) where the exponents of second kind take positive and negative values.
To derive the G-function explicitly, we must compute the action of T θ given in (64) on the series expansion of ϕ 1 (z) in (72), i.e. on even powers of z, where P m ω (z) is a polynomial in z of degree 2m, It follows It is evident that expression (78) becomes most simple at z = 0, and condition (75) would read where the a 2m are recursively determined by (74).However, a Poincaré analysis of the series shows that the convergence radius R(f ) = 1, meaning that f (w) is not defined at w = 1, which is supposed to yield T θ [ϕ 1 ](0).This is easy to understand because the generic form of T θ [ϕ 1 ](z) has a branch-cut at z = 0 and each series expansion of it diverges there.The analysis of (81) gives thus another independent proof of the fact that a non-normalizable but entire solution of ( 54) is mapped onto a solution of (71) which is not analytic at z = 0.
To derive a G-function whose series expansion does converge, we apply I −1 θ to both sides of (75), The expressions in (82) are well defined at z = 0, so it is possible to evaluate the G-function at the invariant point z 0 = 0 of the map T .This has the advantage that φ 2 (z 0 ) − φ 1 (iz 0 ) = 0 is sufficient to determine the regular spectrum in H ++ .If z 0 = 0, one would have two independent conditions φ 2 (z 0 ) = φ 1 (iz 0 ) and φ 2 (−iz 0 ) = φ 1 (z 0 ) [45].Using (63), we find for The Poincaré analysis of the series yields the convergence radius R( f ) = 1+ √ ω 2 − 4/ω > 1, therefore the expression (83) is given by an absolutely converging series.A fortiori the same applies to I −1 θ [ϕ 1 ](0) and we can define a G-function as G ++ (E) = φ 2 (0) − φ 1 (0), that is After adjusting notation, we find G ++ (E) = G 1/4 + (x(E)) which is the G-function for sector H ++ given in ( 16) of [26] and has been originally proposed ten years ago by Chen et al. in [25].By examination of the normalizability properties of formal solutions to the Schrödinger equation in the Bargmann space, we have demonstrated that Chen's G-function yields the complete regular spectrum in H ++ (the other sectors require minor modifications of the derivation).It is shown in Figure 4.
The major qualitative difference of this G-function compared with the others proposed by Zhang [23], Travěnec [43] and Maciejewski/Stachowiak [24] is the presence of simple poles at even integer values of E 1 /ω 1 , that means at energies For most values of the parameters ω and ∆ all eigenenergies are regular and located between the poles given by (86).According to the conjecture formulated in [2], the number of zeros between adjacent poles is restricted to be 0, 1 or 2. This conjecture has been numerically verified in the linear and 2pQRM and can be derived analytically for small ∆.It plays here a similar role as the string hypothesis in systems solvable by the Bethe ansatz [48].As discussed in [26], the conjecture entails that the average distance between eigenenergies in H ++ is given by the distance between the poles, 2 √ ω 2 − 4. The poles are equally spaced along the whole real axis for and become more dense as ω approaches the critical value ω c = 2 from above -the distance between adjacent energy levels goes to zero in this limit.At ω c , the real axis above the threshold value E 0 = −1 becomes densely covered with levels, indicating the transition to a continuous spectrum as derived in section 2 for the critical point itself.In this way the qualitative properties of the spectrum as the system approaches ω c can be deduced from the known pole structure of G ++ (E).There is no "collapse" of the spectrum at the point ω c to the single energy E 0 as falsely inferred from exact diagonalization in a truncated state space but instead the transition to a continuous spectrum at ω c , spanning the whole real axis from E 0 to infinity.The exceptional spectrum in H ++ consists of those energy levels which coincide with a pole energy E (m) in full analogy with the linear QRM.In this case the pole in G ++ (E) is lifted by a zero of the expansion coefficients in (85) and G ++ (E (m) ) is finite.It is easy to see that those eigenstates are exactly the special solutions found in [34].They appear if the asymptotic expansion (58) of the third solution to (54) is convergent because the recurrence (59) breaks off after m steps, leading to another normalizable entire solution, a polynomial.The exceptional level is thus doubly degenerate.We have seen theat the presence of poles in G ++ (E) is instrumental to determine the approach to the collapse point without numerical evaluation of the g-function.On the other hand, it may seem detrimental to have access only to the regular part of the spectrum.As shown in [49,50], it is possible to obtain a G-function for the linear QRM whose zeros give the complete spectrum including the degenerate and non-degenerate exceptional parts by multiplying G ± (E) with appropriate Gamma-functions.The same is possible for G ++ (E).This may be of importance to make connections with the spectral zeta-functions associated with H R and H 2p [9,50].
The degeneracies of Juddian type, between states in H ++ and H +− , are connected to the poles of G ++ and G +− , where E 1 /ω 1 is a positive even integer.Likewise the degeneracies between H −+ and H −− correspond to lifting of a pole in G −± .The poles correspond to E 1 /ω 1 being a positive odd integer.The degeneracies appearing between even and odd eigenstates with E 1 /ω 1 a half-integer are not reflected in these G-functions.They correspond to a factorization property of ( 9) for special values of the parameters ω, ∆ which does not follow from the Z 4 -symmetry alone [36].

Conclusions
Many qualitative aspects of the spectrum of the two-photon quantum Rabi model which define the three different parameter regimes, ω > 2, ω = 2 and ω < 2, can be inferred without numerical computation by examining the singularity structure of the Schrödinger operator in the Bargmann space of entire functions.In these regimes, the spectrum is either pure point (ω > 2) or comprises a continuous and a discrete part (ω = 2).For ω < 2 the spectrum is at least partly continuous with generalized eigenstates resembling asymptotically those of the inverted harmonic oscillator.At the "collapse point" ω = 2, the effective potential following from (32), (33) depends on the energy eigenvalue so that no direct mapping to a one-dimensional problem without spin degree of freedom is possible.Nevertheless qualitative aspects of the spectrum may be derived from the analogy (section 2.3).
The various generalized spectral determinants proposed in the regime ω > 2 have different qualitative properties depending on their derivation.While Zhang's G-function [23], based on a continued fraction, and Chen's G-function [25] based on a Bogoliubov (squeezing) transformation, both possess poles on the real axis located between the zeros (eigenvalues of H 2p ), their position is known analytically only in the latter case.Neither Travěnec's [43] nor Maciejewski's G-functions [24] exhibit poles or other known qualitative features, like the position of maxima or minima, to infer the distribution of the zeros without explicit numerical computation.The average distance of energy levels, while independent from the coupling in the linear QRM, changes dramatically in the 2pQRM if the critical coupling g = ω/2 is approached from below (or ω → 2 from above for fixed g = 1).The levels coalesce on the real axis above the threshold energy E 0 = −1 to form the continuous part of the spectrum at the critical point ω = 2.Because the average distance of levels is given by the pole distance in Chen's G-function, ∆E = 2 √ ω 2 − 4, the exponents of this critical behavior can be computed exactly [26].Moreover, Chen's G-function yields the quasi-exact Juddian points of the spectrum in the same way as the corresponding G-function of the linear QRM.While the validity of this G-function has been confirmed numerically in great detail, a rigorous derivation based on the normalizability of the wave functions was missing.As the numerical checks fail close to the critical point due to truncation of the Hilbert space, the rigorous proof of validity provided here seems not only of formal but also of practical importance.

Figure 1 .
Figure 1.Zhang's G-function in sector H ++ for parameters ω = 2.5 and ∆ = 0.7.The poles of G Z (E) move closer to its zeros En for larger E.

Figure 2 .
Figure 2. Travěnec's G-function for ω = 2.5 and ∆ = 0.7 with z 0 = 5 + 5i.The real (magenta) and imaginary (green) parts of G T (E) share the same zeros which are located at the eigenenergies of H 2p .

Figure 3 .
Figure 3.The global behavior of the third solution of (54) for E 1 /ω 1 / ∈ Z.The branch-cut runs from zero to ∞ within the sector S * .

Figure 4 .
Figure 4. Chen's G-function G++(E) for the same parameters as in Figs.1 and 2. This function is characterized by poles at even integer values of E 1 /ω 1 .The regular zeros are located between the poles which determine the average distance between adjacent zeros.
(12)s easy to see that the integral in(15)is finite for |γ| < 1 in all of C and diverges for |γ| > 1 in the Stokes sectors with Re(γz 2 ) > 0, independent of β and ρ.If |γ| = 1, the convergence of the integral is determined by β = 0.If γ lies on the unit circle and β vanishes,(15)is determined by ρ.From(12)we can infer three regions in the range of ω with different properties of (15).= −γ 2,1 have finite Bargmann norm in all sectors.Solutions with γ 4 < −1 converge in sectors S 0 and S 2 and diverge in sectors S 1 and S 3 while solutions with γ 3 = −γ 4 > 1 diverge in sectors S 0 and S 2 and converge in sectors S 1 and S 3 .