Introduction

Interfaces with inversion symmetry breaking present new avenues for electron pairing, providing a valuable platform for exploring the emergent quantum phenomena associated with unconventional superconductivity1,2,3. The most prominent example is the superconducting interface between two insulating oxides LaAlO3 (LAO) and SrTiO3 (STO)4. In addition to the two-dimensional (2D) superconductivity, the coexisted ferromagnetism at this interface has also been unveiled by high-resolution magnetic torque magnetometry5, scanning superconducting quantum interference device (SQUID)6, and magnetoresistance7. The magnetic exchange is considered to be a key ingredient of high-temperature superconductivity, and the formation of unconventional superconducting pairs is closely related to magnetism8,9,10. Obviously, 2D oxide heterointerfaces provide opportunities for exploring the underlying physics of unconventional superconductivity. However, the extremely low superconducting critical temperature TC (below 250 mK) of the STO-based 2D electron gases (2DEGs)1,11,12 limits a thorough investigation of the superconductivity-magnetism correlation.

The recent discovery of superconducting 2DEGs at the KTaO3 (KTO)-based heterointerfaces has attracted much attention because its TC is up to ~2 K13,14,15,16,17, nearly an order of magnitude higher than that of STO. KTO is in many ways similar to STO. However, the 2DEGs residing in KTO derive from 5d orbitals rather than 3d orbitals, exhibiting a large atomic spin-orbit coupling (SOC) and unique properties compared to the 2DEGs in STO18. For example, the superconductivity of KTO-based 2DEGs is strongly dependent on crystal orientation15,16,17,19,20,21, a feature that is absent in the STO system1,11,12. Although the 2DEGs at KTO interfaces with (001) crystal orientation exhibit no interfacial superconductivity down to 25 mK15,16, higher transition temperatures up to 0.9 and 2.2 K were observed for samples oriented along (110)20 and (111)17, respectively. These results indicate that the superconductivity could be closely related to the interfacial d-orbital reconstruction22. Meanwhile, an admixture of s-wave and p-wave pairing components induced by strong SOC has been revealed in KTO-based 2DEG14. In principle, the integration of strong SOC, superconductivity and long-range ferromagnetic (FM) order could potentially realize the long-sought-after spin-triplet pairing and topological superconductivity23,24,25,26,27,28. We note that the KTO interface has already demonstrated strong SOC and 2D superconductivity, thus it is of great significance to explore the long-range ferromagnetism that coexists with superconductivity. Very recently, magnetization hysteresis loops were observed for the (111)-AM/KTO interfaces over a wide temperature range from 1.8 to 300 K (AM denotes LaAlO3 or YAlO3 amorphous films)29, suggesting coexisted superconductivity and ferromagnetism noting that 1.8 K is slightly lower than the onset temperature of the superconductivity investigated there. However, direct magnetotransport evidences are lacking, especially when the system is in the well-superconducting state at extremely low temperatures. Here, we present magnetotransport evidences for the coexistence of 2D superconductivity and 2D ferromagnetism over a wide temperature range for (111)-oriented a-CZO/KTO interfaces (CZO=CaZrO3 and a-CZO denotes amorphous CZO film). As experimentally demonstrated, butterfly-shaped magnetoresistance (MR) emerges when temperature is lower than the superconducting onset temperature, achieves a peak value around TC, and decreases rapidly with the further decrease of temperature. However, MR remains sizable even at the temperature of 0.05 K which is well below the zero-resistance temperature, indicating the concomitant occurrence of magnetic ordering and superconducting. A rapid increase of the hysteretic MR with the sweep rate of the magnetic field is further observed while fixing the temperature at a constant value, providing a conclusive proof that the magnetization dynamics are at play in the superconducting state. Density functional theory (DFT) calculations show that oxygen vacancy-induced gap states cause an increase in the density of states (DOS) near the Fermi surface, inducing the itinerant FM state according to the Stoner model.

Results and discussion

Amorphous CZO films, 10 nm in thickness, were grown on (111)-oriented KTO single crystal substrates by pulsed laser deposition (PLD). The surface morphology of the film is atomically flat (Fig. S1). The microstructure of the heterointerface was assessed by the technique of high-resolution scanning transmission electron microscopy (STEM). The typical high-angle annular dark-field (HAADF) image of the a-CZO/KTO(111) cross-section, recorded along the [1-10] zone axis of KTO is shown in Fig. 1a. A large-scale HAADF-STEM image is shown in Fig. S2. Structural characterization indicates that the CZO film is amorphous and homogeneous, and the a-CZO/KTO interface is abrupt and smooth. A close inspection of the lattice image indicates a slight blurring in color contrast within an interfacial layer of ~0.5 nm at the a-CZO/KTO interface, which could be attributed to the deficiency of oxygen near the interface. Additionally, X-ray photoelectron spectroscopy (XPS) measurements indicate the presence of Ta4+ ions with a content of 22.8% (see Fig. S3 and Table S1), suggesting the occurrence of oxygen vacancies which could be introduced into the KTO substrate by the PLD deposition of the a-CZO film. Notably, oxygen vacancies act as electron donors, contributing charge carriers to the 2DEGs at the a-CZO/KTO interfaces.

Fig. 1: Two-dimensional superconductivity at a-CZO/KTO(111) interface.
figure 1

a HAADF-STEM image of a-CZO/KTO(111) heterostructure measured along [1-10] zone axis, showing that the CZO film is amorphous. b RS/R4K as a function of temperature for samples with different carrier densities nS, where R4K is the normal-state resistance at 4 K. The direction of the arrow indicates the increase of nS, which grows from 3.8 × 1013 to 9.8 × 1013 cm−2. c nS dependence of TC, where TC is determined at 50% of R4K. The dashed line is a guide for the eyes. d Temperature-dependent I–V measurements for sample with nS = 3.8 × 1013 cm−2. e I–V curves plotted on a logarithmic scale with the same color codes as in (d). The Red dashed line represents V∝I3, which is used to infer the BKT transition temperature. f Temperature dependence of the power-law exponent α, as deduced from the linear fits of the curves in (e). g Sheet resistance dependence of temperature plotted on a [dln(RS)/dT]−2/3 scale. The red solid line is a linear extrapolation from the high-T linear section, which crosses the T-axis at TBKT = 0.706 K.

Full size image

By adjusting the growth parameters, metallic 2DEGs with different carrier densities (nS) were obtained, as detailed in the Methods section and Table S2. The carrier density increases as the oxygen pressure for PLD decreases, lateral evidence for the existence of oxygen vacancies. Figure 1b shows the temperature dependence of the reduced sheet resistance RS/R4K (R4K is the normal-state resistance at 4 K) for the sample with a nS between 3.8 × 1013 and 9.8 × 1013 cm−2, where nS is deduced from the Hall effect of 10 K (Fig. S4). The superconductivity can be clearly seen for all samples, and the superconducting transition temperature TC, defined by the temperature corresponding to R4K/2, shows a strong dependence on carrier density. It exhibits a linear increase with nS with a TC–nS slope of ~2.84 × 10−14 K∙cm2 (Fig. 1c). A similar relation between TC and nS was also found for the 2DEGs at the EuO/KTO(111)15 and AlOx/KTO(111)13 interfaces, displaying a TC–nS slope reduced by ~20%. Remarkably, TC is as high as ~2.4 K in our 2DEG with nS = 9.8 × 1013 cm−2, nearly 12 times higher than that in the LAO/STO system1.

The Berezinskii–Kosterlitz–Thouless (BKT) transition is a typical behavior of the 2D superconducting systems, characterized by a V∝Iα power-law relation with α = 3 at the transition temperature1,30, where I and V are applied current and corresponding voltage, respectively. To get an idea about the BKT transition, Fig. 1d presents the I–V characteristics acquired at different temperatures for the sample with nS = 3.8 × 1013 cm−2. By linearly fitting the logarithmically scaled I–V relation (Fig. 1e, f), we found a TBKT of ~0.618 K, corresponding to the exponent α of 3. In addition, the TBKT can also be estimated from the formula RS(T) = R0exp[−b(T/TBKT−1)−1/2], where R0 and b are material parameters31. Application of such a fit to the RS(T) curve suggests TBKT = 0.706 K for the sample of nS = 3.8 × 1013 cm−2 (Fig. 1g), consistent with the corresponding result obtained from the analysis of the I–V characteristics. A similar analysis for the sample with nS = 9.8 × 1013 cm−2 is shown in Fig. S5, and TBKT of ~2.153 K and 2.078 K are obtained from the analysis of the α exponent and the fit to the RS(T) curve, respectively. These results confirm the 2D nature of the superconductivity at the a-CZO/KTO(111) interfaces.

To reveal the effect of magnetic field on superconductivity, we measured the reduced resistances RS/R0.5T and RS/R8T with out-of-plane field (μ0H⊥) and in-plane magnetic field (μ0H||), respectively, for the sample of nS = 3.8 × 1013 cm−2 at different temperatures, and the corresponding results are shown in Fig. 2a, b. From the data in Fig. 2a, b, the upper critical field can be deduced. Figure 2c illustrates the temperature dependence of the upper critical fields μ0HC2⊥and μ0HC2||, defined by the half normal-state resistances R0.5T and R8T, respectively. We further performed an analysis of the μ0HC2–T relation based on the Ginzburg–Landau theory32:

$${\mu }_{0}{H}_{{\mbox{C}}2\perp }=[{\varphi }_{0}/2\pi {\xi }_{{\mbox{GL}}}^{2}(0)]\left(1-\frac{T}{{T}_{{\mbox{C}}}}\right)$$
(1)
$${\mu }_{0}{H}_{{\mbox{C}}2{\mbox{||}}}=[{\varphi }_{0}\sqrt{12}/2\pi {\xi }_{{\mbox{GL}}}(0){d}_{{\mbox{SC}}}]{\left(1-\frac{T}{{T}_{{\mbox{C}}}}\right)}^{1/2}$$
(2)

where \({\varphi }_{0}\) is the magnetic flux quantum, and dSC is the superconducting layer thickness. Fitting the experiment data to Eqs. (1) and (2) leads to the zero-temperature limit μ0HC2⊥(0) = 0.09 T and μ0HC2||(0) = 3.99 T, indicating a large anisotropic ratio of ~44, and the coherence length ξGL(0) = 61.63 nm more than 13 times larger than dSC = 4.63 nm. These results reconfirm the 2D characteristics of the superconductivity. In addition, the mean free path lmfp of the conducting electrons can be estimated using lmfp = (h/e2)(1/kFRS), where kF = (2πnS)1/2 and h are the Fermi wave number and Planck constant, respectively30. Adopting the nS and RS obtained at T = 10 K, the lmfp of the present 2DEG is found to be ~12.3 nm, which is about 20% of ξGL. Therefore, the superconductivity at the a-CZO/KTO(111) interface occurs in an intermediate range between clean (lmfp ≫ ξGL) and dirty (lmfp ≪ ξGL) limits, which is similar to a-LAO/KTO(110)20 and EuO/KTO(111)15, but different from EuO/KTO(110) for which lmfp > ξGL33. For weakly coupled BCS superconductors, the Pauli paramagnetic limit acts as the theoretical upper bound for the parallel critical field34,35, which is given by \({\mu }_{0}{H}_{{\mbox{C}}2}^{P}\approx 1.76{k}_{{\mbox{B}}}{T}_{{\mbox{C}}}/\sqrt{2}{\mu }_{{\mbox{B}}}\), where μB is Bohr magneton. Taking TC = 0.668 K, we obtain \({\mu }_{0}{H}_{{\mbox{C}}2}^{P}=1.24\,{\mbox{T}}\) (a red solid circle in Fig. 2c), which is only 31% of the measured μ0HC2||(0) = 3.99 T. Obviously, the Pauli paramagnetic limit is drastically exceeded. To clarify the origins of the large μ0HC2||, we systematically investigated the back gating effect for the a-CZO/KTO(111) sample with nS = 3.8 × 1013 cm−2, which corresponds to the nS value of the as-grown sample (Note 1 and Fig. S6). The carrier density nS can be tuned by applying the gate voltage VG. The dependence of μ0HC2|| on nS is obtained, indicating that μ0HC2|| increases from 3.6 to 5.0 T as nS decreases from 3.8 × 1013 cm−2 to 2.9 × 1013 cm−2 (Fig. S6g). Note that a higher μ0HC2|| is observed in the a-CZO/KTO(111) system compared to LAO/STO (111) systems (~4.5 T)40. The analysis of the nS-dependent spin-orbit energy εso reveals that εso increases as nS decreases (Fig. S6i), which is consistent with the variation of μ0HC2|| with nS. Therefore, the high μ0HC2|| value exceeding the Pauli limit could be attributed to the strong SOC of the 5d 2DEG at KTO-based interfaces, which stems from the inversion symmetry breaking at the a-CZO/KTO interface and the presence of relatively heavy tantalum ions18,36,37,38,39.

Fig. 2: MR loops in superconducting state.
figure 2

a RS/R0.5T as a function of the out-of-plane magnetic field at different temperatures ranging from 0.03 to 2 K, where R0.5T is the normal-state resistance when μ0H⊥ = 0.5 T. b In-plane magnetic field dependence of RS/R8T, where R8T is the normal-state resistance when μ0H|| = 8 T. The symbol label is the same as that in (a). c Temperature dependence of the upper critical field μ0HC2, derived from half of the normal-state resistance. Dashed lines are fitting curves based on Ginzburg–Landau theory. The estimated Pauli paramagnetic limit is marked with a red dot. d RS/R8T as a function of μ0H|| at different temperatures. Blue and red arrows indicate the direction of the field sweep. The upper left corner shows the geometry for the measurements. e In-plane magnetoresistance at 0.05 K with the magnetic field parallel or perpendicular to applied current. The upper right corners indicate the measurement geometries. f Temperature dependence of the difference in sheet resistance, denoted as ∆RS = RS(μ0H|| = 0.08 T)−RS(μ0H|| = 0). The marks represent the experimental data. The solid line is fitted using the LAMH theory for TAPS.

Full size image

Strikingly, sharp resistance peaks are observed in the low-field range of the RS/R8T–μ0H|| curves collected at low temperatures (Fig. 2b), suggesting a deviation of the transport behavior from superconductivity. The resistance peaks can be clearly seen in a magnified view of the low-field data (Fig. S7). To get further knowledge about this resistive anomaly, the resistance is re-measured by sweeping μ0H|| in a narrow field range between −1.5 and 1.5 T. Figure 2d shows the reduced resistance recorded with the μ0H|| parallel to applied current I at several typical temperatures (the RS/R8T–μ0H|| curves have been vertically shifted for clarity). The most remarkable observation is the magnetic hysteresis in the superconducting background, as indicated by the appearance of butterfly-shaped RS/R8T–μ0H|| curves. This is a signature of magnetic ordering as will be discussed in detail later. Further analysis shows that the resistive anomaly is robust. Hysteretic MR loops occur no matter μ0H|| is parallel or perpendicular to the applied current (Fig. 2e). Applying current along [11-2] or [1-10] crystal orientations also lead to magnetic hysteresis (Fig. S8). It is therefore a general feature of the magnetotransport behaviors of the a-CZO/KTO interface. Moreover, MR peaks appear over a wide temperature range, covering both the superconducting state (T < TC) and the normal state (T > TC) (Fig. S7). The magnetism in the normal state can also be confirmed by the magnetization dependence on the magnetic field, measured using the SQUID. The ferromagnetic hysteresis loops at various temperatures clearly demonstrate the persistence of normal-state magnetism up to room temperature (Fig. S9), consistent with the findings reported in ref. 29. The undetectability of the MR loop at higher temperatures could be attributed to the precision limitations of the magnetotransport system.

To get a quantitative description about the thermal evolution of the butterfly-shaped RS/R8T–μ0H|| behavior, in Fig. 2f we show the peak resistance, defined by ∆RS = RS(μ0H|| = 0.08 T)-RS(μ0H|| = 0), as a function of temperature. Here μ0H|| = 0.08 T is the position of resistance peaks. As shown in Fig. 2f, ∆RS is low at high temperatures, rapidly increases with the decrease of temperature, and gets the maximum value ~ 814.9 ohms/square (Ω/□) at ~ 0.5 K, a temperature slightly lower than the zero-resistance temperature (~0.6 K). Further cooling results in a drastic decrease in ∆RS. However, ∆RS remains identifiable even at the temperatures well below TC. It is, for example, ~9.7 Ω/□ at 0.05 K where the system is expected to be in the well-superconducting state. The implications of these results are twofold. The first one is the concomitant occurrence of the FM phase and superconducting phase, as suggested by the considerable MR peaks at temperatures well below TC, and the second one is the strong competition between magnetic ordering and superconducting, as implied by the occurrence of the strongest MR effect close to the zero-resistance temperature at which a percolation of the superconducting phase may just take place. An explanation for this peak behavior of ∆RS will be presented later.

Magnetic hysteresis MR has been reported for the superconducting 2DEG at the LAO/STO interface, and ascribed to the effect of vortex moving41. It was also observed in normal-state 2DEG at the EuO/KTO interface33, and attributed to the magnetic proximity effect of the FM EuO over layer on 2DEG42,43. However, this phenomenon was never observed before in the superconducting 2DEG at the KTO-based interfaces.

In general, the occurrence of magnetic hysteresis is a signature of the existence of magnetic domains. As well documented42, a resistance peak emerges when magnetic domains change their polarities as the magnetic field sweeps through coercive field; the structure of the magnetic domains is most disordered at the coercive field and, correspondingly, the resistance is maximal due to enhanced magnetic scattering. It is a typical feature of metallic magnetic materials that MR peaks at coercive field.

The magnetotransport process in superconductors may be somewhat different from that in normal conductors. When the system enters the superconducting state, the onset of superconducting pairing could induce a granular superconductor44. Given the impact of the domain wall of the ferromagnet on superconductivity, the magnetization dynamics in the ferromagnet will generate moving vortices in the superconductor as μ0H|| sweeps41. Notably, the intermediate areas between adjacent superconducting grains can enclose magnetic vortices. In this case, the motion of the magnetic vortices will cause an increase in resistance as the vortices cross the weak links between superconducting grains41,44. When the system is cooled well below TC, the superconducting islands expand and merge with each other, forming percolated superconducting areas that dominate the transport behavior of the 2DEGs. As the intermediate region between the interconnected islands shrinks, the resistance effect induced by the vortex moving across these junctions is quickly weakened, which is in good agreement with the reduction in the MR peak amplitude with decreasing temperature. Note that a finite resistance peak observed at T = 0.05 K clearly verifies that the dissipationless energy transfer in superconductors has been disrupted by vortex motion.

Obviously, the proportion of the weakly linked area between superconducting grains determines the size of the hysteretic MR peaks. With the decrease of this proportion, the MR peak is expected to diminish and, finally, vanish when superconducting areas dominate the transport behavior of the 2DEG. Considering the close relation between superconductivity and carrier density shown in Fig. 1b, c, we further investigated the influence of nS on the RS–μ0H|| loops measured at 0.05 K by applying a magnetic field parallel to current (Figs. 3a and S10). The dependence of the peak resistance ΔRS on carrier density is shown in Fig. S11. As the doping level grows from 3.8 × 1013 to 4.3 × 1013 cm−2, the peak value of the resistance increases from ~9.7 to ~19.6 Ω/□. However, the resistance peak vanishes when nS exceeds 5 × 1013 cm−2, leaving a zero-resistance background. To explore the role of disorder, the mean free path lmfp is plotted as a function of nS (Fig. S11). The results indicate that lmfp does not exhibit a significant dependence on nS. Obviously, the disorder mechanism cannot account for the disappearance of the resistivity peak at densities above 5 × 1013 cm−2. The butterfly-shaped hysteresis could be attributed to the vortex motion across the weak links between superconducting grains. At low nS, superconducting islands are sparsely distributed in a non-superconducting background. In this case, the small number of interconnected islands results in a weak resistance effect when the vortex moves across these junctions, consistent with the small magnitude of ΔRS. As nS increases, the sample enters into a deeper superconducting state due to the higher TC, which means that the superconducting islands expand and form more interconnected islands. The increase in the proportion of the weakly linked area leads to an increase in ΔRS. However, when nS further increases, crossing the percolation threshold, the continued growth and merging of superconducting islands cause the weakly linked regions to be short-circuited by the superconducting phase, resulting in a dropping of ΔRS to zero. Thus, the disappearance of the resistivity peak at carrier densities above 5 × 1013 cm-2 could be attributed to the percolation of the superconducting phases.

Fig. 3: Hysteresis induced by the magnetization dynamics.
figure 3

a In-plane magnetoresistance for samples with different carrier density at 0.05 K. b In-plane magnetoresistance at different VG values measured at 0.05 K for the sample with nS = 3.8 × 1013 cm−2, which corresponds to the nS value of the as-grown sample. c In-plane magnetoresistance at different field sweep rates for nS = 4.3 × 1013 cm−2 at 0.05 K. d RS/R4K as a function of normalized temperature T–TC. The marks represent the experimental data. The solid lines and dashed black lines are fitting curves using the LAMH theory for TAPS and the Golubev–Zaikin model for QPS, respectively. The hollow patterns indicate the values of carrier densities.

Full size image

As mentioned earlier, the back gating effect for the a-CZO/KTO(111) sample with nS = 3.8 × 1013 cm−2, corresponding to the nS value of the as-grown sample, has been investigated (Fig. S6). Figure 3b plots the peak resistance ΔRS as a function of the back gate voltage VG. Notably, as VG changes from +150 V to −100 V, ΔRS decreases from 13.4 to 9.8 Ω/□, and lmfp decreases from 13.7 to 10.2 nm (Fig. S6c). The application of a negative VG will push charge carriers towards the a-CZO/KTO(111) interface, leading to a higher degree of disorder. This reduction in ΔRS cannot be explained by enhanced disorder, as the higher disorder induced by negative VG generally results in increased resistance. The decrease in ΔRS as VG decreases could be attributed to reduced carrier density (Fig. S6c), which results in more sparsely distributed superconducting islands in a non-superconducting background. Consequently, the resistance effect induced by vortex motion across the weak links between superconducting grains is weakened.

To get the knowledge about the dynamic process, the magnetic hysteresis has been systemically explored by varying the sweep rate of μ0H|| for the 2DEG with nS = 4.3 × 1013 cm−2. As shown in Fig. 3c, the magnitude of the MR curve strongly depends on the field sweep rate, and rapidly increases with the increase of the sweep rate. As the sweep rate increases from 0.02 to 0.10 T/min, for example, the peak value grows from ~0.4 to ~19.6 Ω/□. At very slow sweep rates, the peak becomes minor and even negligible. Similar results have been reported for the superconducting LAO/STO interfaces, as the transport evidence for the coexistence of superconductivity and ferromagnetism41,45. As expected from our earlier discussion, the vortex motion generated by the dynamics of the ferromagnet during magnetization reversal leads to a sharp resistance peak. The rate at which the magnetic vortices are generated is proportional to the time dependence of the magnetization of the ferromagnet, and the rate of magnetization change is proportional to the sweep rate of external field. Therefore, the evolution of the RS–μ0H|| loops with the sweep rate of the magnetic field provides conclusive proof that magnetization dynamics are responsible for the resistance peak when magnetic ordering and superconducting take place concomitantly.

Notably, the butterfly-shaped MR hysteresis is absent when measured with out-of-plane magnetic fields (Figs. 2a and S12). This result markedly contrasts with the behaviors of the superconducting 2DEG at the LAO/STO interface41, where MR hysteresis is observed under both in-plane and out-of-plane magnetic fields, indicating the existence of Bloch-type domain walls. Alternatively, Néel domain walls may form at the interface, considering that the easy axis of the 2D magnetism is parallel to the a-CZO/KTO(111) interface. This can be confirmed by comparing the SQUID measurements of the in-plane and out-of-plane ferromagnetic hysteresis loops (Fig. S13). Unlike STO, the strong SOC of 5d Ta could cause magnetic anisotropy21. An in-plane easy magnetic axis alignment has also been observed in the LaAlO3/KTO(111) and YAlO3/KTO(111) systems29. The magnetization lies in the plane of the interface, and the in-plane external magnetic field is able to drive the movement of domain walls, whereas the out-of-plane external magnetic field is not. Therefore, the dependence of MR loops on the magnetic field direction is strongly influenced by the easy-plane anisotropy of magnetism. In our work, the absence of the out-of-plane MR hysteresis corroborates the in-plane easy axis of magnetization. Notably, a similar behavior has also been observed in the normal-state 2DEG at the EuO/KTO (110) interface33. Here, we would like to emphasize that although the magnetization configuration in the actual system during magnetic reversal may be far more complicated than that proposed here, the underlying mechanism is similar.

The vortices will be pinned by domain walls, and their motion is enabled by the domain wall reconfiguration driven by magnetic fields. As vortices move, they can cause the phase of the superconducting order parameter to shift locally. If a vortex moves across a narrow superconducting wire, for example, it can lead to a complete circle of phase change (0–2π or −2π to 0), which is termed a phase slip, i.e., the intrinsic fluctuations of the superconducting order parameter46,47,48,49,50,51,52,53,54. It is therefore necessary to get the information about phase slips in the process of superconducting transition. In general, thermally activated phase slips (TAPS) occur when thermal energy is higher than the energy barrier between two neighboring metastable states, enabling the system to overcome the energy barrier via thermal fluctuations. When the phase slip is thermally activated, it can be described by the formula of the Langer–Ambegaokar–McCumber–Halperin (LAMH) theory47,48,49:

$${R}_{{\mbox{TAPS}}}=\frac{\pi {{{\hslash }}}^{2}\varOmega }{2{e}^{2}{k}_{{\mbox{B}}}T}\exp \left(-\frac{\Delta F(T)}{{k}_{{\mbox{B}}}T}\right)$$
(3)

where \(\Delta F(T)=\frac{8\sqrt{2}}{3}\frac{{{H}_{{\mbox{C}}}}^{2}(T)}{8\pi }S{\xi }_{{\mbox{GL}}}(T)\) is the energy barrier, \(\varOmega=\frac{L}{{\xi }_{{\mbox{GL}}}(T){\tau }_{{\mbox{GL}}}}{[\frac{\Delta F(T)}{{k}_{{\mbox{B}}}T}]}^{1/2}\) is the attempt frequency, \(T_{{{\mbox{C}}}}^{{{\mbox{zero}}}}\) is the zero-resistance temperature, \(\tau_{{\mbox{GL}}}=\pi{\hbar}/8k_{{\mbox{B}}}(T_{{\mbox{C}}}^{{\mbox{zero}}}-T)\) is the characteristic relaxation time in the time-dependent Ginzburg–Landau theory, HC(T) = HC(0) \((1-T/T_{{\mbox{C}}}^{{\mbox{zero}}})\) is the thermodynamic critical field, and ξGL(T) = ξGL(0)\((1-T/T_{{{\mathrm{C}}}}^{{{\mathrm{zero}}}})^{-1/2}\) is the Ginzburg–Landau coherence length. L is the length of the sample in the direction of the current, and S is the cross-sectional area of the sample. kB is the Boltzmann constant and ħ is the reduced Planck constant. Another mechanism for the activation of phase slips is quantum tunneling below the free energy barrier triggered by quantum fluctuations, thus referred to as quantum phase slips (QPS). When thermal energy is insufficient to overcome the free energy barrier, quantum tunneling becomes the dominant mechanism. QPS were theoretically investigated by Golubev–Zaikin51, and the temperature dependence of the resistance is given by refs. 50,51:

$${R}_{{\mbox{QPS}}}(T)={B}_{1}{R}_{{\mbox{Q}}}{S}_{{\mbox{QPS}}}\frac{L}{{\xi }_{{\mbox{GL}}}(T)}\exp (-{S}_{{\mbox{QPS}}})$$
(4)

where \({S}_{{\mbox{QPS}}}={B}_{2}(\frac{{R}_{{\mbox{Q}}}}{{R}_{{\mbox{N}}}})(\frac{L}{{\xi }_{{\mbox{GL}}}(T)})\) is the effective action, B1 and B2 are constants, RQ = h/4e2 is the quantum resistance, and RN is the normal state resistance.

A fingerprint of phase slips is the broadening of the superconducting transition process48. Therefore, we further analyzed the superconducting transition width (ΔT) for the samples with different doping levels, where ΔT is defined as T(RS = 0.8R4K)−T(RS = 0.2R4K). ΔT broadens from 0.12 to 0.32 K when nS increases from 3.8 × 1013 to 9.8 × 1013 cm−2. To verify the fluctuation origin of the broadening transition, we fit Eqs. (3) and (4) to the T–TC dependence of the RS/R4K in the transition region, as depicted in Figs. 3d and S14. The solid lines and dashed lines represent the fits to the TAPS and QPS models, respectively. In fact, TAPS is dominant at the temperature near TC, where thermal fluctuations are significant55. As shown by the solid curves in Fig. 3d, the TAPS model well reproduces the transition-induced rapid resistance drop near TC for all samples. However, when nS > 5 × 1013 cm−2, a weakly temperature-dependent resistance tail appears at lower temperatures, deviating significantly from the TAPS model. As the temperature decreases further, the magnitude of the free energy barrier relative to kBT increases, leading to a rapid reduction in the thermal-activation rate55. As this rate becomes negligible, tunneling through the energy barrier starts to dominate. As shown by the black dashed lines in Fig. 3d, the resistance tail behavior in the lower temperature range is well captured by Eq. (4), indicating that the QPS model performs well where the TAPS model fails. As the temperature decreases, a crossover from TAPS to QPS occurs for highly doped samples. Therefore, all samples follow the TAPS model near TC, while high-density samples undergo a TAPS-to-QPS transition as the temperature decreases. By comparing the transport behaviors under magnetic fields for samples with different nS, we found that the high-nS sample exhibits a shorter ξGL(0). Specifically, ξGL(0) = 61.63 nm and 16.43 nm for nS = 3.8 × 1013 and 8.9 × 1013 cm−2, respectively. The reduced coherence length in high-nS samples could enhance quantum fluctuations53, which becomes more pronounced at lower temperatures, enabling the TAPS-to-QPS transition.

To rule out disorder as the dominant mechanism influencing the change in transition width, we further compare the dependencies of the mean free path lmfp and the superconducting transition width ΔT on carrier density. As mentioned above, lmfp does not exhibit a significant dependence on nS (Fig. S11). However, the relationship of ΔT and carrier density indicates that the magnitude of ΔT is around ~0.12 K for low-nS samples (nS < 5 × 1013 cm−2) but increases to about ~0.29 K for high-nS samples (nS > 5 × 1013 cm−2). Obviously, the disorder mechanism fails to explain the obvious enhancement in ΔT as nS exceeds 5 × 1013 cm−2. However, this behavior is consistent with phase slips fitting results, that is, low-nS samples only follow the TAPS model, while a TAPS-to-QPS transition occurs in high-nS samples thus broadening the superconducting transition width. Therefore, the variation in transition width with carrier density is attributed to phase slips rather than disorder.

Given that the peak resistance ∆RS of the hysteretic MR loop originates from the motion of magnetic vortices cross the weak links between superconducting islands, we further try to adopt the TAPS model to describe the strong-temperature dependence of the ∆RS for the low-nS sample in Fig. 2f. Although the LAMH theory is derived for the case of μ0H = 0, it is applicable to the case with a finite magnetic field because the temperature dependence of both the energy barrier ΔF(T) and the coherence length ξGL(T) has the same form as that of μ0H = 056. Therefore, we conduct a theoretical analysis of the ∆RS–T relationship on the basis of the LAMH formula Eq. (3), adopting ∆RS = RTAPS(μ0H|| = 0.08 T)−RTAPS(μ0H|| = 0). As illustrated in Fig. 2f, the fitting curve reproduces the main features of experimental data around TC, indicating that the TAPS model captures the ∆RS–T peak behavior close to TC. Moreover, the validity of the TAPS model is also confirmed by the analysis of the dependence of peak resistance on both superconducting critical current and μ0H|| sweep rate, as detailed in Note 2 and Fig. S15. Therefore, it can be concluded that vortex motion can induce TAPS, which in turn affects the macroscopic properties of the superconductor and results in finite resistance. Notably, under an external magnetic field, vortex motion generated by magnetization dynamics is strongly influenced by the easy magnetic axis. In our work, the in-plane easy magnetic axis allows an in-plane external magnetic field to effectively drive magnetization dynamics, inducing vortex motion and phase slips, thereby generating a distinct hysteretic MR loop. In contrast, an out-of-plane external field, being orthogonal to the plane of magnetization rotation, is unable to drive magnetization dynamics. Consequently, the absence of the vortex motion implies that no phase slips occur, resulting in a lack of out-of-plane hysteretic behavior. Under a high initial positive in-plane magnetic field, the magnetic moment aligns with the applied field. As the field sweeps from positive to negative, the system retains its original magnetization even as the positive field decreases to zero. When the negative field approaches the coercive field, magnetic moment reversal occurs and magnetization dynamics will generate moving vortices that lead to phase slips, thereby resulting in finite resistance. Similarly, during a negative-to-positive sweep, magnetic moment reversal happens as the positive field reaches the coercive field. This explains the dependence of the resistance peak on the in-plane magnetic field sweep direction.

We summarize previous studies on upper critical field measurements of KTO(111) heterostructures14,15,57,58, as detailed in Table S3. Notably, most of these studies have focused on interfaces with higher carrier densities, while a few studies involving low-carrier-density samples have only measured out-of-plane magnetoresistance. These factors might have impeded the observation of magnetotransport evidence that supports the coexistence of superconductivity and ferromagnetism at KTO(111) interfaces in earlier studies.

To elucidate the origin of the ferromagnetism that coexists with superconductivity, we performed DFT calculations for the KTO(111) layer structure as schematically shown in Fig. 4a (see detailed modeling in Methods). The constructed KTO(111) layer consists of six KTO unit cells with an oxygen vacancy at one KO3 atomic layer (marked by the red dashed circle). Figure 4b shows the calculated orbital-resolved and layer-dependent DOS for the Ta atoms. Combined with the calculated band structures (Fig. S16), it can be concluded that the KTO(111) layer exhibits both metallic and magnetic properties, which coincides with our experimental observations. In Fig. 4b we presented the calculated magnetic moment per Ta atom and the layer-dependent polarization of the Ta 5d states. We can see that the ferromagnetism mainly arises from the first and second Ta atomic layers neighboring the oxygen vacancy, with the magnetic moment of ~0.13 μB/Ta. For Ta atoms located far from the oxygen vacancy, the calculated magnetic moments drop rapidly to ~0.059 μB/Ta. We further analyzed the orbital-resolved DOS for the Ta atoms and oxygen atoms around the oxygen vacancy (Fig. S17), and found that the ferromagnetism of Ta atoms is primarily contributed by the Ta 5dyz orbital which partially hybridizes with the Ta 5dx2−y2 and O 2py orbitals. The layer-by-layer orbital-resolved DOS for the d orbitals from Ta atoms (Fig. S17b) indicates that the magnetism is mainly attributed to the Ta 5dyz orbital in the first two TaO2 layers. The real-space spin density Δρ = ρ↑−ρ↓ is shown in Fig. S17c, and the spin-density isosurface near oxygen vacancy is mainly dominated by the spin-up Ta 5dyz orbitals. In fact, a Ta–Ov–Ta dimer with an oxygen vacancy located between two Ta atoms is formed when an oxygen vacancy is introduced into one KO3 atomic layer. Such absence of local covalent bonding and a reduction in local symmetry will strongly lower the energy of the d orbitals with lobes pointing to the vacancy59. As shown in Fig. S17c, the Ta dyz orbital lies in the yz plane with lobes aligning along the Ta–Ov–Ta direction, which significantly lowers its energy. Consequently, electrons transferred from the oxygen vacancy preferentially occupy the lower-energy state, specifically the Ta dyz orbital. The substantial spin splitting of the Ta dyz orbital near the Fermi level introduces an asymmetry between spin-up and spin-down contributions, leading to spin polarization.

Fig. 4: DFT calculations on KTO(111) layer structure.
figure 4

a Structural model of the KTO(111) layer with an oxygen vacancy at one KO3 atomic layer for DFT calculation. b The orbital-resolved and layer-dependent density of states for the Ta atoms. Positive and negative densities of states represent the spin-up and spin-down states, respectively. Numbers in the figure indicate the net magnetization in the unit of μB/Ta.

Full size image

The DFT study reveals that specific orbitals are responsible for magnetism, which may open the door to exploring the origin of magnetism through orbital-selective spectroscopic techniques. Highly sensitive resonant X-ray techniques, such as X-ray magnetic circular dichroism, X-ray linear dichroism, and X-ray absorption spectroscopy, are particularly well-suited to provide further experimental evidence and shed light on the orbital selectivity of magnetism. We believe that our theoretical findings provide the foundation and guidance for future research employing orbital-selective spectroscopic techniques to verify the contribution of specific orbitals to magnetism in KTO-based 2DEGs.

The physical origin of ferromagnetism for the momentum-confined interacting 2DEG can be understood based on the Stoner model60, where both Stoner parameter IS and the DOS at the Fermi level in the nonmagnetic state N(EF) determine whether the system favors a ferromagnetic or paramagnetic state. IS describes the strength of electron exchange, whereas N(EF) is inversely proportional to the kinetic energy of the electrons. The competition between the exchange and kinetic energy is taken into account by the Stoner criterion61, according to which the system would be unstable at the nonmagnetic state if N(EF)IS > 1. The Stoner FM transition can reduce the total energy of the system and thus stabilize the FM state. As shown in Fig. S18, the total non-spin-polarized DOS at the Fermi level N(EF) is significantly enhanced by the oxygen vacancy-induced gap states, and a sharp DOS peak located at the Fermi level indicates a Stoner-type instability. The Stoner parameter IS can be calculated through IS = ΔμB/M, where Δ is the spin splitting energy and M is the magnetic moment per Ta atom. The calculated N(EF) and IS are 7.0 state/eV and 1.4 eV, respectively. Obviously, the Stoner criterion N(EF)IS = 9.8 > 1 is satisfied. Therefore, the KTO(111) layer structure with oxygen vacancy-doped electrons has an FM ground state. In our spin-polarized calculations for the FM configuration, the DOS peak splits and shifts away from the Fermi level, with the calculated total energy of the FM state being 4 meV lower than that of nonmagnetic states. In the Stoner scenario, the magnetic strength in the normal state is expected to depend on the carrier density. The increase in saturation magnetization as nS grows is confirmed by the magnetic data obtained at 2 K by the SQUID magnetometer (Fig. S19). Given that the carrier density increases with the number of oxygen vacancies, DFT calculations on the KTO(111) layer structure with two oxygen vacancies were performed (Fig. S20). When the number of oxygen vacancies increases from one to two, the DOS at the Fermi level N(EF) in normal state increases from 7 to 9.4 states/eV. Accordingly, the total net magnetization increases from 0.504 to 1.335 μB. Therefore, the experimental results, combined with theoretical calculations, suggest an enhancement in ferromagnetism as the carrier density increases, further validating the Stoner scenario.

Notably, the magnetic moments originate from the itinerant electrons which are relatively localized around the oxygen vacancy, and the conductivity is contributed by itinerant electrons which come from each Ta atomic layer. According to this scenario, both ferromagnetism and superconductivity occur involving the Ta 5d electrons. However, the electron layers that are responsible for magnetic ordering and superconducting have different spatial distributions, as the oxygen vacancies are more localized to the a-CZO/KTO(111) interface. The local imaging measurement using scanning SQUID enables direct visualization of landscapes of ferromagnetism, paramagnetism, and superconductivity6,62. In fact, Julie et al. have directly imaged the coexistence of ferromagnetism and superconductivity at the LAO/STO interface using scanning SQUID, providing evidence for nanoscale spatial phase separation6. Therefore, the coexistence of ferromagnetism and superconductivity at the oxide interface could generally be understood as the simultaneous presence of both magnetic and superconducting phases and does not necessarily imply a strict spatial overlap between the two phases. Nonetheless, this does not diminish the profound interplay between magnetism and superconducting pairing at the oxide interface, which plays a central role in determining the emergence of unconventional superconducting states.

Recent studies have revealed the emergence of superconducting stripes in EuO/KTO systems15,33, induced by the ferromagnetic proximity effect of EuO. In our study, in-plane anisotropic transport behavior, a signature of superconducting stripes, was observed in the sample with nS = 4.3 × 1013 cm−2 but was absent in the sample with nS = 8.9 × 1013 cm−2 (Fig. S21). A similar nS-dependent in-plane anisotropic superconducting behavior has also been observed at the EuO/KTO (110) interface33.

As mentioned earlier, the superconductivity of 2DEGs at KTO interfaces is highly dependent on crystal orientation, with the TC being highest at the (111) interface and lower at the (110) interface. No superconductivity is observed down to 25 mK at the (001) interface. Notably, the crystal orientation dependence of the coexistence of magnetism and superconductivity is also a compelling issue, which could be explored in the (111) and (110) orientations. Although the (111) orientation appears to be the most favorable for superconductivity, this does not necessarily imply that it also favors the coexistence of magnetism and superconductivity, especially given the inherent competition between magnetism and superconductivity. Therefore, due to the complex interplay between magnetism and superconductivity, the dependence of their coexistence on crystal orientation requires further investigations.

In summary, definite signatures for the coexistence of 2D superconductivity and ferromagnetism are observed in the 2DEGs at the nonmagnetic a-CZO/KTO(111) interfaces, as indicated by the appearance of remarkable MR hysteresis in the superconducting background. Butterfly-shaped MR is detected in a wide temperature range, from the onset temperature of the superconductivity down to extremely low temperatures where the system enters deeply into the superconducting state. Specifically, the magnitude of the MR peak is found to be strongly temperature dependent, small at the onset temperature, reaching a maximum value around the superconducting transition temperature, and small again at even lower temperatures. These features are well captured by the model of TAPS. Moreover, the magnitude of the MR strongly depends on the field sweep rate, growing rapidly as the sweep rate increases. This observation reveals the close relation between hysteretic MR and magnetization dynamics. Theoretical analysis shows that the magnetic moments originate from the itinerant electrons, which are relatively localized around the oxygen vacancy and mainly arise from the Ta 5dyz orbital. The present work reveals the strong interplay between ferromagnetism and superconductivity, paving the way towards the exploration of quantum emergent phenomena.

Methods

Sample fabrication and characterization

Amorphous CZO films with the thickness of 10 nm were grown on (111)-oriented KTO single crystalline substrates by the technique of PLD with a ceramic CZO target. A KrF Excimer laser (wavelength is 248 nm) was employed. The repetition rate was 2 Hz and the fluence was ~2 J/cm2. During the deposition process, the substrate temperature was kept at 300 °C, while the O2 partial pressure PO2 varied from 7 × 10−3 to 1 × 10−5 Pa (details in Supporting Information). We obtained the 2DEGs at the a-CZO/KTO(111) interfaces with different doping levels, i.e., carrier densities. After deposition, the temperature of the sample was furnace-cooled to room temperature under the same atmosphere. Film thickness was determined by the number of laser pulses, which have been carefully calibrated by small-angle X-ray reflectivity. The surface morphology of the heterostructure was measured by atomic force microscope (AFM, SPI 3800N, Seiko). Lattice images were recorded by a high-resolution scanning transmission electron microscope (STEM) with double CS correctors (JEOL-ARM200F). The samples have been patterned into the Hall bar configuration using standard optical lithography and argon etching techniques. The central Hall bar bridge is 20 μm in width and 120 μm in length. XPS was measured using an Al Kα source in a Thermo Scientific ESCALAB 250X system. The magnetic properties were measured by a Quantum-designed vibrating sample magnetometer scanning superconducting quantum interference device.

Electrical transport measurements

The transport measurements including RS–T and Hall effect (temperature ranging from 300 to 2 K) were performed on a Quantum-designed physical properties measurement system. An applied current of 1 μA was used. The electrical transport measurements at lower temperature were conducted in a dilution refrigerator equipped with a vector-type magnet by a standard AC lock-in detection method. The AC excitation current of 10 nA at a frequency of 17.7777 Hz was applied using a Keithley 6221 current source. The corresponding AC voltage signals were measured using NF LI5650 lock-in amplifiers. For I–V curves measurements, a DC current was applied by Keithley 2400 and the corresponding voltage was measured using the same instrument. Ultrasonic wire bonding (Al wires of 20 μm in diameter) was used for electric connection. Unless specifically stated, Hall bar configuration was adopted to the measurements, the applied current direction is along [11-2] and the field sweep rate is 0.10 T/min.

Phase slips model

Superconductivity distinguishes itself by its remarkable ability to carry charge current with zero dissipation, which can be characterized by a macroscopic wave function Ψ(r) = |Ψ(r)|eiϕ(r), termed the order parameter. The amplitude and phase coherence of the order parameter vanishes at the transition temperature where the system goes back to the normal phase. However, even below the critical temperature, the phase coherence of the system can be affected by the so-called phase slips, i.e., phase fluctuations of the order parameter, which induce a finite resistance and therefore lead to a destruction of persistent currents. A phase slip is an elementary excitation of the order parameter due to thermal or quantum fluctuations, corresponding to a local suppression of its amplitude and a simultaneous jump of the phase by 2π. A phase slip can be thermally activated when the temperature is higher than the free energy barrier ΔF(T) between two neighboring metastable states and the system may overcome the barrier via thermal fluctuations. Another mechanism for the activation of phase slips is quantum tunneling below the free energy barrier triggered by quantum fluctuations, thus referred to as QPS. According to the Arrhenius law47,49, TAPS is extremely improbable for T ≤ ΔF(T)/kB. The prime fitting parameters used in the simulation of Eqs. (3) and (4) in Figs. 3c and S14, along with a comparison to experimental data, are listed in Table 1.

Table 1 Summary of the prime fitting parameters in the phase slips model
Full size table

Density functional theory calculations

First-principles calculations were performed using the projector-augmented wave method within the DFT63,64, as implemented in the Vienna Ab-initio Simulation Package65,66. To describe the exchange-correlation energy, we used the general gradient approximation with the Perdew-Burke-Ernkzerhof for solids parametrization67,68. The strong on-site Coulomb repulsion among Ta 5d electrons is corrected by using the DFT+U method, where U is the Hubbard parameter. The Ueff = 3 eV is used for Ta 5d states. In our computational investigation, the construction of a KTO(111) layer consists of six KTO unit cells with an oxygen vacancy at one KO3 atomic layer for simulating the electronic and magnetic properties. Our convergence standard requires that the Hellmann-Feynman force on each atom is less than 0.01 eV/Å and the absolute total energy difference between two successive consistent loops is smaller than 10−5 eV. The KTO(111) layer was fully optimized using a Γ-centered 4 × 4 × 1 k-grid, and the plane wave energy cutoff was set to 500 eV. The electronic structure calculations were performed by adopting a Γ-centered 9 × 9 × 1 k-grid.