Density of states in the presence of spin-dependent scattering in SF bilayers: a numerical and analytical approach

Introduction

It is well-known that superconductivity can be induced in a non-superconducting metal in hybrid structures due to the proximity effect [1-7]. For instance, in NS bilayers (where N denotes a normal metal and S denotes a superconductor), the superconducting correlations penetrate into the normal metal layer over a characteristic decay length ξ_n = where D_n is the diffusion constant in the normal metal and T_c is the transition temperature. When a superconductor S is combined with a ferromagnetic layer F, forming an SF bilayer, the superconductivity leaks into the ferromagnetic region over the characteristic length ξ_h = where D_f is the diffusion constant in F layer and h is the exchange field in the ferromagnetic layer [1]. Not only superconductivity is substantially suppressed due to the exchange field, but also Cooper pairs gain a finite center of mass momentum, which leads to the oscillatory behavior of the Cooper pair wave function. These oscillations can be described in the diffusive limit in the framework of the so-called Usadel equations, which are written in terms of the quasiclassical Green’s functions. This approach proved to be very powerful for the description of the proximity effect in diffusive superconducting hybrids [1-3,8-11].

The scientific community has been examining the proximity effect in SF hybrid structures already for a long while. It has been found that the oscillatory behavior of the superconducting wave function can lead to various interesting phenomena that can be observed experimentally [1,2]. For instance, the superconducting transition temperature shows non-monotonous and, in some cases, oscillatory behavior in multilayered SF structures [12-17]. Recently, it has been shown theoretically that similar behavior can be observed in S/TI structures with non-uniform magnetization patterns on the surface of a 3D topological insulator (TI) [18]. The Josephson critical current demonstrates damped oscillatory behavior as a function of the thickness of the ferromagnetic layer in SFS Josephson junctions [19-47]. Similarly, the density of states (DOS) also demonstrates a damped oscillatory dependence as a function of the F layer thickness in SF systems [48-52].

The density of states is one of the crucial spectral characteristics of the proximity effect in superconducting hybrid structures. For example, the DOS calculation is essential for the quasiparticle current computation in SIFS (where I denotes an insulating layer) [29,51,53-57] or SFIFS tunneling Josephson junctions [58]. Therefore, computation of the DOS is also needed for many actively studied areas of research, including the thermospin [59,60] and thermoelectric [61-66] effects, spin and heat valves [67-75], as well as nanoscale refrigerators [76-78]. Presently, the DOS structure at the free edge of a normal metal layer in NS bilayers is well-known [1-3,79]. It has a so-called mini-gap at the subgap energies E < Δ (where Δ is the superconducting gap), whose magnitude depends on the NS interface parameters and the thickness of the N layer [79,80]. Replacing the N layer with a ferromagnetic metal F results in a more sophisticated DOS structure since there is a non-zero exchange field, which causes spin-split densities of states for two spin populations of electrons [1-3,81]. More general considerations should also include possible spin-flips as well as spin–orbit scattering processes in the ferromagnetic region [82].

In this work, we consider a diffusive SF bilayer, assuming a relatively low interface transparency and the presence of magnetic and spin–orbit scattering. For this purpose, the Kupriyanov–Lukichev (KL) boundary conditions at the superconductor/ferromagnet interface are perfectly suitable [83]. We build the model employing the Usadel equations. The quasiclassical Usadel equations are based on the Matsubara Green’s functions and represent a powerful microscopic tool for the treatment of the diffusive superconducting hybrid structures. This approach is valid as long as the elastic scattering length l is much smaller than the superconducting coherence length ξ = where D_s is the diffusion constant in the superconductor. Previously, the DOS in SF bilayers has been studied numerically [84,85]. We revisit this question and propose an analytical model to describe the influence of spin-flip and spin–orbit scattering on the DOS behavior. Then, we provide a comparison with the exact numerical calculation using a self-consistent two-step iterative method. Furthermore, we briefly discuss the consequences of the different kinds of scattering on the current–voltage characteristics in SFIFS junctions. We do not consider any additional effects in the SF boundary such as spin-dependent interfacial phase shifts (SDIPS). The effect of SDIPS on the DOS behavior in SFIFS junctions has been studied both analytically [86] and numerically [87].

The paper is organized as follows. In the section (“Model”) we formulate the theoretical model. In the following sections, the derivation of the analytical results is presented. We discuss the calculations in the section (“Results and Discussion”), and finally we summarize the results in the last section (“Conclusion”).

Model

The theoretical model of the SF structure under consideration is depicted in Figure 1. It consists of a ferromagnetic layer with thickness d_f and a superconducting electrode along the x direction. The SF interface is characterized by the dimensionless parameter γ_B = R_Bσ_n/ξ_f, where R_B is the resistance of the SF interface in units Ω·m², σ_n is the conductivity of the F layer [88,89], ξ_f = D_f is the diffusion coefficient in the ferromagnetic metal, and T_c is the critical temperature of the superconductor [1,2]. We assume ℏ = k_B = 1. We also assume that the SF interface is not magnetically active. We will consider the diffusive limit in this model and neglect the nonequilibrium effects in the structure [90-92].

Our goal is to find the DOS of a single SF bilayer, which can be done by solving the Usadel equations in the ferromagnetic and superconducting layers. We employ the θ parametrization of the normal and anomalous quasiclassical Green’s functions, G = cosθ and F = sinθ, respectively. In the spin space, the anomalous functions are parameterized as

and we can write the Usadel equations in the F layer as [4,51,82,84,93]

where the positive and negative signs correspond to the spin-up (↑) and spin-down (↓) states, respectively. In terms of the electron fermionic operators ψ, the spin-up state corresponds to the anomalous Green’s function F_↑ ∼ ⟨ψ_↑ψ_↓⟩, while the spin-down state corresponds to F_↓ ∼ ⟨ψ_↓ψ_↑⟩. We use the Matsubara Green’s functions, hence, ω = 2 πT(n + 1/2) are the Matsubara frequencies [94]. The exchange field of the ferromagnet is h, and the scattering times are labeled here as τ_z, τ_x and τ_so. The parameter τ_z(x) corresponds to the magnetic scattering parallel (perpendicular) to the quantization axis, and τ_so is the spin–orbit scattering time [82,84,95].

In the model under consideration, we assume a rotational symmetry around the axis, that is, τ_x = τ_y. We consider the possibility of anisotropic magnetic scattering times. These can occur in ferromagnetic superconductors where magnetic disorder can be characterized by two scattering times τ_x = τ_y and τ_z [82]. The Usadel equations employed in this work are the special case of the more general formalism derived by Ivanov and co-workers [96]. Thus, the equations that we use in our calculations are exactly a quasiclassical limit of the model derived in [96] except that in our model we also include spin–orbit scattering. Ivanov et al. describe the spin-flip scattering via a symmetric matrix of scattering rates Depending on the symmetry of the problem, the number of parameters can be different (up to six), for example, in the most symmetric isotropic case, the spin-flip matrix is diagonal and described by one scattering time (τ_x = τ_y = τ_z). In the case of a ferromagnetic ordering in which disorder is anisotropic (one axis is special), the spin-flip rates in the directions along the ferromagnetic axis and perpendicular to it may be different. In this case, as we mentioned earlier, two scattering times would be necessary.

In the S layer, the Usadel equation has the following form [93]:

Here, D_s is the diffusion coefficient in the superconductor, and Δ(x) is the superconducting order parameter (pair potential). From the Usadel equations, it can be shown that there is a symmetry relation between θ_↑ and θ_↓: θ_↑(E) = (−E), where E is the energy (ω_n → − iE) and ^* is the complex conjugation. Equation 2 and Equation 3 should be supplemented with the self-consistency equation for the coordinate dependence of superconducting order parameter Δ,

The resulting system must be complemented by the boundary conditions at the outer boundary of a ferromagnet,

and the Kupriyanov–Lukichev boundary conditions at the FS interface [83],

Here γ = ξ_sσ_n/ξ_fσ_s, σ_s is the conductivity of the S layer, and ξ_s = is the superconducting coherence length. The indices of the θ-parameterized Green’s functions are omitted, since there is no mixing between the components θ_↑ and θ_↓. In other words, we are not considering spin-active interfaces. In the case of spin-active barriers, one should use the boundary conditions introduced in [87,97,98], rather than the standard Kupriyanov–Lukichev boundary conditions in Equation 6 and Equation 7. The parameter γ determines the strength of superconductivity suppression in the S layer by the ferromagnet F (inverse proximity effect). For instance, when γ ≫ 1, the inverse proximity effect is very strong, and the order parameter is heavily suppressed near the SF interface compared to its bulk value. On the contrary, when γ = 0, there is no suppression of the order parameter because there is no inverse proximity effect. In all numerical simulations, we assume that γ ≪ 1, that is, there is almost no superconductivity suppression in the superconductor. The transparency parameter γ_B is proportional to the interface resistance. In the regime of a fully transparent junction, γ_B = 0, the proximity effect is the strongest, and the θ functions are continuous at the SF interface while γ_B ≫ 1 corresponds to the tunnel junction limit, that is, there is no mutual impact between superconducting layer S and ferromagnetic layer F.

To complete the boundary problem, we also set a boundary condition at x = +∞:

where the Green’s functions take the well-known bulk BCS form. The Green’s function method allows us to compute the DOS at the outer F boundary by solving the resulting system of equations above.

The DOS at the outer F boundary N_f(E) is normalized to the DOS in the normal state and can be written as

where N_f↑(↓)(E) are the spin-resolved DOS written in terms of the spectral angle θ,

To calculate Equation 10, we use a self-consistent two-step iterative method. In the first step, we calculate the pair potential coordinate dependence Δ(x) using the self-consistency equation (Equation 4) in the S layer. Then, by proceeding to the analytical continuation in Equation 2 and Equation 3 over the quasiparticle energy iω → E + i0 and using the Δ(x) dependence obtained in the previous step, we find the Green’s functions by repeating the iterations until convergency is reached.

The DOS in the limit of small F layer thickness

In this section, we obtain the analytical result assuming d_f ≪ min(ξ_f, ), which is the case in a thin and weak ferromagnet. Under the condition γ = 0, we can neglect the suppression of superconductivity in the superconductor. Hence, the problem can be reduced to the rigid boundary condition when the order parameter in the S layer is set to its bulk value Δ₀. We will keep all scattering terms in the solution to obtain a more general result. In this case, we can expand the solution of the Usadel equations up to the second order in small spatial gradients. The θ_f functions can be approximated in the following way:

where the coefficients A_↑(↓), B_↑(↓), C_↑(↓) are determined from the boundary conditions.

Inserting the solution in Equation 11 into the Usadel equation in the F layer (Equation 2), we get,

For convenience, we introduced the scattering rate parameters α_z = 1/τ_zΔ, α_x = 1/τ_xΔ, and α_so = 1/τ_soΔ. To find the coefficients, we utilize the boundary conditions in Equation 5 and Equation 7,

From the first equation we obtain B_↑(↓) = 0, while the second equation results in the expression for A_↑(↓).

Now, we will discuss a ferromagnet with strong uniaxial anisotropy, in which case the perpendicular fluctuations of the exchange field are suppressed (α_x ∼ 0). For simplicity, we also assume the ferromagnet with weak spin–orbit interactions and also neglect the spin–orbit scattering time α_so. Finally, assuming α_x = α_so = 0 and keeping the solution to the lowest order, the equation for θ_f takes the form

where sinθ_s = Δ₀ / and cosθ_s = ω_n/ Here, Δ₀ is the bulk value of the pair potential. The equation above can be used for further semi-analytical calculations of the DOS for the case of a thin F layer with the magnetic scattering rate α_z. When α_z = 0, Equation 16 reduces to the well-known result (see, for example, [81] or [86]),

Analytical solution in the low proximity limit and small F layer thickness

In this section, we perform further analytical calculations of the anomalous Green’s function in the F layer based on the results of the previous section. We then analyze the effect of various scattering rates on the superconducting correlations, including the odd-frequency triplet component, which is generated in the adjacent ferromagnet.

The expression for θ-parameterized Green’s functions can be found from the boundary conditions (Equation 14),

In order to simplify the calculation and the final form of the solution θ_f, we consider only positive Matsubara frequencies ω_n and perform further linearization of Equation 18, which is justified in the low-proximity limit. Then, we obtain

Here, Σα_i denotes the sum of all the scattering rates. The above solution is true for thin ferromagnetic layers in the low-proximity limit. A more general analytical solution can be obtained for arbitrary thicknesses using the linearized Usadel equations. In order to find the DOS in the proposed limit, we expand the θ-parametrized normal Green’s function around small values of θ_f. In this case, we have

and, to calculate the total DOS, we need to sum the contributions from two spin populations using Equation 9.

Results and Discussion

In the present section we outline the main results, including both numerical and analytical calculations. The following parameters are fixed throughout the section: T = 0.1T_c, γ = 0.05, and d_f = 0.5ξ_f. First, we discuss general features of the DOS in an SF bilayer in the absence of any scattering. Then the effect of the spin-dependent scattering on the key DOS features in two relevant cases is discussed (see below) and finally, we present the analytical result and compare it with the numerically calculated DOS.

Evolution of the DOS in SF bilayer

It is instructive to discuss the key features of the DOS in an SF bilayer first. That is why, in this section, we briefly discuss the evolution of the DOS for different values of the exchange field h and the barrier transparency γ_B. All the scattering is assumed to be absent for simplicity α_m = α_x = α_so = 0 in this subsection.

In Figure 2, we observe the influence of an increasing exchange field h on the DOS structure calculated for γ_B = 5. In particular, we can see the evolution of the DOS peaks. For h = 0, that is, for the case of an SN bilayer, we see the well-known DOS structure with the characteristic mini-gap at energies E < Δ (Figure 2a, black dotted line) [79]. This proximity-induced mini-gap originates from the effective backscattering of the quasiparticles at the SN interface due to a finite interface resistance [99]. As h increases, the DOS splits for the spin-up and spin-down electrons, which results in the mini-gap peak splitting. For a certain value of h, the mini-gap closes, resulting in the DOS enhancement at zero energy as seen from Figure 2b and Figure 2c. This feature known as a zero-energy peak (ZEP) has been investigated both theoretically [100-103] and experimentally [49]. Another interesting peculiarity of the DOS is the appearance of the characteristic peak at E = h, which arises as the exchange field exceeds the superconducting gap h > Δ. Apparently, this peak arises from the evolution of the second spin-split peak due to a non-zero exchange field. The existence of such an effect offers a method of determining relatively small exchange field values in the F layer via DOS measurements [48,51,104].

In Figure 3, the DOS evolution at increasing interface parameter γ_B is shown. The blue solid line corresponds to the subgap exchange field h = 0.4Δ, whereas the black dotted line corresponds to h = 1.5Δ. From the figure we can notice that an increase of γ_B also has a strong influence on the DOS structure. Sufficiently large interface resistance values can close the mini-gap and lead to the emergence of the ZEP (Figure 3c, blue solid line). However, the peak structure is different for the two exchange fields h as seen from the figure.

One can notice that the peak at E = h in Figure 3a is almost absent. The reason that we observe such a behavior is because we consider relatively small values of F layer thickness d_f and transparency parameter γ_B. More detailed analysis can be made in the limiting case of d_f ≫ ξ_f and γ_B ≫ 1 and γ = 0 [51]. In the absence of any scattering, the analytical DOS expression for E ≥ Δ can be written in the following way:

where ε = One can clearly see the exponential asymptotic of the peak at E = h from Equation 21. We should not forget that Equation 21 is strictly valid only for large d_f/ξ_f, but nevertheless it may qualitatively explain why we do not see the peak at E = h for small a ratio of d_f/ξ_f. If this factor is relatively small, the variation of the exponent {−2(d_f/ξ_f)} near the point E = h is also small. The peak is observable only for h of the order of a few Δ. For larger exchange fields, the peak is very difficult to observe since the energy-dependent prefactor of the exponent in Equation 21 decays as E⁻² for E ≫ Δ.

In what follows, we will examine the effect produced by both spin-flip and spin–orbit scattering on the DOS features, mostly focusing on the mini-gap and the DOS peak at E = h. Unlike previous results on this topic [84,85], we provide both numerical and analytical results for the DOS calculation. Although the analytical expressions have a rather narrow range of applicability such limiting cases are relevant for experiments.

Effect of scattering on the DOS features

Now we discuss the influence of the finite scattering rates on the DOS features mentioned in the previous section. In this paper we consider two cases of the junction transparency: (i) intermediate interface transparency (γ_B ≥ 1) and (ii) low interface transparency (γ_B ≫ 1). In both cases, we fix the thickness of the F layer to d_f = 0.5ξ_f. Focusing on these cases allows us to discuss all major effects on the DOS features utilizing not only numerical solutions of the problem but also some analytical results, which will be presented below.

Intermediate interface transparency (γ_B = 5)

Figure 4 depicts the DOS dependencies in the case of relatively low interface transparency (γ_B = 5) in the presence of spin-flip and spin–orbit scattering. It is clearly seen that the decrease of the parallel magnetic scattering time leads to a smearing of the split peaks with the gradual closing of the induced energy gap (mini-gap) in the F layer (Figure 4a). The influence of the perpendicular magnetic scattering can be observed in Figure 4c. While increasing the scattering rate α_x tends to suppress the split peaks, perpendicular magnetic scattering also moves the peaks towards the Fermi energy destroying the mini-gap. This can be explained in terms of the additional effective exchange field in the system, which is due to non-zero perpendicular scattering rates [84]. Summarizing the results of the calculations, it is obvious that the magnetic scattering tends to destroy the proximity-induced superconductivity in the F layer. Such an effect becomes clear from a more detailed analysis of the linearized Usadel equation (Equation 2) in the low-proximity limit. In this case, the anomalous Green’s function is dependent on the exchange field h and the magnetic scattering rates, which apparently are pair breaking.

Figure 4b shows that a smaller spin–orbit scattering time leads to the vanishing of the peak splitting in the subgap region. In contrast, though the spin–orbit scattering destroys the double peak structure due to an exchange field smearing them into one peak, it does not produce a destructive effect on the mini-gap magnitude (see below Figure 6b). This feature has been reported previously [85].

Low interface transparency (γ_B = 50)

Now we focus on the limit of a highly resistive SF interface and investigate the effects of a spin-dependent scattering. As expected, the influence of the adjacent superconducting layer on the DOS is rather limited due to low transparency of the interface (Figure 5). It can be seen that the mini-gap is hardly recognizable in the cases of both subgap values of h (Figure 5a,b) and h > Δ (Figure 5c). From the plots, we can say that the finite scattering rates suppress the DOS features in the case of low interface transparency as well (Figure 5, solid and dashed lines). Even the DOS peak at E = h is suppressed substantially (Figure 5c). However, a closer examination shows that all the scattering rates slightly differ in the way they modify the DOS structure.

We would like to discuss the effect of scattering on the DOS peak located at the exchange energy in more detail. As we have mentioned above, one of the interesting features in the DOS of the considered system is the peak at E = h (Figure 6a). In Figure 6, we demonstrate the influence of different scattering rates on the DOS peak at E = h using the numerically obtained results. The remaining parameters used for calculations here are γ_B = 50, d_f = 0.5ξ_f, and h = 1.5Δ. In Figure 6b, the plot for different values of α_z is shown. It can be noticed that the uniaxial magnetic scattering not only suppresses the peak but also slightly shifts the DOS peak towards E = 0. The spin–orbit scattering has a similar effect on the peak, though α_so has a stronger effect on the peak height compared to α_z as it can be noticed from Figure 6c. In both cases above, the DOS peak also smears as any of the scattering rates increases. The effect of the perpendicular magnetic scattering α_x is indicated in Figure 6c.

Analytical result for the interfaces with low transparency and qualitative picture

Here, we employ the analytical expression (Equation 20) obtained in the limit of low proximity and thin F layer. Considering the problem in such limit makes it possible to use a simple expression for the qualitative description of the corresponding scattering effects on the DOS structure. It should not be forgotten that the linearized solution of the form in Equation 20 is quite limited in its application. In our case, it is valid when γ_B ≫ 1 and d_f ≪ min(ξ_f,), which is true for γ_B = 50. This tunneling limit is experimentally feasible. Thus, our result could easily be applied.

In Figure 7, the DOS calculated analytically via Equation 20 is illustrated. Here, we focus on the case of zero exchange field h = 0 to investigate the impact of each type of scattering on the mini-gap. We plot the analytically obtained DOS for the SN case in the absence of any scattering for comparison (black dotted line) as well. From the figure, one can see that the spin–orbit scattering does not affect the DOS in any way (Figure 7b). This effect has been shown before numerically in [85]. In contrast, both nonzero types magnetic of scattering, α_x and α_z, have a strong effect on the mini-gap, leading to its complete vanishing at some value of α (Figure 7a,c).

Making comparisons with previous results [84,85], we can say that there is a qualitative agreement in the DOS behavior. In the case of large DOS variations and, especially, singularities, the analytical model introduced above may fail. Nevertheless, we can explain major features of N_f(E) in the presence of a spin-dependent scattering. Examining the linearized Usadel equations, we can analyze the anomalous Green’s functions by studying even-frequency spin-singlet (f_s ∝ (θ_f↑ + θ_f↓)/2) and odd-frequency spin-triplet (f_t ∝ (θ_f↑ − θ_f↓)/2) components. When there is an F layer with a relatively high rate of parallel magnetic scattering, we can simplify the linearized Usadel equation and obtain

From this equation, we can find that θ_f↑ = θ_f↓ = 0, leading to supression of both the singlet and triplet components. In the presence of large in-plane magnetic scattering α_x, we obtain θ_f↑ = −θ_f↓, which leads to f_s = 0, whereas the triplet component f_t is nonzero. This can be understood in a similar way from the linearized Usadel equation

Since the energy gap is defined by the singlet correlations, we observe a detrimental effect of the in-plane scattering on the mini-gap magnitude (Figure 7c). In contrast, in the limit of strong spin–orbit scattering α_so, the Usadel equation in F layer reads

which results in θ_f↑ = θ_f↓, causing strong suppression of the triplet component and not the singlet one, which in turn explains the robustness of the mini-gap (Figure 7b). It can be demonstrated that the suppression of the triplet component due to spin–orbit scattering can actually lead to the appearance of the mini-gap at finite exchange fields h.

One of the possible ways of observing the DOS features is the examination of the current–voltage characteristics. Utilizing the Werthamer expression for the quasiparticle current in tunneling junctions, we can calculate the I–V curves for an SFIFS junction. The current then reads

Here, N_f1,2(E) is the density of states (DOS) in the corresponding ferromagnetic layer at x = 0, f(E) = [1 + e^E/T]⁻¹ is the Fermi–Dirac distribution function, and R = R_B0 is the resistance across the FIF interface. Both densities of states N_f1,2(E) are normalized to their values in the normal state. The abovementioned effects of spin-dependent scattering have a direct influence on the current. Figure 8 demonstrates the current–voltage characteristics of the SFIFS junction calculated in the presence of parallel magnetic (Figure 8a), spin–orbit (Figure 8b), and perpendicular magnetic scattering (Figure 8c). From the plots, we can notice that while a magnetic scattering destroys the mini-gap, the spin–orbit scattering slightly enhances it.

Finally, we compare the analytically derived and numerically calculated DOS in the case of an SF junction with a thin F layer and low-transparency interface. The corresponding result is shown in Figure 9. We can observe a fairly good agreement between the numerical and the analytical calculations. As expected, the analytical expression in Equation 20 cannot describe the features of N_f(E), which are relatively large in scale compared to unity.

Conclusion

We have formulated a model that takes into account magnetic and spin–orbit scattering processes in the framework of the quasiclassical Green’s function approach in the diffusive limit. Based on these equations, the local density of states has been calculated numerically. Applying the developed numerical solution, we have studied some previously overlooked features such as the influence of the scattering rates on the peak at E = h. Moreover, we provide a relatively simple expression to calculate the DOS analytically in the presence of magnetic scattering α_z for thin F layers. In addition, the analytic solution for the anomalous Green’s function has been derived in the limit of low proximity and a thin ferromagnetic layer. Based on this solution, we have been able to present analytical results for the DOS taking into account all spin-dependent scattering. We have demonstrated that the analytical result is in qualitative agreement with the numerical predictions, including previously published findings.