Ultrastrongly dissipative quantum Rabi model

We discuss the equilibrium and out of equilibrium dynamics of cavity QED in presence of dissipation beyond the standard perturbative treatment of losses. Using the dynamical polaron \emph{ansatz} and Matrix Product State simulations, we discuss the case where both light-matter $g$-coupling and system-bath interaction are in the ultrastrong coupling regime. We provide a critical $g$ for the onset of Rabi oscillations. Besides, we demonstrate that the qubit is \emph{dressed} by the cavity and dissipation. That such dressing governs the dynamics and, thus, it can be measured. Finally, we sketch an implementation for our theoretical ideas within circuit QED technology.


I. INTRODUCTION
More than 80 years ago, Rabi studied the interaction of a two level system (TLS) with a classical electromagnetic field [1]. Jaynes and Cummings (JC) quantized this theory [2], focusing on the case of a single mode. This model is the quantum Rabi model ( = 1) In this Hamiltonian, ∆ and Ω are the bare TLS and cavity frequencies, while g denotes the light-matter interaction strength, see Fig. 1b. Considering the dissipation of the two-level system γ and of the cavity κ, we obtain several light-matter regimes [Cf. Fig. 1a]. When coupling outweighs dissipation g {γ, κ}, the qubit and the cavity field exchange excitations in a coherent way. Here (blue zone in Fig. 1a) we distinguish two regimes. If g/∆ 0.1, we are in the strong coupling (SC) regime and light-matter interaction can be simplified into g(σ + a + h.c.) using the Rotating Wave Approximation (RWA). However, if g/∆ 0.1, the RWA fails and the full interaction, i.e. the counter rotating terms -g(σ + a † + h.c.)-are needed. This is the ultrastrong coupling regime (USC) [3]. By analogy, when dissipation dominates, we distinguish between the weak (W) and weak ultrastrong coupling (WUSC), depending on whether we can apply the Markovian approximation or not to describe dissipation (red area in Fig. 1a) .
The goal of this work is to derive a mathematical treatment of the cavity-QED model that provides quantitatively or qualitatively accurate solutions in all coupling and dissipation regimes -WC, SC, WUSC and USC-. There are many ways to solve the cavity-QED model that apply to subsets of these regimes. In absence of dissipation, equation (1) admits an analytical solution [4] and can be solved in the computer for any g-value. If losses are taken into account, they are typically discussed using Markovian master equations [5,6]. They are perturbative in the system-bath interaction [7,8]. Going beyond this perturbative treatment is tricky [9]. Renormalization, path integral expansions or numerical techniques are required [10,11]. In contrast to this zoo of solutions, we will offer a unique method of broad utility with the only restriction that qubit dissipation remains below the quantum phase transition into the strongly correlated regime. The resulting method will be useful in studying all the quantum technologies that are developed around the JC model -single photon emitters, quantum a)

b)
FIG. 1. Cavity QED phase diagram and setup sketch. a) Blue region marks where Rabi oscillations occur. We distinguish between strong (SC) and ultrastrong (USC) coupling regimes. In the red region, losses are big enough and the TLS decays in an overdamped way. Here, we distinguish between weak (WC) and weak ultrastrong coupling regime (WUSC). In the figure α characterizes the TLS losses [Cf. Eq. (3a)]. In b) we draw a cavity QED skectch with the main parameters indicated. computers, spin squeezing [12]-, as well as experiments that exploit the huge dipole moments of superconducting qubits in the ultrastrong coupling regime [13][14][15].
Our method builds on the polaron Hamiltonian [16,17] to develop an effective model that can be analytically or numerically solved. Similar to the Ohmic spin-boson theory [18][19][20], we predict non-Markovian renormalization of the qubit splitting ∆ due to the coupling with the bath, either directly γ or via the cavity. We can also solve the qubit-cavity dynamics, from overdamped decay in the limit of WC or WUSC dissipation, to coherent scenarios that extend well inside the USC both in losses and light-matter coupling, see Fig. 1.
The outline of this work is as follows. In the next section II, we summarize the model, the polaron transformation and the different ways of solving the equilibrium and out of equilibrium dynamics. In section III we announce our results. We discuss the ground state properties of the model, the onset of Rabi oscillations and the noise spectrum. Finally, we give some conclussions and a possible implementation in IV. Several technical details are sent to the appendices.

A. Combined dissipation channels
We consider a qubit and a cavity interacting with each other and coupled to independent baths. The model in the system-bath formalism [21,22] is [cf. Eq, (1)], The qubit and cavity baths have independent modes ω k,i (i = 1, 2), with bosonic quadratures X k,i ≡ b † k,i + b k,i . Both noise channels can be described using spectral density functions, J i (ω) = 2π k c 2 k,i δ(ω − ω k,i ). In this paper we are considering Ohmic noise spectrum for both the cavity and the spin, i.e. J i (ω) ∼ ω. In the Markovian limit, the dissipation strength is determined by the spontaneous emission rates of the qubit (γ) and the cavity κ because [Cf. Fig. 1b]: with α (α cav ) dimensionless parameters characterizing the dissipation strenght for the TLS (cavity). The last term in equation (2) deserves some discussion. This regularization of the bosonic modes arises from a cavity-bath coupling of the form ∼ (a+a † −Φ) 2 , where a+ a † is the cavity quadrature and Φ is the electromagnetic field injected by the bath . This type of coupling -which is very natural in superconducting circuits-ensures that the total energy is bounded from below and leads to the quadratic correction of the bosonic modes. Importantly, the quadratic correction ensures that the resonance of the cavity stays at Ω, irrespective of the dissipation strenght α cav . Note also that we do not find a similar term in the qubit-bath coupling because of saturation: (σ x ) 2 = 1. Further parameter renormalization is associated to quantum many body effects between the bath and the cavity QED system [9,23].
The cavity mode in (1) can be diagonalized together with its environment. In doing so, Hamiltonian (2) is rewritten as a spin-boson model [18] for a two level system coupled to two baths, one of which contains the cavity mode, By joining the cavity modes and the qubit bath, we arrive at the total spectral density J(ω) = 2π k c 2 k δ(ω −ω k ): The second term characterizes the bath containing the cavity mode and is peaked around the cavity frequency Ω [24]. The first term accounts for the intrinsic qubit Ohmic environment. We will correct this spectral function by introducing a hard cut-off, ω c [19]. Finally, note that models, (2) and (4) are completely equivalent. When writing the model as (4), Rabi oscillations can be understood as non-Markovian decaying oscillations comming from the peaked spectral density. Details on the cavity-bath diagonalization and the effective spectral density are given in App. B.

B. Effective RWA models
It has been recently shown that the low energy spectrum of a spin-boson model (5) can be very well approximated by an effective, excitation number conserving Hamiltonian derived from a polaron transformation [16,17]. The basic idea is to construct a unitary transformation that disentangles the TLS from the bath and choosing the displacements f k with the Silbey-Harris prescription that the ground state of H p = U † p HU p be as close as possible to |0 ⊗ |0 , the ground state of the uncoupled TLS |0 and of the bath |0 . Minimization yields the self-consistent relation and the effective Hamiltonian H p is well approximated, within the single-excitation sector, by [25] This model has two important features. First, our TLS appears with the renormalized frequency ∆ r , that determines the dynamics of the qubit in the bath [cf. Sect. III B]. Second, the model develops a conserved quantity [H p , σ z + k b † k b k ] = 0 and becomes tractable with the same techniques as RWA models.

C. Estimating ∆r
The qubit renormalized frequency, ∆ r [Cf. Eqs. (7) and (8)] admits analytical solutions in the continuum limit, where we can write using the UV cutoff ω c . The resulting expression is not tractable, and requires the numerical solution of a transcendental equation for ∆ r in (9). One extreme limit of this equation appears when the qubit decouples from the cavity (g = 0). In that case Eq. (9) can be solved, resulting in the Ohmic spin boson model. Then, ∆ r = ∆(∆/ω c ) α/(1−α) , and the localizationdelocalization transition at α = 1 [18]. When we depart from this limit g = 0, the second summand in (5), decrease the onset of the localization transition to lower values of α.
The predictions of the polaron ansatz can be compared with the adiabatic renormalization (ARG) [18] in the continuum limit. Applied to the spin-boson model, the ARG predicts a qubit frequency at Comparing (9) and (10) we see a difference in the lower limit of the integral. The reason is that the RG flow stops at ω c ∼ ∆ r . If g = 0, both methods yield the same result except corrections of the order of O(∆ r /ω c ). Therefore, they predict the same ∆ r . However, if g = 0 ARG is different from polaron (and becomes less accurate). Our interpretation is that the cavity Ω behaves as an effective cuttoff and for ∆ ∼ Ω the ARG should fail. More generally, we also expect the ARG to fail whenever cavity losses and/or light-matter coupling dominate over the TLS intrinsic dissipation.

D. Modified Wigner-Weisskopf
We have just discussed that (8) conserves the number of excitations. We can solve the dynamics within a single-excitation subspaceà la Wigner-Weisskop. There, the dynamics is fully determined by the wavefunction: Using (8) and (11), the coefficients {ψ, ψ k } satisfy the set of coupled linear equations: From these coefficients we may derive, for instance, the excitation probability of the two-level system In the regime g/∆, α, α cav 1, we can solve for the qubit amplitude ψ applying the Markov approximation on the qubit losses and replacing the second summand in (5) with a Lorentzian centered on the cavity resonance Ω. Then an analytical solution is possible as it is fully developed in our App. C. In the ultrastrong both for losses and light-matter coupling,, analytical advances are possible in the calculation of the qubit noise spectrum S(ω) as we explain in Sect. III C. However, for the time evolution, in general, an analytical solution is no longer possible. Then, we approximate the environment using a finite number of modes, N , as explained in App. A. In doing so, we can solve the set of O(N ) ordinary differential (12a) and (12b) numerically, e.g. using Lanczos, Runge-Kutta or any other available method.

E. Matrix Product States
In order to confirm the predictions of H p , we also run numerical simulations on unapproximated model H P = U † p HU p using Matrix Product State ansatz. While working with H p significantly decreases the amount of entanglement in the MPS simulation, we introduce another optimization and express the H p as a tight-binding model [26]. The simulated model reads The new collective modes c i are constructed from the original ones b i using a Lanczos recursion [26] that also produces the real numbers α i and β i . This new model has the advantage that it can be simulated using both Arnoldi and Trotter-type MPS methods [27].

III. CAVITY-QED BEYOND MARKOVIAN REGIME
A. Ground state We will now show how to apply the previous formalism to study the static and dynamic properties of the cavity-QED setup in all regimes. We begin with the nature and properties of the ground state.
The construction of the polaron Hamiltonian provides a zeroth-order approximation to the ground state, U p |0, 0 , which predicts that the qubit has some probability to be excited. This is consistent with earlier findings in lossless cavities [7,28] and in the spin-boson model without cavity [18], but our new treatment allows us to interpolate between both limits. The equilibrium zmagnetization is proportional to the qubit renormalized frequency, which as we saw before, can be computed from the displacement f k (7). Let us now compare the estimates of σ z and ∆ r , from the adiabatic renormalization group, (10), the polaron method, Eq. (9), and an exact solution of H p with MPS discussed in Sects. II C and II E. Fig. 2 summarizes the ground state properties for different values of the dissipation and coupling strength. Those simulations have been performed with N = 256 modes for the cavity bath, and a similar amount for the qubit bath, ensuring numerical convergence to a quasi-continuum limit. In panel 2a we plot the qubit renormalization ∆ r as a function of the TLS dissipation [α, Cf. Eq. (3a)], for bare parameters Ω = ∆ = 1, an USC coupling strength g = 0.2 and a low cavity spontaneous emission κ = πα cav Ω = π0.01. We compare three methods: polaron, ARG and MPS simulations, as explained in Sect. II. The dependence of ∆ r on α resembles the pure Ohmic spin-boson model. As we have anticipated and explained in Sect. (II C) ARG is not accurate for small TLS intrinsic noise strenght α. The inset of Fig. 2 also shows that the qubit-cavity coupling lowers even further the qubit frequency ∆ r , due to the friction induced by the additional bosonic modes from the cavity and its bath.. Figure 2b probes the ground state for stronger cavity dissipation, entering the WUSC regime. As seen on the main panel, for relatively high cavity losses (κ = πα cav Ω = π0.8) the trend of ∆ r is qualitatively similar to the uncoupled case g = 0. However, if we increase κ (and g) enough both the polaron and ARG models predict a sharp transition, leading to localized solutions ∆ r = 0. Since the MPS are numerically exact simulations and do not exhibit such a transition, we conclude that this is an artifact of the polaron method that constraints its applicability to large values of dissipation and light-matter coupling, g ≥ 0.6 and κ ≥ 1. The remaining of this work will stay well within this regime, in which simulations verify well against MPS.

B. Non perturbative Rabi oscillations
In this section we study the dynamics of the cavity-QED setup, solving numerically and analytically the qubit excitation probability P e (13) with the polaron methods detailed in Sect. II D and in App. C. The first result is the evidence of coherent light-matter (Rabi) oscillations that (i) are resonant around the qubit renormalized frequency ∆ r and (ii) dampen exponentially with a modified spontaneous emission rate determined by the joint spectral function (3a) that we introduced in Sect. II. We find that Rabi oscillations start  approximately at the boundary Above this critical value, the TLS and the cavity exchange excitations coherently; below this boundary, the qubit exhibits overdamped exponential decay without oscillations. It is remarkable that this boundary is formally the same as the one in the RWA and Markovian approximation (g/Ω, α, α cav 1 regime) [29], but extends to the USC and WUSC regimes, which do not admit a perturbative treatment.
In Fig. 3 we show the qubit dynamics for varying g. In 3a we see the appearance of oscillations for g |γ r − κ|/4 ∼ = 0.1. These plots confirm that the TLS dynamical frequency is given by ∆ r ; since this resonance changes with g, we explore on-and off-resonant oscillations as we increase the coupling strength, for fixed Ω. With the parameters used, Ω = ∆ r = 0.68, which is reached at g = 0.3. Fig. 3b shows the qubit dynamics at precised couplings. This includes the resonant case g = 0.3, which shows resonant-like Rabi oscillations, in the middle plot. For lower coupling g = 0.05, losses dominate and the TLS dynamics is overdamped -i.e. exponential decay in the polaron frame-. Since for α = 0.1, ∆ r = ∆, this is an example of WUSC dynamics. The right-hand plot shows a USC coupling dynamics, with g = 0.6, where ∆ r ∼ = 0.4 and the dynamics are non-resonant Rabi-like oscillations.
We have found an analytical approximation that reproduces and explains the TLS-cavity dynamics and the onset of the Rabi oscillations. We take Eq. (5), remove the term that is O(f 2 ) and modify the rest, replacing the effective displacements f k with the original couplings f k → c k . The solution of these simplified equations is formally identical to the one for g/∆, α, α cav 1 (Markov and Lorentzian approximations) but with a renormalized frequency ∆ → ∆ r . We denote it P e (t). Besides, we need to impose that the time converges in the t → ∞ to the correct equilibrium solution given by P eq e = (1−∆ r /∆)/2 [Cf. (13) and (15)]. Notice, that our numerical simula-tions verified thermalization, marked as dotted lines in Fig. 3. To have the correct stationary limit, we use the simplest interpolation for our analytical estimation, P app e (t) ∼ = (1 − P eq e ) P e (t) + P eq e .
In figure, 3b) we show how such a approximation holds for relatively high g (well inside the ultrastrong coupling regime), justifying equation (17). At strong coupling, g = 0.6, the simple approximation P app e stops working (not shown), because we have neglected the O(f 2 ) terms. In any case, the evolution still reflects detuned Rabi oscillations, converging to the expected limit: P e (t → ∞) → 1 2 (∆ r /∆ − 1). We further verified the location of the critical value (17), analyzing numerically the transition from an overdamped dynamics, to the appearance of the first oscillations. For that we compute the maximum of the time We plot the maximum derivative of the qubit excitation probability, max dPe/dt as a function of g for different noise strenght (solid lines). This derivative is negative or zero in the overerdamped regime. Vertical dashed lines mark the predicition for the critical g given by (17). Equal colors mean equal parameters.
derivative of P e on the time interval [0, T ] with T sufficiently large. We denote this quantity as max dP e /dt. If the TLS is overdamped, then dP e /dt < 0 always and the max dP e /dt = 0. However, if some oscillations occur the derivative is sometimes positive. This is represented in figure 4 and compared with the bound (17), exhibiting a very good agreement. Two comments are in order. First, to generate Fig. 4 we have chosen to be approximately at resonance at the critical value of g. Second, since ∆ r can go to zero faster that linearly with α, γ r approaches to zero by increasing α. In the figure, we indeed see that α r (α = 0.2) ∼ = α r (α = 0.3) [Cf. Fig. 4]. Summing up, our simulations justify the use of quantum optics approximations in this non perturbative regime. We find that Rabi oscillations extend qualitatively into a regime where light-matter interaction and dissipation are both non-perturbative. In particular, we have studied a novel regime -denoted Weak Ultrastrong coupling regime (WUSC)-in which g is big enough that we cannot make the RWA, but at same time losses are also large and prevent coherent exchange between light and matter degrees of freedom [cf. Fig. 1].

C. Qubit noise spectrum
The TLS emission spectrum S(ω) is a very useful experimental tool that provides information about the coupling g and TLS line-width. S(ω) is typically computed using the input-output formalism [30]. In this framework, output and input fields are related to the TLS state via a out = a out − i √ Γ TLS X − (t), where Γ TLS is the emission rate into the transmission line in wich the output signal is collected and X − is the negative frequency component of the qubit-TL coupling operators (σ x in our case) [31,32]. The emission spectrum is defined as where R(ω) and Γ(ω) are the real and imaginary part of the self energy of the qubit. The former gives the position of the eigenvalues and the latter the line width. Their explicit expressions ares and K (K ) the real (imaginary) part of We notice that the renormalized frequency, ∆ r , is again explicit. If α, α cav , g/∆ 1, the linewidth reduces to recovering the standard results in cavity QED (using e.g. master equations for dealing with the bath [33]). The expressions (D14) and (D15) evidence important corrections in the response profile, with the most evident fact of an asymmetry between peaks at the dressed resonance Ω = ∆ r [Cf. Fig. 5]. This is a signature of the modified coupling constants c k → f k in the non perturbative polaron Hamiltonian. From Eq. (7) the effective coupling coefficients f k are smaller than c k for ω k < ∆ r and bigger in the region of the spectrum, leading to the asymetric profiles. A similar pheonomenon has been identified in the USC regime, which accounts for large g but only for weak dissipation [34].
In the previous section we tested a simple approximation to the TLS dynamics (13). This consisted in taking the RWA, weak/strong coupling solutions -i.e. α, α cav , g/∆ 1-, but replacing the qubit resonance ∆ → ∆ r , and correcting for the and correct the equilibrium state. If we use this approximation to estimate S(ω) we find that it more or less accounts for the linewidths. However, our qualitative method fails to reproduce the asymmetry of the peaks [cf. Fig. (5)], and also has a minor error in the peak location, due to the the Bloch-Siegert and dissipation-induced shifts.

D. Circuit QED implementation
The model that we have discussed in this work admits a straightforward realization using superconducting circuits. The three elements that we need are (i) a qubit that is ultrastrongly coupled to the cavity [13,14], (ii) a possibly ultrastrong coupling between the same qubit and some external environment, such as a transmission line [20], and (iii) a strong or ultrastrong coupling between the superconducting cavity and its own bath, a regime already achieved in [35].
All these three elements admit full and independent tuneability, probing arbitrary values of g, κ and γ. First of all, the coupling of the qubit to the photons -g and γ in our model-can be tuned using SQUIDs that can be embedded in the design of the qubit itself [19,36], as demonstrated in Ref. [20] for a flux qubit in an open transmission line. Moreover, the coupling between the cavity and its own bath κ, can also be adjusted using in-line dc SQUIDs. This has been demonstrated in the lab, and used for photon trapping and release [37,38] -applications that are much more demanding than the simple, stationary tuning of the parameter κ in our model.
Assuming a superconducting circuit implementations, we can probe the physics of the combined environments in various ways. We have studied the spectral function J(ω), which can be reconstructed from the spontaneous decay of an excited qubit. This requires a protocol in which (i) the cavity is decoupled, κ → 0, (ii) the qubit is excited, (iii) all couplings are switched on for a brief period of time, and (iv) the excited population of the qubit is measured in a non-destructive way [39]. Alternatively, it is possible to relate the spectral function to the qubit spectroscopy, by studying the low-power transmission spectrum of the cavity and relating the total lineshape to the spectral function, using the theory from Ref. [25].

IV. SUMMARY
We have studied a cavity QED model beyond the standard perturbative treatment of losses, using numerical and analytical techniques that apply both in and out of equilibrium. Our study builds on the polaron Hamiltonian [25] and Matrix-Product State simulations. In the former case, we have shown that techniques to solve the RWA and weak noise regime (as Wigner-Weisskopf or S(ω) calculation) can be extended to work with a variety of regimes -USC, weak, strong coupling, etc-.
As concrete applications, we have discussed with detail the case of a two-level system that couples ultrastrongly to both the cavity and the bath -i.e. both g and γ are comparable to the qubit and cavity resonances. Using our techniques, we prove that strong dissipation renormalizes the qubit frequency, leading to a new resonance Ω = ∆ r and changing its decay rate γ r . Our simulations show that this renormalized decay rate can be used to define the onset of Rabi oscillations (g ∼ = |γ r − κ|/4), in a formula that extends beyond RWA for all range of parameters. This suggests a new regime where light and matter are ultrastrongly coupled but losses are large enough to suppress Rabi oscillations. We call this regime the weak ulstrastrong coupling regime (WUSC) [Cf. Fig. (1)].
This work has different possible continuations. On the experimental side, we have shown that all regimes and physics shown in this work can be probed using state-of-the-art circuit QED technology. On the theory side, our numerical methods open the door to extend those experiments to study very challenging cavity QED phenomena, such as transmission-reflection experiments, Dicke physics or nonlinear optics, in the USC and WUSC regimes -enabling ultrasfast, broadband photon sources, opening access to stronger nonlinearities and facilitating the study of non-Markovian open quantum systems, among other new phenomena to be explored. so that, with ∆ω k is the the frequency interval around ω k . When a frequency is degenerate ω k ω k , we have to add up all contributions coming from the different couplings.
A transmission line is a model for an Ohmic bath. Its discrete version is a set of coupled harmonic oscillators The normal modes are known. In the chiral k ≥ 0 case, where the momentum spacing relates to ∆x The dispersion relation is then simply obtained from the speed of light (v = 1) In our model for the TLS ∆ω k = v∆k = ∆k, so that using (A5) and (A3) wich are the coupling used in our numerical simulations. In figure 6 we compare the continuum spectral density and this discretization. The agreement is clear.

The cavity-bath case
The cavity-bath model is the positive defined quadratic Hamiltonian: (A12) This is nothing but the Caldeira-Legget model of dissipation [21]. The spectral density of the cavity-bath is given by: operators and in the spin boson in terms of annihilationcreation operators.
The Caldeira-Legget model can be rewritten as: This model can be further diagonalized using a unitary transformation U and eigenvalues Ω 2 such that We quantized the model as usual: obtaining an expansion for the original cavity mode: In doing so, the qubit-cavity coupling can be written as the spin-boson coupling: (B7) The orthogonal transformation (A16) fulfills the normalization condition: In the second equality we have used (B2). Now, we notice that [Cf. Eq. (B4)] Using Eqs. (B7), (B8) and (B9) we arrive to:

Effective spectral density
We can rewrite the spin-boson coupling (A19) using (B10): Therefore, the effective spectral density for the spinboson reads [Cf. Eq. (A2)]: Using now (B11) for any well behaved function f (ω), we have that Therefore: So far, everything was general (we have not specified the spectral density for the cavity bath). We particularize to an Ohmic spectral density, see (A13) and (3b): . Now, we can compute A defined in (B5) Inserting the last result in the definition of g −1 , taking the imaginary part and using (B14) we get J eff (ω) = 4g 2 κω (Ω 2 − ω 2 ) 2 + (κω) 2 . (B17) which is nothing but the peaked spectral density discussed in the main text and that has been used to test our bath discretization [Cf. Fig. (7)].