Superconductivity near a magnetic domain wall

TL;DR

This paper investigates how magnetic domain walls influence superconducting properties in ferromagnetic insulator/superconductor structures, revealing effects on domain wall size, spin currents, and local density of states.

Contribution

It provides a detailed analysis of the impact of magnetic domain walls on superconductivity and related physical quantities in hybrid structures.

Findings

Domain wall size is affected by superconductivity.

Spin current density is altered near the domain wall.

Local density of states and tunneling conductance are modified by the domain wall.

Abstract

We study the equilibrium properties of a ferromagnetic insulator/superconductor structure near a magnetic domain wall. We show how the domain wall size is affected by the superconductivity in such structures. Moreover, we calculate several physical quantities altered due to the magnetic domain wall, such as the spin current density and local density of states, as well as the resulting tunneling conductance into a structure with a magnetic domain wall.

Figures

Equations

$$\begin{array}{rl}\alpha (x)=& \frac{\pi }{\lambda }\left(x+\frac{\lambda }{2}\right)\Theta \left(x+\frac{\lambda }{2}\right)\Theta \left(\frac{\lambda }{2}-x\right)\\ & +\pi \Theta \left(x-\frac{\lambda }{2}\right),\end{array}$$

$$D\nabla \cdot \left(\check{g}\nabla \check{g}\right)-\left[{\omega }_n{\hat{\tau }}_3+i\mathbf{h}\cdot \mathbf{\sigma }{\hat{\tau }}_3+\check{\Delta }+\check{\Sigma },\check{g}\right]=0,$$

$$\Psi ={\left(\begin{array}{cccc}{\psi }_{\uparrow }(x)& {\psi }_{\downarrow }(x)& -{\psi }_{\downarrow }^{†}(x)& {\psi }_{\uparrow }^{†}(x)\end{array}\right)}^{\top },$$

$$\mathbf{h}=h(\sin \alpha (x),0,\cos \alpha (x)),$$

$$\check{R}=e^{i{\hat{\sigma }}_2\alpha (x)/2},$$

$$\check{g}={\check{R}}^{†}{\check{g}}_0\check{R},$$

$$D{\check{\partial }}_x^A\cdot \left({\check{g}}_0{\check{\partial }}_x^A{\check{g}}_0\right)-\left[{\omega }_n{\hat{\tau }}_3+ih{\hat{\sigma }}_3{\hat{\tau }}_3+\check{\Delta }+\check{\Sigma },{\check{g}}_0\right]=0,$$

$${\check{\partial }}_x^{A}X={\partial }_{x}X-\left[A,X\right],$$

$$A=i{\hat{\sigma }}_2{\alpha }^{\prime }(x)/2.$$

$$D\left[A,{\check{g}}_{0}A{\check{g}}_0\right]=-\frac{1}{2}{\alpha }^{\prime }(x{)}^{2}D\left[{\hat{\sigma }}_2{\check{g}}_0{\hat{\sigma }}_2,{\check{g}}_0\right].$$

$$\check{g}_{0}\partial_{x}\check{g}_{0}\biggr{|}_{x_{b}^{\pm}}-\check{g}_{0}\partial_{x}\check{g}_{0}\biggr{|}_{x_{b}^{\mp}}=\frac{1}{2}\alpha^{\prime}(x)\left[\check{g}_{0}i\hat{\sigma}_{2},\check{g}_{0}\right]\biggr{|}_{x_{b}^{\pm}},$$

$$\begin{array}{rl}{\check{g}}_0=& \cos \theta {\hat{\tau }}_3\left(M_0+i\tan \theta \mathbf{M}\cdot \mathbf{\sigma }\right)\\ & +\sin \theta {\hat{\tau }}_1\left(M_0-i\cot \theta \mathbf{M}\cdot \mathbf{\sigma }\right).\end{array}$$

$$M_0^2-|\mathbf{M}{|}^2=1.$$

$$D{\partial }_x^2\theta +2M_0\left(\Delta \cos \theta -{\omega }_n\sin \theta \right)-2h{M}_3\cos \theta -\frac{1}{4{\tau }_{sf}}\left(2M_0^2+1\right)\sin \left(2\theta \right)=0,$$

$$\begin{array}{rl}D\left(\mathbf{M}{\partial }_x^2{M}_0-M_0{\partial }_x^2\mathbf{M}\right)+& 2D{M}_0{\alpha }^{\prime }(x)\left(\begin{array}{c}-{\partial }_x{M}_3\\ 0\\ {\partial }_x{M}_1\end{array}\right)+D{M}_0{\alpha }^{\prime }(x{)}^2\left(\begin{array}{c}{M}_1\\ 0\\ M_3\end{array}\right)+2\mathbf{M}\left(\Delta \sin \theta +{\omega }_n\cos \theta \right)\\ & -2h{M}_0\sin \theta \left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)+\left[\frac{1}{{\tau }_{so}}+\frac{1}{2{\tau }_{sf}}\cos \left(2\theta \right)\right]M_0\mathbf{M}=0.\end{array}$$

$${\left({\partial }_x{M}_0\right)}^2-{\left({\partial }_x{M}_1\right)}^2-{\left({\partial }_x{M}_3\right)}^2+M_0{\partial }_x^2{M}_0-M_1{\partial }_x^2{M}_1-M_3{\partial }_x^2{M}_3=0.$$

$$\begin{array}{rl}D{\partial }_x^2& M_j-2{\alpha }^{\prime }(x)D{M}_j(M_3{\partial }_x{M}_1-M_1{\partial }_x{M}_3)+D{M}_j\sum _{k=0,1,3}(-1{)}^k{\left({\partial }_x{M}_k\right)}^2-{\alpha }^{\prime }(x{)}^{2}D{M}_j(M_1^2+M_3^2)\\ & -2M_j{M}_0({\omega }_n\cos \theta +\Delta \sin \theta )+2h{M}_j{M}_3\sin \theta -\left[\frac{1}{{\tau }_{so}}+\frac{1}{2{\tau }_{sf}}\cos \left(2\theta \right)\right]M_j\left(M_0^2-{\delta }_{j0}\right)=S_j,\end{array}$$

$$S_j=\{\begin{array}{ccc}-2\left({\omega }_n\cos \theta +\Delta \sin \theta \right)& & j=0\\ -2{\alpha }^{\prime }(x)D{\partial }_x{M}_3+{\alpha }^{\prime }(x{)}^{2}D{M}_1& & j=1\\ 2{\alpha }^{\prime }(x)D{\partial }_x{M}_1+{\alpha }^{\prime }(x{)}^{2}D{M}_3-2h\sin \theta & & j=3\end{array}.$$

$$\displaystyle\partial_{x}\theta\bigg{|}_{x_{b}^{\pm}}=\partial_{x}\theta\bigg{|}_{x_{b}^{\mp}}\begin{split}\left(M_{i}\partial_{x}M_{0}-M_{0}\partial_{x}M_{i}\right)\bigg{|}_{x_{b}^{\pm}}-\alpha^{\prime}(x)M_{0}M_{i^{\prime}}\bigg{|}_{x_{b}^{\pm}}=\left(M_{i}\partial_{x}M_{0}-M_{0}\partial_{x}M_{i}\right)\bigg{|}_{x_{b}^{\mp}},\ i=1,3,\end{split}$$

$$\tan$$

$$M_0=\frac{{\omega }_n+\Delta \tan \theta }{\sqrt{{\omega }_n^2+\left({\Delta }^2-h^2\right){\tan }^2\theta +2{\omega }_n\Delta \tan \theta }}$$

$$M_3=\frac{h\tan \theta }{\sqrt{{\omega }_n^2+\left({\Delta }^2-h^2\right){\tan }^2\theta +2{\omega }_n\Delta \tan \theta }}$$

$$M_1=0.$$

$$\begin{array}{rl}\check{g}& =\cos \theta {\hat{\tau }}_3[M_0+i\tan \theta \cos \alpha (x)\mathbf{M}\cdot \mathbf{\sigma }\\ & +i\tan \theta \sin \alpha (x)\left(M_3{\hat{\sigma }}_1-M_1{\hat{\sigma }}_3\right)]\\ & +\sin \theta {\hat{\tau }}_1[M_0-i\cot \theta \cos \alpha (x)\mathbf{M}\cdot \mathbf{\sigma }\\ & -i\cot \theta \sin \alpha (x)\left(M_3{\hat{\sigma }}_1-M_1{\hat{\sigma }}_3\right)].\end{array}$$

$$\Delta =\frac{1}{2}\pi T\gamma \sum _{{\omega }_n>0}^{{\omega }_D}\text{Tr}\left({\hat{\tau }}_1\check{g}\right)=2\pi T\gamma \sum _{{\omega }_n>0}^{{\omega }_D}{M}_0\sin \theta ,$$

$$\begin{array}{rl}2\pi T\sum _{n=0}^{N_D(T)}\frac{1}{{\omega }_n}& =2\pi T\left(\sum _{n=0}^{N_D(T_c)}+\sum _{N_D(T_c)}^{N_D(T)}\right)\frac{1}{{\omega }_n}\\ & =\frac{1}{\gamma }+\log \left(\frac{T_c}{T}\right)\end{array}$$

$$\Delta \log \left(\frac{T}{T_c}\right)=2\pi T\sum _{{\omega }_n>0}\left(M_0\sin \theta -\frac{\Delta }{{\omega }_n}\right).$$

$$j_{k,i}=\frac{{\sigma }_N}{2e^2}\pi Ti\sum _{{\omega }_n>0}\text{Tr}\left[{\hat{\tau }}_0{\hat{\sigma }}_i\left(\check{g}{\nabla }_k\check{g}\right)\right],$$

$$\begin{array}{rl}{j}_{x,2}=\frac{2{\sigma }_N}{e^2}\pi T\sum _{{\omega }_n>0}& [{\alpha }^{\prime }(x)(M_1^2+M_3^2)\\ & +M_3{\partial }_x{M}_1-M_1{\partial }_x{M}_3]\end{array}$$

$$j_{x,1/3}=0,$$

$$\frac{d{\mathbf{n}}_s}{dt}+\nabla \cdot j_{k,i}=\mathcal{T},$$

$$\mathcal{T}={\partial }_x{j}_{x,2}=\frac{4{\sigma }_N}{e^{2}D}\pi Th\sum _{{\omega }_n>0}\sin \theta M_1.$$

Full text

I Introduction

One of the hallmarks of superconductivity is the Meissner effect, where a superconductor expels an external magnetic field from its interior.Meissner and Ochsenfeld (1933) Increasing the magnetic field, superconductivity gets suppressed and eventually destroyed by fields exceeding a critical field. However, the critical field of thin films in parallel magnetic field is greatly increased than that of the bulk superconductors, since the Meissner effect is negligibly small in this case.Tinkham (2004) Rather, the critical field is determined by the paramagnetic effect.Chandrasekhar (1962); Clogston (1962) This effect suppresses superconductivity by aligning the electrons of Cooper pairs to be separated in energy. The magnetic field in this case then leads to a Zeeman effect.

One important consequence of the Zeeman field is the splitting of the density of states (DoS) in energy.Meservey et al. (1970); Tedrow and Meservey (1971); Meservey et al. (1975); Paraskevopoulos et al. (1977); Meservey et al. (1980) In the absence of the Zeeman field, the DoS of a superconductor shows a singularity at the energy which is equal to the superconducting pair potential. This singularity is separated by the Zeeman field for each spin species. The spin splitting can also be induced by making contact of a superconductor (S) with a ferromagnetic insulator (FI).Tedrow et al. (1986); Hao et al. (1990); Moodera et al. (2007); Xiong et al. (2011); Senapati et al. (2011); Wolf et al. (2014) In this case, the ferromagnetic insulator induces a strong exchange field with a small external magnetic field or even in its absence.

Various properties of superconductors with a spin-splitting density of states have been studied in recent years. An example is the strong thermoelectric effect with a thermopower predicted to exceed $k_{B}/e$ and a possibility of obtaining large values of the figure of merit $ZT\gg 1$ at low temperaturesMachon et al. (2013, 2014); Ozaeta et al. (2014). Indications of this effect were also recently detectedKolenda et al. (2016, 2017). A spin accumulation in such a structure was detected at much longer distances than the relaxation lengths at the normal state.Quay et al. (2013); Hübler et al. (2012); Wolf et al. (2013); Silaev et al. (2015) Such effects are reviewed in Refs. Bergeret et al., 2018, 2017.

In the above studies, the spin-splitting field induced in the superconductors is considered to be homogeneous. However, this is not always the case in a FI/S structure due to the nonuniformity of the ferromagnets. The nonuniformity can be described by a domain structure. Since the typical domain size in ferromagnets is much longer than the superconducting coherence length, the uniform magnetization is a good assumption in many cases. Here we relax this assumption and study the effects of inhomogeneity, especially domain walls.

The inhomogeneous exchange field also brings interesting effects to superconductivity. The existence of spin-singlet superconductivity in a magnetic domain wall in a ferromagnet (F) was studied several decades ago.Buzdin et al. (1984) This effect also has been studied in S/F bilayers,Aladyshkin et al. (2003); Houzet and Buzdin (2006) where a superconducting material is placed on top of a ferromagnet with a domain structure. Decreasing the temperature, superconductivity first appears just above the domain wall. The experimental realization of such domain wall superconductivity was reported in Ref. Yang et al., 2004. The generation of a spin supercurrent in Josephson contacts with a domain wall is discussed in Ref. Shomali et al., 2011. The reconstructed density of states at the end of the superconductorBobkova and Bobkov (2019) and the magnetoelectric effectsRabinovich et al. (2018) in S/F bilayers with a magnetic texture were also studied recently. In the case of FI/S structures, a peculiar tunneling conductance was observed in a recent experiment and the magnetic domains in the ferromagnet was considered responsible for the experimental result.Strambini et al. (2017)

Usually, the inhomogeneity of the exchange field induced by the ferromagnetic insulator in the superconductor can be represented by a multi-domain structure, namely, alternating domains with opposite magnetization directions connected via domain walls. The theoretical model of the tunneling conductance in Ref. Strambini et al., 2017 concentrates on this case. Since the size of the domains is often much longer than the superconducting coherence length, whereas here we consider a single domain wall structure separating two domains with opposite magnetization directions. Moreover, the size of the domain wall in Ref. Strambini et al., 2017 was considered much smaller than the superconducting coherence length, here we also consider larger domain walls and study several physical quantities altered due to the inhomogeneity, such as the spin current density and the local density of states.

This paper is organized as follows. We introduce the FI/S structure containing a magnetic domain wall in Sec. II, and solve the Usadel equation with an extended $\theta$-parameterization. Equilibrium spin current density, as a direct consequence of the inhomogeneous exchange field, and its various properties are discussed in Sec. III. In Sec. IV we discuss the effect of superconductivity on the domain wall size, by considering the contribution of the superconducting free energy to the domain wall energy, and show that the maximal relative effect happens when the domain wall size is of the order of $\xi_{0}$, the superconducting coherence length (here defined at zero temperature and exchange field). For the possible measurement accessing this physics, we discuss the local density of states in Sec. V and the tunneling conductance across a nearby tunneling barrier in Sec. VI.

II Model and Method

II.1 Model

We study the properties of the structure in Fig. 1 in equilibrium. A superconducting wire is placed on top of a FI wire containing a magnetic domain wall. In many thin ferromagnets the domain wall structure energetically prefers a Néel type, in which the rotation of the magnetization happens in the plane of the domain wall. Therefore we consider a Néel domain wall with size $\lambda$. We make a variational ansatz and define the magnetization rotation angle as

[TABLE]

where $\Theta$ is the heaviside step function. The choice of $\alpha(x)$ in Eq. (1) as a linear function of $x$ instead of the typically used hyperbolic functions that lead to a somewhat lower energy brings certain technical advantages. Its derivative is a constant inside the domain wall. This simplifies the Usadel equation, which describes the properties of the superconductor (see Sec. II.2). The non-analytic derivative of $\alpha(x)$ at the boundary of the domain wall, can be transferred to the boundary conditions of the Usadel equation (see Sec. II.3). Besides, the domain wall energy introduced by this choice of $\alpha(x)$ gives only slightly larger energy than the rotation angle constructed with hyperbolic functions (see Sec. IV).

The properties of a superconductor in the diffusive limit, namely, in the case when the elastic relaxation rate is much larger than any other energy scales in the problem, are described by the Usadel equation. In equilibrium, it is convenient to use the Usadel equation in the Matsubara representation ($\hbar=k_{B}=1$)

[TABLE]

where $D$ is the diffusion constant, and $\check{g}$ is the quasiclassical Green’s function satisfying the normalization condition $\check{g}^{2}=\check{1}$. In the commutator, $\omega_{n}=(2n+1)\pi T$ is the Matsubara frequency, $T$ is the temperature and $n$ is an integer, $\bm{h}$ is the exchange field, $\bm{\sigma}=(\hat{\sigma}_{1},\hat{\sigma}_{2},\hat{\sigma}_{3})$ is a vector of Pauli spin matrices, and $\check{\Delta}=\Delta\hat{\tau}_{1}$ is the superconducting pair potential. The self-energy $\check{\Sigma}=\check{\Sigma}_{so}+\check{\Sigma}_{sf}$, where $\check{\Sigma}_{so}=\bm{\sigma}\cdot\check{g}\bm{\sigma}/(8\tau_{so})$ and $\check{\Sigma}_{sf}=\bm{\sigma}\cdot\hat{\tau}_{3}\check{g}\hat{\tau}_{3}\bm{\sigma}/(8\tau_{sf})$ describe spin and charge imbalance relaxation due to the spin-orbit scattering and exchange interaction with magnetic impurities with corresponding relaxation times $\tau_{so}$ and $\tau_{sf}$, and the Pauli matrices $\hat{\tau}_{j}(\hat{\sigma}_{j})$ are in the Nambu (spin) space. We choose the Nambu spinor as

[TABLE]

where $\top$ denotes a transpose.

With the rotation angle described in Eq. (1), the exchange field can be written as

[TABLE]

where $h$ is the exchange field strength. With this choice, the exchange field depends only on $x$, and rotates in the $xz$-plane, namely, $h_{y}=0$ everywhere. Correspondingly, the gradient in Eq. (2) becomes a derivative in the $x$-direction.

The Usadel equation we apply here is based on a lowest-order spherical harmonics expansion of the Green’s function in terms of the momentum direction. Hence it cannot describe the domain wall superconductivity in a d-wave superconductor/ferromagnetic insulator multilayers. Komori et al. (2018) The physics studied in Ref. Komori et al., 2018 is hence outside the scope of the present work.

II.2 Rotation matrix

By introducing a position dependent rotation matrix, we can rotate the spin axis parallel to the local magnetization direction, so that the inhomogeneous exchange field in Eq. (2) can be treated as homogeneous. We define the rotation matrix as

[TABLE]

where $\alpha(x)$ is the rotation angle in Eq. (1), and $\hat{\sigma}_{y}$ is the second Pauli matrix. Considering this rotation matrix, we define a new quasiclassical Green’s function $\check{g}_{0}$

[TABLE]

so that $\check{g}_{0}$ satisfies

[TABLE]

where

[TABLE]

[TABLE]

Thus the problem reduces to solving the Green’s function $\check{g}_{0}$ for a homogeneous exchange field $h$, but with a redefined gradient with an SU(2) type vector potential $A$. Moreover, it is also straightforward to show that $\check{g}_{0}$ satisfies the normalization condition $\check{g}_{0}^{2}=\check{g}^{2}=\check{1}$.

The long derivative in Eq. (4) with the form in Eq. (5) introduces some extra terms with respect to the usual derivative and commutator terms in the Usadel equation. One of them is of the form

[TABLE]

This term has a similar form as the spin-orbit relaxation with a relaxation rate $\alpha^{\prime}(x)^{2}D/2$, but only in one spin direction. It is hence similar to the intrinsic (Rashba or Dresselhaus) spin-orbit coupling.Bergeret and Tokatly (2013, 2014); Gorini et al. (2010) Since $\alpha^{\prime}(x)\neq 0$ inside the domain wall, this term is a direct outcome of the existence of the inhomogeneous exchange field. In other words, inhomogeneous exchange field acts like a spin-orbit relaxation, reducing the effect of the exchange field without affecting the superconducting pair potential. This is shown on some physical quantities in the following sections.

II.3 Boundary condition

The long derivative $\partial^{A}_{x}g$ in the Usadel equation in Eq. (4) has to be continuous, so that a discontinuity in the derivative of $\alpha(x)$ implies a discontinuity in the derivative of $\check{g}$ at $x=\pm\lambda/2$. In order to describe this discontinuity, we integrate Eq. (4) at the boundary, obtaining

[TABLE]

where $x_{b}=\pm\lambda/2$ and $\pm$ refers to the right and left side of the boundary. These boundary conditions together with the solutions in the case of a homogeneous exchange field far from the domain wall form the boundary conditions to the solutions of the Usadel equation in Eq. (4).

II.4 Parameterization

The quasiclassical Green’s function $\check{g}$ and $\check{g}_{0}$ are $4\times 4$ matrices in the Nambu $\otimes$ spin space. Since $\check{g}$ represents the Green’s function for the inhomogeneous exchange field, we parameterize $\check{g}_{0}$ following the parameterization of the quasiclassical Green’s function for inhomogeneous exchange field in Ref. Ivanov and Fominov, 2006,

[TABLE]

The advantage of using this parameterization is that $\theta$ and $M_{0}$ are real scalars and $\bm{M}=(M_{1},M_{2},M_{3})$ is a real vector in the Matsubara representation. The normalization condition $\check{g}_{0}^{2}=1$ adds the constraint

[TABLE]

With the parameterization in Eq. (8), we get a set of differential equations from Eq. (4)

[TABLE]

[TABLE]

We can see directly from the vector differential equation (11) that the parameter $M_{2}$ does not have a source term from the domain wall structure, and therefore $M_{2}=0$ everywhere. Furthermore, one cannot directly solve these equations without the help of the constraint in Eq. (9). Therefore, it is more convenient to use this constraint to transform Eq. (11) to another set of differential equations for each of the component of $M_{i}$ as in Ref. Baker et al., 2016.

Taking the second derivative of the constraint in Eq. (9), we obtain

[TABLE]

Substituting Eq. (12) to each component of Eq. (11), yields

[TABLE]

where $j=0,1,3$, and

[TABLE]

The differential equations in (10) and (13) have to be supplemented by the boundary conditions in Eq. (7) and the solutions at the regions far from the domain wall. With the parameterization in Eq. (8), Eq. (7) becomes

[TABLE]

where $M_{i^{\prime}}=M_{3}$ for $i=1$, and $M_{i^{\prime}}=-M_{1}$ for $i=3$.

II.5 Solutions of the Usadel equation

Without the exchange field $h=0$, the solutions of the differential equations (10) and (13) give the regular results of the $\theta$ parameterization with $M_{0}=1$ and $M_{i}=0$. In the case of a homogeneous exchange field $\bm{h}=h\hat{z}$, the values of $\theta$ and $M_{0}$ are changed, and $M_{3}\neq 0$, due to the odd-frequency spin triplet superconductivity (with zero spin projection on the spin $z$ axis) is induced. In the absence of spin relaxation terms

[TABLE]

These results can be used to describe the solution of the differential equations far from the domain wall structure. For the inhomogeneous exchange field with a domain wall structure, however, the differential equations (10) and (13) cannot be solved analytically, but the numerical solutions are plotted in Fig. 2.

We can see that the domain wall structure brings certain changes to the homogeneous solutions. For $\theta$, $M_{0}$, and $M_{3}$, the effect of the exchange field around the domain wall is reduced. This introduces a Gaussian function like structure for the solutions of these parameters as a function of position. At the center of the domain wall, the effect of the exchange field is reduced the most and the values of these parameters reach their homogeneous solutions gradually away from the domain wall. For a smaller domain wall, the reduction of the effect of the exchange field is more obvious, until for $\lambda\rightarrow 0$, the effect of the exchange field is completely lifted at the center of the domain wall. This is due to the existence of the spin-orbit relaxation like term in Eq. (6) in the rotated Usadel equation with a relaxation rate $\alpha^{\prime}(x)^{2}D/2$ and $\alpha^{\prime}(x)\propto\lambda^{-1}$.

The domain wall structure also introduces a nonzero solution of $M_{1}$ around the domain wall, due to the odd-frequency spin triplet superconductivity (with nonzero spin projection on the spin $z$ axis) induced by the inhomogeneous exchange field. The maximum of $M_{1}$ appears at the boundary of the domain wall and gradually goes to zero away from it. It also changes sign at the two sides of the domain wall center. The amplitude of $M_{1}$ first increases and then decreases with increasing $\lambda$, the maximum taking place at $\lambda\approx 2\xi_{0}$.

Spin relaxation also brings many changes to the solutions of the differential equations. For the regions far from the domain wall, the consequence of spin relaxation is the same as with the case of homogeneous exchange field. The spin-orbit relaxation reduces the effect of the exchange field, therefore, the solutions in these regions approach the ones in the domain wall center. The spin-flip relaxation reduces the superconducting pair potential $\Delta$, and therefore $\theta$ in these regions becomes smaller, and $M_{0/3}$ becomes larger with spin-flip relaxation. These can be seen from Fig. 3(a,b,d).

For the regions around the domain wall, the spin relaxation brings minor changes. This is due to the fact that at the center of the domain wall, the effect of the exchange field is already reduced, and spin relaxation affects superconductivity similarly to the case without exchange field. This can be seen in Fig. 3(b,d). These behaviors can also be revealed in the physical quantities as discussed below.

With these solutions, and using Eq. (3), the unrotated Green’s function is

[TABLE]

In the rest of the paper, we use this Green’s function to calculate various physical quantities.

II.6 Self-consistent calculations

The superconducting pair potential in Eq. (4) has to be determined self-consistently. In the matsubara representation, we have

[TABLE]

where $\check{g}$ is given by Eq. (15), $\gamma$ is the coupling constant, and $\omega_{D}$ is the BCS cutoff energy. The latter gives a temperature dependent cutoff $N_{D}=\omega_{D}/(2\pi T)$ to the sum over $n$. Considering the relations

[TABLE]

we can rewrite the self-consistency equation as

[TABLE]

We can see that $\Delta$ does not explicitly depend on the rotation angle $\alpha(x)$. Its position dependence comes from the parameters $M_{0}$ and $\theta$, whose solutions depend on $\alpha(x)$. This is because we only consider spin-singlet pairing in the self-consistent calculations.

The self-consistent pair potential in Eq. (18) is position dependent, $\Delta=\Delta(x)$. At the regions far from the domain wall, $\Delta(x=\pm\infty)$ is the same as in the case of homogeneous exchange field. A homogeneous exchange field brings many interesting effects on $\Delta$. It suppresses $\Delta$ at a finite temperature and for $T\rightarrow 0$, $\Delta\rightarrow 0$ for a field $h>h_{c}$. Here $h_{c}=\Delta_{0}/\sqrt{2}$ when $T=0$, called the Chandrasekhar-Clogston limitChandrasekhar (1962); Clogston (1962). The superconducting pair potential is reduced by spin relaxation, and $h_{c}$ is enhanced by the spin-orbit relaxation and suppressed by the spin-flip relaxation.Bergeret et al. (2018, 2017)

In the case of an inhomogeneous exchange field, above properties are similar. However, in the domain wall region, the weak effect of the spin relaxation (Fig. 3) brings less changes to $\Delta$, compared to the homogeneous case.

III Equilibrium Spin Current Density

One important consequence of the inhomogeneous exchange field is the equilibrium spin current density. Due to the inhomogenetity of the magnetization, the spin of the quasiparticles rotates following the local magnetization, which creates a flow of spinRalph and Stiles (2008).

In the quasiclassical theory, the spin current can be calculated from

[TABLE]

where $\sigma_{N}=2e^{2}N_{0}D$ is the normal state conductivity and $N_{0}$ is the density of states at the Fermi level. The spin current density $j_{k,i}$ is a tensor, the index $k$ represents the transport direction of spin, and $i$ represents the spin component. With the parameterization in Eq. (15), we get the three components of the spin current density in the $x$ direction,

[TABLE]

[TABLE]

where $\alpha(x)$ is the rotation angle in Eq. (1). Since the rotation of the magnetization happens in the $xz$ plane, the $x$ and $z$ spin components of the spin current density are absent, $j_{x,1/3}=0$. The only nonzero spin component of the spin current density is in the $y$ spin direction.

The spatial dependence of $j_{x,2}$ is shown in Fig. 4(a) for a domain wall size $\lambda=0.5\xi_{0}$. We can see that a nonzero spin current density is created around the domain wall structure. Inside the domain wall, the maximum spin current density $j_{x,2}(0)$ is constant, and smoothly goes to zero outside the domain wall. Both spin-orbit and spin-flip relaxations reduce the spin current density.

The constant spin current density $j_{x,2}(0)$ in the domain wall region is determined from the adiabatic solution of the parameter $M_{3}$ and the domain wall size $\lambda$. From Fig. 2(d) we know that $M_{3}$ describes the effect of the exchange field on superconductivity, and at the center of the domain wall $M_{3}\rightarrow 0$ for a small domain wall size. Therefore, $j_{x,2}(0)\rightarrow 0$ for $\lambda\rightarrow 0$ and becomes large for a stronger exchange field. These are shown in Fig. 4(b,c). We can also see that for a small domain wall, both spin-orbit and spin-flip relaxation reduce $j_{x,2}(0)$. However, for large domain walls, spin-flip relaxation increases $j_{x,2}(0)$ for a fixed domain wall size, but it kills $j_{x,2}(0)$ at the Chandrasekhar-Clogston limit, where the superconductivity is suppressedBergeret et al. (2018, 2017).

The temperature dependence of $j_{x,2}(0)$ is shown in Fig. 4(d). It resembles the temperature dependence of the self-consistent pair potential $\Delta$. This means $j_{x,2}(0)\sim\Delta$ for a fixed domain wall size and exchange field $h$. The effect of spin relaxation also resembles the temperature dependence of $\Delta$.

From the spatial dependence of $j_{x,2}$ we can see that, strictly speaking, the spin current density is not conserved. The position dependence is related with the appearance of a spin-transfer torqueRalph and Stiles (2008). It is exerted on the spins in order to reorient spin flow to follow the direction of local magnetization, namely, it represents the rotation of spins through the domain wall structure. Then in the continuity equation for spin density $\bm{n}_{s}$, we have Sun and Xie (2005)

[TABLE]

where $\mathcal{T}$ is the spin-transfer torque and it is nonzero in the case of an inhomogeneous exchange field and spin-relaxation.

In the static case $\nabla\cdot j_{k,i}=\mathcal{T}$. With the parameterization in Eq. (15),

[TABLE]

It is nothing but the Matsubara sum of $M_{1}$ but constrained by superconductivity ($\sin\theta$). In the normal state $\theta\rightarrow 0$, then the torque is zero in equilibrium.

The expression of $\mathcal{T}$ in Eq. (20) also helps us to understand the properties of the spin current density. The spatial dependence follows Fig. 2(c), but the amplitude is constrained by $\sin\theta\sim\Delta$. That is why the temperature dependence of the spin current density is similar to that of the self-consistent pair potential (Fig. 4(d)). The dependence on $h$ explains the monotonous dependence on $h$ of the spin current density (Fig. 4(c)).

The equilibrium spin-transfer torque $\mathcal{T}$ is related to the superconducting free energy by

[TABLE]

where $V$ is the volume of the superconductor, $\bm{h}$ is the exchange field induced in the superconductor, and $F_{sn}$ is given by

[TABLE]

where $W$ is the cross sectional area of the superconductor and the superconducting free energy densityEilenberger (1968); Usadel (1970); Altland et al. (1998); Taras-Semchuk and Altland (2001); Virtanen et al. (2017) compared to its normal state $f_{sn}$ is given by (see Appendix A)

[TABLE]

In the absence of spin relaxation and exchange field, this agrees with the result in Ref. Virtanen et al., 2017. Here $f_{sn}=f_{sn}[\bm{h}]$, and $\Delta$ and $\check{g}$ are the self-consistent values of the order parameter and the Green function. With the parameterization in Eq. (15), it becomes

[TABLE]

From the relation in Eqs. (19,21), we obtain

[TABLE]

where $f_{sn}^{0}$ is the free-energy density in Eq. (24) with terms that do not directly depend on $\alpha^{\prime}(x)$, and $f_{sn}-f_{sn}^{0}$ is nonzero only in the domain wall region. We can see that the spin current density contributes to the energetics of the system, which in turn influences the formation of the domain wall.

In nonequilibrium spin transport, spin-transfer torque leads to the domain wall motion and influences the orientation of the magnetizationRalph and Stiles (2008). In our model, the equilibrium spin-transfer torque in Eq. (20) does not make the domain wall move, but it contributes to the superconducting free energy via the spin current density (as in Eq. (25)), which in turn affects the domain wall size, as we discuss in the next section.

IV Effect of superconductivity on the domain wall size

In the absence of the external magnetic field, the domain wall size in a ferromagnet is dictated by the competition between exchange and anisotropy energiesLandau et al. (2013). Exchange energy tries to maintain the direction of the magnetization, while the anisotropy energy tends to align the magnetization to a particular direction.

The exchange energy density can be expressed by

[TABLE]

where $m_{i}$ is the component of the magnetization unit vector, and $Q$ is the exchange stiffness constant. Here $\bm{m}=(\cos\alpha,0,\sin\alpha)$, where $\alpha$ is the rotation angle in Eq. (1). Substituting Eq. (1) to $U_{ex}$ we get

[TABLE]

The anisotropy energy density depends on the crystal structure of the system. Most of the ferromagnetic insulators have face centered cubic crystal structureYang et al. (2014). The anisotropy energy density in this case is given by

[TABLE]

where $K_{c1},K_{c2}<0$ for many ferromagnetic insulators, which makes the magnetization lie in one of the easy planes. In the case of thin films, the symmetry is broken in the direction perpendicular to the film plane and this energy density can be expressed by a uniaxial crystal structure asLi et al. (2015)

[TABLE]

where $K_{1}=K_{c1}$, $K_{2}=-7K_{c1}/8+K_{c2}/8$, and $K_{s}$ is the surface anisotropy constant representing the rotation of the easy plane towards an easy axis magnetization with film thickness $t$.

Together with Eqs. (26,27), the domain wall energy is given by

[TABLE]

where $K_{\textrm{eff}}=K_{1}+K_{s}/t+3K_{2}/4$. Minimization of this energy with respect to $\lambda$ gives the domain wall size of the inhomogeneous exchange field in the ferromagnet

[TABLE]

Then the minimized domain wall energy is

[TABLE]

In many studies, the domain wall structure is represented by the hyperbolic functions with Landau et al. (2013)

[TABLE]

where $\lambda^{\prime}=\sqrt{Q/K_{\textrm{eff}}^{\prime}}$ and $K_{\textrm{eff}}^{\prime}=K_{1}+K_{s}/t+2K_{2}/3$. The minimized energy is then given by $4\sqrt{K_{\textrm{eff}}^{\prime}Q}$. Since $K^{\prime}_{\textrm{eff}}\approx K_{\text{eff}}$, this is very close to $F_{1}^{\textrm{min}}$ in Eq. (30). In other words, the domain wall structure in Eq. (1) gives the minimum energy which is only slightly larger with that of more complicated domain wall structures.

The contribution of superconductivity to the domain wall energy is given by the difference of the free energy in the cases of inhomogeneous and homogeneous magnetization

[TABLE]

where $F_{sn}$ is given by Eq. (22). This energy cannot be expressed analytically, but the numerical result is plotted in Fig. 5(a). We can see that $F_{2}$ is negative. This is due to the fact that the existence of the domain wall structure enhances the superconducting condensation energy near the domain wall. This contribution is stronger for a smaller domain wall, since smaller domain wall reduces the effect of the exchange field more in the domain wall region. The effect of spin relaxation is plotted in Fig. 5(b). We can see that spin-orbit relaxation makes $|F_{2}|$ smaller. This is because spin-orbit relaxation reduces the effect of the exchange field, but it has less effect in the domain wall region. This makes the free energy difference $|F_{sn}(\bm{h})-F_{sn}(\bm{h}=h\hat{z})|$ smaller. Spin-flip relaxation on the other hand, makes $|F_{2}|$ larger.

From Fig. 5(a) we can approximate $F_{2}$ as

[TABLE]

where $a,b$ are functions of temperature and exchange field strength and can be determined numerically. For example, for $h=0.3\Delta_{0}$, $a=0.4N_{0}\Delta_{0}^{2}W\xi_{0}^{2}$ and $b=0.085N_{0}\Delta_{0}^{2}W\xi_{0}$.

Minimizing the total energy $F_{1}+F_{2}$ with respect to $\lambda$, we get a compact analytical expression for two extreme limits

(i) $\lambda\ll\xi_{0}$

[TABLE]

(ii) $\lambda\gg\xi_{0}$

[TABLE]

[TABLE]

In both limits superconductivity reduces the domain wall size. (i) For a small domain wall $\lambda\ll\xi_{0}$, superconductivity reduces the domain wall size by effectively increasing the anisotropy constant $K_{\textrm{eff}}$. (From the numerical results $a/(2\xi_{0}^{2})>2b/\xi_{0}$ holds for all $h$.) (ii) For the case of a large domain wall, the situation is more complicated. For a ferromagnet with strong stiffness $Q\pi^{2}>a$, superconductivity reduces the domain wall size by effectively reducing the exchange stiffness constant $Q$. For a ferromagnet with weak stiffness $Q\pi^{2}<a$, which also refers to the case of $K_{\textrm{eff}}\ll 1.0N_{0}\Delta_{0}^{2}W$, the domain wall size is also reduced. However, in this case superconductivity dominates the domain wall energy, and leads to a negative total energy $F_{1}+F_{2}$, which introduces a dense domain structure with a domain length comparable with the domain wall size.Buzdin (2005) Since we are considering a single domain structure, here we only consider domain walls with positive domain wall energy.

In Fig. 5(c), we numerically minimize the domain wall energy $F=F_{1}+F_{2}$ with respect to the domain wall size for the case of $h=0.3\Delta_{0}$ and $T=0.05\Delta_{0}$, and calculate the relative change of the domain wall size. To avoid a negative domain wall energy leading to a transition to a system with many domain walls, we set in each figure a constant $F_{1}\geq 0.12N_{0}\Delta_{0}^{2}W\xi_{0}$ so that $F>0$ for each case considered. In other words, instead of varying $K_{\rm eff}$ and $Q$, we fix $F_{1}$ and vary $\lambda_{0}$ in the figures. The effect of superconductivity on the domain wall size is strongest for the lowest $F_{1}$. If we consider larger values of $F_{1}$, namely, a stronger ferromagnet, the effect of superconductivity on the domain wall size is smaller.

The effect of spin relaxation on the domain wall size is shown in Fig. 5(d). We see that for small domain walls, the spin relaxation brings very little effects on the domain wall size compared to the case without spin-relaxation. However, for larger domain walls, the two types of spin-relaxation mechanisms lead to different effects on the domain wall size. Spin-orbit relaxation makes the effect smaller since it reduces the effect of the exchange field and makes $|F_{2}|$ smaller. Spin-flip relaxation makes the effect stronger, since $|F_{2}|$ is larger.

As we can see, the effect of superconductivity on the domain wall size is pronounced for weak ferromagnets with large domain walls. Domain wall sizes in ferromagnetic insulators are rarely reported, but with Eq. (29) we can evaluate a typical size. For materials with a face centered cubic structure, the exchange stiffness constant $Q$ is related to the exchange integral $J$ by $Q=4S^{2}J/a_{0}$, where $S$ is spin and $a_{0}$ is the lattice constant. Chikazumi and Graham (2009) With the values of the parameters for $Q$ in Ref. Söllinger et al., 2010 and for $K_{\rm{eff}}$ in Ref. Li et al., 2015, we get the domain wall size for a EuS thin film as $50\ \text{nm}$. Considering a typical coherence length of a conventional superconductor in the diffusive superconductors at zero temperature and exchange field between 15 nm and a few hundred nm, we can estimate in this case $\lambda_{0}\sim 0.1\dots 5\xi_{0}$. Therefore, although perhaps in most cases $\lambda_{0}<\xi_{0}$ is the most relevant limit, nothing as such seems to exclude the possibility of the opposite limit as well.

V Density of states

The inhomogeneous exchange field also makes the local density of states (DoS) peculiar. Since the domain wall structure reduces the effect of the exchange field, the local DoS in the domain wall region is different from the one for homogeneous magnetization.

In the quasiclassical theory, the local DoS for each spin species is given by

[TABLE]

where $\pm$ is for spin $\uparrow/\downarrow$. With the parameterization in Eq. (15), it becomes

[TABLE]

In Fig. 6(a) we show the local DoS $N_{+}=N_{\uparrow}+N_{\downarrow}$ at the center of the domain wall for different domain wall sizes in the absence of spin relaxation. At the center of the domain wall $\alpha=\pi/2$ and $M_{1}=0$, then $N_{\uparrow}=N_{\downarrow}$ and $N_{-}=N_{\uparrow}-N_{\downarrow}=0$ due to the symmetry of the model. For $\lambda\rightarrow\infty$, $N_{+}$ has two BCS peaks at $\Delta\pm h$. As the domain wall becomes smaller, the inner peak is shifted towards the outer one and forms a shark-fin structure for a very small domain wall size [$\lambda=0.1\xi_{0}$ in Fig.6(a)]. In this case, although $N_{-}=0$, $N_{+}$ is different from the case with zero exchange field. Due to the spin-orbit relaxation like effect of the inhomogeneous exchange field, the peaks are wider and the superconducting gap is smaller compared to the case of $h=0$.

Since the magnetization direction is opposite on the two sides of the domain wall, $N_{s}$ behaves differently in these two regions, such that $N_{+}(x)=N_{+}(-x)$ and $N_{-}(x)=-N_{-}(-x)$. In Fig. 6(b) we show $N_{\uparrow}$ at different location along the superconducting wire with a small domain wall $\lambda=0.1\xi_{0}$, in the absence of spin relaxation. These results in are similar with those in Ref. Strambini et al., 2017, which concentrates on the limit $\lambda\ll\xi_{0}$. Here we further study the effect of spin relaxation on the local DoS in Fig. 6(c,d). The DoS in the presence of the two kinds of spin relaxation is plotted in Fig. 6(d) at the center of a domain wall with a size $\lambda=0.1\xi_{0}$. For comparison, the DoS at $x=-\infty$, which also refers to the homogeneous exchange field is plotted for the same parameters in Fig. 6(c). In both cases, spin-orbit relaxation broadens the peaks, but keeps $\Delta(x)$ unchanged. Spin-flip relaxation also broadens the peaks, but it suppresses $\Delta(x)$.

These properties of the local DoS are caused by the inhomogeneous magnetization, and are needed for understanding the tunneling conductance as discussed in the next section.

VI Tunneling conductance

The local density of states is visible in measurements of a tunneling conductance through a barrier in contact with the FI/S bilayer. However, the results depend on whether the barrier itself is magnetic or not. Therefore, we consider either a non-magnetic tunnel contact (NISFI) or tunneling through the FI containing the domain wall (SFIN). In either system the tunneling current can be written as

[TABLE]

where $G_{T}$ is the normal-state conductance, $V$ is the applied voltage and $f=\left[1+\exp(\epsilon/k_{B}T)\right]^{-1}$ is the Fermi function. Here the averaged density of states $\bar{N}(\epsilon)$ over the tunneling area is given byBergeret et al. (2012)

[TABLE]

where $P(x)$ is the spin polarization of the junction, $L=x_{b}-x_{a}$ is the size of the tunneling barrier, and $N_{s}$ is given in Eq. (33). Using the definition of $N_{\pm}=N_{\uparrow}\pm N_{\downarrow}$, we can also write

[TABLE]

The tunneling conductance is given by

[TABLE]

In the case of a scanning tunneling microscope measurement with a small normal metal tip ($x_{b}-x_{a}\ll\xi_{0}$ and $P(x)=0$), $\bar{N}=N_{+}$. In this case the tunneling conductance at $T\rightarrow 0$ gives the local DoS $N_{+}$ shown in Fig. 6.

For a planar tunneling spectroscopy measurement with a normal metal electrode having size $x_{b}-x_{a}\geq\xi_{0}$, $\bar{N}$ is the averaged $N_{+}$ since $P(x)=0$. The tunneling conductance in this case is given in Fig. 7(c). We show the comparison of different sizes of tunneling barriers (centered at $x=0$). Since we are considering a single domain wall structure, a large tunneling barrier produces tunneling conductance identical with the case of homogeneous magnetization (green curve). If we choose a smaller tunneling barrier, the effect of the domain wall is more obvious. It reduces the effect of the exchange field, becomes similar to the case for a small homogeneous exchange field strength $h$, and reshapes the peaks in the tunneling conductance (red and blue curves).

If the tunneling is through the ferromagnetic insulator, we need to include the effect of finite spin polarization. Assuming the polarization of the transmission to follow the local magnetization, the position dependent polarization $P$ is related to the rotation angle of the magnetization as $P(x)=P_{0}\cos\alpha$. The tunneling conductance is plotted in Fig. 7(d). The asymmetric dependence in the injection voltage in the conductance is a direct result of the polarization of the ferromagnet. Similar to Fig.7(c), for a large tunneling barrier the tunneling conductance is the same for the case of a homogeneous exchange field $h$ with a finite spin polarization $P_{0}$. The effect of the domain wall is again obvious for smaller tunneling barriers, and the spin-orbit relaxation like effect of the domain wall structure reduces the effect of the exchange field.

VII Conclusion

In conclusion, we have studied various properties of a superconductor in contact to a ferromagnet with a domain wall. We have studied the equilibrium spin current density, which exists due to the inhomogeneity of the exchange field. It also contributes to the superconducting free energy, which in turn affects the domain wall size. We show that the domain wall size is reduced by the contribution of the superconducting free energy to the domain wall energy. We have also studied the peculiar density of states around the domain wall, and the tunneling conductance. Our work can be a precursor to the study of the nonequilibrium effects, in particular domain wall motion.

Acknowledgements.

We thank F. S. Bergeret, M. Silaev and S. Chakraborty for discussions. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 800923 (SUPERTED) and the Academy of Finland projects No. 305256 and 317118.

Appendix A Expression for the free energy

For completeness, we give now a brief argument for the form of the free energy in the superconducting state. There are several ways to arrive at such results, and below we follow a similar procedure as in Ref. Eilenberger, 1968. The free energy density difference between superconducting and normal states can be written as (we set $N(0)=1$ in the following) Maki and Tsuneto (1964); Gor’kov (1959)

[TABLE]

The integral could now in principle be computed numerically from Usadel equations, solving $\Delta$ self-consistently as a function of the coupling constant $\gamma$ in the geometry considered. This is cumbersome, and also unnecessary, as the integral can be evaluated analytically as follows.

Suppose there exists a functional $R=\int\mathrm{d}^{3}x\,\,r$,

[TABLE]

whose saddle point vs. $\vec{z}$ and $\Delta$ defines $\Delta(q)$. In other words, the self-consistency (16) and Usadel equations (2) are defined by variations

[TABLE]

and $\vec{z}$ is some parameterization of $\hat{g}$ under the constraint $\hat{g}^{2}=1$ of the Usadel equation. Now,

[TABLE]

where the last two terms vanish due to the saddle point conditions. Therefore,

[TABLE]

This assumes there is a continuous solution branch connecting the normal and the superconducting states as the coupling constant $\gamma$ (ie. $q$) is changed.

The next step is to find a suitable $R$ that satisfies Eqs. (38). Its form can be guessed (or derived) based on $\sigma$-model results Belitz and Kirkpatrick (1994); Altland et al. (1998); Taras-Semchuk and Altland (2001):

[TABLE]

First, we can note that variation vs $\Delta^{*}$ (or $\Delta$) gives the self-consistency relation (16). For variations with $\hat{g}$, we parametrize $\delta\hat{g}=e^{\delta{}\hat{W}}\hat{g}e^{-\delta{}\hat{W}}-\hat{g}\simeq\delta W\hat{g}-\hat{g}\delta\hat{W}$ to retain the normalization condition. Requiring variation vs. $\hat{W}$ to vanish we find Eq. (2). Hence, the result has the property (38).

The sum defining $R$ is not convergent, and only the difference in Eq. (40) is well defined. There’s also an implicit Matsubara frequency cutoff in the term that appears in the self-consistency equation. We can eliminate this issue by substituting the self-consistency equation back into the $|\Delta|^{2}$ term in Eq. (41). This results to Eq. (23) in the main text. Variation vs $\mathbf{h}$ still gives the correct quasiclassical expression for the magnetization, but $\Delta$ and $\hat{g}$ can no longer be varied and have to be taken at their saddle-point values.