Kerr-free three-wave mixing in superconducting quantum circuits
V. V. Sivak N. E. Frattini V. R. Joshi A. Lingenfelter S. Shankar M. H. Devoret





vladimir.sivak@yale.edu
Department of Applied Physics, Yale University, New Haven, CT 06520, USA

Department of Applied Physics, Yale University, New Haven, CT 06520, USA

Department of Applied Physics, Yale University, New Haven, CT 06520, USA

Department of Applied Physics, Yale University, New Haven, CT 06520, USA

Department of Applied Physics, Yale University, New Haven, CT 06520, USA

michel.devoret@yale.edu
Department of Applied Physics, Yale University, New Haven, CT 06520, USA






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 operating 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.




Kerr-free three-wave mixing in superconducting quantum circuits
    M. H. Devoret
    December 30, 2023
===============================================================





§ 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 <cit.>. Tunable-strength frequency conversion Hamiltonians enter in quantum state transfer between remotely separated modes <cit.>. Bilinear couplings with well-controlled phases are crucial for realizing active nonreciprocity in few-body systems such as parametric circulators and directional amplifiers <cit.>, as well as for simulating many-body physics of topological band structure using photonic systems <cit.> and implementing topological traveling-wave amplifiers <cit.>. Furthermore, bilinear interactions are essential auxiliaries in nonlinear operations, such as the universal two-qubit entangling gate eSWAP <cit.> or Clifford gates on GKP-encoded logical qubits <cit.>. 

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_ DPA=-Δ_b b^† b + g(b^2+b^†2), where b refers to the photon annihilation operator in the rotating frame. It can be constructed with four-wave mixing by pumping n photons at detuning Δ into the Duffing oscillator realized with, for instance, a capacitively shunted Josephson junction or SQUID <cit.>. In such an implementation, the junction four-wave mixing with strength K induces three terms: (i) the desired bilinear squeezing term with g = Kn/2, (ii) the spurious Stark shift Δ_b - Δ = - 2Kn, and (iii) the residual quartic term δ H = Kb^† 2b^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 δ H leads to the distortion of quantum states of light in a microwave cavity at the level of a few photons <cit.>, to “bananization” of squeezing <cit.>, and to saturation of parametric amplifiers <cit.>. In addition, the spurious Stark shift of the type (ii) limits the achievable gain in the multipumping schemes for directionality <cit.> and deteriorates phase matching in traveling-wave amplifiers (TWPA) <cit.>.

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 <cit.>. Alternatively, non-degenerate three-wave mixing is also available via the Josephson ring modulator (JRM) <cit.>, and is utilized in the Josephson parametric converter (JPC) <cit.>. In both of these implementations the residual Kerr has been determined to cause amplifier saturation <cit.>. Recently a two-terminal three-wave mixing element, the Superconducting Nonlinear Asymmetric Inductive eLement (SNAIL), has been introduced <cit.> 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 <cit.> and inductively shunted JRM <cit.>.



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  dB compression power of -102  dBm on par with the best published results <cit.>. 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. 



§ SNAIL PARAMETRIC AMPLIFIER


The device used in this work is the SNAIL parametric amplifier (SPA) <cit.>. Its electrical circuit, shown in Fig. <ref>(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 Φ piercing each loop. This knob enables the tuning of the SPA resonant frequency ω_a in the range 6.2-7.2  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 Ω environment via a finger-capacitor and sets the energy damping rate of the SPA mode κ/2π∈150-240  MHz, while a weakly coupled pump port is used to supply the pump. The Hamiltonian of such system can be approximated as

    H_ SPA/ħ= ω_a(Φ) a^† a + ∑_n=3^∞ g_n(Φ) (a^† +a)^n.

The dependence of the most relevant Hamiltonian parameters on Φ is shown in Fig. <ref>(b) and is calculated and measured using the methods elaborated in Ref. <cit.> and in Appendices <ref> and <ref>.

When a strong off-resonant pump is applied at frequency ω_p=2(ω_a+Δ), this leads to the effective Hamiltonian of a driven oscillator with Kerr nonlinearity

    H_ eff/ħ=-Δ_b b^† b + g(b^2+b^†2) + K/2b^†2b^2,

written in terms of the photon annihilation operator b in the frame rotating at ω_p/2. The parameters Δ_b, g and K of the effective Hamiltonian in Eq. (<ref>) depend on all nonlinearities of the initial Hamiltonian in Eq. (<ref>) 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

    Δ_b(n_p)    ≈Δ - (32/3g_4 - 28g_3^2/ω_a) n_p, 
    
    g(n_p)    ≈ 2g_3√(n_p), 
    
    K(n_p)    = 12g_4^* +O(n_p),  where  g_4^*=g_4 - 5g_3^2/ω_a.
 
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(Φ) can go through zero and change sign while g(Φ) 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 ω_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. <cit.> and Appendix <ref>. The result is given by 

    G = 1 + (2|g|κ)^2/(κ^2/4+Δ_iΔ_s-4|g|^2)^2+κ^2/4(Δ_i-Δ_s)^2,

where Δ_s/i(n_p,n_s,n_i)=Δ_b(n_p)-K(n_p)×(n_s/i+2n_i/s). To obtain large gain, one needs to tune the model described by Eq. (<ref>) close to the parametric instability threshold at n_p=(κ^2+4Δ_b^2)/(8g_3)^2. However, the spurious additional Stark shift K(n_p)×(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. (<ref>), 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. (<ref>) becomes comparable to the n_p-independent term, and the lowest-order perturbative approximation for K(n_p) given by Eq. (<ref>) breaks down. This breakdown is particularly important in the region of suppressed Kerr.



§ 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 ω_p is difficult in practice because it is obscured by the amplification process arising from the second term in Eq. (<ref>). We instead measure the Stark shift caused by n steady-state photons of a strong near-resonant drive placed at ω_d=2π×7.8  GHz, about seven linewidths away from the resonance, which emulates the effect of the pump. The resulting Δ_ Stark(n) is shown in Fig. <ref>(a) for different fluxes in the range (0.24-0.44)Φ_0, where Φ_0=h/2e is the superconducting magnetic flux quantum. To the lowest order in n, this Stark shift is given by 24g_4^*n. The upper bound on n in Fig. <ref>(a) is roughly determined by the onset of chaotic behavior, caused by the excitation of free-particle-like states <cit.>. The inset in Fig. <ref>(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 <cit.>, occurs due to the higher-order terms in the SPA Hamiltonian given by Eq. (<ref>). 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 ω_s1 and ω_s2 are sent to the resonator, and the sidebands at 2ω_s1-ω_s2 and 2ω_s2-ω_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

    IIP_3 = κ/|K|1/(√(G)+1)^3ħω_a κ,

as proved in <cit.>. The measurement result of IIP_3 with the pump off is shown in Fig. <ref>(b) in blue, and the corresponding g_4^* extracted according to Eqs. (<ref>) and (<ref>) is shown in Fig. <ref>(b) in good agreement with that extracted from the linear Stark shift. Eq. (<ref>) 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. (<ref>), and, thus, the peak shape in IIP_3 at Φ=0.40 Φ_0 would remain unchanged with increasing gain. However, in the experiment we find that pumping at ω_p=2ω_a (i.e. Δ=0) leads to the complete disappearance of the peak in IIP_3 at gains G>5  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. (<ref>) is intentionally suppressed. Outside this region, the approximation K(n_p)= const works well up to n_p required for G=20  dB, see Appendix <ref>. 

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. <ref>(a), changing n_p at a fixed Φ can influence the sensitivity of the Kerr-induced terms Δ_s/i(n_p,n_s,n_i) in Eq. (<ref>) to n_s/i, providing us the tool that we need in order to find the effective Kerr-free point. At a fixed Φ, we can vary the pump photon number n_p required to reach G=20  dB by changing the pump detuning Δ, which had been set to zero in the analysis so far. 



§ STABILITY REGIONS


The model described by Eq. (<ref>) exhibits a variety of new effects when Δ 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. <ref>(a) for Φ = 0.30Φ_0. The colors in this diagram represent the three regions in the space of parameters Δ_b and g, in which the classical nonlinear dissipative dynamics of a complex amplitude b has qualitatively different phase portraits <cit.>. 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 <cit.>. This potential can have one, two or three minima depending on the parameters Δ_b and g, and is sketched for illustrative purposes in Fig. <ref>(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  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. <ref>(b) right upper panel and Appendix <ref>. In addition, the saturation curve taken at the signal frequency near ω_p/2 monotonically decreases with increasing signal power, see Fig. <ref>(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 ω_p/2 even with no applied signal power, an effect termed as ”parametric self-oscillation” or ”period-doubling” <cit.>. This region is separated from the rest of the diagram by the parametric instability threshold, where the gain of a parametric amplifier diverges <cit.>. In the K<0 case, due to the Stark shift produced by n_p, there is a maximal positive detuning Δ 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  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 <cit.>.

On the negative detuning side, the G=20  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. <ref>(b) upper middle panel. Nevertheless, as shown in Fig. <ref>(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 (<ref>) for gain G as a function of n_s/i(G), as shown in Appendix <ref>. The characteristic shape of the saturation curve in this regime was called “shark fin” in Ref. <cit.>, 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  dB when the input difference is less than 1  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 <cit.>. Therefore, the NVR is not appreciably degraded at small red detuning inside region (III). Following the 20  dB curve in Fig. <ref>(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 <cit.>. Hence, we attribute the significant increase of NVR at large negative detuning, shown in Fig. <ref>(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 Δ=-500  MHz. In such a state, the amplifier is no longer quantum-limited and cannot be described by the conventional theory that leads to Eq. (<ref>). 

The stability diagram in Fig. <ref>(a) was constructed at Φ=0.30 Φ_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. (<ref>) to fail, we need to perform an empirical search for the dressed Kerr-free point in the space of parameters Φ and Δ.



§ DRESSED KERR-FREE POINT


We show in Fig. <ref> the results of the measurement of gain saturation (upper panels) and NVR (lower panels) for high-resolution sweeps of detuning Δ done at five different fluxes Φ in the region of reduced Kerr. At every detuning, the pump power is adjusted to obtain the small-signal gain of 20  dB. The input signal power is then swept beyond the 1  dB compression point (solid black). The left and right boundaries of the Δ sweep at each Φ correspond to detunings where a small-signal gain G=20  dB could not be achieved at any applied pump power. Moving from smaller to larger Φ, the Δ tuning range decreases from about 900  MHz at Φ=0.19 Φ_0 to about 200  MHz at Φ=0.48 Φ_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 Φ=0.34 Φ_0, Δ=-95  MHz. There, the 1 dB compression power peaks at a value P_ 1 dB=-102  dBm which is an order of magnitude larger than at Δ=0. The tunable bandwidth over which we observe the P_1 dB increase is about 200  MHz. We measure a similar improvement in the IMD properties with peak IIP_3=-94  dBm, see Appendix <ref>. 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.



§ 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_ 1 dB=-102  dB, on par with the best published results, but with room for further improvement. In Appendix <ref> 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. (<ref>). 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 <cit.> or the TWPA with array of RF SQUIDs <cit.> where Kerr cancellation will improve phase matching <cit.>. 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.








§ DEVICE PARAMETERS AND MEASUREMENT SETUP  


The SPA device package is shown in Fig. <ref>. Josephson junctions of the SNAIL are formed by Al/AlO_x/Al layers deposited using the Dolan bridge shadow evaporation process (critical current density j_c=120  A/ cm^2), and have critical currents I_c=8.5  uA and I_c=0.85  uA. The microstrip transmission line resonator is formed by a 2 μ m thick silver film deposited on the back of a 300 μ 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 TMM10 i printed circuit board (PCB), which is soldered to the gold-plated aluminum box, see Fig. <ref>(a). The first resonant mode of the box is at 18  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. <ref>(a). The aluminum box acts as a partial shield against stray external magnetic fields. For better protection, an additional cryogenic μ-metal shield is fitted around the aluminum box. The device is thermally anchored to the base stage of a dilution refrigerator (T=24  mK) through this shield. It is mounted back-to-back with the 4-8 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 <cit.>. 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 ω_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  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 <cit.>, which was not available in our setup. However, our claim that the system is nearly quantum-limited when NVR≈ 8.3  dB is supported by the independent measurement of quantum efficiency η=0.6 in the experiment <cit.> done with the SPA.



§ INTERMODULATION DISTORTION (IMD)   



In the IMD measurement the two signals at ω_s1 and ω_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 ω_p/2 by amount δ_1≡(ω_s1+ω_s2)/2 - ω_p/2=2π×500  kHz to avoid phase-sensitive amplification. The detuning is chosen small compared to the amplifier 3dB-bandwidth B∼2π× 25  MHz. The separation of the two tones is δ_2=ω_s2-ω_s1=2π×100  kHz. Due to the nonlinear intermodulation distortion, the sidebands at 2ω_s2-ω_s1 and  2ω_s1-ω_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. <ref>, and exhibits a peak value IIP_3=-94  dB at the dressed  Kerr-free point similar to the peak in P_1 dB.






§ SATURATION POWER AND SHARK FINS 


Saturation of a good amplifier happens at large average signal and idler photon numbers n_s/i≫1/2, justifying the use of the semiclassical approximation, in which the phase preserving gain G at the signal frequency near ω_p/2 is given by Eq. (<ref>). The small-signal limit of this expression is given by

    G_0=1+4κ^2|g|^2/( Δ_b^2+κ^2/4 -4|g|^2)^2.

The goal of this section is to derive a closed-form expression for the input signal power P_ in at which the gain changes from G_0 to G (which directly gives the closed-form expression for P_1 dB) and to study the phenomenon of shark fins. Using input-output theory, we can establish the relation between the input power P_ in and the intraresonator populations n_s and n_i. For large gain G≫1 we obtain n_s≈ n_i≈ GP_ in/ħω_aκ with a relative error δ n_s/n_s∼1/√(G), which is about 10% for G=20  dB. In this approximation, we can rewrite Eq. (<ref>) in the form

    G=1+4κ^2|g|^2/(Δ_ eff^2+κ^2/4-4|g|^2)^2,
    Δ_ eff=Δ_b-3KGP_ in/ħω_aκ,

where Δ_ eff contains the additional Stark shift created by the signal and idler photons (this formula is given in <cit.>, but with a typo which is corrected here).  After plugging this into Eq. (<ref>) and solving for P_ in,
we obtain

    P_ in=ħω_aκ/3KG(Δ_b±√(4|g|^2-κ^2/4+2κ|g|/√(G))).




The Hamiltonian parameter g can be expressed through G_0 from Eq. (<ref>), which leads to

    P_ in=ħω_aκ^2/3KG{Δ_b/κ±√(Δ_b^2/κ^2+√(G_0)-√(G)/√(G_0G)√(Δ_b^2/κ^2+1/4))}.
 

In this expression both K and Δ_b depend on n_p, and thus implicitly on gain. When we consider a simplified case of small negative K, such that Δ_b≈Δ and K≈ 12g_4^*, we can further simplify Eq. (<ref>) to

    P_ in=ħω_aκ^2/36g_4^*G{Δ/κ±√(Δ^2/κ^2+√(G_0)-√(G)/√(G_0G)√(Δ^2/κ^2+1/4))}.


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. (<ref>) to the measured  P_1 dB at a fixed detuning Δ=0 as a function of flux Φ is shown in Fig. <ref>(a), and comparison to  the measured P_1 dB at a fixed flux Φ=0.30 Φ_0 as a function of detuning Δ is shown in Fig. <ref>(b). For the plotted theory lines, we have assumed that there exists a systematic miscalibration of the resonant frequency (and therefore the detuning Δ) by approximately 25  MHz. The compression power is very sensitive to such miscalibration especially near the detuning Δ=0,  as evident from Fig. <ref>(b). Accounting for this shift aligns the theory and data much better in both Δ and Φ axes, while without this correction we find a 6  dB systematic disagreement between theory and data in Fig. <ref>(a). 

The extrapolation of Eq. (<ref>) to detunings Δ<-300  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. (<ref>) is significantly modified at large negative detuning Δ  due to the secondary high-amplitude minima, which invalidates the derivation of P_1 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. (<ref>) to fluxes Φ≳ 0.30Φ_0 leads to the incorrect prediction of improved compression power at an unpumped Kerr-free point at Φ=0.40 Φ_0, as shown in Fig. <ref>(a). Using Eq. (<ref>) with K extracted from the IMD measurement at gain G=20  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. (<ref>) as a self-consistent equation of state, analogous to the Van der Waals equation for nonideal gas, which exhibits a first-order phase transition <cit.>. In this analogy, the detuning Δ plays the role of temperature, input signal power P_ in plays the role of pressure, and gain G plays the role of volume. As shown in Fig. <ref>(d), for the negative detuning Δ, the equation (<ref>) 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. <ref>(c).

Finally, note that the flux tuning range of the SPA shown in Fig. <ref>(a) corresponds to a wide tunable bandwidth of about 1  GHz. Moreover, at a fixed Φ the tunability can be accomplished with the help of pump detuning Δ, which corresponds to a tuning range of approximately 700  MHz, as shown in Fig.<ref>(b) for Φ=0.30 Φ_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  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. <cit.>.




§ STABILITY DIAGRAM OF THE PARAMETRICALLY DRIVEN KERR OSCILLATOR   


In this section we derive the separation lines between the regions of the stability diagram at Φ=0.30 Φ_0 shown in Fig. <ref>(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 ω_nm = nω_p + mω_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 ω_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:


  ∘ Pump ω_p, signal ω_s and idler ω_i=ω_p-ω_s comprise a minimal set required for quantum-limited amplification.

  ∘ 2ω_p, 2ω_s, 2ω_i, ω_p+ω_s, ω_p+ω_i, ω_s-ω_i and DC are required in the consistent calculation that takes into account linear Stark shift. These harmonics lead to ∼ g_3^2/ω_a corrections to K in Eq. (<ref>) and to pump-induced Stark shift in Eq. (<ref>). 

  ∘ ω_p/2, ω_p/2-ω_s and ω_p/2-ω_i account for a possible period-doubling effect.




Following the steps of <cit.> and Refs. <cit.>, 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) <cit.> 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

    a_-ω^† = ω_a-ω/ω_a+ωa_ω.

The a_-ω^† is usually dropped from the frequency-domain QLE for a_ω if all the signals arriving at the resonator are near the resonance. Eq. (<ref>) can be derived from the relation Φ̇=C_a Q between the canonical coordinate Φ of the mode and its canonical momentum Q, written in second-quantized form.

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

    α(t) = ∑_x (α_x e^-iω_x t+ω_a-ω_x/ω_a+ω_xα_x^* e^iω_x t),

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 ω_s, ω_i and ω_p/2, denoted as α_s, α_i and α_h, respectively, as follows:


    (ω +Δ_b+iκ/2)α_s     = u_s + (4g_3α_p+12g_4^*α_h^2)α_i^* + 12g_4^*(|α_s|^2+2|α_i|^2+2|α_h|^2)α_s,
    
    (-ω+Δ_b+iκ/2)α_i     = u_i + (4g_3α_p+12g_4^*α_h^2)α_s^* + 12g_4^*(|α_i|^2+2|α_s|^2+2|α_h|^2)α_i,
    
    (Δ_b+iκ/2)α_h     = u_h + 4g_3α_pα_h^* + 12g_4^*(|α_h|^2+2|α_s|^2+2|α_i|^2)α_h,


where ω = ω_s-ω_p/2 is the signal detuning, Δ_b is defined in Eq. (<ref>), and u_s, u_i and u_h denote the drive strengths at corresponding frequencies, following the notations introduced in Ref. <cit.>. Note that the pump amplitude α_p=u_p/ω_a here does not include any corrections due to α_s, α_i and α_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 |α_s|^2 and |α_i|^2 can be neglected. Since no input is sent at ω_p/2, we can set the drive strength u_h to zero. Then equation (<ref>) reduces to

    (Δ_b-12g_4^*|α_h|^2+iκ/2 )α_h = 4g_3α_pα_h^*.

After multiplying it by the complex conjugate, we find a trivial solution α_h=0 and possible high-amplitude solutions

    |α_h|^2 = 1/12g_4^*(Δ_b±√((4g_3|α_p|)^2-κ^2/4)).


First, consider Δ_b>0 half-plane. In this case only the solution with the minus sign in Eq. (<ref>) is possible if the additional condition (4g_3|α_p|)^2-κ^2/4>Δ_b^2 is satisfied. This defines the separation line between regions (I) and (II) of the stability diagram. In region (I) only α_h=0 solutions exists and is stable. High-amplitude solution exists in region (II) and is stable, while α_h=0 is unstable.  

Next, consider Δ_b<0 half-plane. The requirement for |α_h|^2 being real leads to the condition (4g_3|α_p|)^2> κ^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. (<ref>) are possible. This yields five stationary points: a stable trivial α_h=0 point, a pair of symmetric unstable points defined by the plus sign in Eq. (<ref>) and a pair of symmetric stable points defined by the minus sign. The two unstable stationary points merge with the stable α_h=0 point on the line which separates regions (II) and (III), defined by the equation (4g_3|α_p|)^2=κ^2/4+Δ_b^2. Above this line in region (II) the α_h=0 point becomes unstable and the two symmetric high-amplitude states remain stable. 

The stability diagram plotted in coordinates Δ and n_p (which directly translate into the experimental control knobs pump frequency ω_p and power P_p) is shown in Fig. <ref>(a). Note that similar stability diagrams can be found in for example Refs. <cit.>. 

Finally we note that by setting α_h=0 (where appropriate, see Appendix <ref>) the system of equations (<ref>)-(<ref>) together with the input-output relation can be solved to yield the expression (<ref>) for gain in the presence of a finite input signal.




§ 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. <cit.>, where steady improvement of P_ 1dB of an SPA was achieved by varying the design parameters in the direction corresponding to the increase of the coupling κ and the decrease of the overall level of nonlinearity. 

Continuing to improve P_1 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

    η_p = GP_ 1 dB/P_p.


As shown in Fig. <ref>, 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  mK for -20  dBm of applied power at the base stage. Typically, an attenuator is used to thermalize the pump to this stage. Assuming 20   dB of attenuation, this puts the limit on allowable pump power at the plane of the device close to -40  dBm. This limit is fuzzy (shown as color gradient in Fig. <ref>), 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 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. <cit.>. These advances are left for future work.





§ 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

    H_ SPA/ħ=ω_aa^†a+g_3(a+a^†)^3+g_4(a+a^†)^4+...


The goal of this section is to obtain the dependence of the parameters ω_a, g_3 and g_4^*≡ g_4-5g_3^2/ω_a on the design parameters of the circuit. The expressions (<ref>), (<ref>) and (<ref>) obtained here were used for fitting and making the design choices in this work and in Ref. <cit.>.



 §.§ 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 ℓ. 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 ϕ(x). For simplicity, we also introduce the notation ϕ^R/L=x→±0limϕ(x). 

The Lagrangian can be written as

    ℒ=∫_-d^-0[c/2(∂_tϕ)^2-1/2ℓ(∂_xϕ)^2]dx+∫_+0^d[c/2(∂_tϕ)^2-1/2ℓ(∂_xϕ)^2]dx-MU_S(ϕ^R-ϕ^L/Mϕ_0),

where ϕ_0=ħ/2e is the reduced flux quantum and U_S(φ) is the potential energy of the SNAIL with a phase drop φ,

    U_S(φ) = -E_J[ αcosφ +3cos(φ_ ext-φ/3) ]

where E_J is the Josephson energy of the large junctions, α is a ratio of the junction inductances in the two arms of the SNAIL, and φ_ ext=2πΦ/Φ_0 is the phase bias corresponding to the external magnetic flux Φ, see <cit.>. 

The action for this Lagrangian is simply given by S=∫_t_1^t_2ℒ(t)dt, and can be split into three separate parts S=S_ L+S_ R+S_ SNAIL corresponding to the separate contributions in the Lagrangian in Eq. <ref>. The variation of action δ S contains three terms

    δ S_ L    =-∫_t_1^t_2dt[1/ℓ∂_xϕ]δϕ|_x=-d^x=-0-∫_t_1^t_2dt∫_-d^-0[c∂_t^2ϕ-1/ℓ∂_x^2ϕ]δϕ dx,
    δ S_ R    =-∫_t_1^t_2dt[1/ℓ∂_xϕ]δϕ|_x=+0^x=d-∫_t_1^t_2dt∫_+0^d[c∂_t^2ϕ-1/ℓ∂_x^2ϕ]δϕ dx,
    δ S_ SNAIL    =-M∫_t_1^t_2dt∂ U_S(φ)/∂φδϕ^R-δϕ^L/Mϕ_0|_φ=ϕ^R-ϕ^L/Mϕ_0.


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


  ∘ In the bulk of the transmission line resonator, for x∈(-d,0) and x∈(0,d), the flux field ϕ(x) obeys the wave equation

    ∂_t^2ϕ-v^2∂_x^2ϕ=0.

where v=1/√(ℓ c) is the phase velocity.

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

    ∂_xϕ=0.


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

    1/ϕ_0∂ U_S(φ)/∂φ|_φ=ϕ^R-ϕ^L/Mϕ_0=1/ℓ∂_xϕ^L/R,

which also implies the current conservation ∂_xϕ^L=∂_xϕ^R across the array. Strictly speaking, the bare SNAIL potential U_S(φ) is 6π-periodic in φ. However, we are interested in a time scale on which the phase slips can be neglected, and, thus, the fact that φ is compact does not matter. In this approximation, we can assume that the phase φ is localized near the bottom of one of the many equivalent potential minima and expand the potential around that point (denoted as φ_ min) 

    U_S(φ)=E_J(c_2/2!(φ-φ_ min)^2+c_3/3!(φ-φ_ min)^3+c_4/4!(φ-φ_ min)^4+...),

where all dimensionless expansion coefficients c_k≡ c_k(Φ,α) depend only on the external magnetic flux Φ and junction inductance ratio α. Using this expansion truncated to the fourth order, we can rewrite the continuity condition at x=±0 as

    E_J/ϕ_0[c_2(ϕ^R-ϕ^L/Mϕ_0-φ_ min)+c_3/2!(ϕ^R-ϕ^L/Mϕ_0-φ_ min)^2+c_4/3!(ϕ^R-ϕ^L/Mϕ_0-φ_ min)^3]=1/ℓ∂_xϕ^L/R.

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

Instead of working with the continuous flux field ϕ(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. (<ref>) in the left and right arms of the transmission line resonator

    ϕ(x,t)     =(a_0^L/R+b_0^L/Rx)+∑_n[a_n^L/Rcos(ω_nx/v)+b_n^L/Rsin(ω_nx/v)],
    ϕ^L/R    =a_0^L/R+∑_na_n^L/R,

where the coefficients a_n^L/R(t) and b_n^L/R(t) should satisfy [∂_t^2+ω_n^2](·)=0. We need to rewrite the boundary condition in Eq. (<ref>) and the continuity condition in Eq. (<ref>) using the eigenmode decomposition. Starting with the boundary condition, by differentiating Eq. (<ref>) we obtain

    b_0^L/R+∑_nω_n/v[± a_n^L/Rsin(ω_nd/v)+b_n^L/Rcos(ω_nd/v)]=0,

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

    ± a_n^L/Rsin(ω_nd/v)+b_n^L/Rcos(ω_nd/v)=0,

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. (<ref>) using Eq. (<ref>) and its derivative

    ∂_xϕ^L/R=b_0^L/R+∑_nω_n/vb_n^L/R.


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 ∂_xϕ^L=∂_xϕ^R instead. With the help of Eq. (<ref>) and applying harmonic balance, this reduces to

    b_n^L=b_n^R.


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

    E_J/ϕ_0[c_2(∑_n=0^∞a_n^R-a_n^L/Mϕ_0-φ_ min)+c_3/2!(∑_n=0^∞a_n^R-a_n^L/Mϕ_0-φ_ min)^2+c_4/3!(∑_n=0^∞a_n^R-a_n^L/Mϕ_0-φ_ min)^3]-1/ℓ∑_nω_n/vb_n^R=0.


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 ω_n. The lowest such eigenmode of the device is used to implement the SPA.



 §.§ Eigenmode decomposition


The linear part of Eq. (<ref>) is simply

    1/L_s(∑_n=0^∞a_n^R-a_n^L/M-ϕ_ min)-1/ℓ∑_nω_n/vb_n^R=0,

where we have introduced the inductance of the large junction of a SNAIL L_J=φ_0/E_J, the flux-dependent SNAIL inductance L_s(Φ)=L_J/c_2(Φ) and the flux offset ϕ_ min=φ_ min ϕ_0. After applying harmonic balance to this equation we obtain

    a_0^R-a_0^L/M-ϕ_ min    =0,
    1/L_sa_n^R-a_n^L/M-1/ℓω_n/vb_n^R    =0.


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

    [[  tan(ω_nd/v)            1            0            0;            0            0 -tan(ω_nd/v)            1;            0            1            0           -1;           -1            0            1 -L_sMω_n/Z_c ]]([ a_n^L; b_n^L; a_n^R; b_n^R ])=0,

where we have introduced the characteristic impedance of the transmission line resonator Z_c=√(ℓ/c). By requiring that the determinant of Eq. (<ref>) 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. (<ref>). The equation for eigenfrequencies is

    2Z_c/ML_s(Φ)=ω_ntan(π/2ω_n/ω_0),

where we have explicitly indicated the flux dependence of the SNAIL inductance and defined the parameter ω_0=π/2v/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 ϕ_n(t)=ϕ_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 ϕ_n we can write the components of the eigenvector of Eq.(<ref>) as

    a_n^R=ϕ_ncos(π/2ω_n/ω_0),     a_n^L=ϕ_n[cos(π/2ω_n/ω_0)-ML_sω_n/Z_csin(π/2ω_n/ω_0)],     b_n^L/R=ϕ_nsin(π/2ω_n/ω_0),

and the flux distribution in each eigenmode as

    ϕ_n(x,t)=ϕ_n[ sgn xcos(ω_nx/v)cos(π/2ω_n/ω_0)+sin(π/2ω_n/ω_0)sin(ω_nx/v)].


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

    ℒ=∑_n(C_nϕ̇_n^2/2-ϕ_n^2/2L_n),

with the mode capacitance and inductance given by

    C_n    =1/2ω_nZ_c[πω_n/ω_0+sin(πω_n/ω_0)],
    
    L_n^-1    =ω_n/2Z_c[πω_n/ω_0+sin(πω_n/ω_0)],

and the mode resonant frequency ω_n=1/√(L_nC_n) and impedance Z_n=√(L_n/C_n). Note that in general both C_n and L_n can depend on Φ. 



 §.§ 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≪ R_Q <cit.>, where Z_n is the mode impedance and R_Q=ħ/(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 ω_n, but also at higher harmonics of this frequency 2ω_n, 3ω_n etc.  The amplitudes of the oscillations at the multiples of the mode frequency ω_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 <cit.>, 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  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. (<ref>), we will denote the frequency of the mode as ω_a and canonical coordinate as ϕ_a. The flux profile of the excited mode should still obey the equation (<ref>) in the left and right arms of the transmission line resonator, and, thus, we can write it as

    ϕ_a(x,t)=u_0^L/R+∑_k[u_k^L/Rcos(kω_ax/v)+v_k^L/Rsin(kω_ax/v)],

instead of Eq. (<ref>). The new coefficients u_k^L/R(t) and v_k^L/R(t) should obey the equation [∂_t^2+(nω_a)^2](·)=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ϕ_ min, and the only nonzero coefficients are u_1^R/L=±ϕ_acos(π/2ω_a/ω_0) and v_1^R/L=ϕ_asin(π/2ω_a/ω_0);
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∝ϕ_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^Rtan(π/2kω_a/ω_0), and thus for each k we only need to use u_k^R. Let us now return to the continuity condition in Eq. (<ref>), 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. (<ref>) reduces to

    ∑_m=1^3c_m+1/ϕ_0^m-1m!(∑_k=1^22u_k^R/M+u_0^R-u_0^L/M-ϕ_ min)^m-ω L_J/Z_c∑_k=1^2ku_k^Rtan(π/2kω_a/ω_0)=0.


By including the time dependence u_1^R(t)=u_1^Rsinω_at and u_2^R(t)=u_2^Rcos2ω_at   and keeping the relevant contributions (in the spirit of perturbation theory) we solve Eq. (<ref>) by harmonic balance. 


  ∘ DC harmonic

    u_0^R-u_0^L/M-ϕ_ min=-c_3/c_21/M^2(u_1^R)^2/ϕ_0.


  ∘ Second harmonic 2ω_a

    u_2^R=c_3/c_2-ω_aL_JM/Z_ctan(πω_a/ω_0)1/2M(u_1^R)^2/ϕ_0.


  ∘ First harmonic ω_a

    [2c_2/M-ω_aL_J/Z_ctan(π/2ω_a/ω_0)]+[c_4-2c_3^2/c_2-c_3^2/c_2-ω_aL_JM/Z_ctan(πω_a/ω_0)]1/M^3(u_1^R)^2/ϕ_0^2=0.

Note that to write the relation (<ref>) we had to use the expressions (<ref>) and (<ref>) to eliminate u_0^R, u_0^L and u_2^R.

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

    Δω=Z_c/2L_Jcos^2(π/2ω_a/ω_0)sin(πω_a/ω_0)/π/ω_0Z_c/ML_s+sin^2(π/2ω_a/ω_0)1/M^3[c_4-2c_3^2/c_2-c_3^2/c_2-ω_aL_JM/Z_ctan(πω_a/ω_0)]ϕ_a^2/ϕ_0^2,

where we have used the relation u_1^R=ϕ_acos(π/2ω_a/ω_0) to express the resulting frequency shift Δω in terms of the canonical coordinate ϕ_a. 

We can also compute the linear shift of the resonant frequency due to the finite average mode population n using the model described by Eq. (<ref>), which results in Δω=12g_4^*n with g_4^*=g_4-5g_3^2/ω_a. To directly compare these two results for Δω, we need to express ϕ_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 ϕ_a=ϕ_ ZPF(a+a^†), where ϕ_ ZPF=√(ħ Z_1/2) and Z_1=√(L_1/C_1) is the impedance of the mode calculated with the help of Eqs. (<ref>) and (<ref>). Comparison of the two expressions for the shift Δω leads to the following

    g_4^*=1/12ω_asin^2(πω_a/ω_0)/c_2M^2tan(π/2ω_a/ω_0)[πω_a/ω_0+sin(πω_a/ω_0)]^2Z_c/R_Q[c_4-c_3^2/c_23+5(ω_aML_s/2Z_c)^2/1+3(ω_aML_s/2Z_c)^2],


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 <cit.>. 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. (<ref>) coming from the array can be written simply as

    ME_Jc_3/3!(ϕ^R-ϕ^L/Mϕ_0-φ_ min)^3→ϕ_a^3/ϕ_08/M^2L_Jc_3/3!cos^3(π/2ω_a/ω_0).


Plugging in the expression ϕ_a=ϕ_ ZPF(a+a^†), we recover the cubic term g_3(a+a^†)^3 with g_3 given by

    g_3=4Z_cc_3/3M^2L_J(cos^2(π/2ω_a/ω_0)/πω_a/ω_0+sin(πω_a/ω_0))^3/2√(Z_c/R_Q).




 §.§ Fitting


To conclude this section, we have calculated the coefficients ω_a( S,Φ), g_3( S,Φ) and g_4( S,Φ) of the SPA Hamiltonian as functions of the design parameters S={M, ω_0,Z_c,L_J,α} of the device, as well as the external magnetic flux Φ. In the experiment, we fit the resonance frequency ω_a(Φ) in order to refine the knowledge of the parameters in S. Those are obtained in the following way: Z_c=45.8 Ω is defined by the microstrip geometry, L_J=38  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 ω_0=2π×16.0  GHz and α=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 S in Fig. <ref>(b). 

apsrev_longbib



