Anomalous currents and spontaneous vortices in spin-orbit coupled superconductors

TL;DR

This paper proposes a new mechanism for generating supercurrents and vortices in superconductors with magnetic inclusions using spin-orbit coupling.

Contribution

The novelty lies in using spin-orbit coupling and magnetic inclusions to generate spontaneous vortices without an external magnetic field.

Findings

Spin-orbit coupling creates an effective gauge field that influences Cooper pairs near magnetic inclusions.

The mechanism can either enhance or reduce magnetization depending on material parameters.

Spontaneous vortices can form in superconductors with the right inclusion distribution.

Abstract

We propose a mechanism which can generate supercurrents in spin-orbit coupled superconductors with charged magnetic inclusions. The basic idea is that through spin-orbit interaction, the in-plane electric field near the edge of each inclusion appears to the electrons as an effective spin-dependent gauge field; if Cooper pairs can be partially spin polarized, then each pair experiences a nonzero net transverse pseudo-gauge field. We explore the phenomenology of our mechanism within a Ginzburg-Landau theory, with parameters determined from a microscopic model. Depending on parameters, our mechanism can either enhance or reduce the total magnetization upon superconducting condensation. Given an appropriate distribution of inclusions, we show how our mechanism can generate superconducting vortices without any applied orbital magnetic field. Our mechanism can produce similar qualitative behavior to the “magnetic memory effect” observed in 4Hb-TaS21. However, the magnitude of the effect in that material seems larger than our model can naturally explain.

Figures

Charged ferromagnetic inclusions in a layered superconductor.a A single inclusion. A vertical electric field (not shown) is assumed to be present, providing Rashba coupling; it may be uniform, or random due to charge disorder above and below the plane. The inclusion exerts an out-of-plane exchange field $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{b}}={b}_{z}{\hat{\bf{z}}}$$\end{document}$ throughout its i

Analytical (curves) and numerical (symbols) results for the current density 〈*J*^*y*^〉 given the Hamiltonian of Eq. (14), as a function of (a) varying vertical electric field $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{z}$$\end{document}$ at fixed superconducting gap Δ, and (b) varying Δ at fixed $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \us

Full text

Introduction

The most dramatic manifestation of the superconducting state is its electromagnetic response^2–4^. As has been known since the early days of the subject, any superconductor exposed to a magnetic flux density will generate dissipationless supercurrents in an effort to screen that flux density. In type-II superconductors, if the applied flux density is sufficiently large, but not so large as to completely destroy the condensate, then flux penetrates through the cores of quantized vortices. These types of superflow, both related to applied flux, are generic features of a charged condensate subject to electromagnetic gauge invariance.

Depending on details, other forms of superflow are possible. In particular, spin-orbit coupling (SOC) intermingles orbital dynamics (such as current flow) with spin dynamics. A variety of superconducting magnetoelectric effects^5–21^ and related phenomena^22–25^ can result, putting new twists on the usual relationship between magnetism and supercurrent. In this paper, we identify a magnetoelectric contribution to the supercurrent near inhomogeneities in superconductors with SOC, which, to the best of our knowledge, has thus far gone unnoticed.

Our basic idea begins with the observation^26^ that spin-orbit coupling of the form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}_{{\rm{SO}}}=\lambda (\boldsymbol{\mathcal{E}} \times {\boldsymbol{\sigma}})\cdot {\bf{p}}$$\end{document} resembles the minimal coupling between the electron momentum operator p and a non-abelian gauge potential. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}$$\end{document} (r) is the local electric field at position r, and σ is the vector of Pauli matrices for spin. Considering a layered quasi-two-dimensional system and focusing on the term proportional to the in-plane electric field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}_{\|}$$\end{document} , we point out that the non-abelian structure reduces to a σ^z^ matrix, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({\boldsymbol{\mathcal{E}}}_{\parallel }\times {\boldsymbol{\sigma }})\cdot {\bf{p}}={\sigma }^{z}({\boldsymbol{\mathcal{E}}}_{\parallel }\times \hat{{\bf{z}}})\cdot {\bf{p}}$$\end{document} . This term looks precisely like the coupling between p and an ordinary U(1) electromagnetic vector potential, but with opposite effective vector potentials for opposite spin projections. If Cooper pairs can be at least partially spin polarized, then each pair experiences a nonzero net effective vector potential \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{{\bf{A}}}$$\end{document} . If \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\boldsymbol{\nabla}}\times {\tilde{\bf{A}}}\,\ne\, {\mathbf{0}}$$\end{document} , then the condensate will respond to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\boldsymbol{\nabla}}\times {\tilde{\bf{A}}}$$\end{document} just as it would to a real magnetic flux, either by screening à la Meissner, or by forming vortices. In this paper, we show how this mechanism can play out when Rashba SOC [proportional to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{z}$$\end{document} ] and exchange fields applied by frozen local moments [modeled as an effective spatially-varying Zeeman splitting b(r) ⋅ σ] combine to facilitate the required Cooper pair spin polarization. We note that Holst et al. studied related physics near the edge of a two-dimensional superconductor in ref. ^27^.

In ref. ^11^, Pershoguba et al. studied a related situation of a superconductor subjected to a local exchange field b(r), but with zero in-plane electric field ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}_{\|}$$\end{document} = 0). They identified two magnetoelectric supercurrent contributions resulting from Rashba SOC:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{j}}}_{{\rm{extra}}}=\alpha {\hat{\bf{z}}}\times {\bf{b}}+\beta\,{\boldsymbol{\nabla}} \times ({b}_{z}{\hat{\bf{z}}}).$$\end{document}

The α term is a well-known magnetoelectric term^5–12^, and the β term was pointed out by Pershoguba et al.^11^. In this work, we identify an additional contribution,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{j}}}_{\gamma }=\gamma {\boldsymbol{\mathcal{E}}}_{\parallel }\times ({b}_{z}\hat{{\bf{z}}}),$$\end{document}

which arises when the in-plane electric field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}_{\|}$$\end{document} does not vanish, corresponding to local breaking of inversion and rotation symmetries; in homogeneous systems, the γ term is forbidden by any point-group symmetry acting nontrivially in the plane (such as inversion or rotation)^28^. Like the α term, and unlike the β term, the γ term is forbidden in the normal state (see Methods section IV C for details). We note that the α term may be thought of as analogous to the γ term, with the electric field oriented out-of-plane and proportional to α. Unlike the α term, and like the β term, the γ term is even under the mirror Mz, so is not forbidden by that symmetry. We analyze the γ term in the context of Ginzburg-Landau theory in the Results section. We propose a scenario where the material contains inclusions which are both charged and magnetic, giving rise to the fields b and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}$$\end{document} mentioned above; see Fig. 1a. In the superconducting state, supercurrents associated with our γ term contribute to the total magnetization. We show that, in the presence of a finite density of such magnetic inclusions and assuming that their magnetizations are aligned, the superconductor may respond by producing vortices located in between the inclusions, as shown in Fig. 1b. Following the Ginzburg-Landau discussion, we provide a microscopic derivation of the coefficient γ in the Results section.Fig. 1. Charged ferromagnetic inclusions in a layered superconductor.a A single inclusion. A vertical electric field (not shown) is assumed to be present, providing Rashba coupling; it may be uniform, or random due to charge disorder above and below the plane. The inclusion exerts an out-of-plane exchange field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{b}}={b}_{z}{\hat{\bf{z}}}$$\end{document} throughout its interior, and a radial electric field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}_{\|}$$\end{document} (light blue arrows) in a region of width d around its edge; see Methods. A current density jγ (dark blue) flows in the edge region. b Given many inclusions of the type shown in panel (a), vortices can form in between the inclusions; see Results.

Our study was motivated by anomalous observations of 4Hb-TaS_2_^1,29–33^, especially those implying a hidden magnetic memory^1^. The memory manifests as an unexplained increase in total magnetization, in the form of spontaneous vortices, upon cooling into the superconducting state. Crucially, the observed vortex density is far in excess of any magnetic flux, internal or external, present during superconducting condensation. We discuss the observations in question in Results section II C. We find that our mechanism can indeed produce a significant change in total magnetization upon condensation, albeit with a magnitude likely too small to explain the experiments.

Results

Ginzburg-Landau theory

We now provide a Ginzburg-Landau description of our effect. We consider a layered 3D superconductor with spin-orbit coupling, containing charged magnetic inclusions which produce an inhomogenous exchange field b(r) and electric field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}$$\end{document} (r). We assume that the magnetic moments of the inclusions are polarized out-of-plane for simplicity, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{b}}={b}_{z}{\hat{\bf{z}}}$$\end{document} ; the α term then vanishes, and we ignore it from now on. Deep in the superconducting phase, up to constant terms (which do not contribute to the supercurrent), the free-energy density is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{F}}[\theta ,{\bf{A}};{\tilde{\bf{A}}}]=\frac{\kappa }{2}{\left({\boldsymbol{\nabla }}\theta +2e{\bf{A}}+2e{\tilde{\bf{A}}}\right)}^{2}+\frac{1}{2{\mu }_{0}}{\left({\boldsymbol{\nabla}}\times {\bf{A}}\right)}^{2}$$\end{document}

where − e < 0 is the electron charge, κ is the superfluid stiffness, θ is the phase of the order parameter, A is the electromagnetic vector potential, and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\rm{A}}}=\frac{-\gamma }{\kappa {(2e)}^{2}}{\boldsymbol{\mathcal{E}}}_{\parallel }\times ({b}_{z}{\hat{\bf{z}}})$$\end{document}

is the spin-orbit induced pseudo-gauge field experienced by a Cooper pair^34^. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}_{\|}$$\end{document} is the in-plane component of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}$$\end{document} . At this stage, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\bf{A}}}$$\end{document} is motivated by the physical intuition we described in the introduction; we will justify it from a microscopic model below. Note that we treat \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\bf{A}}}$$\end{document} as externally imposed, while A is determined by minimizing the free energy. The supercurrent density corresponding to the free energy in Eq. (3) is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{j}}=-2e\kappa \left({\boldsymbol{\nabla }}\theta +2e{\bf{A}}+2e{\tilde{\bf{A}}}\right).$$\end{document}

We now consider the geometry depicted in Fig. 1(a). We suppose that a charged magnetic inclusion constitutes a disk of radius R with some vertical thickness. The inclusion produces a net radial electric field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}_{\|}$$\end{document} in a region of width d at the edge of the inclusion (we justify this field configuration in the Methods section). The electric field corresponds to a potential difference V∥ between the inside and outside of the inclusion, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\boldsymbol{\mathcal{E}}}_{\parallel } \sim -({V}_{\parallel }/d){\hat{\bf{r}}}$$\end{document} within the edge region. Neglecting variation of b across the edge for simplicity, the resulting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\bf{A}}}$$\end{document} in the edge region is uniform and purely azimuthal: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\bf{A}}}=-\gamma {b}_{z}{V}_{\parallel }{\hat{\boldsymbol{\varphi }}}/(4{e}^{2}\kappa d)$$\end{document} . The total γ-term current flowing around the edge of the inclusion, per unit vertical thickness, is then

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}_{\gamma }=d{{\bf{j}}}_{\gamma }\cdot {\hat{\boldsymbol{\varphi }}}=\gamma {b}_{z}{V}_{\parallel }.$$\end{document}

The corresponding magnetic moment per unit vertical thickness is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${M}_{\gamma }=\pi {R}^{2}{I}_{\gamma }=\pi {R}^{2}\gamma {b}_{z}{V}_{\parallel }.$$\end{document}

Note that this current can flow in either direction around the inclusion, depending on the sign of V∥ (i.e., the charge of the inclusion), and that of γ; Mγ can thus either add to or cancel against the spin magnetic moment, which we consider next.

We compare Mγ to the paramagnetic response of the electronic spins in the magnetic inclusion (present in the normal state as well),

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${M}_{{\rm{n}}}=\pi {R}^{2}\chi {b}_{z}$$\end{document}

where χ is the spin susceptibility. To estimate the scale of Mn, we use the Pauli value of χ, yielding

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${M}_{{\rm{n}}}=3{\pi }^{2}{R}^{2}m{\mu }_{{\rm{B}}}{b}_{z}$$\end{document}

where m is the electron mass and μB is the Bohr magneton. There is also the current due to the β term in Eq. (1). This term is not required to vanish in the normal state, but it does so in the microscopic model we consider below. Furthermore, within our microscopic model, the current contribution from the β term in the superconductor is suppressed relative to the one from our γ term by a factor of the characteristic spin-orbit energy over the chemical potential; see Methods. Since this factor is likely to be small, we neglect the β term in our present estimate.

The extra contribution to the total magnetic moment in the superconducting state relative to that present already in the normal state is therefore

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{M}_{\gamma }}{{M}_{{\rm{n}}}}=\frac{\gamma {V}_{\parallel }}{3\pi m{\mu }_{{\rm{B}}}}.$$\end{document}

To estimate this ratio, we need an estimate for the coefficient γ, which depends on a choice of microscopic model. We derive γ for a simple model below, starting from Eq. (14). Before doing so, we address the question of spontaneous vortex formation within the Ginzburg-Landau picture.

Starting from the free-energy density of Eq. (3), we may assume that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\bf{A}}}={\tilde{\bf{A}}}_{T}$$\end{document} is purely transverse; we can always absorb its longitudinal part \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\bf{A}}}_{L}$$\end{document} with an appropriate gauge transformation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\boldsymbol{\nabla }}\theta \to {\boldsymbol{\nabla }}\theta -2e{\tilde{\bf{A}}}_{L}$$\end{document} . With that assumption, we may write \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\bf{A}}}$$\end{document} as the curl of some field. Assuming κ is uniform, we can write

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\boldsymbol{\nabla }}\times {\tilde{\bf{h}}}=-4{e}^{2}{\mu }_{0}\kappa {\tilde{\bf{A}}}={\mu }_{0}{{\bf{j}}}_{\gamma }.$$\end{document}

Note that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\bf{h}}}/{\mu }_{0}$$\end{document} is the magnetization corresponding to the γ-term current, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{j}}}_{\gamma }=-4{e}^{2}\kappa {\tilde{\bf{A}}}$$\end{document} . We may choose \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\bf{h}}}$$\end{document} to be nonzero within the magnetic inclusions and zero outside the inclusions. Neglecting boundary contributions and constant terms (see Methods), the total bulk free energy can be rewritten as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{ll}F[\theta ,{\bf{A}};{\tilde{\bf{A}}}]=\int\,{{\rm{d}}}^{2}{\bf{r}}{\mathcal{F}}\\=\int{{\rm{d}}}^{2}{\bf{r}}\left[\displaystyle\frac{\kappa }{2}{\left({\boldsymbol{\nabla }}\theta +2e{\bf{A}}\right)}^{2}+\displaystyle\frac{\pi }{e{\mu }_{0}}{{\bf{n}}}_{v}\cdot {\tilde{\bf{h}}}+\displaystyle\frac{1}{2{\mu }_{0}}{({\boldsymbol{\nabla }}\times {\bf{A}}-{\tilde{\bf{h}}})}^{2}\right],\end{array}$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{n}}}_{v}=-\frac{1}{2\pi }{\boldsymbol{\nabla }}\times {\boldsymbol{\nabla }}\theta$$\end{document}

is the density of superconducting vortices. The physical configuration is determined by minimizing the free energy over A and nv. Eq. (12) has the form of the standard free energy of a type-II superconductor with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\bf{h}}}$$\end{document} playing the role of the external magnetic field (neglecting variations of the amplitude of the superconducting order parameter), with an additional term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\propto {{\bf{n}}}_{v}\cdot {\tilde{\bf{h}}}$$\end{document} . Here, we have not included the contribution of vortex cores, where the order parameter amplitude is suppressed over scales of the coherence length, giving rise to a positive core energy (and hence a nonzero lower critical field for vortex formation).

It is useful to first consider the situation where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\bf{h}}}$$\end{document} is uniform in the bulk. By Eq. (11), this corresponds to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\bf{A}}}\propto {\boldsymbol{\nabla}}\times {\tilde{\bf{h}}}$$\end{document} (and its associated current) taking nonzero values only near the edge; this situation corresponds to bz being uniform in the bulk, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}_{\|}$$\end{document} ≠ 0 near the edge and zero in the bulk. For a state with no vortices (nv = 0), the free energy is exactly that of a superconductor in a genuine uniform applied field, which is then Meissner screened by the condensate. The bulk free energy of this state is zero. On the other hand, unlike in a usual superconductor subject to a magnetic field, the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{n}}}_{v}\cdot {\tilde{\bf{h}}}$$\end{document} term prevents the formation of vortices in the bulk. To see this, we assume for simplicity that the penetration depth is extremely large, such that B = ∇ × A is essentially uniform. In the usual case of a superconductor in an external field, the superconductor can lower its energy by developing a vortex lattice with spatially-averaged vortex density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\overline{{\bf{n}}}}_{v}=-\overline{{\boldsymbol{\nabla}}\times {\boldsymbol{\nabla }}\theta }/(2\pi )=2e{\boldsymbol{\nabla}}\times {\bf{A}}/(2\pi )=2e{\tilde{\bf{h}}}/(2\pi )$$\end{document} , where overlines denote spatial averages. However, in the present case, that configuration would yield a positive bulk free-energy density through the term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi {\overline{{\bf{n}}}}_{v}\cdot {\tilde{\bf{h}}}/(e{\mu }_{0})={\tilde{\bf{h}}}^{2}/{\mu }_{0}\, > \,0$$\end{document} .

In an inhomogeneous system, however, we now show that spontaneous vortices can form. Intriguingly, the vortices appear away from the magnetic inclusions responsible for their formation. We consider the situation depicted in Fig. 1(b), with many charged magnetic inclusions embedded into a superconductor, whose magnetizations are all aligned along \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+{\hat{\bf{z}}}$$\end{document} . We assume that the inclusions are distributed randomly in 3D. As before, we assume that the penetration depth is sufficiently large (in particular, larger than the distance between the inclusions) so that the magnetic field can be considered to be uniform. Since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\bf{h}}}$$\end{document} is only nonzero within the inclusions and the penetration depth is large, the system can lower its energy by developing vortices in between the inclusions. At the location of the vortices \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\bf{h}}}=0$$\end{document} , so the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{n}}}_{v}\cdot {\tilde{\bf{h}}}$$\end{document} term does not impose any energetic cost on such a vortex^35^. However, the large penetration depth means that the flux associated with each vortex penetrates through a large region, including some nearby inclusions: this arrangement hence saves energy by reducing the average value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${({\boldsymbol{\nabla}}\times {\bf{A}}-{\tilde{\bf{h}}})}^{2}$$\end{document} , while keeping the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{n}}}_{v}\cdot {\tilde{\bf{h}}}$$\end{document} term zero.

Recall now our motivating question: Can this mechanism explain the spontaneous formation of vortices in 4Hb-TaS_2_? (We address this question in a subsection below). This question reduces to the issue of magnetization enhancement we considered previously. Magnetization \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathfrak{m}}$$\end{document} acts as an applied field, and the superconducting condensate forms vortices in order to reduce a free energy term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\propto {({\boldsymbol{\nabla}}\times {\bf{A}}-{\mu }_{0}{\mathfrak{m}})}^{2}$$\end{document} . Each inclusion contributes a magnetic moment Mn + Mγ estimated above in Eqs. (7),(8), and (9), so the average magnetization density would be (Mn + Mγ)ρinc, where ρinc is the density of inclusions. Hence, the ratio of novel (i.e., γ term related) to conventional vortex densities is Mγ/Mn, given by Eq. (10). We now provide a microscopic calculation of the coefficient γ appearing in that equation.

Microscopic model

As discussed in the previous sections, the basic ingredients for the mechanism we propose are superconductivity, spin-orbit coupling, and an exchange field. We consider a simple 2D model containing all three:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{H}}={{\mathcal{H}}}_{0}+{{\mathcal{H}}}_{b}+{{\mathcal{H}}}_{{\rm{SO}}},$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}_{0}=\left(\frac{1}{2m}{{\bf{p}}}^{2}{\tau }^{z}-\mu {\tau }^{z}-\Delta {\tau }^{x}\right){\sigma }^{0}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}_{b}=-{b}_{z}{\sigma }^{z}\,{\tau }^{0}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}_{{\rm{SO}}}=-\lambda e({\boldsymbol{\sigma }}\times {\boldsymbol{\mathcal{E}}} )\cdot {\bf{p}}\,{\tau }^{z}.$$\end{document}

p = ∇/i is the momentum operator (we use ℏ = 1). We work in the Nambu basis \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Psi }^{\dagger }=\left({\psi }_{\uparrow }^{\dagger }{\psi }_{\downarrow }^{\dagger }{\psi }_{\downarrow }-{\psi }_{\uparrow }\right)$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\psi }_{\uparrow /\downarrow }^{\dagger }({\bf{r}})$$\end{document} creates an electron with the indicated spin at position r. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\boldsymbol{\mathcal{E}}} ={\boldsymbol{\mathcal{E}}}_{\parallel }+{{\mathcal{E}}}^{z}{\hat{\bf{z}}}$$\end{document} is the local electric field, and bz is the local vertical exchange field. Note that e \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}$$\end{document} has units of [energy] [length^−1^], and λ has units of [length^2^]. Pauli matrices σ^j^ act on spin, and τ^j^ act on particle/hole. We assume a gauge where Δ > 0 for convenience. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{{\rm{F}}}=\sqrt{2m\mu }$$\end{document} is the Fermi momentum at chemical potential μ absent spin-orbit coupling, with corresponding velocity vF = kF/m. We will treat \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}_{\|}$$\end{document} (r) perturbatively, so the largest characteristic energy scale of the spin-orbit coupling is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$| \lambda e{{\mathcal{E}}}^{z}| {k}_{{\rm{F}}}$$\end{document} .

Before describing our calculation, let us spend a few words on the electric field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}$$\end{document} . The out-of-plane component \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{z}$$\end{document} is only required in order to mix triplet components into the singlet Cooper pairs produced by Δ, so the effect is independent of the sign of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{z}$$\end{document} . Then, since the dependence of our effect on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{z}$$\end{document} is even, the effect could result from charge disorder corresponding to local fields with significantly nonzero average \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\left({{\mathcal{E}}}^{z}\right)}^{2}$$\end{document} , even if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{z}$$\end{document} averages to zero.

With respect to the in-plane electric field, in our mechanism, we assume that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}_{\|}$$\end{document} is produced by charged inclusions, as we discussed above in the Ginzburg-Landau context. Since our goal in the present section is to determine the magnitude of the current, but not its precise distribution, we temporarily disregard the inclusion geometry and instead assume uniform fields b^z^ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}$$\end{document} for simplicity^36^.

Now for the calculation. We provide here the outline and key results, leaving the detailed steps to the Methods section. We begin with the Hamiltonian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{H}}$$\end{document} given in Eq. (14). Exploiting translation invariance, we then express \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{H}}$$\end{document} in momentum space, i.e, p = ∇/i → k. The free-electron dispersion, relative to the chemical potential, is ξk = k^2^/(2m) − μ. We take \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}} =({{\mathcal{E}}}^{x},0,{{\mathcal{E}}}^{z})$$\end{document} without loss of generality—this just amounts to choosing coordinates so that the x axis points along \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}_{\|}$$\end{document} . The q → 0 current-density operator is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{J}}=\frac{(-e)}{2{L}^{2}}\sum _{{\bf{k}}}{\Psi }_{{\bf{k}}}^{\dagger }\left[\frac{{\bf{k}}}{m}{\sigma }^{0}-\lambda e{({\boldsymbol{\sigma }}\times {\boldsymbol{\mathcal{E}}} )}_{x,y}\right]{\tau }^{0}{\Psi }_{{\bf{k}}}\equiv \frac{1}{{L}^{2}}\sum _{{\bf{k}}}{\Psi }_{{\bf{k}}}^{\dagger }{{\bf{J}}}_{{\bf{k}}}{\Psi }_{{\bf{k}}}.$$\end{document}

The subscript x,y indicates that only the in-plane components of the cross product are considered. The capital J distinguishes this current operator from the Ginzburg-Landau current j.

We compute the current perturbatively in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{x}$$\end{document} and bz. We start by perturbing around bz = 0, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{x}$$\end{document} temporarily kept non-perturbative, and then expand the resulting expression in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{x}$$\end{document} . Writing k = (i**ωn, k) with ωn = (2n + 1)π**T a fermionic Matsubara frequency at temperature T, the linear current response to bz is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {\bf{J}}\rangle =\frac{{b}_{z}}{2}\frac{T}{{L}^{2}}\sum _{k}{\rm{tr}}\{{{\bf{J}}}_{{\bf{k}}}{{\mathcal{G}}}_{k}{\sigma }^{z}{\tau }^{0}{{\mathcal{G}}}_{k}\}.$$\end{document}

The Green function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{G}}$$\end{document} here includes spin-orbit coupling, but not the exchange field bz—see Methods. We focus on T → 0 for simplicity. Assuming that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0 \,< \,| \lambda e{{\mathcal{E}}}^{z}| {k}_{{\rm{F}}}\,\ll\, \mu$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{{\rm{F}}}=\sqrt{2m\mu }$$\end{document} is the Fermi momentum at zero spin-orbit coupling, a straightforward calculation (provided in the Methods section) yields

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {\bf{J}}\rangle =\frac{{e}^{2}\lambda m}{4\pi }f\left[{\left(\frac{\Delta }{\lambda e{{\mathcal{E}}}^{z}{k}_{{\rm{F}}}}\right)}^{2}\right]{\boldsymbol{\mathcal{E}}}_{\parallel }\times ({b}_{z}{\hat{\bf{z}}})$$\end{document}

with

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f[{y}^{2}]=1-\frac{{y}^{2}}{\sqrt{1+{y}^{2}}}{\rm{arccoth}}(\sqrt{1+{y}^{2}}).$$\end{document}

In the limit of small Δ (so that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0\, < \,\Delta \,\ll\, | \lambda e{{\mathcal{E}}}^{z}| {k}_{{\rm{F}}}\,\ll\, \mu$$\end{document} ), one has \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f[\left({\Delta }^{2}/{(\lambda e{{\mathcal{E}}}^{z}{k}_{{\rm{F}}})}^{2}\right.]\to 1$$\end{document} . This yields the simple expression

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {\bf{J}}\rangle =\frac{{e}^{2}\lambda m}{4\pi }{\boldsymbol{\mathcal{E}}}_{\parallel }\times ({b}_{z}{\hat{\bf{z}}}).$$\end{document}

In spite of the simplicity of Eq. (19), our result is in fact highly singular. For small but nonzero \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{z}$$\end{document} and Δ (in the order of limits given above), the current is independent of both those parameters, but it vanishes if either one is zero. See the Methods section for details. We now match our microscopic calculation to the Ginzburg-Landau theory to determine the coefficient γ appearing in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\bf{A}}}$$\end{document} :

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{j}}}_{\gamma }=-4{e}^{2}\kappa {\tilde{\bf{A}}}=\gamma {\boldsymbol{\mathcal{E}}}_{\parallel }\times ({b}_{z}{\hat{\bf{z}}})=\langle {\bf{J}}\rangle ,$$\end{document}

i.e.,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma =-\frac{{e}^{2}\lambda m}{4\pi }f\left[{\left(\frac{\Delta }{\lambda e{{\mathcal{E}}}^{z}{k}_{{\rm{F}}}}\right)}^{2}\right]\mathop{\longrightarrow }\limits^{\Delta \ll | \lambda e{{\mathcal{E}}}^{z}| {k}_{{\rm{F}}}}-\frac{{e}^{2}\lambda m}{4\pi }.$$\end{document}

We conclude that in the given limit,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{M}_{\gamma }}{{M}_{{\rm{n}}}}=-\frac{{e}^{2}\lambda {V}_{\parallel }}{12{\pi }^{2}{\mu }_{{\rm{B}}}}.$$\end{document}

To extend our results beyond perturbation theory, and to nonzero temperature, we numerically compute the average current 〈J〉 by diagonalizing the 4 × 4 Hamiltonian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}_{{\bf{k}}}$$\end{document} . Denoting the energy eigenvalues as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon_{{\bf{k}}}^{(n)}$$\end{document} and eigenstates as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$| {\psi }_{{\bf{k}}}^{(n)}\rangle$$\end{document} , the average current is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {\bf{J}}\rangle =\frac{1}{{L}^{2}}\mathop{\sum}\limits_{{\bf{k}}}\mathop{\sum}\limits_{n}\left\langle {\psi }_{{\bf{k}}}^{(n)}\right\vert {{\bf{J}}}_{{\bf{k}}}\left\vert {\psi }_{{\bf{k}}}^{(n)}\right\rangle {n}_{{\rm{F}}}\left({\mathcal{E}}_{{\bf{k}}}^{(n)}\right).$$\end{document}

nF is the Fermi-Dirac distribution at temperature T. Fig. 2 displays our analytical and numerical results. As shown in Fig. 2b, the singular behavior revealed by our analytical calculation at T = 0 is smoothed out at T > 0.Fig. 2. Analytical (curves) and numerical (symbols) results for the current density 〈J^y^〉 given the Hamiltonian of Eq. (14), as a function of (a) varying vertical electric field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{z}$$\end{document} at fixed superconducting gap Δ, and (b) varying Δ at fixed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{z}$$\end{document} .Both panels use \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda e{{\mathcal{E}}}^{x}/{v}_{{\rm{F}}}=-0.01/\sqrt{20}$$\end{document} and b/μ = 0.001. All analytical results use the T → 0 limit. Numerics in panel (a) use T/μ = 0.0001. Panel (b) uses \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda e{{\mathcal{E}}}^{z}/{v}_{{\rm{F}}}=0.5/\sqrt{20}$$\end{document} .

Comparison to 4Hb-TaS2

As mentioned in the introduction, our study was motivated by recent experiments on 4Hb-TaS_2_^1,29–33^. Most relevant to our study is the observation of hidden magnetic memory^1^. We begin this section by summarizing that observation. We then argue that while our mechanism can explain some qualitative features of the experiments, the magnitude of the effect predicted by our theory seems too small to account for the experimental observations.

Ref. ^1^ reported the appearance of spontaneous vortices in the superconducting state of TaS_2_, despite the lack of any applied magnetic flux during superconducting condensation. In particular, vortices were observed if a training field was applied only above the superconducting transition temperature Tc; the training field was turned off before cooling down through Tc. The effect persists as long as training is performed below a second temperature T^^ > Tc. The most obvious explanations are an inhomogeneous superconducting transition, or straightforward ferromagnetism present above Tc^37,38^. Both are ruled out by scanning SQUID magnetometry at intermediate temperatures Tc < T < T^^, which places an upper bound on the magnetization in the normal state: that magnetization is smaller than would be necessary to explain the density of vortices observed in the superconducting state by at least a factor of 10^3^.

The question demanding explanation is then: Why does the total magnetization of 4Hb-TaS_2_ appear to increase drastically when entering the superconducting state? Several answers have been proposed, based on the assumption of: (i) alternating layers of superconductor and either a chiral^39^ or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document} ^40^ quantum spin liquid; (ii) unconventional interlayer pairing^41^, or; (iii) a type-II heavy Fermi liquid^42^. No explanation yet has unambiguous experimental support, so the question remains open.

Consider again the scenario depicted in Fig. 1b, where the material hosts a finite density of inclusions which are both charged and magnetic, and which hence exert a local electric field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}$$\end{document} and a local exchange field b on the itinerant electrons. In the context of 4Hb-TaS_2_, we suppose that these inclusions are weakly coupled to one another via the metal into which they are embedded, and that they order (either into a ferromagnet or a glass) at T^*^. Then, the training field induces some nonzero average exchange field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{b}}={b}_{z}{\hat{\bf{z}}}$$\end{document} , but the associated magnetization is small enough to evade detection in SQUID magnetometry performed above Tc.

In the superconducting state, the spin-orbit coupling induced by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}$$\end{document} then combines with b to produce supercurrent through the mechanism we have identified, which (for an appropriate sign of γ and the charge of the inclusions) increases the total magnetization below Tc. Moreover, if the distance between the inclusions is smaller than the London penetration depth, the superconductor will respond to the magnetic field generated by these supercurrents by creating vortices that are randomly located in between the inclusions. Since the vortices are not pinned to the location of the inclusions, they may appear in different locations on each cooldown (assuming that the inclusions occupy only a small fraction of the sample area, leaving a large area for vortices to explore). Both of these observations are consistent with the experiments: (i) The spontaneous vortices reflect an enhanced magnetization relative to the normal state, and (ii) the locations of the vortices are random and are not pinned to any particular location.

To adequately explain the experiment, however, the magnetization corresponding to the spontaneous vortices generated by our mechanism would need to be several orders of magnitude larger than that corresponding to the magnetization in the normal state. Otherwise, our picture cannot explain why the magnetic order remained invisible when Tc < T < T^^. In the Methods section, we estimate the ratio M_γ*/Mn in Eq. (10) and find that it is likely to be at most of order unity. Therefore, the experimentally observed magnetization enhancement in 4Hb-TaS_2 seems larger than our mechanism can naturally explain. However, the electronic structure of that material is considerably different from our toy model, for example due to its nonzero Berry curvature, multiple bands crossing the Fermi energy, and strong Ising-type spin-orbit coupling. A full microscopic accounting for these properties lies beyond the scope of this work.

We note that the previously-studied α and β terms, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{j}}}_{\alpha }+{{\bf{j}}}_{\beta }=\alpha {\hat{\bf{z}}}\times {\bf{b}}+\beta {\boldsymbol{\nabla}}\times ({b}_{z}{\hat{\bf{z}}})$$\end{document} , are even more unlikely to explain the experiment in question. In that experiment, the training field is applied along z, so one expects \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\bf{z}}}\times {\bf{b}}=0$$\end{document} . The α current jα then vanishes. The β current, jβ, does not require in-plane b, but is unlikely to explain the experiment for two reasons. First, β is not required to vanish in the normal state^43^. Second, we estimate that under reasonable conditions, the magnetization contribution due to β is parametrically smaller than that due to γ; see Methods.

It is also worth mentioning that ref. ^16^ proposed a mechanism of generating spontaneous vortices bound to magnetic impurities in iron-based superconductors. This mechanism is unlikely to explain the observations in 4Hb-TaS_2_ since, as mentioned above, the spontaneous vortices in that system do not seem to be pinned to any particular location in the sample.

Discussion

We have identified a new mechanism by which charged magnetic inclusions can drive anomalous currents in spin-orbit coupled superconductors. Our results expand the family of known mechanisms by which the interplay of local magnetism and spin-orbit coupling can produce currents in the superconducting state.

Concerning the problem of magnetic memory in 4Hb-TaS_2_, by assuming that the London penetration depth is much larger than both the size of the inclusions and the distance between them, we argued that our mechanism can indeed drive the spontaneous formation of vortices in between the inclusions, even when the external training field is removed before cooling into the superconducting state. However, within the model we studied, we found that the effect is unlikely to be large enough to explain how the magnetic memory in 4Hb-TaS_2_ remains hidden in the normal state, given the observed vortex density in the superconducting state. Future work should study whether more realistic electronic structure could yield a larger effect.

Lastly, we point out that the dependence of our mechanism on inclusions suggests a possible experimental signature. Since the spontaneous vortices we have identified avoid the inclusions which produce them, repeated cooldowns should yield different random vortex distributions when the inclusion density is low. Increasing the density of inclusions progressively restricts the low-energy configuration space available to the vortices, ultimately pinning them to fixed locations even on repeated cooldowns, provided that the exchange field configuration does not change.

Methods

The electric field distribution

In the Results section, we assumed a charged inclusion which exerts a radial electric field near its edge and not elsewhere. In this subsection we justify that assumption. The total screened potential is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${V}_{{\rm{tot}},{\bf{q}}}=\frac{{V}_{{\rm{ext}},{\bf{q}}}}{1-{e}^{2}{U}_{{\bf{q}}}{\chi }_{{\bf{q}}}},$$\end{document}

where Vext,q = Uqρext,q is the unscreened potential exerted by the charged inclusion (corresponding to the charge density ρext,q), Uq = 1/(ϵ0∣q∣) is the 2D Fourier transform of the bare Coulomb interaction, and χq is the static charge response function. In the Thomas-Fermi approximation (valid at small q, and all the way up to ∣q∣ = 2kF in the 2D Fermi gas), the response function χ is independent of q, allowing us to write

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${V}_{{\rm{tot}},{\bf{q}}}=\frac{{\rho }_{{\rm{ext}},{\bf{q}}}/({\mathbf{\epsilon}}_{0}| {\bf{q}}| )}{1+{k}_{{\rm{TF}}}/| {\bf{q}}| }\approx \frac{{\rho }_{{\rm{ext}},{\bf{q}}}}{{{\mathcal{E}}}_{0}{k}_{{\rm{TF}}}}.$$\end{document}

The last approximation is valid for momenta much smaller than the Thomas-Fermi inverse screening length, ∣q∣ ≪ kTF. If we assume that the external charge distribution varies slowly compared to the Thomas-Fermi screening length, such that ρext,q vanishes when ∣q∣ approaches kTF, and we have Vtot ≈ ρext/(ϵ0kTF) in coordinate space as well. Assuming that the external charge distribution is constant inside the inclusion and decays to zero over a region of width d ≫ 1/kTF, we obtain the in-plane field configuration used in the Results section.

The free energy with vortices

The full Ginzburg-Landau free energy is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F=\displaystyle\int{{\rm{d}}}^{2}{\bf{r}}{\mathcal{F}}=\displaystyle\int{{\rm{d}}}^{2}{\bf{r}}\left[\frac{\kappa }{2}{\left({\boldsymbol{\nabla }}\theta +2e{\bf{A}}+2e{\tilde{\bf{A}}}\right)}^{2}+\frac{1}{2{\mu }_{0}}{({\boldsymbol{\nabla}}\times {\bf{A}})}^{2}\right].$$\end{document}

Using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\bf{A}}}=-{\boldsymbol{\nabla}}\times {\tilde{\bf{h}}}/(4{e}^{2}{\mu }_{0}\kappa )$$\end{document} as in the Results section, we then write

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{lll}F=\displaystyle\int{{\rm{d}}}^{2}{\bf{r}}\left[\displaystyle\frac{\kappa }{2}{\left({\boldsymbol{\nabla }}\theta +2e{\bf{A}}\right)}^{2}-\displaystyle\frac{1}{2e{\mu }_{0}}\left({\boldsymbol{\nabla }}\theta +2e{\bf{A}}\right)\cdot ({\boldsymbol{\nabla }}\times {\tilde{\bf{h}}})\right.\\\qquad +\,\left.\displaystyle\frac{1}{8{e}^{2}{\mu }_{0}^{2}\kappa }{({\boldsymbol{\nabla }}\times {\tilde{\bf{h}}})}^{2}+\displaystyle\frac{1}{2{\mu }_{0}}{({\boldsymbol{\nabla}}\times {\bf{A}})}^{2}\right].\end{array}$$\end{document}

We integrate the second term by parts and neglect the resulting surface term (justified since we assume that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\bf{h}}}$$\end{document} vanishes outside the material). Dropping the constant (non-dynamical) third term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim {({\boldsymbol{\nabla}}\times {\tilde{\bf{h}}})}^{2}$$\end{document} , we find

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F=\displaystyle\int{{\rm{d}}}^{2}{\bf{r}}\left[\frac{\kappa }{2}{\left({\boldsymbol{\nabla }}\theta +2e{\bf{A}}\right)}^{2}-\frac{1}{2e{\mu }_{0}}{\tilde{\bf{h}}}\cdot {\boldsymbol{\nabla}}\times \left({\boldsymbol{\nabla }}\theta +2e{\bf{A}}\right)+\frac{1}{2{\mu }_{0}}{({\boldsymbol{\nabla }}\times {\bf{A}})}^{2}\right].$$\end{document}

Writing the vortex density nv = ∇− × ***∇*θ/(2π), collecting terms, and freely adding a constant term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim {\tilde{\bf{h}}}^{2}$$\end{document} , we arrive at Eq. (12):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F=\displaystyle\int{{\rm{d}}}^{2}{\bf{r}}\left[\frac{\kappa }{2}{\left({\boldsymbol{\nabla }}\theta +2e{\bf{A}}\right)}^{2}+\frac{\pi }{e{\mu }_{0}}{{\bf{n}}}_{v}\cdot {\tilde{\bf{h}}}+\frac{1}{2{\mu }_{0}}{({\boldsymbol{\nabla}}\times {\bf{A}}-{\tilde{\bf{h}}})}^{2}\right].$$\end{document}

Current in the q → 0 limit

On general grounds, there can never be a nonzero spatially-averaged current in equilibrium, so computing the current under a uniform field configuration as in the second subsection of the Results appears suspicious. However, in the superconducting state, a loophole allows us to estimate the current response by working at q = 0 in the microscopic model. The argument is the same as for the usual electromagnetic response of a superconductor. First, combine the true vector potential and our effective vector potential into \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{A}}}_{{\rm{net}}}={\bf{A}}+{\tilde{\bf{A}}}$$\end{document} . Within linear response, the current is then given by a response tensor Kq: jq = Kq[Anet,q + iqθq/(2e)]. In a gapped superconducting state, Kq is continuous and nonzero at q → 0, so we may estimate its value at small but nonzero q as approximately equal to Kq=0. The fact that Kq=0 ≠ 0 does not violate the requirement of zero net current at equilibrium: for a uniform arrangement of Anet, the superconducting phase θ will vary linearly in space such that 2eAnet + ***∇***θ = 0. Our microscopic calculation is then justified: we assume a uniform arrangement of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}$$\end{document} and b in the microscopic model, and compute the current neglecting the ***∇***θ contribution (which cancels out the q = 0 current in reality).

The above argument also explains why the α and γ terms, corresponding to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{j}}}_{\alpha }=\alpha {\hat{\bf{z}}}\times {\bf{b}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{j}}}_{\gamma }=\gamma {\boldsymbol{\mathcal{E}} }_{\parallel }\times ({b}_{z}\hat{{\bf{z}}})$$\end{document} , are forbidden in the normal state: without the ***∇***θ contribution to the current, both of these terms would allow for nonzero spatially-averaged current in equilibrium. On the other hand, the β term is permitted in the normal state: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{j}}}_{\beta }=\beta {\boldsymbol{\nabla}}\times ({b}_{z}{\hat{\bf{z}}})$$\end{document} contributes no current at q = 0, due to the presence of the curl operator.

The Green function at zero exchange field

At zero exchange field (bz = 0), one can express the Green function cleanly as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{G}}}_{k}=\sum _{\alpha =\pm }{{{\Pi }}}_{{\bf{k}}\alpha }{{\mathcal{G}}}_{k\alpha },$$\end{document}

with projectors

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\Pi }}}_{{\bf{k}}\alpha }=\frac{1}{2}\left({\sigma }^{0}+\alpha {\hat{\bf{d}}}_{{\bf{k}}}\cdot {\boldsymbol{\sigma }}\right){\tau }^{0}$$\end{document}

and Green function components

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{G}}}_{k\alpha }={\sigma }^{0}\frac{i{\omega }_{n}{\tau }^{0}+{\xi }_{{\bf{k}}\alpha }{\tau }^{z}-\Delta {\tau }^{x}}{{(i{\omega }_{n})}^{2}-{E}_{{\bf{k}}\alpha }^{2}}.$$\end{document}

The spin-orbit pseudovector is dk = λ**e \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal{E}}$$\end{document} × k, the spin-orbit-coupled band energies (relative to μ) are

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\xi }_{{\bf{k}}\alpha }=\frac{{{\bf{k}}}^{2}}{2m}-\mu -\alpha | {{\bf{d}}}_{{\bf{k}}}| ,$$\end{document}

and the Bogoliubov quasiparticle energies are

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E}_{{\bf{k}}\alpha }=\sqrt{{\xi }_{{\bf{k}}\alpha }^{2}+{\Delta }^{2}}.$$\end{document}

α = ± labels the spin-orbit-coupled bands according to their eigenvalues of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\bf{d}}}\cdot {\boldsymbol{\sigma }}$$\end{document} .

Computing the current

The trace in 〈J〉 factorizes over particle/hole and spin indices using the representation of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{G}}$$\end{document} from the preceding subsection:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{ll}\langle {\bf{J}}\rangle \,=\,-e\frac{{b}_{z}}{4}\frac{T}{{L}^{2}}\mathop{\sum}\limits _{k}{\rm{tr}}\left\{\left[\frac{{\bf{k}}}{m}{\sigma }^{0}-\lambda e{({\boldsymbol{\sigma }}\times {\boldsymbol{\mathcal{E}}} )}_{x,y}\right]{\tau }^{0}{{\mathcal{G}}}_{k}{\sigma }^{z}{\tau }^{0}{{\mathcal{G}}}_{k}\right\}\\ \qquad\,=\,-e\frac{{b}_{z}}{4}\frac{T}{{L}^{2}}\mathop{\sum}\limits _{k}\mathop{\sum}\limits _{\alpha \beta }{{\rm{tr}}}_{\sigma }\left\{\left[\frac{{\bf{k}}}{m}{\sigma }^{0}-\lambda e{({\boldsymbol{\sigma }}\times {\boldsymbol{\mathcal{E}}} )}_{x,y}\right]{{{\Pi }}}_{{\bf{k}}\alpha }{\sigma }^{z}{{{\Pi }}}_{{\bf{k}}\beta }\right\}{{\rm{tr}}}_{\tau }\left\{{{\mathcal{G}}}_{k\alpha }{{\mathcal{G}}}_{k\beta }\right\}.\end{array}$$\end{document}

The Matsubara summation can be performed using standard methods^44^. For the α = β terms,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}\mathop{\sum}\limits_{i{\omega }_{n}}{{\rm{tr}}}_{\tau }\left\{{{\mathcal{G}}}_{k\alpha }{{\mathcal{G}}}_{k\alpha }\right\}=n_{\rm{F}}^{\prime} ({E}_{{\bf{k}}\alpha })$$\end{document}

Note that at T → 0, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n^{\prime} ({E}_{{\bf{k}}\alpha })=-\delta ({E}_{{\bf{k}}\alpha })=0$$\end{document} as long as Δ ≠ 0. In the normal state, the δ-function results in a Fermi surface integral which must be included to obtain the necessary cancellation of the total current; see the next subsection. For the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha ={\overline{\beta }}$$\end{document} terms,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}\mathop{\sum}\limits_{i{\omega }_{n}}{{\rm{tr}}}_{\tau }\left\{{{\mathcal{G}}}_{k+}{{\mathcal{G}}}_{k-}\right\}\mathop{\longrightarrow }\limits^{T\to 0}\frac{1}{{E}_{+}+{E}_{-}}\left(\frac{{\xi }_{+}{\xi }_{-}+{\Delta }^{2}}{{E}_{+}{E}_{-}}-1\right).$$\end{document}

For the spin sector,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathop{\sum}\limits_{\alpha =\pm }{{\rm{tr}}}_{\sigma }\left\{\left[\frac{{\bf{k}}}{m}{\sigma }^{0}-\lambda e{({\boldsymbol{\sigma }}\times {\boldsymbol{\mathcal{E}}} )}_{x,y}\right]{{{\Pi }}}_{{\bf{k}}\alpha }{\sigma }^{z}{{{\Pi }}}_{{\bf{k}}{\overline{\alpha}}}\right\}=2e\lambda {({{\mathcal{E}}}^{z})}^{2}{{\mathcal{E}}}^{x}\frac{{k}_{x}{k}_{y}\hat{{\bf{x}}}-{k}_{x}^{2}\hat{{\bf{y}}}}{{({{\mathcal{E}}}^{z}{\bf{k}})}^{2}+{({{\mathcal{E}}}^{x}{k}_{y})}^{2}}.$$\end{document}

The x component is odd in kx and ky, so it vanishes in the momentum integral and we may safely ignore it. For the surviving y component, the overall factor of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{x}$$\end{document} means that to leading order, we can set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{x}=0$$\end{document} inside the integral, and the only angle dependence in the integrand is through the overall factor of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${({k}_{x}/k)}^{2}={\cos }^{2}\theta$$\end{document} . We then have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {\bf{J}}\rangle =-{e}^{2}\frac{{b}_{z}\lambda {{\mathcal{E}}}^{x}}{2}\mathop{\displaystyle\int}\nolimits_{-\pi }^{\pi }\frac{{\rm{d}}\theta }{2\pi }{\cos }^{2}\theta \mathop{\displaystyle\int}\nolimits_{0}^{\infty }\frac{{\rm{d}}kk}{2\pi }\frac{1}{{E}_{+}+{E}_{-}}\left[1-\frac{{\xi }_{+}{\xi }_{-}+{\Delta }^{2}}{{E}_{+}{E}_{-}}\right]{\hat{\bf{y}}}$$\end{document}

where ξ± and E± are at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{x}=0$$\end{document} . Note that if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{z}=0$$\end{document} , then the integrand is identically zero (since ξ+ = ξ−) and the current vanishes, as mentioned in the Results section. Physically, without the Rashba coupling generated by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{z}$$\end{document} , the exchange field cannot polarize the singlet Cooper pairs (unless it is strong enough to break them entirely, in which case the current vanishes anyways), and the mechanism we have described in the Introduction and Results sections cannot act.

To evaluate the integral, we linearize around \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{{\rm{F}}}=\sqrt{2m\mu }$$\end{document} , writing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k={k}_{{\rm{F}}}+k^{\prime}$$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k^{\prime} \,\ll\, {k}_{{\rm{F}}}$$\end{document} . In this case,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\xi }_{\pm }=\frac{1}{2m}{({k}_{{\rm{F}}}+k^{\prime} )}^{2}\mp \lambda e{{\mathcal{E}}}^{z}({k}_{{\rm{F}}}+k^{\prime} )={v}_{{\rm{F}}}k^{\prime} \mp {\tilde{v}}({k}_{{\rm{F}}}+k^{\prime} )\approx ({v}_{{\rm{F}}}\mp {\tilde{v}})k^{\prime}$$\end{document}

where the last approximation used \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{v}}=| \lambda e{{\mathcal{E}}}^{z}| \,\ll\, {v}_{{\rm{F}}}={k}_{{\rm{F}}}/m$$\end{document} . Then,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {\bf{J}}\rangle \approx -\frac{{e}^{2}{b}_{z}\lambda {{\mathcal{E}}}^{x}m}{4\pi }f\left[{\left(\frac{\Delta }{\lambda e{{\mathcal{E}}}^{z}{k}_{{\rm{F}}}}\right)}^{2}\right]{\hat{\bf{y}}}$$\end{document}

with

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{l}f[\,{y}^{2}]=\displaystyle\frac{1}{2}\mathop{\int}\nolimits_{-\infty }^{\infty }{\rm{d}}x\frac{1-{x}^{2}-{y}^{2}-\sqrt{{(x-1)}^{2}+{y}^{2}}\sqrt{{(x+1)}^{2}+{y}^{2}}}{\sqrt{{(x-1)}^{2}+{y}^{2}}\sqrt{{(x+1)}^{2}+{y}^{2}}\left(\sqrt{{(x-1)}^{2}+{y}^{2}}+\sqrt{{(x+1)}^{2}+{y}^{2}}\right)}\\ =1-\frac{{y}^{2}}{\sqrt{1+{y}^{2}}}{\rm{arccoth}}(\sqrt{1+{y}^{2}}).\end{array}$$\end{document}

In the limit \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \,\ll\, | \lambda e{{\mathcal{E}}}^{z}| {k}_{{\rm{F}}}$$\end{document} we find the result given by Eq. (19) in the Results section, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {\bf{J}}\rangle =-{e}^{2}{b}_{z}\lambda {{\mathcal{E}}}^{x}m\hat{{\bf{y}}}/(4\pi )$$\end{document} .

Current cancellation in the normal state

Above, we showed that the T → 0 current also has terms containing the δ functions δ(Ekα), which are relevant when Δ = 0 so that Ekα = ∣ξkα∣. The extra piece is then

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{lll}\langle {{\bf{J}}}_{{\rm{extra}}}\rangle =\frac{e{b}_{z}}{4}\displaystyle\int\frac{{{\rm{d}}}^{2}{\bf{k}}}{{(2\pi )}^{2}}\mathop{\sum}\limits _{\alpha }\delta (| {\xi }_{{\bf{k}}\alpha }| ){{\rm{tr}}}_{\sigma }\left\{\left[\frac{{\bf{k}}}{m}{\sigma }^{0}-\lambda e{({\boldsymbol{\sigma }}\times \boldsymbol{\mathcal{E}} )}_{x,y}\right]{{{\Pi }}}_{{\bf{k}}\alpha }{\sigma }^{z}{{{\Pi }}}_{{\bf{k}}\alpha }\right\}\\\qquad\qquad =\frac{e{b}_{z}}{4}\displaystyle\int\frac{{{\rm{d}}}^{2}{\bf{k}}}{{(2\pi )}^{2}}\mathop{\sum}\limits _{\alpha }\delta (| {\xi }_{{\bf{k}}\alpha }| )\left[\alpha \frac{{k}_{y}}{m}{\hat{d}}^{z}-\lambda e{{\mathcal{E}}}^{x}{({\hat{d}}^{z})}^{2}+\lambda e{{\mathcal{E}}}^{z}{\hat{d}}^{x}{\hat{d}}^{z}\right]{\hat{\bf{y}}}.\end{array}$$\end{document}

As in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha ={\overline{\beta}}$$\end{document} terms, the x component again involves an integrand which is odd in kx and ky, and hence vanishes. To be consistent with the rest of our calculation, we take the leading-order term in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{x}$$\end{document} only, which yields

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {{\bf{J}}}_{{\rm{extra}}}\rangle =-\frac{{e}^{2}{b}_{z}\lambda {{\mathcal{E}}}^{x}}{4}\displaystyle\int\frac{{{\rm{d}}}^{2}{\bf{k}}}{{(2\pi )}^{2}}\sum _{\alpha }\delta (| {\xi }_{{\bf{k}}\alpha }| ){\left(\frac{{k}_{y}}{k}\right)}^{2}\left[1-\alpha \frac{k}{| \lambda e{{\mathcal{E}}}^{z}| m}\right]{\hat{\bf{y}}}.$$\end{document}

The integrand is evaluated at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{x}=0$$\end{document} , so the δ functions pick out the Fermi momenta of the Rashba bands

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{\text{F},\alpha }=\sqrt{2m\mu +{(m\lambda e{{\mathcal{E}}}^{z})}^{2}}+\alpha m| \lambda e{{\mathcal{E}}}^{z}| .$$\end{document}

They also give factors of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$| \partial {\xi }_{{\bf{k}}\alpha }/\partial k{| }^{-1}=1/({v}_{{\rm{F}}}-\alpha | \lambda e{{\mathcal{E}}}^{z}| )$$\end{document} , where we assume \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$| \lambda e{{\mathcal{E}}}^{z}|\, < \,{v}_{{\rm{F}}}$$\end{document} so that we can drop the overall absolute-value sign in the denominator. In the limit \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$| \lambda e{{\mathcal{E}}}^{z}|\, \ll \,{v}_{{\rm{F}}}$$\end{document} —the same limit as considered in the superconducting state—we find that the extra contribution in the normal state is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {{\bf{J}}}_{{\rm{extra}}}\rangle =\frac{{e}^{2}{b}_{z}\lambda {{\mathcal{E}}}^{x}m}{4\pi }{\hat{\bf{y}}},$$\end{document}

which is precisely the negative of the term found before. The total current thus vanishes in the normal state, as it must on general grounds.

Estimating the magnetization enhancement

To estimate the size of the magnetization enhancement we have identified relative to the normal state magnetization, we start from Eq. (22). We relate λ to the characteristic Rashba energy \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Delta }_{{\rm{R}}}=| \lambda e{{\mathcal{E}}}^{z}| {k}_{{\rm{F}}}$$\end{document} , and express the out-of-plane electric field in terms of its potential drop \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}^{z}={V}_{\perp }/{d}_{\perp }$$\end{document} . Then, ignoring signs, λ = Δ_R_d⊥/(kFe**V⊥). Since μB ~ e/m, we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{M}_{\gamma }}{{M}_{{\rm{n}}}} \sim \frac{{V}_{\parallel }}{{V}_{\perp }}\frac{m{\Delta }_{{\rm{R}}}{d}_{\perp }}{{k}_{{\rm{F}}}} \sim \frac{{V}_{\parallel }}{{V}_{\perp }}({d}_{\perp }{k}_{{\rm{F}}})\frac{{\Delta }_{{\rm{R}}}}{\mu }.$$\end{document}

We expect V∥ ≲ V⊥ (screening is better in the conducting plane than out), \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${d}_{\perp } \sim [1\,{\text{to}}\,10]\times {k}_{\,{\text{F}}\,}^{-1}$$\end{document} (fields from impurities are probably not felt very far away) and Δ_R_ < μ, so dramatic amplification is unlikely; at most, we expect the extra magnetic moment in the superconducting state to be of similar order to the moment from Pauli paramagnetism.

What about the β term? In the limit \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Delta }_{{\rm{R}}}\,\gg\, m{(\lambda e{{\mathcal{E}}}^{z})}^{2}\,\gg\, \Delta$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta \sim {(m\lambda e{{\mathcal{E}}}^{z})}^{2}/{k}_{\,{\text{F}}\,}^{2}$$\end{document} ^11^. In the geometry we consider, the corresponding current \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{j}}=\beta {\boldsymbol{\nabla }}{b}_{z}\times {\hat{\bf{z}}}$$\end{document} yields a magnetic moment

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${M}_{\beta }=\pi {R}^{2}\beta {b}_{z},$$\end{document}

so

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{M}_{\gamma }}{{M}_{\beta }} \sim \frac{\gamma {V}_{\parallel }}{\beta } \sim \frac{{V}_{\parallel }}{{V}_{\perp }}({d}_{\perp }{k}_{{\rm{F}}})\frac{\mu }{{\Delta }_{{\rm{R}}}}.$$\end{document}

Note the factor of μ/Δ_R_, which is likely quite large. In the presence of an in-plane electric field (i.e., around impurities which are both charged and magnetic), our γ term then likely plays a larger role than the β term identified previously. The β term, on the other hand, does not require an in-plane electric field; neutral magnetic impurities suffice. It is interesting to note that for our model, the β term vanishes in the normal state, even though this is not required on general grounds. For example, if a term ∝ ∣k∣^4^ is added to the bare dispersion relation, then the β term survives in the normal state.