Manipulating multi-vortex states in superconducting structures

TL;DR

This paper presents a technique for controlling and mapping vortex configurations in superconducting structures using a magnetic tip, enabling operations like vortex braiding relevant for topological quantum computing.

Contribution

It introduces a novel method to manipulate and visualize multi-vortex states in superconductors with high spatial resolution, facilitating quantum information processing.

Findings

High-resolution vortex configuration maps obtained

Localized vortex stabilization achieved via magnetic tip

Demonstrated potential for vortex braiding operations

Abstract

We demonstrate a method for manipulating small ensembles of vortices in multiply-connected superconducting structures. A micron-size magnetic particle attached to the tip of a silicon cantilever is used to locally apply magnetic flux through the superconducting structure. By scanning the tip over the surface of the device, and by utilizing the dynamical coupling between the vortices and the cantilever, a high-resolution spatial map of the different vortex configurations is obtained. Moving the tip to a particular location in the map stabilizes a distinct multi-vortex configuration. Thus, the scanning of the tip over a particular trajectory in space permits non-trivial operations to be performed, such as braiding of individual vortices within a larger vortex ensemble -- a key capability required by many proposals for topological quantum computing.

Tables

Table 1. Table 1: Parameters of structures used in experiments.

Ring	R($\mu$m)	w(nm)	$T_c$(K)
1	0.95	133	$\sim 0.95$
2	1.99	236	1.199
3	0.94	130	$\sim 0.98$
4	0.94	125	$\sim 0.98$

Equations

$$F\propto \frac{1}{1+\beta }\left({\widetilde{{\varphi }_1}}^2+2\beta \widetilde{{\varphi }_1}\widetilde{{\varphi }_2}+{\widetilde{{\varphi }_2}}^2\right)$$

$$\frac{{\nu }_r}{{\omega }_0}=\frac{\Delta \omega }{\Delta \gamma }$$

$$B_c(T)=B_c(0)\sqrt{1-T/T_c},$$

$$B_c(0)=\frac{\sqrt{3}}{\pi }\frac{{\Phi }_0}{w\xi (0)},$$

$$F=F(\hat{\varphi },\hat{n}).$$

$${\gamma }_k\equiv \Delta {\phi }_k-(2\pi /{\Phi }_0){\int }_{\text{k-th wire}}\mathbf{A}\cdot d\mathbf{s}.$$

$$\sum _{j:i\text{-th loop}}{\gamma }_j=2\pi (n_i-{\varphi }_i).$$

$$\sum _{j:p\text{-th node}}{I}_j=0,$$

$$I_j=\frac{1}{{\mu }_0{\lambda }^2}\left(\frac{{\Phi }_0}{2\pi }\right)\frac{{\gamma }_k}{l_j}{S}_j,$$

$$L_j^K={\mu }_0{\lambda }^2\frac{l_j}{S_j}.$$

$$I_j=\frac{{\Phi }_0}{2\pi }\frac{1}{L_j^K}{\gamma }_j.$$

$$\sum _{j:p\text{-th node}}\frac{{\gamma }_j}{L_j^K}=0.$$

$$E_j\approx E_j^K=\frac{L_j^K{I}_j^2}{2}={\left(\frac{{\Phi }_0}{2\pi }\right)}^2\frac{{\gamma }_j^2}{2L_j^K}.$$

$$E=\sum _j^{N_{\mathrm{links}}}{E}_j({\gamma }_j)={\left(\frac{{\Phi }_0}{2\pi }\right)}^2\sum _j^{N_{\mathrm{links}}}\frac{{\gamma }_j^2}{2L^K}.$$

$$\{\begin{array}{rl}& {\gamma }_1+{\gamma }_3=-2\pi \widetilde{{\varphi }_1}\\ & {\gamma }_2-{\gamma }_3=-2\pi \widetilde{{\varphi }_2}\\ & \frac{{\gamma }_1}{L_1^K}=\frac{{\gamma }_3}{L_3^K}+\frac{{\gamma }_2}{L_2^K}\end{array}$$

$$E=\frac{1}{2}{\Phi }_0^2{L}_{\Vert }^K\left[\frac{{\widetilde{{\varphi }_1}}^2}{L_1^K{L}_3^K}+\frac{{\widetilde{{\varphi }_2}}^2}{L_2^K{L}_3^K}+\frac{(\widetilde{{\varphi }_1}+\widetilde{{\varphi }_2}{)}^2}{L_1^K{L}_2^K}\right],$$

$$E=\frac{1}{2}\frac{{\Phi }_0^2}{L^K(\beta +1)}\left[{\widetilde{{\varphi }_1}}^2+2\beta \widetilde{{\varphi }_1}\widetilde{{\varphi }_2}+{\widetilde{{\varphi }_2}}^2\right],$$

$$\{\begin{array}{rl}& \left|\widetilde{{\varphi }_1}+\beta \widetilde{{\varphi }_2}\right|<1/2\\ & \left|\beta \widetilde{{\varphi }_1}+\widetilde{{\varphi }_2}\right|<1/2\\ & |\widetilde{{\varphi }_1}-\widetilde{{\varphi }_2}|<1\end{array}$$

$$\beta =\frac{l_2}{l_1}$$

$$P_a^{eq}(\hat{\varphi })={\Gamma }_b/({\Gamma }_a+{\Gamma }_b).$$

$$d{P}_a/dt$$

$$\ddot{x}+2{\gamma }_0\dot{x}+{\omega }_0^{2}x=\frac{{\omega }_0^2}{k}[f(t)+\zeta (t,x)],$$

$$P_a(t)\simeq P_a^{eq}(\widehat{{\varphi }_0})+\delta P{e}^{i(\omega t-\theta )},$$

$$\delta P=\frac{\partial P_a^{eq}}{\partial \hat{\varphi }}\delta \hat{\varphi }\cos \theta ,$$

$$\theta =\arctan \left(\frac{\omega }{{\nu }_r}\right).$$

$$⟨\hat{I}(t)⟩=\widehat{I_a}(\hat{\varphi })P_a(t)+\widehat{I_b}(\hat{\varphi })(1-P_a(t)).$$

$$⟨\hat{I}(t)⟩=[\widehat{I_a}(\widehat{{\varphi }_0})-\widehat{I_b}(\widehat{{\varphi }_0})]\delta P{e}^{i(\omega t-\theta )}+\frac{1}{2}\frac{\partial (\widehat{I_a}+\widehat{I_b})}{\partial \hat{\varphi }}\cdot \delta \hat{\varphi }.$$

$$\Delta \omega \simeq \frac{{\omega }_0}{2}\frac{{\nu }_r^2}{{\nu }_r^2+{\omega }_0^2}K,$$

$$\Delta \gamma =-\frac{{\omega }_0}{2}\frac{{\omega }_0{\nu }_r}{{\nu }_r^2+{\omega }_0^2}K,$$

$$K=\left(\frac{\hat{\kappa }({\mathbf{r}}_{tip})}{k}\frac{d\hat{\varphi }}{dx}\right)\left(\Delta \hat{I}(\hat{\varphi })\frac{\partial P_a^{eq}}{\partial \hat{\varphi }}\right).$$

$$\frac{{\nu }_r}{{\omega }_0}=\frac{\Delta \omega }{\Delta \gamma }.$$

Full text

The dynamics of superconducting vortices plays an important role in many phenomena. In bulk superconductors the motion of macroscopic number of vortices is responsible for flux flow, pinning and creep Tinkham (1996). In mesoscopic superconducting structures, the energetics and dynamics of small numbers of vortices gives rise to a plethora of unexpected vortex states and mesoscopic effects Bezryadin et al. (1996); Geim et al. (1997, 1998); Morelle et al. (2004); Grigorieva et al. (2006); Timmermans et al. (2016); Embon et al. (2017), the understanding of which often requires numerical simulations Peeters et al. (2000); Baelus et al. (2000, 2001); Baelus and Peeters (2002). Being topological defects of superconducting order parameter, vortices play a key role in such profound effects as Aharonov-Casher effect Aharonov and Casher (1984) and Berezinskii-Kosterlitz-Thouless transition in two-dimensional superconductors Kosterlitz and Thouless (1973). In recent years, superconducting vortices have also been proposed as resource for the implementation of topological quantum computation. It has been suggested that vortices in superconductor-topological-insulator-superconductor junctions may host Majorana bound states (MBS) Fu and Kane (2008); Beenakker (2013); Ren et al. (2019) permitting topological quantum computation. More recently, Abrikosov vortices in iron-based superconductors have emerged as a promising platform for manipulating MBSs. The normal state of these materials possesses the necessary topological properties, hence the cores of Abrikosov vortices could support MBSs; recent experiments provide evidence to support the existence of MBSs in iron-based superconducting vortices Zhang et al. (2019); Wang et al. (2018); Liu et al. (2018); Zhang et al. (2018).

Considerable efforts have been devoted to developing new ways to manipulate individual vortices. The control of vortices has been demonstrated by means of electrical currents Kalisky et al. (2009); Embon et al. (2015); Kalcheim et al. (2017); Ji et al. (2016); Mills et al. (2015); de Souza Silva et al. (2006), focused laser beams Veshchunov et al. (2016), local mechanical stress Kremen et al. (2016), local magnetic fields Gardner et al. (2002); Straver et al. (2008); Auslaender et al. (2009), local heating Ge et al. (2017, 2016), and also was investigated numerically Ma et al. (2018); Milošević and Peeters (2010). In most of these studies, vortices were controlled one at a time, which is an excellent approach for investigating the physics of vortex pinning at the nanoscale. On the other hand, for applications such as vortex braiding, it is essential to be able to control individual vortex states, and to be able to evolve these states along desired trajectories within the space of possible vortex configurations.

In this letter we use a recently developed variant of magnetic force microscopy (MFM), which we refer to as $\Phi_{0}$-MFM Polshyn et al. (2018); Polshyn (2017), to probe and manipulate individual multi-vortex states in multiply-connected mesoscopic superconducting structures. The devices, shown schematically in Fig. 1b, are patterned from a thin aluminum film into rings containing one, two, three and four equal-area sectors. Below the superconducting temperature of aluminum ($T_{c}\approx$1.2 K) and in the presence of an applied magnetic flux, these structures can host discrete multi-vortex states, described by the number of vortices in each sector. We use the spatially inhomogeneous magnetic field produced by a micron-sized magnetic particle attached to a cantilever to access complex multi-vortex states, many of which simply can not be stabilized by a homogeneous external magnetic field. We show that as the particle is scanned over a particular structure, the transitions between different vortex sates are marked by a shift in the resonance frequency and dissipation of the cantilever, arising from the cantilever-driven transitions between states with different vortex configurations. We show that the spatial patterns of the frequency shift vs. position that emerge can be mapped directly to transitions between different multi-vortex states. Furthermore, we demonstrate that this mapping enables deterministic control of the multi-vortex states of the superconducting device.

The schematic diagram of $\Phi_{0}$-MFM microscopy setup is shown in Fig. 1a. The $\Phi_{0}$-MFM images were obtained using an ultra-soft silicon cantilever with a micron-size $\mathrm{SmCo_{5}}$ particle attached to the tip. Both the cantilever and the magnetic moment of the $\mathrm{SmCo_{5}}$ particle are perpendicular to the plane of the carrier chip, as shown in Fig. 1b. The magnetic particle generates a static magnetic field $B_{z}(x,y)$ and modulation $\delta B_{z}(x,y)\sin(\omega_{0}t)$ (caused by cantilever oscillations at its resonant frequency $\omega_{0}$) in the pane of the sample. Two cantilevers used for this work have resonance frequencies $\omega_{0}/2\pi\simeq$ 4146 Hz and 7675 Hz and spring constants $k=$0.11 mN/m and 0.18 mN/m, respectively (for more information and fabrication details see Supporting Information SI ). The motion of the cantilever is detected using a fiber-based laser interferometer operating at wavelength of 1510 nm. MFM measurements are performed in a frequency detection mode Albrecht et al. (1991), in which the cantilever is self-oscillated at its resonant frequency, and the frequency is measured using phase-locked-loop (PLL). A PID feedback loop is used to maintain a constant oscillation amplitude between 2.5-7.5 nm and to measure changes in dissipation. A silicon chip carrying the superconducting structures is mounted below the cantilever on a stack of nano-positioners, that permits approach, coarse positioning and scanning of the superconducting device with respect to the cantilever. All scanning is done at a fixed height $h$ between 0.2-2 $\mu$m above the surface of the sample. All measurements were conducted in vacuum in a continuous flow $\mathrm{{}^{3}He}$ refrigerator operating at a base temperature of 350 mK.

In this work we study narrow-wall aluminum rings with radii between 0.5-2 $\mu$m and wall width $w$ between 100-200 nm. The rings are divided into two, three or four sectors by radial crossbars. The devices were patterned by e-beam lithography of PMMA resist spun on a 500-nm thick $\mathrm{SiO_{2}}$ layer grown on a silicon substrate. After exposure and development, the PMMA was metalized by e-beam evaporation with a 5-nm thick Ti adhesion layer and a 45-nm thick Al device layer. The ring structures were obtained by lift-off of the PMMA. We measured more than fifteen structures; here we report the representative results for four of them.

We start by considering vortex states supported by a structure with a ring geometry divided in two halves by a crossbar (Fig. 2b). If the wall width is sufficiently small compared to the superconducting coherence length, which applies to all the structures reported in this work, then the vortex sates are described by pairs of winding numbers $(n_{1},n_{2})$ of the superconducting phase around the bottom and top halves of the ring. For narrow and thin-wall structures, the free energy of the vortex state depends only on $(n_{1},n_{2})$ and fluxes $(\phi_{1}$, $\phi_{2})$ threading the two halves of the structure, where $\phi$ is the magnetic flux in units of the flux quantum $\Phi_{0}=h/2e$. For a symmetric ring with a crossbar (see Supplementary InformationSI ), the energy can be expressed as

[TABLE]

where $\tilde{\phi_{i}}=\phi_{i}-n_{i}$. The parameter $\beta$ controls the coupling between the two halves of the ring. The limiting cases of $\beta=0$ and $\beta=1$ correspond to either two isolated sectors or a vanishing crossbar, respectively. While in general both the magnetic and kinetic inductances contribute to the coupling parameter $\beta$, structures for which $wt<\lambda^{2}$, where $t$ is the device thickness and $\lambda$ is the superconducting penetration depth, the kinetic contribution dominates. For such devices, which includes the one shown in Fig. 2b, a value of $\beta\approx 0.39$ is expectedSI . The region of stability of each vortex state $(n_{1},n_{2})$ in $(\phi_{1},\phi_{2})$ coordinates has a shape of a hexagon centered on $(n_{1},n_{2})$ point (Fig. 2c). The stability diagram of the vortex states is similar to the honeycomb stability diagram of charging states in a double quantum dot van der Wiel et al. (2002). This similarity arises because both systems share the same form of the effective energy. The coupling parameter $\beta$ controls the shape of the hexagons, and can be determined experimentally from the modulation depth of hexagons $\beta=l_{2}/l_{1}$, as shown in Fig. 2c.

The image of the cantilever frequency vs. tip position allows us to directly map the cantilever position to a particular vortex state. Figure 2a shows an $\Phi_{0}$-MFM scan of Ring 1, which has a radius of $R=0.95$ $\mu$m and wall width $w=133$ nm, divided by a crossbar of the same width into two halves. The scan is taken with tip-surface separation $h=0.8$ $\mu$m, at 0.915 K, close to the superconducting transition temperature $T_{c}\approx$0.95 K. In this regime near $T_{c}$, the transitions between vortex states become reversible. The dark colored contours, shown in Fig. 2a, indicate the shift of the resonance frequency of the cantilever , corresponding to transitions between two tip-induced vortex states – the number of vortices changes by one in at least one of the sectors of the structure. At these points, the energies of two or more vortex states are degenerate and the field modulation $\delta B_{z}$ caused by the tip oscillations drives the vortex transitions. The transitions are accompanied by a switching of the supercurrents in the structure between two distinct configurations, which causes a strong back-action on the cantilever and gives rise to a shift of the resonance frequency Polshyn et al. (2018); SI .

Every position of the tip $(x_{tip},y_{tip},h)$ relative to the superconducting structure corresponds to certain values of $(\phi_{1},\phi_{2})$ induced in the structure. Thus, each $\mathrm{\Phi_{0}}$-MFM image is equivalent to a nonlinear mapping from $(x_{tip},y_{tip},h)$ to ($\phi_{1}$, $\phi_{2}$) coordinates. The mapping is determined by the magnetic flux induced in each sector of the device for a given tip position. We use $\mathrm{\Phi_{0}}$-MFM image shown in Fig. 2a and known dimensions of Ring 1 to reconstruct the mapping $(x_{tip},y_{tip},h)\to(\phi_{1},\phi_{2})$. We start with a model of the magnetic tip and calculate $(\phi_{1},\phi_{2})$ for every point of the scan $(x_{tip},y_{tip},h)$. Then, using equation (1), we find the positions of the vortex transitions for a given tip model and parameter $\beta$. Further, the model of the tip is tuned to obtain a good match between the simulated and observed vortex transitions, as shown in Fig. 2d (see Supporting Information for details). Best fit is obtained with $\beta=0.40\pm 0.02$, which is in close agreement with the expected value of 0.39.

In Figs. 2d-e, we show the correspondence between selected tip positions in $\Phi_{0}$-MFM image of Ring 1 and the points $(\phi_{1},\phi_{2})$ on the vortex configuration diagram. Establishing this correspondence identifies all vortex states, some of which are labeled in Fig 2a. The shaded region in Figs. 2e highlights the vortex states that were accessed during the scan. It is worth noting, that in a homogeneous magnetic field, only points for which $\phi_{1}=\phi_{2}$ can be accessed.

The spatial distribution of the static magnetic field $B_{z}(x,y)$ and modulation $\delta B_{z}(x,y)$, obtained from matching vortex transitions in Fig. 2a, are shown in Figs. 2f and 2g respectively. The field modulation is highly spatially inhomogeneous and has two peaks with the opposite sign – one in front and one behind the cantilever about $\sim 1\,\mathrm{\mu m}$ apart, in the oscillation direction of the cantilever. This spatially separated character of the modulation in $\Phi_{0}$-MFM enables efficiently driving transitions between different multi-vortex states.

Most of the static flux generated by the magnetic tip is confined to an area of size $\sim h^{2}$. Thus, it is possible to vary the total number of vortices in the structure by simply changing the tip-surface separation. Figures 3a-d show a series of $\Phi_{0}$-MFM images of Ring 2 ($R=1.99$ $\mathrm{\mu m}$, $w=230\,\text{nm}$) taken at several tip-surface separations $15$ mK below the superconducting transition $T_{c}=1.199$ K of the structure. As $h$ decreases from 1.5 $\mathrm{\mu}$m to 0.8 $\mathrm{\mu}$m, the flux induced by the tip through the structure grows, and the maximum total number of vortices imaged in a scan increases from 7 to 13. Even at the smallest tip-surface separation, the transitions between vortex states are extremely sharp ($\sim 15$ nm in real space) and manifest clear honeycomb pattern as shown in the fine scan of the selected area of Fig. 3d, shown in Fig. 3e.

Controlling the dynamics of vortex transitions is crucial for realizing metastable vortex states. $\Phi_{0}$-MFM provides a unique way of visualizing the vortex energy barrier between different parts of the device. We measure the vortex dynamics in our superconducting devices using a previously described technique, which relies on the dynamical interaction between thermally activated phase slips (TAPS) Little (1967); Langer and Ambegaokar (1967) in the superconducting device and the cantilever Polshyn et al. (2018). Near $T_{c}$, the field modulation caused by the motion of the tip in the presence of TAPS can give rise to a dynamical force that significantly shifts the frequency and dissipation of the cantilever. In the regime where the vortex transition rate $\nu_{r}$ becomes comparable to the cantilever frequency $\omega_{0}$, a simple relationship exists between $\nu_{r}$ and the variations in the cantilever frequency $\Delta\omega$ and dissipation $\Delta\gamma$  Polshyn et al. (2018); SI .

[TABLE]

Figure 3f shows images of $\nu_{r}$ calculated using (2), for the area marked in Fig. 3e. This region contains three individual vortex transitions: (A)$(5,5)\leftrightarrow(5,4)$, (B)$(5,5)\leftrightarrow(6,4)$ and (C)$(5,4)\leftrightarrow(6,4)$, measured at several temperatures near $T_{c}$. The color range in Fig. 3f represents the measured relaxation rate $\nu_{r}/\omega_{0}=\Delta\omega/\Delta\gamma$, while the brightness represents the magnitude of the signal $\sqrt{\Delta\omega^{2}+\Delta\gamma^{2}}$. As is evident from Fig. 3f, the three transition rates $\nu_{r}$ vary differently with respect to temperature, indicating slight differences in the energy barrier responsible for vortex entry corresponding to the different transitions (for more details see SI ). Transitions A and C correspond to vortices entering/leaving the bottom and top halves of the ring respectively, while transition B corresponds to a vortex moving from one half of the ring to the other (see inset in Fig. 3e). For each type of transition, the corresponding phase slip occurs in a different part of the device. The variations in the transition rates could be result of small variations in the crossectional area of the device or the presence of defects.

We can generate more complex vortex states by fabricating superconducting structures with larger number of partitions. For the purposes of demonstration, we study vortex states hosted by circular structures partitioned into three and four equal area sectors, as shown in Fig. 4a and 4e. $\mathrm{\Phi_{0}}$-MFM images of the vortex states in these structures reveal complicated patters of transitions that evolve with tip-sample separation (Fig. 4b-d). Despite the complexity of the transition patterns however, the majority of vortex states can be easily identified from symmetry considerations, as shown in Fig. 4b and Fig. 4f. In general, the full set of transitions observed in a given image can be calculated using a model for the field distribution from the tip and the analytic form of the free energy of the vortex states. Images Fig. 4b-d,f demonstrate several unique capabilities of $\mathrm{\Phi_{0}}$-MFM approach:

1. The majority of the observed vortex states are stable only in the presence of the inhomogeneous magnetic field from the tip. Therefore, the ability to scan the magnetic particle provides a robust and highly reliable means of accessing a large range of vortex configurations. (2) The $\Phi_{0}$-MFM images provide a deterministic mapping between the position of the cantilever and a particular vortex state. This mapping enables a simple and intuitive means of evolving the system through a desired trajectory of vortex configurations by simply tracing the tip through the corresponding trajectory in space.

As a starting point, we demonstrate a method for braiding two vortices in the three-sector device shown in Fig. 4e. We choose $h$=1.6 $\mathrm{\mu}$m and a uniform external field of $B$=0.9 mT, such that for tip positions close to the center of the ring, only two of the sectors are populated with vortices (see Fig. 4g,h). One complete circular motion of the tip around the center drives the structure through a sequence of vortex states, as shown in Fig. 4i, and exchanges the positions of the vortices. Two full circles of the magnetic tip will accomplish a winding of one vortex around another. While, in our system we control fluxons in a superconducting wire network, the same idea could be applied to Josephson or Abrikosov vortices in appropriate superconducting structures.

In summary, we have shown that $\Phi_{0}$-MFM can be used to probe complicated multi-vortex states in narrow-wall multiply connected superconducting structures. The spatially inhomogeneous magnetic field generated by MFM tip enables stabilization of complex vortex states that are not accessible in a uniform field. Strong interaction of the magnetic tip and supercurrents in superconducting structure at points where tip-induced transitions occur enables the mapping of vortex transitions and the characterization of their dynamic properties. This approach provides a versatile way to induce, identify, and manipulate complex multi-vortex states in superconducting structures. Finally, as a potentially useful application of the technique, we describe how $\Phi_{0}$-MFM may be used to braid a pair of vortices around each other, which might be useful for experiments with Majorana bound states. Unique capability to work with multi-vortex or complex vortex states opens new possibilities in the exploration of vortex interactions and designing new vortex-based superconducting devices.

I acknowledgement

We are grateful to Nadya Mason, Taylor Hughes and Alexey Bezryadin for useful discussions. This work was supported by the DOE Basic Energy Sciences under DE-SC0012649, the Department of Physics and the Frederick Seitz Materials Research Laboratory Central Facilities at the University of Illinois.

II Supporting information

II.1 Experimental procedures

In our measurements we use custom made ultra-soft silicon cantilevers with $\mathrm{SmCo_{5}}$ particles mounted at the tip, shown in Fig. S1. To fabricate them, we first attach a $\mathrm{SmCo_{5}}$ magnetic particle of appropriate size to the tip of the cantilever using a micromanipulator. G1 epoxy (Gatan Inc.) is used to glue the particle. The magnetic moment of the particle is aligned with the axis of the cantilever by applying a magnetic field while epoxy is being cured. Next, the particle is trimmed to a desired regular shape using a focused ion beam (FIB). To avoid the ion damage of the magnetic particle, we use low ion currents: 40-20 pA for the rough cuts and 1 pA for the finishing cuts. The magnetic moment of the particle is characterized using cantilever torque magnetometry.

For measurements reported here we use two MFM tips (Fig. S1). Tip A is used to measure Ring 1,3 and 4, while Tip B is used to measure Ring 2. Tip A has cantilever length $L=110~{}\mathrm{\mu m}$, resonant frequency $f_{0}\simeq 4146$ Hz, spring constant $k=1.1\times 10^{-4}$ N/m and quality factor $Q\simeq 10^{3}$. The $\mathrm{SmCo_{5}}$ particle of tip A has magnetic moment $m_{\mathrm{tip},\parallel}=2.2\cdot 10^{-13}$ J/T. Tip B is mounted on a cantilever with length $L=80~{}\mathrm{\mu m}$, resonant frequency $f_{0}\simeq 7675$ Hz, spring constant $k=1.8\times 10^{-4}$ N/m and quality factor $Q\simeq 31,800$ at 4 K. The $\mathrm{SmCo_{5}}$ particle has magnetic moment with components $m_{\mathrm{tip},\parallel}=7.2\cdot 10^{-13}$ J/T and $m_{\mathrm{tip},\perp}=3\cdot 10^{-14}$ J/T.

The superconducting parameters of the measured structures are measured using the same MFM setup. The superconducting transition temperature $T_{c}$ is measured by parking a MFM tip above the structure and monitoring the resonance frequency shift as a function of temperature (Fig. S2 a). The superconducting coherence length $\xi(0)$ is determined from the suppression of the superconducting transition by external magnetic field (Fig. S2)b. The critical field of a thin wall ring perpendicular to the plane of the ring is given by Polshyn et al. (2018)

[TABLE]

with

[TABLE]

where $w$ is the width of the wall of the ring. We obtain $\xi(0)$ from the known value of $w$ and the value of $B_{c}(0)$ determined from the fit (see Fig. S2b). For Ring 2 with wall width $w$=236 nm we find $\xi(0)=106$ nm. From the temperature dependence $\xi(T)=\xi(0)/\sqrt{1-T/T_{c}}$, we estimate that at temperatures used in this work ($T/T_{c}>0.97$) the superconducting coherence length $\xi>630$ nm, which is much larger than the width and thickness of the wall of the studied structures. The same is also true for superconducting penetration depth, which we estimate to be $\lambda>1$ $\mathrm{\mu m}$ using $\lambda(T)=(1-T/T_{c})^{-0.5}\times 166$ nm obtained for similar aluminium structures Polshyn et al. (2018). Thus, the structures are in the limit of 1D superconducting networks with negligible magnetic screening. The parameters of the superconducting structures used in experiments are shown in Tab. 1.

II.2 Energetics of vortex states in multi-loop superconducting structures

Here we derive the energy of multi-loop superconducting structures in which the vortices can sit only inside the loops. A vortex state of such structure, comprised of $N$ elementary loops, can be described by $N$ winding numbers $\hat{n}=(n_{1},\dots,n_{N})$ – the numbers of vortices hosted by each loop. We assume, that the walls of the wires, comprising the structures, are narrow and thin ($w,t\ll\lambda,\xi$) and the direct effect of the magnetic field of the MFM tip on the wires can be neglected. In this case, the free energy of the vortex states, hosted by the structure, depends only on the fluxes threading the elementary loops $\hat{\phi}=(\phi_{1},\dots,\phi_{N})$, rather than on the full spatial distribution of the magnetic field:

[TABLE]

Let us consider a superconducting wire network with $N_{\mathrm{loops}}$ elementary loops, $N_{\mathrm{links}}$ wires and $N_{\mathrm{nodes}}$ nodes. It can be shown that $N_{\mathrm{links}}=N_{\mathrm{loops}}+N_{\mathrm{nodes}}-1$. The gauge-invariant phase difference for each wire is given by:

[TABLE]

Two types of the equations, similar to Kirchoff’s rules, can be written for the network. The first type follows from the fluxoid quantization conditionTinkham (1996) applied to each elementary $i$-th loop:

[TABLE]

The second type of the equations is obtained from the conservation of current in each node, for $p$-th node:

[TABLE]

where $I_{j}$ is the current in the $j$-th wire. The sum in Eq. (S4) is taken over the wires connected to the $p$-th node, respecting the chosen orientation of the wires in the network. There are $N_{\mathrm{nodes}}-1$ independent equations of the second type.

If the wires forming the network are long in comparison to the superconducting coherence length $\xi$, the pairbreaking effects are small, thus further we assume homogeneous superfluid density described by penetration depth $\lambda$. Under this condition, the current $I_{j}$ can be expressed through $\gamma_{j}$ using London equation as:

[TABLE]

where $l_{j}$ and $S_{j}$ are the length and the cross sectional area of the wire $j$. It is convenient to introduce the kinetic inductance $L_{j}^{K}$ of a wire with length $l_{j}$ and cross-sectional area $S_{j}$ as follows:

[TABLE]

Combining equations (S5) and (S6) yields

[TABLE]

Replacing $I_{j}$’s in equation(S4) we obtain

[TABLE]

Eq. (S8), written for $(N_{\mathrm{nodes}}-1)$ nodes of the network, together with Eq. (S3), written for $N_{\mathrm{loops}}$ loops form a system of $N_{\mathrm{links}}$ linear equations. Solving this system allows us to find the gauge-invariant phase differences $\gamma_{k}$ across all the wires of the network. From equations (S8) and (S3), it is easy to see, that $\gamma$’s depend linearly on $(\hat{n}-\hat{\phi})$.

For a wire of small cross-section ($S\ll\lambda^{2}$), the energy $E^{M}$ due to the magnetic field is negligible in comparison to the kinetic energy of the supercurrent $E^{K}$. Thus, the energy of the $j$-th wire is

[TABLE]

The total energy of the superconducting network is obtained by summing the energies of all links:

[TABLE]

Since the $\gamma$’s depend linearly on $(\hat{n}-\hat{\phi})$, the total energy of the network is a quadratic form of $(\hat{n}-\hat{\phi})$.

A vortex state in a ring with a crossbar (Fig. S3) is described by a pair of winding numbers $\{n_{1},n_{2}\}$ for the top and the bottom halves of the rings. Writing down equations S3, S8 for this structure (using the orientation of wires shown in Fig. S3) yields a system of three equations:

[TABLE]

where $\tilde{\phi_{i}}=\phi_{i}-n_{i}$. The total energy of the structure is obtained from Eq. (S10):

[TABLE]

where, for convenience, we introduce $1/L_{\|}^{K}=1/L^{K}_{1}+1/L^{K}_{2}+1/L_{3}^{K}$. In the case of a symmetric structure, such that $L_{1}^{K}=L_{2}^{K}=L^{K}$, the energy of the structure is

[TABLE]

where $\beta=L^{K}_{3}/(L^{K}_{1}+L^{K}_{3})$ is a parameter that characterises the strength of coupling of the two halves. For a ring with a crossbar made of wires with a homogeneous cross section, shown in Fig. S3, $\beta=2/(2+\pi)\approx 0.39$. It is instructive to consider two special limits: if $\beta=0$, equation (S13) yields $E\propto\tilde{\phi_{1}}^{2}+\tilde{\phi_{2}}^{2}$ – the structure behaves like two isolated rings; if $\beta=1$: $E\propto(\tilde{\phi_{1}}+\tilde{\phi_{2}})^{2}$, which corresponds to a ”vanishing” crossbar, and the structures behaves like a single ring.

The stability region of the vortex state with vortex numbers $\{n_{1},n_{2}\}$ has a hexagonal shape in ($\phi_{1},\phi_{2}$) coordinates, set by the following constrains:

[TABLE]

The examples of honeycomb stability diagram for structures with different values of inter-loop coupling are shown in Fig. S4(a-c). Increasing the kinetic inductance of the crossbar increases the inter-loop coupling, which results in a squeezing of the hexagons in (1,1) direction. The depth of the modulation of the honeycomb cells $\beta$, measured as shown in Fig. 2, can be expressed in terms of the inter-loop coupling constant $\beta$:

[TABLE]

Thus, the coupling $\beta$ between two halves of the structure can be determined from the observed pattern of vortex transitions. In asymmetric structures $(L^{K}_{1}\neq L^{K}_{2})$, the hexagons in the stability diagram become oblique, as shown in Fig. S4(d-f).

II.3 Matching of the vortex transitions by simulating magnetic field of the tip

In order to quantitatively understand the pattern of vortex transition observed in Fig. 2 we simulate them using a model for magnetic tip. As a first approximation, we use the dimensions of the magnetic particle measured with SEM (600 x 600 x 1000 nm3) and assume the homogeneous distribution of the magnetic moment $m_{\|}=2.2\times 10^{-13}$ J/T, that is measured using cantilever magnetometry. Further, we calculate the fluxes $(\phi_{1},\phi_{2})$ as a function of the tip position $(x_{tip},y_{tip},h)$ and use Eqs. (S14) to find the positions of vortex transitions. Next, we use $m_{\|}$ and $h$ as tuning parameters of the model to match the observed transitions. A good match is obtained with $m_{\|}=2.43\times 10^{-13}$ J/T and $h=1.1$ $\mathrm{\mu}m$. Finally, we use the refined model of the tip to calculate the field distribution $B_{z}(x,y)$ and the modulation $\delta B_{z}(x,y)=(\partial B_{z}/\partial y)\cdot\delta y$ due to tip oscillation along $y$-axis with amplitude $x_{0}$=5 nm (see Figs.2 f-g).

II.4 Stochastic resonance imaging of the transition rates

Here we generalize the model of the interaction of MFM cantilever with thermally-activated transitions between vortex states in a ring (Ref. Polshyn et al. (2018)) to a more general case of a small superconducting network composed of narrow superconducting wires. We consider two vortex configurations a and b near the transition point at $\hat{\phi_{0}}$, where $F_{a}(\hat{\phi_{0}})=F_{b}(\hat{\phi_{0}})$, and assume that all other vortex states have much lower or higher energies and hence do not contribute to the dynamics. We assume that the flux modulation $\delta\hat{\phi}$ is sufficiently small so that $|F_{a}(\hat{\phi})-F_{b}(\hat{\phi}))|\lesssim k_{B}T$ and hence both vortex states have a substantial probability of being occupied. We also assume that the energy barrier for transition between states $a$ and $b$ is sufficiently low to permit thermally activated transitions.

In the regime described above, the state of the structure exhibit thermally driven fluctuations between the two lowest-energy vortex states of the ring ($a$ and $b$). Consequently, the supercurrents $\hat{I}(t)$ circulating in the structure also contain a two-level stochastic component. The dynamics of the system is governed by $\hat{\phi}$-dependent transition rates $\Gamma_{a}$ and $\Gamma_{b}$ that correspond to transitions $a\to b$ and $b\to a$. The probability to find the structure in state $a$, when it is in thermal equilibrium, and the cantilever is stationary, is given by

[TABLE]

The dynamics of the probability $P_{a}(t)$ is determined by the relaxation rate $\nu_{r}=\Gamma_{a}+\Gamma_{b}$

[TABLE]

At $\hat{\phi}=\hat{\phi_{0}}$: $\Gamma_{a}=\Gamma_{b}$, so that $P^{eq}_{a}(\hat{\phi_{0}})=0.5$ and $\nu_{r}=2\Gamma_{a}$.

The supercurrent in the structure exerts a force on the magnetic particle given by ${\zeta(t)=\hat{\kappa}(\mathbf{r}_{tip})\cdot\hat{I(t)}}$, where $\hat{\kappa}(\mathbf{r}_{tip})$ represents the coupling between the circulating supercurrent and the cantilever. The $j$-th component of $\hat{\kappa}(\mathbf{r}_{tip})$ represents the coupling to the current $I_{j}$ flowing in the $j$-th wire of the structure. The equation of motion for the cantilever becomes:

[TABLE]

where $x$ is the displacement of the tip from its equilibrium position, $\gamma_{0}$ is the unmodified dissipation of the cantilever, and $f(t)$ is the force applied by the feedback controller, which resonantly drives the cantilever at a fixed amplitude $x_{0}$.

Tip oscillations with amplitude $x_{0}$ generate a small modulation of the flux $\delta\hat{\phi}=(\partial\hat{\phi}/\partial x)\,x_{0}$, which in turn modulates the transition rates $\Gamma_{a}$ and $\Gamma_{b}$. The modulation of the transition rates leads to statistical correlation between $\zeta(t)$ and $x(t)$, which causes the shift of the resonant frequency and dissipation of the cantilever. The motion of cantilever can be represented as a sum of the coherent and stochastic terms: $x(t)=x_{0}e^{i\omega t}+x_{s}(t)$. The stochastic part of the motion has a vanishing time time-averaged Fourier component $\langle\hat{x_{s}}(\omega)\rangle=0$. Since we are mainly interested in the effect of the fluctuating force on the frequency and dissipation of the cantilever $\Delta\omega\equiv\left<\omega-\omega_{0}\right>$ and $\Delta\gamma\equiv\left<\gamma-\gamma_{0}\right>$, emerging due to correlations between $\zeta(t)$ and $x(t)$, we consider only the case of weak stochastic force, i.e., $|\zeta(t,x)|\ll kx_{0}$. In this case $|x_{s}(t)|\ll x_{0}$, and we can approximate $x(t)\simeq x_{0}e^{i\omega t}$ in calculating $\langle\hat{\zeta}(\omega)\rangle$. The crucial approximation of weak $|\zeta(t,x)|$, which enables us to effectively decouple the cantilever dynamics from the dynamics of the phase slips, is justified in our measurements, because the observed resonance frequency shifts are small.

For sufficiently weak flux modulation such that $(\partial P^{eq}_{a}/\partial\hat{\phi})\delta\hat{\phi}\ll 1$, the resulting modulation of $P_{a}(t)$ is linear in $\delta\hat{\phi}$, with

[TABLE]

The ensemble-averaged value of current is given by

[TABLE]

Expanding $\hat{I}$ around $\hat{\phi_{0}}$, $\hat{I}(\hat{\phi})=\hat{I}(\hat{\phi_{0}})+(\partial\hat{I}/\partial\hat{\phi})\cdot\delta\hat{\phi}$ we obtain

[TABLE]

The first term in $\langle I(t)\rangle$ describes the contribution to the current from the thermally-activated transitions between the two states. The second term is not relevant to the effect of interest since it describes the flux dependence of the currents in each state. From Eq. S23, we find the Fourier component of the statistically-synchronized stochastic force due to the cantilever-driven phase slips $\langle\hat{\zeta}(\omega)\rangle=-\hat{\kappa}(\mathbf{r}_{tip})\Delta\hat{I}\delta Pe^{-i\theta}$, where $\Delta\hat{I}=\hat{I_{a}}(\hat{\phi_{0}})-\hat{I_{b}}(\hat{\phi_{0}})$. Finally, $\langle\hat{\zeta}(\omega)\rangle$ enables us to find $\Delta\omega$ and $\Delta\gamma$:

[TABLE]

where $K$ is a coupling constant given by

[TABLE]

First factor in Eq. S26 represents the geometric part of the coupling and depends on the mutual position of the structure and the cantilever, while the second factor describes the dependence of $K$ on the strength of supercurrents in the structure.

Remarkably, the ratio $\Delta\omega/\Delta\gamma$, obtained from Eqs. (S24) and (S25), has a simple form and gives a relaxation rate $\nu_{r}$ in units of cantilever frequency:

[TABLE]

Fig. S5 shows data obtained for vortex transitions between states $\{5,5\}$, $\{5,4\}$ and $\{6,4\}$. The frequency shift data shown in Fig. S5b were measured by taking scans along the lines marked as A, B, C in Fig. S5a. All three transitions show no dependence on the cantilever oscillation amplitude except at the highest amplitude of 12 nm. Moreover, the relative frequency shifts are small ($\Delta\omega/\omega_{0}<0.002$). Thus, the modulation is indeed sufficiently weak to justify the assumptions of our model of driven vortex transitions. Hence Eq. (S27) can be used to extract the relaxation rates of the transitions.

The peak shifts in the resonant frequency and the dissipation, which were extracted from the line scans, are shown in figure S5c as a function of temperature. As can be seen from Fig. S5c, the signal on all three vortex transitions shows behavior consistent with stochastic resonance model given by Eqs. (S24) and (S25): the dissipation $\Delta\gamma$ reaches a maximum (when $\nu_{r}\sim\omega_{0}$), while the resonant frequency shift $\Delta\omega$ grows rapidly; moreover, at this point $-\Delta\omega\simeq\Delta\gamma$. The $\nu_{r}\sim\omega_{0}$ condition for transition C is observed around 1.175 K, for transition A – around 1.173 K, and for transitions B at 1.170 K. The rates $\nu_{r}$, derived for each transition using equation Eq. (S27), are plotted in figure S5d. The right axis shows the relative change of the phase slip barrier heights in units of $k_{B}T$, which was calculated under the assumption that the attempt frequency does not change significantly within the temperature range of the measurements: $-\Delta E/k_{B}T\simeq\ln(\nu_{r})$. Figure S5e provides full images of $\nu_{r}$ for several temperatures in addition to those shown in Fig. 3 of the main text. The images of the phase slip rate demonstrate a unique way $\Phi_{0}$-MFM allows us to characterize the dynamics of individual vortex transitions even in complex structures.