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 non-ideal 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 multilayer 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*) 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.


I. 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, KPUPs, and frequency-tuneable superconducting resonators, this nonlinear kinetic inductance is critical to the operation and performance of the devices [2, 3, 10-13]; in other cases, the nonlinear kinetic inductance results in non-ideal behaviour that is important even in common device operating environments [9]. As such, understanding and calculation of the nonlinear kinetic inductance is 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 BardeenCooperSchrieffer (BCS) theory of superconductivity [14][15][16]. 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 [15]. The theory and experiment demonstrated excellent agreement, lending confidence to the use of the Usadel equations as the foundation of our analysis framework. The paper by Anthore et al. 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,17]. 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 [18], we compute the complex conductivities of the thin-films. We then compute the surface impedances using the transfer matrix method [19]. Finally, we calculate the transmission line inductances from the surface impedances by using the appropriate transmission line theory for the geometry of the device [20], 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.

II. THEORY
A. Usadel equations In this analysis, the multilayers are stacked in the x direction, and the supercurrent flows in the z direction. The Usadel equations in one dimension are [14,15,[21][22][23] hD S 2 and where θ is a complex variable dependent on energy E parametrising the superconducting properties, N S is the electron single spin density of states, V 0,S is the superconductor interaction potential, ∆ is the superconductor 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 [24], 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 self-consistency 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. [15] 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 equation (2) can be removed, simplifying equation (2) into iE sin θ + ∆(x) cos θ − Γ cos θ sin θ = 0, where Γ =hD S /2 * (∂ φ /∂ z) 2 is the depairing factor. The above equation can be solved iteratively with equation (3) to obtain ∆(Γ). Numerically, it is easier to solve equation (5) for sin(θ ) using a polynomial root finder, rather than finding θ directly.
In the case of a multilayer 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 [19].
Instead of calculating nonlinearity with respect to Γ, which is not constant across the multilayer, calculations should be performed with respect to ∂ φ /∂ z. ∂ φ /∂ z cannot vary across the multilayer (in the x direction) due to the absence of net supercurrent (in the x direction). Computationtime-wise, it is beneficial adopt the thin-film approximation scheme that has demonstrated good agreement with experiment for multilayer superconductors. The approximation assumes θ varies slowly, and can be accounted by a second order polynomial expansion. [24,25] B

. Complex Conductivities and Impedances
Nam's equations [26] are a generalization of the Mattis-Bardeen [27] theory into strongcoupling 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 [19].
After calculating σ , the complex surface impedance for a homogeneous single layer can then be obtained using [28] where t is the thickness of the homogeneous superconducting film, and µ 0 is the vacuum permeability.
For multilayers, Z s can be found by dividing the multilayer into thin layers of thickness δ x, and then cascading the resultant transfer matrices along the multilayer. A detailed discussion of the above methodology, as well as an analysis of numerical results for Al-Ti multilayers, can be found in [19].

C. Transmission Line Properties
The series impedance and shunt admittance of a transmission line structure can be calculated from Z s as follows [20,29,30]: where k 0 is the free-space wavenumber, η 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 [31,32]. g 1 and g 2 are geometric factors which can be calculated using appropriate conformal mapping theories [20,30]. against energy E/k B at temperature T = 0.01 K for different values of Γ/∆ 0 , where ∆ 0 ≈ 1.764 k B T c 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 [15].
Previous approximations on the inductance nonlinearity [9,17]  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.
The right figure of 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, 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:

IV. CRITICAL TEMPERATURE EXPERIMENT
Many aspects of our analysis routine have been individually experimentally established by previous studies: the analysis of superconducting multilayers using the Usadel equations has been justified by [25,33]; the analysis of supercurrent using the Usadel equations has been justified by [15]; the computation of complex conductivities using Nam's equations has been justified by [18]; the calculation of transmission line properties using conformal mapping analysis has been justified by [34][35][36].
Despite the above experimental justifications, a caveat exists regarding the analysis of supercurrent using the Usadel equations: the physical dimensions of the devices tested in previous studies are smaller than the dimensions typically used in KIDs or TWPAs. For example, the aluminium strip tested by [15] has width w = 120 nm and thickness t = 40 nm; the aluminium strips tested by [37] has dimensions w = 30 − 61 nm, t = 20 − 89 nm. In particular, [37] has identified that, for devices with w > ξ at high current densities, vortex formation will result in deviations from ideal behaviour. Here ξ refers to the coherence length of the material. For superconducting strips with w on the order of a few microns [2, [38][39][40], it is useful to determine the range of current within which the Usadel equations treatment of the supercurrent is valid. 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 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 thermometer. 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.  [25]. 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. 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, [41][42][43]. 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.

V. 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 appreciated in our discussion around the middle figure of 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. * sz311@cam.ac.uk