Non-Markovian dynamics of a superconducting qubit in an open multimode resonator

# Non-Markovian dynamics of a superconducting qubit in an open multimode resonator

Moein Malekakhlagh    Alexandru Petrescu    Hakan E. Türeci Department of Electrical Engineering, Princeton University, Princeton, New Jersey, 08544
July 5, 2019
###### Abstract

We study the dynamics of a transmon qubit that is capacitively coupled to an open multimode superconducting resonator. Our effective equations are derived by eliminating resonator degrees of freedom while encoding their effect in the Green’s function of the electromagnetic background. We account for the dissipation of the resonator exactly by employing a spectral representation for the Green’s function in terms of a set of non-Hermitian modes and show that it is possible to derive effective Heisenberg-Langevin equations without resorting to the rotating wave, two level, Born or Markov approximations. A well-behaved time domain perturbation theory is derived to systematically account for the nonlinearity of the transmon. We apply this method to the problem of spontaneous emission, capturing accurately the non-Markovian features of the qubit dynamics, valid for any qubit-resonator coupling strength.

## I Introduction

Superconducting circuits are of interest for gate based quantum information processing Devoret and Martinis (2004); Blais et al. (2007); Devoret and Schoelkopf (2013) and for fundamental studies of collective quantum phenomena away from equilibrium Houck et al. (2012); Schmidt and Koch (2013); Hur et al. (2016). In these circuits, Josephson junctions provide the nonlinearity required to define a qubit or a pseudo-spin degree of freedom, and low loss microwave waveguides and resonators provide a convenient linear environment to mediate interactions between Josephson junctions Majer et al. (2007); Sillanpää et al. (2007); Filipp et al. (2011a, b); Van Loo et al. (2013); Shankar et al. (2013); Kimchi-Schwartz et al. (2016), act as Purcell filters Houck et al. (2008); Jeffrey et al. (2014); Bronn et al. (2015) or as suitable access ports for efficient state preparation and readout. Fabrication capabilities have reached a stage where coherent interactions between multiple qubits occur through a waveguide Van Loo et al. (2013), active coupling elements Roushan et al. (2016) or cavity arrays McKay et al. (2015), while allowing manipulation and readout of individual qubits in the circuit. In addition, experiments started deliberately probing regimes featuring very high qubit coupling strengths Niemczyk et al. (2010); Forn-Díaz et al. (2010); Todorov et al. (2010) or setups where multimode effects cannot be avoided Sundaresan et al. (2015). Accurate modeling of these complex circuits has not only become important for designing such circuits, e.g. to avoid cross talk and filter out the electromagnetic environment, but also for the fundamental question of the collective quantum dynamics of qubits Mlynek et al. (2014). In this work, we introduce a first principles Heisenberg-Langevin framework that accounts for such complexity.

The inadequacy of the standard Cavity QED models based on the interaction of a pseudo-spin degree of freedom with a single cavity mode was recognized early on Houck et al. (2008). In principle, the Rabi model could straightforwardly be extended to include many cavity electromagnetic modes and the remaining qubit transitions (See Sec. III of Malekakhlagh and Türeci (2016)), but this does not provide a computationally viable approach for several reasons. Firstly, we do not know of a systematic approach for the truncation of this multimode multilevel system. Secondly, the truncation itself will depend strongly on the spectral range that is being probed in a given experiment (typically around a transition frequency of the qubit), and the effective model for a given frequency would have to accurately describe the resonator loss in a broad frequency range. It is then unclear whether the Markov approximation would be sufficient to describe such losses.

Multimode effects come to the fore in the accurate computation of the effective Purcell decay of a qubit Houck et al. (2008) or the photon-mediated effective exchange interaction between qubits in the dispersive regime Filipp et al. (2011b), where the perturbation theory is divergent. A phenomenological semiclassical approach to the accurate modeling of Purcell loss has been suggested Houck et al. (2008), based on the availability of the effective impedance seen by the qubit. A full quantum model that incorporates the effective impedance of the linear part of the circuit at its core was later presented Nigg et al. (2012). This approach correctly recognizes that a better behaved perturbation theory in the nonlinearity can be developed if the hybridization of the qubit with the linear multimode environment is taken into account at the outset Bourassa et al. (2012). Incorporating the dressing of the modes into the basis that is used to expand the nonlinearity gives then rise to self- and cross-Kerr interactions between hybridized modes. This basis however does not account for the open nature of the resonator. Qubit loss is then extracted from the poles of the linear circuit impedance at the qubit port, Z(\omega). This quantity can in principle be measured or obtained from a simulation of the classical Maxwell equations. Finding the poles of Z(\omega) through Foster’s theorem introduces potential numerical complications Solgun et al. (2014). Moreover, the interplay of the qubit nonlinearity and dissipation is not addressed within Rayleigh-Schrödinger perturbation theory. An exact treatment of dissipation is important for the calculation of multimode Purcell rates of qubits as well as the dynamics of driven dissipative qubit networks Aron et al. (2016).

The difficulty with incorporating dissipation on equal footing with energetics in open systems is symptomatic of more general issues in the quantization of radiation in finite inhomogeneous media. One of the earliest thorough treatments of this problem Glauber and Lewenstein (1991) proposes to use a complete set of states in the unbounded space including the finite body as a scattering object. This “modes of the universe” approach Lang et al. (1973); Ching et al. (1998) is well-defined but has an impractical aspect: one has to deal with a continuum of modes, and as a consequence simple properties characterizing the scatterer itself (e.g. its resonance frequencies and widths) are not effectively utilized. Several methods have been proposed since then to address this shortcoming, which discussed quantization using quasi-modes (resonances) of the finite-sized open resonator Dalton et al. (1999); Lamprecht and Ritsch (1999); Dutra and Nienhuis (2000); Hackenbroich et al. (2002). Usually, these methods treat the atomic degree of freedom as a two-level system and use the rotating wave and the Markov approximations.

In the present work, rather than using a Hamiltonian description, we derive an effective Heisenberg-Langevin equation to describe the dynamics of a transmon qubit Koch et al. (2007) capacitively coupled to an open multimode resonator (See Fig. (a)a). Our treatment illustrates a general framework that does not rely on the Markov, rotating wave or two level approximations. We show that the electromagnetic degrees of freedom of the entire circuit can be integrated out and appear in the equation of motion through the classical electromagnetic Green’s function (GF) corresponding to the Maxwell operator and the associated boundary conditions. A spectral representation of the GF in terms of a complete set of non-Hermitian modes Türeci et al. (2006); Martin Claassen and Hakan E. Tureci (2010) accounts for dissipative effects from first principles. This requires the solution of a boundary-value problem of the Maxwell operator only in the finite domain of the resonator. Our main result is the effective equation of motion (29), which is a Heisenberg-Langevin Scully and Zubairy (1997); Gardiner and Zoller (2004); Carmichael (1998) integro-differential equation for the phase operator of the transmon. Outgoing fields, which may be desired to calculate the homodyne field at the input of an amplifier chain, can be conveniently related through the GF to the qubit phase operator.

As an immediate application, we use the effective Heisenberg-Langevin equation of motion to study spontaneous emission. The spontaneous emission of a two level system in a finite polarizable medium was calculated Dung et al. (2000) in the Schrödinger-picture in the spirit of Wigner-Weisskopf theory Scully and Zubairy (1997). These calculations are based on a radiation field quantization procedure which incorporates continuity and boundary conditions corresponding to the finite dielectric Matloob et al. (1995); Gruner and Welsch (1996), but only focus on separable geometries where the GF can be calculated semianalytically. A generalization of this methodology to an arbitrary geometry Krimer et al. (2014) uses an expansion of the GF in terms of a set of non-Hermitian modes for the appropriate boundary value problem Türeci et al. (2006); Martin Claassen and Hakan E. Tureci (2010). This approach is able to consistently account for multimode effects where the atom-field coupling strength is of the order of the free spectral range of the cavity Meiser and Meystre (2006); Krimer et al. (2014); Sundaresan et al. (2015) for which the atom is found to emit narrow pulses at the cavity roundtrip period Krimer et al. (2014). A drawback of these previous calculations performed in the Schrödinger picture is that without the rotating wave approximation, no truncation scheme has been proposed so far to reduce the infinite hierarchy of equations to a tractable Hilbert space dimension. The employment of the rotating wave approximation breaks this infinite hierarchy through the existence of a conserved excitation number. The Heisenberg-Langevin method introduced here is valid for arbitrary light-matter coupling, and therefore can access the dynamics accurately where the rotating-wave approximation is not valid.

In summary, our microscopic treatment of the openness is one essential difference between our study and previous works on the collective excitations of circuit-QED systems with a localized Josephson nonlinearity Wallquist et al. (2006); Bourassa et al. (2012); Nigg et al. (2012); Leib et al. (2012). In our work, the lifetime of the collective excitations arises from a proper treatment of the resonator boundary conditions Clerk et al. (2010). The harmonic theory of the coupled transmon-resonator system is exactly solvable via Laplace transform. Transmon qubits typically operate in a weakly nonlinear regime, where charge dispersion is negligible Koch et al. (2007). We treat the Josephson anharmonicity on top of the non-Hermitian linear theory (See Fig (b)b) using multi-scale perturbation theory (MSPT) Bender and Orszag (1999); Nayfeh and Mook (2008); Strogatz (2014). First, it resolves the anomaly of secular contributions in conventional time-domain perturbation theories via a resummation Bender and Orszag (1999); Nayfeh and Mook (2008); Strogatz (2014). While this perturbation theory is equivalent to the Rayleigh-Schrödinger perturbation theory when the electromagnetic environment is closed, it allows a systematic expansion even when the environment is open and the dynamics is non-unitary. Second, we account for the self-Kerr and cross-Kerr interactions Drummond and Walls (1980) between the collective non-Hermitian excitations extending Bourassa et al. (2012); Nigg et al. (2012). Third, treating the transmon qubit as a weakly nonlinear bosonic degree of freedom allows us to include the linear coupling to the environment non-perturbatively. This is unlike the dispersive limit treatment of the light-matter coupling as a perturbation Boissonneault et al. (2009). Therefore, the effective equation of motion is valid for all experimentally accessible coupling strengths Thompson et al. (1992); Wallraff et al. (2004); Reithmaier et al. (2004); Anappara et al. (2009); Niemczyk et al. (2010); Forn-Díaz et al. (2010); Todorov et al. (2010); Sundaresan et al. (2015).

We finally present a perturbative procedure to reduce the computational complexity of the solution of Eq. (29), originating from the enormous Hilbert space size, when the qubit is weakly anharmonic. Electromagnetic degrees of freedom can then be perturbatively traced out resulting in an effective equation of motion (63) in the qubit Hilbert space only, which makes its numerical simulation tractable.

The paper is organized as follows: In Sec. II, we introduce a toy model to familiarize the reader with the main ideas and notation. In Sec. III, we present an ab initio effective Heisenberg picture dynamics for the transmon qubit. The derivation for this effective model has been discussed in detail in Apps. A and B. In Sec. IV.1, we study linear theory of spontaneous emission. In Sec. IV.2, we employ quantum multi-scale perturbation theory to investigate the effective dynamics beyond linear approximation. The details of multi-scale calculations are presented in App. D. In Sec. IV.3 we compare these results with the pure numerical simulation. We summarize the main results of this paper in Sec. V.

## II Toy model

In this section, we discuss a toy model that captures the basic elements of the effective equations (Eq. (29)), which we derive in full microscopic detail in Sec. III. This will also allow us to introduce the notation and concepts relevant to the rest of this paper, in the context of a tractable and well-known model. We consider the single-mode Cavity QED model, consisting of a nonlinear quantum oscillator (qubit) that couples linearly to a single bosonic degree of freedom representing the cavity mode (Fig. 1). This mode itself is coupled to a continuum set of bosons playing the role of the waveguide modes. When the nonlinear oscillator is truncated to the lowest two levels, this reduces to the standard open Rabi Model, which is generally studied using Master equation Ridolfo et al. (2012) or stochastic Schrödinger equation Henriet et al. (2014) approaches. Here we will discuss a Heisenberg-picture approach to arrive at an equation of motion for qubit quadratures. The Hamiltonian for the toy model is (\hbar=1)

 \displaystyle\begin{split}\displaystyle\hat{\mathcal{H}}&\displaystyle\equiv% \frac{\omega_{j}}{4}\left(\hat{\mathcal{X}}_{j}^{2}+\hat{\mathcal{Y}}_{j}^{2}% \right)+\frac{\omega_{j}}{2}U(\hat{\mathcal{X}}_{j})\\ &\displaystyle+\frac{\omega_{c}}{4}\left(\hat{\mathcal{X}}_{c}^{2}+\hat{% \mathcal{Y}}_{c}^{2}\right)+g\hat{\mathcal{Y}}_{j}\hat{\mathcal{Y}}_{c}\\ &\displaystyle+\sum\limits_{b}\left[\frac{\omega_{b}}{4}\left(\hat{\mathcal{X}% }_{b}^{2}+\hat{\mathcal{Y}}_{b}^{2}\right)+g_{b}\hat{\mathcal{Y}}_{c}\hat{% \mathcal{Y}}_{b}\right],\end{split} (1)

where \omega_{j}, \omega_{c} and \omega_{b} are bare oscillation frequencies of qubit, the cavity and the bath modes, respectively. We have defined the canonically conjugate variables

 \displaystyle\hat{\mathcal{X}}_{l}\equiv(\hat{a}_{l}+\hat{a}_{l}^{{\dagger}}),% \ \hat{\mathcal{Y}}_{l}\equiv-i(\hat{a}_{l}-\hat{a}_{l}^{{\dagger}}), (2)

where \hat{a}_{l} represent the boson annihilation operator of sector l\equiv j,c,b. Furthermore, g and g_{b} are qubit-cavity and cavity-bath couplings. U(\hat{\mathcal{X}}_{j}) represents the nonlinear part of the potential shown in Fig. (b)b with a blue spider symbol.

The remainder of this section is structured as follows. In Sec. II.1, we eliminate the cavity and bath degrees of freedom to obtain an effective Heisenberg-Langevin equation of motion for the qubit. We dedicate Sec. II.2 to the resulting characteristic function describing the hybridized modes of the linear theory.

### II.1 Effective dynamics of the qubit

In this subsection, we derive the equations of motion for the Hamiltonian (1). We first integrate out the bath degrees of freedom via Markov approximation to obtain an effective dissipation for the cavity. Then, we eliminate the degrees of freedom of the leaky cavity mode to arrive at an effective equation of motion for the qubit, expressed in terms of the GF of the cavity. The Heisenberg equations of motion are found as

 \displaystyle\hat{\dot{\mathcal{X}}}_{j}(t)=\omega_{j}\hat{\mathcal{Y}}_{j}(t)% +2g\hat{\mathcal{Y}}_{c}(t), (3a) \displaystyle\hat{\dot{\mathcal{Y}}}_{j}(t)=-\omega_{j}\left\{\hat{\mathcal{X}% }_{j}(t)+U^{\prime}[\hat{\mathcal{X}}_{j}(t)]\right\}, (3b) \displaystyle\hat{\dot{\mathcal{X}}}_{c}(t)=\omega_{c}\hat{\mathcal{Y}}_{c}(t)% +2g\hat{\mathcal{Y}}_{j}(t)+\sum\limits_{b}2g_{b}\hat{\mathcal{Y}}_{b}(t), (3c) \displaystyle\hat{\dot{\mathcal{Y}}}_{c}(t)=-\omega_{c}\hat{\mathcal{X}}_{c}(t), (3d) \displaystyle\hat{\dot{\mathcal{X}}}_{b}(t)=\omega_{b}\hat{\mathcal{Y}}_{b}(t)% +2g_{b}\hat{\mathcal{Y}}_{c}(t) (3e) \displaystyle\hat{\dot{\mathcal{Y}}}_{b}(t)=-\omega_{b}\hat{\mathcal{X}}_{b}(t), (3f)

where U^{\prime}[\hat{\mathcal{X}}_{j}]\equiv dU/d\hat{\mathcal{X}}_{j}. Eliminating \hat{\mathcal{Y}}_{j,c,b}(t) using Eqs. (3b), (3d) and (3f) first, and integrating out the bath degree of freedom via Markov approximation Scully and Zubairy (1997); Walls and Milburn (2008) we obtain effective equations for the qubit and cavity as

 \displaystyle\hat{\ddot{\mathcal{X}}}_{j}(t)+\omega_{j}^{2}\left\{\hat{% \mathcal{X}}_{j}(t)+U^{\prime}[\hat{\mathcal{X}}_{j}(t)]\right\}=-2g\omega_{c}% \hat{\mathcal{X}}_{c}(t), (4a) \displaystyle\begin{split}&\displaystyle\hat{\ddot{\mathcal{X}}}_{c}(t)+2% \kappa_{c}\hat{\dot{\mathcal{X}}}_{c}(t)+\omega_{c}^{2}\hat{\mathcal{X}}_{c}(t% )\\ &\displaystyle=-2g\omega_{j}\left\{\hat{\mathcal{X}}_{j}(t)+U^{\prime}[\hat{% \mathcal{X}}_{j}(t)]\right\}-\hat{f}_{B}(t),\end{split} (4b)

where 2\kappa_{c} is the effective dissipation Senitzky (1960); Caldeira and Leggett (1981); Clerk et al. (2010) and \hat{f}_{B}(t) is the noise operator of the bath seen by the cavity

 \displaystyle\hat{f}_{B}(t)=\sum\limits_{b}2g_{b}\left[\omega_{b}\hat{\mathcal% {X}}_{b}(0)\cos(\omega_{b}t)+\hat{\dot{\mathcal{X}}}_{b}(0)\sin(\omega_{b}t)% \right]. (5)

Note that Eq. (4b) is a linear non-homogoneous ODE in terms of \hat{\mathcal{X}}_{c}(t). Therefore, it is possible to find its general solution in terms of its impulse response, i.e. the GF of the associated classical cavity oscillator:

 \displaystyle\ddot{G}_{c}(t,t^{\prime})+2\kappa_{c}\dot{G}_{c}(t,t^{\prime})+% \omega_{c}^{2}G_{c}(t,t^{\prime})=-\delta(t-t^{\prime}). (6)

Following the Fourier transform conventions

 \displaystyle\tilde{G}_{c}(\omega)\equiv\int_{-\infty}^{\infty}dtG_{c}(t,t^{% \prime})e^{i\omega(t-t^{\prime})}, (7a) \displaystyle G_{c}(t,t^{\prime})\equiv\int_{-\infty}^{\infty}\frac{d\omega}{2% \pi}\tilde{G}_{c}(\omega)e^{-i\omega(t-t^{\prime})}, (7b)

we obtain an algebraic solution for \tilde{G}_{c}(\omega) as

 \displaystyle\tilde{G}_{c}(\omega)=\frac{1}{(\omega-\omega_{C})(\omega+\omega_% {C}^{*})}, (8)

with \omega_{C}\equiv\nu_{c}-i\kappa_{c} and \nu_{c}\equiv\sqrt{\omega_{c}^{2}-\kappa_{c}^{2}}. Taking the inverse Fourier transform of Eq. (8) we find the single mode GF of the cavity oscillator

 \displaystyle G_{c}(t,t^{\prime})=-\frac{1}{\nu_{c}}\sin\left[\nu_{c}(t-t^{% \prime})\right]e^{-\kappa_{c}(t-t^{\prime})}\Theta(t-t^{\prime}), (9)

where since the poles of \tilde{G}_{c}(\omega) reside in the lower-half of the complex \omega-plane, G_{c}(t,t^{\prime}) is retarded (causal) and \Theta(t) stands for the Heaviside step function Abramowitz and Stegun (1964).

Then, the general solution to Eq. (4b) can be expressed in terms of G_{c}(t,t^{\prime}) as Morse and Feshbach (1953)

 \displaystyle\begin{split}&\displaystyle\hat{\mathcal{X}}_{c}(t)=2g\omega_{j}% \int_{0}^{t}dt^{\prime}G_{c}(t,t^{\prime})\left\{\hat{\mathcal{X}}_{j}(t^{% \prime})+U^{\prime}[\hat{\mathcal{X}}_{j}(t^{\prime})]\right\}\\ &\displaystyle+\left.\left(\partial_{t^{\prime}}+2\kappa_{c}\right)G_{c}(t,t^{% \prime})\right|_{t^{\prime}=0}\hat{\mathcal{X}}_{c}(0)-G_{c}(t,0)\hat{\dot{% \mathcal{X}}}_{c}(0)\\ &\displaystyle+\int_{0}^{t}dt^{\prime}G_{c}(t,t^{\prime})\hat{f}_{B}(t^{\prime% }).\end{split} (10)

Substituting Eq. (10) into the RHS of Eq. (4a) and defining

 \displaystyle\mathcal{K}(t)\equiv 4g^{2}\frac{\omega_{c}}{\omega_{j}}G_{c}(t,0), (11a) \displaystyle\mathcal{D}(t)\equiv-2g\omega_{c}G_{c}(t,0), (11b) \displaystyle\mathcal{I}(\omega)\equiv-2g\omega_{c}\tilde{G}_{c}(\omega), (11c)

we find the effective dynamics of the nonlinear oscillator in terms of \hat{\mathcal{X}}_{j}(t) as

 \displaystyle\begin{split}&\displaystyle\hat{\ddot{\mathcal{X}}}_{j}(t)+\omega% _{j}^{2}\left\{\hat{\mathcal{X}}_{j}(t)+U^{\prime}[\hat{\mathcal{X}}_{j}(t)]% \right\}=\\ &\displaystyle-\int_{0}^{t}dt^{\prime}\mathcal{K}(t-t^{\prime})\omega_{j}^{2}% \left\{\hat{\mathcal{X}}_{j}(t^{\prime})+U^{\prime}[\hat{\mathcal{X}}_{j}(t^{% \prime})]\right\}\\ &\displaystyle+\int_{0}^{t}dt^{\prime}\mathcal{D}(t-t^{\prime})\hat{f}_{B}(t^{% \prime})\\ &\displaystyle+\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\mathcal{I}(\omega)% \left[(i\omega+2\kappa_{c})\hat{\mathcal{X}}_{c}(0)-\hat{\dot{\mathcal{X}}}_{c% }(0)\right]e^{-i\omega t}.\end{split} (12)

The LHS of Eq. (12) is the free dynamics of the qubit. The first term on the RHS includes the memory of all past events encoded in the memory kernel \mathcal{K}(t). The second term incorporates the influence of bath noise on qubit dynamics and plays the role of a drive term. Finally, the last term captures the effect of the initial operator conditions of the cavity. Note that even though Eq. (12) is an effective equation for the qubit, all operators act on the full Hilbert space of the qubit and the cavity.

### II.2 Linear theory

In the absence of the nonlinearity, i.e. U[\hat{\mathcal{X}}_{j}]=0, Eq. (12) is a linear integro-differential equation that can be solved exactly via unilateral Laplace transform

 \displaystyle\tilde{f}(s)\equiv\int_{0}^{\infty}dte^{-st}f(t), (13)

since the memory integral on the RHS appears as a convolution between the kernel \mathcal{K}(t) and earlier values of \hat{\mathcal{X}}_{j}(t^{\prime}) for 0<t^{\prime}<t. Employing the convolution identity

 \displaystyle\mathfrak{L}\left\{\int_{0}^{t}dt^{\prime}\mathcal{K}(t-t^{\prime% })\hat{\mathcal{X}}_{j}(t^{\prime})\right\}=\tilde{\mathcal{K}}(s)\hat{\tilde{% \mathcal{X}}}_{j}(s), (14)

we find that the Laplace solution to Eq. (12) takes the general form

 \displaystyle\hat{\tilde{\mathcal{X}}}_{j}(s)=\frac{\hat{\mathcal{N}}_{j}(s)}{% D_{j}(s)}, (15)

where the numerator

 \displaystyle\begin{split}\displaystyle\hat{\mathcal{N}}_{j}(s)&\displaystyle=% s\hat{\mathcal{X}}_{j}(0)+\hat{\dot{\mathcal{X}}}_{j}(0)\\ &\displaystyle-\frac{2g\omega_{c}\left[(s+2\kappa_{c})\hat{\mathcal{X}}_{c}(0)% +\hat{\dot{\mathcal{X}}}_{c}(0)-\hat{\tilde{f}}_{B}(s)\right]}{s^{2}+2\kappa_{% c}s+\omega_{c}^{2}},\end{split} (16)

contains the information regarding the initial conditions and the noise operator. The characteristic function D_{j}(s) is defined as

 \displaystyle\begin{split}\displaystyle D_{j}(s)\equiv s^{2}+\omega_{j}^{2}% \left[1+\tilde{\mathcal{K}}(s)\right]&\displaystyle=s^{2}+\omega_{j}^{2}\\ &\displaystyle-\frac{4g^{2}\omega_{j}\omega_{c}}{s^{2}+2\kappa_{c}s+\omega_{c}% ^{2}},\end{split} (17)

which is the denominator of the algebraic Laplace solution (15). Therefore, its roots determine the complex resonances of the coupled system. The poles of D_{j}(s) are, on the other hand, the bare complex frequencies of the dissipative cavity oscillator found before, z_{c}\equiv-i\omega_{C}. Therefore, D_{j}(s) can always be represented formally as

 \displaystyle D_{j}(s)=(s-p_{j})(s-p_{j}^{*})\frac{(s-p_{c})(s-p_{c}^{*})}{(s-% z_{c})(s-z_{c}^{*})}, (18)

where p_{j} and p_{c} are the qubit-like and cavity-like poles such that for g\to 0 we get p_{j}\to-i\omega_{j} and p_{c}\to-i\omega_{C}\equiv z_{c}. In writing Eq. (18), we have used the fact that the roots of a polynomial with real coefficients come in complex conjugate pairs.

It is worth emphasizing that our toy model avoids the rotating wave (RW) approximation. This approximation is known to break down in the ultrastrong coupling regime Bourassa et al. (2009); Anappara et al. (2009); Niemczyk et al. (2010); Forn-Díaz et al. (2010); Todorov et al. (2010). In order to understand its consequence and make a quantitative comparison, we have to find how the RW approximation modifies D_{j}(s). Note that by applying the RW approximation, only the coupling Hamiltonian in Eq. (1) transforms as

 \displaystyle\hat{\mathcal{Y}}_{j}\hat{\mathcal{Y}}_{c}\underset{\text{RW}}{% \longrightarrow}\frac{1}{2}\left(\hat{\mathcal{X}}_{j}\hat{\mathcal{X}}_{c}+% \hat{\mathcal{Y}}_{j}\hat{\mathcal{Y}}_{c}\right). (19)

Then, the modified equations of motion for \hat{\mathcal{X}}_{j}(t) and \hat{\mathcal{X}}_{c}(t) read

 \displaystyle\hat{\ddot{\mathcal{X}}}_{j}(t)+\left(\omega_{j}^{2}+g^{2}\right)% \hat{\mathcal{X}}_{j}(t)=-g(\omega_{j}+\omega_{c})\hat{\mathcal{X}}_{c}(t), (20a) \displaystyle\begin{split}\displaystyle\hat{\ddot{\mathcal{X}}}_{c}(t)+2\kappa% _{c}\hat{\dot{\mathcal{X}}}_{c}(t)&\displaystyle+\left(\omega_{c}^{2}+g^{2}% \right)\hat{\mathcal{X}}_{c}(t)\\ &\displaystyle=-g(\omega_{j}+\omega_{c})\hat{\mathcal{X}}_{j}(t)-\hat{f}_{B}(t% ).\end{split} (20b)

Note that the form of Eqs. (20a-20b) is the same as Eqs. (4a-4b) except for the modified parameters. Therefore, following the same calculation as in Sec. II.1 we find a new characteristic function D_{j}^{RW}(s) which reads

 \displaystyle\begin{split}\displaystyle D_{j}^{RW}(s)&\displaystyle=s^{2}+% \left(\omega_{j}^{2}+g^{2}\right)\\ &\displaystyle-\frac{g^{2}(\omega_{j}+\omega_{c})^{2}}{s^{2}+2\kappa_{c}s+% \left(\omega_{c}^{2}+g^{2}\right)}.\end{split} (21)

We compare the complex roots of D_{j}(s) and D_{j}^{RW}(s) in Fig. 2 as a function of g. For g=0, the poles start from their bare values i\omega_{j} and i\nu_{c}-\kappa_{c} and the results with and without RW match exactly. As g increases both theories predict that the dissipative cavity oscillator passes some of its decay rate to the qubit oscillator. This is seen in Fig. (a)a where the poles move towards each other in the s-plane while the oscillation frequency is almost unchanged. As g is increased more, there is an avoided crossing and the poles resolve into two distinct frequencies. After this point, the predictions from D_{j}(s) and D_{j}^{RW}(s) for p_{j} and p_{c} deviate more significantly. This is more visible in Figs. (b)b and (c)c that show the difference between the two solutions in the complex s-plane. In addition, there is a saturation of the decay rates to half of the bare decay rate of the dissipative cavity oscillator.

In summary, we have obtained the effective equation of motion (12) for the quadrature \hat{\mathcal{X}}_{j}(t) of the nonlinear oscillator. This equation incorporates the effects of memory, initial conditions of the cavity and drive. It admits an exact solution via Laplace transform in the absence of nonlinearity. To lowest order, the Josephson nonlinearity is a time-domain perturbation \propto\hat{\mathcal{X}}_{j}^{3}(t) in Eq. (12). This amounts to a quantum Duffing oscillator Bowen and Milburn (2015) coupled to a linear environment. Time-domain perturbation theory consists of an order by order solution of Eq. (12). A naive application leads to the appearance of resonant coupling between the solutions at successive orders. The resulting solution contains secular contributions, i.e. terms that grow unbounded in time. We present the resolution of this problem using multi-scale perturbation theory (MSPT) Bender and Orszag (1999); Nayfeh and Mook (2008); Strogatz (2014) in Sec. IV.2.

## III Effective dynamics of a transmon qubit

In this section, we present a first principles calculation for the problem of a transmon qubit that couples capacitively to an open multimode resonator (see Fig. 3). Like the toy model in Sec. II, this calculation relies on an effective equation of motion for the transmon qubit quadratures, in which the photonic degrees of freedom are integrated out. In contrast to the toy model where the decay rate was obtained via Markov approximation, we use a microscopic model for dissipation Senitzky (1960); Caldeira and Leggett (1981). We model our bath as a pair of semi-infinite waveguides capacitively coupled to each end of a resonator.

As shown in Fig. 3, the transmon qubit is coupled to a superconducting resonator of finite length L by a capacitance C_{g}. The resonator itself is coupled to the two waveguides at its ends by capacitances C_{R} and C_{L}, respectively. For all these elements, the capacitance and inductance per length are equal and given as c and l, correspondingly. The transmon qubit is characterized by its Josephson energy E_{j}, which is tunable by an external flux bias line (FBL) Johnson et al. (2010), and its charging energy E_{c}, which is related to the capacitor C_{j} as E_{c}=e^{2}/(2C_{j}). The explicit circuit quantization is explained in App. A following a standard approach Devoret (1997); Clerk et al. (2010); Bishop (2010); Devoret et al. (2014). We describe the system in terms of flux operator \hat{\Phi}_{j}(t) for transmon and flux fields \hat{\Phi}(x,t) and \hat{\Phi}_{R,L}(x,t) for the resonator and waveguides.

The dynamics for the quantum flux operators of the transmon and each resonator shown in Fig. 3 is derived in App. A. In what follows, we work with unitless variables

where \Phi_{0}\equiv h/(2e) is the flux quantum and 1/\sqrt{lc} is the phase velocity. We also define unitless parameters

The Heisenberg equation of motion for the transmon reads

 \displaystyle\hat{\ddot{\varphi}}_{j}(t)+(1-\gamma)\omega_{j}^{2}\sin{[\hat{% \varphi}_{j}(t)]}=\gamma\partial_{t}^{2}\hat{\varphi}(x_{0},t), (25)

where \gamma\equiv\chi_{g}/(\chi_{g}+\chi_{j}) is a capacitive ratio, \omega_{j}\equiv\sqrt{8\mathcal{E}_{c}\mathcal{E}_{j}} is the unitless bare transmon frequency and x_{0} is the location of transmon. The phase field \hat{\varphi}(x,t) of the resonator satisfies an inhomogeneous wave equation

 \displaystyle\left[\partial_{x}^{2}-\chi(x,x_{0})\partial_{t}^{2}\right]\hat{% \varphi}(x,t)=\chi_{s}\omega_{j}^{2}\sin{[\hat{\varphi}_{j}(t)]}\delta(x-x_{0}), (26)

where \chi(x,x_{0})=1+\chi_{s}\delta(x-x_{0}) is the unitless capacitance per unit length modified due to coupling to the transmon qubit, and \chi_{s}\equiv\chi_{g}\chi_{j}/(\chi_{g}+\chi_{j}) is the unitless series capacitance of C_{j} and C_{g}. The effect of a nonzero \chi_{s} reflects the modification of the cavity modes due to the action of the transmon as a classical scatterer Malekakhlagh and Türeci (2016). We note that this modification is distinct from, and in addition to, the modification of the cavity modes due to the linear part of the transmon potential discussed in Nigg et al. (2012). Table 1 lists the unitless variables and parameters used in the remainder of this paper.

The flux field in each waveguide obeys a homogeneous wave equation

 \displaystyle\left(\partial_{x}^{2}-\partial_{t}^{2}\right)\hat{\varphi}_{R,L}% (x,t)=0. (27)

The boundary conditions (BC) are derived from conservation of current at each end of the resonator as

 \displaystyle\begin{split}\displaystyle-\left.\partial_{x}\hat{\varphi}\right|% _{x=1^{-}}&\displaystyle=-\left.\partial_{x}\hat{\varphi}_{R}\right|_{x=1^{+}}% \\ &\displaystyle=\chi_{R}\partial_{t}^{2}\left[\hat{\varphi}(1^{-},t)-\hat{% \varphi}_{R}(1^{+},t)\right],\end{split} (28a) \displaystyle\begin{split}\displaystyle-\left.\partial_{x}\hat{\varphi}\right|% _{x=0^{+}}&\displaystyle=-\left.\partial_{x}\hat{\varphi}_{L}\right|_{x=0^{-}}% \\ &\displaystyle=\chi_{L}\partial_{t}^{2}\left[\hat{\varphi}_{L}(0^{-},t)-\hat{% \varphi}(0^{+},t)\right].\end{split} (28b)

Equations (25-28b) completely describe the dynamics of a transmon qubit coupled to an open resonator. Note that according to Eq. (25) the bare dynamics of the transmon is modified due to the force term \gamma\partial_{t}^{2}\hat{\varphi}(x_{0},t). Therefore, in order to find the effective dynamics for the transmon, we need to solve for \hat{\varphi}(x,t) first and evaluate it at the point of connection x=x_{0}. This can be done using the classical electromagnetic GF by virtue of the homogeneous part of Eqs. (26,27) being linear in the quantum fields (see App. B.1). Substituting it into the LHS of Eq. (25) and further simplifying leads to the effective dynamics for the transmon phase operator

 \displaystyle\begin{split}&\displaystyle\hat{\ddot{\varphi}}_{j}(t)+(1-\gamma)% \omega_{j}^{2}\sin{\left[\hat{\varphi}_{j}(t)\right]}=\\ \displaystyle+&\displaystyle\frac{d^{2}}{dt^{2}}\int_{0}^{t}dt^{\prime}% \mathcal{K}_{0}(t-t^{\prime})\omega_{j}^{2}\sin{\left[\hat{\varphi}_{j}(t^{% \prime})\right]}\\ \displaystyle+&\displaystyle\int_{-\infty}^{+\infty}\frac{d\omega}{2\pi}% \mathcal{D}_{R}(\omega)\hat{\tilde{\varphi}}_{R}^{inc}(1^{+},\omega)e^{-i% \omega t}\\ \displaystyle+&\displaystyle\int_{-\infty}^{+\infty}\frac{d\omega}{2\pi}% \mathcal{D}_{L}(\omega)\hat{\tilde{\varphi}}_{L}^{inc}(0^{-},\omega)e^{-i% \omega t}\\ \displaystyle+&\displaystyle\int_{0^{-}}^{1^{+}}dx^{\prime}\int_{-\infty}^{+% \infty}\frac{d\omega}{2\pi}\mathcal{I}(x^{\prime},\omega)\left[i\omega\hat{% \varphi}(x^{\prime},0)-\hat{\dot{\varphi}}(x^{\prime},0)\right]e^{-i\omega t}.% \end{split} (29)

The electromagnetic GF is the basic object that appears in the various kernels constituting the above integro-differential equation:

 \displaystyle\mathcal{K}_{n}(\tau)\equiv\gamma\chi_{s}\int_{-\infty}^{+\infty}% \frac{d\omega}{2\pi}\omega^{n}\tilde{G}(x_{0},x_{0},\omega)e^{-i\omega\tau}, (30a) \displaystyle\mathcal{D}_{R}(\omega)\equiv-2i\gamma\omega^{3}\tilde{G}(x_{0},1% ^{+},\omega), (30b) \displaystyle\mathcal{D}_{L}(\omega)\equiv-2i\gamma\omega^{3}\tilde{G}(x_{0},0% ^{-},\omega), (30c) \displaystyle\mathcal{I}(x^{\prime},\omega)\equiv\gamma\omega^{2}\chi(x^{% \prime},x_{0})\tilde{G}(x_{0},x^{\prime},\omega). (30d)

Equation (29) fully describes the effective dynamics of the transmon phase operator. The various terms appearing in this equation have transparent physical interpretation. The first integral on the RHS of Eq. (29) represents the retarded self-interaction of the qubit. It contains the GF in the form \tilde{G}(x_{0},x_{0},\omega) and describes all processes in which the electromagnetic radiation is emitted from the transmon at x=x_{0} and is scattered back again. We will see later on that this term is chiefly responsible for the spontaneous emission of the qubit. The boundary terms include only the incoming part of the waveguide phase fields. They describe the action of the electromagnetic fluctuations in the waveguides on the qubit, as described by the propagators from cavity interfaces to the qubit, \tilde{G}(x_{0},0^{-},\omega) and \tilde{G}(x_{0},1^{+},\omega). The phase fields \hat{\varphi}_{L}(0^{-},t) and \hat{\varphi}_{R}(1^{+},t) may contain a classical (coherent) part as well. Finally, the last integral adds up all contributions of a nonzero initial value for the electromagnetic field inside the resonator that propagates from the point 0<x^{\prime}<1 to the position of transmon x_{0}.

The solution to the effective dynamics (29) requires knowledge of \tilde{G}(x,x^{\prime},\omega). To this end, we employ the spectral representation of the GF in terms of a set of constant flux (CF) modes Türeci et al. (2006, 2008)

 \displaystyle\tilde{G}(x,x^{\prime},\omega)=\sum\limits_{n}\frac{\tilde{\Phi}_% {n}(x,\omega)\bar{\tilde{\Phi}}_{n}^{*}(x^{\prime},\omega)}{\omega^{2}-\omega_% {n}^{2}(\omega)}, (31)

where \tilde{\Phi}_{n}(x,\omega) and \bar{\tilde{\Phi}}_{n}(x,\omega) are the right and left eigenfunctions of the Helmholtz eigenvalue problem with outgoing BC and hence carry a constant flux when x\to\pm\infty. Note that in this representation, both the CF frequencies \omega_{n}(\omega) and the CF modes \tilde{\Phi}_{n}(x,\omega) parametrically depend on the source frequency \omega. The expressions for \omega_{n}(\omega) and \tilde{\Phi}_{n}(x,\omega) are given in App. B.3.

The poles of the GF are the solutions to \omega=\omega_{n}(\omega) that satisfy the transcendental equation

 \displaystyle\begin{split}&\displaystyle\left[e^{2i\omega_{n}}-(1-2i\chi_{L}% \omega_{n})(1-2i\chi_{R}\omega_{n})\right]\\ &\displaystyle+\frac{i}{2}\chi_{s}\omega_{n}[e^{2i\omega_{n}x_{0}}+(1-2i\chi_{% L}\omega_{n})]\\ &\displaystyle\times[e^{2i\omega_{n}(1-x_{0})}+(1-2i\chi_{R}\omega_{n})]=0.% \end{split} (32)

The solutions to Eq. (32) all reside in the lower half of \omega-plane resulting in a finite lifetime for each mode that is characterized by the imaginary part of \omega_{n}\equiv\nu_{n}-i\kappa_{n}. In Fig. 4 we plotted the decay rate \kappa_{n} versus the oscillation frequency \nu_{n} of the first 100 modes for x_{0}=0 and different values of \chi_{R}=\chi_{L} and \chi_{s}. There is a transition from a super-linear Houck et al. (2008) dependence on mode number for smaller opening to a sub-linear dependence for larger openings. Furthermore, increasing \chi_{s} always decreases the decay rate \kappa_{n}. Intuitively, \chi_{s} is the strength of a \delta-function step in the susceptibility at the position of the transmon. An increase in the average refractive index inside the resonator generally tends to redshift the cavity resonances, while decreasing their decay rate.

In summary, we have derived an effective equation of motion, Eq. (29), for the transmon qubit flux operator \hat{\varphi}_{j}, in which the resonator degrees of freedom enter via the electromagnetic GF \tilde{G}(x,x^{\prime},\omega) given in Eq. (31).

## IV Spontaneous emission into a leaky resonator

In this section, we revisit the problem of spontaneous emission Purcell et al. (1946); Kleppner (1981); Goy et al. (1983); Hulet et al. (1985); Jhe et al. (1987); Dung et al. (2000); Houck et al. (2008); Krimer et al. (2014), where the system starts from the initial density matrix

 \displaystyle\hat{\rho}(0)=\hat{\rho}_{j}(0)\otimes\ket{0}_{ph}\bra{0}_{ph}, (33)

such that the initial excitation exists in the transmon sector of Hilbert space with zero photons in the resonator and waveguides. \hat{\rho}_{j}(0) is a general density matrix in the qubit subspace. For our numerical simulation of the spontaneous emission dynamics in terms of quadratures, we will consider \hat{\rho}_{j}(0)=\ket{\Psi_{j}(0)}\bra{\Psi_{j}(0)} with \ket{\Psi_{j}(0)}=(\ket{0}_{j}+\ket{1}_{j})/\sqrt{2}. The spontaneous emission was conventionally studied through the Markov approximation of the memory term which results only in a modification of the qubit-like pole. This is the Purcell modified spontaneous decay where, depending on the density of the states of the environment, the emission rate can be suppressed or enhanced Purcell et al. (1946); Kleppner (1981); Goy et al. (1983); Hulet et al. (1985); Jhe et al. (1987). We extract the spontaneous decay as the real part of transmon-like pole in a full multimode calculation that is accurate for any qubit-resonator coupling strength.

A product initial density matrix like Eq. (33) allows us to reduce the generic dynamics significantly, since the expectation value of any operator \hat{\mathcal{O}}(t) can be expressed as

 \displaystyle\operatorname{Tr}_{j}\operatorname{Tr}_{ph}\left\{\hat{\rho}_{j}(% 0)\otimes\hat{\rho}_{ph}(0)\hat{\mathcal{O}}(t)\right\}=\operatorname{Tr}_{j}% \left\{\hat{\rho}_{j}(0)\hat{O}(t)\right\} (34)

where \hat{O}\equiv\operatorname{Tr}_{ph}\{\hat{\mathcal{O}}\} is the reduced operator in the Hilbert space of the transmon. Therefore, we define a reduced phase operator

 \displaystyle\hat{\phi}_{j}(t)\equiv\operatorname{Tr}_{ph}\{\hat{\rho}_{ph}(0)% \hat{\varphi}_{j}(t)\}. (35)

In the absence of an external drive, the generic effective dynamics in Eq. (29) reduces to

 \displaystyle\begin{split}&\displaystyle\hat{\ddot{\phi}}_{j}(t)+\omega_{j}^{2% }\left[1-\gamma+i\mathcal{K}_{1}(0)\right]\operatorname{Tr}_{ph}\left\{\hat{% \rho}_{ph}(0)\sin{\left[\hat{\varphi}_{j}(t)\right]}\right\}\\ &\displaystyle=-\int_{0}^{t}dt^{\prime}\mathcal{K}_{2}(t-t^{\prime})\omega_{j}% ^{2}\operatorname{Tr}_{ph}\left\{\hat{\rho}_{ph}(0)\sin{\left[\hat{\varphi}_{j% }(t^{\prime})\right]}\right\}.\end{split} (36)

The derivation of Eq. (36) can be found in Apps. B.5 and B.6.

Note that, due to the sine nonlinearity, Eq. (36) is not closed in terms of \hat{\phi}_{j}(t). However, in the transmon regime Koch et al. (2007), where \mathcal{E}_{j}\gg\mathcal{E}_{c}, the nonlinearity in the spectrum of transmon is weak. This becomes apparent when we work with the unitless quadratures

where \phi_{\text{zpf}}\equiv(2\mathcal{E}_{c}/\mathcal{E}_{j})^{1/4} is the zero-point fluctuation (zpf) phase amplitude. Then, we can expand the nonlinearity in both sides of Eq. (36) as

 \displaystyle\begin{split}\displaystyle\frac{\sin{\left[\hat{\varphi}_{j}(t)% \right]}}{\phi_{\text{zpf}}}&\displaystyle=\frac{\hat{\varphi}_{j}(t)}{\phi_{% \text{zpf}}}-\frac{\hat{\varphi}_{j}^{3}(t)}{3!\phi_{\text{zpf}}}+\mathcal{O}% \left[\frac{\hat{\varphi}_{j}^{5}(t)}{\phi_{\text{zpf}}}\right]\\ &\displaystyle=\hat{\mathcal{X}}_{j}(t)-\frac{\sqrt{2}\epsilon}{6}\hat{% \mathcal{X}}_{j}^{3}(t)+\mathcal{O}\left(\epsilon^{2}\right),\end{split} (38)

where \epsilon\equiv(\mathcal{E}_{c}/\mathcal{E}_{j})^{1/2} appears as a measure for the strength of the nonlinearity. In experiment, the Josephson energy \mathcal{E}_{j} can be tuned through the FBL while the charging energy \mathcal{E}_{c} is fixed. Therefore, a higher transmon frequency \omega_{j}=\sqrt{8\mathcal{E}_{c}\mathcal{E}_{j}} is generally associated with a smaller \epsilon and hence weaker nonlinearity.

The remainder of this section is organized as follows. In Sec. IV.1 we study the linear theory. In Sec. IV.2 we develop a perturbation expansion up to leading order in \epsilon. In Sec. IV.3, we compare our analytical results with numerical simulation. Finally, in Sec. IV.4 we discuss the output response of the cQED system that can be probed in experiment.

### IV.1 Linear theory

In this subsection, we solve the linear effective dynamics and discuss hybridization of the transmon and the resonator resonances. We emphasize the importance of off-resonant modes as the coupling \chi_{g} is increased. We next investigate the spontaneous decay rate as a function of transmon frequency \omega_{j} and coupling \chi_{g} and find an asymmetric dependence on \omega_{j} in agreement with a previous experiment Houck et al. (2008).

Neglecting the cubic term in Eq. (38), the partial trace with respect to the resonator modes can be taken directly and we obtain the effective dynamics

 \displaystyle\begin{split}\displaystyle\hat{\ddot{X}}_{j}(t)&\displaystyle+% \omega_{j}^{2}\left[1-\gamma+i\mathcal{K}_{1}(0)\right]\hat{X}_{j}(t)\\ &\displaystyle=-\int_{0}^{t}dt^{\prime}\mathcal{K}_{2}(t-t^{\prime})\omega_{j}% ^{2}\hat{X}_{j}(t^{\prime}).\end{split} (39)

Then, using Laplace transform we can solve Eq. (39) as

 \displaystyle\hat{\tilde{X}}_{j}(s)=\frac{s\hat{X}_{j}(0)+\omega_{j}\hat{Y}_{j% }(0)}{D_{j}(s)}, (40)

with D_{j}(s) defined as

 \displaystyle D_{j}(s)\equiv s^{2}+\omega_{j}^{2}\left[1-\gamma+i\mathcal{K}_{% 1}(0)+\tilde{\mathcal{K}}_{2}(s)\right]. (41)

Equations (40) and (41) contain the solution for the reduced quadrature operator of the transmon qubit in the Laplace domain.

In order to find the time domain solution, it is necessary to study the poles of Eq. (40) and consequently the roots of D_{j}(s). The characteristic function D_{j}(s) can be expressed as (see App. C)

 \displaystyle\begin{split}&\displaystyle D_{j}(s)=s^{2}+\omega_{j}^{2}+\\ &\displaystyle\omega_{j}^{2}\left\{-\gamma+\sum\limits_{n}M_{n}\frac{s\{\cos{[% 2\delta_{n}(x_{0})]}s+\sin{[2\delta_{n}(x_{0})]}\nu_{n}\}}{(s+\kappa_{n})^{2}+% \nu_{n}^{2}}\right\},\end{split} (42)

where \delta_{n}(x) is the phase of the non-Hermitian eigenfunction such that \tilde{\Phi}_{n}(x)=|\tilde{\Phi}_{n}(x)|e^{i\delta_{n}(x)}. We identify the term

 \displaystyle M_{n}\equiv\gamma\chi_{s}|\tilde{\Phi}_{n}(x_{0})|^{2} (43)

as the measure of hybridization with individual resonator modes. The form of M_{n} in Eq. (43) illustrates that the hybridization between the transmon and the resonator is bounded. This strength of hybridization is parameterized by \gamma\chi_{s} rather than \chi_{g}. This implies that as \chi_{g}, the coupling capacitance, is increased, the qubit-resonator hybridization is limited by the internal capacitance of the qubit, \chi_{j}:

 \displaystyle\lim\limits_{\frac{\chi_{g}}{\chi_{j}}\to\infty}\gamma\chi_{s}=% \lim\limits_{\frac{\chi_{g}}{\chi_{j}}\to\infty}\left(\frac{\chi_{g}}{\chi_{g}% +\chi_{j}}\right)^{2}\chi_{j}=\chi_{j}. (44)

For this reason, our numerical results below feature a saturation in hybridization as \chi_{g} is increased.

The roots of D_{j}(s) are the hybridized poles of the entire system. If there is no coupling, i.e. \chi_{g}=0, then D_{j}(s)=s^{2}+\omega_{j}^{2}=(s+i\omega_{j})(s-i\omega_{j}) is the characteristic polynomial that gives the bare transmon resonance. However, for a nonzero \chi_{g}, D_{j}(s) becomes a meromorphic function whose zeros are the hybridized resonances of the entire system, and whose poles are the bare cavity resonances. Therefore, D_{j}(s) can be expressed as

 \displaystyle D_{j}(s)=(s-p_{j})(s-p_{j}^{*})\prod\limits_{m}\frac{(s-p_{m})(s% -p_{m}^{*})}{(s-z_{m})(s-z_{m}^{*})}. (45)

In Eq. (45), p_{j}\equiv-\alpha_{j}-i\beta_{j} and p_{n}\equiv-\alpha_{n}-i\beta_{n} are the zeros of D_{j}(s) that represent the transmon-like and the nth resonator-like poles, accordingly. Furthermore, z_{n}\equiv-i\omega_{n}=-\kappa_{n}-i\nu_{n} stands for the nth bare non-Hermitian resonator resonance. The notation chosen here (p for poles and z for zeroes) reflects the meromorphic structure of 1/D_{j}(s) which enters the solution Eq. (40).

An important question concerns the convergence of D_{j}(s) as a function of the number of the resonator modes included in the calculation. The form of D_{j}(s) given in Eq. (45) is suitable for this discussion. Consider the factor corresponding to the mth resonator mode in 1/D_{j}(s). We reexpress it as

 \displaystyle\begin{split}\displaystyle\frac{(s-z_{m})(s-z_{m}^{*})}{(s-p_{m})% (s-p_{m}^{*})}&\displaystyle=\left(1-\frac{z_{m}-p_{m}}{s-p_{m}}\right)\left(1% -\frac{z_{m}^{*}-p_{m}^{*}}{s-p_{m}^{*}}\right)\\ &\displaystyle=1+\mathcal{O}\left(\left|\frac{z_{m}-p_{m}}{s-p_{m}}\right|% \right).\end{split} (46)

The consequence of a small shift |p_{m}-z_{m}| as compared to the strongly hybridized resonant mode |p_{1}-z_{1}| is that it can be neglected in the expansion for 1/D_{j}(s). The relative size of these contributions is controlled by the coupling \chi_{g}. As rule of thumb, the less hybridized a resonator pole is, the less it contributes to qubit dynamics. Ultimately, the truncation in this work is established by imposing the convergence of the numerics.

A numerical solution for the roots of Eq. (42) at weak coupling \chi_{g} reveals that the mode resonant with the transmon is significantly shifted, with comparatively small shifts |p_{m}-z_{m}| in the other resonator modes (See Fig. 5). At weak coupling, the hybridization of p_{j} and p_{1} is captured by a single resonator mode.

Next, we plot in Fig. 6 the effect of truncation on the response of the multimode system in a band around s=p_{j}. As the coupling \chi_{g} is increased beyond the avoided crossing, which is also captured by the single mode truncation, the effect of off-resonant modes on p_{j} and p_{1} becomes significant. It is important to note that the hybridization occurs in the complex s-plane. On the frequency axis \text{Im}\{s\} an increase in \chi_{g} is associated with a splitting of transmon-like and resonator-like poles. Along the decay rate axis \text{Re}\{s\} we notice that the qubit decay rate is controlled by the resonant mode at weak coupling, with noticeable enhancement of off-resonant mode contribution at strong coupling. If the truncation is not done properly in the strong coupling regime, it may result in spurious unstable roots of D_{j}(s), i.e. \text{Re}\{s\}>0, as seen in Fig. (a)a.

The modification of the decay rate of the transmon-like pole, henceforth identified as \alpha_{j}\equiv-\text{Re}\{p_{j}\}, has an important physical significance. It describes the Purcell modification of the qubit decay (if sources for qubit decay other than the direct coupling to electromagnetic modes can be neglected). The present scheme is able to capture the full multimode modification, that is out of the reach of conventional single-mode theories of spontaneous emission Purcell et al. (1946); Kleppner (1981); Goy et al. (1983); Hulet et al. (1985); Jhe et al. (1987).

At fixed \chi_{g}, we observe an asymmetry of \alpha_{j} when the bare transmon frequency is tuned across the fundamental mode of the resonator, in agreement with a previous experiment Houck et al. (2008), where a semiclassical model was employed for an accurate fit. Figure 7 shows that near the resonator-like resonance the spontaneous decay rate is enhanced, as expected. For positive detunings spontaneous decay is significantly larger than for negative detunings, which can be traced back to an asymmetry in the resonator density of states Houck et al. (2008). We find that this asymmetry grows as \chi_{g} is increased. Note that besides a systematic inclusion of multimode effects, the presented theory of spontaneous emission goes beyond the rotating wave, Markov and two-level approximations as well.

Having studied the hybridized resonances of the entire system, we are now able to provide the time-dependent solution to Eq. (39). By substituting Eq. (45) into Eq. (40) we obtain

 \displaystyle\hat{\tilde{X}}_{j}(s)=\left(\frac{\hat{A}_{j}}{s-p_{j}}+\sum% \limits_{n}\frac{\hat{A}_{n}}{s-p_{n}}\right)+H.c., (47)

from which the inverse Laplace transform is immediate

 \displaystyle\hat{X}_{j}(t)=\left[\left(\hat{A}_{j}e^{p_{j}t}+\sum\limits_{n}% \hat{A}_{n}e^{p_{n}t}\right)+H.c.\right]\Theta(t). (48)

The frequency components have operator-valued amplitudes

 \displaystyle\hat{A}_{j}\equiv A_{j}^{X}\hat{X}_{j}(0)+A_{j}^{Y}\hat{Y}_{j}(0), (49a) \displaystyle\hat{A}_{n}\equiv A_{n}^{X}\hat{X}_{j}(0)+A_{n}^{Y}\hat{Y}_{j}(0), (49b)

with the residues given in terms of D_{j}(s) as

 \displaystyle A_{j,n}^{X}\equiv\left.\left[(s-p_{j,n})\frac{s}{D_{j}(s)}\right% ]\right|_{s=p_{j,n}}, (50a) \displaystyle A_{j,n}^{Y}\equiv\left.\left[(s-p_{j,n})\frac{\omega_{j}}{D_{j}(% s)}\right]\right|_{s=p_{j,n}}. (50b)

The dependence of A_{j,n}^{X} and A_{j,n}^{Y} on coupling \chi_{g} has been studied in Fig. 8. The transmon-like amplitude (blue solid) is always dominant, and further off-resonant modes have smaller amplitudes. By increasing \chi_{g}, the resonator-like amplitude grow significantly first and reach an asymptote as predicted by Eq. (44).

### IV.2 Perturbative corrections

In this section, we develop a well-behaved time-domain perturbative expansion in the transmon qubit nonlinearity as illustrated in Eq. (38). Conventional time-domain perturbation theory is inapplicable due to the appearance of resonant coupling between the successive orders which leads to secular contributions, i.e. terms that grow unbounded in time (For a simple example see App. D.1). A solution to this is multi-scale perturbation theory (MSPT) Bender and Orszag (1999); Nayfeh and Mook (2008); Strogatz (2014), which considers multiple independent time scales and eliminates secular contributions by a resummation of the conventional perturbation series.

The effect of the nonlinearity is to mix the hybridized modes discussed in the previous section, leading to transmon mediated self-Kerr and cross-Kerr interactions. Below, we extend MSPT to treat this problem while consistently accounting for the dissipative effects. This goes beyond the extent of Rayleigh-Schrödinger perturbation theory, as it will allow us to treat the energetic and dissipative scales on equal footing.

The outcome of conventional MSPT analysis in a conservative system is frequency renormalization Bender and Orszag (1999); Bender and Bettencourt (1996). We illustrate this point for a classical Duffing oscillator, which amounts to the classical theory of an isolated transmon qubit up to leading order in the nonlinearity. We outline the main steps here leaving the details to App. D.1. Consider a classical Duffing oscillator

 \displaystyle\ddot{X}(t)+\omega^{2}\left[X(t)-\varepsilon X^{3}(t)\right]=0, (51)

with initial conditions X(0)=X_{0} and \dot{X}(0)=\omega Y_{0}. Equation (51) is solved order by order with the Ansatz

 \displaystyle X(t)=x^{(0)}(t,\tau)+\varepsilon x^{(1)}(t,\tau)+\mathcal{O}(% \varepsilon^{2}), (52a) where \tau\equiv\varepsilon t is assumed to be an independent time scale such that \displaystyle d_{t}\equiv\partial_{t}+\varepsilon\partial_{\tau}+\mathcal{O}(% \varepsilon^{2}). (52b)

This additional time-scale then allows us to remove the secular term that appears in the \mathcal{O}(\varepsilon) equation. This leads to a renormalization in the oscillation frequency of the \mathcal{O}(1) solution as

 \displaystyle X^{(0)}(t)=x^{(0)}(t,\varepsilon t)=\left[a(0)e^{-i\bar{\omega}t% }+c.c.\right], (53a) \displaystyle\bar{\omega}\equiv\left[1-\frac{3\varepsilon}{2}|a(0)|^{2}\right]\omega, (53b)

where a(0)=(X_{0}+iY_{0})/2.

One may wonder how this leading-order correction is modified in the presence of dissipation. Adding a small damping term \kappa\dot{X}(t) to Eq. (51) such that \kappa\ll\omega requires a new time scale \eta\equiv\frac{\kappa}{\omega}t leading to

 \displaystyle X^{(0)}(t)=e^{-\frac{\kappa}{2}t}\left[a(0)e^{-i\bar{\omega}t}+c% .c.\right], (54a) \displaystyle\bar{\omega}\equiv\left[1-\frac{3\varepsilon}{2}|a(0)|^{2}e^{-% \kappa t}\right]\omega. (54b)

Equations (54a-54b) illustrate a more general fact that the dissipation modifies the frequency renormalization by a decaying envelope. This approach can be extended by introducing higher order (slower) time scales \varepsilon^{2}t,\eta^{2}t, \eta\varepsilon t etc. The lowest order calculation above is valid for times short enough such that \omega t\ll\varepsilon^{-2},\eta^{-2},\eta^{-1}\varepsilon^{-1}.

Besides the extra complexity due to non-commuting algebra of quantum mechanics, the principles of MSPT remain the same in the case of a free quantum Duffing oscillator Bender and Bettencourt (1996). The Heisenberg equation of motion is identical to Eq. (51) where we promote X(t)\to\hat{X}(t). We obtain the \mathcal{O}(1) solution (see App. D.2) as

 \displaystyle\begin{split}\displaystyle\hat{X}^{(0)}(t)=e^{-\frac{\kappa}{2}t}% \left[\frac{\hat{a}(0)e^{-i\hat{\bar{\omega}}t}+e^{-i\hat{\bar{\omega}}t}\hat{% a}(0)}{2\cos\left(\frac{3\omega}{4}\varepsilon te^{-\kappa t}\right)}+H.c.% \right]\end{split} (55a) with an operator-valued renormalization of the frequency \displaystyle\hat{\bar{\omega}}=\left[1-\frac{3\varepsilon}{2}\hat{\mathcal{H}% }(0)e^{-\kappa t}\right]\omega, (55b) \displaystyle\hat{\mathcal{H}}(0)\equiv\frac{1}{2}\left[\hat{a}^{{\dagger}}(0)% \hat{a}(0)+\hat{a}(0)\hat{a}^{{\dagger}}(0)\right]. (55c)

The cosine that appears in the denominator of operator solution (55a) cancels when taking the expectation values with respect to the number basis \{\ket{n}\} of \hat{\mathcal{H}}(0):

 \displaystyle\bra{n-1}\hat{X}^{(0)}(t)\ket{n}=\sqrt{n}e^{-\frac{\kappa}{2}t}e^% {-i\left(1-\frac{3n\varepsilon}{2}e^{-\kappa t}\right)\omega t}. (56)

Having learned from these toy problems, we return to the problem of spontaneous emission which can be mapped into a quantum Duffing oscillator with \varepsilon=\frac{\sqrt{2}}{6}\left(\mathcal{E}_{c}/\mathcal{E}_{j}\right)^{1/2}, up to leading order in perturbation, coupled to multiple leaky quantum harmonic oscillators (see Eq. (38)). We are interested in finding an analytic expression for the shift of the hybridized poles, p_{j} and p_{n}, that appear in the reduced dynamics of the transmon.

The hybridized poles p_{j} and p_{n} are the roots of D_{j}(s) and they are associated with the modal decomposition of the linear theory in Sec. IV.1. The modal decomposition can be found from the linear solution \mathcal{X}_{j}(t) that belongs to the full Hilbert space as

 \displaystyle\begin{split}\displaystyle\hat{\mathcal{X}}_{j}(t)&\displaystyle=% \left(\hat{\mathcal{A}}_{j}e^{p_{j}t}+\sum\limits_{n}\hat{\mathcal{A}}_{n}e^{p% _{n}t}\right)+H.c.\\ &\displaystyle\equiv\left(u_{j}\hat{\bar{a}}_{j}e^{p_{j}t}+\sum\limits_{n}u_{n% }\hat{\bar{a}}_{n}e^{p_{n}t}\right)+H.c.\end{split} (57)

This is the full-Hilbert space version of Eq. (48). It represents the unperturbed solution upon which we are building our perturbation theory. We have used bar-notation to distinguish the creation and annihilation operators in the hybridized mode basis. Furthermore, u_{j} and u_{n} represent the hybridization coefficients, where they determine how much the original transmon operator \hat{\mathcal{X}}_{j}(t), is transmon-like and resonator-like. They can be obtained from a diagonalization of the linear Heisenberg-Langevin equations of motion (see App. D.3). The dependence of u_{j} and u_{n} on coupling \chi_{g} is shown in Fig. (9) for the case where the transmon is infinitesimally detuned below the fundamental mode of the resonator. For \chi_{g}=0, u_{j}=1 and u_{n}=0 as expected. As \chi_{g} reaches \chi_{j}, u_{1} is substantially increased and becomes comparable to u_{j}. By increasing \chi_{g} further, u_{n} for the off-resonant modes start to grow as well.

The nonlinearity acting on the transmon mixes all the unperturbed resonances through self- and cross-Kerr contributions Drummond and Walls (1980); Nigg et al. (2012); Bourassa et al. (2012). Kerr shifts can be measured in a multimode cQED system Rehák et al. (2014); Weißl et al. (2015). We therefore solve for the equations of motion of each mode. These are (see App. D.3)

 \displaystyle\begin{split}&\displaystyle\hat{\ddot{\bar{\mathcal{X}}}}_{l}(t)+% 2\alpha_{l}\hat{\dot{\bar{\mathcal{X}}}}_{l}(t)\\ &\displaystyle+\beta_{l}^{2}\left\{\hat{\bar{\mathcal{X}}}_{l}(t)-\varepsilon_% {l}\left[u_{j}\hat{\bar{\mathcal{X}}}_{j}(t)+\sum\limits_{n}u_{n}\hat{\bar{% \mathcal{X}}}_{n}(t)\right]^{3}\right\}=0,\end{split} (58)

where \hat{\bar{X}}_{l}\equiv\hat{\bar{a}}_{l}+\hat{\bar{a}}_{l}^{{\dagger}} is the quadrature of the lth mode, and \alpha_{l} and \beta_{l} are the decay rate and the oscillation frequency, respectively. Equation (58) is the leading order approximation in the inverse Q-factor of the lth mode, 1/Q_{l}\equiv\alpha_{l}/\beta_{l}. Each hybridized mode has a distinct strength of the nonlinearity \varepsilon_{l}\equiv\frac{\omega_{j}}{\beta_{l}}u_{l}\varepsilon for l\equiv j,n. In order to do MSPT, we need to introduce as many new time-scales as the number of hybridized modes, i.e. \tau_{j}\equiv\varepsilon_{j}t and \tau_{n}\equiv\varepsilon_{n}t, and do a perturbative expansion in all of these time scales. The details of this calculation can be found in App. D.3. Up to lowest order in \varepsilon, we find operator-valued correction of p_{j}=-\alpha_{j}-i\beta_{j} as

 \displaystyle\hat{\bar{p}}_{j}=p_{j}+i\frac{3\varepsilon}{2}\omega_{j}\left[u_% {j}^{4}\hat{\bar{\mathcal{H}}}_{j}(0)e^{-2\alpha_{j}t}+\sum\limits_{n}2u_{j}^{% 2}u_{n}^{2}\hat{\bar{\mathcal{H}}}_{n}(0)e^{-2\alpha_{n}t}\right], (59a) while p_{n}=-\alpha_{n}-i\beta_{n} is corrected as \displaystyle\begin{split}\displaystyle\hat{\bar{p}}_{n}=p_{n}+i\frac{3% \varepsilon}{2}\omega_{j}&\displaystyle\left[u_{n}^{4}\hat{\bar{\mathcal{H}}}_% {n}(0)e^{-2\alpha_{n}t}+2u_{n}^{2}u_{j}^{2}\hat{\bar{\mathcal{H}}}_{j}(0)e^{-2% \alpha_{j}t}\right.\\ &\displaystyle+\left.\sum\limits_{m\neq n}2u_{n}^{2}u_{m}^{2}\hat{\bar{% \mathcal{H}}}_{m}(0)e^{-2\alpha_{m}t}\right],\end{split} (59b) where \hat{\bar{\mathcal{H}}}_{j}(0) and \hat{\bar{\mathcal{H}}}_{n}(0) represent the Hamiltonians of each hybridized mode \displaystyle\hat{\bar{\mathcal{H}}}_{l}(0)\equiv\frac{1}{2}\left[\hat{\bar{a}% }_{l}^{{\dagger}}(0)\hat{\bar{a}}_{l}(0)+\hat{\bar{a}}_{l}(0)\hat{\bar{a}}_{l}% ^{{\dagger}}(0)\right],\ l=j,n. (59c) These are the generalizations of the single quantum Duffing results (55b) and (55c) and reduce to them as \chi_{g}\to 0 where u_{j}=1 and u_{n}=0. Each hybdridized mode is corrected due to a self-Kerr term proportional to u_{l}^{4}, and cross-Kerr terms proportional to u_{l}^{2}u_{l^{\prime}}^{2}. Contributions of the form u_{l}^{2}u_{l^{\prime}}u_{l^{\prime\prime}} Bourassa et al. (2012) do not appear up to the lowest order in MSPT. In terms of Eqs. (59a-59b), the MSPT solution reads
 \displaystyle\begin{split}\displaystyle\hat{\mathcal{X}}_{j}^{(0)}(t)&% \displaystyle=\frac{\hat{\mathcal{A}}_{j}(0)e^{\hat{\bar{p}}_{j}t}+e^{\hat{% \bar{p}}_{j}t}\hat{\mathcal{A}}_{j}(0)}{2\cos\left(\frac{3\omega_{j}}{4}u_{j}^% {4}\varepsilon te^{-2\alpha_{j}t}\right)}+H.c.\\ &\displaystyle+\sum\limits_{n}\left[\frac{\hat{\mathcal{A}}_{n}(0)e^{\hat{\bar% {p}}_{n}t}+e^{\hat{\bar{p}}_{n}t}\hat{\mathcal{A}}_{n}(0)}{2\cos\left(\frac{3% \omega_{j}}{4}u_{n}^{4}\varepsilon te^{-2\alpha_{n}t}\right)}+H.c.\right],\end% {split} (60)

where \hat{\mathcal{A}}_{j,n} is defined in Eq. (57). In Fig. 10, we have compared the Fourier transform of \braket{\hat{\mathcal{X}}_{j}(t)} calculated both for the MSPT solution (60) and the linear solution (48) for initial condition \ket{\Psi(0)}=\frac{\ket{0}_{j}+\ket{1}_{j}}{\sqrt{2}}\otimes\ket{0}_{ph} as a function of \chi_{g}. At \chi_{g}=0, we notice the bare \mathcal{O}(\varepsilon) nonlinear shift of a free Duffing oscillator as predicted by Eq. (53b). As \chi_{g} is increased, the predominantly self-Kerr nonlinearity on the qubit is gradually passed as cross-Kerr contributions to the resonator modes, as observed from the frequency renormalizations (59a) and (59b). As a result of this, interestingly, the effective nonlinear shift in the transmon resonance becomes smaller and saturates at stronger couplings. In other words, the transmon mode becomes more linear at stronger coupling \chi_{g}. This counterintuitive result can be understood from Eq. (59a). For initial condition considered here, the last term in Eq. (59a) vanishes, while one can see from Fig. 9 that u_{j}<1 for \chi_{g}>0.

### IV.3 Numerical simulation of reduced equation

The purpose of this section is to compare the results from MSPT and linear theory to a pure numerical solution valid up to \mathcal{O}(\varepsilon^{2}). A full numerical solution of the Heisenberg equation of motion (29) requires matrix representation of the qubit operator \hat{\mathcal{X}}_{j}(t) over the entire Hilbert space, which is impractical due to the exponentially growing dimension. We are therefore led to work with the reduced Eq. (36). While the nonlinear contribution in Eq. (36) cannot be traced exactly, it is possible to make progress perturbatively. We substitute the perturbative expansion Eq. (38) into Eq. (36):

 \displaystyle\begin{split}&\displaystyle\hat{\ddot{X}}_{j}(t)+\omega_{j}^{2}% \left[1-\gamma+i\mathcal{K}_{1}(0)\right]\left[\hat{X}_{j}(t)-\varepsilon% \operatorname{Tr}_{ph}{\{\hat{\rho}_{ph}(0)\hat{\mathcal{X}}_{j}^{3}(t)\}}% \right]\\ &\displaystyle=-\int_{0}^{t}dt^{\prime}\omega_{j}^{2}\mathcal{K}_{2}(t-t^{% \prime})[\hat{X}_{j}(t^{\prime})-\varepsilon\operatorname{Tr}_{ph}{\{\hat{\rho% }_{ph}(0)\hat{\mathcal{X}}_{j}^{3}(t^{\prime})\}}],\end{split} (61)

with \varepsilon\equiv\frac{\sqrt{2}}{6}\epsilon. If we are interested in the numerical results only up to \mathcal{O}(\varepsilon^{2}) then the cubic term can be replaced as

 \displaystyle\begin{split}\displaystyle\varepsilon\hat{\mathcal{X}}_{j}^{3}(t)% =\varepsilon\left[\left.\hat{\mathcal{X}}_{j}(t)\right|_{\varepsilon=0}\right]% ^{3}+\mathcal{O}\left(\varepsilon^{2}\right).\end{split} (62)

Since we know the linear solution (57) for \hat{\mathcal{X}}_{j}(t) analytically, the trace can be performed directly (see App. E). We obtain the reduced equation in the Hilbert of transmon as

 \displaystyle\begin{split}&\displaystyle\hat{\ddot{X}}_{j}(t)+\omega_{j}^{2}% \left[1-\gamma+i\mathcal{K}_{1}(0)\right]\left[\hat{X}_{j}(t)-\varepsilon\hat{% X}_{j}^{3}(t)\right]\\ &\displaystyle=-\int_{0}^{t}dt^{\prime}\omega_{j}^{2}\mathcal{K}_{2}(t-t^{% \prime})\left[\hat{X}_{j}(t^{\prime})-\varepsilon\hat{X}_{j}^{3}(t^{\prime})% \right]+\mathcal{O}(\varepsilon^{2}).\end{split} (63)

Solving the integro-differential Eq. (63) numerically is a challenging task, since the memory integral on the RHS requires the knowledge of all results for t^{\prime}<t. Therefore, simulation time for Eq. (63) grows polynomially with t. The beauty of the Laplace transform in the linear case is that it turns a memory contribution into an algebraic form. However, it is inapplicable to Eq. (63).

In Fig. 11, we compared the numerical results to both linear and MSPT solutions up to 10 resonator round-trip times and for different values of \chi_{g}. For \chi_{g}=0, the transmon is decoupled and behaves as a free Duffing oscillator. This corresponds to the first row in Fig. ((a)a) where there is only one frequency component and MSPT provides the correction given in Eq. (55b). As we observe in Fig. (a)a the MSPT results lie on top of the numerics, while the linear solution shows a visible lag by the 10th round-trip. Increasing \chi_{g} further, brings more frequency components into play. As we observe in Fig. 10, for \chi_{g}=0.01 the most resonant mode of the resonator has a non-negligible u_{1}. Therefore, we expect to observe weak beating in the dynamics between this mode and the dominant transmon-like resonance, which is shown in Fig. (b)b. Figures (c)c and (d)d show stronger couplings where many resonator modes are active and a more complex beating is observed. In all these cases, the MSPT results follow the pure numerical results more closely than the linear solution confirming the improvement provided by perturbation theory.

### IV.4 System output

Up to this point, we studied the dynamics of the spontaneous emission problem in terms of one of the quadratures of the transmon qubit, i.e. \braket{\hat{\mathcal{X}}_{j}(t)}. In a typical experimental setup however, the measuarable quantities are the quadratures of the field outside the resonator Clerk et al. (2010). We devote this section to the computation of these quantities.

The expression of the fields \hat{\varphi}(x,t) can be directly inferred from the solution of the inhomogeneous wave Eq. (26) using the impulse response (GF) defined in Eq. 82. We note that this holds irrespective of whether one is solving for the classical or as is the case here, for the quantum fields. Taking the expectation value of this solution (App. B.4) with respect to the initial density matrix (33) we find

 \displaystyle\braket{\hat{\varphi}(x,t)}=\chi_{s}\omega_{j}^{2}\int_{0}^{t}dt^% {\prime}G(x,t|x_{0},t^{\prime})\braket{\sin[\hat{\varphi}_{j}(t^{\prime})]}. (64)

Dividing both sides by \phi_{\text{zpf}} and keeping the lowest order we obtain the resonator response as

 \displaystyle\braket{\hat{\mathcal{X}}^{(0)}(x,t)}=\chi_{s}\omega_{j}^{2}\int_% {0}^{t}dt^{\prime}G(x,t|x_{0},t^{\prime})\braket{\hat{\mathcal{X}}_{j}^{(0)}(t% ^{\prime})}, (65)

where \hat{\mathcal{X}}_{j}^{(0)}(t) is the lowest order MSPT solution (60), which takes into account the frequency correction to \mathcal{O}(\varepsilon). Taking the Laplace transform decouples the convolution

 \displaystyle\braket{\hat{\tilde{\mathcal{X}}}^{(0)}(x,s)}=\chi_{s}\omega_{j}^% {2}\tilde{G}(x,x_{0},s)\braket{\hat{\tilde{\mathcal{X}}}_{j}^{(0)}(s)}, (66)

which indicates that the resonator response is filtered by the GF.

Figure 12 shows the field outside the right end of the resonator, \braket{\hat{\tilde{\mathcal{X}}}(x=1^{+},s=i\omega)}, in both linear and lowest order MSPT approximations. This quadrature can be measured via heterodyne detection Bishop et al. (2009). Note that the hybridized resonances are the same as those of \braket{\hat{\mathcal{X}}_{j}(t)} shown in Fig. 10. What changes is the relative strength of the residues. The GF has poles at the bare cavity resonances and therefore the more hybridized a pole is, the smaller its residue becomes.

## V Conclusion

In this paper, we introduced a new approach for studying the effective non-Markovian Heisenberg equation of motion of a transmon qubit coupled to an open multimode resonator beyond rotating wave and two level approximations. The main motivation to go beyond a two level representation lies in the fact that a transmon is a weakly nonlinear oscillator. Furthermore, the information regarding the electromagnetic environment is encoded in a single function, i.e. the electromagnetic GF. As a result, the opening of the resonator is taken into account analytically, in contrast to the Lindblad formalism where the decay rates enter only phenomenologically.

We applied this theory to the problem of spontaneous emission as the simplest possible example. The weak nonlinearity of the transmon allowed us to solve for the dynamics perturbatively in terms of (\mathcal{E}_{c}/\mathcal{E}_{j})^{1/2} which appears as a measure of nonlinearity. Neglecting the nonlinearity, the transmon acts as a simple harmonic oscillator and the resulting linear theory is exactly solvable via Laplace transform. By employing Laplace transform, we avoided Markov approximation and therefore accounted for the exact hybridization of transmon and resonator resonances. Up to leading nonzero order, the transmon acts as a quantum Duffing oscillator. Due to the hybridization, the nonlinearity of the transmon introduces both self-Kerr and cross-Kerr corrections to all hybdridized modes of the linear theory. Using MSPT, we were able to obtain closed form solutions in Heisenberg picture that do not suffer from secular behavior. A direct numerical solution confirmed the improvement provided by the perturbation theory over the harmonic theory. Surprisingly, we also learned that the linear theory becomes more accurate for stronger coupling since the nonlinearity is suppressed in the qubit-like resonance due to being shared between many hybdridized modes.

The theory developed here illustrates how far one can go without the concept of photons. Many phenomena in the domain of quantum electrodynamics, such as spontaneous or stimulated emission and resonance fluorescence, have accurate semiclassical explanations in which the electric field is treated classically while the atoms obey the laws of quantum mechanics. For instance, the rate of spontaneous emission can be related to the local density of electromagnetic modes in the weak coupling limit. While it is now well understood that the electromagnetic fluctuations are necessary to start the spontaneous emission process Haake et al. (1981), it is important to ask to what extent a quantized electromagnetic field effects the qubit dynamics Scully and Sargent (1972). We find here that although the electromagnetic degrees of freedom are integrated out and the dynamics can systematically be reduced to the Hilbert space of the transmon, the quantum state of the electromagnetic environment reappears in the initial and boundary conditions when computing observables.

Although we studied only the spontaneous emission problem in terms of quadratures, our theory can be applied to a driven-dissipative problem as well and all the mathematical machinery developed in this work can be used in more generic situations. In order to maintain a reasonable amount of material in this paper, we postpone the results of the driven-dissipative problem, as well as the study of correlation functions to future work.

## VI Acknowledgements

We appreciate helpful discussions with O. Malik on implementing the numerical results of Sec. IV.3. This research was supported by the US Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Award No. DE-SC0016011.

## Appendix A Quantum equations of motion

The classical Lagrangian for the system shown in Fig. 3 can be found as sum of the Lagrangians for each circuit element. In the following, we use the convention of working with flux variables Bishop (2010); Devoret et al. (2014) as the generalized coordinate for our system. For an arbitrary node n in the circuit, the flux variable \Phi_{n}(t) is defined as

 \displaystyle\Phi_{n}(t)\equiv\int_{0}^{t}dt^{\prime}V_{n}(t^{\prime}), (67)

while V_{n}(t^{\prime}) stands for the voltage at node n.

The classical Euler-Lagrange equations of motion can then be found by setting the variation of Lagrangian with respect to each flux variable to zero. For the transmon and the resonator we find

 \displaystyle\ddot{\Phi}_{j}+\frac{1}{C_{g}+C_{j}}\frac{\partial U_{j}(\Phi_{j% })}{\partial\Phi_{j}}=\gamma\partial_{t}^{2}\Phi(x_{0},t), (68)
 \displaystyle\partial_{x}^{2}\Phi(x,t)-lc(x,x_{0})\partial_{t}^{2}\Phi(x,t)=l% \gamma\delta(x-x_{0})\frac{\partial U_{j}(\Phi_{j})}{\partial\Phi_{j}}. (69)

where U_{j}(\Phi_{j}) stands for the Josephson potential as

 \displaystyle U_{j}(\Phi_{j})=-E_{j}\cos{\left(\frac{2\pi}{\Phi_{0}}\Phi_{j}% \right)}, (70)

and \Phi_{0}\equiv\frac{h}{2e} is the superconducting flux quantum. Furthermore, 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}). Moreover, l and c are the inductance and capacitance per length of the resonator and waveguides while c(x,x_{0})\equiv c+C_{s}\delta(x-x_{0}) represents the modified capacitance per length due to coupling to transmon.

In addition, we find two wave equations for the flux field of the left and right waveguides as

 \displaystyle\partial_{x}^{2}\Phi_{R,L}(x,t)-lc\partial_{t}^{2}\Phi_{R,L}(x,t)% =0, (71)

The boundary conditions (BC) are derived from continuity of current at each end as

 \displaystyle\begin{split}\displaystyle-\frac{1}{l}\left.\partial_{x}\Phi% \right|_{x=L^{-}}&\displaystyle=-\frac{1}{l}\left.\partial_{x}\Phi_{R}\right|_% {x=L^{+}}\\ &\displaystyle=C_{R}\partial_{t}^{2}\left[\Phi(L^{-},t)-\hat{\Phi}_{R}(L^{+},t% )\right],\end{split} (72a) \displaystyle-\frac{1}{l}\left.\partial_{x}\Phi\right|_{x=0^{+}} \displaystyle=-\frac{1}{l}\left.\partial_{x}\Phi_{L}\right|_{x=0^{-}} (72b) \displaystyle=C_{L}\partial_{t}^{2}\left[\Phi_{L}(0^{-},t)-\Phi(0^{+},t)\right], (72c)

continuity of flux at x=x_{0}

 \displaystyle\Phi(x=x_{0}^{-},t)=\Phi(x=x_{0}^{+},t), (73)

and conservation of current at x=x_{0} as

 \displaystyle\left.\partial_{x}\Phi\right|_{x=x_{0}^{+}}-\left.\partial_{x}% \Phi\right|_{x=x_{0}^{-}}-lC_{s}\partial_{t}^{2}\Phi(x_{0},t)=l\gamma\frac{% \partial U_{j}(\Phi_{j})}{\partial\Phi_{j}}. (74)

In order to find the quantum equations of motion, we follow the common procedure of canonical quantization Devoret et al. (2014):

• Find the conjugate momenta Q_{n}\equiv\frac{\delta\mathcal{L}}{\delta\dot{\Phi}_{n}}

• Find the classical Hamiltonian via a Legendre transformation as \mathcal{H}=\sum\limits_{n}Q_{n}\dot{\Phi_{n}}-\mathcal{L}

• Find the Hamiltonian operator by promoting the classical conjugate variables to quantum operators such that \{\hat{\Phi}_{m},\hat{Q}_{n}\}=\delta_{mn}\to[\hat{\Phi}_{m},\hat{Q}_{n}]=i% \hbar\delta_{mn}. We use a hat-notation to distinguish operators from classical variables.

The derivation for the quantum Hamiltonian of the the closed version of this system where C_{R,L}\to 0 can be found in Malekakhlagh and Türeci (2016) (see App. A, B, and C). Note that nonzero end capacitors C_{R,L} leave the equations of motion for the resonator and waveguides unchanged, but modify the BC of the problem at x=0,L. The resulting equations of motion for the quantum flux operators \hat{\Phi}_{j}, \hat{\Phi}(x,t) and \hat{\Phi}_{R,L}(x,t) have the exact same form as the classical Euler-Lagrange equations of motion.

Next, we define unitless parameters and variables as

where v_{p}\equiv 1/\sqrt{lc} is the phase velocity of the resonator and waveguides. Furthermore, we define unitless capacitances as

as well as a unitless modified capacitance per length as

 \displaystyle\chi(\bar{x},\bar{x}_{0})\equiv 1+\chi_{s}\delta(\bar{x}-\bar{x}_% {0}). (77)

Then, the unitless equations of motion for our system are found as

 \displaystyle\hat{\ddot{\varphi}}_{j}(\bar{t})+(1-\gamma)\bar{\omega}_{j}^{2}% \sin{[\hat{\varphi}_{j}(\bar{t})]}=\gamma\partial_{\bar{t}}^{2}\hat{\varphi}(% \bar{x}_{0},\bar{t}), (78a) \displaystyle\begin{split}\displaystyle\left[\partial_{\bar{x}}^{2}-\chi(\bar{% x},\bar{x}_{0})\partial_{\bar{t}}^{2}\right]\hat{\varphi}(\bar{x},\bar{t})=% \chi_{s}\bar{\omega}_{j}^{2}\sin{[\varphi_{j}(\bar{t})]}\delta(\bar{x}-\bar{x}% _{0}),\end{split} (78b) \displaystyle\partial_{\bar{x}}^{2}\hat{\varphi}_{R,L}(\bar{x},\bar{t})-% \partial_{\bar{t}}^{2}\hat{\varphi}_{R,L}(\bar{x},\bar{t})=0, (78c)

with the unitless BCs given as

 \displaystyle\begin{split}\displaystyle-\left.\partial_{\bar{x}}\hat{\varphi}% \right|_{\bar{x}=1^{-}}&\displaystyle=-\left.\partial_{\bar{x}}\hat{\varphi}_{% R}\right|_{\bar{x}=1^{+}}\\ &\displaystyle=\chi_{R}\partial_{\bar{t}}^{2}\left[\hat{\varphi}(1^{-},\bar{t}% )-\hat{\varphi}_{R}(1^{+},\bar{t})\right],\end{split} (79a) \displaystyle\begin{split}\displaystyle-\left.\partial_{\bar{x}}\hat{\varphi}% \right|_{\bar{x}=0^{+}}&\displaystyle=-\left.\partial_{\bar{x}}\hat{\varphi}_{% L}\right|_{\bar{x}=0^{-}}\\ &\displaystyle=\chi_{L}\partial_{\bar{t}}^{2}\left[\hat{\varphi}_{L}(0^{-},% \bar{t})-\hat{\varphi}(0^{+},\bar{t})\right],\end{split} (79b) \displaystyle\hat{\varphi}(\bar{x}=\bar{x}_{0}^{-},\bar{t})=\hat{\varphi}(\bar% {x}=\bar{x}_{0}^{+},\bar{t}), (79c) \displaystyle\begin{split}\displaystyle\left.\partial_{\bar{x}}\hat{\varphi}% \right|_{\bar{x}=\bar{x}_{0}^{+}}&\displaystyle-\left.\partial_{\bar{x}}\hat{% \varphi}\right|_{\bar{x}=\bar{x}_{0}^{-}}-\chi_{s}\partial_{\bar{t}}^{2}\hat{% \varphi}(\bar{x}_{0},\bar{t})\\ &\displaystyle=\chi_{s}\bar{\omega}_{j}^{2}\sin{[\varphi_{j}(\bar{t})]}.\end{split} (79d)

In Eqs. (78a) and (78b), we have defined the unitless oscillation frequency \bar{\omega}_{j} as

 \displaystyle\bar{\omega}_{j}^{2}\equiv lcL^{2}\frac{E_{j}}{C_{j}}\left(\frac{% 2\pi}{\Phi_{0}}\right)^{2}=8\mathcal{E}_{c}\mathcal{E}_{j}, (80)

where \mathcal{E}_{c} and \mathcal{E}_{j} stand for the unitless charging and Josephson energy given as

 \displaystyle\mathcal{E}_{j,c}\equiv\sqrt{lc}L\frac{E_{j,c}}{\hbar}, (81)

with E_{c}\equiv\frac{e^{2}}{2C_{j}}.

In what follows, we work with the unitless Eqs. (78a-78c) and BCs (79a-79d) and drop the bars.

## Appendix B Effective dynamics of the transmon via a Heisenberg picture Green’s function method

In order to find the effective dynamics of the transmon qubit, one has to solve for the flux field \hat{\varphi}(x,t) and substitute the result back into the RHS of time evolution of the qubit given by Eq. (78a). It is possible to perform this procedure in terms of the resonator GF. In Sec. B.1 we define the resonator GF. In Sec. B.3 we study the spectral representation of the GF in terms of a suitable set of non-Hermitian modes. In Sec. B.4, we discuss the derivation of the effective dynamics of transmon in terms of the resonator GF. Finally, in Secs. B.5 and B.6 we discuss how the generic dynamics is reduced for the problem of spontaneous emission.

### B.1 Definition of G(x,t|x^{\prime},t^{\prime})

The resonator GF is defined as the response of the linear system of Eqs. (78b-78c) to a \delta-function source in space-time as

 \displaystyle\begin{split}\displaystyle\left[\partial_{x}^{2}-\chi(x,x_{0})% \partial_{t}^{2}\right]&\displaystyle G(x,t|x_{0},t_{0})=\delta(x-x_{0})\delta% (t-t_{0}),\end{split} (82)

with the same BCs as Eqs. (79a-79d). Using the Fourier transform conventions

 \displaystyle\tilde{G}(x,x_{0},\omega)=\int_{-\infty}^{\infty}dtG(x,t|x_{0},t_% {0})e^{+i\omega(t-t_{0})}, (83a) \displaystyle 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})}, (83b)

Eq. (82) transforms into a Helmholtz equation

 \displaystyle\begin{split}\displaystyle\left[\partial_{x}^{2}+\omega^{2}\chi(x% ,x_{0})\right]\tilde{G}(x,x_{0},\omega)=\delta(x-x_{0}).\end{split} (84)

Moreover, the BCs are transformed by replacing \partial_{x}\to\partial_{x} and \partial_{t}\to-i\omega as

 \displaystyle\left.\tilde{G}\right|_{x=x_{0}^{+}}=\left.\tilde{G}\right|_{x=x_% {0}^{-}}, (85a) \displaystyle\left.\partial_{x}\tilde{G}\right|_{x=x_{0}^{+}}-\left.\partial_{% x}\tilde{G}\right|_{x=x_{0}^{-}}+\chi_{s}\omega^{2}\left.\tilde{G}\right|_{x=x% _{0}}=1, (85b) \displaystyle\begin{split}\displaystyle\left.\partial_{x}\tilde{G}\right|_{x=1% ^{-}}&\displaystyle=\left.\partial_{x}\tilde{G}\right|_{x=1^{+}}\\ &\displaystyle=\chi_{R}\omega^{2}\left(\left.\tilde{G}\right|_{x=1^{-}}-\left.% \tilde{G}\right|_{x=1^{+}}\right),\end{split} (85c) \displaystyle\begin{split}\displaystyle\left.\partial_{x}\tilde{G}\right|_{x=0% ^{-}}&\displaystyle=\left.\partial_{x}\tilde{G}\right|_{x=0^{+}}\\ &\displaystyle=\chi_{L}\omega^{2}\left(\left.\tilde{G}\right|_{x=0^{-}}-\left.% \tilde{G}\right|_{x=0^{+}}\right).\end{split} (85d)

Note that BCs (85a-85d) do not specify what happens to \tilde{G}(x,x_{0},\omega) at x\to\pm\infty. We model the baths by imposing outgoing BCs at infinity as

 \displaystyle\left.\partial_{x}\tilde{G}(x,x_{0},\omega)\right|_{x\to\pm\infty% }=\pm i\omega\tilde{G}(x\to\pm\infty,x_{0},\omega), (86)

which precludes any reflections from the waveguides to the resonator.

### B.2 Spectral representation of GF for a closed resonator

It is helpful to revisit spectral representation of GF for the closed version of our system by setting \chi_{R}=\chi_{L}=0. This imposes Neumann BC \partial_{x}\tilde{G}|_{x=0,1}=0 and the resulting differential operator becomes Hermitian. The idea of spectral representation is to expand \tilde{G} in terms of a discrete set of normal modes that obey the homogeneous wave equation

 \displaystyle\partial_{x}^{2}\tilde{\Phi}_{n}(x)+\chi(x,x_{0})\omega_{n}^{2}% \tilde{\Phi}_{n}(x)=0, (87a) \displaystyle\left.\partial_{x}\tilde{\Phi}_{n}(x)\right|_{x=0,1}=0. (87b)

Then, the real valued eigenfrequencies obey the transcendental equation

 \displaystyle\sin{(\omega_{n})}+\chi_{s}\omega_{n}\cos{(\omega_{n}x_{0})}\cos{% \left[\omega_{n}(1-x_{0})\right]}=0. (88)

 \displaystyle\tilde{\Phi}_{n}(x)\propto\begin{cases}\cos{\left[\omega_{n}(1-x_% {0})\right]}\cos{(\omega_{n}x)},&0

where the normalization is fixed by the orthogonality condition

 \displaystyle\int_{0}^{1}dx\chi(x,x_{0})\tilde{\Phi}_{m}(x)\tilde{\Phi}_{n}(x)% =\delta_{mn}. (90)

Note that eigenfunctions of a Hermitian differential operator form a complete orthonormal basis. This allows us to deduce the spectral representation of \tilde{G}(x,x^{\prime},\omega) Morse and Feshbach (1953); Economou (1984); Hassani (2013) as

 \displaystyle\tilde{G}(x,x^{\prime},\omega)=\sum\limits_{n\in\mathbb{N}}\frac{% \tilde{\Phi}_{n}(x)\tilde{\Phi}_{n}(x^{\prime})}{\omega^{2}-\omega_{n}^{2}}=% \sum\limits_{n\in\mathbb{Z}\atop n\neq 0}\frac{1}{2\omega}\frac{\tilde{\Phi}_{% n}(x)\tilde{\Phi}_{n}(x^{\prime})}{\omega-\omega_{n}}, (91)

where the second representation is written due to relations \omega_{-n}=-\omega_{n} and \tilde{\Phi}_{-n}(x)=\tilde{\Phi}_{n}(x).

### B.3 Spectral representation of GF for an open resonator

A spectral representation can also be found for the GF of an open resonator in terms of a discrete set of non-Hermitian modes that carry a constant flux away from the resonator. The Constant Flux (CF) modes Türeci et al. (2006) have allowed a consistent formulation of the semiclassical laser theory for complex media such as random lasers Türeci et al. (2008). The non-Hermiticity originates from the fact that the waveguides are assumed to be infinitely long, hence no radiation that is emitted from the resonator to the waveguides can be reflected back. This results in discrete and complex-valued poles of the GF. The CF modes satisfy the same homogeneous wave equation

 \displaystyle\partial_{x}^{2}\tilde{\Phi}_{n}(x,\omega)+\chi(x,x_{0})\omega_{n% }^{2}(\omega)\tilde{\Phi}_{n}(x,\omega)=0, (92)

but with open BCs the same as Eqs. (85a-86). Note that the resulting CF modes \tilde{\Phi}_{n}(x,\omega) and eigenfrequencies \omega_{n}(\omega) parametrically depend on the source frequency \omega.

Considering only an outgoing plane wave solution for the left and right waveguides based on (86), the general solution for \tilde{\Phi}_{n}(x,\omega) reads

 \displaystyle\tilde{\Phi}_{n}(x,\omega)=\begin{cases}A_{n}^{<}e^{i\omega_{n}(% \omega)x}+B_{n}^{<}e^{-i\omega_{n}(\omega)x},&0}e^{i\omega_{n}(\omega)x}+B_{n}^{>}e^{-i\omega_{n}(\omega)x},&x_{0}1\\ D_{n}e^{-i\omega x},&x<0\\ \end{cases} (93)

Applying BCs (85a-85d) leads to a characteristic equation

 \displaystyle\begin{split}&\displaystyle\sin\left[\omega_{n}(\omega)\right]+(% \chi_{R}+\chi_{L})\omega_{n}(\omega)\left\{\cos[\omega_{n}(\omega)]-\frac{% \omega_{n}(\omega)}{\omega}\sin[\omega_{n}(\omega)]\right\}\\ &\displaystyle-\chi_{R}\chi_{L}\omega_{n}^{2}(\omega)\left\{2i\frac{\omega_{n}% (\omega)}{\omega}\cos[\omega_{n}(\omega)]+\left[1+\frac{\omega_{n}^{2}(\omega)% }{\omega^{2}}\right]\sin[\omega_{n}(\omega)]\right\}\\ &\displaystyle+\chi_{s}\omega_{n}(\omega)\left\{\cos[\omega_{n}(\omega)x_{0}]-% \chi_{L}\frac{\omega_{n}(\omega)}{\omega}\left\{i\omega_{n}(\omega)\cos[\omega% _{n}(\omega)x_{0}]+\omega\sin[\omega_{n}(\omega)x_{0}]\right\}\right\}\\ &\displaystyle\times\left\{\cos[\omega_{n}(\omega)(1-x_{0})]-\chi_{R}\frac{% \omega_{n}(\omega)}{\omega}\left\{i\omega_{n}(\omega)\cos[\omega_{n}(\omega)(1% -x_{0})]+\omega\sin[\omega_{n}(\omega)(1-x_{0})]\right\}\right\}=0,\end{split} (94)

which gives the parametric dependence of CF frequencies on \omega. Then, the CF modes \tilde{\Phi}_{n}(x,\omega) are calculated as

 \displaystyle\tilde{\Phi}_{n}(x,\omega)\propto\begin{cases}e^{-i\omega_{n}(% \omega)(x-x_{0}+1)}\left[e^{2i\omega_{n}(\omega)x}+(1-2i\omega_{n}(\omega)\chi% _{L})\right]\left[e^{2i\omega_{n}(\omega)(1-x_{0})}+\left(1-2i\omega_{n}(% \omega)\chi_{R}\right)\right],&01\\ -2i\chi_{L}\omega_{n}(\omega)e^{-i\omega_{n}(\omega)(1-x_{0})}\left[e^{2i% \omega_{n}(\omega)(1-x_{0})}+(1-2i\chi_{R}\omega_{n}(\omega))\right]e^{-i% \omega x}.&x<0\\ \end{cases} (95)

These modes satisfy the biorthonormality condition

 \displaystyle\begin{split}\displaystyle\int_{0}^{1}dx\chi(x,x_{0})\bar{\tilde{% \Phi}}_{m}^{*}(x,\omega)\tilde{\Phi}_{n}(x,\omega)=\delta_{mn},\end{split} (96)

where \{\bar{\tilde{\Phi}}_{m}(x,\omega)\} satisfy the Hermitian adjoint of eigenvalue problem (92). In other words, \tilde{\Phi}_{n}(x,\omega) and \bar{\tilde{\Phi}}_{n}(x,\omega) are the right and left eigenfunctions and obey \bar{\tilde{\Phi}}_{n}(x,\omega)=\tilde{\Phi}_{n}^{*}(x,\omega). The normalization of Eq. (95) is then fixed by setting m=n.

In terms of the CF modes, the spectral representation of the GF can then be constructed

 \displaystyle\tilde{G}(x,x^{\prime},\omega)=\sum\limits_{n}\frac{\tilde{\Phi}_% {n}(x,\omega)\bar{\tilde{\Phi}}_{n}^{*}(x^{\prime},\omega)}{\omega^{2}-\omega_% {n}^{2}(\omega)}. (97)

Examining Eq. (97), we realize that there are two sets of poles of \tilde{G}(x,x^{\prime},\omega) in the complex \omega plane. First, from setting the denominator of Eq. (97) to zero which gives \omega=\omega_{n}(\omega). These are the quasi-bound eigenfrequencies that satisfy the transcendental characteristic equation

 \displaystyle\begin{split}&\displaystyle\left[e^{2i\omega_{n}}-(1-2i\chi_{L}% \omega_{n})(1-2i\chi_{R}\omega_{n})\right]\\ &\displaystyle+\frac{i}{2}\chi_{s}\omega_{n}[e^{2i\omega_{n}x_{0}}+(1-2i\chi_{% L}\omega_{n})]\\ &\displaystyle\times[e^{2i\omega_{n}(1-x_{0})}+(1-2i\chi_{R}\omega_{n})]=0.% \end{split} (98)

The quasi bound solutions \omega_{n} to Eq. (98) reside in the lower half of complex \omega-plane and come in symmetric pairs with respect to the \text{Im}\{\omega\} axis, i.e. both \omega_{n} and -\omega_{n}^{*} satisfy the transcendental Eq. (98). Therefore, we can label the eigenfrequencies as

 \displaystyle\omega_{n}=\begin{cases}-i\kappa_{0},&n=0\\ +\nu_{n}-i\kappa_{n},&n\in+\mathbb{N}\\ -\nu_{n}-i\kappa_{n},&n\in-\mathbb{N}\end{cases} (99)

where \nu_{n} and \kappa_{n} are positive quantities representing the oscillation frequency and decay rate of each quasi-bound mode. Second, there is an extra pole at \omega=0 which comes from the \omega-dependence of CF states \tilde{\Phi}_{n}(x,\omega). We confirmed these poles by solving for the explicit solution \tilde{G}(x,x^{\prime},\omega) that obeys Eq. (84) with BCs (85a-86) with Mathematica.

### B.4 Effective dynamics of transmon qubit

Note that Eqs. (78b-78c) are linear in terms of \hat{\varphi}(x,t) and \hat{\varphi}_{R,L}(x,t) . Therefore, it is possible to eliminate these linear degrees of freedom and express the formal solution for \hat{\varphi}(x,t) in terms of \hat{\varphi}_{j}(t) and G(x,t|x^{\prime},t^{\prime}). At last, by plugging the result into the RHS of Eq. (78b) we find a closed equation for \hat{\varphi}_{j}(t).

Let us denote the source term that appears on the RHS of Eq. (78b) as

 \displaystyle S\left[\hat{\varphi}_{j}(t)\right]\equiv\chi_{s}\omega_{j}^{2}% \sin{[\hat{\varphi}_{j}(t)]}. (100)

Then, we write two equations for \hat{\varphi}(x,t) and G(x,t|x^{\prime},t^{\prime}) Morse and Feshbach (1953) (See Sec. 7.3) as

 \displaystyle\left[\partial_{x^{\prime}}^{2}-\chi(x^{\prime},x_{0})\partial_{t% ^{\prime}}^{2}\right]\hat{\varphi}(x^{\prime},t^{\prime})=S\left[\hat{\varphi}% _{j}(t^{\prime})\right]\delta(x^{\prime}-x_{0}), (101a) \displaystyle\left[\partial_{x^{\prime}}^{2}-\chi(x,x^{\prime})\partial_{t^{% \prime}}^{2}\right]G(x,t|x^{\prime},t^{\prime})=\delta(x-x^{\prime})\delta(t-t% ^{\prime}). (101b)

In Eq. (101b) we have employed the reciprocity property of the GF

 \displaystyle G(x,t|x^{\prime},t^{\prime})=G(x^{\prime},-t^{\prime}|x,-t), (102)

which holds since Eq. (101b) is invariant under

Multiplying Eq. (101a) by G(x,t|x^{\prime},t^{\prime}) and Eq. (101b) by \hat{\varphi}(x^{\prime},t^{\prime}) and integrating over the dummy variable x^{\prime} in the interval [0^{-},1^{+}] and over t^{\prime} in the interval [0,t^{+}] and finally taking the difference gives

 \displaystyle\begin{split}&\displaystyle\int_{0}^{t^{+}}dt^{\prime}\int_{0^{-}% }^{1^{+}}dx^{\prime}\left\{\underbrace{\left(G\partial_{x^{\prime}}^{2}\hat{% \varphi}-\hat{\varphi}\partial_{x^{\prime}}^{2}G\right)}_{\bf(a)}\right.\\ &\displaystyle+\underbrace{\left[\chi(x,x^{\prime})\hat{\varphi}\partial_{t^{% \prime}}^{2}G-\chi(x^{\prime},x_{0})G\partial_{t^{\prime}}^{2}\hat{\varphi}% \right]}_{\bf(b)}\\ &\displaystyle-\left.\underbrace{GS(\hat{\varphi}_{j})\delta(x^{\prime}-x_{0})% }_{\bf(c)}+\underbrace{\hat{\varphi}\delta(t-t^{\prime})\delta(x-x^{\prime})}_% {\bf(d)}\right\}=0,\end{split} (104)

where we have used the shorthand notation G\equiv G(x,t|x^{\prime},t^{\prime}) and \hat{\varphi}\equiv\hat{\varphi}(x^{\prime},t^{\prime}).

The term labeled as (a) can be simplified further through integration by parts in x^{\prime} as

 \displaystyle\int_{0}^{t^{+}}dt^{\prime}\left.\left(G\partial_{x^{\prime}}\hat% {\varphi}-\hat{\varphi}\partial_{x^{\prime}}G\right)\right|_{x^{\prime}=0^{-}}% ^{x^{\prime}=1^{+}} (105)

There are two contributions from term (b). One comes from the constant capacitance per length in \chi(x,x^{\prime}) and \chi(x,x_{0}) that simplifies to

 \displaystyle\int_{0^{-}}^{1^{+}}\,dx^{\prime}\left.\left(\hat{\varphi}% \partial_{t^{\prime}}G-G\partial_{t^{\prime}}\hat{\varphi}\right)\right|_{t^{% \prime}=0}, (106)

where due to working with the retarded GF

 \displaystyle G(x,t|x^{\prime},t^{+})=0, (107)

hence the upper limit t^{\prime}=t^{+} vanishes. The second contribution comes from the Dirac \delta-functions in \chi(x,x^{\prime}) and \chi(x,x_{0}) which gives

 \displaystyle\begin{split}\displaystyle\chi_{s}\int_{0}^{t^{+}}dt^{\prime}&% \displaystyle\left[\hat{\varphi}(x,t^{\prime})\partial_{t^{\prime}}^{2}G(x,t|x% ,t^{\prime})\right.\\ &\displaystyle\left.-G(x,t|x_{0},t^{\prime})\partial_{t^{\prime}}^{2}\hat{% \varphi}(x_{0},t^{\prime})\right]\end{split} (108)

Terms (c) and (d) get simplified due to Dirac \delta-functions as

 \displaystyle\int_{0}^{t^{+}}\,dt^{\prime}G(x,t|x_{0},t^{\prime})S[\hat{% \varphi}_{j}(t^{\prime})], (109)

and \hat{\varphi}(x,t), respectively.

At the end, we find a generic solution for the flux field \hat{\varphi}(x,t) in the domain [0^{-},1^{+}] as

 \displaystyle\begin{split}&\displaystyle\hat{\varphi}(x,t)=\underbrace{\int_{0% }^{t^{+}}\,dt^{\prime}G(x,t|x_{0},t^{\prime})S[\hat{\varphi}_{j}(t^{\prime})]}% _{Source\ Contribution}+\underbrace{\int_{0}^{t^{+}}\,dt^{\prime}\left.\left[% \hat{\varphi}(x^{\prime},t^{\prime})\partial_{x^{\prime}}G(x,t|x^{\prime},t^{% \prime})-G(x,t|x^{\prime},t^{\prime})\partial_{x^{\prime}}\hat{\varphi}(x^{% \prime},t^{\prime})\right]\right|_{x^{\prime}=0^{-}}^{x^{\prime}=1^{+}}}_{% Boundary\ Contribution}\\ &\displaystyle+\underbrace{\int_{0^{-}}^{1^{+}}\,dx^{\prime}\left.\left[\hat{% \varphi}(x^{\prime},t^{\prime})\partial_{t^{\prime}}G(x,t|x^{\prime},t^{\prime% })-G(x,t|x^{\prime},t^{\prime})\partial_{t^{\prime}}\hat{\varphi}(x^{\prime},t% ^{\prime})\right]\right|_{t^{\prime}=0}}_{Initial\ Condition\ Contribution}\\ &\displaystyle+\underbrace{\chi_{s}\int_{0}^{t^{+}}dt^{\prime}\left[\hat{% \varphi}(x,t^{\prime})\partial_{t^{\prime}}^{2}G(x,t|x,t^{\prime})-G(x,t|x_{0}% ,t^{\prime})\partial_{t^{\prime}}^{2}\hat{\varphi}(x_{0},t^{\prime})\right]}_{% Feedback\ induced\ by\ transmon}.\end{split} (110)

According to Eq. (78a), the transmon is forced by the resonator flux field evaluated at x=x_{0}, i.e. \hat{\varphi}(x_{0},t). In the following, we rewrite the GF in terms of its Fourier representation for each term in Eq. (110) at x=x_{0}. The Fourier representation simplifies the boundary contribution further, while also allowing us to employ the spectral representation of GF discussed in Sec. B.3.

The source contribution can be written as

 \displaystyle\chi_{s}\int_{0}^{t}\,dt^{\prime}\int_{-\infty}^{+\infty}\frac{d% \omega}{2\pi}\tilde{G}(x_{0},x_{0},\omega)\omega_{j}^{2}\sin{[\hat{\varphi}_{j% }(t^{\prime})]}e^{-i\omega(t-t^{\prime})}. (111)

The boundary terms consist of two separate contributions at each end. Assuming that there is no radiation in the waveguides for t<0 we can write

 \displaystyle\hat{\varphi}_{R,L}(x,t)=\hat{\varphi}_{R,L}(x,t)\Theta(t), (112a) \displaystyle\partial_{x}\hat{\varphi}_{R,L}(x,t)=\partial_{x}\hat{\varphi}_{R% ,L}(x,t)\Theta(t). (112b)

Using Eqs. (112a-112b) and causality of the GF, i.e. G(x,t|x^{\prime},t^{\prime})\propto\Theta(t-t^{\prime}), we can extend the integration domain in t^{\prime} from [0,t^{+}] to [-\infty,\infty] without changing the value of integral since for an arbitrary integrable function F(t,t^{\prime}), we have

 \displaystyle\begin{split}\displaystyle\int_{0}^{t^{+}}dt^{\prime}F(t,t^{% \prime})&\displaystyle\theta(t^{\prime})\theta(t-t^{\prime})\\ &\displaystyle=\int_{-\infty}^{+\infty}dt^{\prime}F(t,t^{\prime})\theta(t^{% \prime})\theta(t-t^{\prime}).\end{split} (113)

This extension of integration limits becomes handy when we write both \hat{\varphi}_{R}(x^{\prime},t^{\prime}) and G(x_{0},t|x^{\prime},t^{\prime}) in terms of their Fourier transforms in time. Focusing on the right boundary contribution at x^{\prime}=1^{+} we get

 \displaystyle\begin{split}&\displaystyle\int_{-\infty}^{+\infty}dt^{\prime}% \int_{-\infty}^{+\infty}\frac{d\omega_{1}}{2\pi}\int_{-\infty}^{+\infty}\frac{% d\omega_{2}}{2\pi}\left[\hat{\tilde{\varphi}}_{R}(x^{\prime},\omega_{1})% \partial_{x^{\prime}}\tilde{G}(x_{0},x^{\prime},\omega_{2})\right.\\ &\displaystyle\left.\left.-\tilde{G}(x_{0},x^{\prime},\omega_{2})\partial_{x^{% \prime}}\hat{\tilde{\varphi}}_{R}(x^{\prime},\omega_{1})\right]\right|_{x^{% \prime}=1^{+}}e^{-i\omega_{1}t^{\prime}}e^{-i\omega_{2}(t-t^{\prime})}.\end{split} (114)

Next, we write \hat{\tilde{\varphi}}_{R}(x^{\prime},\omega) as the sum of “incoming” and “outgoing” parts

 \displaystyle\hat{\tilde{\varphi}}_{R}(1^{+},\omega_{1})=\hat{\tilde{\varphi}}% _{R}^{inc}(1^{+},\omega_{1})+\hat{\tilde{\varphi}}_{R}^{out}(1^{+},\omega_{1}), (115)

defined as

 \displaystyle\partial_{x^{\prime}}\hat{\tilde{\varphi}}_{R}^{out}(x^{\prime}=1% ^{+},\omega_{1})=+i\omega_{1}\hat{\tilde{\varphi}}_{R}^{out}(x^{\prime}=1^{+},% \omega_{1}), (116a) \displaystyle\partial_{x^{\prime}}\hat{\tilde{\varphi}}_{R}^{inc}(x^{\prime}=1% ^{+},\omega_{1})=-i\omega_{1}\hat{\tilde{\varphi}}_{R}^{inc}(x^{\prime}=1^{+},% \omega_{1}). (116b)

On the other hand, since we are using a retarded GF with outgoing BC we have

 \displaystyle\partial_{x^{\prime}}\tilde{G}(x_{0},x^{\prime}=1^{+},\omega_{2})% =+i\omega_{2}\tilde{G}(x_{0},x^{\prime}=1^{+},\omega_{2}). (117)

By substituting Eqs. (116a, (116b) and (117) into Eq. (114), the integrand becomes

 \displaystyle\begin{split}&\displaystyle i(\omega_{1}+\omega_{2})\tilde{G}(x_{% 0},1^{+},\omega_{2})\hat{\tilde{\varphi}}_{R}^{inc}(1^{+},\omega_{1})\\ &\displaystyle+i(\omega_{2}-\omega_{1})\tilde{G}(x_{0},1^{+},\omega_{2})\hat{% \tilde{\varphi}}_{R}^{out}(1^{+},\omega_{1})\end{split} (118)

By taking the integral in t^{\prime} as \int_{-\infty}^{\infty}dt^{\prime}e^{i(\omega_{2}-\omega_{1})t^{\prime}}=2\pi% \delta(\omega_{1}-\omega_{2}), Eq. (114) can be simplified as

 \displaystyle\int_{-\infty}^{+\infty}\frac{d\omega}{2\pi}\left[2i\omega\tilde{% G}(x_{0},x^{\prime}=1^{+},\omega)\hat{\tilde{\varphi}}_{R}^{inc}(0^{-},\omega)% \right]e^{-i\omega t}, (119)

which indicates that only the incoming part of the field leads to a non-zero contribution to the field inside the resonator. A similiar expression holds for the left boundary with the difference that the incoming wave at the left waveguide is “right-going” in contrast to the right waveguide

 \displaystyle\int_{-\infty}^{+\infty}\frac{d\omega}{2\pi}\left[2i\omega\tilde{% G}(x_{0},x^{\prime}=0^{-},\omega)\hat{\tilde{\varphi}}_{L}^{inc}(0^{-},\omega)% \right]e^{-i\omega t}. (120)

The initial condition (IC) terms can be expressed in a compact form as

 \displaystyle\begin{split}\displaystyle\int_{x_{1}}^{x_{2}}dx^{\prime}\int_{-% \infty}^{\infty}\frac{d\omega}{2\pi}\left\{\chi(x^{\prime},x_{0})\tilde{G}(x_{% 0},x^{\prime},\omega)\right.\\ \displaystyle\left.\left[\hat{\dot{\varphi}}(x^{\prime},0)-i\omega\hat{\varphi% }(x^{\prime},0)\right]\right\}e^{-i\omega t}.\end{split} (121)

Gathering all the contributions, plugging it in the RHS of Eq. (78a) and defining a family of memory kernels

 \displaystyle\mathcal{K}_{n}(\tau)\equiv\gamma\chi_{s}\int_{-\infty}^{+\infty}% \frac{d\omega}{2\pi}\omega^{n}\tilde{G}(x_{0},x_{0},\omega)e^{-i\omega\tau}, (122a) and transfer functions \displaystyle\mathcal{D}_{R}(\omega)\equiv-2i\gamma\omega^{3}\tilde{G}(x_{0},1% ^{+},\omega), (122b) \displaystyle\mathcal{D}_{L}(\omega)\equiv-2i\gamma\omega^{3}\tilde{G}(x_{0},0% ^{-},\omega), (122c) \displaystyle\mathcal{I}(x^{\prime},\omega)\equiv\gamma\omega^{2}\chi(x^{% \prime},x_{0})\tilde{G}(x_{0},x^{\prime},\omega), (122d)

the effective dynamics of the transmon is found to be

 \displaystyle\begin{split}&\displaystyle\hat{\ddot{\varphi}}_{j}(t)+(1-\gamma)% \omega_{j}^{2}\sin{\left[\hat{\varphi}_{j}(t)\right]}=\\ \displaystyle+&\displaystyle\frac{d^{2}}{dt^{2}}\int_{0}^{t}dt^{\prime}% \mathcal{K}_{0}(t-t^{\prime})\omega_{j}^{2}\sin{\left[\hat{\varphi}_{j}(t^{% \prime})\right]}\\ \displaystyle+&\displaystyle\int_{-\infty}^{+\infty}\frac{d\omega}{2\pi}% \mathcal{D}_{R}(\omega)\hat{\tilde{\varphi}}_{R}^{inc}(1^{+},\omega)e^{-i% \omega t}\\ \displaystyle+&\displaystyle\int_{-\infty}^{+\infty}\frac{d\omega}{2\pi}% \mathcal{D}_{L}(\omega)\hat{\tilde{\varphi}}_{L}^{inc}(0^{-},\omega)e^{-i% \omega t}\\ \displaystyle+&\displaystyle\int_{0^{-}}^{1^{+}}dx^{\prime}\int_{-\infty}^{+% \infty}\frac{d\omega}{2\pi}\mathcal{I}(x^{\prime},\omega)\left[i\omega\hat{% \varphi}(x^{\prime},0)-\hat{\dot{\varphi}}(x^{\prime},0)\right]e^{-i\omega t}.% \end{split} (123)

This is Eq. (29) in Sec. III.

### B.5 Effective dynamics for spontaneous emission

Equation (123) is the most generic effective dynamics of a transmon coupled to an open multimode resonator. In this section, we find the effective dynamics for the problem of spontaneous emission where the system starts from the IC

 \displaystyle\hat{\rho}(0)=\hat{\rho}_{j}(0)\otimes\ket{0}_{ph}\bra{0}_{ph}. (124)

In the absence of external drive and due to the interaction with the leaky modes of the resonator, the system reaches its ground state \hat{\rho}_{g}\equiv\ket{0}_{j}\bra{0}_{j}\otimes\ket{0}_{ph}\bra{0}_{ph} in steady state.

Note that due the specific IC (124), there is no contribution from IC of the resonator in Eq. (123). To show this explicitly, recall that at t=0 the interaction has not turned on and we can represent \hat{\varphi}(x,0) and \hat{\dot{\varphi}}(x,0) in terms of a set of Hermitian modes of the resonator as Malekakhlagh and Türeci (2016)

 \displaystyle\hat{\varphi}(x,0)=\hat{\mathbf{1}}_{j}\otimes\sum\limits_{n}% \left(\frac{\hbar}{2\omega_{n}^{(H)}cL}\right)^{1/2}\left[\hat{a}_{n}(0)+\hat{% a}_{n}^{{\dagger}}(0)\right]\tilde{\Phi}_{n}^{(H)}(x), (125a) \displaystyle\hat{\dot{\varphi}}(x,0)=\hat{\mathbf{1}}_{j}\otimes\sum\limits_{% n}-i\left(\frac{\hbar\omega_{n}^{(H)}}{2cL}\right)^{1/2}\left[\hat{a}_{n}(0)-% \hat{a}_{n}^{{\dagger}}(0)\right]\tilde{\Phi}_{n}^{(H)}(x), (125b)

where we have used superscript notation (H) to distinguish Hermitian from non-Hermitian modes. By taking the partial trace over the photonic sector we find

 \displaystyle\begin{split}&\displaystyle\operatorname{Tr}_{ph}\left\{\hat{\rho% }_{ph}\left[\hat{a}_{n}(0)\pm\hat{a}_{n}^{{\dagger}}(0)\right]\right\}\\ &\displaystyle=\bra{0}_{ph}\left[\hat{a}_{n}(0)\pm\hat{a}_{n}^{{\dagger}}(0)% \right]\ket{0}_{ph}=0.\end{split} (126)

With no external drive, \hat{\tilde{\varphi}}_{R,L}^{inc} do not have a coherent part and their expectation value vanish due to the same reasoning as Eq. (126). Therefore, the effective dynamics for the spontaneous emission problem reduces to

 \displaystyle\begin{split}&\displaystyle\hat{\ddot{\phi}}_{j}(t)+(1-\gamma)% \omega_{j}^{2}\operatorname{Tr}_{ph}\left\{\hat{\rho}_{ph}(0)\sin{\left[\hat{% \varphi}_{j}(t)\right]}\right\}\\ &\displaystyle=\frac{d^{2}}{dt^{2}}\int_{0}^{t}dt^{\prime}\mathcal{K}_{0}(t-t^% {\prime})\omega_{j}^{2}\operatorname{Tr}_{ph}\left\{\hat{\rho}_{ph}(0)\sin{% \left[\hat{\varphi}_{j}(t^{\prime})\right]}\right\}.\end{split} (127)

Taking the second derivative of the RHS using Leibniz integral rule, and bringing the terms evaluated at the integral limits to the LHS gives

 \displaystyle\begin{split}&\displaystyle\hat{\ddot{\phi}}_{j}(t)-\omega_{j}^{2% }\mathcal{K}_{0}(0)\operatorname{Tr}_{ph}\left\{\hat{\rho}_{ph}(0)\cos{\left[% \hat{\varphi}_{j}(t)\right]}\hat{\dot{\varphi}}_{j}(t)\right\}\\ &\displaystyle+\omega_{j}^{2}\left[1-\gamma+i\mathcal{K}_{1}(0)\right]% \operatorname{Tr}_{ph}\left\{\hat{\rho}_{ph}(0)\sin{\left[\hat{\varphi}_{j}(t)% \right]}\right\}\\ &\displaystyle=-\int_{0}^{t}dt^{\prime}\mathcal{K}_{2}(t-t^{\prime})\omega_{j}% ^{2}\operatorname{Tr}_{ph}\left\{\hat{\rho}_{ph}(0)\sin{\left[\hat{\varphi}_{j% }(t^{\prime})\right]}\right\},\end{split} (128)

where we have used Eq. (122a) to rewrite time-derivatives of \mathcal{K}_{0}(\tau) in terms of \mathcal{K}_{n}(\tau).

### B.6 Spectral representation of \mathcal{K}_{0}, \mathcal{K}_{1} and \mathcal{K}_{2}

In this section, we express the contributions from the kernels \mathcal{K}_{0}(0), \mathcal{K}_{1}(0) and \mathcal{K}_{2}(\tau) appearing in Eq. (128) in terms of the spectral representation of the GF. For this purpose, we use the partial fraction expansion of the GF in agreement with Leung et al. (1994, 1997); Severini et al. (2004); Muljarov et al. (2010); Doost et al. (2013); Kristensen et al. (2015) in terms of its simple poles discussed in Sec. B.3 as

 \displaystyle\tilde{G}(x,x^{\prime},\omega)=\sum\limits_{n\in\mathbb{Z}}\frac{% 1}{2\omega}\frac{\tilde{\Phi}_{n}(x)\tilde{\Phi}_{n}(x^{\prime})}{\omega-% \omega_{n}}, (129)

where \tilde{\Phi}_{n}(x)\propto\tilde{\Phi}_{n}(x,\omega=\omega_{n}) is the quasi-bound eigenfunction.

Let us first calculate \mathcal{K}_{2}(\tau). By choosing an integration contour in the complex \omega-plane shown in Fig. (a)a and applying Cauchy’s residue theorem Mitrinovic and Keckic (1984); Hassani (2013) we find

 \displaystyle\begin{split}&\displaystyle\oint_{C}d\omega\omega^{2}\tilde{G}(x_% {0},x_{0},\omega)e^{-i\omega\tau}\\ &\displaystyle=\int_{I}d\omega\omega^{2}\tilde{G}(x_{0},x_{0},\omega)e^{-i% \omega\tau}+\int_{II}d\omega\omega^{2}\tilde{G}(x_{0},x_{0},\omega)e^{-i\omega% \tau}\\ &\displaystyle=-2\pi i\sum\limits_{n=0}^{\infty}\frac{1}{2}\left[\omega_{n}[% \tilde{\Phi}_{n}(x_{0})]^{2}e^{-i\omega_{n}\tau}-\omega_{n}^{*}[\tilde{\Phi}_{% n}^{*}(x_{0})]^{2}e^{+i\omega_{n}^{*}\tau}\right]\\ &\displaystyle=-2\pi\sum\limits_{n=0}^{\infty}|\omega_{n}||\tilde{\Phi}_{n}(x_% {0})|^{2}\sin{[\nu_{n}\tau+\theta_{n}-2\delta_{n}(x_{0})]}e^{-\kappa_{n}\tau},% \end{split} (130)

where due to nonzero opening of the resonator, both \omega_{n} and \tilde{\Phi}_{n}(x) are in general complex valued. Therefore, we have defined

 \displaystyle\theta_{n}\equiv\arctan{\left(\frac{\kappa_{n}}{\nu_{n}}\right)}, (131) \displaystyle\delta_{n}(x)\equiv\arctan{\left(\frac{\text{Im}[{\tilde{\Phi}_{n% }(x)}]}{\text{Re}[{\tilde{\Phi}_{n}(x)}]}\right)}. (132)

As the radius of the half-circle in Fig. (a)a is taken to infinity, \int_{II}d\omega\omega^{2}G(x_{0},x_{0},\omega) approaches zero. This can be checked by a change of variables

Substituting this into \int_{II} and taking the limit R_{II}\to\infty gives

 \displaystyle\begin{split}&\displaystyle\lim\limits_{R_{II}\to\infty}\int_{II}% d\omega\omega^{2}\tilde{G}(x_{0},x_{0},\omega)e^{-i\omega\tau}\\ &\displaystyle=\sum\limits_{n=0}^{\infty}\lim\limits_{R_{II}\to\infty}\int_{II% }d\omega\frac{\omega(\omega+i\kappa_{n})[\tilde{\Phi}_{n}(x_{0})]^{2}}{(\omega% -\omega_{n})(\omega+\omega_{n}^{*})}e^{-i\omega\tau}\\ &\displaystyle\propto\int_{0}^{\pi}d\psi\lim\limits_{R_{II}\to\infty}e^{-iR_{% II}\tau\cos{(\psi)}}R_{II}e^{-R_{II}\tau\sin{(\psi)}}=0,\ \tau>0.\end{split} (134)

On the other hand, \int_{I} in this limit reads

 \displaystyle\begin{split}&\displaystyle\lim\limits_{R_{II}\to\infty}\int_{I}d% \omega\omega^{2}\tilde{G}(x_{0},x_{0},\omega)e^{-i\omega\tau}\\ &\displaystyle=\int_{-\infty}^{\infty}d\omega\omega^{2}\tilde{G}(x_{0},x_{0},% \omega)e^{-i\omega\tau},\end{split} (135)

which is the quantity of interest. Therefore, we find

 \displaystyle\begin{split}&\displaystyle\int_{-\infty}^{\infty}d\omega\omega^{% 2}\tilde{G}(x_{0},x_{0},\omega)e^{-i\omega\tau}\\ &\displaystyle=-2\pi\sum\limits_{n=0}^{\infty}|\omega_{n}||\tilde{\Phi}_{n}(x_% {0})|^{2}\sin{\left[\nu_{n}\tau+\theta_{n}-2\delta_{n}(x_{0})\right]}e^{-% \kappa_{n}\tau}.\end{split} (136)

From this, we obtain the spectral representation of \mathcal{K}_{2}(\tau) as

 \displaystyle\mathcal{K}_{2}(\tau)=-\sum\limits_{n=0}^{\infty}A_{n}\sin{\left[% \nu_{n}\tau+\theta_{n}-2\delta_{n}(x_{0})\right]}e^{-\kappa_{n}\tau}, (137)

with A_{n}\equiv\gamma\chi_{s}\sqrt{\nu_{n}^{2}+\kappa_{n}^{2}}\left|\tilde{\Phi}_{% n}(x_{0})\right|^{2}.

\mathcal{K}_{1}(0) can be found through similar complex integration

 \displaystyle\begin{split}&\displaystyle\oint_{C}d\omega\omega\tilde{G}(x_{0},% x_{0},\omega)\\ &\displaystyle=\int_{I}d\omega\omega\tilde{G}(x_{0},x_{0},\omega)+\int_{II}d% \omega\omega\tilde{G}(x_{0},x_{0},\omega)\\ &\displaystyle=-2\pi i\sum\limits_{n=0}^{\infty}\left[\frac{[\tilde{\Phi}_{n}(% x_{0})]^{2}}{2}+\frac{[\tilde{\Phi}_{n}^{*}(x_{0})]^{2}}{2}\right]\\ &\displaystyle=-2\pi i\sum\limits_{n=0}^{\infty}|\tilde{\Phi}_{n}(x_{0})|^{2}% \cos{[2\delta_{n}(x_{0})]}\end{split} (138)

It can be shown again that \int_{II}\to 0 as R_{II}\to\infty from which we find that

 \displaystyle\begin{split}\displaystyle i\mathcal{K}_{1}(0)&\displaystyle=% \gamma\chi_{s}\sum\limits_{n=0}^{\infty}|\tilde{\Phi}_{n}(x_{0})|^{2}\cos{[2% \delta_{n}(x_{0})]}\\ &\displaystyle=\sum\limits_{n=0}^{\infty}\frac{A_{n}}{\sqrt{\nu_{n}^{2}+\kappa% _{n}^{2}}}\cos{[2\delta_{n}(x_{0})]}.\end{split} (139)

\mathcal{K}_{0}(0) has an extra pole at \omega=0, so the previous contour is not well defined. Therefore, we shift the integration contour as shown in Fig. 13. Then, we have

 \displaystyle\begin{split}&\displaystyle\oint_{C}d\omega\tilde{G}(x_{0},x_{0},% \omega)\\ &\displaystyle=\int_{I}d\omega\tilde{G}(x_{0},x_{0},\omega)+\int_{II}d\omega% \tilde{G}(x_{0},x_{0},\omega)\\ &\displaystyle=-2\pi i\sum\limits_{n=0}^{\infty}\frac{1}{2}\left[\frac{[\tilde% {\Phi}_{n}(x_{0})]^{2}}{\omega_{n}}-\frac{[\tilde{\Phi}_{n}^{*}(x_{0})]^{2}}{% \omega_{n}^{*}}\right]\\ &\displaystyle-2\pi i\sum\limits_{n=0}^{\infty}\frac{1}{2}\left[\frac{[\tilde{% \Phi}_{n}(x_{0})]^{2}}{-\omega_{n}}+\frac{[\tilde{\Phi}_{n}^{*}(x_{0})]^{2}}{% \omega_{n}^{*}}\right]=0,\\ \end{split} (140)

where the first sum comes from the residues at \omega=\omega_{n} and \omega=-\omega_{n}^{*}, while the last sum is the residue at \omega=0 and they completely cancel each other and we get

 \displaystyle\mathcal{K}_{0}(0)=0. (141)

From Eq. (141) we find that the effective dynamics for the spontaneous emission problem simplifies to

 \displaystyle\begin{split}&\displaystyle\hat{\ddot{\phi}}_{j}(t)+\omega_{j}^{2% }\left[1-\gamma+i\mathcal{K}_{1}(0)\right]\operatorname{Tr}_{ph}\left\{\hat{% \rho}_{ph}(0)\sin{\left[\hat{\varphi}_{j}(t)\right]}\right\}\\ &\displaystyle=-\int_{0}^{t}dt^{\prime}\mathcal{K}_{2}(t-t^{\prime})\omega_{j}% ^{2}\operatorname{Tr}_{ph}\left\{\hat{\rho}_{ph}(0)\sin{\left[\hat{\varphi}_{j% }(t^{\prime})\right]}\right\}.\end{split} (142)

## Appendix C Characteristic function D_{j}(s) for the linear equations of motion

Up to linear order, transmon acts as a simple harmonic oscillator and we find we find

 \displaystyle\begin{split}\displaystyle\hat{\ddot{X}}_{j}(t)&\displaystyle+% \omega_{j}^{2}\left[1-\gamma+i\mathcal{K}_{1}(0)\right]\hat{X}_{j}(t)\\ &\displaystyle=-\int_{0}^{t}dt^{\prime}\mathcal{K}_{2}(t-t^{\prime})\omega_{j}% ^{2}\hat{X}_{j}(t^{\prime}).\end{split} (143)

Equation (143) is a linear integro-differential equation with a memory integral on the RHS, appearing as the convolution of the memory kernel \mathcal{K}_{2} with earlier values of \hat{X}_{j}. It can be solved by means of unilateral Laplace transform Abramowitz and Stegun (1964); Korn and Korn (1968); Hassani (2013) defined as

 \displaystyle\tilde{f}(s)\equiv\int_{0}^{\infty}dte^{-st}f(t). (144)

Employing the following properties of Laplace transform:

• Convolution

 \displaystyle\begin{split}&\displaystyle\mathfrak{L}\left\{\int_{0}^{t}dt^{% \prime}f(t^{\prime})g(t-t^{\prime})\right\}=\mathfrak{L}\left\{\int_{0}^{t}dt^% {\prime}f(t-t^{\prime})g(t^{\prime})\right\}\\ &\displaystyle=\mathfrak{L}\left\{f(t)\right\}\cdot\mathfrak{L}\left\{g(t)% \right\}=\tilde{f}(s)\tilde{g}(s),\end{split} (145)
• General derivative

 \displaystyle\mathfrak{L}\left\{\frac{d^{N}}{dt^{N}}f(t)\right\}=s^{N}\tilde{f% }(s)-\sum\limits_{n=1}^{N}s^{N-n}\left.\frac{d^{n-1}}{dt^{n-1}}f(t)\right|_{t=% 0}, (146)

we can transform the integro-differential Eq. (143) into a closed algebraic form in terms of \hat{\tilde{X}}_{j}(s) as

 \displaystyle\hat{\tilde{X}}_{j}(s)=\frac{s\hat{X}_{j}(0)+\hat{\dot{X}}_{j}(0)% }{D_{j}(s)}=\frac{s\hat{X}_{j}(0)+\omega_{j}\hat{Y}_{j}(0)}{D_{j}(s)}, (147)

where we have defined

 \displaystyle D_{j}(s)\equiv s^{2}+\Omega^{2}(s), (148a) \displaystyle\Omega^{2}(s)\equiv\omega_{j}^{2}\left[1-\gamma+i\mathcal{K}_{1}(% 0)+\tilde{\mathcal{K}}_{2}(s)\right]. (148b)

and \hat{Y}_{j} is the normalized charge variable and is canonically conjugate to \hat{X}_{j} such that [\hat{X}_{j}(0),\hat{Y}_{j}(0)]=2i.

Note that in order to solve for \hat{X}_{j}(t) from Eq. (147), one has to take the inverse Laplace transform of the resulting algebraic form in s. This requires studying the denominator first which determines the poles of the entire system up to linear order. Using the expressions for \mathcal{K}_{2}(\tau_{1}) and i\mathcal{K}_{1}(0) given in Eqs. (137) and (139) we find

 \displaystyle\begin{split}&\displaystyle i\mathcal{K}_{1}(0)+\tilde{\mathcal{K% }}_{2}(s)=\sum\limits_{n\in\mathbb{N}}\frac{A_{n}}{\sqrt{\nu_{n}^{2}+\kappa_{n% }^{2}}}\cos{[2\delta_{n}(x_{0})]}\\ &\displaystyle-\sum\limits_{n\in\mathbb{N}}A_{n}\frac{\cos{[\theta_{n}-2\delta% _{n}(x_{0})]}\nu_{n}+\sin{[\theta_{n}-2\delta_{n}(x_{0})]}(s+\kappa_{n})}{(s+% \kappa_{n})^{2}+\nu_{n}^{2}}.\end{split} (149)

Expanding the sine and cosine in the numerator of the second term in Eq. (149) as

 \displaystyle\begin{split}&\displaystyle\cos{[\theta_{n}-2\delta_{n}(x_{0})]}% \nu_{n}+\sin{[\theta_{n}-2\delta_{n}(x_{0})]}(s+\kappa_{n})\\ &\displaystyle=\left\{\cos{(\theta_{n})}\cos{[2\delta_{n}(x_{0})]}+\sin{(% \theta_{n})}\sin{[2\delta_{n}(x_{0})]}\right\}\nu_{n}\\ &\displaystyle+\left\{\sin{(\theta_{n})}\cos{[2\delta_{n}(x_{0})]}-\cos{(% \theta_{n})}\sin{[2\delta_{n}(x_{0})]}\right\}(s+\kappa_{n})\\ &\displaystyle=\frac{\left\{\kappa_{n}\cos{[2\delta_{n}(x_{0})]}-\nu_{n}\sin{[% 2\delta_{n}(x_{0})]}\right\}s}{\sqrt{\nu_{n}^{2}+\kappa_{n}^{2}}}\\ &\displaystyle+\frac{(\nu_{n}^{2}+\kappa_{n}^{2})\cos{[2\delta_{n}(x_{0})]}}{% \sqrt{\nu_{n}^{2}+\kappa_{n}^{2}}},\end{split} (150)

Eq. (149) simplifies to

 \displaystyle\begin{split}&\displaystyle\sum\limits_{n=0}^{\infty}\frac{A_{n}}% {\sqrt{\nu_{n}^{2}+\kappa_{n}^{2}}}\left\{\cos{[2\delta_{n}(x_{0})]}-\frac{(% \nu_{n}^{2}+\kappa_{n}^{2})\cos{[2\delta_{n}(x_{0})]}}{(s+\kappa_{n})^{2}+\nu_% {n}^{2}}\right.\\ &\displaystyle-\left.\frac{\left\{\kappa_{n}\cos{[2\delta_{n}(x_{0})]}-\nu_{n}% \sin{[2\delta_{n}(x_{0})]}\right\}s}{(s+\kappa_{n})^{2}+\nu_{n}^{2}}\right\}\\ &\displaystyle=\sum\limits_{n=0}^{\infty}M_{n}\frac{s\{\cos{[2\delta_{n}(x_{0}% )]}s+\sin{[2\delta_{n}(x_{0})]}\nu_{n}\}}{(s+\kappa_{n})^{2}+\nu_{n}^{2}},\end% {split} (151)

where we have defined

 \displaystyle M_{n}\equiv\frac{A_{n}}{\sqrt{\nu_{n}^{2}+\kappa_{n}^{2}}}=% \gamma\chi_{s}|\tilde{\Phi}_{n}(x_{0})|^{2}. (152)

Therefore, D_{j}(s) simplifies to

 \displaystyle\begin{split}&\displaystyle D_{j}(s)=s^{2}+\omega_{j}^{2}+\\ &\displaystyle\underbrace{\omega_{j}^{2}\left\{-\gamma+\sum\limits_{n=0}^{% \infty}M_{n}\frac{s\{\cos{[2\delta_{n}(x_{0})]}s+\sin{[2\delta_{n}(x_{0})]}\nu% _{n}\}}{(s+\kappa_{n})^{2}+\nu_{n}^{2}}\right\}}_{Modification\ due\ to\ % memory}.\end{split} (153)

## Appendix D Multi-Scale Analysis

In order to understand the application of MSPT on the problem of spontaneous emission, we have broken down its complexity into simpler toy problems, discussing each in a separate subsection. In Sec. D.1, we revisit the classical Duffing oscillator problem Bender and Orszag (1999) in the presence of dissipation, to study the interplay of nonlinearity and dissipation. In Sec. D.2, we discuss the free quantum Duffing oscillator to show how the non-commuting algebra of quantum mechanics alters the classical solution. Finally, in Sec. D.3, we study the full problem and provide the derivation for the MSPT solution (60).

### D.1 Classical Duffing oscillator with dissipation

Consider a classical Duffing oscillator

 \displaystyle\ddot{X}(t)+\delta\,\omega\dot{X}(t)+\omega^{2}\left[X(t)-% \varepsilon X^{3}(t)\right]=0, (154)

with initial condition X(0)=X_{0}, \dot{X}(0)=\omega Y_{0}. In order to have a bound solution, it is sufficient that the initial energy of the system be less than the potential energy evaluated at its local maxima, X_{max}\equiv\pm\sqrt{1/3\varepsilon} , i.e. E_{0}<U(X_{max}) which in terms of the initial conditions X_{0} and Y_{0} reads

 \displaystyle\frac{1}{2}Y_{0}^{2}+\frac{1}{2}\left(X_{0}^{2}-\varepsilon X_{0}% ^{4}\right)<\frac{5}{36\varepsilon}. (155)

Note that a naive use of conventional perturbation theory decomposes the solution into a series X(t)=X^{(0)}(t)+\varepsilon X^{(1)}(t)+\ldots, which leads to unbounded (secular) solutions in time. In order to illustrate this, consider the simple case where \delta=0, X_{0}=1 and Y_{0}=0. Then, we find

 \displaystyle\mathcal{O}(1):\ddot{X}^{(0)}(t)+\omega^{2}X^{(0)}(t)=0, (156a) \displaystyle\mathcal{O}(\varepsilon):\ddot{X}^{(1)}(t)+\omega^{2}X^{(1)}(t)=% \omega^{2}[X^{(0)}(t)]^{3}, (156b)

which leads to X^{(0)}(t)=\cos(\omega t) and X^{(1)}(t)=\frac{1}{32}\cos(\omega t)-\frac{1}{32}\cos(3\omega t)+\frac{3}{8}% \omega t\sin(\omega t). The latter has a secular contribution that grows unbounded in time.

The secular terms can be canceled order by order by introducing multiple time scales, which amounts to a resummation of the conventional perturbation series Bender and Orszag (1999). We assume small dissipation and nonlinearity, i.e. \delta,\varepsilon\ll 1. This allows us to define additional slow time scales \tau\equiv\varepsilon t and \eta\equiv\delta t in terms of which we can perform a multi-scale expansion for X(t) as

 \displaystyle\begin{split}\displaystyle X(t)&\displaystyle=x^{(0)}(t,\tau,\eta% )+\varepsilon x^{(1)}(t,\tau,\eta)\\ &\displaystyle+\delta y^{(1)}(t,\tau,\eta)+\mathcal{O}(\varepsilon^{2},\delta^% {2},\varepsilon\delta).\end{split} (157a) Using the chain rule, the total derivative d/dt is also expanded as \displaystyle d_{t}=\partial_{t}+\varepsilon\partial_{\tau}+\delta\partial_{% \eta}+\mathcal{O}(\varepsilon^{2},\delta^{2},\varepsilon\delta). (157b)

Plugging Eqs. (157a-157b) into Eq. (154) and collecting equal powers of \delta and \epsilon we find

 \displaystyle\mathcal{O}(1):\partial_{t}^{2}x^{(0)}+\omega^{2}x^{(0)}=0, (158a) \displaystyle\mathcal{O}(\delta):\partial_{t}^{2}y^{(1)}+\omega^{2}y^{(1)}=-% \omega\partial_{t}x^{(0)}-2\partial_{t}\partial_{\eta}x^{(0)}, (158b) \displaystyle\mathcal{O}(\varepsilon):\partial_{t}^{2}x^{(1)}+\omega^{2}x^{(1)% }=\omega^{2}\left[x^{(0)}\right]^{3}-2\partial_{t}\partial_{\tau}x^{(0)}. (158c)

The general solution to O(1) Eq. (158a) reads

 \displaystyle x^{(0)}(t,\tau,\eta)=a(\tau,\eta)e^{-i\omega t}+a^{*}(\tau,\eta)% e^{+i\omega t}. (159)

Plugging Eq. (159) into Eq. (158b) we find that in order to remove secular terms a(\tau,\eta) satisfies

 \displaystyle(2\partial_{\eta}+\omega)a(\tau,\eta)=0, (160)

which gives the \eta-dependence of a(\tau,\eta) as

 \displaystyle a(\tau,\eta)=\alpha(\tau)e^{-\frac{\omega}{2}\eta}. (161)

The condition that removes the secular term on the RHS of \mathcal{O}(\varepsilon) Eq. (158c) reads

 \displaystyle 2i\omega\partial_{\tau}a(\tau,\eta)+3\omega^{2}|a(\tau,\eta)|^{2% }a(\tau,\eta)=0. (162)

Multipliying Eq. (162) by a^{*}(\tau,\eta) and its complex conjugate by a(\tau,\eta) and taking the difference gives

 \displaystyle\partial_{\tau}|a(\tau,\eta)|^{2}=0, (163)

which together with Eq. (161) implies that

 \displaystyle|a(\tau,\eta)|^{2}=|\alpha(0)|^{2}e^{-\omega\eta}. (164)

Then, a(\tau,\eta) is found as

 \displaystyle a(\tau,\eta)=\alpha(0)e^{-\frac{\omega}{2}\eta}e^{i\frac{3}{2}% \omega|\alpha(0)|^{2}e^{-\omega\eta}\tau}. (165)

Replacing \tau=\varepsilon t and \eta=\delta t, and, the general solution up to \mathcal{O}(\varepsilon^{2},\delta^{2},\varepsilon\delta) reads

 \displaystyle\begin{split}\displaystyle X^{(0)}(t)=x^{(0)}(t,\varepsilon t,% \delta t)=e^{-\frac{\kappa}{2}t}\left[\alpha(0)e^{-i\bar{\omega}(t)t}+c.c.% \right],\end{split} (166)

where we have defined the decay rate \kappa\equiv\delta.\omega and a normalized frequency \bar{\omega}(t) as

 \displaystyle\bar{\omega}(t)\equiv\left[1-\frac{3\varepsilon}{2}|\alpha(0)|^{2% }e^{-\kappa t}\right]\omega. (167)

Furthermore, \alpha(0) is determined based on initial conditions as \alpha(0)=(X_{0}+iY_{0})/2.

A comparison between the numerical solution (blue), \mathcal{O}(1) MSPT solution (166) (red) and linear solution (black) is made in Fig. 14 for the first ten oscillation periods. The MSPT solution captures the true oscillation frequency better than the linear solution. However, it is only valid for \omega t\ll\varepsilon^{-2},\delta^{-2},\varepsilon^{-1}\delta^{-1} up to this order in perturbation theory.

### D.2 A free quantum Duffing oscillator

Consider a free quantum Duffing oscillator that obeys

 \displaystyle\hat{\ddot{X}}(t)+\omega^{2}\left[\hat{X}(t)-\varepsilon\hat{X}^{% 3}(t)\right]=0, (168)

with operator initial conditions

such that \hat{X}(0) and \hat{Y}(0) are canonically conjugate variables and obey [\hat{X}(0),\hat{Y}(0)]=2i\hat{\mathbf{1}}.

Next, we expand \hat{X}(t) and d/dt up to \mathcal{O}(\varepsilon^{2}) as

 \displaystyle\hat{X}(t)=\hat{x}^{(0)}(t,\tau)+\varepsilon\hat{x}^{(1)}(t,\tau)% +\mathcal{O}(\varepsilon^{2}), (170a) \displaystyle d_{t}=\partial_{t}+\varepsilon\partial_{\tau}+\mathcal{O}(% \varepsilon^{2}). (170b)

Plugging this into Eq. (168) and collecting equal powers of \varepsilon gives

 \displaystyle\mathcal{O}(1):\partial_{t}^{2}\hat{x}^{(0)}+\omega^{2}\hat{x}^{(% 0)}=0, (171a) \displaystyle\mathcal{O}(\varepsilon):\partial_{t}^{2}\hat{x}^{(1)}+\omega^{2}% \hat{x}^{(1)}=\omega^{2}\left[\hat{x}^{(0)}\right]^{3}-2\partial_{t}\partial_{% \tau}\hat{x}^{(0)}. (171b)

Up to \mathcal{O}(1), the general solution reads

 \displaystyle\hat{x}^{(0)}(t,\tau)=\hat{a}(\tau)e^{-i\omega t}+\hat{a}^{{% \dagger}}(\tau)e^{+i\omega t} (172)

Furthermore, from the commutation relation [\hat{x}(t,\tau),\hat{y}(t,\tau)]=2i\hat{\mathbf{1}} we find that [\hat{a}(\tau),\hat{a}^{{\dagger}}(\tau)]=\hat{\mathbf{1}}. Substituting Eq. (172) into the RHS of Eq. (171b) and setting the secular term oscillating at \omega to zero we obtain

 \displaystyle\begin{split}&\displaystyle 2i\omega\frac{d\hat{a}(\tau)}{d\tau}+% \omega^{2}\left[\hat{a}(\tau)\hat{a}(\tau)\hat{a}^{{\dagger}}(\tau)\right.\\ &\displaystyle\left.+\hat{a}(\tau)\hat{a}^{{\dagger}}(\tau)\hat{a}(\tau)+\hat{% a}^{{\dagger}}(\tau)\hat{a}(\tau)\hat{a}(\tau)\right]=0,\end{split} (173)

The condition that removes secular term at -\omega, appears as Hermitian conjugate of Eq. (173).

Using [\hat{a}(\tau),\hat{a}^{{\dagger}}(\tau)]=1, Eq. (173) can be rewritten in a compact form

 \displaystyle\frac{d\hat{a}(\tau)}{d\tau}-i\frac{3\omega}{4}\left[\hat{% \mathcal{H}}(\tau)\hat{a}(\tau)+\hat{a}(\tau)\hat{\mathcal{H}}(\tau)\right]=0, (174)

where

 \displaystyle\hat{\mathcal{H}}(\tau)\equiv\frac{1}{2}\left[\hat{a}^{{\dagger}}% (\tau)\hat{a}(\tau)+\hat{a}(\tau)\hat{a}^{{\dagger}}(\tau)\right]. (175)

Next, we show that \hat{\mathcal{H}}(\tau) is a conserved quantity. Pre- and post-multiplying Eq. (173) by \hat{a}^{{\dagger}}(\tau), pre- and post-multiplying Hermitian conjugate of Eq. (173) by \hat{a}(\tau) and adding all the terms gives

 \displaystyle\frac{d\hat{\mathcal{H}}(\tau)}{d\tau}=0, (176)

which implies that \hat{\mathcal{H}}(\tau)=\hat{\mathcal{H}}(0). Therefore, we find the solution for \hat{a}(\tau) as

 \displaystyle\hat{a}(\tau)=\mathcal{W}\left\{\hat{a}(0)\exp\left[+i\frac{3% \omega}{2}\hat{\mathcal{H}}(0)\tau\right]\right\}, (177)

where \mathcal{W}\{\bullet\} represents Weyl-ordering of operators Schleich (2011). The operator ordering \mathcal{W}\left\{\hat{a}(0)f\left(\hat{\mathcal{H}}(0)\tau\right)\right\} is defined as follows:

1. Expand f\left(\hat{\mathcal{H}}(0)\tau\right) as a Taylor series in powers of operator \hat{\mathcal{H}}(0)\tau,

2. Weyl-order the series term-by-term as \mathcal{W}\left\{\hat{a}(0)\left[\hat{\mathcal{H}}(0)\right]^{n}\right\}% \equiv\frac{1}{2^{n}}\sum\limits_{m=0}^{n}{{n}\choose{m}}\left[\hat{\mathcal{H% }}(0)\right]^{m}\hat{a}(0)\left[\hat{\mathcal{H}}(0)\right]^{n-m}.

The formal solution (177) can be re-expressed in a closed form Bender et al. (1986, 1987); Bender and Dunne (1988); Bender and Bettencourt (1996) using the properties of Euler polynomials Abramowitz and Stegun (1964) as

 \displaystyle\hat{a}(\tau)=\frac{\hat{a}(0)e^{i\frac{3\omega}{2}\hat{\mathcal{% H}}(0)\tau}+e^{i\frac{3\omega}{2}\hat{\mathcal{H}}(0)\tau}\hat{a}(0)}{2\cos% \left(\frac{3\omega\tau}{4}\right)}. (178)

Plugging Eq. (178) into Eq. (172) and substituting \tau=\varepsilon t, we find the solution for \hat{X}(t) up to \mathcal{O}(\varepsilon) as

 \displaystyle\begin{split}\displaystyle\hat{X}^{(0)}(t)=\hat{x}^{(0)}(t,% \varepsilon t)&\displaystyle=\frac{\hat{a}(0)e^{-i\hat{\bar{\omega}}t}+e^{-i% \hat{\bar{\omega}}t}\hat{a}(0)}{2\cos\left(\frac{3\omega}{4}\varepsilon t% \right)}\\ &\displaystyle+\frac{\hat{a}^{{\dagger}}(0)e^{+i\hat{\bar{\omega}}t}+e^{+i\hat% {\bar{\omega}}t}\hat{a}^{{\dagger}}(0)}{2\cos\left(\frac{3\omega}{4}% \varepsilon t\right)},\end{split} (179)

where \hat{\bar{\omega}}\equiv\omega[1-\frac{3\varepsilon}{2}\hat{\mathcal{H}}(0)] appears as a renormalized frequency operator.

The physical quantity of interest is the expectation value of \hat{X}^{(0)}(t) with respect to the initial density matrix \hat{\rho}(0). The number basis of the simple harmonic oscillator is a complete basis for the Hilbert space of the Duffing oscillator such that

 \displaystyle\hat{\rho}(0)=\sum\limits_{mn}c_{mn}\ket{m}\bra{n}. (180)

Therefore, calculation of \braket{\hat{X}^{(0)}(t)} reduces to calculating the matrix element \bra{m}\hat{a}(\varepsilon t)\ket{n}. From Eq. (178) we find that the only nonzero matrix element read

 \displaystyle\begin{split}\displaystyle\bra{n-1}\hat{a}(\varepsilon t)\ket{n}&% \displaystyle=\frac{\bra{n-1}\hat{a}(0)\ket{n}e^{i\frac{3\varepsilon\omega}{2}% \bra{n}\hat{\mathcal{H}}(0)\ket{n}}}{2\cos\left(\frac{3\varepsilon\omega}{4}t% \right)}\\ &\displaystyle+\frac{e^{i\frac{3\varepsilon\omega}{2}\bra{n-1}\hat{\mathcal{H}% }(0)\ket{n-1}}\bra{n-1}\hat{a}(0)\ket{n}}{2\cos\left(\frac{3\varepsilon\omega}% {4}t\right)}\\ &\displaystyle=\bra{n-1}\hat{a}(0)\ket{n}e^{i\frac{3n\varepsilon\omega}{2}t},% \end{split} (181)

where we used that \bra{n}\hat{\mathcal{H}}(0)\ket{n}=n+1/2 is diagonal in the number basis.

### D.3 Quantum Duffing oscillator coupled to a set of quantum harmonic oscillators

Quantum MSPT can also be applied to the problem of a quantum Duffing oscillator coupled to multiple harmonic oscillators. For simplicity, consider the toy Hamiltonian

 \displaystyle\begin{split}\displaystyle\hat{\mathcal{H}}&\displaystyle\equiv% \frac{\omega_{j}}{4}\left(\hat{\mathcal{X}}_{j}^{2}+\hat{\mathcal{Y}}_{j}^{2}-% \frac{\varepsilon}{2}\hat{\mathcal{X}}_{j}^{4}\right)\\ &\displaystyle+\frac{\omega_{c}}{4}\left(\hat{\mathcal{X}}_{c}^{2}+\hat{% \mathcal{Y}}_{c}^{2}\right)+g\hat{\mathcal{Y}}_{j}\hat{\mathcal{Y}}_{c},\end{split} (182)

where the nonlinearity only exists in the Duffing sector of the Hilbert space labeled as j. Due to linear coupling there will be a hybridization of modes up to linear order. Therefore, Hamiltonian (182) can always be rewritten in terms of the normal modes of its quadratic part as

 \displaystyle\begin{split}\displaystyle\hat{\mathcal{H}}&\displaystyle\equiv% \frac{\beta_{j}}{4}\left(\hat{\bar{\mathcal{X}}}_{j}^{2}+\hat{\bar{\mathcal{Y}% }}_{j}^{2}\right)+\frac{\beta_{c}}{4}\left(\hat{\bar{\mathcal{X}}}_{c}^{2}+% \hat{\bar{\mathcal{Y}}}_{c}^{2}\right)\\ &\displaystyle-\frac{\varepsilon\omega_{j}}{8}\left(u_{j}\hat{\bar{\mathcal{X}% }}_{j}+u_{c}\hat{\bar{\mathcal{X}}}_{c}\right)^{4},\end{split} (183)

where u_{j,c} are real hybridization coefficients and \hat{\bar{\mathcal{X}}}_{j,c} and \hat{\bar{\mathcal{Y}}}_{j,c} represent j-like and c-like canonical operators. For g=0, u_{j}\rightarrow 1, u_{c}\rightarrow 0, \hat{\bar{\mathcal{X}}}_{j,c}\rightarrow\hat{\mathcal{X}}_{j,c} and \hat{\bar{\mathcal{Y}}}_{j,c}\rightarrow\hat{\mathcal{Y}}_{j,c}. To find u_{j,c} consider the Heisenberg equations of motion

 \displaystyle\hat{\ddot{\mathcal{X}}}_{j}(t)+\omega_{j}^{2}\hat{\mathcal{X}}_{% j}(t)=-2g\omega_{c}\hat{\mathcal{X}}_{c}(t), (184a) \displaystyle\hat{\ddot{\mathcal{X}}}_{c}(t)+\omega_{c}^{2}\hat{\mathcal{X}}_{% c}(t)=-2g\omega_{j}\hat{\mathcal{X}}_{j}(t). (184b)

Expressing \vec{\mathcal{X}}\equiv(\hat{\mathcal{X}}_{j}\ \hat{\mathcal{X}}_{c})^{T}, the system above can be written as \vec{\ddot{\mathcal{X}}}+V\vec{\mathcal{X}}=0, where V is a 2\times 2 matrix. Plugging an Ansatz \vec{\mathcal{X}}=\vec{\mathcal{X}}_{0}e^{i\lambda t} leads to an eigensystem whose eigenvalues are \beta_{j,c} and whose eigenvectors give the hybridization coefficients u_{j,c}.

The Heisenberg equations of motion for the hybdridized modes \hat{\bar{\mathcal{X}}}_{l}(t), l\equiv j,c, reads

 \displaystyle\begin{split}\displaystyle\hat{\ddot{\bar{\mathcal{X}}}}_{l}(t)+% \beta_{l}^{2}\left\{\hat{\bar{\mathcal{X}}}_{l}(t)-\varepsilon_{l}\left[u_{j}% \hat{\bar{\mathcal{X}}}_{j}(t)+u_{c}\hat{\bar{\mathcal{X}}}_{c}(t)\right]^{3}% \right\}=0,\end{split} (185)

where due to hybridization, each oscillator experiences a distinct effective nonlinearity as \varepsilon_{l}\equiv\frac{\omega_{j}}{\beta_{l}}u_{l}\varepsilon. Therefore, we define two new time scales \tau_{l}\equiv\varepsilon_{l}t in terms of which we can expand

 \displaystyle\begin{split}\displaystyle\hat{\bar{\mathcal{X}}}_{l}(t)&% \displaystyle=\hat{\bar{x}}_{l}^{(0)}(t,\tau_{j},\tau_{c})+\varepsilon_{l}\hat% {\bar{x}}_{l}^{(1)}(t,\tau_{j},\tau_{c})\\ &\displaystyle+\varepsilon_{l^{\prime}}\hat{\bar{y}}_{l}^{(1)}(t,\tau_{j},\tau% _{c})+\mathcal{O}(\varepsilon_{j}^{2},\varepsilon_{c}^{2},\varepsilon_{j}% \varepsilon_{c}),\end{split} (186a) \displaystyle d_{t}=\partial_{t}+\varepsilon_{j}\partial_{\tau_{j}}+% \varepsilon_{c}\partial_{\tau_{c}}+\mathcal{O}(\varepsilon_{j}^{2},\varepsilon% _{c}^{2},\varepsilon_{j}\varepsilon_{c}). (186b)

where we have used the notation that if l=j, l^{\prime}=c and vice versa. Up to \mathcal{O}(1) we find

 \displaystyle\mathcal{O}(1):\partial_{t}^{2}\hat{\bar{x}}_{l}^{(0)}+\beta_{l}^% {2}\hat{\bar{x}}_{l}^{(0)}=0, (187)

 \displaystyle\begin{split}\displaystyle\hat{\bar{x}}_{l}^{(0)}(t,\tau_{j},\tau% _{c})&\displaystyle=\hat{\bar{a}}_{l}(\tau_{j},\tau_{c})e^{-i\beta_{l}t}\\ &\displaystyle+\hat{\bar{a}}_{l}^{{\dagger}}(\tau_{j},\tau_{c})e^{+i\beta_{l}t% }.\end{split} (188)

where

 \displaystyle[\hat{\bar{a}}_{l_{1}},\hat{\bar{a}}_{l_{2}}^{{\dagger}}]=\delta_% {l_{1}l_{2}}\hat{\mathbf{1}},\ [\hat{\bar{a}}_{l_{1}},\hat{\bar{a}}_{l_{2}}]=[% \hat{\bar{a}}_{l_{1}}^{{\dagger}},\hat{\bar{a}}_{l_{2}}^{{\dagger}}]=0. (189)

There are \mathcal{O}(\varepsilon_{l}) and \mathcal{O}(\varepsilon_{l^{\prime}}) equations of for each normal mode as

 \displaystyle\begin{split}\displaystyle\mathcal{O}(\varepsilon_{l})\ of\ l:&% \displaystyle\partial_{t}^{2}\hat{\bar{x}}_{l}^{(1)}+\beta_{l}^{2}\hat{\bar{x}% }_{l}^{(1)}=-2\partial_{t}\partial_{\tau_{l}}\hat{\bar{x}}_{l}^{(0)}\\ &\displaystyle-\beta_{l}^{2}\left[u_{j}\hat{\bar{x}}_{j}^{(0)}+u_{c}\hat{\bar{% x}}_{c}^{(0)}\right]^{3}=0,\end{split} (190a) \displaystyle\mathcal{O}(\varepsilon_{l^{\prime}})\ of\ l:\partial_{t}^{2}\hat% {\bar{y}}_{l}^{(1)}+\beta_{l}^{2}\hat{\bar{y}}_{l}^{(1)}=-2\partial_{t}% \partial_{\tau_{l^{\prime}}}\hat{\bar{x}}_{l}^{(0)}. (190b)

By setting the secular terms on the RHS of Eq. (190b) we find that \partial_{\tau_{l^{\prime}}}\hat{b}_{l}=0 which means that q and c sectors are only modified with their own time scale, i.e. \hat{\bar{a}}_{l}=\hat{\bar{a}}_{l}(\tau_{l}). Applying the same procedure on Eq. (190a) and using commutation relations (189) we find

 \displaystyle\begin{split}\displaystyle\frac{d\hat{\bar{a}}_{l}}{d\tau_{l}}-i% \frac{3\beta_{l}}{4}&\displaystyle\left\{u_{l}^{3}\left[\hat{\bar{\mathcal{H}}% }_{l}\hat{\bar{a}}_{l}+\hat{\bar{a}}_{l}\hat{\bar{\mathcal{H}}}_{l}\right]% \right.\\ &\displaystyle+\left.2u_{l}u_{l^{\prime}}^{2}\left[\hat{\bar{\mathcal{H}}}_{l^% {\prime}}\hat{\bar{a}}_{l}+\hat{\bar{a}}_{l}\hat{\bar{\mathcal{H}}}_{l^{\prime% }}\right]\right\}=0,\end{split} (191)

where

 \displaystyle\hat{\bar{\mathcal{H}}}_{l}(\tau_{l})\equiv\frac{1}{2}\left[\hat{% \bar{a}}_{l}^{{\dagger}}(\tau_{l})\hat{\bar{a}}_{l}(\tau_{l})+\hat{\bar{a}}_{l% }(\tau_{l})\hat{\bar{a}}_{l}^{{\dagger}}(\tau_{l})\right]. (192)

By pre- and post-multiplying Eq. (192) by \hat{\bar{a}}_{l}^{{\dagger}}(\tau_{l}) and its Hermitian conjugate by \hat{\bar{a}}_{l}(\tau_{l}) and adding them we find that

 \displaystyle\frac{d\hat{\bar{\mathcal{H}}}_{l}(\tau_{l})}{d\tau_{l}}=0, (193)

which means that the sub-Hamiltonians of each normal mode remain a constant of motion up to this order in perturbation. Therefore, in terms of effective Hamiltonians

 \displaystyle\hat{\bar{h}}_{l}(0)\equiv u_{l}^{3}\hat{\bar{\mathcal{H}}}_{l}(0% )+2u_{l}u_{l^{\prime}}^{2}\hat{\bar{\mathcal{H}}}_{l^{\prime}}(0), (194)

Eq. (195) simplifies to

 \displaystyle\frac{d\hat{\bar{a}}_{l}}{d\tau_{l}}-i\frac{3\beta_{l}}{4}\left[% \hat{\bar{h}}_{l}(0)\hat{\bar{a}}_{l}+\hat{\bar{a}}_{l}\hat{\bar{h}}_{l}(0)% \right]=0. (195)

Equation (195) has the same form as Eq. (174) while the Hamiltonian \mathcal{H}(0) is replaced by an effective Hamiltonian \hat{\bar{h}}_{l}(0). Therefore, the formal solution is found as the Weyl ordering

 \displaystyle\hat{\bar{a}}_{l}(\tau)=\mathcal{W}\left\{\hat{\bar{a}}_{l}(0)% \exp\left[+i\frac{3\beta_{l}}{2}\hat{\bar{h}}_{l}(0)\tau_{l}\right]\right\}. (196)

Note that since [\hat{\bar{a}}_{l},\hat{\bar{\mathcal{H}}}_{l^{\prime}}(0)]=0, the Weyl ordering only acts partially on the Hilbert space of interest which results in a closed form solution

 \displaystyle\hat{\bar{a}}_{l}(\tau_{l})=\frac{\hat{\bar{a}}_{l}(0)e^{i\frac{3% \beta_{l}}{2}\hat{\bar{h}}_{l}(0)\tau_{l}}+e^{i\frac{3\beta_{l}}{2}\hat{\bar{h% }}_{l}(0)\tau_{l}}\hat{\bar{a}}_{l}(0)}{2\cos\left(\frac{3u_{l}^{3}\beta_{l}% \tau_{l}}{4}\right)}. (197)

At last, \hat{\bar{\mathcal{X}}}_{l}^{(0)}(t) is found by replacing \tau_{l}=\varepsilon_{l}t as

 \displaystyle\begin{split}\displaystyle\hat{\bar{\mathcal{X}}}_{l}^{(0)}(t)=% \hat{\bar{x}}_{l}^{(0)}(t,\varepsilon_{l}t)&\displaystyle=\frac{\hat{\bar{a}}_% {l}(0)e^{-i\hat{\bar{\beta}}_{l}t}+e^{-i\hat{\bar{\beta}}_{l}t}\hat{\bar{a}}_{% l}(0)}{2\cos\left(\frac{3u_{l}^{3}\beta_{l}\varepsilon_{l}}{4}t\right)}\\ &\displaystyle+\frac{\hat{\bar{a}}_{l}^{{\dagger}}(0)e^{+i\hat{\bar{\beta}}_{l% }t}+e^{+i\hat{\bar{\beta}}_{l}t}\hat{\bar{a}}^{{\dagger}}(0)}{2\cos\left(\frac% {3u_{l}^{3}\beta_{l}\varepsilon_{l}}{4}t\right)},\end{split} (198)

where \hat{\bar{\beta}}_{l}\equiv\beta_{l}\left[1-\frac{3\varepsilon_{l}}{2}\hat{% \bar{h}}_{l}(0)\right]. Plugging the expressions for \varepsilon_{l} and \hat{\bar{h}}_{l}(0), we find the explicit operator renormalization of each sector as

 \displaystyle\hat{\bar{\beta}}_{j}=\beta_{j}-\frac{3\varepsilon}{2}\omega_{j}% \left[u_{j}^{4}\hat{\bar{\mathcal{H}}}_{j}(0)+2u_{j}^{2}u_{c}^{2}\hat{\bar{% \mathcal{H}}}_{c}(0)\right], (199a) \displaystyle\hat{\bar{\beta}}_{c}=\beta_{c}-\frac{3\varepsilon}{2}\omega_{j}% \left[u_{c}^{4}\hat{\bar{\mathcal{H}}}_{c}(0)+2u_{c}^{2}u_{j}^{2}\hat{\bar{% \mathcal{H}}}_{j}(0)\right]. (199b)

Equations. (199a-199b) are symmetric under j\leftrightarrow c, implying that in the normal mode picture all modes are renormalized in the same manner. The terms proportional to u_{j,c}^{4} and u_{j,c}^{2}u_{c,j}^{2} are the self-Kerr and cross-Kerr contributions, respectively.

This analysis can be extended to the case of a Duffing oscillator coupled to multiple harmonic oscillators without further complexity, since the Hilbert spaces of the distinct normal modes do not mix to lowest order in MSPT. Consider the full Hamiltonian of our cQED system as

 \displaystyle\begin{split}\displaystyle\hat{\mathcal{H}}&\displaystyle\equiv% \frac{\omega_{j}}{4}\left(\hat{\mathcal{X}}_{j}^{2}+\hat{\mathcal{Y}}_{j}^{2}-% \frac{\varepsilon}{2}\hat{\mathcal{X}}_{j}^{4}\right)\\ &\displaystyle+\sum\limits_{n}\frac{\omega_{n}}{4}\left(\hat{\mathcal{X}}_{n}^% {2}+\hat{\mathcal{Y}}_{n}^{2}\right)+\sum\limits_{n}g_{n}\hat{\mathcal{Y}}_{j}% \hat{\mathcal{Y}}_{n},\end{split} (200)

where here we label transmon operators with j and all modes of the cavity by n. The coupling g_{n} between transmon and the modes is given as Malekakhlagh and Türeci (2016)

 \displaystyle g_{n}=\frac{1}{2}\gamma\sqrt{\chi_{j}}\sqrt{\omega_{j}\omega_{n}% }\tilde{\Phi}_{n}(x_{0}). (201)

Then, the Hamiltonian can be rewritten in a new basis that diagonalizes the quadratic part as

 \displaystyle\begin{split}\displaystyle\hat{\mathcal{H}}&\displaystyle\equiv% \frac{\beta_{j}}{4}\left(\hat{\bar{\mathcal{X}}}_{j}^{2}+\hat{\bar{\mathcal{Y}% }}_{j}^{2}\right)+\sum\limits_{n}\frac{\beta_{n}}{4}\left(\hat{\bar{\mathcal{X% }}}_{n}^{2}+\hat{\bar{\mathcal{Y}}}_{n}^{2}\right)\\ &\displaystyle-\frac{\varepsilon\omega_{j}}{8}\left(u_{j}\hat{\bar{\mathcal{X}% }}_{j}+\sum\limits_{n}u_{n}\hat{\bar{\mathcal{X}}}_{n}\right)^{4}.\end{split} (202)

The procedure to arrive at u_{j,c} and \beta_{j,c} is a generalization of the one presented under Eqs. (184a-184b).

The Heisenberg dynamics of each normal mode is then obtained as

 \displaystyle\begin{split}\displaystyle\hat{\ddot{\bar{\mathcal{X}}}}_{l}(t)+% \beta_{l}^{2}\left\{\hat{\bar{\mathcal{X}}}_{l}(t)-\varepsilon_{l}\left[u_{j}% \hat{\bar{\mathcal{X}}}_{j}(t)+\sum\limits_{n}u_{n}\hat{\bar{\mathcal{X}}}_{n}% (t)\right]^{3}\right\}=0,\end{split} (203)

where \varepsilon_{l}\equiv\frac{\omega_{j}}{\beta_{l}}u_{l}\varepsilon for l\equiv j,n. Up to lowest order in perturbation, the solution for \hat{\bar{\mathcal{X}}}_{l}^{(0)}(t) has the exact same form as Eq. (198) with operator renormalization \hat{\bar{\beta}}_{j}

 \displaystyle\hat{\bar{\beta}}_{j}=\beta_{j}-\frac{3\varepsilon}{2}\omega_{j}% \left[u_{j}^{4}\hat{\bar{\mathcal{H}}}_{j}(0)+\sum\limits_{n}2u_{j}^{2}u_{n}^{% 2}\hat{\bar{\mathcal{H}}}_{n}(0)\right], (204a) and \hat{\bar{\beta}}_{n} as \displaystyle\begin{split}\displaystyle\hat{\bar{\beta}}_{n}=\beta_{n}-\frac{3% \varepsilon}{2}\omega_{j}&\displaystyle\left[u_{n}^{4}\hat{\bar{\mathcal{H}}}_% {n}(0)+2u_{n}^{2}u_{j}^{2}\hat{\bar{\mathcal{H}}}_{j}(0)\right.\\ &\displaystyle+\left.\sum\limits_{m\neq n}2u_{n}^{2}u_{m}^{2}\hat{\bar{% \mathcal{H}}}_{m}(0)\right].\end{split} (204b)

In App. D.1, we showed that adding another time scale for the decay rate and doing MSPT up to leading order resulted in the trivial solution  (166) where the dissipation only appears as a decaying envelope. Therefore, we can immediately generalize the MSPT solutions (204a-204b) to the dissipative case where the complex pole p_{j}=-\alpha_{j}-i\beta_{j} of the transmon-like mode is corrected as

 \displaystyle\hat{\bar{p}}_{j}=p_{j}+i\frac{3\varepsilon}{2}\omega_{j}\left[u_% {j}^{4}\hat{\bar{\mathcal{H}}}_{j}(0)e^{-2\alpha_{j}t}+\sum\limits_{n}2u_{j}^{% 2}u_{n}^{2}\hat{\bar{\mathcal{H}}}_{n}(0)e^{-2\alpha_{n}t}\right], (205a) and resonator-like mode p_{n}=-\alpha_{n}-i\beta_{n} as \displaystyle\begin{split}\displaystyle\hat{\bar{p}}_{n}=p_{n}+i\frac{3% \varepsilon}{2}\omega_{j}&\displaystyle\left[u_{n}^{4}\hat{\bar{\mathcal{H}}}_% {n}(0)e^{-2\alpha_{n}t}+2u_{n}^{2}u_{j}^{2}\hat{\bar{\mathcal{H}}}_{j}(0)e^{-2% \alpha_{j}t}\right.\\ &\displaystyle+\left.\sum\limits_{m\neq n}2u_{n}^{2}u_{m}^{2}\hat{\bar{% \mathcal{H}}}_{m}(0)e^{-2\alpha_{m}t}\right].\end{split} (205b)

Then, the MSPT solution for \hat{\mathcal{X}}_{j}^{(0)}(t) is obtained as

 \displaystyle\begin{split}\displaystyle\hat{\mathcal{X}}_{j}^{(0)}(t)&% \displaystyle=u_{j}\frac{\hat{\bar{a}}_{j}(0)e^{\hat{\bar{p}}_{j}t}+e^{\hat{% \bar{p}}_{j}t}\hat{\bar{a}}_{j}(0)}{2\cos\left(\frac{3\omega_{j}}{4}u_{j}^{4}% \varepsilon te^{-2\alpha_{j}t}\right)}+H.c.\\ &\displaystyle+\sum\limits_{n}\left[u_{n}\frac{\hat{\bar{a}}_{n}(0)e^{\hat{% \bar{p}}_{n}t}+e^{\hat{\bar{p}}_{n}t}\hat{\bar{a}}_{n}(0)}{2\cos\left(\frac{3% \omega_{j}}{4}u_{n}^{4}\varepsilon te^{-2\alpha_{n}t}\right)}+H.c.\right].\end% {split} (206)

Note that if there is no coupling, u_{j}=1 and u_{n}=0 and we retrieve the MSPT solution of a free Duffing oscillator given in Eq. (179).

## Appendix E Reduced equation for the numerical solution

In this appendix, we provide the derivation for Eq. (63) based on which we did the numerical solution for the spontaneous emission problem. Substituting Eq. (38) into Eq. (36) we obtain the effective dynamics up to \mathcal{O}(\varepsilon^{2}) as

 \displaystyle\begin{split}&\displaystyle\hat{\ddot{X}}_{j}(t)+\omega_{j}^{2}% \left[1-\gamma+i\mathcal{K}_{1}(0)\right]\left[\hat{X}_{j}(t)-\varepsilon% \operatorname{Tr}_{ph}{\{\hat{\rho}_{ph}(0)\hat{\mathcal{X}}_{j}^{3}(t)\}}% \right]\\ &\displaystyle=-\int_{0}^{t}dt^{\prime}\omega_{j}^{2}\mathcal{K}_{2}(t-t^{% \prime})[\hat{X}_{j}(t^{\prime})-\varepsilon\operatorname{Tr}_{ph}{\{\hat{\rho% }_{ph}(0)\hat{\mathcal{X}}_{j}^{3}(t^{\prime})\}}].\end{split} (207)

If we are only interested in the numerical results up to linear order in \varepsilon then we can write

 \displaystyle\hat{\mathcal{X}}_{j}(t)=\hat{\mathcal{X}}_{j}^{(0)}(t)+% \varepsilon\hat{\mathcal{X}}_{j}^{(1)}(t)+\mathcal{O}(\varepsilon^{2}), (208)

and we find that

 \displaystyle\begin{split}\displaystyle\varepsilon\operatorname{Tr}_{ph}\left% \{\hat{\rho}_{ph}(0)\hat{\mathcal{X}}_{j}^{3}(t)\right\}&\displaystyle=% \varepsilon\operatorname{Tr}_{ph}\left\{\hat{\rho}_{ph}(0)\left[\hat{\mathcal{% X}}_{j}^{(0)}(t)\right]^{3}\right\}\\ &\displaystyle+\mathcal{O}\left(\varepsilon^{2}\right).\end{split} (209)

Note that in this appendix \hat{\mathcal{X}}_{j}^{(0)}(t) differs the MSPT notation in the main body and represents the linear solution. We know the exact solution for \hat{\mathcal{X}}_{j}^{(0)}(t) via Laplace transform as

 \displaystyle\begin{split}\displaystyle\hat{\mathcal{X}}_{j}^{(0)}(t)&% \displaystyle=\mathfrak{L}^{-1}\left\{\frac{s\hat{\mathcal{X}}_{j}^{(0)}(0)+% \omega_{j}\hat{\mathcal{Y}}_{j}^{(0)}(0)}{D_{j}(s)}\right\}\\ &\displaystyle+\mathfrak{L}^{-1}\left\{\frac{\sum\limits_{n}\left[a_{n}(s)\hat% {\mathcal{X}}_{n}^{(0)}(0)+b_{n}(s)\hat{\mathcal{Y}}_{n}^{(0)}(0)\right]}{D_{j% }(s)}\right\}\\ &\displaystyle=\mathfrak{L}^{-1}\left\{\frac{s\hat{X}_{j}^{(0)}(0)+\omega_{j}% \hat{Y}_{j}^{(0)}(0)}{D_{j}(s)}\right\}\otimes\hat{\mathbf{1}}_{ph}\\ &\displaystyle+\hat{\mathbf{1}}_{j}\otimes\mathfrak{L}^{-1}\left\{\frac{\sum% \limits_{n}\left[a_{n}(s)\hat{X}_{n}^{(0)}(0)+b_{n}(s)\hat{Y}_{n}^{(0)}(0)% \right]}{D_{j}(s)}\right\},\end{split} (210)

where we have employed the fact that at t=0, the Heisenberg and Schrödinger operators are the same and have the following product form

 \displaystyle\hat{\mathcal{X}}_{j}^{(0)}(0)=\hat{X}_{j}^{(0)}(0)\otimes\hat{% \mathbf{1}}_{ph}, (211a) \displaystyle\hat{\mathcal{Y}}_{j}^{(0)}(0)=\hat{Y}_{j}^{(0)}(0)\otimes\hat{% \mathbf{1}}_{ph}, (211b) \displaystyle\hat{\mathcal{Y}}_{n}^{(0)}(0)=\hat{\mathbf{1}}_{j}\otimes\hat{Y}% _{n}^{(0)}(0), (211c) \displaystyle\hat{\mathcal{X}}_{n}^{(0)}(0)=\hat{\mathbf{1}}_{j}\otimes\hat{X}% _{n}^{(0)}(0). (211d)

The coefficients a_{n}(s) and b_{n}(s) can be found from the circuit elements and are proportional to light-matter coupling g_{n}. However, for the argument that we are are trying to make, it is sufficient to keep them in general form.

Note that equation (210) can be written formally as

 \displaystyle\hat{\mathcal{X}}_{j}^{(0)}(t)=\hat{X}_{j}^{(0)}(t)\otimes\hat{% \mathbf{1}}_{ph}+\hat{\mathbf{1}}_{j}\otimes\hat{X}_{j,ph}(t). (212)

Therefore, \left[\hat{\mathcal{X}}_{j}^{(0)}(t)\right]^{3} is found as

 \displaystyle\begin{split}&\displaystyle\left[\hat{\mathcal{X}}_{j}^{(0)}(t)% \right]^{3}=\left[\hat{X}_{j}^{(0)}(t)\right]^{3}\otimes\hat{\mathbf{1}}_{ph}+% \hat{\mathbf{1}}_{j}\otimes\hat{X}_{j,ph}^{3}(t)\\ &\displaystyle+3\left\{\left[\hat{X}_{j}^{(0)}(t)\right]^{2}\otimes\hat{X}_{j,% ph}(t)+\hat{X}_{j}^{(0)}(t)\otimes\hat{X}_{j,ph}^{2}(t)\right\}.\end{split} (213)

Finally, we have to take the partial trace with respect to the photonic sector. For the initial density matrix \hat{\rho}_{ph}(0)=\ket{0}_{ph}\bra{0}_{ph}

 \displaystyle\braket{\hat{X}_{j,ph}(t)}_{ph}=\braket{\hat{X}_{j,ph}^{3}(t)}_{% ph}=0. (214)

The only nonzero expectation values in \braket{\hat{\mathcal{X}}_{j,ph}^{2}(t)}_{ph} are \braket{\hat{X}_{n}^{2}(0)}_{ph}=\braket{\hat{Y}_{n}^{2}(0)}_{ph}=1. Therefore, the partial trace over the cubic nonlinearity takes the form

 \displaystyle\begin{split}&\displaystyle\operatorname{Tr}_{ph}\left\{\hat{\rho% }_{ph}(0)\left[\hat{\mathcal{X}}_{j}^{(0)}(t)\right]^{3}\right\}=\\ &\displaystyle\left[\hat{X}_{j}^{(0)}(t)\right]^{3}+3\mathfrak{L}^{-1}\left\{% \frac{\sum\limits_{n}\left[a_{n}^{2}(s)+b_{n}^{2}(s)\right]}{D_{j}(s)}\right\}% \hat{X}_{j}^{(0)}(t).\end{split} (215)

The first term is the reduced transmon operator cubed. The second term is the sum over vacuum fluctuations of the resonator modes. Neglecting these vacuum expectation values we can write

 \displaystyle\begin{split}\displaystyle\operatorname{Tr}_{ph}\left\{\hat{\rho}% _{ph}(0)\left[\hat{\mathcal{X}}_{j}^{(0)}(t)\right]^{3}\right\}&\displaystyle% \approx\left[\hat{X}_{j}^{(0)}(t)\right]^{3}\\ &\displaystyle=\hat{X}_{j}^{3}(t)+\mathcal{O}(\varepsilon^{2}),\end{split} (216)

Substituting Eq. (216) into Eq. (207) gives

 \displaystyle\begin{split}&\displaystyle\hat{\ddot{X}}_{j}(t)+\omega_{j}^{2}% \left[1-\gamma+i\mathcal{K}_{1}(0)\right]\left[\hat{X}_{j}(t)-\varepsilon\left% [\hat{X}_{j}(t)\right]^{3}\right]\\ &\displaystyle=-\int_{0}^{t}dt^{\prime}\omega_{j}^{2}\mathcal{K}_{2}(t-t^{% \prime})\left[\hat{X}_{j}(t^{\prime})-\varepsilon\left[\hat{X}_{j}(t^{\prime})% \right]^{3}\right]+\mathcal{O}(\varepsilon^{2}).\end{split} (217)

## References

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters