Spectral fingerprints of the non-linear dynamics of driven superconductors with dissipation
H. P. Ojeda Collado Gonzalo Usaj José Lorenzana C. A. Balseiro




 



Centro Atómico Bariloche and Instituto Balseiro,
Comisión Nacional de Energía Atómica (CNEA)–Universidad Nacional de Cuyo (UNCUYO), 8400 Bariloche, Argentina
Instituto de Nanociencia y Nanotecnología (INN), Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET)–CNEA, 8400 Bariloche, Argentina

Centro Atómico Bariloche and Instituto Balseiro,
Comisión Nacional de Energía Atómica (CNEA)–Universidad Nacional de Cuyo (UNCUYO), 8400 Bariloche, Argentina
Instituto de Nanociencia y Nanotecnología (INN), Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET)–CNEA, 8400 Bariloche, Argentina


ISC-CNR and Department of Physics, Sapienza University of Rome, Piazzale Aldo Moro 2, I-00185, Rome, Italy


Centro Atómico Bariloche and Instituto Balseiro,
Comisión Nacional de Energía Atómica (CNEA)–Universidad Nacional de Cuyo (UNCUYO), 8400 Bariloche, Argentina
Instituto de Nanociencia y Nanotecnología (INN), Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET)–CNEA, 8400 Bariloche, Argentina



We study the out-of-equilibrium dynamics of a BCS superconductor in the presence of a periodic drive with a frequency ω_d larger than the equilibrium superconducting gap and an external bath providing dissipation. Similar to the dissipationless case, a subset of quasiparticles, resonant with the drive, synchronizes and produces interesting non-linear phenomena. For small dissipation Rabi-Higgs oscillations [Collado et al. Phys. Rev. B 98,  214519  (2018)] can be observed as a transient effect. At long-times in contrast, dissipation leads to a steady state with a population imbalance that increases quadratically for small drives and saturates for large drives as in non-linear quantum optics. We show also that second harmonic generation is allowed for a drive which acts on the BCS coupling constant. We also compute the intensity of time- and angle-resolved photo-emission spectroscopy (tr-ARPES) and time-resolved tunneling spectra. The excited quasiparticle population appears as a decrease (increase) in the photo-emission intensity at energy -ω_d/2 (+ω_d/2) measured from the chemical potential. The tunneling intensity shows a time-dependent depression at ±ω_d/2 due to the lacking of spectral density and population unbalance causing Pauli blockade. We predict that at short times, compared to the energy relaxation time, both experiments will show oscillations with the Rabi-Higgs frequency.  






Spectral fingerprints of the non-linear dynamics of driven superconductors with dissipation
    C. A. Balseiro
    December 30, 2023
===========================================================================================








§ INTRODUCTION




The recent advances in laser technologies have opened new avenues for the study of collective behavior and emergent phenomena in condensed matter <cit.> and ultracold-atom systems <cit.> far from the linear response paradigm. Theoretical examples of these strongly non-linear phenomena are dynamical phase transitions after a quantum quench <cit.> or under a periodic drive <cit.>.

While these studies are of interest per se, they can also, in principle, shine some light on the nature of the precursor equilibrium phases. In this context, superconducting condensates have drawn increasing attention, both in theoretical and experimental research, as challenging cases for the study of collective out of equilibrium states, especially in the case of superconductors with competing orders <cit.>.
Much of the interest in the field was fueled by theoretical studies of quenched systems <cit.>
and experiments in which the superconductor is excited by a pump pulse without a complete suppression of the superconducting state. In the perturbed system, the superconducting order parameter Δ (Higgs mode) or the charge modes evolve with collective oscillations at frequency 2Δ  which are rapidly damped due to dephasing <cit.>. 


Another interesting way to manipulate a condensate is through periodic drives. In solid-state superconductors, different type of drives have been discussed <cit.>. Among them are impulsive stimulated Raman scattering (ISRS) <cit.>, phonon assisted modulation of the density of states (DOS-driving) as  in the case of 2H-NbSe_2 <cit.> or of the coupling constant (λ-driving) as proposed for FeSe <cit.>, microwave drives and THz drives <cit.> as already realized in Ref. [Matsunaga2014]. Periodic drives can also be achieved in ultra-cold atoms where there are well known techniques to modify Hamiltonian parameters at will: DOS-driving can be achieved modifying the depth of a periodic potential, as it has been done for bosons <cit.>, while λ-driven can be implemented in a variety of ways <cit.>.


Previously, we have shown that within a BCS self-consistent dynamics, a periodic excitation at a frequency ω_d in resonance with quasiparticle excitations (ω_d>2Δ) can produce collective Rabi oscillations (“Rabi-Higgs” mode) of the quasiparticle population with a frequency proportional to the strength of the drive <cit.>. This is due to a subset of quasiparticle excitations with energy E_k satisfying ω_d≈ 2 E_k=2√(ξ_k^2+Δ^2) where E_k is the quasiparticle energy and ξ_k is the energy of the fermions measured from the chemical potential. A family of quasiparticles approximately satisfying the above condition synchronize among themselves and perform collective oscillations at the Rabi frequency.

The observation of a Rabi-Higgs mode requires entering into a highly nonlinear regime where damping and decoherence effects can be overcome. It is thus interesting to discuss under which conditions this non-linear regime can be achieved. It is well known from optical Bloch equations for a periodically driven two-level system <cit.> that a finite energy relaxation time τ leads to a linear response regime at long times. The non-linear regime of  Rabi oscillations can be accessed at short times when the  Rabi frequency is larger than 1/τ. Similarly, the observation of Rabi-Higgs oscillations in a superconductor requires long relaxation times and/or strong drives (which result in Rabi-Higgs frequency Ω_R>1/τ). We mention that a similar phenomena has been proposed for driven graphene <cit.> where, however, the physics is much simpler as quasiparticle interactions have not been considered. 


How long can τ be? Quasiparticle relaxation times in superconductors can be extremely long under favorable conditions. A simple estimate can be obtained through  the Dynes parameter <cit.> in tunneling experiments which, as explicitly shown in Ref. [Collado2019], is directly related to the damping of quasiparticles by the bath ultimately limiting the coherent dynamics. In aluminium samples, an inverse Dynes parameter 1/γ≡τ∼ 10^6/Δ∼μs has been measured  <cit.> which suggest that there is a large time window for coherent dynamics. Indeed, in Ref. [Collado2018] we found that for a moderate drive strength the Rabi period τ_R≡ 2π/Ω_R  is in the scale of tenths or hundreds of 1/Δ≪τ (see also Fig. <ref> below). On the other hand, the above estimate for τ is probably too optimistic as the same out-of-equilibrium quasiparticles will open new relaxation channels <cit.> and hence the whole out-of-equilibrium many-body problem can be considered. 


Here we take a simple approach for the interplay of coherent non-linear phenomena and damping by considering a driven superconductor in the presence of a bath providing energy relaxation. We make the simplest possible assumption for the bath leaving a microscopic description for future work. Thus we treat the coupling with the bath as a phenomenological parameter and   
present the spectral signatures  of the Rabi-Higgs modes varying the coupling and other experimental parameters. Differently from previous works, the  dissipative dynamics is treated in a self-consistent way using the method introduced in Ref. [Collado2019]. Both tr-ARPES and tunneling are analyzed in detail. 





§ THE TIME-DEPENDENT BCS MODEL WITH DISSIPATION



In the absence of dissipation, λ-drive and DOS-drive produce qualitatively similar results <cit.>. As we expect the same to be true in the presence of dissipation, we restrict hereon to study λ-driving. 

In addition, in the dissipationless case, even for weak drives and for ω_d>2Δ the Rabi-Higgs mode appears with  frequency Ω_R proportional to the intensity of the drive. Apparently the Rabi-Higgs mode violates any linear response prescription but, as mention in the introduction, without dissipation the system becomes inherently non-linearly at long times no matter how weak the perturbation is. In other words, the order of limits is important. Taking first the Rabi frequency  (proportional to the drive intensity) Ω_R→ 0 and then energy relaxation times τ→∞ a linear response regime is well defined. Inverting the order of limits it is not. In the following we study numerically the effect of dissipation on the Rabi-Higgs response. 




 §.§ Model and formalism

We consider a single-band s-wave superconductor described by the Hamiltonian


    H_BCS=∑_k,σξ_kc_kσ^†c_kσ^-λ(t)∑_k,k^'c_k↑^†c_-k↓^†c_-k^'↓^c_k^'↑^


where c_kσ (c_kσ^†)
destroys (creates) an electron with momentum k, energy ε_k and spin
σ. Here  ξ_k=ε_k-μ measures the energy from the Fermi level μ and the pairing interaction λ(t) is parameterized as


    λ(t)=λ_0[1+Θ(t)αsin(ω_d t)] ,


where Θ(t) is the Heaviside step function. In most of our calculations we take the parameter α∈[0,0.2], that corresponds to a modulation of up to 20% of the equilibrium pairing interaction λ_0, to keep it within the range of experimental accessibility <cit.>. In the thermodynamic limit, the Hamiltonian (<ref>) is equivalent to the mean-field Hamiltonian,

    H_MF=∑_kψ_k^†H_k(t)ψ_k,

written in the Nambu spinor basis ψ_k=(c_k↑,c_-k↓^†)^T, where


    H_k(t)=([     ξ_k   -Δ(t); -Δ(t)^*    -ξ_k ]) ,


and the instantaneous superconducting order parameter is given by


    Δ(t)=λ(t)∑_k⟨ c_k↑^†(t)c_-k↓^†(t)⟩ .


Here ⟨…⟩ denotes the expectation value on the initial state. 

In order to consider dissipation we couple the  superconductor to a reservoir.
For simplicity we will take the bath to be at zero temperature, which means that if the superconductor is out of equilibrium at some time and it is allowed to evolve in the absence of the drive, it will eventually relax to the ground state with the bath absorbing all the excess energy. The method to treat the bath was explained in detail in Ref. [Collado2019]. Here we summarize the main results.
 
To describe the reservoir effect, the self-consistent solution of the gap equation is written in terms of the Keldysh two-time contour Green’s functions. In the Nambu spinor basis, the retarded and lesser Green functions are 2× 2 matrices with matrix elements given by


    G_k^R(t,t^')_αβ   =   -iΘ(t-t^')⟨{ψ_kα(t),ψ_kβ^†(t^')}⟩,
    G_k^<(t,t^')_αβ   =   i⟨ψ_kα^†(t^')ψ_kβ(t)⟩,


respectively. Thus, the superconducting gap equation [Eq. (<ref>)] can be written as


    Δ(t)=-i λ(t) ∑_kG_k^<(t,t)_12 .


When considering the coupling to a reservoir, the lesser Green function satisfies the Keldysh equation in time domain <cit.>


    G_k^<(t,t^')=∫ dt_1∫ dt_2 G_k^R(t,t_1)Σ_k^<(t_1,t_2)G_k^R(t^',t_2)^†


where the dissipation effects are taken into account via the lesser self-energy Σ_k^<(t_1,t_2) and the retarded Green function which is solution of the corresponding Dyson equation with a retarded self-energy Σ_k^R(t_1,t_2).

Following Refs. <cit.> we consider a mechanism for dissipation that couples each pair of states k↑, -k↓ with a reservoir described by a time-independent one body Hamiltonian  H_b=∑_ℓ,σE_ℓd_ℓσ^†d_ℓσ where d_ℓσ^† creates a state in a single particle bath state with energy E_ℓ. In the limit of a wide-band reservoir with identical coupling V_kℓ=V_ℓ for each k  all details of the bath dropout and its effects can be described by a single frequency independent parameter γ describing the effects of inelastic scattering and producing a level broadening ∼γ and a finite lifetime τ=1/γ.
In this case, the retarded and lesser self-energy become momentum independent and diagonal in Nambu space, i.e. Σ_k^R≡ IΣ^R and  Σ_k^<≡ IΣ^<
with


    Σ^R(t_1,t_2)   =   -iγδ(t_1-t_2)/2 ,
    Σ^<(t_1,t_2)   =    iγ∫dω/2πf(ω)e^-iω(t_1-t_2) .


Here f(ω) is the Fermi function evaluated at the bath temperature. Consequently, the Dyson equation for the retarded Green function can be easily integrated, being given by G_k^R(t,t^')=G_k^R (0)(t,t^')e^-γ(t-t^')/2 where G_k^R (0)(t,t^') is the retarded Green function in the absence of dissipation. The latter can be computed by solving the following differential equations (in matrix notation in the Nambu spinor basis and setting ħ=1),


    G_k^R(0)(t,t)   =   -iI,
    
    i∂_tG_k^R(0)(t,t^')   =   H_k(t)G_k^R(0)(t,t^'),       t>t^',
    
    i∂_t^'G_k^R(0)(t,t^')   =-   G_k^R(0)(t,t^')H_k(t^'),       t>t^'.


Replacing Eq. (<ref>) into Eq. (<ref>) and assuming a reservoir at zero temperature we obtain


    G_k^<(t,t^')=-γ/2π∫_-∞^tdt_1∫_-∞^t^'dt_2        G_k^R(0)(t,t_1) G_k^R(0)(t^',t_2)^†
          ×e^-γ(t-t_1+t^'-t_2)/2/t_1-t_2+i0^+ .


Hence, the time dependence of the order parameter [Eq. (<ref>)] can be obtained after computing the lesser Green function [Eq. (<ref>)] for t^'=t. This equal-time lesser Green function G_k^<(t,t)≡G_k^<(t) satisfies the equation of motion


    ∂_tG_k^<(t)=-γG_k^<(t)+ℐ_k(t)-i[H_k(t),G_k^<(t)] ,


where


    ℐ_k(t)=iγ/2π∫_-∞^tdt^'(G_k^R (0)(t,t^')/t-t^'-i0^++G_k^R(0)(t,t^')^†/t-t^'+i0^+)e^-γ(t-t^')/2.


The initial condition for the differential Eq. (<ref>) is given by the equilibrium value of the lesser Green function

    G_k^<(0)=i/2I-i/π E_karctan(2E_k/γ)([  ξ_k -Δ_0; -Δ_0 -ξ_k ]) ,


where E_k=√(ξ_k^2+Δ_0 ^2), which is time-independent and easily obtained after replacing the equilibrium retarded Green function in Eq. (<ref>) (see Ref.[Collado2019]). 

As a consequence, in the presence of dissipation, the equilibrium order parameter Δ_0 is defined, via Eq. (<ref>), by the gap equation


    1=1/π∑_kλ_0/E_karctan(2E_k/γ) .


In the γ→  0 limit Eq. (<ref>) becomes the standard BCS gap equation. In the presence of inelastic scattering (γ) the superconducting order parameter is reduced.  As already mentioned, another important result is that at  equilibrium the present formalism presents a rigorous justification for the Dynes formula for the density of states <cit.>. As will be shows explicitly below, this provides a simple way to estimate the γ parameter close to equilibrium directly from tunneling experiments <cit.>. 




 §.§ Rabi-Higgs modes in the superconducting response and dissipation effects


We now present the numerical solution of our model for λ-driving. For concreteness we shall show simulations for ω_d=4Δ_0 but qualitative similar behavior is obtained for not too large frequencies above the gap (ω_d>2Δ_0).

At t≤0 the system is in equilibrium and the order parameter Δ_0 is given by Eq. (<ref>). At t>0 the drive switches on according to Eq. (<ref>) and we compute the equal-time lesser Green function via  Eqs. (<ref>), (<ref>) and (<ref>) in order to self-consistently determine the superconducting order parameter evolution through Eq. (<ref>).





  §.§.§ Non-linear effects in the undamped dynamics


We first briefly present the superconducting response in the absence of dissipation as a starting point to compare with those in which the reservoir effects play a role. In this case, the calculation can be made using either the Anderson pseudospin language <cit.> or the present formalism by considering the evolution of the equal-time lesser Green function dictated only by the commutator with the  Hamiltonian [Eq. (<ref>) without the two first terms in the r.h.s.]. 

The dynamics of Δ(t) and the expectation value of the momentum distribution function, 


    n_k(t)=∑_σ⟨ c_kσ^†(t)c_kσ(t)⟩ = 1-i [G_k^<(t)_11-G_k^<(t)_22]


are shown in Fig. <ref> for two different values of the perturbation amplitude α.

It is apparent from the figure that the order parameter oscillates with two fundamental frequencies and, after a short transient, averages to a smaller value respect to equilibrium. The drive frequency corresponds to a fast oscillation that can not be resolved on the scale of the figure and leads to the filled black regions of the gap dynamic. In addition, the amplitude of the gap shows the Rabi-Higgs oscillations with a frequency that increases approximately linearly with increasing α. Indeed, 

    Ω_R/Δ_0≈ A(ω_d/Δ_0) α

with A(ω_d/Δ_0)∼ 2.0 for  ω_d/Δ_0=4 [see Eq. (29) in Ref. <cit.> for an analytic approximation].  


The  Rabi-Higgs mode is associated with a periodic inversion of the population of the quasiparticles resonant with the drive. For ω_d=4Δ_0, this is visible in the momentum distribution function n_k(t)  as a narrow time dependent structure at quasiparticle energy  ξ_k≈± 2Δ_0   (Fig. <ref>, middle panel) with a frequency that  matches the Rabi-Higgs period of Δ(t). Notice that the inversion of colors along the anomaly represent a cyclic inversion of population of the resonant quasiparticles.  
Such time-dependent anomaly represent a clear hallmark of the Rabi-Higgs mode and opens the possibility  to detect it through spectroscopies as we shall demonstrate in the next section. 






Another non-linear effect is the generation of a second-harmonic. In the case of an electromagnetic drive, second harmonic generation is not allowed <cit.>. This follows from the general fact that the current response to the vector potential A is J∼ρ_s(A) A where ρ_s(A) is the superfluid stiffness. In the absence of a steady state current, the free energy and  ρ_s(A) are even in A, so the lowest non-linear contribution to J is order A^3, i.e. third harmonic generation <cit.>. In contrast, not such symmetry exist in our case as for λ>0 terms odd in δλ are allowed in the free energy and second harmonic generation is allowed. Indeed, as can be seen in the bottom panel of Fig. <ref> the superconducting response not only contains the driving frequency ω_d=4Δ_0 but also 2ω_d=8Δ_0. 

This can also been seen in the middle panel of Fig. <ref>, where a less intense Rabi-Higgs mode is developed associated with the second harmonic response, i.e. quasiparticles with ξ_k≈± 4Δ_0 also show a narrow time dependent anomaly in the population but with a much smaller Rabi frequency. 








  §.§.§ Dissipative dynamics


The driven superconducting response (Δ(t)) in the presence of a bath is shown in Fig. <ref> for different values of the bath parameter γ and the perturbation amplitude α. As in non-linear optics, in the presence of dissipation the Rabi oscillation is a transient effect. For  sufficiently long times a steady state is achieved where only the drive frequency is present.

 As γ is increased (from top to bottom in Fig. <ref>) the system reaches a steady state more rapidly with a vanishing of Rabi-Higgs oscillations. By increasing α, for a fixed value of γ, Ω_R increases [cf. Eq. (<ref>)] and several Rabi-Higgs oscillations are visible before they disappear as a consequence of relaxation (see for example Fig. <ref> (c,d)). Thus, in order to detect the Rabi-Higgs modes experimentally it is necessary to ensure that Ω_R ≳γ.  In fact, for the γ values used here, the slower Rabi-Higgs mode, associated with the second harmonic generation, is not visible.

Conversely, for  small α and strong dissipation it is possible to get Ω_R≲γ and only the oscillations synchronous with the drive are visible as expected from  linear response theory (bottom panels of Fig. <ref> and Fig. <ref>). Indeed, in this regime, the amplitude of the oscillation in the order parameter increases linearly with α as  shown in Fig. <ref>. However, it is important to note that by increasing α, and after a very fast transient, the superconducting gap decreases in average which constitutes the first nonlinear effect arising in the dynamics. As a conclusion, in the presence of dissipation, we can distinguish two different regimes: one in which the linear response theory is applicable even at short times (Ω_R≲γ) and the other corresponding with a transient strong nonlinear behavior in which Rabi oscillations can be observed (Ω_R≳γ).






We now discuss in more detail the nonlinear regime, in particular the possibility of detection of the Rabi-Higgs mode in experiments. Previously, we have demonstrated that the Rabi-Higgs mode is associated with oscillations in the occupation values n_k (see Fig. <ref>). This charge fluctuations provide an efficient manner to detect the existence of this nonlinear mode with standard experimental techniques as we shall show in the next section. 

Since for the simple electronic structure we are taking, n_k depends on k only through its distance from the Fermi surface (ξ_k) it is useful to introduce the distribution function n(ξ_k)≡ n_k. When convenient, in the following we will drop the momentum dependence and refer to n(ξ) loosely as the momentum distribution function keeping in mind the above equivalence. 

In Fig. <ref> we show the time-dependent n_k distribution taking into account dissipative effects (middle and right column) in comparison with the dissipationless counterpart (left column) at different times within the first Rabi cycle (panels (d)-(f)) and in the steady state (panels (a)-(c)). We have used the same γ and α parameters as in the right column of Fig. <ref>. 

For γ=0 the Rabi oscillations can be observed as oscillations in the occupation value n_k for ξ_k≈± 2Δ_0, which we will refer to as  n(±ξ_R). Small peaks also appear for ξ_k≈± 4Δ_0, corresponding to the second-harmonic  slower Rabi-Higgs mode that starts to develop in the temporal window used in Fig. <ref>(d). At long-times, the two Rabi-Higgs modes can be seen as is shown in Fig. <ref>(a). 

In the presence of weak dissipation (γ=0.05Δ_0), the transient dynamic only shows oscillations of  n(±ξ_R), corresponding to the fastest Rabi-Higgs mode. Finally, 
for strong dissipation (γ=0.2Δ_0), a full Rabi cycle cannot be completed since the relaxation takes place in a very short time. As a consequence n(ξ_R) (n(-ξ_R)) increases (decreases) during a short period of time and rapidly saturates without exhibiting Rabi oscillations.






The  n_k distribution in the steady-state in the presence of dissipation are shown in Fig. <ref>(b) and Fig. <ref>(c) for different times within a period of the drive. In this case, n_k shows a peak at ξ_k≈± 2Δ_0 such that n(±ξ_R)≈ 1 all the time. Taking into account spin this correspond to half-filled single particle states which, as shown below, corresponds to saturation of the population imbalance due to strong drive. Thus, ultimately the fate of the superconducting dynamics is a gap oscillation with the drive frequency (see Fig. <ref> (e,f)) and a stationary population unbalance at the quasiparticle energy ξ_R. 

The population unbalance changes considerably as a function of the amplitude of the perturbation as can be appreciated in Fig. <ref>. Because the Rabi oscillation is only a transient effect, a full population inversion (e.g. a sign change of n(ξ_R)-n(-ξ_R)) is not possible in the steady state. As a consequence, the excited state population can not exceed the half-filled value (counting spin) n(ξ_R)≈1 by considering a bath in the wide-band limit <cit.>.

We have extracted the n_k distribution in the steady state for several values of α with γ=0.2Δ_0. In accordance with the above statement, the value of n(ξ_R) (n(-ξ_R)) in the steady state, increases (decreases) with α and saturates to 1 in the limit of large α as it is shown in Fig. <ref>. Notice that the population imbalance is quadratic in α so it vanishes in the linear response regime which follows also from the fact that, by symmetry, it should be even in α.  

This saturation of the excited population is a common effect in nonlinear quantum optics. There, one usually considers a two-level system in the presence of a driving force and damping effects via a phenomenological parameter in the Bloch equation of motion. For a drive frequency in resonance with the two level system, the excited population in the steady state increases with the amplitude of the perturbation and in the limit of large intensity, the largest possible excited population is equal to the ground state population <cit.>. Again, the difference in the present context is that the phenomena does not corresponds to a single two-level system but a collection of two-level systems  interacting through the self-consistency that determines the superconducting order parameter. 




 





§ THEORETICAL MODELLING OF TR-ARPES AND TUNNELING EXPERIMENTS: SPECTRAL FINGERPRINTS OF RABI-HIGGS MODE


Our main aim in this section is to identify some spectral fingerprints of the Rabi-Higgs mode, and related non-linear phenomena, that could be experimentally detected in the presence of  dissipation. 
Clearly, an emergent technique to detect time dependent phenomena is tr-ARPES and so we discuss such a case first. Yet,  in addition, we shall demonstrate that also time resolved tunneling experiments could be useful to detect the non-linear mode discussed above.





In all the calculations discussed below, we obtain the spectral signals in terms of the lesser Green function Eq. (<ref>) as a function of α and γ. 



 §.§ tr-ARPES


In our setting the  λ-driving is turned on at time t=0 and the photoemission process is induced by a wave packet of photons centered at time t_0 and with central energy ħω_ q larger than the work function of the solid W. For simplicity we use a Gaussian shape for this probe pulse, s(t)=exp(-(t-t_0)^2/2σ^2)/(σ√(2π)) with standard deviation σ. 

In a tr-ARPES experiment the momentum of the outgoing electrons k_e is measured. Energy conservation determines the excitation energy left in the system after the photoemission process,
ħω=ħω_ q-(ħ k_e)^2/(2m_e)-W and momentum conservation yields information on the momentum of the excitations k. The momentum resolved photocurrent in the detector at time t is due to all electrons photoemitted before that time and is determined by
<cit.>,


    I_k(ω,t)=Im∫_-∞^tdt_1∫_-∞^tdt_2s(t_1)s(t_2)e^iω(t_1-t_2)G_k^<(t_1,t_2)_11 .



In order to probe the Rabi-Higgs modes in real time in the following we use 
different probe times t_0 and take the integration limit from a lower cutoff time t_l=t_0-5σ to t=t_0+5σ. Fig. <ref> shows the tr-ARPES intensity  with and without dissipation for different t_0 during the first period τ_R associated with the Rabi-Higgs mode. For simplicity, we consider a parabolic band for electrons ξ_k∝ (k^2-k_F^2) and use a probe pulse with a standard deviation σ=0.1τ_R. From the photo-emission signal at equilibrium (left column in Fig. <ref>) there is an almost imperceptible broadening of the  spectral line as γ is increased, which leads to a decrease of superconducting order parameter according to Eq. (<ref>).  


The presence of Rabi-Higgs oscillations are clearly visible in the top and middle panels of Fig. <ref> via a spectral weight around ω=±2Δ_0 that increases  (decreases) above (below) the Fermi energy in the first half-period and has the opposite behavior in the second 
half-period of the Rabi-Higgs mode. One can visualize the process as an excitation of quasiparticles from the lower quasiparticle branch to the higher quasiparticle branch in the first half-cycle followed by a deexcitation in the second half-cycle or equivalently as a stimulated absorption phase followed by a stimulated emission phase. 

If γ is large enough, Rabi oscillations are not visible and only the spectral weight corresponding with the steady state is observed after a fast transient (bottom panel of Fig. <ref>). 

Clearly the above dynamics is the photoemission image of the momentum distribution function imbalance discussed above. Indeed, integrating in frequency 

    ∫_-∞^∞dω/2πI_k(ω,t)   =   Im∫_-∞^tdt^'s(t^')^2G_k^<(t^',t^')_11
       =   1/2∫_-∞^tdt^'s(t^')^2 n_k(t')

which is clearly a moving average of  n_k(t). For example, in the large γ case the 
unbalance over time in the tr-ARPES intensity for ω=±2Δ_0 matches the n_k distribution in the steady state shown in Fig. <ref>(d).

These results establish  time-resolved photo-emission experiments as a tool to investigate steady state non-linearities, like population imbalance, and dynamic non-linearities like the Rabi-Higgs mode, in driven superconductors. 



  §.§.§ Floquet Analysis


As already mentioned for t≳τ=1/γ one reaches a steady state and the superconducting gap oscillates only with the drive frequency. Therefore the mean-field Hamiltonian  Eq. (<ref>) is time periodic at long times and we can use Floquet theorem to analyze the spectrum. This theorem guarantees the existence of a set of solutions of the time-dependent Schrodinger equation of the form

    |ψ_ν(t)⟩=exp(-iε_νt/ħ)|ϕ_ν(t)⟩

where |ϕ_ν(t)⟩ has the same periodicity of the Hamiltonian Eq. (<ref>) <cit.>. The Floquet states |ϕ_ν(t)⟩ are the solutions of 

    ℋ_F|ϕ_ν(t)⟩=ε_ν|ϕ_ν(t)⟩

where ℋ_F=H_MF-iħ∂/∂ t is the Floquet Hamiltonian and ε_ν is the quasienergy <cit.>. In the Floquet basis this Hamiltonian becomes a time-independent infinite matrix operator. Since we are interested in the low energy spectra we restrict the Floquet Hamiltonian to a large enough but finite subspace containing many multiphoton processes (finite number of replicas) <cit.>. 

For γ=0.2Δ_0 and α=0.2 the steady state condition is reached already for t∼ 0.75 τ_R= 11.25/Δ_0. The quasi-energy spectrum in the long time limit is well described by  Δ(t)=Δ̅+Δ_1 cos(ω_d t) with Δ̅=0.86Δ_0 and Δ_1=0.42Δ_0  [c.f. Fig. <ref>(f)]. 


In the upper panels of Fig <ref> we compare the tr-ARPES with the Floquet spectrum in the steady state. To analyze the details in the right panels we sacrifice temporal resolution to gain energy resolution by using a wider probe pulse. We see that the tr-ARPES signal nicely matches the Floquet spectrum, depicted by the dashed line. Thus the tr-ARPES signal essentially probes the occupied parts of the Floquet band structure.

Surprisingly, also in the transient dynamics the tr-ARPES intensity fits very well with a Floquet spectrum that is obtained from an effective Δ(t) with a monochromatic dependence (see lower panel of Fig. <ref>), even thought in this regime the superconducting gap shows several incommensurate frequencies (Rabi-Higgs and drive frequency) and Floquet theorem is not strictly applicable. For γ=0, we compute the Hamiltonian spectrum for a Δ(t)=0.63Δ_0+0.35Δ_0cos(ω_d t). Thus the transient response shown in Fig. <ref> and the lower panels in Fig. <ref> can be seen as a  time-dependent change in the occupancy of the Floquet band structure. Clearly, the reason this analysis works in the Rabi-Higgs oscillation regime is the large separation of time scales between the slow Rabi-Higgs dynamics and the fast drive oscillations <cit.>. 









 §.§ Time-resolved tunneling experiment


We now discuss a possible setup to detect the Rabi-Higgs mode via tunneling measurements. In the last decades there have been efforts to add temporal resolution to scanning tunneling spectroscopy (STM) techniques and now sub-picosenconds resolution can be reached in experimental setups <cit.>. In principle one can radiate a sample while performing the measurement  making STM a feasible time resolved technique at the same level that pump-probe measurements <cit.>. However, since the spatial resolution is not a requirement to reveal this nonlinear mode, a planar junction would probably be a more stable setting. An antenna could be used to couple the radiation with the superconductors as is done in detectors for photons with energy matching the superconducting gap <cit.>. In what follows, we will refer to STM signal but our formalism applies also for the case of a planar junction.

To obtain the STM signal, we now consider a tip located close above the out of equilibrium superconductor and tunnel-coupled locally to the k states in the system. The entire problem can be described by a generic Hamiltonian H=H_MF+H_T+H_TS where we have added the tip (or metallic contact) Hamiltonian 


    H_T=∑_p(ϵ_p+eV)a_p^†a_p^ ,


and the coupling between both subsystem via the tunneling Hamiltonian


    H_TS=∑_k,p(T_kp c_k^†a_p^+h.c.).


The tip or the metallic contact is connected to an external voltage source with energy eV, a_p (a_p^†) destroys (creates) an electron with momentum p and energy ϵ_p+eV in the tip and T_kp is the strength of the tunneling. The electron current operator is given by


    Î=ie[N̂,Ĥ]=ie(L̂^†-L̂) ,


where e is the electron charge, N̂=∑_pa_p^†a_p is the number operator and L̂=∑_k,pT_kpc_k^†a_p^. Assuming a weak coupling between tip and the superconductor we calculate the current through the tip in the linear response limit as


    I(t)=ie∫_-∞^∞θ(t-t^')⟨[Î(t),H_TS(t^')]⟩_0


where ⟨ ... ⟩_0 denotes the expectation value at zero order in the tunneling Hamiltonian. Notice that we have not consider the spin degree of freedom so far as the current is a spin conserved quantity (not spin-flip are allowed in the tunneling process). This will be included in the final expression of the time-dependent current as a factor 2. From Eq. (<ref>) we obtain 



    I(t)=2eT^2∫_-∞^tdt^'∑_k,p([e^-iϵ_p(t-t^')⟨ c_k^†(t)c_k^(t^')⟩Θ(ϵ_p-eV)-e^iϵ_p(t-t^')⟨ c_k^(t)c_k^†(t^')⟩Θ(-ϵ_p+eV)]+(t↔ t^')) ,



where we have assumed a momentum-independent tunneling coupling T_kp≡ T and a zero-temperature Fermi distribution for the electrons on the tip via the Heaviside step function Θ(x). Taking the derivative of Eq. (<ref>) with respect to the voltage V and using Eq. (<ref>), the time-dependent differential conductance can be written as


    G(t)=dI/dV∝Im∑_k∫_-∞^tdt^'Tr[G_k^<(t,t^')]e^-ieV(t-t^')


where Tr represents the trace in Nambu space.

In the absence of the drive (that is, at equilibrium) the differential conductance becomes time-independent, being proportional to the well-know phenomenological Dynes density of states


    G∝ρ(eV)=ρ_0 Re [ eV+i Γ/√((eV+i Γ)^2-Δ_0^2) ] ,


where the Dynes parameter Γ=γ/2, Δ_0 is the equilibrium order parameter and ρ_0 is the normal phase density of states. It is worth recalling that the Dynes formula was originally introduced phenomenologicaly <cit.>. There have been several theoretical proposals to provide a formal justification including Eliashberg physics <cit.>, inelastic tunneling <cit.> and magnetic impurities <cit.>. In our previous work we have shown that the coupling with a bath <cit.> provides a mechanism to justify Dynes formula and allows to link directly equilibrium tunneling with the γ parameter. Of course since we are not providing a microscopic theory of the bath our justification is still semi-phenomenological. 


Turning to the driven case, for small γ values the out-of-equilibrium differential conductance is computationally very demanding since the whole dynamics has to be integrated in Eq. (<ref>), differently from the photoemission case where the integrals are cutoff by the probe pulse shape. To make the computations feasible we use γ=0.2Δ_0 and α=0.1. The nominal Rabi-Higgs period corresponds to  τ_RΔ_0≈27 but because of the large damping  a steady state is reached before a full Rabi cycle is completed. Still as shown in Fig. <ref>, the time-dependent differential conductance clearly shows a non trivial transient dynamics. Indeed, a clear depression in the time-resolved conductance around eV=2Δ_0 can be observed as a function of time which can be identified with  the beginning of the Rabi-Higgs oscillation. 

As in the case of photoemission two effects determine the features in the spectrum. The upper and lower Bogoliubov bands hybridize with the first Floquet sideband of the  lower and upper Bogoliubov bands respectively creating a  pseudogap at ±ω_d/2 (see right column of Fig <ref>). Furthermore, the depression of the population at energy -ω_d/2 depresses the probability to extract electrons and the excess of population at energy ω_d/2 prevents injecting electrons due to Pauli blocking. 

Another effect of the drive is that whole shape of the conductance gets modified. Due to the Dynes parameter, at equilibrium the differential conductance does not vanishes sharply for energies lower than Δ_0. Once the time-dependent perturbation is turned on, the superconducting coherence peak decreases according to a decreasing in the average  order parameter. 











§ SUMMARY AND OUTLOOKS




We have studied the dynamic of a BCS superconductor subject to a periodic drive in the presence of dissipation focusing on interesting transient dynamics. For γ≲ω_R the transient show   Rabi-Higgs oscillations which persist forever when γ=0. In the opposite limit γ≳ω_R oscillations are not observed and the response is practically described by linear response at all times for a weak drive and show the phenomena of saturation of population imbalance for large drives. Another interesting non-linearity is second harmonic generation, which we have found it is allowed for λ-drive. 

In all regimes the behavior resembles an isolated  driven two level system. However one should keep in mind that the interactions are fundamental in synchronizing a finite fraction of the quasiparticles and make the mode collective. There are close analogies with phenomena in quantum optics where  Rabi oscillations and saturation of population are common phenomena. This analogy paves the way to explore quantum optics protocols to control and manipulate the superconducting state of materials as already proposed in Ref. <cit.>.

The experimental detection of these highly nonlinear behaviors in driven superconductors could be a major step towards the quantum control and manipulation of quantum phases. We have proposed two techniques to detect  Rabi-Higgs oscillations  taking into account dissipation. A cyclic imbalance of quasiparticle population at specific  locations in energy and momentum has been shown to appear in  tr-ARPES. Due to technical reasons only the rise of the population imbalance could be explicitly computed for tunneling experiments but we expect that, in close analogy with photoemission,  
the anomalies oscillate with the Rabi frequency for small dissipation. Even detection of the steady state non-linearities would be a quite interesting experimental achievement. 

The present formalism can be easily extended to take into account more interesting relaxation mechanisms via more sophisticated self-energies beyond the wide-band approximation used here for the bath. An interesting direction is to analyze how robust the Rabi-Higgs mode is in the presence of dephasing, decoherence and relaxation sources from a more microscopic point of view by considering residual Coulomb and electron-phonon interactions where heating effects could be relevant.



J.L. is in debt for enlightening discussions with L.  Benfatto, C. Castellani, G. Seibold, B. Leridon, N. Bergeal and Jérôme Lesueur. We acknowledge financial support from Italian MAECI and Argentinian MINCYT through bilateral project AR17MO7 and from ANPCyT (grants PICTs 2013-1045 and 2016-0791), from CONICET (grant PIP 11220150100506) and from SeCyT-UNCuyo (grant 06/C603). J.L. acknowledges financial support from  Italian MAECI thought collaborative project SUPERTOP-PGR04879, from Italian MIUR though 
Project No. PRIN 2017Z8TS5B, and from Regione Lazio (L. R. 13/08) under project SIMAP.






apsrev4-1_title





