Superconducting Vortices in Half-Metals

TL;DR

This paper reveals that quantum vortices can form in half-metals with odd-frequency superconducting correlations, accompanied by spin supercurrents, and shows how interface disorder affects vortex nucleation, providing experimental insights.

Contribution

It demonstrates the existence of quantum vortices in half-metals with odd-frequency pairing and analyzes the effects of interface disorder on vortex formation.

Findings

Vortices can emerge in half-metals without conventional superconductivity.

A circulating spin supercurrent accompanies these vortices.

Interface magnetic disorder influences vortex nucleation positions.

Abstract

When the impurity mean free path is short, only spin-polarized Cooper pairs which are non-locally and antisymmetrically correlated in time may exist in a half-metallic ferromagnet. As a consequence, the half-metal acts as an odd-frequency superconducting condensate. We demonstrate both analytically and numerically that quantum vortices can emerge in half-metals despite the complete absence of conventional superconducting correlations. Because these metals are conducting in only one spin band, we show that a circulating spin supercurrent accompanies these vortices. Moreover, we demonstrate that magnetic disorder at the interfaces with the superconductor influences the position at which the vortices nucleate. This insight can be used to help determine the effective interfacial misalignment angles for the magnetization in hybrid structures, since the vortex position is experimentally observable via STM-measurements. We also give a brief discussion regarding which superconducting order parameter to use for odd-frequency triplet Cooper pairs in the quasiclassical theory.

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)$$

$$G_i{L}_i{\mathbf{e}}_n\cdot ({\hat{g}}_i^{\text{textsc}r}\bar{\nabla }{\hat{g}}_i^{\text{textsc}r})=\hat{I}({\hat{g}}_i^{\text{textsc}r},{\hat{g}}_j^{\text{textsc}r}),$$

$$\hat{I}$$

$$+\frac{G_{\varphi 2}^i}{8}\left[{\hat{g}}_i,{\hat{m}}_i{\hat{g}}_i{\hat{m}}_i\right]+\frac{i{G}_{\chi }^i}{8}\left[{\hat{g}}_i,F\left({\hat{g}}_j\right){\hat{g}}_i{\hat{m}}_i+{\hat{m}}_i{\hat{g}}_{i}F\left({\hat{g}}_j\right)\right]$$

$$+\frac{i{G}_{\chi }^j}{8}\left[{\hat{g}}_i,F\left({\hat{g}}_j{\hat{m}}_j{\hat{g}}_j-{\hat{m}}_j\right)\right],$$

$$\begin{array}{rl}{G}_0^i& =G_q\sum _{n=1}^N{T}_n^i,G_{\varphi }^i=2G_q\sum _{n=1}^N{\varphi }_n^i,\\ G_2^i& =G_q\sum _{n=1}^N{\left(T_n^i\right)}^2,G_{\chi }^i=G_q\sum _{n=1}^N{T}_n^i{\varphi }_n^i,\\ G_{\varphi 2}^i& =2G_q\sum _{n=1}^N{\left({\varphi }_n^i\right)}^2,\end{array}$$

$${\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),$$

$${\gamma }_{\text{textsc}hm}=\left(\begin{array}{cc}a& 0\\ 0& 0\end{array}\right).$$

$${\nabla }^{2}a+\frac{2\tilde{a}\nabla a\cdot \nabla a}{1-a\tilde{a}}=\frac{4(1+a\tilde{a})e\mathbf{A}\cdot (ae\mathbf{A}+i\nabla a)}{1-a\tilde{a}}-2i\epsilon a,$$

$$\begin{array}{r}{G}_{\text{textsc}hm}{\mathbf{e}}_n\cdot \nabla a=4G_0^{\text{textsc}hm}BCa-G_2^{\text{textsc}hm}{B}^2{C}^{2}a(a\tilde{a}+3)\\ +2i{G}_{\chi }^{\text{textsc}sc}B{C}^2\sin \theta \left(b{\mathrm{e}}^{-i\alpha }-\tilde{b}{\mathrm{e}}^{i\alpha }{a}^2\right)+2i{G}_{\text{textsc}hm}a{\mathbf{e}}_n\cdot \mathbf{A}e,\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}):=N_{\uparrow }(\epsilon ,\mathbf{r})=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 }_1(\mathbf{r},t):=\text{expectationvalue}{\psi }_{\uparrow }(\mathbf{r},t){\psi }_{\uparrow }(\mathbf{r},0) \\ =\frac{-i{N}_0}{2}{\int }_{-\infty }^{\infty }{f}_{\uparrow }(\mathbf{r},\epsilon )\tanh (\epsilon \beta /2)\sin (\epsilon t)\text{differential}\epsilon . \end{gathered}$$

$${\Psi }_2(\mathbf{r}):={\frac{\partial {\Psi }_1}{\partial t}|}_{t=0}=\frac{-i{N}_0}{2}{\int }_{-\infty }^{\infty }\epsilon f_{\uparrow }(\mathbf{r},\epsilon )\tanh (\epsilon \beta /2)\text{differential}\epsilon .$$

$${\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),$$

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

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

$$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/\lambda }+{\mathrm{e}}^{\left(\text{absolutevalue}x-2L_{\text{textsc}sc}\right)/\lambda }\right],$$

Full text

I Introduction

New physical phenomena can emerge at the interface between materials with different quantum order. One such example is in systems combining ferromagnetism and superconductivity, where it is possible to generate Cooper pairs that are both spin-polarized and correlated non-locally in time. This has become the basis for the field of superconducting spintronics Linder and Robinson (2015), which has as one of its goals to enable new types of devices utilizing spin-polarized supercurrents Eschrig (2015). On a more fundamental level, it is of interest to consider the interplay between different types of spontaneous symmetry breaking in such hybrid structures, since symmetry breaking governs a wide range of physical phenomena, including mass differences of elementary particles and phase transitions.

Half-metallic ferromagnets are $100\text{\,}\mathrm{\char 37\relax}$ spin-polarized, meaning that only one spin-band is conducting. Any supercurrent flowing through such a material, as has been observed experimentally Keizer et al. (2006), is therefore necessarily spin-polarized. This makes them especially interesting to study in order to understand how superconductivity adapts to a fully polarized environment. Much experimental and theoretical work has recently been conducted in order to understand hybrid structures involving superconductors (S) and half-metals (H) Keizer et al. (2006); Kalcheim et al. (2012); Singh et al. (2015); Peña et al. (2004); Eschrig et al. (2007, 2003, 2015); Ouassou et al. (2017); Dahal et al. (2009); Wu and Halterman (2018); Mironov and Buzdin (2015); Anwar et al. (2010); Asano et al. (2007).

One hitherto unsolved problem is whether superconducting quantum vortices can form in half-metallic materials. This is an unusual physical situation since the electrons are correlated exclusively non-locally in time, such that the half-metal in fact mimicks a purely odd-frequency Linder and Balatsky (2017) superconducting state. 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). In addition to being interesting from a fundamental physics point of view, understanding the behaviour of vortices is useful on a practical level. Their motion is a source of non-zero electrical resistance Halperin et al. (2010), and recently it has been proposed that vortices can be used as a means for long-range spin transport Kim et al. (2018).

It is known that vortices can form also inside normal metals that are in the proximity to a superconductor Cuevas and Bergeret (2007); Bergeret and Cuevas (2008); Amundsen et al. (2018); Stolyarov et al. (2018). Cooper pairs can then leak into the normal metal through the process of Andreev reflection Pannetier and Courtois (2000). This is the key mechanism behind the proximity effect which consists of weak superconductivity observed in a material placed in contact with a superconductor.

The proximity effect in half-metals is more complicated because it requires a mechanism which converts the spinless (singlet) Cooper pairs to spin-polarized (triplet) pairs. The theorized mechanism to produce such correlations involve spin mixing and spin-flip scattering at the interface Eschrig et al. (2003). Spin mixing introduces triplet correlations at the superconducting side, and spin-flip scattering mediates these correlations to the half-metallic side.

What allows us to investigate SH-heterostructures in the presence of an external magnetic field is the recent derivation of general spin-active boundary conditions for the quasiclassical theory applied to diffusive systems Eschrig et al. (2015); Ouassou et al. (2017). This means that we can apply the quasiclassical Usadel theory in such a way that the Cooper pair conversion mechanism described above is included.

Here, we apply this theory both analytically and numerically to a two-dimensional SHS-junction as depicted in fig. 1 under a constant perpendicular magnetic field. The constant magnetic field will in general have contributions both from the uniform magnetization in the half-metal and from an uniform applied field. We find that vortices indeed form in the purely triplet odd-frequency superconducting condensate existing in the half-metallic ferromagnet. Their location depends not only on the superconducting phase difference, but also on the effective interfacial magnetization directions characterizing either magnetic disorder or artificially inserted thin ferromagnetic layers Singh et al. (2015).

II Methodology

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

II.1 Quasiclassical theory

The SHS junction depicted in fig. 1 can be treated in the quasiclassical formalism under the assumption that the Fermi wavelength is much shorter than all other relevant length scales. This assumption is seemingly broken in the half-metal, where the exchange field is so large that the associated energy is not negligible compared to the Fermi energy Schlottmann (2003). The redemption comes from the realization that the spin-splitting in such system is so severe that there is effectively no interaction between different spin bands. We can then continue to use the quasiclassical theory if we, instead of including an exchange field, treat the spin-bands as non-interacting. If in addition the system is diffusive, meaning that the scattering time is small, the isotropic part of the quasiclassical Green’s function dominates and solves the Usadel equation Belzig et al. (1999); Chandrasekhar (2004); Rammer (2004); Usadel (1970),

[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. In the superconductors $D$ is a scalar while in the haf-metal it is $\operatorname{diag}(D,0)$, assuming that the conducting band is spin-up. The covariant derivative is $\bar{\nabla}\check{g}=\nabla\check{g}-ie\left[{\hat{\rho}_{3}{\bm{A}}}\,,\,{\check{g}}\right]$, where $\bm{A}$ is the vector potential, and

[TABLE]

is the quasiclassical impurity-averaged Green’s function. $\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 Usadel equation can be made dimensionless by introducing the Thouless energy, $\varepsilon_{\textsc{t}}\coloneqq D/L^{2}$, and measuring length scales relative to $L$ and energies relative to $\varepsilon_{\textsc{t}}$.

In general, the Usadel equation has to be solved together with the Maxwell equation in a self-consistent manner. However, we are interested here in the case where the width $W$ is smaller than the Josephson penetration depth. In this case one can ignore the screening of the magnetic field by the Josephson currents and the magnetic field is equal to the external one Barone and Paternò (1982).

To simplify the numerical and analytical calculations we assume that the magnetic field is zero inside the superconductors. This assumption is widely used Cuevas and Bergeret (2007); Alidoust and Halterman (2015); Bergeret and Cuevas (2008) and has been shown to give good agreement with experimental results Chiodi et al. (2012). Belzig et al. (1998) found that including the vector potential in the superconductors leads to corrections proportional to $\lambda_{s}/L$, where $\lambda_{s}$ is the effective penetration depth in the superconductors and $L$ is the length of the junction. Assuming $\lambda_{s}/L$ is small, we neglect the corrections from the vector potential to the Green’s function in the superconductors. We use the vector potential $e\bm{A}=-n\pi y/W\left[\theta(x)-\theta(x-L)\right]\bm{e}_{x}$, where $n=\Phi/\Phi_{0}$ is the number of flux quanta penetrating the half-metal and $\bm{e}_{x}$ is the unit vector in the $x$-direction.

II.2 Boundary Conditions

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 half-metal and superconductors separately, and the solutions must be connected through boundary conditions. These can be written

[TABLE]

where $\bm{e}_{n}$ is the outward-pointing normal vector for region $i$, $G_{i}$ is the bulk conductance of material $i$ and $L_{i}$ is the length of material $i$ in the direction of $\bm{e}_{n}$. $\hat{I}(\hat{g}^{\textsc{R}}_{i},\hat{g}^{\textsc{R}}_{j})$ is the matrix current from material $i$ to material $j$.

The matrix current for general spin-active boundaries between diffusive materials was found in 2015 by Eschrig et al. (2015). The expression to second order in spin-mixing angles and transmission probabilities was simplified in 2017 by Ouassou et al. (2017), making them easier to implement and more efficient to compute. This is the expression we will use here, and it reads

[TABLE]

where, for a half-metallic ferromagnet $F(\hat{v})=\hat{v}+\left\{{\hat{v}}\,,\,{\hat{m}}\right\}+\hat{m}\hat{v}\hat{m}$ and $\hat{m}_{k}=\operatorname{diag}(\bm{m}_{k}\cdot\bm{\sigma},\bm{m}_{k}\cdot{\bm{\sigma}}^{*})$. Here $\bm{\sigma}=(\sigma_{i},\sigma_{2},\sigma_{3})^{\intercal}$ is the vector of Pauli matrices and $\bm{m}_{k}$ is a unit vector in the direction of the magnetization experienced by a particle being reflected in material $k$. Similarly, $\hat{m}=\operatorname{diag}(\bm{m}\cdot\bm{\sigma},\bm{m}\cdot{\bm{\sigma}}^{*})$ where $\bm{m}$ is the unit vector in the direction of the magnetization being felt by a particle which is transmitted. The interface conductances are Ouassou et al. (2017)

[TABLE]

where $T^{i}_{n}$ and $\phi^{i}_{n}$ are respectively the transmission probability and spin mixing angle for tunneling channel $n$ from material $i$ to material $j$. For boundaries interfacing vacuum at $y=\pm W/2$, the matrix current is $\hat{I}=0$.

II.3 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)$.

There is only one conducting spin band in a half-metal, and as a result $\gamma$ has only one nonzero element,

[TABLE]

Substituting this into eq. 1 we get that $a$ solves the equation

[TABLE]

In section III.1 we will show that the Green’s function in the superconductors can be taken to be taken as the bulk Green’s function. Thus, the Ricatti parameter can be written as $\gamma_{\textsc{s}c}=\operatorname{antidiag}(b,-b)$, where $b$ is a function of $\varepsilon$ and the superconducting gap parameter $\Delta$. Inserting this and sections II.2 and 7 into eq. 3 we get

[TABLE]

where $B=b\tilde{b}-1$, $C=1/\left(1+b\tilde{b}\right)$ and $\theta$ and $\alpha$ are the angles for the magnetization directions on the superconducting side as shown in fig. 1. The corresponding equations for $\tilde{a}$ and $\bm{e}_{n}\cdot\nabla\tilde{a}$ is found by tilde conjugating eqs. 8 and 9.

From eq. 9 it can be seen that in order to have a non-zero solution for $a$ it is necessary that either $\sin\theta_{l}\neq 0$ or $\sin\theta_{r}\neq 0$. This means that having the effective magnetization angles at the interface not parallel with the uniform magnetization of the half-metal is necessary for the occurence of superconducting triplet correlations in the half-metal. This is because the creation of long-range spin-triplet require spin mixing and spin-flip scattering/spin rotation Eschrig et al. (2003); Linder and Robinson (2015). When spin-singlet Cooper pairs in a superconductors encounter an interface with effective magnetization direction $\bm{m}_{l}$, the spin-active boundary will produce spin-triplets with zero spin along $\bm{m}_{l}$. If $\bm{m}_{l}$ is not parallel to the magnetization of the half-metal, $\bm{m}$, then the triplet with zero spin along $\bm{m}_{l}$ will have a non-zero projection onto the spin-triplet state with spin one along $\bm{m}$. Hence, it the interfacial magnetization angles are not parallel with $\bm{m}$, the spin-active boundaries will produce equal spin triplets.

II.4 Observables

As mentioned initially, a vortex is accompanied by a normal-state density of states and a circulating supercurrent. This can be extracted from the quasiclassical Green’s function. In the following it will be useful to write

[TABLE]

In the half-metal, $f$ has only one nonzero component, $f_{\uparrow}$.

II.4.1 Local Density of States

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 half-metal we can write eq. 11 in terms of $a$,

[TABLE]

II.4.2 Supercurrent

Written in terms of the quasiclassical Green’s function, the current density is Belzig et al. (1999)

[TABLE]

Inserting eq. 10, 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. 13 can be rewritten

[TABLE]

The spin current can be found by multiplying the matrix in the integrand of eq. 13 by the Pauli matrix corresponding to the appropriate spin direction before taking the trace. For a half-metal magnetized in the $z$-direction, the $z$-component of the spin supercurrent polarization is proportional to the electric current while the remaining spin current components vanish.

II.4.3 Cooper Pair Correlation Function

The study of vortices in diffusive half-metals naturally raises the question of how to define the superconducting order parameter. In a normal superconductor, the order parameter is $\expectationvalue{\psi_{\uparrow}(\bm{r},0)\psi_{\downarrow}(\bm{r},0)}$ where $\psi_{\sigma}(\bm{r},t)$ is the field operator which destroys an electron with spin $\sigma$ at position $\bm{r}$ and time $t$. The same order parameter is used in a normal metal, but the analogous quantity for the half-metal, $\expectationvalue{\psi_{\uparrow}(\bm{r},0)\psi_{\uparrow}(\bm{r},0)}$ is always zero because of the Pauli exclusion principle. That is, the Cooper pair correlation function in a diffusive half-metal must vanish at equal times and is thus temporally non-local Linder and Balatsky (2017).

One approach of defining an order parameter in odd-frequency superconducting condensates, which is often used in the Bogolioubov-de Gennes formalism Halterman et al. (2007), is to keep the relative time coordinate $t$ finite between the field operators, that is

[TABLE]

Another frequently used strategy Abrahams et al. (1995) is to make the order parameter even in time by differentiation. This yields

[TABLE]

Below, we shall compare these two possible choices for order parameter describing the odd-frequency superconducting condensate to see which of them that correctly captures the vortex behavior.

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 quadratic quadrilateral 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 is present also in the full proximity effect regime. Dimensionless quantities are used in the analytics with length 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

In order to justify eq. 9 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.

First, let

[TABLE]

be the solution of the dimensionfull Usadel equation,

[TABLE]

in the superconductor at $x<0$. This gives an equation for $\delta\hat{g}$,

[TABLE]

where we have used that $\hat{g}_{\textsc{bcs}}$ solves the eq. 19 for a bulk superconductor and assumed that the variations of the gap parameter from the bulk value is negligible. 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}^{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 $\real(\lambda)>0$.

Using the boundary condition

[TABLE]

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

[TABLE]

where $C$ is some a function of $y$ and $\varepsilon$ to be determined by the final boundary condition. If the matrix current across this boundary is $\hat{I}_{\textsc{sc}}$, then

[TABLE]

From eq. 28 we see that $\real(\lambda)$ can be interpreted as the penetration depth of $\delta g$. Note that $\real(\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$. From the definition of $\lambda$ in eq. 26 we see that $\real(\lambda)=\lambda$ when $\varepsilon<\absolutevalue{\Delta}$ and $\real(\lambda)=\absolutevalue{\lambda}/\sqrt{2}$ otherwise.

$G_{\textsc{sc}}$ is the conductance over the whole length $L_{\textsc{sc}}$ of the superconductor and is therefore proportional to $1/L_{\textsc{sc}}$. Hence, $G_{\textsc{sc}}L_{\textsc{sc}}/\real(\lambda)$ is the normal-state conductance of a superconductor of length $\real(\lambda)$. From the definition of $\hat{I}$ in section II.2 we see 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[\real(\lambda)]$. A similar calculation shows that we can use $\hat{g}_{\textsc{bcs}}$ also in the superconductor at $x>L$.

Taking the superconducting coherence length $\xi$ as a measure of the inverse proximity effect penetration depth $\real(\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.

Next we turn to the solution of the Usadel equation in the half-metal, eq. 8, together with the boundary condition (9). In order to solve these equations we must make some simplifying assumptions. If we assume the proximity to be weak, we can keep only terms which are linear in $a$ and $\tilde{a}$ and their gradients. In this case the dimensionless Usadel equation (8) decouples:

[TABLE]

and so does the boundary conditions,

[TABLE]

Equation 31 can be further simplified in the so-called wide junction limit, where $n/W\ll 1$. If $\bm{A}=0$, the solution of eq. 31 is constant in the $y$-direction. Assuming this is approximately true also for small $\bm{A}$, we neglect the term $\partial_{y}^{2}a$. Equation 31 can now be solved exactly, as it is a second order ordinary differential equation with constant coefficients. The solution of eq. 31a is

[TABLE]

where $u=-2\pi iny/W$, $k=\sqrt{-2i\varepsilon}$ and $C_{1}$ and $C_{2}$ are independent of $x$.

Determining $C_{1}$ and $C_{2}$ requires the boundary conditions, which can be written

[TABLE]

where

[TABLE]

After some algebra, we find that the solution can be written

[TABLE]

where

[TABLE]

Note that the wide junction approximation is not applicable at small energies. This is because the solution will be slowly varying in the $x$-direction and therefore $\partial_{y}^{2}a$ is no longer negligible compared to $\partial_{x}^{2}a$.

When $\sin\theta_{l}=\sin\theta_{r}$, $a$ vanishes at $x=1/2$ and

[TABLE]

where $N$ is any integer. This means that $f_{\uparrow}$ and hence also the order parameters $\Psi_{1}$ and $\Psi_{2}$ vanish at these points. From eq. 12 we see that the density of states is equal to the normal state density of states at these points, indicating that these are indeed vortices. By Taylor expanding $a$ to first order around a root located at $(1/2,\tilde{y})$ we find

[TABLE]

where $x-1/2\sim\cos\theta$ and $y-\tilde{y}\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. These approximately $n$ roots are the only ones for $\Psi_{2}$, but for $\Psi_{1}$ there are relative times $t$ for which additional roots exist. Since each vortex is associated with a quantum of magnetic flux, $\Phi_{0}$, there should be at most $n$ vortices when the flux is $n\Phi_{0}$. This suggests that $\Psi_{1}$ is less suited for finding vortices than $\Psi_{2}$ if we identify vortices by the roots of the order parameter. Using $\Psi_{2}$ suggests that when $\sin\theta_{l}=\sin\theta_{r}$ and the magnetic flux is $n\Phi_{0}$, there will in the wide-junction limit be $n$ vortices whose location is determined by the difference in the superconducting phases and the magnetization angles.

The situation is more complicated when $\sin\theta_{l}\neq\sin\theta_{r}$. In this case the roots of eq. 36 depend on $\varepsilon$, and we will leave the discussion for how this affects the order parameter to the numerical investigation. However, some insight can still be had from the analysis. Scaling $\sin\theta$ in the boundary condition (32) is equivalent to scaling the conductance $G_{\chi}^{\textsc{sc}}$. That is, if $\sin\theta_{r}<\sin\theta_{l}$, the proximity effect should be weaker at the right side, meaning that the vortices should be pushed to the right. This is indeed what we find numerically.

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}}$, $G_{\textsc{hm}}=3G_{0}^{\textsc{hm}}$, $G_{\chi}^{\textsc{sc}}=0.01G_{0}^{\textsc{hm}}$, $G_{2}^{\textsc{hm}}=0.002G_{0}^{\textsc{hm}}$ and $\phi_{l}=\alpha_{l}=0$ common for all the numerical calculations. We obtain qualitatively similar results for other choices of the conductance parameters $G_{i}^{\textsc{hm/sc}}$. 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).

III.2.1 Local Density of States

In the symmetric case ($\sin\theta_{r}=1.0$), we suspect from the analysis above that for all energies, the local density of states is equal to that of the normal state in $n$ points along the line $x=0.5$, where $n$ is the number of magnetic flux quanta. This is also what we find numerically, as shown in fig. 2. Figure 2 shows the local density of states at various energies for the symmetric case with $n=2$ with red lines close to where it is equal to the normal state density of states. We see that there are indeed two locations where the difference between the local density of states and the normal state density of states vanish for all energies, and that the locations of these points are exactly those predicted by the analysis.

The situation is slightly different for the asymmetric case ($\sin\theta_{r}\neq 1$), as can be seen in fig. 3. There is no longer a single point where $N=N_{0}$ for all energies. Instead there are points where $N$ is equal to or almost equal to $N_{0}$ for all energies, as can be seen in fig. 3. In sections III.2.2 and III.2.3 we will se that these points also are associated with vortices. Qualitatively, the points where $N$ stays close to $N_{0}$ are different in the asymmetric case. In the symmetric case the points at $x=1/2$ are mostly isolated, but in the asymmetric case the point is part of a line where $N=N_{0}$ which stretches towards the side there $\sin\theta$ is smaller. For $n=2$, as can be seen in fig. 3, the location where $N$ is close to or equal to $N_{0}$ for alle nergies occurs at $y=\pm W/4$, which is also the $y$-values where the vortices are in the symmetric case.

Note that the energy-dependence on the position where $N=N_{0}$ is a not unique for the special case of a SHS-junction. The same phenomenon occurs in normal SNS-junctions if the conductances at the interfaces are unequal. As mentioned above, changing $\sin\theta$ is the same as changing the conductance $G_{\chi}^{\textsc{sc}}$.

III.2.2 Supercurrent

Figure 4 shows the current amplitude and direction for the same two cases as was discussed in section III.2.1. In both cases there are exactly two points where the supercurrents vanish and where the supercurrent circles around. This indicates the existence of two superconducting vortices, which is in accordance with the analysis, local density of states and the fact that the system is experiencing two quanta of magnetic flux. We underline that the supercurrents accompanying the vortices induced in the half-metal are fully spin-polarized and carried by triplet Cooper pairs. This is different from the non-polarized charge supercurrents circulating vortices in previously studied hybrid structures Cuevas and Bergeret (2007); Bergeret and Cuevas (2008); Amundsen et al. (2018).

The $y$-value of the points with circulating supercurrents are the same as what is expected by the analysis, given by eq. 38. In the symmetric case these points are midway between the superconductors, i.e.  at $x=L/2$, while in the asymmetric case they are moved slightly toward the side where $\sin\theta$ is smaller. Hence, the vortex locations as given by the supercurrents agrees with the analysis as well as the results from the local density of states.

One feature of the asymmetric case worth noting is that the supercurrent is suppressed to the right of the vortex, which is towards the side where $\sin\theta$ is smaller. This is in agreement with the fact that this region had a local density of states which was closer to the normal state value, as shown in section III.2.1.

III.2.3 Cooper Pair Correlation Function

Investigating vortices in odd-frequency superconductors gives rise to the problem of choosing what order parameter to use. The Cooper pair correlation function which is normally used is not applicable as it is identically equal to zero. The two most obvious remaining choices are

[TABLE]

and

[TABLE]

as mentioned above. From looking at the local density of states and the supercurrent we already know that we expect vortices and that their location should have $y$-value given by eq. 38 and be at $x=L/2$ in the symmetric case and closer to the side where $\sin\theta$ is smaller in the asymmetric case. Comparing the position of the roots of $\Psi_{1}$ and $\Psi_{2}$ with the position of the vortices as predicted by the local density of states and the supercurrent can be used to give an identification of how good the order parameters are as means to find vortices.

Figures 5 and 6 shows $\Psi_{1}$ with various values of $t$ for the two same cases which was shown in sections III.2.1 and III.2.2. Figure 5 shows the symmetric case with $n=2$, $\phi_{r}=\alpha_{r}=0$ and $\sin\theta_{r}=1$, while fig. 6 shows the asymmetric case with $n=2$, $\phi_{r}=\alpha_{r}=0$ and $\sin\theta_{r}=0.5$. $\Psi_{2}(\bm{r})$ looks identical to $\Psi_{1}(\bm{r},0.01/\varepsilon_{\textsc{t}})$, which is shown in the figures, but is a factor $100$ larger in magnitude. $\Psi_{2}$ is therefore not shown. In both cases, there are exactly two roots in $\Psi_{2}$ and $\Psi_{1}$ for small $t$, and the position are the same as those given for the vortices by the local density of states and supercurrent. Around these two roots are a phase winding of $2\pi$.

For larger values of $t$, $\Psi_{1}(\bm{r},t)$ seems to be less suited for finding vortices. In the symmetric case there are additional roots which appear. These additional roots also has a corresponding phase winding of $2\pi$, but does not correspond to vortices when compared to the local density of states or supercurrent. The situation is even worse for large $t$ in the asymmetric case. In addition to having extra roots, the original roots which correspond to vortices are either moved or not present. This is also the case for other choices of $n$, $\phi_{r}$, $\alpha_{r}$ and $\sin\theta_{r}$.

Thus, we conclude that using $\Psi_{2}$ seems best suited as order parameter for the numerical investigation of quantum vortices in an purely odd-frequency superconducting condensate. Alternatively, one could use $\Psi_{1}$ with $\varepsilon_{\textsc{t}}t\ll 1$ which will give the same result since $\Psi_{1}\sim\Psi_{2}t$ as $t\to 0$.

Figures 7 and 8 shows the phase of $\Psi_{2}$, in addition to the supercurrent and amplitude of $\Psi_{2}$, for a magnetic flux of $4\Phi_{0}$ and with $\sin\theta_{r}=1$, $\phi_{r}=\alpha_{r}=0$ and $\sin\theta_{r}=0.5$, $\phi_{r}-\alpha_{r}=\pi/2$, respectively. Note that there is indeed a phase winding of $2\pi$ around the vortices, as mentioned above.

III.2.4 Vortex Position

So far we have looked only at the case with $\phi_{r}=\alpha_{r}=0$. Choosing a nonzero value for $\phi_{r}-\alpha_{r}$ moves the vortices along the $y$-axis, and the locations correspond to those predicted by the analysis and given in eq. 38. Figures 7 and 8 shows the amplitude and phase of $\Psi_{2}$ as well as the supercurrent $\bm{j}$ for a wide junction subjected to a magnetic flux of $4\Phi_{0}$, with $\sin\theta_{r}=1$, $\phi_{r}-\alpha_{r}=0$ and $\sin\theta_{r}=0.5$, $\phi_{r}-\alpha_{r}=\pi/2$, respectively. In fig. 8 the vortices are shifted $W/16$ down along the $y$-axis, which is exactly what is predicted by eq. 38. Note also that the vortices are moved towards the side where $\sin\theta$ is smaller in fig. 8.

The dependence of the vortex position on both $\alpha$ and $\theta$ suggests an experimental method to determine the effective magnetization angles describing disorder in the form of interfacial misaligned moments or artificially inserted misaligned magnetic layers in half-metallic hybrid structures. For a fixed value of the magnetic flux and phase difference $\phi_{r}-\phi_{l}$ (which is tunable by the applied current), the $y$-coordinates of the vortices gives information about the azimuthal angles $\alpha_{l}$ and $\alpha_{r}$, while the $x$-coordinates gives information about the polar angles $\theta_{l}$ and $\theta_{r}$. This approach could possibly be easier than trying to measure the magnetization angles directly, especially if the non-collinear magnetization angle at the interface is produced by the natural misalignment of magnetic dipoles arising from the conjunction of different atomic structures at the interface.

IV Conclusion

Using the quasiclassical Usadel theory we have found both analytically and numerically that superconducting vortices occur also in purely odd-frequency superconducting condensates that exist in proximized half-metallic ferromagnets. Because half-metals only have one conducting spin-band, the vortcex cores are circulated by fully polarized spin supercurrents. An additional feature of vortex structures in half-metals compared to a normal metal is that the vortex position depend on the magnetization angles. We suggest that this insight can be used to help determine these angles in hybrid structures.

The study of vortices in odd-frequency superconducting condensates naturally raises the question of how to define the superconducting order parameter. Comparing the roots of the superconducting order parameter to the location vorticies we have found that the order parameter that works best is that which is made even in time by differentiation, $\frac{\partial}{\partial t}\expectationvalue{\psi_{\uparrow}(\vec{r},t)\psi_{\uparrow}(\vec{r},0)}\rvert_{t=0}$.

Acknowledgements.

We thank M. Amundsen 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”.