Phonon heat capacity and self-heating normal domains in NbTiN nanostrips

Self-heating normal domains in thin superconducting NbTiN nanostrips with the granular structure were characterized via steady-state hysteretic current–voltage characteristics measured at different substrate temperatures. The temperature dependence and the magnitude of the current, which sustains a domain in equilibrium at different voltages, can only be explained with a phonon heat capacity noticeably less than expected for 3D Debye phonons. This reduced heat capacity coincides with the value obtained earlier from magnetoconductance and photoresponse studies of the same films. The rate of heat flow from electrons at a temperature Te to phonons in the substrate at a temperature TB is proportional to (Tep−TBp) with the exponent p ≈ 3, which differs from the exponents for heat flows mediated by the electron–phonon interaction or by escaping of 3D Debye phonons via the film/substrate interface. We attribute both findings to the effect of grains on the phonon spectrum of thin NbTiN films. Our findings are significant for understanding the thermal transport in superconducting devices exploiting thin granular films.


I. INTRODUCTION
Low-dimensional superconducting structures (e.g., thin films, nanowires, nanotubes, nanoparticles, and superlattices) have become building blocks for various fascinating applications such as single-photon detectors, hotelectron bolometers, kinetic inductance detectors, quantum interference devices, quantum bits (qubits), and other circuit elements [1][2][3][4][5][6].For the design and optimization of these elements, the knowledge of transport and thermodynamic properties is crucial while qualitative understanding of the impact of reduced dimensionality on these properties is of fundamental interest.Since, at sufficiently low temperatures, the phonon wavevectors perpendicular to the film plane become restricted by the reciprocal film thickness, the phonon spectrum undergoes strong modifications that manifest in a considerable change of the thermal conductivity [7] and heat capacity [8] of phonons.Although the size effect in the phonon heat capacity has been studied theoretically (for films [8], nanowires [9], and spherical grains [10,11]), the available models are limited to crystalline specimens.Most superconducting devices, however, utilize amorphous and polycrystalline, granular materials.
We have recently reported [12,13] a reduction of phonon heat capacities in thin polycrystalline granular films at low temperatures that was attributed to the effect of the mean grain size.This result was derived from extensive studies of the magnetoconductance and photoresponse of thin NbTiN, NbN, and WSi films.To support these findings and to verify the consistency of techniques in the steady state and in the time domain, here we analyze hysteretic current-voltage characteristics (CVCs) of superconducting NbTiN strips.Hysteresis appears in the regime of current bias as the difference between the experimental critical (switching) current I C and the return current I r , at which the strip returns to the superconducting state when the bias current in the normal state is gradually decreased.In the voltage-bias regime, I r defines the current plateau which is commonly affiliated to a self-heating normal domain with the length controlled by the applied voltage.The domain remains in equilibrium, which is set by the balance between Joule heating of electrons by I r and their cooling via heat diffusion along the strip and heat flow through its interfaces.
The model of the self-heating normal domain describing hysteresis was first proposed in [14].For small differences between the superconducting transition temperature, T C , and the bath (substrate) temperature, T B , the authors admitted that the heat flow from the strip to the underlying substrate is Q ∝ (T −T B ) where T is the strip temperature.In order to cover larger differences between T C and T B , this model was further modified by an introduction of the heat flow Q = K(T p −T p B ), where K is the effective thermal conductance.The approach with p = 4 was implemented in [15].This value of the exponent p is provided by a microscopic model [16] for the heat flow between two solids across their interface via 3-d Debye phonons.Further modifications of the self-heating normal domain model were made by incorporating the state (normal/superconducting) and temperature dependence of the electron thermal conductivity [17][18][19][20] or by varying the exponent p [21,22].Microscopic models [16,23,24] describing the heat flow from electrons to phonons in the film and further to phonons in the substrate show that the exponent p is not necessary an integer and may have any value from 4 to 6.The value of p in particular film is controlled by the dimensionality of phonons and by the degree of disorder.
In this study, we analyze hysteretic CVCs of granular disordered NbTiN strips with different thicknesses.We apply a modified model of the self-heating normal domain with an arbitrary exponent p in order to account for the effect of disorder and the effect of the mean grain size on the heat flow.Furthermore, we considered the heat transfer between three different systems (electrons and phonons in the film and phonons in the substrate) and address the microscopic meaning of the exponent p and the effective thermal conductance.

II. EXPERIMENT AND RESULTS
The strips were fabricated from two NbTiN films with thicknesses 6 and 9 nm, which were studied earlier in [13].The films were deposited on Si substrates on top of a 270 nm-thick thermally-grown SiO 2 layer.They were shaped into straight strips with a length, L, of 150 µm and a width, w, of 200 nm.In order to reduce current crowding, the strips were terminated by tapered contacts.The shape of our specimens is shown in Fig. 1(a).Details of the fabrication procedure have been reported elsewhere [25].The specimens were mounted in a closed compartment inside a continuous-flow cryostat.The steady-state CVCs were measured at a set of fixed bath temperatures, T B , between 2.5 − 8 K.
Fig. 1(b) shows typical CVCs for our strips measured in two regimes: sweeping current from zero value upwards (current-bias) and sweeping voltage from a value in the resistive state downwards (voltage-bias).The former regime reduces the impact of bias electronics on the switching current I C while the latter reveals the current plateau at the return current I r caused by the presence of an equilibrium normal domain in the strip.Adopting the heat flow from electrons to the substrate in the form Q = K(T p e − T p B ) (T e is the electron temperature, K is the effective thermal conductance, see the next section) and assigning a constant electron thermal conductivity λ = D c e (T C ) (D is the electron diffusivity, c e is the electron heat capacity, Table I) to the normal and superconducting parts of the strip, we solved numerically the system of steady-state heat balance equations (Eqs.(1)).The obtained T e (x) profiles along the strip (Fig. 1(c)) were computed with the best-fit value p = 3.2 (see below) for a set of lengths of the normal domain with the edges at T e = T C .Knowing the length and, hence the resistance of the domain, we further obtained the voltage along the strip.Corresponding discrete CVC points are shown with symbols in Fig. 1(b).
From experimental CVCs, we extracted I C (T B ) and I r (T B ) dependences, which are plotted in Fig. 2(a) and Fig. 2(b), respectively.For each strip, we computed the theoretical depairing current, I dep , with Eq. (A1) using the material parameters from Table I.We fitted the I dep (T ) dependences to the experimental I C (T B ) data in the vicinity of the superconducting transition (Fig. 2(a)) where both currents are expected to be equal.The value of T C was used as a fitting parameter in order to account for an expected reduction of the transition temperature in nanostructures as compared to non-structured films [26].The best-fit values of T C listed in Table I are indeed slightly smaller than those of non-structured films (8.41 and 9.51 K for films with thicknesses 6 and 9 nm, respectively [13]).At the bath temperature of 3 K, the ratios I C /I dep for our strips are 0.56 (d = 6 nm) and 0.63 (d = 9 nm).These values are comparable to those for other disordered thin films [27].From the combined analytical and numerical solution of the heat balance equations for the electron temperature, we computed I r (T B ) dependences (Eq.( 4)) using the best-fit values of T C and the parameters from Table I.The exponent p and the scaling factor A in Eq. ( 4) were used as the only fitting parameters.The results shown in Fig. 2(b) with solid lines were obtained with p = 3.2 and 2.9 for strip thicknesses 6 and 9 nm, respectively, and with A ≈ 1 for both strips.TABLE I. Parameters of the NbTiN strips studied here.Electron and phonon heat capacities and the electron-phonon energy relaxation times at the transition temperatures (TC ) of the strips were obtained via extrapolation of corresponding values of non-structured films [13] according to ce ∝ T (Drude model), c ph ∝ T 3 (Debye model) and τEP ∝ T −n , respectively.The values of the exponent n were reported in [13].R is the normal-state sheet resistance, D is the electron diffusivity, τesc is the phonon escape time, and τ is the relaxation time of the electron energy (see Appendix B).The steady-state heat flow from a system at a temperature T to a thermal bath with a temperature T B via a thermal link is often described as Q = K(T p −T p B ) where K is an effective thermal conductance [14-24, 28, 29].For a small change in the system temperature ∆T T , the rate of change in the heat flow is dQ/dt ≈ KpT p−1 dT /dt.On the other hand, for a small deviation from the equilibrium, the relaxation is exponential, i.e. dQ/dt = c dT /dt = c ∆T /τ where τ is the relaxation time and c is the heat capacitance of the system.Hence, the effective thermal conductance is K = c(T )/(p T p−1 τ (T )).It may depend on the system temperature by virtue of arbitrary temperature dependences c(T ) and τ (T ).When applying this approach to a thin metal film on a dielec- tric substrate, one has to consider two subsystems, i.e., electrons and phonons in the film.Electrons are heated by the current with the rate I 2 R /(w 2 d) per unit volume.They further transfer the energy to phonons, which in turn release it to the substrate as it is schematically depicted in Fig. 3.There are two limiting cases [24] separated by the value of the ratio between the phononelectron energy relaxation time τ P E = τ EP c ph /c e and the phonon escape time τ esc = 4d/(u s ᾱ), where ᾱ is the angle-averaged transmission of the film/substrate interface given by the acoustic mismatch model [30], u s is the sound velocity, τ EP is the electron-phonon energy relaxation time, and c ph and c e are the phonon and electron heat capacities, respectively.(i) In the limiting case of a thick film, τ esc τ P E , reabsorption by electrons thermalizes nonequilibrium phonons at a temperature T ph which is slightly larger but close to the electron temperature T e .Since in the steady-state conditions the net heat flows from electrons to phonons, Q EP , and from phonons to the substrate, Q ecs , are equal, the rate of heat removal from electrons per unit volume of the film can be represented as either of these two fluxes.For crystalline metallic films, the net heat flux per unit volume from electrons to phonons in the film was described in [23] as e τ EP )(T 5 e − T 5 ph ).The net heat flux from the film to the substrate via 3-d Debye phonons was described in [16] as ).Alternatively, for thick films with T e ≈ T ph , the rate of the energy removal from electrons can be described as a one-stage process (Fig. 3) as Q ≈ c ph (T e )/(p T p−1 e τ esc ) (T p e − T p B ).For 3-d Debye phonons in the film and in the substrate p = 4 [16].This approach has been implemented in Ref. [15].
(ii) In the opposite limiting case of a thin film, τ esc τ P E , nonequilibrium phonons escape to the substrate without being thermalized that leads to overheating of electrons with respect to equilibrium phonons.Microscopic analysis of this essentially two-stage process has shown [24], that even in this case the heat flow from electtrons to the substrate can be phenomenologically described as a one-stage process in the general form Q = K(T p e −T p B ) with K = c e (T e )/(pT p−1 e τ (T e )).Here, τ is the electron energy relaxation time, which appears as the response time in photoresponse measurements, and p = 5+γ where γ varies depending on the phonon dimensionality and the degree of disorder of the strip material.For 3-d Debye phonons and strongly disordered films, γ = 1.We shall note here that the microscopic expression suggested in [24] for the effective thermal conductance K is temperature independent.This fact imposes a constraint on the temperature-dependent quantities entering K: c e (T e ), τ (T e ), and T p−1 e .Although our NbTiN films with τ esc /τ P E ≈ 3.5 fall into the intermediate regime between two limiting cases discussed above, we assume that even in this regime the heat flow from electrons to phonons in the substrate can be described as a one-stage process.Successful description of our experimental results validates the correctness of this assumption.We apply the modified heat balance model to describe self-heating normal domain in a superconducting strip sustained in equilibrium due to the balance between Joule heating via the current I r and the cooling via thermal diffusion and heat flow through the film-substrate interface.The electron temperature distribution T e (x) along the strip is given by the solution of two steady-state heat-balance equations written for unit volumes of the normal (subscript N ) and superconducting (subscript S) parts of the strip The strip itself and its normal part (normal domain) are symmetrically centered at x = 0.The coordinates of the edges of the normal domain ±x N D are defined from the condition T e (x N D ) = T C .Both the electron thermal conductivities λ N and λ S as well as the effective thermal conductances K N and K S may be different in the normal and in the superconducting parts of the strip.Additionally, λ's generally depend on temperature.The boundary conditions (BCs) for Eqs.In order to proceed analytically, we first assume that λ's are temperature independent but different while the effective thermal conductances are the same in both parts of the strip K N = K S = K and equal the value suggested in Ref. [24] for normal films.This simplifying assumption allows for a very accurate semi-analytical description of the experimental data.We discuss the impact of temperature dependent heat conductivity below.With the simplifying assumptions, analytical solution for p = 1 [14] and numerical for p = 4 [15] showed that for x N D , (L/2 − x N D ) L T , where L T is the effective thermal length defining the width of the domain edge, one can introduce two additional approximate BCs: (iv) ∂T e /∂x = 0 at the ends of the strip x = ±L/2; and (v) ∂ 2 T e /∂x 2 = 0 at the strip center x = 0.These additional BC's are easy to recognize in numerically computed and plotted in Fig. 1(c .We also confirmed the validity of these additional BC's for temperature dependent λ's (see the discussion below) and p = 1 by solving numerically Eqs.(1).
For temperature independent but different λ's, the first terms in Eqs. ( 1) reduce to λ N,S ∂ 2 T e /∂x 2 .Defining the Joule temperature as T p J = I 2 r R /(Kw 2 d) and applying BC (v) to the first equation in Eqs. ( 1), one gets where T e0 = T e (x = 0).One can further reduce Eqs.( 1) to two first-order differential equations via the substitution (∂T /∂x) 2 = 2 dT (∂ 2 T /∂x 2 ).Further applying remaining BC's (i) -(iv) to these first order equations and using Eq. ( 2) one can analytically relate T J with T B as follows Solving Eq. ( 3) numerically and substituting K = c e (T e0 )/(pT p−1 e0 τ (T e0 )) in the definition of the Joule temperature one obtains T J (T B ) and the return current as This final expression contains an additional scaling factor A which was used as one of the fitting parameters.We found (see the discussion below) that the maximum envisaged difference in λ's relevant to the N/S interface does not noticeably affect the best fit value of the exponent p.We therefore used λ N = λ S = λ N (T C ) in Eqs.(1) and Eq. ( 3) to compute model I r (T B ) curves (Fig. 2(b)) and the profiles of the electron temperature (Fig. 1(c)).Other parameters were taken from Table I.The τ was computed (Appendix B) as a function of the electronphonon energy relaxation time τ EP , the phonon escape time τ esc , and the ratio between electron and phonon heat capacities c e /c ph , which were taken from Table I.For strips with thicknesses 6 and 9 nm, we obtained the best-fit values p = 3.2 and 2.9 and A = 0.91 and 0.94, respectively.With these values of p and with the Drude temperature dependence c e ∝ T e , the constraint imposed by the temperature-independent effective thermal conductance K [24] implies that τ ∝ T −1 e .Such temperature dependence was indeed found at temperatures around T C for the τ computed in the framework of the two-temperature (2-T) model (Appendix B).It is important to stress here that both the magnitude and the temperature dependence of τ around T C are noticeably affected by the magnitude of c ph .We could satisfy the constraint on K only by using c ph value found in [13], which is noticeably less than the Debye heat capacity expected for a given sound velocity.Our best fit values of the exponent p require γ ≈ −2 in the microscopic theory [24].Recalling that for 3-d Debye phonons γ = 1, we attribute the change in γ and the reduced phonon heat capacity to the size effect imposed by grains on the phonon spectrum in our thin granular NbTiN films.0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 4. Dimensionless electron and phonon temperatures in the center of the normal domain vs. dimensioneless bath temperature computed for the 6 nm-thick NbTiN strip.
It's worth to emphasise here that for computing the return current with Eq. ( 4), the heat capacities c e and c ph should be taken at the electron and phonon temperatures at x = 0, T e0 and T ph0 , and the energy relaxation time τ at T e0 .The phonon heat capacity, c ph , enters Eq. ( 4) through τ , which is a function of τ EP (T e0 ), τ esc , and c e (T e0 )/c ph (T ph0 ).In our computation procedure, T ph0 was found from Eqs. (B1) in the steady-state conditions with q = 2 + n where n was taken from Table I and s = 4 (as for 3-d Debye phonons).Fig. 4 shows how the computed temperatures T e0 and T ph0 vary with the bath temperature.
We discuss now how temperature dependent λ N (T ) and λ S (T ) may affect the best fit value of the exponent p.In the normal state, λ N (T ) = Dc e (T ) ∝ T by virtue of the c e (T ) ∝ T dependence.In the superconducting state, λ S (T ) decreases much quicker.Down to the relative temperature 0.3T C , with a good accuracy, it can be approximated as λ S ≈ λ N 4/3 (T /T C − 0.3) [31].Solving Eqs.(1) numerically for p = 1 with temperature dependent λ's and with λ N = λ S = λ N (T C ), we found that introduction of the temperature dependent conductivities cause a change of the temperature in the center of the normal domain T e0 but do not affect additional BCs (iv) and (v).We therefore used Eq. ( 3) to evaluate the effect of different but temperature independent λ N and λ S for p > 0. Setting the ratio λ S /λ N = 0.5, we found a 15% decrease in the best-fit value of the exponent p.Note that λ N and λ S enter Eq. ( 3) via the boundary condition at x = x N D where they are equal.Physically, the heat flow through this N/S interface can be affected by the temperature distribution in the layer with a thickness of the order of the electron-phonon thermal length l EP = √ Dτ EP which is about an order of magnitude less than the effective thermal length L T .Since the temperature change around the domain edge is T e0 − T B , the temperature difference in the layer with the thickness l EP is ∆T = l EP /L T (T e0 −T B ) ≈ 1 K. Corresponding change in λ S is less than 20% which would cause a correction to the best-fit value of the exponent p remaining beyond our experimental accuracy.We, therefore, neglected the temperature variations in the electron thermal conductivity and used λ S = λ N = λ N (T C ) for the description of our experimental data.
In the interpretation of photoresponse data with the 2-T model in [13], a possible contribution of diffusion cooling was neglected by the authors and τ was associated with the measured response time.In order to check the validity of this approximation, we numerically solved the time-dependent heat balance equation including diffusion for small temperature deviations and computed the response time as a function of the strip length.The dependence shown in Fig. 5(b) (Appendix B) supports the assumption that the diffusion cooling was negligible for strip lengths used in [13].

IV. CONCLUSION
We have analyzed the hysteretic current-voltage characteristics of straight nanostrips fabricated from thin granular NbTiN films.We have shown that the results can be quantitatively explained only with the phonon heat capacity c ph , which is drastically reduced compared to the value expected for 3-d Debye phonons.The same reduced c ph was obtained earlier from photoresponse studies.This shows the compatibility of steady-state and time-resolving experimental approaches for the evaluation of heat capacities and spectra of phonons in nanostructured superconducting thin films.
Furthermore, we have shown that the steady-state experimental approach is self-consistent since it yields the same temperature dependence for the relaxation time of the electron energy as the two-temperature model predicts.
We have also observed that the heat flow from electrons to the substrate ∝ (T p e − T p B ) with the exponent p ≈ 3 which differs from both mediated by the electron-phonon interaction and by escaping of 3-d Debye phonons via film/substrate interface.
This finding, along with the reduced c ph , is attributed to the effect of the mean grain size on the phonon spectrum of thin granular films.Our results provide important insights into thermal transport in thin nanostrips exploited in superconducting devices and thus pave the way for improving their performance.
where P is the power dissipated per unit volume in the electron subsystem.Small periodic variations in the dissipated power ∆P e −jωt (∆P P and ω is the circular frequency) cause small periodic oscillations of T e and T ph , i.e. ∆T e e −jωt and ∆T ph e −jωt .Substituting them in Eqs.(B1) and cancelling steady-state parts, one gets the linearized equations for ∆T e and ∆T ph identical to those suggested in [34].The known solution of linearized equations for ∆T e is given by [34], [12] ∆T e (ω) = Re ∆P c e τ 2 τ 3 τ 1 The characteristic times are τ 1 = (Γ 2 + Γ 3 ) −1 and τ   I.At small temperatures the electron energy relaxation time asymptotically approaches τ EP ∝ T −3.5 while at T T C it saturates at the temperature independent value τ esc .Another important observation is that around T C we find τ EP ∝ T −1 .This temperature dependence meets the constraint imposed by temperature-independent microscopic expression for the effective thermal conductance K.

Impact of diffusion
The 2-T model in the form of Eq. (B1) does not account for the heat removal from electrons via electron diffusion, e.g., to the contacts.In order to check when the diffusion becomes important, we find the relaxation time of the total electron energy in the strip τ * as a function of the strip length by numerically solving the linearized onedimensional time-dependent heat balance equation with diffusion cooling τ ∂T e /∂t = L 2 T ∂ 2 T e /∂x 2 − (T e − T B ), and applying a δ-like uniform heat source and Dirichlet boundary conditions (T e = T B ) at the strip edges.Here L T = √ Dτ is the appropriate (p = 1) thermal length.We further compute the weighted energy 0 δT e (t, x) dx, where δT e (t, x) is the solution of the differential equation, and define τ * from the condition δE(τ * ) = δE(0)/exp (1).As seen in Fig. 5(b), diffusion affects τ * when L ≤ 10 L T .For the 6 nm and 9 nm-thick NbTiN films, L T ≈ 60 nm.The τ was obtained in [13] for 1 µm-long NbTiN bridges, for which L > 10 L T .Therefore, diffusive cooling had no impact on the energy relaxation time and, hence, on the phonon heat capacity evaluated in [13].
FIG. 1.(a) Sketch of the studied specimens (not in scale).(b) Curves (solid and dashed) represent CVCs of the NbTiN strip with d = 6 nm measured at the bath temperature 2.9 K.The legend indicates the sweep directions and the corresponding bias regimes.Symbols represent discrete CVC points numerically computed (Eqs.(1), boundary conditions (i-iii)) with actual strip parameters, λS = λN (TC ) and the best-fit value p = 3.2.The dotted curve is to guide the eyes.(c) Profiles of the electron temperature along the strip numerically computed for the discrete CVC points shown in panel (b); x = 0 corresponds to the center of the strip and of the the normal domain.Curve colors in panel (c) correspond to the symbol colors in panel (b).
(1) are: (i) (∂T e /∂x) N = 0 at the center of the strip x = 0; (ii) the temperature and the heat flux are continuous at the N/S interfaces where additionally the temperature equals the transition temperature, i.e. (T e ) N = (T e ) S = T C and λ N (∂T e /∂x) N = λ S (∂T e /∂x) S at x = ±x N D ; (iii) T e = T B at the ends of the strip x = ±L/2.With these three BCs and p = 1, the system Eqs.(1) has no analytical solution.
) T e (x) profiles for a set of different values of x N D .They were obtained for p = 3.2 under the assumption λ S = λ N = D c e (T C ).In this case L T = D c e (T C )/K T p−1 C ph , and Γ 3 = τ −1 esc .For an electron subsystem obeying only one relaxation time τ , the solution would have the form |∆T e (ω)|= 1/ 1 + (ωτ )2 .In this case, |∆T e (1/τ )|= |∆T e (0)/ √ 2|.We use the very same criterion to define τ for ∆T e (ω) given by Eq. (B2).

FIG. 5 .
FIG. 5. Electron energy relaxation time for 6 nm-thick NbTiN strip in the double-logarithmic scale.(a) τ vs. temperature computed with the uniform 2-T model without diffusion cooling.At low (Te TC ) and high (Te TC ) temperatures, τ asymptotically approaches τEP and τesc, respectively.(b) τ * vs. reduced strip length at Te = TC computed with the heat balance equation including diffusion.For L LT , when diffusion cooling can be neglected, both models give the same relaxation time τ = τ * .