Laser modulation of superconductivity in a cryogenic widefield nitrogen-vacancy microscope

The cryogenic widefield NV microscope used for this study houses the sample in a closed-cycle cryostat with a base temperature down to $4.0$ K (see Methods for details), equipped with a superconducting vector magnet allowing application of a uniform magnetic field up to $1$ T along an arbitrary axis, $\bf{B}_{\text{app}}$. The cryostat contains a cage mounted optics column which includes a high numerical aperture objective lens, and is accessible via a window at the top of the cryostat [Fig. 1(a)]. A $532$ nm laser is used to initialise and readout the NV-spin ensemble, and the red photoluminescence (PL) ($650-800$ nm),Chen et al. (2011) is collected via the same optics column and focused onto an sCMOS camera [Fig. 1(a)]. The NV-diamond used in these experiments is a 50-$\mu$m-thick membrane irradiated to form an NV imaging layer extending $200$ nm below the diamond surface (see Methods).

Four $200$-nm-thick Nb devices were fabricated directly on the diamond surface, each device featuring two square bonding pads ($500$ $\mu$m) connected by a narrow channel ($500$ $\mu$m $\times$ 40 $\mu$m) [Fig. 1(b)] (see Methods). The diamond was mounted to a glass cover slip featuring an omega-shaped microwave resonator for NV spin-state driving, which itself was mounted to a printed circuit board (PCB) for electrical contact to the resonator and Nb devices. This PCB was mounted to a stage equipped with a thermistor and heater for control and measurement of the near-sample temperature. Imaging of the near surface NV-layer occurs through the cover slip and bulk of the diamond, with the devices on the underside of the sample [Fig. 1(c)].

Prior to imaging, the devices were characterised electrically to identify their critical temperature. The resistance across each device ($R$) was measured as a function of temperature as read by the thermistor on the sample holder ($T_{\text{sensor}}$). All devices showed a superconducting transition at a critical temperature $T_{c}\approx 9$ K, where, for example, $R$ decreases from the normal state resistance, $R_{\text{N}}=50.2$ $\Omega$, to the superconducting resistance, $R_{\text{SC}}=34.6$ $\Omega$, which corresponds to the resistance of the non-superconducting leads [Fig. 1(d)]. A typical PL image of the NV-layer beneath the Nb film under continuous wave (CW) illumination shows the Gaussian beam profile, with a $43$ $\mu$m beam waist, giving reasonable illumination across a $100$ $\mu$m field of view [Fig. 1(e)]. At a total laser power of $P_{\text{laser}}=1.0$ mW, we note that imaging the sample results in only a slight change of the temperature as measured by the thermistor $\Delta T_{\text{sensor}}=0.05$ K.

We now move to NV magnetic imaging of a Nb device, focusing initially on vortices in the large Nb pad under no applied current. The sample was cooled to $T_{\text{sensor}}=4.3$ K under a uniform magnetic field perpendicular to the film plane ($z$ axis), $B_{\text{app}}=1.5$ G, with the laser off. The net magnetic field was imaged using CW optically detected magnetic resonance (ODMR), at a range of laser powers. The magnetic field images are presented as the field in the $z$ direction, $B_{z}$, which is calculated from projection of the field along the NV axes (see SI, section II). At $P_{\text{laser}}=0.5$ mW, the $B_{z}$ image shows vortices distributed across the field of view [Fig. 2(a)], with a number density that is consistent with the theoretical value for such films, $n=B_{\text{app}}/\Phi_{0}$, where $\Phi_{0}$ is the magnetic flux quantumStan et al. (2004) (see SI, Fig. S8). At $P_{\text{laser}}=1.0$ mW, vortices disappear from pockets near the centre of the image, while those towards the edge remain fixed [Fig. 2(b)]. At $P_{\text{laser}}=2.0$ and $4.0$ mW, we see a disc centred on the laser spot in which the vortices are removed while the vortices towards the edge remain fixed [Figs. 2(c,d)]. Comparing the $B_{z}$ images with ODMR contrast maps from the same measurements [Fig. 2(e)-(h)], we observe that the regions where the vortices are removed correlate with regions in which the contrast is reduced, indicating a local reduction in the MW field strength (see SI, section II).

These observations are explained by local laser heating of the Nb film raising the temperature of some regions within the field of view above $T_{c}$ [Fig. 2(i)]. As $P_{\text{laser}}$ is increased from $0.5$ mW to $1.0$ mW, pockets of normal state Nb are formed where the local laser intensity is largest (and $T>T_{c}$), thereby removing the vortices and attenuating the MW field (the superconducting state Nb is mostly transparent to the MW field given that the frequencies used are less than the superconducting gap). These regions are highlighted by subtracting the broader Gaussian curve from a PL image [Fig. 2(j)]. Increasing $P_{\text{laser}}$ further, these pockets merge to form a normal state disc centred in the field of view.

Additionally, we observe clustering of vortices around intensity maxima in the laser profile when $P_{\text{laser}}$ is reduced from powers giving large areas of normal state Nb ($P_{\text{laser}}=2.0$ and $4.0$ mW) to a less invasive power ($P_{\text{laser}}=0.5$ mW) [Fig. 2(k,l)]. This is because as the normal region shrinks, the vortices re-nucleate in the superconducting region where they are attracted by the nearest hot spot,Vadimov et al. (2018) i.e. here the normal/superconducting boundary. The vortices therefore cluster around local temperature maxima, where they are pinned once the Nb cools further (see SI, section III, for modelling and further discussion). A faster reduction in $P_{\text{laser}}$ reduces this clustering effect (see SI, Fig. S9). Moreover, uniform vortex configurations are recovered by heating the system globally above $T_{c}$, and cooling in the absence of laser (see SI, Fig. S9). Recently, thermal gradients arising from focused laser beams have achieved patterning at the single vortex level,Veshchunov et al. (2016) while ensembles of vortices can be manipulated by nano-patterned current profiles,Kalcheim et al. (2017) temperature patterning,González et al. (2018) and local magnetic fields.Polshyn et al. (2019)

The images presented demonstrate the invasiveness of widefield NV microscopy in this case. Imaging the Nb film with $P_{\text{laser}}=2.0$ mW gives a normal state region nearly the size of the laser spot, indicating heating in this region upwards of $>5$ K given the base temperature $4.3$ K and the Nb $T_{c}\approx 9$ K. Such a laser power corresponds to a relatively modest peak intensity of $40$ W/cm², which is $4$ orders of magnitude lower than the laser intensity needed to saturate the optical cycling of the NV,Manson et al. (2006) as is often employed in single NV experiments. Surprisingly, the temperature at the sensor remains close to the base temp in these measurements, with $T_{\text{sensor}}<4.4$ K at $P_{\text{laser}}=2$ mW. A linear response of $T_{\text{sensor}}$ to $P_{\text{laser}}$ is observed [Fig. 3(a)], with a slope that depends on the location of the laser beam: the thermistor is heated most efficiently when a greater portion of the laser is incident on the transparent diamond. In any case, $T_{\text{sensor}}$ remains below $4.7$ K for $P_{\text{laser}}$ up to $10$ mW, implying that there exist strong temperature gradients across the sample, which may be overlooked if local indicators are not available.

To further quantify the local heating, we measure the resistance, $R$, as a function of $T_{\text{sensor}}$ (controlled by the heater) and $P_{\text{laser}}$, when the laser spot is focused at various locations around the device. When the laser is focused $200$ $\mu$m from the Nb film, the $R$ versus $T_{\text{sensor}}$ curves shift to lower $T_{\text{sensor}}$ as $P_{\text{laser}}$ increases, indicating a global heating of the Nb device [Fig. 3(b)]. When the laser spot is focused on the Nb bonding pad, the curves shift again, but the shape is severely distorted, as $T_{-}$ moves to lower $T_{\text{sensor}}$ at a faster rate than $T_{+}$ as $P_{\text{laser}}$ increases [Fig. 3(c)] ($T_{+/-}$ defined in Fig. 1(d)). This scenario is exacerbated when the laser is focused on the narrow channel, where the full width of the current path is encompassed by the laser spot [Fig. 3(d)].

The minimum and maximum temperature increase along the current path through the device ($\Delta T_{\text{min}}$ and $\Delta T_{\text{max}}$) can be quantified by the shift in $T_{+}$ and $T_{-}$ at a given $P_{\text{laser}}$ from their values when the laser is off [Figs. 3(e,f)]. In the case where the laser is focused on the narrow channel, where the current must pass through the maximally heated part of the device [Fig. 3(g) lower], we find a maximum temperature increase of $>2$ K/mW, sufficient to completely quench the superconductivity of the illuminated section of the strip at $P_{\text{laser}}=2$ mW. This is consistent with Fig. 2(g) where imaging at $P_{\text{laser}}=2.0$ mW results in a normal state Nb disc centred on the laser spot $\sim 50$ $\mu$m in diameter. The minimum temperature increase of $\sim 0.5$ K/mW suggests that the whole Nb device experiences significant heating even $1$ mm away from the laser spot. This is possibly indicative of a relatively poor thermal conductivity of our implanted diamond substrate, and/or a poor thermal contact with the Nb film. When the laser is focused on the bonding pad, the current can avoid the hottest part of the device directly under the laser spot [Fig. 3(g) upper], and so $\Delta T_{\text{max}}$ is reduced and is closer to the global minimum heating of $\sim 0.5$ K/mW as compared to the previous case [Fig. 3(h)]. The heating of the Nb film is reduced but still measurable with the laser off of the Nb device ($\sim 0.2$ mm away), $\sim 0.1$ K/mW, just a factor $3$ shy of the heating measured at the sensor. Heating of the sample by the microwave field was also assessed and found to give a small global temperature change ($\sim 0.1$ K) across the device at powers relevant to most imaging applications (see SI, section I).

Resistance measurements allow us to characterise local heating due to the laser and infer its impact on the current path, however, the impact of laser heating on the local current distribution can be imaged directly by measuring the Ørsted field using ODMR. Direct reconstruction of current paths in superconductors is of particular interest to superconducting single-photon and single-electron detectors, which rely on quenching superconductivity at the site of detection.Bulaevskii et al. (2012); Adami et al. (2013) Here we apply a biasing field, $\bf{B}\rm_{\text{app}}=(47.5,97.4,19.1)$ G, to resolve the different NV orientations and measure the net vector magnetic field, $\mathbb{B}=(B_{x}$,$B_{y}$,$B_{z})$. This allows us to reconstruct the current density in the Nb device with good accuracy, by inverting the Biot-Savart lawTetienne et al. (2017, 2019) (see SI, section IV). The Ørsted field was measured at three different laser powers, $P_{\text{laser}}=0.25$, $0.5$, and $1.0$ mW, with a constant total current, $I=20$ mA, for all measurements. We show only the field projected in the $z$ direction, $B_{z}$, as this was used to reconstruct the two-dimensional current density, $\bf{J}$ (see SI, section IV).

At $P_{\text{laser}}=0.25$ mW and $0.5$ mW, $R=34.6$ $\Omega$ for the duration of measurement (approximately $1$ hour), indicating fully superconducting current pathways across the device (see Methods). The measured $B_{z}$ [Fig. 4(a,b)] and the associated total current density, $|\bf{J}|$ [Fig. 4(d,e)], show a non-uniform current density distribution, that is more laterally confined at the lower laser power. At $P_{\text{laser}}=1.0$ mW, $R=50.2$ $\Omega$, indicating a fully normal state current pathway, as expected from cascade Joule heating from the comparably large current. Consequently, $B_{z}$ [Fig. 4(c)] and $|\bf{J}|$ [Fig. 4(f)] show a current distribution that is uniform across nearly the full width of the Nb channel. The current density tapers at the edge of the field of view due to a reconstruction error (see SI, Fig. S14, where we show that $|\bf{J}|$ reconstructed from $B_{y}$ does not show this tapering effect). Note that imaging with $P_{\text{laser}}=0.5$ mW but using the heater to raise the temperature above $T_{c}$ ($T_{\rm sensor}=12$ K) gave results identical to the $P_{\text{laser}}=1.0$ mW case, for both the $|\bf{J}|$ map and the ODMR contrast.

The non-uniform current density through the superconducting state Nb channel is a direct consequence of the temperature profile imprinted by the laser [Fig. 4(g)]. As the local temperature increases, the superconducting gap, and hence the critical current density, $J_{c}$, is reduced. The measured current density is therefore larger where the local temperature is lower. Increasing the total laser power reduces $J_{c}$ across the channel, and the current density distribution broadens to maintain the same total current. Line cuts of the current density in the $x$ direction highlight this effect [Fig. 4(h)]. We note that the imaged $|\bf{J}|$ under fixed laser conditions retains a consistent shape across the Nb channel as the field of view is translated, indicating that the non-uniformity observed arises from the excitation laser, rather than local variation in the film (see SI, Fig. S15).

In this work, we demonstrated a cryogenic widefield NV microscope featuring a submicrometer spatial resolution and a field of view of $100~{}\mu$m, and applied it to the imaging of vortices and transport currents in superconducting Nb devices. The demonstrated field of view is five times larger than in the previous demonstration of widefield NV imaging at cryogenic temperatures, which reported a field of view of $20~{}\mu$m.Schlussel et al. (2018) This new capability is ideal for real-space investigations of mesoscopic phenomena in a variety of materials and devices, and also enables imaging of several samples in parallel. For instance, atomically thin samples of van der Waals materials prepared by mechanical exfoliation typically come in the form of multiple micrometre-sized flakes with different properties (thickness, shape), and so widefield imaging of such samples would allow simultaneous studies of many of them, greatly speeding up the characterisation process. This could be particularly useful to investigate, for example, the magnetic properties of ferromagnetic van der Waals materials and heterostructures.

Our work also highlighted a limitation of NV microscopy for low temperature measurements, where the laser illumination required to optically interrogate the NV centres can lead to significant heating of the sample under study. Using the Nb superconducting film ($T_{c}\approx 9$ K) as a local temperature probe, we found that even modest illumination powers (2 mW, corresponding to a peak intensity of $40$ W/cm²) can locally quench the superconductivity of the film, implying that the sample temperature exceeds 9 K even when the temperature measured with a nearby sensor remains below 5 K, close to the base temperature of the cryostat. This work thus demonstrates the need for caution in NV sensing experiments at low temperature, setting a limit to the laser power that can be used for minimally-invasive imaging. In our experiments, a total power of 0.5 mW (peak intensity of $10$ W/cm²) was sufficiently low to keep the Nb devices fully superconducting, allowing an array of frozen superconducting vortices to be imaged. Although the acceptable illumination conditions will depend on the details of the experimental setup and sample under study, it is likely that the most sensitive samples will require strategies to mitigate laser-induced heating. These include the introduction of a thin high-reflectance metallic film between the diamond and sample of interest, the use of better thermal conductors between the sample and the cooling elements of the cryostat, or the implementation of an optimised illumination geometry to reduce the required laser power. These precautions may be necessary even in single NV experiments which can have similar laser power densities at the point of imaging despite using less total power. Such steps will unlock the potential of NV microscopy for a broader range of low temperature condensed matter phenomena.

This work was supported by the Australian Research Council (ARC) through Grants No. DE170100129, CE170100012, LE180100037, FT180100211, DP190101506 and DP200101118. We acknowledge the AFAiiR node of the NCRIS Heavy Ion Capability for access to ion-implantation facilities. D.A.B. and S.E.L. are supported by an Australian Government Research Training Program Scholarship. We acknowledge the contributions of D. Creedon and S. Yianni in developing the fabrication of superconducting Nb films, and A. Martin for useful discussions.

Most NV experiments utilise a microwave resonator or strip-line to drive NV spin-state transitions. The required microwave field power varies greatly depending on the measurement, and the microwave-source-sample configuration. Here, we assess heating of the Nb sample due to the microwave resonator by measuring the device resistance and $T_{\text{sensor}}$ across a range of microwave powers, $P_{\text{MW}}$. The $R$ versus $T_{\text{sensor}}$ curves of the device imaged in the main text at various $P_{\text{MW}}$ show a translation of the critical temperature to lower $T_{\text{sensor}}$ as $P_{\text{MW}}$ is increased [Fig. S4(a)]. The shape of the curve is preserved across the range of $P_{\text{MW}}$, as shown by matching changes in $\Delta T_{\text{min}}$ and $\Delta T_{\text{max}}$ (defined in main text) with $P_{\text{MW}}$ ($\sim 32$ mK/mW) [Fig. S4(b)], indicating negligible temperature gradients across the length of the Nb film. The change in temperature as measured by the thermistor, $\Delta T_{\text{sensor}}$ is comparably small ($\sim 5$ mK/mW), indicating a significant temperature gradient between resonator-sample plane, and the thermistor, which are separated by $>2$ mm. The ODMR based imaging presented in the main text uses $P_{\text{MW}}=2$ mW, and hence, the microwave heating of the Nb film is negligible as compared to the laser heating in the same measurement.

The magnetic images of vortices presented in Fig. 2 of the main text originate from continuous wave (CW) ODMR measurements, at low field, such that the vortex density is compatible with our imaging resolution limit, $\sim 700$ nm. Here, we describe these measurements in detail and the subsequent analysis. We show additional vortex images that demonstrate the dependence of vortex density on the magnetic field, and the dependence of vortex clustering on the cooling protocol used.

CW ODMR images of vortices were produced by PL accumulation over a ($100$ $\mu$m)² field of view [Fig. S5(a)], within the large bonding pad of the Nb film [Fig. S5(b) upper], synchronised with the continuous laser illumination and microwave (MW) field at a given MW frequency [Fig. S5(b) lower]. ODMR spectra in this region, acquired under a field $B_{0}=1.5$ G in the $z$-direction, exhibit two broad resonances corresponding to the two electron spin transitions of the NV centres, $f_{+}$ and $f_{-}$ [Fig. S5(c)]. In this low field regime the exact line shape of these two resonances is non-trivial due to the ensemble nature of the measurement (ensemble of NVs with different symmetry axes and local charge/spin/strain environments) and due to fluctuations of each NV’s environment (charge and spin). In particular, the narrow feature separating the two resonances has been explained by Mittiga et al.,Mittiga et al. (2018) as arising from random distributions of both local electric field (caused by charge fluctuations in nearby defects) and local magnetic field (spin fluctuations). While in principle it is possible to use the model from Ref. Mittiga et al. (2018) to fit our data and extract the mean value of the magnetic field, we found the convergence of the fit to be very sensitive to the initial guess, and so were not able to obtain satisfying fits across the whole set of data ($>10,000$ spectra per image, each having slightly different microwave broadening conditions, different local environments etc.). Instead, we fitted our data with a simple sum of two Lorentzian functions with independent frequencies, amplitudes and widths (solid lines in Fig. S5(c)). While this fit function does not capture the narrow central feature and as a result may lead to systematic errors in the estimation of the magnetic field, we found this to be the most robust and sensitive method and so is well suited to the imaging of small magnetic features such as vortices.

The fit parameters for each resonance, obtained by fitting spectra across an entire image, show long-range variations that are somewhat correlated to the laser intensity distribution [Fig. S5(d) - (i)]. These variations are thus attributed to fitting errors due to the simplified model as explained above. For instance, the fits to the two spectra in Fig. S5(c) indicates frequencies ($f_{-}=2867.8$ MHz, $f_{+}=2886.4$ MHz) for Area 1, against ($f_{-}=2866.8$ MHz, $f_{+}=2888.0$ MHz) for Area 2, but the difference can be explained by a change in the width of the resonances, caused mainly by a difference in laser intensity between the two areas, rather than by an actual difference in magnetic field. The resulting magnetic field is calculated from the fit frequencies, by the approximate formula $B_{\text{NV}}=(f_{+}-f_{-})/(2\gamma_{e})$,Rondin et al. (2014) [Fig. S5(j)]. The non-physical long-range variations that carry over from the frequency maps have been removed by applying a smoothing filter to the image, which fits features in the raw $B_{\text{NV}}$ map that vary over more than $20$ pixels [Fig. S5(k)] and subtracts them, leaving only the short range features (i.e. the vortices). All the data in Fig. 2 of the main text have been obtained by applying the fitting and post-processing method outline above. Note that the $B_{\text{NV}}$ thus extracted corresponds to the average projection of the magnetic field along the four different NV axes, since the different NV families are not resolved in the ODMR spectra. For a field pointing in the $z$-direction (i.e. normal to the diamond surface, which is the case at the centre of the vortices), the projection is identical for all four NV families in our $\langle 100\rangle$-oriented diamond, given by $B_{\text{NV}}=B_{z}/\sqrt{3}$. The $B_{z}$ images in Fig. 2 of the main text have been scaled accordingly.

The magnetic images produced by the method outlined above, allow us to determine the size and magnitude of the magnetic field features associated with the vortices. Examining our $B_{\text{NV}}$ images over smaller areas better highlights the features of individual vortices [Fig. S6(a)]. A line profile through one such vortex shows the peak measured field projected along the NV axis is approximately $\sim 0.3$ G, with a full width at half maximum (FWHM) of $1.4$ $\mu$m [Fig. S6(b)], comfortably above our optical resolution limit, $\sim 700$ nm. To compare this result with theory, we calculated the expected vortex size following Carneiro et al. Carneiro and Brandt (2000). We consider a straight vortex in a $200$ nm thick film with an isotropic London response, with London penetration depth, $\lambda=400$ nm, and superconducting coherence length, $\xi=8$ nm, at temperature $T/T_{c}=0.75$, as inferred from Fig. 3 in the main text. The NV response is accounted for by considering this field at a distance $130$ nm below the film, which corresponds to the mean NV depths, as shown in Fig. S1. The field as seen by the NV is taken as the average projection of this field along each NV axis, which we can compare directly to the measured distribution [Fig. S6(b)]. The simulated vortex field has a FWHM of $600$ nm, less than half the width measured, and a peak field strength of $4$ G, an order of magnitude above our measured value. This discrepancy in size is likely due to disorder within the film, which will alter the straightness of the vortex within the film, and hence its size, and explain the non-uniform vortex sizes across our field of view. The peak field discrepancy is partially explained by the observed broadening, but is likely dominated by our simplified fitting model, which neglects residual splitting of the NV spin-resonance lines, and hence reduces $B_{\text{NV}}$ Broadway et al. (2018, 2019).

Local reductions in ODMR contrast, as seen in Fig. 2 of the main text, are attributed to a reduction in the local microwave field strength due to the normal state Nb or superconducting Nb with a superconducting gap close to the microwave frequency. However, PL spin contrast can be reduced by a number of effects including the intensity of the excitation laser, large magnetic field gradients,Tetienne et al. (2018b) and interactions of NVs with nearby defectsBluvstein et al. (2019) and materials.Lillie et al. (2018) To validate our interpretation, we compare ODMR contrast and Rabi curve measurements of the same region, imaged at base temperature with $P_{\text{laser}}=1.0$ mW. An ODMR contrast map of an area within the Nb bonding pad taken under these conditions shows localised pockets of reduced PL contrast towards the centre of the image [Fig. S7(a)], similar to those seen in Fig. 2(f). Rabi curves from these pockets show slower Rabi oscillations as compared to neighbouring regions with greater ODMR contrast [Fig. S7(b)], i.e. for a given interaction time with the microwave field, the spin-state dependent PL evolves less rapidly in these regions, giving up to $0.2\%$ more PL for an interaction time of $0.15$ ms. Quantitative comparisons with ODMR contrast measurements are non-trivial given the pulsed nature of the Rabi measurement versus the CW ODMR measurement, however, this result demonstrates that locally reduced microwave field power at least partially explains the reductions in ODMR contrast. We note that the laser heating in each measurement should be comparable given that a $1$ s initialising laser pulse was used prior to camera exposure ($0.5$ s), which has a minimum laser duty cycle $0.2$ for the microwave pulse durations explored.

The number of vortices within a given area is such that the sum of their fluxes is equal to the total flux passing through the same area in the normal state. Each vortex is threaded by a magnetic flux quantum, $\Phi_{0}$, and hence $N\,\Phi_{0}=\bf{B}_{0}\cdot\bf{A}$, where $N$ is the number of vortices within an area, $\bf{A}$, at a magnetic field, $\bf{B}_{0}$. The vortex images presented in Fig. 2 in the main text were taken at a magnetic field $B_{0}$ $\approx 1.5$ G in the $z$-direction. This field is actually due to some background field and residual magnetisation within the unshielded cryostat, as there was no field applied by the superconducting coils for these measurements, i.e. $B_{\text{app}}=0$ G. Here we demonstrate the effect of varying our applied magnetic field strength, $B_{\text{app}}$, which differs slightly to the field seen at the sample $B_{0}$, given the un-shielded cryostat and remanent magnetisation.

Varying $B_{\text{app}}$ from $2.0$ to $-4.0$ G in the $z$-direction, the number of vortices within the field of view changes considerably [Fig. S8(a)-(e)]. Counting the number of vortices within a fixed area, ($33.81$ $\mu$m)², and fitting the number as a function of the applied field, shows a strong linear relation with an effective zero field at $B_{\text{app}}=-1.5(1)$ G. Accounting for this offset in field strength, $B_{0}=B_{\text{app}}+1.5$ G, the number of vortices compares well with theory [Fig. S8(f)], with a slight discrepancy likely due an uncalibrated magnification from a non-ideal optics setup. We note that the $B_{0}$ values quoted in the main text for current imaging do not include this correction to the field at the sample, as it is small compared to the applied field.

Vortex images shown in the main text demonstrate clustering of vortices around local hot-spots when the laser power is reduced from powers which quench superconductivity, $P_{\text{laser}}\geq 1.0$ mW, to a less invasive imaging power, $P_{\text{laser}}=0.5$ mW. The final arrangement of vortices after cooling is sensitive to both the initial laser power, as shown in the main text in Figs. 2(k) and (l), but also the rate at which the laser power is reduced. Reducing the laser power from $P_{\text{laser}}=4.0$ mW [Fig. S9(a)] to $P_{\text{laser}}=0.5$ mW by turning down the power over a few seconds and then imaging at $P_{\text{laser}}=0.5$ mW [Fig. S9(b)] results in greater clustering than if the beam is cut immediately (in $\approx 10$ ns using the acousto-optic modulator) and then imaged at $P_{\text{laser}}=0.5$ mW once $T_{\text{sensor}}$ has stabilised [Fig. S9(c)]. Similar clustering is observed when the cryostat is cooled to base temperature from $T_{\text{sensor}}>T_{c}$ with the laser on, at $P_{\text{laser}}=0.5$ mW [Fig. S9(d)]. A uniform vortex configuration is recovered when the cryostat is heated above $T_{c}$ and cooled to base temperature in the absence of laser [Fig. S9(e)], indicating that the vortex clustering is rewritable.

In the following, we describe how to model superconducting vortices in a temperature gradient, induced for instance by heating a sample locally with a laser.

We consider a rectangular geometry with periodic boundary conditions, where we define a periodic temperature gradient (shown in Fig. S10(a,b)) given by

	$\displaystyle\tau(x,y)$	$\displaystyle=a_{0}+\bigg{[}\cos\Big(x\frac{2\pi}{L_{x}}\Big{missing})+\cos\Big(y\frac{2\pi}{L_{y}}\Big{missing})-2\bigg{]}\frac{a_{0}-a_{1}}{4},$		(1)
	$\displaystyle\tau(x,y)$	$\displaystyle=1-\frac{T(x,y)}{T_{c}},$		(2)

where $L_{x}$ and $L_{y}$ denote the length of the system in $x$- and $y$-direction, respectively. $a_{0}$ ($a_{1}$) denotes the highest (lowest) appearing value of $\tau$. Eq. (1) is constructed such that $\tau=a_{0}$ in the corners of the rectangle and $\tau=a_{1}$ in its middle, i.e.

	$\displaystyle\tau(0,0)=\tau(0,L_{y})=\tau(L_{x},0)=\tau(L_{x},L_{y})$	$\displaystyle=a_{0},$
	$\displaystyle\tau(L_{x}/2,L_{y}/2)$	$\displaystyle=a_{1}.$

Using the London equationLondon et al. (1935),

$$\vec{h}(x,y)+\lambda^{2}\nabla\times\nabla\times\vec{h}(x,y)=0,$$

(3)

we can calculate the magnetic field, $\vec{h}$, and ultimately the energy of a given configuration using the London energy functional which is given by

$$E=\int\big{[}\vec{h}^{2}+(\lambda\nabla\times\vec{h})^{2}\big{]}{\rm d}\vec{A}.$$

(4)

With the Maxwell equation $\nabla\cdot\vec{h}=0$, the temperature dependance of the penetration length, $\lambda=\lambda(0)/\sqrt{\tau}$, and $\vec{h}=h\hat{e}_{z}$, Eq. (3) becomes

$$\tau h-\lambda(0)^{2}\nabla^{2}h=0.$$

(5)

Implementing an arbitrary magnetic flux configuration $\Phi(x,y)$ into the system is simply done by adding it as a inhomogeneity to the right hand side of Eq. (5), which leads to

$$\tau h-\lambda(0)^{2}\nabla^{2}h=\Phi(x,y).$$

(6)

The integral of $\Phi$ over the whole space yields the total flux introduced. For a single vortex at position $\vec{r}_{0}$ we add the flux in form of a delta distribution, $\Phi(\vec{r}_{0})=\Phi_{0}\delta(\vec{r}_{0}-\vec{r})$, with the flux quantum $\Phi_{0}$. Multiple vortices at positions $\vec{r}_{i}$ are introduced by simply summing up their respective inhomogeneities, i.e.

$$\Phi(x,y)=\Phi_{0}\sum_{i}\delta(\vec{r}_{i}-\vec{r}).$$

(7)

Since we calculate the energy of a given vortex configuration from the total magnetic field, using Eq. (4), vortex-vortex interactions are automatically accounted for. This energy calculation can then be used for a Monte Carlo simulation to find the ground state for an arbitrary temperature gradient and for any number of vortices.

In order to image transport currents, as shown in Fig. 4 of the main text, we need to accurately measure the magnetic field projection (including the sign) of at least one NV family out of the four available. This requires to apply a bias magnetic field of at least $40$ G typically (for the diamond sample used here) to separate one NV family from the other three in the ODMR spectrum. Moreover, to improve the accuracy of the reconstruction it is preferable to measure all four NV families simultaneously rather than just one.Tetienne et al. (2017, 2019) This can be achieved by aligning the bias magnetic field such that all four NV families can be resolved in the ODMR spectrum. Here we apply a magnetic field $\bf{B}_{\text{app}}=\rm{(47.5,97.4,19.1)}$ G expressed in lab frame coordinates (see definition of the $xyz$ axes in Fig. 1(e) of the main text). In this frame the NV axes have unit vectors $\bf{u}_{\text{NV}}=\rm{(\pm 1,\pm 1,1)/\sqrt{3}}$. Example ODMR spectra recorded under this field are shown in Fig. S12(a).

To analyse the ODMR data, we first fit the spectrum at each pixel with a sum of eight Lorentzian functions with free frequencies, amplitudes and widths (solid lines in Fig. S12(a)). The eight resulting frequencies $\{f_{i}\}_{i=1\dots 8}$ are then used to infer the total magnetic field, $\bf{B}_{\rm tot}$, by minimising the root-mean-square error function

$$\varepsilon(D,{\bf{B}}_{\text{tot}})=\sqrt{\frac{1}{8}\sum_{i=1}^{8}\left[f_{i}-f_{i}^{\rm calc}(D,\bf{B}\rm_{tot})\right]^{2}}$$

(8)

where $\{f_{i}^{\rm calc}(D,\bf{B}\rm_{tot})\}_{i=1\dots 8}$ are the calculated frequencies obtained by numerically computing the eigenvalues of the spin Hamiltonian for each NV orientation,

$\displaystyle{\cal H}$

$\displaystyle=DS_{Z}^{2}+\gamma_{\rm NV}\bf{S}\cdot\bf{B}~{},$

(9)

and deducing the electron spin transition frequencies. Here $\mathbb{S}=(S_{X},S_{Y},S_{Z})$ are the spin-1 operators, $D$ is the temperature-dependent zero-field splitting, $\gamma_{\text{NV}}=28.035(3)$ GHz/T is the isotropic gyromagnetic ratio, and $XYZ$ is the reference frame specific to each NV orientation, $Z$ being the symmetry axis of the defect Doherty et al. (2012, 2013) with unit vector $\bf{u}_{\rm NV}$ defined previously. The resulting total magnetic field $\bf{B}_{\rm tot}$ contains both the applied bias field and the field from the Nb film due to the magnetic response as well as the transport current (Ørsted field). Since the applied field is uniform over the field of view, we simply subtract a constant offset to $\bf{B}_{\rm tot}$, taken to be the field measured far from the Nb wire, to obtain the Nb-induced field. Illustrative results of this analysis are shown in Fig. S12, corresponding to the case studied in Fig. 4(b) and (e) of the main text ($I=20$ mA, $P_{\rm laser}=0.5$ mW): the PL image is shown in (b), the background-subtracted $B_{z}$ map in (c), the contrast and width of the lowest ODMR resonance in (d) and (e), the $D$ parameter in (f) and the residue $\varepsilon$ in (g). The fit residue seems to correlate with the contrast: $\varepsilon<1$ MHz beside the Nb wire but it reaches $2$ MHz under the Nb where the ODMR contrast is lower as can be seen in the ODMR spectra of Fig. S12(a). This suggests that $\varepsilon$ is dominated by the noise in the data and that the model Eq. (9) is adequate. In comparison, the shifts induced by the Oersted field are much larger, up to $10$ MHz. The zero-field splitting parameter $D\approx 2877$ MHz is roughly as expected at this temperature,Doherty et al. (2014) with a step of $\sim 0.5$ MHz under the Nb wire attributed to strain.Broadway et al. (2019)

The three components $(B_{x},B_{y},B_{z})$ of the magnetic field are shown in Fig. S13 for the three cases studied in Fig. 4 of the main text. To separate the Ørsted field from the magnetic response of the Nb film, we performed the same measurements prior to applying the transport current (just after cooling the sample below $T_{c}$) and after turning off the current following the measurement with current presented before. The $B_{z}$ maps for those three cases are shown in Fig. S14, along with corresponding maps of the ODMR contrast. Upon cooling in the applied field of $\bf{B}\rm_{app}=(47.5,97.4,19.1)$ G, the field exhibits a small expulsion as apparent from a reduction of $B_{z}$ by $\approx 0.4$ G under the Nb wire (Fig. S14(a)). The expected density of vortices due the applied $B_{z}$ is about $1$ vortex$/\mu$m²,Stan et al. (2004) hence the weak expulsion is in broad agreement with the London penetration depth estimated from the device resistance ($\lambda\approx 400$ nm). Upon turning the current on, the field becomes dominated by the Ørsted contribution, reaching $B_{z}\pm 5$ G (Fig. S14(b)). Turning off the current, the small field expulsion is recovered but some additional features have appeared, indicating a rearrangement of the vortices during the application of the current. Because this response field may be different with the current on, and because in any case it remains much smaller than the Ørsted field, we will simply ignore it when reconstructing the current density. That is, the reconstructed current density will also include a small contribution from the supercurrent associated with the response to the applied field. This supercurrent is localised near the edges of the wire and has opposite signs for each edge, so that the net current integrated across the width of the wire will be only due to the transport current.

Comparing the ODMR contrast maps in Fig. S14, a drop in contrast can be seen at the centre of the wire with the current on, whereas the contrast is relatively uniform in the absence of current. From the observations made in Fig. 2 of the main text, we tentatively ascribe this drop in contrast to a portion of the Nb wire turning normal, but not extending through the entire thickness of the film since the resistance remains $R=R_{\text{SC}}$. In the vortex imaging presented in Fig. 2 of the main text, there was no significant change in local contrast at $P_{\text{laser}}=0.5$ mW, however here the presence of a relatively large current density may lower $T_{c}$ sufficiently to allow the local heating due to the laser to locally quench the superconductivity.

We now describe the methods for reconstructing the current density in the Nb film. In principle, the three magnetic field components $(B_{x},B_{y},B_{z})$ are related to each other via Ampère’s law, $\nabla\crossproduct{\bf B}=0$,Lima and Weiss (2009) and there are several ways (nominally equivalent) to use them to reconstruct the current density. For instance, the projected current density $\tilde{\bf J}=\int\bf J\,\rm dz$ can be obtained from the in-plane components $B_{x}$ and $B_{y}$ using Tetienne et al. (2019)

	$\displaystyle\tilde{J}_{x}$	$\displaystyle=$	$\displaystyle-\frac{2}{\mu_{0}}B_{y}$		(10)
	$\displaystyle\tilde{J}_{y}$	$\displaystyle=$	$\displaystyle\frac{2}{\mu_{0}}B_{x},$		(11)

where we utilised the fact that the Nb-NV distance (between $0$ - $400$ nm) is small compared to the lateral spatial resolution of the measurements ($\approx 1\,\mu$m given by the optical resolution).Tetienne et al. (2019) Below, this method will be referred to as the first method. Alternatively, we can use the $B_{z}$ component together with the continuity of current ($\nabla\cdot{\bf J}=0$) to obtain relations in the Fourier plane,

	$\displaystyle\tilde{\mathcal{J}}_{x}$	$\displaystyle=$	$\displaystyle-\frac{2}{\mu_{0}}\frac{ik_{y}}{k}\mathcal{B}_{z}$		(12)
	$\displaystyle\tilde{\mathcal{J}}_{y}$	$\displaystyle=$	$\displaystyle\frac{2}{\mu_{0}}\frac{ik_{x}}{k}\mathcal{B}_{z}.$		(13)

Here $\mathcal{F}(k_{x},k_{y})$ denotes the two-dimensional Fourier transform of $F(x,y)$, where ${\bf k}=(k_{x},k_{y})$ is the spatial frequency vector and $k=\sqrt{k_{x}^{2}+k_{y}^{2}}$. This is the second method.

The projected current density reconstructed from these two methods are shown in Fig. S15 for the three cases presented in Fig. 4 of the main text. We note that the current densities shown in the main text are in fact projected current densities, labeled as $\bf{J}$ (without the tilde) for simplicity. There are several differences between the two methods. Considering the case where the Nb wire is in the normal state (Fig. S15(right column)), we see that the first method gives a $|\tilde{\bf{J}}|$ map that looks sharper and more uniform along the Nb wire than the second method, for which the current appears tapered near the left and right boundaries of the image. This is because the latter is prone to truncation artefacts when calculating the Fourier transform of $B_{z}$, and those are more pronounced near the boundaries where the $B_{z}$ data has more noise. In contrast, there is no Fourier transform involved in the first method and so no truncation artefact; the noise in $|\tilde{\bf{J}}|$ directly mirrors the noise in the $B_{x}$, $B_{y}$ data. While this observation would normally motivate the use of the first method, examination of the maps in Fig. S15(middle column) reveals another effect. With the first method, there appears to be a “missing” current near the centre of the image. Indeed, integrating the $x$-component of the current density across the width of the wire, i.e. $I_{\rm int.}=\int\tilde{J}_{x}\,dy$, we approximately recover the electrically measured current of $I=20$ mA near the boundaries of the image but the value of $I_{\rm int.}$ drops to $14$ mA (a $30$% drop) near the centre. This corresponds to a drop in $B_{y}$ seen in Fig. S13(b) and also correlates with a drop in ODMR contrast observed in this region (third column in Fig. S15) indicating a local quenching of superconductivity. At lower laser power (Fig. S15(left column)), the contrast is mostly uniform again along the Nb wire, indicating that it is fully superconducting with no normal regions, and the missing current is recovered. In previous work, a similar effect of missing current was observed in normal metals at room temperature,Tetienne et al. (2019) where it was interpreted as an apparent delocalisation of the current towards the diamond substrate, which decreases $|\tilde{\bf{J}}|$ in the first method but not in the second method. A possible explanation for this effect may involve a magnetic response in the diamond itself which can be significant for nitrogen-doped diamond at low temperatures,Barzola-Quiquia et al. (2019) however this requires further investigation (see Ref.Tetienne et al. (2019) for a more complete discussion). For this reason, we chose to use the second method in the main text, assuming the resulting $|\tilde{\bf{J}}|$ is a good indicator for the transport in the Nb film. With the second method, the integrated current $I_{\rm int.}$ remains within $10$% of the electrically measured value ($20$ mA) in all cases at any point along the wire.

So far we have presented results from a single location along the Nb wire. In Fig. S16, we show the current density $|\tilde{\bf{J}}|$ (obtained with the second method) for two other locations [Fig. S16(a,d) and (c,f)] and compare with the original location [Fig. S16(b,e)]. The current pattern is broadly similar in all cases, supporting the interpretation that it is mostly dictated by the laser profile which imprints a temperature profile. Moving the sample relative to the laser does not change the temperature profile in the field of view and therefore, the current density pattern is essentially unchanged. There are small differences still, and those somewhat correlate with features seen in the ODMR contrast [Fig. S16(d-f)]. For instance, at the left most location [Fig. S16(a,d)] a hot spot of current can be seen which corresponds to an island of normal ODMR contrast (i.e. a fully superconducting region) in the middle of a patch of reduced ODMR contrast (i.e. a region of reduced superconductivity).