Cutoff-free Circuit Quantum Electrodynamics

TL;DR

This paper demonstrates that respecting gauge invariance is essential in circuit QED calculations to avoid divergences and accurately compute spontaneous emission rates and energy shifts.

Contribution

The authors identify the source of divergence in multimode circuit QED calculations and develop a gauge-invariant method to obtain finite, cutoff-free results.

Findings

Divergences occur when gauge invariance is not maintained.

A gauge-invariant approach yields finite spontaneous emission rates.

The method accurately captures the influence of all cavity modes.

Abstract

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, where the coupling strengths can be substantial, 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.

Equations

$${\gamma }_P=\sum _n(g_n/{\delta }_n{)}^2{\kappa }_n,$$

$${\hat{\ddot{\phi }}}_j(t)+(1-\gamma ){\omega }_j^2\sin [{\hat{\phi }}_j(t)]=\gamma {\partial }_t^2\hat{\phi }(x_0,t),$$

$$\left[{\partial }_x^2-\chi (x,x_0){\partial }_t^2\right]\hat{\phi }(x,t)={\chi }_s{\omega }_j^2\sin [{\hat{\phi }}_j(t)]\delta (x-x_0),$$

$$-{{\partial }_x{\tilde{\phi }}_n(x)]}_{x_0^{-}}^{x_0^{+}}={\chi }_s{\omega }_n^2{\tilde{\phi }}_n(x_0),$$

$${\hat{\mathcal{H}}}_{\text{JC}}=\frac{{\omega }_j}{2}{\hat{\sigma }}_z+\sum _n{\nu }_n{\hat{a}}_n^{†}{\hat{a}}_n+\sum _n{g}_n\left({\hat{\sigma }}^{+}{\hat{a}}_n+{\hat{\sigma }}^{-}{\hat{a}}_n^{†}\right)$$

$$\begin{array}{r}{\hat{\ddot{X}}}_j(t)+{\omega }_j^2\left[1-\gamma +i{\mathcal{K}}_1(0)\right]{\hat{X}}_j(t)\\ =-{\omega }_j^2{\int }_0^{t}d{t}^{\prime }{\mathcal{K}}_2(t-t^{\prime }){\hat{X}}_j(t^{\prime }),\end{array}$$

$${\hat{\tilde{X}}}_j(s)=\frac{s{\hat{X}}_j(0)+{\hat{\dot{X}}}_j(0)}{D_j(s)},$$

$$D_j(s)\equiv s^2+{\omega }_j^2\left[1-\gamma +i{\mathcal{K}}_1(0)+{\tilde{\mathcal{K}}}_2(s)\right].$$

$$D_j(s)=(s-p_j)(s-p_j^{\ast })\prod _m\frac{(s-p_m)(s-p_m^{\ast })}{(s-z_m)(s-z_m^{\ast })}.$$

$${\hat{\beta }}_j={\beta }_j-\frac{\sqrt{2}ϵ}{4}{\omega }_j\left[u_j^4{\hat{\mathcal{H}}}_j(0)+\sum _{n}2u_j^2{u}_n^2{\hat{\mathcal{H}}}_n(0)\right]$$

$${\Phi }_n(t)\equiv {\int }_0^{t}d{t}^{\prime }{V}_n(t^{\prime }).$$

$$\begin{array}{rl}\mathcal{L}& =\underset{{\mathcal{L}}_j}{\underbrace{\frac{1}{2}{C}_j{\dot{\Phi }}_j(t{)}^2-U_j({\Phi }_j(t))}}\\ & +\underset{{\mathcal{L}}_{\text{Res}}}{\underbrace{{\int }_{0^{+}}^{L^{-}}dx\left[\frac{1}{2}c({\partial }_t\Phi {)}^2-\frac{1}{2l}({\partial }_t\Phi {)}^2\right]}}\\ & +\underset{{\mathcal{L}}_{\text{RW}}}{\underbrace{{\int }_{L^{+}}^{\infty }dx\left[\frac{1}{2}c({\partial }_t{\Phi }_R{)}^2-\frac{1}{2l}({\partial }_x{\Phi }_R{)}^2\right]}}\\ & +\underset{{\mathcal{L}}_{\text{LW}}}{\underbrace{{\int }_{-\infty }^{0^{-}}dx\left[\frac{1}{2}c({\partial }_t{\Phi }_L{)}^2-\frac{1}{2l}({\partial }_x{\Phi }_L{)}^2\right]}}\\ & +\underset{{\mathcal{L}}_{C_L}}{\underbrace{\frac{1}{2}{C}_L{\left[{\dot{\Phi }}_L(0^{-},t)-\dot{\Phi }(0^{+},t)\right]}^2}}\\ & +\underset{{\mathcal{L}}_{C_R}}{\underbrace{\frac{1}{2}{C}_R{\left[{\dot{\Phi }}_R(L^{+},t)-\dot{\Phi }(L^{-},t)\right]}^2}}\\ & +\underset{{\mathcal{L}}_{C_g}}{\underbrace{\frac{1}{2}{C}_g{\left[{\dot{\Phi }}_j(t)-\dot{\Phi }(x_0,t)\right]}^2}},\end{array}$$

$$c(x,x_0)\equiv c+C_s\delta (x-x_0).$$

$$\begin{array}{rl}& x\to \frac{x}{L},t\to \frac{t}{\frac{L}{v_p}},\omega \to \frac{\omega }{v_p}L,\\ & \hat{\phi }\equiv 2\pi \frac{\hat{\Phi }}{{\Phi }_0},\hat{n}\equiv \frac{\hat{Q}}{2e}\end{array}$$

$$\chi (x,x_0)\equiv 1+{\chi }_s\delta (x-x_0).$$

$${\hat{\ddot{\phi }}}_j(t)+(1-\gamma ){\omega }_j^2\sin [{\hat{\phi }}_j(t)]=\gamma {\partial }_t^2\hat{\phi }(x_0,t),$$

$$\begin{array}{r}\left[{\partial }_x^2-\chi (x,x_0){\partial }_t^2\right]\hat{\phi }(x,t)={\chi }_s{\omega }_j^2\sin [{\phi }_j(t)]\delta (x-x_0),\end{array}$$

$${\partial }_x^2{\hat{\phi }}_{R,L}(x,t)-{\partial }_t^2{\hat{\phi }}_{R,L}(x,t)=0,$$

$$\begin{array}{rl}-{{\partial }_x\hat{\phi }|}_{x=1^{-}}& =-{{\partial }_x{\hat{\phi }}_R|}_{x=1^{+}}\\ & ={\chi }_R{\partial }_t^2\left[\hat{\phi }(1^{-},t)-{\hat{\phi }}_R(1^{+},t)\right],\end{array}$$

$$\begin{array}{rl}-{{\partial }_x\hat{\phi }|}_{x=0^{+}}& =-{{\partial }_x{\hat{\phi }}_L|}_{x=0^{-}}\\ & ={\chi }_L{\partial }_t^2\left[{\hat{\phi }}_L(0^{-},t)-\hat{\phi }(0^{+},t)\right],\end{array}$$

$$\hat{\phi }(x=x_0^{-},t)=\hat{\phi }(x=x_0^{+},t),$$

$$\begin{array}{rl}{{\partial }_x\hat{\phi }|}_{x=x_0^{+}}& -{{\partial }_x\hat{\phi }|}_{x=x_0^{-}}-{\chi }_s{\partial }_t^2\hat{\phi }(x_0,t)\\ & ={\chi }_s{\omega }_j^2\sin [{\phi }_j(t)].\end{array}$$

$${\mathcal{E}}_{j,c}\equiv \sqrt{lc}L\frac{E_{j,c}}{\hslash },E_c\equiv \frac{e^2}{2C_j}.$$

$$\begin{array}{rl}\left[{\partial }_x^2-\chi (x,x_0){\partial }_t^2\right]& G(x,t|x_0,t_0)=\delta (x-x_0)\delta (t-t_0),\end{array}$$

$$\tilde{G}(x,x_0,\omega )={\int }_{-\infty }^{\infty }dtG(x,t|x_0,t_0)e^{+i\omega (t-t_0)},$$

$$G(x,t|x_0,t_0)={\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }\tilde{G}(x,x_0,\omega )e^{-i\omega (t-t_0)},$$

$$\begin{array}{r}\left[{\partial }_x^2+{\omega }^2\chi (x,x_0)\right]\tilde{G}(x,x_0,\omega )=\delta (x-x_0).\end{array}$$

$$\begin{array}{rl}{{\partial }_x\tilde{G}|}_{x=1^{-}}& ={{\partial }_x\tilde{G}|}_{x=1^{+}}\\ & ={\chi }_R{\omega }^2\left({\tilde{G}|}_{x=1^{-}}-{\tilde{G}|}_{x=1^{+}}\right),\end{array}$$

$$\begin{array}{rl}{{\partial }_x\tilde{G}|}_{x=0^{-}}& ={{\partial }_x\tilde{G}|}_{x=0^{+}}\\ & ={\chi }_L{\omega }^2\left({\tilde{G}|}_{x=0^{-}}-{\tilde{G}|}_{x=0^{+}}\right).\end{array}$$

$${\tilde{G}|}_{x=x_0^{+}}={\tilde{G}|}_{x=x_0^{-}},$$

$${{\partial }_x\tilde{G}|}_{x=x_0^{+}}-{{\partial }_x\tilde{G}|}_{x=x_0^{-}}+{\chi }_s{\omega }^2{\tilde{G}|}_{x=x_0}=1,$$

$${{\partial }_x\tilde{G}(x,x_0,\omega )|}_{x\to \pm \infty }=\pm i\omega \tilde{G}(x\to \pm \infty ,x_0,\omega ).$$

$${\partial }_x^2{\tilde{\phi }}_n(x)+\chi (x,x_0){\omega }_n^2{\tilde{\phi }}_n(x)=0,$$

$${{\partial }_x{\tilde{\phi }}_n(x)|}_{x=0,1}=0.$$

$$-{{\partial }_x{\tilde{\phi }}_n(x)|}_{x_0^{-}}^{x_0^{+}}={\chi }_s{\omega }_n^2{\tilde{\phi }}_n(x_0),$$

$$\sin ({\omega }_n)+{\chi }_s{\omega }_n\cos ({\omega }_n{x}_0)\cos \left[{\omega }_n(1-x_0)\right]=0,$$

Full text

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 Kleppner (1981); Goy et al. (1983); Hulet et al. (1985); Jhe et al. (1987), which is called the Purcell effect Purcell et al. (1946). An accompanying effect is the Lamb shift, a radiative level shift first observed in the microwave spectroscopy of the hydrogen ${}^{2}P_{1/2}-~{}^{2}S_{1/2}$ transition Lamb Jr and Retherford (1947). These quantities have been experimentally accurately characterized for superconducting Josephson junction (JJ) based qubits coupled to coplanar transmission lines Fragner et al. (2008); Houck et al. (2008) and three-dimensional resonators Nigg et al. (2012). In the dispersive regime where a qubit with transition frequency $\omega_{j}$ is far-detuned from the nearest resonant cavity mode (frequency $\nu_{r}$, loss $\kappa_{r}$), single mode expressions exist for the Purcell decay rate, $\gamma_{P}=(g/\delta)^{2}\kappa_{r}$ and the Lamb shift, $\Delta_{L}=g^{2}/\delta$. Here $g$ denotes the coupling between the qubit and the cavity mode and $\delta=\omega_{j}-\nu_{r}$ denotes their detuning Boissonneault et al. (2009). However, for large couplings accessible in circuit QED, the single mode approximation is often inaccurate Fragner et al. (2008); Houck et al. (2008). 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 Houck et al. (2008). 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 Fragner et al. (2008). The Purcell rate has been generalized to account for all modes

[TABLE]

where $g_{n}$ and $\delta_{n}=\omega_{j}-\nu_{n}$ are coupling to and detuning from resonator mode $n$ with frequency $\nu_{n}$ and decay rate $\kappa_{n}$. Expression (1) is divergent without imposing a high-frequency cutoff Houck et al. (2008). Divergences appear as well in the Lamb shift and other vacuum-induced phenomena, e.g. photon-mediated qubit-qubit interactions Filipp et al. (2011). 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 Lamb Jr and Retherford (1947), but less noted for the spontaneous emission rate. Indeed, free space spontaneous emission rate diverges as well, as we show in Malekakhlagh et al. . The finite result by Wigner and Weisskopf Weisskopf and Wigner (1930); Scully and Zubairy (1997) 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 Krimer et al. (2014). So far, no satisfactory theoretical explanation has been given for these divergences. Here we address this issue within the framework of circuit quantum electrodynamics Devoret (1997) (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. Lamb et al. (1987), 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 Koch et al. (2007)Schreier et al. (2008)) capacitively coupled to a transmission-line resonator that in turn is coupled at both ends to semi-infinite waveguides (Fig. 1a). We assign flux variables to nodes, $\Phi_{n}(t)=\int^{t}d\tau\,V_{n}(\tau)$, with $V_{n}(t)$ being the instantaneous voltage at node $n$ with respect to the ground node Devoret (1997); Devoret et al. (2014). Fixing the ground amounts to a particular gauge choice Devoret (1997). For the connection geometry in Fig 1a, the light-matter interaction derives from the energy on the coupling capacitor in the dipole approximation, $T_{\text{int}}=\frac{1}{2}C_{g}[\dot{\Phi}(x_{0})-\dot{\Phi}_{j}]^{2}$ Malekakhlagh et al. , with $x_{0}$ the qubit position. If from the three terms in its expansion, $T_{\text{EM}}=\frac{1}{2}C_{g}\dot{\Phi}(x_{0})^{2}$, $T_{\text{EM-JJ}}=-C_{g}\dot{\Phi}(x_{0})\cdot\dot{\Phi}_{j}$ and $T_{\text{JJ}}=\frac{1}{2}C_{g}\dot{\Phi}_{j}^{2}$, only the direct interaction $T_{\text{EM-JJ}}$ is kept, a multimode JC model in terms of circuit parameters can be derived 111See 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. Scully and Zubairy (1997); Bishop (2010); Devoret et al. (2014); Malekakhlagh et al. (2016); Malekakhlagh and Türeci (2016); Clerk et al. (2010); Morse and Feshbach (1953); Economou (1984); Hassani (2013); Türeci et al. (2006); Boissonneault et al. (2009); Senitzky (1960), but gives rise to a diverging Purcell rate using Eq. (1). 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_{\text{EM-JJ}}$ violates gauge invariance. We find that inclusion of all terms, in particular $T_{\text{EM}}$, equivalent to the diamagnetic $A^{2}$ term in the minimal coupling Hamiltonian $(p-eA)^{2}/2m$ Malekakhlagh and Türeci (2016), 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 Milonni (2013). 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. 1a, extending from $x=-\infty$ to $x=\infty$, are Malekakhlagh et al. (2016, )

[TABLE]

Here $\hat{\varphi}_{j}(t)$ and $\hat{\varphi}(x,t)$ are dimensionless flux operators for the JJ and the resonator-waveguide system, respectively, $\gamma\equiv C_{g}/(C_{g}+C_{j})$ is a capacitive ratio, $\chi_{s}=\gamma C_{j}/cL$ is the dimensionless series capacitance of $C_{g}$ and $C_{j}$, $\omega_{j}$ is the dimensionless transmon frequency, and $\chi_{i}\equiv C_{i}/(cL)$ for $i=g,j,R,L$ Malekakhlagh et al. . These two inhomogeneous equations show that the flux field at $x_{0}$ drives the dynamics of the JJ [Eq. (2)], while the JJ acts as a source driving the EM fields [Eq. (3)]. 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_{\text{int}}$ in Eqs. (2-3). $T_{\text{JJ}}$ modifies the qubit frequency, renormalizing $\gamma$ from $C_{g}/C_{j}$ to $C_{g}/(C_{g}+C_{j})$, while the direct interaction term $T_{\text{EM-JJ}}$ gives source terms in both equations. Most importantly, $T_{\text{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 $\chi(x,x_{0})=1+\chi_{s}\delta(x-x_{0})$. Consequently, these equations are consistent Malekakhlagh and Türeci (2016) with Kirchhoff’s law of current conservation. In particular, at $x=x_{0}$, Eq. (3) yields $\left.\partial_{x}\hat{\varphi}(x,t)\right]_{x_{0}^{-}}^{x_{0}^{+}}=\chi_{s}\partial_{t}^{2}\hat{\varphi}(x_{0},t)+\chi_{s}\omega_{j}^{2}\sin[\hat{\varphi}_{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 Nigg et al. (2012); Bourassa et al. (2012); Malekakhlagh and Türeci (2016).

Equation 3 can be solved in the Fourier domain, where $\hat{\tilde{\varphi}}(x,\omega)=\int_{-\infty}^{\infty}dt\,\hat{\varphi}(x,t)e^{-i\omega t}$ can be expanded in the basis $\tilde{\varphi}_{n}(x,\omega)$ that solves the generalized eigenvalue problem $\left[\partial_{x}^{2}+\chi(x,x_{0})\,\omega^{2}\right]\tilde{\varphi}_{n}(x,\omega)=0$, subject to continuity conditions at the ends of the resonator, i.e. $\partial_{x}\tilde{\varphi}_{n}(1^{-},\omega)=\chi_{R}\omega^{2}[\tilde{\varphi}_{n}(1^{-},\omega)-\tilde{\varphi}_{n}(1^{+},\omega)]$ and $\partial_{x}\tilde{\varphi}_{n}(0^{+},\omega)=\chi_{L}\omega^{2}[\tilde{\varphi}_{n}(0^{-},\omega)-\tilde{\varphi}_{n}(0^{+},\omega)]$, which models the coupling to the waveguides and associated loss. The Dirac $\delta$-function in $\chi(x,x_{0})$ leads to the discontinuity

[TABLE]

resulting in a modified current-conserving (CC) basis Malekakhlagh and Türeci (2016). 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 Eq. (4) 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 Boissonneault et al. (2009). An effective multimode JC model

[TABLE]

can be obtained from our first principles model Malekakhlagh et al. , 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 $\hat{\dot{\rho}}=-i[\hat{\mathcal{H}}_{\text{JC}},\hat{\rho}]+\kappa_{n}\left(2\hat{a}_{n}\hat{\rho}\hat{a}_{n}^{\dagger}-\{\hat{\rho},\hat{a}_{n}^{\dagger}\hat{a}_{n}\}\right)$. The expressions of cavity frequencies $\nu_{n}$, associated losses $\kappa_{n}$ and modal interaction strengths $g_{n}$ are given in the Supplementary Material Malekakhlagh et al. . All these quantities are functions of $\chi_{s}$, the strength of the modification of the capacitance per unit length. In particular, the light-matter coupling is found as $g_{n}=\frac{1}{2}\gamma\sqrt{\chi_{j}}\sqrt{\omega_{j}\nu_{n}}\tilde{\varphi}_{n}(x_{0})$. We show in Fig. 2a that $g_{n}$ is non-monotonic Malekakhlagh and Türeci (2016) for any $\chi_{s}\neq 0$, first increasing, then turning over at a critical $\chi_{s}$-dependent mode $n$, decreasing as $g_{n}\sim 1/\sqrt{n}$ in the large-$n$ limit Malekakhlagh et al. . This high frequency behavior of $g_{n}$ renders the multimode Purcell rate finite, without an imposed cutoff 222We 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. 1a. The underlying physics is the conservation of current at the position $x_{0}$ of the qubit. At high frequency, the series capacitance $\chi_{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\rightarrow\infty$ (Fig. 2a). Moreover, eliminating the continuum degrees of freedom of the waveguides gives an effective decay rate for each mode, $\kappa_{n}$, which increases monotonically as $\kappa_{n}\sim n^{0.3}$ (Fig. 2b). In the Supplementary Material, we show that for $\chi_{s}=0$ the resulting series Eq. (1) diverges Malekakhlagh et al. , as pointed out in previous studies Houck et al. (2008); Filipp et al. (2011). For any nonzero $\chi_{s}$, individual terms in the sum (1) display a universal power law $\sim n^{-2.7}$ (Fig. 2c), which guarantees convergence 333The power law dependence of $\kappa_{n}$ and $g_{n}$, though universal with respect to $\chi_{s}$, are specific to the chosen circuit topology..

Solution of the Heisenberg-Langevin equations. Although we showed that the expression (1) for the Purcell decay rate converges, it is only valid in the dispersive regime $g_{n}\ll\delta_{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. 3). 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 Eqs. (2-3) perturbatively in the transmon’s weak nonlinearity. EM degrees of freedom can be integrated out by solving Eq. (3) exactly, plugging into Eq. (2) and tracing over the photonic Hilbert space. To lowest order in the transmon nonlinearity $\epsilon=(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 Malekakhlagh et al. (2016)

[TABLE]

where $\hat{X}_{j}(t)=\operatorname{Tr}_{ph}\{\hat{\rho}_{ph}(0)\hat{\varphi}_{j}(t)\}/\phi_{\text{zpf}}$ is the reduced flux operator traced over the photonic degrees of freedom and $\phi_{\text{zpf}}\equiv(\sqrt{2}\epsilon)^{1/2}$ is the magnitude of the zero-point phase fluctuations. This delay equation features the memory kernels $\mathcal{K}_{n}(\tau)\equiv\gamma\chi_{s}\int_{-\infty}^{+\infty}\frac{d\omega}{2\pi}\,\omega^{n}\,G(x_{0},x_{0},\omega)e^{-i\omega\tau}$, where $G(x,x^{\prime},\omega)$ is the classical EM Green’s function defined by $\left[\partial_{x}^{2}-\chi(x,x_{0})\partial_{t}^{2}\right]G(x,x^{\prime},\omega)e^{-i\omega t}=e^{-i\omega t}\delta(x-x^{\prime})$ implying that $G(x,x^{\prime},\omega)$ is the amplitude of the flux field created at $x$ by a transmon oscillating with a frequency $\omega$ at $x^{\prime}$ Malekakhlagh et al. (2016). The term on the right hand side of 6 is therefore proportional to the fluctuating current driving the qubit at time $t$, that was excited by itself at an earlier time $t^{\prime}$. This Green’s function correctly encodes the modification of the capacitance per length. Equation (6) can be solved exactly in the Laplace domain

[TABLE]

where $\tilde{h}(s)\equiv\int_{0}^{\infty}dt\,h(t)\,e^{-st}$, with $D_{j}(s)$ defined as Malekakhlagh et al. (2016)

[TABLE]

We express the characteristic function $D_{j}(s)$ in meromorphic form

[TABLE]

The poles of $1/D_{j}(s)$ are the hybridized qubit-like and resonator-like complex-valued excitation frequencies, $p_{j}\equiv-\alpha_{j}-i\beta_{j}$ and $p_{n}\equiv-\alpha_{n}-i\beta_{n}$, respectively, of the qubit-resonator system, while its zeroes $z_{n}\equiv-i\omega_{n}=-\kappa_{n}-i\nu_{n}$ correspond to bare non-Hermitian Malekakhlagh et al. cavity resonances. The real part of the qubit-like pole, $\alpha_{j}$, is the Purcell loss rate, while $\beta_{j}-\omega_{j}$ is the Lamb shift, as shown in Fig. 1b. In the Supplementary Material, we show that $D_{j}(s)$ is convergent, and hence so are all hybridized frequencies, for any nonzero $\chi_{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 Caldeira and Leggett (1983); Leggett (1984). This problem has also been discussed in the context of the quantum theory of laser radiation Schwabl and Thirring (1964).

Perturbative corrections. The transmon nonlinearity neglected in Eq. (6) can be reintroduced as a weak perturbation. The leading order correction to the hybridized resonances amounts to self- and cross-Kerr interactions Nigg et al. (2012); Bourassa et al. (2012). Using multi-scale perturbation theory Bender and Orszag (1999); Malekakhlagh et al. (2016), the correction to the transmon qubit-like resonance $\beta_{j}$ is given by

[TABLE]

where the coefficients $u_{j,n}$ define the transformation from the hybridized to the unhybridized modes and $\hat{\mathcal{H}}_{j,n}(0)$ are the free Hamiltonians of the transmon and mode $n$, respectively. For $\chi_{g}\to 0$, we find $u_{j}\to 1$, $u_{n}=0$ and $\beta_{j}\to\omega_{j}$ such that we recover the frequency correction of free quantum Duffing oscillator $\hat{\bar{\omega}}_{j}=\omega_{j}[1-\frac{\sqrt{2}\epsilon}{4}\hat{\mathcal{H}}_{j}(0)]$ Bender and Bettencourt (1996). 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 Malekakhlagh et al. . 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 $\gamma_{P}$ in Eq. (1) (black dashed) as the transmon frequency is tuned across the fundamental mode in Figs. 3a-3b. First, the spontaneous decay is asymmetric, since there are (in)finitely many modes with frequency (larger) smaller than $\omega_{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}/\delta_{n}$ and hence yields a divergent result (fake kink) on resonance regardless of coupling constant, contrary to our result 9 which predicts a finite value even at ultrastrong coupling (Fig. 3b and caption).

In Figs. 3c-3d 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) $\sum_{n}g_{n}^{2}/\delta_{n}$ Boissonneault et al. (2009). 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 et al., 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

I Heisenberg Equations Of Motion

In this section, we present the Heisenberg equations of motion in terms of flux variables Bishop (2010); Devoret et al. (2014). These equations were derived before by the authors Malekakhlagh et al. (2016) (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

[TABLE]

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}(\Phi_{j})$ be the nonlinear Josephson potential):

[TABLE]

From Eq. (12) one can derive, via a Legendre transformation followed by quantization Devoret (1997); Clerk et al. (2010); Devoret et al. (2014), the Hamiltonian operator associated with the quantum circuit. The quantum Hamiltonian for $C_{R,L}\to 0$ is in Ref. Malekakhlagh and Türeci, 2016. $C_{R,L}\neq 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 $\hat{\Phi}_{j}$, $\hat{\Phi}(x,t)$ and $\hat{\Phi}_{R,L}(x,t)$ turn out to be formally identical to Euler-Lagrange equations for (12) with classical fields promoted to operators.

To express the Heisenberg equations of motion in a compact way, we introduce the following notations. $\Phi_{0}\equiv\frac{h}{2e}$ is the superconducting flux quantum and $E_{j}$ is the Josephson energy. $C_{s}\equiv C_{g}C_{j}/(C_{g}+C_{j})$ is the series capacitance of $C_{j}$ and $C_{g}$ and $\gamma\equiv 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}$:

[TABLE]

$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}\equiv 1/\sqrt{lc}$)

[TABLE]

The newly introduced operators $\hat{\varphi}$ and $\hat{n}$ represent phase and number and are canonically conjugate: $[\hat{\varphi}_{j},\hat{n}_{j}]=i$ and $[\hat{\varphi}(x,t),\hat{n}(x^{\prime},t^{\prime})]=i\delta(x-x^{\prime})\delta(t-t^{\prime})$. Below we use unitless capacitances $\chi_{i}\equiv C_{i}/(cL)$, $i=R,L,j,g,s$, and the unitless capacitance per unit length becomes

[TABLE]

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

[TABLE]

with boundary conditions

[TABLE]

In Eqs. (16a) and (16b), the oscillation frequency is unitless $\omega_{j}^{2}=8\mathcal{E}_{c}\mathcal{E}_{j}$, in terms of unitless Josephson and charging energies

[TABLE]

Equations (16a-16b) are Eqs. ($2$-$3$) in the main text.

II 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 et al., 2016 for a complete derivation). The resonator Green’s function appears if one follows this aim in Eqs. (16a) and (16b): one has to solve for $\hat{\varphi}(x,t)$, which is driven by the qubit in Eq. (16b), and substitute into (16a). The resonator Green’s function is defined as the response of the resonator fields, described by the left hand sides of Eqs. (16b-16c), to a $\delta$-function source in space-time

[TABLE]

obeying BCs (17a-17c) with $\hat{\varphi}(x,t)$ replaced by $G(x,t|x_{0},t_{0})$. Introducing Fourier transforms

[TABLE]

Equation (19) becomes a Helmholtz equation

[TABLE]

while the BCs (17a-17c) are transformed by replacing $\partial_{t}\to-i\omega$ to

[TABLE]

Lastly, outgoing BCs at infinity model the baths:

[TABLE]

Excitations leaving the resonator never reflect back towards it.

II.1 Spectral representation of Green’s function for $\chi_{R,L}=0$

Setting $\chi_{R}=\chi_{L}=0$ (amounting to a closed resonator) imposes Neumann BC $\partial_{x}\tilde{G}|_{x=0,1}=0$ and the problem for $\tilde{G}$ is Hermitian. $\tilde{G}$ can be expanded in terms of a discrete set of normal modes satisfying

[TABLE]

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

[TABLE]

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

[TABLE]

while the eigenfunctions read

[TABLE]

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

[TABLE]

Equation (26) can be solved numerically or asymptotically as $n\to\infty$, as we do in Sec. V.

The spectral representation of $\tilde{G}(x,x^{\prime},\omega)$ Morse and Feshbach (1953); Economou (1984); Hassani (2013) is

[TABLE]

since $\omega_{-n}=-\omega_{n}$ and $\tilde{\varphi}_{-n}(x)=\tilde{\varphi}_{n}(x)$.

II.2 Spectral representation of Green’s function for $\chi_{R,L}\neq 0$

If the resonator is open, $\chi_{L,R}\neq 0$, we resort to a spectral representation in terms of a discrete set of non-Hermitian modes Malekakhlagh et al. (2016) that carry constant flux away from the resonator, Constant Flux (CF) modes Türeci et al. (2006); *Tureci_Strong_2008. CF modes satisfy the homogeneous wave equation

[TABLE]

with BCs (22a)-(22c) and (23). Both the modes $\tilde{\varphi}_{n}(x,\omega)$ and their frequencies $\omega_{n}(\omega)$ depend on the source frequency $\omega$.

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

[TABLE]

Applying BCs (22c-22b) leads to a transcendental equation analogous to the closed case which fixes the parametric dependence $\omega_{n}(\omega)$ Malekakhlagh et al. (2016).

The CF modes satisfy now a biorthonormality Türeci et al. (2006) condition

[TABLE]

where $\{\bar{\tilde{\varphi}}_{m}(x,\omega)\}$ obey the Hermitian adjoint of (30). $\tilde{\varphi}_{n}(x,\omega)$ and $\bar{\tilde{\varphi}}_{n}(x,\omega)$ are right and left eigenfunctions and obey $\bar{\tilde{\varphi}}_{n}(x,\omega)=\tilde{\varphi}_{n}^{*}(x,\omega)$.

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

[TABLE]

There are two sets of poles of $\tilde{G}(x,x^{\prime},\omega)$ in the complex plane. When the denominator of (33) vanishes, $\omega=\omega_{n}(\omega)$, which corresponds to quasi-bound eigenfrequencies that obey

[TABLE]

The solutions reside in the lower half of complex $\omega$-plane and come in symmetric pairs with respect to the $\Im\{\omega\}$ axis, i.e. if $\omega_{n}$ satisfies (34), so does $-\omega_{n}^{*}$. Therefore the eigenfrequencies are

[TABLE]

where $\nu_{n}>0$ and $\kappa_{n}>0$ are the oscillation frequency and decay rates of quasi-bound mode $n$, respectively. The dependence of $\kappa_{n}$ on mode number $n$ is plotted in Fig. $2$ of the main letter. Note the existence of a pole at $\omega=0$, which comes from the $\omega$-dependence of CF states $\tilde{\varphi}_{n}(x,\omega)$ Türeci et al. (2006); *Tureci_Strong_2008.

III Multimode Jaynes-Cummings Hamiltonian

The classical Hamiltonian for the cQED system can be found from the circuit Lagrangian (12) Malekakhlagh and Türeci (2016)

[TABLE]

where $z\equiv Z/R_{Q}$ where $Z\equiv\sqrt{l/c}$ is the characteristic impedance of the resonator and $R_{Q}\equiv 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 $\chi_{g}[\dot{\varphi}_{j}(t)-\dot{\varphi}(x_{0},t)]^{2}/2$. In contrast, a phenomenological product coupling $\chi_{g}\dot{\varphi}_{j}(t)\dot{\varphi}(x_{0},t)$ would yield a $\mathcal{H}_{sys}$ with $\chi_{s}=0$ which results in bare resonator modes.

For the purpose of quantizing $\mathcal{H}_{sys}$, we find the spectrum of the resonator by solving the corresponding Helmholtz eigenvalue problem that has been discussed in Sec. (II.1). We find

[TABLE]

where have defined the canonically conjugate variables $\hat{\mathcal{X}}_{l}\equiv(\hat{a}_{l}+\hat{a}_{l}^{{\dagger}})$ and $\hat{\mathcal{Y}}_{l}\equiv-i(\hat{a}_{l}-\hat{a}_{l}^{{\dagger}})$, where $\hat{a}_{l}$ represent the boson annihilation operator of sector $l\equiv j,c$. Moreover, $\omega_{j}\equiv\sqrt{8\mathcal{E}_{j}\mathcal{E}_{c}}$ and $\epsilon\equiv\sqrt{\mathcal{E}_{c}/\mathcal{E}_{j}}$ is a measure for the strength of transmon nonlinearity. For $\epsilon=0$, we recover $\omega_{j}(\hat{\mathcal{X}}_{j}^{2}+\hat{\mathcal{Y}}_{j}^{2})/4$, the Hamiltonian of a simple harmonic oscillator. In the transmon regime where $\epsilon\ll 1$, the leading contribution is $-\sqrt{2}\epsilon\omega_{j}\hat{\mathcal{X}}_{j}^{4}/48$. The coupling between qubit and the $n$th CC mode of the resonator is

[TABLE]

There are typically two approaches to diagonalize Eq. (37). 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. (37) by projecting the quadratures to Pauli sigma matrices, $\hat{\mathcal{X}}_{j}\rightarrow\hat{\sigma}^{x}$ and $\hat{\mathcal{Y}}_{j}\rightarrow\hat{\sigma}^{y}$, which yields

[TABLE]

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

[TABLE]

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}\ll|\omega_{j}-\omega_{n}|$ Boissonneault et al. (2009). 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 $\hat{\dot{\rho}}=-i[\hat{\mathcal{H}}_{\text{JC}},\hat{\rho}]+\frac{\kappa_{n}}{2}\left(2\hat{a}_{n}\hat{\rho}\hat{a}_{n}^{\dagger}-\{\hat{\rho},\hat{a}_{n}^{\dagger}\hat{a}_{n}\}\right)$, where $\kappa_{n}$ can be replaced from the solutions to Eq. (34). The second approach treats the nonlinearity as a weak perturbation and is explained in the next section.

IV 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 et al., 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

[TABLE]

where $\varepsilon\equiv\sqrt{2}\epsilon/6$, $\beta_{j,n}$ are hybridized frequencies and $u_{j,n}$ are hybridization coefficients: $\hat{\mathcal{X}}_{j}=u_{j}\hat{\bar{\mathcal{X}}}_{j}+\sum\limits_{n}u_{n}\hat{\bar{\mathcal{X}}}_{n}$.

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

[TABLE]

where $\varepsilon_{l}\equiv\frac{\omega_{j}}{\beta_{l}}u_{l}\varepsilon$ for $l\equiv j,n$. Up to lowest order in the perturbation Malekakhlagh et al. (2016), we find an operator valued correction of the linear theory qubit-like frequency $\beta_{j}$:

[TABLE]

where $\hat{\bar{\mathcal{H}}}_{l}(0)\equiv\frac{1}{4}[\hat{\bar{\mathcal{X}}}_{l}^{2}(0)+\hat{\bar{\mathcal{Y}}}_{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 $\hat{{\beta}}_{j,n}$, Malekakhlagh et al. (2016)

[TABLE]

This equation takes into account corrections up to $\mathcal{O}(\varepsilon)$ in frequencies. To extract these corrections, we must evaluate the expectation value of Eq. (44) with respect to the initial density matrix. We chose $\hat{\rho}=\ket{\Psi}_{j}\bra{\Psi}_{j}\otimes\ket{0}_{\text{ph}}\bra{0}_{\text{ph}}$ with $\ket{\Psi}_{j}=(\ket{0}_{j}+\ket{1}_{j})/\sqrt{2}$. The correction to the transmon like frequency is obtained from the Fourier components of $\left<\hat{\mathcal{X}}_{j}(t)\right>$. This is the correction plotted in Fig. $3$ of the main text.

V Asymptotic Behavior of Light-Matter coupling

In this section we find the asymptotic behavior of the eigenfrequencies $\omega_{n}$ and eigenmodes $\tilde{\varphi}_{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 $\chi_{s}$, let us consider the closed resonator ($\chi_{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. (26) simplifies to

[TABLE]

which can be rewritten as

[TABLE]

The large $\omega_{n}$ solution for $\chi_{s}\neq 0$ is then obtained

[TABLE]

which is independent of the value for $\chi_{s}$. This implies that the effect of a nonzero $\chi_{s}$ on $\omega_{n}$ is a total shift $\pi/2$ (half of the free spectral range) in comparison with the case $\chi_{s}=0$. Substituting $x_{0}=0^{+}$ in Eq. (27), the normalization factor $\mathcal{N}_{n}$ is found via Eq. (28) as

[TABLE]

which gives

[TABLE]

Therefore the eigenmode is found as

[TABLE]

Using the trigonometric identity

[TABLE]

and Eq. (46) we can rewrite Eq. (50) as

[TABLE]

which now provides the algebraic dependence of $\tilde{\varphi}_{n}(x_{0})$ on $\omega_{n}$. According to Eq. (52), for large enough $\omega_{n}$ ($\chi_{s}\omega_{n}\gg 1+\chi_{s}$), we find

[TABLE]

where the symbol $\thicksim$ represents asymptotic equivalence. This imposes a natural cut-off on the light matter coupling for $n\to\infty$, since

[TABLE]

Next, we would like to find the asymptotic behavior of $\omega_{n}$ and $\tilde{\varphi}_{n}(x_{0})$ for a general $x_{0}$. In order to bring Eq. (26) into a similar form to Eq. (46), we first replace $\sin(\omega_{n})=\sin[\omega_{n}x_{0}+\omega_{n}(1-x_{0})]$ and then divide by $\cos(\omega_{n}x_{0})\cos[\omega_{n}(1-x_{0})]$ to obtain

[TABLE]

Next, the normalization factor $\mathcal{N}_{n}$ is found from Eq. (28) as

[TABLE]

Plugging this into Eq. (27) we find

[TABLE]

Equations (55) and (57) provide the asymptotic behavior of $\omega_{n}$, $\tilde{\varphi}_{n}(x_{0})$ and $g_{n}$ for a general $x_{0}$.

VI Characteristic function $\large 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 (37) in the linear regime ($\epsilon=0$) for $\hat{\mathcal{X}}_{j,n}(t)$ as

[TABLE]

where $\kappa_{n}$ and $\hat{f}_{n}$ are the decay rate and noise operator coming from coupling to the waveguide degrees of freedom Senitzky (1960).

Equations (58a-58b) are linear constant coefficient ODEs and can be solved exactly via the unilateral Laplace transform

[TABLE]

Taking the Laplace transform of Eqs. (58a-58b) we obtain

[TABLE]

The solution for $\hat{\tilde{\mathcal{X}}}_{j}(s)$ then reads

[TABLE]

where the numerator

[TABLE]

contains the operator initial conditions and the denominator

[TABLE]

is the characteristic function whose roots give the hybridized poles of the full system. Therefore, we can represent $D_{j}(s)$ as

[TABLE]

where $p_{j,n}\equiv-\alpha_{j,n}-i\beta_{j,n}$ stand for the transmon-like and the $n$th resonator-like poles, respectively. Furthermore, $z_{n}\equiv-\kappa_{n}-i\sqrt{\omega_{n}^{2}-\kappa_{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 (61).

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 $\chi_{s}$ this sum is divergent, since $g_{n}\thicksim\sqrt{\omega}_{n}\thicksim\sqrt{n}$ and for a fixed s we obtain

[TABLE]

Hence, the series in divergent. On the other hand, we found that for a non-zero $\chi_{s}$, $g_{n}\thicksim 1/\sqrt{\omega_{n}}\thicksim 1/\sqrt{n}$. Therefore we find

[TABLE]

and the series becomes convergent. In writing Eq. (66), we used the fact that $\omega_{n}\thicksim n$ and $\kappa_{n}$ has a sublinear asymptotic behavior found numerically.

VII 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 Scully and Zubairy (1997)). 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

[TABLE]

which conserves the total number of excitations

[TABLE]

As a result, a number conserving Ansatz for the wavefunction can be written as

[TABLE]

where there is either no photon in the cavity and the qubit is in excited state $\ket{e}$, or there is a photon at frequency $\omega_{k}$ with qubit in the ground state $\ket{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

[TABLE]

where the memory Kernel $\mathcal{K}(\tau)$ is given by

[TABLE]

Next, we replace the expression for $g_{k}(x_{0})$, derived in Sec. III, as

[TABLE]

Note that without respecting the gauge symmetry of interaction $|\tilde{\varphi}_{k}(x_{0})|=\mathcal{N}(x_{0})$ is $k$-independent. Moreover, the sum over $k$ can be replaced as

[TABLE]

for a continuum of modes, where $v_{p}$ is the phase velocity of the medium. Inserting Eqs. (72) and (73) into the effective Eq. (70) we obtain

[TABLE]

Importantly, the integral over $\omega_{k}$ in Eq. (74) does not converge since the integrand grows unbounded as $\omega_{k}\to\infty$. 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 $\omega_{k}$ can be replaced by $\omega_{j}$ and by extending the lower limit of integral over $\omega_{k}$ to $-\infty$ we can use the identity

[TABLE]

to arrive at a finite value for the spontaneous decay as

[TABLE]

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. (70) with the memory kernel

[TABLE]

has an exact solution in Laplace domain as

[TABLE]

where $\tilde{\mathcal{K}}(s)\equiv\int_{0}^{\infty}d\tau\mathcal{K}(\tau)e^{-s\tau}$ is the Laplace transform and is found as

[TABLE]

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

[TABLE]

Replacing Eq. (80) into expression (79) for $\tilde{\mathcal{K}}(s)$ we obtain

[TABLE]

Interestingly, with the corrected expression for the eigenmodes, the integrand behaves like $1/\omega_{k}^{2}$ at $\omega_{k}\to\infty$, and as a result the integral converges. Otherwise, the integrand behaves like a constant at $\omega_{k}\to\infty$ and the result is divergent.