Dynamical Onset of Light-Induced Unconventional Superconductivity – a Yukawa-Sachdev-Ye-Kitaev study

,

Strong electronic correlations, especially in cuprates and iron based superconductors, make systematic theoretical investigations challenging.In these systems, superconductivity is believed to emerge out of a strongly correlated non-Fermi liquid (nFl) normal state [46][47][48][49] and studies of the real time emergence of the superconducting condensate typically have to rely on perturbative approximations or purely phenomenological models.It is hence highly desirable to investigate the dynamical onset of light-induced superconductivity in a strongly correlated, yet theoretically controllable system, to gain intuition and assess the validity of common approximations.
First, we consider (quasi-static) light-induced modifications of Hamiltonian parameters, that we model as interaction quenches across the superconducting phase transition.We find, that in contrast to BCS theory [83] and previous studies in the Hubbard model [84][85][86][87][88], it is not possible to induce superconductivity via interaction quenches in the YSYK model due to excessive heating of the nFl normal state.The light-induced cooling on the other hand, which we model as the relaxation of undercooled normal states to the superconducting ground state, is found to universally lead to superconductivity at late times.Providing a full characterization of the non-equilibrium relaxation, we identify an overdamped and oscillatory relaxation regime, and demonstrate that a coarse grained Ginzburg Landau theory can be effectively determined from microscopics to describe the emerging order parameter dynamic despite the strong correlations.
Model -The YSYK model describes N fermions randomly interacting with N dispersionless bosons and is defined by the Hamiltonian [50][51][52]55] Here c ( †) iα are fermionic ladder operators, while the canonical boson fields φ k and their conjugate momentum π k satisfy [φ k , π k ] = iδ kk .The electron-boson vertices g ij,k are uncorrelated, random all-to-all couplings which are sampled from a Gaussian-orthogonal ensemble (GOE) [89] with mean g ij,k = 0 and variance g ij,k g * ij,k = g 2 .The normal state at low temperatures, is a strongly interacting quantum critical non-Fermi liquid (nFl) with power law fermionic G(ω) ∼ |ω| 2∆−1 and bosonic D(ω) ∼ |ω| 1−4∆ spectral functions and universal ∆ ≈ 0.42 [50,53].Allowing for pairing (U(1) symmetry breaking), superconductivity emerges for all coupling strengths g.It supersedes the quantum critical nFl and is the final low temperature state.The superconducting condensate is a strongly interacting Cooper-pair fluid, unlike the noninteracting pair fluid in BCS theory [90].The YSYK model hence provides a paradigmatic example for strong coupling unconventional superconductivity.
Non-equilibrium formalism -In the thermodynamic limit N → ∞ the model becomes self-averaging [91] and exactly solvable due to the vanishing of vertex corrections.This allows the derivation of a closed system of self-consistent equations for fermion, boson and anomalous Green's functions (GF) on the Keldysh contour.Focusing on s-wave superconductivity, we introduce Nambu vectors ψ i = (c i↑ , c † i↓ ) T and the disorderaveraged Nambu-Keldysh G(t, t ) and bosonic D(t, t ) GF's defined by where T C denotes contour-time ordering along the Keldysh contour C [92].The time evolution of the Greens functions is governed by the Kadanoff-Baym equations (KBe) with exact fermionic and bosonic self-energies [93] Equations ( 4) to (7) represent the exact solution of the YSYK model, and they are structurally similar to the famous Migdal-Eliashberg equations of superconductivity [94], but on the Keldysh contour.Providing initial conditions in the lower quadrant of the two-time plane (t, t < 0), we start the non-equilibrium protocol at t = 0 and follow the exact evolution of the system by integrating the resulting equations numerically; see Supplementary Material for details [93].
After the full integration of the KBe's, the dynamics of the superconducting order parameter can be directly extracted from the GF's and provides a measure for the pairing strength of the system.Using generalized fluctuation dissipation relations (gFDR), we can further extract effective temperatures and nonequilibrium occupation functions.To that end, we transform the GF's to center of mass coordinate T = (t+t )/2 and relative coordinate τ = t − t and Fourier transform the latter (Wigner transformation).Writing the gFDR with G R/K the normal state electronic GF's, defines the electronic non-equilibrium occupation function n(T , ω), which in thermal equilibrium is a Fermi-distribution n f = (e βω + 1) −1 .In non-equilibrium, we can extract an effective temperature T eff (T ) from Eq. ( 9), provided that n(T , ω) resembles n f at low energies [73].To study the dynamical transition from a normal into a superconducting state, we provide an explicit U(1) symmetry breaking at t = 0, which is implemented as The symmetry breaking strength α, provides an infinitesimal seed to start the relaxation process and is chosen as to not affect the final state after relaxation [93].
Dressed Hamiltonians -First, we consider lightinduced modifications of Hamiltonian parameters, that we model as interaction quenches across the superconducting phase transition.The quenches are implemented as a time dependent fermion-boson coupling that modifies the self-energies Eqs. ( 6) and (7) via g 2 → g 2 f (t)f (t ).Starting from a normal state at g i with temperature T i > T c (g i ) we quench into the superconducting phase g i → g f such that T i < T c (g f ).The exact dynamics following the interaction quench is illustrated in Fig. 2. Shortly after the quench, the system is non-thermal and observables show an oscillatory dynamics [Fig.2(c)], with amplitudes that increase in the strong coupling regime, similar to observations made for quenches in the Holstein model [26].At larger times, the system universally relaxes towards a thermalized equilibrium state.The final temperature is determined from matching the instant energy density (t) = H YSYK (t) with that of a thermal ensemble . . .Tinst by solving the self-consistency relation (t) = H(t) Tinst(t) , hence defining the instant temperature T inst (t).For the current protocol T inst (t) is constant for t ≷ 0 respectively and agrees with T eff in the long time limit [Fig.2(c)], indicating thermalization [73].This we confirm by also analyzing other thermodynamic observables and non-equilibrium occupation functions.The final temperature after the quench with weakly g i dependent exponent.We universally find showing that no (thermal) superconductivity remains in the long time limit.
Indeed, we generally find that it is not possible to induce superconductivity by interaction quenches in the YSYK model, neither transiently nor at late times.The dynamics of the superconducting order parameter scales as ∆(t) ∼ O(α) for all times [Fig.2(b)], so that for the physical limit α → 0, no superconducting fluctuations remain.Instead, the non-Fermi liquid normal state of the YSYK models rapidly heats following the interaction quench, prohibiting the formation of cooper pairs.Our surprising results remain true, when considering quenches with finite, linear interaction ramp [93].This is distinct from BCS theory, where one finds a coherently oscillating, yet finite order parameter following interaction quenches [83].Further, it deviates from observations made for similar protocols in Hubbard models, treated perturbatively or with DMFT, where a transiently ordered state can emerge following an interaction quench [84][85][86][87][88].
Dressed Statistics -Next, we consider the emergence of superconductivity due to light-induced cooling.This we model as the dynamical relaxation from an undercooled YSYK normal state, into the superconducting ground state, following an explicit U(1) symmetry breaking; see Eq. ( 10).An undercooled state, that is an unordered ∆ = 0 state with T < T c , is unstable to such perturbations and its exact relaxation dynamics is shown in Fig. 3. Similar to the quench protocol above, we find that the system universally thermalizes in the long time limit, but now with T eff (∞) < T c , leading to emergent superconducting states.
The exact dynamics of the superconducting order parameter in shown in Fig. 3(a).While the thermalization time t * (circles) shows a strong dependence on α, the overall dynamics as well as the final state ∆(t > t * ) are independent of the symmetry breaking strength [93].The latter agrees with the thermal gap function evaluated for T eff (t → ∞) (color matched arrows), further illustrating thermalization at late times.
The early time dynamics (t < t * ; pentagon), is uni- (a) Order parameter ∆(t) [Eq.( 8)] for different symmetry breaking strengths α [Eq.( 10)] and temperatures.Close to Tc (blue) the dynamics is overdamped while for T Tc (red) it becomes underdamped and oscillatory.Both regimes are well described by Ginzburg Landau theory (dashed lines).The color matched arrows mark the thermal ∆ evaluated at T eff (∞), indicating thermalization.(b) Non-equilibrium occupation function n(T , ω) [Eq.( 9)] for different center mass times T indicated as dots in the effective temperature inset.Close to Tc the system is locally thermal and follows a hydrodynamic time evolution.For T Tc the early times are non-thermal and only at late times does the system thermalize.
versally determined by an exponentially growing order parameter with α independent realxation rate Γ(T ), while at later times (t > t * ) the dynamics in the weak (T T c ) and strong (T T c ) undercooling regime is qualitatively different.Close to T c it is underdamped, in the sense that δ(t) = ∆(t) − ∆(∞) has no zero crossing, while for T T c we observe an oscillatory dynamics (Fig. 3(a); inset).This difference can be understood intuitively.Close to T c , the superconducting and normal state do not deviate strongly such that the 'effective force', driving the relaxation towards the ordered state is small compared to internal relaxation rates.Locally the system is in thermal equilibrium and the time evolution is hydrody- namic, with non-equilibrium occupation functions resembling Fermi-distributions (Fig. 3(b); upper panel).For T T c the 'effective force' becomes large compared to the internal relaxation rate so that the system is nonthermal at short and intermediate times and only thermalizes in the long time limit (Fig. 3(b); lower panel).
The superconducting order parameter dynamics can be reproduced by time dependent Ginzburg-Landau theory (GL) with F GL = c 1 ∆ 2 + c 2 ∆ 4 and dynamical equation For weak undercooling, the dynamics is described by overdamped GL (η → 0), while the oscillatory case requires both finite η, δ to reproduce the exact data [dashed lines in Fig. 3(a)].BCS theory corresponds to the limit δ → 0, implying that the interactions effectively introduce a damping for the order parameter dynamics [93].
The temperature dependence of the relaxation rate Γ(T ) is a defining characteristic of the relaxation process.For T T c it agrees well with the overdamped GL prediction Γ(T ) ∼ |T − T c |, while it saturates for T T c (Fig. 4).The reason for the different regimes is again intuitively understood via the 'effective force' argument.Close to T c the speed of the relaxation dynamics is limited by the force, while for T T c internal relaxation processes determine the relaxation timescale.
Generically, the relaxation rate vanishes as Γ(T ) ∼ |T − T c | γ with critical exponent γ, a phenomenon known as critical slowing down [90,95].From our exact dynamics we find γ = 1, which despite the strong correlations and nFl normal state, agrees with time dependent GL theory.
While the general structure of the relaxation is similar for different couplings g, we find that the dynamics is faster at strong coupling ∂ g Γ(T ) > 0, but like the critical temperature saturates for g 2 [50,93].Further, the oscillatory order parameter dynamics is most pronounced for g 1, while the amplitude of the oscillations becomes strongly suppressed at stronger interactions [93].
Discussion -We studied the dynamical onset of superconductivity in a paradigmatic, solvable toy model for strongly correlated superconductivity by analyzing two dynamical protocols.These are motivated by the mechanism of light induced superconductivity; that is the dressing of Hamiltonian parameters and light-induced cooling.
Modeling the (quasi-static) dressing of Hamiltonian parameters, as interaction quenches across the superconducting phase transition we found that no superconductivity is induced due to excessive heating of the non-Fermi liquid normal state.This contrasts results obtained for BCS theory and DMFT studies in the Hubbard model.
Imitating light-induced cooling by analyzing the dynamical relaxation from an undercooled YSYK normal state, we universally observed superconductivity at late times.Depending on the strength of the undercooling, we found an overdamped or oscillatory relaxation dynamics towards these states.The resulting order parameter dynamics could be reproduced by time dependent Ginzburg Landau theory.
An exciting extension of our work is the investigation of periodically and parametrically driven YSYK models, as this protocol has been suggested to be able to induce superconductivity [19][20][21]27].Our studies of the driven YSYK model revealed that driving leads to excessive heating, prohibiting the formation of superconductivity, but the additional coupling to a thermal bath [66,[75][76][77]81] might mitigate this effect and lead to Floquet steady-states [96,97] that could host non-thermal light-induced superconductivity.This promises to further enhance our understanding of the non-equilibrium superconductivity in strongly correlated systems.
Note Added: After the completion of this work we learned about an independent study of the nonequilibrium dynamics of the YSYK normal state, following a quench coupling to an external reservoir [98].
Supplementary Material for: Dynamical Onset of Light-Induced Unconventional Superconductivitya Yukawa-Sachdev-Ye-Kitaev study Lukas Grunwald, 1, 2, * Giacomo Passetti, 1 and Dante M. Kennes We derive the exact solution of the YSYK model on the Keldysh contour at large N .For completeness, we consider a slight generalization to the model introduced in the main text by considering spin non-diagonal interactions as well as a finite chemical potential.The Hamiltonian reads [1][2][3] where σ x αβ , x = 0, . . ., 3 denotes the Pauli-matrices and f (t) a generic time dependent function.We will derive the general solution of Eq. ( 1) and restrict to our specific model with x, µ = 0 in the end.The starting point is the path-integral on the Keldysh-contour [4,5].For the YSYK model, the associated action appearing in the generating functional We imply a summation over repeated Greek indices and use the shorthand C ( ) (. . . ) = C dt ( ) (. . . ) for integrals over the Keldysh-contour C = C + ∪ C − , which consists of forward branch C + , running from t 0 → ∞ and backward branch C − running back from ∞ → t 0 .We assume no initial correlations and set t 0 → −∞ [4,5], which is sufficient for the non-equilibrium protocols we analyze in this work (see below).
arXiv:2307.09935v1 [cond-mat.supr-con]19 Jul 2023 From the non-interacting action S 0 , we can directly extract the free fermionic and bosonic contour propagators with C δ(t, t ) = 1.Next we execute the disorder average over the interaction vertices g ij,k ∈ GOE (Gaussian-Orthogonal-Ensemble [6]), where for each bosonic index k, g ij,k represents a random matrix in the indices i, j.The GOE is an ensemble of symmetric random matrices g ij (suppress bosonic index k), which is defined over the mean g ij = 0 and the variance g ij g km = (δ ik δ jm + δ im δ jk )g 2 of the matrices.The probability density and measure read so that the evaluation of the disorder average reduces to calculating Gaussian integrals (no replica trick is needed in the Keldysh formalism [4]).The resulting disorder averaged action reads Note the appearance of anomalous terms ∼ ψ α (t)ψ β (t ), which give rise to superconductivity and which are entirely due to g ij,k ∈ GOE being real.For complex g ij,k ∈ GUE (Gaussian-Unitary-Ensemble), like in the SYK model [7], we would obtain Eq. ( 7), but without the anomalous contributions.Next, we introduce the site averaged contour Greens functions and corresponding self-energies (Lagrange-multipliers) as dynamical mean-fields via fat-unities The signs and prefactors are chosen in anticipation of obtaining Dyson equations (see below).Using the dynamical mean-fields, one finds for the interaction part (∼ g 2 ) of the averaged action where the tr[. ..] is to be taken over the spin indices.The ψ, ψ and φ fields are now decoupled and can be integrated out separately after taking into account the self-energy contributions from the fat unities.This process is straight forward in the bosonic case In the fermionic sector, the action contains both normal ∼ ψα ψ β and anomal ∼ ψ α ψ β contributions, the latter stemming from the fat unities Eqs. ( 10) and (11).We introduce Nambu vectors ψ i = (ψ i,↑ , ψ i,↓ , ψi,↑ , ψi,↓ ) to rewrite the remaining fermionic action as a Grassman-Gaussian integral with Here we introduced the Nambu-Keldysh self-energy Σ(t, t ) and free-propagator G0 (t, t ) and understand the transpose as A αβ (t, t ) T = A βα (t , t).Note that the integral-measure only contains ψ but not ψ † , since the latter is just a rearrangement of the former.Hence, Eq. ( 14) is not an ordinary Grassman-Gaussian integral.For M ∈ Skew n and Grassman fields θ = (θ 1 , . . ., θ n ) integrals of this type equate to [8] Equation ( 14) can be brought into exactly this form by permuting elements in [. . .], and we always have even n because of the 4-component Nambu vectors.Using Eq. ( 16) and reversing the permutation in the log Pf[. . .] term, we find The resulting GΣ-action only depends on the dynamical mean fields G, Σ, F, φ, F , φ and reads where the Tr(. ..) is to be taken over all indices and times and where we arranged fermionic propagators G, F, F into the Nambu-Keldysh Greens function The shorthand (. .)abbreviates convolutions over the support of the functions involved, in this case the Keldyshcontour C as well as the matrix structure of the GF's.
The action is globally multiplied by N , such that in the large N limit, the exact solution can be obtained by evaluating the variation δS GΣ .This yields equations of motion for the dynamical mean-fields.Similar to derivatives of log det(. ..) terms, (functional) derivatives of Pfaffians can be calculated as [8] 1 Pf(A) Varying the action with respect to the self-energies Σ, φ, φ and Π, we obtain the contour Dyson equations [4,5] for the Numbu-Keldysh GF's and the bosonic propagator respectively During the evaluation of the former we used [ G−1 The equations for the self-energies follow by varying S GΣ with respect to the Greens functions G, F, F , D, resulting in Equations ( 22) to (25) determine the exact self-energies of the YSYK model on the Keldysh contour (t, t ∈ C).We additionally obtained them using a diagrammatic derivation as a consistency check.Note that this system of equations only depends on the site averaged GF's and self-energies.This illustrates the self-averaging property, since the exact solution has no reference to the site indices of the original theory.The hence derived equations contain both singlet and triplet pairing components.As in previous works on the YSYK model, we focus on the x = 0 case, where the model has a SU(2) symmetry.The presence of triplet superconductivity would spontaneously break this symmetry.Following [1,2], we assume no spontaneous SU(2) symmetry breaking and will thus focus on singlet anomalous propagators.On the Keldysh contour this implies F . Further, due to the SU(2) symmetry, the normal GF and self-energy become spin independent, i.e.
Henceforth, we can work with reduced Nambu-vectors ψ i = (c i↑ , c † i↓ ) T and the associated contour GF's and self-energies (t, t ∈ C) Equations ( 22) to ( 25) then simplify to the equations from the main text (t, t ∈ C) where In order to solve Eqs. ( 28) to (31), we need to analytically continue them to real-times, and physical Greens functions G ≷ , D ≷ or G R/A/K , D R/A/K .The latter are obtained by fixing the position of the contour times t, t ∈ C on the Keldysh contour t, t ∈ C ± , viz.(similar for D Using the Langreth theorems [5,9], Eqs. ( 28) to (31) translate to a system of non-linear Volterra integro-differential equations for G ≷ and D ≷ -the so called Kadanoff-Baym equations (KBe's) where (A B)(t, t ) = ∞ t0 A(t, τ )B(τ, t ) dτ and matrix multiplications are implied.Due to the causality structure of the GF's, these equations can be integrated as an initial value problem (see next section).
In this work the initial conditions are fixed in the lower quadrant of the two-time plane t, t < 0 as the equilibrium solution of Eqs.(37) to (42), i.e. for f (t) ≡ const.In equilibrium the GF's only depend on time differences, allowing us to map the differential equations' [Eqs.(37) to (40)] onto algebraic equations in frequency space.There self-energies and Greens functions are connected by Further G R (ω), G K (ω) and D R , D K are related by fluctuation dissipation theorems respectively [5].Hence, Eq. ( 43) together with Eqs. ( 40) and ( 41) provides a closed, algebraic system of equations for the equilibrium GF's.It can be solved numerically using a fixed point iteration in combination with numerical Fourier transforms, which we evaluate using an interpolation based FFT-scheme [10].The fixed point iteration is analog to the SYK model and previous studies of the YSYK model; see [1,7,11,12] and references therein.

INTEGRATION OF KADANOFF-BAYM EQUATIONS
General Idea -In this section we discuss the general procedure -time stepping scheme -for solving Kadanoff-Baym equations.Since the general time-stepping is the same for bosons and fermions, we explain our scheme using the generic KBe's where h(t) denotes the non-interacting Hamiltonian and the matrix structure is left implicit.The bosonic equations encountered in the last section are equivalent, apart from ∂ t → ∂ 2 t .We will comment on the implications of this below.
Because of the causality structure of G R/A and Σ R/A , the r.h.s of the KBe's for G ≷ (t, t ) only depends on times smaller than max(t, t ).Thus, providing initial conditions for G ≷ (t, t ) in some interval t, t ∈ (t 0 , t 1 ), the KBe's can be integrated as an initial value problem [13].Due to them being integro-differential equations, each time-step depends on all previous times via the history integrals A time-step is defined as advancing δt forward in time, naturally enforcing a time discretization t n = t 0 + nδt.The different KBe's correspond to different directions in the (t, t )-plane and in order to perform a time-step, they need to be used in combination.This we illustrate in Fig. 1.Equation (44) propagates the Greens functions down [green, G ≷ (t 1 + δt, t ∈ (t 0 , t 1 ))] and Eq. ( 45) propagates to the right [blue, G ≷ (t ∈ (t 0 , t 1 ), t 1 + δt)].For the diagonal propagation [orange, G ≷ (t 1 + δt, t 1 + δt)] we subtract Eqs. ( 44) and ( 45) to obtain Equations ( 44) to (46) allow the calculation of all points needed for further time steps (dotted line in Fig. 1).They can thus be used to propagate the GF's in the entire (t, t )-plane, given some initial conditions.The full integration of the KBe's yields G ≷ (t, t ) for t, t ∈ (t 0 , t max ).
The bosonic equations [Eqs.(39) and (40)] have exactly the same analytical structure as Eqs.( 44) and ( 45), and we can hence use the same scheme as outlined above with one caveat.Because of ∂ 2 t it is not possible to generate an equation for the diagonal, comparable to Eq. ( 46), which only contains a second order derivative.Instead, the diagonal is determined as a combination of propagation to the right and down [Eqs.(39)   KBe time stepping scheme.Given G > (t, t ), G < (t, t ) on an initial interval t, t ∈ (t0, t1), here illustrated as gray dots, we can obtain the solution at all times by time-stepping with the KBe's.One time step comprises the calculation of the points indicated by the green G ≷ (t1 + δt, t ∈ (t0, t1)), blue G ≷ (t ∈ (t0, t1), t1 +δt), and orange G ≷ (t1 +δt, t1 +δt) arrows, corresponding to the propagation with ( 44), ( 45), (46) respectively.

QUENCH WITH FINITE RAMPING TIME
The inability to induce superconductivity following a quench across the phase boundary also persists when considering a finite ramp for the interaction quench.We study this for quenches from g i = 0.75 to g f = 1.0,where the finite ramping time τ is implemented as a linear ramp in g ij,k f (t) with The exact dynamics for ramping times τ = 0, . . ., 7, where τ = 0 corresponds to instant quenches, is shown in Fig. 3. Analog to the instant quenches, the system relaxes to the new equilibrium state with an oscillatory dynamics.While it is possible to decrease the final temperature after the quench by increasing the ramp time τ , it does not seem to be possible to achieve T eff < T c , at least for the coupling considered here.The final temperature saturates for τ > 6 [orange and red curves in Fig. 3(b)].We universally find ∆(t) ∼ α [saturated and desaturated; Fig. 3(a)], implying that no transient superconductivity is induced.
At the largest ramping times τ ≥ 6, the system is 'locally thermal' during and after the quench, as we illustrate in Fig. 4, which shows the non-equilibrium occupation function n(T , ω) for early and intermediate times.While for instant quenches and short ramps, the occupation function n(T , ω) shows strong non-eq features (n > 1 or n < 0), these gradually decrease and vanish as the ramping times τ increases.At τ ≥ 6, the occupation function resembles as Fermi-distribution with time dependent temperature during the time evolution.

THERMALIZATION TIME AND CRITICAL SLOWING DOWN
In the main text we remarked on the strong α dependence of the thermalization time t * .Its functional dependence given by with temperature dependent constants C, t * ,0 > 0, as we illustrate in Fig. 5(a) (note the logarithmic x-axis).The turning-point time t * , up to which the order parameter shows an exponential behavior, has exactly the same slope, but a different constant t * ,0 → t * ,0 .The dynamical onset of superconductivity is hence exponentially activated in the strength of the symmetry breaking seed.This is readily understood via the exponential dynamics at early times ∆(t) ∼ exp(Γt).The system inherently produces an exponentially growing seed, a short time after the initial symmetry breaking occurs.Hence, only providing an exponentially bigger initial seed can modify this self-propelled dynamics in a meaningful way.Another interesting question is the parametric dependence of t * on the initial temperature T of the normal state, Close to Tc (normal state; orange) we observe a critical slowing down, i.e. t * (T ) ∼ (T − Tc) −y for which we extract the α independent exponent y by fitting T ≥ 0.03.For temperatures T Tc, the exact t * (T ) deviates from this power-law form and saturates.Close to the superconducting transition temperature (Tc ≈ 1/15.1) the relaxed ∆ depends strongly on α, while deep within the superconducting phase it is much less sensitive.Since we want the final state to be independent of α on the scale of the plots, we must decrease it as we approach Tc.This, combined with the critical slowing down (see above), makes the numerical study of undercooled quenches close to Tc challenging.
which we analyze in Fig. 5(b).Here, we observe a critical slowing down of the relaxation dynamics as we approach the critical temperature T c (orange region), i.e. lim T →Tc t * (T ) = ∞, which is simply a different manifestation of lim T →Tc Γ(t) = 0 discussed in the main text [18].This feature is very intuitive, since the closer we get to T c , the smaller the 'effective force' becomes, that pushes the system into a superconducting ground state.Right at T c , the system is critical so that the restoring force vanishes and the thermalization time diverges.
While the dynamics is exponentially activated in the symmetry breaking strength, the thermalized value ∆(t > t * ) as well as the relaxation rate Γ(T ) are α-independent, for small enough α.That is desirable, since α is supposed to represent an infinitesimally small perturbation to the system that kick-starts our relaxation process, but does not affect the final state after relaxation (at least on the scale of the plots).This in turn means, that as we approach T c , α needs to be become smaller for the final state to be α-independent, since the thermal equilibrium value of ∆ continuously decreases to zero for T → T c .We illustrate this in Fig. 6.This in combination with the critical slowing down makes the analysis for T ∼ T c numerically challenging.

UNDERCOOLING AT DIFFERENT COUPLING STRENGTHS
Here we comment on the undercooling dynamics at different coupling strengths g.The relaxation rate Γ(T ) extracted from the early time dynamics, is shown in Fig. 7 for a variety of different couplings (note the rescaled y-axis for better visibility).Close to T c the prediction from overdamped Ginzburg-Landau theory (dashed) agrees well with the exact results, while the curves universally deviate from this behavior for stronger undercooling, where they show a saturation behavior.We observe that the dynamics is generally faster at stronger coupling.Further, we find that the oscillatory order parameter dynamics is most pronounced at weak and intermediate coupling g 1, while the amplitude of the oscillations becomes strongly suppressed at strong coupling.

GINZBURG-LANDAU PARAMETER FITS
We analyze the temperature and coupling dependence of the Ginzburg Landau theory (GL) parameters, mentioned in the main text.We obtained these from fits of the Ginzburg-Landau ∆ GL (t) to the exact order parameter dynamics.Our Ginzburg Landau theory is defined by the Free energy F GL = a∆ 2 GL + b∆ 4 GL and dynamical equation which for η = 0 has to be solved numerically due to the non-linearity.Generally we find that the fitted parameters are very sensitive to small deviations in ∆(t), making an accurate numerical extraction difficult.The extracted values lead to ∆ GL (t)-curves, that on the scale of the plots, coincide with the exact ∆(t) dynamics, but one shouldn't put too much emphasis on the exact numerical values obtained here.This procedure can nonetheless be used to determine the general trends for the GL parameters.The fits are obtained using the BFGS routine, implemented in Optim.jl[19], on ∆(t) and ∂ t ∆(t) together with DifferentialEquations.jl[20] for the numerical solution of the ODE's.
The extracted parameters are shown in Fig. 8 for different couplings g.Note that generally a/η → 0 as T → T c as expected [Fig.8(a)].In this limit δ/η → ∞ [Fig.8(c)], implying that the dynamics becomes overdamped.We independently verified this by fitting GL theory with η = 0 to the exact numerical results close to T c .This leads to excellent agreement with the exact data.
We generally find larger damping δ in the stronger coupling regime, with the rough scaling δ ∼ O(g 2 ); curves in the figure become of the same order of magnitude when applying this rescaling.Hence, in order to observe the oscillatory order parameter dynamics at strong coupling, a stronger undercooling is necessary.

Figure 1 .
Figure 1.Cartoon of light irradiated YSYK model.Irradiation of the YSYK model consisting of randomly interacting fermions (blue) and bosons (red) could dynamically induce superconductivity due to two effective mechanisms.(a) Dressing of coupling constants; here modeled as interaction quenches.(b) Cooling of the system; here modeled as the relaxation of an undercooled normal state into the superconducting ground state.

T c /ω 0 ≈ 1 / 15 . 1 T /ω 0 = 1/ 20 T /ω 0 = 1 Figure 3 .
Figure 3. Dynamical relaxation from undercooled normal state into superconductivity at g = 1.(a) Order parameter ∆(t) [Eq.(8)] for different symmetry breaking strengths α [Eq.(10)] and temperatures.Close to Tc (blue) the dynamics is overdamped while for T Tc (red) it becomes underdamped and oscillatory.Both regimes are well described by Ginzburg Landau theory (dashed lines).The color matched arrows mark the thermal ∆ evaluated at T eff (∞), indicating thermalization.(b) Non-equilibrium occupation function n(T , ω) [Eq.(9)] for different center mass times T indicated as dots in the effective temperature inset.Close to Tc the system is locally thermal and follows a hydrodynamic time evolution.For T Tc the early times are non-thermal and only at late times does the system thermalize.

Figure 4 .
Figure 4. Relaxation rate Γ(T ) from the early time dynamics ∆(t) ∼ e Γ(T )t at g = 1.Γ extracted at different symmetry breaking strengths α (different shapes) becomes indistinguishable for α < 10 −5 illustrating the α independence.Close to Tc we find Γ ∼ |T − Tc| in agreement with the overdamped Ginzburg-Landau prediction (dashed line).Deep within the superconducting phase T Tc the rate saturates.Note limT →Tc Γ(T ) = 0, indicating the critical slowing down.

Figure 3 .
Figure 3. Quenches across superconducting transition with finite ramping time for g = 0.75 → 1.0 with α = 10 −5 (full lines) and α = 10 −6 (desaturated lines; only distinguishable in the left panel).(a) The superconducting order parameter ∆(t) stays on the order of the symmetry breaking ∆(t) ∼ α, implying that no superconductivity is induced.(b) After initial oscillations at early times, the effective temperature T eff agrees with the instant temperature Tinst (color matched arrows), indicating thermalization at late times.Universally T eff (t → ∞) > Tc.(c) The bosonic ocuupation n boson (t) follows an analog behavior to the effective temperature, viz.oscillations at early times that relax in the long time limit.

Figure 3 (
a) plots the gap function ∆(t), Fig. 3(b) the effective temperature and Fig. 3(c) the bosonic occupation.A further increase in ramping time does not change the dynamics in any meaningful way, nor decrease the final effective temperature.

Figure 4 .
Figure 4. Non-equilibrium ocuupation function n(T , ω) after quenches across superconducting transition at fixed α = 10 −5 for instant (blue) quenches and finite ramps.While instant quenches lead to strong non-equilibrium features in the non-equilibrium occupation function (n > 1, n < 0), these features gradually vanish as the ramping time increases.

Figure 6 .
Figure 6.Relaxed gap function ∆(α)|t * at thermalization time t * , for g = 1.0 for different temperatures.Close to the superconducting transition temperature (Tc ≈ 1/15.1) the relaxed ∆ depends strongly on α, while deep within the superconducting phase it is much less sensitive.Since we want the final state to be independent of α on the scale of the plots, we must decrease it as we approach Tc.This, combined with the critical slowing down (see above), makes the numerical study of undercooled quenches close to Tc challenging.

Figure 7 .
Figure 7. Relaxation rate Γ(T ) extracted from early time dynamics ∆(t) ∼ e Γ(T )t at α = 10 −5 and different couplings g.The main panel shows g • Γ(T ) (for better visibility) together with linear fits (dashed lines) Γ(T ) ∼ c1|T − Tc|/Tc + c2 to the exact data (dots).The inset illustrates the g-dependence for the slope of the linear fit c1 = ∂T Γ(T ) T →Tc .Note that c2 ≈ 0, indicating the critical slowing down.The dynamics is generally faster at stronger coupling ∂gc1(g) > 0, but like the critical temperature, saturates for g 2.

Figure 8 .
Figure 8. Ginzburg-Landau parameters for undercooling protcol at different couplings g.Axis are rescaled with the coupling strength g and critical temperatures Tc for better visualization.While the exact numerical values are sensitive to numerical details, the general trends appear to be universal.
t t (t max , t max ) (t min , t min ) (t min , t max ) (t max , t min ) δt δt (t n , t m )