Intrinsic Paramagnetic Meissner Effect due to s-wave Odd-Frequency Superconductivity

In 1933, Meissner and Ochsenfeld reported the expulsion of magnetic flux, the diamagnetic Meissner effect, from the interior of superconducting lead. This discovery was crucial in formulating the Bardeen-Cooper-Schrieffer (BCS) theory of superconductivity. In exotic superconducting systems BCS theory does not strictly apply. A classical example is a superconductor-magnet hybrid system where magnetic ordering breaks time-reversal symmetry of the superconducting condensate and results in the stabilisation of an odd-frequency superconducting state. It has been predicted that under appropriate conditions, odd-frequency superconductivity should manifest in the Meissner state as fluctuations in the sign of the magnetic susceptibility meaning that the superconductivity can either repel (diamagnetic) or attract (paramagnetic) external magnetic flux. Here we report local probe measurements of faint magnetic fields in a Au/Ho/Nb trilayer system using low energy muons, where antiferromagnetic Ho (4.5 nm) breaks time-reversal symmetry of the proximity induced pair correlations in Au. From depth-resolved measurements below the superconducting transition of Nb we observe a local enhancement of the magnetic field in Au that exceeds the externally applied field, thus proving the existence of an intrinsic paramagnetic Meissner effect arising from an odd-frequency superconducting state.


I. INTRODUCTION
Below the superconducting transition of a conventional (s-wave) Bardeen-Cooper-Schrieffer (BCS) superconductor such as Nb, the electrons stabilise into Cooper pairs in a spin-singlet state meaning that the electrons of a pair have oppositely aligned spins. The screening supercurrent density (J) that is generated by a superconductor in response to a weak magnetic field is linearly proportional to the vector potential (A) via the density of pairs present ns (J=-e 2 nsA/mc, where c, e and m are the speed of light, the electron charge and the electron rest mass, respectively). Consequently, the amplitude of the screening supercurrent density (J) is negative and a diamagnetic Meissner effect is observed [1].
At the surface of an s-wave superconductor proximity coupled to a magnetic metal, the exchange field of the magnet can induce an odd-frequency superconducting state in which the Cooper pairs are in a spintriplet state with a density (nt) that is a mixture of spin-zero and spin-one pair projections [11][12][13]. This means that J is dependent on the magnitude and sign of ns-nt (i.e. J=-e 2 (ns-nt)A/mc) [14] and so odd frequency triplets should act to reduce the screening current [6]. Since ns and nt have different decay envelopes in an exchange field, J should reverse in sign as a function of magnetic layer thickness when nt exceeds ns at which point the magnetic susceptibility is positive and an inverse -paramagnetic -Meissner effect [6][7][8][9] prevails.
Evidence for spin-triplet pairing has recently been demonstrated in experimental studies involving magnetically inhomogeneous superconductor/ferromagnet (S/F) hybrids, such as via transition temperature measurements of S/F1/F2 spin valves [15], long-ranged supercurrents in S/F/S Josephson junctions [16], and various spectroscopy measurements of F/S systems [17].
To investigate the Meissner effect in a superconductor-magnet system, we measure the depth profile of the local magnetic susceptibility of a Au(27.5nm)/Ho(4.5nm)/Nb(150nm) trilayer by low-energy muon spin spectroscopy (LE-µSR). The antiferromagnetic rare earth metal Ho breaks time-reversal symmetry of the pair correlations in Au and has a thickness that is comparable to the known coherence length for singlet pairs in Ho [18] to ensure pair transmission into Au. The Au layer is necessary since a Meissner state cannot be probed by muons directly in a magnetic material due to their rapid depolarization in a strong magnetic field. Here, we report the discovery of the paramagnetic Meissner effect in Au, which is found to be an intrinsic property of the odd-frequency superconducting state that is generated via the superconductor proximity effect.

II. RESULTS
LE-µSR offers extreme sensitivity to magnetic fluctuations and spontaneous fields of less than 0.1 Gauss with a depth-resolved sensitivity of a few nanometers [19][20][21][22][23]. To probe the depth dependence [z coordinate in Fig. 1(a)] of the Meissner response in Au/Ho/Nb by LE-µSR, an external field (Bext) is applied parallel to the sample plane [along the y coordinate in Fig. 1(a)] and perpendicular to the muon initial spin polarization (oriented in the x-z plane) as sketched in Fig. 1(a). The muon spin polarization is proportional to the asymmetry of decay positrons from the implanted muons as shown in Fig. 1(b), which is experimentally determined as a difference in the number of counting events of the two detectors, as discussed in the Supplemental Material [24].
In this transverse-field (TF) configuration, a muon's spin polarization precesses on average at a frequency ̅ = ̅ about the average local field ̅ sensed by the implanted muons, with γµ = 2π * 135.5 MHz T -1 being the muon gyromagnetic ratio. Assuming a local field profile Bloc(z) within the sample, muons implanted with energy E and a corresponding stopping distribution p(z,E), precess at an average frequency ̅ = ̅ = ∫ ( ) ( , ) [25]. The asymmetry spectrum As(t,E) can be approximated as ∝ − ̅ cos[γµ ̅ t+φ0(E)] for each implantation energy E ( ̅ and φ0(E) being the mean depolarization rate and starting phase of the muon precession, respectively). The experimental Bloc(z) profile is therefore sampled as series of mean field values ̅ of the magnetic field distributions p(Bloc) [Fig. 1(c)] determined as fits of the corresponding asymmetry functions As(t,E) measured at different energies E [24].
To investigate the paramagnetic Meissner effect, implantation energies in the 3-6 keV range were used to determine the Bloc(z) profile in the Au layer. At the lowest energy of 3 keV, the muons contributing to the asymmetry stop within the Au, while for increasing energy, an increasing fraction stops within the Ho and Nb layers [ Fig. 1(a) and Fig. S1 in the Supplemental Material [24]]. The implantation profiles were calculated using the Monte Carlo algorithm TrimSP [26]. To minimize the contribution from backscattered muons [26] to the measured signal, implantation energies below 3 keV were not used. For muon energies above ~7 keV, the contribution of the Nb becomes dominant and therefore not relevant for probing the Meissner state in the Au (Fig. S1 [24]). However, energies above 7 keV are important to confirm the emergence of a conventional (diamagnetic) Meissner response in Nb in the superconducting state. Muons stopping in the Ho layer depolarize almost immediately and do not contribute to the measured asymmetry.  Fig. S1(b) of the Supplemental Material [24]. Since the two energy scans at 3 K and 5 K were performed at different stages, the normal-state (10 K) energy scan is reacquired to avoid any influence on the measurement data from the specific magnetic configuration reached by Ho after cooling through its magnetic transitions. to the absence of a magnetic interface [27]. Although the measured increase in the local field, ΔBloc(z), in the superconducting state relative to the normal state is small, it exceeds the statistical and systematic measurement error in Bloc(z) (error bars in Fig. 2). Furthermore, ΔBloc(z) increases at lower temperatures, which implies that the magnitude of ΔBloc(z) is related to the amplitude of the superconducting order parameter which is consistent with theory [7,8].
Although the asymmetry fits for single implantation energy in Fig. 2 show a paramagnetic Meissner effect in Au below the Nb superconducting transition, the Bloc(z) profiles obtained with this approach include depth averaging due to the width of muon stopping distributions. To obtain an accurate Bloc(z) profile, a global fit for all implantation energies with a common field profile is used [19][20][21][22][23]. The common field profile in the Au layer is modelled as Bloc(z) = Bext + M(z), where we set the magnetization term to M(z) = Ba sin(z/κ), which is a parameterization of the theoretical magnetization profile calculated for the Au/Ho/Nb heterostructure.
The theoretical magnetization is computed from the vector potential A determined as a solution of the Maxwell equation where the supercurrent J is assumed proportional to the vector potential A via the term Jx(z) including the dependence on the anomalous Green's function (see Supplemental Material [24]). In this expression, Jx(z) also represents the component of the supercurrent density J along the x-axis in Fig. 1(a). Both odd-frequency and even-frequency pairing correlations contribute to J, which is calculated using the quasiclassical theory of superconductivity under the assumption that time-reversal symmetry is spontaneously broken by the spatially-dependent exchange field of the Ho which forms a conical pattern along the z coordinate in Fig. 1(a). We also take into account the spin-selective scattering taking place at the interface between Nb and Ho by using spin-dependent boundary conditions [24]. Our model excludes the presence of a Fulde-Ferrel-Larkin-Ovchinnokov (FFLO) state which can theoretically compete with the paramagnetic Meissner state, but only if the superconducting layer is thinner than the magnetic screening length [6]. In Nb the magnetic screening length is about 90 nm, which is much shorter than the thickness of the Nb used here of 150 nm, and so contributions from the FFLO state can be ignored, as stated in Ref. [6]. In the particular case of Au/Ho/Nb, using realistic values for the physical parameters involved in the description of the proximity effect occurring in Au and Ho, the expected theoretical profile for Bloc(z) shows a single oscillation reaching a maximum inside Au (blue curve in Fig. 3). Therefore, making also the realistic assumption that Bloc(z) matches the applied external field Bext at the Au/vacuum interface, it is clear that Bloc(z) = Bext + Ba sin(z/κ) represents an appropriate parameterization of the oscillatory local magnetic field profile in Au to use in the global energy fit. This parametrization is also in agreement with the experimental profiles determined at 3 and 5 K by sampling Bloc(z) at different energies, which can be approximated by half-period sine functions (blue curves for Au in Fig. 2). The global energy fit was implemented on the measurement data at 3 K which show the most significant paramagnetic Meissner response in Au. An exponential depolarization function G(t,E) = − ̅ was used for the fit, with ̅ being a fitting parameter common for all energies [24]. In the analytic expression for Bloc(z), Bext was set equal to the normal-state field obtained from the global fit of the measurement data at 10 K under the assumption that Bloc(z) can be modelled as a constant field at this temperature [magenta curve in

III. DISCUSSION
While the fits to the data are in close agreement with the predicted paramagnetic Meissner effect, we first rule out other possibilities that could lead to an increase in magnetic flux in Au. However, as we discuss below, none of these should show temperature dependence in the measured range of 3-10 K. One mechanism that can result in a magnetization enhancement in Au is Ruderman-Kittel-Kasuya-Yosida (RKKY)-type oscillations in the spin-polarization of Au induced via an interaction with Ho. The largest period predicted [28] and reported experimentally [29] for RKKY oscillations in epitaxial Au (001) is 7-8 monolayers or ~ 1.6 nm, which is much too small to explain the behaviour observed in Fig. 2, where the oscillation period of the magnetization exceeds several tens of nanometers. Furthermore, RKKY oscillations should also lead to an additional broadening of the magnetic field distribution experienced by muons (other than to the measured shift in average field) [30] and also be present in the normal-state data at 10 K, which we do not observe. A second possibility is an enhancement of the magnetization of Ho for decreasing temperature. Similarly to the case of RKKY oscillations, however, such an enhancement should only result in a broadening of the magnetic field distribution rather than the observed shift because the magnetic domains have a finite size and a random orientation giving a random dipolar field profile in Au.
Oscillations in the magnetic susceptibility induced by unconventional (odd-frequency) superconductivity are therefore the most likely explanation for the paramagnetic Meissner effect in Au due to the presence of Ho [6][7][8][9]. In Fig. 2 it is shown that a conventional Meissner effect is measured in Nb up to the interface with Ho, where the contribution of spin-singlet Cooper pairs to the screening supercurrent Jx(z) is larger than that due to the long-ranged spin-triplet pairs (Fig. 3). Nevertheless, while spin-singlet pairs are rapidly filtered out by the exchange field in Ho, the spin-triplet pairs rotate into each other within the Ho layer Finally, the results in Fig. 3 demonstrate that the paramagnetic effect is most strong in Au where the odd frequency state dominates over the singlet state. This indicates that the paramagnetic response is an intrinsic property of the odd frequency superconducting state and that the superconductivity must therefore carry a net magnetization. Future experiments should explore ways to harness the magnetization generated by odd frequency superconductivity in order to explore the potential for driving magnetization-reversal process in the superconducting state [13].

I. SAMPLES PREPARATION AND LE-μSR MEASUREMENT SETUP
The Au(27.5nm)/Ho(4.5nm)/Nb(150 nm) thin film multilayers were grown onto unheated oxidised Si substrates by DC magnetron sputtering in Ar plasma at 1.5 Pa. Before and during the film growth, the walls of the deposition chamber were cooled via a liquid nitrogen jacket to achieve a base pressure of below 10 -8 Pa (verified using an in-situ residual gas analyser). The substrates were placed on a circular table and rotated below stationary sputtering targets. Deposition rates were pre-determined for each material using an atomic force microscope to measure pre-deposited step edges on calibration samples. The calibrated deposition rates were used to set the rotation speeds and deposition powers needed to achieve the desired thickness for each metallic layer.
In contrast to bulk µSR where ~4 MeV muon beams are used, in the LE-µSR apparatus at Paul Scherrer Institute (PSI) lower muon implantation energies can be obtained using a few-hundred-nanometer-thick Ar solid gas moderator capped with ~ 10 nm N2 and grown on top of a silver foil [31] (about 100 µm in thickness).
Choosing appropriate beam transport and sample voltages, the muon implantation energy can be varied between 0.5 keV and 30 keV, which allows tuning of the mean muon stopping depth in the range 1-200 nm with an accuracy of a few nanometers. The energies needed to probe the local magnetic field at specific depths inside the Au/Ho/Nb samples were chosen on the basis of the stopping profiles calculated using a Monte Carlo algorithm TRIM.SP [32]. A transverse-field geometry [33] was adopted with the external field Bext ~ 101 Gauss applied in the film plane and perpendicular to the initial polarization of the muons in order to bring the sample into a Meissner state. This configuration, usually adopted when is larger than any internal field, is extremely sensitive to sample magnetic inhomogeneities [34], which would result in a spreading of the Larmor frequencies for the individual muon decay events. LE-µSR measurements were carried out at two different temperatures (3 K and 5 K) below the Nb superconducting transition, with 3 K being the lowest temperature achievable with the measurement apparatus at PSI.

II. THEORY AND FITTING OF THE LE-μSR MEASUREMENTS
A schematic diagram of the apparatus used to probe the local magnetic field of our Au/Ho/Nb heterostructure by LE-µSR has been reported in Fig. 1(a). In a transverse-field (TF) configuration, where the initial polarization of the muon spin is perpendicular to the applied field, the starting point for the analysis of the LE- Significant improvement to the fit can be made by taking the contribution of the muon stopping profile p(E,z) and the magnetic field profile into consideration when calculating the asymmetry signal. In order to perform this improved fit, the following more general expression is used for As(t,E) instead of equation ( where the integral is now extended over the entire sample depth range probed by the muon beam at a given implantation energy E. The local magnetic field Bloc(z) is not treated as a constant field for a given muon implantation energy like in (1.4), and the goal is to determine a functional form for Bloc(z) = Bext + M(z) which is consistent with the measurement data simultaneously for all the energies considered. A good convergence of the chi-square minimization algorithm run by the software musrfit [35] (chi-square/number of degrees of freedom ratio close to 1) for all the energies is normally only achieved for an appropriate choice of the functional form for Bloc(z).
Since the theoretical profile calculated for Bloc(z) in Au can be parameterized by a half-cycle sinusoid with Bloc(z) = at the Au/vacuum interface (blue curve in Fig. 3), we have used the following functional form for Bloc(z) in (1.5) to perform the global energy fit: where the amplitude Ba and the angular frequency κ -1 of the sinusoid are the model parameters. This common field profile is also consistent with the trend of the experimental ̅ (z) values determined from single-energy asymmetry fits in the superconducting state at 3 K and 5 K (blue curves in Fig. 2).
Since the expression (1.6) does not apply to the local field profile in the Nb layer, only energies up to 6 keV were taken into account to implement the global fit. This is consistent with the simulated muon stopping fractions in Fig. S1(c), which show that the contribution of muons stopping in Nb to the asymmetry signal becomes non-negligible at energies higher than 6 keV. Bext in (1.6) was determined from the global fit of the measurement data in the normal state at 10 K assuming that Bloc(z) = constant = Bext through the entire Au/Ho/Nb multilayer.
Using ̅ and as fitting parameters and fixing κ = 13.58 nmwhich corresponds to a sine function having the same peak position inside Au as the theoretical ( ) profile (blue curve in Fig. 3

III. ANALYTICAL DESCRIPTION OF THE ANOMALOUS MEISSNER RESPONSE
Normalizing the z-axis coordinate by L=dN+dF, where dN(dF) is the thickness of the N(F) layer, so that z=0 corresponds to the Nb/Ho interface and z=1 to the Au/vacuum interface, the Maxwell equation that determines the magnetization response reads: where A is the vector potential and Jx(z) is the normalized screening current density along the x-axis in Fig. 1.
Using the linear response theory formalism, the screening supercurrent density J of Maxwell equation has been written in equation (1.7) as proportional to the vector potential A via a factor Jx(z) containing the dependence on the anomalous Green's functions [36].
The induced magnetization M in the N layer and the vector potential A are related through the following expression: Here 2) is normally difficult to solve analytically [37].
Nevertheless, to understand the basic physics underlying the anomalous Meissner response using this equation, it is possible to consider some limiting cases.
Firstly, it can be proven [38] that the sign of the screening current in (1. For our Au/Ho/Nb heterostructure, we expect that even-and odd-frequency pairings are mixed together in Au. Therefore, any signature of an anomalous Meissner effect in N implies that the odd-frequency component is dominant herein [38].

IV. SOLUTION OF THE USADEL AND MAXWELL EQUATIONS
To quantify the proximity effect that generates superconducting correlations in the Ho and Au layers, we use the quasiclassical theory of superconductivity. The correlations are then described by a Green's function matrix g that satisfies the Usadel equation in the non-superconducting region: Here, D is the diffusion constant, is the quasiparticle energy, 3 =diag(1,1,-1,-1), while R(z) is a 4x4 matrix that contains the exchange field profile {for further details see [39]}. To compute g, one also needs to use boundary conditions at the interfaces of the junction. At the Au/vacuum interface, the derivative of the Green's function must be zero, while at the Nb/Ho interface we use the standard Kupryianov-Lukichev boundary conditions valid for a non-ideal interface (as this is experimentally realistic): where ς is the ratio between the interface resistance and the resistance of the non-superconducting layer, L = dN + dF is the length of the non-superconducting layer, gBCS is the bulk superconducting Green's functions, Gϕ is a phenomenological parameter capturing the spin-dependent phase-shifts occurring at the interface, while β is a matrix describing the orientation of the interface magnetic moment [39].
In the superconductor, we use the bulk superconducting Green's function since the superconductor is much larger than the other regions and it acts as a reservoir. To model the Ho/Au region, we use a spatiallydependent exchange field with a magnitude constant in Ho (rotating in direction) but gradually dropping to zero over a region of few nanometers centred at the Ho/Au interface. As long as this drop is not too sharp (i.e. occurring over 1 nm or less), this assumption is allowed by the quasiclassical formalism, which demands that all length-scales involved must be much larger than the Fermi wavelength (typically of the order of 1 Å).
Within the linear response theory, the amplitude J of the shielding current density J generated in the junction reads: where N0 is the normal-state density of states at the Fermi level, e is the electron charge, ̆ is the 8x8 Green The amplitudes of the anomalous Green's functions have been derived solving the exact Usadel equation in the non-linearized regime (absence of a weak proximity effect). Each anomalous Green's function amplitude is defined as sum over a fine mesh of energies between zero and ten times the superconducting gap. The quantization axis is taken to be along the interface magnetic moment. This figure is meant to illustrate the different behaviour of the spatial evolution of the triplet and singlet components by defining |f| 2 as a measure of the strength of the proximity effect. We emphasize that even if the total triplet amplitude defined in this way is larger than the singlet even near the Nb interface, this does not mean that the triplet contribution to the actual screening current necessarily is larger in this region, since this quantity depends on the anomalous Green's functions in a different way than |f| 2 .