Transient dynamics of a magnetic impurity coupled to superconducting electrodes: exact numerics versus perturbation theory

Impurities coupled to superconductors offer a controlled platform to understand the interplay between superconductivity, many-body interactions, and non-equilibrium physics. In the equilibrium situation, local interactions at the impurity induce a transition between the spin-singlet to the spin-doublet ground state, resulting in a supercurrent sign reversal ($0-\pi$ transition). In this work, we apply the exact time-dependent density matrix renormalization group method to simulate the transient dynamics of such superconducting systems. We also use a perturbative approximation to analyze their properties at longer times. These two methods agree for a wide range of parameters. In a phase-biased situation, the system gets trapped in a metastable state characterized by a lower supercurrent compared to the equilibrium case. We show that local Coulomb interactions do not provide an effective relaxation mechanism for the initially trapped quasiparticles. In contrast, other relaxation mechanisms, such as coupling to a third normal lead, make the impurity spin relax for parameter values corresponding to the equilibrium $0$ phase. For parameters corresponding to the equilibrium $\pi$ phase the impurity converges to a spin-polarized stationary state. Similar qualitative behavior is found for a voltage-biased junction, which provides an effective relaxation mechanism for the trapped quasiparticles in the junction.


I. INTRODUCTION
Superconductors are macroscopic coherent materials, described by a complex order parameter ∆, measuring the pairing amplitude between electrons. They are characterized by a non-dissipative electric current through the material for sufficiently low bias voltages [1]. In conventional superconductors, the electron pairing induces a gap in the density of states at the Fermi level with a value |∆|. However, a spatial modulations of ∆ can give rise to discrete excitation energies inside the gap due to the multiple Andreev reflections [2][3][4][5]. These states have been observed in superconductors coupled to other materials, including semiconductors, normal metals, and other superconductors, for recent reviews see [6,7].
The coupling between superconductors and magnetic impurities or small quantum dots offers the possibility to investigate the competition between superconductivity and many-body interactions [8][9][10][11][12][13][14][15][16]. Coulomb blockade tends to fix the number of electrons in the system, suppressing charge fluctuations and blocking the current through the system [17]. In the strong Coulomb blockade limit, the ground state of an impurity with an odd number of electrons is doubly degenerate, hosting an unpaired spin [9]. In this regime, spin fluctuations increase the low-bias conductance, known as the Kondo effect [18,19], characterized by the Kondo temperature, T K .
The time dynamics of impurities coupled to superconductors can answer fundamental questions, like the time onset of electron correlations, formation of subgap states, dependence on the initial conditions, and relaxation of quasiparticles. These questions are inaccessible from stationary calculations. They are essential for controlling the internal state of the system, relevant for quantum information applications [55], and to study the interplay between superconductivity and many-body interactions under non-equilibrium conditions, including the possible excitation of the Higgs mode [56].
In the absence of electron-electron interactions, the transient dynamics of a superconducting nanojunction has been studied in Refs. [50,[57][58][59][60][61][62][63]. In this case, the system is not able to thermalize due to the trapping of quasiparticles in the subgap states [58,59,64]. These quasiparticles are long-lived, being an obstacle towards the efficient control and manipulation of superconducting devices [65,66]. The effect of electron correlations has been much less studied in the non-equilibrium case [11,[67][68][69][70][71] and the transient regime [72,73]. In this work, we develop an exact numerical method that simulates the real-time evolution of the system, based on the time-dependent density matrix renormalization group (tDMRG). This method is a powerful numerical technique to solve the timedependent Schrödinger equation in finite systems, accounting for many-body interactions [74][75][76][77][78][79]. To study the long-time properties, we use a second-order perturbation expansion on the Coulomb strength, describing the system dynamics for a wide range of parameters.
We show that local interactions do not provide an effective relaxation mechanism for the trapped quasiparticles in the subgap states. In contrast, the inclusion of a third normal lead or a finite bias allows for quasiparticle relaxation due to transitions to excited states with energies above |∆| [80]. In this last case, the superconducting phase difference evolves linearly in time, giving rise to a time oscillating current due to multiple Andreev reflections (MAR) [59,[81][82][83][84][85][86][87][88][89][90]. We show that the impurity converges to a state with spin-rotation symmetry for parameters in the 0 phase. In contrast, it converges to a spin-polarized state for parameters in the π phase and an initially trapped spin in the impurity.
The rest of the manuscript is organized as follows: In Sec. II we introduce the model and the two methods used in our study: tDMRG and perturbative Green functions. In Sec. III we present the main results of our work. We first focus on the equilibrium situation, Sec. III A, where we compare both methods and analyze the convergence to the steady-state. In Sec. III B we analyze the question of thermalization in the case where a metallic lead couples weakly to the impurity. The voltage-biased situation is discussed in Sec. III C, where we show that it provides an effective relaxation mechanism for the initially trapped quasiparticles. Additional results complementing the ones shown in the main text can be found in the Appendix. Finally, we present the main conclusions and a brief outlook in Sec. IV.

A. Model
We consider an impurity coupled to superconducting electrodes, described by the Hamiltonian The impurity is described by a spin-degenerate level where Ψ 0 = (d 0↑ , d † 0↓ ) T is the Nambu spinor, with d 0σ being the annihilation operator for an electron with spin σ at the impurity, andτ are the Pauli matrices in the Nambu space. Here, ε 0 is the impurity energy level, U the charging energy, and n σ = d † 0σ d 0σ , with σ =↑, ↓. The leads are described as a 1-dimensional chain of electronic sites, given by where t 0 is the inter-site chain hopping parameter, h νj = ξ ντz + ∆ ντx , with ξ being the onsite energy of electron at site j, and ∆ ν the superconducting gap with ν = L, R.
The hopping between the impurity and the leads is described by Here,V ν (t) = θ(t)V ντz e iτzφν (t)/2 , where we have considered the tunneling amplitude, |V ν |, to be momentumindependent. The tunneling rate to the leads is given by Γ ν = π|V ν | 2 ρ ν , where ρ ν = 1/πt 0 is the normal density of states at the Fermi level, and Γ = Γ L + Γ R . We consider the impurity to be suddenly connected to the leads at where µ ν is the chemical potential of ν electrode. In the following we consider a symmetric voltage drop in the junction, µ L = −µ R = V /2 with V being the applied bias voltage.
In this work we focus on the time-evolution of single particle observables. The spin-polarization of the impurity is given by is the occupation per spin and n 0 (t) = n ↑ (t) + n ↓ (t) the impurity charge. The current at ν interface is given by where the trace is taken in the Nambu space. We define the symmetrized current as I = (I L − I R ) /2. Finally, the electron pair amplitude at the impurity is given by B. Time-dependent density matrix renormalization group The method tDMRG is a powerful numerical technique to solve the time-dependent Schrödinger equation in finite systems with a great degree of control and numerical precision. In the same way as in its ground state version, tDMRG relies on a matrix product state representation of the wave function plus a Suzuki-Trotter decomposition of the evolution operator to evolve the state [74][75][76][77][78][79]. Among other applications, it has been used to study time-dependent transport in a number of one-dimensional setups [91][92][93][94][95][96][97][98][99][100][101]. The formulation is extensively discussed in the cited literature, so we limit ourselves to a brief summary of the protocol used in this work. Typically, in all the situations discussed here, we are interested in the evolution of the the system after it is perturbed at time t = 0 by some external potential (a quench of the tunneling amplitudes, for instance). We usually observe that the wave function develops a wavefront that propagates away from the center, translating into a transient regime. After some characteristic time that depends on the model parameters, observables such as the current reach a pseudo-steady-state, equilibrating around a constant value. Eventually, the wave front hits the boundaries of the open system, the wave packet bounces back, and the current is reversed. In order to extend the duration of the steady-state time window, we delay the wave packet at the edges by introducing "damped boundary conditions" with exponentially decaying hopping amplitudes near the boundaries of the system, slowing down the electrons near the edges. [102,103]. We use a hoping that decays exponentially with the distance to the impurity as t j = t 0 Λ −j/2 between j and j + 1 sites in the left and right leads. For our calculations, we use a chain with N = 120 sites and Λ = 20, which we found optimal for convergence.
In this work, we extend the tDMRG method to the case where the leads are superconducting, including a local pairing term ∆d n↑ d n↓ +H.c. at every site of the chain, except at the impurity. A particular consideration to take into account about the included pairing term is that the particle number is conserved modulo two [104] (i.e., number parity and spin are good quantum numbers). To preserve unitarity and prevent loss of accuracy, we pick a bond dimension (the dimension of the DMRG basis) large enough to keep the truncation error under a predetermined tolerance, and we stop the simulations either when the accuracy grows beyond this value, or when the effects of the boundaries become dominant. We typically use a bond dimension of 200 and a tolerance of 10 −6 . The observables we represent are given by Eqs. (5-7) for the impurity population, current, and pair amplitude.
Our protocol is implemented as follows: before calculating the evolution, we prepare the system in its global ground state. In this work, we focus on the situation where the impurity is initially detached from the leads. The initial impurity occupation is controlled by the initial level energy, ε 0 (t = 0 − ). We also introduce an initial Zeeman splitting in the impurity Hamiltonian (2) of the form At the initial time, t = 0, the impurity is connected to the superconducting leads and the impurity Hamiltonian takes the form of Eq. (2), with h z = 0 and ε 0 assuming its final value. We then evolve the wave function using tDMRG. Finally, we include the bias voltage as a timedependent hoping between the impurity and the leads, using the expression in Eq. (4).

C. Non-equilibrium Green functions
Perturbation theory on the interaction strength, U , provides accurate results in the U/Γ 1 regime. This approximation has been shown to describe accurately physical observables in the stationary regime and the nonmagnetic phase [22,53,54]. In this work, we extend this approximation to the time domain, allowing us to describe the non-stationary and non-equilibrium situations. For this task, we extend to the superconducting case, the self-consistent formalism introduced in Refs. [52,105], describing the formation of the Kondo resonance at the Fermi level for metallic electrodes. The method is based on a time discretization of the non-equilibrium Green functions (NEGFs) [106][107][108][109]. The impurity Green function is given bŷ We solve the problem in a discrete time mesh, where the inverse of the isolated impurity Green function is given byĝ The electron and hole parts of the inverse Green function in the discrete mesh are described in Refs. [58,59]. The leads self-energy is given by [15] Σ αβ where α, β = ± denote the Keldysh branch. The leads boundary Green functions are given bŷ The bare Green functions for a superconducting chain are given by and with W = 4t 2 0 + ∆ 2 being the bandwidth and 0 + an infinitesimal. The Keldysh components are given by To calculate the interaction self-energy, we perform a perturbation in U/Γ. Up to second order in the perturbative parameter, the self-energy is given byΣ int = Σ HF +Σ (2) . The linear term in U is the mean-field approximation, given by where n = 1, 2 denote the Nambu component andn the opposite component to n. Here, n nn denotes the charge per spin, given by n ↑ = n 11 (t) = −iG +− 11 (t, t) and n ↓ = n 22 (t) = −iG −+ 22 (t, t). The other two components, off-diagonal in Nambu space, describe the induced pairing at the impurity at the mean field level, n nn (t) = −iG +− nn (t, t). The pairing amplitude in the impurity, ∆ i is given by the n 12 (t) component.
Electron correlation effects are described by higher order terms in the expansion. They contain information about the electron-hole pair creation due to the local Coulomb repulsion. Up to second order in U/Γ, the selfenergy is given by where we have perturbed over the mean-field Green function, given bŷ This approximation has been shown to provide accurate results in the normal case, describing the onset of the Kondo peak at temperatures lower than T K [51,52,110,111].
The system observables can be obtained from the Green's function (9). The average charge and spin at the impurity are given by n 0 = n ↑ +n ↓ and S Z = (n ↑ −n ↓ )/2. The current at a given interface can be calculated as [22] where the trace is taken over the Nambu space.

A. Phase-biased junction
In the phase-biased situation (V = 0), a supercurrent can flow through the junction if the two leads have different superconducting phases. In the stationary regime, the current at zero temperature is given by the phase derivative of the ground state energy as a function of the phase difference [112]. In Fig. 1 we show results for the current and the spin after a sudden connection of the impurity to the leads. We compare results for the two methods developed in our work, tDMRG and perturbative NEGFs, for an initial spin in the impurity, i.e. (n ↑ (0), n ↓ (0)) = (0, 1). In the absence of interactions (U = 0), both methods provide similar results. The small differences found, mostly visible at long times, are due to the renormalization of the electron pairing amplitude in the leads close to the impurity. This effect is taken into account within tDMRG but absent in the perturbative NEGFs approach. Both methods agree for a relatively wide range of U/Γ values. In Fig. 1 we show results for S Z and the current through the system. In App. B 1 we also display the time evolution of the pairing potential at the impurity. The results of both methods start to deviate for U/Γ ≈ 6, predicting different stationary values for S Z (not shown).
In the limit ∆ T K , the two methods describe the trapping of a magnetic moment in the impurity at long times [58,59]. We have checked that no spin relaxation happens for times t 50/∆ (the steady-state value has been reached at t ∼ 10/∆). This is an indication that localized Coulomb interactions do not provide an effective relaxation mechanism for the initially trapped spin in the impurity. On the contrary, localized Coulomb interactions tend to quench the short-time spin relaxation, leading to a higher S Z stationary value when increasing U , Fig. 1 (a). The absence of thermalization causes a supercurrent reduction with respect to the expected value at equilibrium [58].
In the case of a fully occupied or empty impurity level, the charge exhibits undamped oscillations [59,60], shown in Fig. 2. The oscillation frequency is given by the energy of the subgap states formed at the interface between the impurity and the superconductors (analytic expressions for the oscillations are given in Ref. [59] for U = 0). The amplitude of the oscillations depend on the interaction strength, which tends to decrease its height. This is due to the increased energy difference between the ground state, with n 0 = 1, and the excited states with two and zero electrons. In the limit U/∆ 1, the oscillations quenched and n 0 converges to the expected stationary state.

B. Thermalization
The results in Fig. 1 and 2 suggest the trapping of long-lived quasiparticles in the system. Their effects are better illustrated in the long-time current, shown in Fig.  3 for the exact tDMRG method. The supercurrent converges to a lower steady-state, Fig. 3 (a), compared to the equilibrium result, Fig. 3 (b). The difference is better observed before the supercurrent sign reversal at U c /Γ ∼ 3.7. The lower supercurrent value is due to the trapping of quasiparticles in the subgap states, which have an infinite lifetime within this model. Therefore, the long-time state is dominated by the short-time dynamics, when the Andreev bound states are formed [58].
The absence of thermalization is consistent with the observed long relaxation times for trapped quasiparticles observed in experiments [65,66]. The finite relaxation time in these systems could be due to photon or phonon emission/absorption processes [113] which are not included in the present analysis. In addition, a small normal density of states inside the gap can provide a finite quasiparticle lifetime. Finally, transitions to the continuum states forced by an applied bias can allow quasiparticles to diffuse into the bulk of the superconductor, making the junction relax. The obtained results suggest that local Coulomb interactions do not provide by themselves an effective relaxation mechanism for trapped quasiparticles in the junction. This might not be the case if one considers extended interactions in the leads [114].
To simulate the mechanisms commented above for quasiparticle relaxation, we attach an additional normal lead to the impurity [58,60,61], Fig. 4. In this case, the supercurrent exhibits a sign change when increasing U/Γ, characteristic of the 0−π transition, Fig. 3 (b). We find a smaller stationary value with respect to the exact calculations in Fig. 3 (b), due to the finite tunneling coupling to the normal electrode, η. The relaxation rate is given by η when the system is far away from the 0 − π transition. In contrast, the system relaxation close to the transition (U c /Γ ≈ 3.7 for this choice of parameters) becomes slower, with a relaxation time exceeding 1/η. This behavior is more clearly seen in panel (a), where we represent S Z . It indicates that perturbations take longer to relax close to the transition. This phenomenology is similar to the critical slowing down close to a phase transition. However, in our case, the 0−π transition becomes smooth due to the finite η chosen. Therefore, it is not possible to observe the expected divergence of the relaxation timescale at the transition.
The additional relaxation mechanism considered makes the spin in the impurity converge to different stationary values, depending on the spin of the ground state, Fig. 4 (a). For the smallest U/Γ values shown, S Z converges to zero, consistent with expected singlet ground state. The situation is different in the π phase, where S Z converges to a non-zero value for an initial condition with broken spin-rotation symmetry. Similar stationary results have been obtained for different η values. In the next subsection, we show that a finite bias voltage provides a similar qualitative behavior, which acts as an effective relaxation mechanism for quasiparticles trapped in the subgap states.

C. Voltage-biased junction
A bias voltage applied to the junction produces a linear increase of superconducting phase difference in time as φ(t) = φ(0) + V t. For V < 2∆, multiple Andreev reflections lead to the transference of multiple charges between the superconductors, resulting in a time-oscillating as well as a dc steady-state current. In the low-bias situation (V ∆, Γ), the subgap states can adapt instantaneously to the phase change and an adiabatic picture is valid. The situation is more complex in the case where V is comparable to ∆ and Γ. In this situation, non-adiabatic effects are important and higher harmonics have to be considered [84].
In Fig. 5 we compare the results from tDMRG and the perturbative NEGFs for V = ∆. Similarly to the phasebiased situation studied before, both methods agree for U/Γ 1. They start to deviate for the largest U/Γ values shown, where the differences are more evident in the current evolution. We show additional calculations for different V values in App. B 2. We additionally compare the pairing amplitude, ∆ i in App. B 1, exhibiting time-dependent oscillations with the same period as S Z .
An applied bias voltage provides a relaxation mechanism for the trapped quasiparticles, making the initial condition relax, as shown in the upper panel of Fig. 5   FIG. 4. (a) Evolution of the impurity spin polarization using NEGFs. We tunnel couple the impurity to a normal lead, with a tunneling amplitude η/∆ = 0.1. For parameters in the 0 phase, U/Γ =0 (black), 1 (red), and 2 (green), the impurity spin converges to an unpolarized stationary state. In contrast, SZ converges to a finite value for parameters in the π phase, U/Γ =4 (blue) and 6 (yellow). (b) Current evolution for the same cases shown in panel (a). The remaining parameters are the same as in Fig. 1: Γ = ∆ = 0.2t0, φ = 0.6π, and ε0 = −U/2. for the smallest U/Γ values. This behavior is better illustrated at longer times. In Fig. 6 we show the time evolution of S Z for an initially trapped spin for the perturbative NEGFs method described in Sec. II C. In this case, the local magnetic moment of the impurity converges to the expected stationary value (S Z = 0 for a system with spin-rotation symmetry) at long times, Fig.  5 (a).
In this case, relaxation is due to transitions between the subgap states and the continuum of states at energies above |∆|. In the limit where the bandwidth W ∆, Γ, the transition rate is given by [80] where J n are the Bessel functions of order n,∆ the induced gap at the impurity, and Σ(ω) = ω √ ∆ 2 − ω 2 /Γ. Eq. (19) predicts a relaxation rate that decreases with the energy gap between the subgap and the continuum of states,∆. It increases with the voltage, showing steps at V = 2|∆|/(2n + 1), Fig. 6 (a). At these voltages, a new MAR channel involving the excitation of a single electron to the continuum of states opens. The applied bias voltage leads to an exponential relaxation of the initially trapped spin in the impurity, 6 (b). We represent the function Ae −γ T t with dashed lines, describing the longtime relaxation rate of the impurity. Here, γ T is given by Eq. (19) and A is taken as the only free parameter in the fit.
For Γ/∆ 1, the subgap states appear close to zero energy exhibiting a small modulation with φ. In this case, the relaxation rate is suppressed for V < ∆. In the opposite regime, Γ/∆ 1, known as the quantum point contact limit, the energy difference between the subgap and the continuum of states tends to zero. In this situation, the relaxation happens in a much faster timescale, determined the ac Josephson period for any voltage [59].
A weak electron-electron interaction does not modify the qualitative picture described above, Fig. 7 (a). We note that Coulomb interactions tend to increase the relaxation rate with respect to the non-interacting situation in the 0 phase, as shown by the green curve. The qualitative picture changes for strong Coulomb interactions, where S Z converges to a non-zero long-time value, blue curve in Fig. 7. This change in the qualitative picture is related to the ground state transition between spin-singlet to spin-doublet [20,27]. It indicates that the state with broken spin-rotation symmetry is a good stationary state of the system. As shown in App. B 2, this state is reached for an initially broken spin-rotation symmetry and V /∆ < 1, illustrating its robustness against perturbations. In contrast, for V /∆ 2 the S Z converges to a zero value, regardless of the initial condition.
In the case where the level is initially fully occupied, a finite bias voltage provides a relaxation mechanism for the impurity charge oscillations, Fig. 7 (b). The charge converges to the expected steady-state in a time scale determined by 1/γ T (19).
Finally, we analyze the non-equilibrium current through the impurity. We note that the symmetrized current, I = (I L − I R )/2, shows a short time dependence on the initial conditions, Fig. 8. This is in contrast to the Γ/∆ 1 limit, where the dependence on the initial condition is lost at very short times [59]. At sufficiently long times (t 1/Γ T ), the current converges to the same value, independently from the impurity initial condition. The current exhibits the characteristic time oscillations of the ac-Josephson effect, with a period given by π/V . At lower voltages, the dependence of the current on the initial conditions is visible at times much larger than any inverse energy in the model, illustrating the exponential dependence of the relaxation rate on the bias voltage, Fig. 6 (a).

IV. CONCLUSIONS
In this work, we have studied the non-equilibrium transient dynamics of an impurity tunnel-coupled to superconducting leads in the presence of local Coulomb interactions. We have developed two theoretical tools to study the time evolution of single-particle observables after a quench of the tunneling amplitudes. We have extended the exact time-dependent density matrix renormalization group (tDMRG) to the case where the leads feature superconducting correlations. We have further-  more developed an approximate method based on the second-order perturbation expansion in the interaction strength to study how the system approaches the steadystate. These two methods agree in the limit of weak to intermediate Coulomb interaction strengths.
We have shown that the system gets trapped in an excited state for generic initial conditions, characterized by a lower supercurrent, similarly to the non-interacting case [58] for the non-interacting case. Local Coulomb interactions are not sufficient to make the system relax to the expected stationary state. We have also considered two additional relaxation mechanisms: (i) coupling the impurity to a normal electrode, where excited quasiparticles can escape, and (ii) a finite bias voltage, which allows quasiparticles to tunnel to the continuum of states where they can diffuse. Both mechanisms predict a spindegenerate steady-state in the 0 phase. In the π phase, the system converges to a steady-state with finite spinpolarization for an initial condition with broken spinrotation symmetry.
The methods developed in this work open the possibility of studying several interesting problems. The numerically unbiased tDMRG method is suitable to analyze the role of extended interactions in the leads and several coupled impurities under non-equilibrium conditions. This last situation is of current interest for the development 1-D topological devices [115]. On the other hand, the perturbative method is suitable for studying the long-time dynamics under generic non-equilibrium conditions. The method can be also extended to including topological superconducting leads [116,117] or additional impurities and leads, which can be either normal or superconducting.  In this section, we show the current-phase relation for some U/Γ values using exact DMRG calculations, Fig.  9. In the non-interacting situation, the current is always positive. It exhibits a non-sinusoidal shape due to the high transmission chosen. In contrast, for 2 < U/Γ < 4 we find a sudden transition of the current from a positive to a negative value, associated to the ground state change from singlet to doublet. Finally, for larger values of the interaction strength, the current is reversed with respect to the U/Γ = 0 case, showing a sinusoidal dependence on φ.

Appendix B: Additional comparisons between tDMRG and NEGFs
In this section, we show some results comparing the two develop methods. We include results for the evolution of the pair amplitude and the evolution of S Z with different bias voltages.

Pair amplitude
In Fig. 10 we show the time evolution of the pair amplitude at the impurity after the connection for the two initial charge configurations: (0,1) and (0,0). In both cases, the exact tDMRG method and the approximated perturbative expansion using NEGFs agree remarkably. Coulomb interaction tends to damp the pair amplitude in the system, becoming nearly zero for U/Γ 0. This is due to the increased energy cost of having two electrons in the system due to the interaction. We note some differences between both methods for U/Γ = 0, generated by the renormalization of the order parameter in tDMRG and absent in the NEGFs perturbative calculation.
The situation becomes richer for an initially occupied or empty impurity level, where the imaginary part of the FIG. 10. Pairing amplitude at the impurity for V = 0 and φ = 0 and an initial occupation (0, 1) (a) and (0, 0) (b). We compare tDMRG (solid lines) with perturbative NEGFs (dashed lines). The remaining parameters are the same as in Fig. 1: Γ = ∆ = 0.2t0, ε0 = −U/2, and U/Γ = 0 (black), 1 (red), 2 (green), and 4 (blue). order parameter oscillates in time around zero value. The real part converges to a stationary value, similarly to Fig.  10 (a). The oscillation frequency is determined by the energy of the subgap state formed at the interface between the impurity and the superconductor [58,59]. These oscillations are undamped for any value of U/Γ. However, their amplitude decrease when increasing Coulomb interaction as the energy difference between the ground state, with an electron in the impurity, and the excited states with a fully occupied or empty level increases.
In the bias voltage situation, the pair amplitude oscillates with a period given by the Josephson frequency, V /π, Fig. 11. This is consistent with the fact that the phase difference in the leads increases linearly in time. We note that the shape of ∆ i is not sinusoidal due to non-adiabatic effects on the time evolution of the system. Similar to the phase-biased situation, ∆ i decreases with the interaction due to the enhanced energy cost of the doubly occupied (or vacant) state in the impurity.

Finite bias voltage
In this section, we show additional results comparing tDMRG with the perturbative NEGF method for a voltage-biased junction. We show that both methods agree for a relatively wide range of parameters. They start to deviate for bias voltages larger than the gap and U/Γ = 4, Fig. 12 (b). However, both methods seem to converge to similar stationary values, approaching S Z = 0 for V > ∆. FIG. 11. Pairing amplitude at the impurity for V = ∆. The remaining parameters are the same as in Fig. 1: Γ = ∆ = 0.2t0, ε0 = −U/2, and U/Γ = 0 (black), 1 (red), 2 (green), and 4 (blue).