Kerr-free three-wave mixing in superconducting quantum circuits

TL;DR

This paper demonstrates a method to achieve Kerr-free three-wave mixing in superconducting circuits using a SNAIL, significantly enhancing the dynamic range of quantum amplifiers by reducing signal-induced detuning.

Contribution

The authors introduce a flux and pump tuning technique to operate a SNAIL-based amplifier at a Kerr-free point, increasing photon processing capacity.

Findings

10-fold increase in photon processing capacity

Successful in situ Kerr-free operation

Method applicable to other quantum circuits

Abstract

Quantum-limited Josephson parametric amplifiers are crucial components in circuit QED readout chains. The dynamic range of state-of-the-art parametric amplifiers is limited by signal-induced Stark shifts that detune the amplifier from its operating point. Using a Superconducting Nonlinear Asymmetric Inductive eLement (SNAIL) as an active component, we show the ability to in situ tune the device flux and pump to a dressed Kerr-free operating point, which provides a 10-fold increase in the number of photons that can be processed by our amplifier, compared to the nominal working point. Our proposed and experimentally verified methodology of Kerr-free three-wave mixing can be extended to improve the dynamic range of other pumped operations in quantum superconducting circuits.

Equations

$$H_{\mathrm{SPA}}/\hslash ={\omega }_a(\Phi )a^{†}a+\sum _{n=3}^{\infty }{g}_n(\Phi )(a^{†}+a{)}^n.$$

$$H_{\mathrm{eff}}/\hslash =-{\Delta }_b{b}^{†}b+g(b^2+b^{†2})+\frac{K}{2}{b}^{†2}{b}^2,$$

$${\Delta }_b(n_p)$$

$$g(n_p)$$

$$K(n_p)$$

$$G=1+\frac{(2|g|\kappa {)}^2}{(\frac{{\kappa }^2}{4}+{\Delta }_i{\Delta }_s-4|g{|}^2{)}^2+\frac{{\kappa }^2}{4}({\Delta }_i-{\Delta }_s{)}^2},$$

$$II{P}_3=\frac{\kappa }{|K|}\frac{1}{(\sqrt{G}+1{)}^3}\hslash {\omega }_a\kappa ,$$

$$G_0=1+\frac{4{\kappa }^2|g{|}^2}{({\Delta }_b^2+\frac{{\kappa }^2}{4}-4|g{|}^2{)}^2}.$$

$$G=1+\frac{4{\kappa }^2|g{|}^2}{({\Delta }_{\mathrm{eff}}^2+\frac{{\kappa }^2}{4}-4|g{|}^2{)}^2},$$

$${\Delta }_{\mathrm{eff}}={\Delta }_b-3K\frac{G{P}_{\mathrm{in}}}{\hslash {\omega }_a\kappa },$$

$$\displaystyle P_{{\rm in}}=\frac{\hbar\omega_{a}\kappa}{3KG}\bigg{(}\Delta_{b}\pm\sqrt{4|g|^{2}-\frac{\kappa^{2}}{4}+\frac{2\kappa|g|}{\sqrt{G}}}\bigg{)}.$$

$$\displaystyle P_{{\rm in}}=\frac{\hbar\omega_{a}\kappa^{2}}{3KG}\bigg{\{}\frac{\Delta_{b}}{\kappa}\pm\sqrt{\frac{\Delta_{b}^{2}}{\kappa^{2}}+\frac{\sqrt{G_{0}}-\sqrt{G}}{\sqrt{G_{0}G}}\sqrt{\frac{\Delta_{b}^{2}}{\kappa^{2}}+\frac{1}{4}}}\bigg{\}}.$$

$$\displaystyle P_{{\rm in}}=\frac{\hbar\omega_{a}\kappa^{2}}{36g_{4}^{*}G}\bigg{\{}\frac{\Delta}{\kappa}\pm\sqrt{\frac{\Delta^{2}}{\kappa^{2}}+\frac{\sqrt{G_{0}}-\sqrt{G}}{\sqrt{G_{0}G}}\sqrt{\frac{\Delta^{2}}{\kappa^{2}}+\frac{1}{4}}}\bigg{\}}.$$

$$a_{-\omega }^{†}=\frac{{\omega }_a-\omega }{{\omega }_a+\omega }{a}_{\omega }.$$

$$\alpha (t)=\sum _x\left({\alpha }_x{e}^{-i{\omega }_{x}t}+\frac{{\omega }_a-{\omega }_x}{{\omega }_a+{\omega }_x}{\alpha }_x^{\ast }{e}^{i{\omega }_{x}t}\right),$$

$$(\omega +{\Delta }_b+i\kappa /2){\alpha }_s$$

$$(-\omega +{\Delta }_b+i\kappa /2){\alpha }_i$$

$$({\Delta }_b+i\kappa /2){\alpha }_h$$

$$\displaystyle\big{(}\Delta_{b}-12g_{4}^{*}|\alpha_{h}|^{2}+i\kappa/2\big{)}\alpha_{h}=4g_{3}\alpha_{p}\alpha_{h}^{*}.$$

$$|{\alpha }_h{|}^2=\frac{1}{12g_4^{\ast }}\left({\Delta }_b\pm \sqrt{(4g_3|{\alpha }_p|{)}^2-\frac{{\kappa }^2}{4}}\right).$$

$${\eta }_p=\frac{G{P}_{1dB}}{P_p}.$$

$$H_{\mathrm{SPA}}/\hslash ={\omega }_a{a}^{†}a+g_3(a+a^{†}{)}^3+g_4(a+a^{†}{)}^4+...$$

$$\mathcal{L}=\int_{-d}^{-0}\bigg{[}\frac{c}{2}(\partial_{t}\phi)^{2}-\frac{1}{2\ell}(\partial_{x}\phi)^{2}\bigg{]}dx+\int_{+0}^{d}\bigg{[}\frac{c}{2}(\partial_{t}\phi)^{2}-\frac{1}{2\ell}(\partial_{x}\phi)^{2}\bigg{]}dx-MU_{S}\bigg{(}\frac{\phi^{R}-\phi^{L}}{M\phi_{0}}\bigg{)},$$

$$U_S(\phi )=-E_J\left[\alpha \cos \phi +3\cos \left(\frac{{\phi }_{\mathrm{ext}}-\phi }{3}\right)\right]$$

$$\delta S_{\mathrm{L}}$$

$$\delta S_{\mathrm{R}}$$

$$\delta S_{\mathrm{SNAIL}}$$

$${\partial }_t^2\varphi -v^2{\partial }_x^2\varphi =0.$$

$${\partial }_x\varphi =0.$$

$$\frac{1}{\phi_{0}}\frac{\partial U_{S}(\varphi)}{\partial\varphi}\bigg{|}_{\varphi=\frac{\phi^{R}-\phi^{L}}{M\phi_{0}}}=\frac{1}{\ell}\partial_{x}\phi^{L/R},$$

$$U_{S}(\varphi)=E_{J}\bigg{(}\frac{c_{2}}{2!}(\varphi-\varphi_{{\rm min}})^{2}+\frac{c_{3}}{3!}(\varphi-\varphi_{{\rm min}})^{3}+\frac{c_{4}}{4!}(\varphi-\varphi_{{\rm min}})^{4}+...\bigg{)},$$

$$\frac{E_{J}}{\phi_{0}}\bigg{[}c_{2}\left(\frac{\phi^{R}-\phi^{L}}{M\phi_{0}}-\varphi_{{\rm min}}\right)+\frac{c_{3}}{2!}\left(\frac{\phi^{R}-\phi^{L}}{M\phi_{0}}-\varphi_{{\rm min}}\right)^{2}+\frac{c_{4}}{3!}\left(\frac{\phi^{R}-\phi^{L}}{M\phi_{0}}-\varphi_{{\rm min}}\right)^{3}\bigg{]}=\frac{1}{\ell}\partial_{x}\phi^{L/R}.$$

$$\varphi (x,t)$$

$${\varphi }^{L/R}$$

Full text

I Introduction

Superconducting circuits offer the attractive possibility to synthesize systems with tailor-made Hamiltonians that display a variety of physical phenomena. Among the vast diversity of the Hamiltonians that can be produced, bilinear multimode bosonic Hamiltonians are the simplest, and yet can be wielded to generate nontrivial operations for quantum information processing. For example, single- and two-mode squeezing Hamiltonians are used to construct quantum-limited parametric amplifiers Clerk et al. (2010); Roy and Devoret (2016). Tunable-strength frequency conversion Hamiltonians enter in quantum state transfer between remotely separated modes Axline et al. (2018); Campagne-Ibarcq et al. (2018); Kurpiers et al. (2018). Bilinear couplings with well-controlled phases are crucial for realizing active nonreciprocity in few-body systems such as parametric circulators and directional amplifiers Sliwa et al. (2015); Lecocq et al. (2017), as well as for simulating many-body physics of topological band structure using photonic systems Fang et al. (2012) and implementing topological traveling-wave amplifiers Peano et al. (2016). Furthermore, bilinear interactions are essential auxiliaries in nonlinear operations, such as the universal two-qubit entangling gate eSWAP Gao et al. (2018, 2019) or Clifford gates on GKP-encoded logical qubits Gottesman et al. (2001); Terhal and Weigand (2016).

For such applications, it is important to master the implementation of bilinear Hamiltonians with high dynamic range. A common recipe for realizing a bilinear Hamiltonian is to pump the four-wave mixing nonlinearity of the Josephson junction with an external microwave drive. While this procedure synthesizes the target Hamiltonian, it unavoidably introduces additional spurious terms of higher order. To illustrate this inherent problem encountered in the task of Hamiltonian design in the realm of superconducting quantum circuits, consider the particular task of realizing the Hamiltonian of a degenerate parametric amplifier $H_{\rm DPA}=-\Delta_{b}b^{\dagger}b+g(b^{2}+b^{\dagger 2})$, where $b$ refers to the photon annihilation operator in the rotating frame. It can be constructed with four-wave mixing by pumping $\overline{n}$ photons at detuning $\Delta$ into the Duffing oscillator realized with, for instance, a capacitively shunted Josephson junction or SQUID Eichler and Wallraff (2014); Roy and Devoret (2016). In such an implementation, the junction four-wave mixing with strength $K$ induces three terms: (i) the desired bilinear squeezing term with $g=K\overline{n}/2$, (ii) the spurious Stark shift $\Delta_{b}-\Delta=-2K\overline{n}$, and (iii) the residual quartic term $\delta H=Kb^{\dagger 2}b^{2}/2$. Kerr terms of the type (iii) have in general serious consequences for the dynamic range of various pumped processes. For example, the photon-number-dependent rotation of the phase space caused by $\delta H$ leads to the distortion of quantum states of light in a microwave cavity at the level of a few photons Vlastakis et al. (2013); Kirchmair et al. (2013), to “bananization” of squeezing Boutin et al. (2017); Malnou et al. (2018), and to saturation of parametric amplifiers Eichler and Wallraff (2014); Liu et al. (2017); Frattini et al. (2018); Planat et al. (2019). In addition, the spurious Stark shift of the type (ii) limits the achievable gain in the multipumping schemes for directionality Sliwa et al. (2015); Lecocq et al. (2017) and deteriorates phase matching in traveling-wave amplifiers (TWPA) O’Brien et al. (2014); Macklin et al. (2015); Vissers et al. (2016); Ranzani et al. (2018).

Is it possible to generate pumped bilinear interactions without the detrimental sideeffects provided by Kerr? A necessary ingredient toward this goal is three-wave mixing, which is possible in superconducting circuits biased with external magnetic flux or DC current. For amplification, the third-order nonlinearity can be directly pumped to generate the squeezing term $g$ without the need for the Kerr term $K$. One example of such a three-wave mixing implementation of the DPA is based on flux-modulating the SQUID Yamamoto et al. (2008); Zhou et al. (2014); Simoen et al. (2015). Alternatively, non-degenerate three-wave mixing is also available via the Josephson ring modulator (JRM) Abdo et al. (2013), and is utilized in the Josephson parametric converter (JPC) Bergeal et al. (2010). In both of these implementations the residual Kerr has been determined to cause amplifier saturation Zhou et al. (2014); Liu et al. (2017). Recently a two-terminal three-wave mixing element, the Superconducting Nonlinear Asymmetric Inductive eLement (SNAIL), has been introduced Frattini et al. (2017, 2018) as a tool to fight Kerr in situ by tuning the magnetic flux. In addition, three-wave mixing with Kerr suppression capability is also possible using RF SQUID Yurke et al. (1989); Zorin (2016) and inductively shunted JRM Chien et al. (2019).

However, as will be explained below, the Kerr constant $K$ can be significantly dressed by the presence of the pump. Therefore, the strategy to suppress Kerr must take into account this dressing. In this work, we demonstrate that it is indeed possible to realize the bilinear Hamiltonian corresponding to the degenerate parametric amplifier with three-wave mixing, while suppressing the effect of Kerr in the presence of the pump. We confirm the practicality of this Kerr cancellation by experimentally achieving an order of magnitude improvement in the saturation power and intermodulation distortion (IMD) properties of the parametric amplifier at a Kerr-free sweet spot, reaching the $1\rm\,dB$ compression power of $-102\rm\,dBm$ on par with the best published results Macklin et al. (2015). We further show that this is possible without sacrificing nearly quantum-limited noise performance. More broadly, these results lead us to argue that the optimal strategy to engineer bilinear Hamiltonians should combine three ingredients: a) presence of a three-wave mixing capability, b) classical pump tone variable in strength and frequency and c) sufficiently versatile nonlinearity that allows for Kerr cancellation.

II SNAIL Parametric Amplifier

The device used in this work is the SNAIL parametric amplifier (SPA) Frattini et al. (2018). Its electrical circuit, shown in Fig. 1(a), consists of an array of $M=20$ SNAILs embedded into a microstrip transmission line resonator. The array is flux-biased, with a magnetic flux $\Phi$ piercing each loop. This knob enables the tuning of the SPA resonant frequency $\omega_{a}$ in the range $6.2-7.2\rm\,GHz$ and gives access to the cubic nonlinearity of the SNAIL’s flux-dependent potential energy. The signal port of the device is strongly coupled to a $50\>\rm\Omega$ environment via a finger-capacitor and sets the energy damping rate of the SPA mode $\kappa/2\pi\in 150-240\rm\,MHz$, while a weakly coupled pump port is used to supply the pump. The Hamiltonian of such system can be approximated as

[TABLE]

The dependence of the most relevant Hamiltonian parameters on $\Phi$ is shown in Fig. 1(b) and is calculated and measured using the methods elaborated in Ref. Frattini et al. (2018) and in Appendices A and F.

When a strong off-resonant pump is applied at frequency $\omega_{p}=2(\omega_{a}+\Delta)$, this leads to the effective Hamiltonian of a driven oscillator with Kerr nonlinearity

[TABLE]

written in terms of the photon annihilation operator $b$ in the frame rotating at $\omega_{p}/2$. The parameters $\Delta_{b}$, $g$ and $K$ of the effective Hamiltonian in Eq. (2) depend on all nonlinearities of the initial Hamiltonian in Eq. (1) and, importantly, also on the pumping condition. The lowest order dependence of these parameters on the average number of pump photons in the resonator $n_{p}$ is given by

[TABLE]

These relations translate into the remarkable capability of the SPA to harness Kerr-free three-wave mixing: the parameters $g$ and $K$ are defined by different order nonlinearities, and their flux dependence is such that $K(\Phi)$ can go through zero and change sign while $g(\Phi)$ remains large enough to enable the operation of the SPA as an amplifier.

The phase-preserving gain of such an amplifier at the signal frequency near $\omega_{p}/2$ in the presence of $n_{s}$ signal and $n_{i}$ idler intraresonator steady-state photons can be calculated in the input-output formalism using a semiclassical harmonic balance approximation as in Refs. Zhou et al. (2014); Sundqvist and Delsing (2014); Frattini et al. (2018) and Appendix D. The result is given by

[TABLE]

where $\Delta_{s/i}(n_{p},n_{s},n_{i})=\Delta_{b}(n_{p})-K(n_{p})\times(n_{s/i}+2n_{i/s})$. To obtain large gain, one needs to tune the model described by Eq. (2) close to the parametric instability threshold at $n_{p}=(\kappa^{2}+4\Delta_{b}^{2})/(8g_{3})^{2}$. However, the spurious additional Stark shift $K(n_{p})\times(n_{s/i}+2n_{i/s})$, created by the signal and idler photons, detunes the system from its operating condition, causing amplifier saturation. In order to improve the saturation power, we need to find an operating point at which $G$, given by Eq. (6), is highly insensitive to $n_{s/i}$. To first order, this sensitivity is determined by $K$ dressed by the pump. As we will show next, in the presence of $10^{3}-10^{4}$ pump photons the $O(n_{p})$ term in Eq. (5) becomes comparable to the $n_{p}$-independent term, and the lowest-order perturbative approximation for $K(n_{p})$ given by Eq. (5) breaks down. This breakdown is particularly important in the region of suppressed Kerr.

III Pump-induced dressing of Kerr

Experimentally, we can acquire insight into the behavior of $K(n_{p})$ at large numbers of steady-state pump photons by measuring the Stark shift and intermodulation distortion. Directly measuring the Stark shift coming from photons in the pump at $\omega_{p}$ is difficult in practice because it is obscured by the amplification process arising from the second term in Eq. (2). We instead measure the Stark shift caused by $\overline{n}$ steady-state photons of a strong near-resonant drive placed at $\omega_{d}=2\pi\times 7.8\rm\,GHz$, about seven linewidths away from the resonance, which emulates the effect of the pump. The resulting $\Delta_{\rm Stark}(\overline{n})$ is shown in Fig. 2(a) for different fluxes in the range $(0.24-0.44)\Phi_{0}$, where $\Phi_{0}=h/2e$ is the superconducting magnetic flux quantum. To the lowest order in $\overline{n}$, this Stark shift is given by $24g_{4}^{*}\overline{n}$. The upper bound on $\overline{n}$ in Fig. 2(a) is roughly determined by the onset of chaotic behavior, caused by the excitation of free-particle-like states Verney et al. (2019). The inset in Fig. 2(a) shows that in a certain flux range we are able to suppress the linear contribution to the Stark shift and observe the beginning of an oscillation. We emphasize that this oscillation, similar to that theoretically predicted for a Josephson junction Kochetov and Fedorov (2015); Verney et al. (2019), occurs due to the higher-order terms in the SPA Hamiltonian given by Eq. (1). This phenomenon corresponds to the dressing of the Kerr constant, which becomes especially relevant in the flux region of suppressed Kerr.

We further confirm the pump-induced dressing of $K$ by investigating the IMD properties of the SPA. In this experiment, two near-resonant signals at $\omega_{s1}$ and $\omega_{s2}$ are sent to the resonator, and the sidebands at $2\omega_{s1}-\omega_{s2}$ and $2\omega_{s2}-\omega_{s1}$ appear at the output due to four-wave mixing. The input-referred third-order intercept point $IIP_{3}$ is then calculated as the extrapolated signal power for which the power in the sidebands is equal to the signal power. Here, it is related to the Kerr constant via

[TABLE]

as proved in Frattini et al. (2018). The measurement result of $IIP_{3}$ with the pump off is shown in Fig. 2(b) in blue, and the corresponding $g_{4}^{*}$ extracted according to Eqs. (7) and (5) is shown in Fig. 1(b) in good agreement with that extracted from the linear Stark shift. Eq. (7) suggests that $IIP_{3}$ decreases as $G^{-3/2}$ for large $G$ under the assumption that $K$ is independent of $n_{p}$ as in Eq. (5), and, thus, the peak shape in $IIP_{3}$ at $\Phi=0.40\,\Phi_{0}$ would remain unchanged with increasing gain. However, in the experiment we find that pumping at $\omega_{p}=2\omega_{a}$ (i.e. $\Delta=0$) leads to the complete disappearance of the peak in $IIP_{3}$ at gains $G>5\rm\,dB$, meaning that $K(n_{p}(G))$ is no longer small. This dressing is most pronounced in the region of reduced Kerr where the contribution given by Eq. (5) is intentionally suppressed. Outside this region, the approximation $K(n_{p})=\rm const$ works well up to $n_{p}$ required for $G=20\rm\,dB$, see Appendix C.

We next formulate a procedure to correct for the dressing of $K$ by the pump. As suggested by the measurement of nonlinear Stark shift shown in Fig. 2(a), changing $n_{p}$ at a fixed $\Phi$ can influence the sensitivity of the Kerr-induced terms $\Delta_{s/i}(n_{p},n_{s},n_{i})$ in Eq. (6) to $n_{s/i}$, providing us the tool that we need in order to find the effective Kerr-free point. At a fixed $\Phi$, we can vary the pump photon number $n_{p}$ required to reach $G=20\rm\,dB$ by changing the pump detuning $\Delta$, which had been set to zero in the analysis so far.

IV Stability regions

The model described by Eq. (2) exhibits a variety of new effects when $\Delta$ is varied away from zero. We first discuss these effects, which can be understood by considering the stability diagram of the Kerr oscillator driven parametrically near twice its resonant frequency, shown in Fig. 3(a) for $\Phi=0.30\Phi_{0}$. The colors in this diagram represent the three regions in the space of parameters $\Delta_{b}$ and $g$, in which the classical nonlinear dissipative dynamics of a complex amplitude $b$ has qualitatively different phase portraits Dykman et al. (1998); Wustmann and Shumeiko (2013). In general, the motion of $b$ is not governed by a Hamiltonian that has separate potential and kinetic energy terms. However, under certain conditions, due to the separation of time scales one can integrate out the fast component of the motion and reduce the problem to the slow motion in an effective 1D potential Zorin and Makhlin (2011); Lin et al. (2014). This potential can have one, two or three minima depending on the parameters $\Delta_{b}$ and $g$, and is sketched for illustrative purposes in Fig. 3(a). Such a 1D potential helps to qualitatively think about the problem, but the partition of the stability diagram into three separate single-, bi- and tristable regions can be done regardless, because it is based on the full 2D phase portrait of the system.

Region (I) of the stability diagram has a trivial ground state, in which $b$ fluctuates near zero. At $G=20\rm\,dB$ operating points in this region, the device behaves as a nearly quantum-limited amplifier which is seen from the measurement of the noise visibility ratio (NVR), defined as the ratio of noise power spectral densities with the SPA pump on and off, see Fig. 3(b) right upper panel and Appendix A. In addition, the saturation curve taken at the signal frequency near $\omega_{p}/2$ monotonically decreases with increasing signal power, see Fig. 3(b) right lower panel.

Region (II) has two stable points and, in the resulting ground state, $b$ has a finite mean amplitude, implying the resonator oscillates at frequency $\omega_{p}/2$ even with no applied signal power, an effect termed as ”parametric self-oscillation” or ”period-doubling” Dykman (2012). This region is separated from the rest of the diagram by the parametric instability threshold, where the gain of a parametric amplifier diverges Wustmann and Shumeiko (2013). In the $K<0$ case, due to the Stark shift produced by $n_{p}$, there is a maximal positive detuning $\Delta$ beyond which large gain can no longer be obtained. Immediately next to the threshold, the black line denotes the pump photon number and detuning for $G=20\rm\,dB$. In our experiment, we always operate on this line below the threshold and thus do not enter region (II), although note that it is an interesting region for operating the device as a parametric oscillator Krantz et al. (2016); Lin et al. (2014).

On the negative detuning side, the $G=20\rm\,dB$ line enters the region (III) of the oscillator stability diagram, which has three stable points. As long as the device operates at the single, global minimum at zero amplitude, it remains nearly quantum-limited as seen by the NVR trace in Fig. 3(b) upper middle panel. Nevertheless, as shown in Fig. 3(b) lower middle panel, the response of the amplifier to large signals changes qualitatively. When the signal power is increased above a certain threshold, the system undergoes the equivalent of a first-order phase transition and adopts a new working point with increased gain. This phenomenon is captured by the self-consistent equation (6) for gain $G$ as a function of $n_{s/i}(G)$, as shown in Appendix C. The characteristic shape of the saturation curve in this regime was called “shark fin” in Ref. Liu et al. (2017), where a similar phenomenon was reported for the JPC. Interestingly, shark fins can be extremely sensitive to the level of input signal power, changing the output by $15\rm\,dB$ when the input difference is less than $1\rm\,dB$.

The presence of additional stable points in the region (III) implies that the period-doubling effect is also possible here. The associated instability where the system dynamically switches between the global minimum and the secondary ones brings the danger of increased noise temperature of the amplifier. However, the high-amplitude minima are separated from the trivial minimum by potential barriers that prevent the system from switching Dykman et al. (1998); Svensson et al. (2017). Therefore, the NVR is not appreciably degraded at small red detuning inside region (III). Following the $20\rm\,dB$ curve in Fig. 3(a), it is only beyond the detuning where all three coexisting steady states have equal populations that the system prefers to stay in the high-amplitude state Lin et al. (2015). Hence, we attribute the significant increase of NVR at large negative detuning, shown in Fig. 3(b) upper left panel, to the enhanced switching between the high- and low-amplitude stable points. This effect is quantified by the increase of the noise temperature of the SPA by a factor of $2.9$ at $\Delta=-500\rm\,MHz$. In such a state, the amplifier is no longer quantum-limited and cannot be described by the conventional theory that leads to Eq. (6).

The stability diagram in Fig. 3(a) was constructed at $\Phi=0.30\,\Phi_{0}$, away from the region of reduced Kerr. Since we cannot predict the stability regions at the fluxes where we expect the approximation in Eq. (5) to fail, we need to perform an empirical search for the dressed Kerr-free point in the space of parameters $\Phi$ and $\Delta$.

V Dressed Kerr-free point

We show in Fig. 4 the results of the measurement of gain saturation (upper panels) and NVR (lower panels) for high-resolution sweeps of detuning $\Delta$ done at five different fluxes $\Phi$ in the region of reduced Kerr. At every detuning, the pump power is adjusted to obtain the small-signal gain of $20\rm\,dB$. The input signal power is then swept beyond the $1\rm\,dB$ compression point (solid black). The left and right boundaries of the $\Delta$ sweep at each $\Phi$ correspond to detunings where a small-signal gain $G=20\rm\,dB$ could not be achieved at any applied pump power. Moving from smaller to larger $\Phi$, the $\Delta$ tuning range decreases from about $900\rm\,MHz$ at $\Phi=0.19\,\Phi_{0}$ to about $200\rm\,MHz$ at $\Phi=0.48\,\Phi_{0}$ (not shown) and its median shifts from positive to negative values.

We find the Kerr-free operating point in the presence of the pump to be at $\Phi=0.34\,\Phi_{0},\,\Delta=-95\rm\,MHz$. There, the 1 dB compression power peaks at a value $P_{\rm 1dB}=-102\rm\,dBm$ which is an order of magnitude larger than at $\Delta=0$. The tunable bandwidth over which we observe the $P_{1\rm dB}$ increase is about $200\;\rm MHz$. We measure a similar improvement in the IMD properties with peak $IIP_{3}=-94\rm\,dBm$, see Appendix B. Importantly, as indicated by the NVR, the amplifier remains nearly quantum-limited at this operating point. This result demonstrates that it is indeed possible to realize the degenerate parametric amplifier Hamiltonian with three-wave mixing, while suppressing the effect of Kerr in the presence of the pump.

VI Conclusion

To conclude, through theory and experiments we have formulated the methodology for realizing a quantum-limited degenerate parametric amplifier with suppressed spurious Kerr effect, and hence improved dynamic range. This work raises several directions for future investigation. Our device achieved $P_{\rm 1dB}=-102\rm\,dB$, on par with the best published results, but with room for further improvement. In Appendix E we provide an outlook for how such an optimization would proceed and introduce a new concept of “pump power efficiency” that would need to be addressed to achieve further gains in dynamic range.

Another important direction would be to investigate the limits on the tunable “Kerr-free bandwidth” and on the compression power at the dressed Kerr-free point, which are possibly due to even higher-order terms not included in the effective Hamiltonian in Eq. (2). This naturally raises the question of whether it is possible to mitigate the detrimental effect of these higher order terms by circuit design and choice of pumping condition.

More broadly, we expect that the strategy of operating at a Kerr-free point while accounting for the dressing by the pump will be crucial for improved performance of other flux-biased three-wave mixing amplifiers that are, in principle, capable of operating at a Kerr-free point, such as the JPC with an inductively shunted JRM Chien et al. (2019) or the TWPA with array of RF SQUIDs Zorin (2016) where Kerr cancellation will improve phase matching Zhang et al. (2017). Finally, our methodology for achieving Kerr-free three-wave mixing can be applied for engineering other pumped bilinear Hamiltonians with superconducting circuits while suppressing the effect of spurious processes.

We acknowledge helpful discussions with S. Mundhada, P. D. Kurilovich, V. D. Kurilovich, G. Liu and T. Roy. We also acknowledge the Yale Quantum Institute. Facilities use was supported by the Yale SEAS clean room and YINQE. This research was supported by AFOSR under Grant No. FA9550-15-1-0029, and by ARO under Grants No. W911NF- 18-1-0212, W911NF-18-1-0020 and W911NF-16-1-0349.

Appendix A Device parameters and measurement setup

The SPA device package is shown in Fig. 5. Josephson junctions of the SNAIL are formed by $\rm Al/AlO_{x}/Al$ layers deposited using the Dolan bridge shadow evaporation process (critical current density $j_{c}=120\,{\rm A}/{\rm cm}^{2}$), and have critical currents $I_{c}=8.5\;\rm uA$ and $I_{c}=0.85\;\rm uA$. The microstrip transmission line resonator is formed by a $2\,\mu m$ thick silver film deposited on the back of a $300\,\rm\mu m$ thick silicon wafer to act as a ground plane, and by the aluminum traces on the chip. The silver back plane of the chip is glued using conducting silver paste to the copper back plane of the ${\rm TMM}10{\rm i}$ printed circuit board (PCB), which is soldered to the gold-plated aluminum box, see Fig. 5(a). The first resonant mode of the box is at $18\,{\rm GHz}$, well above relevant frequencies for the operation of the SPA. The pump and signal aluminum transmission line traces on the chip are wire bonded to the copper transmission line traces on the PCB, which are soldered to the edge-mount SMA connectors.

A superconducting NbTi coil that applies magnetic flux to the SNAILs is mounted under the PCB inside of the box, as shown in Fig. 5(a). The aluminum box acts as a partial shield against stray external magnetic fields. For better protection, an additional cryogenic $\mu$-metal shield is fitted around the aluminum box. The device is thermally anchored to the base stage of a dilution refrigerator ($T=24\rm\,mK$) through this shield. It is mounted back-to-back with the $4-8{\rm\,GHz}$ cryogenic circulator, which is connected to the signal port of the SPA with a short 2 inch cable. This configuration allowed us to significantly suppress the small-scale ripples in the data resulting from the impedance mismatches Mutus et al. (2014); Roy et al. (2015); Frattini et al. (2018). We found that the close proximity of the commercial ferrite circulator does not influence the SPA performance.

The SPA pump tone is generated by an Agilent PSG E8257D generator. The DC current is supplied to the coil by a Yokogawa GS200 current source. All measurements, including the IMD and NVR, are done using the corresponding measurement classes of a Keysight PNA-X N5242A network analyzer.

We found that PNA-X “Noise Figure” measurement class sometimes displays a sharp spike at $\omega_{p}/2$, which is absent in a separate measurement done with an Agilent EXA N9010A Spectrum Analyzer, and likely results from the lower tolerance of the PNA-X to the pump leaking through the isolators of the output chain at frequencies close to their cutoff at $12\,\rm GHz$. This spike does not affect the rest of the data, as verified by the independent measurements of NVR, and is filtered out in software.

Extraction of the amplifier noise temperature from NVR requires the calibration of the gain of the output chain, for example with the shot noise thermometer Spietz et al. (2003, 2010), which was not available in our setup. However, our claim that the system is nearly quantum-limited when ${\rm NVR}\approx 8.3\rm\,dB$ is supported by the independent measurement of quantum efficiency $\eta=0.6$ in the experiment Touzard et al. (2019) done with the SPA.

Appendix B Intermodulation distortion (IMD)

In the IMD measurement the two signals at $\omega_{s1}$ and $\omega_{s2}$ of equal power are sent to the input of the SPA. The median of the signals is detuned from the center of the Lorenzian gain at $\omega_{p}/2$ by amount $\delta_{1}\equiv(\omega_{s1}+\omega_{s2})/2-\omega_{p}/2=2\pi\times 500\rm\,kHz$ to avoid phase-sensitive amplification. The detuning is chosen small compared to the amplifier 3dB-bandwidth $B\sim 2\pi\times 25\rm\,MHz$. The separation of the two tones is $\delta_{2}=\omega_{s2}-\omega_{s1}=2\pi\times 100\rm\,kHz$. Due to the nonlinear intermodulation distortion, the sidebands at $2\omega_{s2}-\omega_{s1}$ and $2\omega_{s1}-\omega_{s2}$ appear at the output and are measured using the IMD measurement class of the PNA-X. The $IIP_{3}$ extracted from the raw data is shown in Fig. 6, and exhibits a peak value $IIP_{3}=-94\rm\,dB$ at the dressed Kerr-free point similar to the peak in $P_{1\rm dB}$.

Appendix C Saturation power and shark fins

Saturation of a good amplifier happens at large average signal and idler photon numbers $n_{s/i}\gg 1/2$, justifying the use of the semiclassical approximation, in which the phase preserving gain $G$ at the signal frequency near $\omega_{p}/2$ is given by Eq. (6). The small-signal limit of this expression is given by

[TABLE]

The goal of this section is to derive a closed-form expression for the input signal power $P_{\rm in}$ at which the gain changes from $G_{0}$ to $G$ (which directly gives the closed-form expression for $P_{1\rm\,dB}$) and to study the phenomenon of shark fins. Using input-output theory, we can establish the relation between the input power $P_{{\rm in}}$ and the intraresonator populations $n_{s}$ and $n_{i}$. For large gain $G\gg 1$ we obtain $n_{s}\approx n_{i}\approx GP_{{\rm in}}/\hbar\omega_{a}\kappa$ with a relative error $\delta n_{s}/n_{s}\sim 1/\sqrt{G}$, which is about $10\%$ for $G=20\,\rm dB$. In this approximation, we can rewrite Eq. (6) in the form

[TABLE]

where $\Delta_{\rm eff}$ contains the additional Stark shift created by the signal and idler photons (this formula is given in Frattini et al. (2018), but with a typo which is corrected here). After plugging this into Eq. (10) and solving for $P_{{\rm in}}$, we obtain

[TABLE]

The Hamiltonian parameter $g$ can be expressed through $G_{0}$ from Eq. (8), which leads to

[TABLE]

In this expression both $K$ and $\Delta_{b}$ depend on $n_{p}$, and thus implicitly on gain. When we consider a simplified case of small negative $K$, such that $\Delta_{b}\approx\Delta$ and $K\approx 12g_{4}^{*}$, we can further simplify Eq. (12) to

[TABLE]

At positive detuning, only one solution is possible for $G<G_{0}$ and no solutions for $G>G_{0}$, while at negative detuning two solutions are possible for $G>G_{0}$ and one for $G<G_{0}$, which corresponds to the “shark fin” shape of the saturation curve.

The comparison of Eq. (13) to the measured $P_{1{\rm dB}}$ at a fixed detuning $\Delta=0$ as a function of flux $\Phi$ is shown in Fig. 7(a), and comparison to the measured $P_{1{\rm dB}}$ at a fixed flux $\Phi=0.30\,\Phi_{0}$ as a function of detuning $\Delta$ is shown in Fig. 7(b). For the plotted theory lines, we have assumed that there exists a systematic miscalibration of the resonant frequency (and therefore the detuning $\Delta$) by approximately $25\rm\,MHz$. The compression power is very sensitive to such miscalibration especially near the detuning $\Delta=0$, as evident from Fig. 7(b). Accounting for this shift aligns the theory and data much better in both $\Delta$ and $\Phi$ axes, while without this correction we find a $6\rm\,dB$ systematic disagreement between theory and data in Fig. 7(a).

The extrapolation of Eq. (13) to detunings $\Delta<-300\rm\,MHz$ predicts persisting shark fins, which are not observed in data. We understand this effect to arise from the fact that the ground state of the metapotential in Eq. (2) is significantly modified at large negative detuning $\Delta$ due to the secondary high-amplitude minima, which invalidates the derivation of $P_{1\rm dB}$ presented in this section. Empirically, we locate this region by the increased NVR at large negative detuning, as discussed in the main text. Similarly, the extrapolation of Eq. (13) to fluxes $\Phi\gtrsim~{}0.30\Phi_{0}$ leads to the incorrect prediction of improved compression power at an unpumped Kerr-free point at $\Phi=0.40\,\Phi_{0}$, as shown in Fig. 7(a). Using Eq. (13) with $K$ extracted from the IMD measurement at gain $G=20\rm\,dB$ leads to a much better functional agreement, indicating that the IMD experiment better captures the dressing of the Kerr constant by the pump.

We can view Eq. (10) as a self-consistent equation of state, analogous to the Van der Waals equation for nonideal gas, which exhibits a first-order phase transition Landau and Lifshitz (1969). In this analogy, the detuning $\Delta$ plays the role of temperature, input signal power $P_{\rm in}$ plays the role of pressure, and gain $G$ plays the role of volume. As shown in Fig. 7(d), for the negative detuning $\Delta$, the equation (10) predicts the possible existence of a metastable state with high gain. In the sweep of input signal power, the system switches to this state along the line analogous to the isobar in the Maxwell construction, as shown in the experimental data in Fig. 7(c).

Finally, note that the flux tuning range of the SPA shown in Fig. 7(a) corresponds to a wide tunable bandwidth of about $1\rm\,GHz$. Moreover, at a fixed $\Phi$ the tunability can be accomplished with the help of pump detuning $\Delta$, which corresponds to a tuning range of approximately $700\rm\,MHz$, as shown in Fig.7(b) for $\Phi=0.30\,\Phi_{0}$. However, not all of this range corresponds to the quantum-limited performance as discussed in the main text, reducing the usable tunability at a fixed flux to about $400\rm\,MHz$. Due to this property, such an SPA is a good candidate for time-multiplexed readout schemes, in which the pump frequency can switch to match the required readout channel, similarly to the proposal of flux switching in Ref. Abdo et al. (2017).

Appendix D Stability diagram of the parametrically driven Kerr oscillator

In this section we derive the separation lines between the regions of the stability diagram at $\Phi=0.30\,\Phi_{0}$ shown in Fig. 3(a), using the semiclassical harmonic balance method. When the pump and signal tones are incident on the SPA nonlinear resonator, the mixing between them results in the amplification of the signal, creation of the idler and creation of other intermodulation products at frequencies that can be parametrized as $\omega_{nm}=n\omega_{p}+m\omega_{s}$ where $n$ and $m$ are integers. However, under certain conditions, generation of radiation which does not belong to this family of intermodulation products is possible. Of particular interest to us is the coherent period-doubling tone at frequency $\omega_{p}/2$ generated solely by the pump.

To grasp the major processes that happen in the pumped SPA, we perform a second-order harmonic balance calculation which accounts for the following harmonics:

• $\circ$

Pump $\omega_{p}$, signal $\omega_{s}$ and idler $\omega_{i}=\omega_{p}-\omega_{s}$ comprise a minimal set required for quantum-limited amplification.

• $\circ$

$2\omega_{p}$, $2\omega_{s}$, $2\omega_{i}$, $\omega_{p}+\omega_{s}$, $\omega_{p}+\omega_{i}$, $\omega_{s}-\omega_{i}$ and $\rm DC$ are required in the consistent calculation that takes into account linear Stark shift. These harmonics lead to $\sim g_{3}^{2}/\omega_{a}$ corrections to $K$ in Eq. (5) and to pump-induced Stark shift in Eq. (3).

• $\circ$

$\omega_{p}/2$, $\omega_{p}/2-\omega_{s}$ and $\omega_{p}/2-\omega_{i}$ account for a possible period-doubling effect.

Following the steps of Frattini et al. (2018) and Refs. Zhou et al. (2014); Sundqvist and Delsing (2014), we derive the self-consistent system of equations that links the amplitudes of all chosen harmonics. This is done by solving the quantum Langevin equation (QLE) Gardiner and Zoller (2004) without RWA approximation (because the out-of-band harmonics do not satisfy it). In order to do this, we need to take into account that the negative frequency Fourier component of the annihilation operator $a$ is linked to the positive frequency component via the relation

[TABLE]

The $a_{-\omega}^{\dagger}$ is usually dropped from the frequency-domain QLE for $a_{\omega}$ if all the signals arriving at the resonator are near the resonance. Eq. (14) can be derived from the relation $\dot{\Phi}=C_{a}Q$ between the canonical coordinate $\Phi$ of the mode and its canonical momentum $Q$, written in second-quantized form.

Taking into account Eq. (14), we seek the semiclassical harmonic balance solution to QLE in the form

[TABLE]

where $x$ runs over all harmonics described previously.

Equations for the out-of-band harmonics can be partially solved, and the system further reduced to three complex equations for the amplitudes at $\omega_{s}$, $\omega_{i}$ and $\omega_{p}/2$, denoted as $\alpha_{s}$, $\alpha_{i}$ and $\alpha_{h}$, respectively, as follows:

[TABLE]

where $\omega=\omega_{s}-\omega_{p}/2$ is the signal detuning, $\Delta_{b}$ is defined in Eq. (3), and $u_{s}$, $u_{i}$ and $u_{h}$ denote the drive strengths at corresponding frequencies, following the notations introduced in Ref. Frattini et al. (2018). Note that the pump amplitude $\alpha_{p}=u_{p}/\omega_{a}$ here does not include any corrections due to $\alpha_{s}$, $\alpha_{i}$ and $\alpha_{h}$ – they have been reabsorbed into the above equations.

To derive the boundaries in the stability diagram, we solve these equations in the small-signal approximation, in which the Stark shift due to $|\alpha_{s}|^{2}$ and $|\alpha_{i}|^{2}$ can be neglected. Since no input is sent at $\omega_{p}/2$, we can set the drive strength $u_{h}$ to zero. Then equation (18) reduces to

[TABLE]

After multiplying it by the complex conjugate, we find a trivial solution $\alpha_{h}=0$ and possible high-amplitude solutions

[TABLE]

First, consider $\Delta_{b}>0$ half-plane. In this case only the solution with the minus sign in Eq. (20) is possible if the additional condition $(4g_{3}|\alpha_{p}|)^{2}-\frac{\kappa^{2}}{4}>\Delta_{b}^{2}$ is satisfied. This defines the separation line between regions (I) and (II) of the stability diagram. In region (I) only $\alpha_{h}=0$ solutions exists and is stable. High-amplitude solution exists in region (II) and is stable, while $\alpha_{h}=0$ is unstable.

Next, consider $\Delta_{b}<0$ half-plane. The requirement for $|\alpha_{h}|^{2}$ being real leads to the condition $(4g_{3}|\alpha_{p}|)^{2}>~{}\frac{\kappa^{2}}{4}$. This condition defines the separation line between region (I) and (III) of the stability diagram. Below the separation line only the trivial solution is possible, while above this line both signs in Eq. (20) are possible. This yields five stationary points: a stable trivial $\alpha_{h}=0$ point, a pair of symmetric unstable points defined by the plus sign in Eq. (20) and a pair of symmetric stable points defined by the minus sign. The two unstable stationary points merge with the stable $\alpha_{h}=0$ point on the line which separates regions (II) and (III), defined by the equation $(4g_{3}|\alpha_{p}|)^{2}=\frac{\kappa^{2}}{4}+\Delta_{b}^{2}$. Above this line in region (II) the $\alpha_{h}=0$ point becomes unstable and the two symmetric high-amplitude states remain stable.

The stability diagram plotted in coordinates $\Delta$ and $n_{p}$ (which directly translate into the experimental control knobs pump frequency $\omega_{p}$ and power $P_{p}$) is shown in Fig. 3(a). Note that similar stability diagrams can be found in for example Refs. Zorin and Makhlin (2011); Wustmann and Shumeiko (2013).

Finally we note that by setting $\alpha_{h}=0$ (where appropriate, see Appendix C) the system of equations (16)-(17) together with the input-output relation can be solved to yield the expression (6) for gain in the presence of a finite input signal.

Appendix E Further optimization of the SPA: an outlook

Given the findings of this work, one can raise the question of how to improve the compression power of resonant parametric amplifiers beyond what has been achieved at a Kerr-free point in our particular device. Although sweet spots are in general useful, it is desirable to have a systematic approach to the problem of saturation of quantum-limited amplifiers. We would like to emphasize that such an approach was suggested in Ref. Frattini et al. (2018), where steady improvement of $P_{\rm 1dB}$ of an SPA was achieved by varying the design parameters in the direction corresponding to the increase of the coupling $\kappa$ and the decrease of the overall level of nonlinearity.

Continuing to improve $P_{1\rm dB}$ in this fashion will require more pump power to be delivered to the device at the base of the dilution refrigerator. As such, it is important to compare state-of-the-art quantum-limited parametric amplifiers on the metric of power efficiency, here defined as the ratio of the signal output power at the 1 dB compression point to the required pump power $P_{p}$

[TABLE]

As shown in Fig. 8, the amplifiers with off-resonant pumps, such as three-wave mixing SPAs and JPCs and flux-pumped JPAs, suffer from poor power efficiency, compared to the four-wave mixing amplifiers with near-resonant pumps and traveling-wave amplifiers.

This is an important challenge on the way towards larger dynamic range. Indeed, as the required pump power increases, undesired effects such as increased noise temperature or fridge heating, may result. An Oxford Triton dilution refrigerator will heat up by about $10\;\rm mK$ for $-20\rm\;dBm$ of applied power at the base stage. Typically, an attenuator is used to thermalize the pump to this stage. Assuming $20\;\rm dB$ of attenuation, this puts the limit on allowable pump power at the plane of the device close to $-40\rm\;dBm$. This limit is fuzzy (shown as color gradient in Fig. 8), because it depends on the exact details of the attenuation and filtering chosen for the pump line. Nevertheless, it is clear that further increase of $P_{1\rm dB}$ at the expense of pump power is not a viable approach for optimization.

Instead, this issue needs to be addressed with better design of the pump coupling, which will allow to obtain the required $n_{p}$ with less pump power $P_{p}$ supplied to the plane of the device. Note that the large spectral separation of signal and pump in three-wave mixing amplifiers can be utilized to separately engineer the environmental properties seen by them. For example, the pump could be coupled through a separate resonant filter mode, which would also protect against the leakage of the signal through this port. At the same time, the signal port of the SPA can be optimized for bandwidth using the methods employed in Refs. Mutus et al. (2014); Roy et al. (2015); Naaman et al. (2017). These advances are left for future work.

Appendix F Derivation of the SPA Hamiltonian

The truncated Hamiltonian of the SPA, limited to only the fundamental mode $a$, further submitted to pumping, is given by

[TABLE]

The goal of this section is to obtain the dependence of the parameters $\omega_{a}$, $g_{3}$ and $g_{4}^{*}\equiv g_{4}-5g_{3}^{2}/\omega_{a}$ on the design parameters of the circuit. The expressions (43), (55) and (57) obtained here were used for fitting and making the design choices in this work and in Ref. Frattini et al. (2018).

Equations of motion in the distributed-element model of the SPA

Let the length of each arm of the transmission line resonator be $d$, the capacitance to ground per unit length be $c$, and the inductance per unit length be $\ell$. We treat the array of $M$ SNAILs as a lumped circuit element located at $x=0$ with the phase distributed equally among all SNAILs in the array. We will denote the generalized flux field in the transmission line as $\phi(x)$. For simplicity, we also introduce the notation $\phi^{R/L}=\underset{x\to\pm 0}{\lim}\phi(x)$.

The Lagrangian can be written as

[TABLE]

where $\phi_{0}=\hbar/2e$ is the reduced flux quantum and $U_{S}(\varphi)$ is the potential energy of the SNAIL with a phase drop $\varphi$,

[TABLE]

where $E_{J}$ is the Josephson energy of the large junctions, $\alpha$ is a ratio of the junction inductances in the two arms of the SNAIL, and $\varphi_{\rm ext}=2\pi\Phi/\Phi_{0}$ is the phase bias corresponding to the external magnetic flux $\Phi$, see Frattini et al. (2018).

The action for this Lagrangian is simply given by $S=\int_{t_{1}}^{t_{2}}\mathcal{L}(t)dt$, and can be split into three separate parts $S=S_{{\rm L}}+S_{{\rm R}}+S_{{\rm SNAIL}}$ corresponding to the separate contributions in the Lagrangian in Eq. 23. The variation of action $\delta S$ contains three terms

[TABLE]

According to the variational principle, by requiring $\delta S=0$ we obtain the equations of motion for the flux $\phi(x,t)$ and the boundary conditions.

• $\circ$

In the bulk of the transmission line resonator, for $x\in(-d,0)$ and $x\in(0,d)$, the flux field $\phi(x)$ obeys the wave equation

[TABLE]

where $v=1/\sqrt{\ell c}$ is the phase velocity.

• $\circ$

At the boundaries $x=-d$ and $x=d$ the flux field has zero current boundary condition

[TABLE]

• $\circ$

In the center of the resonator, where it is interrupted by the SNAIL array, the left and right continuity conditions should hold

[TABLE]

which also implies the current conservation $\partial_{x}\phi^{L}=\partial_{x}\phi^{R}$ across the array. Strictly speaking, the bare SNAIL potential $U_{S}(\varphi)$ is $6\pi$-periodic in $\varphi$. However, we are interested in a time scale on which the phase slips can be neglected, and, thus, the fact that $\varphi$ is compact does not matter. In this approximation, we can assume that the phase $\varphi$ is localized near the bottom of one of the many equivalent potential minima and expand the potential around that point (denoted as $\varphi_{\rm min}$)

[TABLE]

where all dimensionless expansion coefficients $c_{k}\equiv c_{k}(\Phi,\alpha)$ depend only on the external magnetic flux $\Phi$ and junction inductance ratio $\alpha$. Using this expansion truncated to the fourth order, we can rewrite the continuity condition at $x=\pm 0$ as

[TABLE]

This condition will allow us to match the solutions of the wave equation (28) in the left and right arms of the microstrip transmission line.

Instead of working with the continuous flux field $\phi(x)$, it is convenient to decompose it using the eigenmode excitations of the circuit. This is done by using as a basis the general solution to Eq. (28) in the left and right arms of the transmission line resonator

[TABLE]

where the coefficients $a_{n}^{L/R}(t)$ and $b_{n}^{L/R}(t)$ should satisfy $[\partial_{t}^{2}+\omega_{n}^{2}](\cdot)=0$. We need to rewrite the boundary condition in Eq. (29) and the continuity condition in Eq. (32) using the eigenmode decomposition. Starting with the boundary condition, by differentiating Eq. (33) we obtain

[TABLE]

and since this has to be satisfied at all times, applying harmonic balance leads to

[TABLE]

and also $b_{0}^{L/R}=0$, meaning that there is no uniform static flux gradient in the resonator. Now we can similarly treat the continuity condition in Eq. (32) using Eq. (34) and its derivative

[TABLE]

The continuity condition contains two equations for the left and right side of the array. Equivalently, we can use one of them and the current conservation condition $\partial_{x}\phi^{L}=\partial_{x}\phi^{R}$ instead. With the help of Eq. (37) and applying harmonic balance, this reduces to

[TABLE]

As the second continuity equation we choose Eq. (32) on the right side of the array, which with the help of Eqs. (33) and (37) can be transformed to

[TABLE]

Note that this equation is nonlinear, which leads to the coupling of various eigenmodes by the SNAIL array mixer. In the limit of weak excitations, we can linearize the system and decompose it into a collection of harmonic oscillators with eigenfrequencies $\omega_{n}$. The lowest such eigenmode of the device is used to implement the SPA.

Eigenmode decomposition

The linear part of Eq. (39) is simply

[TABLE]

where we have introduced the inductance of the large junction of a SNAIL $L_{J}=\varphi_{0}/E_{J}$, the flux-dependent SNAIL inductance $L_{s}(\Phi)=L_{J}/c_{2}(\Phi)$ and the flux offset $\phi_{\rm min}=\varphi_{\rm min}\,\phi_{0}$. After applying harmonic balance to this equation we obtain

[TABLE]

Now the eigenmodes of the circuit can be found using the equations (36), (38) and (41). It is convenient to combine them in a matrix form

[TABLE]

where we have introduced the characteristic impedance of the transmission line resonator $Z_{c}=\sqrt{\ell/c}$. By requiring that the determinant of Eq. (42) be zero we obtain the equation for eigenfrequencies, while the eigenvectors define the coefficients needed to calculate the flux profile of the eigenmodes in Eq. (33). The equation for eigenfrequencies is

[TABLE]

where we have explicitly indicated the flux dependence of the SNAIL inductance and defined the parameter $\omega_{0}=\frac{\pi}{2}\frac{v}{d}$, whose physical meaning is the frequency of the fundamental harmonic of the resonator in the absence of the array.

We will use the value of the flux field at the right boundary $\phi_{n}(t)=\phi_{n}(d,t)$ as the canonical coordinate for each mode. It is a natural choice, since this variable is the one to which we couple the input transmission line via a weak coupling capacitor. In terms of $\phi_{n}$ we can write the components of the eigenvector of Eq.(42) as

[TABLE]

and the flux distribution in each eigenmode as

[TABLE]

Note also that the constant flux offset between the two arms of the resonator is given by Eq. (40). After substituting the flux distribution given by Eq. (33) into the Lagrangian and calculating the $x$-integral, we can diagonalize it and reduce to the collection of independent harmonic oscillators

[TABLE]

with the mode capacitance and inductance given by

[TABLE]

and the mode resonant frequency $\omega_{n}=1/\sqrt{L_{n}C_{n}}$ and impedance $Z_{n}=\sqrt{L_{n}/C_{n}}$. Note that in general both $C_{n}$ and $L_{n}$ can depend on $\Phi$.

Nonlinearities of the SPA

Having considered the linear properties of the SPA resonator, we will now focus on deriving the dependence of nonlinearities $g_{3}$ and $g_{4}^{*}$ on the design parameters of the circuit. Such a perturbative approach, with the first step consisting of solving the linearized problem, is valid if $Z_{n}\ll R_{Q}$ Manucharyan (2012), where $Z_{n}$ is the mode impedance and $R_{Q}=\hbar/(2e)^{2}$ is the reduced resistance quantum.

Because of the nonlinearity of the array, exciting a particular nonlinear mode $n$ leads to the oscillations not only at the frequency $\omega_{n}$, but also at higher harmonics of this frequency $2\omega_{n}$, $3\omega_{n}$ etc. The amplitudes of the oscillations at the multiples of the mode frequency $\omega_{n}$ are small if the mode is not significantly excited, and if the system is weakly nonlinear. This allows to solve the problem in a perturbative manner, similar to the calculation in Wallquist et al. (2006), using the eigenmodes found in the previous section as the first approximation.

In this section, we focus on the lowest mode $n=1$ of the circuit. This single-mode approximation is justified in the SPA resonator for two reasons. First, the SPA mode is in the nearly “lumped” regime, where almost all inductance comes from the array and all capacitance comes from the microstrip pads. We have experimentally verified using two-tone spectroscopy that the next mode of the SPA is at $25\rm\;GHz$ (data not shown). Since the detuning of the higher modes is much larger than the lowest mode frequency, we expect the effect of these modes on the nonlinearities of the lowest mode to be small. Second, all the driving signals in our experiments are significantly below these higher modes, and, therefore, do not excite them.

To be consistent with Eq. (22), we will denote the frequency of the mode as $\omega_{a}$ and canonical coordinate as $\phi_{a}$. The flux profile of the excited mode should still obey the equation (28) in the left and right arms of the transmission line resonator, and, thus, we can write it as

[TABLE]

instead of Eq. (45). The new coefficients $u_{k}^{L/R}(t)$ and $v_{k}^{L/R}(t)$ should obey the equation $[\partial_{t}^{2}+(n\omega_{a})^{2}](\cdot)=0$. In the linear approximation, we have already found the coefficients of this expansion in the previous section. The DC part obeys $u_{0}^{R}-u_{0}^{L}=M\phi_{{\rm min}}$, and the only nonzero coefficients are u_{1}^{R/L}=\pm\phi_{a}\cos\big{(}\frac{\pi}{2}\frac{\omega_{a}}{\omega_{0}}\big{)} and v_{1}^{R/L}=\phi_{a}\sin\big{(}\frac{\pi}{2}\frac{\omega_{a}}{\omega_{0}}\big{)}; the rest of the harmonics with $k>1$ are not excited. Using the perturbative approach, we will find that the amplitudes of these higher harmonics are suppressed, for example we will find that $u_{2}^{R/L}\propto\phi_{a}^{2}$. But first, by using the boundary conditions and harmonic balance we can establish a useful relation between $u_{k}^{L/R}$ and $v_{k}^{L/R}$, which will allow us to reduce the number of independent coefficients. This way we find that $u_{k}^{L}=-u_{k}^{R}$ and v_{k}^{R}=u_{k}^{R}\tan\big{(}\frac{\pi}{2}\frac{k\omega_{a}}{\omega_{0}}\big{)}, and thus for each $k$ we only need to use $u_{k}^{R}$. Let us now return to the continuity condition in Eq. (39), which we have previously linearized, but this time treat it as nonlinear. We will set the amplitudes of all higher modes with $n>1$ to zero, but include the higher harmonics of the first mode, truncated at the second order. Then Eq. (39) reduces to

[TABLE]

By including the time dependence $u_{1}^{R}(t)=u_{1}^{R}\sin\omega_{a}t$ and $u_{2}^{R}(t)=u_{2}^{R}\cos 2\omega_{a}t$ and keeping the relevant contributions (in the spirit of perturbation theory) we solve Eq. (50) by harmonic balance.

• $\circ$

DC harmonic

[TABLE]

• $\circ$

Second harmonic $2\omega_{a}$

[TABLE]

• $\circ$

First harmonic $\omega_{a}$

[TABLE]

Note that to write the relation (53) we had to use the expressions (51) and (52) to eliminate $u_{0}^{R}$, $u_{0}^{L}$ and $u_{2}^{R}$.

The equation (53) is essentially equivalent to the dispersion relation Eq. (43), corrected by an additional term proportional to the intensity of the excitation – a result of the Kerr effect. By subtracting Eq. (43) from Eq. (53) and replacing the finite difference with the derivative we obtain

[TABLE]

where we have used the relation u_{1}^{R}=\phi_{a}\cos\big{(}\frac{\pi}{2}\frac{\omega_{a}}{\omega_{0}}\big{)} to express the resulting frequency shift $\Delta\omega$ in terms of the canonical coordinate $\phi_{a}$.

We can also compute the linear shift of the resonant frequency due to the finite average mode population $\overline{n}$ using the model described by Eq. (22), which results in $\Delta\omega=12g_{4}^{*}\overline{n}$ with $g_{4}^{*}=g_{4}-5g_{3}^{2}/\omega_{a}$. To directly compare these two results for $\Delta\omega$, we need to express $\phi_{a}^{2}$ in terms of the average photon number as well, which can be done using the conventional method of canonical quantization. This leads to $\phi_{a}=\phi_{{\rm ZPF}}(a+a^{\dagger})$, where $\phi_{{\rm ZPF}}=\sqrt{\hbar Z_{1}/2}$ and $Z_{1}=\sqrt{L_{1}/C_{1}}$ is the impedance of the mode calculated with the help of Eqs. (47) and (48). Comparison of the two expressions for the shift $\Delta\omega$ leads to the following

[TABLE]

The derivation of $g_{4}^{*}$ presented here relies crucially on the careful treatment of higher harmonics of the excited mode. Such treatment is necessary for any higher-order nonlinearity of the system. However, to find the lowest-order nonlinearity $g_{3}$, we can employ the traditional black-box quantization method Nigg et al. (2012). In this method, the circuit is first linearized, and the nonlinearity is then applied to the phase drop across the nonlinear elements of the original circuit. Therefore, the cubic contribution to the Hamiltonian in Eq. (22) coming from the array can be written simply as

[TABLE]

Plugging in the expression $\phi_{a}=\phi_{{\rm ZPF}}(a+a^{\dagger})$, we recover the cubic term $g_{3}(a+a^{\dagger})^{3}$ with $g_{3}$ given by

[TABLE]

Fitting

To conclude this section, we have calculated the coefficients $\omega_{a}({\cal S},\Phi)$, $g_{3}({\cal S},\Phi)$ and $g_{4}({\cal S},\Phi)$ of the SPA Hamiltonian as functions of the design parameters ${\cal S}=\{M,\omega_{0},Z_{c},L_{J},\alpha\}$ of the device, as well as the external magnetic flux $\Phi$. In the experiment, we fit the resonance frequency $\omega_{a}(\Phi)$ in order to refine the knowledge of the parameters in $\cal S$. Those are obtained in the following way: $Z_{c}=45.8{\rm\>\Omega}$ is defined by the microstrip geometry, $L_{J}=38\,{\rm pH}$ is calculated using the room-temperature measurements of junction’s resistance and the Ambegaokar-Baratoff formula, $M=20$ is the number of SNAILs in the array, and $\omega_{0}=2\pi\times 16.0\,{\rm GHz}$ and $\alpha=0.065$ are fitted and are found to be close to their design values. The measurement results of $g_{4}^{*}$ and $g_{3}$ are compared to the predicted values using the refined parameters $\cal S$ in Fig. 1(b).