Observation and stabilization of photonic Fock states in a hot radio-frequency resonator

TL;DR

This paper demonstrates the direct observation, control, and stabilization of photonic Fock states in a hot radio-frequency resonator using a superconducting qubit, advancing quantum thermodynamics and quantum interface technologies.

Contribution

It introduces a method to observe and stabilize quantum states in a high-temperature RF resonator, extending circuit QED to new regimes.

Findings

Successful stabilization of photonic Fock states.

Real-time observation of re-thermalization dynamics.

Quantum control of thermal photons at megahertz frequencies.

Abstract

Detecting weak radio-frequency electromagnetic fields plays a crucial role in wide range of fields, from radio astronomy to nuclear magnetic resonance imaging. In quantum mechanics, the ultimate limit of a weak field is a single-photon. Detecting and manipulating single-photons at megahertz frequencies presents a challenge as, even at cryogenic temperatures, thermal fluctuations are significant. Here, we use a gigahertz superconducting qubit to directly observe the quantization of a megahertz radio-frequency electromagnetic field. Using the qubit, we achieve quantum control over thermal photons, cooling to the ground-state and stabilizing photonic Fock states. Releasing the resonator from our control, we directly observe its re-thermalization dynamics with the bath with nanosecond resolution. Extending circuit quantum electrodynamics to a new regime, we enable the exploration of thermodynamics at the quantum scale and allow interfacing quantum circuits with megahertz systems such as spin systems or macroscopic mechanical oscillators.

Tables

Table 1. Table S1: Four-wave mixing terms Only half of terms are shown, the other half can be obtained by taking the hermitian conjugate of all these terms. Terms become approximately time-independent around the frequency ω p subscript 𝜔 p \omega_{\text{p}} given in the left column.

${\omega }_{\text{p}}\simeq$	prefactor	interaction
Stark shift
	$-2A_{\text{H}}{\left|{\xi }_{\text{p}}\right|}^2$	${\hat{a}}^{†}\hat{a}$
	$-\chi {\left|{\xi }_{\text{p}}\right|}^2$	${\hat{b}}^{†}\hat{b}$
Heating interactions
$\left({\omega }_{\text{H}}+{\omega }_{\text{L}}\right)/2$	$-A_{\text{H}}^{\frac{3}{4}}{A}_{\text{L}}^{\frac{1}{4}}{\left({\tilde{\xi }}_{\text{p}}^{\ast }\right)}^2$	$\hat{a}\hat{b}$
${\omega }_{\text{H}}+2{\omega }_{\text{L}}$	$-\chi {\tilde{\xi }}_{\text{p}}^{\ast }/2$	$\hat{a}{\hat{b}}^2$
$2{\omega }_{\text{H}}+{\omega }_{\text{L}}$	$-A_{\text{H}}^{\frac{3}{4}}{A}_{\text{L}}^{\frac{1}{4}}{\tilde{\xi }}_{\text{p}}^{\ast }$	${\hat{a}}^2\hat{b}$
Cooling interactions
$\left({\omega }_{\text{H}}-{\omega }_{\text{L}}\right)/2$	$-A_{\text{H}}^{\frac{3}{4}}{A}_{\text{L}}^{\frac{1}{4}}{\left({\tilde{\xi }}_{\text{p}}^{\ast }\right)}^2$	$\hat{a}{\hat{b}}^{†}$
${\omega }_{\text{H}}-2{\omega }_{\text{L}}$	$-\chi {\tilde{\xi }}_{\text{p}}^{\ast }/2$	$\hat{a}{\left({\hat{b}}^{†}\right)}^2$
$2{\omega }_{\text{H}}-{\omega }_{\text{L}}$	$-A_{\text{H}}^{\frac{3}{4}}{A}_{\text{L}}^{\frac{1}{4}}{\tilde{\xi }}_{\text{p}}^{\ast }$	${\hat{a}}^2{\hat{b}}^{†}$
Unused interactions
$3{\omega }_{\text{H}}$	$-A_{\text{H}}{\tilde{\xi }}_{\text{p}}/3$	${\hat{a}}^3$
${\omega }_{\text{H}}/3$	$-A_{\text{H}}{\left({\tilde{\xi }}_{\text{p}}^{\ast }\right)}^3/3$	$\hat{a}$
${\omega }_{\text{H}}$	$-A_{\text{H}}{\left({\tilde{\xi }}_{\text{p}}^{\ast }\right)}^2/2$	${\hat{a}}^2$
${\omega }_{\text{H}}$	$-\chi {\tilde{\xi }}_{\text{p}}$	$\hat{a}{\hat{b}}^{†}\hat{b}$
${\omega }_{\text{H}}$	$-A_{\text{H}}{\tilde{\xi }}_{\text{p}}$	${\hat{a}}^{†}{\hat{a}}^2$
${\omega }_{\text{H}}$	$-A_{\text{H}}{\left({\tilde{\xi }}_{\text{p}}^{\ast }\right)}^3$	$\hat{a}$
${\omega }_{\text{H}}$	$-A_{\text{H}}{\tilde{\xi }}_{\text{p}}$	$\hat{a}$
${\omega }_{\text{H}}$	$-\chi {\tilde{\xi }}_{\text{p}}/2$	$\hat{a}$
$3{\omega }_{\text{L}}$	$-A_{\text{H}}^{\frac{1}{4}}{A}_{\text{L}}^{\frac{3}{4}}{\tilde{\xi }}_{\text{p}}/3$	${\hat{b}}^3$
${\omega }_{\text{L}}/3$	$-A_{\text{H}}^{\frac{3}{4}}{A}_{\text{L}}^{\frac{1}{4}}{\left({\tilde{\xi }}_{\text{p}}^{\ast }\right)}^3/3$	$\hat{b}$
${\omega }_{\text{L}}$	$-\chi {\tilde{\xi }}_{\text{p}}^2/4$	${\hat{b}}^2$
${\omega }_{\text{L}}$	$-2A_{\text{H}}^{\frac{3}{4}}{A}_{\text{L}}^{\frac{1}{4}}{\tilde{\xi }}_{\text{p}}$	${\hat{a}}^{†}\hat{a}\hat{b}$
${\omega }_{\text{L}}$	$-A_{\text{H}}^{\frac{1}{4}}{A}_{\text{L}}^{\frac{3}{4}}{\tilde{\xi }}_{\text{p}}$	${\hat{b}}^{†}{\hat{b}}^2$
${\omega }_{\text{L}}$	$-A_{\text{H}}^{\frac{3}{4}}{A}_{\text{L}}^{\frac{1}{4}}{\left({\tilde{\xi }}_{\text{p}}^{\ast }\right)}^3$	$\hat{b}$
${\omega }_{\text{L}}$	$-A_{\text{H}}^{\frac{1}{4}}{A}_{\text{L}}^{\frac{3}{4}}{\tilde{\xi }}_{\text{p}}$	$\hat{b}$
${\omega }_{\text{L}}$	$-A_{\text{H}}^{\frac{3}{4}}{A}_{\text{L}}^{\frac{1}{4}}{\tilde{\xi }}_{\text{p}}$	$\hat{b}$

Table 2. Table S2: Fitted system parameters

Quantity	Symbol	Value	Equation
Hamiltonian parameters
Dressed high-mode frequency ($|g,0⟩\to |e,0⟩$)	${\omega }_H$	$2\pi \times$ 5.911 GHz	${\overline{\omega }}_H+A_H/\hslash +\chi /2\hslash$
Dressed low-mode frequency ($|g,0⟩\to |g,1⟩$)	${\omega }_L$	$2\pi \times$ 173 MHz	${\overline{\omega }}_L+A_L/\hslash +\chi /2\hslash$
Bare high-mode frequency	${\overline{\omega }}_H$	$2\pi \times$ 6.113 GHz	$\sqrt{\frac{L+L_J}{L{L}_J{C}_L}}$
Bare low-mode frequency	${\overline{\omega }}_L$	$2\pi \times$ 182 MHz	$\frac{1}{\sqrt{\left(L+L_J\right)C_L}}$
High-mode anharmonicity	$A_H$	$h\times$ 192 MHz	$\frac{e^2}{2C_H}\left(\frac{L}{L+L_J}\right)$,
Low-mode anharmonicity	$A_L$	$h\times$ 495 kHz	$\frac{e^2}{2C_L}{\left(\frac{L_J}{L+L_J}\right)}^3$,
Cross-Kerr	$\chi$	$h\times$ 21.29 MHz	$2\sqrt{A_L{A}_H}$
Dissipation rates
High-mode dissipation rate	$\kappa$	$2\pi \times$ 3.70 MHz
External coupling rate	${\kappa }_{\text{ext}}$	$2\pi \times$ 1.63 MHz
Low-mode dissipation rate	$\gamma$	$2\pi \times$ 23.50 kHz
Low-mode external coupling rate	${\gamma }_{\text{ext}}$	$2\pi \times$ 1.99 Hz
High-mode quality factor	$Q_H$	1599	${\omega }_H/\kappa$
High-mode external quality factor	$Q_H^{\left(\text{ext}\right)}$	3617	${\omega }_H/{\kappa }_{\text{ext}}$
Low-mode quality factor	$Q_L$	7348	${\omega }_L/\gamma$
Low-mode external quality factor	$Q_L^{\left(\text{ext}\right)}$	87 $\times {10}^6$	$Z_0\sqrt{\frac{C_L}{L+L_J}}{\left(\frac{C_c}{C_L}\right)}^2$
Thermal parameters
High-mode temperature	$T_H$	112 mK
Low-mode temperature	$T_L$	17 mK
High-mode occupation number	$n_{\text{th}}^{\left(H\right)}$	0.09	$\frac{1}{e^{\frac{\hslash {\omega }_H}{k_B{T}_H}}-1}$
Low-mode occupation number	$n_{\text{th}}$	1.62	$\frac{1}{e^{\frac{\hslash {\omega }_L}{k_B{T}_L}}-1}$
Circuit parameters
Josephson energy	$E_J$	$h\times$ 4.01 GHz	$\frac{{\hslash }^2{\overline{\omega }}_H^2{\overline{\omega }}_L^2}{8{\left({\overline{\omega }}_H\sqrt{A_L}+{\overline{\omega }}_L\sqrt{A_H}\right)}^2}$
Josephson inductance	$L_J$	41 nH	$\frac{{\hslash }^2}{4e^2{E}_J}$
Low-mode capacitance	$C_L$	11.1 pF	$\frac{e^2\sqrt{A_L}{\overline{\omega }}_H^3}{2{\left({\overline{\omega }}_L\sqrt{A_H}+{\overline{\omega }}_H\sqrt{A_L}\right)}^3}$
High-mode capacitance	$C_H$	40.7 fF	$\frac{e^2{\overline{\omega }}_L}{2\left(\sqrt{A_H{A}_L}{\overline{\omega }}_H+A_H{\overline{\omega }}_L\right)}$
High-mode inductance	$L$	28.2 nH	$\frac{2\sqrt{A_H}{\left({\overline{\omega }}_L\sqrt{A_H}+{\overline{\omega }}_H\sqrt{A_L}\right)}^2}{e^2\sqrt{A_L}{\overline{\omega }}_H^3{\overline{\omega }}_L}$
Coupling capacitor	$C_c$	0.95 fF	$C_H\sqrt{\frac{{\kappa }_{\text{ext}}L{L}_J}{Z_0(L+L_J))}}$
Feedline impedance	$Z_0$	50$\Omega$

Equations

$$A_{\text{L}}=-\frac{e^2}{2C_{\text{L}}}{\left(\frac{L_{\text{J}}}{L+L_{\text{J}}}\right)}^3,A_{\text{H}}=-\frac{e^2}{2C_{\text{H}}}\left(\frac{L}{L+L_{\text{J}}}\right).$$

$$\begin{array}{rl}\hat{H}& =\hslash {\omega }_{\text{H}}{\hat{a}}^{†}\hat{a}+\hslash {\omega }_{\text{L}}{\hat{b}}^{†}\hat{b}\\ & -\frac{A_{\text{H}}}{2}{\hat{a}}^{†}{\hat{a}}^{†}\hat{a}\hat{a}-\frac{A_{\text{L}}}{2}{\hat{b}}^{†}{\hat{b}}^{†}\hat{b}\hat{b}\\ & -\chi {\hat{a}}^{†}\hat{a}{\hat{b}}^{†}\hat{b},\end{array}$$

$$P_n{\kappa }_{\text{ext}}/{\kappa }_n,$$

$$\begin{array}{rl}{\hat{H}}_{\text{int}}=& -\hslash g\sqrt{n+1}|f,n⟩⟨g,n+1|+h.c.,\end{array}$$

$$\frac{P_1}{P_0}\simeq \frac{n_{\text{th}}}{1+n_{\text{th}}+C}.$$

$$\begin{array}{rl}{\dot{P}}_n& =-n\gamma (n_{\text{th}}+1)P_n+n\gamma n_{\text{th}}{P}_{n-1}\\ & -(n+1)P_n\gamma n_{\text{th}}+(n+1)P_{n+1}\gamma (n_{\text{th}}+1).\end{array}$$

$$A_m=-\frac{2e^2}{L_{\text{J}}{\bar{\omega }}_m^2(\text{Im}{Y}^{\prime }({\bar{\omega }}_m){)}^2}.$$

$$\begin{array}{rl}\hat{H}& =\sum _m\hslash {\bar{\omega }}_m{\hat{a}}_m^{†}{\hat{a}}_m+\underset{\text{junction non-linearity}}{\underbrace{E_{\text{J}}[1-\cos \hat{\phi }]-E_{\text{J}}\frac{{\hat{\phi }}^2}{2}}},\\ \text{where }\hat{\phi }& =\sum _m{\left(2A_m/E_{\text{J}}\right)}^{1/4}({\hat{a}}_m^{†}+{\hat{a}}_m).\end{array}$$

$$\text{Im}Y(\omega )=\frac{1}{\omega L_J}\frac{{\left(\frac{\omega }{{\omega }_L}\right)}^2-1}{1-{\left(\omega \sqrt{L{C}_L}\right)}^2},$$

$${\omega }_L=\frac{1}{\sqrt{(L+L_J)C_L}}.$$

$$\text{Im}\frac{\partial Y}{\partial \omega }({\omega }_L)=2C_L{\left(\frac{L+L_J}{L_J}\right)}^2$$

$$A_L=-\frac{e^2}{2C_L}{\left(\frac{L_J}{L+L_J}\right)}^3.$$

$$\text{Im}Y(\omega )=C_H\omega \left(1-\frac{{\omega }_H^2}{{\omega }^2}\right),$$

$${\omega }_H=\sqrt{\frac{L+L_J}{L{L}_J{C}_H}}.$$

$$\text{Im}\frac{\partial Y}{\partial \omega }({\omega }_H)=2C_H$$

$$A_H=-\frac{e^2}{2C_H}\left(\frac{L}{L+L_J}\right).$$

$$\begin{array}{rl}{\hat{H}}_{4,\text{diag}}=& \hslash {\omega }_{\text{H}}{\hat{a}}^{†}\hat{a}-\frac{A_{\text{H}}}{2}{\hat{a}}^{†}{\hat{a}}^{†}\hat{a}\hat{a}\\ & +\hslash {\omega }_{\text{L}}{\hat{b}}^{†}\hat{b}-\frac{A_{\text{L}}}{2}{\hat{b}}^{†}{\hat{b}}^{†}\hat{b}\hat{b}\\ & -\chi {\hat{a}}^{†}\hat{a}{\hat{b}}^{†}\hat{b},\end{array}$$

$$\begin{array}{rl}\hslash {\omega }_{\text{H}}=\hslash {\bar{\omega }}_{\text{H}}-A_{\text{H}}-\frac{\chi }{2},& \\ \hslash {\omega }_{\text{L}}=\hslash {\bar{\omega }}_{\text{L}}-A_{\text{L}}-\frac{\chi }{2}.& \end{array}$$

$$\frac{\text{d}}{\text{d}t}\hat{a}(t)=\frac{i}{\hslash }[{\hat{H}}_{\text{undr}},\hat{a}(t)]-\frac{\kappa }{2}\hat{a}(t)+\sqrt{{\kappa }_{\text{ext}}}{\tilde{a}}_{\text{in}}(t).$$

$$S_{\text{11}}(t)=\frac{{\tilde{a}}_{\text{out}}(t)}{{\tilde{a}}_{\text{in}}(t)}=1-\sqrt{{\kappa }_{\text{ext}}}\frac{\hat{a}(t)}{{\tilde{a}}_{\text{in}}(t)},$$

$${\tilde{a}}_{\text{in}}(t)=\sqrt{\frac{P_{\text{in}}}{\hslash {\omega }_{\text{d}}}}{e}^{-i{\omega }_{\text{d}}t},$$

$$\begin{array}{rl}\frac{\text{d}}{\text{d}t}\hat{a}(t)& =\frac{i}{\hslash }[{\hat{H}}_{\text{undr}}+{\hat{H}}_{\text{drive}},\hat{a}(t)]-\frac{{\kappa }_{\text{ext}}}{2}\hat{a}(t),\\ \text{where }{\hat{H}}_{\text{drive}}& =i\hslash {ϵ}_{\text{d}}\left(e^{-i{\omega }_{\text{d}}t}{\hat{a}}^{†}(t)-e^{i{\omega }_{\text{d}}t}\hat{a}(t)\right),\\ {ϵ}_{\text{d}}& =\sqrt{\frac{{\kappa }_{\text{ext}}{P}_{in}}{\hslash {\omega }_{\text{d}}}}.\end{array}$$

$$\frac{\text{d}}{\text{d}t}\hat{a}=\frac{i}{\hslash }[U_{\text{probe}}^{†}{\hat{H}}_{\text{undr}}{U}_{\text{probe}}+{\tilde{H}}_{\text{drive}},\hat{a}]-\frac{{\kappa }_{\text{ext}}}{2}\hat{a},$$

$$\begin{gathered} {\tilde{H}}_{\text{drive}}=-\hslash {\omega }_{\text{d}}{\hat{a}}^{†}\hat{a} \\ +i\hslash {ϵ}_{\text{d}}\left({\hat{a}}^{†}-\hat{a}\right). \end{gathered}$$

$${\hat{S}}_{\text{11}}({\omega }_{\text{d}})=1-\frac{{\kappa }_{\text{ext}}}{{ϵ}_{\text{d}}}\hat{a},$$

$$\begin{array}{rl}\frac{\text{d}}{\text{d}t}\hat{a}=& \frac{i}{\hslash }[U_{\text{probe}}^{†}{\hat{H}}_{\text{undr}}{U}_{\text{probe}}+{\hat{H}}_{\text{drive}},\hat{a}]\\ +& \left(L^{†}\hat{a}L-\frac{1}{2}\left(\hat{a}{L}^{†}L+L^{†}L\hat{a}\right)\right),\\ \text{where }L=& {\kappa }_{\text{ext}}\hat{a},\end{array}$$

$$\begin{array}{rl}\frac{\text{d}}{\text{d}t}\rho =& -\frac{i}{\hslash }[U_{\text{probe}}^{†}{\hat{H}}_{\text{undr}}{U}_{\text{probe}}+{\hat{H}}_{\text{drive}},\rho ]\\ +& \left(L\rho L^{†}-\frac{1}{2}\left(\rho L^{†}L+L^{†}L\rho \right)\right)\\ =& -\frac{i}{\hslash }[U_{\text{probe}}^{†}{\hat{H}}_{\text{undr}}{U}_{\text{probe}}+{\hat{H}}_{\text{drive}},\rho ]+{\kappa }_{\text{ext}}\mathcal{L}[\hat{a}],\end{array}$$

$$\begin{array}{rl}\frac{\text{d}}{\text{d}t}\rho =& -\frac{i}{\hslash }[H_{\text{4,diag}}+{\hat{H}}_{\text{drive}},\rho ]\\ & +\kappa (n_{\text{th}}^{(H)}+1)\mathcal{L}[\hat{a}]+\kappa n_{\text{th}}^{(H)}\mathcal{L}[{\hat{a}}^{†}]\\ & +\gamma (n_{\text{th}}+1)\mathcal{L}[\hat{b}]+\gamma n_{\text{th}}\mathcal{L}[{\hat{b}}^{†}].\end{array}$$

$${\rho }_{\text{ss}}={\left[\begin{array}{cccc}{P}_g& 0& 0& \\ 0& P_e& 0& \\ 0& 0& P_f& \\ & & & \ddots \end{array}\right]}_H\otimes {\left[\begin{array}{cccc}{P}_0& 0& 0& \\ 0& P_1& 0& \\ 0& 0& P_2& \\ & & & \ddots \end{array}\right]}_L,$$

$$\rho ={\rho }_{\text{ss}}+{ϵ}_{\text{d}}{\rho }_{\text{pert}},$$

Full text

Detecting and manipulating single photons becomes more difficult at lower frequencies because of thermal fluctuations. A hot environment randomly creates and annihilates photons causing decoherence in addition to creating statistical mixtures of states from which quantum state preparation is challenging. This can be mitigated by using a colder system as a heat sink, to extract the entropy created by the environment. Such a scheme, known as reservoir engineering, was first developed in trapped ions Diedrich et al. (1989), where hot degrees of freedom are cooled via the atomic transitions of ions.

Using superconducting electronics, circuit quantum electrodynamics (cQED) has made extensive use of reservoir engineering to cool, but also manipulate electromagnetic fields at the quantum level. With the prospect of building a quantum computer, or to demonstrate fundamental phenomena, experiments have shown the cooling or reset of qubits to their ground-state Valenzuela et al. (2006); Geerlings et al. (2013); Magnard et al. (2018), also in the megahertz regime Vool et al. (2018), quantum state stabilization Murch et al. (2012); Shankar et al. (2013); Leghtas et al. (2015), and quantum error correction Ofek et al. (2016). Using superconducting circuits, reservoir engineering is commonplace in electromechanical systems Teufel et al. (2011a); Eichler and Petta (2018), but with weak non-linearity, such schemes have only limited quantum control Ockeloen-Korppi et al. (2018); Viennot et al. (2018) compared to typical cQED systems. Despite the many applications of quantum state engineering in cQED, obtaining control over the quantum state of a hot resonator, where the environment temperature is a dominant energy scale, remains a largely unexplored and challenging task.

We directly observe the quantization of radio-frequency electromagnetic fields in a thermally-excited megahertz photonic resonator, and manipulate its quantum state using reservoir engineering. Specifically, we cool the 173 MHz resonator to $90\%$ ground-state occupation, and stabilize one- and two-photon Fock states. Releasing the resonator from our control, we observe its re-thermalization with photon-number resolution.

We use the paradigm of cQED, where a resonator can be read out and controlled by dispersively coupling it to a superconducting qubit. Achieving this with a GHz qubit and MHz photons is challenging, since in a conventional cQED architecture the coupling would be far too weak Gely et al. (2018). To overcome this, we present a circuit enabling a very strong interaction, resulting in a cross-Kerr coupling larger than the qubit and resonator dissipation rates, despite an order of magnitude difference in their resonance frequencies.

The circuit (Fig. 1A) comprises of a Josephson junction ($L_{\text{J}}=41$ nH) connected in series to a capacitor ($C_{\text{L}}=11$ pF) and a spiral inductor ($L=28$ nH). At low frequencies, the parasitic capacitance of the spiral inductor is negligible, and the equivalent circuit (Fig. 1B) has a first transition frequency $\omega_{\text{L}}=2\pi\times 173\ \text{MHz}$. At gigahertz frequencies, $C_{L}$ behaves as a short, and the capacitance of the spiral inductor $C_{\text{H}}=40$ fF becomes relevant instead. The resulting parallel connection of $L_{\text{J}}$, $L$ and $C_{\text{H}}$ (Fig. 1C) has a first transition frequency $\omega_{\text{H}}=2\pi\times 5.91\ \text{GHz}$. The two modes share the Josephson junction. The junction has an inductance that varies with the current fluctuations traversing it, and consequently the resonance frequency of the high-frequency (HF) mode shifts as a function of the number of excitations in the low-frequency (LF) mode and vice versa. This cross-Kerr interaction is quantified by the shift per photon $\chi=2\sqrt{A_{\text{H}}A_{\text{L}}}$, where the anharmonicity of the LF and HF mode $A_{\text{L}}=h\times 495\ \text{kHz}$ and $A_{\text{H}}=h\times 192\ \text{MHz}$ are given by SI

[TABLE]

The system is described by the Hamiltonian SI

[TABLE]

where $\hat{a}$ ($\hat{b}$) is the annihilation operator for photons in the HF (LF) mode. The second line describes the anharmonicity or Kerr non-linearity of each mode. The last term describes the cross-Kerr interaction. By combining it with the first term as $(\hbar\omega_{\text{H}}-\chi\hat{b}^{\dagger}\hat{b})\hat{a}^{\dagger}\hat{a}$, the dependence of the HF mode resonance on the number of photons in the LF mode becomes apparent.

The cross-Kerr interaction manifests as photon-number splitting Schuster et al. (2007) in the measured microwave reflection $S_{11}$ (Fig. 1D). Distinct peaks correspond to the first transition frequency of the HF mode $|g,n\rangle\leftrightarrow|e,n\rangle$, with frequencies $\omega_{\text{H}}-n\chi/\hbar$ where $\chi/h=21$ MHz. We label the eigenstates of the system $|j,n\rangle$, with $j=g,e,f,...$ ($n=0,1,2,...$ ) corresponding to excitations of the HF (LF) mode. The amplitude of peak $n$ is proportional to

[TABLE]

where $P_{n}$ is the occupation of photon-number level $|n\rangle$ in the LF mode and $\kappa_{\text{ext}}/\kappa_{n}$ is the ratio of external coupling $\kappa_{\text{ext}}/2\pi=1.6\cdot 10^{6}s^{-1}$ to the total line-width $\kappa_{n}$ of peak $n$. From the Bose-Einstein distribution of peak heights $P_{n}$, we extract the average photon occupation $n_{\text{th}}=1.6$ corresponding to a mode temperature of $17$ mK.

The resolution of individual photon peaks is due to the condition $\kappa_{n}\ll\chi/\hbar$. The peak line-widths increase with $n$ following $\kappa_{n}=\kappa(1+4n_{\text{th}}^{(H)})+2\gamma(n+(1+2n)n_{\text{th}})$, where $\kappa/2\pi=3.7\cdot 10^{6}s^{-1}$ is the dissipation rate of the HF mode, $n_{\text{th}}^{(H)}\simeq 0.09$ its thermal occupation (see Fig. S10), and $\gamma/2\pi=23\cdot 10^{3}s^{-1}$ is the dissipation rate of the LF mode (obtained through time-domain measurement Fig. 4). The condition $\kappa_{n}\ll A_{\text{H}}/\hbar$ makes the HF mode an inductively-shunted transmon qubit Koch et al. (2007), making it possible to selectively drive the $|g,n\rangle\leftrightarrow|e,n\rangle$ and $|e,n\rangle\leftrightarrow|f,n\rangle$ transitions. Despite its low dissipation rate $\gamma$, the LF mode has a line-width of a few MHz (measured with two-tone spectroscopy, Fig. S15) which originates in thermal processes such as $|g,n\rangle\rightarrow|e,n\rangle$ occurring at rates $\sim\kappa n_{\text{th}}^{(H)}$ larger than $\gamma$ SI . The LF mode line-width is then an order of magnitude larger than $A_{\text{L}}$, making it essentially a harmonic oscillator that we will refer to as the resonator.

The junction non-linearity enables transfer of population between states by coherently pumping the circuit at a frequency $\omega_{\text{p}}$. The cosine potential of the junction imposes four-wave mixing selection rules, only allowing interactions that involve 4 photons. One such interaction is

[TABLE]

activated when driving at the energy difference between the two coupled states $\omega_{\text{p}}=2\omega_{\text{H}}-\omega_{\text{L}}-\left(2n\chi+A_{\text{H}}\right)/\hbar$. This process, enabled by a pump photon, annihilates a photon in the resonator and creates two in the transmon. The number of photons involved in the interaction is four, making it an allowed four-wave mixing process. The induced coupling rate is $g=A_{\text{H}}^{\frac{3}{4}}A_{\text{L}}^{\frac{1}{4}}\xi_{\text{p}}$, where $|\xi_{\text{p}}|^{2}$ is the amplitude of the coherent pump tone measured in number of photons SI .

We use this pump tone in combination with the large difference in mode relaxation rates to cool the megahertz resonator to its ground-state (Fig. 2A). The pump drives transitions between $|g,1\rangle$ and $|f,0\rangle$ at a rate $g$. The population of $|g,1\rangle$, transfered to $|f,0\rangle$, subsequently decays at a rate $2\kappa$ to the ground-state $|g,0\rangle$. Cooling occurs when the thermalization rate of the resonator $n_{\text{th}}\gamma$ is slower than the rate $C\gamma$ at which excitations are transfered from $|g,1\rangle$ to $|g,0\rangle$, where $C=2g^{2}/\kappa\gamma$ is the cooperativity (proportional to cooling-pump power SI ).

For different cooling pump strengths, we measure $S_{11}$ (Fig. 2B). The pump frequency is adapted at each power since the AC-stark effect increasingly shifts the qubit frequency as a function of power (see Fig. S9). The data is fitted to a sum of complex Lorentzians, with amplitudes given by Eq. (3) and line-widths $\kappa_{n}$, from which $P_{n}$ is extracted. Thermal effects lead to the ratio $P_{n+1}/P_{n}=n_{\text{th}}/(1+n_{\text{th}})$ between neighboring photon-number states for $n\geq 1$, and the cooling pump changes the ratio of occupation of the first two states

[TABLE]

The ground-state occupation hence increases with cooperativity and we attain a maximum $P_{0}=0.82$. At higher cooperativity, $P_{0}$ diminishes due to the off-resonant driving of other four-wave mixing processes such as $|f,n+1\rangle\langle g,n|+h.c.$ which tend to raise the photon number of the resonator. This effect is simulated using an adaptive rotating-wave approximation Baker et al. (2018) (Fig. 2C and S6).

Neighbouring four-wave mixing processes are measured by sweeping the pump frequency whilst monitoring the spectrum (Fig. 3A). When cooling with a single pump they eventually limit performance, but can be resonantly driven to our advantage. By driving multiple cooling interactions $|g,n\rangle\leftrightarrow|f,n-1\rangle$, less total pump power is required to reach a given ground-state occupation, hence minimizing off-resonant driving. By maximizing the ground-state peak amplitude as a function of the power and frequency of four cooling tones, we achieve $P_{0}=0.90$ (Fig. 3B).

By combining cooling $|g,n\rangle\leftrightarrow|f,n-1\rangle$ and raising $|g,n\rangle\leftrightarrow|f,n+1\rangle$ tones (inset of Fig. 3C), we demonstrate stabilization of higher Fock states, non-Gaussian states commonly considered as non-classical phenomena Rips et al. (2012). The optimum frequencies for the raising and cooling tones adjacent to the stabilized state were detuned by a few MHz from the transition frequency (see dashed lines in the inset of Fig. 3C), otherwise one pump tone would populate the $|f\rangle$ level, diminishing the effectiveness of the other.

Finally we investigate dynamics in a photon resolved manner (Fig. 4). Whilst probing $S_{11}$ at a given frequency, we switch the cooling or single photon stabilization pumps on and off for intervals of $50\ \upmu$s. We perform this for a sequence of probe frequencies, resulting in $S_{11}$ as a function of both frequency and time (see full spectrum in SI ). The spectrum is fitted at each time to extract $P_{n}$ as a function of time. After reaching the steady state, the pumps are turned off and we observe the thermalization process which follows the semi-classical master equation

[TABLE]

Our cQED architecture enables the readout and manipulation of a radio-frequency resonator at the quantum level. Utilizing the fast readout methods of cQED, single-shot readout or the tracking of quantum trajectories could enable even finer resolution of thermodynamic effects at the quantum scale. Coupling many of these devices together could enable the exploration of many-body effects in Bose-Hubbard systems with dynamically tunable temperatures Rigol et al. (2008); Sorg et al. (2014). This circuit architecture could also be used to interface circuit quantum electrodynamics with different physical systems in the MHz frequency range, such as spin systems Ares et al. (2016) or macroscopic mechanical oscillatorsTeufel et al. (2011a). Finally, this circuit could enable sensing of radio frequency radiation with quantum resolution, a critical frequency range for a number of applications, from nuclear magnetic resonance imaging to radio astronomy.

**Acknowledgments: ** We acknowledge Ya. M. Blanter, S. M. Girvin, J. D. P. Machado for useful discussions. **Funding: ** This work was supported by the European Research Council under the European Union’s H2020 program under grant agreements 681476 - QOM3D and 785219 - GrapheneCore2, by the Dutch Foundation for Scientific Research (NWO) through the Casimir Research School, and by the Army Research Office through Grant No. W911NF-15-1-0421. **Author contributions: ** MFG and RV developed the theoretical description of the experiment. MFG designed the device. MFG fabricated the device with help from JD and MK. MFG, MK, CD, JD and MJ participated in the measurements. MFG and CD analyzed the data. BB provided the software and input for the adaptive rotating-wave-approximation simulation. MFG wrote the manuscript with input from MK, CD and GAS. All co-authors reviewed the manuscript and provided feedback. GAS supervised the project. **Competing interests: ** The authors declare that they have no competing interests. **Data and materials availability: ** Raw data as well as all measurement, data-analysis and simulation code used in the generation of main and supplementary figures is available in Zenodo with the identifier 10.5281/zenodo.2551258

1 Fabrication

The circuit is fabricated on a high resistivity silicon substrate using aluminum as a superconductor. Using shadow evaporation, we first pattern Al/AlOx/Al Josephson junctions, the bottom plate of the capacitor and an underpass (a line connecting the center of the spiral inductor to the capacitor). We then deposit $260$ nm of hydrogenated amorphous silicon (a-Si:H) as a dielectric layer, motivated by its expected low dielectric loss O’Connell et al. (2008), using PECVD (plasma enhanced chemical vapor deposition) The a-Si:H is patterned to form a dielectric layer for the parallel plate capacitor, a bridge over the spirals underpass, and a protection layer above the junctions. Finally we sputter-deposit and pattern aluminum to form the rest of the circuitry, after an argon-milling step to ensure a galvanic connection to the first aluminum layer. The resulting circuit is shown in detail in Fig. S1.

2 Experimental setup

3 Data filtering

In Figs 1A, 3A, S4B, S10B,C, S14, S13A,B,C we applied a Gaussian filter with a standard deviation of one increment in the x-axis (and y-axis when applicable). The filtering was used in the construction of the figure for clarity. No filtering was applied before fitting the data.

4 Theory

4.1 Circuit Hamiltonian

In this section, we derive the Hamiltonian for the circuit shown in Fig. S2A using the black-box quantization method Nigg et al. (2012). This method allows the systematic derivation of the resonance frequency $\bar{\omega}_{m}$ and anharmonicity $A_{m}$ of the different modes $m$ of a circuit from the admittance $Y(\omega)=1/Z(\omega)$ across the Josephson junction if we replace the latter by a linear inductor $L_{\text{J}}=\hbar^{2}/4e^{2}E_{\text{J}}$. The resonance frequencies $\bar{\omega}_{m}$ are the zeros of the admittance $Y(\bar{\omega}_{m})=0$, and the anharmonicities are given by

[TABLE]

The idea is to quantify through $A_{m}$ the amount of current traversing the Josephson junction for an excitation in mode $m$. The Hamiltonian of the circuit is then

[TABLE]

In the circuit of Fig. S2A , there are two modes, a high-frequency one and a low-frequency one. By comparing to a black-box quantization of the full circuit, we find that taking the approximation of $C_{H}\simeq 0$, $C_{c}\simeq 0$ for the low-frequency mode and $C_{L}\simeq\infty$ for the high-frequency mode results in corrections of only $0.2$, $1.2$, $0.3$ and $2.1$ % in the value of $\omega_{L}$, $\omega_{H}$, $A_{L}$ and $A_{H}$ respectively. It is therefore a good approximation, which has the additional advantage of producing simple analytical equations for the frequencies and anharmonicities of the circuit. Starting with the low-frequency mode shown in Fig. S2B, we find the (imaginary part of the) admittance across the linearized junction to be

[TABLE]

yielding the resonance frequency

[TABLE]

Taking the derivative of the imaginary part of the admittance at $\omega=\omega_{L}$ yields:

[TABLE]

Substituting this into Eq. (S1) yields

[TABLE]

Turning to the high-frequency mode shown in Fig. S2C, we find the (imaginary part of the) admittance across the linearized junction to be

[TABLE]

yielding the resonance frequency

[TABLE]

Taking the derivative of the imaginary part of the admittance at $\omega=\omega_{H}$ yields:

[TABLE]

Substituting this into Eq. (S1) yields

[TABLE]

A Taylor expansion of the junctions cosine potential is justified if the anharmonicities are weak and only a few photons populate the circuit. Whilst numerical calculations in this work consider the 8-th order expansion, much understanding can be gleaned by stopping the expansion at the fourth-order

[TABLE]

where $\chi$ is the cross-Kerr coupling: the amount by which the high-mode transition shifts as a result of adding an excitation in the low mode and vice versa. We defined the first transition frequencies of both modes

[TABLE]

In Eq. (LABEL:eq:H4_diag), we have neglected terms in the expansion which are off-diagonal in the Fock basis and do not modify the eigenergies to leading order perturbation theory. The eigenfrequencies of the system are summarized in the energy diagram of Fig. S3

4.2 Comparison to the typical circuit QED architecture

We now compare our circuit architecture with a more conventional circuit QED coupling scheme Koch et al. (2007), where the transmon qubit with a frequency $\omega_{\text{H}}$ couples capacitively at a rate $G$ to an $LC$-oscillator ($\omega_{\text{L}}$). In this architecture, the cross-Kerr coupling would be $4A_{\text{H}}(\bar{g}\omega_{\text{L}}/\omega_{\text{H}}^{2})^{2}$, to first order in $\bar{G}$ and $A_{\text{H}}$ Gely et al. (2018). If we would want a cross-Kerr coupling to exceed $\kappa$, for the large frequency difference $\omega_{\text{L}}\ll\omega_{\text{H}}$ demonstrated in this work, the couplings $G$ would have to be very large. As is well known from ultra-strong coupling circuit QED architectures, this translates to both high impedance resonators Bosman et al. (2017a) and large coupling capacitors Bosman et al. (2017b). These elements are all present in this circuit if we consider $C_{L}$ as a coupling capacitor between the high impedance $L_{\text{H}}C_{\text{H}}$-oscillator and the qubit constituted of the junction and the junctions capacitance (that we have neglected in the previous Hamiltonian derivation). However, the basis of normal modes presented in the previous section present a much more convenient framework to understand the system.

4.3 Translating the measured $S_{11}$ to a quantum operator

We now introduce a driving term in the Hamiltonian and consider losses to both the environment and the measurement port. Following input-output theory Vool and Devoret (2017); Clerk et al. (2010), the quantum Langevin equation for $\hat{a}(t)$ is

[TABLE]

Where the undriven Hamiltonian $\hat{H}_{\text{undr}}$ corresponds to that of Eq. (S2), where the degree of expansion of the non-linearity is yet unspecified. The microwave reflection measured in spectroscopy (here in the time-domain) is given by

[TABLE]

where $\tilde{a}_{\text{in}}(t)$ ($\tilde{a}_{\text{out}}(t)$) is the incoming (outgoing) field amplitude, $\kappa_{\text{ext}}$ ($\kappa$) is the external (total) coupling rate of the high-frequency mode. The coupling of the low mode to the feedline $\gamma_{\text{ext}}/2\pi=2s^{-1}$ is much smaller than coupling of the high mode to the feedline $\kappa_{\text{ext}}/2\pi=1.63\cdot 10^{6}s^{-1}$, we therefore assume that a drive tone only affects the high-frequency mode. For a coherent drive, characterized by a drive frequency $\omega_{\text{d}}$ and an incoming power $P_{\text{in}}$ (equal to the average power $\langle P(t)\rangle$ of the oscillating input signal), the wave amplitude is

[TABLE]

and the drive term can be incorporated in the Hamiltonian of the system

[TABLE]

Additionally, we also remove the time-dependence in the drive Hamiltonian by moving to a frame rotating at $\omega_{\text{d}}$ with the unitary transformation $U_{\text{probe}}=e^{i\omega_{\text{d}}t\hat{a}^{\dagger}\hat{a}}$,

[TABLE]

where $\hat{a}e^{i\omega_{d}t}=\hat{a}(t)$ and

[TABLE]

In this rotating frame, the reflection coefficient becomes

[TABLE]

of which we measure the expectation value when probing the system. From now on, and in the main text we use the shorthand $S_{11}(\omega_{\text{d}})=\langle\hat{S}_{\text{11}}(\omega_{\text{d}})\rangle$.

Note that by casting the quantum Langevin Eq. (S17) in the form

[TABLE]

it can be readily transformed to a Lindblad equation

[TABLE]

better suited to numerical calculations using QuTiP Johansson et al. (2012).

4.4 Derivation of $S_{11}[\omega]$: probing the system

In this section we derive the spectrum of the high mode for arbitrary states of the low mode. We append the Lindblad equation of Eq. (S21) to take into account additional interactions of the system with the environment. Internal dissipation of the high mode $\kappa_{\text{int}}$, is added to the external dissipation rate to constitute its total dissipation rate $\kappa=\kappa_{\text{int}}+\kappa_{\text{ext}}$. The low mode is attributed a dissipation rate $\gamma$. The average thermal occupation of the two modes are denoted by $n_{\text{th}}^{(H)}$ and $n_{\text{th}}$ for the high and low mode respectively. We can estimate the response function $S_{11}(\omega_{\text{d}})$ analytically using the Hamiltonian of Eq. (LABEL:eq:H4_diag). The unitary $U_{\text{probe}}$ leaves this Hamiltonian unchanged and the complete Lindblad equation is then

[TABLE]

In the un-driven case $\epsilon_{\text{d}}=0$, we assume the steady-state solution to be a diagonal density matrix $\rho_{\text{ss}}$ as a consequence of thermal effects

[TABLE]

where $P_{g},P_{e},P_{f}$ ($P_{0},P_{1},P_{2}$) corresponds to the occupation of high (low) mode levels. Note that when we pump the system, effectively coupling levels of the high and low mode, this approximation breaks down and that particular limit is discussed below. We now look for a perturbative correction to this matrix at a small driving rate $\epsilon_{\text{d}}$

[TABLE]

where $\rho_{\text{pert}}$ has the unit time. The objective is to determine the expectation value of the reflection coefficient

[TABLE]

We substitute the perturbative expansion of $\rho$ into Eq. (S22) and keep only terms to first order in $\epsilon_{\text{d}}$. This equation is solved analytically in reduced Hilbert-space sizes using the software Wolfram Mathematica. The largest Hilbert-space sizes for which Mathematica could provide an analytical solution in a reasonable amount of time were: (4,0), (3,2), (2,5) where the first (second) number designates the number of levels included in the high (low) mode. We extrapolate the obtained results to construct the reflection coefficient

[TABLE]

which corresponds to a sum of Lorentzian functions, each associated to high-mode level $i$ and a low-mode level $n$, with line-width $\kappa_{i,n}$ centered around $\Delta_{i,n}=0$, where

[TABLE]

Note that in the main text we use the notation $\kappa_{g,n}=\kappa_{n}$.

By numerically computing $S_{11}$ as described in Sec. 5.2, we find that the expression for the line-widths $\kappa_{i,n}$ remains accurate, whilst the center of the Lorentzians $\Delta_{i,n}$ will slightly shift from Eq. (S27), as shown in Fig. S4B. When fitting data, we hence use the Eqs. (S26) whilst fixing $\Delta_{i,n}$ with a diagonalization of the Hamiltonian Eq. (S2) Taylor expanded to the 8-th order. In Fig. S10(C), we show that Eq. (S26) is in excellent agreement with both data and numerics.

4.4.1 The impact of a pump tone on $S_{11}(\omega)$

Pump tones can invalidate Eq. (S26) in different ways. As an example let us take the cooling scheme where a pump tone couples the levels $\ket{g,1}$ and $\ket{f,0}$ at a rate $g$. This is simulated by numerically finding the steady state of the Hamiltonian

[TABLE]

written in a frame rotating at the probe frequency $\Delta=\omega_{H}-\omega_{d}$ and where the levels $\ket{g,1}$ and $\ket{f,0}$ are made resonant. As shown in Fig. S4A, a peak corresponding to a transition to or from a level which is being pumped will be broadened in line-width and eventually will split into two peaks with increasing $g$. This is not an issue in the cooling scheme since we do not use the driven $n=1$ peak to extract Fock-state fidelity, only the $n=0$ peak.

We do however off-resonantly pump $\ket{g,0}\leftrightarrow\ket{f,1}$ for example, along with many other transitions involving either state $\ket{g,0}$ or $\ket{e,0}$. Off-resonant pumping should also lead to line-width broadening, this time of a peak used in extracting a Fock state fidelity. To mitigate this issue we extract $P_{n}$ – when stabilizing the $n$-th Fock state – by using a fixed line-width $\kappa_{n}$ defined in Eq. (S26). This means that we always give a lower bound to $P_{n}$. By comparing the pumped and un-pumped line-width of $n=0$ peak (see Fig. S9(B)), we notice no change in line-width with increasing pump power, indicating that our underestimation is certainly not drastic.

Finally, pump tones could drive the steady-state away from our assumption of a purely diagonal density matrix Eq. (S23). However we find that in the cooling experiment of Fig. 2, the adaptive rotating-wave simulation suggests that at maximum $P_{0}$, all off-diagonal terms of the density matrix are below $2.3\times 10^{-3}$. This issue can safely be disregarded.

4.5 Four wave mixing

4.5.1 Analytical derivation of the pump-induced coupling rates

In this section we will consider the probe tone to be very weak and hence negligible. Following Refs. Leghtas et al. (2015), we add a pump tone driving the high mode with frequency $\omega_{\text{p}}$ and strength $\epsilon_{\text{p}}$ to the system Hamiltonian

[TABLE]

We move to the displaced frame of the pump through the unitary transformation

[TABLE]

Where $\tilde{\xi}_{\text{p}}$ is defined by the differential equation

[TABLE]

For $t\gg 1/\kappa$, and for far detuned drives $|\omega_{\text{H}}-\omega_{\text{p}}|\gg\kappa$, this equation is solved by

[TABLE]

In this frame, the Hamiltonian becomes

[TABLE]

We now Taylor expand the cosine non-linearity to fourth-order, neglecting terms which are off-diagonal in the Fock basis except when they depend on $\tilde{\xi}_{\text{p}}$ . The latter can be made relevant depending on our choice of $\omega_{\text{p}}$.

[TABLE]

Where $\hat{H}_{4,\text{diag}}$ was given in Eq. (LABEL:eq:H4_diag). The terms dependent on the pump power and frequency are assembled in the term $\hat{H}_{\text{p}}$ and written in Table S1, along with the approximate pumping frequency $\omega_{\text{p}}$ necessary to eliminate their time-dependence. As shown in the next paragraph, this occurs when the pump frequency matches the transition frequency between the two states coupled by the interaction term.

We now move to the interaction picture through the unitary transformation

[TABLE]

$\hat{H}_{4,\text{diag}}$ is diagonal in the Fock state basis $\left\{\ket{j,n}\right\}_{\begin{subarray}{c}n=0,1,2,..\\ j=g,e,f,..\end{subarray}}$

[TABLE]

To determine $\hat{H}_{\text{p}}$ in this frame, it suffices to know the expression of annihilation operators in this frame. We will take as an example the term we use for cooling, which reads in the interaction picture

[TABLE]

Since $\hat{H}_{0}$ is diagonal, exponentiating it only requires exponentiating each of the diagonal elements, and the annihilation operators in the interaction picture are

[TABLE]

Note that if the system were harmonic, these expressions would simplify to $e^{-i\omega_{\text{H}}t}\hat{a}$ and $e^{-i\omega_{\text{L}}t}\hat{b}$. If we substitute Eqs. (S38) into Eq. (S37), one of the terms we obtain is

[TABLE]

where we defined the interaction strength

[TABLE]

By choosing the pump frequency $\omega_{\text{p}}=2\omega_{\text{H}}-A_{\text{H}}/\hbar-\omega_{\text{L}}$, the term becomes time-independent, making it more relevant than the other terms of $\hat{H}_{P}$ as we will derive next. More generally, we can engineer the cooling interactions

[TABLE]

by choosing the pump frequencies

[TABLE]

This is the interaction used in all expriments presented in the last three figures of the main text. Cooling by driving the $|g,1\rangle\leftrightarrow|e,0\rangle$ transition may seem like a more natural choice, but it is a two pump-photon process (due to four-wave mixing selection rules), and hence requires higher pumping power. Additionally, due to its higher energy, the $|f,0\rangle$ state has a lower thermal occupation than $|e,0\rangle$. As discussed below, high pump powers and thermal occupation of the qubit place strong limitations on the cooling efficiency.

Rather than lowered, the number of excitations in the low mode can also be raised using interactions of the form

[TABLE]

which are realized by choosing the pump frequencies

[TABLE]

4.5.2 Derivation of cooling rate

In this section we focus on the cooling interaction of Eq. (S41), however the methodology described is generalizable to all interaction terms. The objective of this section is to translate the interaction term derived previously into a cooling rate for the low mode. We assume that this interaction is sufficiently weak to enable us to perform first-order perturbation theory, considering the high mode as a fluctuating quantum noise source $\hat{F}_{H}$ perturbing the low mode following App. B.1 of Ref. Clerk et al. (2010). An initial state of the low mode $|n\rangle$ will evolve following

[TABLE]

where $\hat{F}_{H}(\tau)=\left(\ket{f}\bra{g}\right)(\tau)$ is treated as an independent noise source acting on the Hilbert space of the high mode. We consider the transition is off-resonantly driven such that the time-dependence in the interaction picture is not completely eliminated and the interaction term rotates at

[TABLE]

The probability amplitude of finding the low mode in $|n-1\rangle$ is

[TABLE]

leading to a probability

[TABLE]

Note that $|\langle n-1|\psi(t)\rangle|^{2}$ is still a quantum operator acting on the high-mode Hilbert space. To obtain a classical probability, we now calculate its expectation value $\langle.\rangle_{H}$, provided that the high mode evolves in steady-state under thermal effects and dissipation

[TABLE]

As in Appendix A.2 of Clerk et al. (2010), we transform the double integral $S$ to

[TABLE]

For time-scales larger than the decay rate of the high mode $\tau\gg 1/\kappa$, the two time-dependent high-mode operators are not correlated and the integrand will vanish (see Appendix A.2 of Clerk et al. (2010)). We can therefore extend the range of the inner integral to $\pm\infty$ in estimating the probability at a time $t\gg 1/\kappa$.

[TABLE]

Using time-translation invariance, we can remove the dependence on $T$

[TABLE]

such that the rate becomes time-independent

[TABLE]

Using time-translation invariance, we find that for negative values of $\tau$,

[TABLE]

leading to

[TABLE]

In the steady state of the system, the quantum regression theorem can be shown to reduce the expression to

[TABLE]

where $\hat{\rho}$ is the steady-state density matrix of the high mode and $e^{\mathcal{L}t}$ its propagator, a function which takes a density matrix as an input and evolves it up to a time $t$ following the Lindblad equation. Reducing the high mode to a three-level system and considering dissipation and thermal effects, this trace can be calculated analytically using the QuantumUtils Mathematica library

[TABLE]

By only considering dissipation and thermalization, we made the assumption that an excitation could not be driven back from $\ket{f,n-1}$ to $\ket{g,n}$ under the effect of pumping, *i.e. *we assume $2\kappa\gg\sqrt{n}g$, that we are far from the strong coupling regime. After integration, we obtain

[TABLE]

Following the same method, we also obtain for the hermitian conjugate of this interaction term

[TABLE]

if the $\ket{f}$ level is populated, we find that there is a probability for the pump to raise the number of excitations in the low mode rather than lower it. We refer to the steady state population of the ground and second-excited state of the high mode as $P_{g}$ and $P_{f}$ respectively. The same calculation can be performed for the raising interaction, which yields identical rates only with $P_{g}$ and $P_{f}$ interchanged. A good figure of merit of the cooling efficiency is then to compare this rate with $\gamma$, yielding the cooperativity

[TABLE]

4.5.3 Semi-classical description of the cooling process

With the cooling rate above, we can construct a semi-classical set of rate equations describing the competition between thermalization and cooling. They would correspond to the diagonal part of a Lindblad equation, and equates the population leaving and arriving to a given state of the low mode. We restrict ourselves to the driving of $\ket{f,0}\leftrightarrow\ket{g,1}$ as in the experiment of Fig. 2, where these equations can be written as

[TABLE]

[TABLE]

and, for $n\geq 2$

[TABLE]

In steady state ($\dot{P}=0$), the solution is a function of $P_{0}$

[TABLE]

We reach a unique solution by imposing $\sum_{n}P_{n}=1$, which yields an expression for $P_{0}$

[TABLE]

This expression is used in Fig. S5 to show the temperature limited evolution of $P_{0}$ as a function of cooperativity. A more accurate description of the cooling process at high cooperativities comes from a numerical simulation taking the strong coupling limit and off-resonant driving of other four-wave mixing processes into account.

4.6 Limiting factors to cooling

Here, we discuss three limiting factors to the cooling experiment (Fig. 2), ending with some notes on how to improve the device cooling performance The first limiting factor is the thermal occupation of the high mode. The pump tone drives the population from $|g,1\rangle$ to $|f,0\rangle$, but the reverse process also occurs since the $f$ level has a small thermal population $P_{f}\simeq 0.006$ (see Sec. 7.2). This leads to the limit $P_{1}/P_{0}>P_{f}/P_{g}$ (dashed line in Fig. S5) for which we have derived an exact analytical expression (Eq. (S65))

The second limiting factor is that of strong coupling (similar to in optomechanical cooling Teufel et al. (2011b)), where the pump hybridizes the $|g,1\rangle$ and $|f,0\rangle$ states. If $g$ exceeds the decay rate $2\kappa$, the population of state $|g,1\rangle$ will be driven to $|f,0\rangle$ and then transfered back to $|g,1\rangle$ without having the time to decay to $|e,0\rangle$. To simulate this effect, we compute the steady state of the system by solving a Lindblad equation numerically (see  5.3). The result is shown as a dotted line in Fig. S5, which additionally takes into account the population of the high mode. As with the thermal effect, the strong coupling limit only imposes an upper bound on $P_{0}$, rather than predicting its decrease at high $C$.

When the cooling tone is detuned by $\Delta$ from its transition frequency, the cooperativity acquires a factor $1/\left(1+\Delta^{2}/\kappa^{2}\right)$ (Eq. (S59)). A similar formula applies to all other four-wave mixing processes, including raising interactions (Eq. (S43)). If the latter are far-detuned, their off-resonant driving will have little impact on the system. However, as the cooling process starts to saturate due to the previously discussed limiting factors, the driving of other transitions is still far from saturation and can overpower the cooling effect. What ensues is a competition between off-resonantly driven transitions that cool and raise the photon occupation. We simulate this by following the bootstrap step of the adaptive rotating-wave approximation method of Ref. Baker et al. (2018), which offers a way to include the most relevant off-resonantly driven transitions to the system Hamiltonian (see Sec. 5.3). The result is shown as the solid curve of Fig. S5 which predicts the maximum $P_{0}$ and the strong cooperativity behavior. We emphasize that, except for a small shift on the calibrated cooperativity-axis, the theoretical curves do not correspond to a fit to the data, but rather constitute a prediction based on the independently determined dissipation rates, thermal occupations and circuit parameters. From this simulation we extract that, at maximum $P_{0}=0.82$, the average photon number in the cooled resonator is $\bar{n}=0.65$. Note that

[TABLE]

The first 10 most populated levels are: $P_{g0}=0.736$, $P_{e0}=0.067$, $P_{g1}=0.036$, $P_{g2}=0.028$, $P_{g3}=0.028$, $P_{g4}=0.022$, $P_{g5}=0.017$, $P_{g6}=0.014$, $P_{g7}=0.011$, $P_{f0}=0.009$, where $P_{j,n}$ refers to the occupation of state $\ket{j,n}$. Taking only the contribution of these states into account in the above formula already leads to $\bar{n}=0.51$, and including the occupation of all 50 simulated levels leads to $\bar{n}=0.65$.

Determining the ideal system parameters to improve cooling (and Fock-state stabilization fidelity) is not straightforward. One path to improvement could lie in determining values of $A_{H}$ and $\chi$ which minimize the effect of off-resonant driving by moving the most problematic transitions away from the cooling frequency. Another is to reach a higher ground-state occupation before being limited by strong coupling, which can only be achieved by reducing the resonators dissipation $\gamma$. Decreasing the high mode dissipation $\kappa$ is not necessarily beneficial: it diminishes off-resonant driving, but strong coupling would occur at smaller pump powers. For our system, decreasing $\kappa$ in the simulation of Fig. 2C results in a lower ground-state occupation.

5 Numerical procedures

5.1 Spectrum

The eigenfrequencies of the system are determined by diagonalizing the system Hamiltonian. Unless specified otherwise, we diagonalize the Hamiltonian of Eq. S2 with the junction non-linearity Taylor expanded to 8-th order. We consider 10 excitations in the high mode and 20 in the low mode, and have verified that extending the Hilbert space further only leads to negligeable changes in the obtained spectrum. This diagonalization also provides the dressed eigenstates $\ket{j,n}$, which are to be distinguished from the bare eigenstates $\widetilde{\ket{j,n}}_{\begin{subarray}{c}n=0,1,2,..\\ j=g,e,f,..\end{subarray}}$.

5.2 Microwave reflection

In order to compute the microwave reflection of the device, we solve a Lindblad equation using Qutip Johansson et al. (2012). The Hamiltonian is written in the dressed basis defined above, it is hence diagonal with entries corresponding to the eigenfrequencies obtained in the diagonalization. We consider 5 high-mode excitations and 10 low-mode excitations. We add the drive term $i\hbar\epsilon_{d}(\hat{a}^{\dagger}-\hat{a})$ defined in the dressed basis, and move to the frame rotating at the drive frequency $\omega_{d}$ by adding $-\hbar\omega_{d}\hat{a}^{\dagger}\hat{a}$. We add jump operators defined in the dressed basis by

[TABLE]

to describe dissipation and thermal effects. Finally, we compute the expectation value of $\hat{S}_{11}=1-\frac{\kappa_{\text{ext}}}{\epsilon_{\text{d}}}\hat{a}$ for different drive frequencies. As shown in Fig. S10, this computation is in excellent agreement with the sum of Lorentzian formula of Eq. (S26).

5.3 Cooling simulation

We use a similar method for the adaptive rotating-wave approximation (aRWA) simulation of Fig. 2. We start with the same diagonal Hamiltonian. We denote by $\omega_{j,n}$ the eigenfrequency of the dressed eigenstates $\ket{j,n}$. As a result of the collapse operators of Eq. (S67), a dressed state of the system $\ket{j,n}$ will have a total decay rate to other states of the system

[TABLE]

Following Ref. Baker et al. (2018), we can then estimate the impact of a pump tone at a frequency $\omega_{p}$ and driving rate $\epsilon_{p}$ on the steady state of the system. Two states $\ket{k}=\ket{j,n}$ and $\ket{k^{\prime}}=\ket{j^{\prime},n^{\prime}}$ will be coupled by this pump. And to first order in $\epsilon_{p}$, the only change in the steady state density matrix will be in its off-diagonal element

[TABLE]

where $P_{k}$ is the occupation of state $\ket{k}$ under the collapse operators of Eq. (S67). The dipole moment $V_{kk^{\prime}}=\bra{k}\epsilon_{p}(\tilde{a}+\tilde{a}^{\dagger})\ket{k^{\prime}}$ is computed using annihilation $\tilde{a}$ and creation operators $\tilde{a}^{\dagger}$ defined in the bare basis. The transitions between all the states are then ranked with decreasing $|\rho_{kk^{\prime}}|$ (i.e. decreasing relevance). The most relevant terms are added in the form $\hbar V_{kk^{\prime}}\ket{k}\bra{k^{\prime}}$ to the Hamiltonian which is moved to the rotating frame in which states $\ket{k}$ and $\ket{k^{\prime}}$ are resonant.

In Fig. S5, we perform this calculation for $\omega_{p}=\omega_{f,0}-\omega_{g,1}$. We show both the result of including a maximum number of transitions (465) and a single transition. It was only possible to include 465 transitions out of the 650 transitions which have a non zero dipole moment. This is due to limitations in the construction of the rotating frame, for more details see Ref. Baker et al. (2018).

In Fig. S6, we study how each transition affects the steady-state of the system. Ranking using $|\rho_{kk^{\prime}}|$ does not take into account that multiple transitions may interact. To rank the relevance of the transitions in a more realistic way, we further rank the transitions following their impact on $P_{0}$. We add transitions one by one in the simulation, recording for each transition the change $\Delta P_{0}$ that ensues. We then rank the transitions with decreasing $|\Delta P_{0}|$, and repeat: we add the transitions in the new order one by one, rank them and start again until reaching convergence.

This simulation is in good agreement with the data except at the very highest powers (see Fig. S6D). There are four possible limitations in our aRWA simulation that could be the cause of this discrepancy. First, our implementation of aRWA does not take into account the AC-Stark shift of each level. Present only at high powers, these AC-Stark shifts could bring certain transitions in or out of resonance with each-other, modifying the final steady-state of the system. Secondly, we work with a Hilbert space of only 10 excitations in the low mode. At the highest power, the simulation indicates an average low-mode occupation of $~{}5$ and a larger Hilbert space may be needed to reach more accurate results. Thirdly, only first order transitions were considered in the ranking of the transitions, so no higher order processes, such as those shown in Fig. S14C, are taken into account. Fourthly, we rank transitions with Eq. (S69) using $P_{k}$ the occupation of states $\ket{k}$ under the collapse operators of Eq. (S67). However $P_{k}$ may change under the effect of the driving, modifying the relevance of a given transition. This can be taken into account as described in Ref. Baker et al. (2018), but is too computationally expensive with the Hilbert-space size used here.

6 Background subtraction

6.1 Network analysis

Most of our data analysis relies on fitting a sum of complex Lorentzians (see Eq.( S26)), to the measured microwave reflection $S_{11}$ in both phase and amplitude. The signal we acquire is affected by the imperfections of the microwave equipment used to carry the signals to and from the device.

These can be modeled by a two port network with $s$ parameters $s_{11},s_{22}$, corresponding to the reflections at the VNA ports (reflected back to the VNA) and at the device (reflected back to the device) respectively, and $s_{21},s_{12}$, corresponding to the attenuation chain from the VNA to the device and the amplification chain from the device to the VNA respectively.

We hence measure with our VNA the effective microwave reflection

[TABLE]

Note that these $s$ parameters are generally frequency dependent. We make the approximation $s_{11},s_{22}\ll s_{12},s_{21}$, meaning we attribute most of the measured microwave background to the frequency dependent transmission of the attenuation and amplification chain. The signal we want to measure is now proportional to a so-called “microwave background”

[TABLE]

which we have to experimentally measure.

6.2 Measuring the microwave background

As shown in Fig. S8, when probing the system at high power the device response is $S_{11}=1$, allowing us to extract the microwave background $s_{12}s_{21}$. This phenomenon is a consequence of super-splitting as explained in Bishop et al. (2009), which we will briefly summarize here.

To understand super-splitting, we have to truncate the high mode to a two-level system constituted of its two first levels $\ket{g}$ and $\ket{e}$. In the Bloch sphere, the probe tone will cause rotations around the y-axis and $1-S_{11}$ corresponds to the projection of the state vector on the x-axis. For driving rates faster than $\kappa$, the state vector will rapidly rotate around the y-axis yielding a zero projection on the x-axis hence $S_{11}=1$ and no peak. For driving rates slower than $\kappa$, random decays of the state vector will be very likely to occur before the state vector can rotate around the y-axis, yielding a non zero projection on the x-axis and a dip in the microwave reflection. A signature of this effect is the splitting of the absorption peak in two for large probe powers. Whilst our signal to noise does not allow the resolution of this feature, it is present in the fitted simulation, supporting this explanation.

At even higher power, the system starts to resonate at a different frequency, corresponding to the junction being replaced by an open circuit when the current traversing the junction exceeds the critical current. This effect is shown in the inset, Fig. S8B.

We use the disappearance of peaks at a high power indicated by the arrow “calibrating power” to acquire a microwave background that is subtracted (divided) in phase (amplitude) to all datasets.

[TABLE]

7 Fitting

Here, we summarize our fitting routine. We start by extracting $\gamma$ from the time-domain data, which will be used in the formula for the linewidth $\kappa_{n}$ in all subsequent fits. By fitting the microwave reflection $S_{11}$ to a sum of Lorentzians (see Eq. (S26)), we get access to the peaks linewidths and amplitudes which allows us to determine $\kappa$, $\kappa_{\text{ext}}$ and $n_{\text{th}}^{(H)}$. By fitting $S_{11}$ to the eigenfrequencies obtained from a diagonalization of the Hamiltonian of Eq. (S2), we determine the values of the circuit elements. The occupation of the low mode $n_{\text{th}}$ is determined separately for each individual experiment. Each step is detailed in the subsections below.

7.1 Low-frequency mode dissipation

We start by fitting the thermalization from the ground-state measured in time-domain (Fig. S13A) to determine $\gamma$. Since the line-width of the $S_{11}(t)$ peaks is a function of $\gamma$ and $n_{\text{th}}$, we start by postulating these two values to extract a first estimate of the time evolution of $P_{n}$. By fitting the evolution of $P_{n}$ to the rate equation of Eq. (6), we extract a new value for $\gamma$ and $n_{\text{th}}$. We then repeat this process many times, each time using the new values $\gamma$ and $n_{\text{th}}$ to fit $S_{11}(t)$, until we converge to $\gamma/2\pi=23.5\cdot 10^{3}s^{-1}$.

The low-frequency mode dissipation can also be measured without recourse to time-domain experiments. The knowledge of the power dependent AC-stark shift and the cooperativity, measured in a single tone cooling experiment, is sufficient to extract $\gamma$. We use this method to confirm our time-domain results, as well as verify the theory developed in Sec.  4.5.2. First we measure the AC-stark shift of the $n=0$ peak, from which we extract the the proportionality factor $\xi^{2}/P$, between pumping rate $\xi$ and pump power $P$ (Fig. S9A). Secondly we determine the power at which the strong coupling regime arises (Fig. S9B). Above this power, the line-width of the $n=1$ peak will rise as the state $\ket{g,1}$ hybridizes with $\ket{f,0}$ under the effect of the cooling pump. Below this power, the line-width of the $n=1$ peak is approximatively constant, and its height provides an accurate measure of $P_{1}$. In this regime, we thirdly extract the ratio of probabilities $P_{2}/P_{1}$ and $P_{1}/P_{0}$ Following Eqs. (S64), the former should remain constant $P_{2}/P_{1}=n_{\text{th}}/(1+n_{\text{th}})$. The latter, however, decreases with power, $P_{1}/P_{0}=(n_{\text{th}}+CP_{f})/(1+n_{\text{th}}+CP_{g})$, and fitting this curve provides the conversion factor between cooperativity $C$ and power. If we also know the anharmonicity $A_{\text{H}}$, cross-Kerr $\chi$, the high-mode occupation $n_{\text{th}}^{(H)}$ and dissipation rate $\kappa$, we can estimate the low-mode dissipation $\gamma=2(\xi^{2}/P)/(C/P)A_{\text{H}}\chi/\kappa/(1+3n_{\text{th}}^{(H)}/2)\simeq 2\pi\cdot 16\cdot 10^{3}s^{-1}$ close to the value obtained in time-domain. The discrepancy is due to the inaccuracy of the relation $P_{1}/P_{0}=(n_{\text{th}}+CP_{f})/(1+n_{\text{th}}+CP_{g})$, arising from the off-resonant driving of other four-wave mixing transitions.

7.2 High-frequency mode dissipation and device temperature

Using $\gamma/2\pi=23\cdot 10^{3}s^{-1}$, we fit the spectra shown in Fig. S10 to fix $\kappa$, $\kappa_{\text{ext}}$ and $n_{\text{th}}^{(H)}$. Here, the fridge temperature is varied, and from a fit of Eq. (S26) we extract $\kappa,\kappa_{\text{ext}}$, $n_{\text{th}}$ and $n_{\text{th}}^{(H)}$ at each temperature. We took care to let the system thermalize for $\sim 10$ minutes at each temperature before starting measurements. The linear scaling of low-mode temperature with fridge temperature, shown in Fig. S10B, confirms that we can extract a realistic mode temperature from the Bose-Einstein distribution. A large difference in temperature is measured between low and high mode, which could be explained by the difference in external coupling to the feedline of the two modes. We fix the values of $\kappa$, $\kappa_{\text{ext}}$ and $n_{\text{th}}^{(H)}$ to the lowest fridge temperature fit (Fig. S10C).

We leave $n_{\text{th}}$ as free parameters in the other experiments as it was found to vary by 10 to 20 percent on a time-scale of hours. In the main text we quote the value of $n_{\text{th}}$ of the lowest point in the temperature sweep ($n_{\text{th}}=1.62$), but in the Fock state stabilization measurement, we measured $n_{\text{th}}=1.40$, in the cooling experiment $n_{\text{th}}=1.81$ and in the time-domain $n_{\text{th}}=1.37$. These fluctuations are much smaller than the uncertainty in fitting $n_{\text{th}}$. In both the cooling and Fock state stabilization experiments, $n_{\text{th}}$ was extracted from an initial measurement of $S_{11}$ in absence of pump tones.

7.3 Circuit parameters

The frequency of the system transitions (and hence the circuit parameters) is determined by fitting a numerical steady-state calculation of $S_{11}$ to the lowest temperature data (Fig. S9C). This simulation, described in 5.2, starts with a diagonalization of the Hamiltonian of Eq. (S2). In this fit we additionally impose that the transition frequency $\ket{g,0}\leftrightarrow\ket{g,1}$ match the value measured in two-tone spectroscopy (Fig. S15A).

We further verify the values of $C_{\text{H}}$, $C_{\text{L}}$, $L_{\text{J}}$ and $L$, as well as the black-box circuit analysis of Sec. 4.1, by extracting $A_{\text{H}}$, $\chi$ and $\omega_{\text{H}}$ for a varying $L_{\text{J}}$. The junction or rather SQUID inductance is modified by sweeping the flux traversing it. This is done by current-biasing a coil situated beneath our sample. We show in Figs. S11B,C,D the result of fitting a sum of Lorentzians to the flux-dependent spectrum (Fig. S11A). For each extracted parameter, we plot the theoretical evolution with flux obtained through a numerical diagonalization of the Hamiltonian of Eq. (S2) (Taylor expanded to the 8-th order), as well as the analytical expressions obtained from black-box quantization (Eqs. (S4,S6,S8,S10)). The only discrepancy is between the numerical and analytical estimation of $A_{\text{L}}$ and $\chi$. It arises due to a term obtained from the quartic non-linearity of the junction proportional to: $(\hat{a}^{\dagger}\hat{a}+1)(\hat{a}^{\dagger}+\hat{a})(\hat{b}^{\dagger}+\hat{b})+h.c.$. This term resembles a beam-splitter interaction which typically makes an oscillator more anharmonic when coupled to an oscillator more non-linear than itself.

The asymmetry of the SQUID dictating the dependence of $L_{\text{J}}$ on flux was a fit parameter in the construction of this figure and was found to be $20\%$. This experiment also suffered from a number of flux jumps, where the transition frequency of the circuit suddenly jumped to a different value. The flux was then swept until we recovered the same frequency before continuing the scan. This data-set is thus assembled from 6 different measurements. Therefore, an entire flux periodicity was not successfully measured, making the conversion between the current fed into a coil under the sample and flux a free parameter.

8 Supplementary experimental data

8.1 Flux dependence of thermal and dissipation parameters

The flux sweep shown in Fig. S11 also gives access to the temperature and dissipation rates of the modes as a function of flux which are shown in Fig. S12. These are extracted from a fit of the sum of Lorentzians (Eq. (S26)), i.e. from the line-widths and amplitudes of the measured peaks. This relies on our estimation of $\gamma$, which is assumed to be a constant, and $n_{\text{th}}^{(H)}$, which is difficult to extract due to the low signal-to-noise ratio as well as the $\ket{e,0}\leftrightarrow\ket{f,0}$ peak crossing the $\ket{g,n}\leftrightarrow\ket{e,n}$ peaks. To investigate the accuracy of these fits, we plot in Fig. S12C the external quality factor of the high mode as a function of its frequency (Fig. S12C) . We find a clear mismatch with the behavior expected from our circuit analysis. This indicates that we cannot confidently state that the temperature of the low mode and dissipation of the high mode fluctuate with flux as shown here, or even provide meaningful error bars. Further analysis could take the form of time-domain measurements at each flux points to determine $\gamma$, or higher signal-to-noise measurements of the $\ket{e,0}\leftrightarrow\ket{f,0}$ to fix $n_{\text{th}}^{(H)}$.

8.3 Four-wave mixing spectrum

By measuring the spectrum whilst sweeping the frequency of a pump tone, we show in Fig. S14 the multitude of four-wave mixing processes possible in this system. Panel A is particularly relevant to the cooling experiment, and is shown in Fig. 3A, as one can see the relevant transitions lying next to the cooling transition $\ket{g,0}\leftrightarrow\ket{f,1}$. We tested different combinations of raising and cooling four-wave mixing processes (panels B and C) for cooling and Fock-state stabilization, but these alternatives consistently produced lower state occupations than the results shown in the main text.

Two transitions in panel A are unexpected from a simple four-wave mixing approach to the system: $\ket{g,n}\leftrightarrow\ket{f,n+3}$ and $\ket{e,0}\leftrightarrow\ket{h,n+3}$. These are six-wave mixing processes and one could expect them to have very weak effects. However in this system the cross-Kerr $\chi$ is a considerable fraction of $\omega_{L}$. The usually neglected term of the quartic non-linearity of the junction proportional to $\chi(2\hat{a}^{\dagger}\hat{a}+1)(\hat{b}\hat{b}+\hat{b}^{\dagger}\hat{b}^{\dagger})$ then leads to the dressed low Fock state $\ket{n}$ having a significant overlap with the bare states $\ket{n\pm 2k}$ where k is a positive integer. The transition $\ket{g,0}\leftrightarrow\ket{f,3}$ is thus visible since $\ket{f,3}$ has a large overlap with $\ket{f,1}$ and $\ket{g,0}\leftrightarrow\ket{f,1}$ is an easily drivable four-wave mixing transition.

8.4 Low-frequency spectrum

We monitor the height of the $\ket{g,0}\leftrightarrow\ket{e,0}$, whilst sweeping the frequency of a secondary pump tone. As shown in Fig. S15A,B, this allows us to easily measure the anharmonicity of the high mode and the frequency of the low mode.

The line-width of the low-mode peak is considerably larger than the previously determined low-mode dissipation rate $\gamma/2\pi=23\cdot 10^{3}s^{-1}$. If the line-width was equal to $\gamma$, we would expect to see photon number splitting, distinct peaks separated by the low-mode anharmonicity $A_{\text{L}}$, corresponding to the transitions $\ket{g,n}\leftrightarrow\ket{g,n+1}$. To understand why this is not the case, we fit a steady-state numerical computation of a pumped and probed Hamiltonian

[TABLE]

with the collapse operators of Eq. (S22). The only free parameter is the pumping strength $\epsilon_{\text{p}}\sim 16\times\gamma$, the probe strength was taken to be negligibly small with respect to all other rates in the model. By varying simulation parameters, we can then explore the origin of this broad line-width. These results are summarized in Fig. S15C. Reducing the pumping strength $\epsilon_{\text{p}}$ will suppress what is usually referred to as ‘power broadening’, at the expense of the signal-to-noise ratio, but does not reveal photon-number splitting. By reducing $\gamma$ to a negligibly small rate, photon number splitting can only be glimpsed behind a line-width broadening induced by the process $\ket{g,n}\rightarrow\ket{e,n}$ which occurs at a rate $\kappa n_{\text{th}}^{(H)}$. This becomes clear if we instead keep $\gamma/2\pi=23\cdot 10^{3}s^{-1}$ and take the limit $n_{\text{th}}^{(H)}=0$, making the first two peaks apparent. As derived in Eq. (S26), the line-width of a thermally populated anharmonic oscillator broadens significantly with its thermal occupation, which is responsible in this case for the disappearance (broadening) of peaks $n\geq 2$. By reducing both $\gamma$ and $n_{\text{th}}^{(H)}$, photon-number resolution would become visible.