Spectral Properties and Quantum Phase Transitions in Superconducting Junctions with a Ferromagnetic Link

TL;DR

This paper theoretically investigates the spectral and transport properties of superconducting junctions with magnetic links, revealing quantum phase transitions characterized by changes in total spin and supercurrent sign.

Contribution

It derives a general equation for the subgap spectrum in magnetic Josephson junctions and demonstrates the occurrence of quantum phase transitions in finite-length systems.

Findings

Quantum phase transition occurs when bound states cross the Fermi level.

Total spin changes discontinuously by integer jumps during the transition.

Phase transitions are reflected in the sign change of the supercurrent.

Abstract

We study theoretically the spectral and transport properties of a superconducting wire with a magnetic defect. We start by modelling the system as a one dimensional magnetic Josephson junction and derive the equation determining the full subgap spectrum in terms of the normal-state transfer matrix for arbitrary length and exchange field of the magnetic region. We demonstrate that the quantum phase transition predicted for a short-range magnetic impurity, and associated with a change of the total spin of the system, also occurs in junctions of finite length. Specifically, we find that the total spin changes discontinuously by integer jumps when bounds states cross the Fermi level. The spin can be calculated by using a generalization of Friedel sum rule for the superconducting state, which we also derive. With these tools, we analyze the subgap spectrum of a junction with the length of the magnetic region smaller than the superconducting coherence length and demonstrate how phase transitions also manifest as change of the sign of the supercurrent.

Figures

b). Such coupling lifts the degeneracy at $\varphi=0$ and "pushes" one of the states to energies closer to the continuum. By solving Eq ([6](#S2.E6)) in this limit for a general value of $\varphi$, one obtains four bound states Costa *et al.* ([2018](#bib.bib15)):

Equations

$$\hat{H}(x)=\left(\begin{array}[]{cc}-\frac{\hbar^{2}\partial_{x}^{2}}{2m}-\mu-V(x)&\Delta(x)\\ \Delta^{*}(x)&\frac{\hbar^{2}\partial_{x}^{2}}{2m}+\mu+\bar{V}(x)\end{array}\right)\;.$$

$$\begin{array}[]{ll}\Psi^{\sigma}_{L}(x<\tfrac{-L}{2})=&e^{x/\xi}\Bigg{[}\Bigg{(}\begin{array}[]{c}A_{L}^{\sigma}\\ A_{L}^{\sigma}e^{i\alpha}e^{-i\varphi/2}\end{array}\Bigg{)}e^{ik_{F}x}\\[10.00002pt] &+\left(\begin{array}[]{c}B_{L}^{\sigma}\\ B_{L}^{\sigma}e^{-i\alpha}e^{-i\varphi/2}\end{array}\right)e^{-ik_{F}x}\Bigg{]},\end{array}$$

$$\begin{array}[]{ll}\Psi^{\sigma}_{R}(x>\tfrac{L}{2})=&e^{-x/\xi}\Bigg{[}\Bigg{(}\begin{array}[]{c}A_{R}^{\sigma}\\ A_{R}^{\sigma}e^{-i\alpha}e^{i\varphi/2}\end{array}\Bigg{)}e^{ik_{F}x}\\[10.00002pt] &+\left(\begin{array}[]{c}B_{R}^{\sigma}\\ B_{R}^{\sigma}e^{i\alpha}e^{i\varphi/2}\end{array}\right)e^{-ik_{F}x}\Bigg{]},\end{array}$$

$${\text{C}}_R=\check{T}{\text{C}}_L$$

$${\text{C}}_R=e^{-i\phi }{e}^{i\hat{\alpha }}\check{\overline{T}}{e}^{i\check{\alpha }}{\text{C}}_L$$

$$\det\bigg{(}e^{i\varphi}-\check{T}^{-1}e^{i\hat{\alpha}}\check{\overline{T}}e^{i\hat{\alpha}}\bigg{)}=0\;.$$

$${\hat{G}}^R(ϵ)=(ϵ-{\hat{H}}_0-\hat{V}+i{0}^{+}{)}^{-1},$$

$$S_{i}=-\frac{\hbar}{4\pi}\int_{-\infty}^{\infty}d\epsilon f_{F}(\epsilon)\text{Im}\bigg{[}\text{Tr}\Big{\{}{\hat{\sigma}}_{i}\hat{G}^{R}(\epsilon)\Big{\}}\bigg{]},$$

$${\hat{G}}^R={\hat{G}}_0^R+{\hat{G}}_0^R\hat{V}(I-{\hat{G}}_0^R\hat{V}{)}^{-1}{\hat{G}}_0^R.$$

$$S=\frac{\hbar}{2}\int_{-\infty}^{\infty}\frac{d\epsilon}{2\pi}f_{F}(\epsilon)\frac{d}{d\epsilon}\big{[}\delta_{-}(\epsilon)-\delta_{+}(\epsilon)\big{]}\;,$$

$$2S/\hbar=\frac{1}{2\pi}\big{[}\delta_{-}(0)-\delta_{+}(0)\big{]}.$$

$$\check{T}=\left(\begin{array}[]{cc}\hat{T}^{+}&0\\ 0&\hat{T}^{-}\end{array}\right)=\left(\begin{array}[]{cccc}T^{+}_{++}&T^{+}_{+-}&0&0\\ T^{+}_{-+}&T^{+}_{--}&0&0\\ 0&0&T^{-}_{++}&T^{-}_{+-}\\ 0&0&T^{-}_{+-}&T^{-}_{--}\end{array}\right)\;.$$

$$\displaystyle T_{++}^{\sigma}=\bigg{[}\cos\Big{(}q_{\sigma}L\Big{)}+\frac{i}{2}\frac{{q_{\sigma}}^{2}+{q_{0}}^{2}}{q_{\sigma}q_{0}}\sin\Big{(}q_{\sigma}L\Big{)}\bigg{]}e^{-iq_{0}L},$$

$$\displaystyle T_{+-}^{\sigma}=\frac{i}{2}\frac{{q_{\sigma}}^{2}-{q_{0}}^{2}}{q_{\sigma}q_{0}}\sin\Big{(}q_{\sigma}L\Big{)}\;,$$

$$\frac{ϵL}{\hslash v_F}+\sigma \Phi -\arccos \frac{ϵ}{\Delta }\pm \frac{\phi }{2}=\pi n,$$

$$\displaystyle\epsilon=\pm\frac{\Delta}{\Phi^{2}+1}\bigg{[}$$

$$\displaystyle\pm\Phi\sqrt{2\Phi^{2}(1+\cos\varphi)+\sin^{2}\varphi}\ \bigg{]}^{1/2}\;.$$

$$ϵ=\pm \Delta \frac{1-{\Phi }^2}{1+{\Phi }^2}.$$

$$\psi(x)=\left\{\begin{array}[]{lcl}A_{L}^{\sigma}e^{iq_{0}x}+B_{L}^{\sigma}e^{-iq_{0}x}&\text{if}&x<-L/2\\ C^{\sigma}e^{iq_{\sigma}x}+D^{\sigma}e^{-iq_{\sigma}x}&\text{if}&-L/2<x<L/2\\ A_{R}^{\sigma}e^{iq_{0}x}+B_{R}^{\sigma}e^{-iq_{0}x}&\text{if}&x>L/2\end{array}\right.\ ,$$

$$\left(\begin{array}[]{c}A_{R}^{\sigma}\\ B_{R}^{\sigma}\end{array}\right)=\left(\begin{array}[]{cc}T_{11}^{\sigma}&T_{12}^{\sigma}\\ T_{21}^{\sigma}&T_{22}^{\sigma}\end{array}\right)\left(\begin{array}[]{c}A_{L}^{\sigma}\\ B_{L}^{\sigma}\end{array}\right),$$

$$\displaystyle T_{11}^{\sigma}=\bigg{[}\cos\Big{(}q_{\sigma}L\Big{)}+\frac{i}{2}\frac{{q_{\sigma}}^{2}+{q_{0}}^{2}}{q_{\sigma}q_{0}}\sin\Big{(}q_{\sigma}L\Big{)}\bigg{]}e^{-iq_{0}L},$$

$$\displaystyle T_{12}^{\sigma}=\frac{i}{2}\frac{{q_{\sigma}}^{2}-{q_{0}}^{2}}{q_{\sigma}q_{0}}\sin\Big{(}q_{\sigma}L\Big{)}\;,$$

$$\cos\varphi-\text{Re}\bigg{[}T_{12}^{\sigma}\overline{T}_{12}^{\sigma}+e^{2i\alpha}T_{11}^{\sigma^{*}}\;\overline{T}_{11}^{\sigma}\bigg{]}=0,$$

$$2\cos \phi$$

$$-\frac{{q_{\sigma }}^2+k_F^2}{q_{\sigma }{k}_F}\sin (2\alpha )\sin (q_{\sigma }L)\cos (\overline{q_{\sigma }}L)+\frac{{\overline{q_{\sigma }}}^2+k_F^2}{\overline{q_{\sigma }}{k}_F}\sin (2\alpha )\cos (q_{\sigma }L)\sin (\overline{q_{\sigma }}L)=0.$$

Full text

I Introduction

The study of Josephson magnetic junctions and magnetic impurities in superconductors has attracted a great deal of attention in the past years. The research is mainly motivated by the search of a topological superconducting state in magnetic impurity chains and clusters embedded in a superconductor Yazdani et al. (1997); Franke et al. (2011); Nadj-Perge et al. (2014); Heinrich et al. (2018); Choi et al. (2018). In this context it is essential to understand the spectral properties around the magnetic region. In a quasi one-dimensional setup this is equivalent to study the spectrum of a superconductor-ferromagnet-superconductor (SFS) junction.

Ballistic SFS junctions have been widely explored in the past, mainly in two limiting cases. One of them is the semiclassical limit, in which the Fermi energy, $\mu$, is assummed to be much larger than any other energy involved in the system, including the superconducting gap, $\Delta$, and the Zeeman splitting, $h$ Golubov et al. (2004); Buzdin (2005); Bergeret et al. (2005); Konschelle et al. (2016a, b). In this limit, one can directly apply the Bohr-Sommerfeld semiclassical quantization condition Duncan and Györffy (2002) and demonstrate that, in the absence of interface barriers, the spectrum consists of two double-degenerate Andreev bound states with opposite energies. This degeneracy of the bound states reflects the degeneracy of the $\pm k_{F}$ Fermi momentum valleys, which remain uncoupled in the absence of normal reflection, as schematically showed in Fig 1b.

The second widely studied limiting case is when the spin-splitting field is very large, $h\gg\mu$, and concentrated in a region much smaller than $k_{F}^{-1}$ Yu (1965); Shiba (1968); Rusinov (1968); Costa et al. (2018). This has been described as $\delta$-like magnetic impurity that strongly couples both propagation directions to form two nondegenerate bound states within the gap with opposite energies. These states, known as the Yu-Shiba-Rusinov (YSR) states, may cross the Fermi level at a certain strength of the exchange energy. At this crossing, the system undergoes a quantum phase transition (QPT) Sakurai (1970); Balatsky et al. (2006) that has been widely studied within the $\delta$-like impurity model. However, the discussion of whether such a QPT may take place beyond the impurity model is an open question. To address it, one needs to understand how these two known limiting cases are connected.

The goal of this work is twofold. On the one hand, we derive a general equation, Eq (6), that determines the subgap spectrum of a one-dimensional junction in terms of the normal-state transfer matrix for an arbitrary spin-dependent potential describing the F region, assuming that $\Delta=0$ within F. For the particular case of a collinear (unidirectional) magnetization in the F region, we derive a generalized Friedel sum rule, Eq (11), adapted for the superconducting state. This rule states that every time a bound state crosses the Fermi energy, the total electronic spin changes by the amount of $\hbar/2$. Importantly, this sum rule is valid not only for one-dimensional systems, but applies universally to any dimension, size and shape of a localized magnetic region.

On the other hand, in Sec. III, we use these findings to provide a complete analysis of the subgap spectrum of a ballistic one-dimensional SFS junction for and arbitrary homogenoeus exchange field $h$. We focus on the short junction regime, where the ferromagnetic region is shorter than the superconducting coherence length, $\xi$. In this case, the presence of a superconducting gap in the ferromagnet due to the proximity effect has no effect on the subgap spectral properties of the junction, so we set $\Delta=0$ in F. For this system we obtain the normal-state transfer matrix and, from it, we determine all spectral properties of the system from the central expression Eq (6). We recover the well-established limiting cases, i.e. delta-like and semiclassical magnetic region, but also the subgap spectrum for all intermediate situations. We identify the values of $h$ and $L$ at which zero-energy crossings of bound states occur. As in the YSR case, these crossings are associated with a QPT, which manifests as a change of the total electronic spin of the system, in accordance with the sum rule derived in section II.1. We finally demonstrate that this change of the total spin at the QPT is associated with the change of sign of the supercurrent in the SFS junction or, equivalently, to a change of the ground state phase difference between the superconductors from [math] to $\pi$.

II Model and General Properties

We consider a one-dimensional geometry consisting of a superconducting wire interrupted by a ferromagnetic region, as sketched in Fig 1a. The Bogoliubov-De Gennes (BdG) Hamiltonian of the system reads

[TABLE]

Here $\mu$ is the chemical potential, $\Delta(x)$ is the superconducting gap that is only finite on the S electrodes, $\Delta(|x|>L/2)=|\Delta|e^{\pm i\varphi/2}$, the length of the F region is labeled by $L$ and $\pm\varphi/2$ is the superconducting phase, where the plus (minus) sign stands for the right (left) superconductor. The potential $V(x)=V_{0}(x)+\bm{h}(x)\cdot\bm{\sigma}$ is only finite, but arbitrary, within the region $|x|<L/2$ and it consists of an scalar component $V_{0}$ and a spin-dependent one $\bm{h}(x)\cdot\bm{\sigma}$. The "bar" denotes time-reverse conjugation such that $\bar{V}=\hat{\sigma}_{y}V^{*}\hat{\sigma}_{y}$.

We focus on the subgap spectra, $\epsilon<|\Delta|$, which determines the main transport features at zero voltage and low temperatures. For such energies, the decaying wavefunctions into the left (L) and right (R) superconducting leads for each spin component read

[TABLE]

[TABLE]

where the upper (lower) element of the Nambu spinors stand for electrons (holes), the index $\sigma=\pm$ labels components of the spin spinor, $\xi=\hbar v_{F}/\sqrt{\Delta^{2}-\epsilon^{2}}$ is the decaying length of the wavefunction into the superconductor and $k_{F}$ and $v_{F}$ stand for the Fermi wavenumber and the Fermi velocity respectively. The quantity $\alpha$ is the phase associated with each Andreev reflection at the S/F interface and it is given by $\cos\alpha=\frac{\epsilon}{\Delta}$.

The coefficients $A_{L(R)}^{\sigma}$ and $B_{L(R)}^{\sigma}$ in Eqs (2) and (3) are the constants of integration at the left (right) superconductor for the quasiparticles consisting of right moving (those multiplied by $e^{ik_{F}x}$) and left moving (those multiplied by $e^{-ik_{F}x}$) electrons respectively. At this stage it is convenient to define the four vectors $\textbf{C}_{L(R)}\equiv(A_{L(R)}^{+},B_{L(R)}^{+},A_{L(R)}^{-},B_{L(R)}^{-})^{T}$ for the left (right) superconductor.

The wave functions on opposite sides of the F region are connected via the normal state electronic T-matrix, $\check{T}$

[TABLE]

for the electrons and

[TABLE]

for the holes. In Eqs (4) and (5), $e^{i\hat{\alpha}}$ is a diagonal matrix with elements $[e^{i\alpha},e^{i\alpha},e^{-i\alpha},e^{-i\alpha}]$. Notice that time conjugation also implies to change the sign of the quasiparticle energy ($\epsilon\rightarrow-\epsilon$), so that $\check{\overline{T}}(\epsilon)=\hat{\sigma}_{y}\check{T}^{*}(-\epsilon)\hat{\sigma}_{y}$.

After substitution of $\mathbf{C}_{R}$ from Eq (4) into Eq (5) and multiplication by $\check{T}^{-1}$ from the left one obtains a homogeneous equation for $\mathbf{C}_{L}$ that leads to the condition determining the bound states:

[TABLE]

This expression is a generalization of Beenakker’s equation for the Andreev spectrum of a SNS junctions derived from the scattering matrix Beenakker (1991). The second term inside the determinant describes an “Andreev loop”. Namely, from right to left, first an electron from F is Andreev reflected as a hole at one F/S interface. The hole propagates to the opposite interface and it is converted again into an electron via the Andreev reflection. The electron is finally transferred back to the origin. After this cycle, the wavefunction accumulates a phase equal to $\varphi$.

II.1 Sum rule for spin in gapped systems

Before using Eq (6) to calculate the subgap spectrum in a one-dimensional geometry, we can anticipate changes of the total spin of the system associated with bound states crossing the Fermi level. Note that the derivation presented here is valid for any dimension, so that the result that we obtain is not restricted to the one dimensional problem described by Eq (1).

We start by considering the retarded Green’s function (GF) for the BdG equations,

[TABLE]

where $\hat{H}_{0}$ is the unperturbed BdG Hamiltonian of the system and $\hat{V}$ is a general perturbation potential operator, as introduced in Eq (6) . The component $i=\{x,y,z\}$ of the total electronic spin is given by

[TABLE]

where $f_{F}(\epsilon)=(e^{\epsilon/k_{B}T}+1)^{-1}$ is the Fermi distribution function, ${\hat{\sigma}}_{i}$ is the $i$-th Pauli matrix, and the trace runs over the whole coordinate$\times$Nambu$\times$spin space.

The full GF in Eq (8) can be also written in terms of the unperturbed GF, $\hat{G}_{0}$, and the potential $\hat{V}$ via Dyson’s equation, $\hat{G}^{R}=\hat{G}_{0}^{R}+\hat{G}_{0}^{R}\hat{V}\hat{G}^{R}$. Solving it for $\hat{G}^{R}$ and substituting it back into the right hand side, we obtain the expression determining the exact $\hat{G}^{R}$

[TABLE]

As the total spin of the unperturbed system is zero, only the second term of $\hat{G}^{R}$ in Eq (9) contributes to the trace in Eq (8).

Let us now assume that $\hat{V}$ is an energy independent local perturbation, and its magnetic part is collinear with $z$-axis (i.e. it commutes with $\hat{\sigma}_{z}$). Noticing that \big{(}\hat{G}_{0}^{R}\big{)}^{2}=-\frac{d\hat{G}_{0}^{R}}{d\epsilon}, one can use the cyclic property of trace to obtain from Eq (8) the $z$-component of the total spin:

[TABLE]

where \delta_{\sigma}(\epsilon)=\text{Im}\ln\big{[}\det(I-\hat{G}_{0}^{R}(\epsilon)V_{\sigma})\big{]} is a generalized phase shift. Notice that $\hat{V}$ can have any spatial distribution and that the determinant inside the logarithm is the quantization condition coming from the Lippmann-Schwinger equation. In particular, zeros of this determinant determine the spectrum of the bound states. Therefore, in a one-dimensional case, it has to be proportional to the left hand side of Eq (6). At zero temperature ($T=0$), Eq (10) becomes especially simple,

[TABLE]

This result is analogous to the well-known Friedel sum rule that relates the charge/spin induced by a local perturbation to the phase shifts at the Fermi level.

The important feature of the superconducting state is its gap at the Fermi level ($\epsilon=0$), where the unperturbed Green’s function is real (and, therefore, $\det(I-\sigma\hat{G}_{0}^{R}(0)V)$ is real too). Thus, $\delta_{\sigma}(0)/\pi$ can only take integer values, which will only change discontinuously by $\pm 1$ when a spin polarized bound state crosses the middle of the gap, as the determinant changes its sign. The electron-hole symmetry requires that the spin-up/down polarized states cross zero simultaneously while moving in opposite directions. As a result, at every crossing event the normalized spin $2S/\hbar$ jumps by one 111The stepwise process of the spin polarization that follows from our phase-shift arguments agree with the picture of Ref. Balatsky et al. (2006) based of the analysis of the spin structure of the many-body BCS wavefunction in the $\delta$-like impurity case. However, it should be note that the result of this section is valid for any energy independent local perturbation potential $\hat{V}$ acting on a system with a gap at the Fermi energy and Green’s function $\hat{G}_{0}$, as long as $\hat{V}$ commutes both with $\hat{G}_{0}$ and $\hat{\sigma}_{z}$..

In the above derivation we only assume that the perturbation $\hat{V}$ is localized in space and has a collinear magnetic structure. Therefore our sum rule relating the total induced spin to the in-gap spectrum applies to any dimension and any size and shape of a finite magnetic region. For example, it can be directly used to analyze the behavior of the total spin in a magnetic chain on top of a superconductor, as the one studied in Ref Björnson et al., 2017.

III One-dimensional SFS junction

We now apply the results of previous section to compute the spectral properties of a one-dimensional SFS junction. We assume that the scattering F region is described by the potential $\hat{V}(x)=h\hat{\sigma}_{z}$ for $|x|<L/2$. In such a case, the T-matrix in Eq (6) has a block-diagonal structure in spin-space,

[TABLE]

This considerably simplifies the problem, since we only need to calculate the normal state transfer for each spin orientation, $\sigma=\pm$, separately (see Appendix A for details). The elements of $\check{T}$ in Eq (12) read

[TABLE]

where $q_{\sigma}(\epsilon)=k_{F}\sqrt{1+\tfrac{\epsilon+\sigma h}{\mu}}$ and $q_{0}(\epsilon)=k_{F}\sqrt{1+\tfrac{\epsilon}{\mu}}$ are the momentum of the electron in the ferromagnet and the normal metal respectively and $\mu$ stands for the Fermi energy. Due to the symmetry of the problem, one can verify that other components of $\hat{T}^{\sigma}$ are related to the ones in Eq (III) by complex conjugation, T_{--}^{\sigma}=\big{(}T_{++}^{\sigma}\big{)}^{*} and T_{+-}^{\sigma}=\big{(}T_{-+}^{\sigma}\big{)}^{*}. The diagonal terms, ${T}_{aa}^{\sigma}$, describe a direct transmission (forward scattering) within one valley, whereas the off-diagonal terms represent backscattering events that couple the opposite valleys at $\pm k_{F}$ (see Fig 1b).

The solution of Eq (6), after substitution of Eq (III) in it, determines the full subgap spectrum of the SFS junction. For analytic results, we focus on the semiclassical limit where $\mu$ is the largest energy, so that $\epsilon,\Delta,h\ll\mu$. In this case the quassiparticle momenta in the F and S regions are approximated by $q_{\sigma}(\epsilon)\approx k_{F}+\tfrac{\epsilon+\sigma h}{\hbar v_{F}}$ and $q_{0}(\epsilon)\approx k_{F}+\tfrac{\epsilon}{\hbar v_{F}}$, respectively. To the leading order in the semiclassical approximation, the off-diagonal elements of the T-matrix in Eq (III) are negligible and the diagonal terms are given by $T_{++}^{\sigma}\approx e^{\sigma i\Phi}$, where $\Phi\equiv\tfrac{hL}{\hbar v_{F}}$ is referred to as the magnetic phase. This expression for the T-matrix has a simple physical interpretation: within the semiclassical approach the incoming electrons have an energy of the order of $\mu$, much larger than the scattering potential height, $h$. Hence, incoming particles have a unit probability to be transmitted through the F region. Propagation through the F region results only in the additional phase $\Phi$. Clearly, the spectrum obtained from Eq (6) in this limit coincides with the result of the Bohr-Sommerfeld quantization condition for the spectrum:

[TABLE]

where $n$ is an integer. Eq. 14 determines the spectrum of Andreev bound states (ABS) Konschelle et al. (2016b, a). In a short junction, $L\ll\xi_{0}$, where $\xi_{0}\equiv\hbar v_{F}/\Delta$ is the superconducting coherence length, one obtains \epsilon_{\sigma}=\pm\Delta\cos\big{(}\sigma\Phi+\tfrac{\varphi}{2}\big{)}. It follows that by changing the magnetic phase, the energy of the ABS can be tuned between $\pm\Delta$. In particular, zero energy single states can be created by proper choice of $\Phi$ and $\varphi$.

The other widely studied limiting case is the YSR limit in which the F region is described by a $\delta$-like potential, i.e. its length tends to zero, $L\rightarrow 0$, while $\Phi$ is kept finite. One can read directly from Eqs (III) that, within this limit, $T_{++}^{\sigma}\approx 1+\sigma i\Phi$ and the off-diagonal elements are non-zero, $T_{+-}^{\sigma}\approx\sigma i\Phi$. This means that, in the presence of a $\delta$-like potential, the backscattering probability is finite. The latter leads to a coupling between the $\pm k_{F}$ valleys (see sketch in Fig 1b). Such coupling lifts the degeneracy at $\varphi=0$ and "pushes" one of the states to energies closer to the continuum. By solving Eq (6) in this limit for a general value of $\varphi$, one obtains four bound states Costa et al. (2018):

[TABLE]

Here the $\pm$ signs are mutually independent and the bound states have to appear inside the gap, $|\epsilon|\leq\Delta$. For a zero phase-difference, $\varphi=0$, there are only two states inside the gap, which correspond to the well-known YSR solution:

[TABLE]

The other two states remain at the gap edges, $\epsilon=\pm\Delta$, independently of the value of $\Phi$. Whereas the YSR are nondegenerate, ABS states, Eq (14), are double degenerate. Moreover, with increasing $\Phi$ the ABS cross zero energy every time $\Phi=(2n+1)\pi/2$. In contrast, YSR states cross the zero only once at $\Phi=1$, where, as explained below, a quantum phase transition takes place Sakurai (1970); Balatsky et al. (2006); Franke et al. (2011).

III.1 Spectrum in an intermediate range of parameters

We address now the question about the spectrum in an intermediate case, between the semiclassical and the YSR limits. This may correspond to a cluster of magnetic atoms or a small ferromagnetic island with a large but finite exchange field. The expression determining the bound states can be obtained from Eqs (6) and (III) and it is explicitly shown in the appendix, Eq (22). In Fig 2 we show with solid lines the subgap spectrum of the SFS structure as a function of the normalized length of the magnetic region, $k_{F}L$, for $\varphi=0$ . Different panels correspond to different values of the magnetic phase $\Phi$. For small $k_{F}L\lesssim 1$ there are only two nondegenerate states within the gap. These are the YSR states. Figures 2a and 2b correspond, respectively, to the situations before and after the YSR states cross at zero energy. Further increase of $\Phi$ pushes the states towards the gap edges. In contrast, for longer junctions, $k_{F}L\gg 1$, two pairs of bound states can be found within the gap. These pairs of states are non-degenerate (except at certain values of $k_{F}L$) and their energy oscillates with a period $2\pi/k_{F}$ around the semiclassical value determined by Eq (14). The oscillations stem from interference effects that are ignored in the semiclassical limit. Further increase of the junction length towards $L\sim\xi_{0}$ will bring additional bound states into the gap, which are not considered here.

It is worth noticing that Figs 2b-c show zero-energy crossings for finite length junctions at $\varphi=0$. At each crossing the total spin of the system change by one, as calculated from Eq (10) and shown by dashed black lines in Fig 2. In other words, Fig 2 demonstrates that a QPT also takes place beyond the YSR limit. Moreover, a sequence of QPTs with a stepwise change of the total spin may exist in a finite length junction.

The number of zero-energy crossings as a function of $L$ grows with increasing $\Phi$ 222See multimedia material in the supplementary information. As it follows from Eq (14), in a short junction within the semiclassical limit, $k_{F}^{-1}\ll L\ll\xi_{0}$, the ABS cross zero periodically at values of the magnetic phase $\Phi=(2n+1)\pi/2$. Each of these “asymptotic” crossings should be accompanied with, at least, two additional zero-energy crossings at intermediate values of $k_{F}L$ (Fig 2b-c). Fast oscillations of the bound state energies as a function of $L$ may increase further the number of zero-energy crossings by an even number (Fig 2c).

III.2 Josephson current

The subgap spectrum can be measured by means of tunneling spectroscopy Yazdani et al. (1997); Ji et al. (2008); Pillet et al. (2010); Ruby et al. (2015); Heinrich et al. (2018). In addition, measurements of the Josephson current in SFS junctions can also shed light on the spectral properties Goffman et al. (2000); Agraıt et al. (2003), in particular on the ground state of the junction. In conventional SNS junctions the Josephson energy is minimized when the phase difference vanishes, $\varphi=0$. However, it is known that, in SFS junctions, this minimum can also be found at $\varphi=\pi$ by tuning the exchange field or the length of the F region Bulaevskii et al. (1977); Buzdin et al. (1977); Lazar et al. (2000); Ryazanov et al. (2001); Chtchelkatchev et al. (2001). In the context of a delta-like magnetic impurity the connection between the zero energy YSR state and the [math]-$\pi$ transition has been recently discussed in Ref Costa et al. (2018). As we discuss next, the transition between the [math]- and $\pi$-junction behaviour is closely related to the QPTs described above for arbitrary junctions.

For this sake, we compute the ground state energy of the junction as a function of the phase difference $\varphi$. If the energy has a unique minimum at $\varphi=0$ or $\varphi=\pi$, one says that the junction is in the 0- or $\pi$- phase respectively. If the Josephson energy has minima both at $\varphi=0$ and at $\varphi=\pi$, then the ground state is denoted as $0^{\prime}$ or $\pi^{\prime}$ depending on the location of the global minimum Rozhkov and Arovas (1999); Vecino et al. (2003); Bergeret et al. (2006).

In Fig 3a we show the phase diagram in the $L$-$\Phi$ plane. This diagram provides an interesting connection: the QPTs associated with the zero-energy crossings shown in Figs 2(b-c) (horizontal dashed lines in Fig 3a), correspond to transitions between the ([math], $0^{\prime}$, $\pi^{\prime}$) states and the $\pi$ state.

Finally, in Fig 3b we show the dependence of the total spin of the system on $\Phi$ for a junction with $k_{F}L=10$. The dashed light line (dashed dark line) shows the spin if the junction is forced to stay in the math-state. The solid red line shows the spin of the system if the junctions always stays in the true ground state, i.e., if it follows the global energy minimum when the parameters are changed. Notice that, whenever the ground state corresponds to $\varphi=0$ ($\varphi=\pi$), the total electronic spin of the system is even (odd).

IV Conclusions

In conclusion, we present a complete study of equilibrium properties of a superconducting wire with a magnetic defect. We derive a general expression, Eq (6), that determines the full subgap spectrum provided that the T-matrix of the F region in the normal state is known. We also demonstrate in Eqs (10) and (11) that the total spin of a SFS junction undergoes integer jumps in units of $\hbar/2$ associated with zero-energy crossings of the bound states. Specifically, we analyze the spectrum of a one-dimensional ballistic SFS Josephson with a F region smaller than the superconducting coherence length but arbitrary strength of the exchange field. Our theoretical analysis bridges nicely two previously disconnected limiting cases: the YSR and the semiclassical ones. We demonstrate that the QPT predicted by the YSR model can be also found for SFS junctions of finite length $L$. Such phase transitions are associated not only to the integer jumps of the total spin described by our generalized Friedel sum rule, but also to a change of the sign of the supercurrent.

Acknowledgments

We acknowledge funding by the Spanish Ministerio de Economía y Competitividad (MINECO) (Projects No. FIS2014-55987-P, FIS2016-79464-P and FIS2017-82804-P). I.V.T. acknowledges support by Grupos Consolidados UPV/EHU del Gobierno Vasco (Grant No. IT578-13). M.R and F.S.B. acknowledge funding from the EU’s Horizon 2020 research and innovation programme under grant agreement No. 800923 (SUPERTED).

Appendix A T-matrix of the F region

Here we derive the normal state T-matrix of a ferromagnetic region of length $L$ and Zeeman splitting $h$ centered at the origin between two metallic electrodes. This matrix enters Eq (6) and hence it is pivotal to obtain the the bound states. In the normal state electrons and holes are decoupled, so we will only focus on the electrons. The wavefunction reads

[TABLE]

where $q_{\sigma}=k_{F}\sqrt{1+\frac{\epsilon+\sigma h}{\mu}}$ and $q_{0}=k_{F}\sqrt{1+\frac{\epsilon}{\mu}}$ are the wavenumbers at the ferromagnet and the normal metal, $k_{F}$ is the Fermi wavenumber and $\sigma=\pm$ stands for the spin orientation. From the continuity of the wavefunction in Eq (17) and its first derivative, we obtain a set of four equations that we have to solve. First writing $C^{\sigma}$ and $D^{\sigma}$ in terms of $A_{L}^{\sigma}$ and $B_{L}^{\sigma}$, and substituting them into the expressions for write $A_{R}$ and $B_{R}$, we finally get a connection between the wavefunction at the left and right superconductor,

[TABLE]

where

[TABLE]

and the remaining two components are related to these ones by complex conjugation, $T_{22}^{\sigma}=(T_{11}^{\sigma})^{*}$ and $T_{21}^{\sigma}=(T_{12}^{\sigma})^{*}$. As it is defined in Eq (4), the matrix in Eq (18) is the normal state transfer matrix of the ferromagnetic region.

Appendix B The S-F-S subgap spectra

Here we obtain the spectrum of a homogeneous S-F-S junction, whatever the values of the width and the exchange field strength of the magnetic region are. We start from the secular equation (6) and assume that $\hat{T}^{\sigma}$ is a generic $2\times 2$ matrix, like the one in Eq (18). After some algebra and exploiting the relations between the elements of the transfer matrix, we get a rather simple equation

[TABLE]

from which, substituting the expressions for the elements of the T-matrix in Eqs (19) and (20), we obtain

[TABLE]

In Eq (22), $\overline{q_{\sigma}}=k_{F}\sqrt{1-\frac{\epsilon+\sigma h}{\mu}}$ is the time conjugate of the electron wavenumber in F and we have approximated $q_{0}\approx k_{F}$, which is totally justified by the fact that $\Delta\ll\mu$ is fulfilled in any superconductor and that $q_{0}$ did not appear in any trigonometric function (where the accumulated phases along long distances would eventually be non-negligible, $\epsilon/\mu\cdot k_{F}L\sim 2\pi$).