Cutoff-free Circuit Quantum Electrodynamics
Moein Malekakhlagh Alexandru Petrescu Hakan E. Türeci





Department of Electrical Engineering, Princeton University, Princeton, New Jersey, 08544


Any quantum-confined electronic system coupled to the electromagnetic continuum is subject to radiative decay and renormalization of its energy levels. When coupled to a cavity, these quantities can be strongly modified with respect to their values in vacuum. Generally, this modification can be accurately captured by including only the closest resonant mode of the cavity. In the circuit quantum electrodynamics architecture, it is however found that the radiative decay rates are strongly influenced by far off-resonant modes. A multimode calculation accounting for the infinite set of cavity modes leads to divergences unless a cutoff is imposed. It has so far not been identified what the source of divergence is. We show here that unless gauge invariance is respected, any attempt at the calculation of circuit QED quantities is bound to diverge. We then present a theoretical approach to the calculation of a finite spontaneous emission rate and the Lamb shift that is free of cutoff.

Cutoff-free Circuit Quantum Electrodynamics
    Hakan E. Türeci
    December 30, 2023
===========================================



Introduction. An atom-like degree of freedom coupled to continuum of electromagnetic (EM) modes spontaneously decays. When the atom is confined in a resonator, the emission rate can be modified compared with its value in free space, depending on the EM local density of states at the atomic position <cit.>, which is called the Purcell effect <cit.>. An accompanying effect is the Lamb shift, a radiative level shift first observed in the microwave spectroscopy of the hydrogen ^2P_1/2- ^2S_1/2 transition <cit.>. These quantities have been experimentally accurately characterized for superconducting Josephson junction (JJ) based qubits coupled to coplanar transmission lines <cit.> and three-dimensional resonators <cit.>. In the dispersive regime where a qubit with transition frequency ω_j is far-detuned from the nearest resonant cavity mode (frequency ν_r, loss κ_r), single mode expressions exist for the Purcell decay rate, γ_P = (g/δ)^2 κ_r and the Lamb shift, Δ_L = g^2 / δ. Here g denotes the coupling between the qubit and the cavity mode and δ=ω_j - ν_r denotes their detuning <cit.>. However, for large couplings accessible in circuit QED, the single mode approximation is often inaccurate <cit.>. In addition, due to particular boundary conditions imposed by the capacitive coupling of a resonator to external waveguides, the qubit relaxation time is limited by the EM modes that are far-detuned from the qubit frequency <cit.>. Similarly the measured Lamb shift in the dispersive regime can only be accurately fit with an extended Jaynes-Cummings (JC) model including several modes and qubit levels <cit.>. The Purcell rate has been generalized to account for all modes




    γ_P=∑_n(g_n/δ_n)^2  κ_n,

where g_n and δ_n =ω_j - ν_n are coupling to and detuning from resonator mode n with frequency ν_n and decay rate κ_n. Expression (<ref>) is divergent without imposing a high-frequency cutoff <cit.>. Divergences appear as well in the Lamb shift and other vacuum-induced phenomena, e.g. photon-mediated qubit-qubit interactions <cit.>. These divergences are neither specific to the dispersive limit nor to the calculational scheme used to compute QED quantities. This issue is well-known for the Lamb shift <cit.>, but less noted for the spontaneous emission rate. Indeed, free space spontaneous emission rate diverges as well, as we show in <cit.>. The finite result by Wigner and Weisskopf <cit.> is due to Markov approximation which filters out the ultraviolet divergence. Recent generalizations of the Wigner-Weisskopf approach impose an artificial cut-off to obtain a finite result <cit.>. So far, no satisfactory theoretical explanation has been given for these divergences. Here we address this issue within the framework of circuit quantum electrodynamics <cit.> (QED) and show that finite expressions can be obtained when gauge invariance is respected. We focus here on a superconducting artificial atom coupled to an open transmission-line resonator, but our results should be valid for other types of one-dimensional open EM environments as well. 



Gauge invariance in circuit QED. The role of gauge invariance in accounting for light-matter interaction has been a vexing question since the beginnings of QED (see Ref. <cit.>, and references therein). Hence, we first discuss gauge invariance in  superconducting electrical circuits, and its impact on QED observables.  

We consider a weakly nonlinear charge qubit (e.g. transmon <cit.><cit.>) capacitively coupled to a transmission-line resonator that in turn is coupled at both ends to semi-infinite waveguides (Fig:cQED-open). We assign flux variables to nodes, Φ_n (t) = ∫^t dτ  V_n(τ), with V_n(t) being the instantaneous voltage at node n with respect to the ground node <cit.>. Fixing the ground amounts to a particular gauge choice <cit.>. For the connection geometry in Fig <ref>, the light-matter interaction derives from the energy on the coupling capacitor in the dipole approximation, T_int = 1/2 C_g [Φ̇(x_0) - Φ̇_j]^2 <cit.>, with x_0 the qubit position. If from the three terms in its expansion, T_EM =  1/2 C_g  Φ̇(x_0)^2, T_ EM-JJ = - C_g Φ̇(x_0) ·Φ̇_j and T_JJ = 1/2 C_g  Φ̇_j^2, only the direct interaction T_ EM-JJ is kept, a multimode JC model in terms of circuit parameters can be derived [See also the Supplementary Material for a brief derivation of the Heisenberg-Langevin equation of motion and a discussion of the multimode convergence of its characteristic function, which includes Refs. <cit.>], but gives rise to a diverging Purcell rate using Eq. (<ref>). This open JC Model involves a two level approximation (TLA) of the JJ Hilbert space, the rotating wave approximation (RWA) to drop nonresonant contributions, and the Born and Markov approximations leading to a Master equation accounting for losses due to resonator-waveguide coupling. It is unclear which approximation underlies the divergence, or whether the divergence can be resolved within the effective subgap circuit QED field theory.

We first note that keeping only the direct interaction T_ EM-JJ violates gauge invariance. We find that inclusion of all terms, in particular T_EM, equivalent to the diamagnetic A^2 term in the minimal coupling Hamiltonian (p - e A)^2/2m <cit.>, is essential to make all studied QED observables finite. 

The A^2-term is thought to have no impact on transition frequencies in vacuum-induced effects such as the Lamb shift. Because it does not involve atomic operators, it is expected to make the same perturbative contribution to every atomic energy level, precluding observable shifts in transition frequencies <cit.>. This argument relies on perturbation theory in the A^2-term. We show that the diamagnetic term does have an impact when accounted for exactly to all orders.  

Heisenberg equations of motion describing the infinite network in Fig:cQED-open, extending from x=-∞ to x=∞, are <cit.>

    φ̂̈̂_j(t)+ ( 1- γ) ω_j^2sin[φ̂_j(t)]  = γ∂_t^2φ̂(x_0,t),  
       [∂_x^2-χ(x,x_0)∂_t^2]φ̂(x,t)  = χ_sω_j^2 sin[φ̂_j(t)]δ(x-x_0),

Here φ̂_j(t) and φ̂(x,t) are dimensionless flux operators for the JJ and the resonator-waveguide system, respectively, γ≡ C_g/(C_g+C_j) is a capacitive  ratio, χ_s = γ C_j / cL is the dimensionless series capacitance of C_g and C_j, ω_j is the dimensionless transmon frequency, and χ_i≡ C_i/(cL) for i=g,j,R,L <cit.>. These two inhomogeneous equations show that the flux field at x_0 drives the dynamics of the JJ [eqn:TransDyn], while the JJ acts as a source driving the EM fields [eqn:ResDyn]. In addition, the fields are subject to continuity conditions at the ends of the resonator x=0,1 (in units of L).

It is instructive to trace the individual terms of T_int in Eqs. (<ref>-<ref>). T_JJ modifies the qubit frequency, renormalizing γ from C_g/C_j to C_g/(C_g+C_j), while the direct interaction term  T_ EM-JJ gives source terms in both equations. Most importantly, T_EM introduces an effective scattering term in the wave equation describing the fields in the transmission line, by modifying the unitless capacitance per length from 1 to χ(x,x_0) = 1 + χ_s δ(x-x_0). Consequently, these equations are consistent <cit.> with Kirchhoff's law of current conservation. 
In particular, at x=x_0, Eq. (<ref>) yields .∂_xφ̂(x,t)]_x_0^-^x_0^+=χ_s∂_t^2φ̂(x_0,t)+χ_sω_j^2sin[φ̂_j(t)], where the discontinuity in the resonator current is equal to the total current through the capacitive and Josephson branches of the transmon. Similar modification of resonator dynamics has been pointed out before for JJ-based qubits <cit.>.

Equation <ref> can be solved in the Fourier domain, where φ̂̃̂(x,ω) = ∫_-∞^∞ dt  φ̂(x,t) e^-iω t can be expanded in the basis φ̃_n (x,ω) that solves the generalized eigenvalue problem [∂_x^2 + χ(x,x_0)  ω^2 ] φ̃_n (x,ω)=0, subject to continuity conditions at the ends of the resonator, i.e. ∂_x φ̃_n (1^-,ω)= χ_Rω^2[ φ̃_n (1^-,ω)-φ̃_n (1^+,ω)] and  ∂_x φ̃_n (0^+,ω)= χ_Lω^2 [φ̃_n (0^-,ω)-φ̃_n (0^+,ω)], which models the coupling to the waveguides and associated loss. The Dirac δ-function in χ(x,x_0) leads to the discontinuity

    -.∂_x φ̃_n(x)]_x_0^-^x_0^+ = χ_s ω_n^2 φ̃_n(x_0),

resulting in a modified current-conserving (CC) basis <cit.>. These modifications in the spectrum of the transmission line resonator impact the qubit dynamics that is driven by resonator fluctuations.

The role of modal modification in eqn:Current-Conservation can be illustrated with a phenomenological model. Previously, the Purcell rate and the Lamb shift have been calculated using the Lindblad formalism in the dispersive limit <cit.>. An effective multimode JC model

    ℋ̂_JC = ω_j/2σ̂_z  + ∑_n ν_n â_n^†â_n + ∑_n g_n (σ̂^+â_n+σ̂^- â_n^†)

can be obtained from our first principles model <cit.>, which incorporates the modifications to the resonator modes and the qubit dynamics. Resonator losses are included through a Bloch–Redfield equivalent zero-temperature master equation for the reduced density matrix of the resonator and qubit ρ̂̇̂ = -i [ ℋ̂_JC, ρ̂ ] + κ_n ( 2 â_n ρ̂â_n^† - {ρ̂, â_n^†â_n }). The expressions of cavity frequencies ν_n, associated losses κ_n and modal interaction strengths g_n are given in the Supplementary Material <cit.>. All these quantities are functions of χ_s, the strength of the modification of the capacitance per unit length. In particular, the light-matter coupling is found as g_n=1/2γ√(χ_j)√(ω_jν_n)φ̃_n(x_0). We show in Fig. <ref> that g_n is non-monotonic <cit.> for any χ_s≠ 0, first increasing, then turning over at a critical χ_s-dependent mode n, decreasing as g_n ∼ 1/√(n) in the large-n limit <cit.>. This high frequency behavior of g_n renders the multimode Purcell rate finite, without an imposed cutoff [We note that this result is valid in the dispersive limit i.e. away from cavity resonances. In that limit, we expect this result to be fairly accurate when compared to the rate extracted from the exact time evolution of the Master equation for the multimode JC model.].  

This phenomenon is not specific to the resonator geometry in Fig:cQED-open. The underlying physics is the conservation of current at the position x_0 of the qubit. At high frequency, the series capacitance χ_s becomes a short-circuit to ground, acting as a low-pass filter and suppressing mode amplitude at x_0. This is the cause of the power law drop of g_n as n →∞ (subfig:gXrXl1Em3). Moreover, eliminating the continuum degrees of freedom of the waveguides gives an effective decay rate for each mode, κ_n, which increases monotonically as κ_n ∼ n^0.3 (subfig:KappaXrXl1Em3). In the Supplementary Material, we show that for χ_s=0 the resulting series eqn:Multimode Purcell Rate diverges <cit.>, as pointed out in previous studies <cit.>. For any nonzero χ_s, individual terms in the sum (<ref>) display a universal power law ∼ n^-2.7 (subfig:LogLogSpEmRateXrXl1Em3), which guarantees convergence [The power law dependence of κ_n and g_n, though universal with respect to χ_s, are specific to the chosen circuit topology.].




Solution of the Heisenberg-Langevin equations. Although we showed that the expression (<ref>) for the Purcell decay rate converges, it is only valid in the dispersive regime g_n ≪δ_n. This estimate for the Purcell decay rate and the Lamb shift will deviate substantially from the exact result for a range of order g_n around each cavity resonance, diverging as the qubit frequency approaches the resonance (see Fig. <ref>). This fictitious divergence can in principle be cured by solving the full multimode Master equation. Even if computational challenges relating to the long-time dynamics in such a large Hilbert space can be addressed, the resulting rate would still be subject to the TLA, RWA, Born and Markov approximations, casting a priori an uncertainty on its reliability.

An improved analytic result that is uniformly valid in the transmon frequency, and is not limited by the aforementioned approximations can be found by solving eqn:TransDyneqn:ResDyn perturbatively in the transmon's weak nonlinearity. EM degrees of freedom can be integrated out by solving eqn:ResDyn exactly, plugging into eqn:TransDyn and tracing over the photonic Hilbert space. To lowest order in the transmon nonlinearity ϵ = (E_c/E_j)^1/2, where E_c and E_j are the charging and Josephson energy, respectively, the effective equation for the qubit is <cit.>

    Ẍ̂_j(t)+ω_j^2[1-γ+i𝒦_1(0)]X̂_j(t)
    
    =- ω_j^2 ∫_0^tdt' 𝒦_2(t-t') X̂_j(t'),

where X̂_j (t) =  _ph{ρ̂_ph(0)φ̂_j(t)}/ϕ_zpf is the reduced flux operator traced over the photonic degrees of freedom and ϕ_zpf≡(√(2)ϵ)^1/2 is the magnitude of the zero-point phase fluctuations. This delay equation features the memory kernels  𝒦_n(τ)≡γχ_s∫_-∞^+∞dω/2π ω^n   G (x_0,x_0,ω)e^-iωτ, where G(x,x',ω) is the classical EM Green's function defined by [∂_x^2-χ(x,x_0)∂_t^2] G(x,x',ω) e^-i ω t = e^-i ω tδ(x-x') implying that G (x,x',ω) is the amplitude of the flux field created at x by a transmon oscillating with a frequency ω at x' <cit.>. The term on the right hand side of <ref> is therefore proportional to the fluctuating current driving the qubit at time t, that was excited by itself at an earlier time t'. This Green's function correctly encodes the modification of the capacitance per length. Equation (<ref>) can be solved exactly in the Laplace domain

    X̂̃̂_j(s)=sX̂_j(0)+ Ẋ̂_j(0)/D_j(s),

where h̃(s)≡∫_0^∞ dt   h(t)   e^-st, with D_j(s) defined as <cit.>

    D_j(s)≡ s^2+ω_j^2[1-γ+i𝒦_1(0)+𝒦̃_2(s)].

We express the characteristic function D_j(s) in meromorphic form 

    D_j(s)=(s-p_j)(s-p_j^*)∏_m(s-p_m)(s-p_m^*)/(s-z_m)(s-z_m^*).

The poles of 1/D_j(s) are the hybridized qubit-like and resonator-like complex-valued excitation frequencies, p_j≡ -α_j-iβ_j and p_n≡ -α_n-iβ_n, respectively, of the qubit-resonator system, while its zeroes z_n ≡ -iω_n=-κ_n-iν_n correspond to bare non-Hermitian <cit.> cavity resonances. The real part of the qubit-like pole, α_j, is the Purcell loss rate, while β_j - ω_j is the Lamb shift, as shown in Fig:poles_sf_labeled. In the Supplementary Material, we show that  D_j(s) is convergent, and hence so are all hybridized frequencies, for any nonzero χ_s.

The A^2-term kept in our calculation to enforce gauge invariance plays the role of the “counterterm" discussed by Caldeira and Leggett to cancel infinite frequency renormalization <cit.>. This problem has also been discussed in the context of the quantum theory of laser radiation <cit.>.




Perturbative corrections. The transmon nonlinearity neglected in Eq. (<ref>) can be reintroduced as a weak perturbation. The leading order correction to the hybridized resonances amounts to self- and cross-Kerr interactions <cit.>. Using multi-scale perturbation theory <cit.>, the correction to the transmon qubit-like resonance β_j is given by

    β̂_j=β_j-√(2)ϵ/4ω_j[u_j^4ℋ̂_j(0)+∑_n 2u_j^2u_n^2ℋ̂_n(0)]

where the coefficients u_j,n define the transformation from the hybridized to the unhybridized modes and ℋ̂_j,n(0) are the free Hamiltonians of the transmon and mode n, respectively. For χ_g→0, we find u_j→1, u_n=0 and β_j→ω_j such that we recover the frequency correction of free quantum Duffing oscillator ω̂̅̂_j=ω_j[1-√(2)ϵ/4ℋ̂_j(0)] <cit.>. We note three features of this result. Firstly, the correction is an operator and that expresses the fact that transmon levels are anharmonic. The anharmonicity can be calculated from the expectation value of a corrected quadrature operator <cit.>. Secondly, by virtue of the lowest order result being convergent without a cutoff, the perturbative corrections are also convergent in the number of modes included. Finally, this result is not limited by the qubit-resonator coupling strength or the openness of the cavity. The final result is finite for all qubit frequencies, as opposed to the dispersive-limit result. The correction to the Purcell decay is higher order and forms the subject of future work.

We compared the spontaneous decay from the linear theory (blue solid) to the dispersive limit estimate γ_P in Eq. (<ref>) (black dashed) as the transmon frequency is tuned across the fundamental mode in Figs. <ref>-<ref>. 
First, the spontaneous decay is asymmetric, since there are (in)finitely many modes with frequency (larger) smaller than ω_j. This feature is captured by both theories. Second, the spontaneous decay is enhanced as the qubit frequency approaches the fundamental resonator frequency. However, the dispersive limit estimate is perturbative in g_n/δ_n and hence yields a divergent result (fake kink) on resonance regardless of coupling constant, contrary to our result <ref> which predicts a finite value even at ultrastrong coupling (Fig. <ref> and caption).

In Figs. <ref>-<ref> we compare the Lamb shift from the linear theory (blue solid) and the leading order perturbation theory (red dotted) to the dispersive multimode estimate (black dashed) ∑_n g_n^2/δ_n <cit.>. Below the fundamental mode, the Lamb shift is negative due to the collective influence of all higher modes that redshifts the qubit frequency. Above the fundamental mode, there appears a competition between the hybridization with the fundamental mode and all higher modes. Close enough to the fundamental mode, the Lamb shift is positive until it changes sign, as predicted by all three curves.
 
Conclusion. We have presented a framework to calculate the spontaneous decay and the Lamb shift of a transmon qubit, convergent in the number of resonator modes without the need for rotating-wave, two-level, Born or Markov approximations, or a high frequency cutoff. This is achieved by an ab initio treatment of the quantum circuit equations of motion containing the A^2-term to enforce gauge invariance. Therefore, the modes of the resonator are modified such that the light-matter coupling is suppressed at high frequencies. Formulating the cavity resonances in terms of non-Hermitian modes provides access to the spontaneous decay, the Lamb shift, and any other QED observables in a unified way.

Acknowledgements. We acknowledge helpful discussions with Zlatko Minev and S. M. Girvin. This work was supported by the US Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering, under Award No. DE-SC0016011.

Note. While finishing this manuscript we became aware of Ref. Gely_Convergence_2017, which arrives at a similar conclusion for the Lamb shift in the dispersive regime through a different approach.  



Supplementary Material: Cutoff-free Circuit Quantum Electrodynamics




§ HEISENBERG EQUATIONS OF MOTION


In this section, we present the Heisenberg equations of motion in terms of flux variables <cit.>. These equations were derived before by the authors <cit.> (see App. A), but the main steps are summarized below for clarity. The flux variable is defined at any node n in terms of the voltage at that node with respect to a fixed ground node

    Φ_n(t)≡∫_0^tdt' V_n(t').

The classical Lagrangian is the sum of the Lagrangians for the Josephson junction, resonator, right and left waveguides, capacitive coupling between the resonator and the waveguides and the transmon-resonator capacitive coupling, respectively (let U_j(Φ_j) be the nonlinear Josephson potential):  

    ℒ   =1/2C_jΦ̇_j(t)^2-U_j(Φ_j(t))_ℒ_j
       +∫_0^+^L^- dx[1/2 c(∂_tΦ)^2-1/2l(∂_tΦ)^2]_ℒ_Res
       +∫_L^+^∞ dx[1/2 c(∂_tΦ_R)^2-1/2l(∂_xΦ_R)^2]_ℒ_RW
       +∫_-∞^0^- dx[1/2 c(∂_tΦ_L)^2-1/2l(∂_xΦ_L)^2]_ℒ_LW
       +1/2C_L[Φ̇_L(0^-,t)-Φ̇(0^+,t)]^2_ℒ_C_L
       +1/2C_R[Φ̇_R(L^+,t)-Φ̇(L^-,t)]^2_ℒ_C_R
       +1/2C_g[Φ̇_j(t)-Φ̇(x_0,t)]^2_ℒ_C_g,

From Eq. (<ref>) one can derive, via a Legendre transformation followed by quantization <cit.>, the Hamiltonian operator associated with the quantum circuit. The quantum Hamiltonian for C_R,L→ 0 is in Ref. Malekakhlagh_Origin_2016. C_R,L≠ 0 leave equations of motion unchanged, but change boundary conditions (BCs) at x=0,L. Importantly, Heisenberg equations of motion for the quantum flux operators Φ̂_j, Φ̂(x,t) and Φ̂_R,L(x,t) turn out to be formally identical to Euler-Lagrange equations for (<ref>) with classical fields promoted to operators.


To express the Heisenberg equations of motion in a compact way, we introduce the following notations. Φ_0≡h/2e is the superconducting flux quantum and E_j is the Josephson energy. C_s ≡ C_gC_j/(C_g+C_j) is the series capacitance of C_j and C_g and  γ≡ C_g/(C_g+C_j). There is a modified capacitance per unit length in the resonator due to the coupling to the transmon qubit at position x_0:

    c(x,x_0)≡ c+C_sδ(x-x_0).

c and l are the capacitance and inductance per unit length in the resonator and the waveguides.
 
We pass to unitless coordinates and operators (v_p≡ 1/√(lc))

    x→x/L,  
    t→t/L/v_p,  ω→ω/v_pL,
       φ̂≡ 2πΦ̂/Φ_0,   n̂≡Q̂/2e

The newly introduced operators φ̂ and n̂ represent phase and number and are canonically conjugate: [φ̂_j,n̂_j]=i and [φ̂(x,t),n̂(x',t')]=iδ(x-x')δ(t-t'). Below we use unitless capacitances χ_i≡ C_i/(cL), i=R,L,j,g,s, and the unitless capacitance per unit length becomes

    χ(x,x_0)≡ 1+χ_s δ(x-x_0).


In terms of the quantities introduced, the Heisenberg equations of motion for the superconducting phase operators are:


    φ̂̈̂_j(t)+(1-γ)ω_j^2sin[φ̂_j(t)]=γ∂_t^2φ̂(x_0,t),


    [∂_x^2-χ(x,x_0)∂_t^2]φ̂(x,t)=χ_sω_j^2 sin[φ_j(t)]δ(x-x_0),


    ∂_x^2φ̂_R,L(x,t)-∂_t^2φ̂_R,L(x,t)=0,


with boundary conditions


    -.∂_xφ̂|_x=1^-   =-.∂_xφ̂_R|_x=1^+
       =χ_R∂_t^2[φ̂(1^-,t)-φ̂_R(1^+,t)],


    -.∂_xφ̂|_x=0^+   =-.∂_xφ̂_L|_x=0^-
       =χ_L∂_t^2[φ̂_L(0^-,t)-φ̂(0^+,t)],


    φ̂(x=x_0^-,t)=φ̂(x=x_0^+,t),


    . ∂_xφ̂|_x=x_0^+   - .∂_xφ̂|_x=x_0^-
    -χ_s∂_t^2 φ̂(x_0,t)
       =χ_s ω_j^2 sin[φ_j(t)].



In Eqs. (<ref>) and (<ref>), the oscillation frequency is unitless ω_j^2 = 8ℰ_cℰ_j, in terms of unitless Josephson and charging energies

    ℰ_j,c≡√(lc) LE_j,c/ħ,   E_c≡e^2/2C_j.

Equations (<ref>-<ref>) are Eqs. (2-3) in the main text.



§ SPECTRAL REPRESENTATION OF THE GREEN'S FUNCTION

In this section we introduce a spectral representation of the Green's function. The Green's function enters the effective Heisenberg equation of motion for the superconducting phase of the transmon qubit (see Ref. Malekakhlagh_NonMarkovian_2016 for a complete derivation). The resonator Green's function appears if one follows this aim in Eqs. (<ref>) and (<ref>): one has to solve for φ̂(x,t), which is driven by the qubit in Eq. (<ref>), and substitute into (<ref>). The resonator Green's function is defined as the response of the resonator fields, described by the left hand sides of Eqs. (<ref>-<ref>), to a δ-function source in space-time

    [∂_x^2 -χ(x,x_0)∂_t^2]   G(x,t|x_0,t_0)=δ(x-x_0)δ(t-t_0),

obeying BCs (<ref>-<ref>) with φ̂(x,t) replaced by G(x,t|x_0,t_0). Introducing Fourier transforms


    G̃(x,x_0,ω)=∫_-∞^∞dt G(x,t|x_0,t_0) e^+iω(t-t_0), 
       G(x,t|x_0,t_0)=∫_-∞^∞dω/2πG̃(x,x_0,ω)e^-iω(t-t_0),


Equation (<ref>) becomes a Helmholtz equation 

    [∂_x^2 +ω^2χ(x,x_0) ]G̃(x,x_0,ω)=δ(x-x_0).

while the BCs (<ref>-<ref>) are transformed by replacing ∂_t → -iω to


    .∂_xG̃|_x=1^-   =.∂_x G̃|_x=1^+
       =χ_R ω^2 (.G̃|_x=1^--.G̃|_x=1^+),


    .∂_x G̃|_x=0^-   = .∂_x G̃|_x=0^+
       =χ_L ω^2 (.G̃|_x=0^--.G̃|_x=0^+).


    .G̃|_x=x_0^+=.G̃|_x=x_0^-, 
    
       .∂_x G̃|_x=x_0^+-.∂_x G̃|_x=x_0^-+χ_sω^2.G̃|_x=x_0=1,


Lastly, outgoing BCs at infinity model the baths:

    . ∂_x G̃(x,x_0,ω)|_x→±∞=± iωG̃(x→±∞,x_0,ω).

Excitations leaving the resonator never reflect back towards it.



 §.§ Spectral representation of Green's function for χ_R,L = 0


Setting χ_R=χ_L=0 (amounting to a closed resonator) imposes Neumann BC ∂_x G̃|_x=0,1=0 and the problem for G̃ is Hermitian. G̃ can be expanded in terms of a discrete set of normal modes satisfying


    ∂_x^2φ̃_n(x)+χ(x,x_0)ω_n^2φ̃_n(x)=0,
    
       .∂_x φ̃_n(x)|_x=0,1=0.


An important feature of the modes is that their derivative is discontinuous

    -.∂_x φ̃_n(x)|_x_0^-^x_0^+=χ_sω_n^2 φ̃_n(x_0),

Physically, this is the continuity equation at x_0, or current conservation. The mode amplitude at x_0 is suppressed. These observations lead us to name this set of resonator eigenmodes the current-conserving (CC) basis. 

The CC basis eigenfrequencies obey a transcendental equation

    sin(ω_n)+χ_sω_ncos(ω_n x_0)cos[ω_n (1-x_0)]=0,

while the eigenfunctions read

    φ̃_n(x)∝cos[ω_n (1-x_0)]cos(ω_n x),   0<x<x_0
    cos(ω_n x_0)cos[ω_n (1-x)],   x_0<x<1

and the basis is orthonormal over [0,1]: 

    ∫_0^1dxχ(x,x_0)φ̃_m(x)φ̃_n(x)=δ_mn.

Equation (<ref>) can be solved numerically or asymptotically as n→∞, as we do in Sec. <ref>. 

The spectral representation of G̃(x,x',ω) <cit.> is

    G̃(x,x',ω)=∑_n∈ℕφ̃_n(x)φ̃_n(x')/ω^2-ω_n^2=∑_n∈ℤ n≠ 01/2ωφ̃_n(x)φ̃_n(x')/ω-ω_n,

since ω_-n=-ω_n and φ̃_-n(x)=φ̃_n(x).


 §.§ Spectral representation of Green's function for χ_R,L≠ 0

 
If the resonator is open, χ_L,R≠ 0, we resort to a spectral representation in terms of a discrete set of non-Hermitian modes <cit.> that carry constant flux away from the resonator, Constant Flux (CF) modes <cit.>. CF modes satisfy the homogeneous wave equation 

    ∂_x^2φ̃_n(x,ω)+χ(x,x_0)ω_n^2(ω)φ̃_n(x,ω)=0,

with BCs (<ref>)-(<ref>) and (<ref>). Both the modes φ̃_n(x,ω) and their frequencies ω_n(ω) depend on the source frequency ω. 

An outgoing plane wave solution for the left and right waveguides obeying (<ref>), is

    φ̃_n(x,ω)=
    
    A_n^<e^iω_n(ω) x+B_n^< e^-iω_n(ω) x,   0<x<x_0
    
    A_n^>e^iω_n(ω) x+B_n^> e^-iω_n(ω) x,   x_0<x<1
    
    C_ne^iω x,   x>1 
    
    D_ne^-iω x,   x<0

Applying BCs (<ref>-<ref>) leads to a transcendental equation analogous to the closed case which fixes the parametric dependence ω_n(ω) <cit.>.  

The CF modes satisfy now a biorthonormality <cit.> condition

    ∫_0^1dxχ(x,x_0)φ̅̃̅_m^*(x,ω)φ̃_n(x,ω)=δ_mn,

where {φ̅̃̅_m(x,ω)} obey the Hermitian adjoint of (<ref>). φ̃_n(x,ω) and φ̅̃̅_n(x,ω) are right and left eigenfunctions and obey φ̅̃̅_n(x,ω)=φ̃_n^*(x,ω).

The CF mode spectral representation of the Green's function of the open resonator is

    G̃(x,x',ω)=∑_nφ̃_n(x,ω)φ̅̃̅_n^*(x',ω)/ω^2-ω_n^2(ω).

There are two sets of poles of G̃(x,x',ω) in the complex plane. When the denominator of (<ref>) vanishes, ω=ω_n(ω), which corresponds to quasi-bound eigenfrequencies that obey  

    [e^2iω_n-(1-2iχ_Lω_n)(1-2iχ_Rω_n)]
       +i/2χ_sω_n[e^2iω_n x_0+(1-2iχ_Lω_n)]
       ×[e^2iω_n (1-x_0)+(1-2iχ_Rω_n)]=0.

The solutions reside in the lower half of complex ω-plane and come in symmetric pairs with respect to the {ω} axis, i.e. if ω_n  satisfies (<ref>), so does -ω_n^*. Therefore the eigenfrequencies are

    ω_n=
    -iκ_0,       n=0
    
    +ν_n-iκ_n,      n∈+ℕ
    
    -ν_n-iκ_n,      n∈-ℕ

where ν_n > 0 and κ_n >0 are the oscillation frequency and decay rates of quasi-bound mode n, respectively. The dependence of κ_n on mode number n is plotted in Fig. 2 of the main letter. Note the existence of a pole at ω=0, which comes from the ω-dependence of CF states φ̃_n(x,ω) <cit.>. 


§ MULTIMODE JAYNES-CUMMINGS HAMILTONIAN


The classical Hamiltonian for the cQED system can be found from the circuit Lagrangian (<ref>) <cit.>

    ℋ_sys   =4ℰ_c n_j^2(t)-ℰ_jcos[φ_j(t)]
       +∫_0^1dx {n^2(x,t)/2χ(x,x_0)+1/2[∂_x φ(x,t)]^2}
       +2πγ z n_j(t)∫_0^1dxn(x,t)/χ(x,x_0)δ(x-x_0),

where z≡ Z/R_Q where Z≡√(l/c) is the characteristic impedance of the resonator and R_Q≡ h/(2e)^2 is the superconducting resistance quantum. The modification in capacitance per length originates from the system Lagrangian that contains the gauge-invariant qubit-resonator coupling χ_g[φ̇_j(t)-φ̇(x_0,t)]^2/2. In contrast, a phenomenological product coupling χ_gφ̇_j(t)φ̇(x_0,t) would yield a ℋ_sys with χ_s=0 which results in bare resonator modes.

For the purpose of quantizing ℋ_sys, we find the spectrum of the resonator by solving the corresponding Helmholtz eigenvalue problem that has been discussed in Sec. (<ref>). We find


    ℋ̂_sys≡ω_j/4{𝒴̂_j^2-√(2)/ϵcos[(2ϵ^2)^1/4𝒳̂_j]}
       +∑_n{ν_n/4[𝒳̂_n^2+𝒴̂_n^2]+g_n𝒴̂_j𝒴̂_n},

where have defined the canonically conjugate variables 𝒳̂_l≡(â_l+â_l^†) and 𝒴̂_l≡ -i(â_l-â_l^†), where â_l represent the boson annihilation operator of sector l≡ j,c. Moreover, ω_j≡√(8ℰ_jℰ_c) and ϵ≡√(ℰ_c/ℰ_j) is a measure for the strength of transmon nonlinearity. For ϵ=0, we recover ω_j(𝒳̂_j^2+𝒴̂_j^2)/4, the Hamiltonian of a simple harmonic oscillator. In the transmon regime where ϵ≪ 1, the leading contribution is -√(2)ϵω_j𝒳̂_j^4/48. The coupling between qubit and the nth CC mode of the resonator is

    g_n=1/2γ√(χ_j)√(ω_jν_n)φ̃_n(x_0).


There are typically two approaches to diagonalize Eq. (<ref>). In the first approach, assuming that the qubit nonlinearity is strong, one performs a two level reduction. Then, the multimode Rabi Hamiltonian can be derived from Eq. (<ref>) by projecting the quadratures to Pauli sigma matrices, 𝒳̂_j→σ̂^x and 𝒴̂_j→σ̂^y, which yields

    ℋ̂_Rabi   =ω_j/2σ̂^z+∑_nν_nâ_n^†â_n
       -∑_ng_n(â_n-â_n^†)(σ̂^–σ̂^+).

In the rotating wave approximation, Eq. (<ref>) transforms into the multimode Jaynes-Cummings Hamiltonian

    ℋ̂_JC   =ω_j/2σ̂^z+∑_nν_nâ_n^†â_n+∑_ng_n(σ̂^+â_n+σ̂^-â_n^†)

used in the main text. Analytic results can be found for the Purcell decay rate and the Lamb shift in the dispersive limit where g_n≪|ω_j-ω_n| <cit.>. In a Lindblad calculation, resonator losses are included by a Bloch-Redfield approach through the Master equation for the reduced density matrix of the resonator and qubit degrees of freedom ρ̂̇̂ = -i [ ℋ̂_JC, ρ̂ ] + κ_n/2( 2 â_n ρ̂â_n^† - {ρ̂, â_n^†â_n }), where κ_n can be replaced from the solutions to Eq. (<ref>). The second approach treats the nonlinearity as a weak perturbation and is explained in the next section.


§ WEAKLY NONLINEAR TRANSMON

In this section we summarize the steps necessary to derive Eq. (10) of the main text. The full development of multi scale perturbation theory is in Ref. Malekakhlagh_NonMarkovian_2016. By keeping the lowest order nonlinearity (Kerr terms which are quartic in the transmon quadrature), the Hamiltonian can be rewritten in a new basis that diagonalizes the quadratic part

    ℋ̂_sys   ≡β_j/4(𝒳̂̅̂_j^2+𝒴̂̅̂_j^2)+∑_nβ_n/4(𝒳̂̅̂_n^2+𝒴̂̅̂_n^2)
       -εω_j/8(u_j 𝒳̂̅̂_j+∑_nu_n 𝒳̂̅̂_n)^4,

where ε≡√(2)ϵ/6, β_j,n are hybridized frequencies and u_j,n are hybridization coefficients: 𝒳̂_j=u_j𝒳̂̅̂_j+∑_n u_n𝒳̂̅̂_n.

The Heisenberg equations of motion for quadratures become a set of quantum Duffing equations coupled via the quartic terms

    𝒳̂̈̂̅̂̈̂_l(t)+β_l^2{𝒳̂̅̂_l(t)-ε_l[u_j𝒳̂̅̂_j(t)+∑_n u_n𝒳̂̅̂_n(t)]^3}=0,

where ε_l≡ω_j/β_lu_lε for l≡ j,n. Up to lowest order in the perturbation <cit.>, we find an operator valued correction of the linear theory qubit-like frequency β_j: 


    β̂_j=β_j-√(2)ϵ/4ω_j[u_j^4ℋ̂̅̂_j(0)+∑_n2u_j^2u_n^2ℋ̂̅̂_n(0)],

and an analogous correction of the resonator like frequency β_n as

    β̂_n=β_n-√(2)ϵ/4ω_j   [u_n^4ℋ̂̅̂_n(0)+2u_n^2u_j^2ℋ̂̅̂_j(0).
       +.∑_m≠ n2u_n^2u_m^2ℋ̂̅̂_m(0)],


where ℋ̂̅̂_l(0)≡1/4[𝒳̂̅̂_l^2(0)+𝒴̂̅̂_l^2(0)] for l=j,n. In the main text, Eq. (10), the bar notation is dropped. The lowest order MSPT solution for the qubit quadrature becomes, in terms of renormalized frequencies β̂_j,n, <cit.>

    𝒳̂_j^(0)(t)   =u_jâ_j(0)e^-iβ̂_j t+e^-iβ̂_j tâ_j(0)/2cos(3ω_j/4u_j^4ε t)+H.c.
       +∑_n[u_nâ_n(0)e^-iβ̂_n t+e^-iβ̂_n tâ_n(0)/2cos(3ω_j/4u_n^4ε t )+H.c.].

This equation takes into account corrections up to 𝒪(ε) in frequencies. To extract these corrections, we must evaluate the expectation value of Eq. (<ref>) with respect to the initial density matrix. We chose ρ̂=|Ψ⟩_j⟨Ψ|_j⊗|0⟩_ph⟨0|_ph with |Ψ⟩_j=(|0⟩_j+|1⟩_j)/√(2). The correction to the transmon like frequency is obtained from the Fourier components of <𝒳̂_j(t)>. This is the correction plotted in Fig. 3 of the main text.


§ ASYMPTOTIC BEHAVIOR OF LIGHT-MATTER COUPLING


In this section we find the asymptotic behavior of the eigenfrequencies ω_n and eigenmodes φ̃_n(x) of the resonator discussed in the main text. This provides an analytical understanding of the high frequency suppression in the light-matter coupling g_n. 

To point out the origin of the suppression that arise from a nonzero χ_s, let us consider the closed resonator (χ_R,L=0) case. Consider the special case of x_0=0^+ first. This is of experimental interest in order to achieve the maximum coupling to all modes of a resonator. Then, the transcendental Eq. (<ref>) simplifies to

    sin(ω_n)+χ_sω_ncos(ω_n)=0,

which can be rewritten as 

    tan(ω_n)=-χ_sω_n.

The large ω_n solution for χ_s≠0 is then obtained 

    lim_n→∞ω_n=nπ-π/2,
 
which is independent of the value for χ_s. This implies that the effect of a nonzero χ_s on ω_n is a total shift π/2 (half of the free spectral range) in comparison with the case χ_s=0. Substituting x_0=0^+ in Eq. (<ref>), the normalization factor 𝒩_n is found via Eq. (<ref>) as

    ∫_0^1dxcos^2[ω_n(1-x)]+χ_scos^2(ω_n)=1/𝒩_n^2,

which gives

    𝒩_n=√(2)/√(1+χ_scos^2(ω_n)).

Therefore the eigenmode is found as

    φ̃_n(x_0=0^+)=√(2)cos(ω_n)/√(1+χ_scos^2(ω_n)).

Using the trigonometric identity 

    cos^2(ω_n)=1/1+tan^2(ω_n)

and Eq. (<ref>) we can rewrite Eq. (<ref>) as

    φ̃_n(x_0=0^+)=√(2)/√(1+χ_s+χ_s^2ω_n^2),

which now provides the algebraic dependence of φ̃_n(x_0) on ω_n. According to Eq. (<ref>), for large enough ω_n (χ_sω_n ≫ 1+χ_s), we find

    φ̃_n(x_0)1/ω_n,

where the symbol  represents asymptotic equivalence. This imposes a natural cut-off on the light matter coupling for n→∞, since

    g_n∝√(ω_n)φ̃_n(x_0)1/√(ω_n).


Next, we would like to find the asymptotic behavior of ω_n and φ̃_n(x_0) for a general x_0. In order to bring Eq. (<ref>) into a similar form to Eq. (<ref>), we first replace sin(ω_n)=sin[ω_nx_0+ω_n(1-x_0)] and then divide by cos(ω_n x_0)cos[ω_n (1-x_0)] to obtain

    tan(ω_n x_0)+tan[ω_n (1-x_0)]=-χ_sω_n.

Next, the normalization factor 𝒩_n is found from Eq. (<ref>) as


    𝒩_n=√(2)/√(x_0cos^2[ω_n(1-x_0)]+(1-x_0)cos^2(ω_n x_0)+χ_scos^2[ω_n(1-x_0)]cos^2(ω_nx_0)),

Plugging this into Eq. (<ref>) we find

    φ̃_n(x_0)=√(2)/√(1+χ_s+x_0tan^2(ω_nx_0)+(1-x_0)tan^2[ω_n(1-x_0)])


Equations (<ref>) and (<ref>) provide the asymptotic behavior of ω_n, φ̃_n(x_0) and g_n for a general x_0. 


§ CHARACTERISTIC FUNCTION D_J(S) AND ITS CONVERGENCE


In this section we derive the expression for the characteristic function D_j(s) and compare its convergence in number of resonator modes with and without the modification we found for g_n. 

Consider the Heisenberg-Langevin equations of motion corresponding to Hamiltonian (<ref>) in the linear regime (ϵ=0) for 𝒳̂_j,n(t) as


    (d_t^2+ω_j^2)𝒳̂_j(t)=-∑_n2g_nω_n𝒳̂_n(t),
    
       (d_t^2+2κ_n d_t+ω_n^2)𝒳̂_n(t)=-2g_nω_j𝒳̂_j(t)-f̂_n(t),


where κ_n and f̂_n are the decay rate and noise operator coming from coupling to the waveguide degrees of freedom <cit.>.

Equations (<ref>-<ref>) are linear constant coefficient ODEs and can be solved exactly via the unilateral Laplace transform 

    h̃(s)=∫_0^∞dt h(t)e^-st.


Taking the Laplace transform of Eqs. (<ref>-<ref>) we obtain


    (s^2+ω_j^2)𝒳̂̃̂_j(s)+∑_n2g_nω_n 𝒳̂̃̂_n(s)=
    
    s𝒳̂_j(0)+𝒳̂̇̂_j(0),


    (s^2+2κ_ns+ω_n^2)𝒳̂̃̂_n(s)+2g_nω_j𝒳̂̃̂_j(s)=
    
    (s+2κ_n)𝒳̂_n(0)+𝒳̂̇̂_n(0)+f̂̃̂(s).


The solution for 𝒳̂̃̂_j(s) then reads

    𝒳̂̃̂_j(s)=N̂_j(s)/D_j(s),

where the numerator

    N̂_j(s)   =s𝒳̂_j(0)+𝒳̂̇̂_j(0)
       -∑_n2g_nω_n[(s+2κ_n)𝒳̂_n(0)+𝒳̂̇̂_n(0)-f̂̃̂_n(s)]/s^2+2κ_n s+ω_n^2,

contains the operator initial conditions and the denominator

    D_j(s)≡ s^2+ω_j^2-∑_n4g_n^2ω_jω_n/s^2+2κ_n s+ω_n^2.

is the characteristic function whose roots give the hybridized poles of the full system. Therefore, we can represent D_j(s) as

    D_j(s)=(s-p_j)(s-p_j^*)∏_n(s-p_n)(s-p_n^*)/(s-z_n)(s-z_n^*),

where p_j,n≡ -α_j,n-iβ_j,n stand for the transmon-like and the nth resonator-like poles, respectively. Furthermore, z_n ≡-κ_n-i√(ω_n^2-κ_n^2) is the nth bare non-Hermitian resonator mode. The notation (p for poles and z for zeros) is chosen based on 1/D_j(s) that appears in the Laplace solution (<ref>). 

In order to compute the hybridized poles in practice, we need to truncate the number of resonator modes in D_j(s). This truncation is only justified if the function D_j(s) converges as we include more and more modes. First, note that without the correction give by χ_s this sum is divergent, since g_n√(ω)_n√(n) and for a fixed s we obtain

    4g_n^2ω_jω_n/s^2+2κ_n s+ω_n^2ω_n^2/ω_n^2 1.

Hence, the series in divergent. On the other hand, we found that for a non-zero χ_s, g_n 1/√(ω_n) 1/√(n). Therefore we find

    4g_n^2ω_jω_n/s^2+2κ_n s+ω_n^21/ω_n^21/n^2,

and the series becomes convergent. In writing Eq. (<ref>), we used the fact that ω_n n and κ_n has a sublinear asymptotic behavior found numerically.


§ DIVERGENCE IN THE WIGNER-WEISSKOPF THEORY OF SPONTANEOUS EMISSION

Divergence of the Purcell decay rate appears in other frameworks besides the dispersive limit Jaynes-Cummings model as well. In this appendix, we show that the spontaneous decay rate of a qubit coupled to continuum of modes is also divergent, unless the gauge invariance of the interaction is incorporated as presented in this manuscript. The impression of an (erroneous) finite decay rate in free space goes back to Wigner and Weisskopf's original work on spontaneous atomic decay, which implicitly makes a Markov approximation (See Sec. 6.3 of <cit.>). We emphasize that employing the Markov approximation always yields a finite value for the decay rate regardless of the form of spectral function for electromagnetic background. 

To see this explicitly, we go over the Wigner-Weisskopf theory of spontaneous emission for a two-level system coupled to a continuum of modes inside an infinitely long 1D medium.  In interaction picture, the Hamiltonian reads 

    ℋ̂_I=∑_kħ[g_k^*(x_0)σ̂^+â_ke^i(ω_j-ω_k)t+H.c.],

which conserves the total number of excitations 

    N̂≡σ̂^+σ̂^-+∑_k⃗â_k⃗^†â_k⃗.
 
As a result, a number conserving Ansatz for the wavefunction can be written as

    |Ψ(t)⟩=c_e(t)|e,0⟩+∑_kc_g,k(t)|g,1_k⟩,

where there is either no photon in the cavity and the qubit is in excited state |e⟩, or there is a photon at frequency ω_k with qubit in the ground state |g⟩. By solving the Schrodinger equation we obtain the time evolution of the unknown probability amplitudes c_e(t) and c_g,k(t). Combining these equations yields an effective equation for c_e(t) as

    ċ_e(t)=-∫_0^tdt'𝒦(t-t')c_e(t'),

where the memory Kernel 𝒦(τ) is given by

    𝒦(τ)≡∑_k|g_k(x_0)|^2e^i(ω_j-ω_k)t.

Next, we replace the expression for g_k(x_0), derived in Sec. <ref>, as

    |g_k(x_0)|^2=γχ_s/4ω_jω_k|φ̃_k(x_0)|^2.

Note that without respecting the gauge symmetry of interaction |φ̃_k(x_0)|=𝒩(x_0) is k-independent. Moreover, the sum over k can be replaced as

    ∑_k→L/2π∫_0^∞dk=L/2π v_p∫_0^∞dω_k,

for a continuum of modes, where v_p is the phase velocity of the medium. Inserting Eqs. (<ref>) and (<ref>) into the effective Eq. (<ref>) we obtain

    ċ_e(t)   =-1/2πγχ_sω_j 𝒩^2(x_0) L/4v_p
       ×∫_0^∞dω_k ω_k ∫_0^tdt' e^i(ω_j-ω_k)(t-t')c_e(t')

Importantly, the integral over ω_k in Eq. (<ref>) does not converge since the integrand grows unbounded as ω_k→∞. To resolve this, Wigner and Weisskopf assumed that the dominant contribution comes from those modes of continuum whose frequency are close to the qubit frequency. Therefore, the factor ω_k can be replaced by ω_j and by extending the lower limit of integral over ω_k to -∞ we can use the identity

    ∫_-∞^+∞dω_ke^i(ω_j-ω_k)(t-t')=2πδ(t-t'),

to arrive at a finite value for the spontaneous decay as


    ċ_e(t)≈-Γ_sp/2c_e(t),
       Γ_sp≡γχ_sω_j^2 𝒩^2(x_0) L/2v_p.


It is worth mentioning that using Markov approximation, one always obtains a finite expression for the spontaneous decay rate regardless of the form for the spectral function. This happens because instead of integrating over the entire frequency span, the Markov approximation picks a small window around qubit frequency.

Next, we show how our natural high frequency cut-off for light-matter coupling resolves the divergence of Wigner-Weisskopf theory. First, note that applying Markov approximation is indeed unnecessary, since the Volterra Eq. (<ref>) with the memory kernel 

    𝒦(τ)=1/2πγχ_sω_jL/4v_p∫_0^∞dω_k ω_k|φ̃_k(x_0)|^2e^i(ω_j-ω_k)τ,

has an exact solution in Laplace domain as

    c̃_e(s)=c_e(0)/s+𝒦̃(s),

where 𝒦̃(s)≡∫_0^∞dτ𝒦(τ)e^-sτ is the Laplace transform and is found as

    𝒦̃(s)=1/2πγχ_sω_jL/4v_p∫_0^∞dω_k ω_k|φ̃_k(x_0)|^2/s+i(ω_k-ω_j).

Second, when the gauge-invariance of the interaction is incorporated, the mode amplitude is frequency dependent that experiences a high frequency suppression as

    |φ̃_k(x_0)|∼1/ω_k.

Replacing Eq. (<ref>) into expression (<ref>) for 𝒦̃(s) we obtain

    𝒦̃(s)∝∫ dω_k1/ω_k[s+i(ω_k-ω_j)].

Interestingly, with the corrected expression for the eigenmodes, the integrand behaves like 1/ω_k^2 at ω_k→∞, and as a result the integral converges. Otherwise, the integrand behaves like a constant at ω_k→∞ and the result is divergent.

apsrev4-1

