The heterostructures consist of ferromagnets and superconductors (FS) have been a focus of nanoelectronic research in recent years because of their relevance to spintronics applications as well as their potential impact on fundamental research^[1-2]. There are many theoretical works for various structures, such as diffusive FS^[3-5], clean FS^[6-8], F/unconventional-superconductor^[9], and nonequilibrium FS junctions^[10]. Many of these studies have been the generation and detection of long-range superconducting correlations in the ferromagnet. For a simple ferromagnet/spin-singlet superconductor junction, the proximity effect is often limited to be short-range due to the ferromagnetic exchange field acting differently on the spins of the Cooper pairs. This limitation would not arise if the Cooper pairs in the superconductor had spin-triplet symmetry. Cooper pairs in the equal-spin state are immune to the exchange field, hence they can propagate over a long distance into the ferromagnet.

It was predicted by Bergeret et al.^[2] and Kadigrobov et al.^[11] that spin-triplet pair correlations can be generated in an FS system involving inhomogeneous magnetizations, even when all superconductors present in the system have conventional spin-singlet symmetry. By now the predictions of these authors have been verified in a wide variety of different experiments; for reviews, see Refs.[1, 12-13]. In practice, instead of considering inhomogeneous profiles within a ferromagnetic layer as sources of long-range spin-triplet correlations, one can obtain the same effect in a multilayered geometry with non-collinear arrangements of the magnetizations between adjacent layers. An F₁F₂S trilayer is one of the simplest geometries of these structures, which is usually studied as a spin valve^[14], utilizing both charge and spin degrees of freedom. When long-range equal-spin pairs are achieved in this junction, the Fermi-Dirac statistics require them to be odd in time^[6] or frequency^[15]. The signatures of the DOS in the F layers of these hybrid structures can reveal the existence and type of superconducting correlations in that regions^[16-17]. There are many theories that stated odd-frequency pairs would cause an enhancement of the zero-energy DOS^[18-23], which is usually manifested as the ZEP.

The authors of Ref. [20] have discussed the link between the equal-spin triplet correlations and the ZEPs in the DOS, assuming F₁=F₂ in an F₁F₂S structure. However, in the half-metallic limit of the ferromagnet, they failed to find the ZEP, although this limit is most favorable for the spin-triplet correlation. In Ref. [21] the ZEP is calculated confined to the weak ferromagnetic limit; the height of the ZEP seems too low to be probed experimentally. In this work, we use the well-established Green's function method to analyze an F₁F₂S junction in the clean limit, with no restriction of the magnitudes or directions of the ferromagnets. Our primary aim is to find the interdependence of the long-range triplet correlation and ZEP. We first calculate the pair correlations in the F₁ layer and discuss how to create pure spin-triplet correlations. Then we calculate the ZEP in DOS, due to the zero-energy Andreev bound state. Finally, we discuss the connection between this ZEP and the long-range triplet correlation.

1 Formulation

Let us consider a ferromagnet/ferromagnet/s-wave superconductor (F₁F₂S) heterojunction as shown in the clean limit in Fig. 1. The transport direction is along the x-axis. The two interfaces locating at x=0 and x=l are described by δ-type potentials. The direction of the exchange ferromagnetic field in each ferromagnet layer is specified by their polar and azimuthal angles. The whole junction can be described by the Bogoliubov-de Gennes (BdG) equation $\hat{{\boldsymbol{H}}} {\boldsymbol{\varPsi}}=E {\boldsymbol{\varPsi}}$,

where

$ \hat{{\boldsymbol{H}}}=\left(\begin{array}{cc} \hat{{\boldsymbol{H}}}_0({\boldsymbol{r}}) & \hat{{\boldsymbol{\varDelta}}}({\boldsymbol{r}}) \\ \hat{{\boldsymbol{\varDelta}}}^{\dagger}({\boldsymbol{r}}) & -\hat{{\boldsymbol{H}}}_0^{{\mathrm{T}}}({\boldsymbol{r}}) \end{array}\right), $

(1)

with

$ \begin{aligned} \hat{{\boldsymbol{H}}}_0({\boldsymbol{r}})= & {\left[-\frac{\hbar^2 \nabla^2}{2 m}-E_F+U_1 {\mathtt{δ}}(x)+U_2 {\mathtt{δ}}(x-l)\right] \hat{1}-\hat{{\boldsymbol{\sigma}}} \cdot } \\ & \left({\boldsymbol{h}}_1 {\boldsymbol{\varTheta}}(-x)+{\boldsymbol{h}}_2 {\boldsymbol{\varTheta}}(x) {\boldsymbol{\varTheta}}(l-x)\right), \end{aligned} $

(2)

$ \begin{gathered} {\boldsymbol{h}}_i=h_i\left[\sin \alpha_i \cos \beta_i {\boldsymbol{e}}_x+\sin \alpha_i \sin \beta_i {\boldsymbol{e}}_y+\cos \alpha_i {\boldsymbol{e}}_z\right], \\ (i=1, 2) \end{gathered} $

(3)

and

$ \hat{{\boldsymbol{\Delta}}}({\boldsymbol{r}})=i \varDelta(T) \hat{{\boldsymbol{\sigma}}}_y {\boldsymbol{\varTheta}}(x-l). $

(4)

Here $\hat{{\boldsymbol{\sigma}}}_y=\left(\begin{array}{cc} 0 & -i \\ i & 0 \end{array}\right)$, Θ(x) is the Heavyside step function, and the temperature dependence of the energy gap is given by Δ(T)=Δ₀tanh ($1.74 \sqrt{T_c / T-1}$) with T_c the superconducting critical temperature.

After solving the BdG equation in each region, we can write down the right-moving particle (quasiparticle) wavefunctions as follows:

$ {\boldsymbol{\varphi}}_{F_{j}^e \uparrow, \downarrow}(x)={\mathrm{e}}^{{\mathrm{i}} {\boldsymbol{p}} \cdot {\boldsymbol{\rho}}} {\mathrm{e}}^{{\mathrm{i}} k_{j \uparrow, \downarrow}^e x} \hat{{\boldsymbol{F}}}_{j(1, 2)}, \\{\boldsymbol{\varphi}}_{F_{j}^ h \uparrow, \downarrow}(x)={\mathrm{e}}^{{\mathrm{i}} {\boldsymbol{p}} \cdot {\boldsymbol{\rho}}} {\mathrm{e}}^{{\mathrm{i}} k_{j \uparrow, \downarrow}^h (-x)} \hat{{\boldsymbol{F}}}_{j(3, 4)}, \\{\boldsymbol{\varphi}}_{S e \uparrow, \downarrow}(x)={\mathrm{e}}^{{\mathrm{i}} {\boldsymbol{p}} \cdot {\boldsymbol{\rho}}} {\mathrm{e}}^{{\mathrm{i}} q_{\uparrow, \downarrow}^e x} \hat{{\boldsymbol{S}}}_{1, 2}, \\{\boldsymbol{\varphi}}_{S h \uparrow, \downarrow}(x)={\mathrm{e}}^{{\mathrm{i}} {\boldsymbol{p}} \cdot {\boldsymbol{\rho}}} {\mathrm{e}}^{{\mathrm{i}} q_{\uparrow, \downarrow}^h(-x)} \hat{{\boldsymbol{S}}}_{3, 4} . $

(5)

Here k_{j↑, ↓}^e= $\sqrt{\frac{2 m}{\hbar^2}\left(E_F-{\boldsymbol{\epsilon}}_p+E \pm h_j\right)}$, k_{j↑, ↓}^h= $\sqrt{\frac{2 m}{\hbar^2}\left(E_F-{\boldsymbol{\epsilon}}_p-E \pm h_j\right)}$, (j=1, 2), are the wave vectors of the electron-like and hole-like quasiparticles in the ferromagnets, and q_{↑, ↓}^e= $\sqrt{\frac{2 m}{\hbar^2}\left(E_F-{\boldsymbol{\epsilon}}_p+\varOmega\right)}=q^{+}$, q_{↑, ↓}^h= $\sqrt{\frac{2 m}{\hbar^2}\left(E_F-{\boldsymbol{\epsilon}}_p-\varOmega\right)}$ =q^- with $\varOmega=\sqrt{E^2-\Delta^2}$ are the wave vectors in the superconductors. The quantities p =(k_y, k_z), ρ =(y, z), and ${\boldsymbol{\epsilon}}_p=\frac{{\boldsymbol{p}}^2}{2 m}$ are the components parallel to the interfaces.

The basis vectors of the wavefunctions for ferromagnets ($\hat{{\boldsymbol{F}}}_{i j}$) are given by:

$ \hat{{\boldsymbol{F}}}_{i 1}=\left(\begin{array}{c} \cos \frac{\alpha_i}{2} {\mathrm{e}}^{-\frac{{{{\rm{i}}}} \beta_i}{2}} \\ \sin \frac{\alpha_i}{2} {\mathrm{e}}^{\frac{{{{\rm{i}}}} \beta_i}{2}} \\ 0 \\ 0 \end{array}\right), \hat{{\boldsymbol{F}}}_{i 2}=\left(\begin{array}{c} -\sin \frac{\alpha_i}{2} {\mathrm{e}}^{-\frac{{{{\rm{i}}}} \beta_i}{2}} \\ \cos \frac{\alpha_i}{2} {\mathrm{e}}^{\frac{{{{\rm{i}}}} \beta_i}{2}} \\ 0 \\ 0 \end{array}\right), \\\hat{{\boldsymbol{F}}}_{i 3}=\left(\begin{array}{c} 0 \\ 0 \\ \cos \frac{\alpha_i}{2} {\mathrm{e}}^{\frac{{{{\rm{i}}}} \beta_i}{2}} \\ {{{\rm{si}}}}n \frac{\alpha_i}{2} {\mathrm{e}}^{-\frac{{{{\rm{i}}}} \beta_i}{2}} \end{array}\right), \hat{{\boldsymbol{F}}}_{i 4}=\left(\begin{array}{c} 0 \\ 0 \\ -\sin \frac{\alpha_i}{2} {\mathrm{e}}^{\frac{{{{\rm{i}}}} \beta_i}{2}} \\ \cos \frac{\alpha_i}{2} {\mathrm{e}}^{-\frac{{{{\rm{i}}}} \beta_i}{2}} \end{array}\right), \\(i=1, 2) . $

(6)

and those for superconductors ($\hat{{\boldsymbol{S}}}_i$) are given by:

$ \hat{{\boldsymbol{S}}}_1=\left(\begin{array}{l} u \\ 0 \\ 0 \\ v \end{array}\right), \hat{{\boldsymbol{S}}}_2=\left(\begin{array}{c} 0 \\ u \\ -v \\ 0 \end{array}\right), \\\hat{{\boldsymbol{S}}}_3=\left(\begin{array}{c} 0 \\ -v \\ u \\ 0 \end{array}\right), \hat{{\boldsymbol{S}}}_4=\left(\begin{array}{l} v \\ 0 \\ 0 \\ u \end{array}\right) . $

(7)

Here $u=\sqrt{(1+\varOmega / E) / 2}$, $v=\sqrt{(1-\varOmega / E) / 2}$. To construct the Green's functions of the two outside superconductor leads, we need to consider eight types of elementary processes: the scattering of a spin up (down) electron (hole)-like quasiparticle incident from the left (right) superconductor. For example, as shown in Fig. 2, the wavefunction (in each region) describing the scattering of a spin-up electron-like quasiparticle incident from the left is given by:

$ \begin{gathered} {\boldsymbol{\varPsi}}_{1 F_1}(x)={\boldsymbol{\varphi}}_{F_1 e^{\uparrow}}(x)+a_1 {\boldsymbol{\varphi}}_{F_1 h \downarrow}(-x)+ \\ b_1 {\boldsymbol{\varphi}}_{F_{1}e \uparrow}(-x)+a_1^{\prime} {\boldsymbol{\varphi}}_{F_1 h \uparrow}(-x)+b_1^{\prime} {\boldsymbol{\varphi}}_{F_1 e \downarrow}(-x), \\ {\boldsymbol{\varPsi}}_{1 F_2}(x)=e_1 {\boldsymbol{\varphi}}_{F_2 e \uparrow}(x)+e_1^{\prime} {\boldsymbol{\varphi}}_{F_2 e \downarrow}(x)+ \\ f_1 {\boldsymbol{\varphi}}_{F_{2} e\uparrow}(-x)+f_1^{\prime} {\boldsymbol{\varphi}}_{F_2 e \downarrow}(-x)+g_1 {\boldsymbol{\varphi}}_{F_2 h \uparrow}(x)+ \\ g_1^{\prime} {\boldsymbol{\varphi}}_{F_2 h \downarrow}(x)+h_1 {\boldsymbol{\varphi}}_{F_2 h \uparrow}(-x)+h_1^{\prime} {\boldsymbol{\varphi}}_{F_2 h \downarrow}(-x), \\ {\boldsymbol{\varPsi}}_{1 S}(x)=c_1 {\boldsymbol{\varphi}}_{S e}(x)+d_1 {\boldsymbol{\varphi}}_{S h \downarrow}(x)+c_1^{\prime} {\boldsymbol{\varphi}}_{S e \downarrow}(x)+ \\ d_1^{\prime} {\boldsymbol{\varphi}}_{S h \uparrow}(x) . \end{gathered} $

(8)

Here, the coefficients a₁, b₁, a′₁, b′₁ describe the usual Andreev reflection, normal reflection, novel Andreev reflection, and the normal reflection with spin-flip, respectively. c₁, d₁, c′₁, d′₁, e₁, e′₁, f₁, f′₁, g₁, g′₁, h₁, h′₁ are the transmission coefficients, corresponding to the processes described in Fig. 2.

All the scattering coefficients can be determined by the matching conditions:

$ \begin{aligned} & {\boldsymbol{\varPsi}}_i\left(0^{-}\right)={\boldsymbol{\varPsi}}_i\left(0^{+}\right), {\boldsymbol{\varPsi}}_i\left(l^{-}\right)={\boldsymbol{\varPsi}}_i\left(l^{+}\right), \\ & \frac{\partial {\boldsymbol{\varPsi}}_i\left(0^{+}\right)}{\partial x}-\frac{\partial {\boldsymbol{\varPsi}}_i\left(0^{-}\right)}{\partial x}=k_F Z_1 {\boldsymbol{\varPsi}}_i(0), \\ & \frac{\partial {\boldsymbol{\varPsi}}_i\left(l^{+}\right)}{\partial x}-\frac{\partial {\boldsymbol{\varPsi}}_i\left(l^{-}\right)}{\partial x}=k_F Z_2 {\boldsymbol{\varPsi}}_i(l), i=1, 2, \cdots, 8 . \end{aligned} $

(9)

where Z_i=2mU_i/ħ²k_F, i=1, 2 are dimensionless parameters describing the magnitude of the interfacial resistance, and $k_F=\sqrt{\frac{2 m E_F}{\hbar^2}}$. The wavefunctions of the other seven types of processes can be obtained similarly.

The 4×4 matrix retarded Green's function is constructed from a linear combination of outer the outer products of the wavefunctions Ψ_i and their conjugates Ψ_i^†[24-28]:

$ \begin{array}{l} {\boldsymbol{G}}^r\left(x<x^{\prime} ; E, {\boldsymbol{\epsilon}}_p\right)=\sum\limits_{i=5}^8 \sum\limits_{j=1}^4 \alpha_{i j} {\boldsymbol{\varPsi}}_i(x) {\boldsymbol{\varPsi}}_j^{\dagger}\left(x^{\prime}\right), \\ {\boldsymbol{G}}^r\left(x<x^{\prime} ; E, {\boldsymbol{\epsilon}}_p\right)=\sum\limits_{j=1}^4 \sum\limits_{i=5}^8 \beta_{j i} {\boldsymbol{\varPsi}}_j(x) {\boldsymbol{\varPsi}}_i^{\dagger}\left(x^{\prime}\right). \end{array} $

(10)

The coefficients α and β are obtained from the conditions below:

$ \begin{gathered} {\boldsymbol{G}}^r\left(x, x+0^{+} ; E, {\boldsymbol{\epsilon}}_p\right)={\boldsymbol{G}}^r\left(x, x-0^{+} ; E, {\boldsymbol{\epsilon}}_p\right), \\ \left.\frac{\partial {\boldsymbol{G}}^r\left(x, x^{\prime} ; E, {\boldsymbol{\epsilon}}_p\right)}{\partial x}\right|_{x=x^{\prime}+0^{+}}-\left.\frac{\partial {\boldsymbol{G}}^r\left(x, x^{\prime} ; E, {\boldsymbol{\epsilon}}_p\right)}{\partial x}\right|_{x=x^{\prime}-0^{+}} \\ =\frac{2 m}{\hbar^2}\left(\begin{array}{cc} I & 0 \\ 0 & -{\boldsymbol{I}} \end{array}\right) . \end{gathered} $

(11)

We finally get the retarded Green's functions in F₁ and S:

$ \begin{array}{l} & {\boldsymbol{G}}_{F_1}^r\left(x, x^{\prime}\right)=\sum\limits_{i=1}^4 {\boldsymbol{G}}_i\left(x, x^{\prime}\right), \\ & {\boldsymbol{G}}_S^r\left(x, x^{\prime}\right)=\sum\limits_{i=5}^8 {\boldsymbol{G}}_i\left(x, x^{\prime}\right), \end{array} $

(12)

where

$ \begin{aligned} & {\boldsymbol{G}}_1\left(x, x^{\prime}\right)=-\frac{{\mathrm{im}} E}{k_{11} \varOmega_{F_1}} {\boldsymbol{\varPsi}}_{1 F_1}(x)\left({\boldsymbol{\varphi}}_{F_{1}e \uparrow}\left(x^{\prime}\right)\right)^{\dagger}, \\ & {\boldsymbol{G}}_2\left(x, x^{\prime}\right)=-\frac{{\mathrm{im}} E}{k_{12} \varOmega_{F_1}} {\boldsymbol{\varPsi}}_{2 F_1}(x)\left({\boldsymbol{\varphi}}_{F_{1}e \downarrow}\left(x^{\prime}\right)\right)^{\dagger}, \\ & {\boldsymbol{G}}_3\left(x, x^{\prime}\right)=-\frac{{\mathrm{im}} E}{k_{13} \varOmega_{F_1}} {\boldsymbol{\varPsi}}_{3 F_1}(x)\left({\boldsymbol{\varphi}}_{F_{1}h \uparrow}\left(x^{\prime}\right)\right)^{\dagger}, \\ & {\boldsymbol{G}}_4\left(x, x^{\prime}\right)=-\frac{{\mathrm{im}} E}{k_{14} \varOmega_{F_1}} {\boldsymbol{\varPsi}}_{4 F_1}(x)\left({\boldsymbol{\varphi}}_{F_{1}h \downarrow}\left(x^{\prime}\right)\right)^{\dagger}, \\& {\boldsymbol{G}}_5\left(x, x^{\prime}\right)=-\frac{\operatorname{im} E}{q^{+} \varOmega_S} {\boldsymbol{\varPsi}}_{5 S}(x)\left({\boldsymbol{\varphi}}_{S e\uparrow} \left(-x^{\prime}\right)\right)^{\dagger}, \\ & {\boldsymbol{G}}_6\left(x, x^{\prime}\right)=-\frac{\operatorname{im} E}{q^{+} \varOmega_S} {\boldsymbol{\varPsi}}_{6 S}(x)\left({\boldsymbol{\varphi}}_{S e \downarrow}\left(-x^{\prime}\right)\right)^{\dagger}, \\ & {\boldsymbol{G}}_7\left(x, x^{\prime}\right)=-\frac{\operatorname{im} E}{q^{-} \varOmega_S} {\boldsymbol{\varPsi}}_{7 S}(x)\left({\boldsymbol{\varphi}}_{S h \uparrow}\left(-x^{\prime}\right)\right)^{\dagger}, \\ & {\boldsymbol{G}}_8\left(x, x^{\prime}\right)=-\frac{\operatorname{im} E}{q^{-} \varOmega_S} {\boldsymbol{\varPsi}}_{8 S}(x)\left({\boldsymbol{\varphi}}_{S h \downarrow}\left(-x^{\prime}\right)\right)^{\dagger} . \end{aligned} $

(13)

and Ω_F1=E.

It can be checked that our Green's functions satisfy the Gor'kov equations $\{E-\hat{{\boldsymbol{H}}}(x)\}$ G^r(x, x′; E, ${\boldsymbol{\epsilon}}_p$)=δ(x-x′)$\hat{1}$. After an analytic continuation E→iω_n, we can also get the 4×4 matrix Matsubara Green's function $\mathcal{G}\left(x, x^{\prime} ; {\mathrm{i}} \omega_n, {\boldsymbol{\epsilon}}_p\right)$, where ω_n=(2n+1)πk_BT.

Here we introduce the spin-singlet pair correlation f_S and the spin-triplet pair correlations f₀, f₁, and f_-1 in terms of the anomalous Matsubara Green's functions as follows:

$ \begin{aligned} & f_S= \frac{1}{\sqrt{2}} k_B T \sum\limits_{{\boldsymbol{\epsilon}}_p} \sum\limits_n\left({\boldsymbol{\mathcal{G}}}_{14}\left(x, x ; {\mathrm{i}} \omega_n, {\boldsymbol{\epsilon}}_p\right)-\right. \\ &\left.{\boldsymbol{\mathcal{G}}}_{23}\left(x, x ; {\mathrm{i}} \omega_n, {\boldsymbol{\epsilon}}_p\right)\right), \\ & f_0=\frac{1}{\sqrt{2}} k_B T \sum\limits_{{\boldsymbol{\epsilon}}_p} \sum\limits_n\left({\boldsymbol{\mathcal{G}}}_{14}\left(x, x ; {\mathrm{i}} \omega_n, {\boldsymbol{\epsilon}}_p\right)+\right. \\ &\left.{\boldsymbol{\mathcal{G}}}_{23}\left(x, x ; {\mathrm{i}} \omega_n, {\boldsymbol{\epsilon}}_p\right)\right), \\ & f_1= k_B T \sum\limits_{{\boldsymbol{\epsilon}}_p} \sum\limits_n\left({\boldsymbol{\mathcal{G}}}_{13}\left(x, x ; {\mathrm{i}} \omega_n, {\boldsymbol{\epsilon}}_p\right)\right), \\ & f_{-1}=k_B T \sum\limits_{{\boldsymbol{\epsilon}}_p} \sum\limits_n\left({\boldsymbol{\mathcal{G}}}_{24}\left(x, x ; {\mathrm{i}} \omega_n, {\boldsymbol{\epsilon}}_p\right)\right) . \end{aligned} $

(14)

where f₀, f₁, and f_-1 are related to the probabilities for the spin projection being S_z=0, S_z=1, and S_z=-1, respectively. The local density of states is defined as:

$ \begin{gathered} {\boldsymbol{\rho}}(x, E)=-\sum\limits_{{\boldsymbol{\epsilon}}_p} \frac{1}{{\mathtt{π}}} \operatorname{Im}\left[{\boldsymbol{G}}_{11}^r\left(x, x ; E, {\boldsymbol{\epsilon}}_p\right)+\right. \\ \left.{\boldsymbol{G}}_{22}^r\left(x, x ; E, {\boldsymbol{\epsilon}}_p\right)\right]. \end{gathered} $

(15)

Throughout in the following numerical calculations, we take Δ(0)=0.005E_F, T/T_c=0.1, and Z₁=0.1, Z₂=0 assuming rather transparent interfaces.

2 Results and discussions

In this section, we describe our results obtained according to the formula given in the previous section. We shall confine ourselves to the properties of the left ferromagnet. We start with the discussion about the induced pair correlations and then present the results for the local density of states.

2.1 Pair correlations

In this subsection, we calculate the properties of various pair correlations induced in the F₁ region by proximity effect using the equations (14). The correlations in this work are all normalized to the value of the bulk spin-singlet superconductor. To begin with, let us consider the standard magnetization configuration for our F₁F₂S system where the exchange filed in the F₁ layer h₁ is fixed in the z-direction (α₁=0), while the exchange filed in the F₂ layer h₂ is allowed to rotate within the x-z plane (β₁, β₂=0). The relative angle between h₁ and h₂ (|α₂-α₁|) is designated by α in the following. When α =0, it is found that only the opposite-spin correlations Re(f_S) and Im(f₀) are induced in the F₁ layer. The spatial dependence of Re(f_S) and Im(f₀) for h₁/E_F=0.3, h₂/E_F=0.1 and lk_F=20 are depicted in Fig. 3(a). Correlations Re(f_S) and Im(f₀) exhibit a damped oscillatory behavior with an oscillating length ħv_F/h₁ where v_F denotes the Fermi velocity. Because they decay rapidly with increasing exchange fields, they are said to be short-ranged. Note that the correlations Re(f_S) and Im(f₀) can also be generated even if the F₁ layer is replaced by a normal metal one (NF₂S junction) as shown in the inset; they can penetrate deeply into the F₁ layer and therefore are called long-ranged.

Next, we turn to the case of α≠0. We take, as an example, α=0.5π, and h₁/E_F=0.3, h₂/E_F= 0.1, lk_F=20. Presented in Fig. 3 (b) are spatial variations for the induced pairing correlations. As we see, two new spin-triplet correlations Im(f₁) and Im(f_-1) are generated in addition to the spin-singlet correlation Re(f_S) and the spin-triplet correlation Im(f₀). The spin-triplet correlations Im(f₁) and Im(f_-1) decay rather slowly and monotonically similar to the NF₂S case discussed above. Such types of spatially long-range spin-triplet correlations can be created in the F₁ region when the exchange field h₁ is in the direction of the spin quantization axis and not collinear with h₂ at the same time. Plotted in Fig. 3(c) are the spatial behaviors of our correlations for h₁=0.98E_F with other system parameters taken the same as that in Fig. 3(b). Due to the strong pair-breaking effect of ferromagnetic exchange splitting, the amplitudes of f_S, f₀ as well as that of f_-1 are significantly reduced in comparison with the case of h₁/E_F= 0.3 (see Fig. 3(b)). As we see, however, the effect on the equal-spin component Im(f₁) is rather limited. We further consider the half-metallic case h₁/E_F=1, corresponding to the situation where only the up-spin electron energy band exists in this region. In this case, though the three correlations Re(f_S), Im(f₀), and Im(f_-1) are completely destroyed, the spin-triplet correlation Im(f₁) is still survived with its magnitude reduced slightly. Here we obtain a special F₁F₂S system where only the long-range and highly spin-polarized triplet pair correlation Im(f₁) is generated in F₁ region by superconducting proximity effect.

To further explore the influence of the exchange field h₁ on various pair correlations, we consider the correlations fixed at xk_F=-100 as functions of h₁ (see Fig. 4(a)). When h₁/E_F is increased from 0 to 1, correlation Im(f₁) changes slightly, and the correlations Re(f_S) and Im(f₀) undergo similar quickly damping oscillation. It is noted however that, despite the presence of the pair breaking effect, the f_-1 amplitude is increased up to h₁/E_F≈0.5, and its magnitude is comparable with that of f₁ over a wide range of exchange fields, except in the very vicinity of h₁/E_F=1. In the previous discussions, the relative exchange field angle α is mostly taken as 0.5π. As mentioned above, the equal-spin correlations Im(f₁) and Im(f_-1) can be generated in F₁ region when the exchange fields h₁ and h₂ are not collinear. Now we study their α dependencies more closely. Illustrated in Fig. 4(b) are the four correlations as functions of α at xk_F= -100, with h₁/E_F=0.8 and lk_F=20. As we see, the amplitudes of f_-1 and f₁ depend sensitively on α; they are zero when α=0 or α=π, and largest near α=0.5π. On the contrary, the amplitudes of f_S and f₀ undergo little changes as angle α varies from 0 to π. We obtain that when the exchange field reaches the half-metallic limit h₁/E_F=1, all pair correlations in Fig. 4(b) but Im(f₁) disappear (see inset).

Before closing this subsection, let us discuss an interesting problem that how the results obtained above are altered if the spin quantization axis is taken differently or if the exchange field h₁ is not collinear with the spin quantization axis with a relative angle β (α remains fixed). We present in Fig. 5(a) the situation where α=0, and β=0.25π. As we see, two more pair components Im(f_-1) and Im(f₁) are induced in the F₁ region in comparison to the case of α =0, and β =0 (see Fig. 3(a)), and moreover, all triplet correlations are short-ranged. As another example we show in Fig. 5(b) the result for h₁/E_F=1, α=0.5π, and β=0.5π which is compared with the result for h₁/E_F=1, α=0.5π, β=0 discussed above (see Fig. 3(d)). As can be seen, the triplet pair correlations f₀ and f_-1 have become long-ranged and survived in the F₁ region along with the f₁ correlation.

It is found that, however such seemingly quite different results can be boiled down to a unified description. We introduce a new set of pair correlations:

$ \begin{array}{l} \tilde{f}_S=\frac{1}{2} k_B T \sum\limits_{{\boldsymbol{\epsilon}}_p} \sum\limits_n \mid {\boldsymbol{\mathcal{G}}}_{14}\left(x, x ; {\mathrm{i}} \omega_n, {\boldsymbol{\epsilon}}_p\right)- \\ \left.{\boldsymbol{\mathcal{G}}}_{23}\left(x, x ; {\mathrm{i}} \omega_n, {\boldsymbol{\epsilon}}_p\right)\right|^2, \end{array} $

(16)

for spin-singlet correlation, and

$ \begin{array}{l} \tilde{f}_T=k_B T \sum\limits_{{\boldsymbol{\epsilon}}_p} \sum\limits_n\left(\frac{1}{2} \mid {\boldsymbol{\mathcal{G}}}_{14}\left(x, x ; {\mathrm{i}} \omega_n, {\boldsymbol{\epsilon}}_p\right)+\right. \\ \left.{\boldsymbol{\mathcal{G}}}_{23}\left(x, x ; {\mathrm{i}} \omega_n, {\boldsymbol{\epsilon}}_p\right)\right|^2+\left|{\boldsymbol{\mathcal{G}}}_{13}\left(x, x ; {\mathrm{i}} \omega_n, {\boldsymbol{\epsilon}}_p\right)\right|^2+ \\ \left.\left|{\boldsymbol{\mathcal{G}}}_{24}\left(x, x ; {\mathrm{i}} \omega_n, {\boldsymbol{\epsilon}}_p\right)\right|^2\right), \end{array} $

(17)

for spin-triplet correlation. We can verify that correlation $\tilde{f}_S$ and $\tilde{f}_T$ are independent of the relative angle β as well as the choice of spin quantization axis^[29]; their angular dependencies are only through the angle α, just as in the case of the DOS (see below).

2.2 Local density of states

After the discussion of the prominent features of various pair correlations induced within the F₁ region, we now turn to another topic of this paper: the local density of states measured in the F₁ region. The local density of states in our calculation is normalized by DOS(x, E)=ρ(x, E)/ρ(x, E)|Δ=0. The signatures of the proximity-induced electronic density of states in the F₁ layer can reveal the existence and type of superconducting correlations in this region. Practically, the presence of a zero-energy peak in the single-particle DOS spectrum is considered to be strong evidence for the existence of spin-triplet correlation in the F₁F₂S system^[20-21].

Here we calculate first the DOS for this salient case keeping the possible ZEP in mind. As a first step, we calculate the zero-energy DOS (DOS(E=0)) as a function of the width of the middle ferromagnet lk_F (0 < lk_F < 120). The results for α=0.5π, at xk_F=-20, with various values of h₁ and h₂ are illustrated in Fig. 6. The DOS(E=0)s display oscillatory behavior with a certain period, implying the formation of the zero-energy Andreev bound states in our F₁F₂S quantum well system^[30-31]. The magnitude of the peaks of DOS(E=0), which appears at some critical values of lk_F, as well as the periods of the oscillations, depend sensitively on the exchange fields h₁/E_F and h₂/E_F. The results in Fig. 6(a) and 6(b) show that for a given value of h₁/E_F both the peak and the period of the oscillation increase with decreasing the exchange field h₂/E_F. On the other hand, the peak reaches its maximum when h₁/E_F=1 for the systems with the same value of h₂/E_F (see Fig. 6(a), 6(c), and 6(d)).

To understand the origin of these zero-energy DOS patterns obtained numerically, further analytical study is required. Following the standard approach of the Green's function method, we seek analytically the poles of our one-particle retarded Green's function in the F₂ region for E=0. We focus on the α=0.5π case and assume Z₁=Z₂=0 for the sake of simplicity. In this situation, we need to solve the scattering problem described in Fig. 2 in the absence of the left incident process. Upon putting the determinant of the 16×16 matrix for the coefficients to be zero, we arrive at the following conditions which should be satisfied for the formation of the zero-energy Andreev bound state:

$ \begin{gathered} l k_F\left(\sqrt{\frac{\left(h_1+h_2\right)}{E_F}}-\sqrt{\left.\frac{\left(h_1-h_2\right)}{E_F}\right)}=\frac{{\mathtt{π}}}{2}(1+2 n), \right. \\ n=0, 1, 2, \cdots \end{gathered} $

(18)

In deriving the above result, we have used the limit h₁/h₂≫1 for the wave vectors, except for that in the exponentials^[32]. The critical widths lk_F calculated from the equation (18) are marked on the horizontal axis in Fig. 6. As can be seen, the agreement between the numerical and analytical results is excellent. In Fig. 7 we present, as a supplement to Fig. 6, the DOS(E=0) versus lk_F curves for some moderate values of h₂/E_F. For h₂/E_F=0.2, the DOS(E=0) displays eight resonance peaks with large average height (~ 6) (see Fig. 7(a)). In the case of h₂/E_F=0.4, we encounter 16 sharp peaks; although the average height is diminished drastically in comparison with the former case, it is still considerably large (~2) (see Fig. 7(b)). The critical widths calculated from equation (18) are found to be in good agreement with that obtained in both cases.

In Fig. 6 and Fig. 7 the peculiar peak structures of the DOS(E=0) are discussed only for α=0.5π. The equation (18) is also derived assuming α=0.5π. Now we turn to the α-dependence of the DOS(E=0). In Fig. 8(a) we show the DOS(E=0) as a function of the α for h₁/E_F=1, h₂/E_F=0.1, and lk_F=15.7 at xk_F=20; it is tightly packed in the vicinity of 0.5π. Shown in Fig. 8(b) are the DOS versus E curves for α=0, 0.4π and 0.5π. We note in particular that no zero-energy peak appears for α=0; the same feature is found for the cases of h₁/E_F < 1 in general. As our results show, both the correlation Im(f₁) and the ZEP get their maximum values when h₁/E_F=1, α~0.5π, and vanish together when α=0 irrespective of the value of the h₁/E_F, indicating an intimate correlation between the equal-spin triplet pair amplitude and the zero-energy Andreev bound state. By the way, we would like to point out that, in the study of Ref. [20] no zero-energy Andreev bound state is considered and the ZEP in their DOS spectrum gets its minimum the half-metallic limit is taken.

Finally, let us discuss the long-range nature of the ZEP. We return to the case of α=0.5π, h₁/E_F=1, and h₂/E_F=0.1.

The dependence of the DOS(E=0) on the width of F₂ layer lk_F at xk_F=-20 is already depicted in Fig. 6 (b), but here we consider the xk_F=-200 case (see Fig. 8(c)). As can be seen, although the heights of peaks are declined substantially the pattern as a whole is almost unchanged. Illustrated in Fig. 8(d) is the evolution of the amplitude of the first peak (corresponding to lk_F≈15.7) as a function of xk_F. The DOS at xk_F=-500 is plotted in the inset; the ZEP is still prominent. Clearly, similar to the triplet correlation f₁, the ZEP structure of DOS is long-ranged.

3 Summary

In summary, we have systematically studied the pair correlations and the local DOS induced an F₁F₂S heterostructure in the clean limit. The analytical and numerical results are obtained within the framework of the standard Green's function method. We found that different pair correlations may be induced in the F₁ region depending on the value of the relative angle α and the magnitude of the exchange field in the F₁ region h₁/E_F. In particular, the long-range triplet pair correlation Im(f₁) is generated only when α≠0. This amplitude of Im(f₁) is maximal when α~0.5π, and more importantly, it may become the only pair correlation present if h₁/E_F=1. In the numerical calculation of the one-particle DOS(E=0), we found an oscillatory multi-peaked behavior which is attributed to the formation of the zero-energy Andreev bound state. The condition for the occurrence of the zero-energy Andreev bound state is derived analytically for α=0.5π; this condition may be useful in the superconducting electronic device design and application in the future. Due to the zero-energy Andreev bound state, the DOS shows sharp ZEP. It is found that the ZEP structure is pronounced mostly when h₁/E_F=1 and α~0.5π and also of long-range. Therefore, a ZEP structure that appears deep in the strong ferromagnetic F₁ layer may be considered to be an unequivocal signature of the emergence of the long-range triplet pair correlation in our F₁F₂S system.