Colloquium: Nonequilibrium effects in superconductors with a spin-splitting field

We review the recent progress in understanding the properties of spin-split superconductors under non-equilibrium conditions. Recent experiments and theories demonstrate a rich variety of transport phenomena occurring in devices based on such materials that suggest direct applications in thermoelectricity, low-dissipative spintronics, radiation detection and sensing. We discuss different experimental situations and present a theoretical framework based on quantum kinetic equations. Within this framework we provide an accurate description of the non-equilibrium distribution of charge, spin and energy, which are the relevant non-equilibrium modes, in different hybrid structures. We also review experiments on spin-split superconductors and show how transport measurements reveal the properties of the non-equilibrium modes and their mutual coupling. We discuss in detail spin injection and diffusion and very large thermoelectric effects in spin-split superconductors.

We review the recent progress in understanding the transport properties of spin-split superconductors under non-equilibrium conditions. Recent experiments and theories demonstrate a rich variety of transport phenomena occurring in devices based on such materials that suggest direct applications in thermoelectricity, low-dissipative spintronics, radiation detection and sensing. We discuss different experimental situations and present a complete theoretical framework based on the quantum kinetic equations. Within this framework we provide an accurate description of the non-equilibrium distribution of charge, spin and energy, which are the relevant non-equilibrium modes, in different hybrid structures. We also review experiments on spin-split superconductors and show how transport measurements reveal the properties of the non-equilibrium modes and their mutual coupling. We discuss in detail phenomena such as spin injection and diffusion, relaxation of the different electronic degrees of freedom, very large thermoelectric effects and the non-equilibrium dynamics of spin-split superconductors. Ferromagnetism and spin-singlet superconductivity are antagonist orders and hardly coexist in bulk materials. However, hybrid nanostructures allow for the possibility of combining the two phenomena via mutual proximity effects. The combination leads to the presence of novel features not present in either system alone. Such features are discussed in this review. We can make a distinction among those characteristics affecting the spectral properties of the materials, showing up when the probed systems are in equilibrium, and those related to nonequilibrium phenomena.

CONTENTS
In the first group of features we highlight the spatial oscillation of the superconducting pair correlations in a superconductor-ferromagnet (S/F) structure that leads to the 0-π transition in S/F/S Josephson junctions, and the appearance of superconducting triplet correlations that may allow for a coexistence of ferromagnetic and superconducting ordering in S/F structures. The equilibrium properties of hybrid superconductor-ferromagnet structures have been studied in great detail and are reviewed for example by Buzdin (2005) and Bergeret et al. (2005). On the superconducting side, the magnetic proximity effect typically resides within a coherence length ξ s from a contact to the ferromagnet, i.e., in typical structures at low temperatures from tens up to a couple of hundred nanometers. An example of the magnetic proximity effect is the spin-splitting field created in a superconducting film in contact with a ferromagnetic insulator. Such a spin-splitting field is of central importance in this article.
A spin splitting can be also created by applying an inplane field to a thin superconducting film. In this case the orbital (diamagnetic) effect is weak and the field induces a paramagnetic response. This was first studied by Chandrasekhar (1962) and Clogston (1962), who found the critical paramagnetic field needed to suppress superconductivity (see Sec. II.F.2). Later on, this problem was also addressed by Fulde and Ferrell (1964) and Larkin and Ovchinnikov (1965). They predicted that around the region of the critical fields a spontaneous transition to an inhomogeneous superconducting state may take place. This FFLO state is characterized by an oscillatory superconducting order parameter. Much later, an analogue of this effect was measured in the oscillations of the critical current through a magnetic weak link (Kontos et al., 2001;Ryazanov et al., 2001).
Another consequences of the spin-splitting field is the spin splitting of the superconducting density of states (DOS), which was first explored experimentally by Meservey et al. (1975Meservey et al. ( , 1970. In those experiments the field was applied in the plane of a superconducting film of thickness much smaller than the London penetration length, and also smaller than the superconducting coherence length. In this case the magnetic field penetrates the film almost uniformly and the screening, Meissner, currents are negligibly small, leading to the dominance of the paramagnetic spin effects. The spin-split DOS was utilized to determine the spin polarization of ferromagnets attached to such superconductors (Meservey et al., 1980;Meservey and Tedrow, 1994;Paraskevopoulos et al., 1977;Meservey, 1971, 1973). These experiments lead to more refined studies of the superconducting state, including Fermi-liquid effects, and spin-flip and spin-orbit scattering (Alexander et al., 1985). As anticipated above, a similar effect can also arise in thin superconducting films from a magnetic proximity effect to a nearby magnetic material (Tedrow et al., 1986). In such a case spin splitting of the density of states can be observed by applying very small magnetic fields or even at zero field. These situations are discussed in Sec. II.
Recently, the combination of spin-splitting fields with strong spin-orbit interaction in superconducting nanowires has raised considerable interest as a platform for realizing topological phases and Majorana fermions, with possible applications in topological quantum computing (Aasen et al., 2016). These effects are outside the context of this review, but also there inducing large Zeeman splitting via the magnetic proximity effect can be a great advantage.
In the present review we focus on the second group of effects mentioned in the first paragraph, i.e., those related to nonequilibrium properties of superconductors in the presence of a spin-splitting (exchange or Zeeman) field h. These nonequilibrium effects can survive to much higher distances than ξ s , as their decay scales are determined via the various inelastic and spin-flip scattering lengths. Moreover, they can be studied at a weak tunneling contact to ferromagnets, making the analysis in some cases more straightforward than in proximity experiments. Nonequilibrium properties are related to the deviation of the electron distribution function from its equilibrium form, which leads to a nonequilibrium distribution (imbalance) of charge, energy or spin degrees of freedom. We refer to these different types of deviations from equilibrium as nonequilibrium modes. 1 We explore the coupling between these modes in superconductors with a spin-splitting field, and discuss unusually strong thermoelectric response and long-range spin signals.
The above mentioned ability to characterize the spin polarized Fermi surface of metallic magnets with the help of spin-split superconductors has a direct connection with 4 FIG. 1 Schematic pictures of various nonequilibrium phenomena occurring at normal metal/superconductor (NS) and ferromagnetic/superconductor (FS) interfaces discussed in this review. (a) Spin-resolved tunneling from a ferromagnetic metal to a spin-split superconductor that leads to the spin valve effect. (b) Creation of spin accumulation in the voltage biased FS junction. (c,d) Schematic picture of thermally excited currents in NS and FS junctions with a spin-split superconductor. (c) Spin Seebeck effect in NS junction: A pure spin current is generated by the temperature bias between a superconductor at temperature TS and a normal metal at temperature TN > TS. (d) Thermoelectric effect in a FIS junction: Here the spin current is partially converted to the charge current due to the spin-dependent density of states in the ferromagnet. spintronics, and in particular with the search for spin valves with larger efficiencies than in the structures exhibiting giant magnetoresistance (Baibich et al., 1988;Binasch et al., 1989). Indeed, a superconductor with a spin-splitting field has also an intrinsic energy dependent spin polarization around the Fermi level. This allows for studying different spintronic effects in a setting of a controllable non-linearity arising from the superconducting gap. Some of these effects are schematically shown in Fig. 1. This review explains those phenomena in detail.
The basis of the spin-valve effect, extensively studied by Meservey and Tedrow (1994), is the spin-resolved tunnelling into the superconductor with Zeeman splitting, shown in Fig. 1(a). This schematic picture illustrates how by properly tuning the voltage across the junction, the electronic transport is dominated by only one of the spin species. That results in the peculiar asymmetric differential conductance curves dI/dV (V ) = dI/dV (−V ) (see Fig. 2) which were observed in experiments. This idea has been recently used to probe the spatially resolved spin polarization of different magnetic materials by means of scanning tunnelling microscopy with Zeeman split superconducting tips (Eltschka et al., 2014(Eltschka et al., , 2015. In normal metals and superconductors a spin accumulation, or spin imbalance, can be created by injection of a charge current from a ferromagnetic electrode (Gu et al., 2002;Jedema et al., 2001;Johnson, 1994;Johnson and Silsbee, 1985;Poli et al., 2008;Shin et al., 2005;van Son et al., 1987;Takahashi and Maekawa, 2003). This state is characterized by the excess population in one of the spin subbands, determined by the balance between spin injection and relaxation or spin diffusion rates. In normal metals the nonequilibrium spin imbalance decays due to spin-flip scattering at typical distances of several hundreds of nanometers. In the superconducting state, at low temperatures k B T ∆ the injection of any amount of carriers just above the energy gap shifts the chemical potential of quasiparticles rather strongly due to the large amount of quasiparticles at the gap edge [ Fig. 1(b)]. This leads to a strong spin signal in SF junctions (Poli et al., 2008;Takahashi and Maekawa, 2003).
The spin relaxation length in normal metals depends only weakly on temperature. In the superconducting state, however, the scattering length is drastically modified. According to the first theory and experiments on spin injection in superconductors, the spin relaxation length was found to be reduced compared to the normal state (Morten et al., 2004;Poli et al., 2008). However, subsequent experiments showed, contrary to expectations, an increase of the spin decay length (Hübler et al., 2012;Quay et al., 2013;Yang et al., 2010). It is now understood that these findings can only be explained by taking into account the spin-splitting field inside the superconductor Bobkov, 2015, 2016;Kr-ishtop et al., 2015;Silaev et al., 2015a). Due to this field, as shown in Sec. III.A, it is necessary to take into account four types of nonequilibrium modes describing spin, charge, energy, and spin-energy imbalances. These modes provide the natural generalization of the charge and energy imbalances introduced by Schmid and Schön (1975). In Sec. IV we show how the spin-splitting field couples pairwise these modes: Charge to spin energy and spin to energy. Such a coupling leads to striking effects. For example, the coupling between the spin and energy modes leads to the long-range spin-accumulation observed in the experiments by Hübler et al. (2012) and Quay et al. (2013). As we show in Sec. IV this long-range effect is related to the fact that the energy mode can only relax via inelastic processes which at low temperatures are rare.
The coupling between different modes shows up also in tunnel contacts with spin-split superconductors. Because the spin-splitting field shifts the spin-resolved DOS away from the chemical potential of the superconductor, the system exhibits a strong spin-dependent electron-hole asymmetry. The spin-averaged density of states is still electron-hole symmetric, and therefore does not violate fundamental symmetries of the (quasiclassical) superconducting state. This spin-resolved electron-hole asymmetry leads to a large spin Seebeck effect shown schematically in Fig. 1(c) and discussed in Sec. V.F. A temperature difference across a tunneling interface between a normal metal and a spin-split superconductor leads to a pure spin current between the electrodes, without transport of charge. If one of the electrodes is small so that the spin injection rate is large or comparable to the rate for spin relaxation, a spin accumulation forms in this electrode.
However, it was noticed in several recent works (Kalenkov et al., 2012;Machon et al., 2013Machon et al., , 2014Ozaeta et al., 2014) that in certain situations the relevant observables are not spin-averaged, resulting in an effective electron-hole asymmetry showing up also in the charge current. The spin components are weighted differently in a setup consisting of the spin-filter junction connected to the spin-split superconductor (Ozaeta et al., 2014), shown schematically in Fig. 1(d). As a result of this effective electron-hole symmetry breaking, the system exhibits a very large thermoelectric effect. This is discussed in Sec. V.
The main body of the review is organized as follows. In Sec. II we discuss different mechanisms causing the spinsplitting field in superconductors, explain the effects of spin-flip and spin-orbit scattering, and describe shortly different types of equilibrium and linear-response effects resulting from the spin-splitting field. In Sec. III we introduce the quasiclassical formalism to describe superconducting systems driven out of equilibrium. In particular we introduce there the different nonequilibrium modes, the different currents related to them, and the transport through hybrid interfaces. We also discuss the relaxation mechanisms of different nonequilibrium modes. Section IV focuses on the spin injection and diffusion in superconducting systems. We review different experiments on the detection of spin and charge imbalance, and discuss the long-range spin accumulation and the predicted behavior of the spin Hanle effect in superconductors. In Sec. V we first give an overall review of thermoelectric effects in superconductors, and then concentrate on describing the giant thermoelectric response of a system exhibiting spin-polarized tunneling into a superconductor with a spin-splitting field. In all those sections we focus on stationary phenomena. In contrast, in Sec. VI we consider the case of time-dependent fields and currents, which is another method of driving the systems out of equilibrium, besides stationary voltages and currents. In particular, we discuss magnetic resonance and current rectification effects, and the microwave-induced modifications of the superconducting state, especially in the presence of the spin-splitting fields. Finally, we present our conclusions and an outlook on possible future developments in the field in Sec. VII. This review is written by four theorists, and therefore it has a strong focus on the proper theoretical description of the various effects. However, we have attempted to serve also those readers wishing to get an overall picture of the involved physics and a historical outline of the different developments. Therefore, each section starts with such an outline, without dwelling on the theoretical details. We have also aimed at separating the formal descriptions from the overall picture of the physical effects. Although this article is a review, we have included several new results, especially with regard to the thermoelectric effects in Sec. V, and reformulated results in the theoretical framework used in the review, as in Sec. VI. Simultaneously we have tried to give all necessary details needed to carry out the theoretical calculations presented in the review. We hope that they serve as a starting point for those who want to expand the theory to completely new realms.

II. SUPERCONDUCTOR WITH AN EXCHANGE FIELD
In this section we describe the main features of superconductors with a spin-split density of states. This is one of the main focus points of this review. Such a splitting can originate either by an adjacent ferromagnetic insulator or by an external magnetic field. We start by discussing these two situations.

A. Spin splitting caused by an external magnetic field
The response of a superconductor to a magnetic field is twofold: On the one hand the field couples to the or-6 bital motion of the electrons and creates circulating currents (Meissner effect) (Meissner and Ochsenfeld, 1933) that try to expel the field from the bulk (diamagnetic response). By increasing the amplitude of the applied field superconductivity is gradually reduced. We denote below this mechanism of suppression of superconductivity by the orbital depairing effect. This mechanism dominates in bulk samples or in thin films with the field applied perpendicular to the plane of the film. If the amplitude of the applied field reaches a critical value H c , the created currents increase the free energy such that the system undergoes a transition into the normal state (Tinkham, 1996).
On the other hand the magnetic field may also couple to the spin of the conduction electrons (Zeeman effect). This effect tries to align the spin of the electrons and acts as a paramagnetic depairing mechanism that also leads to the suppression of superconductivity. The spin paramagnetic effect dominates with respect to the orbital one in superconducting films with thickness smaller than the London penetration length. The critical value H c for a magnetic field applied in the plane of the film largely exceeds the critical value H c⊥ of a perpendicular field Tinkham, 1996). In this case the magnetic field penetrates uniformly the film, screening currents are relatively small and therefore H c is limited by the spin paramagnetic effect that tries to align the spin of the original singlet Cooper pairs, as demonstrated by Clogston (1962) and Chandrasekhar (1962). At T = 0, and in the absence of spin-orbit coupling and magnetic impurities, the critical field due to the paramagnetic effect is given by H cP = ∆ 0 /( √ 2µ B ), where ∆ 0 is the superconducting gap at zero field and temperature and µ B is the Bohr magneton. For conventional BCS superconductors with critical temperature ranging from 1 to 10 K, H cP can be of the order of several Tesla and hence superconductivity exists within a large range of field strengths.
The paramagnetic effect in thin superconducting films leads to a Zeeman splitting of the density of states (DOS), equal to 2µ B |H| as first observed by Meservey et al. (1970) in Al films, and later in several other works (Hao et al., 1990;Meservey et al., 1980;Paraskevopoulos et al., 1977;Tedrow and Meservey, 1971;Xiong et al., 2011). The splitting can be detected by measuring the tunneling conductance of the superconductor with a normal or a ferromagnetic probe. For example, if one uses a normal metal (N) tunnel coupled to a superconducting film S, the current through the normal-insulator-superconductor (NIS) junction is determined by the tunneling expression (Giaever and Megerle, 1961)  where N and N s are the reduced density of states of the normal metal and the superconductor, respectively, n F the Fermi distribution function and G T is the conductance of the tunneling barrier I which we first assume to be spin-independent. Electronic transport occurs mainly at energies close to the Fermi level where the density of states of the normal probe N (ε) can be accurately approximated by a function independent of the energy. In such a case, and at low temperatures the measured differential conductance G = dI/dV is proportional to the density of the states of the superconductor. If the normal metal probe is non-magnetic, i.e., electrons are not spinpolarized, the measured G shows the spin-split peaks and is symmetric with respect to the voltage as sketched in Fig. 2(a). If the tunneling barrier has a spin-dependent transmission, i.e., a spin filter, or a ferromagnetic electrode with a non-vanishing spin polarization, the observed peaks are asymmetric (Paraskevopoulos et al., 1977) [see Fig. 2 (b)]. The asymmetry is proportional to the spin polarization of the conduction electrons of the electrode or the spin-filter efficiency of the barrier. Indeed spinsplit superconductor-ferromagnet bilayers have been used for determining the spin polarization of magnetic metals (Meservey et al., 1980;Paraskevopoulos et al., 1977;Tedrow and Meservey, 1971).
The normalized total density of states of the spin-split superconductor, sketched in Fig. 2 (a), can be written as the sum of the DOS of each spin species, 2 2 As conventional, here and below ε 2 − |∆| 2 means sgn(Re ε) ε 2 − |∆| 2 in terms of the principal branch Re √ · ≥ 0 of the square root.
where ±h eff is the effective spin-splitting field. Equation (2) for the DOS is a simplified description of spin-split superconductors, which is accurate, as we see below, only in certain limiting cases. The expression does not take into account the effect of magnetic impurities or spin-orbit coupling (SOC). For example, magnetic impurities suppress superconductivity and eventually lead to a gap-less situation (Abrikosov and Gor'kov, 1960a). On the other hand the SOC counteracts the magnetic impurities but it may also lead to a broadening of the coherence peaks in the spectrum (Fulde, 1973;Maki, 1966;Meservey et al., 1975;Werthamer et al., 1966). Moreover, if the spin-splitting field h eff in Eq. (2) is due to an external magnetic field, also orbital depairing contributes to the pair breaking.
We present a more quantitative analysis of these effects below, in Sec. II.D. But first we discuss another way of inducing a spin-split DOS in a superconductor by means of the magnetic proximity effect.

B. Superconductor-ferromagnetic insulator structures
An alternative way of creating a spin-split density of states in a superconductor is via the exchange field induced by an adjacent ferromagnetic insulator (FI). This may lead to a large spin splitting, of the order of the superconducting gap, even without an externally applied magnetic field. The first evidence of an exchange field induced in a superconductor via the magnetic proximity effect was observed in a thin Al film in contact with an EuO film (Tedrow et al., 1986) and with EuS (Hao et al., 1990;Wolf et al., 2014b;Xiong et al., 2011). The equivalent internal field of the Zeeman splitting in Al-EuS systems has been reported to be as large as 5 T (Moodera et al., 2007). Another S-FI combination explored so far is GdN-NbN (Senapati et al., 2011). This system shows, however, a weaker and less clear spin splitting than the Al-europium chalcogenide combination, due to the influence of the larger spin-orbit interaction in Nb compared to the Al-based devices.
The advantage of using a FI instead of an external magnetic field is that one avoids the depairing effects and all complications caused by the need to apply magnetic fields in superconducting devices. Moreover, because the electrons of the superconductor cannot penetrate into the insulating FI, the superconducting properties are only modified by the induced spin-splitting field at the S/FI interface, and not by the leakage of Cooper pairs into the FI. The latter would be relevant in the case of metallic ferromagnets. Finally, FIs can also be used as spinfilter barriers (Moodera et al., 2007), in some cases with a very high spin-filtering efficiency, and therefore they will play a crucial role for several of the applications discussed in subsequent sections which require high spin-filter efficiency. For all these reasons the emphasis of this review FIG. 3 Color plot of measured differential conductance, dI/dV of a EuS/Al/Al2O3/Al junction as a function of the applied voltage and external magnetic field. Hco denotes the coercive field of the EuS layer when the magnetization switches. Courtesy of Elia Strambini, Francesco Giazotto (SNS Pisa), and Jagedeesh Moodera (MIT).
is on S-FI structures.
In Fig. 3 we show an example of the measured differential conductance of an EuS/Al/Al 2 O 3 /Al junction. The Al layer adjacent to the EuS has a spin-split density of states that shows up as the splitting peaks (bright stripes in the figure) in dI/dV . Even at zero applied magnetic field the splitting is nonzero. The magnetization reversal of EuS at H c ≈ −18.5 mT manifests as a discontinuity of the conductance peaks.
The amplitude of the spin splitting in a FI/S structure depends on both the intrinsic properties of the superconductor, such as the amount of magnetic impurities and the strength of spin-orbit coupling, and on the quality of the S/FI interface. The latter is crucial for obtaining a large splitting, as shown for example by Hao et al. (1990), on EuS/Al systems. In Table I we show a list of FI/S combinations and the reported induced exchange splittings and spin-filter efficiencies (barrier spin polarizations).
The spin splitting caused by a ferromagnetic insulator on an adjacent superconductor can be understood with the help of the following model (Izyumov et al., 2002;Khusainov, 1996). The effective Hamiltonian describing a FI/S bilayer consists of the ferromagnetic coupling within the FI layer, Middle column shows the spin-filter efficiency characterized by the polarization P = (G ↑ − G ↓ )/(G ↑ + G ↓ ) of the FI barrier (red) with normal-state conductance Gσ for spin σ. The exchange splittings measured in the superconductor (blue) are listed in the right column. The data is extracted from 1 (Tedrow et al., 1986); 2 (Moodera et al., 1988); 3 (Hao et al., 1990); 4 (Moodera et al., 1993); 5 (Senapati et al., 2011); 6 (Pal and Blamire, 2015). Note that µB·1 T=58 µeV ∼ = 670 mK. where J r,r is the exchange coupling between the localized spins S r . The term H an describes the magnetic anisotropy that in films is usually in the plane of the film.
The S layer we describe by the usual BCS Hamiltonian, whereas for the FI/S interface we consider a model of an ensemble of localized magnetic moments that interact with the spin of the conduction electrons of the superconductor via exchange interaction, Here the symbolΣ means that we consider only the magnetic moments S r localized at the interface. J ex is an effective parameter which for example describes the s-d or s-f exchange interaction. We consider a ferromagnetic insulator with a Curie temperature much larger than the superconducting transition temperature. Thus we can assume that the magnetization of the FI is only determined by the Hamiltonian in Eq. (3) and not affected by the superconducting state (Bergeret et al., 2000;Buzdin and Bulaevskii, 1988). With this assumption Hamiltonian (4) describes conduction electrons interacting with an effective exchange field proportional to the average local spin S r at the FI/S interface. This average can be computed by solving independently the magnetic Hamiltonian (3) (Strambini et al., 2017). In the superconducting state the localized exchange field leads to a modified density of states characterized by the splitting of the coherent BCS peaks. This modification of the spectral properties of the superconductor is non-local and survives over distances away from the FI/S interface of the order of the coherence length ξ s (Bergeret et al., 2004;Tokuyasu et al., 1988). If the thickness d of the S film is much smaller than ξ s , the spin splitting can be assumed as homogeneous across the film. Thus the density of states is given by Eq. (2) with an effective exchange field h ef f ≈ J ex S r a/d (de Gennes, 1966a;Khusainov, 1996;Tokuyasu et al., 1988), where a is the characteristic distance between the localized spins.
In this review we mainly focus on thin S films and hence we adopt this approximation for the description of the uniform magnetic proximity effect when the film is adjacent to a FI layer. For an inhomogeneous magnetic configuration of the FI layer the spatial scale that determines the effective splitting is the superconducting coherence length. Thus the splitting becomes observable if the ferromagnet consists of magnetic domains with sizes larger than the superconducting coherence length. This is for example the case of EuS films with magnetic domains of micrometer size (Tischer, 1973) that explain the spin splitting observed in experiments (Hao et al., 1990) on Al-EuS structures, where the Al layers are a few nanometers thick and with a coherence length of the order of 100 nm. In contrast, a spatially fast changing magnetization averages to zero and results in a vanishing effective splitting (Strambini et al., 2017).
It also follows from this description that the strength of the spin splitting depends crucially on the quality of the FI/S interface. The growth of any non-magnetic oxide between the FI and S layers drastically suppresses the effective exchange interaction J ex and hence reduces the spin splitting (Hao et al., 1990). In addition to the quality of the interfaces, the spin splitting also depends on intrinsic properties of the superconducting film. For example, magnetic disorder may lead to a strong suppression of superconductivity in the S film (Abrikosov and Gor'kov, 1960a). But even in the absence of magnetic disorder and magnetic impurities, spin-orbit coupling may lead to strong modifications in the DOS. In particular, large spin-orbit scattering rates lead to a rounding of the BCS peaks and to a less sharp spin splitting. This explains why splitting has been observed in materials with relative small atomic number Z, such as Al, Be, V, but not in heavier materials such as Pb. For a detailed description of the intrinsic properties of different superconductors in this respect, we refer the reader to the excellent review by Meservey and Tedrow (1994).
In order to move forward and describe theoretically the transport properties of diffusive hybrid structures with spin-split superconductors we introduce in the next section the quasiclassical Green's function formalism for superconductors in the presence of spin-dependent fields and spin-polarised interfaces.
C. Diffusive-like equation for superconductors in the presence of spin-dependent fields We mainly focus on diffusive superconductors, for which the elastic mean free path l, due to scattering at non-magnetic impurities, is much smaller than any other length involved in the problem. 3 This condition is equivalent to the assumption that the inverse of the elastic relaxation time 1/τ is much larger than all characteristic energies. 4 It was shown by Usadel (1970) that in such a limiting case the superconducting properties can be described by a diffusion-like equation, now called the Usadel equation, which describes the behavior of the Keldysh Green's function (GF) Here D is the diffusion coefficient,Σ =Σ(r, ε) is the general self-energy term andǧ(r, ε) is the (momentum isotropic) quasiclassical Green's function. The latter two depend on the spatial coordinate and the energy. Their structure in the Keldysh space can be represented via three components. For example, where the retarded (R), advanced (A) and Keldysh (K) components are 4×4 matrices in the Nambu ⊗ spin space. The quasiclassical Green's functions are obtained by integrating over the quasiparticle energy the microscopic Green's functions constructed from the Nambu bi-spinor For instance, in the absence of gradients in the superconducting order parameter we can fix the general form of R and A functions tô The matrices τ j , j = 1, 2, 3 in Eqs. (5) and (8), are Pauli matrices in the Nambu (or electron-hole) space. 5 The offdiagonal components, i.e., those proportional to τ 1 , are characteristic for the superconducting state and describe the anomalous GFs. The matrices σ j (j = 0, 1, 2, 3), are the usual Pauli matrices in spin space, with σ 0 denoting the unit matrix. The component of the anomalous GFs proportional to σ 0 describes the usual singlet correlations in the BCS theory. On the other hand, the terms proportional to the Pauli matrices σ j , j = 1, 2, 3, describe the three triplet components of the condensate that appear in the presence of spin-dependent fields (Bergeret et al., 2001b(Bergeret et al., , 2005. The real and imaginary parts of the different components of the retarded Green's function are plotted in Fig. 4 in the case of homogeneous magnetization pointing in the z direction. In that case, the GF components proportional to σ 1,2 vanish. In Eq. (5), h denotes the spin-splitting field either generated by an external field or generated by the magnetic proximity effect in a FI/S junction, and∆ is the order parameter matrix∆ = ∆τ ↑ + ∆ * τ ↓ , diagonal in Keldysh 10 and spin space. ∆ is the superconducting gap that has to be determined self-consistently from where λ is the effective coupling constant and Ω D is the Debye cutoff frequency. The self-energy in Eq. (5) can be written as the sum of elastic and inelastic contributions,Σ =Σ el +Σ in . The elastic partΣ el =Σ so +Σ sf +Σ orb consists of three terms due to scattering at impurities with spin-orbit couplinǧ Σ so , the spin flips at magnetic impuritiesΣ sf and the orbital depairing due to the applied magnetic field Σ orb . Within the Born approximation these components are given byΣ The characteristic relaxation times are defined as τ so , τ sf and τ orb . The orbital depairing leads to the suppression of superconductivity in thin films by the in-plane magnetic field (de Gennes, 1999). It does not enter explicitly the kinetic equation for the distribution function, but modifies indirectly its energy dependence by affecting the spectral properties of the superconductor. In the case of a thin magnetic film with thickness d, and an in-plane applied magnetic field the orbital depairing relaxation time is given by (Belzig et al., 1996) 1/τ orb = De 2 B 2 d 2 /6. In contrast, the spin relaxation self-energies,Σ so ,Σ sf , not only change the coefficients in kinetic equations by affecting the spectral properties, but also provide additional relaxation terms for the non-equilibrium modes as discussed in detail below. We describe these relaxation mechanism in terms of two parameters: the total spin relaxation time, τ −1 sn = τ −1 so + τ −1 sf and the relative strength of the two scattering mechanisms β = (τ so − τ sf )/(τ so + τ sf ). Spin-flip scattering dominates for β ≈ 1, and spin-orbit scattering for β ≈ −1. These parameters are material and to some extent sample specific. For example, in thin Al films, the reported values from a set of spin injection experiments are τ sn ≈ 100 ps (Jedema et al., 2002;Poli et al., 2008) and β ≈ 0.5 indicating the dominance of spin-flip relaxation over spin-orbit scattering. According to (Abrikosov and Gor'kov, 1962;Meservey and Tedrow, 1978), the spinorbit scattering rate is the momentum relaxation rate times (Zα) 4 , where Z is the atomic number and α is the fine-structure constant. Therefore, the spin-orbit scattering rate grows rapidly as the atomic number grows. For example, the reported value of τ sn in Nb is only 0.2 ps, and is strongly dominated by spin-orbit scattering (Wakamura et al., 2014). Relevant inelastic processes, particle-phonon and particle-particle collisions, are described by the inelastic self-energyΣ in =Σ eph +Σ ee . These processes do not conserve the energies of colliding quasiparticles, but conserve the total spin. A model self-energy for inelastic relaxation due to electron-phonon scattering is given by (Eliashberg, 1972) are parts of Fermi surface averages of the free phonon propagators, T ph is the phonon temperature and we denote X RA = X R − X A for X = D,ĝ. Inelastic particle-particle collision self-energiesΣ ee within the quasiclassical theory have been discussed by Eliashberg (1972) and Serene and Rainer (1983). Equation (5) is supplemented by the normalization conditionǧ We do not present here the derivation of the Usadel equation Eq. (5), that can be found in several previous works (Belzig et al., 1999;Bergeret et al., 2005;Eilenberger, 1968;Usadel, 1970). For the readers not familiar with the quasiclassical method it is sufficient to understand that the knowledge of the GF components, Eq. (6), satisfying Eq. (5) shed light on the spectral properties via the retarded and advanced components and on transport properties via the Keldysh component. In particular, the self-energies provide the collision integrals of the transport theory for the different scattering processes, as described in more detail in Sec. III.
As an example solution of the Usadel equation, we analyze next the spectral properties of a superconducting film with an intrinsic exchange field.

D. Superconducting film in a homogeneous spin-splitting field
By combining the retarded part of the Usadel equation, Eq. (5), with the self-consistency equation, Eq. (9), one can find the equilibrium value of the superconducting order parameter ∆ as a function of the parameters of the model, such as temperature, exchange field or spinorbit/spin-flip scattering rates. This analysis also reveals the critical parameter values above which the material ceases to be superconducting. The results of such an analysis are presented in Sec. II.F. Additionally, superconductivity can also be suppressed by driving the system out of equilibrium as discussed in Sec. III.H and VI.B. We denote by ∆ the self-consistent gap, whereas ∆ 0 is the pair potential at T = h = 0, in the absence of extra scattering, i.e.,Σ el =Σ in = 0. The corresponding zero-exchange field critical temperature is denoted by T c0 = ∆ 0 /1.76.
In order to determine the spin-averaged density of states (DOS) of a superconductor in equilibrium within the quasiclassical approach, we only need to compute the R and A components of the GF, Eq. (6), which have the general form Eq. (8). From the knowledge of these functions we can determine the reduced density of states relative to the normal state, is the difference of the DOS of the spin-split bands.
If we assume that the exchange or Zeeman field induced in S is homogeneous and points in z-direction then the spectral GF, Eq. (8), reduces tô The advanced GF can be obtain from g A = −τ 3 g R † τ 3 . In a homogeneous situation one can neglect the derivatives in Eq. (5) and g R is determined from the algebraic matrix equation and the normalization condition Eq. (14). We have neglected the self-energy due to orbital effects by assuming, for example, that spin splitting h is induced by the proximity with a FI. If we also neglect the relaxation terms, σ so and σ sf , it is straightforward to show that the total DOS coincides with Eq. (2). Inclusion of magnetic impurities and spin-orbit scattering does not allow for an analytical closed-form expression of the DOS. The different GF components calculated numerically are shown in Fig. 4. In Fig. 5 we show the density of states calculated numerically for T → 0 and h = 0.4∆ 0 . In Fig. 5(a) we assume no spin-orbit coupling, β = 1, and plot the DOS for different values of τ sn . The effect of magnetic impurities is to suppress drastically the superconducting properties. By increasing the number of magnetic impurities, i.e., by decreasing the value of 1/τ sf , one clearly sees that the superconducting gap decreases and the coherent peaks at energies close to the gap are broadened. For a relatively short spin-flip time τ sf < τ sf,cr (h, T ) the energy gap vanishes (Abrikosov and Gor'kov, 1960a), e.g. τ sf,cr ∆(h = 0, T = 0) = 3/4 and τ sf,cr ∆(h = 0.4∆ 0 , T = 0) ≈ 0.44 In Fig. 5(b) we show the DOS for different values of the spin-orbit scattering time τ so ∆ 0 , when there is no spin-flip scattering, β = −1. The overall consequence of the spin-orbit coupling is to counteract the effect of the exchange field in the density of states. As one sees in the figure, by decreasing τ so ∆ 0 the density of states tends to the DOS in the h = 0 case. Notice, however, that in general the spin-orbit coupling also causes a broadening of the coherent peaks in the case that h = 0, while it has no effect on the DOS when h = 0.

E. Static spin susceptibility of conventional superconductors
Besides broadening of the DOS singularities, the spindependent scattering has an important effect on the paramagnetic spin susceptibility of the superconductor as it becomes non-vanishing even in the zero-temperature limit (Abrikosov and Gor'kov, 1960b;Yosida, 1958). In general the spin susceptibility χ can be found by calculating the spin polarization S = −χh generated by a small spin-splitting field h. In the superconducting state the susceptibility normalized with respect to that of the normal state χ n can be expressed in terms of the difference between spin-up and spin-down DOS N − = N ↑ − N ↓ and the equilibrium distribution function f eq = tanh(ε/2T ) (Abrikosov and Gor'kov, 1960b;Bergeret et al., 2005), Using Eq. (2) for the spin-split DOS in the absence of magnetic depairing processes one can see that Therefore for T → 0 and in the absence of spindependent scattering χ(T = 0) = 0. This is a con-sequence of the lack of polarizability of the condensate consisting of spin-singlet Cooper pairs. However, if either the spin-orbit and/or the spin-flip scattering rate is non-zero, the superconducting condensate may exhibit a finite paramagnetic susceptibility for T → 0. For a pure spin-orbit relaxation this result was obtained first by Abrikosov and Gor'kov (1960b).
A generalization of this result can be obtained by solving Eqs. (5,10,11) to find the DOS difference related to the GF component N − = Reg R 33 and substituting into Eq. (18). In Fig. 6 we show the low-temperature dependence of the spin susceptibility as a function of the normal-state spin scattering rate 1/τ sn and the parameter β determining the relative strength of spin-flip and spin-orbit scattering.
From these plots one can see that both scattering mechanisms result in a finite susceptibility at low temperatures. Notice that only the spin-flip scattering generates a strong suppression of the order parameter (see Fig. 7). Therefore, at T T c0 the growth of χ as a function of the rate 1/τ sn towards the normal-state value is much faster when the spin-flip scattering dominates the spinrelaxation in the normal state, i.e. when β > 0. The nonzero susceptibility explains the observed Knight shifts in ordinary superconductors (Meservey and Tedrow, 1978).
The static spin susceptibility (18) characterizes the paramagnetic response of the superconductor to the external magnetic field. This expression cannot be used for calculating the non-equilibrium spin accumulation induced by the spin-dependent chemical potential shift ±δµ when the distribution functions in different spin subbands are given by f ↑ (E) = f 0 (E + δµ) and f ↓ (E) = f 0 (E − δµ). The theory of non-equilibrium spin states in superconductors is explained in detail in Sec. IV. In particular, there we explicitly show [see the discussion after Eqs. (59,60)] that the non-equilibrium spin accumulation generated by δµ is always exponentially suppressed at low temperatures provided that there is a non-zero gap in the quasiparticle spectrum. This is the case even in the presence of the strong spin-orbit interaction [cf. Fig. 5(b)] although the static spin susceptibility (18) is non-zero at T = 0.

F. Other effects related to the Zeeman and exchange fields
The results obtained in the previous section are used throughout this review for the description of superconductors with spin-split density of states. It is however important to mention other effects that we are not considering but that may appear in the presence of Zeeman and exchange fields. FIG. 6 Paramagnetic susceptibility of a thin superconducting film for different values of β as a function of (a) the normalstate spin scattering rate 1/τsnTc0 and (b) the temperature. In (a) T = 0.1Tc0 and in (b) 1/τsn = 0.5Tc0, where Tc0 is the critical temperature in the absence of spin scattering.

Fermi liquid effects
The equations presented in Sec. II.C-II.D assume noninteracting electrons, which is a rather good approximation in metals. Corrections due to interactions can be incorporated by renormalizing certain parameters of the theory with the help of the powerful Landau's Fermiliquid theory (Landau, 1956). In the context of the quasiclassical theory such effects have been studied by Alexander et al. (1985). In particular, and important for our work, is the renormalization of the spin-splitting field amplitude h ≡ |h| entering Eq. (5) which according to the Fermi liquid theory is given by (Catelani et al., 2008;Xiong et al., 2011) : where H ext is the external magnetic field, H F I is the field induced by a FI, and G 0 is the effective antisymmetric Fermi-liquid parameter of Landau's theory. The value of G 0 in the superconducting state depends on temperature, is negligibly small at low temperatures (T T c ) and finite for temperatures close to the superconducting critical temperature T c and above it. In what follows we neglect the Fermi-liquid corrections by assuming that G 0 1. This is justified at the low temperatures where majority of our work concentrates. However, when describing equilibrium properties of spin-split superconductors at temperatures close to T c , the Fermiliquid corrections to the field h in Eq. (5) should also be included for some superconductors. For example, for Al films G 0 ∼ 0.16−0.3 (Catelani et al., 2008;Tedrow et al., 1984) close to T c and above it, whereas in V it is negligibly small at all temperatures (Gibson and Meservey, 1989). 13 2. Paramagnetic depairing mechanisms and the FFLO state When we calculate the DOS of a superconductor in an exchange field in Sec. II.D, we assume that the ground state of the S layer corresponds to a spatially homogeneous order parameter ∆. The value of the order parameter is determined from the self-consistency equation, written in terms of the Matsubara frequencies ω n = (2n + 1)πT and a free-energy functional discussed in Appendix A.
In the absence of spin relaxation the results coincide with those of Maki and Tsuneto (1964): There is a critical value of temperature, T 0 ≈ 0.556T c0 , above which the transition, driven by the exchange field, between the normal and superconducting states is of the second order. For T < T 0 , the transition is of the first order. At low temperatures the critical field of paramagnetic depairing is known as the Chandrasekhar-Clogston limit (Chandrasekhar, 1962;Clogston, 1962) and is given by h = ∆ 0 / √ 2, where ∆ 0 = 1.76T c0 is the order parameter at T = 0 without depairing.
A non-zero spin relaxation time leads to a quantitative modification of this behavior (Bruno and Schwartz, 1973). In Fig. 7 we show the computed ∆ as a function of temperature and exchange field for a normal-state spin relaxation rate 1/T c0 τ sn = 0.96, where T c0 is the critical temperature in the absence of gap-suppressing effects. Panels (a) and (b) in Fig. 7 show the effects of spin-flip β = 1 and spin-orbit β = −1 relaxation, respectively. The phase-transition line, ∆(T, h) = 0, is shown in Fig. 7 by the solid curve in the (T, h) plane. The red part on this curve corresponds to the second-order transition where ∆(h) continuously goes to zero with increasing h. This behavior changes to the abrupt first-order transition at T < T 0 . The first-order transition line is shown by circles in the (T, h) plane and it is different from the ∆(T, h) = 0 curve shown by the green line at T < T 0 . Both relaxation mechanisms reduce the range of temperatures for which the first-order transition takes place. In other words, the threshold temperature T 0 between firstand second-order phase transitions is suppressed as compared to the case without relaxation.
In other respects the modification of the superconducting state strongly depends on the spin relaxation mechanism. Spin-flip scattering breaks the time-reversal symmetry and therefore leads to a strong suppression of the superconducting gap and T c . On the other hand, spinorbit interaction is time-reversal invariant and keeps the T c intact.
The striking effect of spin-orbit scattering is to enhance the critical field of the Chandrasekhar-Clogston limit (Bruno and Schwartz, 1973). This tendency is shown in Fig. 8, where we compare the dependencies of the superconducting gap on h and relaxation rates 1/τ sn for (a) spin-flip and (b) spin-orbit relaxation, respectively. These two cases are characterized by the opposite behaviors of the critical field h c as a function of 1/τ sn -it is suppressed by spin-flip scattering and enhanced by the spin-orbit one. This effect can be understood from the comparison to the changes in the density of states caused by the two spin-relaxation mechanisms (see Fig. 5): spin-flip scattering primarily lifts the gap in the density of states, whereas spin-orbit scattering acts to nullify the spin splitting without affecting the gap. In both cases there is a threshold value of the scattering rate τ −1 0 when the superconducting phase transition changes from the second order at τ −1 sn > τ −1 0 to the first order at τ −1 sn < τ −1 0 . The assumption of a homogeneous order parameter does, however, not always correspond to the lowestenergy state of the system. Fulde and Ferrell (1964) and Larkin and Ovchinnikov (1965) predicted that the exchange or Zeeman field can cause a transition to an inhomogeneous superconducting state, with a spatial periodical modulation of the order parameter at the scale of the coherence length. Such an FFLO state exists for temperatures below T 0 and only for exchange fields that satisfy 0.71∆ 0 < h < 0.755∆ 0 when T → 0.
The FFLO state has not been observed in conventional superconductors. One of the main reasons for this is that it is very sensitive to disorder (Aslamazov, 1969) and spin-orbit coupling and hence it is expected to occur only in extremely clean samples. In the present review we mainly focus on diffusive systems and therefore we do not pay further attention to the FFLO state. Moreover, the effective exchange fields we mostly consider are far below the range for the FFLO state to occur.

Long-range triplet correlations
In this review we mainly focus on the transport properties of the superconductor in S/F structures. For completeness, however, we briefly discuss in this section the appearance of triplet correlations in ferromagnetic metals induced by the superconducting proximity effect (Bergeret et al., 2001b(Bergeret et al., , 2005. Such correlations may induce non-dissipative spin-polarized currents that could be particularly useful to lower the energy consumption in spintronic devices (Eschrig, 2011(Eschrig, , 2015Linder and Robinson, 2015).
Superconducting correlations can be induced in nonsuperconducting (N) metals by means of the proximity effect (de Gennes, 1964). This can be seen as the "leakage" of Cooper pairs from the superconductor into the normal region. The microscopic origin of this "leakage" is the Andreev reflection (Andreev, 1964), where an incident electron from N is retro-reflected as a hole forming a Cooper pair in S. An important consequence of the Andreev reflection is that the electron-hole pair induced in N leads to a formation of a non-vanishing pair amplitude in N within a characteristic distance ξ N away from the S-N interface. This distance depends on the intrinsic properties of N, such as the degree of disorder, spin-dependent fields, etc. For example, in a diffusive normal metal without magnetic impurities this length is of the order of the thermal length D/T , where D is the diffusion coefficient.
In a superconductor-ferromagnet (S-F) junction the situation is very different. In usual ferromagnets the intrinsic exchange field h is much larger than the superconducting characteristic energies. The proximity-induced superconducting correlations decay and oscillate over the length ξ F = Re D/h The oscillation of the condensate leads to the well-established 0-π transition in S-F-S Josephson junctions predicted by Buzdin et al. (1982) and observed much later in experiments by Ryazanov et al. (2001) and Kontos et al. (2001) on SFS structures with ferromagnetic alloys with a weak exchange field. In conventional ferromagnets like Fe, Ni or Co, the effective exchange fields is large (∼ 0.1 − 1 eV) and hence ξ F is very short (Verbanck et al., 1994). This explains why the Josephson effect in S-F-S junctions with transition metals has been only observed for thickness of the F layer smaller than 10 nm (Blum et al., 2002;Robinson et al., 2006).
The situation is very different if the exchange field is spatially non-homogeneous, as for example in a magnetic domain wall. In such a case, as shown by Bergeret et al. (2001b), triplet components with equal spin components can generated. These are insensitive to the local exchange field and hence can penetrate the F metal over the thermal length. Such long-range triplet component of the superconducting condensate explains the Josephson long-range currents observed in S-F-S junctions made for example of half-metal (Anwar et al., 2012;Keizer et al., 2006b), ferromagnetic multilayers (Khaire et al., 2010), and ferromagnets with intrinsic inhomogeneous magnetization (Robinson et al., 2010).
We emphasize that the physics of the triplet component addressed in this section describes the equilibrium state of ferromagnets in contact with superconducting electrodes. The main focus of the present review is however the non-equilibrium properties of superconductors in spin-splitting fields, either generated by an adjacent FI or by an external field. As discussed before, in order to preserve the superconducting state, those fields have to be smaller than the superconducting gap, and hence according to Eq. (5), all components of the condensate, singlet and triplets, vary over similar length scales. This in particular means that no distinction between shortand long-range has to be made.

Cryptoferromagnetic state
The way we model the exchange field generated in a superconducting thin film by an adjacent ferromagneticinsulating layer, see Sec. II.B, is based on the assumption that the magnetic ordering of the FI is the same in both the normal and superconducting states of the S layer. This requires that the effective exchange constants J and J ex in Eqs. (3) and (4), respectively, satisfy J J ex . In this case the ferromagnetic order is robust against the exchange interaction with the conduction electrons and the superconducting condensate. Roughly speaking, this corresponds to the case when the Curie temperature T Curie of the FI layer is larger than the effective exchange field, whose maximum value can be set equal to the superconducting order parameter ∆. For example in EuS/Al junctions T Curie ∼ 16.6 K and ∆ 0 ∼ 0.25 meV, thus k B T Curie /∆ 0 ∼ 5.3. In other FIs like EuO or NbN, the Curie temperature is even larger and hence this ratio increases. Thus the approximation made in Sec. II.B is well justified.
In the case of a weak magnetic stiffness of the FI compared to the characteristic superconducting energy ∆, one should take into account the effect of the superconducting condensate on the magnetic ordering mediated by the Hamiltonian in Eq. (4). The competition between superconducting and magnetic ordering in such a case has been first considered by Anderson and Suhl (1959). They considered a ferromagnetic state mediated via the Ruderman-Kittel-Kasuya-Yosida (RKKY) indirect exchange interaction. Above the superconducting critical temperature T c the ground state of the system is a homogeneous magnetic state. When the temperature is lowered below T c the ground state of the system corresponds to a spatially inhomogeneous magnetic structure characterized by a modulation with a wave vector of the order of (λ 2 F ξ s ) 1/3 , where λ F is the Fermi wavelength and ξ s is the superconducting coherence length. This state was called by Anderson and Suhl the cryptoferromagnetic state and may describe the situation of certain ternary rare earth compounds (Bulaevskii et al., 1985). Such a magnetic modulation was also predicted for thin metallic (Bergeret et al., 2000;Buzdin and Bulaevskii, 1988) and insulating (Izyumov et al., 2002) ferromagnetic films on top of bulk superconductors.
Although in this review we assume that the condition k B T Curie /∆ 0 1 is satisfied, and hence no cryptoferromagnetic state is induced, it is important to keep in mind that such a state may exist in the case of FIs with weak magnetic stiffness.

III. NONEQUILIBRIUM MODES AND CURRENTS IN SUPERCONDUCTORS
In the subsequent sections we study transport properties of superconducting systems in an out-of-equilibrium state. Formally this means that we analyze the full Keldysh structure of the equations presented in Sec. II.C. Such a state can be characterized by the presence of nonequilibrium modes. For example, injection of an electric current from a normal electrode into a superconductor generates a charge imbalance mode (Clarke, 1972;Hübler et al., 2010;Tinkham, 1972;Tinkham and Clarke, 1972;Yagi, 2006) that diffuses into the S region. This nonequilibrium mode reflects an imbalance of the quasiparticle population between the electron-like and holelike spectrum branches. The charge imbalance measurements made in the 1970s were to our knowledge the first to study such nonequilibrium modes in non-local multiterminal settings. This technique was later adapted to spintronics, to study the nonequilibrium spin accumulation induced by spin-polarized electrodes (Johnson and Silsbee, 1985).
In the absence of spin-dependent fields there is one more nonequilibrium mode that can be excited in a superconductor besides charge imbalance. It describes the excess energy stemming from an equal change in the quasiparticle populations of the electron-like and holelike branches. This energy mode affects charge transport properties indirectly via the self-consistency equation for ∆. This mechanism explains, for example, the enhancement of the superconducting transition temperature in the presence of a microwave field (Ivlev et al., 1973;Klapwijk et al., 1977).
In this section we generalize the description in terms of nonequilibrium modes to account for superconductors with spin-split density of states. For this it is useful first to analyze in detail the full Keldysh structure of the Usadel equation, Eq. (5). We show that in the presence of a Zeeman or exchange field, it is possible to excite the spin-dependent components of the nonequilibrium distribution function in the superconductor. In a second stage we give a quantitative description of the nonequilibrium modes associated to each of these components. In Sec. IV we apply the formalism to describe the spin injection in superconductors.

A. Quasiclassical description of currents, spin and charge accumulation
In Sec. II.C we present the full Usadel equation, Eq. (5) but only focus on the spectral properties of a mesoscopic superconductor, i.e., on the retarded and advanced components of Green's function. For a quantitative description of transport properties, we need also to consider the Keldysh componentĝ K of the quasiclassical GF [Eq. (6)]. By exploiting the normalization condition, Eq. (14),ĝ K can be expressed in terms of the retarded and advanced components and a generalized matrix distribution functionf (Langenberg and Larkin, 1986) In the most general case, the 4×4 matrix distribution functionf can be written aŝ Here we use the labeling introduced by Schmid and Schön (1975), generalized for the spin-dependent case. They pointed out that a change in the L/T function causes a deviation of the superconducting order parameter parallel/perpendicular to the equilibrium value in the complex plane. As shown in Eq. (54) below, our notation of the spin-dependent functions does not follow this scheme. The L-labeled functions describe longitudinal modes, the (spin) energy degrees of freedom, and are antisymmetric in energy with respect to the Fermi level, ε = 0. The T -labeled functions describe transverse modes and are symmetric in energy. In equilibrium, the distribution function is proportional to the unit matrix in Nambu and spin space, and given bŷ Each component off in Eq. (22) can be associated with a nonequilibrium mode as discussed in the next subsection. The charge, spin and heat current densities can be related to the Keldysh components of the matrixǧ∇ǧ by taking different traces in the Nambu ⊗ spin space, as follows (for the k = x, y, z-components) (Belzig et al., 1999).
(i) The charge current density (ii) the spin current density polarized in a-direction (iii) the energy current density and (iv) the spin energy current density Here σ N = 2e 2 Dν F is the normal state conductivity and ν F is the normal-state density of states per spin at the Fermi level. Compactly, we introduce the spectral current tensor J a kb , The particular current density components can be obtained from Eq (28) by multiplying with σ N /(2e 2 ) and integrating with energy; for the energy currents after multiplying with the energy ε.
The structure of the matrix current (28) in Nambu space is determined by the choice of spin-Nambu basis which in this review is the one in Eq. (7). The advantage of such a choice with respect to others used in the literature (see for example (Bergeret et al., 2005)) is that the Nambu structure is the same for all spin components.
Typical experiments on the transport properties of superconductors use multi-terminal devices and non-local measurements, in the sense that charge currents are injected from an electrode whereas the voltage is measured in a different part of the superconductor where no current is flowing (Beckmann et al., 2004;Clarke, 1972;Hübler et al., 2012;Quay et al., 2013;Wolf et al., 2014aWolf et al., , 2013Wolf et al., , 2014bYagi, 2006). Alternatively, the current induced in a detector electrode is measured (Hübler et al., 2010). By using normal electrodes one can determine in such setups the dependence of charge imbalance on various parameters, such as the distance between the injector and the detector, the injection voltage or the magnetic field.
In the same setup, but with ferromagnetic electrodes, one can probe both the charge and spin accumulations, which in terms of the Keldysh component, Eq. (21), of the GF are given by The local energy and spin-energy accumulations are defined similary, Above, a = 1, 2, 3 denotes the polarization direction of the nonequilibrium spins. Spin accumulation is a standard observable in spintronics (Jedema et al., 2002;Johnson and Silsbee, 1985). The local energy accumulation is typically measured via electron thermometry (Giazotto et al., 2006). The spinenergy accumulation was measured recently in normalstate nanopillar spin valves (Dejene et al., 2013). To our knowledge this quantity has not been directly studied experimentally in superconducting systems.
The nonequilibrium modes, i.e., the different parts of the matrix distribution functionf , are obtained from a kinetic equation forf . In a diffusive conductor in the steady-state limit, the kinetic equation can be obtained by inserting Eq. (21) in the Usadel equation (5) and evaluating its different matrix components. This results into a set of diffusion-type equations of the generic form where I ab coll = Trτ b σ a [Σ,ǧ] K /8 is a collision integral describing a given scattering process with self-energyΣ. The term H ab = Trτ b σ a [−ih · στ 3 ,ĝ K ]/8 describes the Hanle precession of spin caused by the exchange field, and R ab = Trτ b σ a [∆,ĝ K ]/8 the conversion between quasiparticles and the superconducting condensate. The detailed forms of the different components of the kinetic equations in spin-split superconductors are discussed in Sec. IV.A. Note that the definition of the collision integrals here depends on the form of the diffusion equation adopted; another definition is discussed in Sec. VI.A.
Equation (33) neglects the charge-spin coupling induced by the spin-orbit interaction that leads to the spin Hall and other magnetoelectric effects (Sinova et al., 2015). In this review we only consider conventional superconductors, i.e., metals with a centrosymmetric crystal structure. In this case, the spin-orbit coupling originates from the random potential of the impurities. After averaging over it, a bulk superconductor remains inversion symmetric and hence the only effect of the spin-orbit interaction is to relax the spin and the triplet component of the condensate (Bergeret and Tokatly, 2016). In contrast, currents flowing through superconducting films of finite dimensions may induce a spin accumulation at the film edges via the spin Hall effect (Bergeret and Tokatly, 2016;Espedal et al., 2017). The spin accumulation is proportional to the spin Hall angle θ, which is a material dependent parameter. In Al films this contribution is negligibly small, and the effect of the SOC is mainly to relax the spin. This is the limit we adopt in Sec. IV, when studying the transport properties on lateral spin valves (see . Neglecting terms of the order of θ in Eq. (33) reduces the problem considerably to a 1D boundary value problem.
In a non-superconducting conductor, Eq. (33) reduces to a diffusive equation of the form where f ab = Tr[τ b σ af ]/4 are the distribution functions characterizing the different nonequilibrium modes. In this case, it is unnecessary to separate the energy-even (transverse, b = 3) and odd (longitudinal, b = 0) modes, and rather consider the spin-dependent full distribution In the normal state, solutions of this kinetic equation are discussed for example in (Brataas et al., 2006).

B. Nonequilibrium modes in superconductors with spin splitting
The nonequilibrium modes can be understood in a framework similar to that often used in the discussion of charge imbalance in the absence of spin dependence [see, for example, (Tinkham, 1996)]. Schematically, nonequilibrium modes can be represented in terms of the electron/hole branches of the spectrum of the superconductor, illustrated in Fig. 9. The spin splitting 2h = 2µ B B of spin states gives rise to four distinct quasiparticle branches, electron/hole and spin up/down. If we consider a unique spin polarization direction parallel to the z-axis, and neglect noncollinear spin components off , there are only four components of the electron distribution function Eq. (22) corresponding to the modes illustrated in Figs. 9(b)-(e). Two of these modes have electron-hole branch imbalance, f T and f L3 , while f T 3 and f L are particle-hole symmetric.
The filled circles in Fig. 9 represent the occupied states. As a reference, panel (a) corresponds to the equilibrium distribution functionf = f 0 L1 = tanh(ε/2T )1. In order to excite the nonequilibrium modes, f T , f T 3 and f L3 , one only needs to move the populated states (filled circles) between the different spectral branches in an elastic process, i.e., between equal-energy states (marked by horizontal dashed arrows). These modes can also relax back to equilibrium due to elastic scattering processes. The relaxation mechanisms depend on intrinsic material properties, degree and type of disorder, and also on the superconducting spectrum, and are discussed in more detail below.
The last nonequilibrium mode, the deviation of f L from f 0 L , is characterized by a change in the total quasiparticle number and energy content, corresponding to an increase or decrease of the effective temperature. It can be excited by increasing the number of occupied states to higher energies, and its relaxation requires inelastic processes.
In the absence of spin splitting, the charge imbalance is determined by f T , and the energy imbalance by f L . The spin splitting changes the system properties, mixing the coupling between spin-dependent modes and physical observables. Qualitatively, the outcome can be seen by counting the number of occupied states on the different branches in Fig. 9. For example, the charge imbalance µ is determined by the difference between the number of occupied states in the electron and hole branches. Both f T and f L3 components contribute to it, as seen in Figs. 9(b) and (d).
On the other hand, a nonzero spin accumulation µ z can be induced by exciting the modes f T 3 or f L [Figs. 9(c),(e)]. These two contributions to the total spin accumulation have important differences: The mode f T 3 contributes to spin imbalance also in the absence of spin splitting. Spin imbalance in this mode can be induced for example by a spin-polarized injection from a ferromagnetic electrode, in both the normal and the superconducting state. The relaxation of the spin accumulation created in this way is determined by elastic scattering processes [ Fig. 9(d)]. The second mechanism of inducing spin accumulation is by exciting the longitudinal mode f L , in the presence of spin splitting [see Fig. 9(e)]. Since energyconserving transitions do not result in the relaxation of the f L mode, this component of the spin imbalance is not suppressed by elastic scattering. In other words, its relaxation can be only provided by inelastic processes, e.g., electron-phonon and electron-electron scattering. This result, obtained here on a phenomenological level, is crucial in understanding the long-range spin signal observed in superconductors, for example by Hübler et al. (2012). A quantitative analysis based on the Usadel equation is provided in Sec. IV.B.

C. Description of hybrid interfaces
The currents given by Eqs. (24-28) are defined in the bulk of a superconductor (S) or a normal metal (N). For the description of transport in hybrid structures, we need in addition a description of interfaces between different materials in the form of boundary conditions. Such interfaces usually are described by sharp changes of the potential and material parameters over atomic distances, and thus cannot be included straightforwardly in the quasiclassical equations which describe properties over distances much larger than λ F . The description of hybrid interfaces requires then the derivation of suitable boundary conditions, first done in the quasiclassical approach by Zaitsev (1984).
Boundary conditions for the Usadel equation trace back to the work of Kupriyanov and Lukichev (1988). These boundary conditions are applicable for nonmagnetic N-N, S-S and S-N interfaces with low transmissivity (Lambert et al., 1997). Later Nazarov general-ized these boundary conditions for an arbitrary interface transparency (Nazarov, 1999). Tokuyasu et al. (1988) derived the boundary condition for an interface between a superconductor and a ferromagnetic insulator and introduced the concept of the spin-mixing angle, which describes the spin-dependent phase shifts acquired by the electrons after being scattered at the FI/S interface. Later Cottet et al. (2009) and Zhao et al. (2004) extended these boundary conditions to magnetic metallic structures, such as F-S or S-F-S systems, though with low polarization. Boundary conditions for large polarization and low transmission have been presented by Bergeret et al. (2012a) and Machon et al. (2013). General boundary conditions for arbitrary spin polarization and transmission have been extensively discussed by Eschrig et al. (2015).
Here we mainly deal with low transmissive barriers between a mesoscopic superconductor and normal and magnetic leads and use the description presented by Bergeret et al. (2012a). When the interface is polarized in the zdirection, the transmission through it can be described by the tunneling matrixΓ = tτ 3 + uσ 3 . It is defined through the normalized transparencies satisfying t 2 + u 2 = 1 and determining the interface polarization via 2ut = P . The possible non-collinearity between the interface and the bulk magnetizations can be described by introducing rotation matrices in spin space. We assume that this inter-face connects a normal-metal electrode with bulk Green's functionǧ N with a spin-split superconductor with GFǧ. The boundary condition gives an expression for the matrix current through the interface in terms ofΓ,ǧ N anď g:ǧ where R characterizes the spin-averaged barrier resistance per unit area. Its precise microscopic definition is given in (Bergeret et al., 2012a;Eschrig et al., 2015). In experiments it is usually used as a fitting parameter. Equation (35) is valid for an arbitrary spin polarization P ∈ [−1, 1] of the interface. From this condition and Eqs. (24-27) we obtain the spectral currents through spin-polarized barriers as The normal electrode Green's functionǧ N is defined bŷ g is the distribution function of the normal metal. The boundary conditions for both spectral GF and distribution functions are given by the continuity of the matrix current n k J a kb = I ab , where the tunneling current I ab and bulk current J a kb are provided by Eqs. (36) and (28), respectively, and n k are the components of the unit vector normal to the interface.
Low effective transmissivity of interfaces which corresponds to a large interface resistance (tunneling limit), allows for simplified approximations to be made in the description of nonequilibrium effects. In such a case, changes in the spectrum due to the coupling can be neglected, and the probability for Andreev reflection (Andreev, 1964) vanishes. In the lowest order in transmissivity, we can also assume the electrodes to be in local thermal equilibrium. If, in addition, P = 0, the charge current between two electrodes is given by the usual tunneling expression Eq. (1) with G T = A/R , where A is the cross-sectional area of the contact. Below we describe the more general limit of an arbitrary P , and the different components of the spectral current.
It is worth emphasizing that the boundary condition, Eq. (35), when combined with the quasiclassical equations, allows for the description of effects which are beyond the quasiclassical limit and that involve strong ferromagnets with large spin polarization (Bergeret et al., 2012b;Silaev et al., 2017).

D. Linear response and generalized Onsager relations
A way of exciting the nonequilibrium modes in superconductors is by injecting electrons from either a voltage or temperature biased metallic lead. In this case and in the language of Eq. (22), the distribution function in the normal metal is given byf L/T = [f eq (ε−eV )±f eq (ε+eV )]/2 is the voltage-biased distribution function in the normal-metal electrode and f eq (ε) = tanh(ε/2T N ) is the equilibrium distribution corresponding to the normal-metal temperature T N .
In a more general situation, all components of the distribution function, Eq. (22), have to be taken into account. The different components of the spectral current through a tunneling contact between a normal electrode and a superconductor can be derived from the boundary condition in Eq. (36). In the case of collinear magnetizations the relevant components are the spectral charge j c ≡ I 30 , energy j e ≡ I 00 , spin j s ≡ I 03 and spin energy j se ≡ I 33 currents. They satisfy where the parameter κ = 1/(R e 2 ) describes the inter- , between the superconductor and normal-metal electrodes. The response matrix is here described by the spin polarization P and the energysymmetric and energy-antisymmetric parts of the density of states, N ± defined in Eq. (15) and below it.
Expression (37) is valid for arbitrarily large deviations from the equilibrium state and describes several intriguing effects which come into play in the nonequilibrium situation in junctions combining spin-split superconductors (N − = 0) and a spin-polarized injector, modeled by a nonzero spin polarization P = 0. If both N − = 0 and P = 0 simultaneously, all the nonequilibrium modes are coupled. In turn, this coupling results in a rich variety of cross-couplings between different potentials and currents.
To illustrate the basic phenomena which can be expected in such a tunneling contact, let us consider the linear response limit. The different components of the distribution function can be described by spindependent temperatures 6 and voltages introduced according to Bauer et al. (2012) where V and ∆T are the voltage and the temperature bias and V s and ∆T s the spin-dependent biases. Integrating Eq. (37) over energy allows us to calculate the total charge I = e ∞ −∞ j c dε, energyQ = ∞ −∞ εj e dε, spin I s = ∞ −∞ j s dε, and spin energyQ s = ∞ −∞ εj e dε currents. They are related to the generalized potentials via a 4 × 4 Onsager matrix 7 where the conductance, heat conductance and thermoelectric coefficient are given by respectively. 8 Equation (39) is consistent with the time-reversal transformation properties expected for a linear-response matrix (Jacquod et al., 2012;Machon et al., 2013;Onsager, 1931). Note that V, −∆T /T, I andQ are all odd in the time-reversal transformation, whereas the corresponding spin-resolved quantities V s , −∆T s /T, I s andQ s are even. Therefore, the four 2 × 2 quadrants of the response matrix have different time-reversal symmetries: the quadrants on the diagonal are symmetric, whereas the off-diagonal quadrants coupling spin and charge are antisymmetric. This result follows from the fact that G and G th are symmetric whereas α and P are antisymmetric in time reversal.
From Eq. (39) we can, for example, study the spin and charge currents created at the SF interface in response to a temperature bias between the normal lead and the superconductor. These are the thermospin and thermoelectric effects depicted schematically in Fig. 1(c-d). Interestingly, Eq. (39) demonstrates also that the charge current can be induced by the spin-dependent temperature bias ∆T s even without spin-filtering (P = 0) but in the presence of a spin-split DOS, N − = 0. Qualitatively the corresponding non-equilibrium mode f L3 [ Fig. 9 can be interpreted as the spin-dependent particle-hole imbalance which produces electric signal due to N − = 0.
Let us focus on the case where the nonequilibrium state is generated by applying a voltage and a temperature bias to the injecting normal electrode in the absence of spin-dependent potentials, ∆T s = 0 and V s = 0. If the interface lacks spin polarization, P = 0, then from Eq. (39) it follows that all currents are decoupled from each other, recovering standard expressions for charge and heat currents in terms of the local electrical and thermal conductances I = GV andQ = −G th ∆T /T . In addition, the remaining block of the Onsager matrix leads to nonzero spin current I s = −α∆T /T and spin heat currentQ s = αV . From this linear response analysis we can conclude that in a setup consisting of a normal and a spin-split superconducting electrode, a spin current can be generated even in the absence of ferromagnetic electrodes. It is worth mentioning that this thermospin effect leads to a pure spin current I s that does not carry electric charge and therefore does not produce Joule heating and Ohmic losses. Clearly this situation is very different from the one occurring in normal metal systems, where spin currents can be generated only by injecting a charge current from a ferromagnet -a process corresponding to the relation I s = P GV obtained from Eq. (39) and schematized for a spin-polarized contact to a superconductor in Fig. 1(a).
If P = 0 all four types of currents are injected into the superconducting region just by applying the temperature bias ∆T . In particular a very large thermoelectric effect can be generated, described by I = −P α∆T /T . This effect is discussed in detail in Sec. V.

E. Spin relaxation mechanisms
In metals and semiconductors there are mainly two types of mechanisms that relax the spin: scattering from magnetic impurities and spin-orbit coupling. As shown in Sec. II.D, these processes can modify substantially the spectral properties of a superconductor. Here we consider the additional influence on the relaxation of the nonequilibrium modes.
In the case of spin-orbit coupling we can distinguish two types of relaxation mechanisms according to the origin, intrinsic or extrinsic. Intrinsic spin-orbit coupling occurs in systems with the lack of inversion symmetry, generated either by geometry constraints or by the crystal potential associated with the electronic band structure. Momentum-dependent spin precession together with random momentum relaxation due to impurities leads to the Dyakonov-Perel relaxation (Dyakonov and Perel, 1972). For example, for a diffusive 2D system with a Rashba spin-orbit coupling α R , the Dyakonov-Perel relaxation time equals to τ = 1/D(2mα R ) 2 for a in-plane spin-polarization, and τ ⊥ = τ /2 for spins perpendicular FIG. 10 Spin relaxation rates for (a) h = 0 and h = 0.3∆0 and for spin-orbit (β = −1) and spin-flip (β = 1) relaxation mechanisms. We assume that the normal-state spin relaxation rate is τ −1 sn = 0.2Tc0. to the 2D system. In contrast, the extrinsic spin-orbit effect originates from a random impurity potential and hence leads to an isotropic spin relaxation known as the Elliott-Yafet relaxation mechanism. It is described by the self-energy Eq. (10). A detailed discussion of spin relaxation mechanisms originated from the spin-orbit coupling can be found in the review byŽutić et al. (2004) and references therein.
In this review we focus on centrosymmetric materials, and therefore we only consider the extrinsic relaxation mechanism. Moreover, we do not consider effects related to the spin-charge coupling, such as the spin Hall effect (Sinova et al., 2015) and the spin galvanic effects (Edelstein, 1990). Such effects enter the kinetic equations in a higher order of the gradient expansion and can be neglected in the leading order. The spin-charge coupling in superconductors leads to non-dissipative magnetoelectric effects and the appearance of an anomalous phase in Josephson junctions, as discussed by Konschelle et al. (2015) and Bergeret and Tokatly (2016).
In the case of the extrinsic relaxation, the collision integrals entering the kinetic equation due to spin-dependent scattering are given by where the components a, b refer to the different nonequilibrium modes and spin projections. The self-energies are defined in (10) and (11). In the superconducting state, the form of the spin-flip and spin-orbit contributions to the collision integrals differ due to their different properties under time reversal.
Let us now discuss the relaxation of the nonequilibrium modes related with the spin and the spin energy imbalances, i.e., the components f Specifically, the spin relaxation, shown schematically by the dotted lines in Fig. 9(c), is described by the colli- In the normal state these processes are energy independent (within the quasiclassical approximation), and have the rate τ −1 sn . In the superconducting state, the spinrelaxation is described by Eq. (44), with the energydependent rate τ −1 s = S T 3 (ε)/N + . Here the collision integral is normalized with the spin-averaged density of states N + to counterbalance effects due to the mere changes of the quasiparticle spectrum in the superconducting state. The typical behavior of τ −1 s (ε) is shown in Fig. 10 for dominating spin-flip and spin-orbit mechanisms for (a) zero and (b) finite spin splitting.
If the spin-orbit mechanism dominates, β = −1, superconductivity suppresses spin relaxation, which can be qualitatively understood as the decrease in cross section of the quasiparticle momentum scattering at the energies near the gap edge ε ∼ ∆. The suppression of τ −1 s is however compensated by the decrease in the quasiparticle group velocity v g ∼ v F 1 − |∆| 2 /ε 2 so that the spin relaxation length λ so ∼ v g τ s remains almost unchanged in the superconducting state (Morten et al., 2004). On the contrary, if the spin-flip mechanism dominates, β = 1, the spin relaxation is not related to the momentum scattering because the interaction with magnetic impurities does not depend on the propagation direction and the quasiparticle spin does not depend on energy. This results in an increase of τ −1 s which is equivalent to a decrease of the spin-relaxation length in the superconducting state (Morten et al., 2004).
In the presence of a Zeeman or exchange field, h = 0, this situation changes drastically. In this case the nonequilibrium spin accumulation couples to the energy mode, f L . This mode is robust with respect to elastic spin relaxation processes, resulting to a qualitative change in the observed spin relaxation. This effect, as well as the detailed energy dependence of the spin relaxation lengths, are discussed in detail in Sec. IV.B.
For the relaxation of the spin-energy mode f L3 , the situation is more complicated. As sketched in Fig. 9(d), the relaxation of this mode involves transitions between electron and hole branches of the quasiparticle spectrum. As in the case of charge imbalance relaxation, such interbranch transitions involve the formation of Cooper pairs via Andreev reflection. This results in two contributions to the relaxation of the spin-energy mode: One stems from the (j = 3, k = 3) component of the collision integral (43) that describes transitions between spin-up and spin-down subbands S s The second contribution originates from the particle-hole relaxation described by the non-diagonal self-energy in the Keldysh-Usadel equation (5), The relaxation of the spin-energy mode in the superconductor is given by the sum R T + S L3 . Note that the particle-hole conversion processes responsible for the spin energy relaxation also lead to the generation of charge imbalance producing a coupling between spin and charge degrees of freedom. These mechanisms are discussed in more detail in Sec. IV.B. At large energies ε ∆ the charge and spin energy nonequilibrium modes are decoupled. In this limit the charge relaxation length tends to infinity because the coherence factor Reg 01 → 0. At the same time, the spin-energy relaxation length approaches its normal-state value equal to that of the spin relaxation length Λ sn = √ Dτ sn .

F. Electron-phonon heat flow in a superconductor with spin splitting
One of the main features limiting some of the effects discussed in the following sections arises from the coupling between the quasiparticles and phonons, residing at different temperatures T qp and T ph , respectively. Due to the strong energy dependence of the phonon density of states, the effect of this coupling decreases rapidly towards low temperatures, and eventually phonons decouple from electrons, and the main heat relaxation occurs via other processes such as quasiparticle diffusion. The presence of superconductivity further modifies the electron-phonon heat conduction (Eliashberg, 1972;Kaplan et al., 1976;Kopnin, 2001) as also the electronic spectrum is energy dependent. Here we discuss the main features of the electron-phonon collisions in spin-split superconductors.
The collision integral for electron-phonon process in a spin-split superconductor was discussed in Grimaldi and Fulde (1997) and can be obtained from the self-energies described in Eq. (13) . The overall collision integral entering the kinetic equations [Keldysh part of Eq. (5)] is obtained from a commutator between the self-energy and Green's function,Ǐ eph = [Σ eph ,ǧ] K . This is still a Nambu-spin matrix, its different matrix components describe the relaxation of the different nonequilibrium modes. Electron-phonon interaction causes inelastic scattering of quasiparticles, which means that energy flow between quasiparticles and phonons becomes possible. That process is described by the collision integral I L ≡ Tr[Î eph ]/8. Also the other matrix components of the collision integral are in general nonvanishing; for example the component Tr[τ 3Îeph ] affecting the energy-even distribution function is detailed in (Heikkilä and Giazotto, 2009). However, since electronphonon interaction conserves charge and spin, the energy integrals of Tr[τ 3Îeph ] and Tr[σ 3Îeph ] vanish (see Appendix B). The remaining component Tr[τ 3 σ 3Îeph ] describes the relaxation of the difference in the thermal energies of electrons with opposite spins due to electronphonon interaction (Heikkilä et al., 2010). In the following, we only concentrate on energy relaxation. Substituting the self-energy (13) and assuming thatf = f L1 we can write Typically this collision integral is used to calculate the total heat currentQ eph = 2ν F Ω dεεI L (ε) between the electron and phonon systems in a volume Ω in a metal with the Fermi-level density of states ν F per spin. Assuming that the electron system resides in temperature T qp , we have f L = tanh[ε/(2k B T qp )] and thus the heat current is (Maisi et al., 2013;Timofeev et al., 2009), Here Σ = 3072ζ(5)ν F k 5 B g eph describes the materials dependent magnitude of the electron-phonon coupling [see tabulated values in (Giazotto et al., 2006)] typically used for heat flow, and ζ(x) is the Riemann zeta function. The kernel L E,E depends on the spin-splitting field, as it is where N σ (ε) = Re(g 03 + σg 33 )/2 is the density of states for electrons with spin σ and F σ = Re(g 01 +σg 31 )/2 is the anomalous function for spin σ (see Sec. II.D). The latter k B T =" formula applies in the absence of spin-flip or spin-orbit scattering.
Strictly speaking, Eq. (49) requires two conditions to be met: (i) A well-defined local quasiparticle temperature that may deviate from the phonon temperature. In other words, the quasiparticles should first equilibrate between themselves. (ii) In the presence of a spin-splitting field, also spin imbalance affects the quasiparticle-phonon heat flow . Neglecting it therefore requires that the spin relaxation rate exceeds the injection rate causing the out-of-equilibrium situation. 9 In the general case Eq. (49) needs to be evaluated numerically. However, we may consider two limiting cases. In a normal metal ∆ = 0 implies L ε,ε = 1. In that case the heat current between the two systems becomes (Wellstood et al., 1994) Another tractable limit is that of linear response, where 9 Note that strictly speaking one should then include this spin relaxation into the kernel L E,E as spin relaxation affects the spectral functions, see Fig. 5. The relevant scales for the two effects differ, however, as the kernel is modified when /τsn is not very much weaker than ∆. This is further quantified in Sec. V.F. where∆ = ∆/k B T ,h = h/k B T , and f 1 (x) = 3 n=0 C n x n , f 2 (x) = 2 n=0 B n x n , C 0 ≈ 440, C 1 ≈ 500, C 2 ≈ 1400, C 3 ≈ 4700, B 0 = 64, B 1 = 144, B 2 = 258. The two terms in Eq. (52) describe scattering and recombination processes, respectively. Equation (52) is compared to the exact result in Fig. 11. The recombination process, whose heat conductance is almost independent of the spin splitting (as long as h ∆), dominates at high temperatures, whereas the two become of the same order of magnitude for T ≈ 0.1∆. Note that the exponential suppression of the electron-phonon heat conductance ∼ e −∆/T at low temperatures does not directly make it irrelevant in a superconductor, because any effect related to quasiparticles contains such exponential terms due to the gap in the superconducting density of states. However, as lowering the temperature reduces the phase space available for acoustic phonons, also the prefactor of the exponential becomes small low temperatures. Thus, at very low temperatures quasiparticles decouple from phonons, and other heat conduction mechanisms become relevant.

G. Particle-particle collisions
The collision integralÎ ee related to the particle-particle collisions in superconductors and superfluids is discussed in (Eliashberg, 1972;Kopnin, 2001;Serene and Rainer, 1983). These concentrate on a contact interaction model and disregard self-consistent screening effects in dirty superconductors (Feigel'man et al., 2000;Kamenev and Levchenko, 2009;Narozhny et al., 1999). The spin structure of collision integrals in the normal state with spin accumulation is discussed in (Chtchelkatchev and Burmistrov, 2008;Dimitrova and Kravtsov, 2008). Below, we only point out conservation laws that a consistent formulation of the collision integrals must imply in a system involving collinear magnetizations. These are the conservation of charge, spin, and overall energy by spinrotation symmetric scattering within the particle system: The spin energy, involving the component ετ 3 σ 3 of the collision integral, is generally not conserved in electronelectron collisions, because electrons with opposite spins can exchange energy (Heikkilä et al., 2010). The above conservation laws are relevant in terms of studying what happens in the limit of strong electronelectron scattering. In that case, the usually assumed model is that the quasiparticle distribution functions tend to equilibrium forms, but retaining some effective temperature and a spin-dependent chemical potential. This quasiequilibrium limit (without spin accumulation) is often used in analyzing heat transport in superconductors driven out of equilibrium (for example, see (Muhonen et al., 2012)).
The far-from-equilibrium results discussed in Sec. IV disregard the particle-particle collisions as their full inclusion would complicate the theory quite much and the simpler theory already describes effects not very far from the measured ones. On the other hand, Sec. V mostly concentrates on the quasiequilibrium limit, where also spin accumulation is lost due to a strong spin relaxation.

H. Nonequilibrium effects on the superconducting order parameter
Inserting the parametrization (21,22) of the nonequilibrium Green's function in Eq. (9), the gap equation becomes (54) In equilibrium, only the first term contributes with f L = f eq = tanh[ε/(2T )]. The imaginary terms affect the phase of ∆ and ensure charge current conservation in situations where quasiparticle current is converted to supercurrent (described by the R ab term in Eq. (33)). In the absence of spin splitting, the nonequilibrium modifications in the size of the gap, described by the first term, are the most relevant. Often such modifications are studied within the phenomenological model of Rothwarf and Taylor (1967) that neglects the energy dependence of the distribution functions, and rather concentrates on the overall number of quasiparticles. This model was derived under some simplifiying assumptions from the full energy dependent kinetic equations in (Chang and Scalapino, 1977). This model has been used to study the gap suppression due to nonequilibrium injection, see for example (Yeh and Langenberg, 1978). According to them, for small changes around equilibrium, such models agree with the predictions of the phenomenological models by (Owen and Scalapino, 1972;Parker, 1975). For large changes, the full kinetic equations with energy dependent distribution functions should be used. An extreme example is given in (Keizer et al., 2006a) describing a voltage-driven transition of a superconducting wire to the normal state.
In the spin-dependent case we need to consider both the effects of spin splitting and the spin-dependent distribution functions describing spin accumulation. Disregarding the terms affecting mostly only the phase of the order parameter, we can also write Eq. (54) as where F σ = (g R 01 + σg R 31 )/2 and f σ = (f L + σf T 3 )/2; the different signs of σ representing the two spin directions. Takahashi et al. (1999) went a step further, assuming that the spin accumulation in a superconductor can be described in terms of a simple spin-dependent chemical potential shift µ s , neglecting all other nonequilibrium effects. In that case f σ = f eq (ε − σµ s ). As a result, by a simple shift of the energy (ε → ε + σµ s ) in Eq. (55), µ s shows up as a shift of the energy of the anomalous function F σ and a small shift of the cutoff Ω D . When Ω D µ s , the latter effect can be disregarded, and the net effect of the spin accumulation is the same as that of the spin splitting field, explained in Sec. II.D, and eventually leading to the suppression of superconductivity. Moreover, in the presence of both spin splitting and spin accumulation, this model shows how in the special case h = µ s spin accumulation can actually lead to the recovery of superconductivity suppressed by spin splitting (Bobkova and Bobkov, 2011).
This physics can be probed in a FISIF system where a superconducting island or layer is placed between two ferromagnetic electrodes. In this case, the current induced suppression of the superconducting gap should be larger when the magnetizations of the two ferromagnets are antiparallel than when they are parallel (Takahashi et al., 1999). Only in the previous case the spin accumulation builds up. This effect was measured in (Yang et al., 2010). They indeed found a stronger suppression of superconductivity in the antiparallel configuration. However, they found that the required spin relaxation time for fitting the results to the above theory is much longer than that expected from the normal-state measurements. This discrepancy may result from the somewhat simplified model for the spin accumulation described above.
The suppression of superconducting properties due to spin injection has been measured in the case of hightemperature superconductors (Gim et al., 2001;Koller et al., 1998;Vas'ko et al., 1997;Yeh et al., 1999). However, these features are typically attributed to the non-conventional character of superconductivity, and are therefore outside the scope of this review.
Besides the superconducting gap, in principle also the induced spin-splitting field can obtain nonequilibrium corrections in the presence of spin injection. The theory for such effects was outlined already by (Alexander et al., 1985), but to our knowledge such effects have not been thoroughly examined in spin-split superconductors.
I. Overall strategy to explore the transport and spectral properties of hybrid superconducting systems Spectral and transport properties of a diffusive metal or superconductor of mesoscopic size attached to electrodes can be fully described by solving the boundary problem defined by the Keldysh-Usadel equation (5), the normalization condition forǧ, Eq. (14), and the boundary conditions at the interfaces with the electrodes, Eq. (35). Once the Green's functions are determined one can compute the currents and potentials from Eqs. (24-32).
If the system under consideration is in the normal state, the equations for the retarded and advanced GFs are decoupled from each other and from the Keldysh one. In such a case one can solve first the spectral part and then plug the retarded and advanced GFs in the equation for the Keldysh component.
If the system is a superconductor attached to leads, the equations for the retarded and advanced GFs are coupled to the Keldysh component through the selfconsistency equation for the superconducting order parameter, Eq. (9). This complication can be overcome if the superconductor and the electrodes are coupled via tunneling junctions. In such a case the spectral properties of the superconductor remain unchanged within the leading order in the interface transmission. This means that retarded and advanced GFs coincide with those of the homogeneous superconductor. In Sec. IV, we concentrate on the case where the normal-state tunnel junction conductance satisfies G T R E 1, where R E is the normal-state resistance of a wire with length equal to the energy relaxation length (due to electron-phonon or electron-electron scattering) E at the energy around the amplitude of the superconducting gap. On the other hand, in Sec. V, we mostly consider the cooling of the superconducting island with a starting temperature much below T c . In these cases we can disregard the nonequilibrium effects on the energy gap, and can rather use the self-consistent gap that is calculated from the equilibrium version of self-consistency equation, i.e., Eq. (54), with In most of the discussed examples of this review we follow this approach and use the self-consistent ∆ calculated in the decoupled spin-split superconductor. From this we obtain the spectral functions of the spin-split superconductor that yield the coefficients of the kinetic equation for the nonequilibrium distribution functions. The latter is obtained from the Keldysh component of Eq. (5). Only in Sec. VI.B, where we study the effect of an ac field on the spectrum of a spin-split superconductor, we compute the self-consistent gap using the nonequilibrium distribution function (see Fig. 36).
Inspection of the different components in Keldysh space of Eq. (5) provide a first insight about the characteristic lengths involved in the different situations. In a superconductor changes of the spectral properties due, for example, to the magnetic proximity effect with a ferromagnetic insulator, occur over the superconducting coherence length ξ s . On the other hand, this ξ s depends on the temperature, and on the concentration of magnetic and spin-orbit impurities via the self-consistent ∆.
Different length scales govern the decay of nonequilibrium components of the distribution function generated, for example, at the contact with an electrode by injection of a current. In this case the characteristic lengths depend on the nature of the excited mode (charge, energy or spin) and the type of relaxation process in the system. The calculation of these characteristic lengths is one of the main goals of Sec. IV. As we show there, the spin splitting plays a crucial role in determining these length scales, since it couples the different nonequilibrium modes, and hence changes the range over which non-local spin and charge signals can be detected.

IV. SPIN INJECTION AND DIFFUSION IN SUPERCONDUCTORS
Electrical injection of spins into a superconductor was first studied by Aronov (1976). The injection of an electric current from a F electrode into a superconductor creates not only spin, but also charge imbalance. Moreover, such imbalances relax away from the injection point over distances different from those in the normal state.
In the absence of a spin-splitting field in the superconductor the spin imbalance is decoupled from the other modes. However, in spin-split superconductors charge, spin and energy modes couple with each other. The goal of this section is to describe this coupling, the characteristic relaxation lengths of the different nonequilibrium modes and effects that occur in superconducting structures as a consequence of such coupling.
We start by reviewing experiments on charge and spin injection in superconductors and by summarizing the main theoretical works on spin injection in superconductors with no spin splitting. Further on, we take into account the spin splitting induced by external fields, and describe different experimental situations with the help of the kinetic equations derived in Sec. III.

A. Detection of spin and charge imbalance: Non-local transport measurements
Studies of the nonequilibrium modes started with the pioneering experiment of Clarke (1972), who realized a way of detecting the charge imbalance in a superconductor. The geometry used in that experiment was a vertical structure as the one sketched in Fig. 12(a). Injection of a current from a normal metal (injector) into the superconductor (blue stripe) leads to a charge imbalance that in turns leads to a shift of the chemical potential of the quasiparticles with respect to the one of the condensate. This shift of the chemical potential can be detected by a second electrode [detector in Fig. 12(a)], that probes the voltage between the superconductor and the detector. First experiments demonstrating the charge imbalance induced in superconductors by current injection were performed at temperatures close to the critical temperature T c by using vertical samples with large cross sections, as the one depicted in Fig. 12(a). More recent experiments have allowed for non-local transport measurements in superconducting lateral structures and lower temperatures, as the one sketched in Fig. 12(b) (Beckmann et al., 2004;Hübler et al., 2010;Poli et al., 2008;Quay et al., 2013;Wolf et al., 2014aWolf et al., , 2013. A scanning electron microscopy image of a lateral structure is shown in Fig. 13. Such structures allow for an accurate measurement of the charge imbalance and nonequilibrium energy mode and their spatial dependence by placing contacts (detectors) at different distances from the injector. For example, in an aluminum wire at 100 mK the charge and energy modes decay over 5 and 10 µm, respectively, as reported by Arutyunov et al. (2011). A detailed overview of the experiments on charge and energy imbalance can be found in the recent topical reviews by Beckmann (2016) and Quay and Aprili (2017). Substituting the normal injector by a ferromagnet, the injected current becomes spin polarized and the charge injection is accompanied by a spin injection into the superconductor. The nonequilibrium charge and spin relax over different lengths. When the blue wire in Fig. 13 is in the normal state, charge accumulation is negligible at the detector and only a spin imbalance contributes to the nonlocal signal. The nonequilibrium spin density induced at the interface is polarized in the direction of the injector's magnetization and it diffuses into the normal wire over the spin-relaxation length, which can be several hundreds of nanometers. The spin accumulation can be detected by measuring the non-local voltage between the (ferromagnetic) detector and the N wire. First measurement of the nonequilibrium spin accumulation was reported by Johnson and Silsbee (1985) on a single crystal aluminum bar at temperatures below 77 K. This pioneering experiment did not only demonstrate the spin accumulation associated with the injection from a ferromagnetic electrode, but also the coherent spin precession, the Hanle effect, that occurs when an external magnetic field non-collinear with the injector magnetization is applied. Later these effects have been reported at room temperature in Al-based nanostructures by Jedema et al. (2001Jedema et al. ( , 2003. Since then, electrical injection of spins have been used in several experiments on metallic spintronics devices, for example for the direct detection of the spin Hall effect in metallic structures (Valenzuela and Tinkham, 2006), electronic spin transport in graphene (Tombros et al., 2007), and modulation of the spin density in metal/ferromagnetic insulator bilayers (Dejene et al., 2015;Villamor et al., 2015).
In the superconducting state the situation drastically changes, as reported in several experiments on lateral superconducting structures with ferromagnetic electrodes. In contrast to the the charge imbalance, that has been well understood since the 1980s, spin injection and accumulation in superconductors is a more recent research line. In 1990 Kivelson and Rokhsar (1990) showed theoretically that charge and spin should exhibit different relaxation times in superconductors, leading to the possibility of separating charge and spin transport. Also Zhao and Hershfield (1995) studied the spin injection from a ferromagnet to a bulk superconductor and found theoretically that, while the charge imbalance survives only within the field penetration length from the surface, the spin imbalance may also exist in the bulk. Experiments at that time on F-S-F layered structures, however, did not show any evidence of such a spin-charge separation (Johnson, 1994) and the different relaxation times for spin and charge accumulation in superconductors remained an open question.
First measurement of the spin relaxation length was reported by Gu et al. (2002) in a multilayer vertical F-S-F structure. The experiment showed a decrease of the spin-relaxation length in the superconducting state. Later, clearer insight into the spin and charge modes has been obtained in experiments using lateral nanostructures with ferromagnetic injectors and detectors (Beckmann et al., 2004;Cadden-Zimansky et al., 2007;Hübler et al., 2012;Kolenda et al., 2016b;Poli et al., 2008;Quay et al., 2013;Shin et al., 2005;Wolf et al., 2013Wolf et al., , 2014bYang et al., 2010). By switching the magnetization of the injector and the detector between parallel and antiparallel configurations it has been possible to prove the charge-spin separation.
In order to describe theoretically the different modes excited in experiments on lateral S/F structures, it is convenient to use the generalized quasiclassical model introduced in Sec. III. First theoretical works on spin injection into mesoscopic superconductors (Morten et al., 2004(Morten et al., , 2005 addressed the question of how spin relaxation changes in the superconducting state and how these changes depend on the source of spin relaxation. These works pointed out that in the presence of spin-polarized injected currents additional spin-resolved components of the distribution function may appear. The main conclusion by Morten et al. (2004Morten et al. ( , 2005 was that the spinrelaxation length changes drastically in the superconducting state and strongly depends on the energy of the injected quasiparticles. In particular, for electrons with energy close to the superconducting gap the spin-flip rate may increase. Although these works provided an explanation to some experiments, two important features observed in subsequent works could not be explained in terms of that theory. One intriguing observation was that the spin accumulation can be detected at distances from the injector much larger than the spin-relaxation length measured in the normal state (Hübler et al., 2012;Quay et al., 2013;Wolf et al., 2013). In addition, an unexpected spin accumulation was observed even if the current was injected from a non-magnetic electrode (Wolf et al., 2013).
As we discuss in this section, in order to explain these two observations we need to take into account the spin splitting of the superconductor caused by the applied magnetic field in those experiments. Recent works showed that such a modification of the spectrum of the superconductor leads to an intriguing coupling between the nonequilibrium modes in a superconductor Bobkov, 2015, 2016;Krishtop et al., 2015;Silaev et al., 2015a,b;Virtanen et al., 2016). How this coupling between the modes affects the transport properties of a multi-terminal superconducting structure can be theoretically explored with the quasiclassical formalism developed in Sec. III. Here we review these theory works and provide a quantitative explanation of the long-range spin accumulation detected in multi-terminal superconducting devices.
FIG. 14 Schematic view of the setup for nonlocal transport measurements. (a) Collinear configuration where the polarizations of the magnetic contacts are collinear with the magnetic field, Pinj P det B. The non-local conductance g nl = dI det /dVinj is measured. (b) General non-collinear case with the magnetic field having finite transverse component with respect to the polarizations Pinj, P det . The nonlocal voltage V det is measured in the absence of the current in the detector electrode, I det = 0. The latter detection scheme enhances the spin signal in the superconducting state, see discussion in the text.

B. Nonequilibrium properties of a superconductor with spin splitting
Spin accumulation and transport in superconductors has been considered in a number of theory papers Bobkov, 2015, 2016;Chevallier et al., 2014;Krishtop et al., 2015;Morten et al., 2004Morten et al., , 2005Silaev et al., 2015a,b;Takahashi and Maekawa, 2003;Virtanen et al., 2016). Here we follow the approach developed in (Silaev et al., 2015a) that enables the description of all non-equilibrium modes, the effect of spin-splitting field and the various relaxation mechanisms. As discussed in Sec. III, nonequilibrium conditions lead to the excitation of charge, energy and spin modes in superconducting systems. These modes are related to the different terms of the distribution function introduced in Eq. (22), A prototypical experimental geometry to reveal the different nonequilibrium modes in superconductors is the lateral structure shown in Fig. 14. It consists of a mesoscopic superconductor attached to four electrodes. One of them serves as an injector for the current, whereas another one is used as a detector. We assume first that both electrodes, injector and detector, are ferromagnetic and that their magnetizations are collinear P inj B P det .

28
The tunnelling current at the detector is then given by where R det = R /A is the detector interface resistance in the normal state, A is the cross-sectional area of the detector, µ is the charge imbalance and µ z the spin imbalance defined in Eqs. (29-30). If the current through the detector is measured at zero bias V det = 0, then the nonlocal differential conductance is defined as This quantity provides information about the modes that being excited at the injector via the injected current, reach the detector by diffusing in the superconducting wire. The nonlocal conductance, g nl , depends on the distance L det between the injector and the detector. Such a dependence reveals the characteristic length scale over which the modes relax.
To quantify this dependence, we write the quasiparticle electrostatic and spin potentials µ and µ z , Eqs. (29-30), in terms of the different components of the distribution functions where N + = N ↑ +N ↓ is the total density of states (DOS), N − = N ↑ − N ↓ is the DOS difference between the spin subbands, and f eq (ε) = tanh(ε/2T ) is the equilibrium distribution function. We consider collinear magnetizations along the z-axis and therefore only the spin components f T 3 and f L3 enter these expressions. A noncollinear case is discussed in the next subsection. According to Eqs. (59,60) the full description of the non-local current requires all four non-equilibrium modes. When the distance between the injector and detector electrodes is much larger than the charge imbalance relaxation length the contribution of µ to the detected current can be neglected (Bobkova and Bobkov, 2015;Krishtop et al., 2015). In the absence of spin splitting, the spin-dependent part of I det is determined only by the f T 3 mode related to the spin-dependent chemical potential shift (Morten et al., 2004(Morten et al., , 2005Takahashi and Maekawa, 2003) (see Eq. (38)).
Here the quantities µ and µ z characterize arbitrary nonequilibrium states of the superconducting wire, including those that cannot be described by the effective potentials V , ∆T , V s and ∆T s introduced in Sec. III.D. However, if we restrict the analysis to linear response then Eq. (57) reduces to the expression for the current obtained from Eq. (39) with P = P det , κ = A/(e 2 R det ), µ = R det (α∆T s /2T − GV ) and µ z = R det (α∆T /T − GV s /2).
In the case without spin splitting, h = 0 (N − = 0), The spin accumulation generated by a spin-dependent chemical potential shift δµ = eV s /2 is exponentially small. This can be seen from the first term on the r.h.s of Eq. (60), by noticing that if the energy gap ∆ is nonzero then N + (|ε| < ∆) = 0. In this case the nonequilibrium spin state is described by f T 3 = (∂f eq /∂ε)eV s /2 and hence the spin accumulation is given by µ z = −(eV s /2) ∞ ∆ dεN + (∂f eq /∂ε) so that µ z ∝ e −∆/T at T ∆ . More interesting is the contribution from the second term in the r.h.s of Eq. (60). It is nonzero when the spin splitting is finite and it provides a qualitative explanation of the experiments in (Hübler et al., 2012;Quay et al., 2013;Wolf et al., 2013): The spin imbalance µ z , being related to the energy nonequilibrium mode f L , once excited can only relax via inelastic processes, especially mediated by the electron-phonon interaction. At low temperatures the corresponding decay length can be much larger than the spin decay length in normal metals. This explains the long-range non-local signal observed in the experiments. It can be determined by substituting Eq. (60) into Eqs. (57,58). The observed long-range spin accumulation can be understood to result from the spin accumulation generated by the effective heating of the superconducting wire caused by the injection of nonequilibrium quasiparticles with energies larger than the superconducting gap Bobkov, 2015, 2016;Krishtop et al., 2015;Silaev et al., 2015a,b;Virtanen et al., 2016). Such a heating can originate, for example, by an injected current even from the non-ferromagnetic electrode. This can be seen from the second term in the r.h.s. of Eq. (60), which shows that µ z can be finite even if the injection occurs from a non-magnetic lead, provided that N − = 0, i.e., if the superconductor shows a spin-split spectrum (Wolf et al., 2013). Finally, the heating is not sensitive to the sign of the bias voltage at the injector and hence the generated spin imbalance must be an even function of the voltage, µ z (V inj ) = µ z (−V inj ). This leads to an antisymmetric shape of the non-local spin signal in g nl with respect to V inj . This is also what has been observed in the experiments (Wolf et al., 2014a). All these features occur only if the superconductor has a spin-split density of states induced either by an external magnetic field or by the proximity to a ferromagnetic insulator.

Quantitative evaluation for the decay of the nonequilibrium modes
To give a more quantitative description of these effects we now calculate the potentials µ, µ z and the non-local conductance by using the kinetic equations for superconductors with spin-split subbands introduced in Sec. III.A.
Our starting point is the general Usadel equation (5). By assuming that the transparencies of the interfaces be-tween the superconductor and the electrodes are small we can neglect the changes of the spectral properties of the superconductor due to the proximity effect. This means that the retarded and advanced GFs correspond to the homogeneous superconducting state (cf. Sec. II.D), and we only have to focus on the calculation of the components of the distribution function, Eq. (56). In this discussion we disregard inelastic relaxation, assuming the corresponding scattering length to far exceed the spin relaxation length. We thus setΣ in = 0.
In the collinear situation considered here only the four components entering Eqs. (59-60) are finite. They are pairwise coupled. The first pair, f T and f L3 , determines the charge j c and spin-heat j se spectral currents The spectral coefficients appearing here are defined in terms of the spectral GFs discussed in Sec. II.D, and the diffusion coefficient D According to Eq. (5) the modes f T and f L3 satisfy the diffusion equations: where R T = 2∆Reg 01 , R L3 = 2∆Reg 31 , and the collision integral S L3 is given by Eq. (45). The kinetic equations (65,66) are supplemented by the boundary conditions at the injector electrode given by the second and the third rows on the matrix of Eq. (37).
On the other hand, the second pair of modes, f L , f T 3 , determines the energy and spin currents and satisfy the diffusion equations ∇ · j e = 0 (69) Here and the spin-relaxation term S T 3 is given by Eq. (44). Typical energy dependencies of diffusion coefficients in the presence of spin splitting and relaxation are shown in Appendix D.
Charge and spin-energy modes. The two systems of diffusion equations (65-66,69-70) describe the coupled transport of spin, charge and heat in superconductors. We start by solving the system (65,66) describing the coupled charge and spin energy modes. The solution for f T and f L3 is a superposition of two exponentially decaying functions where the coefficients A i have to be determined from the boundary conditions. The inverse lengths k 1,2 are energy dependent. This dependence is shown in Figs. 15(a-d).
Notice that both, the charge and spin-heat imbalance relaxation, is nonvanishing for all energies, below and above the gap due to the magnetic pair breaking effects (Nielsen et al., 1982;Schmid and Schön, 1975). In the absence of spin splitting, h = 0, R L3 = D L3 = 0 and the charge and spin-energy modes are decoupled [cf. Eqs. (65,66)]. Then k 1 and k 2 can be ascribed to the spin energy and charge imbalance relaxation, respectively. The non-zero Zeeman splitting h = 0 leads to the coupling of f T and f L3 , at low energies. At high energies ε ∆ though, when the superconducting correlations become negligible, the decoupled behavior is restored. This can be seen in the asymptotic behavior for → ∞ shown in Figs. 15(b,d). In this limit k 2 → 0 corresponding to a vanishing charge relaxation, whereas k 1 → λ −1 sn which is the spin-energy relaxation length in the normal state.
Spin and energy modes. We now analyze the system of equations (69-70) that describes the coupling between the spin and energy modes, i.e., the components f L and f T 3 of the distribution function. The solution of the system can be written as the sum of two qualitatively different terms with coefficients B j determined by the boundary conditions.
The first term in (74) describes a decay of the (spectral) spin imbalance with a characteristic length scale k s = S T 3 D L /(D 2 L − D 2 T 3 ). In the absence of a spinsplitting field and depairing mechanisms the expression for k se yields the energy dependent spin-relaxation length in superconductors obtained by Morten et al. (2005): k −1 s = λ sn √ ε 2 − ∆ 2 / ε 2 + β∆ 2 . This length is strongly renormalized by superconductivity if the spinflip scattering rate is nonzero, β = −1. A more accurate calculation of k se (ε), including the modification of spectral functions in the superconductor due to the spin In each plot β = −1; 1 corresponding to dominating spin-orbit and spin-flip scattering, respectively. The temperature is T = 0.1Tc0, normal-state spin relaxation rate (τsn∆0) −1 = 0.1. Fig. 15(c). We see that k s is constant for β = −1 and has a pronounced peak near the self-consistent gap edge for β = 1. Note that the gap edge is suppressed by spin-flip scattering to the values smaller than ∆ 0 .

relaxation, is shown in
Comparing Figs. 15(c,e) one can see that the relaxation lengths of charge imbalance and spin accumulation are very different even in the absence of a spin-splitting field. For example at subgap energies k s < k 2 . In contrast, for quasiparticles with energies above the gap the charge relaxation is suppressed at low temperatures whereas the spin relaxation remains finite. Hence particle injection above the gap leads at large distances from the injector mainly to a charge imbalance.
It is worth noticing that in the presence of a non-zero spin-splitting field the behavior of k s (ε) is very similar for spin-orbit and spin-flip relaxation mechanisms as shown FIG. 16 Nonlocal conductance as a function of the injecting voltage, g nl (Vinj) for α orb = 1.33, spin relaxation rate (τsnTc0) −1 = 0.2 and T = 0.05 Tc0, effective inelastic relaxation length L = 20λsn, L det = 5λsn and detector polarization P det = −0.5. Spin relaxation mechanism is (a) spinflip dominated β = 0.5; (b) spin-orbit dominated β = −0.9. The conductance is normalized to g0 = R ξ /(RinjR det ), where R ξ = ξ/(AsσN ) is the normal-state resistance of the wire with length ξ and cross section As.
We now focus on the second term in the r.h.s. of Eq. (74). It describes an effective increase of quasiparticle temperature associated with the f L nonequilibrium mode. This mode can only decay via inelastic scattering which is disregarded in the above analysis. In the following we assume the presence of an electrode at a distance L much larger than the spin relaxation length, and assume that the energy mode relaxes there.

Non-local conductance measurements
In a typical experiment quasiparticles are injected from a normal or a ferromagnetic metal electrode by applying a voltage V inj between the electrode and the superconductor. In this way the nonequilibrium states described by Eqs. (73,74) can be excited. We determine the coefficients A 1,2 and B 1,2 using the general bound-ary conditions (37) for the currents, where the voltagebiased normal electrodes are described by the distribution functions f (N ) We plot the non-local conductance in three different cases (the three curves in Fig. 16).
The blue curve is for a vanishing spin-splitting field and injector polarization. In that case the signal is solely due to charge imbalance. The black curve is with a nonzero spin-splitting field, but vanishing injector polarization, and in the red curve both the spin-splitting field and the injector polarization are non-zero. As discussed below, the former shows the thermally created spin accumulation, whereas the latter shows a combination of the thermally and electrically created spin accumulations.
In an experiment where the spin-splitting field is controlled by an external magnetic field, we also need to take into account the orbital depairing effect [Eq. (12)] of the magnetic field in addition to the Zeeman effect. Whereas the latter leads to the coupling of the energy and spin modes, the former affects especially the relaxation of charge imbalance. For the results in Fig. 16, we describe the relative strength of the orbital depairing and the spinsplitting field by setting (τ orb T c0 ) −1 = α orb (µ B B/T c0 ) 2 and choosing the dimensionless quantity α orb to a value close to the experiments in (Hübler et al., 2012).
First, in the absence of a Zeeman field, N − = 0, and according to Eqs. (58-60) only the modes f T and f T 3 contribute to g nl . In such a case the contribution stemming from the spin accumulation is finite only if P inj = 0, which is the condition to obtain a finite f T 3 . However, this function decays over the spin diffusion length and therefore is negligibly small at the distances L det > λ sn = √ Dτ sn from the injector. Thus, the detected signal in this case is mostly determined by the charge imbalance µ. This explains the approximate symmetry with respect to the injecting voltage: g nl (V inj ) = g nl (−V inj ) (blue curves in Fig. 16). The charge imbalance contribution to g nl grows monotonically when |eV inj | > ∆ g . This behavior is determined by the increase of the charge relaxation scale at high energies, k −1 2 → ∞ shown in Fig. 15(c). In the presence of an applied magnetic field the charge relaxation is strongly enhanced due to the orbital depairing. As a result, the charge imbalance background signal is strongly suppressed by an increased h. On the other hand the spin imbalance contribution stemming from the f L mode is large. This term describes heat injection in the presence of a finite spin-splitting field h and has a long-range behavior. This contribution leads to the large peaks in g nl (V inj ) shown in Fig. 16. The non-linear heating produced by quasiparticles injected at voltages exceeding the energy gap explains the large electric signal observed in the experiments (Hübler et al., 2012;Quay et al., 2013;Wolf et al., 2014a). Notice that the peaks do not have exactly the same form so that g nl (V inj ) = −g nl (−V inj ). The deviation from the perfect antisymmetric form is due to the small but finite injector polarization P inj , which modifies the boundary conditions for the spin current.
One important feature shown in Fig. 16 is that the spin polarization peaks exist always in the presence of spin splitting, h = 0, even if the injector electrode is not ferromagnetic, P inj = 0. This explains the non-local conductance measurements performed by Wolf et al. (2013) by using a normal metal as injector.
The results summarized in Fig. 16 can have direct applications in the field of spintronics. On the one hand spin accumulation can be created without using ferromagnets as injectors. On the other hand, due to the coupling between the spin and energy mode, the spin signal can be controlled over a very long range. The only requirement for these two features to occur is to use a superconductor with a spin-split DOS due to a Zeeman or an exchange field.

C. Spin Hanle effect in normal metals and superconductors
In the previous section we study the spin injection and spin accumulation in a superconductor assuming that the magnetization of the injector and the applied field are collinear. We now lift this assumption and consider non-collinear field and magnetizations. This situation has been widely studied in the normal state. The noncollinearity between the external field and the spin of the injected electrons leads to a precession of the latter. This is the spin Hanle effect (Jedema et al., 2001(Jedema et al., , 2003Johnson and Silsbee, 1985;Villamor et al., 2015;Yang et al., 2008). This precession can be measured via the non-local conductance in a multi-terminal sample as a function of the applied field.
In order to understand the spin Hanle effect in superconductors it is instructive to first discuss the spin precession in a normal metal within the formalism described above. We then demonstrate that in the superconducting state the spin precession can be either enhanced or suppressed, depending on the spin relaxation mechanism.

Normal-metal non-local spin valve
In the normal state the Usadel equation can be drastically simplified because the retarded and advanced Green's functions do not depend on energy and have the simple formǧ R(A) = ±τ 3 . Thus the reduced density of states equals to unity according to Eq. (15). It follows from Eq. (21) that the Keldysh component is then proportional to the distribution matrix:ǧ K = 2τ 3f .
It is then useful to integrate Eq. (5) over energies and take the trace after multiplication with the Pauli matrices. This results in a diffusion equation for the nonequilibrium spin density S given by This is the spin diffusion equation (Dyakonov and Khaetskii, 2008;van Son et al., 1987;Zhang, 2000) widely used in spintronics. Generalization of this equation in the presence of spin-orbit coupling has been done in several works (Duckheim et al., 2009;Mal'shukov and Chao, 2005;Mishchenko et al., 2004;Stanescu and Galitski, 2007). Its right hand side describes the effective spin relaxation time defined as τ −1 sn = τ −1 so + τ −1 sf . The second term on the left hand side describes the torque induced by the external field via the spin splitting field. This torque leads to the spin Hanle effect studied in the language of quasiclassics by Hernando et al. (2000).
Let us consider a lateral spin valve as the one sketched in Fig. 14(b). We assume that a field is applied in the x-direction whereas the injector and detectors are polarized in the z-direction. If the width and thickness of the N wire are smaller than the spin relaxation length the problem reduces to a quasi 1D geometry along the length of the wire. In such a case Eq. (75) leads to two linear coupled equations for the components S y,z , We have defined the spin diffusion length λ 2 sn = Dτ sn , the magnetic length l 2 h = D/h and κ = λ −1 sn 1 + i(λ sn /l h ) 2 . The coefficients A and B have to be determined from the boundary condition, Eq. (35), which in this particular case has the simple form Here V I is the voltage across the injector/N interface and R I the resistance per unit area of the barrier. The resulting spin accumulation in the N wire is where V I and R I are the voltage drop and the resistance of the injector, respectively. A ferromagnetic voltage detector at a distance x = L from the injector detects the spin potential µ z = S z (L)/ν F . The resulting µ z (B) dependence, described by Eq. (76), coincides with the Hanle-shaped curves measured experimentally (Jedema et al., 2001(Jedema et al., , 2003.

Hanle effect in spin-split superconductors
If the wire in Fig. 14 is in the superconducting state the non-local signal can change drastically as discussed by Silaev et al. (2015b). In contrast to the normal-metal case, the retarded and advanced Green's functions have a non-trivial energy dependence and therefore one cannot integrate straightforwardly the diffusion equation over the energy.
As in previous sections, we assume that the retarded and advanced GFs are position independent and the density of states is the one shown in Fig. 5. In contrast to the normal case the non-local Hanle signal in a superconducting wire depends on the dominating spin relaxation mechanism. We consider again the two mechanisms: magnetic impurities, Eq. (11), and extrinsic spin-orbit coupling, Eq. (10), and introduce the parameter β describing the relative strength between them β = (τ so − τ sf )/(τ so + τ sf ).
In principle the Hanle signal in the superconductor would include contributions from the long-range energy mode discussed in Sec. IV.B for the case of collinear field and magnetizations. To separate that contribution from the bare Hanle effect, we concentrate on the linear response regime of low injector and detector voltages, where the effect of the energy mode can be disregarded (cf. Fig. 16). In such a case the non-local conductance is exponentially small and therefore we focus the non-local resistance V det /I inj , where V det is the voltage at the detector in the absence of a current and I inj is the current at the injector. To uncover the Hanle effect, we furthermore study the difference V S between the voltages for the parallel and anti-parallel orientations of the detector and injector polarizations. The non-local resistance of interest is thus R S = V S /I inj .
If we assume that the external field is applied in zdirection and disregard the orbital effects, the retarded quasiclassical Green's function is described by the four components in Eq. (16). The energy dependence of these components is shown in Fig. 4.
Whereas the spectral terms remain constant in the S field, the spin-polarized current I inj that flows through the interface with the injector causes a spin accumulation characterized by the vector µ s . The current measured at the detector is then given by (Silaev et al., 2015b) where R det = R /A is the normal-state barrier resistance between the superconducting wire and the detector with cross sectional area A and spin polarization P det , and V det is the voltage measured at the detector with respect to the superconductor. Y = ∞ 0 dεN + ∂ ε f eq describes the effect of the density of states in Eq. (40). Its lowtemperature analytical estimate for weak spin relaxation is given in Eq. (117) below.
To find the voltage induced in the electrically open detector circuit we set I det = 0 in Eq. (77). The spin-dependent part of the non-local resistance R S = [V det (P det ) − V det (−P det )]/I inj can be found using the linear response relation for the injector current I inj = V inj Y /R inj with the normal-state resistance R inj of the injector: where V inj is the voltage applied at the injector in order to inject the current. We then need to determine the spin accumulation µ s (x) induced along the S wire which can be written in terms of the Keldysh component of the GF as with m(ε, x) = Tr(τ 3 Sg K )/8. For simplicity we assume that P inj ⊥ h so that the f L mode is not induced by the thermoelectric coupling at linear response, see Sec. V.D.2. Therefore the spin accumulation does not contain the long-range contribution discussed in Sec. IV.B. The remaining part of m in Eq. (79) is given by: where the vector f T = (f T 1 , f T 2 , f T 3 ) contains the spin components of the general distribution function (22).
The first term in this expression is the usual quasiparticle contribution which also appears in the normal state. It is determined by the nonequilibrium distribution function at energies larger than the superconducting gap, that is where the total density of states is non-zero. In contrast, the second term in Eq. (80) only appears in the superconducting case and is finite for subgap energies due to the prefactor Img 33 (cf. Fig. 4). This term can be related to the difference between the DOS in spin subbands using the Kramers-Kronig relation for the retarded GF Let's assume that the spin-splitting field in the superconductor points in z-direction, h = he z . Then the transverse components f T 1 , f T 2 satisfy the kinetic equations obtained from Eq. (5), In contrast to the normal case, Eq. (75), the diffusion coefficient is now a tensor with components that depend on energy, The terms on the right hand side of Eq. (81) are X 1 = (S T 1 − H 1 ), X 2 = (S T 2 + H 2 ), with S T 2 = 2τ −1 sn (Img 33 Reg 03 − βImg 01 Reg 31 ).
The H 1,2 terms are the "Hanle" terms that describe the coherent spin rotation and relaxation due to the action of the external field. In the normal case only H 2 and the first term of S T 1 are non-zero and we recover Eq. (75).
We describe the spin-relaxation, Eqs. (85,86), by two parameters: the normal-state spin relaxation time τ sn and the relative strength parameter β as in Sec. II.C.
In order to solve Eqs. (81a-81b) we need the boundary conditions at the interface with the injector. They can be obtained from Eq. (35). By keeping terms to leading order in the interface parameter κ I = 1/(R I σ N ) we obtain (Silaev et al., 2015b): These are evaluated at the position of the interface. Equations (81a,81b) can be re-written for the spectral density of spin polarization, Eq. (80), in a more familiar form, similar to the Bloch-Torrey transport equation (Torrey, 1956) for the magnetic moment: Here γ = −2 is the electron gyromagnetic ratio and j sk is the spin current density in the k-th spatial direction.
The diffusion coefficients are defined as Equations (88-89) are written in the steady-state limit, which is enough for the description of the spin Hanle effect. In the dynamical case, the magnetization depends explicitly on time m = m(ε, t), and the l.h.s. of Eq. (88) contains additional terms. When the spin dynamics is slow enough, with a characteristic frequency much smaller than the superconducting energy gap, ω ∆, the Bloch-Torrey equation in the superconducting state can be written as 10 In Eq. (88), τ S and h are the transverse spin relaxation time and a correction to the effective Zeeman field appearing in the presence of spin-relaxation processes. The latter determines the (Larmor) precession frequency of the spins. Specifically, Both the transverse relaxation rate τ −1 S and the effective field shift h s = |h s | are proportional to the normal-state spin relaxation rate τ −1 sn . In the superconducting state τ −1 S and h s depend differently on energy for different spin-relaxation mechanisms. The typical dependencies are shown in Fig. 17 for β = ±1. One can see that τ −1 S has a step-wise behavior as a function of energy when the spin-orbit relaxation dominates [ Fig. 17(a)]. In contrast, τ −1 S (ε) shows peaks near the spin-split gap edges for spinflip relaxation [ Fig. 17(b)]. The effective field shift h s has different signs for β = ∓1 as shown in Fig. 17(c,d). Interestingly, the transverse spin relaxation τ −1 S vanishes at the subgap energies. Equations (81a-81b) imply that the transverse components of f T have the general form where A is an integration constant determined by the boundary conditions, and with Rek T > 0. In contrast to the normal metal case considered in Sec. IV.C.1, there is no straightforward relation between k T and the parameters τ S and h s . This is because the diffusion coefficients in the superconducting state entering Eq. (89) are energy dependent and hence the spectral spin current has a complicated expression that couples the different components of the spin polarization. The real part Rek T determines the spin relaxation length and the imaginary part Imk T gives the quasiparticle spin precession. In the superconducting state both these scales are energy dependent.
In Fig. 18 we show how the precession and relaxation scales depend on energy and on the parameter β. It follows that at low energies Re k T is larger than in the normal case which corresponds to the limit ε ∆. At ε = 0 the characteristic relaxation length k −1 T (0) is determined by the coherence length ξ = D/∆ 0 rather than by the normal-state spin relaxation length λ sn . In contrast, the imaginary part of k T vanishes at the energies below the spectral gap where the DOS shown by the dashed lined in Fig. 18 is absent N + = 0. In the case of dominating spin-orbit scattering (β < 0), the precession Imk T has a shallow peak at ε ≈ ∆ + h [see Fig. 18(a)]. If the spin-flip mechanism dominates (β > 0), this peak of Im k T is suppressed, Fig. 18(b).
Putting all this together allows us to evaluate the nonlocal resistance R S exhibiting the Hanle curves in the superconducting state as shown in Fig. 19. The data at T = T c corresponds to the normal-metal Hanle signal. The result depends on the relative directions of the injector and detector polarizations. Panels a and c show the non-local signal R Sy obtained from Eq. (78) when the detector polarization P det = P det y so that P det P inj ⊥ h. One can see that both, the precession and decay of the nonlocal signal, disappear at T → 0, whereas the shape of the curves at intermediate temperatures depends on the type of spin relaxation.
If instead of the above configuration we assume that the three vectors (P det , P inj , h) are perpendicular to each other (e.g., P det = P det x, P inj = P inj y, h = hz), the subgap current is absent in the detector circuit and the corresponding spin signal R sx has a strong dependence on h even at the temperatures well below T c . This is shown in Figs. 19(c,d).
We are not aware of the measurement of the spin Hanle signals in the superconducting state. We emphasize that the above results are obtained in the linear response regime of small voltages. They disregard the thermal effects coming into play at higher voltages, especially at  those of the order of the energy gap ∆ ± h. In the linear response regime the transport quantities in superconductors are exponentially suppressed due to the gap in the density of states, and hence the above theory is rather made for the thermally activated quasiparticles. Because of this, the best way to uncover the superconducting effects on the Hanle response would be to study it rather close but below the critical temperature T c .
To summarize: The kinetic equations for the nonequilibrium modes constitute a rather generic theoretical tool that we have used in this section to provide a quantitative description of some transport properties, such as spin injection and precession, occurring in superconducting structures in the presence of a spin-splitting field. The predictive capability of the derived kinetic equations has been proven in this section by contrasting the results with existing experiments especially in the case of collinear magnetizations. This makes these equations an ideal tool to study further effects that involve the coupling between different nonequilibrium modes in superconductors, as for example, the possible thermal spin Hanle and other non-linear effects taking place in systems with non-collinear magnetization.
In the next section we focus on thermoelectric effects that are also direct manifestations of the spin-splitting fields in the superconductor. In contrast to this section where the focus is on the diffusion of the nonequilibrium modes within the superconducting wire, there we concentrate especially on interface effects, and disregard the position dependence of the nonequilibrium modes within the superconductors.

V. THERMOELECTRIC EFFECTS
Thermoelectric effects relate temperature differences to charge currents, and electrical potentials to heat currents. The traditional view of thermoelectric effects in superconductors is that if they exist, they must be very weak. In bulk superconductors, this is partially because any thermoelectrically generated quasiparticle current is easily screened by a supercurrent (Meissner, 1927). Alternatively, one could then measure this supercurrent via an additional constraint to the phase of the superconducting order parameter in multiply connected structures (Ginzburg, 1944). However, even this thermally created phase gradient tends to be weak, owing to the near-complete electron-hole symmetry in superconductors (Galperin et al., 1974).
In fact, superconductors do contain some ingredients for strong thermoelectric effects, because the latter typically require strongly energy dependent density of states of the charge carriers. This is provided by the BCS density of states. Hence, if one can break the electron-hole symmetry of the transport process via some mechanism, superconductors can become very strong thermoelectrics. This is precisely what happens in spin-split superconductors, as an exchange field breaks the symmetry in each spin sector, but so that the overall spin-summed energy spectrum remains electron-hole symmetric. As discussed in Sec. I and III.D, transport through a spin filter to a spin-split superconductor then can provide large thermoelectric effects because the two spins are weighed differently (Machon et al., 2013(Machon et al., , 2014Ozaeta et al., 2014). This prediction was also recently confirmed experimentally (Kolenda et al., 2017(Kolenda et al., , 2016b. In this section we give an overview of the different types of thermoelectric effects discussed for superconductors, and then concentrate on the ones obtained in superconductors with a spin-splitting field. We discuss both the linear response regime and beyond it, and consider also the limiting features such as the electron-phonon coupling. We show that under suitable conditions, in particular for close-to-optimal spin filters, the efficiency of thermoelectric conversion can become very large and exceed that obtained for best thermoelectric devices operated at or above room temperature. Besides the regular quasiparticle current, in devices coupling two superconductors, one with a spin-splitting field and one without, the thermoelectric effect can also be converted to a phase gradient (Giazotto et al., 2015a). We discuss a prediction for such a large thermophase in Sec. V.E.
As these effects require low operating temperatures, they obviously cannot be directly used to improve the efficiency of various everyday devices. However, this effect may become relevant for other types of applications, such as sensors, where the measured (wide-band) signal consists of heating one part of the system . The thermoelectric conversion can then be used to convert the resulting temperature difference to a charge current.
A. Thermoelectric effects and heat engines Thermoelectric effects are typically described via the linear response relation between charge and heat currents I,Q and bias voltage and temperature difference V and ∆T across a contact The 2 × 2 conductance matrix in Eq. (97) is a part of the generalized 4 × 4 Onsager matrix (39) connecting different interface currents and potentials as discussed in Sec. III.B. Here we consider a conventional situation when the spin-dependent potentials V s and T s are negligibly small. This is relevant in the case of large tunneling barriers. This approximation allows the reduction of the general boundary condition (39) to the simpler one (97) describing thermoelectric response.
With a non-zeroα, electrical energy may be converted to heat or cooling, or reciprocally a temperature difference may be converted to electrical power. The efficiency of this conversion can be described by constructing a model for a generic electrical heat engine (see Fig. 20). There, a load with resistance R L is driven by the power drawn from the thermoelectric element across which there is a temperature difference ∆T . The power dissipated on the load is P = IV = R L I 2 , whereas the voltage across the thermoelectric element is −IR L . Plug-ging this into Eq. (97) yields The extracted power hence is On the other hand, the thermoelectric element extracts heat from the temperature difference, i.e., trying to balance it, with the poweṙ The efficiency of thermoelectric conversion is hence where N =α 2 /(G th GT ) ≤ 1. Alternatively, we can write N in terms of the usual thermoelectric figure of merit, where S =α/(GT ) is the thermopower (Seebeck coefficient) andG th = G th −α 2 /(GT ) is the thermal conductance at a vanishing current. Now we should choose R L to optimize the device somehow. For example, the maximum efficiency is obtained This result is consistent with that obtained by Snyder and Ursell (2003) in the linear response limit T cold /T hot ≈ 1. On the other hand, optimizing the device to yield a maximum power output requires GR L = 1, corresponding to the limit (Curzon and Ahlborn, 1975;Novikov, 1957) Here η CA ≡ 1 − T cold /T hot ≈ ∆T /(2T ) is the maximum efficiency obtained when ZT → ∞. Both of these efficiencies are maximized when the thermoelectric figure of merit becomes large. At or above room temperature, the record-high figures of merit are obtained in certain strongly doped semiconductor structures (Kim et al., 2015;Zhao et al., 2016). A typical record value for those cases is ZT 1 . . . 2. This value and the optimal temperature where it is obtained results from a competition of two generic temperature dependencies: that of phonon heat conductance, and that of the (typically activated) process yielding the thermopower. Both of these decrease towards low temperatures, but the previous decreases as a power law, whereas the latter decreases exponentially. In an electron system the heat conductance is typically close to the Wiedemann-Franz limit G el th ≈ LGT , where L = π 2 k 2 B /(3e 2 ). The total heat conductance is obtained from the sum of the electronic and phononic contributions, G th = G el th + G ph th . Since ZT ∝ GT /G th = 1/[L + G ph th /(GT )], the only way to improve this is to minimize the phononic contribution. Generally the latter decreases towards low temperatures faster than linearly. However, the thermopower is typically an exponential function of temperature, S ∝ exp(−∆/k B T ) as the best thermoelectrics contain a gapped dispersion with the gap ∆. Therefore, the optimal ZT takes place at a temperature which is some fraction of ∆. As shown in Sec. V.D (see in particular Fig. 29), the same is true for a normal metal island coupled to spin-splitted superconductors. In that case, the magnitude of ∆ is just orders of magnitude lower than in semiconductors structures.

B. Thermoelectric effects in superconductors
Research on thermoelectric effects in superconductors dates back to the 1920's, when Meissner (1927) concluded them to be absent because any thermoelectric current is cancelled by a counterflowing supercurrent. Ginzburg (1944) showed that this is no longer the case in multiply connected bimetallic superconducting structures. This is because the presence of supercurrent is linked to the gradient of the phase of the order parameter, and in multiply connected structures uniqueness of the order parameter imposes a constraint on the relation between the flux through the ring and the circulating supercurrent. This situation can be visualized as in Fig. 21. For definiteness, let us consider a bimetallic superconducting loop, where the thermoelectric effects mostly take place at the two contacts. The current through contact i = 1, 2 consists of a sum of the supercurrent and the thermoelectric current, i.e., In a closed circuit this current must equal for both junctions, and thereby it yields a relation between the phases ϕ 1,2 . A second relation fixing the phases is obtained in the presence of a flux Φ through the loop with inductance L = L T + L B . Without loss of generality the total loop inductance can be described with a single quantity. Namely, the phases must be fixed so that they minimize the energy where n ∈ Z, Φ is the magnetic flux through the junction and Φ 0 = h/(2e) is the flux quantum. In the following we denote Φ/Φ 0 = f +2π(m−n) with f ∈ [0, 2π[ and m ∈ Z. For f 1 and I th,i I c,i , the induced phases are small, and we may linearize the current-phase relations. As a result, we obtain a circulating current where L tot = L + (I −1 c1 + I −1 c2 )/2e is the total inductance of the superconducting loop.
Let us assume heating the lower superconductor of the bimetallic loop so that the temperature increases from T to T + ∆T . Such a heating induces thermoelectric currents I th,1 =α 1 ∆T /T, I th,2 = −α 2 ∆T /T.
This produces a circulating thermoelectric current It can be measured for example by placing a SQUID on top of the bimetallic loop, and measuring the induced flux Φ ind = M I circ , where M is the mutual inductance between the two systems. The size of the coefficientα was calculated by Galperin et al. (1974) where the latter form is due to the reduction of the quasiparticle density in the superconducting state, andα N is the size of the thermoelectric coefficient in the normal state. The precise value ofα N depends on the exact electronic spectrum of the metals in question and its calculation needs extending the theory beyond the quasiclassical approximation employed in this review, as within that approximationα N = 0. For example, for a simple quadratic dispersionα N = , where E F is the Fermi energy. At temperatures T ∆/k B ,α is thus expected to be a product of two small coefficients, α N ∝ k B T /E F , and G(∆/T ). This is very small and not easy to measure quantitatively. Because of many spurious effects in such measurements, it is not simple to make the experiments agree quantitatively with this theory.
Nevertheless, at least close to the critical temperature the thermoelectric flux should be observable, and the first experiments to measure it were done in the 1970s and early 1980s. One set of experiments (Katsovnik et al., 1981;Ryazanov and Schmidt, 1981) followed the idea of Aronov and Gal'perin (1974) and measured the thermoelectric voltage across a superconductor-normalmetal-superconductor contact as the thermoelectric current exceeded its Josephson critical current. The results of these experiments are in line with the expected magnitude of the thermoelectric signal. On the other hand, Van Harlingen et al. (1980) measured the thermoelectrically generated flux in a bimetallic loop formed from superconducting Pb and Ti, close to the critical temperature of the latter. These experiments are also discussed by Galperin et al. (2002). Surprisingly, the experiments demonstrated fluxes five orders of magnitude larger than predicted by theory. This discrepancy annoyed Ginzburg so much that he devoted an entire chapter on the topic in his Nobel colloquium (Ginzburg, 2004).
Very recently, Shelly et al. (2016) claim to have solved this discrepancy with new experiments performed on much smaller superconducting loops than what was possible in the early 1980s. According to them, the discrepancy originated from the temperature dependence of the inductances L tot and M as well as the flux f (via the effective area of the loop), and hence they change as one end of the bimetallic loop is being heated. This produces additional contributions to the thermoelectric flux, The additional contributions are in practice much larger than the pure thermoelectric effect. According to Shelly et al. (2016), whereas the effects from f can be accounted for by measuring the period of oscillations as the external flux is altered, the effects due to the trapped flux, accounted for by m, is much harder to deduce, and was likely the reason for the discrepancy, as in the early measurements a geomagnetic field amounted to m ∼ 10 6 . Using much smaller loops, Shelly et al. (2016) were able to control the number of trapped flux quanta, and hence get rid of the spurious effects. The remaining thermoelectric flux that they observed is more or less in accord with the value obtained from the above theory. However, the temperature dependence of the measured flux depends, besides G(x) above, also on the temperature dependent inductance, and on the temperature dependence of the main heat contact between the heated electrons and the phonon bath, so the measurement could not directly deduce G(x). The conclusion from these theory and experimental works is that in conventional superconductors thermoelectric effects can be finite, but they are extremely weak, and therefore difficult to access. In the following section we show how this situation can be completely reversed in the presence of the spin-splitting field, provided we add a second ingredient into the theory: spin filtering. In this case such superconductor/ferromagnet hybrids can become almost ideal thermoelectric devices.
We note here that other types of superconducting heterostructures besides the bimetallic loop do contain relatively strong thermoelectric effects. In particular, a supercurrent flowing along a temperature gradient leads to an appearance of charge imbalance (Clarke and Tinkham, 1980;Pethick and Smith, 1979;Schmid and Schön, 1979), which can be measured (Clarke et al., 1979) for example via a non-local geometry similar to those discussed in Sec. IV. A similar type of an effect was found in Andreev interferometers (Chandrasekhar, 2009;Eom et al., 1998;Parsons et al., 2003) composed of a hybrid multiterminal geometry of a normal metal in contact with a superconducting loop. These effects were considered theoretically by Seviour and Volkov (2000); Titov (2008); Virtanen and Heikkilä (2004); and Virtanen and Heikkilä (2007). Similar effects are also predicted in ballistic systems as discussed by Jacquod and Whitney (2010) and Kalenkov and Zaikin (2017).
Quite recently Zaikin (2014, 2015a); Machon et al. (2013Machon et al. ( , 2014; and Ozaeta et al. (2014) showed another mechanism for creating large thermoelectric effects in superconducting proximity structures via coupling such structures to ferromagnets. In particular, Machon et al. (2013Machon et al. ( , 2014 showed how threeterminal proximity-coupled superconductor-ferromagnet devices can show non-local thermoelectric effects: in this case the density of states (DOS) in a normal metal coupled both to a superconductor and a ferromagnet becomes spin dependent, and the spin-resolved DOS is also electron-hole asymmetric, resulting in the strong thermoelectric effect. On the other hand, Zaikin (2014, 2015a) showed how a metallic bilayer consisting of two superconductors or a superconductor and a normal metal, separated by a spin-active interface (i.e., an interface whose transmission properties are characterized by a spin-dependent scattering matrix), can exhibit large thermoelectric response. These mechanisms are closely related to the one where a superconductor with a spinsplitting field is tunnel coupled to a normal metal via a spin-polarized interface (Ozaeta et al., 2014). We explain this mechanism in more detail below.
Another way to affect the thermoelectric response via magnetism was discussed by Zaitsev (1986) and later by Kalenkov et al. (2012), who argued how magnetic impurities inside a superconductor enhance the thermoelectric coefficient by a large factor k F compared toα in Eq. (110). Here is the elastic mean free path and k F is the Fermi wavenumber. This effect results from the electron-hole asymmetric Andreev states forming in the vicinity of the magnetic impurities.
The above discussion concerns metallic structures. In semiconductor quantum dots the thermoelectric effects can be large when the spectrum of the quantum dot is electron-hole asymmetric (Beenakker and Staring, 1992), even without the presence of superconductivity or magnetism. Theoretical studies showed that the combination of the latter affect the symmetries and might increase the thermoelectric response in quantum dots coupled to superconducting and magnetic electrodes (Hwang et al., 2016a,b).

C. Thermoelectric effects at a spin-polarized interface to a spin-split superconductor
Consider a tunnel contact to a superconductor in a spin-splitting field as described in Sec. III.C. Let us assume that the tunnel contact is magnetic, so that the conductance through it is spin-polarized. Denoting the spin-dependent conductances (in the normal state) by G ↑ , G ↓ we can parameterize them by the total conductance G T = G ↑ + G ↓ and the spin polarization P = (G ↑ − G ↓ )/G T . The total tunneling quasiparticle charge and heat currents are now expressed as a sum over spindependent contributions, each of the form The formula for the current is thus the spin-resolved version of Eq. (1). Here f L/R = n F (E − µ L/R ; T L/R ), n F (E; T ) = {exp[E/(k B T )] + 1} −1 are the (Fermi) functions of the reservoirs biased at potentials µ L/R and temperatures T L/R . We assume N R (ε) spin independent for simplicity; the possible splitting of the Fermi sea in electrode R is included in G σ . The reduced density of states in the superconductor for spin σ is N σ (ε). The heat cur-rentQ i σ is calculated separately for i = R, S, using the potential µ R/S , because the two heat currents differ by the Joule power I(µ R − µ S )/e. Let us denote the spindependent reduced density of states via N + = N ↑ + N ↓ and N − = N ↑ −N ↓ . In that case the spin-averaged tunnel currents are In the following, we first analyze these currents in detail, and then discuss a method of building a near-optimal heat engine based on such junctions. In particular, we use the densities of states described in Sec. II.D. In most of the analysis below, we disregard the spin relaxation effects on the density of states, because this assumption allows for some analytically treatable limits. This assumption is lifted in Fig. 31.

Heat current outside linear response
Let us consider first the case when the reservoir R is a normal metal and therefore N R (ε) = 1. The heat current from reservoir R is a non-monotonous function of voltage even in the absence of spin polarization or temperature difference. In particular, for voltage V = (µ R − µ S ) ≈ ∆/e, it is positive, i.e., reservoir R cools down (Leivo et al., 1996;Nahum et al., 1994;Pekola et al., 2004). However, that heat current is quadratic in the voltage, and therefore it does not result from the usual Peltier effect (Eq. (97) forQ) where the cooling power is linear in voltage. In the presence of spin polarization P and with a non-zero spin-splitting field h in the superconductor, the cooling power can be slightly improved.
The cooling power from reservoir R as a function of voltage is shown in Fig. 22 for various values of h, assuming the ideal case of unit spin polarization P = 1.
In an electron refrigerator, the cooled element is an island coupled to two electrodes. Let us consider a normalmetal island playing the role of the reservoir R, and coupled to two spin-split superconductors via ferromagnetic insulators with polarizations P L and P R , respectively. The optimal situation is realized when P L = −P R = ±1 or when the exchange fields in the two superconductors are reversed. The total cooling power from the island iṡ Q N =Q R (V /2, P L ) +Q R (−V /2, P R ). It works against other relaxation mechanisms, so that the stationary temperature T of reservoir R is determined from heat balance (Giazotto et al., 2006), Here we assume that the dominant heat relaxation mechanism on the reservoir R with volume Ω N is due to electron-phonon interaction with strength Σ N , described by Eq. (51). In addition, we assume that the spin accumulation on the island, produced by the nonequilibrium eV =" driving, is negligibly small due to spin relaxation. The assumptions relevant for this limit are discussed in Sec. V.F below. The resulting island temperatures are shown in Fig. 23 for different parameters of the S-FI-N-FI-S junction. We have chosen the parameters so that for h = P = 0 they correspond to a typical T (V ) curve found experimentally [see for example (Leivo et al., 1996)]. In the absence of spin filtering, P = 0, temperature obtains a minimum at eV ≈ 2(∆ − h), and T (V ) is a symmetric function of voltage. Generally, spin splitting in this case makes the cooling worse, so that the minimum reached temperature is higher than in its absence.
The behavior of T (V ) changes in the presence of spin filtering, P = 0. In particular, the curve becomes nonsymmetric, and there is cooling even in the linear response regime, i.e., low voltages. This Peltier effect is discussed more below. In addition, the minimum reached electron temperature is generally lowered by an increasing P , and for a large P ≈ 1, the lowest temperature may be obtained at a non-zero exchange field. However, this effect seems rather weak for the considered parameters.
One possibly relevant aspect of such a magnetic cooler is the fact that the optimum temperature is obtained at a lower absolute value of the voltage (Kolenda et al., 2016a). This translates into a somewhat lowered Joule power injected to the device. As this heat is dumped into the superconductor, at the lowest temperatures the heating of the superconductor becomes the dominant limiting obstacle for cooling instead of the electron-phonon coupling (Kauppila et al., 2013;Nguyen et al., 2014). The consequences of this are analyzed in (Rouco et al., 2017). Contrary to the spin-independent case, the magnetic element can be used to refrigerate the superconductor (Rouco et al., 2017). We illustrate this by plotting the cooling power of the superconductor in the case of an N-FI-S junction in Fig. 24. Without spin splitting and a spin-polarized interface, the cooling power is always negative, i.e., corresponding to heating. However, since the eV =" Peltier cooling with a non-zero P is linear in voltage, it has to result to a non-vanishing cooling power also from the superconductor. The maximum cooling power is less than in the case of cooling the normal metal (Fig. 22).
However, as indicated in Fig. 25, the resulting temperature drop is nevertheless appreciable because the energy gap also weakens the electron-phonon coupling (see Eq. (52)).
D. Linear response: heat engine based on a superconductor/ferromagnet structure As can be seen in Figs. 22 and 24, the simultaneous presence of the non-vanishing spin polarization P and a spin-splitting field h lead to a heat current that has a linear component in the voltage V . This component is nothing else than the Peltier effect, and the whole linear response can be described within the scheme in Sec. III.D. In the limit k B T ∆ − h the linear-response coefficients evaluate to (Ozaeta et al., 2014) At low temperatures the thermopower is maximized for It can hence become much larger than the "natural scale" k B /e, and even diverge towards low temperatures. However, such a divergence comes together with the vanishing of the conductance, Eq. (117), and therefore is in practice either cut off by circuit effects, where the impedance to the voltmeter becomes lower than the contact impedance, due to spin relaxation neglected above, or alternatively by additional contributions beyond the BCS model. The latter ones are described in more detail in (Ozaeta et al., 2014). Nevertheless, with proper circuit design one should be able to measure a thermopower much exceeding k B /e in this setup. The above theoretical predictions in the linear response regime were recently confirmed experimentally by the group of Detlef Beckmann (Kolenda et al., 2017(Kolenda et al., , 2016b. In particular, they prepared a sample containing a crossing of three types of metals, a normal-metallic Cu, ferromagnetic Fe, and superconducting Al (see Fig. 26). The measured configuration is also sketched in Fig. 27(a). The electrons in the ferromagnetic wire were heated with FIG. 26 Scanning electron micrograph of the sample used in (Kolenda et al., 2016b) for the measurement of the thermoelectric effects. From Kolenda et al. (2016b).
FIG. 27 a) Schematic setup for measuring the thermoelectrically induced current, used in (Kolenda et al., 2016b). S, F, and N stand for a superconductor, ferromagnet and a normal metal, whereas FI is a ferromagnetic insulator. b) Setup used for a direct measurement of the Seebeck effect. c) Heat engine realized in a lateral setup with "n-doped" and "p-doped" junctions using a FNF trilayer with antiparallel magnetization directions. To disregard spin accumulation, the island has to be large compared to the spin relaxation length. d) Heat engine with a spin-split superconducting island. The ferromagnets can also be replaced by a normal metal if the interfaces to the superconductor contain a ferromagnetic insulator. In (c) and (d), the heating power P heat is partially converted to "useful" work P work dissipated on the load. the heater current I heat , producing a temperature difference between the ferromagnet and the superconductor. The contact between the ferromagnet and the normal metal is ohmic and therefore the temperature difference between them is negligible small. Then the thermoelectric current was measured as a function of the magnetic field B applied parallel with the ferromagnetic 43 FIG. 28 Thermoelectric current as a function of the applied magnetic field, measured in (Kolenda et al., 2016b). The circles show the measurement values, the solid lines show a comparison to Eq. (114). The three solid lines correspond to slightly different temperature differences; for further details, see (Kolenda et al., 2016b). From Kolenda et al. (2016b). wire. The agreement between the experimental results and the above described tunneling theory was excellent (see Fig. 28). The temperature difference between the ferromagnet and the superconductor was a fitting parameter, whereas the polarization P was fitted from the conductance spectrum. In the experiment it was fitted to the value P = 0.08, a modest value attributed to the thin oxide barrier between the Fe and the Al layers. In principle larger values of P can be obtained by increasing the thickness of the oxide barrier (Münzenberg and Moodera, 2004), but this of course would reduce the amplitude of the thermoelectric current.
In the experiment, the thermoelectric current was measured rather than the voltage. In that case the impedance of the sample dominated that of the measurement lines. This is why the measurement yielded the exponentially low thermoelectric current, which nevertheless was large enough to be measured. The measurement configuration in Fig. 27(b) would have directly measured the generated voltage drop (i.e., Seebeck effect) instead of the current. This voltage results from the ratio of two exponentially small functions, the thermoelectric coefficientα and the conductance G, and itself is not small. Such a measurement would then tell about spurious effects, for example due to spin relaxation, or due to the presence of fluctuations or states inside the gap, limiting the diverging Seebeck coefficient at low temperatures (Ozaeta et al., 2014). Better still, replacing the normal metal with another superconductor with an inverse spin-splitting field, would have resulted to twice as large signal (corresponding to a series of p-and n-doped thermoelectric elements), but would not be possible to create as such with a magnetic field. The solution would be furthermore to replace the ferromagnetic wire by an FNF heterostructure [Fig. 27(c), where the ferromagnets have antiparallel magnetizations, for example due to different coercive fields, and the normal metal N would serve as a spacer between them]. To reach high figures of merit, the ferromagnetic metals should also replaced by ferromagnetic insulators, which can reach very high values of spin polarization (see Table I), with P exceeding 0.9. The island setup in Figs. 27(c) and (d) also realizes a thermally isolated structure, in contrast to those in (a) and (b). This allows realizing a heat engine, where the voltage measurement is replaced by the "device" to be powered with the engine, with resistance matched to the thermoelectric element as discussed in the beginning of the section. In the following, let us assume that it is only the electrons of the ferromagnetic island that are heated, instead of both the electrons and phonons. 11 In this case the spurious heat conduction mechanisms would arise from the electron-phonon coupling, with heat conductance G e−ph = 5ΣΩT 4 [see Eq. (51)], and from the heat conductance to the extra superconducting electrodes S, with an exponentially suppressed heat conductance, provided the heating current is small. Let us also assume the electron-phonon coupling dominates the other spurious contributions. In this case the heat engine has an efficiency described by a figure of merit that for k B T h and P = 1 tends to where g = 5k 5 B √ 2πe 2 ΣΩ∆ 3 /(2G T ) is a dimensionless quantity characterizing the strength of electron-phonon coupling. Here we take into account the fact that the thermoelectric effect takes place across both contacts, and hence all quantities except the spurious heat conduction channel are doubled. If we do not describe halfmetals, i.e., 1−P 2 (h+∆) 2 /(h−∆) 2 exp[−2h/(k B T )], we can simplify this more by dropping the last term in the denominator. We plot this as a function of temperature in Fig. 29 for some representative parameters. For example, for Ω = 0.005 µm 3 , Σ = 10 9 W µm −3 K −5 and 1/G T = 3 kΩ, g = 100. The figure of merit exhibits a peak at around k B T ≈ 0.1∆, and its peak value depends strongly on the relative strength of the spurious heat conduction channels. As seen in the figure, reaching high values for the figure of merit is quite challenging in this way, as besides a high polarization close to one, it requires quite low normal-state resistance of the junctions, a combination not very easy to reach with ferromagnetic insulators.
Note that it is really the presence of the spurious electron-phonon heat conduction that limits the highest available values of ZT . Often such spurious mechanisms are disregarded from the theoretical analysis, for example in the case of quantum dots (Hwang et al., 2016a).
Another possible setting for the heat engine is the one where the island is the spin-split superconductor, and it is connected to normal-metal or ferromagnetic electrodes via spin filters as in Fig. 27(d). In the linear response regime this is otherwise similar to the previous case, except that the electron-phonon coupling inside the island has an exponentially suppressed heat conductance, Eq. (52). In this case, as long as the pure tunneling limit remains valid, and the island volume is not overly large, the low-temperature behavior of ZT is dictated mostly by the deviation of the polarization from unity as shown in Fig. 30. However, because of the exponentially decaying heat conductance, also in this case other spurious processes besides electron-phonon coupling start to limit ZT , making it vanish at T → 0. In the case of Fig. 30, the simplest model for them is due to the small but nonvanishing density of states assumed for the superconductor (Ozaeta et al., 2014), and described by the Dynes param- The figure has been calculated with h = 0.5∆ and g = 1000, without calculating ∆ self-consistently. The figure of merit at low temperatures reaches very close to P 2 /(1 − P 2 ) unless P is very close to unity, but the exact temperature scale where this happens depends on the value of polarization. At the lowest temperatures ZT is limited by another spurious heat conduction process, due to nonzero density of states inside the gap, described here by the Dynes Γ parameter eter Γ (Dynes et al., 1984). 12 Rather than the Dynes parameter, a more relevant limitation of the figure of merit in spin-split superconductors is due to the spin mixing caused by spin-orbit and spinflip scattering that are disregarded in the above results. This becomes an issue especially in heavy-metal based superconductors where the spin-orbit scattering is strong. As an example, Fig. 31 shows how the figure of merit decays as the spin-orbit relaxation becomes stronger. The high efficiencies are obtained only when the spin relaxation rates are small compared to the gap energy ∆. Note that for some materials the stronger spin-orbit relaxation due to the increase in the atomic number Z (say, compared to Al) is partially compensated by the increased critical temperature and therefore the energy gap.
(Color online) Temperature dependencies of the figure of merit ZT (T ) = α 2 /(GthG − α 2 ) calculated according to the ns Eq. (39,40,41) in the Review draft in the presence of spin-orbital scattering β = −1. Curves from top to bottom d to the spin relaxation times 1/(τsnTc0) = 0; 0.5; 1.0; 2.0; 4.0; 6.0; 10; 25. Even if the true figure of merit of the type of heat engine discussed above can be made high, these systems cannot obviously be used to replace room-temperature thermoelectric devices to be applied for example in energy harvesting. However, there are other applications where the large figure of merit may turn out to be essential. For example, this type of thermoelectric heat engine can be used for thermal radiation sensing at low temperatures (Giazotto et al., 2006;Heikkilä et al., 2017). Another possible use of the thermoelectric effects would be in non-invasive low-temperature thermometry, where the temperature (difference) profiles could be read from the thermopower, without having to apply currents. In a scanning mode this would hence be a low-temperature version of the method used in (Menges et al., 2016).

Nonlinear response heat engine
The above calculations on the heat engine efficiency have been done in the strict linear response limit ∆T T . In this limit both the output power and the Carnot efficiency, proportional to ∆T /T , are vanishingly small. However, the device can have a rather large output power and efficiency also at nonlinear response. To illustrate this, we plot the maximum extractable output poweṙ W max = max V (−IV ) from a N-FI-S junction in Fig. 32 at large differences T S − T N between the spin-split superconductor and a normal reservoir. The output power is somewhat larger for T S > T N than the opposite situation with the same temperature difference. The plotted quantity corresponds to the maximum power that is possible to extract from the junction after optimizing the load resistance -note, however, that this optimized resistance value depends on the exact temperature difference and other parameters. The corresponding efficiency η =Ẇ max /Q is plotted in Fig. 33. There we take into account only the heat current flowing through the junction k B T =" itself. The efficiency is reduced if other heat relaxation mechanisms become relevant.

Effect of non-collinear magnetizations
Above we assume that the polarization P = Pû P of the tunneling contacts is collinear with the exchange field h = hû h in the spin-split superconductor. This is the limit that guarantees the largest thermoelectric current, but it is also interesting to consider how non-collinearity would show up in the thermoelectric response, in particular in the coefficientα of Eq. (119). This coefficient is by definition a scalar quantity. Therefore, it can only depend on a scalar combination of the vectors P and h, or rather,û P andû h . This combination is the inner product,û P ·û h . Therefore, the generalization of Eq. (119) to the non-collinear case is whereα c is the thermoelectric coefficient in the collinear case, for example given by Eq. (119).

E. Thermophase in a S(FI)S contact
We begin the discussion of thermoelecric effects in superconductors above by considering the case of thermally induced phase gradients, or thermally induced circulating currents in bimetallic multiply connected structures. Similar type of setting in spin-split superconductors was theoretically investigated in (Giazotto et al., 2015a). The total current in this case consists of the sum of a thermoelectric current I th and the supercurrent, where I th is obtained from (114) and I c is the critical current for the junction with a phase difference ϕ of the order parameters across the contact. The critical current is obtained from (Bergeret and Giazotto, 2014) where M L (ε) = [F L (ε + h) + F L (ε − h)]/2, F L,R (ε) = ∆ L,R / (ε + i0 + ) 2 − ∆ 2 L,R , and f L,R (ε) = tanh[ε/(2k B T L,R )]. As shown in (Bergeret and Giazotto, 2014), the exchange field affects the precise value of the critical current.
In an electrically open configuration, the two currents must cancel, and instead a thermophase ϕ th develops across the junction. This is obtained from The thermophase can be detected as described in Subs. V.B, using a bimetallic loop with two contacts, characterized by critical currents I c1,2 and thermophases ϕ th 1,2 . For non-zero exchange field and spin polarization P , the resulting thermophases can be much larger than in ordinary bulk superconductors. Hence the temperature dependence of the inductances play a more minor role. For junctions with non-equal thermophases and for negligible loop inductance (in practice, 2eLI c1,2 / 1) in the absence of an external flux the circulating current is In the case of symmetric junctions both thermophases are the same and the circulating current in the absence of an external flux vanishes. However, as discussed in (Giazotto et al., 2015a), the thermoelectric current affects the response of the circulating current to the external flux, allowing for their measurement also in that case. Equation (126) requires that both sides of the equation have an absolute value of at most unity, i.e., |I th | < I c . For a very large thermoelectric current, its cancellation with a supercurrent is no longer possible, and instead a voltage across the contact forms. In this case the direct current response of the junction is more similar to the case discussed above in the linear response limit for a N-FI-S junction. This regime was investigated in detail by Linder and Bathen (2016). Moreover, the nonvanishing dc voltage across the superconducting junction leads to Josephson oscillations at the frequency 2eV /h, where h is the Planck constant. Hence, the device can be used as a temperature (difference) to frequency converter as discussed in more detail in (Giazotto et al., 2015b).

F. Thermally induced spin currents
Besides the large thermoelectric effect discussed above, the contact between spin-split superconductors with other conducting materials can exhibit a large (longitudinal) spin Seebeck effect, where a temperature difference drives spin currents to/from the spin-split superconductor (Ozaeta et al., 2014). The spin current I S = I ↑ − I ↓ can be obtained from Eq. (112). In linear response, Eq. (39), it is described by the same coefficient α/P = α as the thermoelectric coefficient, Eq. (119), but without the spin polarization P of the interface. When either of the two materials realizes an island, the spin current can be converted into a spin accumulation µ z that is determined from the balance between thermally induced spin currents and spin relaxation within the island. For weak spin relaxation, the resulting spin Seebeck coefficient µ z /∆T can be large [i.e., not proportional to the number of thermally excited quasiparticles ∼ e −(∆−h)/(k B T ) as the terms in Eqs. (117)-(119)], similar to the thermopower in Eq. (121). The spin currents induced in the case of two spin-split superconductors, and the additional effects of Josephson coupling, magnetization texture and spin-orbit effects are discussed in (Bathen and Linder, 2017;Linder and Bathen, 2016).
The above discussion of the thermoelectric cooling and heat engines disregards the spin accumulation caused by the nonequilibrium biases. This assumption is valid if the overall spin relaxation within the island, described by a spin-flip conductance G sf = e 2 ν F Ω/τ sn is much larger than the conductances of the tunnel junctions. Here ν F is the normal-state density of states in the Fermi level, Ω is the volume of the island, and τ −1 sn is the spin relaxation rate within the island. In the case of a spin-split superconducting island, both the conductance and the spin relaxation rate contain the same exponential term ∼ e −(∆−h)/(k B T ) , and for an estimate of the overall magnitude of the effect, it is hence enough to compare G sf in the normal state to the normal-state tunnel conductance G T . In other words, disregarding the spin accumulation is justified provided e 2 ν F Ω/(τ sn G T ) 1. On the other hand, most of the discussion disregards the effect of spin relaxation on the density of states of the spin-split superconductor. As shown in Fig. 31, spin relaxation starts to affect the results when /τ sn k B T c . These two constraints can be simultaneously satisfied depending on the precise value of τ sn , provided that G T (k B T c /δ I )e 2 /h, where δ I = 1/ν F Ω is the average energy level spacing on the island. This constraint is always satisfied in tunnel barriers between metallic systems.
To our knowledge, the spin Seebeck effect has not been directly measured in spin-split superconductors driven by a temperature difference. However, outside linear response, this is the effect that causes the long-range spin signal in spin injection experiments as discussed in Sec. IV.B. Besides affecting the nonequilibrium dynamics, the presence of spin accumulation in a superconductor with an exchange field couples to the self-consistency equation for the superconducting order parameter ∆ as in Sec. III.H. Bobkova and Bobkov (2017) utilized this fact and predicted that the thermally induced spin accumulation in a spin-split superconductor can result to changes in the critical temperature: besides only suppressing it at low temperatures, it can enhance the critical temperature at some intermediate temperatures, or even lead to a situation where superconductivity shows up only between two critical temperatures.
This spin Seebeck effect should be contrasted to the analogous phenomenon discussed in nonsuperconducting materials (Uchida et al., 2014). There, the major contribution to the spin Seebeck signal is typically due to the thermally induced spin pumping from ferromagnetic insulators (Hoffman et al., 2013).

VI. NONEQUILIBRIUM DYNAMICS OF SPIN-SPLIT SUPERCONDUCTORS
The conduction electrons of a superconductor in an oscillating electromagnetic field can absorb energy via excitation of the existing quasiparticles, and by creation of new quasiparticles from breaking of Cooper pairs. In a stationary state, the excitation is balanced by corre- sponding relaxation processes. These are sensitive to the spectrum of states available in a spin-split superconductor.
Electrons couple to the electromagnetic field via the Zeeman and the orbital terms. The Zeeman coupling is the source for the conduction electron spin resonance (CESR), which was theoretically considered in early works (Aoi and Swihart, 1970;de Gennes, 1966b;Kaplan, 1965;Tsallis, 1977;Yafet, 1983). The corresponding dynamical susceptibility for the spin-split case was discussed by Kosov and Kochelaev (1978); Maki (1973); and Tagirov and Trutnev (1987). The linewidth of the resonant absorption peak is determined by the spin relaxation time in the normal state, τ sn , and is not affected by the slow relaxation of the thermal long-range spin imbalance discussed in Sec. IV.B. The magnetic screening complicates the experimental observation of the resonance. Results however exist in the mixed state of type II superconductors as reported by Nemes et al. (2000); Vier and Schultz (1983); and Yafet et al. (1984), and more recently, also in spin-split Al thin films by Quay et al. (2015) [see Fig. 34]. For the Al films, linewidths of 100 ps were observed, consistent with the expected magnitude of τ sn in aluminum.
The CESR physics is also related to pumping effects, where an externally driven spin precession drives a nonequilibrium state or currents in the superconductor, or currents across a junction. Effects of this type have been considered in different ferromagnet/superconductor structures in several works, in adiabatic (Brataas and Tserkovnyak, 2004;Skadsem et al., 2011) and nonadiabatic cases (Holmqvist et al., 2011;Houzet, 2008;Richard et al., 2012;Trif and Tserkovnyak, 2013).
As mentioned above, a rf field also couples to the orbital degree of freedom of the electrons. This coupling drives ac currents that also excite quasiparticles. In the absence of spin splitting, the linear and nonlinear response of superconductors has been extensively studied (Tinkham, 1996). In the linear response regime, the complex impedance for spin-split superconducting films was considered theoretically and experimentally by van Bentum and Wyder (1986) adding Zeeman energy shifts to BCS theory results (Abrikosov et al., 1959;Mattis and Bardeen, 1958). In the nonlinear regime, strong driving can modify the observed quasiparticle spectrum . Moreover, excited quasiparticles can lead to an increase of the superconducting gap ∆ (Eliashberg, 1970;Ivlev et al., 1973;Klapwijk et al., 1977). In spinsplit superconductors, this effect is modified by the spinsplit density of states and spin-flip scattering .
In this section we discuss the dynamic response of spinsplit superconductors in terms of the time-dependent quasiclassical equations, both in the linear and nonlinear regimes.
In addition to the effects discussed in this section, oscillating fields can excite magnetic impurities and nuclear moments. Discussion of these extrinsic resonance effects can be found in the articles by Baberschke (1976); Barnes (1981); Garifullin (2015); and Taylor (1975), and are beyond the scope of this review.
A. Linear response of a spin-split superconductor to a rf field The dynamic nonequilibrium response of spin-split dirty superconducting thin films can be studied on the basis of the Usadel equation, Eq. (5), which in timedependent representation we write in a compact form Nowǧ(ω, ω ) = ∞ −∞ dt dt e iωt−iω t ǧ(t, t ) depends on both time indices and matrix products involve energy convolutions, (a • b)(ω, ω ) = dω1 2π a(ω, ω 1 )b(ω 1 , ω ). Above, (ω, ω ) = 2πδ(ω − ω )ω, and∇X = ∇X − i[Aτ 3 • ,X]. External classical electromagnetic fields couple to the electrons via the exchange field h(ω − ω ) and the vector potential A(ω−ω ). We choose the gauge such that the local charge neutrality (Artemenko and Volkov, 1979;Gor'kov and Kopnin, 1975) fixes the scalar potential to φ(t) = π 4 Trǧ K (t, t), and in the spatially uniform cases discussed below φ = 0. As above, we consider thin films, and neglect magnetic screening. A review of the quasiclassical formulation for time-dependent problems can be found in a book by Kopnin (2001).
The kinetic equation for the distribution functionf follows from Eq. (128) and the above assumptions, where the collision integral is defined aŝ andẐ Note that this definition of collision integrals differs from Eq. (33) due to the presence of the vector potential. In a spatially uniform case with A = 0, the Eq. l.h.s of130 vanishes and the kinetic equation reduces to the condi-tionÎ = 0. The CESR emerges when considering a Zeeman field of the form h(ω) = h 0ẑ + h 1 (ω)x, where an ac field h 1 perpendicular to a static exchange field h 0 excites the electron spins. The spin dynamics follows a Bloch-like equation including spin relaxation, similar to the spin Hanle effect (see Sec. IV.C and Eq. (88)). Such equations can be derived from the Usadel equation in linear response to h 1 , and from other approaches (Inoue et al., 2017;Kosov and Kochelaev, 1978;Maki, 1973;Tagirov and Trutnev, 1987). Here we briefly discuss the physics within the Usadel framework.
It is useful to separate the corrections to the Keldysh Green's function in Eq. (128 due to the distribution function from those due to the modification of the spectral functions: δĝ K = δĝ K reg + δĝ K an , defining δĝ K reg = δĝ R • f eq − f eq • δĝ A . We split δX K similarly. The equation [X • ,ǧ] K = 0 then reduces in linear response tô The solution to Eq. (133) can be obtained with the Ansatz δĝ K an =ĝ R eq • δf − δf •ĝ A eq , where δf = f T x σ x + f T y σ y . To determine the solution components, we can take the trace 1 8 Tr[(. . .)σ] of the (ε, ε−ω) frequency component of Eq. (133). This results in where m = 1 8 Tr[τ 3 σδĝ K an (ε, ε−ω)] describes the nonequilibrium magnetization (spin accumulation), m = 1 8i Tr[τ 1 σδĝ K an (ε, ε − ω)] describes a correction to spin scattering by superconductivity, and N + = [g 03 (ε) + g 03 (ε − ω) * ]/2, F + = [g 01 (ε) − g 01 (ε − ω) * ]/(2i) are finitefrequency generalizations of the spin-averaged density of states and anomalous functions, respectively. The term I = ih 1 (ω)[N +x − N zŷ ][f eq (ε) − f eq (ε − ω)] describes the exciting field. There is redundancy in the representation, which enables writing m in terms of m, resulting in the final Bloch equation. Analogously to the spin Hanle effect discussed in Sec. IV.C.2, the spin relaxation time τ S is renormalized with respect to the one in the normal state, τ sn , and there is also a correction h s to the Zeeman field: where N z = [g 33 (ε)−g 33 (ε−ω) * ]/(2i) and F z = [g 31 (ε)+ g 31 (ε − ω) * ]/2. For slow driving, ω → 0, the spin relaxation time τ S and the Zeeman field correction h s coincide with those visible in the Hanle effect (cf. Eq. (93) and Fig. 17), which also involves precession of the transverse spin component. The result (135) describes resonant excitation of the transversal modes, f T j , which correspond to a nonequilibrium contribution to the spin accumulation (30), where δµ reg s arises from the modification of the spectral functions, δg R/A . The result contains the conduction electron spin resonance peak, δµ s ∝ h 1 (ω)/[4(h + h s ) 2 − (ω + i/τ S ) 2 ] at frequency ω 2|h 0 |. As discussed in section IV.B, the f T j modes can relax due to elastic spin-flip scattering, which determines the peak absorption linewidth 1 τ S ∝ τ −1 sn . The result by Yafet (1983) can be obtained from Eq. (135) in the quasiequilibrium approximation δf (∂ E f eq )µ s · σ and in the limit h 0 ∆. Equation (135) also coincides with the result by Maki (1973).
If we focus on the orbital effects, the complex impedance of the superconductor, i.e., the linear response to an oscillating electric field described by A(ω), can also be obtained within the Usadel framework. In the above chosen (London) gauge, the perturbation to time-averagedǧ is of the order D|A(ω)| 2 /ω (Eliashberg, 1970) and can be neglected in linear response at nonzero frequency. The charge current response is then given by where R σ (ε, ε ) = N σ (ε)N σ (ε ) + Im g σ,1 (ε) Im g σ,1 (ε ) , (142) R σ (ε, ε ) = Im g σ,3 (ε)N σ (ε ) − Re g σ,1 (ε) Im g σ,1 (ε ) , and g ↑/↓,1 = (g 01 ± g 31 )/2, N ↑/↓ = Re[g 03 ± g 33 ]/2. The kernels R σ have the BCS form (Abrikosov et al., 1959;Mattis and Bardeen, 1958), and are decoupled for each spin species. In other words, the electric field couples to the orbital motion of the electrons and, in our approach, it conserves spin. Notice that the derivation of the Usadel diffusion equation neglects spin-orbital effects when dealing with the vector potential, and assumes the elastic scattering time τ el as the shortest time scale. The prefactors in Eqs. (139)-(141) are given by the normalstate conductivity, σ N = 2e 2 ν F D. Phenomenological modification of the above by allowing spin flips in the coupling matrix element was considered by van Bentum and Wyder (1986), who obtain a reduction in the pairbreaking threshold frequency from 2∆ to 2∆ − 2h (c.f. Fig. 35). Spin-orbit impurity scattering with τ so ∼ ∆ −1 in the model used above also reduces the threshold, but by an amount that depends on τ so .

B. Nonlinear spin imbalance of a spin-split superconductor in a rf field
We now go beyond the linear response regime and focus on the stationary quasiparticle distribution which enters the observables. This is determined by disregarding the gradients, in which case the collision term in Eq. (130) is equal to the contribution from the time-dependent vector potential. This results in This rate equation, similar to that used by Grimaldi and Fulde (1997), describes the balance between excitation of quasiparticles induced by ac fields (Î ac ), the spin relaxation by the elastic spin-flip scattering (Î sn ), and the inelastic relaxation (Î in ). Each term in Eq. (144) can be described within the quasiclassical approach. The orbital ac term of the collision integral can be obtained by assuming a time-dependent spatially uniform vector potential A(t) = A 0 cos(ωt). Here, it is useful to simplify the problem by expanding in DA 2 0 /ω 1, and reduce it to a form that only involves the dc component of the distribution function. In this case (Eliashberg, 1970) I ac,σ (ε) = DA 2 for the I ↑/↓ = Tr[(1 ± σ z )Î/2] components. We assume that all spin dependent fields are collinear and hence we can define f σ=↑/↓ = f L ± f T 3 , cf. Eq. (22). Notice that the charged modes f T , f L3 are not excited by a uniform A(t).
The absorption kernel R is the one appearing in the real part of the conductivity, Eq. (142). Qualitatively, this collision integral results in the depletion (accumulation) of quasiparticles in an energy band of width ω above (below) the gap edges |E| > |∆ ± h| (see Fig. 35). At higher energies the electron distribution corresponds to an increased temperature.
The collision integrals for the relaxation processes can be written as, cf. Eq. (44), where we describe inelastic scattering in a relaxation time approximation. The spin-conserving microwave absorption term does not directly generate spin imbalance in the superconductor. However, as illustrated in Fig. 35, spin relaxation converts the accumulation of quasiparticles above the gap to an imbalance in the number of spin-up and spin-down quasiparticles. From the kinetic equation (144) we obtain in the limit 1/τ sf 1/τ in , The first term contributes zero total spin imbalance (60), but the second gives a nonzero contribution of the order of τ in /τ sn . A related effect in quasiparticle injection was discussed by Grimaldi and Fulde (1997). The generated spin imbalance can be in principle measured by a spin-polarized probe junction, as discussed in Sec. IV.A. Note, however, that in this type of experiments other processes may need to be also considered. On the one hand, photoelectric signals can also occur in the absence of spin splitting as shown by Kalenkov and Zaikin (2015b) and Zaitsev (1986). These processes, however, generally scale with ∝ τ el rather than ∝ τ in . On the other hand an ac bias over the S/F junction may also result in rectification, as discuss in the next section.
It is interesting to notice that the nonequilibrium quasiparticle accumulation generated by the microwave absorption affects the magnitude of the the superconducting gap ∆, leading to deviations from the results by Eliashberg (1970). Results for the self-consistent ∆(T ) from numerical calculations  are shown in Fig. 36 for h = 0 and h > 0. In the absence of spin-splitting, h = 0, a gap enhancement (Eliashberg, 1970) occurs. For h > ω, an additional instability develops at ∆(T ) = h, corresponding to a coexistence of two solutions where the spin-averaged DOS is either gapless (|∆| < h) or gapped (|∆| > h). The instability requires the presence of strong spin-flip scattering: for τ sn τ in the exchange field does not cause significant qualitative changes, because the quasiparticle accumulations above the gap edges that generate the enhancement, are approximatively described as two independent copies of the effect at h = 0.

C. S/F tunnel junction dynamics
For adiabatic excitation, ω ∆ dynamics of S/F tunnel junctions can be described by the dc relations discussed in Sec. III and V. That this limit can be reached at time scales shorter than the spin-relaxation and inelastic times, was used in a recent experiment (Quay et al., 2016) to probe the relaxation dynamics in spin-split thinfilm superconductors.
At higher frequencies, photoassisted tunneling breaks the adiabatic description. This can be taken into account via a standard tunneling Hamiltonian approach (Larkin and Ovchinnikov, 1967;Werthamer, 1966). The result for the tunneling current is where T is the tunneling matrix element, G >(<) = 1 2 G K ± 1 2 (G R − G A ) are (non-quasiclassical) Green functions for the superconducting (k) and nonsuperconducting (q) sides, and the phase φ(t) = e t dt V (t ) is related to the time-dependent voltage V (t) across the junction. The above result applies for spin-conserving tunneling with collinear magnetizations. If the terminals are in an internal equilibrium (no spin or charge imbalance), standard approximations and changing integration variables yields the time-dependent generalization of Eq. (1), − e iV (t)(t−t ) ] , where I dc is given by Eq.
The resulting current-voltage relation is asymmetric, I dc (−V ) = − I dc (V ) for h = 0. This implies that such junctions rectify ac signals (Quay et al., 2016). For an ac signal V (t) = V dc + V ac cos(ωt), from the above we have the average dc current (Tucker and Feldman, 1985) I dc = At V dc = 0 and small V dc at uniform temperature T h, ω, the rectified current is I dc ∝ P F V 2 ac h/(R T ∆ 2 ),  proportional to the exchange field in S. Measurements attempting to use F probes for the nonequilibrium ac effects discussed in the previous sections need to take into account such rectification, or try to suppress V ac e.g. via a suitable microwave circuit design (Horstman and Wolter, 1981).
While dynamical effects in superconductors have been studied extensively experimentally, we are aware of relatively few studies on spin-split superconductors. Tunnel junction rectification effects were observed for example by Quay et al. (2016). However, to the best of our knowledge, higher-frequency experiments probing gap enhancement or photoassisted tunneling have not been reported so far.

VII. SUMMARY AND OUTLOOK
This review focuses on transport and thermal properties of superconducting hybrid structures with a spinsplit density of states. Such a splitting can be achieved either by an external magnetic field, or, more interestingly, by placing a ferromagnetic insulator (FI) adjacent to a superconducting layer (S) (Sec. II). We discuss several experimental situations with the help of a theoretical framework (see Sec. II.C and III.A) based on the quasiclassical formalism, with which one can account for both thermodynamical and nonequilibrium properties of such hybrid structures. In order to account for effects beyond quasiclassics, as for example strong spin polarization, we combine the quasiclassical equations with effective boundary conditions.
Out-of equilibrium superconductivity by itself leads to a decoupling between the charge and energy degrees of freedom of the electronic transport. In this review we show that the combination between superconductivity and magnetism leads on one hand to additional nonequilibrium modes, spin and spin-energy, and on the other hand to a coupling between them all. This leads to novel and intriguing phenomena discussed in this review with direct impact in latest research activities and proposed future technologies based on superconductors and spin dependent-fields (Eschrig, 2011(Eschrig, , 2015Linder and Robinson, 2015). By using the theoretical formalism presented in this review one can predict and explain phenomena such as the spin injection and relaxation (Sec. III.E and IV) in superconductors with an intrinsic exchange field and their consequence in the transport properties (Sec. IV.B-IV.C). We also discuss a number of striking thermoelectric effects in superconductors with a spinsplitting field (Sec. V).
The best scenario for the phenomena and applications discussed here, and in particular for the thermoelectric effects, are FI-S systems where the spin splitting can be achieved without the need of an applied magnetic field. Hence it becomes important to look for ideal FI-S material combinations. So far europium chalcogenides (EuO, EuS and EuSe) together with Aluminum films have shown large splittings and hence these are the best combination. In addition, thin films of EuO or EuS can be used as almost perfect spin filters (see Table I) and hence they are good candidates for realizing the near-optimal heat engines proposed in Sec. V. One of the main challenges from this perspective is to find FI-S combinations with large superconducting critical temperature and simultaneously a large spin splitting. Superconductors like Nb or Pb on the one hand increase T C with respect to Al-based structures, but on the other the spinorbit coupling may spoil the sharp splitting as discussed in Sec. II.D. Recent experiments on GdN-NbN suggest large splittings (Pal and Blamire, 2015) but further research in this direction is needed.
In Sec. VI we discuss the dynamics of spin-split superconductors in rf fields, and review how nonequilibrium effects induced by them manifest in the quasiclassical theory. Historically, magnetic resonance effects in superconductors are well studied, but fewer experiments have probed spin-split thin films.
Besides the effects discussed in this review, several theoretical studies made striking predictions in mesoscopic systems with spin-split superconductors, such as the creation of highly polarized spin currents (Giazotto and Bergeret, 2013b;Giazotto and Taddei, 2008;Huertas-Hernando et al., 2002), large supercurrents in FI-S-I-S-FI junctions (Bergeret et al., 2001a), junctions with switchable current-phase relations (Strambini et al., 2015), and an almost ideal heat valve based on S-FI elements (Giazotto and Bergeret, 2013a).
Although many of the transport phenomena in spinsplit superconductors are now well-understood, we foresee a number of exciting avenues for future research. We outline such research directions below.
One further perspective of the present work is the extension of the Keldysh quasiclassical formalism used in this review to include magneto-electric effects associated with the spin-orbit coupling (SOC). For a linear in momentum SOC the generalization of this can be done by introducing an effective SU(2) gauge potential. The quasiclassical equations in this case have been derived by Bergeret and Tokatly (2013, 2016. Effects such as the spin-Hall and spin-galvanic effect in superconductors have been studied in the equilibrium case (Konschelle et al., 2015). Extending these results to a nonequilibrium situation, and also to time-dependent fields, would be an interesting further development and would allow for a detailed study of the well-controlled non-linearities associated to these effects in superconductors. First steps in this direction have been taken in (Espedal et al., 2017).
Recent discoveries of skyrmionic states in chiral magnets (Nagaosa and Tokura, 2012) have attracted a lot of attention due to the effects resulting from the interplay of magnetism and SOC (Soumyanarayanan et al., 2016) which can induce chiral Dzyaloshinskii-Moriya interactions between magnetic moments. Currently it is very interesting to study these effects in the presence of the additional component -superconductivity, when the exchange interaction is mediated by the Cooper pairs (de Gennes, 1966c). One can expect that in such systems superconductivity can induce a non-trivial magnetic ordering and dynamics. These effects can show up in various systems including ferromagnet/superconductor bilayers, surface magnetic adatoms and bulk magnetic impurities inducing the localized Yu-Shiba-Rusinov states modified by the SOC (Pershoguba et al., 2015).
Superconducting structures with strong spin-orbit coupling and exchange fields are also of high interest in view of engineering a platform for realization of topological phases and Majorana bound states (Alicea, 2012;Beenakker, 2013;Hasan and Kane, 2010;Qi and Zhang, 2011). Understanding and controlling the behavior and relaxation of nonequilibrium quasiparticles in these systems is also of importance, not least because of their influence on the prospects of solid-state topological quantum computation (Nayak et al., 2008).
This review focuses exclusively on the nonequilibrium properties of superconductors in proximity to magnets. We expect the inclusion of the magnetization dynamics and its coupling to the electronic degrees of freedom via the reciprocal effects of spin transfer torque and spin pumping (Tserkovnyak et al., 2005) in the farfrom equilibrium regime to lead to completely new type of physics, as the two types of order parameters affect each other. The coupling of supercurrent on magnetization dynamics and texture has been studied during the past decade (Houzet, 2008;Richard et al., 2012;Waintal and Brouwer, 2002), but the work where both systems are out of equilibrium has been mainly concentrated on Josephson junctions (Hikino et al., 2011;Holmqvist et al., 2014Holmqvist et al., , 2011Kulagina and Linder, 2014;Mai et al., 2011) and much less attention has been paid to quasiparticle effects (Linder et al., 2012;Skadsem et al., 2011;Trif and Tserkovnyak, 2013).
Besides the rich physics offered by spin-split superconductors, they have been long used as tools to characterize equilibrium properties of magnets, especially their spin polarization. In this review (see end of Sec. V.D) we outline two further possibilities related to their large thermoelectric response: accurate radiation sensing and non-invasive scanning thermometry. We believe there are also many other avenues to be uncovered, opened by the possibility for realizing a controlled combination of magnetism and superconductivity.