Spectroscopic Signatures of Gate-Controlled Superconducting Phases

We investigate the tunneling conductance of superconductor-insulator-normal metal (SIN) and superconductor-insulator-superconductor (SIS) heterostructures with one superconducting side of the junction that is electrically driven and can exhibit $\pi$-pairing through a modification of the surface inversion asymmetric couplings. In SIN tunneling we find that the variation of the electrically driven interactions generally brings an increase of quasi-particles in the gap due to orbitally polarized depaired states, irrespective of the inter-band phase rearrangement. The peak of SIN conductance at the gap edge varies with a trend that depends both on the strength of the surface interactions as well as on the character of the gate-induced superconducting state. While this shift can be also associated with thermal effects in the SIN configuration, for the SIS geometry at low temperature the electric field does not yield the characteristic matching peak at voltages related with the difference between the gaps of the superconducting electrodes. This observation sets out a distinctive mark for spectroscopically distinguishing the thermal population effects from the quantum gate-driven signatures. In SIS the electrostatic gating yields a variety of features with asymmetric peaks and broadening of the conductance spectral weight. These findings indicate general qualitative trends for both SIN and SIS tunneling spectroscopy which could serve to evaluate the impact of gate-control on superconductors and the occurrence of non-centrosymmetric orbital antiphase pairing.


I. INTRODUCTION
Understanding the interplay of electricity and superconductivity is a fundamental problem that stands out for its general relevance, the great impact it can have for accessing, controlling or driving new phases of quantum matter, and the enticing perspectives for the development of future quantum technologies. Electric field effects have been successfully employed to drive or control the superconducting phase in materials with low-to-moderate carrier density [1][2][3][4], in thin films [5,6,8] and interfaces [9], down to the two-dimensional limit [7,10,11]. In this context, to achieve a control of the chargecarrier density by gating the materials have to be in a lowdensity regime and thin enough to avoid shielding of the electric field. In fact, for superconductors with large electron density the strong electric screening typically allows small changes in the superconducting critical temperature [12][13][14][15].
Recently, the discovery of unconventional gating effects in metallic superconductors [16] has challenged the current view of how electric field can affect superconductivity. In particular, the overall experimental observations which have been so far accumulated, provide a variety of unexpected fingerprints of the consequences of electrostatic gating on superconductors. The major findings demonstrate the following effects: reduction or full suppression of the critical supercurrent in nanowires [17][18][19] and Dayem bridges [20][21][22][23], manipulation of superconducting phase in interferometric setups [22,24], enhancement of phase fluctuations [25], increase of in-gap quasiparticle population [18], and weak interrelation between the critical magnetic field and the critical voltage associated with the vanishing supercurrent transition [26]. Remarkably, a suppression of the critical supercurrent has been also achieved in fully suspended gate-controlled Ti nano-transistors [27] and with ionic-gating [28], thus posing bounds on the role of charge injections or electron leakage within the observed phenomena.
The emerging experimental scenario is definitely opening fundamental challenges about the coupling of static electric field and superconductivity. Indeed, the impact of electrostatic gating on the metallic superconductor points to different channels of interaction. One possibility is that the application of large electrostatic gating may cause the injection of highly energetic quasiparticles. In this framework, the induced excitations in the superconductor, with energies that are greater than the superconducting gap, would mainly lead to thermally driven population unbalances and non-equilibrium effects which substantially may end up into a suppression of the amplitude of the superconducting order parameter. Along this line, the increase of the quasiparticle in-gap spectral weight in the tunneling conductance has been indeed recently interpreted as an evidence of high-energy injection of quasiparticles into the superconductor [18,19].
Another scenario places its roots in the fundamental interaction between an electrostatic potential and the electronic states on the surface of the superconductor through the modification of the strength of the inversion asymmetric interaction [29]. Electric field breaks inversion symmetry on the surface and can lead to orbital polarization of the electronic states through orbital Rashba couplings that, however, are active only in the spatial region where the electrostatic field can penetrate the metal, i.e. of the order of the Thomas-Fermi length [30]. In this context, it has been indeed recently recognized that an analog of the spin Rashba coupling [31,32] at the surface of low-dimensional electron systems or noncentrosymmetric materials arises due to the coupling between the atomic orbital angular momentum L and the crystal wavevector k [29,[33][34][35]. For materials having configurations close to the Fermi level with non-vanishing atomic angular momentum, an intrinsic crystalline potential or an applied electric field, breaking spatial inversion or mirror, yield nonlocal odd-parity matrix elements among atomic orbitals contributing to the Bloch states. The resulting orbital Rashba 0-SC π-SC π z -SC N α OR (intra-layer Orbital Rashba) Electric Field Electric Field λ (inter-layer Orbital Rashba) 1. (a) Illustration of the superconducting slab embedded into an electrostatic field generated by gates with voltages V G 1 and V G 2 . For the multilayered superconducting thin-film, due to the screening on the Thomas-Fermi atomic length, one can distinguish between the surface layers (blue), where the electric field fully penetrates and affects the electronic structure, from the inner layers (yellow) that are not directly influenced by electric field induced interactions. On the top surface layer we show a representative pattern of the chiral textures in the Brillouin zone due to the electric field induced orbital-Rashba couplings. (b) Schematic illustration of the superconducting phase diagram which is obtained by applying an electric field through the modification of the induced intra-(α OR ) and inter-(λ ) orbital polarizing interactions at the surface layers of the superconducting thin film. Within the phase diagram, there occurs a conventional spin-singlet superconducting state with orbitally polarized surface states, that is labeled as 0-SC (c). Otherwise, the electric field generates π-pairing with the superconducting order parameter having a π-phase shift among different bands at the Fermi level. Here, we consider a three-band modeling with ∆ a ,∆ b , and ∆ c the corresponding spin singlet order parameters. π z −SC stands for an inhomogeneous superconducting phase along the confining direction of the thin film (z) with layers having zero interband phase shift alternating with layers marked by π-pairing (e). Finally, N indicates the normal metal state configuration with strong orbitally polarized depaired states at the Fermi level. An overall schematic illustration of the 0-, π-and π z -SC phases is reported in (c-e) with the layer and band dependent superconducting order parameters determined by iterative computation on the lowest energy state. Schematic of the SIN (f) and SIS (g) junction with one side subjected to electrostatic gating. The coloured surface on the gated superconductor sketches the penetration of the electrostatic field.
The essential consequence of such orbital polarizing effects is that the electric field is able to break the inter-band superconducting phase rigidity or fully suppress the amplitude of the order parameter [29]. Then, two main consequences are in place as schematically depicted in the phase diagram of Fig.  1. The gating can rearrange the band dependent superconducting phases with a π-shift (i.e. resulting into a π-or π zphase depending on the spatial homogeneity of the π-pairing, Fig. 1 (d)-(e)), or it can completely suppress the amplitude of the order parameter by increasing the population of depaired orbitally polarized quasi-particles until reaching the normal phase N. Both phases (i.e. the π-states and the normal metal with strong orbitally polarized Fermi surfaces) can be in principle marked by a vanishing supercurrent, however the underlying mechanism that leads to the supercurrent suppression is fundamentally different. In fact, for the π-states the supercurrent suppression is due to the superconducting phase frustration of different bands (i.e. it mainly arises from the interference of the phases of the superconducting order parameters at the Fermi level), while for the electrically induced normal metal phase it is essentially due to a suppression of the amplitude of the superconducting order parameter. Remarkably, the application of a magnetic field can help to discern between the two scenarios concerning the way the electric field impacts on the amplitude and phase of the superconductor [26]. In this context, the weak dependence of the critical magnetic fields on the critical voltages and viceversa seems to be mostly compatible with a scenario with the electrostatic gating driving the transition from a conventional superconducting phase to a superconducting state with π-pairing [29] where the inter-band phase frustration is responsible of the vanishing of the supercurrent.
It is interesting to point out that π-pairing is also expected to occur in unconventional superconductors in the absence of external perturbations as indicated, for instance, by the recent observations in iron-based materials [44] or theoretically proposed at oxides interface [45,46]. Hence, the addressed problem has a wider framework and our results can have a more general application. Indeed, they can provide insight into the tunneling spectroscopy of superconducting thin films where the inter-band phase reconstruction already manifests in the absence of an applied electric field and the application of other drives can affect the orbital polarization of the electronic state on the surface and further reconstruct the superconducting phase.
Starting from this physical outlook it is relevant to consider whether spectroscopically one can assess the nature of the gate-driven superconducting phases and search for fingerprints which can be employed to understand the way the electric field affects the superconductivity. For this purpose, we explore the tunneling conductance of superconductorinsulator-normal metal (SIN) and superconductor-insulatorsuperconductor (SIS) heterostructures with one superconducting side of the junction being subjected to the electrostatic gating. Our strategy is to evaluate the tunneling conductance for SIN and SIS configurations in the phase space spanned by the electric field through the variation of the interactions at the surface layers of the superconductor moving within different types of π-paired configurations. The key target is to track the evolution of the tunneling conductance in all phases that are obtained within the scenario of an electrostatically triggered orbital polarization at the surface ( Fig. 1(b)) and extract the main spectral signatures. In this way one can get a significant insight into the spectroscopic response that can result when considering the gate-driven superconducting transitions. Starting from the SIN tunneling spectroscopy the electric field is shown to yield an increase of spectral weight in the gap, due to the orbitally polarized depaired states, and this trend is observed independently of the allowed inter-band phase rearrangements. The voltage position of conductance peak nearby the gap edge varies with a tendency that is sensitive to the amplitude of the orbital Rashba interactions as well as the character of the gate-driven superconducting state. While this shift can be also associated with thermal effects, for the SIS geometry at low temperature the electric field does not yield the characteristic matching peak for voltages that correspond to the difference between the gaps of the superconductors. This observation thus represents a distinctive mark for spectroscopically disentangling the thermal population effects from signatures that are mainly related to a variation of the electric field. Apart from such feature, in SIS the electric field yields a variety of marks with asymmetric redistribution of the spectral weight of the main conductance peak. These characteristics when tracked as a function of the electric field indicate distinct trends for both SIN and SIS configurations with characteristic tunneling behavior in gate-controlled superconducting phases.
The paper is organized as follows. In Sect. II we present the model and the methodology for the tunneling conductance analysis. Sect. III is devoted to the main results by focusing on the SIN and SIS tunneling spectroscopy. In Sect. IV we provide the discussion on the resulting effects and the concluding remarks. In the Appendix we provide details of the profile of the electrostatic potential close to the surface of the superconductor. Additionally, we present the basic elements for the derivation of the tight-binding model in the presence of an electrostatic potential at the surface and the temperature dependence of orbital dependent superconducting order parameters obtained by self-consistently solving the gap equations.

II. MODEL AND METHODOLOGY
In this section we present the model and the methodogy which have been employed in order to evaluate the tunneling conductance for the examined SIN and SIS geometries. To achieve the goal, we need determine the density of states of the superconducting side of the junction that is subjected to the electrostatic gating. For this part of the heterostructure, the superconductor is considered to have a slab geometry with n z layers ( Fig. 1 (a)) and we assume a conventional s−wave spin-singlet pairing.
Before going into the details of the results for the SIN and SIS conductance we deepen the discussion about the main as-pects of the model that has been used to describe the effects of the applied electric field. Concerning the inclusion of the gate voltage in the employed model Hamiltonian, there are three relevant working hypotheses. Firstly, we are considering a metallic superconductor (i.e. with a high density of charge carriers and thus large Fermi surface). In this context, the electrostatic potential in the inner layers of the superconductor is vanishing due to screening effects for distances from the surface that are greater than the Thomas-Fermi length λ TF (typically of few unit cells along the out-of-planeẑ-direction). The second aspect concerns the character of the gate voltage. We assume to insert the layered superconducting system in a capacitor like structure (open circuit configuration), that fixes the value of the gate voltage to be opposite at the two sides of the layered superconductor (see Fig. 1(a)). Then, a final point refers to the investigated electrical regime. We deal with stationary conditions, i.e. constant in time electrostatic gating. Hence, the reference Maxwell equation to be considered is that one for the scalar electrostatic potential. Basically, within the capacitor configuration we have ∇ 2 V (r) = 0 in the region between the plates and the surface of the superconductor, while ∇ 2 V (r) − k 2 TF V (r) = ρ(r)/(4πε 0 ) inside the superconductor, with k T F = λ −1 TF being the inverse of Thomas-Fermi length, and ρ(r) the induced charge at the surface.
Taking into account the above assumptions, one can determine the profile of the electrostatic potential nearby the surface of the superconductor. Since the screening effects in a metal leads to an exponential decaying of the electric field within a distance of the order of few unit cells, the effects of the electrostatic potential V (z) are confined at the surface layers of the slab. Hence, without loss of generality and due to the geometry (i.e. the thin film has a dimension in the xy plane that is larger than the thickness along the z-direction) and symmetry of the problem, the electrostatic potential can be taken as only dependent on the z coordinate and exponentially vanishing inside the superconductor. The solution of the electrostatic equation for this configuration can be handled and is reported in Appendix A. To proceed further and construct the tight-binding model for the layered superconductor, one has to evaluate the matrix elements of V (z) for the selected multiorbitals Wannier basis having non-vanishing atomic angular momentum L. For this purpose, we take a 3×3 sector which is suitable for p− or d− (e.g. t 2g multiplet of the d-manifold in cubic/tetragonal symmetry) states with L = 1, being a description that can be effectively applied to a large class of materials. In this manifold the electrostatic potential has diagonal and off-diagonal elements. The diagonal terms renormalize the chemical potential while the off-diagonal ones lead to the orbital Rashba coupling. In the employed model, we are neglecting the modification of the chemical potential at the surface layer because the correction to the electron density is substantially affecting only the surface layer and is negligible for a metallic system (see Appendix A), with also a minor impact on the superconducting order parameter [30]. In turn, the off-diagonal terms lead to orbital Rashba couplings. In Appendix B we provide the main steps for the derivation of the orbital Rashba interaction at the surface.
Hence, since the applied electric field on the surface of the superconductor is parallel to the out-of-planeẑ-direction and thus the electrostatic potential close to the surface, at the linear order in z, can be described by a potential V s = −E s z with E s being constant in amplitude (assuming the electric charge e is unit). Taking into account the steps for the derivation of the tight binding model (see Appendix B) [26,29,[33][34][35], the matrix elements of V s in the Bloch basis can yield an intra-(α OR ∼ E s ) and inter-layer (λ ∼ E s ) inversion asymmetric interactions, whose ratio depends on the inter-atomic distances and distortions occurring at the surface layers [29]. For convenience one can indicate as (a, b, c) the (yz, xz, xy) d−orbitals. Then, after introducing the creation d † α,σ (k, i z ) and annihilation d α,σ (k, i z ) operators with momentum k, spin (σ = [↑, ↓]), orbital (α = (a, b, c)), and layer i z , one can construct a spinorial basis . In this representation, the Hamiltonian can be generally expressed in a compact way as [26,29]: witĥ where the orbital angular momentum operatorsL have com- within the (yz, xz, xy) subspace, τ i (i = x, y, z) are the Pauli matrices for the electron-hole sector, and δ i, j the Kronecker delta function. Due to the anisotropy and directionality of the d-orbitals, the kinetic energy for the in-plane electron itinerancy is expressed by , with η being a term that takes into account deviations from the ideal cubic symmetry. Other long-range terms or processes involving inter-orbital hoppings that are activated by distortions do not change the quality of the addressed phenomenology [29]. We assume that the layer dependent spin-singlet order parameter (OP) is non-vanishing only for electrons belonging to the same band. Here, we have neglected the pair-hopping term of the form In the early works on two-band superconductivity [47][48][49], the pairhopping term has been included to remove the degeneracy among the configurations with different phase difference between the intra-orbital superconducting order parameters in systems without single particle orbital hybridization. For the examined model, we have inter-orbital mixing and the 0-and π-phases are separated in energy (see for instance Fig. 3 in Ref. [29]), with the 0-configuration being the groundstate in absence of an applied electric field even without the pair hopping V αβ terms. In this respect, the inclusion of the term V αβ would further stabilize the 0-state in the phase diagram and one would need a slightly larger inter-layer orbital Rashba coupling to induce the transition from the 0-to the πphase. Thus, since for the study of the tunneling conductance we have investigated all the regimes in the phase diagram, we do not expect qualitative changes in the results.
Concerning the evaluation of the superconducting order parameters, , we performed it by computing the trace of the pairing operator P over all the eigenstates |n, k of the Hamiltonian associated to negative energies E n,k < 0 at zero temperature (or to all energy configurations at finite temperature weighted by the Fermi function). Since the eigenstates |n, k depend on ∆ α (i z ) and the orbital Rashba interactions couple the momentum with the local angular momentum, the gap equations of the orbital dependent order parameters are strongly coupled between each other. To make evident this dependence one can express the gap equation in the following way the Fermi-Dirac distribution at given temperature T . In this framework, we have computed self-consistently the order parameters for each orbital character as a function of the temperature (see Appendix B). Since the deviation from the canonical BCS behavior does not alter the qualitative outcome of the analysis, in order to have a uniform comparison of the various regimes, we have determined the conductances by taking the BCS profile for the temperature dependence.
Here, N = n x × n y sets the dimension of the layer in terms of the linear lengths n x and n y , while translation invariance is taken in the xy-plane and n z is the number of layers along the z−axis (Fig. 1). We notice that g is not modified by the electric field. This is physically plausible because due to screening effects the electric field cannot induce an inversion asymmetric potential inside the thin film on distances from the surface that goes beyond the Thomas-Fermi length.
For completeness, we point out that the ground state of the investigated model Hamiltonian cannot sustain a time dependent phase dynamics under the applied voltage. In order to get a non-trivial time dependence of the Josephson phase it is crucial to have two superconducting condensates separated by an insulator that is thin enough to keep an energy difference of the Cooper pairs as given by the applied voltage. In turn, for the examined superconducting multi-layered system, we have only one superconducting condensate. The voltage difference across the surface layers does not lead to an energy difference among the Cooper pairs in neighbor layers as they are not isolated and they are in good electrical contact since we are dealing with a metallic system. Although there is an inhomogeneous electrostatic potential at the surface com-pared to the inner layers, both the single particle states and the Cooper pairs are not localized in the corresponding regions of the superconductor. Additionally, there is no phase difference between the superconducting order parameters at the top and bottom surface layers (or among those at the surface layers and the neighbor inner layers) because such phase difference would be associated to a non-equilibrium state with a current flowing across the slab. Such flow would then lead to charge imbalance across the superconductor which is energetically unfavourable and prevented by the electrical screening in the metallic superconductor. According to these observations, it is useful to remark that this regime for the phase in the superconducting system can be broken if the superconducting materials are marked by inhomogeneities in the region where the electric field penetrates that yield isolated islands or clusters separated by insulating barriers, i.e. spatially distributed effective weak-links near the surface of the superconductor. Concerning the model Hamiltonian in Eq. 2, we also observe that it has only one term that is not compatible with the charge conjugation transformation involved in the time reversal symmetry operation. This term is the inter-layer orbital Rashba interaction which changes sign under a complex conjugation transformation. Such aspect, however, has an impact only on the phase of the superconducting order parameter. This consequence has been taken fully into account by allowing a complex value for the orbital-dependent order parameters in the self-consistent simulation of the superconducting order parameter. Moreover, since this is an inter-layer orbital dependent charge transfer that is activated by the electrostatic gating only at the surface of the superconductor, it does not affect the structure of the s-wave pairing in the inner layers which keeps its form of zero momentum with electron pairs having opposite spin at k and −k close to the Fermi surface. We observe that the electric field is a surface perturbation and does not influence at all the pairing interaction in the inner layers of the superconductor, thus the employed conventional form of the superconducting order parameter is a physically valid assumption for the examined model. For completeness, we also point out that a kind of pair density wave along the z− direction (i.e. π z phase) can be achieved in the strong coupling regime of orbital Rashba interactions [29]. The occurrence of this type of pair density reflects the fact that, due to the electrostatic potential, a reconstruction of the order parameter can occur along the zdirection. This outcome does not imply a variation in the structure of the in-plane pairing.
To proceed further, we consider representative profiles of the superconducting order parameters for the various bands including the states with π-pairing ( Fig. 1 (c)-(e)) and correspondingly determine the layer dependent density of states by computing the energy spectra and the eigenvectors of the Hamiltonian.
Tunneling spectroscopy in junctions or in scanning tunneling experiments is an important probe of the density of states (DOS) of target materials. The current across the junction at a finite applied bias V is generally expressed as a convolution of the DOS of the normal metal electrode, N n (E), that one of the material upon examination, N s (E), and the Fermi-Dirac FIG. 2. Thermal evolution of the normalized tunneling conductance, G sn (V ), for the SIN junction evaluated in the configuration with vanishing external electrostatic gating (i.e. for α OR = λ = 0 within the implemented modeling) for several temperatures T below the superconducting critical temperature T c . ∆ 0 is is the amplitude of the lowest energy excitation for the zero electric field state assuming a uniform profile over the layers of the superconducting order parameter. The computation has been performed for a superconducting electrode with n z = 6 layers. For the temperature dependence of the superconducting order parameter we employ the phenomenological BCS expression ∆(T ) = ∆ 0 cos(πT /2T c ) [51].
3. Normalized tunneling current I ss (V ) (a) and conductance G ss (V ) (b) for an SIS 2 junction in the absence of an applied electric field (i.e. for α OR = λ = 0) as a function of temperature. For S 2 we are considering a conventional BCS superconductor with ∆ S 2 = 2.5∆ 0 . The vertical gray lines mark the position of the bias voltage V − = ∆ S 2 − ∆ S and V + = ∆ S + ∆ S 2 at zero temperature, with ∆ S = ∆ 0 . In the insets we show a zoom in the voltage bias region where the matching peaks occur. function f (E), while the conductance is given by dV . For our purposes we assume that the normal metal in the SIN has a constant density of states in the range of investigated voltages while the non-gated superconducting electrode in the SIS is described by a conventional BCS-type density of states.
The presence of the Fermi-Dirac distribution in the expression of the SIN tunneling conductance limits the energy resolution (i.e. ∼ k B T ), with the thermal smearing that often precludes the detection of features with low intensity and small energy separation in the density of states. Indeed, the inspection of the DOS for the superconducting electrode shows that there are more fringes and spectral features at energies that are smaller than the thermal smearing. In Fig. 2 we report the thermal evolution of the SIN tunneling conductance by assuming that the superconducting electrode is described by the tight-binding model Hamiltonian of Eq. (2). Here, we take a representative case with n z = 6, the analysis with a larger number of layers does not affect the overall spectroscopic outcome. Additionally, we assume that the superconducting order parameter is spatially uniform and independent of the orbital index. This position has two main motivations. It is substantially consistent with the typical profiles that are obtained by fully minimizing the free energy with layer dependent superconducting order parameters, as explicitly depicted in Fig.  1(c)-(e). Additionally, in order to single out the spectroscopic fingerprints that uniquely arise from the modification of the electric field, through the orbital Rashba interactions, it is convenient to neglect the small amplitude differences between the superconducting order parameters in different bands. This strategy is quite neat because it helps in clarifying the impact of the electric field without mixing with other effects related to pairing amplitude variation. Indeed, the presence of multigap structure of superconductors leads to other features in the tunneling conductance and its spectroscopy detection is an intricate problem that stands alone even without the application of the electric field [50].
Furthermore, for the junction's electrode subjected to the electrostatic gating, we explicitly compute the density of states for a multilayered superconductor as described by the model Hamiltonian in Eq. 1. Here, the analysis of the electronic energy states and the corresponding eigenfunctions is performed within the whole Brillouin zone and, then, by evaluating the momentum integration of the Bogoliubov energies taking into account the band and layer dependence of the superconducting order parameters (Fig. 1 (c)-(e)).
The outcome in Fig. 2 sets the reference for our analysis for the configuration with vanishing electric field. As expected the conductance has a peak at voltages V max ∼ ∆ 0 and the position of V max moves to values of ∼ 1.5∆ 0 while approaching T c . Furthermore, there is a full gap at low temperature that fills up uniformly from zero to one as the temperature reaches the critical temperature.
Replacing the normal electrode with a superconducting one in the SIS geometry, one can overcome the thermal broadening issues by exploiting the non-linearities of the density of states of the probing superconducting electrode. Indeed, it is well known that the SIS spectroscopy can lead to a significant enhancement of the resolution of the tunneling spectroscopy. To set the reference we consider an SIS heterostructure where one side of the junction, labeled as S 2 , is described by a conventional BCS-type density of states including the Dynes phenomenological parameter Γ, The parameter Γ is usually employed to quantify the consequences of the pair breaking processes, while N 0 is the value of the density of states at the Fermi level in the normal metal phase.
In Fig. 3 we show the SIS tunnel conductance for the case of zero applied gating (i.e. α OR = λ = 0). The outcome describes the typical profile of the SIS tunnel current and conductance. Here, we assume that the gap amplitude of the electrically gated superconductor (S) is ∆ S = ∆ 0 while that one of the conventional electrode (S 2 ) is ∆ S 2 = 2.5∆ 0 . The current profile exhibits the typical peak when the applied voltage matches the difference of the gaps, V − = ∆ S 2 − ∆ S (1.5∆ 0 ), and the rapid increase at the sum of the gaps, V + = ∆ S + ∆ S 2 (3.5∆ 0 ). The peak at V − becomes visible and more pronounced with reduced temperatures T /T c getting above about 0.5. This structure in the current profile gives a nonmonotonous behavior for the conductance with a change of sign close to V − . Furthermore, the rapid upturn of the current corresponds to a peak in the conductance that moves to lower voltages while increasing the temperature upon reaching T c . We notice that for voltages above V + there are small amplitude oscillations in the conductance that reflect the intrinsic features of the electronic structure in the layered superconducting thin film described by the Hamiltonian in Eq. 1.

III. RESULTS
In this Section we present the evolution of the SIN and SIS tunneling conductance due to the application of an electrostatic gating onto one superconducting electrode constituting the junction. The gating can drive the superconducting state, resulting into 0-, π-, or π z -configurations with orbitally polarized surface states (Fig. 1). Since the amplitude of the microscopic parameters directly affected by the electric field is generally proportional to E s with a factor that depends on materials characteristics, we have taken different trajectories in the phase space spanned by α OR and λ . The focus is on the low temperature regime because it can allow to single-out the consequences of the electrostatic gating and distinguish from the canonical thermal effects in conventional junctions.

A. SIN conductance
Let us start by considering the SIN tunneling conductance for two representative trajectories in the phase space assuming that a variation of E s allows us to move along a path with λ that is smaller (Fig. 4(a)-(c)) or greater (Fig. 4(d)-(f)) than α OR , respectively. The range of variation for α OR and λ is taken from 0 to 5 ∆ 0 , with ∆ 0 setting the superconducting gap in absence of external perturbations. This is an energy range for the electrostatic gating that is suitable for the experimental observations and also microscopically consistent for covering the allowed phase transitions from 0-to π-paired states. We remind that λ and α OR depend on E s in such a way that larger amplitudes for E s would typically lead to a normal phase with significant orbital polarization at the Fermi level. The strategy we follow is to explore the response of the three different phases (i.e. 0-, π-, and π z -SC) that can be achieved in a multi-band superconductor upon the application of an electrostatic gating independently of their stability in the phase diagram. This general approach allows us to get a complete view of the behavior of the tunneling conductance for each superconducting state in the whole phase space. Beginning from the 0-SC state we find that the application of the electric field generally leads to a filling up of the gap (Fig. 4(a),(d)). Depending on the ratio between λ and α OR there can be a high (λ > α OR ) or low (λ < α OR ) rate of increase of the quasiparticle in-gap population. This result clearly indicates that the surface inter-layer coupling λ is a source of depairing by orbitally polarizing the electronic configurations at the Fermi level. This is expected indeed because the λ term acts as an effective orbital current and can locally break the time reversal symmetry. Here, we notice that in the regime of large λ ( Fig. 4(d)) the filling of the gap is accompanied by the formation of a bump at voltages that are below the gap edge of the order of ∼ 0.5∆ 0 . Since the impact of the electrically driven couplings is primary on the surface, we argue that such feature mostly originates from significant surface reconstruction of the in-gap electronic states.
Another remarkable spectroscopic aspect of the electrically driven SIN junction refers to the evolution of the maximum of conductance. Starting from the case at zero applied electric field we observe that the peak tends to move to voltages larger than ∆ 0 while approaching the boundaries of the explored phase space. This shift is generally accompanied by a suppression of the amplitude of the conductance peak. A close inspection of the evolution indicates a non-monotonous behavior of the voltage position (V max ) of the conductance peak. In fact, in the regime of λ and α OR being comparable to ∆ 0 , corresponding to the A-C path, with λ > α OR in Fig.4(d), V max tends to decrease to lower voltage intensity. The overall trend keeps being non-monotonous also along the D-E path. Although the inward shift is typically small, this is a qualitative distinctive mark that can be also experimentally detectable.
Let us now move to the π-paired phase. In this context, we deal with two possible configurations with uniform (π-SC) or spatially inhomogeneous (π z -SC) order parameters along the z-direction of the superconducting thin film (Fig. 1 (d),(e)). For the uniform π-SC the tunneling conductance exhibits similar features of those found for the 0-SC case. Indeed, the increase of the electric field amplitude along the two identified trajectories in the (λ , α OR ) phase space generally yields a growing in-gap conductance contributions ( Fig. 4(b),(e)). Furthermore, the in-gap bump obtained approaching the point F in the λ > α OR path (Fig. 4(e)) is robust as it occurs in the π-SC state too. When considering the evolution of the maximum of conductance, some differences are encountered with respect to the 0-SC configuration. Indeed, we find that the peak moves monotonously to high voltages irrespective of the type of trajectory in the phase space. In particular, we find that for the path with λ > α OR the maximum can shift up to V ∼ 1.5∆ 0 . Thus, in the uniform π-SC phase the conductance peak is suppressed and can move to large voltage amplitudes. The shift is particular pronounced for the trajectory with λ > α OR .
Finally, we study the SIN conductance for a superconducting electrode that is marked by π z -pairing. In this configuration, we have that 0-and π-SC coexist in the thin film as they occur in different layers along the confining direction. More specifically, for the analyzed layered superconductor the top/bottom surface layers have no interband phase shift, while the remaining layers are marked by π-pairing. Hence, apart for the intrinsic inter-band phase shift at the Fermi level, the superconducting order parameter exhibits an amplitude spatial modulation with a sign change when moving along the z-direction. This is similar to unconventional pairing with spatial modulation of the amplitude as for the case of the Larkin-Ovchinnikov state [52], such variation is expected to result into nodal excitations or to generally bring an effective suppression of the superconducting gap. This is indeed obtained when considering the limit of vanishing electric field. The zero bias conductance in that case is large indicating a substantial contribution of in-gap electronic states nearby the Fermi level. Then, a variation of the orbital Rashba couplings drives an increase of the in-gap conductance which however is of the same order of magnitude of that observed in the uniform 0-and π-SC phases. The evolution of the peak in the conductance also shows a trend with an increase of the voltage and a suppression of the amplitude as one moves along the two selected trajectories. In particular, the shift of the conductance maximum is amplified when proceeding along the path corresponding to λ > α OR reaching a value of V max ∼ 1.5∆ 0 .
At this point it is instructive to compare the profile of the conductance as a function of the applied electrostatic gating (Fig. 4) with that one obtained by varying the temperature (Fig. 2). In Fig. 2 we have shown the canonical temperature dependence of an SIN junction for a conventional BCS superconductor. The thermal effects induce a filling up of the ingap conductance with a suppression of the maximum amplitude and a shift of its peak position at high voltages while approaching the superconducting transition temperature T c . We notice that the increase of the conductance is quite uniform within the gap while the zero bias amplitude approaches its normal value close to T c . On the other hand, it is unexpected that at given temperature the application of an electrostatic gating leads to a similar qualitative trend when compared to the thermal drive. This implies that the SIN tunneling conductance with the superconducting electrode subjected to an external electric field can yield outcomes which in principle are difficult to get distinguished from thermal effects. In particular, the analogies are evident for the π-phases. Instead, the non-monotonous evolution of the maximum of the conductance upon the application of an electrostatic gating for the 0-SC represents a mark of unique electric effects in SIN tunneling detection, although such consequence is expected only in the regime of small electric field amplitudes.

B. SIS conductance
Having shown that the SIN tunneling conductance cannot directly provide distinctive features to disentangle the role of electrostatic gating from thermal effects, we discuss the major differences that can be extracted once considering the SIS spectroscopic response.
As for the SIN tunneling, the strategy is to search for clearcut fingerprints which might be used to assess the nature of electrically driven superconducting phases by exploring the 0-, π-and π z -SC states. In Sect. II we have shown that for SIS 2 heterostructures with inequivalent superconducting electrodes, a maximum in the conductance is expected to occur at voltage positions that correspond to the difference and sum of the superconducting gaps, with the former being pronounced only above an effective temperature. Hence, it is relevant to ask about the impact of the electric field in the energy range where the matching peak and the characteristic conductance features manifest. Furthermore, taking a temperature where the lowest voltage matching peak is vanishing and not detectable, it is also worth evaluating whether the electric field can induce it, thus mimicking thermal effects as in SIN tunneling spectroscopy, or bring substantially different consequences.
In Figs. 5(a) and (d) we present the evolution of the SIS 2 conductance for the 0-SC phase assuming that ∆ S = ∆ 0 and ∆ S 2 = 2.5∆ 0 . This choice is convenient for having well separated characteristic voltages at the sum and difference of the superconducting gaps. At zero applied electric field the SIS tunneling conductance is basically marked by a main peak at V + = 3.5∆ 0 with satellite features at higher voltages oscillating around the normal state value. The variation of the surface orbital Rashba couplings along the paths with λ greater or smaller than α OR leads to distinct trends for the main peak. Indeed, for the trajectory with λ < α OR (Fig. 5(a)) we observe a reduction of the intensity of the main conductance peak at V + that evolves into a dip accompanied by two peaks at voltages below and above V + . On the other hand, when considering the λ > α OR (Fig. 5(d)) we have that the reduction of the spectral weight of the peak at V + is accompanied by a downward shift of the conductance structures at low voltage bias. The result is a changeover from a single peak at V + into two structures, staying at V + and another one developing at V * ∼ 3∆ 0 . The presence of the conductance peak at V * is connected with the bump structure developing in the DOS of the superconductor subjected to the electric field, at voltages inside the gap of the order of 0.5 ∆ 0 (Fig. 6). We argue that the substantial filling up of the gap when λ is larger than α OR is the source of the peak formation in the SIS conductance at lower voltages with respect to V + . Furthermore, for this regime of coupling there is also a significant spectral weight redistribution at high voltages above V + .
Let us now move to the π-paired phases. As reported in Fig. 5(b),(e), for the uniform π-SC, we observe that the SIS conductance exhibits a profile having strong similarity with that of the 0-SC phase. Indeed, the variation of the electrically driven surface couplings suppresses the main peak at V + which evolves into multiple structures with a broad distribution with respect to the case without an applied electrostatic gating. The width of the conductance spectra is generally larger for the λ > α OR trajectory, thus confirming the dominant role of the λ coupling in yielding non-standard spectroscopic features. However, in this regime, differently from the 0-SC (Fig. 5(d)) we find that the distribution of the spectral weight keeps being symmetric around V + .
Finally, we discuss the SIS spectroscopic behavior assuming that one superconducting electrode is in the π z superconducting phase (Fig. 5(c),(f)). As we have understood from the SIN tunneling conductance, the π z -SC is associated with a significant spectral weight in the gap. Then, already in the regime of small λ and α OR , a non-vanishing conductance structure is obtained at voltage bias of V 2 = ∆ S 2 . The resulting shoulder in the conductance spectra extends down to about 2∆ 0 while The gated superconductor S has a layered geometry with n z = 6 layers.
moving along the trajectories in the phase space from A to F (Fig. 5(c),(f)). We also notice that the changeover of the conductance is not equivalent when comparing the paths with λ < α OR or λ > α OR . In the former, the maximum of the conductance stays at V + while for the case with λ > α OR the peak at V + becomes broad and flat.
A relevant outcome of the above analysis is that although the main peak of the SIS conductance gets modified there are no traces of other structures emerging at a voltage bias that matches the difference of the superconducting gaps in the two electrodes. Furthermore, by a comparative inspection of the SIS tunneling conductance we can state that the main conductance peak gets generally suppressed by the application of an electrostatic gating resulting into a spectral redistribution with asymmetric features and broad structures whose profile is sensitive to microscopic details and to the character of the electrically induced superconducting phases.
The observed features are directly linked to the modification of the DOS in the superconductor within the various induced phases (Fig. 6). The quasiparticle peak at the gap edge typically shifts to high energies and gets suppressed in amplitude. The trajectory with λ < α OR leads to an ingap structure below the gap edge that moves inwards within the gap accompanied with a slight increase of spectral weight as the strength of the electrostatic gating grows. In this regime, the main differences between the 0-,π-and π z phases occur close to zero energy. In fact, the π z configuration shows an enhanced filling up of the spectral weight in the gap due to the spatial sign change of the order parameter along the z-direction of the su-perconducting slab. As a consequence, one observes a non vanishing conductance already around V 2 = ∆ S 2 (Figs.5(c)-(f)) even for small values of α OR and λ . On the other hand, when considering the trajectory with λ > α OR the suppression of the peak at the gap edge is generally more pronounced as well as the increase of the in gap spectral weight. When considering the π-and π z -phases the peak at E ∼ ∆ 0 loses its spectral weight and develops into multiple structures both at higher energies and inside the gap.

IV. DISCUSSION AND CONCLUSIONS
In conclusion we have unveiled the main spectroscopic features that would arise in SIN and SIS junctions with one electrode being a conventional superconductor subjected to an electrostatic gating that can drive the formation of orbital antiphase paired phases by means of strong orbital Rashba effects at the surface. Our theoretical analysis reveals a complete set of distinctive marks for spectroscopically accessing the electrically induced superconducting phases and suggests an experimental way to disentangle thermal population unbalance from effects that are mainly due to the a variation of the electric field strength.
The analysis has been performed in a way to track the general evolution of the driven superconducting phases moving from the regimes of weak to strong gating amplitude. This approach is particularly useful in a system where the electrostatic field can lead to different types of transitions of the ) is scaled to the normal state one (N n (E)). The DOS from B to F are shifted along the vertical axis by 2. We investigate two representative paths in the [λ , α OR ] parameters space with (a) λ < α OR (tr 1 ) and (e) λ > α OR (tr 2 ), respectively. (b),(c) and (d) present the evolution of the density of states for the trajectory tr 1 in the 0-, π-and π z -phases, respectively. The behavior of the density of states for the trajectory tr 2 is depicted in (f),(g) and (h) for the 0-,π-and π z -phase, respectively. In-gap structures are more pronounced for the path tr 2 . The suppression and shift of the peak at the gap edge is generally observed with an amplitude that depends on the character of the superconducting phase. type 0-π or 0-π-π z [26,29]. Here, the transitions are substantially accompanied by a phase rearrangement rather than an amplitude reconstruction of the band dependent order parameters. In this scenario, our results demonstrate that, since both SIN and SIS conductances share similar qualitative characteristics in the spectra while moving in the electric parameters phase space, we expect to observe smooth changeover across the transitions which would substantially manifest by a variation of the spectral weight in the conductance. This is an important outcome of the study because it tells that there will not be dramatic reconstructions of the conductance. On the contrary, one needs to track changes of the states in the gap or close to the conductance peak to identify the character of the electrically reconstructed superconducting states in the superconducting thin film. In particular, depending on whether the inter-layer orbital current processes are more relevant than the intra-layer ones the conductance both for the SIN and SIS can exhibit more pronounced features in the conductance. This is for instance the case of the shift to high voltages of the conductance maximum in the SIN configuration for the π z phase compared with 0-SC configuration which instead has a nonmonotonic behavior.
Another relevant concluding observation refers to the role of the thermal effects. We have demonstrated that the SIN conductance exhibits variations associated to a change of the electric field amplitude that can share wide similarities to those due to an increase of the temperature in the absence of an external electrical perturbation. This implies a difficulty to disentangle the two effects especially when the gating can also lead to heating or thermal gradient in the device. For this reason, we propose to combine the SIN with SIS tunneling probes. Our analysis indeed shows that at low temperature the electric field would not lead to the characteristic matching peak which is instead observed due to thermal excitations in conventional SIS junctions with different gap amplitudes. Finally, for the SIS configuration we have found that the application of an orbital polarizing electric field leads to a peculiar reconstruction of the main conductance peak with multiple structures whose profile is intricately tied to the character of the gate driven superconducting states. Concerning the charge reconstruction, in order to get an estimate of the induced electron density at the surface of the metal, we assume the electric field to be uniform in the xy plane and the induced charge to have a peaked distribution, delta function, at the surface along the zdirection, i.e. ρ(r) = ρ S δ (z). Then, for the electric field generated by the capacitor plates at a distance L C and the thin superconducting film of thickness d placed inside the capacitor at equal distances from the two plates, there are two regions along the z direction with inequivalent profile of the electrostatic potential V (z). Region While for the region II inside the superconductor we have for 0 < z < d/2. By integrating the Poisson equation close to the surface, i.e. at z = −d/2, one can obtain the expression for the induced in-plane charge density ρ S = 4πε 0 V G /L C + λ T F sech(2λ T F d) 2 . For typical amplitudes of the applied voltages for the devices realized in Ref. [17] and related works, i.e. V G ∼ 50V with split gates placed at a distance L C ∼ 100 nm, we have that the variation of the electron density in the unit cell is of the order of δ n = 0.05 (assuming that the thickness d = 50 nm, the in-plane unit cell size is A c = 4Å 2 , and λ T F = 4Å ). Such induced charge for a metal having an average density at the Fermi energy of n e ∼ 1.0 per unit cell, implies that the overall charge reconstruction due to the electrical gating is typically of the order of few percents. Since this variation is confined at the surface, in turn, it has also a small impact on the superconducting order parameter both at the surface and in the inner layers. Similar conclusions on the negligible impact of the charge reconstruction in BCS metallic superconductors have been demonstrated in Refs. [30,53].

B. Electrostatic potential and orbital Rashba coupling at the surface
The external electric field on the surface of the superconductor is parallel to the out-of-planeẑ-direction and thus the previously derived electrostatic potential close to the surface, at the linear order in z, can be described by a potential V s = −E s z with E s being constant in amplitude (assuming the electric charge e is unit). Hence, to construct the model Hamiltonian one has to consider a Bloch state representation and explicitly evaluate the matrix elements of the electrostatic potential V s . Since the translational symmetry is broken along theẑ-direction due to the finite thickness of the thin film and the presence of the electric field, the out-of-plane momentum, k z , is not a good quantum number. Thus, a representation with a Bloch wave function associated to each layer is the most appropriate one to evaluate the effects of the electric field and the way it enters in the tight-binding model. One can use the index i z to label different Bloch wave functions along theẑdirection as follows with the Bravais vector R ν,i z identifying the position of the atoms in the x − y plane for the layer labeled by i z , β indicating the atomic Wannier orbitals, and N the total number of atomic sites. A central aspect in the derivation is that the atomic Wannier functions φ β (r − R ν,i z ) span a manifold with non-vanishing angular momentum L, e.g. p or d states.
The intra-layer matrix elements of the electrostatic potential can be then expressed as with l and m spanning the orbital space, and c ψ the normalization factor of the Bloch state. Since the Wannier orbitals are significantly localized around the atomic position, the dominant terms are for nearest neighbor atoms. When evaluating those contributions in Eq. 5 for l = m we obtain matrix elements that involve the L x and L y orbital angular momentum components and that are odd parity under in-plane spatial inversion [29]. The resulting term yields the intra-layer orbital Rashba coupling as in Eq. 2. A similar derivation can be made for the inter-layer orbital Rashba term. The form of inter-layer term is due to the structure of the matrix elements of the electrostatic potential between Wannier functions in neighbor layers along theẑ−direction. By evaluating these matrix elements [29], it turns out that the electric field can induce an orbital polarization on nearest neighbors atoms in adjacent layers only if one allows for displacements/distortions of the atoms in the plane with respect to the high-symmetry positions. A detailed derivation of the interaction is presented in Ref. [29].
C. Temperature dependence of the superconducting order parameters We have computed self-consistently the order parameters for each orbital component as a function of the temperature in order to compare it with the canonical BCS behavior and assess the possible consequences or deviations concerning the presented results. The outcome is reported in Fig. 7 for various representative cases in the phase space. As one can observe, the orbital dependence of the order parameters manifest in two aspects. The crystal field splitting and the orbital anisotropy of the kinetic energy can introduce an amplitude imbalance of the superconducting order parameters in the 0-SC at zero temperature with ∆ c being smaller than ∆ a and ∆ b (Fig. 7 (a),(b),(c). This asymmetry evolves in temperature      Behavior of the orbital dependent superconducting order parameter ∆ α (α = a, b, c), calculated with an iteratively self-consistent approach, for a system with n z = 6 layers, t ⊥ = 1.5t, orbital Rashba coupling α OR = 0.2t. In (a)-(c) we report the superconducting order parameter at T = 0 as a function of the inter-layer orbital Rashba coupling λ focusing on the evolution of the the first three layers from the top surface (i.e. i z = 1, 2, 3). The other corresponding layers starting from the bottom surface have the same amplitude. For the selected parameters we have a transition from the conventional (0-SC) to the π−phase (π-SC) for λ = 0.2t. The dots indicate the values of λ for which the orbital dependent superconducting OP have been evaluated as a function of the temperature in the panels (d)-(i). In (d)-(f) we show the behavior of ∆(T /T c ) for the layers i z = 1, 2, 3, also marked in the insets, at λ = 0.1t. In this regime, the superconductor exhibits a conventional SC phase (0-SC), with all ∆ α having the same sign. In (g)-(h) the evolution ∆(T /T c ) is plotted for λ = 0.3t, being an interaction amplitude that stabilizes an orbital anti-phase π−pairing state (π-SC), with a superconducting order parameter of a given orbital character having a π-shift in the phase respect to the other ones. In panels (f) and (h) the gray lines indicate the temperature dependence of the canonical BCS gap as given by the phenomenological expression ∆(T ) = ∆ 0 cos( πT /T c ). Here, ∆ 0 is the value of ∆ α (i z , λ ) at T = 0. For the c band in (f) to mimic the BCS profile for the smaller gap we assume an effective lower critical temperature. In all panels δ 0 is the value of ∆ a (λ = 0) at the surface layer and it is used as a reference scale.
with the smaller gap typically having a suppression at temperatures in proximity of the superconducting transition T c (Fig.  7 (d),(e),(f). As shown in Fig. 7, this behavior is substantially independent of the layer position. The resulting behavior in the 0-SC phase indicates that all the superconducting order parameters deviate from the canonical BCS profile (gray dotted line in Fig. 7(f)). Moreover, the inter-layer orbital Rashba coupling modifies the tail of the small superconducting order parameter close to the critical temperature when considering the π-phase (Fig. 7(g),(h),(i)). In turn, for λ values stabilizing the π-phase, we have that the temperature dependence of the superconducting order parameters follows quite well the BCS behavior ( Fig. 7(i)). A similar trend is also obtained for the π z -phase (Fig. 8) where a good degree of matching with the BCS profile is obtained (Fig. 8(d)). On the basis of these results, we conclude that the tunneling conductance evaluated by means of the BCS profile are suitable in the πand π z phases for fully reproducing the temperature dependence of the multi-orbital superconducting order parameters. On the other hand, for the 0-SC configuration one can expect that quantitative changes might occur, especially close to the critical transition. We argue that, although these outcomes are directly related to a specific microscopic model, the inequivalent temperature dependence of the 0-and π-phases can man- (1) (a) π z -SC T = 0 λ/t=0.6; αOR=2.0t (1) (c) layer iz=2 π z -SC λ/t=0.6; αOR=2.0t π z -SC λ/t=0.6; αOR=2.0t Behavior of the orbital dependent superconducting order parameter ∆ α (α = a, b, c), calculated with an iteratively self-consistent approach, for a system with n z = 6 layers and for a set of parameters such that the superconductor is in the π z -phase, namely t ⊥ = 0.9t, orbital Rashba coupling α OR = 2.0t and inter-layer OR coupling λ = 0.6t. In (a) we report the superconducting order parameter at T = 0 as a function of the layer index i z . In (b)-(d) we show the behavior of ∆(T /T c ) for the layers i z = 1, 2, 3, also marked in the insets. In fig. (d), the gray lines indicate the temperature dependence of the canonical BCS gap as given by the phenomenological expression ∆(T ) = ∆ 0 cos( πT /T c ), where ∆ 0 is the value of ∆ α (i z , λ ) at T = 0. In all panels ∆ 00 a (1) is the value of ∆ a (α OR = 0, λ = 0) at the surface layers and it is used as a reference scale in this figure. ifest in a different evolution of the main conductance peaks for the SIN and SIS spectra in the corresponding phases.