Charge and Spin Currents in Ferromagnetic Josephson junctions

# Charge and Spin Currents in Ferromagnetic Josephson junctions

Klaus Halterman Michelson Lab, Physics Division, Naval Air Warfare Center, China Lake, California 93555    Oriol T. Valls Also at Minnesota Supercomputer Institute, University of Minnesota, Minneapolis, Minnesota 55455 School of Physics and Astronomy, University of Minnesota, Minneapolis, Minnesota 55455    Chien-Te Wu School of Physics and Astronomy, University of Minnesota, Minneapolis, Minnesota 55455 The James Franck Institute, The University of Chicago, Chicago, Illinois, 60637
July 20, 2019
###### Abstract

We determine, using a self consistent method, the charge and spin currents in ballistic Josephson junctions consisting of several ferromagnetic () layers sandwiched between superconducting () electrodes (-type junctions). When there are two layers, we also consider the experimentally relevant configuration where a normal () nonmagnetic spacer layer separates them. We study the current-phase relationships as functions of geometrical parameters that are accessible experimentally including particularly the angles that characterize the relative orientation of the magnetization in the layers. Our self-consistent method ensures that the proper charge conservation laws are satisfied, and that important proximity effects are fully and properly accounted for. We find that as we vary the phase difference between the two outer electrodes, multiple harmonics in the current phase relations emerge, the extent of which depends on the interface scattering strength and on the relative layer widths and magnetization orientations. By manipulating the relative layer magnetization orientations, we find that the charge supercurrent can reverse directions or vanish altogether. These findings are discussed in the context of the generation and long-range nature of triplet pair correlation within these structures. We also investigate the spin currents and associated spin transfer torques throughout the entire junction regions. For noncollinear relative magnetizations, the non-conserved spin currents in a given region gives rise to net torques that can switch directions at particular magnetic configurations or values. The details of the spin current behavior are shown to depend strongly on the degree of magnetic inhomogeneity in the system, including the number of layers and the relative widths of the and layers.

## I Introduction

When a phase difference, , exists between two superconductor () electrodes separated by a non-superconducting material in a Josephson junction, the corresponding charge supercurrent is directly controllable via . Motivated by the interplay between ferromagnetism and superconductivity, researchers are also interested in the dc Josephson effect in superconducting junctions that contain a central ferromagnet () region, which in turn can give rise to an additional spin degree of freedom. More specifically, this kind of Josephson effect provides a venue for the study of spin currents that can be manipulated in cryogenic spintronic systems. eschrig (); hikino (); grein () Besides numerous practical applicationsgolubov () involving these -based Josephson junctions, it is found that novel and interesting phenomena can arise. For example, the realization of a state, ryazanov (); golubov (); ah (); hv4 (); hv5 () where the ground state of the system corresponds to across the junction. Moreover, if the ferromagnet region consists of at least two layers that each have a uniform magnetizations (e.g., a structure), manipulation of the angle between the magnetization vectors can serve to generate long range triplet supercurrents bergeret (); trifunovic (); pajovic (); Shomali (); baker (); hbvprl () in addition to the ordinary singlet ones. Additional control of the magnetic state can also occur from the spatially varying spin current within the layers of the junction, causing mutual torques to act on their respective magnetic moments. Therefore, -based junctions that contains multiple layers, present many opportunities for controlling the charge and spin currents, and their influence on the magnetization in terms of the torque they produce.

Although interest in the study of multilayer structures has recently increased considerably, work on Josephson junctions actually started long ago. The Josephson and critical current oscillations were found to occur as a function of the ferromagnet exchange fieldbuz (); buz2 () and the thickness of the magnetbuz2 (). An essential principle behind many of these important phenomena is the damped oscillatory nature of the singlet Cooper pairs in the ferromagnetic regions, and the associated phase shift in the superconducting order parameter. Due to the intrinsic exchange field in the ferromagnet regions, electrons of Cooper pairs with spin-up (down) decrease (increase) their kinetic energy and the Cooper pairs acquire a nonzero center-of-mass momentum. It further leads to an oscillatory order parameter in the regions demler (). Owing to this oscillatory nature, not only the Josephson critical current but also the superconducting critical temperatures oscillate as a function of exchange field and thickness of magnets. Josephson junctions of the type can also be realized by using this principle by adjusting either the exchange field, the magnet width, or both hv5 (); robinson2 (); linder2 (). The proximity effects between the and regions thus give rise to phenomena buzdin1 (); demler (); hv1 (); klaus () that subsequently play crucial roles in the charge and spin currents that may be manipulated in low temperature nanoscale devices, including nonvolatile memory elements, where the dissipationless nature of the supercurrent flow offers reduced energy loss and Joule heating.

In equilibrium, singlet Cooper pairs carry no net spin, hence any spin current in the system either flows only within the ferromagnets due to their exchange interaction, or it flows by means of induced equal-spin triplet correlations, where the Cooper pairs have a net spin of on the spin quantization axis and they can reside in both and regions. The generation of long-range triplet proximity effects in superconducting heterostructures with magnetic inhomogeneity has been both theoretically predicted and experimentally confirmed: by introducing magnetic inhomogeneity, e.g., inclusion of an additional magnet with misaligned exchange field, the Hamiltonian no longer commutes with the total spin operator and equal-spin triplet correlations can then be induced. bergeret2 (); hvb08 (); hbvprl () Due to the imbalance between majority and minority spins in a ferromagnet, conventional singlet superconducting correlations do not survive for long once inside the magnetic region. However, Cooper pairs with electrons that carry the same spin are not subject to paramagnetic pair-breaking and can in principle propagate for large distances inside the ferromagnet, limited only by coherence breaking processes. bergeret (); bergeret2 (); hvb08 (); hbvprl () Such equal-spin triplet correlations thus play an important role in Josephson junctions with inhomogeneous ferromagnets. bergeret (); alidoust_spn () Indeed, it has been reported experimentally that with the presence of magnetic inhomogeneity, the Josephson critical current decays much more slowly with increases in the F layer thicknesses, robinson (); khaire () as compared to junctions with homogeneous magnetization. One of the simplest ways to introduce magnetic inhomogeneity in a Josephson junction, is through the insertion of bilayer or trilayer of uniformly magnetized ferromagnets. Experimentally, this has the advantages of reproducibility and easy manipulation of the relative exchange field orientation. These structures also provide direct evidence of triplet correlations. ilya (); wvh14 (); volkov (); birge () Recently, long-range coherent transport of triplet pairs was studied in double magnet junctions, richard (); hikino () further demonstrating that it is not always necessary to have a trilayerbuzdin () ferromagnet structure of misaligned ferromagnets to generate equal-spin triplet components to the supercurrent. This finding was in consistent with the long-range phenomena found in a similar structure trifunovic () with asymmetric widths, trifunovic () and orthogonal exchange fields. The behavior of the triplet amplitude is often anticorrelatedwvh14 () to that of the critical temperature and associated singlet correlations, suggesting singlet to triplet conversion.

The oscillatory and long-range pair correlations also lead to new behaviors in the current-phase relations (CPRs) in juncitons pajovic (); golubov (); gold () with a nontrivial magnetic structure. The CPR can in general contain not only the first harmonic but also higher order harmonics, i.e. . The appearance of additional harmonics in the current phase relation has also been discussed in the diffusive and clean regimes for ferromagnetic Josephson junction structures. buzdin1 (); 2nd_hrmnc_3 (); zareyan (); alidoust_spn () In conventional Josephson junctions without any magnetic interaction, the magnitude of is much smaller than that of . However, in Josephson junctions, near a to phase transition, the roles of the first harmonic and the second harmonic can be reversed, and the CPR can be largely dominated by the second harmonic 2nd_hrmnc_3 (). In this regime, both and states can be stable or metastable, and they can coexist. The physical origins of the second and higher order harmonics are believed to lie in the long ranged triplet component of the supercurrent. In this regard, within the vicinity of the - transitions, the triplet correlations can be tuned accordingly. Since the first harmonic is suppressed due to the supercurrent flow reversing direction, the higher order harmonics are revealed at the 0- transition point. The influence of interface scattering on the higher order harmonics were investigated in the quasiclassical clean limit2nd_hrmnc_1 () and experimentally detected. 2nd_hrmnc_2 () The measured supercurrent at the 0- transition point ryazanov () was attributed to the presence of higher harmonics. Subsequent work with ferromagnetic Josephson junctions demonstrated that the higher harmonics can naturally arise when varying the location of domain walls dw (), and also in ballistic double magnetic Josephson junctions, provided that the thicknesses of the magnetic layers are unequal trifunovic (). Recently, evidence of higher harmonics has been observed in Josephson junctions with spin dependent tunneling barriers. pal ()

The interaction of the spin current with the magnetization in layered ferromagnetic junctions with multiple ferromagnets has important consequences for memory technologies. Indeed, storage of information bits depends on the precise relative orientation of the magnetizations in two layers, where nonconserved spin currents reflect the mutual torque acting on the magnetic moments. The corresponding spin transfer torque (STT) can also switch magnetizations when a spin-polarized electrical current flows perpendicular to the layers. Spin-transfer torque is known to occur in a very broad variety of materials, making it an attractive switching mechanism. For equilibrium spin currents, governed by spin-polarized Andreev bound statessauls (), tuning the supercurrent (via ) directly influences the STT when varying the relative in-plane magnetization angle. brou () The direction of supercurrent flow however is not simply related to the direction of the induced torque that tends to align the magnetic moments. The triplet correlations generated in these types of Josephson junctions (with noncollinear relative magnetizations), can also induce spatial variations in the spin currents responsible for the mutual torques acting on the ferromagnets. Shomali ()

When considering superconducting proximity effects, it is important to make sure that the self-consistency condition for the pair potential is fulfilled wvh14 () in order to obtain the correct physical picture. The self-consistency condition is often neglected in the literature liu () mainly because it is difficult to properly implement in the theoretical studies. As we will show in the Sec. II, a source term in the continuity equation for the charge currents usually arises when a non self-consistent superconducting order parameter is used. More importantly, when the self-consistent condition is achieved, the free energies of these proximity-coupled systems are properly minimized. This concept is crucial especially when studying Josephson junctions, since the superconducting proximity effects are the fundamental mechanism behind the nontrivial charge and spin currents that flow within these structures. Furthermore, when the solutions are self-consistent, the charge conservation law is satisfied by properly accounting for the proximity effects and transport properties for the charge current riedel (). When magnetic inhomogeneity is present, the spin density is not conserved and STT can arise and interact with the charge dependent quantities. Although there is no continuity equation for the spin density, since its gradients are nonzero, there is still a corresponding, and fundamental, conservation law that balances current gradients and the STT wvh14 (); linder ().

Above, we have discussed various physical phenomenon associated with the magnetic inhomogeneity such as triplet correlations and the generations of STT. Thus understanding the interplay between spin and charge transport with the long-range proximity effects is an important topic and constitutes the main goal of this work. In this work, we therefore consider nanoscale Josephson junctions, where the region contains multiple layers: spin valves consisting of two metallic ferromagnetic layers separated by a non-magnetic normal metal spacer, and trilayer junctions. In each scenario, a supercurrent is established via a phase difference between the terminals, and different relative magnetization orientations are considered. The paper is organized as follows: We describe our method in Sec II, and derive conservation laws for both spin and charge currents. Our method is based on solving the microscopic Bogoliubov-de Gennes (BdG) equations self-consistently. All important physical quantities, e.g., the magnetization, can be extracted or constructed from the self-consistent solutions. In Sec III, we present a detailed study of the transport properties. We conclude with our main findings in Sec IV.

## Ii Methods

The general method that we use in this paper is that of numerical diagonalization of the self-consistent Bogoliubov-de Gennes (BdG) equations. Since many aspects of this method have been extensively discussed in previous workhv1 (); dw (); wvh12 (); wvh14 () we will only include here a brief review of these points, as needed to make this paper understandable. We will discuss in more detail additional aspects needed for the transport calculations described in this work.

The derivation of the BdG equations for general magnetization configurations begins with the effective BCS Hamiltonian, ,

 H=∫d3r{ψ†(r)[He− h⋅σ]ψ(r)+Δ(r)ψ†↑(r)ψ†↓(r) +Δ∗(r)ψ↓(r)ψ↑(r)}, (1)

where are the usual fermionic operators, , and denote the set of Pauli matrices. We describe the magnetism of the layers by effective Stoner exchange fields which in our case have components in all directions (see Fig. 1). The spin independent scattering potential is denoted by , and is the pair potential.

To diagonalize the effective Hamiltonian, the field operators and are expandedbdg () by means of a Bogoliubov transformation:

 ψ↑(r) =∑n(un↑(r)γn−v∗n↑(r)γ†n), (2a) ψ↓(r) =∑n(un↓(r)γn+v∗n↓(r)γ†n), (2b)

where and are the quasiparticle and quasihole amplitudes, which are chosen so that the Hamiltonian is diagonalized in terms of the fermionic operators. Therefore, and . Also, the thermal expectation values involving and are given by the usual Fermi functions song () The anticommutation relations for and yield,

 [ψ↑(r),H] =(He−hz)ψ↑(r)−[hx−ihy]ψ↓(r)+Δ(r)ψ†↓(r), (3a) =(He+hz)ψ↓(r)−[hx+ihy]ψ↑(r)−Δ(r)ψ†↑(r). (3b)

It is convenient at this point to simplify to the quasi one-dimensional geometry (Fig. 1) of our problem. Then, by inserting (2) into (3) and using the commutation relations, we obtain the general spin-dependent BdG equations for this geometry,

 ⎛⎜ ⎜ ⎜ ⎜⎝H0−hz−hx+ihy0Δ−hx−ihyH0+hzΔ00Δ∗−(H0−hz)−hx−ihyΔ∗0−hx+ihy−(H0+hz)⎞⎟ ⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜ ⎜⎝un↑(x)un↓(x)vn↑(x)vn↓(x)⎞⎟ ⎟ ⎟ ⎟⎠ =ϵn⎛⎜ ⎜ ⎜ ⎜⎝un↑(x)un↓(x)vn↑(x)vn↓(x)⎞⎟ ⎟ ⎟ ⎟⎠, (4)

where are the quasiparticle energies, is normal to the layers, which lie in the plane (see Fig. 1), and is the pair potential, to be found self consistently as explained below. Here the single particle Hamiltonian is written,

 H0=12m(−∂2∂x2+k2y+k2z)−EF+U(x). (5)

The components of the exchange field in each of the layers take the form:

 hi=h(cosθi,sinθisinϕi,sinθicosϕi), (6)

where denotes one of the magnetic layers. We will assume that the magnitude of the exchange field is the same in all magnetic layers, and that it vanishes elsewhere. The angles and will in general be taken to vary from layer to layer.

One obtains in the usual waybdg () the self-consistency condition for the pair potential, using . Here is the superconducting coupling constant, which vanishes in the non- layers. In our geometry we find, after using the Bogoliubov transformation and making use of the appropriate averages such as , the pair potential can be expressed in terms of the quasiparticle amplitudes as an appropriate sum over states:

 Δ(x)=g2∑n′[un↑(x)v∗n↓(x)+un↓(x)v∗n↑(x)]tanh(ϵn/2T), (7)

where the prime on the sum indicates that only those states that have energies within a “Debye energy”, , are included.

The problem is then solved iteratively: the potential is initially taken to be in the first layer and in the second layer, where is the initial guess of the magnitude for the pair amplitudes. The Hamiltonian is then numerically diagonalized and the new pair potential is found via Eq. (7). Iteration is continued until convergence. The details of the procedure, including the way to ensure, in the Josephson calculations, that the phase difference between the right and left ends of the sample remains , is explained in Appendix A.

Equilibrium and transport properties in the structures considered are strongly influenced by the existence of “odd” triplet pairs. The existence of such pairs is allowed by conservation laws since, unless all of the layers have magnetizations along the same direction, the total spin of the Cooper pairs is not a conserved quantity. Because of the Pauli principle, these -wave triplet pairs must have wavefunctions odd in frequencybrz () or equivalentlyhbvprl (); hvb08 () in time. Within the BdG framework, the time formulation is much more convenient. Accordingly, we will describe the triplet pair correlations via the following amplitude functions, in terms of the field operators:

 f0(r,t)= 12[⟨ψ↑(r,t)ψ↓(r,0)⟩+⟨ψ↓(r,t)ψ↑(r,0)⟩], (8a) f1(r,t)= 12[⟨ψ↑(r,t)ψ↑(r,0)⟩−⟨ψ↓(r,t)ψ↓(r,0)⟩]. (8b)

Taking the quantization axis along the direction, the triplet amplitudes, and , can be rewrittenhbvprl (); hvb08 () in terms of the quasiparticle amplitudes:

 f0 =1/2∑n(g↑↓n−g↓↑n)ζn(t), (9) f1 =1/2∑n(g↑↑n+g↓↓n)ζn(t), (10)

where , and we define . It is sometimes necessary to evaluate the triplet amplitudes along a different spin axis. For example, one may wish to use the direction of the local magnetization (defined below) as the axis of quantization. To do so one rotates the quantization axis so that it is aligned with the local magnetization direction using the spin rotation matrices discussed in Appendix B).

We will consider here spin currents, as well as charge currents. In our structures spin transport is influenced by the leakage of magnetism out of the layers and into the superconductors. This can be characterized by the local magnetization ,

 m(r)=−μB⟨η(r)⟩, (11)

where is the spin density operator,

 η(r)=ψ†(r)\boldmathσψ(r), (12)

and the Bohr magneton. For our quasi-1D geometry, we can rewrite the components of in terms of the quasiparticle amplitudes:

 mx(x) =−μB∑n{[u∗n↑(x)un↓(x)+u∗n↓(x)un↑(x)]fn −[vn↑(x)v∗n↓(x)+vn↓(x)v∗n↑(x)](1−fn)}. (13) my(x) =−iμB∑n{[un↑(x)u∗n↓(x)−un↓(x)u∗n↑(x)]fn +[vn↑(x)v∗n↓(x)−vn↓(x)v∗n↑(x)](1−fn)}. (14) mz(x) =−μB∑n{[|un↑(x)|2−|un↓(x)|2]fn +[|vn↑(x)|2−|vn↓(x)|2](1−fn)}. (15)

We now turn to the appropriate expressions for the currents. As stressed in the Sec. I, one needs to carefully establish proper conservation laws when discussing the transport properties of the system kadanoff (). We first discuss the charge supercurrent. In our geometry the charge current has only one component, which depends on the coordinate. In the absence of an external magnetic field, the total charge current, , is found from the standard quantum mechanical expression, . This leads to the result:

 Jxσ(x)=e2m⟨−iψ†σ∂∂xψσ+i(∂∂xψ†σ)ψσ⟩. (16)

This expression for the current ensures, together with the self consistency condition, that charge conservation is satisfied, that is, in the steady state. wvh14 (); dw (); sols () It is of course convenient numerically to rewrite the expression for the current in terms of the calculated quasiparticle amplitudes and energies. After inserting the Bogoliubov transformations in Eq. (2), we can write the total charge current, as given by Eq. (16) summed over spins as:

 Jx(x) =2em∑n,σIm[unσ∂u∗nσ∂xfn+vnσ∂v∗nσ∂x(1−fn)]. (17)

One can verify once again the conservation law by taking the divergence of the current in Eq. (17) and using the BdG equations (4), to find: dw (); wvh14 ()

 ∂Jx(x)∂x=2eIm{Δ(x)∑n[u∗n↑vn↓+u∗n↓vn↑]tanh(ϵn2T)}. (18)

When the self-consistency condition is satisfied, the right hand side vanishes, and charge is properly conserved. If the self-consistency condition is not strictly satisfied, the terms on the right side act effectively as sources or sinks of current.sols (); wvh14 (); dw () We will consider large contacts with the amplitude and phase of the order parameter determined self-consistently except near the sample edges (see Appendix A) where sources and sinks of charge exist via the implicit external electrodes. This gives the necessary charge conservation condition in the region of interest. We emphasize here that with self-consistent solutions, we are able to correctly determine the effect of triplet correlations on both the charge and spin transport.

The extension of the above considerations to spin transport is relatively straightforward.brou (); wvh14 () As in the case of the charge density, the Heisenberg picture is utilized to determine the time evolution of the spin density, ,

 ∂∂t⟨η(r,t)⟩=i⟨[H,η(r,t)]⟩, (19)

where is given in Eq. (12). The associated continuity equation now reads,

 ∂∂t⟨η(r,t)⟩+∂S∂x =τ+JS, (20)

where is the spin current which in our geometry is a vector (in general it is a tensor). The spin-transfer torque, , is given by:

 τ=−i⟨ψ†(r)[h⋅σ,σ]ψ(r)⟩=2⟨ψ†(r)[σ×h]ψ(r)⟩. (21)

The term has components,

 JSx =2Im{Δ⟨ψ†↓(r)ψ†↓(r)−ψ†↑(r)ψ†↑(r)⟩}, (22) JSy =2Re{Δ⟨ψ†↓(r)ψ†↓(r)+ψ†↑(r)ψ†↑(r)⟩}. (23)

For the -wave superconductors considered in this paper, we have by virtue of the Pauli principle, since only equal time correlations are involved. Thus in the absence spin transfer torque, we have . However, in generalwvh14 () spin transfer torque is present and, in the steady state, the derivatives of with respect to do not vanish.

The expression for the spin-current, , is found from taking the commutator in Eq. (19) and using Eq. (II):

 S=−i2m⟨ψ†(r)σ∂ψ(r)∂x−∂ψ†(r)∂xσψ(r)⟩, (24)

where we recall that the vector represent spin current flowing along the direction for our quasi-one-dimensional systems. We can now expand each spin component of the spin current in terms of the quasiparticle amplitudes to obtain:

 Sx=−i2m∑n{fn[u∗n↑∂un↓∂x+u∗n↓∂un↑∂x−un↓∂u∗n↑∂x−un↑∂u∗n↓∂x] −(1−fn)[vn↑∂v∗n↓∂x+vn↓∂v∗n↑∂x−v∗n↑∂vn↓∂x−v∗n↓∂vn↑∂x]}, (25) Sy=−12m∑n{fn[u∗n↑∂un↓∂x−u∗n↓∂un↑∂x−un↓∂u∗n↑∂x+un↑∂u∗n↓∂x] −(1−fn)[vn↑∂v∗n↓∂x−vn↓∂v∗n↑∂x+v∗n↑∂vn↓∂x−v∗n↓∂vn↑∂x]}, (26) Sz=−i2m∑n{fn[u∗n↑∂un↑∂x−un↑∂u∗n↑∂x−u∗n↓∂un↓∂x+un↓∂u∗n↓∂x] −(1−fn)[−vn↑∂v∗n↑∂x+v∗n↑∂vn↑∂x+vn↓∂v∗n↓∂x−v∗n↓∂vn↓∂x]}. (27)

In the case of layers with uniform magnetization, there is no net spin current. The introduction of an inhomogeneous magnetization texture however results in a net spin current imbalance that is finitewvh14 () even in the absence of a Josephson current. This will be discussed in greater detail below.

To compute the spin transfer torque, it is useful to express it in terms of the quasiparticle amplitudes. A convenient approach involves directly taking the expectation values of Eq. (21):

 τ=2⟨ψ†(r)σψ(r)⟩×h=−2μBm×h, (28)

where we have used Eq. (11). The magnetization components are given in Eqs. (13)-(15). Since the exchange field is prescribed, it is the self consistently calculated magnetization that determines the torque acting on the ferromagnet layers. Equivalently, one can use the continuity equation in the steady state to determine the torque transfer by evaluating the derivative of the spin current as a function of position:

 τi=∂Si∂x. (29)

It is however safer to evaluate both sides of Eq. (29) independently and use this equation as a consistency check. We have performed extensive numerical checks using this procedure. In most of the results presented we have calculated the torques using Eq. (28), thus avoiding the numerical derivatives that arise when using the right side of Eq. (29). Additional physical insight can be gained by integrating Eq. (29) over a particular region, e.g., :

 Sx(b)−Sx(a)=∫F1dxτx=τx,tot, (30)

which means that the change in spin current through (from to ) is equivalent to the net torque acting within those boundaries.

## Iii Results

The results of our systematic investigations are presented below in terms of convenient dimensionless quantities. Our choices are as follows: all length scales, including the position , and widths () are normalized by the Fermi wavevector, . For the superconducting correlation length we choose the value , and the computational region occupied by the electrodes corresponds to a width of (see Appendix B for numerical details). All temperatures are measured in units of , the transition temperature of bulk material, and we consider the low temperature regime, . Energy scales are normalized by the Fermi energy, , including the Stoner field interaction and the energy cutoff, , used in the self-consistency condition, Eq. (7). The latter is set at : results are independent of this cutoff choice. As mentioned above, the strength of the magnetic exchange fields, is taken to be the same in both magnets: we set its dimensionless value to a representative . We vary the orientation angles of the magnetic exchange field in each of the regions, depending on the quantity being studied. The magnetization is normalized by , where is the electron density, . The normalization for the torque follows from the normalizations for and and Eq. (28). When presenting results for the currents, we normalize the charge current densities by , where , and is the Fermi velocity. All three components of the spin current are normalized similarly, by the quantity , where involves the normalization of and a factor of . The interface scattering is represented by delta functions of strength at all the interfaces. The corresponding dimensionless parameter is . The self-consistency of the pair potential that characterizes an accurate representation of the Cooper pair correlations throughout the system, is associated with the proximity effects and depends to varying degrees on the parameters outlined above. In some cases the dependence is rather obvious: for instance, large results in weaker proximity effects. In other cases it is more intricate and will be analyzed more carefully.

### iii.1 Current-Phase relations

We begin by showing our results for the self-consistent current phase relations in a simple double layer ferromagnet Josephson junction. At the interfaces between the and regions, quasiparticles undergo Andreev and conventional reflections.radovic2 (); radovic1 (); beenaker2 (); beenaker1 () The superposition of these waves in the regions results in subgap bound states that contribute, together with the continuum states, to the total current flow. In Fig. 2, we show the supercurrent as a function of phase difference (current phase relation CPR) for two ferromagnets of unequal width: and . This asymmetric choice of widths helps ensuretrifunovic (); ilya () that equal-spin triplet correlations are generated in the system. The angular parameters in this figure are fixed at , and , corresponding to in-plane magnetization orientations. The layer has its magnetization aligned in the direction (). The first two panels, (a) and (b), display three different relative in-plane magnetization configurations in the layers: parallel (P) (), antiparallel (AP) (), and normal (N)convention () (). Two different strengths of the interface scattering parameter are considered. In (a) there is no interface scattering (), while in (b), a rather high rate of scattering is present, with . The CPR for the collinear configurations (P or AP) possesses the conventional periodicity, and the supercurrent flows oppositely for the two alignments. When the relative magnetizations are orthogonal to each another, the CPR becomes periodic as revealed by the sawtooth-like pattern in (a) or the more sinusoidal behavior in (b), both of which change sign at . This is a consequence of the emergence of equal-spin triplet correlations2nd_hrmnc_1 (); 2nd_hrmnc_2 (); 2nd_hrmnc_3 (); alidoust_spn (); buzdin1 (); zareyan () that are absent when the exchange fields in the ferromagnets point along the same direction. When strong interface scattering is present, Fig. 2(b) shows that the periodic CPR (N case) is substantially diminished, relative to the P or AP cases. This is because the proximity effect is weakened, with a resulting reduction of the associated equal-spin triplet correlations. To further examine the effects that interface scattering has on this -periodic supercurrent, we consider in Fig. 2(c), the same junction with varying degrees of scattering strengths and with the relative in-plane magnetizations fixed and orthogonal to each another (). Increasing clearly leads to a crossover in the CPR from a sawtooth to sinusoidal form and to a marked reduction of the supercurrent flow. As this occurs, the phase difference yielding the critical current density also declines.

We next consider, in Fig. 3, the effect on the CPR when changing and in a structure. Panels (a)-(e) label the various widths considered. We keep the total junction length in which supercurrent flows constant, i.e., , is fixed to a representative value of . In these panels, we consider the same three orthogonal magnetization orientations as in Fig. 2. One sees that, in the N configuration, the structure with the greatest geometric asymmetry, as given by the ratio , tends to have a more pronounced superharmonic CPR relative to the P and AP collinear configurations. Upon increasing the width , the jagged sawtooth peaks at become smoothed and the -periodic CPR is reduced substantially. Eventually in panel (e), where both widths are equal, the current no longer undergoes a sign change and the additional harmonics in the CPR that previously reflected the existence of equal-spin triplet correlations are now drastically modified. Remarkably, for this situation, over most of the range, the linear variations of the currents in the orthogonal and antiparallel configurations tend to overlap. Moreover, it is also evident that the currents in the N and AP configurations that flow opposite to the currents in the magnetically uniform P case, are substantially greater. The observed reversal of current direction for a given magnetic configuration when changing the widths is a direct consequence of the damped oscillatory behavior of the singlet and triplet correlations, as will be discussed below (see also Fig. 8). By comparing Fig. 3(a) with Fig. 2(a) and noting that the only difference between them is the thickness of the second magnet, one finds that the magnitude of the supercurrent for the N case with thicker does not drop significantly, as it does in the P and AP configurations. Such behavior is a signature of the long-range nature of the equal-spin triplet correlations. To show further the effects of increased width on the supercurrent, we consider in Fig. 3(f) the symmetric geometry configuration for several equal layer widths . As shown in the legend, a broad range of are considered. It is clear from the figure that the equal-spin triplet correlations are strongly impacted, and only the -periodic supercurrent arises. By increasing the ferromagnet widths, the rate at which the current changes with tends to decline for wider junctions, and the “critical” phase difference where the current is suddenly reduced becomes smaller. In this panel the ferromagnetic region is no longer constrained to have the same total width, hence increasing reduces the overall current flow.

After these examples of two-magnet Josephson junctions, we now consider a trilayer junction, where a nonmagnetic normal metal “spacer” separates the two ferromagnets. Such spacers are often needed experimentally when it is wishedilya () to rotate the magnetization in one magnet only. To focus on the case where the CPRs have important additional harmonics, we keep the in-plane mutual magnetizations orthogonal (with , and ). In Figs. 4(a) and (b), the ferromagnet widths are set at and , while the width of the normal metal spacer, , varies from zero to 500 (corresponding to ) as indicated in the legend. The interface scattering parameter has the value for the curves shown in (a), whereas in (b) we have . In either case, the effect of increasing the spacer width is to cause an overall reduction in the supercurrent flow. For large interface scattering, the CPR becomes less sensitive to , as seen by comparing (a) and (b). A wider junction with is shown in (c), with ferromagnet widths and . The details of the usual sawtooth CPR reveals that in these cases, the supercurrent flow can be quite sensitive to the spacer width, as it abruptly reverses direction at considerably different phase differences, for incremental changes in . In this panel, results for an additional small value of spacer thickness, are also shown. Since the -periodic CPR is closely related to the generation of the equal-spin triplet correlation, one can infer from Fig. 4 that the introduction of an additional non-magnet metallic layer can quantitatively changes the transport properties.

To examine the effects of increased magnetic inhomogeneity, we show in Fig. 5 results for a pentalayer junction. This structure is complementary to that studied in Fig. 3, the main difference being an additional ferromagnet layer with an out-of-plane magnetic exchange field oriented in the direction (corresponding to in Fig. 1). The relative magnetic orientations are labeled in the legend by the directions of the axes along with the magnetizations are aligned in each layer. For example, denotes a sample in which and are normal to each other (the configuration labeled N in the previous figures) with an additional out-of plane magnetization in . The width of the layer is identical to that of the layer: in (a), and in (b). The thicker ferromagnet has width . The -periodic CPR that arises in double magnet junctions with orthogonal in-plane magnetizations remains relatively unchanged by the addition of the additional out-of-plane intermediate ferromagnet. However, it was shown in Figs. 3(a)-(e) that the collinear P or AP magnetic states with still periodicities behaved in an approximately piecewise linear fashion and that the current maintained its direction when varying . Now, on the other hand, the insertion of an additional layer between the collinear ferromagnets, makes it possible for equal-spin triplet correlations to be generated and the supercurrent becomes drastically modified. For phase differences in the vicinity of or , Figs. 5(a) and (b) show that the current is approximately linear in the phase difference, but there is a broad intermediate range of , where the supercurrent flow is relatively uniform. Moreover, for either the P or AP configuration, varying the phase can result in the Josephson current switching direction. These trends are the same for each of the cases presented in (a) or (b).

### iii.2 Magnetic orientation and CPR

Having discussed the current phase relations for a few different ferromagnetic Josephson junction configurations, we now will study in more detail the effect of varying magnetization orientations on the CPR. We will set the macroscopic phase difference to a prescribed value and study the supercurrent response within the junctions for a range of relative magnetization orientations.

Beginning with a basic Josephson junction whose phase difference is set at , we consider in Fig. 6(a), the normalized supercurrent density as a function of the in-plane magnetization angle (), where we introduce , with being the angle shown in Fig. 1. In terms of , the P and AP states then correspond to , and respectively (since the magnetization in is fixed along ). Four interface scattering strengths are considered in panel (a), as indicated in its legend. In all cases, by tuning the relative alignment angles, supercurrent switching occurs when the mutual magnetization orientations are approximately orthogonal. As expected, the supercurrent flow is greatest for transparent interfaces (), and decreases with increasing , as the sensitivity to becomes weaker. It is notable that the maximum current flow occurs at different values, depending on the interface scattering strength. In Fig. 6(b), an additional layer, of variable strength as indicated, is inserted between the two ferromagnets and the interface scattering parameter is set to . The solid curve corresponds to the curve in panel (a). The figure shows that increasing the layer thickness tends to generally dampen the current through the junction. The current flow peaks when the two ferromagnets are aligned in the P state and it vanishes altogether near the N configuration, before reversing direction as the relative magnetizations approach the AP state. In Fig. 6(c) we investigate the supercurrent flow in a wide junction, with , and three different widths (see legend). The macroscopic phase difference is set at . It is seen that for this geometry, the current has some contrasting features compared to the previous cases involving thinner layers. For instance, when tuning , now the weakest current flow occurs when in the P state, and the maximum occurs in the orthogonal configuration. For all width considered, the current undergoes rapid changes in the vicinity of the middle value between N and P or N and AP configurations, with the case changing the most, and then remains nearly constant at angles near the P and AP configurations. Since the prescribed phase difference is near where the maximum currents occur in -periodic Josephson junctions, it is evident that the results in Fig. 6(c) are a direct consequence of the induced equal-spin triplet correlations. We will see below in Sec. III.3 that this correlates with the triplet generation behavior.

We next move to the study of the magnetization orientation role in the supercurrent for more complicated junctions. As an example, we consider a scenario where the and layers have magnetizations pinned along the and directions respectively, that is, a relative N configuration, but the magnetization in the central layer rotates on the plane () from along the -axis to the -axis. In Fig. 7(a), we show the normalized current density as a function of for few different layer widths, . We set , and consider values of that lead to both symmetric and asymmetric structures. In each case, , and interface scattering is present with , except for the widest junction with where the interfaces are transparent (). In the case , where all thicknesses are identical, the current flow is approximately antisymmetric around . The current flow is suppressed for the asymmetric situation by decoherence effects arising from the larger width. While for the simpler junctions in Fig. 6 was periodic in , variations in for trilayer ferromagnetic junctions are in general periodic, as seen in the figure. The asymmetric case also demonstrates a more intricate structure as the magnetization angle is swept. For the narrower junction with , the orientations in which the current switches direction are near and . When , the current never changes its direction. Although the rotating out-of-plane exchange field of does affect the strength of current with periodicity. For the highly asymmetric case, the current again maintains its direction of flow over the full angular range of . Also, the current is approximately constant except for orientations when the exchange field in points along : or . For orientations near , we find that the current is strongly suppressed. In Fig. 7(b), the in-plane magnetization in the first ferromagnetic layer is now varied, while the other two are kept fixed: for , it is along , and for , it is along , an N configuration. Two types of structures are considered, with the total width of the three ferromagnetic regions being constant. In the first case, , and , while in the second one all three layers have width, . Both cases exhibit similar behavior as a function of . The current vanishes when is antiparallel to and is highest when the magnetization lies nearly in between the those of and . Thus, the charge supercurrent which flows oppositely in the two structures, can be effectively switched on or off by manipulating the in-plane magnetization angle of the first ferromagnet.

It is known that in heterostructures, including bilayersklaus () and Josephson junctions, variations in the magnetic exchange field and ferromagnet thickness induce damped oscillations in the spatial behavior of the Cooper pair amplitudes, resulting in modulation of physical quantities as a function of either or . The damped oscillations in the clean limit have a wavelength that goes as the inverse of the exchange field. To investigate this phenomenon in an junction, we show in Fig. 8, the supercurrent that flows through the junction as a function of the width, . For the adjacent layer, the width is , and with interface scattering strength . Two exchange fields are considered: and . The period of oscillations in each case is seen to be approximately , where , and , respectively. The current for both cases is maximal when the two regions are equal, and then slowly dampens out with increasing . This decay length is inversely proportional to the magnitude of the exchange field. Moreover, the charge current in each case periodically changes sign. Since the oscillations in the current as a function of increase with larger exchange fields, the current direction can be very sensitive to fabrication tolerances for strong magnets.

### iii.3 Induced triplet pairing

We now discuss the induced triplet pairing correlations in ferromagnetic Josephson junctions. The presence of multiple misaligned ferromagnets yields both the (Eq. (9)) and the (Eq. (10)) triplet pair amplitudes as permitted by conservation laws and the Pauli principle. To gain an overall view of the opposite-spin triplet amplitudes, , and the equal-spin amplitudes, , we illustrate in Fig. 9 the spatial behavior of these correlations in the and regions of an junction. We focus on the real parts of these complex quantities, keeping in mind that the imaginary components obey similar trends. The geometrical parameters in this figure are , , and (top panels) or (bottom panels). Thus, in all panels the region is occupied by while occupies the region in the top panels and in the bottom panels (c) and (d), where the vertical dotted line denotes the spacer boundary. Hence, different horizontal scales are used in each case. The scattering parameter is set to , and each curve corresponds to a different phase difference, as indicated by the legend. The exchange fields are in-plane and normal to each other. Within the region, panels (a) and (c) reveal that the magnitude of the pair correlations are approximately of the same magnitude, decreasing in the vicinity of the interface located at . The system with the wider normal metal layer, is slightly more sensitive to phase variations. Within the ferromagnet , the same panels (a) (c) show the oscillatory nature of , which behaves similarly to the singlet pair amplitude, the periodicity arising from the difference in spin-up and spin-down wavevectors. For the chosen exchange field, the oscillations are limited in due to the confined width. Turning now to the equal-spin triplet correlations , panels (b) and (d) display behavior which contrasts with the results. In particular, within the narrow region, the triplets are negligibly small, and the correlations clearly dominate. In the region, the correlations nearly vanish, while the equal-spin triplets peak near the interface (at ), before dropping within the normal metal (see Fig. 9(d)). Finally, within , the triplets assume a slow, long-range variation compared to the damped oscillatory behavior of the curves.

In addition to investigating the spatial behavior of the triplet amplitudes, it is instructive to also examine the spatially averaged triplet and singlet correlations as functions of the relevant system parameters. For example, by tuning the relative magnetization angle, (varying at fixed ), important overall features can be revealed. The top panels (a) (b) (c) of Fig. 10 show the -dependence of the magnitudes of the triplet and singlet amplitudes averaged over the region for the structure studied initially in Fig. 6(b). Four representative layer widths are considered as indicated in the legend. The proximity effects and hence coupling of the two ferromagnets diminish with increasing , and therefore the pair correlations become less sensitive to variations in , as observed for the largest case. Other than the diminished magnitudes, the overall trends and behavior however do not depend strongly on the presence of the normal metal spacer. This may be important in experiment design, where spacers are often needed. The opposite-spin triplet correlations, and the singlet pair amplitude () in (a) and (c) behave in rather similar ways, but is more symmetric about the orthogonal direction . When the relative magnetization orientation varies in inhomogeneous systems, the process of singlet-triplet conversion plays a role in the transport and thermodynamic properties of such systems. It is apparent from panels (a) and (b) that the orientation that leads to a minimum in and , corresponds to that where is largest. These occurrences arise when the exchange fields in and are nearly orthogonal. This angle for slightly shifts with variations in .

In the bottom panels (d) and (e) of Fig. 10 we display the spatial dependence of the real components of the triplet correlations throughout each of the three junction regions discussed in the top panels. Results for five different relative magnetic orientations are presented, (see legend) including the (), AP (), and N () configurations. Considering first Fig. 10(d) in the region, we see behavior similar to that found for the wider junction case in Fig. 9, including a relatively weak dependence on the orientation angle (as opposed to the phase difference ). The central region is most affected by variations in : the real part of changes sign when is swept from the P to AP state, and nearly vanishes altogether at . In , a series of oscillations emerge with a periodicity similar to that found in Fig. 9, since an exchange field of the same strength was used. In Fig. 10(e), the equal-spin amplitudes exhibit considerably different behavior. First, only three of the considered yield nonzero results, since the P and AP configurations cannot generate equal-spin triplet correlations. Interestingly, within the region does not exhibit a slow decay, but rather oscillates with a period that is much shorter than the opposite spin pairs governed by the difference in spin-up and spin-down wave vectors. We find that within the normal metal layer, there is nearly a complete absence of equal-spin correlations, this is accompanied by the appearance of opposite-spin correlations (see panel (d)). The amplitudes are largest in the layer, for the relative orientation of , in agreement with the averaged results in Fig. 10(b). For the relative magnetization angles of , and , the ’ amplitudes are identical due to the symmetry about . One can now correlate the features Fig. 10 and 11 to Fig. 4. From Fig. 4, we learn that the N magnetic configuration often leads to the appearance of -Josephson junctions. Here we are able to give concrete proof that in the N cases the equal-spin triplet correlations are maximized and are insensitive to . Therefore, the CPRs for different magnetic configurations are essentially characterized by their detailed singlet/triplet nature.

We next examine, in Fig. 11, the behavior of the averaged triplet and singlet amplitudes, as the magnetic orientation angles, (top panels), and (bottom panels) are changed, in more complicated Josephson junctions. This study is therefore complementary to the results shown in Fig. 7 involving the charge supercurrent. The geometric parameters are