Quantum versus classical switching dynamics of driven-dissipative Kerr resonators

TL;DR

This paper investigates the driven-dissipative dynamics of Kerr resonators in the mesoscopic regime, highlighting the importance of quantum effects and the limitations of classical models in describing switching behavior.

Contribution

It provides a first-principles quantum analysis of switching dynamics in Kerr resonators, revealing the inadequacy of classical escape models and introducing a universal crossover parameter.

Findings

Significant deviations from steady-state predictions in nonlinear resonance lineshapes.

Emergence of a slow switching time scale longer than linear decay rates.

Quantum treatment necessary for accurate switching rate predictions.

Abstract

We report a first-principles study of the driven dissipative dynamics for Kerr oscillators in the mesoscopic regime. This regime is characterized by large Kerr nonlinearity, realized here using the nonlinear kinetic inductance of a large array of Josephson junctions. The experimentally measured nonlinear resonance lineshapes of the junction array modes show significant deviations from steady-state numerical predictions, and necessitate time-dependent numerical simulations indicative of strong measurement-induced dephasing in the system arising from the large cross-Kerr effect between array modes. Analytical and numerical calculations of switching rate corroborate this by showing the emergence of a slow time scale, which is much longer than the linear decay rate and is set by fluctuation-induced switching times in the bistable regime. Furthermore, our analysis shows that the usual quantum-activated escape treatment is inadequate for prediction of the switching rates at large frequency shifts caused by strong nonlinearities, necessitating a quantum treatment that utilizes the full system Liouvillian. Based on our analysis, we identify a universal crossover parameter that delineates the regimes of validity of semiclassical and quantum descriptions, respectively. Our work shows how dynamical switching effects in strongly nonlinear systems provide a platform to study quantum-to-classical transitions.

Equations

$$H_k={\omega }_k{a}_k^{†}{a}_k-K_k(a_k^{†}{a}_k{)}^2+{\epsilon }_k(t)a_k^{†}+{\epsilon }_k^{\ast }(t)a_k,$$

$$\displaystyle\dot{\rho}=-i[H_{k},\rho]+\kappa_{k}\Big{[}a_{k}\rho a_{k}^{\dagger}-\frac{1}{2}(a_{k}^{\dagger}a_{k}\rho+\rho a_{k}^{\dagger}a_{k})\Big{]}.$$

$${\omega }_k\to {\omega }_k-\sum _l{K}_{kl},$$

$${\Delta }_0=-K_{01}⟨a_1^{†}{a}_1{⟩}_s,$$

$$\begin{split}d{\rho_{Z}}&=-idt\,[H_{1}+H_{0}+H_{c},\rho_{Z}]\\ &\quad+dt\,\kappa_{1}\big{[}a_{1}\rho_{Z}a_{1}^{\dagger}-\frac{1}{2}(a_{1}^{\dagger}a\rho_{Z}+\rho_{Z}a_{1}^{\dagger}a_{1})\big{]}\\ &\quad+dt\,\kappa_{0}\big{[}a_{0}\rho_{Z}a_{0}^{\dagger}-\frac{1}{2}(a_{0}^{\dagger}a_{0}\rho_{Z}+\rho_{Z}a_{0}^{\dagger}a_{0})\big{]}\\ &\quad+\sqrt{\kappa_{0}}\Big{[}dZa_{0}\rho_{Z}+dZ^{*}\rho_{Z}a_{0}^{\dagger}\\ &\quad-\text{Tr}(dZa_{0}\rho_{Z}+dZ^{*}\rho_{Z}a_{0}^{\dagger})\,\rho_{Z}\Big{]}\end{split}$$

$$H_c=-K_{01}{a}_1^{†}{a}_1{a}_0^{†}{a}_0^{}$$

$$J_T(t)=$$

$$\displaystyle+\frac{1}{2}\Big{\{}i\sqrt{\kappa_{0}}\,\text{Tr}[\rho_{Z}(a_{0}-a_{0}^{\dagger})]+dW_{b}/dt\Big{\}}^{2},$$

$$\displaystyle\dot{a}_{1}=-i\big{[}a_{1}\,,\,H_{1}\big{]}-\frac{\kappa_{1}}{2}a_{1}+\sqrt{\kappa_{1}}a_{1}^{\rm in}(t),$$

$$\dot{\alpha }=-i\Delta \alpha +2iK|\alpha {|}^2\alpha -i\epsilon -\frac{\kappa }{2}\alpha +\sqrt{\kappa }\zeta (t),$$

$$\rho (t)=\sum _{\lambda }{c}_{\lambda }{e}^{\lambda t}{\rho }_{\lambda }$$

$${\tilde{\rho }}_{\chi }={\rho }_{\chi }-\frac{⟨{\rho }_s,{\rho }_{\chi }⟩}{⟨{\rho }_s,{\rho }_s⟩}{\rho }_s,$$

$$\rho (t)={\beta }_0(t){\rho }_s+{\beta }_1(t)\tilde{\rho }+\sigma (t),$$

$${\dot{\beta }}_0(t)=⟨{\rho }_s,\mathcal{L}[\rho (t)]⟩.$$

$$\displaystyle\lambda_{e}={\langle\,\rho_{s},\mathcal{L}[\rho(t=0)]\,\rangle}\Bigg{(}{1-\frac{\langle\rho_{s},\rho(t=0)\rangle}{\langle\rho_{s},\rho_{s}\rangle}}\Bigg{)}^{-1}.$$

$$F(z)=$$

$$-2iK(z+e^{-i\phi }{\alpha }_0^{\ast })(z+e^{i\phi }{\alpha }_0{)}^2.$$

$$\text{Im}[F(x)]=\dot{p}=-\frac{dU(x)}{dx}.$$

$$U(x)=-\int \text{Im}[F(x)]dx,$$

$$\displaystyle\lambda_{e}=\gamma_{0}\exp\Big{(}-\frac{\Delta U}{\kappa}\,\Big{)},$$

$$T_{\gamma }=\frac{\hslash {\gamma }_0}{2\pi k_B}\left(\sqrt{\frac{{\kappa }^2}{4{\gamma }_0^2}+1}-\frac{\kappa }{2{\gamma }_0}\right),$$

$$\xi \equiv \frac{T_{\gamma }}{T_{\kappa }},$$

$${\mathcal{L}}_{array}$$

$$-\frac{1}{2L_J}({\varphi }_n-{\varphi }_{n+1}{)}^2.$$

$${\mathcal{L}}_{end}=\frac{C_s}{2}{\dot{\varphi }}_1^2+\frac{C_g}{2}{\dot{\varphi }}_1^2+\frac{C_e}{2}{\dot{\varphi }}_{N+1}^2,$$

$$\mathcal{L}$$

$$={\dot{\vec{\varphi }}}^T\frac{\mathbb{C}}{2}\dot{\vec{\varphi }}-{\vec{\varphi }}^T\frac{\mathbb{L}}{2}\vec{\varphi },$$

$$C_k={\vec{v}}_k^T\mathbb{C}{\vec{v}}_k,$$

$$\mathcal{L}=\sum _k\frac{C_k}{2}{\dot{\varphi }}_k(t{)}^2-\frac{1}{2L_k}{\varphi }_k(t{)}^2.$$

$$q_k=C_k{\dot{\varphi }}_k^2,$$

$$H=\sum _k\frac{{\hat{q}}_k^2}{2C_k}+\frac{{\hat{\varphi }}_k^2}{2L_k},$$

$$H=\sum _k{\omega }_k{a}_k^{†}{a}_k^{}.$$

Full text

††thanks: chanders@phys.ethz.ch††thanks: Present address: IBM T.J. Watson Research Center, Yorktown Heights, NY 10598, USA††thanks: Present address: Physikalisches Institut, Karlsruhe Institute of Technology, 76131 Karlsruhe, Germany

Nonlinear optics spans a broad class of phenomena that involve light-induced variation of optical properties of a system. Interestingly, the nonlinearity in optical systems often originates from an inherently quantum mechanical process, but the description of the resulting output can be either classical or quantum, depending on the system under consideration and intensity of the light fields. For instance, nonlinear crystals implementing frequency mixing and stimulated scattering Boyd (2003) are described using classical descriptions rooted in average nonlinear material susceptibilities. On the other hand, a new class of effects and non-classical states emerge on quantizing the light fields, which are dealt with in the framework of quantum optics Carmichael (2009). In order to distinguish between these two descriptions, it is important to understand the transition from quantum to classical dynamics, especially given the burgeoning presence of nonlinear systems in applications triggered by quantum information processing Nielsen and Chuang (2010); Braunstein and van Loock (2005).

An example of a nonlinear phenomenon that is of importance in both quantum and classical optics is the Kerr effect, which changes the optical properties in proportion to the intensity of the incident field. Under a coherent drive and single photon loss, the Hamiltonian of the $k$th mode of a physical system subject to the Kerr effect is given by ($\hbar=1$)

[TABLE]

and the dynamics of the mode, under the assumption of Markovian decay, is described by the master equation

[TABLE]

Here, $\omega_{k}$ is the bare resonance frequency mode, $K_{k}$ is the Kerr coefficient responsible for shifting the mode frequency, $\varepsilon_{k}(t)=\varepsilon_{k}e^{-i\omega_{d}t}$ is the amplitude of a coherent drive with angular frequency $\omega_{d}$ and $\kappa_{k}$ is the single-photon decay rate of the mode. We note that the Hamiltonian above assumes that all nonlinear modes are independent; in principle, there can be additional cross-Kerr coupling between different $k$ modes of the system. This driven-dissipative Kerr resonator is ubiquitous in physical systems spanning fiber optics Agarwal (2006), superconducting quantum circuits with Josephson junctions (JJs) Yurke (1987); Castellanos-Beltran et al. (2008); Zueco et al. (2012); Bourassa et al. (2012); Leib et al. (2012); Eichler and Wallraff (2014); Bhupathi et al. (2016), optomechanical systems Bose et al. (1997); Ludwig et al. (2012) and atomic ensembles Gupta et al. (2007).

In particular, JJ-based superconducting circuits naturally realize strong Kerr nonlinearities at the single photon level, which have been extensively utilized both for fundamental quantum optical studies Bozyigit et al. (2011); Kirchmair et al. (2013) and quantum information-inspired applications Blais et al. (2007); Córcoles et al. (2015); Kelly et al. (2015). This breadth of applications is enabled by the highly flexible nature of the nonlinearity realized in superconducting circuits, which can be tuned either in-situ through the flux-tunable SQUID-based designs Bell et al. (2012) or ex-situ through appropriate selection of junction parameters; the latter is usually accomplished by designing junctions with different ratios of Josephson energy and charging energy or, alternatively, by using arrays of junctions Eichler and Wallraff (2014); Masluk et al. (2012); Weißl et al. (2015).

Figure 1 maps the different nonlinearity regimes readily accessible with JJ circuits. For nonlinearities much smaller than the mode linewidth ($K_{k}\ll\kappa_{k}$), the dynamics can be entirely captured by that of a driven classical Duffing oscillator Dykman (2012). This is the domain of the Josephson parametric amplifiers (JPA) Ristè et al. (2012); Johnson et al. (2012); Lin et al. (2013) and Josephson bifurcation amplifiers (JBA) Mallet et al. (2009); Vijay et al. (2009) that are routinely employed for quantum-limited measurements Murch et al. (2013); Sun et al. (2014). Of particular interest in this regime is the bistable behavior occurring when a lossy Kerr resonator is strongly driven. In this bistable regime, the Kerr resonator may switch between two stable states corresponding, respectively, to small and large photon populations of the driven resonator Muppalla et al. (2018). For weak nonlinearity, this switching behavior is well-described by a semi-classical quantum activation theory Vijay et al. (2009); Dykman and Krivoglaz (1980); Dykman (2012). On the other hand, for large nonlinearity ($K_{k}\gg\kappa_{k}$), the system enters the photon blockade regime in which the nonlinearity-induced frequency shift limits the number of photons in the oscillator under a drive of fixed frequency. The nonlinear oscillator then becomes an effective two- or ‘few’-level system; this is the so-called transmon regime Koch et al. (2007). In contrast to the semi-classical picture, oscillators in this regime do not exhibit any bistability or hysteresis and a full quantum description is necessary to describe their dynamics.

A thorough understanding of fluctuation-induced switching rates in the cross-over region is important for understanding and optimizing Kerr oscillators, both from fundamental physics and application points of view. For instance, nonlinear oscillators exhibit scale invariant behavior near bistability; also, the time needed to switch from one metastable state to the other limits the qubit measurement time in JBAs Mallet et al. (2009); Vijay et al. (2009). However, the crossover of driven-dissipative dynamics from (semi)-classical to the quantum regime remains poorly understood. In this work, we perform a detailed investigation of this intermediate mesoscopic regime between semi-classical bistability and the transmon regimes for driven Kerr oscillators. The workhorse of our studies is a Kerr oscillator in the strongly nonlinear mesoscopic regime, based on a superinductance formed from an array of Josephson junctions Manucharyan (2012). The eigenmodes of the array form highly nonlinear oscillator modes, where the nonlinear shift per photon is larger than the natural oscillator linewidths. Specifically, we use the time scale, here labelled $\tau$, associated with decay into the steady state of the Kerr resonator as a benchmark to delineate semi-classical and quantum dynamics. We report the results of an experiment with a nonlinear resonator realized with a superconducting quantum circuit. The results are simulated using a stochastic master equation and we find signatures of oscillator relaxation time $\tau$ much longer than the intrinsic decay time $1/\kappa$, which motivates further theoretical investigations. We go beyond a linear treatment and present both numerical and analytical calculations of switching rates in this system, considering the situation where the oscillation amplitude is locked in one of the two metastable dynamical states in the presence of a strong drive. Our theoretical studies indicate a breakdown of usual semi-classical treatments that describe oscillator decay primarily as quantum activation in a metapotential or, alternatively, by thermal activation. Instead, in the mesoscopic regime, a quantum treatment is essential to capture the relaxation timsescale of the oscillator. We characterize this transition from a semi-classical to quantum description by introducing a crossover parameter, $\xi=T_{\gamma}/T_{\kappa}$, as a ratio of a temperature associated with tunneling-induced escape $T_{\gamma}$ and effective temperature associated with fluctuations seen by the oscillator $T_{\kappa}$. When $\xi>1$, quantum effects introduce a new decay channel and fluctuation-induced activation is inadequate to describe the switching dynamics in this regime.

This paper is organized as follows: Section I introduces our experimental implementation of a Kerr resonator in the mesoscopic regime realized by an array of Josephson junctions Masluk et al. (2012). In Sec. II, we present the theoretical model describing the ground state properties and Kerr coupling between the distributed modes of the array. By doing a first principle quantization of the array, we confirm that the Kerr coefficients for its distributed modes lie in the regime of interest, $K_{k}/\kappa_{k}\approx 2-5$. We then describe how stochastic master equation simulations accurately predict the experimentally observed nonlinear resonance lineshape for an eigenmode of the array in the presence of an external drive. In Secs. III and IV, we perform detailed numerical and analytical investigations of fluctuation-induced switching rates in order to understand the nonlinear effects observed in the experiment. Sec. V concludes the paper with a discussion of the main results and provides an outlook for future theoretical and experimental studies. Additional numerical and analytical results are included in Appendicies A and B.

I Experiment with Josephson junction arrays

Figure 2 depicts the experimental system realizing a mesoscopic Kerr oscillator, which consists of an array of 80 Josephson junctions capacitively coupled to a transmission line in a hanger geometry Masluk et al. (2012). The sample is mounted inside a copper box in a dilution fridge at a temperature of 15 mK. Both coherent driving and measurement of the resonator are performed through the transmission line. The transmission signal at the output port 2 is amplified by a HEMT amplifier before being demodulated and recorded.

The resonator has a fundamental frequency of $\omega_{0}/2\pi=4.357$ GHz, with an internal quality factor of 37,000 and external quality factor of 5,000. Larger internal quality factors were observed for other samples Masluk et al. (2012) but, in all cases, the quality factor was dominated by the external coupling. The experimental setup is constrained by circulators and a low-pass filter to measure signals with frequencies ${<}\,12$ GHz. As a result, only the fundamental mode of the Kerr resonator can be directly probed. Leveraging the nonlinearity-induced mode-mode coupling, it is however possible to probe higher frequency modes. Indeed, a continuous drive applied to the array at the frequencies of its higher modes leads to a shift of the fundamental mode frequency. This shift arises due to the cross-Kerr coupling represented by the interaction term $\sum_{k\neq l}-K_{kl}a_{k}^{\dagger}a_{k}a_{l}^{\dagger}a_{l}$ where $K_{kl}$ is the cross-Kerr coupling between array modes $k$ and $l$ while $a_{k}^{(\dagger)}$ destroys (creates) an excitation in mode $k$. This mode-mode interaction leads to a dispersive readout mechanism for the modes that are above the experimental frequency cutoff (see Appendix A). This technique was already used in Ref. Masluk et al. (2012) to map the full dispersion relation of the distributed modes of the array. In this way, the first mode was identified to be at $\omega_{1}/2\pi=11.9$ GHz.

To probe the photon number in this first mode as a function of the drive detuning, a continuous pump tone is applied close to $\omega_{1}$. Simultaneously, the transmission spectrum of the fundamental mode is measured by sweeping the frequency of a weak continuous probe tone. The frequency at which the probe is maximally reflected indicates the resonance of the fundamental mode. Figure 3 shows the experimentally measured frequency shift (dots) obtained using this method for different pump powers and as a function of the detuning of the pump drive. The driven-dissipative Kerr resonator always shows a single resonance peak. Similarly, we note that the measured frequency shifts show abrupt changes between a low and a high value, with no values in between the two. The lack of intermediate values for the frequency shift indicates that no switching occurs between two bistable states during the time of the experiment set by the the measurement time-scale, $t_{m}=1\,\mu$s and averaged 1000 times. We present a detailed theoretical description of the experimentally observed resonance lineshapes in the next section.

II Theoretical model

Since only the fundamental mode is directly measurable experimentally, the Kerr coefficients are inferred from a theoretical model of the distributed junction array. To this end, we write the Lagrangian for the linear part of the array in terms of the system capacitances and Josephson inductance ($L_{J}$), as shown in Fig. 2(c) Devoret (1995). Following a quantization of these modes, we perturbatively include the nonlinear contributions from the Josephson junctions (see Appendix A for details) Bourassa et al. (2012); Weißl et al. (2015). From room-temperature resistance measurements, we infer an inductance per junction of $L_{J}=1.9$ nH and an effective array plasma frequency (highest mode frequency) of 18.2 GHz. The effective array plasma frequency is not equal to the single junction plasma frequency $\omega_{p}=1/\sqrt{C_{J}L_{J}}$, as would be the case if the junctions were purely linear elements. We moreover infer that the system capacitances illustrated in to have values $C_{0}=0.066$ fF, $C_{J}=26.54$ fF, $C_{g}=10.4$ fF, $C_{s}=3$ fF and $C_{e}=10.84$ fF close to the design parameters found by simulations Masluk et al. (2012). With these parameters, we estimate the self-Kerr coefficient of the fundamental mode to be $K_{0}/2\pi=0.5$ MHz, and that for the first mode to be $K_{1}/2\pi=5.7$ MHz. The latter is almost a factor of two bigger than the corresponding linewidth $\kappa_{1}/2\pi=2.9$ MHz, and hence lies in the intermediate (or mesoscopic) regime identified in Fig. 1.

Reintroducing the junction nonlinearity in the analysis, we find that the mode frequencies get shifted by self-Kerr and cross-Kerr couplings. The $k$th mode’s frequency, $\omega_{k}$, is shifted such that the new effective frequency becomes

[TABLE]

where $K_{kl}$ is the cross-Kerr contribution between modes $k$ and $l$ (see Appendix A). Although each $K_{kl}$ is a factor of $10^{2}$–$10^{3}$ smaller than $\omega_{k}$ (on the order of a few MHz), the cumulative effect of all the modes results in significant shifts of the eigenfrequencies of the array by up to a GHz compared to the bare mode frequencies. Moreover, as a result of the cross-Kerr couplings, when driving the Kerr-resonator, the average photon number, $\langle a_{1}^{\dagger}a_{1}\rangle$, increases leading to a measurable frequency shift of the fundamental mode.

As a first, simple, approach to reproducing the data of Fig. 3, the steady-state photon number, $\langle a_{1}^{\dagger}a_{1}\rangle_{s}=\text{Tr}(\rho_{s}a_{1}^{\dagger}a_{1})$, is numerically computed using the single-mode master equation Eq. (2). From this, the frequency shift of the fundamental mode due to photon population in the first mode is then estimated to be

[TABLE]

where $K_{01}=4\sqrt{K_{0}K_{1}}$ denotes the cross-Kerr coupling between the two modes Bourassa et al. (2012). The result of this calculation corresponds to the dashed lines in Fig. 3. Evidently, this steady-state result does not match the experimental data (dots). Indeed only some features of the experiment are captured by this single-mode steady-state calculation. For instance, the detunings are approximately correct. On the other hand, the shape of the frequency shift as a function of the detuning is very poorly reproduced. In particular, the shoulders to the left of the maxima of the shift are not observed experimentally.

In order to theoretically reproduce the data, we now turn to an approach that resembles the experimental situation more closely. First, rather than computing the steady-state photon population, we compute the full time-dependent response of the transmitted power along the transmission line coupled to the array. As in the experiment, this response is integrated over a finite measurement time $t_{m}$ and averaged over 1000 realizations. Any dynamics that occur on a time-scale much slower than $t_{m}$, such as slow switching between two bistable states, are therefore neither resolved in the measurement nor in the simulations. Second, while the cross-Kerr interaction was used in computing the expected frequency shift in Eq. (4), we now add this interaction directly to the system Hamiltonian. This allows for capturing the measurement-backaction in the form of measurement-induced dephasing of the first mode by the probe tone. The addition of measurement-induced dephasing increases the effective linewidth of the nonlinear mode, thereby leading to a reduction in the number of photons in the mode and, consequently, to a smaller observed frequency shift.

In order to obtain a description as close to the actual experiment, we numerically solve the full two-mode Hamiltonian including both the mesoscopic Kerr mode ($a_{1}$) and probe mode ($a_{0}$) with the stochastic master equation,

[TABLE]

valid for heterodyne detection of the fundamental mode $a_{0}$ Wiseman and Milburn (2009). In this expression, $H_{k}$ is Hamiltonian given by Eq. (1), while

[TABLE]

is the relevant cross-Kerr interaction. In Eq. (5), the density matrix $\rho_{Z}$ is conditioned on the result of a heterodyne measurement where $dZ=dW_{a}+idW_{b}$, with $dW_{i}$ a stochastic Wiener processes with $\langle dW_{i}\rangle=0$ and $\langle dW_{i}^{2}\rangle=dt$111The subscripts of $dW_{i}$ indicates solely that these are different independent Wiener processes.. The last term of Eq. (5) is associated to the homodyne measurement record which is used to update our description of the state of the system $\rho_{Z}$. From the numerical integration of Eq. (5), the transmitted power can then be calculated as

[TABLE]

which, averaged over time $t_{m}$, yields the measured signal. The value $\Delta_{0}$ corresponding to the maximum transmission signal, $J_{T}$, gives the frequency shift plotted in Fig. 3. The simulation were performed using standard numerical integration for the deterministic part of the equation and after each integration step, $\delta t$, a random number is drawn from a Gaussian distribution with variance $\delta t$ for each stochastic term $dW_{i}$ such that the stochastic part of $d\rho_{Z}$ can be readily calculated.

In Fig. 3, the numerical simulations using the parameters of the experiment are presented as solid lines, which are in excellent agreement with the experimental data. The difference between the steady-state values used in Eq. (4) and the stochastic master equation simulations is mainly due to the facts that the latter (i) includes the full two-mode interaction and, as a result, the measurement backaction, and (ii) explicitly takes into account the finite measurement time $t_{m}$. Therefore, dynamics on time-scales much longer than the characteristic measurement time, which are not probed by the experiment, are also not observed in the stochastic master equation simulation. The most predominant feature observed from the stochastic simulations, and not the simpler single-mode steady-state calculations, is the lack of intermediate values in the abrupt transition from low photon number to high photon number state when changing the drive detuning, i.e., the “shoulders” in the steady-state photon number. For strong Kerr nonlinearities the energy levels of a Kerr oscillator are $\omega_{k}-K_{k}n$ and can be well-resolved. The “shoulders” at large detunings are, thus, signatures of switching events at corresponding photon-number-selective transitions. These intermediate values are not resolved in the experiment due to the finite measurement time and this fact indicates a time-scale of potential switching dynamics much longer than the measurement time $t_{m}$ and comparably or even longer than the total experimental time of $t_{m}$ times the number of averages. To fully capture the dynamics of the experiment, a two-mode model is needed, however, the contributions from the two-mode interaction do not significantly contribute to the switching time-scale. Rather, the effects of the cross-Kerr interaction is to account for the measurement-induced backaction. This backaction affects the detuning and drive strength for which the increased time-scales appear. Thus, the physics of the long time-scales can be understood from only a single mode model. In the next two sections, we confirm this by performing a detailed investigation of the switching rates for a single mode Kerr resonator with parameters in the mesoscopic regime, and explicitly comparing the predictions from semi-classical and quantum treatments.

III Numerical treatment of the switching rates

Since both the experiment and stochastic master equation simulations show a qualitatively different behavior than the single-mode steady-state master equation, a more complete quantum description is necessary to quantitatively understand the switching dynamics of Kerr resonator in large-nonlinearity regimes. In particular, as evident from Fig. 3, the lack of intermediate photon number values in the bistable regime indicates that there is no dynamical switching between the two bistable states within the time-scale of the experiment. In this section, we perform numerical simulations of the steady-state photon number and switching rates of the Kerr oscillator introduced in the previous section. Specifically, in order to identify the appropriate description for the nonlinear oscillator in mesoscopic regime, we compare the results obtained from semi-classical simulations with that from quantum master equation simulations.

We begin with a generic description of the steady-state response of a nonlinear oscillator. It is convenient to delineate the time-averaged response into three regions, as sketched in Fig. 4 for an oscillator with $K=\kappa$ (this corresponds to the mesoscopic regime depicted as the blue region in Fig. 1). The thick full line is a representative curve obtained using a master equation simulation, while the narrow dashed line is obtained from a semi-classical analysis with the same parameters. In region I, the resonator relaxes to a state with a low photon number on a time scale of $T_{\kappa}=1/\kappa$. There may be two stable classical solutions in this region, as indicated by the dashed line in Fig. 4, but the fluctuations associated with the low photon-number state are not sufficient to bring the system to the high photon number state. For the lower branch solution of this region, the semi-classical and the full quantum treatment give identical results for both the relaxation time into the steady-state and for the steady-state photon number. Indeed, in this regime, only few photons are present and, as a consequence, the nonlinear effects play only a minor role. Region III is conceptually comparable to region I, except that the system relaxes into a high photon number state. In region II, the dynamics are such that the photon number initially latches to a low photon-number state, but after some time it jumps to, and continues to fluctuate around, a high photon-number state Siddiqi et al. (2005); Santori et al. (2014). This switching between the high and low photon number states continues and we refer to the time-scale of this dynamics as the switching time, $\tau$. We henceforth focus on region II, and study switching times as a function of the ratio $K/\kappa_{1}$. In particular we are interested in quantifying the difference between a weakly nonlinearly Kerr resonator, which we expect to behave classically, and a highly nonlinear Kerr resonator where quantum fluctuations are expected to play a larger role.

We first consider a semi-classical description which we obtain from the Heisenberg-Langevin equation satisfied by the first mode,

[TABLE]

where $H_{1}$ is given by Eq. (1). Here the input field, $a_{1}^{\rm in}(t)$, satisfying the commutation relation $[a_{k}^{\rm in}(t),a_{k^{\prime}}^{\rm in\dagger}(t^{\prime})]=\delta_{kk^{\prime}}\delta(t-t^{\prime})$, accounts only for the quantum fluctuations induced by the environment. The driven-dissipative Kerr resonator is most conveniently treated in a frame rotating at the drive frequency, $\omega_{d}$, such that we make the replacements $\varepsilon_{1}(t)\rightarrow\varepsilon_{1}$ and $\omega_{1}\rightarrow\Delta_{1}=\omega_{1}-\omega_{d}$ in the Hamiltonian.

We first attempt to solve the time-evolution using a semi-classical trajectory approach Armen and Mabuchi (2006), which replaces the system field operators with complex numbers, $\langle a\rangle\rightarrow\alpha$. This replacements reduces Eq. (8) to an equation of motion for the phase-space variable $\alpha$ (see also Appendix B for a deterministic classical treatment). However, when replacing the input field with a stochastic variable to account for the quantum noise associated with $a_{1}^{\rm in}(t)$, this approach becomes only approximate in the presence of nonlinear mixing terms in the Hamiltonian, since it does not account for up/down-converted quantum noise due to mode mixing. The resulting semi-classical stochastic equation of motion can then be written as,

[TABLE]

corresponding to an equation for the coherent state amplitude of the system, and where $\zeta(t)$ is stochastic Wiener process that models the quantum vacuum noise 222Here, $\zeta(t)=(dW_{a}(t)/dt+i\,dW_{b}(t)/dt)/\sqrt{2}$ represents a stochastic Wiener process with $\langle dW_{i}\rangle=0$ and $\langle dW_{i}^{2}\rangle=dt$ that, on average, corresponds to the input field being the quantum vacuum state.. Here, we have suppressed the subscript 1 both for brevity and to underscore the generality of this treatment. A particular realisation of $\zeta(t)$ is referred to as a semi-classical trajectory. In the bistable regime of the oscillator, a trajectory initially latches to a low photon-number state, and after some time jumps to a high photon-number state Siddiqi et al. (2005); Santori et al. (2014). As the switching continues with a waiting time between the successive switches approximately Poisson distributed, we use an exponential function to fit the average over many trajectories to extract the time scale to reach steady-state.

To understand the quantum effects that arise in the switching dynamics of the Kerr resonator, we next perform numerical master simulations, Eq. (2), of the driven system with the Hamiltonian, Eq. (1), and analyze the relaxation towards the steady-state. From these simulations, we observe an oscillatory behavior of the photon number that relaxes on a time scale $T_{\kappa}$ followed by an exponential relaxation towards the steady state. The relaxation time scale, $\tau$, can therefore be readily extracted from the exponent of exponential decay obtained using master equation simulations.

In Fig. 5, we compare the results obtained from both simulations described above, for a strongly nonlinear Kerr oscillator with $K=2\kappa$ (solid lines) and for a weakly nonlinear Kerr oscillator with $K=0.2\kappa$ (dashed lines). The results of the master equation simulation are shown in the main panel, while those of the semi-classical trajectories are shown in the inset. In both cases, the relaxation time is shown in blue (left axis) and the average photon number in red (right axis). First, it is clear from both the classical and quantum numerical treatments that the switching time (blue lines) can significantly exceed the linear decay time $T_{\kappa}=1/\kappa$ near the bistability. Further, for weak nonlinearity $K=0.2\kappa$, both the quantum and the semi-classical approaches show a sharp rise in $\tau$ close to the detuning where the steady state maximum photon number reaches its maximum (compare dashed lines between the main figure and the inset). This behavior is a generic feature of weakly nonlinear systems. Semi-classically, it can be understood as being caused by the switching rates between two stable states becoming equal.

For larger nonlinearity, $K=2\kappa$, the discrepancy between the full quantum simulations and semi-classical numerics is much more pronounced (compare solid lines between the main figure and the inset). Specifically, a much wider distribution of time scales is predicted by the quantum treatment in the strong nonlinearity regime, which the semi-classical approach completely fails to capture. This observation is not surprising and is reinforced by the fact that the steady-state in this regime shows a negative $Q$-parameter, $Q=(\langle(\Delta a_{1}^{\dagger}a_{1})^{2}\rangle-\langle a_{1}^{\dagger}a_{1}\rangle)/\langle a_{1}^{\dagger}a_{1}\rangle\approx{-}0.4$, which indicates sub-Poissonian statistics. Similarly, the steady-state Wigner function displays negative values (not shown), a clear sign of the non-classical nature of the state of the system. Numerical simulations with varying $K/\kappa$ further confirm that the average behavior observed in Fig. 5 gradually transitions from a situation where the semi-classical treatment matches well with full quantum simulation for small $K$, towards a larger discrepancy when $K$ is increased (not shown).

The numerical results presented here can now be compared to the results of Fig. 3. There, the steady-state behavior showed intermediate values in the bistable regime, while the experimental data was only correctly reproduced using a two-mode stochastic quantum analysis that take into account the finite measurement time. This is consistent with the emergence of a very slow time-scale, which can be significantly longer than the linear decay rate $\kappa$ of the oscillator; since dynamics on very long time-scales were not probed by the experiment, it invalidates a steady-state analysis of the data. Moreover, for sufficiently large nonlinearities ($K>\kappa$) as relevant to the experiment with JJ arrays, this slow time-scale cannot be attributed only to semi-classical switching dynamics. In the next section, we present analytical studies highlighting the differences between the semi-classical picture and the full quantum calculations of the switching rate in the large nonlinearity regime ($K/\kappa\geq 1$), while comparing them against the numerical results obtained in this section.

IV Analytical treatment of the switching rates

IV.1 Quantum calculation

In the last section, we saw that the increase in the switching time in this regime is not captured by a semi-classical treatment. Therefore, we now consider an analytical description that takes into account the full nonlinear quantum dynamics.

To this end, we consider the Liouvillian describing the time evolution of the driven-dissipative single-mode Kerr-oscillator in the mesoscopic regime [Eqs. (1) and (2)]. In a Hilbert-space of dimension $N$ ($N$ being the number of Fock states included in the calculation), the Liouvillian has $N^{2}$ eigenvalues. The time-evolved density matrix can be expressed in terms of these eigenvalues as

[TABLE]

where $\rho_{\lambda}$ are the eigenstates of the Liouvillian, $\mathcal{L}[\rho_{\lambda}]=\lambda\,\rho_{\lambda}$. In the presence of damping, all the eigenvalues have negative real parts, with the exception of the steady state $\rho_{s}$ for which $\lambda=0$. This also represents the only physical state of the system that survives at long times. To preserve the trace at all times, including in steady state, $\rho_{s}$ always appears in this decomposition with a constant prefactor of $c_{s}=1$. Since the time evolution is trace-preserving this implies that $\text{Tr}(\rho_{\lambda\neq 0})=0$ and that $\rho_{\lambda\neq 0}$ are not valid density matrices with direct physical significance Peres (1995).

Here, we are interested in the time-dependent fraction of population in the steady state and therefore we re-express the expansion as a linear combination of orthogonal unit-trace density matrices. To accomplish this, we consider the eigenvalue, $\chi\neq 0$ of the Liouvillian, $\mathcal{L}[\cdot]$, defined as the eigenvalue with the smallest real part Risken et al. (1987). We observe that $\langle\rho_{s},\rho_{\chi}\rangle\neq 0$, with $\langle A,B\rangle=\text{Tr}(A^{\dagger}B)$ respresenting the matrix inner product. To ensure orthogonality, we use a Gram-Schmidt construction to define

[TABLE]

with $\langle\tilde{\rho}_{\chi},\rho_{s}\rangle=0$. The normalized density matrix $\tilde{\rho}=\tilde{\rho}_{\chi}/\text{Tr}(\tilde{\rho}_{\chi})$ has unit trace, but it is not ensured to be positive definite. Nevertheless, we can express the density matrix $\rho(t)$ as

[TABLE]

where $\sigma(t)$ only ensures the positive definiteness of the full density matrix, while the factors $\beta_{i}$ denote the decay rates.

In order to mimic the exponential relaxation observed in the numerical analysis of Sec. III, we surmise that $\beta_{0}(t)=(1-\tilde{\beta}_{0}e^{-\lambda_{e}t})$ with $\lambda_{e}$ represents the “escape rate”, thereby explicitly assuming Markovian switching dynamics. To estimate the rate at which the density matrix approaches the steady state in the long-time limit, we can calculate the time scale $\lambda_{e}$ by inserting Eq. (12) into the master equation, Eq. (2), and taking the inner product with $\rho_{s}$ to obtain,

[TABLE]

The quantum-induced switching rate can now be expressed as $\lambda_{e}=\dot{\beta}_{0}(t)/(1-\beta_{0}(t))$. Since $\sigma(t)$ only plays a role in preserving positive definiteness, the dominant time scale is not affected by it and we can omit $\sigma(t)$ in the analysis of the time scale without introducing unphysical behavior. Without $\sigma(t)$, the expression for $\lambda_{e}$ only accounts for the dynamics induced by the steady state and the smallest eigenvalue of $\mathcal{L}[\rho(t)]$, the two dominant contributions in the long time limit. While $\lambda_{e}$ is explicitly time-dependent, in the regime of interest it is approximately constant and we can therefore neglect the time dependence and simply evaluate for $t=0$ where $\rho(t=0)$ is the initial vacuum state. This leads to the following expression for $\lambda_{e}$:

[TABLE]

The solid lines in Fig. 6 (a) and (b) show the results of the analytical calculation of the switching times $\lambda_{e}^{-1}$ for $K/\kappa=2$ and $K/\kappa=0.2$, respectively. The dashed lines represent the time-scales extracted from the numerical master equation simulation also shown in Fig. 5. Though Eq. (14) predicts a larger value of $\tau=\lambda_{e}^{-1}$ than that observed numerically, the analytical calculations qualitatively match the numerical results. We can understand the larger value for $\tau$ obtained analytically from the fact that we neglected $\sigma(t)$. Indeed, in general $\sigma(t)$ has contributions from all eigenvalues of the Liouvillian which includes contributions with a larger negative real part than $\chi$. Therefore, the dynamics associated with $\sigma(t)$ must be strictly faster than $\tau$. In particular, the Liouvillian spectrum may show a two-fold degeneracy in the real part of the eigenvalues, which may speed up the dynamics by up to a factor of two. The prolonged time-scales observed in both experiment and in numerical simulation can, therefore, directly be qualified using a simple analytical quantum calculation.

IV.2 Semi-classical calculation: Quantum activation

The quantum calculation based on the Liouvillian qualitatively reproduces the dynamics of the full master equation. However, for weak nonlinearity, we expect a semi-classical treatment to be sufficient. Here, we compare the switching rates predicted by the numerical simulations with analytical results for a fluctuation-induced escape from a metapotential Dykman (2012). A metapotential is an effective potential that corresponds to the same equation of motion as the semi-classical model. In this treatment, the fluctuations associated with the noise $\zeta(t)$ in Eq. (9) are transformed into an effective temperature and the switching rate $\tau$ can then be evaluated as a thermal escape rate from a local minimum of the metapotential Dykman and Krivoglaz (1980); Vijay (2008). The increased time-scale in region II in this picture corresponds to a higher effective barrier in the metapotential.

To estimate the semi-classical escape rate, we begin by first considering the semi-classical equation, Eq. (9), with $\zeta(t)=0$ (See Appendix B), and later reinstate the effect of the noise. The general approach is to rephrase the complex equation into a real equation for a generalized position variable, $x$, that changes slowly in time. To this end, we solve Eq. (9) with $\zeta(t)=0$ and denote the low-amplitude solution as $\alpha_{0}$ and the unstable solution as $\alpha_{u}$. Next, we define a rotation angle $\varphi=\text{tan}^{-1}[-\text{Im}(\alpha_{u}-\alpha_{0})/\text{Re}(\alpha_{u}-\alpha_{0})]$ such that the quantity $x_{0}=e^{i\varphi}(\alpha_{u}-\alpha_{0})$ is a real number and the axis on the line from 0 to $x_{0}$ constitutes our position variable $x$. In order to determine the metapotential $U(x)$, we make a substitution $\alpha(t)=\alpha_{0}+e^{-i\varphi}z(t)$ and rewrite the equation of motion, Eq. (9), in the form $\dot{z}=F(z)$, where

[TABLE]

represents the effective force on the particle. The complex variable $z$ can be represented in terms of two real variables, playing the role of the coordinate and momentum, $z=x+ip$. Since the two states $\alpha_{0}$ and $\alpha_{u}$ are steady states, it follows that $F(0)=F(x_{0})=0$. Further, since the imaginary part accounts for the dynamics of the momentum $p$, we can write the derivative of the metapotential $U(x)$ as

[TABLE]

Integrating this equation results in a one-dimensional metapotential,

[TABLE]

with a minimum at 0 and maximum at $x_{0}$ as illustrated in the inset of Fig. 6 (c). Following Kramer’s escape law Kramers (1940); Hänggi et al. (1990), the escape rate over the barrier can be written as Dykman and Krivoglaz (1980); Vijay (2008)

[TABLE]

with $\Delta U=U(x_{0})-U(0)$ denoting the activation energy. Note that effective temperature is set by the mode linewidth, $\kappa$, since it enters as the prefactor for $\zeta(t)$ in Eq. (9). The rate $\gamma_{0}$ is the attempt frequency and it is extracted from the quadratic term of the potential, which can be expressed as $\frac{\gamma_{0}^{2}}{2\kappa}x^{2}$. When increasing the drive amplitude, $\epsilon$, the barrier height is decreased leading to a faster escape rate Dykman (2007).

In Figs. 6 (c) and (d), we show $\lambda_{e}$ obtained using this quantum-activation approach (solid lines) and compare it with the switching rates extracted from the semi-classical simulations (dashed lines, also shown in the inset of Fig. 5). We observe that, for small detunings, the simulations match the escape rate calculation of Eq. (18) quite well. However, a large discrepancy is observed for larger detunings when the switching time dramatically increases near the bifurcation point, where the rates between the two stable states equilibrates as already explained in the context of the numerical simulations. In the escape rate calculation, we are however only calculating a one-way rate and, thus, we do not capture the same increase in switching time. More interestingly, and as should be expected, we observe that while the escape rate calculation captures the behavior of the semi-classical simulations well, it does not capture the quantum corrections relevant in the strongly nonlinear regime ($K/\kappa=2$), as observed from its deviations from predictions of quantum calculation and master equation simulations [c.f. Figs 6 (a) and (c)].

To summarize this section, we have considered the time scales relevant for the relaxation of a Kerr-resonator towards the steady state. Going back to the regions introduced in Fig. 4, we noticed that, most significantly, the time scale for relaxation towards the steady state becomes very large in region II. This is consistent with experimental data presented in Sec. I, where we found a steady-state description to be inadequate due to the finite time-scale set by the measurement rate (see also the experiments presented in Muppalla et al. (2018)). We also found that while a semi-classical calculation adequately predicts the dynamics of a weakly nonlinear system, it fails to capture the structure of the time scales for relaxation, $\tau$, seen for large nonlinearities. The time scale for relaxation can be accurately predicted by a simplified but fully quantum model. On the other hand, an escape rate calculation matches the semi-classical trajectory simulations only partly, and completely fails to describe the dynamics for large nonlinearities.

IV.3 Crossover parameter

The observed breakdown of the semi-classical theory can be understood in a phenomenological manner by introducing a damping-dependent crossover temperature $T_{\gamma}$ Ingold (1995),

[TABLE]

where $\gamma_{0}$ is the attempt frequency in the metapotential introduced in the semi-classical calculation. Using this definition, we can define a crossover parameter $\xi$ as

[TABLE]

where ${T_{\kappa}=\hbar\kappa/k_{B}}$ denotes the effective temperature of the quantum fluctuations coupled to the oscillator. The semi-classical to quantum crossover boundary is set by $\xi=1$ Ingold (1995).

Figure 7 shows $\xi$ as a function of detuning of the oscillator. For $\xi<1$ or $T_{\kappa}>T_{\gamma}$, the decay dynamics can be described largely using the semi-classical activation treatment. However, as the detuning of the oscillator increases, and ${T_{\gamma}>T_{\kappa}}$, quantum tunneling effects can become essential to describe the nonlinear decay dynamics. We find that the attempt frequency $\gamma_{0}$, and consequently $\xi$, does not depend on the strength of the nonlinearity explicitly. This makes $\xi$ a universal quantity for a dissipative decay of a metastable state; nonetheless, strong nonlinearity $K/\kappa>1$ is essential for the oscillator to bifurcate at large enough detunings and enter the quantum regime as shown in Fig. 6. This intuition is consistent with significant deviations of both numerical simulations and Liouvillian-based analytical estimates from quantum activation results at large nonlinearity. Nonlinear oscillators with large $K/\kappa$ thus access a qualitatively different regime from that described by (thermal or quantum) activation-based models, in which rate of escape via tunneling is exponentially small Dykman (2007).

V Discussion and Outlook

In conclusion we have investigated the switching dynamics of a Kerr resonator in a combined experimental, numerical and analytical framework. We used experimental microwave spectroscopy data of a distributed mode of a Josephson junction array with $K/\kappa\approx 2$ and observed that the relaxation into the quantum steady state is not resolved in the experiment due to switching-times much longer than the probe time. We also simulate the experiment using a two-mode stochastic master equation and find that the numerical simulations match very well with the experimental data, thus confirming our interpretation. A recent experiment Muppalla et al. (2018) has directly measured the very slow switching dynamics in a weakly nonlinear Kerr resonator.

To further analyze the slow switching rates observed in the experiment we performed both semi-classical trajectory simulations and quantum master-equation simulations. We find that, for a range of parameters, the time scale to relax into the steady state is increased significantly beyond the natural decay time of the resonator, especially for strong nonlinearity ($K/\kappa\geqslant 1$) as compared to that for weak nonlinearity ($K/\kappa\ll 1$). We find that a semi-classical trajectory method is able to describe this slowdown only in the weak nonlinearity regime, while a full quantum master equation treatment is essential to calculate the switching rates for strongly nonlinear oscillators. To analytically estimate the time scales to reach the steady state, we compared a simplified quantum model and a semi-classical metapotential model. In the semi-classical metapotential treatment, the fluctuation-induced switching between bistable states is modelled as an activation over barrier in a metapotential, with the fluctuation-intensity determined by thermal or quantum noise. We find significant deviations from this model, especially when switching rates are small; moreover, these deviations are especially pronounced for large $K/\kappa$ and persist even far from bifurcation. In contrast, we obtain good qualitative agreement between quantum calculations and full master equation simulation in this regime. This is not entirely surprising since thermal or quantum activation models necessarily assume weak nonlinearity or weak driving conditions for the oscillator Dykman (2007). Our results indicate that switching dynamics in mesoscopic oscillators i.e. when $K/\kappa\gtrsim 1$ may be dominated by some other mechanism, such as dynamical quantum tunneling.

The conclusions drawn in this work regarding the time scales associated with switching dynamics in nonlinear Kerr-resonators are highly relevant for the characterization of state-of-the-art applications of Josephson devices especially in the strong nonlinearity regimes where quantum effects are more pronounced. In particular, we describe how the interplay of switching rates and the repetition rate of the experiment is essential to explain the experimentally measured nonlinear response of the superinductance presented in Fig. 3. We expect our results to guide the design of applications that similarly aim to observe non-linear spectroscopic signatures in Kerr resonators. In addition, our results can directly be used for assessing the performance of a bifurcation readout scheme for superconducting qubits Mallet et al. (2009); Vijay et al. (2009). In such a scheme, the required measurement time has to be a few multiples of the characteristic switching time. Thus, Eqs. (14) and (18) can be directly used for finding the required measurement time and, consequently, the expected qubit readout fidelity. These concerns become increasingly important as devices based on arrays of Josephson junctions, similar to the application in this work, continue to play an important role in experiments and theory proposals including various readout schemes Castellanos-Beltran et al. (2008); Schmitt et al. (2014); Eichler et al. (2014), quantum controllers Andersen et al. (2016); Sliwa et al. (2015); Kerckhoff et al. (2015) and even qubit architectures Bell et al. (2012); Vool et al. (2014); Dempster et al. (2014).

Furthermore, our calculations for the Kerr coefficients for the higher distributed modes of the Josephson array presented here (see Appendix A) indicate that these modes should be in the ‘mesoscopic’ regime (defined as $K\gtrsim\kappa$) investigated here. This regime is also optimal for direct observation of dissipative quantum tunneling Serban and Wilhelm (2007), which is usually obscured by the activation-dominated switching observed in the JBA regime Siddiqi et al. (2005). Moreover, since the switching rate does not follow activation dependence on noise intensity, this suggests that mesoscopic Kerr oscillators can be useful platforms to test dynamics resulting from non-Gaussian noise. More generally, our study provides a framework to explore multi-photon quantum effects Gevorgyan et al. (2012), quantum noise properties, and quantum-to-classical transitions in strongly nonlinear systems.

Acknowledgements

The authors acknowledge useful and insightful discussions with E. Doucet, S. Boutin and G. Kirchmair. CKA acknowledges fruitful discussions with K. Mølmer, the hospitality of Université de Sherbrooke and financial support from the Villum Foundation and the Danish Ministry of Higher Education and Science. The authors acknowledge support by the U.S. Army Research Office under the grant number W911NF-18-1-0212, by the U.S. Department of Energy under grant number DE-SC0019515, by NSERC, and in part by funding from the Canada First Research Excellence Fund.

Appendix A Kerr-resonators with arrays of Josephson junctions

We consider a Kerr resonator, consisting of a linear 1D array of $N$ Josephson junctions (JJ) forming a nonlinear inductance, with first and last junctions capacitively shunted to ground. We will, at first, assume that each of these junctions are sufficiently linear such that we can neglect the nonlinearity of each junctions. Each junction is described by its effective inductance $L_{J}$ and capacitance $C_{J}$. Furthermore, we include a parasitic capacitance to ground for each junction, $C_{0}$. This gives us the (linearized) Lagrangian for the array:

[TABLE]

The terms in Lagrangian corresponding to the end capacitances can be written as,

[TABLE]

where $C_{s}$ is the capacitance to the transmission line which controls the external quality factor of the array resonances, $C_{g}$ is the shunt capacitance on the first junction, and $C_{e}$ is the shunt capacitance on the last junction. Including the terms due to shunt capacitances, the full Lagrangian can then be witten as

[TABLE]

where we have introduced the symmetric matrices $\mathbb{C}$ and $\mathbb{L}$ along with the vector \phi=\big{\{}\phi_{1},\phi_{2},\ldots,\phi_{N+1}\big{\}}^{T}

Standing waves across the array now constitute a set of normal modes for the array. Formally, we find these standing mode by diagonalizing the matrix ${\Omega^{2}=\mathbb{C}^{-1}\mathbb{L}^{-1}}$, such that the Euler-Lagrange equations for these modes decouple. The eigenvalues of the matrix, $\Omega^{2}$ are the squares of the normal-mode frequencies of the Kerr resonator. We can similarly define the effective capacitances and inductances for these modes as

[TABLE]

where $\vec{v}_{k}$ are the corresponding eigenvectors of $\Omega^{2}$. The eigenfrequencies are now, by construction, given by $\omega_{k}=1/\sqrt{L_{k}C_{k}}$. Keep in mind that the actual modes are described by the physical phase variable $\phi(t)=\sum_{k}\phi_{k}(t)\vec{v}_{k}$, where $\phi_{k}(t)$ is a function oscillating with $\omega_{k}$; thus, the physical amplitude of the phase is in the $\phi_{k}(t)$ variables and no physical quantify depends on the normalization of $\vec{v}$.

A.1 Quantization of the modes

Having found the normal modes we can apply a canonical quantization scheme to develop a quantum model for the JJ array Devoret (1995). This is done by first rewriting the Lagrangian into diagonal form

[TABLE]

This Lagrangian now yields the conjugate variables

[TABLE]

such that we can introduce the quantum operators $\hat{\phi}_{k}$ and $\hat{q}_{k}$ satisfying the commutation relation $[\hat{\phi}_{k},\hat{q}_{k}]=i\hbar$. The Hamiltonian for the system is then readily obtained as

[TABLE]

which can be recast into the form

[TABLE]

by introducing the ladder operators $(a_{k}^{\dagger},a_{k})$ as

[TABLE]

The real physical phase variable is now described by the quantum variables, $\hat{\phi}=\sum_{k}\hat{\phi}_{k}\vec{\phi}_{k}$ and we confirm that for a given wave-function of the system the phase-variable is independent of the normalization of $\vec{\phi}_{k}$.

A.2 Reintroduction of the nonlinearity

The potential energy quadratic in phase variables, as used in Eq. (21) implicitly assumes small phase excursions $\phi\ll\varphi_{0}$ (with $\varphi_{0}=\hbar/(2e)$) for which Josephson junction potential $-E_{J}\cos(\phi/\varphi_{0})$ is approximated as $\cos x\approx 1-x^{2}$. To study the effect of nonlinear contributions of the junctions in the array, we include the next term in the expansion that gives rise to a nonlinear term in the potential, $\mathcal{U}_{nl}$, seen by the phase $\phi$

[TABLE]

Introducing the variable

[TABLE]

where ${v}_{k}[n]$ denotes the $n$th entry in the vector $\vec{v}_{k}$, we can rewrite the potential in terms of the quantum mode operators as,

[TABLE]

Generally, such a term will lead to Kerr terms as well as beam-splitter term. The beam-splitter terms are, however, only important between modes with small frequency difference, or if the system is strongly pumped by some external field of appropriate frequency that makes these terms resonant. Since we are primarily interested in the few photon regime of the lowest modes, which are all well-separated in frequency, we only consider the self-Kerr and the cross-Kerr terms. These can be expressed with the Hamiltonian

[TABLE]

where the Kerr-coefficients are found by inserting Eq. (29) into Eq. (32) and rearranging the terms such that we obtain Eq. (33) as

[TABLE]

with $\delta_{kl}$ being the Kronecker-delta. Figure 8 shows the Kerr shifts for the first eight array modes, calculated using the parameters for the array presented in Sec. I. We note that the presence of end capacitances loads the array and decreases the mode frequencies significantly Masluk et al. (2012). As evident from Eq. (34), the mode frequencies and concomitant Kerr shifts per photon can be much larger for an unloaded Josephson array.

As a last important detail, we should mention that the nonlinear coupling between the modes also drags the eigenfrequencies down such that the real mode frequencies become

[TABLE]

This extra contribution is very important as this can shift the frequency of each mode by as much as 2 GHz for parameters similar to the experiment of Sec. I. The Hamiltonian for the array is, therefore, expressed as

[TABLE]

Appendix B Classical solution for Kerr oscillator

We adopt a simple classical description of the driven-dissipative Kerr resonator by replacing the quantum operator $a$ with the classical complex variable $\alpha$ and ignore vacuum fluctuations. This yields the equation of motion (for convenience we change the phase of the drive)

[TABLE]

In steady state, this leads to the following expression for $n=|\alpha|^{2}$,

[TABLE]

which can be used to find the drive power needed for a given target photon number. We may take the inverse to obtain $n(\varepsilon)$, however, we are not ensured that this function will be single-valued. As a matter of fact, there will always be a set of parameters, $\Delta$ and $K$, for a given $\kappa$ where $n(\varepsilon)$ is multivalued and we can write the highest and lowest $n$’s in this bistable regime as

[TABLE]

This expression for $n_{c,\pm}$ can now be inserted into Eq. (38) to get the critical drive power to be in the bistable regime. An alternative approach would be to numerically propagate the dynamical equation, Eq. (37) until the steady state is reached. In contrast to the multi-valued solution of Eq. (38), this would only give a single value, but if the parameters are chosen to be in bistable regime the steady state will depend on the initial condition. Experimentally this appears as hysteresis in the system, where, say, a continuous adiabatic change of $\Delta$ will yield different results for increasing or decreasing $\Delta$.

It is worth noting that when $\kappa\gg K$, the dissipation in the resonator dominates the dynamics, destroying all the quantum coherent effects. In this respect, for a Kerr-resonator with $K/\kappa\ll 1$ a classical description of dynamics should suffice. To investigate this further, we compared the results of master equation and classical field equation simulations, shown in Fig. 9. The first thing that we observe in Fig. 9 is that for an initial coherent state, $\alpha(0)=\sqrt{30}$, which is very close to the steady state, we have initially a small amount of oscillations in both the quantum [Fig. 9 (a)] and classical [Fig. 9 (b)] simulations before both converge to the same steady state value for the photon number. Interestingly, the initial phase of the oscillations is very different between the classical and quantum simulations.

The second observation is that if we change the initial phase of the coherent state to $\alpha(0)=-\sqrt{30}$, the difference between the quantum and the classical simulation becomes pronounced. The quantum simulation shows a small drop in photon number followed by a convergence towards a steady state value different from the one with opposite phase. The classical simulation shows oscillations, since we are far away from the steady state, but eventually the photon number converges to the same steady-state value that we found for the simulation with opposite sign on $\alpha(0)$. Now one might be tempted to interpret this behavior as the mean photon number being very different in the classical case and in the quantum case. Nonetheless, on taking a closer look at the steady state photon number distribution obtained from the quantum simulation [Fig. 9(c)-(d)], we find that this is not the case. We see that the origin of the intermediate photon number comes from a dual-peaked photon number distribution; this can be interpreted as the steady state being a mixed state between the two peaks, each very close to a Poissonian distribution. The classical simulation is only single-valued, so it always selects only one of the solutions, depending on the size of the initial coherent state. We observe a similar behavior, when we start in the vacuum state. Here, there is always a very few number of photons, so the nonlinearity plays a smaller role. The quantum and the classical simulations are very similar, but, due to the quantum fluctuations, there is a small probability for the oscillator to end up in the high photon number state even when it starts in vacuum [see Fig. 9 (d)].

Finally, we see that if we change the sign of the detuning, $\Delta$, the quantum and classical simulation are identical. This is because the sign change moves us away from the bistable regime. So in both cases, the dynamics is deterministic and, since we are always at a low photon number, the nonlinearity again plays a small role.