Long-range effect of a Zeeman field on the electric current through the helical metal-superconductor interface in Andreev interferometer
A. G. Mal'shukov
cond-mat.supr-con cond-mat.mes-hall



Institute of Spectroscopy, Russian Academy of Sciences, Troitsk, Moscow, 108840, Russia

It is shown that the spin-orbit and Zeeman interactions result in phase shifts of Andreev-reflected holes propagating at the surface of  a topological insulator, or in Rashba spin-orbit-coupled two dimensional normal metals, which are in a contact with an s-wave superconductor.  Due to interference of  holes reflected through different paths of Andreev interferometer the electric current through external contacts varies depending on the strength and direction of the Zeeman field. It also depends on mutual orientations of Zeeman fields in different shoulders of the interferometer.  Such a nonlocal effect is a result of the long-range coherency caused by the superconducting proximity effect. This current  has been calculated within the semiclassical theory for Green functions in the diffusive regime,  by assuming a strong disorder due to elastic scattering of electrons.

74.45.+c, 74.78.-w, 74.25.Ha
Long-range  effect of a Zeeman field on the electric current through the helical metal-superconductor interface in Andreev interferometer.
    A. G. Mal'shukov
    Accepted, Received; in original form 
==========================================================================================================================================




§ INTRODUCTION


Due to a combined effect of a Zeeman field and  the spin-orbit coupling (SOC) the wave functions of Cooper pairs in s-wave superconductors acquire a  phase dependent factor. This phase is responsible for the magnetoelectric effect <cit.>, which leads to a spontaneous supercurrent in the presence of  a nonuniform static Zeeman field <cit.>, so that the spatial distribution of this current depends in a peculiar way on coordinate variations of the field. A similar phase also characterizes the electron wave function of a normal metal placed in a contact with a superconductor, if the strong enough spin-orbit and Zeeman interactions are presented in this metal. For example, it results in a spontaneous current through a  superconductor-normal metal-superconductor Josephson junction, the so called φ-junction <cit.> which has been observed experimentally in Ref.Szombati. These physical phenomena provide important building blocks for low dissipative spintronic applications based on interaction of magnetic and superconducting systems.

It is natural to expect that in superconductor-normal metal proximity systems the phase shift, which is induced by the Zeeman field and SOC, may be observed in the Andreev reflection <cit.>, where an electron scatters from a normal metal-superconductor interface as a hole.  Interesting possibilities for studying  the phase coherent phenomena are provided by Andreev interferometers <cit.>. These devices have several alternative paths for incident electrons and backscattered holes. In the previous studies a phase shift between interfering  scattered waves  has been provided by a magnetic flux. On the other hand, it is important to understand, if the Zeeman field can produce the phase shift that is strong enough to result in measurable effects on   the electric current through the Andreev interferometer. This problem has not been addressed  yet.

 It is clear that a strong enough SOC is needed to produce a  magnetoelectric effect which may be effective in a system of a micron size. Indeed, some two-dimensional (2D) systems have a strong intrinsic SOC <cit.>, which results in a considerable spin splitting of electron bands. In 2D systems these spin-split bands are characterized by opposite spin helicities. However, in the practically important semiclassical regime, when the Fermi energy (chemical potential) μ is larger than SOC, the magnetoelectric effect is reduced by a competition of bands with opposite helicities which cancel each other up to the terms ∼ h_F/μ, <cit.> where h_F is the spin orbit splitting at the Fermi energy. On the other hand, this cancelation does not occur in Dirac systems, such as surface electrons in a three dimensional topological insulator (TI), because in TI only the odd number of surface helical bands cross the Fermi energy. Therefore, it is reasonable to take  a TI wire as a basic component of the device. At the same time, it will be demonstrated that the results obtained for TI may  also be extended to a conventional 2D wire having a very strong SOC  h_F∼μ.

A simple interferometer is shown at Fig.1. Due to interference of paths through the upper and lower branches of the TI wire the electric current between the normal and superconducting leads can be varied by changing magnitudes or directions of Zeeman fields in the branches. For example,  the current might be changed by flipping a magnetization direction in one of the branches. Such a nonlocal dependence of the conductance would demonstrate a long-range phase coherence created in the the TI wire by the proximity effect at low enough temperatures. The Zeeman field in TI is assumed to be directed parallel to the x,y surface of the rectangular wire. It may be created by a ferromagnetic (antiferromagnetic) insulator deposited on top of TI, or by magnetic doping. Instead of fabricating TI wires, one could deposit superconducting and normal leads, as well as magnetic films on a TI flake. We will consider in detail the former setup, although qualitative results will be valid for both.

The electric current through the interferometer will be studied within the semiclassical theory for electron Green functions <cit.>. A strong elastic scattering on impurities will be assumed in the TI wire, so that the corresponding mean free path is much smaller than its dimensions. Also, the elastic scattering rate is much larger than the Zeeman splitting, but much less than  the chemical potential. At the same time, for sufficiently short wires  in the micrometer range, the low-temperature inelastic scattering of electrons will be ignored.

The article is organized in the following way. In Sec.II the Usadel equation and boundary conditions for the semiclassical Green function are derived for a TI wire. In Sec.III linearized  Usadel equations are derived for the case of a weak proximity effect and the analytic expression for the current is found in the low-bias regime. A summary of the results is presented in Sec.IV.








§ USADEL EQUATIONS


The effective one-particle Hamiltonian of  electrons on the  surface of TI can be written in the form <cit.>

    H=τ_3v𝐞(𝐤̂×σ)-τ_3μ+𝐙(𝐫)σ+V(𝐫),

where 𝐤̂=-i∂/∂𝐫 and the Pauli matrices  τ_1,τ_2,τ_3 operate in the  Nambu space, so that the electron destruction operators in the chosen basis have the form ψ_↑, ψ_↓, ψ^+_↓, -ψ^+_↑ with the arrows denoting spin directions. The third term in Eq.(<ref>) represents the Zeeman interaction, where 𝐙(𝐫) is parallel to the xy plane (the coordinate axes are shown in Fig.1), and the last term is a random impurity potential. 𝐞 is the unit vector which is parallel to the external normal to the wire surface. It is assumed that the wire width in the y-direction is much larger than its thickness in the z-direction. Therefore, electrons spend a relatively short time on flank surfaces. For this reason these surfaces are not taken into account in Eq.(<ref>).

The semiclassical Eilenberger equations for electron Green functions are obtained by expanding the Dyson equation with respect to small Fermi wavelengths, in comparison with other characteristic lengths. These equations serve for calculation of the so called semiclassical Green functions. The latter are obtained from initial Green functions by integration over the particle energy at a fixed momentum direction, which is represented by the unit vector 𝐧. These functions are  combined into the 2×2 matrix ĝ_𝐧(𝐫) whose components are g_11=g^r, g_22=g^a, g_12=g^K and g_21=0, where g^r, g^a and g^K are the retarded, advanced and Keldysh functions, respectively. These functions, in turn, are matrices in spin and Nambu spaces. The  procedure for the derivation of the Eilenberger equations is well described in literature <cit.>.  As long as all characteristic energies are much less than the Fermi energy, transitions between bands with opposite helicities can be neglected within the semiclassical approximation. In this case the spin dependence of the Green functions is locked to a momentum direction. Therefore, the initial Eilenberger equations can be projected onto the electron or hole helical bands, depending on a location of the Fermi level. The semiclassical Green function, in turn, takes the form ĝ_𝐧=ĝ_𝐧0(1±𝐧×σ)/2, where at μ>0 the "+" sign must be chosen and vice versa. The function ĝ_𝐧0 does not depend on spin and satisfies the normalization condition ĝ_𝐧0^2=1. For a dirty system, where the mean free path is smaller than other lengths, the Eilenberger equations can be transformed into diffusive Usadel equations <cit.> for  the matrix ĝ_0(𝐫), which is obtained from ĝ_𝐧0(𝐫) by averaging over 𝐧. By this way the Usadel equation has been obtained in Ref. <cit.> for Dirac electrons and in Ref.<cit.> for a superconductor with  Rashba SOC, which is larger than the elastic scattering rate. For the TI wire this equation can be written in the form

    D_t(b)∇̃(ĝ_0∇̃ĝ_0)+i[ωτ_3,ĝ_0]=0  ,

where ∇̃*=∇*+i[τ_3𝐅,*] and the gauge-field vector components are 𝐅=𝐙(𝐫)×𝐞_z/v. The parameters D_t and D_b denote electron diffusion coefficients on the top and bottom surfaces of the wire, respectively. In general these coefficients are different, because environments and surface potentials vary at these interfaces. It is interesting to note that the  Zeeman field enters Eq.(<ref>) in the same way as the vector potential of the magnetic field. An important difference is, however, that one can not change 𝐅 by a gauge transformation. Therefore, it is impossible to eliminate the "longitudinal" part of 𝐅 by such a transformation. In superconductors this part results in the so called helix phase with a spatially dependent order parameter, <cit.> as well as to spontaneous supercurrents around ferromagnetic islands. <cit.>


 When the wire length is much larger than its width w and  ∇_x g_0 is much smaller than w^-1, the Green function will tend to distribute uniformly over the wire width (in y-direction). If, in addition, g_0 is continuous on the wire flanks, it becomes constant around its perimeter.  Let us consider the case when  𝐅 is zero on the bottom surface. As shown in Appendix A, by averaging Eq.(<ref>) over y it can be reduced to the one-dimensional equation

    D∇̃_x(ĝ_0∇̃_xĝ_0)+i[ωτ_3,ĝ_0]-
          D(γ_x F_x^2+γ_y F_y^2)(τ_3ĝ_0τ_3ĝ_0-ĝ_0τ_3ĝ_0τ_3)=0 ,

where ∇̃_x*=∇_x* + i(D_t/2D)[τ_3F_x,*], D=(D_t+D_b)/2,  γ_x=D_tD_b/4D^2 and γ_y=D_t/2D. It should be noted that an equation of the same form may be obtained for a Rashba 2D electron gas with large SOC, such that h_F∼μ,  by formal replacing the constants γ and D_t/D_b with parameters from Ref.[Houzet], which depend on the ratio between the Rashba constant and the Fermi velocity.

Let us consider a weak coupling of the TI  wire to the superconducting lead through tunneling barriers, which are shown in Fig.1 at contact points 2 and 3. Therefore, Eq.(<ref>) has to be supplemented by boundary conditions (BC) at these interfaces. For  a 2D Dirac system the usual semiclassical BC <cit.> must be modified, as shown in Ref. [Zyuzin]. The modified BC  has the  form

    Dĝ_0∇̃_xĝ_0=Γ_S[ĝ_0,ĝ_s],

where ĝ_s is the Green function in the superconducting lead and Γ_S is a tunneling parameter on the interface of TI with the superconducting lead. This parameter can be written in terms of the barrier resistance R_b=ρ_TID/2Γ, where ρ_TI is the wire resistance per unit length. <cit.> The Green functions and 𝐅 in Eq.(<ref>) should be taken near barriers. If the Zeeman interaction vanishes near these interfaces, then 𝐅=0  and Eq.(<ref>) coincides with a conventional expression from Ref. <cit.>. Since it is assumed that the Zeeman interaction is induced by  magnetic layers on top of TI, it may vanish or not at the contacts, depending on sample preparation. It is expected that magnetization directions of the magnetic islands in the two interferometer arms may be varied independently of each other. Therefore, these islands must be separated to some extent in branching point 1.

A tunneling contact will be also assumed at the interface of the TI wire with the normal lead at point 1. At this point the Green functions of electrons in both TI branches coincide. One more BC is an evident generalization of  Eq. (<ref>) that takes into account two branches which make a contact with the normal lead. We apply here the ideas of Refs. [Zaitsev,Stoof] on how to write BC in branching points. By assuming that 𝐅=0 at  contact  point 1, this BC can be written as

    Dĝ_0∇_x_2ĝ_0 +Dĝ_0∇_x_3ĝ_0=-Γ_N[ĝ_0,ĝ_N],

where x_2 and x_3 are coordinates in the branches. They are chosen so, that x_2 and x_3  are directed from contact 1 towards respective contacts 2 and 3 with the superconductor. The tunneling parameter Γ_N may be expressed through the barrier resistance R_b1=ρ_TID/2Γ_N, in the same way as for the TI-S contact. For the massive normal lead one may assume that its Green function is unperturbed by a contact with the TI wire. Therefore, ĝ_N^r/a=±τ_3





§ ANDREEV REFLECTION AND ELECTRIC CURRENT


We consider the case of the low temperature T and small bias voltage V, which are much less than the superconducting gap. Therefore, the electric current between the normal and superconducting leads is determined by the Andreev reflection. This current may be expressed via the conductance G(ω), according to the well known expression <cit.>

    j=1/e∫ dω[tanhω+eV/2k_BT-tanhω-eV/2k_BT]G(ω) .

Let us focus on the high barrier regime, when the barrier resistance R_b at TI-S interface  is much larger than the resistance of the TI wire and the barrier resistance R_b1 at the TI-N interface. In this case G(ω) is given by <cit.>

    G(ω)=1/8R_b(M_2+M_3) ,

where

    M_2(3)=Tr[(g_0^rτ_3-τ_3g_0^a)(g_s^rτ_3-τ_3g_s^a)]|_x_2=L_2(x_3=L_3) .

The functions g_0  are taken in TI wires near contacts 2 and 3. L_2 and L_3 are the lengths of the wires between contact 1 and contacts 2 and 3, respectively.  We assume a massive superconducting lead whose Green function is not perturbed by a proximity to TI wires. Therefore, at both contacts these functions have the form g^r_s=g^a_s=(-iτ_3ω + τ_2 Δ)/√(Δ^2-ω^2) for Δ>ω. At high R_b the Green functions in TI are weakly perturbed by the superconductor, so that  they can be represented as sums of unperturbed functions and small corrections δ g^r(a)_j, namely

    g^r(a)_0(x_j)=±τ_3+δ g^r(a)_j ,

where j=2,3. The functions δ g^r(a)_j≪ 1 are the anomalous Green functions which are nondiagonal in the Nambu variables. By linearizing Eq.(<ref>) with respect to δ g^r(a)_j it can be transformed to

    D((-1)^j∇_x_j+2iτ_3F̃_xj)^2δ g^r(a)_j± 2 iωδ g^r(a)_j-
          4D(γ_x F_x^2+γ_y F_y^2)δ g^r(a)_j=0 ,

where F̃_xj=F_x (x_j)D_t/(D_t+D_b). In turn, boundary conditions Eq.(<ref>) take the linearized form

    D(∇_x_j+2iτ_3F̃_xj)δ g^r(a)_j|_x_j=L_j=Γ_Sτ_3 [τ_3,g_s^r(a)]

At the same time, M_2 and M_3 become

    M_j=2Δ/√(Δ^2-ω^2)Tr[(δ g^r_j+δ g^a_j)τ_2]|_x_j=L_j .




The solutions of Eq.(<ref>) contain the phase factors exp(± 2i∫ dx_jF̃_xj) which result in spatial oscillations of Green functions. Besides these oscillations, the Zeeman interaction leads to a suppression of the superconductor proximity effect. For instance, due to the third term in Eq.(<ref>), δ g_2 and δ g_3 decrease with increasing  distances from contacts 2 and 3, respectively. Therefore, the length L_Z of the region where 𝐙≠ 0 should not be too long. The corresponding condition is 2L_Z(γ_x F_x^2+γ_y F_y^2)^1/2≲ 1. By choosing the direction of 𝐙 perpendicular to x (F_y=0), the suppression effect can be reduced in samples having the  smaller ratio D_b/D_t of the diffusion constants, as follows from the definition of γ_x. It is also possible to construct appropriate barriers at the flanks of the wire to guarantee a weak Klein tunneling between the top and bottom surfaces. By making the angular averaged tunneling rate much less than the Thouless energy E_T=D/L^2, where L=max[L_2,L_3], the bottom surface of TI may almost completely be turned off, that will result in the small damping effect.  It should be noted that the third term in Eq.(<ref>) vanishes completely if the Zeeman fields are finite on both surfaces and  are equal in magnitude and antiparallel (both are perpendicular to x). However, such a situation is probably difficult to realize in practice.






 §.§ Short wires, low bias regime


A simple analytic result may be obtained in the case of V≪ k_BT  at small enough L_2 and L_3, so that k_BT≪ E_T. In this case one may set ω=0 in G(ω) in Eq.(<ref>). Let us assume that  L_2=L_3=L and L_Z is slightly less than L <cit.>. Hence, the phase Φ(x_j)≡2(-1)^j∫_0^x_j dx_jF̃_xj≃ 2(-1)^jF̃_xjx_j. The solutions of Eq.(<ref>) in both TI branches have the form δ g_j=exp(iτ_3Φ_j)[A_jexp(κ_j x_j)+B_jexp(-κ_j x_j)], where κ_j^2=4γ_x F_xj^2± 2iω at F_y=0 (± for retarded and advanced functions, respectively). In a symmetric device, that will be assumed below for simplicity, |F_x2|=|F_x3|. The 2×2 matrices A and B can be obtained from boundary conditions Eq.(<ref>), Eq.(<ref>) and the continuity of Green functions of the wire branches in contact point 1 . By substituting the so calculated δ g_2 and δ g_3 into Eq.(<ref>) we obtain the current from Eq.(<ref>) in the form


    j=ρ_TI/R_b^2Re[β(α + cosΔΦ)]V ,

where ΔΦ=Φ_2|_x_2=L-Φ_3|_x_3=L,

    α=   2(1+Λ)sinh^2κ_0 L+1 ,  
       β=   2/κ_0(1+Λ)sinh2κ_0 L ,

Λ=ρ_TIκ_0 L/2κ_0 R_b1 and κ_0=κ|_ω=0. For more details of the calculation, see Appendix B. It follows from these expressions that the oscillating part of the current may be of the same order as the constant term, if κ L ≲ 1 and Λ≲ 1.  As can be seen from Fig.2, the current's oscillations are strongest at D_t/D_b=0 and they are strongly damped already at D_t/D_b=0.1. The oscillations  almost vanish at D_t/D_b=0.5. In the considered symmetric device the phase-dependent part of the current and the oscillations turn to zero when the Zeeman fields at two branches are antiparallel, so that in Eq.(<ref>) ΔΦ=0. The difference of conductances Δ G for the parallel and antiparallel alignments is shown in Fig.3 at various ratios D_t/D_b and the zero temperature (T≪ E_T).  An alignment switch can be performed by changing a magnetization in one of the magnetic islands. For example, one may adjust their hysteretic characteristics in such a way that an external  magnetic field of a definite strength flips the magnetization of one of them, while the other island stays in its initial state.

It is important that the considered in this subsection short wire regime is valid at low enough temperatures which provide the sufficiently large coherence length ξ=√(D/k_BT), such that  ξ≫ L. Otherwise, one can not simply set ω=0 in G(ω). Instead of that, the integral over ω in  Eq.(<ref>) must be taken.





 §.§ High temperatures


In this subsection the numerical results are presented beyond the short wire regime, at k_BT ≳ E_T. G(ω) can be obtained from Eqs.(<ref>) and (<ref>). In turn, the Green functions, that enter in Eqs.(<ref>), are calculated in Appendix B. Fig.4 shows differences of conductances for parallel and antiparallel alignments of Zeeman fields, at various ratios of k_B T and E_T. These plots are normalized by the temperature dependent conductance G_0T in the absence of the Zeeman field. Fig.4 shows that so normalized Δ G decreases in the considered temperature interval , but not dramatically, that makes it possible to observe the phase shift produced by the Zeeman field even at relatively high temperatures. Note, that the absolute reduction of Δ G is larger, considering almost  a threefold decrease of G_0T  in the same temperature interval. In order to evaluate E_T, let us take the mean free path l=10 nm, as in Bi_1.5Sb_0.5Te_1.7Se_1.3 <cit.> and a typical Dirac velocity v=5 10^5 m/s, that gives the diffusion constant of a 2D gas D=vl/2=25 cm^2/s. With this constant the Thouless energies are   E_T=80mK and 20mK  for the interferometer shoulders L=500nm and 1000nm, respectively.


Now let us evaluate typical values of Z which can provide the strong enough phase shift Φ =2ZLD_t/v(D_t+D_b). As can be seen from Figs.3 and 4, the maximum effect on Δ G is observed for 1≲ΔΦ≲ 2. For D_t/(D_t+D_b) ≃ 1, v=5· 10^5m/sec and  L= 1μm the phase  ΔΦ=2Φ reaches 1.5 at Z≃ 0.1meV. Such a field is  well below the Fermi energy, that is in agreement with the semiclassical approximation used in this work. Note, that the above evaluation of the Zeeman field is valid only for a special case of the magnetic island which covers almost the entire TI wire. Therefore, the field must be stronger for smaller sizes of the islands.







§ CONCLUSION


In conclusion, it is shown that due to a quantum interference of Andreev-scattered waves in  wires made of three dimensional TI,  the electric current through a TI-superconductor system can be varied by changing the mutual orientations of Zeeman fields in distant parts of the TI wire. This effect is a direct consequence of the long-range Cooper correlations created by the superconducting proximity effect and the Zeeman-field-induced phase shifts of the pairing functions. This effect is damped at strong Zeeman fields. The damping can be reduced by a special design of the interferometer. On the other hand, it is shown that the discussed interference effects may be observed  even at weak fields, due to the strong spin-orbit coupling of TI surface states. With some modification of parameters the theory may be extended  to ordinary two-dimensional electron systems with sufficiently strong Rashba interaction.

Acknowledgement. The work was supported by RAS Program "Actual problems of low-temperature physics".



99
Edelstein
V. M. Edelstein, Sov. Phys. JETP 68, 1244 (1989)
Malsh island
A.G. Mal'shukov, Phys. Rev. B 93, 054511 (2016).
Pershoguba
S. S. Pershoguba, K. Björnson, A. M. Black-Schaffer, and
A. V. Balatsky, Phys. Rev. Lett. 115, 116602 (2015).
Hals
Kjetil M. D. Hals, Phys. Rev. B 95, 134504 (2017)
Krive
 I. V. Krive, A. M. Kadigrobov, R. I. Shekhter and M. Jonson, Phys. Rev. B 71, 214516 (2005)
Reinoso
A. Reynoso, G.Usaj, C.A. Balseiro, D. Feinberg, M.Avignon, Phys. Rev. Lett. 101, 107001 (2008)
Zazunov
 A. Zazunov, R. Egger, T. Martin, and T. Jonckheere, Phys.Rev. Lett. 103, 147004 (2009)
ISHE
A. G. Mal'shukov, S. Sadjina, and A. Brataas, Phys. Rev. B 81, 060502 (2010)
Liu
 J.-F. Liu and K. Chan, Phys. Rev. B 82, 125305 (2010)
Yokoyama
 T. Yokoyama, M. Eto, Y. V. Nazarov, Phys. Rev. B 89, 195407 (2014)
Konschelle
 F. Konschelle,  I. V. Tokatly and F. S. Bergeret, Phys. Rev. B 92,125443 (2015)
Szombati
D. B. Szombati, S. Nadj-Perge, D. Car, S. R. Plissard, E. P. A. M. Bakkers, L. P. Kouwenhoven, Nat. Phys. 2,
568 (2016).
Sakano
M. Sakano, M. S. Bahramy, A. Katayama, T. Shimojima, H. Murakawa, Y. Kaneko, W. Malaeb, S. Shin,
K. Ono, H. Kumigashira, R. Arita, N. Nagaosa, H.Y. Hwang, Y. Tokura, and K. Ishizaka, Phys.Rev. Lett. 110, 107204 (2013)
Ast
C. R. Ast, J. Henk, A. Ernst, L. Moreschini, M. C. Falub,
D. Pacilé, P. Bruno, K. Kern, and M. Grioni, Phys. Rev.
Lett. 98, 186807 (2007).
Lesne
E. Lesne,Y. Fu, S. Oyarzun, J.C. Rojas-Sanchez, D.C.
Vaz, H. Naganuma, G. Sicoli, J.-P. Attane, M. Jamet,
E. Jacquet, J.-M. George, A. Barthelemy, H. Jaffres, A.
Fert, M. Bibes and L. Vila, Nature Materials (2016);
doi:10.1038/nmat4726
Song
Qi Song, Hongrui Zhang, Tang Su, Wei Yuan, Yangyang
Chen, Wenyu Xing, Jing Shi, Ji Rong Sun, and Wei Han,
arXiv:1609.06207
Andreev
A. F. Andreev, Zh. Eksp. Teor. Fiz. 46, 1823 (1964) [JETP 19, 1228 (1964)].
Zaitsev
A. V. Zaitsev, Physica B 203, 274 (1994)
Stoof
T. H. Stoof, Yu. V. Nazarov, Phys. Rev. B 53, 14496 (1996); Yu. V. Nazarov and T. H. Stoof, Phys. Rev. Lett 76, 823 (1996)
Golubov
A. A. Golubov, F. K. Wilhelm, A. D. Zaikin, Phys. Rev. B 55, 1123 (1997)
Lambert
C.J. Lambert and R. Raimondi, J. Phys. Condens. Matter 10, 901 (1998)
Eilenberger
G. Eilenberger, Z.Phys. 214, 195 (1968)
Larkin semiclass
A. I. Larkin, and Y. N. Ovchinnikov, Zh. Eksp. Teor. Fiz. 55, 2262 (1968) [Sov. Phys. JETP 28, 1200 (1965)].

Usadel
K.D. Usadel, Phys. Rev. Lett. 25, 507 (1970)
Qi RMP
X. L. Qi and S. C. Zhang, Rev. Mod. Phys. 83, 1057

Rammer
J. Rammer, H. Smith, Rev. Mod. Phys. 58, 323 (1985)
Kopnin
N. Kopnin, Theory of Nonequilibrium Superconductivity (Oxford Science, London, 2001).


Zyuzin
A. Zyuzin, M. Alidoust, and D. Loss, Phys. Rev. B 93, 214502 (2016).
Bobkova
I. V. Bobkova, A. M. Bobkov, A. A. Zyuzin, and M.Alidoust, Phys. Rev. B 94, 134506 (2016)
Linder
Henning G. Hugdal, Jacob Linder, and Sol H. Jacobsen, Phys. Rev. B 95, 235403 (2017)
Houzet
Manuel Houzet and Julia S. Meyer, Phys. Rev. B 92, 014509 (2015).

Samokhin
V. P. Mineev and K. V. Samokhin, Zh. Eksp. Teor. Fiz. 105, 747 (1994) [Sov. Phys. JETP 78, 401 (1994)]
Kaur
 R. P. Kaur, D. F. Agterberg, and M. Sigrist, Phys. Rev. Lett. 94, 137002 (2005)
 Agterberg
 D. F. Agterberg, Physica C 387, 13 (2003)
Dimitrova
O. Dimitrova and M.V. Feigel'man, Phys. Rev. B 76, 014522 (2007)
Barzykin
V. Barzykin and L. P. Gor'kov, Phys. Rev. Lett. 89, 227002 (2002)

Zaitsev BC
A. V. Zaitsev, Sov. Phys. JETP 59, 1015 (1984).
Kupriyanov
M. Y. Kuprianov and V. F. Lukichev, Sov. Phys. JETP 67, 1163 (1988).
Volkov
A.F.Volkov, A.V.Zaitsev, and T.M.Klapwijk, Physica C  210 , 21 (1993).
comment
As long as characteristic lengths of Z gradients near magnetic island edges, as well as distances from the island edges to contacts are much less than the  coherence length √(D/k_BT) and L, one may represent the spacial dependence of Z as a stepwise  function.
TIflake
AE.S. Tikhonov, D.V. Shovkun, V.S. Khrapai, M. Snelder, M.P. Stehno, A. Brinkman, Y. Huang, M.S. Golden, A.A. Golubov, Phys. Rev. Lett. 117, 147001 (2016)






§ AVERAGING OVER A WIRE PERIMETER


Let us assume that the wire has a rectangular cross section and the coordinate l runs along its perimeter, so that it coincides with y and -y on the top and bottom surfaces, respectively. Then, the part of Eq.(<ref>), which is associated with the derivative over l, can be represented in the form:

    ∇_l(D(l)ĝ_0∇̃_̃l̃ĝ_0)+ iD(l)[F_lτ_3,ĝ_0 ∇̃_̃l̃ĝ_0]                ,

where  ∇̃_̃l̃*=∇_l*+i[τ_3F_l,*] and F_l represents a projection of the field 𝐅 onto the l coordinate. For completeness, the lateral surfaces of the wire (l ∥ z) are also taken into account in Eq.(<ref>). D(l) denotes the l-dependent diffusion constant.

 We will denote the average over the perimeter as (...)=∮ (...)dl/∮ dl and assume that ĝ_0 is constant as a function of l. For a diffusive transport the latter assumption is valid if the perimeter is  much smaller than the characteristic lengths which characterize variations of  Green functions  along the wire. Therefore, ∇_lĝ_0=0. Hence, ∇̃_̃l̃ĝ_0=i[τ_3F_l,ĝ_0] in Eq.(<ref>). Further, since F_l and D(l) are periodic functions of l, the average of the first term in Eq.(<ref>) is 0. Therefore, the averaging of  Eq.(<ref>) yields

    -D(l)F_l^2[τ_3,ĝ_0[τ_3,ĝ_0]]   .


By averaging the remaining terms in Eq.(<ref>) over the perimeter we arrive to the one-dimensional equation

    D∇̃_x(ĝ_0∇̃_xĝ_0)+i[ωτ_3,ĝ_0]+
    (DF_x^2/D-DF_x^2-DF_l^2)[τ_3,ĝ_0[τ_3,ĝ_0]] =0        .

By assuming that the thickness of the wire is much smaller than its width one may neglect the contribution of the lateral surfaces into the average. If F_x and F_y are finite only on the top surface, Eq.(<ref>) reduces to Eq.(<ref>)



§ DERIVATION OF EQ.(<REF>)


In each shoulder j the substitution δ g_j=e^iΦ(x_j)f_j , where Φ(x_j)≡2(-1)^j∫_0^x_j dx_jF̃_xj,  allows to transform Eq.(<ref>) to the form

    ∇_x_j^2 f^r(a)_j +κ_j^2f^r(a)_j=0 ,

where κ^2_j=4γ_xF^2_xj± 2i(ω/D) at F_y=0. The "±" signs in κ^2 correspond to retarded and advanced functions, respectively. Due to coordinate dependence of F^2_xj the parameter κ varies with x_j. If F^2_xj is a step-function, the wire can be divided into several parts, so that  in each of them κ^2 is a constant. The function f and its derivative must be continuous at boundaries between these parts, as it follows from Eq.(<ref>). Let us consider a simple case where the homogeneous magnetic islands in each shoulder occupy almost the entire wire, except for small regions near contacts with the leads. When the lengths of these regions is much smaller than the coherence length √(D/2|ω|)∼√(D/ k_BT), the function f and its derivative are almost constant there. Therefore, by neglecting their weak spacial variation one may replace in BC (<ref>) and (<ref>) the function f and ∇_x f with corresponding values in an adjacent magnetic domain. By this way it is possible to skip the small nonmagnetic regions of the wire. 

The solutions δ g_2 and δ g_3 of   in TI branches 12 and 13, respectively, have the form

    δ g_2   =   e^iτ_3Φ(x_2)[A_2e^κ x_2+B_2e^-κ x_2], 
    δ g_3   =   e^iτ_3Φ(x_3)[A_3e^κ x_3+B_3e^-κ x_3] ,

where the labels r and a are skipped for a wile. They will be restored later, where necessary. In the assumed symmetric case there is a common factor κ in both branches. The four coefficients A_j and B_j  anticommute with τ_3 and  can be found from the boundary conditions. According to the definition of the phase Φ,  we have Φ(0)=0 and Φ(x_j)|_x_j=L≡Φ_j at contacts j=2 and j=3, respectively, where Φ_j≃ 2(-1)^jF̃_jL. From the boundary conditions Eq.(<ref>), Eq.(<ref>) and the continuity of the Green functions in branches 2 and 3 at contact 1, it is easy to obtain the following equations near contact 1

    (A_2+B_2)   -   (A_3+B_3)=0 
    
    (A_2-B_2)   +   (A_3-B_3)=Λ_Nδ g(0) ,

where Λ_N=ρ_TI/R_b1κ and δ g(0)=A_2+B_2=A_3+B_3 is the Green function at contact 1.
At contacts 2 and 3 the boundary conditions have the form:

    A_2e^κ L-B_2e^-κ L   =   e^-iτ_3Φ_2τ_2Λ_S, 
    
    A_3e^κ L-B_3e^-κ L   =   e^-iτ_3Φ_3τ_2Λ_S ,

where Λ_S=(Δ/√(Δ^2-ω^2))(ρ_TI/R_bκ). From equations  Eq.(<ref>), Eq.(<ref>) the factors A and B can be expressed as

    A_j=(Φ_+/sinhκ L-(-1)^jΦ_-/coshκ L)τ_2Λ_S/4- 
          Λ_Nδ g(0) e^-κ L/4sinhκ L ,
          B_j=(Φ_+/sinhκ L+(-1)^jΦ_-/coshκ L)τ_2Λ_S/4- 
          Λ_Nδ g(0) e^κ L/4sinhκ L ,

where Φ_±=exp(-iτ_3Φ_3)±exp(-iτ_3Φ_2).

By calculating δ g(0)=A_2+B_2 from Eq.<ref>  we obtain the expression for δ g(0) in the form

    δ g(0)=Φ_+ τ_2Λ_S/2sinhκ L+Λ_N coshκ L .

According to Eq.(<ref>), the spectral conductance G(ω) is proportional to (M_2+M_3). The latter may be expressed from Eq.(<ref>) through the sum δ g^r(a)_2(L)+δ g^r(a)_3(L). By substituting coefficients A and B, that are given by Eq.(<ref>), into Eq.(<ref>) at x_2=x_3=L we obtain

    δ g^r_2(L)+δ g^r_3(L)+δ g^a_2(L)+δ g^a_3(L)=τ_2Λ_S ×
          Re[|Φ_+|^2(κ L+Λ_N/2)/1+Λ_N/2κ L+|Φ_-|^2tanhκ L] .

It is easy to see that in this expression only the phase difference Φ_2-Φ_3 enters, as it should be. Eq.(<ref>) finally gives the result Eq.(<ref>) in the low-bias regime where ω may be set to zero.

It is instructive to see how the phase dependence of the current vanishes in the case when only one of the two interferometer arms is conducting. Let us, for example, turn off branch 3. In this case only the second lines should be left in BC Eqs.(<ref>) and (<ref>), where A_3=B_3=0. It is easy to see that the solutions of these equations at x_2=L have the form A_2=exp(-iτ_3Φ_2)f_a(κ) and B_2=exp(-iτ_3Φ_2)f_b(κ), where the functions f_a/b(κ) do not depend on the phase Φ_2. They depend only on κ. Therefore, the function δ g_2(L), which is given by Eq.(<ref>), does not depend on Φ_2, as well as the conductance G, as can be seen from Eqs.(<ref>) and (<ref>) at δ g_3=0.
Therefore, the only effect of the Zeeman field is a suppression of the proximity effect by the damping factor κ. It produces only a monotonous decreasing of the current at higher Zeeman fields and does not depend on its sign.


