Nonlinear Properties of Supercurrent-Carrying Single- and Multi-Layer Thin-Film Superconductors

Superconducting thin films are central to the operation of many kinds of quantum sensors and quantum computing devices: kinetic inductance detectors (KIDs), travelling-wave parametric amplifiers (TWPAs), qubits, and spin-based quantum memory elements. In all cases, the nonlinearity resulting from the supercurrent is a critical aspect of behaviour, either because it is central to the operation of the device (TWPA), or because it results in nonideal second-order effects (KID). Here, we present an analysis of supercurrent-carrying superconducting thin films that is based on the generalized Usadel equations. Our analysis framework is suitable for both homogeneous and multi-layer thin films, and can be used to calculate the resulting density of states, superconducting transition temperature, superconducting critical current, complex conductivities, complex surface impedances, transmission line propagation constants, and nonlinear kinetic inductances in the presence of supercurrent. Our analysis gives the scale of kinetic inductance nonlinearity (I∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_*$$\end{document}) for a given material combination and geometry, and is important in optimizing the design of detectors and amplifiers in terms of materials, geometries, and dimensions. To investigate the validity of our analysis across a wide range of supercurrent, we have measured the transition temperatures of superconducting thin films as a function of DC supercurrent. These measurements show good agreement with our theoretical predictions in the experimentally relevant range of current values.


Introduction
Owing to their low-loss, high-quality-factor characteristics below their superconducting transition temperatures (T c ), superconducting thin films are important to the operation of many kinds of quantum sensors and quantum computing devices, such as kinetic inductance detectors (KIDs) [1], travelling-wave parametric amplifiers (TWPAs) [2], kinetic inductance parametric up-converters (KPUPs) [3], superconducting qubits [4], and spin-based quantum memory elements [5,6]. When designing these superconducting devices, an important consideration is the nonlinearity in superconducting kinetic inductance with respect to supercurrent [7,8]. The nonlinear inductance of a superconducting device is expected to have the form [9] L = L 0 1 + where L is the inductance of the device, L 0 is the inductance in the absence of supercurrent, I is the supercurrent, and I * is the scale of the quadratic inductance nonlinearity.
In the case of TWPAs and KPUPs, this nonlinear kinetic inductance is critical to the operation and performance of the devices [2,3,[10][11][12]; in other cases, the nonlinear kinetic inductance results in nonideal behaviour that is important even in common device operation power regimes, which often involve high readout power in order to improve noise performance [9]. As such, understanding and calculation of the nonlinear kinetic inductance are important to the quantitative design processes of these thin-film devices.
Analyses of supercurrent in superconducting thin films can be based on the Usadel equations, which is a set of diffusive limit equations derived from the Bardeen-Cooper-Schrieffer (BCS) theory of superconductivity [19][20][21]. Anthore et al. have calculated and experimentally measured the resultant density of states in a superconducting thin film due to supercurrent using the Usadel equations [20]. The theory and experiment demonstrated excellent agreement, lending confidence to the use of the Usadel equations as the foundation of our analysis framework. Further, the work by Clem et al. [21] based on the Usadel equations has been applied experimentally to estimate the depairing current of superconducting nanowires to good agreement [22]. The paper by Anthore et al. in particular presents a series expansion of the superconducting order parameter (Δ) with respect to supercurrent for single-layer superconducting thin films. This series expansion has been used by other studies to estimate the superconductor complex conductivities and kinetic inductances [9,23]. As we shall demonstrate in this study, this approximate approach does not account for the change in the shape of the density of states and underestimates the impact of supercurrent.
Using the full density of states as an input to Nam's equations [24], we compute the complex conductivities of the thin films. We then compute the surface impedances using the transfer matrix method [25]. Finally, we calculate the transmission line induc-tances from the surface impedances by using the appropriate transmission line theory for the geometry of the device [26], such microstrip transmission line or coplanar waveguide.
We have also measured the supercurrent dependence of the superconducting transition temperatures for single-layer titanium (Ti) and multi-layer aluminium-titanium (Al-Ti) thin films. Our results confirm the validity of the Usadel theory approach for experimentally realistic device dimensions and current regimes.

Usadel Equations
In this analysis, the multi-layers are stacked in the x direction, and the supercurrent flows in the z direction. The Usadel equations in one dimension are [19,20,[27][28][29] and where θ is a complex variable dependent on energy E parametrizing the superconducting properties, N S is the electron single spin density of states, V 0,S is the superconductor interaction potential, Δ is the superconducting order parameter, k B Θ D,S is the Debye energy, k B is the Boltzmann constant, T is the temperature of the superconducting film, D S is the diffusivity constant, given by D S = σ N /(N S e 2 ) [30], e is the elementary charge, Im(x) takes the imaginary part of x, and finally, σ N is the normal state conductivity, at T just above T c . Equation (3) is the selfconsistency equation for order parameter Δ. We have introduced the superfluid velocity where φ is the superconducting phase and #» A is the magnetic vector potential. We assume that the effect due to the induced field is negligible compared to that of supercurrent [20].
The supercurrent density #» j is given by For supercurrent flowing in the z direction, #» v 2 s = D 2 S (∂φ/∂z) 2 . For the case of a homogeneous BCS superconductor, the first term of Eq. (2) can be removed, simplifying Eq. (2) into where Γ = D S /2 * (∂φ/∂z) 2 is the depairing factor. The above equation can be solved iteratively with Eq. (3) to obtain Δ(Γ ). Numerically, it is easier to solve Eq. (5) for sin(θ ) using a polynomial root finder, rather than finding θ directly.
In the case of a multi-layer superconductor, the boundary conditions (BCs) between the layers need to be taken into account. The BCs suitable for the Usadel equations can be found in [25]. Instead of calculating nonlinearity with respect to Γ , which is not constant across the multi-layer, calculations should be performed with respect to ∂φ/∂z. ∂φ/∂z cannot vary across the multi-layer (in the x direction) due to the absence of net supercurrent (in the x direction). Computation time-wise, it is beneficial to adopt the thin-film approximation scheme that has demonstrated good agreement with experiment for multi-layer superconductors. The approximation assumes θ varies slowly and can be accounted by a second-order polynomial expansion [30,31].

Complex Conductivities and Impedances
Nam's equations [32] are a generalization of the Mattis-Bardeen [33] theory into strong coupling and impure superconductors. Nam's equations compute the complex conductivity σ = σ 1 − jσ 2 using a pair of integrals of θ across energy E. The integrals as well as their evaluations for Al-Ti bilayers can be can be found in [25].
After calculating σ , the complex surface impedance for a homogeneous single layer can then be obtained using [34] where t is the thickness of the homogeneous superconducting film and μ 0 is the vacuum permeability. For multi-layers, Z s can be found by dividing the multi-layer into thin layers of thickness δx and then cascading the resultant transfer matrices along the multi-layer. A detailed discussion of the above methodology, as well as an analysis of numerical results for Al-Ti multi-layers, can be found in [25].

Transmission Line Properties
The series impedance and shunt admittance of a transmission line structure can be calculated from Z s as follows [26,35,36]: where k 0 is the free-space wave number, η 0 is the impedance of free space, subscript n identifies superconductor surfaces, which are upper, lower, and ground surfaces, denoted by subscripts u, l, and g, respectively, f m is the effective modal dielectric constant, which is given by existing normal conductor transmission line theories, for example [37,38]. g 1 and g 2 are geometric factors which can be calculated using appropriate conformal mapping theories [26,36]. After obtaining the series impedance and the shunt admittance, other properties of the superconducting transmission line can be calculated straightforwardly. The characteristic impedance is given by η = (Z /Y ) 1/2 . The propagation constant is given by γ = α + jβ = (ZY ) 1/2 , where α is the attenuation constant and β is the phase constant. The inductance per unit length L can finally be calculated using L = Im(Z )/ω. The calculation can then be iterated for different values of I to obtain L(I ), which allows the extraction of I * using a polynomial fit.

Results and Discussion
The left figure in Fig. 1 shows Ti superconducting density of states (DoS) is the superconducting energy gap of Ti in the absence of supercurrent. The presence of supercurrent broadens the DoS. This is a real effect and it has been experimentally observed by [20]. Previous approximations on the inductance nonlinearity [9,23] assume the effect of this new DoS on conductivity can be approximated by a single parameter Δ. Effectively, these studies have assumed that the new DoS can be approximated by a zero-supercurrent DoS shifted to an altered DoS gap at Δ. For convenience, we label this simplified DoS function as n j=0 [Δ(Γ )]. As we see in the middle figure, this assumption leads to underestimation on the impact of the supercurrent.
The middle figure in Fig. 1 shows a plot of normalized reactive conductivity σ 2 /σ N against Γ /k B for Ti at T = 0.01 K, frequency f = 10 GHz. The red line shows calculation performed by solving Nam's equations using the full densities of states shown in left figure. The blue line shows calculation performed using n j=0 [Δ(Γ )]. The black line shows calculation performed using n j=0 [Δ g (Γ )], where Δ g is the energy at which the broadened DoS becomes nonzero. Since ω Δ 0 , the blue and black lines have approximate forms σ 2 /σ N = (π Δ)/( ω) and σ 2 /σ N = (π Δ g )/( ω), respectively. Comparing the red line with the blue line, we notice that approximation using n j=0 [Δ(Γ )] underestimates the effect of supercurrent. This shows that the broadened DoS in the presence of supercurrent cannot be approximated well using a single energy parameter Δ(Γ ). Comparing the red line with the black line, approximation using n j=0 [Δ g (Γ )] overestimates the effect of supercurrent. This is because, in the presence of supercurrent, the DoS is broadened. As a result, Δ g shifts further than the overall DoS. The above results highlight the need to perform the full calculation as detailed in this manuscript. For practical purposes, we give the approximation of σ 2 /σ N as a function of Γ /Δ 0 , valid for k B T Δ 0 , ω 2Δ 0 , and Γ /Δ 0 < 0.2: This expansion can be used in conjunction with previous results from [20], which states that the supercurrent is given by , S is the cross-sectional area, ξ is the superconducting coherence length, e is the electron charge, and U S can be approximated by: The right figure in Fig. 1 shows a plot of inductance per unit length L against squared supercurrent I 2 for a Ti microstrip line with thickness t = 100 nm, width w = 5 µm, dielectric height h = 250 nm, and ground plane Ti thickness t g = 200 nm. We see from the figure that L can be approximated well by a quadratic expansion on I at small current values. At larger values, an additional quartic term is needed to encapsulate the superconductor response: where I * ,4 is the scale of the quartic order of inductance nonlinearity. The Ti microstrip studied here has I * = 8.5 mA and I * ,4 = 5.5 mA. Inset of the right figure in Fig. 1 shows a plot of I * against Al thickness t Al for a bilayer Al-Ti microstrip with fixed Ti thickness t Ti = 100 nm, width w = 5 µm, dielectric height h = 250 nm, and ground plane Ti thickness t g = 200 nm. As t Al increases, the nonlinear behaviour of the microstrip decreases in significance: this is reflected in the higher I * values. This trend agrees with our expectations: the presence of an Al layer decreases the resistivity of the multi-layer. This lower resistivity in turn results in smaller nonlinearity [2,23].

Critical Temperature Experiment
Many aspects of our analysis routine have been individually experimentally established by previous studies: the analysis of superconducting multi-layers using the Usadel equations has been justified by [31,39]; the analysis of supercurrent using the Usadel equations has been justified by [20]; the computation of complex conductivities using Nam's equations has been justified by [24]; and the calculation of transmission line properties using conformal mapping analysis has been justified by [40][41][42].
Despite the above experimental justifications, a caveat exists regarding the analysis of real superconducting devices using the Usadel equations: the physical dimensions of the devices tested in previous studies have physical dimensions smaller than the relevant length scales of the material system, i.e. the perpendicular field penetration depth λ ⊥ and the superconducting coherence length ξ as identified by [43]. To illustrate, the aluminium strip tested by [20] has width w = 120 nm and thickness t = 40 nm; the aluminium strips tested by [43] has dimensions w = 30−61 nm and t = 20−89 nm. These dimensions are much smaller than those typically used in the design of KIDs and TWPAs, and many real devices have dimensions comparable to, or exceeding, one or both of the length scales. For a thin film with w > λ ⊥ , the supercurrent distribution becomes nonuniform as current piles up near the edges [44]; for a thin film with w > 4.4ξ , at high current densities, vortex formation will result in deviations from ideal behaviour [45]. The thin-film parallel field penetration depth is given by λ ≈ λ L √ ξ 0 /l, the thin-film perpendicular field penetration depth is given by λ ⊥ = λ 2 /t, and the low-temperature coherence length is given by ξ ≈ √ ξ 0 l where ξ 0 = v F /(π Δ 0 ) is the bulk coherence length, l is the mean free path, v F is the Fermi velocity, λ L = m e /(μ 0 ne 2 ) is the London depth, m e is the electron mass, and n is the electron density [20,43,44]. For our 25 nm aluminium thin films, λ ⊥,Al = 0.40 µm and ξ Al = 0.20 µm; for our 100 nm titanium thin films, λ ⊥,Ti = 0.12 µm and ξ Ti = 0.56 µm. Here, we have used data for Al and Ti properties from [31], supplemented by mean free path data from [46,47]. For superconducting strips with w on the order of a few microns [2,14,17,[48][49][50] on the border of the relevant length scales, it is useful to determine the range of current within which the 1D Usadel equation treatment of the supercurrent provides a good prediction of device behaviour. To this end, we have performed an experiment measuring the T c of a superconducting strip for a given supercurrent I .
Ti and Al-Ti films were deposited by DC magnetron sputtering at a base pressure of 2 × 10 −10 Torr or below. For bilayer films, Al layers were deposited after Ti layers without breaking the vacuum. The films were patterned to achieve four-terminal sensing geometry and connected to electronics via Al wire bonds. The samples were mounted to the cold stage of a dilution refrigerator inside a niobium magnetic shield. Temperature monitoring was achieved using a calibrated ruthenium oxide thermome- ter. For each set of measurement, a fixed current was first injected to the mounted superconducting film. The temperature of the cold stage was then slowly raised until transition from superconducting to normal state had occurred. The potential difference across the film was continuously measured throughout this transition process.
The left figure in Fig. 2 shows a plot of scaled current density j/ j 0 and scaled superconducting order parameter Δ/Δ 0 against scaled depairing factor Γ /Δ 0 for temperature T = 0.01 K. There exist a maximum current density j c . The pair of values ( j c , T ) marks out a curve on the phase diagram within which the material is in the superconducting state, and beyond which the material is in the normal state.
At T ≈ 0 K, j c ≈ 0.746 j 0 , where j 0 = N S σ N Δ 3 0 / . It is worth noting that when j = 0, transition to normal state happens when Δ is nonzero. Computationally, this means that the small Δ, θ approximation technique, commonly used to compute the T c of j = 0 transitions [30,31], cannot be applied for these j = 0 transitions. The middle figure in Fig. 2 shows a plot of critical current in reduced units [I /I 0 ] 2/3 against critical temperature in reduced units T c /T c,0 for a Ti strip with thickness t = 100 nm. The right figure in Fig. 2 shows a similar plot for Al-Ti bilayers with Al thickness t Al = 25 nm and Ti thickness t Ti = 100 nm. The y-axis is chosen to reflect the Ginzburg-Landau result in the small supercurrent limit [43] which states that I /I 0 ∝ 1 − T c /T c,0 3/2 . For both plots, the dotted line shows the values obtained from theoretical calculations using the Usadel equations; the scattered markers show the experimentally measured values for different widths of superconducting lines. The physical parameters used to generate the theoretical lines are the same as those used in [31]. To convert from j to I , we have used I = jtw, where the thickness t is deduced from calibrated deposition time and the width w is part of the design of the deposition mask. As expected from the above analysis of length scales, within each plot, wider superconducting lines result in earlier deviation from the ideal theoretical calculations. Denote I c,0 as the actual critical current of a device at close to 0 K (not to be confused with I 0 which is the theoretical critical current). For most devices, the experimental data demonstrate good agreement with the theoretical prediction at I < I c,0 /2. For the widest bilayer device, good agreement is still obtain at I < I c,0 /3. This range encapsulates the common operating current values for typical TWPAs and KIDs systems: current much smaller than I c,0 is usually chosen to avoid the onset of high current dissipation, or to avoid resonator bifurcation [2,[51][52][53]. In this study, we have chosen a conservative thickness of 100 nm. We expect an even bigger range of agreement for thinner devices such as the coplanar waveguides studied in [2], which have thickness t = 35 nm.

Conclusions
We have presented a numerical routine for analysing the inductance nonlinearity of thin-film superconductors with respect to supercurrent. Our analysis routine is based on the Usadel equations, Nam's equations for complex conductivity, transfer matrix calculation for complex surface impedances, and transmission line models. As mentioned in our discussion around the middle figure in Fig. 1, our analysis takes into account the full shape of superconducting densities of states and avoids an underestimation made by previous analyses on this subject. We have measured the superconducting transition temperature as a function of supercurrent for Ti single layers and Al-Ti bilayers.
Our results show that the theory is in agreement with the experimental data in the current range that most thin-film superconductor devices are operated at and therefore allows this analysis to be integrated in the design and optimization of future thin-film superconducting devices. Care needs to be taken when applying the numerical routine to AC applications, as both AC current distribution [26,36] and field quantization (in particular coherent excited states) [40,54] effects are important at frequencies comparable to the superconductor pair-breaking frequency, and are likely to result deviation from the DC treatment in the Usadel equation formalism. Future studies should be conducted to investigate the extent of applicability as well as techniques to adapt the routine to AC applications.