A generalization of the quantum Rabi model: exact solution and spectral structure

TL;DR

This paper presents an exact solution to a generalized quantum Rabi model with an added Stark-like term, revealing unique spectral features in ultrastrong coupling regimes.

Contribution

It extends the quantum Rabi model by including a Stark-like coupling and provides an exact analytical solution using Braak's method.

Findings

Avoided level crossings for same parity states

Rapid onset of near-degenerate levels with increasing coupling

Spectral features absent in the original model

Abstract

We consider a generalization of the quantum Rabi model where the two-level system and the single-mode cavity oscillator are coupled by an additional Stark-like term. By adapting a method recently introduced by Braak [Phys. Rev. Lett. {\bf 107}, 100401 (2011)], we solve the model exactly. The low-lying spectrum in the experimentally relevant ultrastrong and deep strong regimes of the Rabi coupling is found to exhibit two striking features absent from the original quantum Rabi model: avoided level crossings for states of the same parity and an anomalously rapid onset of two-fold near-degenerate levels as the Rabi coupling increases.

Equations

$${\mathcal{H}}_{\mathrm{Rabi}}$$

$$=$$

$$\mathcal{H}$$

$$=$$

$$\mathcal{H}=\omega a^{†}a+\Delta {\sigma }^x+g{\sigma }^z(a+a^{†})+\gamma {\sigma }^x{a}^{†}a.$$

$$⟨u|v⟩=\frac{1}{\pi }{\int }_{\mathbb{C}}d\Re (z)d\Im (z)\overline{u(z)}v(z)e^{-z\bar{z}},$$

$${\mathcal{H}}=\left(\begin{array}[]{cc}z\frac{d}{dz}+g\left(z+\frac{d}{dz}\right)&\gamma z\frac{d}{dz}+\Delta\\ \gamma z\frac{d}{dz}+\Delta&z\frac{d}{dz}-g\left(z+\frac{d}{dz}\right)\end{array}\right).$$

$${\mathcal{U}}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}1&1\\ {\mathcal{T}}&-{\mathcal{T}}\end{array}\right),$$

$${\mathcal{U}}^{-1}{\mathcal{H}}\,{\mathcal{U}}=\left(\begin{array}[]{cc}{\mathcal{H}}_{+}&0\\ 0&{\mathcal{H}}_{-}\end{array}\right)$$

$${\mathcal{H}}_{\pm }=z\frac{d}{dz}+g\left(z+\frac{d}{dz}\right)\pm \left(\gamma z\frac{d}{dz}+\Delta \right)\mathcal{T}$$

$${\mathcal{H}}_{\pm }{\psi }^{(\pm )}(z)=E_{\pm }{\psi }^{(\pm )}(z)$$

$$z\frac{d}{dz}{\psi }^{(\pm )}(z)+g\left(z+\frac{d}{dz}\right){\psi }^{(\pm )}(z)\pm \left(\gamma z\frac{d}{dz}+\Delta \right){\psi }^{(\pm )}(-z)=E_{\pm }{\psi }^{(\pm )}(z).$$

$$\mathcal{P}=(-1{)}^{a^{†}a}{\sigma }^z.$$

$$|\psi ⟩=|E,p⟩.$$

$$\varphi (z)\equiv \psi (z)\text{mbox}and\bar{\varphi }(z)\equiv \psi (-z),$$

$$(z+g)\frac{d}{dz}\varphi (z)+\left(gz-E\right)\varphi (z)+\gamma z\frac{d}{dz}\bar{\varphi }(z)+\Delta \bar{\varphi }(z)$$

$$(z-g)\frac{d}{dz}\bar{\varphi }(z)-\left(gz+E\right)\bar{\varphi }(z)+\gamma z\frac{d}{dz}\varphi (z)+\Delta \varphi (z)$$

$$\Gamma (z){\varphi }^{\prime }$$

$$\Gamma (z){\bar{\varphi }}^{\prime }$$

$$\Gamma (z)=(1-{\gamma }^2)(z-w)(z+w)$$

$$\Lambda (z)$$

$$\mathcal{E}(z)$$

$$\varphi (z)=\sum _{n=0}^{\infty }{A}_n(z-z_0{)}^{n+r}$$

$$r=r_1=\frac{E+g^2+\Delta \gamma }{1-{\gamma }^2}\equiv x_{\gamma }\ge 0$$

$$r=r_2=0,$$

$$\bar{\varphi }(z)=\sum _{n=0}^{\infty }{\bar{A}}_n(z-z_0{)}^{n+r}$$

$$\varphi (z)$$

$$\bar{\varphi }(z)$$

$$\frac{d\varphi (z)}{dz}$$

$$\frac{d\bar{\varphi }(z)}{dz}$$

$$(1-{\gamma }^2)(y-2w)y{\rho }^{\prime }$$

$$(1-{\gamma }^2)(y-2w)y{\bar{\rho }}^{\prime }$$

$$K_2$$

$$K_1$$

$$K_0$$

Full text

See pages 26-27 for an added corrigendum with corrected equations (39) and (40), and new figures replacing figures 4 and 6. The general conclusions drawn in the paper remain unchanged and valid.

Keywords: quantum optics, quantum Rabi-Stark model, regular and exceptional spectrum

1 Introduction

The integration of coherent nanoscale systems with quantum resonators is a focal point of current quantum engineering of states and devices. Examples range from trapped ions interacting with a cavity field [1] to superconducting charge qubits in circuit QED architectures [2]. The paradigmatic model for these systems is the Rabi model [3] which was first introduced 80 years ago to discuss the phenomenon of nuclear magnetic resonance in a semi–classical way. While Rabi treated the atom quantum mechanically, he still construed the rapidly varying weak magnetic field as a rotating classical field [4].

In the course of investigating the relationship between the quantum theory of radiation and the corresponding semi-classical theory, Jaynes and Cummings [5] discussed a model similar to Rabi’s. However, their model of an idealized atom consisting of only two levels coupled to a single quantised oscillator mode in an optical cavity was now a fully quantum mechanical model, the quantum Rabi model (sometimes also designated as the quantum electrodynamic Rabi problem [4]). Jaynes and Cummings [5] also introduced an important approximation to the quantum Rabi model, the so–called rotating wave approximation (RWA), leading to a model which can be solved exactly by elementary means and which now bears the name quantum Jaynes–Cummings model.

The quantum Rabi model, on the other hand, although still describing the interaction between matter and light in one of the simplest ways, only recently yielded to an exact and complete analytical solution [2] when Braak found an ingenious way to exploit the underlying $\mathbb{Z}_{2}$ parity symmetry of the model to derive its energy spectrum. While the quantum Jaynes–Cummings model has sufficed for a long time to describe experiments in quantum optics, recently it has become more and more necessary to go beyond the RWA as the larger Rabi coupling strengths of the ultrastrong and deep strong regimes come within experimental reach [7, 8].

In connection with his investigation of the exact solvability of the quantum Rabi model, Braak also developed a new proposal for quantum integrability [2]. This proposal is of considerable importance in view of the ongoing quest for a consistent notion of quantum integrability [9, 10, 11, 12].

Concurrent with this theoretical breakthrough, and motivated mostly by novel experimental setups, there has been an avalanche of studies of the quantum Rabi model and its many generalizations, revealing a plethora of intriguing and intrinsically nonclassical effects (for a recent review, see [13]).

A particularly interesting generalization of the model was proposed by Grimsmo and Parkins in 2013 [14]. These authors inquired about the possibility to realize the quantum Rabi model with a single atom coupled to a high–finesse optical cavity mode. They arrived at a scheme where two hyperfine ground states of a multilevel atom emulate an effective two–level system, with resonant Raman transitions between the two states induced by the cavity field and two auxiliary laser fields. Importantly, this scheme allows for a realization of the quantum Rabi model where coupling constants and effective frequencies can be freely and independently tuned, opening an experimental inroad to systematically probe also the ultrastrong and deep–strong coupling regimes. These are the regimes where the Rabi model comes into its own, while the time–honoured RWA $-$ which allowed the Rabi model to be replaced by the much simpler Jaynes–Cummings model [5] $-$ breaks down.

For generic values of the parameters of the model, however, the Grimsmo–Parkins scheme requires the addition of a new term to the quantum Rabi Hamiltonian, a nonlinear coupling term between the two–level system and the quantum oscillator. Such a coupling term has been discussed in the quantum optics literature under the name of dynamical Stark shift, a quantum version of the Bloch–Siegert shift [15]. Accordingly, we shall call the quantum Rabi model augmented by a nonlinear term of the kind discussed by Grimsmo and Parkins the quantum Rabi–Stark model.

Note, however, that in the usual dynamical Stark shift the corresponding nonlinear coupling strength is determined by the parameters of the underlying quantum Rabi model. In the scheme proposed by Grimsmo and Parkins [14] also the Stark coupling can be adjusted freely and independently.

Grimsmo and Parkins conjecture [16] that the Rabi–Stark model may undergo a superradiant transition in the deep strong coupling regime of the Rabi coupling when the Stark coupling strength becomes equal to the frequency of the cavity mode. The additional nonlinear term in the Hamiltonian, the Stark term, may therefore give rise to new physics. It will therefore be of considerable importance to thoroughly investigate the spectral properties of the quantum Rabi–Stark model.

The exact solvability of the model has been elegantly demonstrated in recent work by Maciejewski et al. [17, 18], using a Bargmann representation. The ensuing coupled set of differential equations were then solved by a technique involving Wronskian determinants in the general case and an analysis based on the Stokes phenomenon [19] for the special case when the Stark coupling becomes equal to the quantum oscillator frequency.

In this paper we take a different route to obtain the exact solution of the quantum Rabi–Stark model, adapting Braak’s method from 2011 [2] developed for the original quantum Rabi model. This alternative approach has the virtue of laying bare certain structural similarities between the two models, and highlights the importance of the underlying $\mathbb{Z}_{2}$ parity symmetry which is present also in the quantum Rabi–Stark model. In particular $-$ according to Braak’s criterion for quantum integrability [2] $-$ the retaining of the $\mathbb{Z}_{2}$ symmetry implies that also the Rabi-Stark model is integrable.

Almost all energy eigenvalues are determined by the zeros of two transcendental functions, obtained from a Frobenius analysis of the coupled singular differential equations which define the eigenvalue problem in the Bargmann representation. Provided that the model parameters are chosen so that these transcendental functions are reasonably well–behaved, this allows for numerical access to large portions of the spectrum. Fortunately, the parameter regimes where this property holds cover the most interesting cases for current experiments: the ultrastrong and opening deep strong regimes of the Rabi coupling.

There also exist, again like in the original quantum Rabi model, exceptional spectral points which do not correspond to zeros of these transcendental functions, but to points in parameter space where the singularities of the transcendental functions are lifted. As for the original quantum Rabi model, the exceptional solutions may define level crossings in the spectrum between energy levels of different parity. By increasing the magnitude of the Stark coupling we find that these level crossing points become less and less frequent. Instead there is a stronger tendency $-$ as compared to the original quantum Rabi model $-$ for neighboring levels to coalesce and eventually become two-fold degenerate. This surprising effect comes about from a “reshuffling” of energy levels caused by the added nonlinear Stark coupling, yielding a compressed spectrum which favors pairwise degenerate levels as the two–level system gets coupled to the quantum oscillator more strongly.

The layout of the paper is as follows: In the next section, section 2, we introduce the model, with reference to [14], and discuss some of its key properties. Section 3 contains the analytical solution of the model, leading up to the construction of the transcendental functions, the zeros (lifted singularities) of which determine the regular (two-fold degenerate exceptional) part of the exact spectrum (which becomes complete when adding also the non-degenerate exceptional part of the spectrum, as discussed in section 3). In section 4, the spectral structure in the ultrastrong and opening deep strong coupling regimes is extracted numerically from the exact solution, and the novel features $-$ as compared to that of the original quantum Rabi model $-$ are highlighted and discussed. Section 5, finally, contains a summary and outlook.

2 The quantum Rabi–Stark model

As we have expounded in the introduction, the quantum Rabi model describes the interaction between light and matter, next to the Jaynes–Cummings model, in the simplest possible way and is used as a basic model in many fields of physics [20].

The simplest generic experimental set–up to realize the quantum Rabi model, a cavity quantum electrodynamics (cavity QED) system, consists of a single atom put into a single–mode photon field which is enclosed by mirrors in a cavity. The frequency of the single–mode photon is chosen in such a way as to interact predominantly only with two levels of the atom [21]

In an experiment, there will inevitably be processes which lead to dissipative losses. In a cavity QED experiment, such processes include the dissipative loss of photons from the cavity (at rate $\kappa$) and the emission of the atom into other modes than the single cavity mode (at rate $\tau$). If such losses can be made small compared to the interaction strength between the single photon mode and the atom, described now as a two–level system, the experimental situation can be described by the quantum Rabi model Hamiltonian

[TABLE]

where $a^{\dagger}$ and $a$ are the creation and annihilation operators of the quantum oscillator mode with frequency $\omega$. The two–level atom is described by the Pauli matrices $\sigma^{x}$ and $\sigma^{z}$ with the splitting between the two levels given by $\Delta$. The interaction strength between the single photon mode and the two–level system is $g$ which we call the Rabi coupling to distinguish it from the Stark coupling which will be introduced below.

As already mentioned in the introduction, the Rabi model was originally introduced as the basis to understand nuclear magnetic resonance [3] and has since been applied to physical systems ranging from quantum optics to condensed matter physics, e.g. cavity and circuit quantum electrodynamics, quantum dots, trapped ions, and superconducting qubits. Moreover, it is used to describe nanoelectromechanical devices where the role of the photons is taken by phonons (see, for instance, [22] and [23]). These physical systems are also under investigation as candidates for the physical realization of quantum information processing.

Grimsmo and Parkins [14] propose an experimental arrangement where the two relevant levels of a 87Rb atom in the single–mode cavity is subjected to two auxiliary laser beams. Under conditions equivalent to the ones described above where losses can be neglected, Grimsmo and Parkins can describe their proposed experimental arrangement by an effective Hamiltonian

[TABLE]

where an additional term, $\gamma\sigma^{z}a^{\dagger}a$, appears compared to the original quantum Rabi Hamiltonian ${\mathcal{H}}_{\mathrm{Rabi}}$. This additional term models a nonlinear coupling between the two–level atom and the single–mode cavity oscillator. In the introduction, we gave an argument for naming this Hamilton and the corresponding model the quantum Rabi–Stark Hamiltonian and model, respectively, with the coupling constant $\gamma$, the Stark coupling.

The Hamiltonian (2) of the original quantum Rabi model is solvable by elementary means, employing the RWA (see for example [4] where also the classical and semi–classical versions of the Rabi model are discussed). The resulting model, the quantum Jaynes–Cummings model, emerges through the RWA by neglecting the terms $a^{\dagger}\sigma^{+}$ and $a\sigma^{-}$ in the Hamiltonian (2).

The Jaynes–Cummings model can also be investigated with an analogous nonlinear Stark term added. Interestingly, this variant of the Jaynes–Cummings model sheds light on the Bethe ansatz solution of the original Jaynes–Cummings model. The former can be solved by a standard algebraic Bethe ansatz procedure which allows to extract the algebraic solution of the latter in the limit when the Stark term vanishes [24, 1].

3 Exact solution of the Rabi-Stark model

In this section, we shall outline the exact solution of the quantum Rabi–Stark model represented by the Hamiltonian (4). In doing so, we shall generalize the method introduced by Braak in [2] for the solution of the original quantum Rabi model, described by the Hamiltonian (1), and especially highlight those aspects where the two models differ.

3.1 Bargmann space representation of the eigenvalue problem

It will prove advantageous to rewrite the quantum Rabi–Stark Hamiltonian (4) in the spin–Boson representation, achieved through a unitary rotation of the Hamiltonian by the operator $e^{i\pi\sigma^{y}/4}$. The Hamiltonian (4) then becomes

[TABLE]

In order to calculate the eigenvalues of this Hamiltonian exactly, we employ the Bargmann space representation [26] (for a recent summary, with a view on its application to the quantum Rabi model, of the properties of the Bargmann space representation, which is isomorphic to the space of square integrable functions ${\mathscr{L}}^{2}(\mathbb{R})$, see [27]). In the Bargmann space representation, the quantum oscillator creation operator is replaced by a complex variable $z$, i.e. $a^{\dagger}\rightarrow z$, and the quantum oscillator annihilation operator by the derivative with respect to the complex variable, i.e. $a\rightarrow d/dz$. The state vector $|\psi\rangle$ is represented in the Bargmann space representation by a wave function $u(z)$ depending on the complex variable $z$.

We briefly state the two requirements a function $u(z)$ needs to satisfy in order to be an admissible function of the Bargmann space $\mathscr{B}$, i.e. to be a physically allowed wave function. These requirements were carried over by Bargmann from the corresponding requirements which wave functions have to satisfy in the space of square integrable functions ${\mathscr{L}}^{2}(\mathbb{R})$. The first requirement is that the function must have a finite norm $\langle u|u\rangle<\infty$, where the scalar product is defined by

[TABLE]

and the second requirement that it be holomorphic everywhere in $\mathbb{C}$, i.e. be an entire function [26].

Measuring energy in units of the quantum oscillator frequency, i.e. formally putting $\omega=1$, the Rabi–Stark Hamiltonian (5) becomes in the Bargmann representation

[TABLE]

The canonical Fulton–Gouterman transformation [28]

[TABLE]

employing the parity operator ${\mathcal{T}}[u(z)]=u(-z)$, transforms the Hamiltonian (7) onto diagonal form

[TABLE]

with the Hamiltonians

[TABLE]

in the parity Hilbert spaces $\mathscr{H}_{\pm}$. The corresponding Schrödinger equations in the positive and negative parity sectors, respectively,

[TABLE]

become, written explicitly, non–local functional differential equations

[TABLE]

These two differential equations are converted into each other by the simultaneous replacements $\gamma\rightarrow-\gamma$ and $\Delta\rightarrow-\Delta$. It is therefore sufficient, and we shall do this in the following, to concentrate on one differential equation, here chosen as the one in the positive parity sector.

The reducibility of the Bargmann representation (7), into two blocks ${\mathcal{H}}_{\pm}$ with definite parities $\pm 1$, reflects that the Rabi-Stark Hamiltonian (3) is invariant under the $\mathbb{Z}_{2}$ parity transformation

[TABLE]

Hence, the eigenstates $|\psi\rangle$ can be labeled by the energy eigenvalue $E$ and the parity eigenvalue $p=\pm 1$,

[TABLE]

The $\mathbb{Z}_{2}$ parity symmetry is crucial for both, the exact solution of the model, and also its quantum integrability according to the quantum integrability criterion proposed by Braak [2].

Returning to (12), in order to deal with the non–locality of the differential equation for $\psi^{(+)}$, we define the two new functions (dropping the upper index $(+)$ for the time being)

[TABLE]

thus obtaining a set of two local differential equations. Note that this definition means that we now have two representations of the same function $\psi(z)$ which are to be determined from the two coupled local differential equations. With these definitions and rearranging terms, this set of two coupled local differential equations becomes explicitly

[TABLE]

Note that these two first–order complex differential equations are coupled in both, the unknown functions $\phi(z)$ and $\bar{\phi}(z)$ and their derivatives $d\phi(z)/dz$ and $d\bar{\phi}(z)/dz$. This is an important difference and complication compared to the original quantum Rabi model and is due to the nonlinear term proportional to the Stark coupling strength $\gamma$.

The two coupled first–order differential equations can be partially decoupled with respect to the coupling of the derivatives. In compact notation, we obtain the set of two first–order ordinary differential equations

[TABLE]

where we introduced the functions

[TABLE]

with $w=g/\sqrt{1-\gamma^{2}}$, and

[TABLE]

From these functions, especially (20), we observe that the differential equations are singular with regular singularities at $z=\pm w$ (see figure 1). Note that the regular singularities in the Rabi–Stark model depend on both, the Rabi coupling $g$ and the Stark coupling $\gamma$.

Furthermore, the equations have an irregular singularity at $z=\infty$ of s–rank $R(\infty)=2$ [29] which can be demonstrated by transforming the equations into second–order equations outside of a sufficiently large disk of radius $|z|=R$ which includes all singularities lying in a finite region of the complex plane. The s–rank $R(\infty)=2$ of the differential equations guarantees that the solutions have a finite norm asymptotically for $z\rightarrow\infty$ and are thus members of the Bargmann space [29, 30, 27].

3.2 Frobenius analysis of the singular differential equations

An indicial analysis [30] of the Frobenius ansatz around the regular singular points $z_{0}=\pm w$

[TABLE]

of the decoupled second–order differential equation for $\phi$ $-$ obtained from the coupled first-order differential equations in (18) and (19) $-$ reveals that there is one indicial exponent

[TABLE]

at each of the regular singular points $z_{0}=\pm w$.

The other indicial exponent is given by

[TABLE]

again at both regular singular points $z_{0}=\pm w$. The same indicial exponents are also obtained for the Frobenius ansatz

[TABLE]

from an indicial analysis of the second–order differential equation for $\bar{\phi}$, again at both regular singular points $z_{0}=\pm w$.

There is a subtle point to note about the indicial analysis. The limit $\gamma\rightarrow 0$ does not in general reproduce the indicial exponents of the differential equations for the original quantum Rabi model [27]. The reason for this is that the indical analysis requires a limit $z\rightarrow\pm w$ which cannot be interchanged with the limit $\gamma\rightarrow 0$.

It is important to stress that the indicial exponents determine whether the series solutions of the differential equations are also physically acceptable solutions, i.e. wave functions, belonging to the Bargmann space $\mathscr{B}$. If $r\in\mathbb{N}_{0}$, this is the case. However, solutions for generic $r$, i.e. for values of $r\notin\mathbb{N}_{0}$, although mathematically valid, are not members of the Bargmann space of physical wave functions.

For the quantum Rabi model where $\gamma=0$, the further analysis of our differential equations (18) and (19) can proceed directly [27] or after transforming them into second–order equations [31]. In the present case, the transformation to second–order differential equations generates further singularities not present in the first–order equations which make the analysis difficult. It is therefore preferable to directly solve the first–order equations as we shall do in the following.

For generic values of the parameters $\{\Delta,\gamma,g\}$ and the energy eigenvalue $E$, the indicial exponent $r_{1}$ will be a positive non–integer real number and, hence, the corresponding Frobenius solution, exhibiting a branch cut, will not be a member of the Bargmann space $\mathscr{B}$, i.e. will not be a physical solution. In these cases only the indicial exponents $r_{2}=0$ correspond to physical solutions $\phi(z)$ and $\bar{\phi}(z)$ belonging to the Bargmann space. The corresponding energy eigenvalues constitute the regular spectrum [32, 33, 34] of the Rabi–Stark model Hamiltonian.

However, for special combinations of the parameters $\{\Delta,\gamma,g\}$ and the energy eigenvalue $E$, the indicial exponent $r_{1}$ may become a non–negative integer. Such combinations give rise to the exceptional spectrum of the model, in close analogy with how exceptional spectra emerge in Jahn-Teller-like systems, first discussed by Judd [32].

In the following section 3.3, we concentrate our attention on the regular spectrum, while we shall discuss the exceptional spectrum in section 3.4.

3.3 Regular spectrum

Through the solution of the set of coupled differential equations (18) and (19) for the case $r_{2}=0$, we obtain the regular part of the spectrum. We focus on the singularity at $z_{0}=-w$ and introduce the new complex variable $y=z-z_{0}=z+w$ to perform a transformation of the functions $\phi(z)$ and $\bar{\phi}(z)$ according to

[TABLE]

which implies for the first derivatives

[TABLE]

such that the two first–order differential equations become

[TABLE]

with the constants $K_{2},\ldots,C_{0}$ depending on the parameters $\{\Delta,\gamma,g\}$ and the energy eigenvalue $E$:

[TABLE]

and [35],

[TABLE]

Writing $\rho(y)$ and $\bar{\rho}(y)$ as a power series

[TABLE]

where the expansion coefficients $\alpha_{n}$ and $\bar{\alpha}_{n}$ depend on the parameters $\{\Delta,\gamma,g\}$ and the energy eigenvalue $E$, we obtain a set of two coupled recursion relations for $n\geq 2$,

[TABLE]

The recursion relations for $n=0$ and $n=1$ can be obtained directly but also by the formal requirement that the expansion coefficients $\alpha_{n}$ and $\bar{\alpha}_{n}$ with index $n=-2$ and $n=-1$ vanish in the recursion relations (47) and (48).

For $n=0$, we obtain

[TABLE]

i.e. a set of two homogeneous algebraic equations for $\alpha_{0}$ and $\bar{\alpha}_{0}$. These algebraic equations have a non–trivial solution only if the coefficient determinant vanishes,

[TABLE]

This determinant indeed vanishes identically for all values of the parameters $\{\Delta,\gamma,g\}$ and all values of the energy eigenvalue $E$. The solutions of (49) and (50),

[TABLE]

can therefore be used as initial values for the coupled recursion relations (47) and (48).

With the procedure described above, we have now obtained the holomorphic solutions $\phi(z)$ and $\bar{\phi}(z)$ at the regular singular point $z_{0}=-w$ of the coupled set of the two first–order ordinary differential equations (18) and (19). These solutions are valid in a disk of convergence of radius $2w$ around the regular singular point $z_{0}=-w$ (see figure 1). They will, however, in general, i.e. for arbitrary values of the energy eigenvalue $E$ not be holomorphic at the other regular singular point, $z_{0}=w$, but will develop branch cuts at this singular point.

On the other hand, by a corresponding analysis we can find the holomorphic solutions $\phi(z)$ and $\bar{\phi}(z)$ to (18) and (19) which are valid in a disk of convergence of radius $2w$ around the regular singular point $z_{0}=w$. Again, these expansions, holomorphic at the regular singular point $z_{0}=w$, will in general not be holomorphic at the other regular singular point $z_{0}=-w$.

The symmetry of the differential equations (18) and (19) under reflection $z\rightarrow-z$ reveals that the two combinations, written in vector notation as $\left(\phi(z),\bar{\phi}(z)\right)^{T}$ and $\left(\bar{\phi}(-z),\phi(-z)\right)^{T}$, satisfy the set of differential equations (18) and (19). This property implies that, having obtained a holomorphic solution at one regular singularity through the procedure outlined above, say at $z_{0}=-w$, we also have one at the other regular singularity, i.e. at $z_{0}=w$. However, they represent one and the same function, as required in (15), only if the corresponding energy eigenvalue $E$ belongs to the discrete spectrum of the Hamiltonian (4). Then these solutions can serve as analytic continuations of each other. Together with the s–rank $R(\infty)=2$ for the irregular singularity at $z\rightarrow\infty$, this guarantees that we can find solutions of (18) and (19) which satisfy the requirements for physical solutions of the Bargmann space $\mathscr{B}$.

In practice, the coupled recursion relations can only be solved numerically. Assuming that we have obtained the expansion coefficients, at least to a sufficient degree of numerical accuracy, we can extract the energy eigenvalue $E$ from the solutions of the first–order differential equations, i.e. the wave functions in the Bargmann space representation. This is done by adapting the $G_{\pm}$ function formalism developed by Braak [2] for the quantum Rabi model to our purposes of the generalization of the Rabi model, the quantum Rabi–Stark model. Reintroducing the parity label $(\pm)$ for the wave functions $\phi$ and $\bar{\phi}$, we accordingly introduce the $G_{\pm}$ functions which are functions of the energy eigenvalue $E$, the parameters of the Hamiltonian $\{\Delta,\gamma,g\}$, measured in units of the quantum oscillator frequency $\omega$, and the complex variable $z$

[TABLE]

These functions must vanish for $E$ being an eigenvalue of the Hamiltonian (4), i.e. their zeros at, e.g. $z=0$, $G_{\pm}(E;0)=0$, determine the energy eigenvalues $E$ of the regular spectrum.

3.4 Exceptional spectrum

We have seen in the previous section that the zeros of the functions $G_{\pm}$ determine the energy eigenvalues of the regular spectrum of the quantum Rabi–Stark model.

However, the functions $G_{\pm}$ have poles at certain discrete values of the energy $E$. Thus, while almost all eigenvalues belong to the regular spectrum, in order to determine the complete spectrum, one has to investigate also the values of $E$ where at least one of the $G_{\pm}$ functions diverges. These values of $E$ cannot belong to the regular spectrum, as this is determined by the set of zeros of the $G_{\pm}$ functions.

Instead, these values appear as candidates for the exceptional eigenvalues, which, together with particular combinations of the model parameters $\{\Delta,\gamma,g\}$, turn the indicial exponent $r_{1}=(E+g^{2}+\gamma\Delta)/(1-\gamma^{2})=x_{\gamma}$ into a non–negative integer. Thus, in addition to the Frobenius solutions (23) and (26) corresponding to the indicial exponent $r_{2}=0$ which always belongs to the (physical) Bargmann space, now also the Frobenius solutions corresponding to an indicial exponent $r_{1}=x_{\gamma}\in\mathbb{N}_{0}$ in (23) and (26) become members of the Bargmann space $\mathscr{B}$.

Similar to the case of the original quantum Rabi model [27], we expect two possibilities for the exceptional spectrum. This expectation is borne out by our numerical exploration of our exact solution of the Rabi–Stark model which we report on in the next section 4.

4 Spectral structure

In this section, we report on our numerical procedure to extract the spectrum of the Rabi–Stark model and present our numerical findings.

4.1 Numerical procedure for the regular spectrum

Given the formal solution of the quantum Rabi–Stark model, as derived in section 3.3, the recipe to numerically extract the regular part of the energy spectrum can be summarized as follows:

In order to access the regular part of the spectrum in the positive parity sector for generic values of the model parameters $\{\Delta,\gamma,g\}$ (as before, always having set $\omega=1$), determine the expansion coefficients $\alpha_{n}$ and $\bar{\alpha}_{n}$ for $n=1,2,\ldots,N,$ from the recursion relations (47) and (48) with initial conditions as given in (49) and (50), supplemented by the definitions $\alpha_{-2}=\bar{\alpha}_{-2}=\alpha_{-1}=\bar{\alpha}_{-1}=0$; 2. 2.

Insert the expressions for $\alpha_{n}$ and $\bar{\alpha}_{n}$ from (i) into (23) and (26) (with $r=r_{2}=0$) via (45) and (46) as well as (27) and (28) and sum the first $N+1$ terms to obtain truncated series representations of $\phi^{(+)}(z)$ and $\bar{\phi}^{(+)}(z)$ (for book keeping purposes, now labeled as belonging to the positive parity sector); 3. 3.

Refer to (54) to construct the corresponding $G_{+}$ function; 4. 4.

Locate the zeros (a.k.a. energy eigenvalues) $E_{1},E_{2},\ldots$ of $G_{+}(\Delta,\gamma,g|E;0)$.

The regular spectrum of the negative parity sector is obtained by repeating the steps (i)-(iv) above, but with the replacements $\Delta\rightarrow-\Delta$ and $\gamma\rightarrow-\gamma$ (and with $\phi(z)$ and $\bar{\phi}(z)$ in (ii) now labeled as $\phi^{(-)}(z)$ and $\bar{\phi}^{(-)}(z)$ respectively, and with the energy eigenvalues obtained as the zeros of the corresponding function $G_{-}(-\Delta,-\gamma,g|E;0)$ in (54)).

As long as the $G_{\pm}$ functions are reasonably well-behaved (as they are, if one does not venture too far into the deep strong coupling regime $g>1$), the numerical root–finding can be carried out expeditiously, with stable results already for a truncation of the series in (23) and (26) to $N=12$ terms.

It is worth pointing out that the essential difference from the analogous protocol for obtaining the regular spectrum of the original quantum Rabi model [2] is that the expansion coefficients $\alpha_{n}$ and $\bar{\alpha}_{n}$ now have to be derived from two coupled recursion relations, (47) and (48). As discussed in sections 3.1 and 3.2, this reflects the fact that the differential equations (18) and (19) which determine the eigenfunctions $\phi^{(+)}(z)$ and $\bar{\phi}^{(+)}(z)$ (and $\phi^{(-)}(z)$ and $\bar{\phi}^{(-)}(z)$, respectively) of the quantum Rabi–Stark model have a more complex structure as compared to the case of the original quantum Rabi model.

4.2 Numerical procedure for the exceptional spectrum

Let us now turn to the exceptional part of the spectrum which can be obtained by the following route:

For fixed model parameters $\{\Delta,\gamma,g\}$, rewrite the recursion relations (47) and (48) in matrix form, i.e.

[TABLE]

where the vector $\mathbf{V}_{n}(E)$ is defined as

[TABLE]

with

[TABLE]

and where we have defined the determinant and then used (24),

[TABLE] 2. 2.

Find the zeros $E_{1},E_{2},...$ of the determinant $D_{n}(E)$. These zeros locate the common singularities of the functions $G_{+}$ and $G_{-}$ since they cause a divergence of the corresponding $\alpha_{n}$ and $\bar{\alpha}_{n}$ coefficients in (55). 3. 3.

For each $E_{j}$ thus identified, determine whether it is also a zero of the vector $\mathbf{V}_{n}(E)$ defined in (56). If this is the case, the singularity is lifted in both parity sectors (since the zeros of the vector $\mathbf{V}_{n}(E)$ are invariant under $\Delta\rightarrow-\Delta$ and $\gamma\rightarrow-\gamma$), and $E_{j}$ becomes a two–fold degenerate exceptional energy eigenvalue, determining a crossing between a positive and a negative parity energy level. 4. 4.

If $E_{j}$ is not a zero of $\mathbf{V}_{n}(E)$, the vanishing of $D_{n}(E_{j})$ still makes room for $E_{j}$ to become an exceptional solution. This is because the vanishing of $D_{n}(E_{j})$ corresponds to the indicial exponent $r_{1}=(E_{j}+g^{2}+\gamma\Delta)/(1-\gamma^{2})$ becoming a positive integer, i.e. $r_{1}=n\in\mathbb{N}$. As a consequence, and as explained at the end of section 3.4, $E_{j}$ becomes a nondegenerate exceptional energy eigenvalue in one of the parity sectors, corresponding to the Frobenius solution now turned into a new physical Bargmann wave function at this particular juncture of parameters which turns $r_{1}$ into a positive integer.

As we have seen from the discussion in this section, the eigenvalue spectrum consists of a continuous part, the regular spectrum, which is interrupted or punctured by isolated points of the exceptional spectrum. These latter punctures are characterized by zeros of the determinant $D_{n}(E)$ which cause divergences of the $G_{+}$ or $G_{-}$ function. The degenerate exceptional points occur simultaneously in both parity sectors and, thus, determine the level crossing points (as will be studied in an example in section 4.3).

As for the nondegenerate exceptional points, the continuity of the energy levels as functions of any of the model parameters $\Delta,\gamma$ or $g$ implies that also these points can only “fill out” some isolated punctures in the energy levels of either one or the other parity sector. Their locations are, thus, not immediately visible in a numerical plot of the spectrum, but must be calculated analytically. Since the nondegenerate exceptional points carry no particular significance for the interpretation of the spectrum, and also, since their detailed analytical determination is quite involved, we shall henceforth not elaborate upon these solutions.

For a discussion of the exceptional spectrum in the case of the original quantum Rabi model, see [36, 27]; for a discussion of the exceptional spectrum of a different generalization of the quantum Rabi model, obtained by adding an asymmetric term $\epsilon\sigma^{x}$, see [37, 38]. A detailed mathematical symmetry analysis using Lie algebra representations of $\mathfrak{sl}_{2}(\mathbb{R})$ is given for the spectrum of the original quantum Rabi model in [39] and of the asymmetric quantum Rabi model in [40].

4.3 Level crossings

It is instructive to witness in detail how a level crossing emerges by the lifting of a singularity in the $G_{\pm}$ functions. This is but one of the advantages of the $G$ function approach pioneered by Braak [2]: It allows for a compact encoding of the key features of the energy spectrum.

Figure 2 exhibits a case study, where $G_{+}$ ($G_{-}$) is shown in red (blue) versus $x=E+g^{2}$ in the interval $[-1,2]$ (with $E$ a running parameter which takes energy eigenvalues when $x$ becomes a zero of the corresponding $G_{\pm}$ function). The different panels correspond to different values of $g$, all with $\Delta=0.4$ and $\gamma=0.5$. In all panels, the two zeros closest to the singularity at $x=x_{s}\approx 0.55$ are marked with black circles. In the upper left panel a), the red (blue) zero is seen to be to the right (left) of $x_{s}$. As $g$ decreases, the two zeros creep closer to $x_{s}$, panel b), to eventually coalesce and annihilate at $x_{s}$ for a value of $g=g_{s}$ at which the singularity gets lifted, panel c). By further decreasing $g$, the zeros move away from $x_{s}$, which has now regained its role as a locus of a singularity in $G_{\pm}$. As seen in panel d), the zeros have traded their relative positions.

To sum up, the zeros of the $G_{+}$ and $G_{-}$ functions trade places as $g$ is varied across a common singularity of the two functions by lifting the singularity. As a consequence, a crossing between the positive and negative parity energy levels develops at $E_{\mathrm{cross}}=x_{s}-g_{s}^{2}$. We should add that while the loci of the $G_{\pm}$ singularities in the original quantum Rabi model appear at integer values of $x$, the loci for the quantum Rabi–Stark model now depend on the Stark coupling $\gamma$, with their presence being conditioned by the vanishing of the determinant (59).

4.4 Spectral structure of the quantum Rabi model

Before we present our numerical results for the spectrum of the quantum Rabi–Stark model, let us set the stage by recalling the main characteristics of the original quantum Rabi spectrum [2, 27]. A low–lying part of the spectrum with the levels as function of the Rabi coupling $g$ is depicted in figure 3, here with $g$ ranging continuously from the Jaynes-Cumming limit, $0<g\ll\Delta<1$, into the opening deep strong coupling regime, $1<g<1.6$, with the splitting of the two–level system $\Delta=0.4$.

The most notable feature in figure 3 is the absence of crossings between energy levels of the same parity. This allows for a unique labeling of the corresponding eigenstates, using the pair of quantum numbers $p$ and $n$, with $p=\pm 1$ denoting the eigenvalues of the parity operator ${\mathcal{P}}$, (13), and with $n=0,1,2,...$ indexing the progression of levels of increasing energy, identified as the zeros of $G_{\pm}$. According to the criterion proposed by Braak [2], the quantum Rabi model is quantum integrable because the eigenstates can be uniquely identified by using two quantum numbers ($p$ and $n$), equal to the number of degrees of freedom of the system (one two–dimensional degree of freedom characterizing the states of the two–level system, one infinite-dimensional degree of freedom for the quantum oscillator).

Since crossings, corresponding to the two–fold degenerate exceptional solutions (cf. figure 3), appear only between levels of different parity, one may find the resulting non–violation of the Wigner-von Neumann non-crossing rule [41] surprising: Quantum integrable systems are believed to violate the non-crossing rule [42, 43]. However, as expounded in [44], crossings between levels belonging to the same invariant subspace of a symmetry group (here: $\mathbb{Z}_{2}$ with positive and negative parity subspaces) are inevitable only for quantum integrable Hamiltonians where the number of local conserved quantities which depend linearly on the control parameter (here: the Rabi coupling $g$) is maximal, i.e. equal to the total number of constants of motion. Given that the quantum Rabi model does not belong to this class, there is no contradiction with the criterion suggested by Braak [2].

As seen in figure 3, levels of different parity with the same $n$ cross $n$ times before coalescing into near-degenerate levels for large Rabi–coupling $g$. This feature, present when $0<\Delta<1$, is also known from an analysis of the two-fold degenerate exceptional solutions, being of “Juddian” type [32] and accessible analytically [45]. In contrast, when $\Delta>1$, levels of opposite parities disentangle for small and intermediary values of $g$, with at most avoided level crossings remaining [46].

Given the analytical solution for the two-fold degenerate exceptional levels when $0<\Delta<1$ [45], one may further infer how two neighboring levels of different parity coalesce into a near-degenerate band [33]: For a given $n$ and for large $g\gg\Delta$, the levels will tend to the curve $E=n-g^{2}$, corresponding to one of the two–fold degenerate levels of the quantum Rabi model with $\Delta=0$. This behaviour is also easily read off from figure 3. It has a simple explanation: The two–fold degeneracy at $\Delta=0$ reflects the presence of a parity–flip symmetry: When $\Delta=0$, the quantum Rabi Hamiltonian (2) commutes with the parity-flip operator $\sigma_{x}$. This symmetry is destroyed when turning on $\Delta$, and thus, the two-fold degenerate levels get split. However, as the Rabi term $\sim g$ starts to dominate the level splitting $\Delta$ of the two–level system, there is a smooth crossover to the two–fold degenerate level with an emergent “approximate” parity-flip symmetry for very large $g$ (“approximate” in the sense that the residual terms which remain after commuting the Rabi Hamiltonian with the parity–flip operator $\sigma_{x}$ are small).

4.5 Spectral structure of the quantum Rabi–Stark model

With the description of the original quantum Rabi spectrum as a backdrop, we now turn to the quantum Rabi–Stark model, defined by the Hamiltonian (4). Its fourteen lowest energy levels for $\Delta=0.4$ are shown in figure 4 [47] as functions of $g$ in the interval $0\leq g\leq 1.6$ for $\gamma=0.5$ and $\Delta=0.4$.

Higher energy levels as well as spectra for larger values of $\Delta,\gamma$ or $g$ can also be extracted from the series representations of the $G_{\pm}$ functions in (54). However, the proliferation of singularities in the $G_{\pm}$ functions and the slowdown of the convergence of the series in (23) and (26) in these cases make the numerics more costly. We here confine our attention to the chosen parameter and energy regime in figure 4 [47].

Inspection of the spectrum in figure 4 [47] shows that crossings of energy levels of the same parity remain absent in the presence of the added nonlinear Stark coupling term. What may first appear as equal–parity level crossings (e.g. between the fifth and sixth blue curves close to $g=0.5$ in figure 4 [47]), at close scrutiny are revealed to be avoided level crossings, cf. figure 5.

Intriguingly, by tuning the splitting $\Delta$ of the two-level system, the avoided level crossings can be made progressively sharper, suggesting the possibility of a nonanalyticity for a critical value of $\Delta$, cf. figure 4 [47].

While it is tempting to speculate that this incipient nonanalyticity may be a precursor of an “excited state quantum phase transition” [48, 49, 50], this would be premature. In order to present support for such a transition, one must first and foremost establish a critical energy below which there is a symmetry breaking, with the critical energy accompanied by a singularity in the density of states. Let us note in passing that Puebla et al., using an effective Hamiltonian, have recently conjectured that such a transition may actually be present in the original quantum Rabi model [51] (see also [52]). Their approach was very recently generalized [53] for an anisotropic quantum Rabi model where the rotating and counterrotating parts of the Rabi coupling term acquire different coupling strengths: $g\left(\sigma^{+}+\sigma^{-}\right)\left(a+a^{\dagger}\right)\rightarrow g_{r}\left(\sigma^{+}a+\sigma^{-}a^{\dagger}\right)+g_{cr}\left(\sigma^{+}a^{\dagger}+\sigma^{-}a\right)$. We further note, again in passing, that this anisotropic model is also within the reach of the experimental proposal of Grimsmo and Parkins [14], that it admits an exact solution and that it can be used in a variety of physical situations [54, 34, 55], including a proposed realization of supersymmetry [56]. In view of these developments, it will be very interesting to examine the anisotropic generalization of the Rabi–Stark model.

As for the avoided level crossings in the quantum Rabi–Stark model, we expect that they rather reflect the model’s integrability (in the sense of Braak [2]): In order to uphold integrability as the original quantum Rabi levels get reshuffled by the nonlinear Stark term with strength $\gamma$, same–parity avoided level crossings appear in various parts of the spectrum. If same-parity level crossings had developed, this would have required the model to be ”superintegrable” [10] for nonzero values of $\gamma$, supporting an additional ”good quantum number” by which the energy levels could be uniquely labeled. This $-$ by itself quite unlikely $-$ scenario is made the more improbable by the presence of the avoided level crossings in figure 4 [47].

As is evident from figure 4 [47], the reshuffling of levels as $\gamma$ increases also leads to more densely spaced levels. This latter “compression” effect is anticipated from the trivial solution of the Rabi-Stark model with Rabi coupling $g=0$ which is solvable by elementary means since for this case the Hamiltonian is diagonal. Explicitly, the eigenvalues for $g=0$ are

[TABLE]

from which it can be seen that more and more levels accumulate at $E^{-}_{n}(g=0,\gamma\rightarrow 1)\rightarrow-\Delta$ as $\gamma\rightarrow 1$. In the same limit $\gamma\rightarrow 1$, the other levels become equally spaced, $E^{+}_{n}(g=0,\gamma\rightarrow 1)\rightarrow\Delta+2n$, starting from $+\Delta$. Please note that the label $\pm$ in (60) does not refer to parity. The parity of these levels is determined by the parity operator (13).

In this context: We have already dispelled a possible concern about nonviolations of the Wigner–von Neumann no–crossing rule. But what about the Berry–Tabor criterion [57] that an integrable model exhibits a Poissonian distribution of energy levels? Similar to the original quantum Rabi model, the levels for the quantum Rabi–Stark model when $\gamma\neq 0$ appear to be distributed fairly regularly (cf. figure 4 [47]) and not Poissonian. Thus both, the quantum Rabi and Rabi–Stark spectral distributions fail this test of integrability. However, it is important to be precise about the range of applicability of the Berry–Tabor criterion: It has been proved only in the semiclassical limit, and moreover assumes that the theory supports only continuous degrees of freedom [57]. None of this applies to the quantum Rabi or Rabi–Stark model.

Figure 4 [47] reveals that the level crossings of opposite–parity levels for the Rabi–Stark model ($\gamma\neq 0$) no longer follow the simple “braiding rule” of the Rabi model where two neighbouring levels with quantum number $n$ cross $n$ times. As a case in point, when $\gamma=0.5$ (figure 4 [47]), the first four pairs of opposite-parity levels cross at most two times before coalescing into a near-degenerate level. The implied reduction of two–fold degenerate exceptional solutions of the differential equations (18) and (19) when $\gamma\neq 0$, $-$ underlying the reduction of opposite–parity level crossings $-$ should have an explanation in terms of the $\gamma$-dependent loci of the singularities in the $G_{\pm}$ functions (cf. the discussion in section 4.2). However, to pinpoint the resulting structure of level crossings in the spectrum goes beyond the aim of this work. Indeed, a closer analytic examination of the $G_{\pm}$ functions remains a challenge for the future.

As already mentioned, there occurs a compression of the spacings of the energy levels as the coupling $\gamma$ of the nonlinear Stark term increases. This effect may facilitate $-$ but does not explain $-$ that for large $\gamma$, neighbouring levels with opposite parity coalesce into near–degenerate levels already for quite small values of $g$. For an example, see figure 6 [58], where the two lowest pairs of opposite–parity levels are shown as functions of $g$ when $\gamma=0.95$ (with the two lowest pairs of levels for $\gamma=0$ shown for comparison in the inset). Already for $g=0.6$ in the figure, the difference between the two lowest levels is $<0.00004$ and then rapidly decreases as $g$ increases further.

The rapid approach in figure 6 [58] to near–degeneracy of the energy levels for $\gamma=0.95$, as compared to the case with $\gamma=0$, in fact is surprising in view of the Hamiltonian ${\cal H}$ in (4). In order to make ${\cal H}$ approximately invariant under a parity flip with parity–flip operator $\sigma_{x}$ $-$ which guarantees near–degeneracy of neighbouring levels of opposite parity $-$ the Rabi term must now dominate both the term of the two–level system (with splitting $\Delta$) and the nonlinear Stark term (with coupling parameter $\gamma$). When $\gamma>\Delta$, as in figure 6 [58], one would expect that the approximate parity–flip symmetry, and thus, the concurrent near–degeneracy, would set in for values of $g$ larger than what is required when $\gamma=0$. However, as revealed by the same figure (with its inset), the opposite is the case!

These numerical observations suggest that the spectral compression dramatically enhances the effect of a reduced parity-flip symmetry breaking (so as to boost the early onset of near-degenerate energy levels), or else, that some hidden symmetry is at play for values of $g$ for which the parity–flip symmetry is still manifestly broken. While we find this latter alternative to be rather unlikely, we should alert the reader that there are indeed claims that already the original quantum Rabi model has a hidden symmetry, with implications for its dynamics [59]. In any event, a further analysis of the compression of the spectrum for increasing $\gamma$, with the concurrent rapid emergence of near–degeneracy, seems to be called for.

5 Discussion and summary

The generalized quantum Rabi model described by the Hamiltonian (4), the quantum Rabi–Stark model, is particularly interesting from the point of view that it offers a further tunable parameter, the Stark coupling $\gamma$, which can be used to investigate various regimes of the model which may be less accessible for the original quantum Rabi model where the Stark coupling vanishes. This is an especially intriguing aspect of the model since an experimental realization has been proposed with all the energy parameters $\{\omega,\Delta,\gamma,g\}$ freely and independently variable [14].

In the investigation reported here, the quantum Rabi–Stark model has been shown to be exactly solvable and also quantum integrable in the sense of quantum integrability introduced by Braak in his seminal work on the original quantum Rabi model [2]. Furthermore, we have obtained the exact analytical solution of the quantum Rabi–Stark model Hamiltonian adapting the methods devised for the original quantum Rabi model in [2]. In particular, we highlighted the differences created by the nonlinear Stark term $\gamma\sigma^{z}a^{\dagger}a$ in the generalized model. One of these differences concerns the reproduction of the results for the original quantum Rabi model from those of the generalized model. The naive limit $\gamma\rightarrow 0$ fails for the non–zero indicial exponents of the Frobenius solution of the generalized model. This observation will be crucial for a study of the exceptional points in the spectrum of the quantum Rabi–Stark model.

From the exact solution, we constructed functions $G_{\pm}(E;z)$ which can be used to numerically extract the regular part of the spectrum of the model. The exact solution also allows for a classification of the exceptional part of the spectrum which consists in its turn of a degenerate and a nondegenerate part. The spectrum and its properties, especially its dependence on the parameter $\gamma$ of the non–linear Stark term in the Hamiltonian has been the major aim of the investigation we have reported on in this paper. As detailed in the last section, the low-lying Rabi-Stark spectrum in the experimentally most relevant ultrastrong and opening deep strong regimes of the Rabi coupling exhibits two striking features absent from the original quantum Rabi spectrum: Distinctive avoided level crossings within each parity sector, and the onset of two-fold near-degenerate levels already in the ultrastong regime when $\gamma$ becomes sufficiently large. While the same-parity avoided level crossings most likely reflect the integrabiltiy of the model $-$ as bolstered by the underlying $\mathbb{Z}_{2}$ parity symmetry [2] $-$ the rapid onset of near-degeneracy remains more intriguing. To provide for its interpretation or explanation remains an open problem.

We thank Murray Batchelor, Daniel Braak, Mang Feng, and Michael Tomka for illuminating discussions, and Arne Grimsmo and Jonas Larson for valuable correspondence. Furthermore we thank Elinor Irish for pointing out an erroneous statement in an early version of the manuscript. This work was supported by STINT (Grant No. IG2011-2028) and the Swedish Research Council (Grant No. 621-2014-5972).

6 Corrigendum (published in J. Phys. A: Math. Theor.

56, 345302 (2023))

In our paper (pages 1-25; published in J. Phys. A: Math. Theor. 50 294004 (2017)), we solved exactly the quantum Rabi-Stark model, a generalization of the quantum Rabi model (for a recent elementary discussion of the quantum Rabi model, see [1]), and discussed its spectral structure. Doing this required a generalization of Braak’s method [2] developed for the exact solution of the quantum Rabi model. This method crucially uses the $\mathbb{Z}_{2}$ symmetry of the models. Working in the Bargmann space representation, the ensuing coupled first-order singular ordinary differential equations are solved by a modified power series ansatz, the so-called Frobenius method. The coefficients in this power series are then given, for the quantum Rabi-Stark model, by coupled three-term recurrence relations. Once these recurrence relations are solved numerically, the energy spectra are obtained as the zeros of the so-called $G$ functions.

The recurrence relations, equations (47) and (48), for the power series coefficients contain twelve coefficient functions, equations (33)-(44) in our original paper, which depend on the model parameters (i.e. the oscillator frequency ($\omega\equiv 1$), the qubit splitting ($\Delta$), the Rabi ($g$), and Stark ($\gamma$) coupling) and the energy $E$.

Two of these coefficients, $\bar{C}_{2}$ in equation (39) and $\bar{C}_{1}$ in equation (40) were miscalculated and should be replaced by

[TABLE]

where $w=g/\sqrt{1-\gamma^{2}}$.

Using the corrected coefficients (61) and (62), one obtains the energy spectra in figures 7 and 8 (replacing figures 4 and 6, respectively, in our original paper). Also shown in figures 7 and 8 are comparisons with energy spectra obtained by exact numerical diagonalization with a truncation number of $N=400$. The numerically exact diagonalization of the quantum Rabi-Stark model is completely independent of the procedure carried out in our original paper. However, the fit of the spectra obtained by these two independent methods is excellent.

While the energy spectra in figures 4, 5 and 6 of our paper are numerically incorrect and should not be used for extracting quantitative information, the general conclusions drawn in the paper remain unchanged and valid. This includes the scenario for the lifting of degeneracies in the $G$ functions illustrated in figure 2 of our original paper.

Our special thanks go to Lei Cong for performing the numerical calculations for this corrigendum, in particular the exact numerical diagonalization results. Furthermore, we thank Daniel Braak, Lei Cong, and Elinor Irish for many useful discussions.