Theory of magnetic response in finite two-dimensional superconductors

We present a theory of magnetic response in a finite-size two-dimensional superconductors with Rashba spin-orbit coupling. The interplay between the latter and an in-plane Zeeman field leads on the one hand to an out-of-plane spin polarization which accumulates at the edges of the sample over the superconducting coherence length, and on the other hand, to circulating supercurrents decaying away from the edge over a macroscopic scale. In a long finite stripe of width W both, the spin polarization and the currents, contribute to the total magnetic moment induced at the stripe ends. These two contributions scale with W and W2 respectively, such that for sufficiently large samples it can be detected by current magnetometry techniques.

Superconductivity in two-dimensional (2D) and quasi-2D systems has been attracting a great deal of attention over past decades [1,2]. Examples of such systems range from ultra-thin metallic films, heavy fermion superlattices, and interfacial superconductors to atomic layers of metal dichalcogenides, and organic conductors.
Most 2D superconductors exhibit large spin-orbit coupling (SOC) because of broken space inversion symmetry. In this regard two types of 2D superconductors can be distinguished: Those exhibiting SOC of Rashba-type due to a broken up-down (out-of-plane) mirror symmetry, denoted here as Rashba superconductors, and those, in which a 2D in-plane inversion symmetry is broken due to a non-centrosymmetric crystal structure. The latter are exemplified by 2D transition metal dichalcogenides [3,4]. To the first group, on which we focus here, belong for example ultra-thin superconducting metallic films [5][6][7].
In this letter we demonstrate that finite size effects drastically modify the magnetic response of Rashba superconductors leading to hitherto unknown phenomena. Our main findings are the following: (i) In a response to a Zeeman field the system exhibits a spin texture (Fig.  1a) with a transverse component of the spin localized near the edge on the scale of superconducting coherence length. (ii) Because of the spin-charge coupling mediated by the SOC, a non-homogenous charge current appears in the system with a spatial distribution that depends on the direction of the applied field and geometry of the system (Fig. 1b); (iii) In particular, for a finite stripe oriented along the field, macroscopic currents loops appear at the stripe ends (Fig. 1c). Both, the transverse spin and the edge currents contribute to the total magnetic moment which can be detected by state-of-the-art magnetometry techniques.
These findings can be qualitatively understood recalling the concepts of spin currents and spin galvanic effect (see Fig. 1). The key feature of 2D materials without up-down mirror symmetry is the Rashba SOC, Here e z is a vector normal to the 2D plane, v = p/2m is the quasiparticle velocity, m its effective mass, σ is the vector of Pauli matrices, and α is the SOC constant [22]. The SOC acts as an effective p-dependent spin splitting field. Let us assume that the system is subject to an external Zeeman field B, and for some reason the induced spin polarization S differs from the equilibrium Pauli response χ P B, where χ P is the Pauli paramagnetic polarizability. Then the excess spin δS = S − χ P B will experience an inhomogeneous precession in the effective Rashba field, generating a momentum anisotropy of the density matrix. In the presence of disorder the precession rateR = i[H R , δS · σ] is balanced by the momentum relaxation, which results in a steady spin current in the bulk of the system J a bulk,k = −τ tr v k σ aR = αD(δS z δ ka − δS k δ az ), where τ is the momentum relaxation time, and D = τ v 2 F /2 is the diffusion coefficient. Under equilibrium conditions δS = 0 in normal systems, but in superconductors pair correlations modify the Pauli response leading to a finite δS [12,23]. This leads to finite equilibrium spin-currents in Rashba superconductors generated by the Zeeman field. For example, a field applied in x-direction in a bulk superconductor produces a spin-current with an out-of-plane polarization, J z bulk,x = −αDδS x . Due to the spin-Hall magneto-electric coupling in Rashba materials the bulk spin-current generates a transverse charge current according to j bulk,y ∝ α 2 J z bulk,x = α 3 DδS x , which is nothing, but the anomalous supercurrent well known for bulk superconductors with SOC [8,11,17]. In a finite system currents must vanish at the edges of the sample. This condition can be fulfilled only if the distribution of the excess spin δS(r) is inhomogeneous near the edge, so that the diffusion spin-current J a dif f,k = −D∂ k δS a compensates the bulk contribution. For concreteness, if we assume a boundary with vacuum at x = 0, the zero spin-current condition for a field applied in x-direction reads: D∂ x δS z = J z bulk,x , which implies that a finite component δS z (x) transverse to the field is induced at the edges of the sample. In this case the spin density exhibits a texture as sketched in Fig. 1a [24]. In the presence of SOC both the edge and the bulk spin-currents are converted into a charge current flowing parallel to the boundary, via the spin-galvanic effect, Fig. 1b. In a realistic finite system currents must vanish at all edges. The anomalous charge currents at the boundaries should then be compensated by supercurrents which stem from a gradient of the superconducting phase. As a consequence, in a stripe geometry, an in-plane field induces current loops at the edges as shown in Fig. 1c. The magnetic moment induced by this currents and by the transverse spin can in principle be measured to directly detect the effects we predict here. In the rest of the paper we provide a quantitative derivation of these effects, calculate the induced magnetic moment, and propose materials in which our predictions can be verified.
Specifically, we consider a 2D disordered superconductor with Rashba SOC. We assume that the Fermi energy corresponds to the largest energy scale, so that spectral and transport properties can be accurately described by the quasiclassical Green's functions (GFs) [25,26]. In the diffusive limit these functions are isotropic in momentum and they obey the Usadel equation which in the presence of a Zeeman field and Rashba SOC reads [27][28][29]: Here h = µ B B, σ = (σ x , σ y , σ z ) and τ 2,3 are Pauli matrices spanning spin and Nambu space, respectively, ω is the Matsubara frequency, ∆ is the superconducting order parameter, and SOC enters via the covariant derivativẽ summation over repeated indices is implied, and k = x, y [30]. The quasiclassical GFǧ in Eq. (1) is a 4×4 matrix in the Nambu-spin space, which satisfies the normalization conditionǧ 2 = 1. In the absence of spindependent fields it readsǧ 0 = (ω/E) It is easy to check by substitution into Eq. (1), that in the absence of Zeeman fieldǧ 0 is also the solution of the Usadel equation for arbitraryÂ k .
To compute the response to an external magnetic field we linearize Eq. (1) with respect to h and write the solution asǧ ≈ǧ 0 + δǧ. It is convenient to define δǧ ≡ iǧ 0 [τ 3 ,ǧ 0 ]Q, whereQ is a matrix in spin space that satisfies the following equation [31]: The excess spin density δS a is then determined by [32]: For a homogeneous infinite 2D superconductor, the solu- Equations (3)- (5) reproduce the bulk spin response of Rashba superconductor [12,14], which is finite even at T = 0 and depends on the direction of the applied field. This situation changes drastically in a finite system. First, we assume that the system is infinite in y-direction, and bounded to the region |x| < L/2 in the x-direction. The solution to Eq. (2) can be written as the sum of the bulk contribution and a contribution from the sample edges,Q =Q b + δQ(x). According to Eq. (2) the latter satisfies: The last terms in these equations describe precession of the excess spin, caused by SOC. Importantly, the precession terms are finite only for inhomogeneous systems. The boundary conditions to the above equations are obtained by imposing zero-current at the edges, x = ±L/2 [28,33]: Here the left hand side is proportional to the inhomogeneous spectral spin-current which cancels the bulk one in the right hand side. The boundary problem of Eqs. (7)-(8) has a nontrivial solution only if the right-hand-side in Eq. (8), that is the bulk spin-current, is finite. According to Eqs. (4)-(5), this is the case when the magnetic field has either z-or x-components. How to obtain the solution for δQ a is discussed in the Supplementary Material.
Here we present the spatial dependence of the induced spin obtained from Eq. magnetic response for Rashba superconductors [12,14] to the case of finite samples. In addition to the finite spin response at T = 0, the SOC in superconductors also leads to the spin-galvanic effect, that is, a creation of charge currents by a Zeeman field [11,17,20,34]. In the stripe geometry (see middle panels of Fig. 2) the so called anomalous charge current is induced in y-direction, j an y = (θ/m)(∂ x δS z − 2αδS x ) [33], where θ = 2Dτ α 2 is a dimensionless parameter which in normal systems describes the spin-charge conversion [35]. Within the diffusive approximation it is a small parameter which we treat perturbatively. Here δS a is obtained by substituting the solution of Eqs. (6)-(7) into Eq. (3). This results in In the second line we identify two contributions to the anomalous current: the bulk contribution j an b , widely studied in homogeneous superconductors [9,11,17,20] and given by the last term in the brackets in the first line and (red arrows in middle panels in Fig. 2), and the boundary contribution j an edge , determined by the first two terms. The latter are localized at the edges of the sample within the scale of superconducting coherence length (green arrows in middle panels of Fig. 2). In the geometry under consideration, the "bulk" contribution to the current is finite only for fields applied across the stripe (x-direction). The spatial dependence of the charge current density is shown on Fig. 2 b and d for fields in x-and zdirection, respectively. Because of zero spin-current condition, Eq. (8), the charge current of Eq. (9) also vanishes at the boundaries. When the field is applied in x-direction, Fig. 2b, both, the bulk and edge contributions, are finite. The maximal value of the total current is the "bulk" value reached deeply inside the sample, away from the edges. The spatial distribution of the current is symmetric and the net current through the stripe is nonzero. In contrast, if the field is applied in z-direction, Fig. 2d, there is no bulk contribution, because Q z b does not contribute to the current, see Eq. (9). Only edge currents, opposite on opposite sides, appear. Clearly in this case the total charge current vanishes.
The above results apply for an infinite superconducting stripe, and whether the obtained currents may exist in real finite systems depends on transverse boundary conditions. For example, if the finite stripe is wrapped in a cylinder, the periodic boundary conditions indeed allow for the above current patterns when the field is applied along the the cylinder axis [17]. Here we consider a more experimentally relevant situation: a finite 2D superconductor of rectangular shape, which occupies the region |x| < L/2 and |y| < W/2 (see Fig. 3). Obviously, the charge current through all boundaries must vanish. For out-of-plane field this condition is trivially satisfied by closing the boundary streamlines, which generates a circulating edge current. However, in this situation, orbital effects also create circulating diamagnetic supercurrents and the detection of a pure anomalous current is rather difficult. Therefore, here we focus on the case of in-plane field. As shown in Fig. 2(b), in this case the anomalous current has the same direction at both edges. Hence, in the presence of impenetrable boundaries at y = ±W/2, one naturally expects generation of closed streamlines at each edge as shown schematically in Fig. 1(c). In what follows we study this situation quantitatively.
To satisfy the zero-current condition at y = ±W/2, the anomalous current of Eq. (9) needs to be compensated by a supercurrent generated by a gradient of the superconducting phase ϕ(r). The total charge current in a superconductor reads, where n s is the superfluid density. Since the anomalous current is divergenceless, the continuity equation for the current reduces to the Laplace equation for the phase, which has to be solved by imposing a zero current condition at the edges of the sample, n k ∂ k ϕ| edge = − 2m n s n k j an k edge (12) where n is a unit vector normal to the edge. The boundary problem, Eqs. (11)-(12), can be solved following the procedure described in Refs. [35,36]. Here we present the results for a narrow stripe oriented along x-axis with ξ s L W . In this case the currents near opposite edges at x = ±L/2 are independent from each-other, and one can treat each edge separately. For the current near the left edge at x = −L/2 we find, where I edge y is the net anomalous edge current, obtained by integrating j an edge (x) in Eq. (9) from the boundary (x = −L/2) to the bulk of the sample (x = 0). Fig. 3 shows the current streamlines in the whole stripe obtained from these equations. One identifies two different scales: (i) a mesoscopic scale of the order of the coherence length ξ s , over which the anomalous current j an edge (x) decays, and (ii) a macroscopic scale related to the size of the sample. In the present limiting case, ξ s L W , there are two clearly defined contributions to the y-component of the current : the edge anomalous current localized at the boundary (first term in the integrand of Eq. (13) represented by a green arrow in Fig. 3), and the counterflow of supercurrent given by the second term in Eq. (13), which decays over the scale ∼ W and compensates the anomalous edge component such that dxj y (x, y) = 0.
The current stream generates a finite orbital angular momentum L z , see Fig. 1c. The latter is computed from its general definition L z = m dxdy(xj y − yj x )/2, and scales as W 2 : Similarly the total spin angular momentum accumulated at the edge is obtained by integrating the z-component of the spin, Eq. (3), S z = 0 −L/2 δS z (x)dx. In the Supplementary Material we derive the following analytical expressions for both the spin and orbital angular momenta at T = 0 in two limiting cases: and L z ∝ −N F h x θW 2 ξ 0 α, for ξ 0 α 1 (ξ 0 α) −4 ln(ξ 0 α), for ξ 0 α 1 (17) with ξ 0 = D/2∆. The measurable magnetic moment is then given by M z = µ B (L z / + S z ), where µ B is the Bohr magneton [37]. Both contributions have the same sign. However, the spin angular momentum scales with W , while L z scales with W 2 and therefore dominates in macroscopic samples.
In conclusion, we presented a complete theory of the magnetic response of finite size Rashba superconductors. When the field is applied in-plane our theory predicts on the one hand a finite out-of-plane spin polarization localized at the edge of the sample on the scale of superconducting coherence length. On the other hand, the SOC also leads to supercurrents circulating in the sample. Both the spin and the orbital momentum of supercurrents contribute to the total magnetic moment, which is induced at the edges and can be measured by stateof-the-art magnetic sensors [38,39]. Whereas the contribution from the spin angular momentum scales with the width W of a rectangular stripe, the contribution from the currents scales with W 2 and therefore dominates in large samples. There are several superconducting materials with Rashba SOC and whihc our findings can be verified. These range from Pb and Tl-Pb monolayers [6,[40][41][42], to thin MoS 2 , NbRe, β-Bi 2 Pd films [43][44][45], and 2D superconductivity at the LaAlO 3 /SrTiO 3 interface [46][47][48][49].
In the main text all observables, currents, angular momenta, etc, are expressed in terms of δS, which is the deviation of the spin-density from the Pauli paramagnetic spin. According to Eq. (3), the components of δS are determined from the matrixQ. there are two types of contributions: One stemming from the bulk of the sample,Q b and another one from the boundary,δQ. The latter is obtained by solving the boundary problem Eqs. (6)(7)(8). Here we describes the main features of the solution procedure.
The precession term in Eqs. (6)(7) couples the x and z components, δQ x,z . These components are also coupled via the boundary condition, Eq. (8), which written explicitly reads: Clearly, a nontrivial solution exists only if the righthand-side in these equation is non-zero. This is the case when the magnetic field has either z-or x-components.
It is straightforward to check that the general solution for the boundary quantity, δQ a , reads where κ ± = κ 2 E − 2α 2 ± 2iα 7α 2 + 4κ 2 E , κ 2 E = 2E/D, and Reκ ± > 0, is the characteristic length over which δQ decays away from the edge. Then, the transverse spin induces at the edge of the sample decays into bulk over a distance that we denote ξ s . At low temperatures: where ξ 0 = D/2∆. For sufficiently large L, the solution near each boundary can be found independently. For example the solution at the x = −L/2 boundary consists only on the decaying exponents with with The coefficients B z ± are then determined by Eqs. (S1). The solution, Eq. (S4), describes the spin density accumulated at the edge of the sample,x = −L/2, and decays into the bulk over the superconducting coherence length ,and eventually oscillates. In Fig. 2(a,c) of the main text we show the total spin-density obtained by substituting the solution Eq. (S4) into Eq. (3) of the main text.
Appendix B: Calculation of the angular momenta, Eqs. (15)(16)(17) In this section we calculate the total angular momentum at the left edge, x = −L/2 of the stripe shown in Fig.  3 when ξ 0 W L. Because both, the anomalous current and spin densities generated at the boundary, decay into the bulk over the coherence length ξ 0 the integrals in Eqs. (16)(17) of the main text is can be taken in the interval −L/2 < x < 0. We first determined the total orbital angular momentum L z given in Eq. (15) which is proportional to the contribution to the anomalous current stemming from the boundary of the sample, denoted in the main text as I edge y , and defined in Eq. (16). After substituting the expression j an edge from Eq. (9) we obtain: where δQ x =