Heat propagation models for superconducting nanobridges at millikelvin temperatures

Nanoscale superconducting quantum interference devices (nanoSQUIDs) most commonly use Dayem bridges as Josephson elements to reduce the loop size and achieve high spin sensitivity. Except at temperatures close to the critical temperature Tc, the electrical characteristics of these bridges exhibit undesirable thermal hysteresis which complicates device operation. This makes proper thermal analysis an essential design consideration for optimising nanoSQUID performance at ultralow temperatures. However the existing theoretical models for this hysteresis were developed for micron-scale devices operating close to liquid helium temperatures, and are not fully applicable to a new generation of much smaller devices operating at significantly lower temperatures. We have therefore developed a new analytic heat model which enables a more accurate prediction of the thermal behaviour in such circumstances. We demonstrate that this model is in good agreement with experimental results measured down to 100 mK and discuss its validity for different nanoSQUID geometries.

(Some figures may appear in colour only in the online journal) Nanoscale superconducting quantum interference devices (nanoSQUIDs) are probably the most promising sensors for performing highly sensitive magnetic characterisation of solid state systems with very small spin populations, e.g. magnetic nanoparticles [1], magnetic molecules [2] and related entities such as flux qubits [3]. The spin sensitivity of a nanoSQUID can be improved in several ways. Firstly, operating at temperatures well below 4.2 K reduces the thermal noise floor, and allows the investigation of magnetic systems that show quantum effects at these temperatures, for example some oxide interfaces [4] or topological insulators [5]. Secondly, decreasing the loop size reduces the geometric inductance of the SQUID and improves the flux-to-voltage transfer function, ( ) = F F V V d d max , and the flux sensitivity. Finally, decreasing the track width allows for the flux coupling between the sample and the SQUID to be improved if the sample can be attached to the structure [6]. Various techniques have been used to fabricate the weak links required in nanoSQUIDs, but the most commonly used technique is to form nanobridge constrictions (Dayem bridges) in superconducting thin films patterned using either electron beam lithography [7] or focussed ion beam lithography [8,9].
However at very low temperatures the current−voltage characteristics of most nanobridges are affected by thermal hysteresis which makes reading out the SQUID challenging [10]. This thermal hysteresis has a different origin from the electrical hysteresis which can arise for SQUIDs based on tunnel junctions due to their much larger capacitance.
One common approach to help reduce thermal hysteresis is to add a normal metal capping layer on top of the superconductor to help remove heat. This can be effective for Nb devices close to T c , but at lower temperatures the critical current will be higher, and hysteresis generally reappears. Another approach is to use a lower T c superconductor such as titanium to reduce the critical current, but even then the SQUID geometry still has to be optimised to overcome the excess heat extraction problem [11]. Thus it is desirable to use a thermal model when designing any type of nanoSQUID for sub-kelvin operation. The first widely accepted model for this hysteresis was suggested by Skocpol, Beasley and Tinkham (SBT) [12] who analysed tin microbridges operating close to 4.2 K. In their model, once the bias current is increased above the critical current I c , a normal 'hotspot' region is created in the bridge and adjoining banks, and is self-sustained by the energy dissipated by the Joule effect. They showed that in the hysteretic regime, the bias current has to be reduced well below I c , to a level called the retrapping current I r , for the device to return to its superconducting state. The non-hysteretic regime occurs when > I I r c . SBT considered that heat propagates one-dimensionally along the bridge then radially (in-plane) in the banks. By solving the heat equation for this system, they obtained an analytical expression for the radius r 0 of the normal region in the banks which is proportional to a thermal length scale they called the thermal healing length, h k a = d , where d is the thickness of the film, α is the heat transfer coefficient per unit area to the substrate, and κ is the thermal conductivity of the film, which they assumed to be a constant and similar for the superconducting and normal states close to T c .
Recently Hazra, Pascal, Courtois and Gupta (HPCG) [13] carried the SBT analysis further by including the decrease of the thermal conductivity of the superconductor k s below T c due to the formation of Cooper pairs which cannot carry heat. In the HPCG model, this effect is approximated as a decrease of the form k k = T T , s n c where k n in the thermal conductivity in the normal state just above T c estimated using the Wiedemann−Franz law. The typical temperature variation close to a nanobridge in the regime considered by HPCG is shown in figure 1(a).
These models allowed the authors to successfully analyse their devices operating close to 4.2 K, but in general they are not strictly applicable to smaller devices, or the case that we are especially interested in, i.e. devices operating at much lower temperatures. This is because they assume that the banks have sufficient surface area to give a physically possible solution for r 0 at equilibrium. This is reasonable for say micron-sized Nb devices at 4. give h m 5 m for a 100 nmthick film. However it clearly becomes an issue in modelling nanoscale devices as the dimensions are shrunk below η, or for modelling devices with even reasonably large banks at lower temperatures since the heat transfer to the substrate [14] will become poorer, and thus η will become larger depending on the materials used. For instance, at 100mK the Nb device mentioned above would have k~--0.2 W m K 1 1 and a~-100 W m K 2 , which would give h m 14 m. This means the existing models tend to overestimate the heat extraction in the region close to the nanobridge. This in turn means a larger current is needed to sustain the hotspot in equilibrium, leading to an overestimated I r well below T c . For that reason we have developed a simple new analytical model for our devices operating in that regime.

A simple thermal model for millikelvin temperatures
To develop a model for our devices in the millikelvin regime, we make the simplifiying assumption that the thermal physics is dominated by one-dimensional heat flow along the extended length of the banks rather than radial heat flow from the ends of the nanobridge into the banks as assumed by SBT and HPCG. The thermal length scale will be such that we consider the banks to include both the arms of the nanoSQUID loop and a stretch of the electrical connecting track to the nanoS-QUID. For the moment we will consider the bank as an idealised long thin structure extending away from the nanobridge and ignore the likely curved geometry of the real SQUID loop. We will re-consider the effect of a curved geometry in section 4.
As shown in figure 1(b), we consider a nanobridge of b , of normal state resistivity r b connected to infinitely long banks of width W, thickness d and The temperature of the bridge including a small semicircular region at each end is assumed to be uniform. In the semi-infinite banks the temperature decreases radially to reach T c at a critical radius = r r 0 . The bridge and the region < r r 0 are in the normal state and dissipating heat. The region > r r 0 is in the superconducting state. The substrate is assumed to be at the bath temperature ¥ T . (b) Schematic of the heat distribution in the regime considered in the present paper. The increased thermal length scale is such that the heat flow is assumed to be one dimensional in the banks reaching T c at a critical distance = x x 0 . Unlike the HPCG model, the banks are considered to have a finite width W.
normal state resistivity ρ. We assume that the bridge is small enough to have a uniform temperature T 1 . The system is completely symmetrical with respect to the centre of the nanobridge and the problem is equivalent to half a bridge transmitting heat to a single bank. We assume perfect materials with uniform thickness and only consider the system in the steady-state since our current−voltage measurements are typically quasi-dc. We also assume the power generated by the Joule effect in the bridge is uniformly spread over the whole width of the bank. This approximation is accurate if  h W 1, which is likely to be the case at millikelvin temperatures considering the typical nanoscale dimensions and the large η due to poor heat coupling to the substrate. In the opposite limit (  h W 1), corresponding to larger devices at liquid helium temperatures, the SBT and HPCG models are more physically accurate.
At any position x along the axis parallel to the banks, in the steady state we have where the terms on the left hand side represent the heat flow per second along the film, the heat transferred from the banks to the substrate, and the heat transferred from the bridge to the substrate respectively. The terms on the right hand side represent the power generated by Joule effect in the bank As the heat transfer to the substrate is greatly reduced at low temperatures and the nanobridge to substrate contact area is very small, we can neglect the heat flow from the bridge to the substrate to simplify the equations without loss of generality. We assume ( ) k T has the same dependency on temperature as used by HPCG, and that the resistivity is a constant over the temperatures of interest.
By considering the derivative of equation (1) with respect to x with these assumptions, we obtain the system of equations: The differential equation for the normal region ( < x x 0 ) has a general solution of the form ( ) x T sinh p , where A and B are constants determined from boundary conditions, and T p is a particular solution which is a temperature defined by: To solve the differential equation for the superconducting region, we apply the change of variable = y T T , it is possible to offer an analytical solution. This approximation is fairly accurate over a wide range of temperatures and becomes more precise as y approaches ¥ t 2 . Equation (3) can then be rewritten which has a general solution corresponding to a temperature variation in the superconducting region, where C and D are constants.
To determine the constants in the above two solutions, we use the boundary conditions. Considering the heat flow at the end of the nanobridge, equation (1) at x=0 gives Considering the temperature as we approach the normal/superconducting boundary from the normal side, Then by considering that ( ) Finally, using the continuity of T x d d at the normal/superconducting boundary we have Combining equations (7)-(9) and writing we can factor out the current I as: The retrapping current can be obtained by finding the minimum of this function either graphically as shown in figure 2, or numerically. In practice, this is a weak minimum so that the value tends to be relatively close to ( ) = I x 0 0 in most cases, and then Figure 2 shows the position of x 0 along the x-axis for a typical nanobridge made of titanium (nominal = T 400 mK c ). As expected, a higher bath temperature means a smaller I r and a hotspot extending further along the film. Unlike the HPCG model, the width of the bank is taken into account and strongly affects the value of I r . We can finally plot the retrapping current as a function of temperature for a typical bridge as shown in figure 3. As expected, I r is well overestimated by the previous models in this regime.

Comparison with experiment
In order to test the predictions of the new model in the millikelvin regime we made a series of Ti and Ti/Au nanobridge SQUIDs by e-beam lithography and liftoff in the geometry shown in figure 4. The Ti devices were made from a 90 nm thick Ti thin film deposited by e-beam evaporation. The Ti/Au nanobridges were made from a Ti(90 nm)/Au(15 nm) bilayer. Our Ti films have a higher T c than the nominal value for the pure material. In the Ti/Au devices the Au layer reduces the T c due to the proximity effect, and is consistent with there being a good interface quality. The devices were measured at different temperatures down to 100 mK using an Adiabatic Demagnetization Refrigerator. The SQUID loops were made fairly large with a track width half that of the connecting tracks to make it easy to model the combined heat flow generated in the two weak-links. From separate measurements of the SQUID modulation, we infer that the two weak links have critical currents matched to better than 10%.
For most SQUID applications it is desirable for devices to be non-hysteretic, but in order to be able to measure I r it is necessary for the devices to be hysteretic, i.e. < I I r c . This is almost always the case for Ti devices due to the poor heat extraction below T c . For the Ti/Au devices we deliberately kept the layer thicknesses thin enough to give hysteretic behaviour below some temperature in our measurement range. The measured and modelled values of I c and I r are shown in figures 5 and 6 for Ti and Ti/Au devices respectively. To model the temperature dependence of ( ) I T c , we fit the Kulik −Omelyanchuk (KO-1) theory [15] for the I R c n  product of short metallic weak links in the dirty limit to the experimental data. Since the Au layer will contribute a normal shunt resistance in the Ti/Au devices, we just treat R n as a fitting parameter (which we term R n,0 ) together with T c . We then fit our model for I r to the measured data and extract estimates of ρ and α. When applying our model to the Ti/Au bilayer device we assume that the superconducting/normal layers are thin compared to the relevant coherence length x x Ti Au and we can approximate the bilayer as a single material of averaged characteristics, so that ( )( ) r r r r Ti Au Au Ti Ti Au Au Ti where r d , Au Au and r d ,

Ti
Ti are the resistivity and thickness of respectively Ti and Au. The thermal conductivity is then determined just above T c using the Wiedemann−Franz law. For the data in figure 5 the fitted value for ρ is consistent with the estimated resistivity of our Ti film at these temperatures. As can be seen in figure 3 the HPCG model predicts a reasonably similar temperature dependence to our model, however, when we carried out (not shown) a fit to the measured data, it yields a  similar value of α but a resistivity of W 1711 n . m, which is an order of magnitude higher than the known resistivity of the film in this range of temperature, and not physically plausible. As explained above, this is a consequence of the HPCG model overestimating the heat extraction in our measurement regime.
The success of fitting the model to the Ti/Au devices shows it could be further applied to determine the optimal thickness of normal metal to deposit on a superconductor to get non-hysteretic I −V characteristics ( < I I c r ) at a given temperature. It has to be sufficiently thick to carry enough heat away, but not excessively thick to suppress the superconductivity more than necessary, which would reduce I R c n and as a consequence reduce the flux sensitivity of the device.

Heat saturation current
As we discussed previously, as we bias a SQUID above I c the heat generated diffuses into the arms of the loop. As the bias current increases the hotspot can rapidly reach such an extent that not only the bridges turn normal but the superconductivity of the entire loop is lost. We call this current the heat saturation current I sat which is an important consideration for the bias operation range of the device. Having sufficient margin between I sat and I c is desirable for avoiding unwanted thermal runaway. Our model enables the direct determination of I sat by using equation (10) which links the extent of the hotspot to the current that would allow thermal equilibrium. We find I sat by considering the limit  ¥ x 0 in equation (10). I sat is the largest current the system can sustain while still being able to return to the superconducting state at some point. The HPCG model cannot be used to predict this behaviour since it does not take into account the finite dimensions of the system: in HPCG the ( ) I r r relationship diverges as  ¥ r , which means the hotspot can spread indefinitely.
Our model shows that the equation for this asymptotic behavior is given by: This provides a useful estimate of the bias current range available to operate the device, but also enables the extraction of thermal parameters even from devices that are non-hysteretic down to the base temperature of the system. Figure 7 shows I c and I sat for a non-hysteretic Ti/Au device. Here the nanobridge dimensions were identical to device in figure 6   but the Ti/Au bilayer was inverted to allow more Ti to be selectively added on the banks using a shadow mask technique [11]. This increases the value of the thickness d in equation (12), increasing the usable bias range of the device and also reducing the kinetic inductance of the SQUID loop. The latter improves the flux sensitivity without introducing the adverse effect of a larger I c  that would arise should we increase the bridge thickness as well.

Comparison with finite element simulations
To justify the simplifying assumptions we made in our analytical model, we also used finite element methods (using COMSOL Multiphysics) to model the heat flow produced by Joule heating when a given bias current passes through the nanobridge and banks. In the finite element model we were able to include heat flow into the subtrate under the nanobridge bridge and the lateral temperature variation across the facing edges of the two banks close to the nanobridge, both of which are ignored in our analytic model. For typical conditions encountered at millikelvin temperatures, we can observe that the radius of curvature of the isotherms is larger than the width of the track within 300 nm of the nanobridge, i.e. within h 0.23 , where η is the thermal healing length in these conditions. Compared to our model, the difference between the predicted position of the superconductor-normal interface, x 0 , is less than 100 nm (or h 0.07 ), which indicates that our model is consistent with reality.

Conclusion
In this paper we have described a new simple thermal model more adapted to describe the behaviour of nanobridges at millikelvin temperatures. By taking into account the geometry of the tracks, it can be successfully applied to study even non-hysteretic nanobridges, therefore making possible the extraction of thermal parameters from functional devices such as nanoSQUIDs.