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

R. Seoane Souto  Division of Solid State Physics and NanoLund, Lund
University, S-22100 Lund, Sweden  Departamento de Física Teórica de la Materia
Condensada, Condensed Matter Physics Center (IFIMAC) and Instituto Nicolás
Cabrera, Universidad Autónoma de Madrid, E-28049 Madrid, Spain  A. E. Feiguin
Department of Physics, Northeastern University, Boston, Massachusetts 02115,
USA  A. Martín-Rodero  Departamento de Física Teórica de la Materia
Condensada, Condensed Matter Physics Center (IFIMAC) and Instituto Nicolás
Cabrera, Universidad Autónoma de Madrid, E-28049 Madrid, Spain  A. Levy Yeyati
Departamento de Física Teórica de la Materia Condensada, Condensed Matter
Physics Center (IFIMAC) and Instituto Nicolás Cabrera, Universidad Autónoma de
Madrid, E-28049 Madrid, Spain

###### Abstract

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 $\Delta$, 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 $\absolutevalue{\Delta}$. However, a spatial
modulations of $\Delta$ 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 leading to a blockade of 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 lead to an
increase of the low-bias conductance, known as the Kondo effect [18, 19].

The competition between Kondo and pairing correlations lead to a rich phase
diagram [20, 21, 22, 23]. At the transition between the doublet and the
singlet ground state, the subgap states (also known as Yu-Shiba-Rusinov
states) cross the Fermi level [24, 25, 26] leading to a supercurrent sign
reversal, the so-called $0-\pi$ transition [27, 28, 29, 30, 31]. At the ground
state transition, BCS pair correlations are suppressed and the electron
pairing is predominantly spin-triplet close to the impurity [32, 33].

In equilibrium at zero bias voltage, this problem has been studied using
numerically exact methods such as numerical renormalization group [34, 35, 36,
37, 38, 39, 40], quantum Monte Carlo [41, 42, 43, 44], and also real-time
diagrammatic methods [45]. In the regime of weak Coulomb interaction (in
comparison to the gap and the tunnel coupling to the leads), the problem is
well-described by the lowest-order expansion in the Coulomb strength [46, 20,
38, 47]. The resulting mean-field approximation breaks the spin-rotation
symmetry for sufficiently strong Coulomb interaction. Higher-order corrections
in the expansion, describing quantum fluctuations, cures this shortcoming
[48]. A second-order description succeeds in describing the $0-\pi$ transition
and observables in the $0$ phase [22, 49, 50].

Transient dynamics of a superconducting nanojunction has been studied in the
absence of electron correlations [51, 52, 53, 54, 55, 47, 56, 57]. In this
case, the system is not able to thermalize due to the trapping of
quasiparticles in the subgap states [58, 52, 53]. These quasiparticles are
long-lived, being an obstacle towards the efficient control and manipulation
of superconducting devices [59, 60].

The effect of electron correlations has been much less studied in the non-
equilibrium case [61, 11, 62, 63, 64] and in the transient regime. In this
work, we develop an exact 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 time-
dependent Schrödinger equation in finite systems, accounting for many-body
interactions [65, 66, 67, 68, 69]. 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 above $\absolutevalue{\Delta}$
[70]. 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) [71, 72, 73, 74, 75, 76, 77, 78, 79, 53, 80]. 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 $\pi$ 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.1, where we compare both methods and
analyze the convergence to the steady-state. In Sec. III.2 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.3, 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.

##  II Model and formalism

###  II.1 Model

We consider an impurity coupled to superconducting electrodes, described by
the Hamiltonian

$H=H_{i}+H_{leads}+H_{T}(t)\,.$ (1)

The impurity is described by a spin-degenerate level

$H_{i}=\epsilon_{0}\hat{\Psi}^{\dagger}_{0}\hat{\tau}_{z}\hat{\Psi}_{0}+Un_{0\uparrow}n_{0\downarrow}\,,$
(2)

where $\Psi_{0}=(d_{0\uparrow},d^{\dagger}_{0\downarrow})^{T}$ is the Nambu
spinor, with $d_{0\sigma}$ being the annihilation operator for an electron
with spin $\sigma$ at the impurity, and $\hat{\tau}$ are the Pauli matrices in
the Nambu space. Here, $\epsilon_{0}$ is the impurity energy level, $U$ the
charging energy and $n_{0\sigma}=d^{\dagger}_{0\sigma}d_{0\sigma}$, with
$\sigma=\uparrow,\,\downarrow$. The leads are described as a 1-dimensional
chain of electronic sites, given by

$H_{leads}=\sum_{\nu,j}\left(\hat{\Psi}^{\dagger}_{\nu j}\hat{h}_{\nu
j}\hat{\Psi}_{\nu j}+t_{0}\hat{\Psi}^{\dagger}_{\nu
j}\hat{\tau}_{z}\hat{\Psi}_{\nu j-1}\right)\,,$ (3)

where $t_{0}$ is the inter-site chain hopping parameter, $h_{\nu
j}=\xi_{\nu}\hat{\tau}_{z}+\Delta_{\nu}\hat{\tau}_{x}$, with $\xi$ being the
onsite energy of electron at site $j$, and $\Delta_{\nu}$ the superconducting
gap.

The hoping between the impurity and the leads is described by

$H_{T}(t)=\hat{\Psi}^{\dagger}_{\nu
0}\hat{V}_{\nu}(t)\hat{\Psi}_{0}+\mbox{H.c.}\,.$ (4)

Here,
$\hat{V}_{\nu}(t)=\theta(t)V_{\nu}\hat{\tau}_{z}e^{i\hat{\tau}_{z}\phi_{\nu}(t)/2}$,
where we have considered the tunneling amplitude, $|V_{\nu}|$, to be momentum-
independent. The tunneling rate to the leads is given by
$\Gamma_{\nu}=\pi|V_{\nu}|^{2}\rho_{\nu}$, where $\rho_{\nu}=1/\pi t_{0}$ is
the normal density of states at the Fermi level, and
$\Gamma=\Gamma_{L}+\Gamma_{R}$. We consider the impurity to be suddenly
connected to the leads at $t=0$. The superconducting phase is
$\phi_{\nu}(t)=\phi_{\nu}(0)+\mu_{\nu}t$, where $\mu_{\nu}$ is the chemical
potential of $\nu$ electrode.

In this work we focus on the time-evolution of single particle observables.
The spin-polarization of the impurity is given by
$S_{Z}(t)=[n_{\uparrow}(t)-n_{\downarrow}(t)]/2$, where

$n_{\sigma}(t)=\left\langle
d_{0\sigma}^{\dagger}(t)d_{0\sigma}(t)\right\rangle\,,$ (5)

is the occupation per spin and $n_{0}(t)=n_{\uparrow}(t)+n_{\downarrow}(t)$
the impurity charge. The current at $\nu=L,R$ interface is given by

$I_{\nu}(t)=i\,\mbox{Tr}\left[\,\hat{V}_{\nu}(t)\left\langle\hat{\Psi}^{\dagger}_{\nu
0}(t)\hat{\Psi}_{0}(t)\right\rangle\right]+\mbox{H.c.},$ (6)

where the trace is taken in the Nambu space. We define the symmetrized current
as $\left\langle I\right\rangle=\left(I_{L}-I_{R}\right)/2$. Finally, the
electron pair amplitude at the impurity is given by

$\Delta_{i}(t)=\left\langle d_{0\uparrow}(t)d_{0\downarrow}(t)\right\rangle.$
(7)

###  II.2 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
[65, 66, 67, 68, 69]. Among other applications, it has been used to study
time-dependent transport in a number of one-dimensional setups [81, 82, 83,
84, 85, 86, 87, 88, 89, 90]. 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 wave-front
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. [91, 92].

In this work, we extend the tDMRG method to the case where the leads are
superconducting. A particular consideration to take into account in the
presence of superconductivity is that the particle number is conserved modulo
two [93] (i.e., number parity and spin are used as 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
pre-determined tolerance (typically $10^{-6}$), and we stop the simulations
either when the accuracy grows beyond this value, or when the effects of the
boundaries become dominant. The observables we represent are given by Eqs.
(5-7) for the impurity population, the current, and the pair amplitude.

###  II.3 Non-equilibrium Green functions

Perturbation theory on the interaction strength, $U$, provides accurate
results in the $U/\Gamma\lesssim 1$ regime. This approximation has been shown
to describe accurately physical observables in the stationary regime and the
non-magnetic phase [22, 49, 50]. 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. [94, 95], 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) [96, 97, 98, 99]. The impurity Green function is given by

$\hat{G}_{0}=\left(\hat{g}^{-1}_{0}-\hat{\Sigma}_{leads}-\hat{\Sigma}_{int}\right)^{-1}\,.$
(8)

We solve the problem in a discrete time mesh, where the inverse of the
isolated impurity Green function is given by

$\hat{g}^{-1}_{0}=\begin{pmatrix}g^{-1}_{e}&0\\\ 0&g^{-1}_{h}\end{pmatrix}\,.$
(9)

The electron and hole parts of the inverse Green function in the discrete mesh
are described in Refs. [52, 53].

The leads self-energy is given by [15]

$\hat{\Sigma}^{\alpha\beta}_{leads}(t,t^{\prime})=\alpha\beta\sum_{\nu}\hat{V}_{\nu}(t)\hat{g}^{\alpha\beta}_{\nu}\hat{V}^{\ast}_{\nu}(t^{\prime})\,,$
(10)

where $\alpha$, $\beta=\pm$ denote the Keldysh branch. The leads boundary
Green functions are given by

$\hat{g}^{\alpha\beta}_{\nu}(\omega)=\begin{pmatrix}\mathfrak{g}^{\alpha\beta}_{\nu}(\omega)&\mathfrak{f}^{\alpha\beta}_{\nu}(\omega)\\\
\mathfrak{f}^{\alpha\beta}_{\nu}(\omega)&\mathfrak{g}^{\alpha\beta}_{\nu}(\omega)\end{pmatrix}\,.$
(11)

The bare Green functions for a superconducting chain are given by

$\mathfrak{g}^{r/a}_{\nu}(\omega)=\frac{\omega\pm
i0^{+}}{2t^{2}_{0}}\left[1-\frac{\sqrt{W^{2}-(\omega\pm
i0^{+})^{2}}}{\sqrt{\Delta^{2}-(\omega\pm i0^{+})^{2}}}\right]$ (12)

and

$\mathfrak{f}^{r/a}_{\nu}(\omega)=-\Delta/(\omega\pm
i0^{+})\mathfrak{g}^{r/a}_{\nu}(\omega)$ (13)

with $W=4t_{0}^{2}+\Delta^{2}$ being the bandwidth and $0^{+}$ an
infinitesimal. The Keldysh components are given by
$\hat{g}^{+-}_{\nu}(\omega)=[\hat{g}^{a}_{\nu}(\omega)-\hat{g}^{r}_{\nu}(\omega)]n_{F}(\omega)$
and
$\hat{g}^{-+}_{\nu}(\omega)=[\hat{g}^{a}_{\nu}(\omega)-\hat{g}^{r}_{\nu}(\omega)][n_{F}(\omega)-1]$,
where $n_{F}$ is the Fermi distribution function.

To calculate the interaction self-energy we perform a perturbation in
$U/\Gamma$. Up to second order in the perturbative parameter, the self-energy
is given by $\hat{\Sigma}_{\rm int}=\hat{\Sigma}_{\rm\rm
HF}+\hat{\Sigma}^{(2)}$. The linear term in $U$ is the mean-field
approximation, given by

$\begin{split}\left[\Sigma_{HF}(t,t^{\prime})\right]^{\alpha\beta}_{nn}=U\alpha\,n_{\bar{n}\bar{n}}(t)\,\delta(t,t^{\prime})\delta_{\alpha\beta}\,,\\\
\left[\Sigma_{HF}(t,t^{\prime})\right]^{\alpha\beta}_{n\bar{n}}=U\alpha\,n_{n\bar{n}}(t)\,\delta(t,t^{\prime})\delta_{\alpha\beta}\,,\end{split}$
(14)

where $n=1,2$ denote the Nambu component and $\bar{n}$ the opposite component
to $n$. Here, $n_{nn}$ denotes the charge per spin, given by
$n_{\uparrow}=n_{11}(t)=-iG^{+-}_{11}(t,t)$ and
$n_{\downarrow}=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_{n\bar{n}}(t)=-iG^{+-}_{n\bar{n}}(t,t)$. The pairing
amplitude in the impurity, $\Delta_{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/\Gamma$, the self-
energy is given by

$\begin{split}\left[\Sigma^{(2)}(t,t^{\prime})\right]^{\alpha\beta}_{nn}&=-U^{2}\sum_{j=1,2}(-1)^{j-n}\left[G_{HF}(t,t^{\prime})\right]_{nj}^{\alpha\beta}\left[G_{HF}(t,t^{\prime})\right]_{\bar{n}j}^{\alpha\beta}\left[G_{HF}(t^{\prime},t)\right]_{\bar{n}\bar{n}}^{\beta\alpha}\,,\\\
\left[\Sigma^{(2)}(t,t^{\prime})\right]^{\alpha\beta}_{n\bar{n}}&=U^{2}\sum_{j=1,2}(-1)^{j-n}\left[G_{HF}(t,t^{\prime})\right]_{nj}^{\alpha\beta}\left[G_{HF}(t,t^{\prime})\right]_{\bar{n}j}^{\alpha\beta}\left[G_{HF}(t^{\prime},t)\right]_{n\bar{n}}^{\beta\alpha}\,,\end{split}$
(15)

where we have perturbed over the mean-field Green function, given by

$\hat{G}_{HF}=\left(\hat{g}^{-1}_{0}-\hat{\Sigma}_{leads}-\hat{\Sigma}_{HF}\right)^{-1}\,.$
(16)

This approximation has been shown to provide accurate results in the normal
case, describing the onset of the Kondo peak at temperatures lower than the
Kondo temperature [100, 94].

The system observables can be obtained from the Green’s function (8). The
average charge and spin at the impurity are given by $\left\langle
n_{d}\right\rangle=n_{\uparrow}+n_{\downarrow}$ and $\left\langle
S_{Z}\right\rangle=(n_{\uparrow}-n_{\downarrow}$)/2. The current at a given
interface can be calculated as [22]

$I_{\nu}=\mbox{Tr}\left\\{\tau_{z}\left[\hat{V}_{\nu}(t)\hat{G}^{+-}_{\nu
0}(t,t)-\hat{V}^{\ast}_{\nu}(t)\hat{G}^{+-}_{0\nu}(t,t)\right]\right\\},\ $
(17)

where the trace is taken over the Nambu space.

##  III Results

###  III.1 Phase-biased junction

Figure 1: Comparison between tDMRG (solid line) and perturbative NEGFs (dashed
line) methods. We show the time evolution of the impurity spin, $S_{Z}$ (a)
and the symmetrized current (b) for $U/\Gamma=$0, 1, 2 and 4 from bottom to
top. Current curves have been shifted upwards by $U/40$ for clarity. The
remaining parameters are $\Gamma=\Delta=0.2t_{0}$, $\phi=0.6\pi$,
$\epsilon_{0}=-U/2$ and the initial impurity occupation
$(n_{\uparrow},n_{\downarrow})=(1,0)$.

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
[101]. 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_{\uparrow}(0),n_{\downarrow}(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/\Gamma$ values. In Fig. 1 we show results for $S_{Z}$ and the current
through the system. In App. B.1 we display the time evolution of the pairing
potential at the impurity. The results of both methods start to deviate for
$U/\Gamma\approx 6$, predicting different stationary values for $S_{Z}$ (not
shown).

The two methods describe the trapping of a magnetic moment in the impurity at
long times [52, 53]. We have checked that no spin relaxation happens for times
$t\lesssim 50/\Delta$ (the steady-state value has been reached at $t\sim
10/\Delta$). This is an indication that localized Coulomb interactions do not
provide an effective relaxation mechanism for the initially trapped spin in
the system. On the contrary, localized Coulomb interactions tend to quench the
short-time spin relaxation, leading to a higher stationary value of $S_{Z}$
when increasing $U$, Fig. 1 (a). The absence of thermalization causes a
supercurrent reduction with respect to the expected value at equilibrium [52].

Figure 2: Evolution of the impurity charge for an initially occupied level,
$(n_{\uparrow}(0),n_{\downarrow}(0))=(1,1)$. The remaining parameters are the
same as in Fig. 1

In the case of a fully occupied or empty impurity level, the charge exhibits
undamped oscillations [53, 54], see 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. [53] 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/\Delta\gg 1$,
the oscillations become negligible and $n_{0}$ converges to the expected
stationary state.

###  III.2 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, (a) panel, compared to the equilibrium result, (b) panel.
The difference is better observed before the supercurrent sign reversal at
$U_{c}/\Gamma\sim 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 [52].

Figure 3: (a) Time-dependent current calculated using tDMRG for the same
parameters shown in Fig. 1. (b) Stationary current, where the colored dots
indicate the value of $U$ in the left panel.

The absence of thermalization is consistent with the observed long relaxation
times for trapped quasiparticles observed in experiments [59, 60]. The finite
relaxation time in these systems could be due photon or phonon
emission/absorption processes [102] 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.

To simulate the mechanisms commented above for quasiparticle relaxation, we
attach an additional normal lead to the impurity [52, 54, 55], Fig. 4. In this
case, the supercurrent exhibits a sign change when increasing $U/\Gamma$,
characteristic of the $0-\pi$ 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, $\eta$. The relaxation
rate is given by $\eta$ when the system is far away from the $0-\pi$
transition. In contrast, the system relaxation close to the transition
($U_{c}/\Gamma\approx 3.7$ for this choice of parameters) becomes slower, with
a relaxation time exceeding $1/\eta$. 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-\pi$ transition becomes smooth due to the finite $\eta$ 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/\Gamma$ values shown, $S_{Z}$ converges
to zero, consistent with the convergence of the system to the expected singlet
ground state. The situation is different in the $\pi$ 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
$\eta$ 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.

Figure 4: (a) Evolution of the impurity spin polarization using NEGFs. We
tunnel couple the impurity to a normal lead, with a tunneling amplitude
$\eta/\Delta=0.1$. We show results for $U/\Gamma=0$, 1, 2, 4 and 6 from bottom
to top (b) Current evolution for the same cases shown in panel (a). The
remaining parameters are the same as in Fig. 1.

###  III.3 Voltage-biased junction

A bias voltage applied to the junction produces a linear increase of
superconducting phase difference in time as $\phi(t)=\phi(0)+V\,t$. For
$V<2\Delta$, multiple Andreev reflections lead to the transference of multiple
charges between the superconductors, resulting in a time-oscillating as well
as a dc current in the steady-state. In the low-bias situation
($V\ll\Delta,\Gamma$), 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 $\Delta$ and $\Gamma$. In this
situation, non-adiabatic effects are important and higher harmonics have to be
considered [74].

In Fig. 5 we compare the results from tDMRG and the perturbative NEGFs for
$V=\Delta$. Similarly to the phase-biased situation studied before, both
methods agree for $U/\Gamma\lesssim 1$. They start to deviate for the largest
$U/\Gamma$ 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, $\Delta_{i}$ in B.1,
exhibiting time-dependent oscillations with the same period as $S_{Z}$.

Figure 5: Spin polarization (a) and current (b) for a voltage-biased junction,
$V=\Delta$. The remaining parameters are the same as in Fig. 1. Curves in
panel (b) have been shifted down by $U/20$ for clarity.

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 for the smallest $U/\Gamma$ 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.3. 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 $\absolutevalue{\Delta}$. In the
limit where the bandwidth $W\gg\Delta,\Gamma$, the transition rate is given by
[70]

$\Gamma_{T}=2\sum_{n\geq
0}J_{n}^{2}(\tilde{\Delta}/V)\mbox{Im}\left[\Sigma(nV+V/2)\right]\,,$ (18)

where $J_{n}$ are the Bessel functions of order $n$, $\tilde{\Delta}$ the
induced gap at the impurity, and
$\Sigma(\omega)=\omega\sqrt{\Delta^{2}-\omega^{2}}/\Gamma$. Eq. (18) predicts
a relaxation rate that decreases with the energy gap between the subgap and
the continuum of states, $\tilde{\Delta}$. It increases with the voltage,
showing steps at $V=2\absolutevalue{\Delta}/(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^{-\gamma_{T}t}$ with dashed lines, describing
the long-time relaxation rate of the impurity. Here, $\gamma_{T}$ is given by
Eq. (18) and $A$ is taken as the only free parameter in the fit.

For $\Gamma/\Delta\ll 1$, the subgap states appear close to zero energy
exhibiting a small modulation with $\phi$. In this case, the relaxation rate
is suppressed for $V<\Delta$. In the opposite regime, $\Gamma/\Delta\gg 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 [53].

Figure 6: (a) Analytic relaxation rate of a voltage-biased junction, Eq. (18).
(b) Relaxation of the initially trapped spin, comparing numerical results
(solid lines) and the analytic rates (dashed lines). We show results for
$U/\Gamma=0$ and $V/\Delta=0.2$, $0.3$, $0.5$ and $0.8$ from top to bottom.
The remaining parameters are the same as in Fig. 1.

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
[27, 20]. 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/\Delta<1$,
illustrating its robustness against perturbations. In contrast, for
$V/\Delta\gtrsim 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/\gamma_{T}$ (18).

Figure 7: (a) Mean value of $S_{Z}$ for the $(0,1)$ initial condition. (b)
Impurity charge for the $(1,1)$ initial condition. We show results for the
NEGFs with $V=\Delta/2$ and $U/\Gamma=0$ (black), $1$ (red), $2$ (green) and
$4$ (blue). The remaining parameters are the same as in Fig. 5.

Finally, we analyze the non-equilibrium current through the impurity. We note
that the symmetrized current, $\left\langle I\right\rangle=(I_{L}-I_{R})/2$,
shows a short time dependence on the initial conditions, Fig. 8. This is in
contrast to the $\Gamma/\Delta\gg 1$ limit, where the dependence on the
initial condition is lost at very short times [53]. At sufficiently long times
($t\gg 1/\Gamma_{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 $\pi/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).

Figure 8: Symmetrized current for $V=0.2\Delta$ (black lines) and
$V=0.5\Delta$ (red lines), for $U/\Delta=1$. We show results for $(0,1)$
(solid line) and $(1,1)$ (dashed lines) initial occupations. The remaining
parameters are the same as in Fig. 1.

##  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 furthermore developed an approximate
method based on the second-order perturbation expansion in the interaction
strength to study how the system approaches the steady-state. 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 as shown in Ref.
[52] for the non-interacting case. In this work, we have found that localized
Coulomb interactions are not sufficient to make the system relax to the
expected equilibrium 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 spin-degenerate steady-state in the $0$ phase. In
the $\pi$ phase, the system converges to a steady-state with finite spin-
polarization for an initial condition with broken spin-rotation 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 [103]. 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 [104] or additional impurities and
leads, which can be either normal or superconducting.

##  V Acknowledgments

RSS acknowledges funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme under Grant
Agreement No. 856526 and from QuantERA project “2D hybrid materials as a
platform for topological quantum computing”. AEF is supported by the US
Department of Energy (DOE), Office of Science, Basic Energy Sciences grant
number DE-SC0019275. ALY and AMR acknowledge support from the Spanish MICINN
through grants FIS2017-84860-R and through the “María de Maeztu” Programme for
Units of Excellence in R&D (Grant No. MDM-2014-0377) and by EU through grant
no. 828948 (AndQC).

## References

  * Bardeen  _et al._ [1957] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 106, 162 (1957).
  * de Gennes and Saint-James [1963] P. de Gennes and D. Saint-James, Physics Letters 4, 151 (1963).
  * Andreev [1964] A. F. Andreev, Sov. Phys. JETP 19, 1228 (1964).
  * Andreev [1966] A. F. Andreev, Sov. Phys. JETP 22, 18 (1966).
  * Kulik [1970] I. O. Kulik, Sov. Phys. JETP 30, 944 (1970).
  * Sauls [2018] J. A. Sauls, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 376, 20180140 (2018).
  * Prada  _et al._ [2020] E. Prada, P. San-Jose, M. W. A. de Moor, A. Geresdi, E. J. H. Lee, J. Klinovaja, D. Loss, J. Nygård, R. Aguado, and L. P. Kouwenhoven, Nature Reviews Physics 2, 575 (2020).
  * Cuevas  _et al._ [2001] J. C. Cuevas, A. Levy Yeyati, and A. Martín-Rodero, Phys. Rev. B 63, 094515 (2001).
  * Balatsky  _et al._ [2006] A. V. Balatsky, I. Vekhter, and J.-X. Zhu, Rev. Mod. Phys. 78, 373 (2006).
  * Jarillo-Herrero  _et al._ [2006] P. Jarillo-Herrero, J. A. van Dam, and L. P. Kouwenhoven, Nature 439, 953 (2006).
  * Eichler  _et al._ [2007] A. Eichler, M. Weiss, S. Oberholzer, C. Schönenberger, A. Levy Yeyati, J. C. Cuevas, and A. Martín-Rodero, Phys. Rev. Lett. 99, 126602 (2007).
  * Grove-Rasmussen  _et al._ [2007] K. Grove-Rasmussen, H. I. Jørgensen, and P. E. Lindelof, New J. Phys. 9, 124 (2007).
  * Eichler  _et al._ [2009] A. Eichler, R. Deblock, M. Weiss, C. Karrasch, V. Meden, C. Schönenberger, and H. Bouchiat, Phys. Rev. B 79, 161407 (2009).
  * De Franceschi  _et al._ [2010] S. De Franceschi, L. Kouwenhoven, C. Schönenberger, and W. Wernsdorfer, Nature Nanotechnology 5, 703 (2010).
  * Martín-Rodero and Yeyati [2011] A. Martín-Rodero and A. L. Yeyati, Advances in Physics 60, 899 (2011).
  * Meden [2019] V. Meden, Journal of Physics: Condensed Matter 31, 163001 (2019).
  * Beenakker and van Houten [1991] C. Beenakker and H. van Houten, in  _Semiconductor Heterostructures and Nanostructures_, Solid State Physics, Vol. 44, edited by H. Ehrenreich and D. Turnbull (Academic Press, 1991) pp. 1–228.
  * Goldhaber-Gordon  _et al._ [1998] D. Goldhaber-Gordon, H. Shtrikman, D. Mahalu, D. Abusch-Magder, U. Meirav, and M. A. Kastner, Nature 391, 156 (1998).
  * Kouwenhoven and Glazman [2001] L. Kouwenhoven and L. Glazman, Physics World 14, 33 (2001).
  * Rozhkov and Arovas [1999] A. V. Rozhkov and D. P. Arovas, Phys. Rev. Lett. 82, 2788 (1999).
  * Rozhkov and Arovas [2000] A. V. Rozhkov and D. P. Arovas, Phys. Rev. B 62, 6687 (2000).
  * Vecino  _et al._ [2003] E. Vecino, A. Martín-Rodero, and A. L. Yeyati, Phys. Rev. B 68, 035105 (2003).
  * Tanaka  _et al._ [2007] Y. Tanaka, A. Oguri, and A. C. Hewson, New J. Phys 9, 115 (2007).
  * Pillet  _et al._ [2010] J.-D. Pillet, C. H. L. Quay, P. Morfin, C. Bena, A. L. Yeyati, and P. Joyez, Nature Physics 6, 965 (2010).
  * Dirks  _et al._ [2011] T. Dirks, T. L. Hughes, S. Lal, B. Uchoa, Y.-F. Chen, C. Chialvo, P. M. Goldbart, and N. Mason, Nature Physics 7, 386 (2011).
  * Jellinggaard  _et al._ [2016] A. Jellinggaard, K. Grove-Rasmussen, M. H. Madsen, and J. Nygård, Phys. Rev. B 94, 064520 (2016).
  * Glazman and Matveev [1989] L. I. Glazman and K. A. Matveev, JETP Lett. 49, 659 (1989).
  * van Dam  _et al._ [2006] J. A. van Dam, Y. V. Nazarov, E. P. A. M. Bakkers, S. De Franceschi, and L. P. Kouwenhoven, Nature 442, 667 (2006).
  * Cleuziou  _et al._ [2006] J.-P. Cleuziou, W. Wernsdorfer, V. Bouchiat, T. Ondarçuhu, and M. Monthioux, Nature Nanotechnology 1, 53 (2006).
  * Jørgensen  _et al._ [2007] H. I. Jørgensen, T. Novotný, K. Grove-Rasmussen, K. Flensberg, and P. E. Lindelof, Nano Letters 7, 2441 (2007).
  * Delagrange  _et al._ [2018] R. Delagrange, R. Weil, A. Kasumov, M. Ferrier, H. Bouchiat, and R. Deblock, Physica B: Condensed Matter 536, 211 (2018).
  * Kuzmanovski  _et al._ [2020] D. Kuzmanovski, R. S. Souto, and A. V. Balatsky, Phys. Rev. B 101, 094505 (2020).
  * Perrin  _et al._ [2020] V. Perrin, F. L. N. Santos, G. C. Ménard, C. Brun, T. Cren, M. Civelli, and P. Simon, Phys. Rev. Lett. 125, 117003 (2020).
  * Yoshioka and Ohashi [2000] T. Yoshioka and Y. Ohashi, Journal of the Physical Society of Japan 69, 1812 (2000).
  * Choi  _et al._ [2004] M.-S. Choi, M. Lee, K. Kang, and W. Belzig, Phys. Rev. B 70, 020502 (2004).
  * Oguri  _et al._ [2004] A. Oguri, Y. Tanaka, and A. C. Hewson, Journal of the Physical Society of Japan 73, 2494 (2004).
  * Karrasch  _et al._ [2008] C. Karrasch, A. Oguri, and V. Meden, Phys. Rev. B 77, 024517 (2008).
  * Martín-Rodero and Yeyati [2012] A. Martín-Rodero and A. L. Yeyati, Journal of Physics: Condensed Matter 24, 385303 (2012).
  * Domański  _et al._ [2016] T. Domański, I. Weymann, M. Barańska, and G. Górski, Scientific Reports 6, 23336 (2016).
  * Grove-Rasmussen  _et al._ [2018] K. Grove-Rasmussen, G. Steffensen, A. Jellinggaard, M. H. Madsen, R. Žitko, J. Paaske, and J. Nygård, Nature Communications 9, 2376 (2018).
  * Siano and Egger [2004] F. Siano and R. Egger, Phys. Rev. Lett. 93, 047002 (2004).
  * Luitz and Assaad [2010] D. J. Luitz and F. F. Assaad, Phys. Rev. B 81, 024509 (2010).
  * Luitz  _et al._ [2012] D. J. Luitz, F. F. Assaad, T. Novotný, C. Karrasch, and V. Meden, Phys. Rev. Lett. 108, 227001 (2012).
  * Pokorný and Novotný [2021] V. Pokorný and T. c. v. Novotný, Phys. Rev. Research 3, 023013 (2021).
  * Governale  _et al._ [2008] M. Governale, M. G. Pala, and J. König, Phys. Rev. B 77, 134513 (2008).
  * Shiba [1973] H. Shiba, Progress of Theoretical Physics 50, 50 (1973).
  * Wrześniewski  _et al._ [2020] K. Wrześniewski, B. Baran, R. Taranko, T. Domański, and I. Weymann, arXiv:2007.10747 (2020).
  * Clerk and Ambegaokar [2000] A. A. Clerk and V. Ambegaokar, Phys. Rev. B 61, 9109 (2000).
  * Žonda  _et al._ [2015] M. Žonda, V. Pokorný, V. Janiš, and T. Novotný, Scientific Reports 5, 8821 (2015).
  * Žonda  _et al._ [2016] M. Žonda, V. Pokorný, V. Janiš, and T. Novotný, Phys. Rev. B 93, 024523 (2016).
  * Perfetto  _et al._ [2009] E. Perfetto, G. Stefanucci, and M. Cini, Phys. Rev. B 80, 205408 (2009).
  * Souto  _et al._ [2016] R. S. Souto, A. Martín-Rodero, and A. L. Yeyati, Phys. Rev. Lett. 117, 267701 (2016).
  * Souto  _et al._ [2017] R. S. Souto, A. Martín-Rodero, and A. L. Yeyati, Phys. Rev. B 96, 165444 (2017).
  * Taranko and Domański [2018] R. Taranko and T. Domański, Phys. Rev. B 98, 075420 (2018).
  * Taranko  _et al._ [2019] R. Taranko, T. Kwapiński, and T. Domański, Phys. Rev. B 99, 165419 (2019).
  * Averin  _et al._ [2020] D. V. Averin, G. Wang, and A. S. Vasenko, Phys. Rev. B 102, 144516 (2020).
  * Taranko  _et al._ [2021] R. Taranko, K. Wrześniewski, B. Baran, I. Weymann, and T. Domański, Phys. Rev. B 103, 165430 (2021).
  * Stefanucci  _et al._ [2010] G. Stefanucci, E. Perfetto, and M. Cini, Phys. Rev. B 81, 115446 (2010).
  * Zgirski  _et al._ [2011] M. Zgirski, L. Bretheau, Q. Le Masne, H. Pothier, D. Esteve, and C. Urbina, Phys. Rev. Lett. 106, 257003 (2011).
  * Janvier  _et al._ [2015] C. Janvier, L. Tosi, L. Bretheau, Ç. Ö. Girit, M. Stern, P. Bertet, P. Joyez, D. Vion, D. Esteve, M. F. Goffman, H. Pothier, and C. Urbina, Science 349, 1199 (2015).
  * Avishai  _et al._ [2003] Y. Avishai, A. Golub, and A. D. Zaikin, Phys. Rev. B 67, 041301 (2003).
  * Dell’Anna  _et al._ [2008] L. Dell’Anna, A. Zazunov, and R. Egger, Phys. Rev. B 77, 104525 (2008).
  * Andersen  _et al._ [2011] B. M. Andersen, K. Flensberg, V. Koerting, and J. Paaske, Phys. Rev. Lett. 107, 256802 (2011).
  * Avriller and Pistolesi [2015] R. Avriller and F. Pistolesi, Phys. Rev. Lett. 114, 037003 (2015).
  * White and Feiguin [2004] S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004).
  * Feiguin and White [2005] A. E. Feiguin and S. R. White, Phys. Rev. B 72, 220401 (2005).
  * Feiguin [2011] A. E. Feiguin, in  _XV Training Course in the Physics of Strongly Correlated Systems_ , Vol. 1419 (AIP Proceedings, 2011) p. 5.
  * Feiguin [2013] A. E. Feiguin, in  _Strongly correlated systems: Numerical methods_ , edited by A. Avella and F. Mancini (Springer, Heidelberg, Berlin, 2013).
  * Paeckel  _et al._ [2019] S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig, Annals of Physics 411, 167998 (2019).
  * Yeyati  _et al._ [2003] A. L. Yeyati, A. Martín-Rodero, and E. Vecino, Phys. Rev. Lett. 91, 266802 (2003).
  * Arnold [1987] G. B. Arnold, Journal of Low Temperature Physics 68, 1 (1987).
  * Averin and Bardas [1995] D. Averin and A. Bardas, Phys. Rev. Lett. 75, 1831 (1995).
  * Bratus’  _et al._ [1995] E. N. Bratus’, V. S. Shumeiko, and G. Wendin, Phys. Rev. Lett. 74, 2110 (1995).
  * Cuevas  _et al._ [1996] J. C. Cuevas, A. Martín-Rodero, and A. L. Yeyati, Phys. Rev. B 54, 7366 (1996).
  * Scheer  _et al._ [1997] E. Scheer, P. Joyez, D. Esteve, C. Urbina, and M. H. Devoret, Phys. Rev. Lett. 78, 3535 (1997).
  * Scheer  _et al._ [1998] E. Scheer, N. Agraït, J. C. Cuevas, A. L. Yeyati, B. Ludoph, A. Martín-Rodero, G. R. Bollinger, J. M. van Ruitenbeek, and C. Urbina, Nature 394, 154 (1998).
  * Ludoph  _et al._ [2000] B. Ludoph, N. van der Post, E. N. Bratus’, E. V. Bezuglyi, V. S. Shumeiko, G. Wendin, and J. M. van Ruitenbeek, Phys. Rev. B 61, 8561 (2000).
  * Goffman  _et al._ [2000] M. F. Goffman, R. Cron, A. Levy Yeyati, P. Joyez, M. H. Devoret, D. Esteve, and C. Urbina, Phys. Rev. Lett. 85, 170 (2000).
  * Yeyati  _et al._ [1997] A. L. Yeyati, J. C. Cuevas, A. López-Dávalos, and A. Martín-Rodero, Phys. Rev. B 55, R6137 (1997).
  * Villas  _et al._ [2020] A. Villas, R. L. Klees, H. Huang, C. R. Ast, G. Rastelli, W. Belzig, and J. C. Cuevas, Phys. Rev. B 101, 235445 (2020).
  * Schmitteckert [2004] P. Schmitteckert, Phys. Rev. B 70, 121302 (2004).
  * Al-Hassanieh  _et al._ [2006] K. A. Al-Hassanieh, A. E. Feiguin, J. A. Riera, C. A. Büsser, and E. Dagotto, Phys. Rev. B 73, 195304 (2006).
  * Al-Hassanieh  _et al._ [2008] K. A. Al-Hassanieh, F. A. Reboredo, A. E. Feiguin, I. González, and E. Dagotto, Phys. Rev. Lett. 100, 166403 (2008).
  * Boulat  _et al._ [2008] E. Boulat, H. Saleur, and P. Schmitteckert, Phys. Rev. Lett. 101, 140601 (2008).
  * Dias da Silva  _et al._ [2010] L. G. G. V. Dias da Silva, K. A. Al-Hassanieh, A. E. Feiguin, F. A. Reboredo, and E. Dagotto, Phys. Rev. B 81, 125113 (2010).
  * Heidrich-Meisner  _et al._ [2010] F. Heidrich-Meisner, I. González, K. A. Al-Hassanieh, A. E. Feiguin, M. J. Rozenberg, and E. Dagotto, Phys. Rev. B 82, 205110 (2010).
  * Rincón  _et al._ [2014] J. Rincón, K. A. Al-Hassanieh, A. E. Feiguin, and E. Dagotto, Phys. Rev. B 90, 155112 (2014).
  * Petrović  _et al._ [2021a] M. D. Petrović, P. Mondal, A. E. Feiguin, and B. K. Nikolić, Phys. Rev. Lett. 126, 197202 (2021a).
  * Petrović  _et al._ [2021b] M. D. Petrović, P. Mondal, A. E. Feiguin, P. Plecháč, and B. K. Nikolić, Spintronics meets density matrix renormalization group: Quantum spin-torque-driven nonclassical magnetization reversal and dynamical buildup of long-range entanglement (2021b).
  * Bragança  _et al._ [2021] H. Bragança, M. F. Cavalcante, R. G. Pereira, and M. C. O. Aguiar, Phys. Rev. B 103, 125152 (2021).
  * Dias da Silva  _et al._ [2008] L. G. G. V. Dias da Silva, F. Heidrich-Meisner, A. E. Feiguin, C. A. Büsser, G. B. Martins, E. V. Anda, and E. Dagotto, Phys. Rev. B 78, 195317 (2008).
  * Feiguin  _et al._ [2008] A. Feiguin, P. Fendley, M. P. A. Fisher, and C. Nayak, Phys. Rev. Lett. 101, 236801 (2008).
  * Feiguin  _et al._ [2007] A. E. Feiguin, S. R. White, and D. J. Scalapino, Phys. Rev. B 75, 024505 (2007).
  * Seoane Souto  _et al._ [2018] R. Seoane Souto, R. Avriller, A. Levy Yeyati, and A. Martín-Rodero, New Journal of Physics 20, 083039 (2018).
  * Seoane Souto [2020] R. Seoane Souto,  _Quench Dynamics in Interacting and Superconducting Nanojunctions_ (Springer Nature, 2020).
  * Tang  _et al._ [2014] G.-M. Tang, F. Xu, and J. Wang, Phys. Rev. B 89, 205310 (2014).
  * Tang and Wang [2014] G.-M. Tang and J. Wang, Phys. Rev. B 90, 195422 (2014).
  * Seoane Souto  _et al._ [2015] R. Seoane Souto, R. Avriller, R. C. Monreal, A. Martín-Rodero, and A. Levy Yeyati, Phys. Rev. B 92, 125435 (2015).
  * Avriller  _et al._ [2019] R. Avriller, R. S. Souto, A. Martín-Rodero, and A. L. Yeyati, Phys. Rev. B 99, 121403 (2019).
  * White [1992] J. A. White, Phys. Rev. B 45, 1100 (1992).
  * Annett [2004] J. F. Annett,  _Superconductivity, Superfluids and Condensates_, Oxford Master Series in Physics (OUP Oxford, 2004).
  * Olivares  _et al._ [2014] D. G. Olivares, A. L. Yeyati, L. Bretheau, i. m. c. O. Girit, H. Pothier, and C. Urbina, Phys. Rev. B 89, 104504 (2014).
  * Feldman  _et al._ [2017] B. E. Feldman, M. T. Randeria, J. Li, S. Jeon, Y. Xie, Z. Wang, I. K. Drozdov, B. Andrei Bernevig, and A. Yazdani, Nature Physics 13, 286 (2017).
  * Seoane Souto  _et al._ [2020] R. Seoane Souto, K. Flensberg, and M. Leijnse, Phys. Rev. B 101, 081407 (2020).

##  Appendix A Current-phase relation

Figure 9: Current-phase relation of a phase-biased junction for the same
parameters as in Fig. 5 using tDMRG.

In this section, we show the current-phase relation for some $U/\Gamma$ 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\Gamma<4$ we note a sudden
transition of the current from a positive to a negative value. Finally, for
larger values of the interaction strength, the current is reversed with
respect to the $U/\Gamma=0$ case, showing a sinusoidal dependence on $\phi$.

##  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.

###  B.1 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/\Gamma\gg 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/\Gamma=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 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 [52, 53]. These oscillations are undamped for any value of
$U/\Gamma$. 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.

Figure 10: Pairing amplitude at the impurity for $V=0$ and $\phi=0$ and an
initial occupation $(0,1)$ (a) and $(0,0)$ (b).

In the bias voltage situation, the pair amplitude oscillates with a period
given by the Josephson frequency, $V/\pi$, 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 $\Delta_{i}$ is not sinusoidal due to non-adiabatic
effects on the time evolution of the system. Similar to the phase-biased
situation, $\Delta_{i}$ decreases with the interaction due to the enhanced
energy cost of the doubly occupied (or vacant) state in the impurity.

Figure 11: Pairing amplitude at the impurity for $V=\Delta$.

###  B.2 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/\Gamma=4$, Fig. 12 (b). However,
both methods seem to converge to similar stationary values, approaching
$S_{Z}=0$ for $V>\Delta$.

Figure 12: Magnetization evolution for $U/\Delta=2$ (a) and $4$ (b), and
several voltage values. We compare tDMRG (solid lines) with the second-order
approximation (dashed). The remaining parameters are the same as in Fig. 7.