A general characterization method for nonlinearities in superconducting circuits

Detailed knowledge of nonlinearity in superconducting microwave circuits is required for the optimal control of their quantum state. We present a general method to precisely characterize this nonlinearity to very high order. Our method is based on intermodulation spectroscopy at microwave frequencies and does not require DC-connection or DC-measurement of an on-chip reference structure. We give a theoretical derivation of the method and we validate it by reconstructing a known nonlinearity from simulated data. We experimentally demonstrate the reconstruction of the unknown nonlinear current-phase relation of a microwave resonator with superconducting nanowires.


Introduction
Superconducting microwave circuits with Josephson junctions are a promising platform for quantum technologies [1][2][3]. With their well-controlled and lossless nonlinearity they are useful for the manipulation of single quanta of microwave radiation [4][5][6], creation of non-classical states of the electromagnetic field, and amplification of signals at the quantum limit of added noise [7][8][9][10][11]. Recent proposals and experimental realizations with Schrödinger-cat states of the harmonic oscillator have demonstrated optimal control in the quantum regime [12][13][14][15]. The complexity of the circuits and their control schemes is rapidly growing, providing incentive for accurate determination of the parameters of their Hamiltonian [3]. In particular, advanced pulseshaping and optimal control algorithms [16][17][18] can suffer from loss of accuracy due to inaccurate approximation of the nonlinear terms (i.e. truncation to second order in a Taylor expansion). Higher order nonlinear terms also affect the dynamic range of parametric amplifiers [11,19] and oscillators [20,21]. In principle, the high order nonlinear terms can be calculated if the individual circuit elements are accurately known [22][23][24][25] e.g. each Josephson junction has a known critical current. However, while it may be safe to assume that an individual Josephson tunnel junction has a sinusoidal current-phase relation, a multi-junction circuit [9,[26][27][28] or high kinetic inductance material display nonlinearities significantly different from that of the single tunnel junction [10,[29][30][31][32][33]. In this work we reconstruct the coefficients of a polynomial model of the nonlinearity from a measurement of intermodulation in the nonlinear circuit. A multifrequency lock-in technique is used to capture high order frequency mixing, or intermodulation between carefully chosen drive tones. Through analysis of the intermodulation spectrum we determine the model coefficients using methods of linear algebra.
Intermodulation is a nonlinear effect which moves signal power between the Fourier components of a system's response, something that does not happen in linear systems. This usually unwanted phenomenon causes signals to appear in frequency bands where they are not wanted. Intermodulation can distort control pulse sequences in a quantum computer and might become a serious obstacle to the parallelization of quantum gates with superconducting circuits. Thus, methods to accurately characterize intermodulation distortion are needed, and strategies to circumvent its effect must be developed. Our work is a step in this general direction.
The paper is structured as follows: in section 2 we introduce the theory of our method for reconstructing the nonlinear current-phase relation of a Josephson junction circuit based on the numerical pseudo-inverse of a Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. matrix constructed from the measured intermodulation spectrum. In section 3 we demonstrate the performance and limitation of the matrix inversion technique by reconstructing a known nonlinearity from simulated data. Reconstruction of an unknown nonlinearity from measured data on an array of superconducting nanowires is then demonstrated in section 4.

Method to reconstruct the current-phase relation from an intermodulation spectrum
Our reconstruction method is based on previous methods developed for nonlinear mechanical systems [34][35][36]. It can be summarized in the following way: we Fourier transform the equation of motion arriving at a set of equations in the frequency domain that can be written in matrix form. Numerical pseudo-inverse of a matrix then allows one to find both linear and nonlinear circuit parameters.
Consider a nonlinear circuit element with unknown energy-phase relation E(j) embedded in a microwave resonator as depicted in the schematic figure 1. We require that the energy-phase relation E(j) has a global minimum at j=0. This will insure stable orbits of the systems response for any applied drives and exclude instabilities like Josephson junctions switching into the 'voltage state'. It is convenient to express the circuit energy in terms of the flux Φ, which is related to the phase via j/2π=Φ/Φ 0 ≡ X.
We consider a microwave drive signal V drive (t) applied to the left of the input capacitor C in . Both the drive and the response are further assumed to be periodic with the same period T. This assumption necessarily restricts our reconstruction method to regime where the response is not chaotic. Furthermore with appropriate choice of the drive frequencies and integration time T int =nT with integer n, all mixing products in the nonlinear response are contained in the set of discrete frequencies ω j =jΔω with Δω=2π/T int . In principle, sub harmonic response can also be captured, e.g. period doubling when n=2, but here we will restrict ourselves to periodic response.
The equation of motion for the node flux is where the total damping rate γ, containing both internal losses as well as losses through the input and output capacitors, and the input coupling constant κ in , are defined in appendix A. The nonlinearity in equation (1) is the current-flux relation I(Φ)=∂E(Φ)/∂Φ. Fourier transforming equation (1) gives the following set of equations valid at the discrete frequencies ω j =jΔω.
Here j  denotes the discrete Fourier transform and the discrete frequency spectrum of the dimensionless node flux is denoted with the complex numbers X ĵ . Similarly, the drive spectrum is denoted as D V j j The resonance frequency ω 0 is defined as  where the jth row of the matrix H jl is given by and p l are the components of a column vector of circuit parameters All that is left to do is to invert the matrix H to find the circuit parameters including the coefficients of the polynomial current-flux relation In the experiment we measure the response in the frequency domain obtaining a comb of M intermodulation products (IMPs) X ĵ near ω 0 . Using the inverse Fourier transform we obtain X t X , where the ≅ sign is due to the fact that we use a partial spectrum, neglecting response far from resonance. We then Fourier transform products of X(t) to build the matrix H. Due to the Fast Fourier Transform algorithm, this procedure is computationally more efficient than performing successive convolutions of the signal X ĵ with itself in the frequency domain. To reconstruct N polynomial coefficients, the numerical pseudo-inverse H −1 requires MN. In practice we find M2N is necessary for stable inversion.
The above procedure gives the coefficients of an optimal polynomial approximation, in a least-square sense, of the current-phase relation over the interval of oscillation probed by the given drive. Here we treat the linear circuit parameters ω 0 and γ as unknown. If they are known a priori the corresponding terms in equation (2) can be brought over to the other side and taken out of the vector p. One might question whether a sum of monomials is the best choice for this approximation. In appendix C we discuss reconstruction using orthogonal polynomials, showing that they do not provide a better approximation.
In this work we focus on a system fully described by two degrees of freedom, the phase j and its derivative j . Systems containing additional (external or internal) degree's of freedom, such as resonators coupled to a qubit, require an extension of the model presented here and are beyond the scope of this work.

Performance of reconstruction based on matrix inversion
In order to demonstrate the utility and capacity of the method we simulate spectral data from a known nonlinear current-flux relation and use this data to reconstruct the nonlinearity. We assume a Josephson tunnel junction with a sinusoidal current-phase relation I I sin 2 in parallel with a linear inductor L  . The resulting potential E(Φ) is a parabola superposed with a cosine. This element is embedded in the circuit figure 1 with parameters listed in the table. Equation (2) can be written as two coupled first order differential equations is defined through the total inductance which for this circuit is We numerically integrate equation (2) using scipy.integrate. odeint from the SciPy library [37], a python wrapper for LSODA from ODEPACK [38].
We apply a two-tone drive with equal drive amplitudes at frequencies f 1 =0.996 ω 0 /(2π) and f 2 =0.998 ω 0 /(2π). The resulting X(t) is plotted in figure 2(b) (orange) together with the applied drive (dark green). The node flux undergoes a fast oscillation enveloped by a slow amplitude and phase modulation or 'beat'. The Fourier spectrum which is depicted in figure 2(a) shows response not only at the drive frequencies but also at integer linear combinations of the two drive tones. We take a section of the intermodulation spectrum around resonance to build up the matrix H defined in equation (5). Using a total of M=30 IMPs and a polynomial expansion up to N=15 coefficients, we reconstruct the energy-phase relation. The reconstructed potential (blue) and the original potential (black dashed) are plotted in figure 2(c). The area highlighted in orange shows the region of X which is actually probed by the applied drive. Within this region we find excellent agreement with the original potential, but extrapolating outside this region shows substantial deviation.
In order to test the reconstruction method for a strongly nonlinear system, we use the same potential as above but we increase the drive voltage by a factor of 1000. The simulated time evolution of the strongly nonlinear response is plotted in figure 3(b). The beat pattern is distorted by the nonlinearity and the corresponding frequency spectrum shows a large number of IMPs. Using 500 IMPs, centered around the drives to reconstruct the simulated potential, we see in figure 3(c) that even for this strongly nonlinear case the reconstructed potential (blue) agrees well with the simulated potential (black dashed) in the probed region.
In order to quantify the accuracy of reconstruction and find a criterion for the number of IMPs needed, we take the simulated data from figure 3 and perform the reconstruction using a different number of IMPs around  resonance. The inverse Fourier transform used to obtain the time evolution is performed on a truncated, or windowed intermodulation spectrum. A square window function causes ripples in the time domain (Gibbs phenomenon), which in turn cause errors in the reconstruction. These errors can be reduced by letting the IMPs smoothly decay to zero at the edges of the frequency window. We use a window function given by sinc where f c is the center frequency, Δ f the width of the window and p controls the strength of the smoothing. Because larger p suppresses higher order IMPs, the reconstruction improves only if the error introduced by smoothing is small compared to the error caused by truncation.
In order to quantify this error, we integrate the mean-square deviation between the reconstructed and original potentials over the region probed by the beating oscillation In figure 4(a) we plot the error ò as a function of the number M of IMPs used for the reconstruction. We show three different values of the strength p of smoothing. The corresponding window functions are depicted in figure 4(b). In all three cases the error falls sharply with only a few IMPs before saturating for larger M. With a square window (blue solid line) the mean square deviation shows large oscillation with M. These oscillations are suppressed by smoothing the window function (increasing p), while at the same time increasing the error. Figure 4(c) shows the mean relative error as a function of the response power fraction (rpf) or the fraction of the total response power which is in the IMPs used for the reconstruction. We see that the best reconstruction has rpf 1 dB =or approx 80 % of the total response power. This can be used to estimate the spectral region needed for a reliable reconstruction. The dependence of the rpf on the number of IMPs M is plotted in figure 4(d).
In the case demonstrated in figure 4 the Josephson junction is driven above its critical current, resulting in extremely strong intermodulation and making it difficult to capture all the response power. However, even for much lower drive powers where only a few IMPs are observed, the rpf never reaches 0 dB because the nonlinearity not only transfers power from the drives to the IMPs, but also to response at higher harmonics. Taking only the response around resonance limits the resolution of our reconstruction method. Including the response at higher harmonics or subharmonics is in principle possible but difficult to realize in experiment as it would require large bandwidth or additional up-and down-conversion circuits. Using the effect of parametric down-conversion it should be possible to also reconstruct even order terms in the expansion of the currentphase relation by applying the two-tone drive at twice the resonance frequency and measuring the response around resonance.

Reconstruction from measured data
As an example of reconstruction with experimental data we determine the current-phase relation of an array of superconducting niobium nanowires incorporated in the center conductor of a coplanar waveguide (CPW) resonator. The sample and measurement are described in appendix C. When measuring transmission, the Fourier amplitudes of the response B ĵ and drive amplitudes E ĵ are related to the amplitude inside the resonator as where G is the total power gain of the amplification and attenuation in the output line, κ out is the coupling strength through the output coupling capacitor, α is the attenuation of the input line and δ and β are phases picked up by the signal while propagating through the input and output lines, respectively. Equation We now build the matrix H according to equation (5) using B ĵ instead of X ĵ , the matrix inversion equation (7) yields the following parameter vector The phase factor in front of the damping term can be determined by numerically rotating the drive by 90°and then taking the quotient of the first entry of p T for the rotated and unrotated drive. The nonlinear coefficient however contain powers of the gain and the output coupling constant which cannot be eliminated from the nonlinear terms. Thus a precise quantitative characterization of an unknown nonlinear system requires prior knowledge of the κ out and G. In this experimental example we use a narrow frequency band of 200 kHz for the reconstruction. For such a narrow band we can consider the gain and the phase delay in the cables to be frequency independent. If wider bands are used for reconstruction a frequency dependent calibration of gain and phases might be necessary. For our experiment we use thermal noise to calibrate the gain and the output coupling coefficient estimated from the sample geometry. The values we extract from the reconstruction for α agree with the estimated attenuation of the input line. The reconstructed current-phase relation for the nonlinear superconducting resonator is shown in figure 5(a). In this experiment the current-phase relation is only weakly nonlinear in the region probed by the drive. In order to make the nonlinearity clearly visible we subtract the linear part of the current-phase relation in 5(b).

Conclusions
We presented a general method to reconstruct nonlinearities in microwave circuits from an intermodulation spectrum. The method to finds the coefficients of a polynomial approximation of the nonlinear potential by numerical pseudo-inverse of an intermodulation matrix. Using numerically simulated spectra we showed that the method faithfully reconstructs the original potential in the weak and strongly nonlinear regimes, provided that sufficient spectral data is available. We found that the truncation of the intermodulation spectrum causes ripples in the reconstruction fidelity.
In addition we reconstructed the current-phase relation of a coplanar niobium microwave resonator rendered nonlinear by the introduction of nanowire constrictions in the center conductor. The intermodulation spectrum was measured using a multifrequency lock-in measurement setup. A quantitative reconstruction from measured data requires prior knowledge of the gain of the amplification chain used to measure the intermodulation spectrum. The method we present here is general and easy to implement experimentally, with the potential to become a versatile tool for the characterization of nonlinearity in superconducting microwave circuits. The matrix inversion method can be in principle be extended to also include even order nonlinear coefficients as well as nonlinear damping mechanisms.
) respectively. Using the Euler-Lagrange formalism one derives the following equations of motion for the node-fluxes In the experiment we will apply a periodic drive with a period of T=(2π)/(Δω) and so we can assume that the system's response will also be periodic on the time interval T. In this case we can express the drive and system response as a discrete sum of Fourier components with frequencies ω j =jΔω. X e A . 5 where D j is the complex drive amplitude at frequency ω j . B ĵ is the output amplitude. The relation between A B , j ĵâ nd D ĵ and the amplitude X j inside the resonator are yet to be determined. Inserting into equations Around resonance we can neglect the frequency dependence of the coupling so that In practice, the phases j in and j out can be included in the input and output phases β and δ and be set to zero in A. 11

Appendix B. Approximation of the current-phase relation with orthogonal polynomials
In the main body of the paper we use a sum of monomials to approximate the nonlinear current-phase relation I (Φ). Monomials have the drawback that the different orders are nonorthogonal. That means that adding or removing higher order terms can alter the coefficients of the lower order terms To circumvent this problem, one can instead use a set of orthogonal polynomials for the matrix inversion.
As an example we show in figure B1 the reconstruction of the simulated potential from section 3 using Hermite and Chebyshev polynomials. Similar to the reconstruction in section 3 we use only the odd order polynomials for the current-phase relation I(Φ). Including the even order polynomials introduces large errors as the spectrum around resonance, which contains only odd order IMPs, does not contain any information about the even order nonlinear components.
In figure B1 the mean relative deviation between reconstructed and simulated potential E(Φ) is plotted as a function of the polynomial order. The three curves correspond to the sum of monomials (blue solid), Hermite polynomials (orange dashed) and Chebyshev polynomials (green dashed-dotted). For approximations up to 13th order the three polynomial approximations show identical reconstruction fidelity, demonstrating that the reconstruction is in fact independent of the choice of basis functions. For higher polynomial orders the reconstruction with Hermite polynomials becomes numerical unstable due to the large prefactors of the higher order 'physicist's' Hermite polynomials. nanowires is shown in figure C1. Nineteen wires with a length of 200 nm and a diameter of 50 nm were inserted. The intermodulation spectrum was measured using a multifrequency lock-in amplifier (MLA) in combination with an up-and down-conversion circuit schematically depicted in figure C1. The multifrequency drive from the MLA is upconverted using an IQ-mixer and a local oscillator signal at around 5.5 GHz. The drive amplitudes at the I and Q port of the upconversion IQ-mixer are adjusted such that only the upper sidebands appear at the output of the mixer while the lower sidebands are canceled. This upper sidebands are placed near resonance of the CPW resonator. The signal transmitted through the sample passes a cryogenic isolator, is amplified with a cryogenic amplifier and then down-converted using the same local oscillator. We demodulate the signal at 21 frequencies centered around an IF-frequency of 6 MHz and we use a digital image cancelation technique to capture only response from the high sidebands.
The gain of the amplification chain has been previously calibrated using the Johnson-Nyquist noise of the 50 Ω termination of the cryogenic circulator. Off resonance the transmission line resonator shows almost  perfect reflection, so that the thermal noise of the resistor is reflected into the amplification chain. Heating the dilution refrigerator from base-temperature up to 800 mK the temperature dependence of the noise can be fitted and the gain and added noise of the amplification chain can be extracted.