Controllable Vortex Loops in Superconducting Proximity Systems

TL;DR

This paper proposes a method to create and control superconducting vortex loops in a normal metal via the proximity effect, enabling their experimental detection through scanning tunneling microscopy.

Contribution

It introduces a novel technique to generate and tune vortex loops in superconducting proximity systems, with analytical and numerical validation.

Findings

Vortex loops can be generated by piercing the junction with a current-carrying wire.

The radii of vortex loops can be analytically expressed and tuned.

Vortex loops can be made large enough to be observed directly with microscopy.

Abstract

Superconducting vortex loops have so far avoided experimental detection despite being the focus of much theoretical work. We here propose a method of creating controllable vortex loops in the superconducting condensate arising in a normal metal through the proximity effect. We demonstrate both analytically and numerically that superconducting vortex loops emerge when the junction is pierced by a current-carrying insulated wire and give an analytical expression for their radii. The vortex loops can readily be tuned big enough to hit the sample surface, making them directly observable through scanning tunneling microscopy.

Equations

$$D\bar{\nabla }\cdot \left(\check{g}\bar{\nabla }\check{g}\right)+i\left[\epsilon {\hat{\rho }}_3+\hat{\Delta },\check{g}\right]=0.$$

$$\check{g}=\left(\begin{array}{cc}{\hat{g}}^{\text{textsc}r}& {\hat{g}}^{\text{textsc}k}\\ 0& {\hat{g}}^{\text{textsc}a}\end{array}\right)$$

$${\zeta }_i{L}_i{\mathbf{e}}_n\cdot ({\hat{g}}_i^{\text{textsc}r}\bar{\nabla }{\hat{g}}_i^{\text{textsc}r})=\frac{1}{2}\left[{\hat{g}}_i^{\text{textsc}r},{\hat{g}}_j^{\text{textsc}r}\right],$$

$$\mathbf{B}=\frac{\mu }{4\pi }\int \frac{\mathbf{J}({\mathbf{r}}^{\prime })\times (\mathbf{r}-{\mathbf{r}}^{\prime })}{\text{absolutevalue}{\mathbf{r}-{\mathbf{r}}^{\prime }}^3}\text{differential}[3]r^{\prime },$$

$$\mathbf{J}=I\delta (y)\delta (z)\left[\theta (x+L/2)-\theta (x-L/2)\right],$$

$$\mathbf{B}=\frac{\mu I}{4\pi \rho }\left(\frac{L/2+x}{\sqrt{(L/2+x{)}^2+{\rho }^2}}+\frac{L/2-x}{\sqrt{(L/2-x{)}^2+{\rho }^2}}\right){\mathbf{e}}_{\varphi }$$

$$\bm{A}=\frac{\mu I}{4\pi}\ln\left(\frac{\sqrt{(L/2-x)^{2}+\rho^{2}}+L/2-x}{\sqrt{(L/2+x)^{2}+\rho^{2}}-L/2-x}\right)\Bigl{[}\theta(x+L/2)-\theta(x-L/2)\Bigr{]}\bm{e}_{x},$$

$${\hat{g}}^{\text{textsc}r}=\left(\begin{array}{cc}N& 0\\ 0& -\tilde{N}\end{array}\right)\left(\begin{array}{cc}1+\gamma \tilde{\gamma }& 2\gamma \\ 2\tilde{\gamma }& 1+\tilde{\gamma }\gamma \end{array}\right),$$

$$\begin{array}{r}{\nabla }^{2}a=\frac{2\tilde{a}\nabla a\cdot \nabla a}{1+a\tilde{a}}+\frac{4(1-a\tilde{a})Le\mathbf{A}\cdot (aLe\mathbf{A}+i\nabla a)}{1+a\tilde{a}}\\ +2iLe(\nabla \cdot \mathbf{A})a-2i\epsilon a,\end{array}$$

$$\begin{array}{r}{\mathbf{e}}_n\cdot \nabla a=\frac{(1+a\tilde{b})(b-a)}{\zeta (b\tilde{b}+1)}+2ia{\mathbf{e}}_n\cdot \mathbf{A}eL,\end{array}$$

$${\hat{g}}^{\text{textsc}r}=\left(\begin{array}{cc}g& f\\ -\tilde{f}& -\tilde{g}\end{array}\right).$$

$$N_{\sigma }(\epsilon ,\mathbf{r})=N_0\Re \{g_{\sigma \sigma }(\epsilon ,\mathbf{r})\},$$

$$N(\epsilon ,\mathbf{r}):=\frac{N_{\uparrow }(\epsilon ,\mathbf{r})+N_{\downarrow }(\epsilon ,\mathbf{r})}{2}=N_0\frac{1-a\tilde{a}}{1+a\tilde{a}}.$$

$$\mathbf{j}=\frac{N_{0}eD}{4}{\int }_{-\infty }^{\infty }\text{Tr}\left({\hat{\rho }}_3{\left[\check{g}\bar{\nabla }\check{g}\right]}^{\text{textsc}k}\right)\text{differential}\epsilon .$$

$$\displaystyle\begin{split}\bm{j}=\frac{N_{0}eD}{2}\int_{-\infty}^{\infty}\tanh\left(\frac{\beta\varepsilon}{2}\right)\Tr\Bigl(\operatorname{Re}\left[{\tilde{f}}^{\dagger}\nabla{f}^{\dagger}-f\nabla\tilde{f}\right]{\\ }+2e\bm{A}\operatorname{Im}\left[f\tilde{f}-{\tilde{f}}^{\dagger}{f}^{\dagger}\right]\Bigr{missing})\differential{\varepsilon}.\end{split}$$

$$\begin{gathered} \Psi (\mathbf{r}):=\text{expectationvalue}{\psi }_{\uparrow }(\mathbf{r},0){\psi }_{\downarrow }(\mathbf{r},0) \\ =\frac{N_0}{2}{\int }_{-\infty }^{\infty }{f}_s(\mathbf{r},\epsilon )\tanh (\epsilon \beta /2)\text{differential}\epsilon . \end{gathered}$$

$${\hat{g}}_{\text{textsc}bcs}=\left[\frac{\theta \left({\epsilon }^2-\text{absolutevalue}{\Delta }^2\right)}{\sqrt{{\epsilon }^2-\text{absolutevalue}{\Delta }^2}}\mathrm{sgn}(\epsilon )-\frac{\theta \left(\text{absolutevalue}{\Delta }^2-{\epsilon }^2\right)}{\sqrt{\text{absolutevalue}{\Delta }^2-{\epsilon }^2}}i\right]\left(\epsilon {\hat{\rho }}_3+\hat{\Delta }\right),$$

$$D_{\text{textsc}sc}\nabla \cdot \left({\hat{g}}^{\text{textsc}r}\nabla {\hat{g}}^{\text{textsc}r}\right)+i\left[\epsilon {\hat{\rho }}_3+\hat{\Delta },{\hat{g}}^{\text{textsc}r}\right]=0$$

$${\hat{g}}^{\text{textsc}r}={\hat{g}}_{\text{textsc}bcs}+\delta \hat{g}$$

$$D_{\text{textsc}sc}\nabla \cdot \left(\left[{\hat{g}}_{\text{textsc}bcs}+\delta \hat{g}\right]\nabla \delta \hat{g}\right)+i\left[\epsilon {\hat{\rho }}_3+\hat{\Delta },\delta \hat{g}\right]=0,$$

$$D_{\text{textsc}sc}{\nabla }^2\delta \hat{g}+i{\hat{g}}_{\text{textsc}bcs}\left[\epsilon {\hat{\rho }}_3+\hat{\Delta },\delta \hat{g}\right]=0.$$

$${\left({\hat{g}}_{\text{textsc}bcs}+\delta \hat{g}\right)}^2=1⟹\left\{{\hat{g}}_{\text{textsc}bcs},\delta \hat{g}\right\}=0.$$

$${\hat{g}}_{\text{textsc}bcs}\left[\epsilon {\hat{\rho }}_3+\hat{\Delta },\delta \hat{g}\right]$$

$$=\left\{\delta \hat{g},(\epsilon {\hat{\rho }}_3+\hat{\Delta }){\hat{g}}_{\text{textsc}bcs}\right\}.$$

$${\left(\epsilon {\hat{\rho }}_3+\hat{\Delta }\right)}^2={\epsilon }^2-{\Delta }^2$$

$${\nabla }^2\delta \hat{g}={\lambda }^{-2}\delta \hat{g}$$

$$\displaystyle\lambda^{-2}=-\frac{2i}{D_{\textsc{sc}}}\Biggl{[}\operatorname{sgn}(\varepsilon)\sqrt{\varepsilon^{2}-\absolutevalue{\Delta}^{2}}\theta\left(\varepsilon^{2}-\absolutevalue{\Delta}^{2}\right)$$

$$\displaystyle+i\sqrt{\absolutevalue{\Delta}^{2}-\varepsilon^{2}}\theta\left(\absolutevalue{\Delta}^{2}-\varepsilon^{2}\right)\Biggr{]}.$$

$$\nabla\delta\hat{g}\bigr{\rvert}_{\bm{r}\in\Omega}=0,$$

$$\delta \hat{g}(\epsilon ,x,y)=C\left[{\mathrm{e}}^{-\text{absolutevalue}x+1/2/\lambda }+{\mathrm{e}}^{-2L_{\text{textsc}sc}/\lambda +\text{absolutevalue}x+1/2/\lambda }\right],$$

$$C=\frac{\lambda {\hat{g}}_{\text{textsc}bcs}\left[{\hat{g}}_{\text{textsc}bcs}+\delta \hat{g},{\hat{g}}_{\text{textsc}n}\right]}{2\left(1-{\mathrm{e}}^{-2L_{\text{textsc}sc}/\lambda }\right){\zeta }_{\text{textsc}sc}{L}_{\text{textsc}sc}}.$$

Full text

I Introduction

Many key properties of physical systems are determined by topological defects such as dislocations in solids, domain walls in ferroics, vortices in superfluids, magnetic skyrmions in condensed matter systems and cosmic strings in quantum field theories. In superconductors, the topological entities are vortex lines of quantized magnetic flux. The topological nature of these vortices makes them stable, which is important for potential applications such as superconducting qubits Fedorov et al. (2014); Devoret and Schoelkopf (2013); Herr et al. (2007), digital memory and long-range spin transport Kim et al. (2018). Vortices have non-superconducting cores and a phase winding of an integer multiple of $2\pi$ in the superconducting order parameter, leading to circulating supercurrents Kwok et al. (2016).

The formation of superconducting vortex loops is topologically allowed, and has theoretically been predicted to form around strong magnetic inclusions inside the superconductor Doria, M. M. and Romaguera, A. R.De C. and Milošević, M. V. and Peeters, F. M. (2007), in cylindrically shaped current-carrying superconductors Samokhvalov (1998, 1997); Kozlov and Samokhvalov (1993) or through vortex cutting and recombination Glatz et al. (2016); Berdiyorov et al. (2018). However, no observation of vortex loops in superconducting systems has been found to date. One challenging aspect is that vortex loops are typically small in conventional superconductors and difficult to stabilize for an extended period of time Schönenberger et al. (1996). Recently it has been shown that vortex loops can be formed in proximity systems by inserting physical barriers, around which the vortices can wrap Berdiyorov et al. (2018).

In this manuscript, we present a way to create controllable vortices in mesoscopic proximity systems in a manner which makes them experimentally detectable through scanning tunneling microscopy. The system considered is a three-dimensional SNS junction pierced by a current-carrying wire which creates the inhomogeneous field responsible for the vortex loops. In planar SNS-junctions with uniform applied magnetic field, changing the superconducting phase difference between the two superconductors shifts the vortex lines in the vertical direction Cuevas and Bergeret (2007). We here show that the corresponding effect on vortex loops in three dimensions is to change their size. Thus, these vortex loops are easily tunable. This makes it possible to make the vortices touch the surface, leaving distinct traces which are directly observable by scanning tunneling spectroscopy Stolyarov et al. (2018).

Vortex loops in superconducting systems has previously been predicted using the phenomenological Ginzburg-Landau theory Berdiyorov et al. (2018); Glatz et al. (2016); Doria, M. M. and Romaguera, A. R.De C. and Milošević, M. V. and Peeters, F. M. (2007). Here we use a fully microscopic framework known as quasiclassical theory and solve the Usadel equation relevant for diffusive systems Usadel (1970). By showing that vortex loop formation occurs in a microscopic theory, we give valuable support to the earlier proposed mechanisms for superconducting vortex loops. Finally, we discuss how the proposed setup can be realized experimentally.

II Methodology

In this section we discuss the quasiclassical Usadel theory and how it may be used to analyse the SNS-junction depicted in fig. 1. We first present the mathematical tools and end with the numerical implementation.

II.1 Quasiclassical Theory

In the Usadel theory, the system is described by a quasiclassical Green’s function from which physical properties can be extracted. The SNS junction depicted in fig. 1 can be treated in the quasiclassical formalism under the assumptions that the Fermi wavelength is much shorter than all other relevant length scales. If the system is diffusive, meaning that the scattering time is small, the isotropic part dominates and solves the Usadel equation Zaikin et al. (2002); Chandrasekhar (2004); Rammer (2004); Usadel (1970), which in the normal metal can be written

[TABLE]

Here, $D$ is a diffusion constant, $\hat{\rho}_{3}=\operatorname{diag}(1,1,-1,-1)$ and $\hat{\Delta}=\operatorname{antidiag}(+\Delta,-\Delta,+\Delta^{*},-\Delta^{*})$ where $\Delta$ is the superconducting gap parameter. The covariant derivative is $\bar{\nabla}\check{g}=\nabla\check{g}-ie\left[{\hat{\rho}_{3}{\bm{A}}}\,,\,{\check{g}}\right]$, where $e=-\absolutevalue{e}$ is the electron charge, $\bm{A}$ is the vector potential and

[TABLE]

is the quasiclassical impurity-averaged Green’s function. Finally, $(x,y,z)\in[-L/2,L/2]\times[-W/2,W/2]\times[-H/2,H/2]$ in the normal metal.

$\check{g}$ is normalized such that $\check{g}\check{g}=1$. We use the convention that when two matrices of different dimensionality is multiplied, the smaller matrix is elevated to the dimensionality of the larger matrix by the tensor product of an identity matrix of the appropriate size. In equilibrium, the components of the $8\times 8$ Green’s function in eq. 2 are related by the identities $\hat{g}^{\textsc{k}}=\left(\hat{g}^{\textsc{r}}-\hat{g}^{\textsc{a}}\right)\tanh(\varepsilon\beta/2)$ and $\hat{g}^{\textsc{a}}=-\hat{\rho}_{3}\hat{g}^{\textsc{r}\dagger}\hat{\rho}_{3}$, which means that in this case it is sufficient to solve for the retarded component $\hat{g}^{\textsc{r}}$.

The quasiclassical formalism is not applicable across boundaries because the associated length scale is too short. The Usadel equation must therefore be solved in the normal metal and superconductors separately, and the solutions must be connected through boundary conditions. If we assume a low-transparency interface, we may use the Kupriyanov-Lukichev boundary condition

[TABLE]

where $\bm{e}_{n}$ is the outward-pointing normal vector for region $i$, $\zeta_{i}$ is the ratio of the bulk and interface conductances of material $i$ and $L_{i}$ is the length of material $i$ in the direction of $\bm{e}_{n}$. For the boundaries interfacing vacuum, $\bm{e}_{n}\cdot\bar{\nabla}\hat{g}^{\textsc{r}}=0$.

The Usadel equation can be made dimensionless by introducing the Thouless energy, $\varepsilon_{\textsc{t}}\coloneqq D/L^{2}$. The Usadel equation then becomes dimensionless by doing the substitutions $(x,y,z)\to(x/L,y/L,z/L)$, $\varepsilon\to\varepsilon/\varepsilon_{\textsc{t}}$, $\Delta\to\Delta/\varepsilon_{\textsc{t}}$ and $\bar{\nabla}\to L\bar{\nabla}$.

II.2 Electromagnetic vector potential

The magnetic field should satisfy Biot-Savart’s law,

[TABLE]

where $\mu$ is the permeability and $\bm{J}$ is the electric current density. In general, $\bm{J}$ includes the contribution the induced currents in the normal metal and superconductors in addition to that from the insulated current-carrying wire along the $x$-axis. However, we will make some assumptions in order to simplify the analytical and numerical calculations. Firstly, we will assume that the width $W$ and height $H$ is smaller than the Josepshon penetration depth. In this case we can ignore the screening of the magnetic field by currents inside the normal metal Barone and Paternò (1982). Secondly, we will neglect the magnetic field produced by the supercurrents produced inside the normal metal. Thirdly, we will assume that the magnetic field vanish inside the superconductors due to the screening currents. These last two assumptions is widely used in the context of hybrid structures with constant applied magnetic fields Cuevas and Bergeret (2007); Alidoust and Halterman (2015); Bergeret and Cuevas (2008), and has in such conditions been shown to give good agreement with experimental results Chiodi et al. (2012). Finally, we will assume that the part of the wire which is inside the superconductors does not contribute to the magnetic field in the normal metal.

The last two assumptions are inaccurate very close to the wire. Close to the wire the details of screening currents will be important for the magnetic field, but far away we assume that the total contribution from the currents inside the superconductor is zero. A more precise model could be developed by taking into account screening currents inside the superconductors and solving the Usadel equation self-consistently with Maxwell’s equation and the superconducting gap-equation inside the superconductors. However, we are here interested in the solution far away from the wire, and the details of the magnetic field near the wire should not significantly alter the results. For this reason we also model the wire as being infinitely thin.

With the assumptions presented above we get a current density which is

[TABLE]

where $\theta$ is the Heaviside step function. Inserting eq. 5 into eq. 4 we get

[TABLE]

for $x\in(-L/2,L/2)$, where $\rho=\sqrt{y^{2}+z^{2}}$ and $\bm{e}_{\phi}=(y\bm{e}_{z}-z\bm{e}_{y})/\rho$. $\bm{B}=0$ for $x<-L/2$ and $x>L/2$. A vector potential which satisfies $\bm{B}=\nabla\times\bm{A}$ is

[TABLE]

as can be seen from insertion or calculated directly from Biot-Savarts law by using that $\nabla\times(\bm{J}(\bm{r}^{\prime})/\absolutevalue{\bm{r}-\bm{r}^{\prime}})=\bm{J}(\bm{r}^{\prime})\times(\bm{r}^{\prime}-\bm{r})/\absolutevalue{\bm{r}^{\prime}-\bm{r}}^{3}$.

II.3 The Ricatti Parametrization

In the Ricatti parametrization Schopohl (1998) of $\hat{g}^{\textsc{r}}$, the parameter is the $2\times 2$ matrix $\gamma$ and the retarded Green’s function is written

[TABLE]

where $N\coloneqq\left(1-\gamma\tilde{\gamma}\right)^{-1}$ and tilde conjugation is $\tilde{\gamma}(\varepsilon)={\gamma}^{*}(-\varepsilon)$.

Since the superconducting correlations in our system are spin-singlet, we may write $\gamma_{\textsc{n}}=\operatorname{antidiag}(a,-a)$ and $\gamma_{\textsc{bcs}}=\operatorname{antidiag}(b,-b)$, where $\gamma_{\textsc{n}}$ and $\gamma_{\textsc{bcs}}$ are the Ricatti parameters in the normal metal and superconductors, respectively. Substituting this into eqs. 1 and 3 we obtain the dimensionless equations

[TABLE]

and

[TABLE]

where $L$ is the length which is used to define the Thouless energy, $\varepsilon_{\textsc{t}}$. The corresponding equations for $\tilde{a}$ and $\bm{e}_{n}\cdot\nabla\tilde{a}$ is found by tilde conjugating eqs. 9 and 10.

II.4 Observables

As mentioned initially, a vortex is accompanied by a non-superconducting core and a circulating supercurrent. Both the superconducting order parameter and the supercurrent can be extracted from the quasiclassical Green’s function. In the following it will be useful to write

[TABLE]

There are only singlet correlations in the SNS system, so $f=\operatorname{antidiag}(f_{s},-f_{s})$.

The local density of states for spin-band $\sigma$ at energy $\varepsilon$ and location $\bm{r}$ can be written

[TABLE]

where $N_{0}$ is the normal state density of state at the Fermi surface. In the normal metal we can write eq. 12 in terms of $a$,

[TABLE]

In the cores of vortices we expect $N=N_{0}$ for all energies, which happens when $a(\varepsilon)\equiv 0$.

The current density is Zaikin et al. (2002)

[TABLE]

Inserting eq. 11, using the relations $\hat{g}^{\textsc{a}}=-\hat{\rho}_{3}\hat{g}^{\textsc{r}\dagger}\hat{\rho}_{3}$, $\hat{g}^{\textsc{k}}=\left(\hat{g}^{\textsc{r}}-\hat{g}^{\textsc{a}}\right)\tanh(\varepsilon\beta/2)$, eq. 14 can be rewritten

[TABLE]

Written in terms of the quasiclassical Green’s function, the superconducting order parameter is

[TABLE]

where $\psi_{\sigma}(\bm{r},t)$ is the field operator which destroys an electron with spin $\sigma$ at position $\bm{r}$ and time $t$, $N_{0}$ is the normal state density of states and $\beta=1/k_{\textsc{b}}T$.

II.5 Numerics

The Usadel equation was solved numerically using a finite element scheme. See for instance Amundsen and Linder (2016) to see how to set up solve the nonlinear Usadel equations in a finite element scheme by the use of the Newton-Rhapson method. The program was written in Julia Bezanson et al. (2017), we used linear hexehedral elements and JuAFEM.jl Carlsson (2019) was used to iterate through the cells. Gauss-Legandre quadrature rules of fourth order was used to integrate through the cells and Romberg integration was used to integrate over energy. See for instance Sauer (2013). Finally, forward-mode automatic differentiation Revels et al. (2016) was used to calculate the Jacobian.

III Results and Discussion

Here we present first an analytical solution of the Usadel equation in the weak proximity effect regime, then we show numerically that the findings are also present also in the full proximity effect regime. Dimensionless quantities are used in the analytics with distances being measured relative to the length of the half-metal, $L$, and energies being measures relative to the Thouless energy $\varepsilon_{\textsc{t}}=D/L^{2}$, where $D$ is the diffusion constant in the half-metal.

III.1 Analytics

Before solving the Usadel equation we must determine the solution in the superconducors. We will show that it suffices to use the bulk solution

[TABLE]

in the superconductors when a certain condition is fulfilled. Let $\lambda$ (to be defined quantitatively below) be the length-scale over which the Green function recovers its bulk value in the superconductor. The criterion for neglecting the inverse proximity effect in the superconductors is then that the normal-state conductance of the superconductors for a sample of length $\lambda$ is much larger than the interface conductance and that the length of each superconductor is not small compared to $\lambda$. We now proceed to prove this.

The vector potential, (7), is zero inside the superconductors, so the Usadel equation simplifies to

[TABLE]

in the superconductor at $x<-1/2$. To show that we can use the bulk solution in the limit $L_{\textsc{sc}}\to\infty$, let

[TABLE]

This gives an equation for $\delta\hat{g}$,

[TABLE]

where we have used that $\hat{g}_{\textsc{bcs}}$ solves the eq. 18 for a bulk superconductor. Next, assume the inverse proximity effect to be weak, such that $\delta\hat{g}\ll\hat{g}_{\textsc{bcs}}$. Using that $\hat{g}_{\textsc{bcs}}\hat{g}_{\textsc{bcs}}=1$, this yields

[TABLE]

$\hat{g}_{\textsc{bcs}}+\delta\hat{g}$ must also satisfy the normalization condition $(\hat{g}^{\textsc{r}})^{2}=1$, so

[TABLE]

Hence, using that $\left[{\varepsilon\hat{\rho}_{3}+\hat{\Delta}}\,,\,{\hat{g}_{\textsc{bcs}}}\right]=0$,

[TABLE]

Finally, from

[TABLE]

we get that $\delta\hat{g}$ is an eigenfunction of the Laplacian,

[TABLE]

where

[TABLE]

We can choose the sign of $\lambda$ to be such that $\Re(\lambda)>0$.

Let $L_{\textsc{sc}}$ be the length of the superconductor in multiples of the length of the normal metal. Using the boundary condition

[TABLE]

where $\Omega$ is the boundary not interfacing the normal metal, we get

[TABLE]

where $C$ is some a function of $y$ and $\varepsilon$ to be determined by the final boundary condition. From the remaining boundary condition, eq. 3, we get

[TABLE]

From eq. 28 we see that $\Re(\lambda)$ can be interpreted as the penetration depth of $\delta g$. Note that $\Re(\lambda)$ is bounded by including the effect of inelastic scattering, which is done by the substitution $\varepsilon\to\varepsilon+i\delta$ for some positive scattering rate $\delta$ Dynes et al. (1984). This ensures that $1/\left(1-\mathrm{e}^{-2L_{\textsc{sc}}/\lambda}\right)$ remains finite as $\varepsilon\to\Delta$. Thus, we see from eq. 29 that $C$, and therefore $\delta g$, becomes negligble when

[TABLE]

provided that the length of the superconductor $L_{\textsc{sc}}$ is not small compared to the maximal penetration depth, $\max[\Re(\lambda)]$.

$\zeta_{\textsc{sc}}$ is proportional to the conductance of the whole superconductor and therefore also with $1/L_{\textsc{sc}}$. Therefore, $\zeta_{\textsc{sc}}L_{\textsc{sc}}/\Re(\lambda)$ is the ratio of the normal-state conductance of a superconductor of length $\Re(\lambda)$ to the interface conductance. Taking the superconducting coherence length $\xi$ as a measure of the inverse proximity effect penetration depth $\Re(\lambda)$, we see that the criterion eq. 30 is indeed experimentally feasible. The equation is fulfilled for a low-transparency interface and for a superconductor that is larger than the coherence length. A similar calculation shows that we can use $\hat{g}_{\textsc{bcs}}$ also in the superconductor at $x>1/2$.

Solving for the Ricatti parameter in the superconductors we get that $\gamma_{\textsc{bcs}}=\operatorname{antidiag}(b,-b)$ with

[TABLE]

The non-linear Usadel equation does not have a general analytical solution, but it can be solved analytically in an approximate manner far away from the wire. If we assume the proximity effect to be weak, we can keep only terms which are linear in $a$, $\tilde{a}$ and their gradients. In this case the Usadel equation (9) decouples:

[TABLE]

Equation 32 can be further simplified when we only consider regions where $\rho\gg 1$, with $\rho=\sqrt{y^{2}+z^{2}}$. The solution of eq. 32 is constant in $y$ and $z$ when $\bm{A}=\bm{0}$, and by assuming this is also approximately true when $\absolutevalue{eL\bm{A}}\ll 1$, we can neglect the terms $\partial_{y}^{2}a$ and $\partial_{z}^{2}a$. Finally, we can simplify the calculations further by Taylor expanding the vector potential,

[TABLE]

where

[TABLE]

We keep only the first term in the Taylor expansion.

Equation 32 can now be solved exactly, and by applying the linearized boundary conditions,

[TABLE]

the solution can be written on the form

[TABLE]

where

[TABLE]

From eq. 36 we see that $a$ vanishes at $x=0$ and $i\delta\phi-u=i(2N+1)\pi$, where $N$ is any integer. This happens at

[TABLE]

This means that $f$ and hence also $\Psi$ vanish at these points. By Taylor expanding $a$ to first order around a root located at $(0,\tilde{\rho})$ we find

[TABLE]

where $x\sim\cos\theta$ and $\rho-\tilde{\rho}\sim\sin\theta$, $B_{1}^{2}=5\absolutevalue{k}^{2}/4-\absolutevalue{k}d+2d^{2}$, $B_{2}^{2}=\absolutevalue{k}^{2}/4+d^{2}$, $\alpha_{1}=\tan^{-1}[(\absolutevalue{k}/2+d)/(\absolutevalue{k}-d)]$ and $\alpha_{2}=\tan^{-1}(\absolutevalue{k}/2d)$. Hence, these roots have a phase winding of $2\pi$, as is characteristic for vortices. Equation 40 is our main analytical result as it predicts how the radius of the vortex loops depends on the tunable parameters of the system: the current through the wire and the applied phase difference. Although it was obtained using approximations, we demonstrate below that it matches the full numerical solution of the exact Usadel equation very well.

Note that the radius, $\rho$, of the largest vortex loop given eq. 40 can be made arbitrary large by letting $\phi_{r}-\phi_{l}$ approach $\pi$. Thus, for a given sample size $L\times W\times H$ and current $I$, there is a superconducting phase difference for which the vortex loop hits the surface and can be directly detected experimentally.

It is expected that a change in the superconducting phase difference will change the radii of vortex loops. This is because changing the phase difference is equivalent to changing the applied supercurrent through the junction. The applied current will be deflected by the circulating currents associated with the vortices and hence produce a reactionary force on the vortices. See for instance Sonin (1997). What is more surprising, however, is that changing the superconducting phase difference can make the vortices arbitrarily large so that they can always be made to hit the surface. If this feature is gerenally true for other systems with vortex loops it could prove useful for the study of systems containing vortex loops which are less obviously controllable than the one considered in the present manuscript, but which are easier to design in a lab. For instance, one possibility is to grow the normal metal around a magnetic dipole. Ref. Doria, M. M. and Romaguera, A. R.De C. and Milošević, M. V. and Peeters, F. M. (2007) found that vortex loops can form around magnetic dipole inclusions in superconductors if the magnetic field is strong enough, so thre are reasons to believe that vortex loops can also form around magnetic dipoles embedded in a SNS-junction. The magnetic field from a dipole can, unlike the magnetic field from a wire, not be altered in strength. Nevertheless, if the field is strong enough to produce vortices, altering the superconducting phase difference could be a way to increase the size of the vortex to the point where it touches the surface and becomes directly observable.

III.2 Numerics

We now proceed to show numerical results in the full (non-linear) proximity effect regime. We have set the parameters $\absolutevalue{\Delta}=4\varepsilon_{\textsc{t}}$, $\zeta=3$, $W=H=6L$ and $\phi_{l}=0$ common for all the numerical calculations. We include the effect of inelastic scattering by doing the substitution $\varepsilon\to\varepsilon+i\delta$ where $\delta=0.001\absolutevalue{\Delta}$ in order to avoid the divergence of $\hat{g}_{\textsc{bcs}}$ at $\varepsilon=\absolutevalue{\Delta}$ Dynes et al. (1984).

Numerically we find that vortex loops do indeed form at the locations predicted by the analysis. There are circular paths around the origin where the superconducting order parameter vanish and the local density of states is equal to that of the normal state. This can be seen in figs. 2 and 3 which shows the local density of states and the amplitude of the Cooper pair correlation function $\Psi$, respectively. Around these loops there are a circulating supercurrent, as can be seen in fig. 4, and a phase winding in the order parameter of $2\pi$. Figure 5 shows a contour plot of $\absolutevalue{\Psi}$, which shows the location of the vortices, together with the circulating supercurrent $\bm{j}$ as well as the phase of $\Psi$, which shows that there is indeed a phase winding of $2\pi$ around the vortices.

We find that the positions of the vortex loops match with eq. 40 for vortices with radius that are between $2L$ and $3L$. Figures 6 and 7 shows how the sizes of the vortex loops depend in superconducting phase difference $\phi$ and magnetic field strength $n$, respectively. We find that increasing $\phi$ can make the vortices arbitrary large, but does not increase the number of vortices. Increasing $n$, on the other hand, also increase the number of vortices, but the sizes grow only linearly with $n$. Note that as the vortex loops hit the surface, they curve so as to hit normally to the surface. This is consistent with previous results Doria, M. M. and Romaguera, A. R.De C. and Milošević, M. V. and Peeters, F. M. (2007); de C. Romaguera et al. (2007), and can be understood from the circulating currents. There should be no current component normal to the surface, and the only way for the current circulating the vortices to adhere to this is if the vortices hit the surfaces at a right angle.

IV Experimental Realization

Normal SNS junctions are created by vertically growing first a superconducting material such as niobium, then a normal metal such as copper and finally the same superconducting material. The layers are grown for instance by a sputter deposition technique such as direct current sputtering Miller et al. (1993) or radio frequency sputtering Krasnov et al. (1994). The setup presented here adds an extra complication by requiring an isolated conducting nanowire to penetrate the system. One possible way to achieve this could be to first grow a vertical insulated nanowire and then grow the superconductor and normal metal around it in a layerwise fashion.

Growing a wire is more complicated than growing a plane because one must localize the growth to happen at the tip at the wire, even though most of the surface area will be on the sides. Nevertheless, growing vertical nanowires has successfully been done by methods such as the vapor-liquid-solid method Ng et al. (2004); Schmidt et al. (2006); Tomioka et al. (2012) and template-directed synthesis Xia et al. (2003). The vapor-liquid-solid method works by using droplets of, for instance, gold which are a few angstroms in width to localize the growth Schmidt et al. (2006), and temple-directed synthesis works by having the wire grow inside a premade template which can later be removed Xia et al. (2003). The vapor-liquid-solid method has already been used to produce vertical surround-gate field-effect transistors with a precision exceeding what should be necessary for the system presented here Schmidt et al. (2006).

Schmidt et al. (2006) made nanowires using the vapor-liquid-solid method which were $40\text{\,}\mathrm{nm}$ is diameter and $400\text{\,}\mathrm{nm}$ in length. This should be on the same length scale as would be necessary for the system considered in this manuscript. The superconducting energy gap of niobium is $\absolutevalue{\Delta}=$30.5\text{\times}{10}^{-4}\text{,}\mathrm{eV} Kittel *et al.* ([1976](#bib.bib38)), which is equivalent to about $2.46\text{\,}{\mathrm{mm}}^{-1}$ in natural units. The Fermi velocity and scattering time for copper are about $v_{F}=$3.70\text{\times}{10}^{-3}\text{\,} and $\tau=$10.8\text{,}\mathrm{\SIUnitSymbolMicro m}$$, respectively Gall (2016). The diffusion coefficient is defined as

[TABLE]

so the diffusion coefficient for copper is about $D=$49.2\text{,}\mathrm{pm}$$. In the numerics we have used

[TABLE]

so

[TABLE]

which is on the same scale as what has been made with the vapor-liquid-solid method. Of course, other metals and superconductors could be used, giving different physical lengths corresponding to the values being used in the numerics here. Moreover, from the analysis it seems vortex loops would form also for other values of $\absolutevalue{\Delta}/\varepsilon_{\textsc{t}}$. The calculation above is merely to show that the length scales used here is not unreasonable compared to what has already been experimentally achieved.

V Conclusion

We have used quasiclassical Usadel theory to demonstrate that controllable superconducting vortex loops can emerge in a Josephson junction pierced by an insulated current-carrying wire. The size and number of vortices depend on the phase difference between the superconducting order parameter in the superconductors, $\phi_{r}-\phi_{l}$, as well as the strength of the magnetic field. The radius of the vortices can be made arbitrarily large by tuning of the superconducting phase difference, which means that they can always be manipulated so that they intersect the surface. This makes them directly observable by scanning tunelling microscopy, which has already been used to detect normal vortices in proximized metals Stolyarov et al. (2018). If this ability of the superconducting phase difference to expand vortex loops to arbitrary sizes is a general feature of SNS-junctions, it could be used to detect vortex loops in systems where controlling the magnetic field strength is not an option, such as in system with a magnetic dipole inclusion.

Acknowledgements.

We thank M. Amundsen and V. Risinggård for helpful discussions. This work was supported by the Research Council of Norway through grant 240806, and its Centres of Excellence funding scheme grant 262633 “QuSpin”.