Anomalous Josephson Coupling and High-Harmonics in Non-Centrosymmetric Superconductors with $S$-wave Spin-Triplet Pairing

We study the Josephson effects arising in junctions made of non-centrosymmetric superconductors with spin-triplet pairing having $s$-wave orbital-singlet symmetry. We demonstrate that the orbital dependent character of the spin-triplet order parameter determines its non-trivial texture in the momentum space due to the inversion symmetry breaking and spin-orbit interactions. The emergence of this pattern is responsible for the occurrence of an anomalous Josephson coupling and a dominance of high-harmonics in the current phase relation. Remarkably, due to the spin-orbital couplings, variations in the electronic structure across the heterostructure can generally turn the ground state of the junction from 0- to a generic value of the Josephson phase, thus realizing the so-called $\varphi$-junction. Hallmarks of the resulting Josephson behavior, apart from non-standard current-phase relation, are provided by an unconventional temperature and magnetic field dependence of the critical current. These findings indicate the path for the design of superconducting orbitronics devices and account for several observed anomalies of the supercurrent in oxide interface superconductors.


INTRODUCTION
In the presence of conventional Cooper pairing the current-phase relation (CPR) of a superconductorinsulator-superconductor junction is given by I J = I c sin(ϕ) [1], with I c being the critical current and ϕ the phase difference. While the sinusoidal shape is not always preserved [2], as in extended superconductornormal metal-superconductor junction or narrow ballistic weak links [3][4][5][6][7], a vanishing supercurrent and nondegenerate minimum of the Josephson energy at ϕ = 0 are robust marks of the CPR in conventional Josephson junctions (JJ). In this context, progress in materials science and nanofabrication have led to several physical cases with an unconventional CPR. Deviations from standard CPR indeed can manifest as a Josephson energy offset of a fractional flux quantum, leading to the so called φ 0 -junction [8][9][10][11][12][13][14][15][16][17][18] that violates time-reversal symmetry. On the other hand, for an offset of half-integer flux quantum a π-junction [19] is realized. Apart from anomalous phase shifts, the energy of the Josephson junction can keep the symmetry of phase-inversion but exhibits a minimum at values of the phase which is different from 0 or π, setting out a φ-Josephson junction. Typical requirements to achieve a φ-junction are i) combination of 0and π-Josephson couplings [20][21][22][23][24], ii) specific parameters range and geometric configuration of the junction [25][26][27][28], or iii) higher harmonics [20].
Observations of CPR in form of φ 0 or φ-junctions as well as 0-π transitions have been reported in a variety of devices based on InSb nanowires [29], ferromagnetic heterostructures [23,30], superconducting spin valves [31,32], and junctions with superconducting materials having non-trivial gap symmetry (e.g. iron picnitides [33], oxides interface [34,35], and cuprates [25,26,36]). Here, competing 0-and π-Josephson channels lead to a vanishing first harmonic with a consequent dominant role of the second harmonic in determining the Josephson CPR. Having a CPR with non-negligible harmonics higher than the second one is however quite unusual and difficult to achieve without fine tuning.
Recently, the combination of inversion symmetry breaking and multiple orbital degrees of freedom has emerged as an innovative route to tailor unconventional Josephson effects. This is mostly due to the expectation in acentric materials of a superconducting order parameter that goes beyond the canonical singlet-triplet mixed parity [37], as for the inter-band anti-phase pairing (e.g., s +− and s +−− ) [38,39] or pure even-parity inter-orbital spin-triplet pairs [40]. Striking experimental evidences of anomalous Josephson effects and supercurrents have been indeed reported in noncentrosymmetric superconductors based on oxides interface [34,35,41] which have been ascribed to the occurrence of competing 0-and π-Josephson channels as well as to second harmonics in the CPR. The microscopic origin of the observed effects is however not yet fully settled.
Motivated by this challenge and, in general, by the origin of anomalous Josephson effects in low-dimensional non-centrosymmetric superconductors (NCSs), we demonstrate how to achieve 0-, π-and φ-Josephson couplings together with robust highharmonics in Josephson junction by exploiting multiorbital degrees of freedom. In spin-triplet superconductors the Cooper pairs are typically described by the so-called d-vector whose components express the amplitude of the zero spin projection of the spin-triplet pairing states [42]. Here, we find that the Josephson effects are a striking manifestation of the non-trivial orientations of the d-vector in momentum space arising from the intertwining of its orbital dependent character with the inversion symmetry breaking and spin-orbit couplings. We demonstrate that 0 to π transitions and φ-Josephson phase can be generally realized and manipulated by varying the character of the electronic structure and the strength of inversion symmetry breaking interaction, thus being highly tunable through gating or electric field. Hallmarks of a Josephson junction based on this type of NCSs are: i) CPR with φ-phase and dominant high-harmonics, ii) anomalous temperature dependence of the critical current with a linear upturn at low temperature and iii) maximum of the critical current at finite applied magnetic field in the Fraunhofer pattern. We discuss how this type of pairing can account for the recent observations on supercurrents behavior in oxides interface superconductors [34,35].

ORBITAL-SINGLET SPIN-TRIPLET PAIRS AND d-VECTOR PROFILE ALONG THE FERMI LINES
We consider a multi-orbital 2D electronic system with spin-triplet s-wave pairing. In the normal state we have three bands arising from atomic orbitals spanning an L = 1 angular momentum subspace, such as d a orbitals with a = (yz, zx, xy). Here, we refer to dorbitals localized at the site of a square lattice assuming a C 4v point group symmetry. The breaking of mirror symmetry, in the plane of the junction, sets out a polar axis z leading to an orbital Rashba interaction (α OR ) that couples the atomic angular momentum L with the crystal wave-vector k in the standard form ∼ [L x sin(k y ) −L y sin(k x )] [39,[43][44][45][46]. The atomic spinorbit coupling (λ SO ) expresses the interaction between the spin and angular momentum at each site. Taking the basis of the local creation operator of electrons of dorbitals,Ĉ † k = [c † yz,↑k , c † zx,↑k , c † xy,↑k , c † yz,↓k , c † zx,↓k , c † xy,↓k ], the Hamiltonian can be generally expressed aŝ whose details of the matrix structure and the electronic dispersion with nearest-neighbor hoppings are provided in the Methods. In the superconducting state, the Bogoliubov-de Gennes (BdG) Hamiltonian is then directly constructed by including the pair potential∆(k). In the present study, our focus is on the symmetry allowed local (swave) spin-triplet pairing with orbital-singlet character and B 1 symmetry in the C 4v group. This type of pairing is energetically favorable when considering that interorbital interactions are dominant with respect to the intra-orbital ones [40]. The B 1 order parameter is described [40,47,48]  . A schematic illustration of the spin and orbital structure of the Cooper pair is shown in Fig. 1a. This type of pairing exhibits nodal points along the diagonal of the Brillouin zone ([110] direction) with nonzero topological number that is due to the chiral symmetry of the BdG Hamiltonian [40,[47][48][49][50][51][52]. Before analysing the Josephson effects in the junction formed by interfacing NCS [Fig. 1b], it is useful to consider the structure of the spin-triplet pairing in momentum space along the Fermi lines. This is evaluated by considering the anomalous components of the Green's function (see Methods for details). To this aim we choose two representative configurations whose electron density yields four [ Fig. 1c] and six [ Fig. 1d] Fermi lines around the Γ point. The inner Fermi lines are more isotropic, while the outer ones exhibit a more pronounced anisotropy. For the examined 2D tetragonal configuration, the xy orbital at k = 0 is lower in energy as compared to the (zx, yz) bands. We start by observing how the orientation of the d-vector at a given k-point is modified along the Fermi lines. For those bands having large Fermi momentum, i.e. bands 1 and 2 in Fig. 1e and Fig. 1f, we have obtained that the d-vector is mostly pointing along the x (y) direction for crystal wave vectors above (below) the diagonal. This is because the electronic configurations at the Fermi level have dominant xy character and mixing with zx and yz states through the spin-orbit and the orbital Rashba couplings. When we consider the bands closer to the center of the Brillouin zone [Figs. 1g-j], we observe that the d-vector exhibits a completely different pattern which is marked by an orientation that is mostly along the [110] direction [Figs. 1g,h] or perpendicular to it [Figs. 1i,j]. Such behavior is mostly due to the dominant (zx, yz) character of the electronic states and the fact that the B 1 pairing does not involve a direct coupling between such orbital states, thus implying an equal weight for d x and d y components. The results in Figs. 1cj provide evidence of the Cooper pairs along the Fermi lines having a spin-triplet configuration not uniform in orientation and amplitude that strongly depends on the orbital character of the corresponding electronic states.

MISALIGNMENT OF d-VECTORS ACROSS THE INTERFACE AND 0-φ JOSEPHSON COUPLING
Let us now consider the Josephson current in the NCS-NCS junction. For the computation of the Josephson current, we adopt the recursive Green's function method [53,54] assuming two semi-infinite superconductors and setting the pair potential in the Josephson junction aŝ ∆ L =∆ B1 and∆ R =∆ B1 e iϕ , with ϕ being the phase difference between the superconductors forming the junction [ Fig. 1b]. Below the critical temperature T c , Josephson current I c (ϕ) at the temperature T is then deter- < l a t e x i t s h a 1 _ b a s e 6 4 = " S F N V i n j f E 8 8 D N m J K a r R K 2 / g d z X 0 = " > A A A C h n i c h V G 7 T g J B F D 2 s L 8 Q H q I 2 J D Z F g r M h g U I g V i Y 0 l D 3 k k S M j u O u K G f W V 3 I S L x B 0 x s p b D S x M L 4 A X 6 A j T 9 g w S c Y S 0 x s L L w s m x g l 4 p 3 M z J k z 9 9 y Z k y u Z q m I 7 j P V 8 w s T k 1 P S M f z Y w N 7 + w G A w t L R d t o 2 n J v C A b q m G V J d H m q q L z g q M 4 K i + b F h c 1 S e U l q b E 3 u C + 1 u G U r h n 7 g t E 1 e 1 c S 6 r h w r s u g Q l W / U T m u h C I s a t q u W 6 R W V p k X K M K L s h d 2 z P n t m D + y V f f 5 Z q + P W G P y l T b s 0 1 H K z F r x Y z X / 8 q 9 J o d 3 D y r R q j k C h 7 v C c H x 0 i 5 X h T y Z r r M w K U 8 r N 8 6 6 / b z u 7 l o Z 4 P d s j f y d 8 N 6 7 I k c 6 q 1 3 + S 7 L c 9 c I U I P i v 9 s x C o p b s f h O L J F N R N I p r 1 V + r G E d m 9 S P J N L Y R w Y F e r e O S 1 y h K / i F m L A t J I e p g s / T r O B H C O k v v T u Q z A = = < / l a t e x i t > k x < l a t e x i t s h a 1 _ b a s e 6 4 = " x H n G 1 F e + 5 + k R U g g w S O x l a g j x J H w = " > A A A C h n i c h V G 7 T g J B F D 2 s L 8 Q H q I 2 J D Z F g r M h g U I g V i Y 0 l D 3 k k S M j u O u C G f W V 3 I U H i D 5 j Y S m G l i Y X x A / w A G 3 / A g k 8 w l p j Y W H h Z N j F K x D u Z m T N n 7 r k z J 1 c y V c V 2 G O v 7 h K n p m d k 5 / 3 x g Y X F p O R h a W S 3 a R s u S e U E 2 V M M q S 6 L N V U X n B U d x V F 4 2 L S 5 q k s p L U v N g e F 9 q c 8 t W D P 3 I 6 Z i 8 q o k N X a k r s u g Q l W / W O r V Q h M W Y G + F x E P d A B F 5 k j N A j j n E C A z J a 0 M C h w y G s Q o R N o 4 I 4 G E z i q u g S Z x F S 3 H u O c w R I 2 6 I s T h k i s U 1 a G 3 S q e K x O 5 2 F N 2 1 X L 9 I p K 0 y J l G F H 2 w u 7 Z g D 2 z B / b K P v + s 1 X V r D P / S o V 0 a a b l Z C 1 6 s 5 z / + V W m 0 O z j 9 V k 1 Q S J Q 9 2 Z O D O l K u F 4 W 8 m S 4 z d C m P 6 r f P e o P 8 f i 7 a 3 W K 3 7 I 3 8  mined [48,55] as witht N being the nearest-neighbor hopping matrix in the normal layer, the non-local Green's functionŝ G 01 (k y , iε n , ϕ) andĜ 10 (k y , iε n , ϕ) [ Fig. 1b and Methods], and the fermionic Matsubara frequency iε n = i(2n + 1)πk B T . Here, Tr ′ means that the trace is only in the electron space.
Having demonstrated that the d-vector has a nontrivial texture along the Fermi lines, we expect that the relative orientation of d-vectors in the two sides of the junction has a key role in setting out the Josephson effect. To this aim we employ two representative d-vectors configurations with a distribution of misalignment angles γ that can be ∼ 0 or close to π/2.  Figures  1,2,3). Here, we focus on the region of momentum space where the left-and right-side Fermi lines mostly overlap. We start considering electronic configurations that result into d-vectors that are about collinear on the two sides of the junction [ Fig. 2c]. Hence, taking into account those spin-triplet tunneling processes, we find that the derivative of Josephson current is positive at low ϕ and the CPR has a standard profile with maximum at about ϕ ∼ π/2 [ Fig. 2d]. Next, as reported in Figs. 2e-g, the texture of d-vector in the region nearby the diagonal of the Brillouin zone at the crossings of the left-and rightside Fermi lines indicate that the misalignment angle is about π/2. On the other hand, for momentum along the  [100] or [010] directions, the d-vectors are mostly aligned. In this case, the resulting Josephson current [ Fig. 2h] is marked by a change of sign in the derivative of Josephson current at low ϕ together with vanishing critical current for ϕ different from the time-reversal points at ϕ = 0 or π. This implies that the ground state realizes the so called φ-Josephson configuration. By tuning the amplitude of the orbital Rashba coupling or of the spin-orbit interaction it is possible to induce 0 to π Josephson phases and observe CPR dominated by harmonics higher than the second one close to the critical points (Figs. S2, S3 in the Supplementary Information).
To understand these Josephson effects, we recall that as shown in Ref. [56] for a specific geometrical design of the junction, the Josephson current J(ϕ) between spintriplet superconductors marked by d-vectors that are misaligned by an angle γ is expressed as with ϕ being the applied phase difference and Z t the transmission amplitude across the junction. This ex-pression indicates that in the tunneling process the spintriplet pairs with opposite spin polarization undergo an antiphase shift to keep the time-reversal symmetry with an amount that is proportional to the misalignment angle γ. Here, it is immediate to deduce that for a value of γ that is about π/2, the current J(ϕ) at small phase difference can change sign, thus turning the Josephson phase behavior from 0 to π. Since for the examined NCS-NCS junction the d-vector has a variation of the orientation in the momentum space, it is useful to consider the Josephson behavior resulting from the superposition of different misaligned d-vectors. Taking the representative case of a pair of Josephson channels with different configurations of d-vector misalignment, as in Figs. 3a,b, one can find that when the angles γ and γ 1 are both close to π/2, the resulting current is marked by non-vanishing and comparable amplitude of the first four harmonics (i.e I 1 ,I 2 ,I 3 ,I 4 ). In this case, the CPR has a profile that yields a φ-Josephson coupling.
Another peculiar hallmark of the Josephson current in the examined NCS-NCS junction is represented by the occurrence of high-harmonics in the CPR with nonvanishing and comparable amplitude. We find that this behavior can be also ascribed to the nontrivial misalign-  Fig. 3c]. Within this effective description, we find that high-harmonics generally occur in the CPR [Fig. 3d]. The CPR has an anomalous profile with a positive derivative at low ϕ and multiple oscillations due to the competing high harmonics. Since the d-vector texture is not homogeneous along the Fermi lines, it is natural to expect that various channels with inequivalent mismatch angles and transmission amplitudes cooperate to build up such behavior for the supercurrent. In Figs. 3e,f we provide evidence of the CPR for the NCS-NCS junction, whose profile is significantly marked by high-harmonics at low temperature which in turn tend to get suppressed by increasing the temperature towards the transition into the normal state. Finally, we evaluate the response to an applied magnetic field that acts to modulate the critical current amplitude (details in the Supplementary Information). We find that the resulting Fraunhofer pattern [Figs. 3i,j], as a consequence of the intrinsic competing 0 and π Josephson couplings, does not have a maximum at zero magnetic flux, as expected in conventional spin-singlet JJ junctions. Interestingly, even for a CPR with a 0-Josephson coupling, as in Fig. 3i, the critical current has a local minimum rather than a maximum at zero applied magnetic field [Figs. 3i,j]. This demonstrates that the intrinsic tendency of having competing Josephson channels in the case of NCS superconductors with s-wave orbital-singlet and spin-triplet pairing has clearcut and detectable signatures in the magnetic field response.

DISCUSSION
We now discuss the impact of our findings in SrTiO 3 based hetero-structures, such as LaAlO 3 /SrTiO 3 (LAO-STO) [57,58]. The LAO-STO is an ideal 2D electron system with non-centrosymmetric multi-orbital superconductivity [38,40,59] exhibiting a remarkable control of the superconducting critical temperature by electrostatic gating [60][61][62] together with Rashba spin-orbit coupling [63,64] and the occupation of the Ti 3d orbitals (d xy ,d zx ,d yz ) [65,66]. The superconducting phase exhibits several anomalous properties that cannot be easily addressed within a conventional spin-singlet scenario. Clearcut unconventional observations are provided by the superconducting gap suppression nearby the Lifshitz transition [67,68], the anomalous magnetic field dependence of critical current in weak links [34,41] and uniform nanowires [69], and several in-gap bound states probed by tunneling spectroscopy [70]. The significant role of inhomogeneities also poses fundamental questions on the nature of the superconducting state in LAO-STO interface, excluding p-wave spin-triplet pairing as a candidate, and pointing to an even-parity (s-wave) multi-band superconductivity which is robust to disorder.
Recently, superconducting transport measurements in nano-devices [34,35] have provided direct and significant evidences of an orbital dependent unconventional pairing. The central experimental findings demonstrate: i) an anomalous enhancement of the critical current at weak applied magnetic field, ii) an asymmetric response with respect to the magnetic field direction, and iii) the supercurrent anomalies are gate dependent, getting pronounced when reaching the Lifshitz transition at the onset of the occupation of the (d zx ,d yz ) bands. In this context, our results can account for the central observations of the superconductng transport measurements. Our analysis is based on planar 2D junctions as for the experimental configuration and the unveiled Josephson effect, due to the inter-orbital spin-triplet pairing, leads to π-phase and high-harmonics in the Josephson current. As demonstrated by the analysis of the Fraunhofer pattern, the combination of π-phase and high-harmonics provides the enhancement of the critical current at weak applied magnetic field and we find it to be robust with respect to variations of the electronic parameters in the two sides of the junction. This result can be particularly relevant for the LAO-STO where the gating is related to an inhomogeneous distribution of the electron density in the 2D electron gas and of the inversion symmetry breaking at the interface. Concerning the asymmetric response with respect to the magnetic field direction, our results indicate that the combination of π-phase and high-harmonics is not sufficient to yield the effect. One has to consider it in a spatially inhomogeneous array of Josephson junctions [35]. Hence, we expect that the observed asymmetry of the Fraunhofer pattern with respect to the magnetic field can be dependent on the spatially homogeneity of the superconductor in terms of granularity and inclusion of nanometric sized islands.
Our findings highlight the relevant role of s-wave interorbital spin-triplet pairing in noncentrosymmetric junctions for achieving competing 0-and π-Josephson couplings together with high-harmonics in the current phase relation. Differently to single-band superconductors, the resulting CPR is tied to the character of the electronic structure and thus can be potentially tuned in a controlled way by changing the electron density and the Rashba coupling through gating. This aspect has a direct impact on the transport properties of LAO-STO junctions and points to distinct design of Josephson devices. We also point out that the distribution of the spin moment of the Cooper pairs in momentum space is a key quantum resource for the generation of CPR with dominance of high-harmonics. Remarkably, in all these configurations the magnetic field response reveals a Fraunhofer pattern with a minimum of the critical current at van-ishing field. While nonstandard magneto-electric effects [71][72][73][74] are based on the magnetic field tunability of the d-vector orientation [52,71,75], our findings highlight the role of the orbital rather than the spin degree of freedom. Hence, in perspective, due to the demonstrated intrinsic orbital tunability of the d-vector texture by orbital Rashba coupling, our findings set out innovative routes to design orbitally driven magneto-electric effects and Josephson devices for superconducting orbitronics with orbital control of the supercurrent. M.C., P.G. and Y. Correspondence should be addressed to Mario Cuoco, email: mario.cuoco@spin.cnr.it.
H SO andĤ is (k) stand for the atomic spin-orbit coupling and inversion symmetry breaking terms, respectively. Here,σ i=0,x,y,z denote the Pauli matrices in spin space andL j=0,x,y,z the t 2g -orbital angular momentum operators projected onto L = 2, in the [d yz , d zx , d xy ] basis. For the (100) oriented surface, the local termũ(k y ) and nearest neighbor hopping matrixt(k y ) can be explicitly derived. The local termũ(k y ) is given bỹ u(k y ) =û 0 (k y ) +û SO +û is (k y ), and the nearest neighbor hoppingt(k y ),

Recursive Green function method
In the superconducting state the Green's function is expressed as with the anomalous Green's function for the orbital indices (α, β) which is given bŷ Since we consider the interorbital even-frequency/spintriplet/orbital-singlet/s-wave pair potential, we focus on the d-vector in this pair amplitude. We define the d-vector for the interorbital even-frequency/spintriplet/orbital-singlet/s-wave pair amplitude, Then, we also define with j = L, R superconductors, respectively. In the numerical calculation, both d (L) and d (R) are real number. Thus, the misalignment of d-vectors between two superconductors is obtained by with θ j (a) = arg[d (j) x + id Here, we report the main steps for the determination of the Josephson current by the recursive Green's function method [54].ũ with the chemical potentials in left (right)-side superconductors µ L (µ R ), and the normal layer µ N , respectively.
with the transparency t int . In this calculation, we fix the transparency as t int = 1.0. First, we calculate the semi-infinite surface Green's functionĜ L andĜ R in the left and right-side superconductors, respectively [54]. When we add a normal layer, we obtain the surface Green's functionsĜ L0 (k y , iε n ) and G R1 (k y , iε n , ϕ), with the local term in the normal layerũ N , the tunnel Hamiltoniant L0 =t 1R , and the fermionic Matsubara frequency iε n = iπ(2n + 1)k B T . For the connection of two superconductors with a normal layer, we calculate the local Green's functions, and the nonlocal Green's functions,

JOSEPHSON CURRENT
We report about the evolution of the Josephson current as a function of the electronic parameters. We start by considering a configuration with an amplitude mismatch of the chemical potentials between the left (µ L ) and right (µ R ) side of the junction (Fig. 1). We track the evolution of the first four harmonics of the Josephson current as a function of µ R at a given value of µ L = 0.0 ( Fig. 1(a) as well of the CPR for two representative cases of µ L (( Fig. 1(b),(c)). For these configurations, corresponding to equal amplitude of the spin-orbit (λ SO ) and orbital Rashba couplings on the two sides of the junction, the CPR is marked by a 0-type Josephson coupling. It is more interesting to observe the evolution as a function of the spin-orbit and orbital Rashba interactions assuming a fixed amplitude mismatch of the chemical potentials. In Figs. 2,3 we show the evolution of the first four harmonics and the corresponding CPR as a function of the orbital Rashba coupling (α OR ) and the spin-orbit coupling (λ SO ). We find that the first harmonic can change sign at a critical value of α OR or λ SO thus signalling a 0-to π-phase transition. In proximity of the transition the CPR is dominated by high harmonics as manifested by the presence of double maximum and asymmetric amplitudes of the maximum and minimum in the current phase relations.
To evaluate the dependence of the critical current on an applied magnetic field H we follow the standard approach for a short junction with length L and width W that is thread by a magnetic flux Φ = LW H. Then, we can expand the current phase relation in all the harmonics and consider that the spatial dependent vector potential associated to the magnetic field modifies the phase difference ϕ across the junction with a factor 2πΦx whereΦ = Φ/Φ 0 , with Φ 0 being the unit of magnetic flux, andx indicating the ratio x/W , i.e. the scaled length with respect to the lateral size of the junction. The general expression of the supercurrent expressed via an expansion in harmonics is given by I(Φ, ϕ) = n w 0 I n sin nϕ + 2πΦx dx . (1) The resulting critical current at a given applied magnetic flux is obtained by maximizing I(Φ, ϕ) with respect to ϕ.