Bistability in a Mesoscopic Josephson Junction Array Resonator

TL;DR

This paper experimentally studies stochastic switching in a mesoscopic Josephson junction array resonator, revealing how bistability and switching rates depend on pump strength, temperature, and readout, with implications for quantum measurement and control.

Contribution

It demonstrates the observation of stochastic bistable switching in a Josephson junction array resonator at GHz frequencies with low photon numbers, highlighting the role of nonlinearity and coupling.

Findings

Switching rates follow Kramer's law and vary from Hz to mHz.

Bistability appears at very low photon numbers due to anharmonicity.

Interplay of nonlinearity and coupling offers potential for quantum measurement applications.

Abstract

We present an experimental investigation of stochastic switching of a bistable Josephson junctions array resonator with a resonance frequency in the GHz range. As the device is in the regime where the anharmonicity is on the order of the linewidth, the bistability appears for a pump strength of only a few photons. We measure the dynamics of the bistability by continuously observing the jumps between the two metastable states, which occur with a rate ranging from a few Hz down to a few mHz. The switching rate strongly depends on the pump strength, readout strength and the temperature, following Kramer's law. The interplay between nonlinearity and coupling, in this little explored regime, could provide a new resource for nondemolition measurements, single photon switches or even elements for autonomous quantum error correction.

Figures

b of the main paper, where the maximum switching rate and the associated detuning are plotted as a function of $\bar{n}_{P}$ at the point of the maximum switch. In order to compare to the above model, by varying $\eta_{P}$ in the bistable region $\left[\eta_{+},\eta_{-}\right]$ (with $K_{PP}$ and $\kappa_{i}$ chosen as the values in the experiment) we first numerically locate the detuning $\tilde{\Delta}_{P}^{\mathrm{Max}}$ for a given value of $\bar{n}_{P}$ satisfying the symmetric potential condition namely $U(\bar{n}_{+})=U(\bar{n}_{-})=U(\bar{n}_{P})$. After adding a constant shift of $K_{

Tables

Table 1. Table 1: Parameters for the three array modes that can be directly measured with the VNA. f r subscript 𝑓 𝑟 f_{r} and Q tot subscript 𝑄 tot Q_{\text{tot}} were extracted from data. The Kerr and cross-Kerr coefficients are calculated from a fit to the dispersion relation Fig.  S3 d.

Mode	$f_r$	$Q_{\text{tot}}$	${\kappa }_i$	$K_i$	$K_{i7}$
	(GHz)		(kHz)	(kHz)	(kHz)
5	4.816(1)	26000	181	66	187
7	7.1058(2)	950	7500	133	-
9	9.278(1)	3375	2750	218	343

Equations

$$H/\hslash =\sum _{i=1}^N({\omega }_i{a}_i^{†}{a}_i+\frac{K_i}{2}{a}_i^{†}{a}_i{a}_i^{†}{a}_i)+\sum _{\begin{array}{c}i,j=1\\ i\ne j\end{array}}^N{K}_{i,j}{a}_i^{†}{a}_i{a}_j^{†}{a}_j.$$

$$S_{21}=1-\frac{\frac{1}{Q_{ext}}-2ı\frac{\delta f}{f_R}}{\frac{1}{Q_{tot}}+2ı\frac{f-f_R}{f_R}}.$$

$${\Delta }_i(P)=K_{i}P+\frac{K_i^{\prime }}{2}{P}^2+\frac{K_i^{\prime \prime }}{3}{P}^3$$

$$H/\hslash =({\tilde{\Delta }}_P+K_{PR}{\bar{n}}_R){\hat{a}}_P^{†}{\hat{a}}_P+\frac{K_{PP}}{2}{\left({\hat{a}}_P^{†}{\hat{a}}_P\right)}^2+i{\eta }_P({\hat{a}}_P^{†}-{\hat{a}}_P).$$

$$\frac{d{\alpha }_P}{dt}=\left[-i({\tilde{\Delta }}_P+K_{PP}/2+K_{PR}{\bar{n}}_R)-i{K}_{PP}|{\alpha }_P{|}^2-\frac{{\kappa }_P}{2}\right]{\alpha }_P+{\eta }_P$$

$${\bar{n}}_P=\frac{{\eta }_P^2}{{\left[\delta +K_{PP}{\bar{n}}_P\right]}^2+{\kappa }_P^2/4}$$

$$n_{c,\pm }=\frac{-2\delta }{3K_{PP}}\left[1\mp \sqrt{1-\frac{3}{4}\left(1+\frac{{\kappa }_P^2}{4{\delta }^2}\right)}\right].$$

$$U({\bar{n}}_P)=\frac{K_{PP}^2}{4{\eta }_P^2}{\bar{n}}_P^4-\frac{2K_{PP}\delta }{3{\eta }_P^2}{\bar{n}}_P^3+\frac{1}{2{\eta }_P^2}\left({\delta }^2+{\kappa }_P^2/4\right){\bar{n}}_P^2-{\bar{n}}_P.$$

$$P_L=\exp (-{\beta }_{eff}U({\bar{n}}_L))/\left[\exp (-{\beta }_{eff}U({\bar{n}}_L))+\exp (-{\beta }_{eff}U({\bar{n}}_R))\right]$$

Full text

The non-linearity provided by atoms and Josephson junctions is a necessary ingredient to observe quantum mechanical effects in cavity quantum-electro-dynamics (QED) and circuit QED (cQED) systems. Strong nonlinearites, much larger than the linewidth of the transition, are required to realize qubits Paik et al. (2011), implement quantum information protocols Kelly et al. (2015); Córcoles et al. (2015) and realize textbook quantum optics experiments Sun et al. (2014); Wang et al. (2016). Non-linearities much smaller than the linewidth of the transition are typically exploited for parametric processes Siddiqi et al. (2004); Bergeal et al. (2010); Roy and Devoret (2016) like amplification or frequency conversion at the quantum level.

Besides quantum information applications, there has been a growing interest to exploit cavity QED for ultralow-power classical logic elements Mabuchi (2009); Kerckhoff et al. (2010); Bose et al. (2012). This interest has been sparked by the ever growing all optical communication networks. Remarkably, a single photon transistor Chen et al. (2013), reminiscent of an electronic transistor, has been implemented for the optical domain. In this device a single photon can switch a large optical field. Realizing such devices has been a challenging endeavour as the required non-linearity is hard to realize, due to the weak interaction of optical light with atoms.

Much stronger light matter interactions can be achieved in the microwave regime using the cQED platform. In this context Josephson junction arrays (JJAs) have proven to be an ideal circuit element to build superconducting qubits with excellent coherence properties and unique tuning capabilities Manucharyan et al. (2009); Pop et al. (2014); Bell et al. (2016). Similarly, JJAs have also been used to build quantum limited parametric amplifiers Macklin et al. (2015); Castellanos-Beltran et al. (2008); Roy and Devoret (2016); Eichler and Wallraff (2014). Recently, the coherence properties of the self resonances of JJAs Hutter et al. (2011); Masluk et al. (2012), as well as their self-Kerr and cross-Kerr coefficients have been measured Weißl et al. (2015). The measured Kerr coefficient showed good agreement with a model based on a second order expansion of the Josephson potential Bourassa et al. (2012).

A regime of particular interest arises when the self-Kerr $K_{i}$ and cross-Kerr $K_{ij}$ nonlinear coefficients are on the order of the linewidth $\kappa$ of the system. In this regime the system will show a pronounced bistability Marthaler and Dykman (2006); Guo et al. (2011) at the single to few photon level. Bistability is a phenomenon which is relevant in many fields, ranging from chemistry García-Müller et al. (2008) and biology  McDonnell and Ward (2011); Eguíluz et al. (2000) to Josephson junction physics Silvestrini et al. (1988); Turlot et al. (1989) and cQED Mavrogordatos et al. . Very recently, an optically levitated nanoparticle has been shown to exhibit a stochastic bistability Ricci et al. (2017) and Kramers turnover Rondin et al. (2017).

In this Letter, we report on the realization of a JJA resonator with multiple modes and strong self-Kerr and cross-Kerr coefficients in the regime $K_{i},K_{ij}\approx\kappa$. We investigate the bistability of one mode of the JJA and characterize the dependence of the switching rate on the pump strength, readout strength and temperature.

The JJA consists of $10^{3}$ cascaded Josephson junctions, with a small capacitance to ground $C_{0}$, coupled to a 6 mm long microwave antenna and a shunt capacitance $C_{s}$, as shown in Fig. 1. The junctions are fabricated on a sapphire substrate using electron beam lithography and bridge-free double-angle evaporation Lecocq et al. (2011). An electron beam image of the junctions can be found in Fig. 1d. The junctions were designed to have a large ratio ${E_{J}}/{E_{C}}\approx 200$, in order to suppress coherent quantum phase slips (CQPS) Pop et al. (2010); Manucharyan et al. (2012). Here $E_{J}$ is the Josephson junction energy and $E_{C}$ is the charging energy. The parameters of the JJA were designed such that the fundamental resonance of the JJA combined with the shunt capacitance is around 1 GHz. The mode spacing for the first 10 modes is about 1.2 GHz and progressively becomes smaller for higher resonances Masluk et al. (2012); Weißl et al. (2015).

The JJA is placed inside a copper waveguide Kou et al. with a 6 GHz cutoff, as shown in Fig. 1a. Due to the capacitive coupling of the JJA to the waveguide, we can characterize the sample by performing microwave transmission measurements using a vector network analyzer (VNA). Due to the relative symmetry of the electric field of the waveguide and the antenna, the even modes of the JJA will couple poorly to the waveguide and not be visible in transmission measurements. The waveguide with the sample is mounted on the mixing chamber stage (10 mK) of a cryogen free dilution refrigerator. The sample is enclosed in a double layer cryoperm shield inside a completely closed copper can. The stainless steel input lines are attenuated with 20 dB at 4 K and 30 dB at base temperature. They are filtered with a combination of a 12 GHz low pass and an Eccosorb filter. The output stage consists of a 12 GHz low pass filter, two 4-12 GHz isolators and a 4-8 GHz high electron mobility transistor amplifier. The effective measurement bandwidth for direct transmission measurements using a VNA is limited to about 4-9 GHz due to the cutoff of the waveguide and the combined bandwidth of other microwave components.

Within the accessible measurement bandwidth we can characterize three resonances by fitting their transmission data to a notch type response function Probst et al. (2015). From these measurements we can extract the resonance frequencies ($\omega_{i}/2\pi$), the internal and the coupling quality factors (see supplementary material).

To observe the resonances of the modes outside the measurement bandwidth, we exploit the cross-Kerr interaction, which is induced by the junction non-linearity. The Hamiltonian for the JJA, up to second order, is given by

[TABLE]

The Hamiltonian consists of a self-Kerr term $K_{i}$ which leads to a photon number $n_{i}=a_{i}^{\dagger}a_{i}$ dependent frequency shift of mode $i$ and a cross-Kerr interaction $K_{ij}$, which leads to a frequency shift of mode $i$ depending on the photon number in all other modes $j$.

In Fig.2 we show the results of a two-tone spectroscopy where we sweep the frequency of a pump tone around mode five while weakly probing mode seven with the VNA. When the pump tone is resonant with mode five, we see a shift of the resonance frequency of mode seven. The measured frequency of mode five, using two-tone spectroscopy, matches the direct VNA measurement. We repeat this procedure, using mode seven as a readout, for frequencies between 900 MHz and 16 GHz. By comparing the measured frequencies of the entire JJA with the re-normalized mode frequencies calculated by diagonalizing the capacitance matrix (see supplementary material), we extract the JJA parameters $C_{0}$ = 0.152 fF, $C_{J}$ = 34 fF, $L_{J}$ = 1.25 nH, $C_{S}$ = 18 fF with a confidence range of about 20 %. These parameters match well to the expected design values and room temperature resistance measurements of the junctions.

With these parameters we can calculate $K_{i}$ and $K_{i,j}$ similar to Ref. Weißl et al. (2015). The results are summarized in the supplementary material for modes five, seven and nine. Using these coefficients we can calibrate the photon number by measuring the resonance frequency shift as a function of the applied power (see supplementary material). The attenuation extracted from these Kerr measurements matches within 4% to the attenuation in the cryostat, determined by independent transmission measurements.

Upon closer inspection of the two-tone scan in Fig.2a one can observe a bistable region (see Fig.2b) for a detuning of $\Delta_{P}=-2.58$ MHz from the bare resonance frequency of mode five. Around this frequency, two different cross-Kerr shifts of mode seven can be observed: A shift of $\Delta_{R}=-7.32$ MHz corresponding to a photon number of about 115 photons in mode five and a shift of 133 kHz corresponding to about one photon. In the two-tone scan we observed residence times exceeding ten seconds in either the high or low photon number states.

To precisely characterize the residence time, we implemented a readout scheme similar to the dispersive state detection of a superconducting qubit in a circuit QED architecture Blais et al. (2004). We monitor mode seven continuously on resonance with a readout power $\bar{n}_{R}$ of about half a photon such that we can not see an apparent shift or broadening of mode seven. When we pump mode five with $\bar{n}_{P}$ photons, we observe a change in the transmitted readout signal corresponding to jumps between the high and the low amplitude states. For each data point we average 500 measurements each lasting 19.9 $\mu$s to get a good signal to noise ratio. In addition, the analog to digital converter needs another 5 ms to transfer the data. Thus, it takes about 15 ms to acquire one data point.

In Fig.3a we show a typical time trace for a measurement time of 15 s and $\bar{n}_{P}$=9 photons in mode five. One can clearly observe two distinct amplitudes in transmission corresponding to two distinct photon numbers in mode five. We define the switching rate $\Gamma$ as the inverse of the average time between two transitions from the low to the high power state. The transient time between these two states is much faster than the data acquisition rate.

A histogram of the data shows a well separated bi-modal distribution for the two amplitude states. Furthermore, we extract the residence times in the high amplitude state, the resulting histogram shows an exponential behaviour typical for a Poissonian statisics (see supplementary material).When we sweep the detuning of the pump frequency $\Delta_{P}$ across the bistable point we see a change of the relative height of the peaks in the amplitude distribution $\rho(A)$, Fig.3b-e. There is a $\Delta_{P}^{\mathrm{Max}}$ for a given pump strength where the heights of the peaks are equal as shown for example in Fig.3b.

We repeat the switching rate measurements for a total of ten different pump powers ranging from 2.6 to 115 photons in mode five. In Fig.4b we plot the maximally achieved switching rate $\Gamma_{\mathrm{Max}}$ and the corresponding pump detuning $\Delta_{P}^{\mathrm{Max}}$ versus photon number. The lowest power we can observe switching for is $\bar{n}_{P}=2.6$ photons at a detuning of -169 kHz from the low power resonance. This detuning matches well to the prediction Ong et al. (2011) of $\Delta_{P}=\nicefrac{{\sqrt{3}}}{{2}}\,\kappa_{5}=$ 160 kHz with $\kappa_{5}=\leavevmode\nobreak\ 181$ kHz.

Typically, stochastic switching in a bistable system is described by Kramers law Hänggi et al. (1990). There, the switching rate is determined by the potential landscape and the fluctuations. For a symmetric potential it is given by $\Gamma_{\mathrm{Max}}=\Gamma_{0}\,\mathrm{exp}(-E_{b}/k_{b}T)$, where $E_{b}$ is the barrier height between the two stable solutions and the prefactor $\Gamma_{0}$ which depends on the relative strength of the dissipation and the potential Hänggi et al. (1990). Here, the activation of the switching between the two stable solutions likely originates from the dispersive shift of the resonator frequency due to photon number fluctuations in the readout mode and thermal fluctuations of the photon number in all modes.

In our case, the potential landscape is created by the interplay between the pump, the self-Kerr effect and the damping of the mode. An intuitive choice for this potential Marthaler and Dykman (2006); Guo et al. (2011) is provided by integrating the equation for the photon number in the steady state of a damped Kerr oscillator Agrawal and Carmichael (1979). From this model we extract the scaling of $E_{b}$ Luchinsky et al. (1999) and $\Gamma_{0}$ for the maximum switching rate as a function of the pump photon number $\bar{n}_{P}$. We use this scaling in a fit function (see Fig. 4b) with two free parameters to match our data (see supplementary material). We also add a $\Gamma_{\mathrm{res}}$ to take into account the finite switching rate for high pump powers. This finite rate could be limited by phase slips on the junctions of the JJA, which we estimated to be in the range of a few mHz.

Additionally, we observe a change of the switching rate with respect to $\Delta_{P}$ following a lorentzian curve. The point of maximum switching rate $\Delta_{P}^{\mathrm{Max}}$ also corresponds to a symmetric amplitude distribution. In Fig.4a one can see the change of the switching rate with respect to $\Delta_{P}$ for three different pump powers. The width of this lorentzian is about 72 kHz and the center shifts with increasing photon number in the pump. From our model we can also extract the shift of $\Delta_{P}^{\mathrm{Max}}$ with $\bar{n}_{P}$ which agrees well with the experiment as shown in Fig. 4b. The deviation at high photon numbers can be explained by higher order Kerr effects which are not taken into account in the model.

To better identify the origin of the switching we also varied the power in the readout tone. The measured results are plotted in Fig.4c for constant $\bar{n}_{P}=6$ photons. We can observe two effects for an increasing photon number in the readout mode: I) due to the cross-Kerr effect the bistable point moves to lower frequencies. II) in contrast to lowering the switching rate with the pump photon number, the readout photon number increases the switching rate. For $\bar{n}_{R}=2.5$ photons we see an increase in the switching rate by about a factor of three. This can be understood as a form of measurement induced dephasing. As the photon number in the readout resonator fluctuates, the position of the bistable point moves in frequency due to the cross-Kerr effect. This moves the pump in and out of the bistable region into the regime where either the low or the high power state are more likely, as the cross-Kerr frequency is on the order of the linewidth. This photon number fluctuation happens at a rate $\kappa_{7}=7.5$ MHz and is thus much faster than our acquisition time. For $\bar{n}_{P}>2.5$ photons, the switching becomes much faster and we cannot observe it any more due to the limited measurement bandwidth and signal to noise. Remarkably, changing $\bar{n}_{R}$ by only one photon for a constant $\Delta_{P}$ we can switch the state of mode five from the low to the high photon number occupation.

Lastly, to study the influence of thermal noise on the switching rate, we increased the cryostat base plate temperature from 10 mK to 50 mK ($\bar{n}_{P}=6$ photons, $\bar{n}_{R}=0.5$ photons). Increasing the temperature increases the average thermal population as well as the fluctuations of the photon number in all of the modes of the JJA, most notably for the lower frequency modes. As a consequence, depicted in Fig. 4d, we observe a shift in the magnitude and location of the maximum switching rate, similar to Fig. 4c. This is again due to the cross-Kerr interactions of mode five with all other modes. The observed shift can be explained, by an increase in the JJA temperature up to about 100-130 mK depending on the initial temperature. We can get an upper bound for the JJA temperature at the base temperature of the fridge, if we assume that the linewidth broadening we observe for mode five (see supplementary material) is due to the thermal population in the other modes of the JJA. This broadening can be explained by a minimal temperature of the JJA of $\approx 50$ mK. This is consistent with other cQED experiments Geerlings et al. (2013); Jin et al. (2015), where the devices are well above the fridge base temperature.

In conclusion, we have measured a stochastic bistability in a $10^{3}$ Josephson junction array which appears at a pump strength of only a few photons. This switching at low power is achieved by engineering the Kerr interaction strength to be comparable to the linewidth. We have seen an exponential decrease of the maximal switching rate for increasing pump strength as expected from Kramers law. For an increase in readout strength or temperature, the switching rate increases, likely due to photon induced dephasing through cross-Kerr interactions.

This proof of principle device demonstrates that it is possible for a few microwave photons in the readout mode to switch the photon occupation number by two orders of magnitude in the pump mode. As such it is a promising system to implement novel types of nondemolition measurements Kumar and DiVincenzo (2010), single photon microwave transistors Neumeier et al. (2013); Manzoni et al. (2014), single photon microwave switches Liao et al. (2009), Flip-Flop memories Andersen and Mølmer (2015) or even elements for autonomous quantum error correction Kerckhoff et al. (2010).

We gratefully acknowledge fruitful discussions with Christian Kraglund Andersen, Alexandre Blais, Michel Devoret, Linghzen Guo, Archana Kamal, Zlatko Minev, Michael Marthaler, and Peter Talkner. We also thank Christian Schneider for valuable comments on the manuscript. We thank Luigi Frunzio for help with the fabrication of a first generation junction array resonator. IMP and LG acknowledge funding provided by the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award endowed by the German Federal Ministry of Education and Research. Facilities use was supported by the KIT Nanostructure Service Laboratory (NSL). MP, OG, PV, GK is supported by the Austrian Federal Ministry of Science, Research and Economy (BMWFW). MLJ is supported by European Council (EUP0252).

I Experimental setup and transmission measurements

Fig. S1 shows the experimental apparatus used in the experiments. The waveguide with the sample is mounted on the mixing chamber stage (10 mK) of a cryogen free dilution refrigerator. The sample is enclosed in a double layer cryoperm shield inside a completely closed copper can. The stainless steel input lines are attenuated with 20 dB at 4 K and 30 dB at base temperature. They are filtered with a combination of a 12 GHz low pass and Eccosorb filters. The output stage consists of a 12 GHz low pass filter, two 4-12 GHz isolators and a 4-8 GHz high electron mobility transistor (HEMT) amplifier. We additionally use DC-blocks at room temperature.

Using a vector network analyzer (VNA), we measure the transmission $S_{21}$ of the waveguide. Since the sample is strongly capacitively coupled to the waveguide, we are able to directly detect some of the resonances on the network analyzer. Three resonant modes out of the 1000 resonances of the array are shown in Fig. S2a, c, d. We fit the measured $S_{21}$ parameter using a notch type response function to determine the quality factors and the resonance frequencies:

[TABLE]

The extracted parameters are summarized in Table 1.

Fig. S2b shows a measurement of the decay time $T_{1}$ of mode five for a pump strength of about $\bar{n}_{P}=$ 115 photons. To excite the resonator the pump was detuned by about -2.54 MHz from the bare frequency $\omega_{P}/2\pi=$ 4.8156 GHz of mode five. Similar to a two-tone measurement we use mode seven as a readout, with a readout strength of $\bar{n}_{R}=$ 0.5 photons. To perform the decay time measurement we excite mode five for a few $\mu$s with a pulse before performing a readout on the pump mode again using a pulse of a few $\mu s$. By varying the delay between the two pulses we measure how the excitation decays over time. From the measurements we find $T_{1}\approx 3\mu s$. Comparing this to the direct VNA measurements we see that we have an additional dephasing mechanism which broadens the resonances. As explained in the main text we believe this is due to a fluctuating cross Kerr shift induced by the thermal population in the other modes.

II Two-tone measurements and Dispersion relation

Since the other resonance frequencies of the JJA’s are out of the HEMT amplification bandwidth or below the cutoff frequency of the waveguide, we utilize a two-tone measurement technique to indirectly measure them. We use mode seven as read-out and apply a second pump tone. We sweep the frequency of this pump tone from  960 MHz to  20 GHz. Such a two-tone spectroscopy measurement can be seen as a pump probe experiment. The VNA is used to monitor one mode while a signal generator is used to excite a second one. We then observe the shift in the read-out mode due to populating other resonances through the pump tone. From this we can observe e.g. the fundamental mode of the array at  963 MHz and other higher resonant modes up to  20 GHz which is the limit of our signal generator. Two tone measurements for our array are shown in Fig.  S3.

In Fig. S3d, we show the dispersion relation of the array. The frequencies outside the HEMT bandwidth are extracted by two-tone spectroscopy measurements as shown in Fig.S3a, b, c. The mode spacing for the first 10 modes is about 1.2 GHz. For higher resonant modes the spacing between modes becomes smaller. The • are obtained by diagonalizing the capacitance matrix which includes the shunt capacitance ($C_{s}$), junction capacitance ($C_{J}$), the ground capacitance ($C_{0}$) and the Josephson inductance ($L_{J}$) for the entire array structure. Most of the even modes can not be observed from two tone measurements as the electric field distribution of these modes has a symmetry that does not couple to the waveguide. For higher mode numbers this symmetry is somewhat broken due to inhomogeneities of the junctions and we can again excite these modes.

We also note a discrepancy between the model and the measurements for modes 9 and 10. This is most likely due to the capacitive cross coupling between the parallel segments of the chain of Josephson junctions (see Fig. 1 in main manuscript). This effect can be accounted for by introducing additional capacitances with a minimal impact on the values of the self-Kerr and cross-Kerr coefficients for the lower modes.

III Kerr coefficients and Photon number calibration.

The Kerr-nonlinearity manifests itself as a frequency shift that depends on the circulating power in each of the resonant modes Tweissl2015 ; Jbourassa2012 . By varying the input power on a low signal level, we are able to measure this frequency shift. The Self-Kerr and cross-Kerr shift can be used to calibrate the photon number in a resonant mode by measuring the frequency response as a function of drive power.

III.1 Self-Kerr Measurements.

For the self-Kerr measurements only one mode $i$ of the chain is excited. We record the frequency shift $\Delta_{P}$ with respect to the input power using direct transmission measurements. We fit each data-set with a notch type response function to extract the frequency $f_{r}$ and $Q_{\text{tot}}$. From the shift relative to the bare resonance frequency we can extract the self-Kerr coefficient $K_{i}$.

Fig. S4a shows the self-Kerr measurements on mode five. For low input power $P$ we observe that $\Delta_{5}$ changes linearly with power. For high input power, higher order Kerr terms ($K_{i}^{\prime},K_{i}^{\prime\prime}$) start to play a role. To take this into account we fit the following dependence

[TABLE]

to the data. With this we can extract the Kerr coefficients in units of Hz/W.

Fig. S4b shows the self-Kerr measurements on mode seven. Here we observe that $\Delta_{7}$ changes linearly with the power due to the higher linewidth ($\kappa_{7}$) of this mode.

III.2 Cross-Kerr Measurements $K_{75}$ utilizing two-tone spectroscopy

For determining the cross-Kerr coefficients $K_{ij}$ we utilize two-tone spectroscopy measurements. One mode $i$ of the array is used as the readout with a power of about one photon and the other mode $j$ is excited with varying power. From the shift $\Delta_{R}$ of the resonance frequency of the readout mode $\omega_{R}/2\pi$ upon application of a pump tone, we can extract the cross-Kerr shift $K_{ij}$.

Fig. S5a shows a typical two-tone spectroscopy measurements to extract the cross-Kerr coefficient $K_{75}$ by pumping mode five and observing the frequency shift on mode seven. We drive mode seven with a constant readout strength of $\bar{n}_{R}=$ 0.5 photons. To extract the maximal frequency shift $\Delta_{R}$, we fit each pump frequency to a notch type response function and extract the resonance frequency(solid black line). The maximal shift of the readout resonator for a given pump power results in one datapoint in Fig. S4c. We then repeat this procedure for different pump powers e.g. Fig. S5b-c. From the measurement in Fig. S5 one can also clearly observe that the mode becomes bistable as the power is increased.

Fig. S4c shows the result of all two tone measurement to determine the cross-Kerr shift $K_{75}$ when driving mode five and using mode seven as the readout. For low input power we observe that $\Delta_{7}$ changes linearly with the power but then rapidly higher order terms come into play. In this case, even a third order polynomial fit does not agree with the measured $K_{75}$. We choose instead to characterize cross-Kerr coefficient by the slope at low powers in Fig. S4c. By using the dispersion relation fit in Fig. S3d and the resulting diagonalized capicitance matrix, we extract the ratio between $K_{5}$ and $K_{75}$. This ratio, together with $K_{5}$ determined from the linear part of the fit function Fig. S4a allows us to compute $K_{75}$. A linear fit with the computed slope $K_{75}$, shown by solid red line in Fig.S4c, shows good agreement with the measurements at low powers.

Fig. S4d shows the cross-Kerr measurements $K_{57}$ on mode five using mode seven as the readout. Here the third order polynomial fit agrees well with the data. When we take the ratio of the first order self and cross-Kerr coefficients extracted from fitting the data in Fig. S4b and Fig. S4d and compare it to the Kerr coefficients extracted from fitting the dispersion relation we find an agreement within 99.8%.

Furthermore, if we convert the Kerr coefficients extracted from Fig. S4a-d from Hz/W to Hz/photon we find that the required conversion factor matches within 4% to the attenuation in the cryostat, determined by independent transmission measurements.

For 0.1 nW input power we get about 10 photons in mode five. We can apply this same conversion factor to the whole polynomial Eq. S3 and then use the observed frequency shift for a given input power to calculate the circulating photon number.

III.3 Power dependence of the line-width

Fig. S6 shows the linewidth of mode five and seven with respect to the input power. The line-widths are extracted from direct transmission measurements. Each measurement is fit with a notch type response function to extract the resonance frequency $f_{r}$ and $Q_{\text{tot}}$ for a given drive power. The increase of the line-width $\kappa_{i}$ as a function of the circulating power is expected from the self-Kerr effect. For a circulating power of less than one photon in mode five we find $\kappa_{5}(\bar{n}_{P}\rightarrow 0)=$ 181 kHz, and for mode seven $\kappa_{7}(\bar{n}_{P}\rightarrow 0)=$ 7.5 MHz.

IV Width of bistable region for varying photon number

Fig. S7 shows the width of the bistable region obtained from the switching rate measurements (see main paper) for the ten different pump strengths shown in Fig. 4b in the main text. For very low pump strengths up to 10 photons the width of the bistable region is approximately constant $\approx$ 72 kHz. For photon numbers greater than 10 photons the bistable region becomes wider.

V Residence time and state population inversion

Fig. S8a, b, c show the exponential dependence of the residence time for a pump strength of \bar{n}_{P}=9\leavevmode\nobreak\ photons, \bar{n}_{P}=6\leavevmode\nobreak\ photons and \bar{n}_{P}=4\leavevmode\nobreak\ photons at a detuning of $\Delta_{P}=\Delta_{P}^{Max}$ MHz from mode five. This indicates that the transitions are random and follow a Poissonian statistics. From the exponential fit to the data we extract a mean residence time in the high state $\left<T_{UP}\right>$.

Fig.S8d, e, f show the probability for being either in the low or high photon state as we scan $\Delta_{P}$ across the bistable region for three different pump strengths \bar{n}_{P}=9\leavevmode\nobreak\ photons, \bar{n}_{P}=6\leavevmode\nobreak\ photons and \bar{n}_{P}=6\leavevmode\nobreak\ photons. It shows the state population inversion between the low photon and the high photon state, following a sigmoid behaviour $f(x)=1/(1+\textnormal{e}^{-x})$.

We note that the statistics of our results for different $\bar{n}_{P}$, for both the residence time and the probability distributions, are consistent with bistable systems well described by the Kramers model SHaenggi90 .

VI Theoretical Model for the Switching Rates

In this section we provide details regarding the theoretical model, based on Kramers theory of switching SHaenggi90 , used to obtain the fits in Fig. 4b of the main paper. To this end, we restrict to a simplified version of the full hamiltonian (1) of the main paper focusing only on the pump mode $5$ (indexed by $P$ in what follows) and include in addition a constant shift due to the cross-Kerr interaction $K_{PR}\bar{n}_{R}$ with the readout mode $7$ (indexed by $R$):

[TABLE]

Here, we have written the hamiltonian in a rotating frame with respect to the driving of strength $\eta_{P}$ and with frequency $\omega_{d}$ and the detuning parameter $\tilde{\Delta}_{P}=\omega_{P}-\omega_{d}$ with respect to the bare frequency of the pump mode $\omega_{P}$. Taking the classical limit of the Heisenberg-Langevin equation for the mode $\hat{a}_{P}$, we can obtain the following equation for the complex amplitude $\alpha_{P}=\langle\hat{a}_{P}\rangle$ mode (ignoring noise terms):

[TABLE]

The steady state solution for this equation, defining $\delta=\tilde{\Delta}_{P}+K_{PP}/2+K_{PR}\bar{n}_{R}$ written in terms of the average photon number $\bar{n}_{P}=|\alpha_{P}|^{2}$ is:

[TABLE]

This equation can have either one or three real roots. When it has three real roots, we have bistability. For a given value of $\delta$ and $K_{PP}$, it has three real solutions when $|\delta|>\sqrt{3}\kappa_{P}/2$ SAgarwal79 ; SOng11 , and the driving strength falls in the range $\eta_{-}\leq\eta_{P}\leq\eta_{+}$, where $\eta_{\pm}^{2}=n_{c,\pm}\left(\left[\delta+K_{PP}n_{c,\pm}\right]^{2}+\kappa_{P}^{2}/4\right)$, with the extremal photon numbers given by

[TABLE]

As we commented in the main paper, one key ingredient to analyze switching rates in a bistable system within the Kramers framework is a potential landscape in which the two stable solutions occur as local minima. A simple choice for such a potential is obtained by integrating Eq. S6 SAgarwal79 with respect to $\bar{n}_{P}$ (and dividing by $\eta_{P}^{2}$ to make it dimensionless) giving

[TABLE]

Note that this is not a real potential in the Hamiltonian sense, but a fictitious one for the average photon number $\bar{n}_{P}$ treated as an independent degree of freedom. The critical points of $U(\bar{n}_{P})$, namely points where $dU/d\bar{n}_{P}\equiv 0$ precisely satisfy Eq. S6 and the solutions to the equations identify the extrema of the potential landscape. From the quartic form of the potential, as shown in Fig. S9, we can anticipate that in the bistable region the potential has a double well shape with the two local minima at $\bar{n}_{P}=\bar{n}_{L},\bar{n}_{H}$ and the local maxima (top of the barrier between the two wells) at $\bar{n}_{P}=\bar{n}_{0}$ (with $\bar{n}_{L}<\bar{n}_{0}<\bar{n}_{H}$). Let us denote the oscillation frequencies at the bottom of the wells (top of the barrier) as $\omega_{L},\omega_{H}$ ($\omega_{0}$). The barrier height between $\bar{n}_{L}$ ($\bar{n}_{H}$) and $\bar{n}_{H}$ ($\bar{n}_{L}$) is given by $E_{b,L}=U(\bar{n}_{0})-U(\bar{n}_{L})$ ($E_{b,H}=U(\bar{n}_{0})-U(\bar{n}_{H})$). Given these parameters, the Kramers formula SHaenggi90 predicts that the rate of transition out of the wells are of the form $\Gamma_{L\rightarrow H}=\Gamma_{0,L\rightarrow H}\exp(-\beta_{eff}E_{b,L})$ ($\Gamma_{H\rightarrow L}=\Gamma_{0,H\rightarrow L}\exp(-\beta_{eff}E_{b,H})$) with the functional form of the pre-factor decided by the relative strengths of the damping rate and the well frequencies $\omega_{L,H},\omega_{0}$. In our treatment note that $\beta_{eff}$ is an effective dimensionless temperature. In addition in the limit of a two state model with localized states at $\bar{n}_{L}$ and $\bar{n}_{H}$, the average state population will be given by

[TABLE]

and $P_{H}=1-P_{L}$. The total switching rate, $\Gamma_{H\rightarrow L}+\Gamma_{L\rightarrow H}$, is maximised when the energy barriers are the smallest and we find that this happens for a symmetric configuration with the same barrier height for both directions which we denote as $E_{b}$.

In order to compare the experimental results with the theoretical approach above, we first note that the parameter $\eta_{P}$ is not directly known but has to be inferred from the circulating photon number $\bar{n}_{P}$. Let us consider the situation in Fig. 4b of the main paper, where the maximum switching rate and the associated detuning are plotted as a function of $\bar{n}_{P}$ at the point of the maximum switch. In order to compare to the above model, by varying $\eta_{P}$ in the bistable region $\left[\eta_{+},\eta_{-}\right]$ (with $K_{PP}$ and $\kappa_{i}$ chosen as the values in the experiment) we first numerically locate the detuning $\tilde{\Delta}_{P}^{\mathrm{Max}}$ for a given value of $\bar{n}_{P}$ satisfying the symmetric potential condition namely $U(\bar{n}_{+})=U(\bar{n}_{-})=U(\bar{n}_{P})$. After adding a constant shift of $K_{PP}/2+K_{PR}\bar{n}_{R}$ to our numerically determined $\tilde{\Delta}_{P}^{\mathrm{Max}}$, as shown by the solid red line in Fig. 4b of the main paper, we find good agreement with the experimental result. From Eq.S8, we can extract the barrier height $E_{b}$, and frequencies $\omega_{L}$ and $\omega_{0}$ for this symmetric point as a function of $\bar{n}_{P}$. In order to fit the switching rate to Kramers equation, we need to choose a form for the pre-factor $\Gamma_{0}$. We found that in general the form $\Gamma_{0}\propto\omega_{L}\omega_{0}$ valid for the over-damped regime fits best to the experimental data. Fitting the measured rates $\Gamma_{max}$ to the functional form $A\omega_{L}\omega_{0}\exp(-\beta_{fit}E_{b})+\Gamma_{res}$ on a log-log scale using least squares procedure, we find the fit parameters $\Gamma_{res}=(0.006\pm 0.002)$ Hz, $A=(54.0\pm 15.0)$ Hz and the effective temperature $\beta_{fit}=(2.4\pm 0.2)$. The errors we have quoted here are the standard deviation on the estimated best fit parameters. The fitted curve was depicted by the black line in Fig. 4b in the main paper. We reiterate that in the fitting procedure the parameters $E_{b}$, $\omega_{L}$ and $\omega_{0}$ were calculated numerically from the potential Eq. S8. In addition for the cases of $\bar{n}_{P}=4,6,9$, we also fitted the form Eq. S9 (with $U(\bar{n}_{L,R})$ calculated numerically) for the average state population as a function of $\Delta_{P}$ and found the effective temperatures $\beta_{fit}=\{3.2\pm 0.3,1.4\pm 0.3,1.3\pm 0.1\}$ respectively from the fits depicted in Fig. S8 d, e, f. We can see that this temperature range is similar to the value we obtained from the fitting to Fig.4b. From Fig. S8 d, e, f we can see that while the numerical model helps to locate $\Delta_{P}^{\mathrm{Max}}$ accurately, there are some discrepancies in terms of the width of the bistable region, especially regarding the $\Delta_{P}$ for the onset of bistability once we fix $\eta_{P}$.