Transient superconductivity in three-dimensional Hubbard systems by combining matrix product states and self-consistent mean-field theory

We combine matrix-product-state (MPS) and mean-field (MF) methods to model the real-time evolution of a three-dimensional (3D) extended Hubbard system formed from one-dimensional (1D) chains arrayed in parallel with weak coupling in-between them. This approach allows us to treat much larger 3D systems of correlated fermions out-of-equilibrium over a much more extended real-time domain than previous numerical approaches. We deploy this technique to study the evolution of the system as its parameters are tuned from a charge-density wave phase into the superconducting regime, which allows us to investigate the formation of transient non-equilibrium superconductivity. In our ansatz, we use MPS solutions for chains as input for a self-consistent time-dependent MF scheme. In this way, the 3D problem is mapped onto an effective 1D Hamiltonian that allows us to use the MPS efficiently to perform the time evolution, and to measure the BCS order parameter as a function of time. Our results confirm previous findings for purely 1D systems that for such a scenario a transient superconducting state can occur.


Introduction
Superconductivity (SC) has remained a phenomenon of great interest to researchers ever since its discovery in 1911 by H. K. Onnes.Explaining SC in metals at low-temperature equilibrium was already a challenge, taking more than 40 years until the Bardeen-Cooper-Schrieffer (BCS) framework could explain it via a suitable MF theory.In the 1980s, SC at high critical temperatures T c [1][2][3][4] was discovered, which seemed not to be described by BCS theory.In fact, its theoretical description presents a still-ongoing challenge.It is believed that strongly correlated electron motion is the underlying reason for this type of SC state.Many-body models such as the  or the t-J-model [4,[10][11][12][13] have been investigated to study this question.In more recent developments, experiments claim to have detected transient, 1 light-induced SC states after pushing layered high-T c materials out-of-equilibrium in pump-probe setups, even above the equilibrium T c [14][15][16][17][18]. One interpretation of these results is via the concept of pre-formed pairs, namely that above T c up-and down-spin electrons are still paired, up to the so-called pseudogap (PG) temperature T PG .In this view, at temperature T with T c < T < T PG , the system cannot achieve SC order in equilibrium because the pre-formed pairs within each layer lack inter-layer phase coherence.Thus, the pump-pulse is deliberately designed to increase inter-layer coupling, and results in the observation of a seemingly SC response in the optical conductivity, but only as long as T < T PG [16,17].
There are two possible objections to this particular interpretation of the experimental literature: (1) The reliance of experiments on the time-dependent optical conductivity as a probe for nonequilibrium SC has been questioned.Paeckel et al. [19] recently showed that this measure lacks specificity for SC order and proposed alternative measurements, which would be better suited to detecting the onset of the SC state in the dynamically evolving system.The setup studied in the article of Paeckel et al. consists of a quench performed on a 1D extended Hubbard system at T = 0 using a MPS description.This MPS approach, while unbiased and highly accurate, is effectively restricted to 1D systems, especially when treating out-of-equilibrium dynamics.The question is thus if the findings of Paeckel et al. are specific to 1D, with its strong quantum and thermal fluctuations, or whether they also apply to higher dimensional systems.
(2) This leads to the second objection: quantitative or even qualitative theoretical understanding of the solid state high-T c materials evolving under a pump-pulse probe, or even a generalization of the setup of Paeckel et al. to a 3D system, has been lacking (even understanding the high-T c materials in equilibrium still presents an enormous challenge).That is, no one has been able to show even in principle whether the proposed scenario for interpreting many solid-state experiments -initially phase-incoherent pairs of fermions dynamically acquiring macroscopic coherence as inter-layer coupling is rapidly increased -is theoretically possible, and, if so, under which conditions and on which time scales.While some numerical work has been carried out towards this, [19][20][21], many basic questions about the mechanisms that could lead to dynamically induced SC remain open.To a substantial degree, this is due to the significant challenge of accurately capturing those system sizes and time scales required to observe any emerging SC states.The central question then becomes which class of algorithms has any chance of reaching those regimes, thereby delivering the currently-lacking theoretical insight into how dynamically induced SC develops and which microscopic conditions would be fundamentally required for this.On their own, even in 1D, MPS methods may require exponentially increasing resources as simulation time grows in order to maintain a set accuracy.This is due to the strong growth in bipartite entanglement in these systems with time: for MPS approaches to be efficient, this entanglement should not be too large.Furthermore, already for equilibrium calculations long-range interactions, which are needed to represent two-dimensional (2D) and 3D systems in 1D, increase the entanglement dramatically.Hence, the time evolution of generic 2D and 3D systems are entirely out of reach for any brute-force MPS-based approach, as they could not capture even the initial equilibrium state.
For such higher-dimensional systems, real-time non-equilibrium dynamical mean-field theory (DMFT) could be a powerful alternative approach [22,23].For this technique, one or a few lattice sites -the impurity or, respectively, the cluster -are retained explicitly, including all interactions of the original, infinitely-large lattice.In DMFT, the effect of this remainderlattice on the cluster is mimicked via a free-electron bath that is coupling to it.The parameters of this bath are fixed via self-consistency conditions.Solving these cluster-bath systems within this self-consistency constraint is typically achieved by applying quantum Monte Carlo (QMC) techniques in the real-time domain.These techniques suffer from a strong sign-problem, i.e., their numerical error grows exponentially as the cluster-size and the real-time domain, over which the simulation runs, are increased.In practice, a few sites and time scales on the order of the electron tunneling are accessible.Alternatively, MPS solvers can be used within such real-time non-equlibrium DMFT; however, due to the long-range tunneling in these systems between bath and cluster sites, and the strong growth of entanglement with time, these will also be limited to a few sites and short times.
This leads us to the scope of the present paper: with current methods it seems practically impossible to perform meaningful simulations of dynamically-induced SC in a 3D system.For MPS methods, the growth of entanglement with system size and simulation time is immediately prohibitive; and for non-equilibrium real-time DMFT, the large clusters and long times required to resolve the onset of a potentially weak SC order appear to be out of reach.
However, as we demonstrate in the following, it is possible to make a specific category of 3D fermionic systems amenable to real-time evolution via MPS techniques using a static MF ansatz, by exploiting certain gaps in the excitation spectrum of these cases.In this way, it is possible to capture strong correlations by using a MPS, and treat the full 3D system more accurately than by applying a pure MF treatment.We show explicitly that our approach can describe the dynamical emergence of SC from a state that is not SC to begin with in a PGlike system, i.e., a system that has singlet-pairing of spinful fermions but with no initial phase coherence of these pairs.In this, our method can address exactly the setup thought to be at the heart of the solid state experiments, and for which we present here the first efficient many-body numerics.Several of the authors have previously developed related approaches for systems in equilibrium [24], reproducing physical behavior correctly at zero and finite temperature compared to appropriate QMC simulations [25,26], with the overestimation of SC properties due to the MF-component of our technique a constant one, and modest at that in 3D systems.In these approaches, weakly coupled chains or ladders are stacked up into 3D cubic systems, which thus have anisotropic tunneling -much stronger inside the 1D systems than in-between them in the two orthogonal directions.For the case of fermions, the MF approximation can be introduced if each of the constituent 1D systems has a gapped energy sector, such as a spin gap, and thus single-fermion tunneling in-between 1D systems is suppressed in this weak-coupling regime [25].Just as for the equilibrium case [25], it is this crucial ingredient that allows us to perform real-time evolution for a much higher number of correlated sites than non-equilibrium real-time DMFT, as well as extending the real-time domain enough to perform a meaningful simulation of the dynamically-induced SC in a 3D system.Within this wellbehaved domain, we apply our real-time MPS+MF technique to study the time evolution of the BCS order parameter after a fast ramp of the system from an insulating starting state into a parameter regime where the system would be SC in equilibrium.As a consequence, we observe the onset of a non-equilibrium SC state.
We note that the present work is focused purely on the algorithmic challenge of numerically studying such state, which entails simulating many-body non-equilibrium systems of fermions for system sizes and time scales that were previously inaccessible.As to how the interrelated questions of validating our approach with alternative methods and applying it to concrete experiments could be answered, we provide an outlook to that in the Conclusions.Further, as shown in Sec. 3 and Fig. 2, the basic algorithm is completely generic, i.e., it is not predicated on the dynamics of the system being of a specific type.Thus, our approach can be immediately generalized to quasi 1D systems with, e.g., explicit periodic driving (i.e., Floquettype dynamics [27,28]), or non-unitary dynamics as encountered in open quantum systems (quantum trajectories [29] or full master-equation [30] for Markovian baths, hierarchy of pure states (HOPS) [31] or projected purified MPS (ppMPS) enabled quantum trajectories [32] for non-Markovian ones); MPS-approaches have extensive track records in all these domains, but each of them would represent separate projects of their own and are thus not the subject of the present work.
The paper is structured as follows: in Sec. 2, we recapitulate the MF ansatz for weakly coupled Hubbard chains used in equilibrium, developed originally in [25].In Sec. 3, we introduce the extension to a self-consistent time-dependent MPS+MF scheme to study the time evolution of a 3D extended Hubbard system, which consists of weakly coupled chains.In Sec. 4, we present our results for the BCS order parameter and a detailed discussion of the convergence behavior of the method when treating 3D arrays formed from chains, each up to L = 30 lattice sites long.The time evolution of the SC order parameter shows indeed that in both finite systems as well as the thermodynamic limit a transient SC state can be entered.We further analyze the dependence of our results on the parameters of the simulations.In Sec. 5 we conclude and provide an outlook as to how the technique established in the present work could be validated and applied to experiments.The appendices discuss further details on the method at equilibrium, as well as further details of the simulations out-of-equilibrium.

Mapping of the 3D system onto a 1D self-consistent chain
As we aim to describe a 3D model system with a method that is mainly suitable for 1D, namely MPS, we first need to identify a class of 3D models amenable to mapping onto an effective 1D description.Following the work of Bollmark et al. [25,26], we focus on 3D systems constructed out of gapped 1D fermions.We arrange these 1D systems, which extend in the x-direction, in parallel into a square array in the ŷ − ẑ-plane, forming effectively a cubic lattice.We choose fermion tunneling to be anisotropic in this lattice, denoted by t ⊥ in the ŷand ẑ-directions.Adapting from Bollmark et al. [25], we choose an extended Hubbard chain as the 1D building block.The Hamiltonian constructed in this manner is illustrated in Fig. 1 and is given by Figure 1: Two-dimensional cross-section of the three dimensional model.For the sake of clarity, the 3D extension of the system out of the plane is not shown here.
Each box denotes a lattice site.The sites are coupled into chains in x-direction, as indicated by the thick lines between the boxes.Furthermore, all chains are weakly coupled by the transverse hopping t ⊥ .This way, we obtain an extension in ŷ and ẑ-direction. and Here, ĉ † n,R i ,σ and ĉn,R i ,σ denote the fermionic creation and annihilation operators on site n and for spin σ on a chain that is labeled by the 2D vector R i in the ŷ − ẑ-plane.They obey the anticommutation relations The indices i and j stand for different combinations of n, R i , and σ.The operator nn,R i ,σ = ĉ † n,R i ,σ ĉn,R i ,σ is the particle number operator for the corresponding site, chain, and spin.We use open boundary conditions and include a term for the chemical potential µ.The latter allows us to control the number of particles in the system.
The only non-1D term is the transverse hopping Ĥ⊥ .We are able to eliminate the beyond-1D nature of this term through a combination of perturbation theory on the transverse hopping and a MF decoupling of adjacent 1D systems.In the following we briefly recap the key steps, a detailed derivation of this approach can be found in the publication of Bollmark et al. [25].
Since we are interested in a model system for SC, we specify U < 0 in the chain-Hamiltonian Eq. ( 3).This negative-U term gives rise to pairing of opposite-spin fermions already in isolated systems at t ⊥ = 0.This is expressed by the finite spin gap ∆E s and a finite pairing energy ∆E p of these isolated chains, defined as follows: Here, E 0 (S z , N ) denotes the ground-state energy of Hamiltonian Ĥ0 for a single chain-index at total spin S z and total number of fermions N .Thus, ∆E s and ∆E p represent the minimal energy required for flipping a spin inside a chain and for breaking up a pair on a chain by moving one constituent to another chain in the full 3D system, respectively.From the definitions, it is easy to see that ∆E s ≤ ∆E p , and for our specific choice of 1D systems ∆E s = ∆E p .As outlined in the following, ∆E p becomes important in the actual numerical routine, directly entering the effective Hamiltonian Eq. ( 13).In practice, we can determine ∆E p from a single chain via an extrapolation in the system size L → ∞.
To carry out the second-order perturbation theory in Ĥ⊥ -specifically in t ⊥ /∆E p -we follow [33].We sort the eigenenergies E i,α of Ĥ0 ( Ĥ0 |i, α〉 = E i,α |i, α〉) into manifolds.The lowest-energy manifold with E i,α=0 , corresponds to those states in which each 1D system in the array is balanced between up-and down-spins and thus has S z = 0, and i indexes the states within this manifold.That is, in this manifold there are no broken pairs.The highenergy manifold E i,α=1 is at least ∆E p above the low-energy manifold, corresponding to excited states with at least one broken pair, i.e., where the pair-constituents have moved onto separate chains.In the perturbative regime, we thus assume to hold.We therefore target a small transverse hopping strength t ⊥ with respect to ∆E s and ∆E p .Introducing the projector onto the lowest-energy manifold P0 = i |E i,0 〉〈E i,0 |, the secondorder perturbation theory for Hamiltonian Eq. ( 1) yields: = Ĥpair + Ĥexc .(10) Within Eq. ( 10), we identify two contributions, namely a pairing term Ĥpair , which denotes the hopping of electron-electron pairs of opposite spin between neighboring chains and an exchange term Ĥexc , denoting the exchange of particles of the same spin between neighboring chains.
In the following we use MF theory to eliminate the non-1D nature of Ĥ2 ⊥ .Here, we make use of the relation and assume c ( †) j 〉 to be small.We, moreover, assume which means that all the chains are exact copies of each other.We end up with an effectively 1D expression for a Hamiltonian describing a higher-dimensional model, namely with and thus identify α n,m with the MF-approximated pairing part of Eq. ( 10) and β n,r,σ with its exchange part.Here, we introduced the coordination number z c , which denotes the number of neighboring chains.In our case z c = 4, as the chains are assembled into a 2D square grid in the ŷ − ẑ-plane.The parameters α n,m and β n,r,σ are the so-called MF parameters, meaning they need to be calculated self-consistently for all times.The work in [25] explains this for the ground state and for the finite-temperature equilibrium of the 3D system.Since the present work aims to test and benchmark the real-time dynamical version of MPS+MF itself, in the following we are working with the simplest possible version of the Hamiltonian Eq. ( 13).We neglect the exchange term β n,r,σ and allow only for site-independent onsite pairing, meaning α n,m ≡ α n,n ≡ α.This leads to with In this last expression we are adapting the evaluation of the order parameter α to the open boundary conditions.Obtaining α from an average across the entire system removes the spatial variation that is solely due to these open boundaries.
Regarding the reliability of the partial MF decoupling in the two perpendicular weak-tunneling directions, we expect that properties like the order parameter will inevitably be overestimated, as in any MF theory.For the equilibrium case, several of the authors demonstrated that the MPS+MF approach produces the correct physics compared against QMC, in regimes in which the latter approach is quasi-exact, in a negative-U Hubbard model on a 2D square lattice with anisotropic tunneling [25].That work also shows that the error in T c for the SC within the MPS+MF framework is a quasi-constant one in t ⊥ over a significant range.Moreover, at zero temperature, the overestimation of the SC order parameter becomes systematically better as t ⊥ decreases.We also point out that the degree of overestimation decreases strongly as the dimensionality of the system grows, as expected from the concurrent decrease of quantum and thermal fluctuations, which our MF-treatment partially neglects.In equilibrium, we know this due to work by some of the present authors [26], applying the MPS+MF framework to 3D lattice-bosons, which yields much more modest overestimations of key quantities such as T c for superfluidity compared to the 2D case. 2 Naturally, the good performance of the MPS+MF framework in equilibrium systems and especially in 3D by itself does not guarantee comparable performance to the non-equilibrium systems studied here.The question is thus how the framework could be validated independently.Unfortunately, while the performance for equilibrium systems could be checked in-depth for bosons, and in limiting cases for fermions, using QMC-techniques, there are no other classical computational methods that can currently reach the system sizes and time scales of the framework shown in this work.But this would be essential for any meaningful comparison to our approach, as is clear from the data shown in Sec. 4. There, it is plain that dynamically induced SC from a non-SC starting state will only build up over time scales that are well over an order of magnitude larger than t −1 , and that for 3D systems could only set in at O( 102 ) sites, given the SC healing length will be on the scale of several sites at least in x-direction, and well larger than that in the two other directions.If the linear dimension of the cluster dropped below these scales, the small size could easily preclude any SC.As discussed in Sec. 1, the existing alternatives could not reach such time and length scales due either to the growth of entanglement entropy (for techniques based solely on MPS), or the fermionic sign problem (for QMC-based ones), even for special or limiting cases.Fortunately, there is a powerful alternative for checking the performance of the scheme free from these constraints, discussed in detail in Sec. 5, which applies to that subclass of models limited to purely on-site interactions.For these, existing experiments on ultracold atomic lattice gases already offer all the features that would be required to independently verify the predictions of the dynamical MPS+MF framework.For the scope of the present work, that leaves internal consistency checks, which are done as part of the case study in Sec. 4.

MPS+MF-Algorithm for self-consistent time evolution
The expectation values needed to compute the MF parameter α in Eq. ( 17) are calculated using a self-consistent scheme for both the time evolution and for the ground-state search of our model system.In this section a schematic description of the time-evolution routine is presented, which is one of our main results.The algorithm is based on the work of H. Strand et al. published in [34], where a non-equilibrium version of real-time DMFT for bosons is introduced.Our work incorporates this real-time scheme into a MPS framework and adapts it to 3D lattices of correlated fermions built from weakly coupled 1D systems.All results obtained in the following were generated with Ian McCulloch's matrix product toolkit [35] using its time-evolving block decimation (TEBD) implementation.Note that the described algorithm is not limited to TEBD, but also other MPS based time-evolution methods [36] can be used instead.The initial ground states from which the time evolution proceeds were generated from a self-consistent scheme introduced by Bollmark et al. in [26], which is also briefly described in App. A.
At the beginning of each time step, we start with a state |ψ(t 1 )〉 at time t 1 , which we already have obtained before (either as a previous step or as initial state).From this state, we measure the value of the MF parameter α(t 1 ).Now, we guess which value α might take after one discrete time step dt.In this work, at the start of the self-consistency iterations for each time step, we just assume that the α value does not change at all.In any case, the guess for α at t 2 = t 1 + dt, is labeled α guess (t 2 ).Then, we evolve the system from t 1 to t 2 using the mean of α(t 1 ) and α guess .From the resulting tentative |ψ(t 2 )〉 we can once again measure the MF parameter α new (t 2 ).Next we calculate the distance between the measured and the guessed value and compare it to a chosen precision ϵ, For all data shown in this work, we have set ϵ = 10 −12 .If Eq. ( 18) is fulfilled, we keep the Our scheme achieves this at each discrete time step the algorithm advances, by trying to evolve with an equally weighted average of the current value of α(t) and an updated α-value, which, in the first attempt is just a heuristic guess.If the measured α-value of the new state thus evolved does not match the updated α-value (very likely in the early loops) up to some pre-defined precision ϵ, the wavefunction is discarded, and a new attempt is made, this time with the just-measured α-value as the new guess for the updated α.
state |ψ(t 2 )〉 and proceed with the next time step.Otherwise, we discard |ψ(t 2 )〉 and repeat the time step using the mean of α(t 1 ) and α new (t 2 ).The loop is repeated until Eq. ( 18) is fulfilled.A schematic of the algorithm is depicted in Fig. 2.

Transient superconductivity after a fast ramp of the nearestneighbor interaction
In this section, we present our results using the self-consistent MPS+MF scheme and find that in the extended Hubbard model Eq. ( 1) the BCS order parameter for SC grows in time and begins to oscillate around a finite value on the treated time scales.This indicates the formation of transient SC, which is the second main result of this paper.In the following, all parameters are measured in units of the hopping parameter t ≡ 1.
More specifically, we follow Paeckel et al. [19] and tune the system's parameters from an insulating charge-density wave (CDW) phase into a SC phase.However, we find that the sudden quench performed in [19] is numerically more challenging3 within the self-consistent scheme (see App. B), so we instead perform a fast ramp.
In order to check the equilibrium phases of the 3D model we use the self-consistent MPS+MF approach to compute the ground states using the routine introduced by Bollmark et al. [26] for different parameters and measure the expectation value of the MF parameter α.
We find that for t ⊥ = 0.2, U = −4 and V = 0.25 the system possesses the main properties of a CDW phase relevant for us, i.e., we find alternating occupation of the lattice sites by the electrons and a vanishing value of α.For U = −4 and V = −0.25 instead, the system is SC, as here α ∼ 10 −1 becomes finite and density oscillations less pronounced.These are the same parameters treated by Paeckel et al. in [19] for the purely 1D system.Hence, we perform a fast ramp by tuning the values of the nearest-neighbor interaction from V = 0.25 to V = −0.25 as further detailed below.
Since the effective Hamiltonian Eq. ( 16) depends on the MF parameter α(t) the question of how to choose α ini := α(t = 0) arises.For the CDW system α = 0 and it is hence difficult for it to grow with the method outlined in Fig. 2. Because of this, unless otherwise noted, our default value for this work is α ini = 10 −4 /dt, where dt is the size of the discretized time step of the simulation.Such a small yet finite value is justified by the fact that any system will either have a microscopic fraction of pairs in the center-of-mass zero-momentum state to begin with, or such a fraction is generated during the ramp or quench.Scaling α ini inversely in dt ensures that simulations with different dt agree over long times, see Fig. 3.
The MF term of the Hamiltonian causes the effective model to be no longer particle-number conserving, hence, we need to adjust the value of the chemical potential µ corresponding to the system size and to the onsite repulsion U in order to fix the average density of the total system.From the ground-state calculations we find the values of µ that are listed in table 1.We keep the values of µ, determined in this manner, fixed throughout the whole time evolution in order to keep our algorithm simple and stable.However, we still need to keep track of the overall density of our system during the time evolution to check if this assumption of a time-independent chemical potential is justified.Indeed, for our simulations, the value of the density is preserved to a good accuracy over the time scales treated by us (see Figs. 3 and 4).In general, however, it might be necessary to also include a variation of µ into the self-consistency scheme.

Time evolution of the BCS order parameter and of the total energy
In the following, we investigate the time evolution of the BCS order parameter α(t) (see Eq. ( 17)) and of the total energy E(t) of the system.The latter cannot be expected to remain  1.We compare the ramp scenario (solid violet and dashed green) with an evolution during which we keep the nearest neighbor interaction at V = 0.25 constant (dotted blue).
Even though this state is given the same initial α-value as a seed as in the ramp scenario, α(t) rapidly decays and remaining at near-constant and near-zero values.This shows that no SC state develops without a quench to parameters corresponding to SC in equilibrium.For this calculation we chose a time step of dt = 0.01.constant as the MF term changes the Hamiltonian Eq. ( 16) during evolution.In addition, we monitor the total density of the system, which should stay at a value of ρ = 1 (half filling) during the whole time evolution.
Since we find fast ramps to have lower errors over the simulated time windows than instantaneous quenches, we linearly decrease the value of the nearest-neighbor interaction V from V = 0.25 to V = −0.25 within a time window of ∆t ramp = 3.0.A more detailed discussion of the effect of the size of the time window ∆t ramp can be found in App.B. In Fig. 3 we see the results for a 30-site system for an evolution up to time t end = 50.Since α(t) is complex-valued we show the evolution of the magnitude |α(t)| and of the phase ϕ(t) of the order parameter in Figs. 3 to 5. We find that |α(t)| grows up to time t ∼ 45 to a value of approximately |α| ≈ 0.06, which is clearly non vanishing and hence indicates the formation of a non-equilibrium SC state.In contrast, if we consider a time evolution without a quench or ramp, i.e., V = 0.25 during the whole evolution, the value of α stays unchanged at an order of magnitude of 10 −5 throughout the whole time evolution as can be seen by the dotted blue Both the linear and the quadratic fit suggest a finite and comparable value of t SC in the limit L → ∞.Evolution of the total energy of the system and the total density are shown in (c) and (d).All data were obtained with a bond dimension of χ = 500, an initial guess of the MF parameter of α ini = 10 −4 /dt, a ramp time window ∆t ramp = 3, and the chemical potential was taken from table 1. lines in Fig. 3.The phase ϕ(t) decreases as long as V is decreasing, then oscillates around a value of approximately ϕ(α)/π ≈ −0.8 and seems to increase again slightly when |α| has reached its maximum.We interpret this behavior as an expression of a Josephson effect in-between 1D chains to the extent it can be captured by a single 1D system with time-evolving MF amplitudes.As a kernel of SC order manifests itself in the different chains of the 2D array the macroscopic phases of SC states, within each chain, will be initially uncorrelated, then start aligning via the Josephson effect.With density fluctuating within each individual chain the Josephson effect will keep the phase fluctuating while the system finds a new equilibrium after the rapid ramp, as Fig. 3b shows.
In Fig. 3d we show the evolution of the total energy per site E(t)/L and in Fig. 3e the deviation of the total density ρ(t) from the desired value ρ target = 1.We find that this deviation is of the order of 3 • 10 −5 or smaller for all the times treated, indicating that keeping the chemical potential µ fixed leads only to small errors.The total energy per site E/L behaves as expected during the ramp and decreases almost linearly for the duration of the ramp.Afterwards, we first observe a nearly constant behavior, then a strong decrease until a minimum at time t ≈ 45, shown in the inset of Fig. 3a.We read the behavior of E(t)/L, especially at long times, as the system starting to further lower its energy through condensing Cooper pairs, as the drop in E(t)/L coincides markedly with the onset of a finite value of α(t).We also study the effect of system size, to make certain the dynamical onset of SC would survive in the thermodynamic limit.In Fig. 4 we compare the results for different chain lengths L. From these, we extract the instant t SC , at which |α(t)| reaches its first maximum.The data of the 12-site system shows the onset of oscillation for |α(t)| around a finite value, indicating a dynamically induced SC phase (longer-time simulations for L = 12 further confirm this, as shown in Figs. 5 to 7 for times up to t max = 100).The inset of Fig. 5a displays an extrapolation in inverse chain Figure 5: Evolution of the MF parameter α split up into its magnitude (a) and its phase (b) for different initial guesses α ini in a 12-site system.We see that the reduction of α ini induces merely a shift in the data, at least up to time t SC at which the first maximum of |α| occurs.The inset in (a) shows t SC vs. α ini and a linear fit on a semilogarithmic scale.This shows that t SC grows merely logarithmically with α ini .
The data shown was obtained with χ = 500 and dt = 0.01, a ramp time window ∆t ramp = 3, and µ was taken from table 1.We stress that the evolution of |α| up to the first maximum at t SC is nearly identical for different α ini up to a shift, as highlighted by the logarithmic scale.Beyond t SC , curves for different α ini oscillate around closely similar averages, and any deviations of the α(t)-curves from each other are due to the inevitable finite-size effects, which will naturally depend on the values of α ini .
length 1/L of t SC .In order to see whether t SC diverges we performed a quadratic and a linear fit, both indicating a finite value in the limit L → ∞.Since for the larger system sizes |α(t)| starts to oscillate at around the maximal time reached by us, it is difficult to obtain a finite-size extrapolation of the value of the SC order parameter.In order to do so, one needs to extend the simulations for the larger systems to substantially longer times, which is beyond the scope of this paper.

Accuracy of results and sensitivity to simulation parameters
The results so far were all obtained using the same parameters for the self-consistency cycle.
In the following we study how sensitive the results are on parameters like the initial guess of the MF parameter α ini (see Sec. 3), the bond dimension of the MPS calculations, or the discrete time step dt.To study these effects, we focus on the 12-site system in order to reach the longest time scales.Figure 5 shows the evolution of the magnitude and phase of α(t) for different initial values α ini .Decreasing the value of α ini induces a shift of t SC to later times.In order to further analyze this, we plot the value of t SC against the value of α ini in the inset of Fig. 5a.Speaking to the soundness of our MF approximation, we find that t SC increases only very weakly with α ini , i.e., logarithmically.While this indicates a diverging time for the onset of SC order in the limit α ini → 0, this is merely consistent with α ini = 0 being an unstable fix point of the dynamic MF algorithm in the regime we ramp into.But any finite value, even a microscopic one, will yield dynamically induced SC order in finite time when ramping into the SC parameter regime.Furthermore, as argued at the outset of Sec. 4, on general physical grounds there will always be some electron pairs whose center-of-mass momentum is zero.
In Fig. 3 we compare two different discretized time steps, dt = 0.005 and dt = 0.01, respectively.The results are nearly identical, only a small deviation of the total density, which agrees up to ∼ 10 −5 , can be seen in Fig. 3e.Regarding the discarded weight of our simulations, we find that even for the smallest bond dimension these values stay below 10 −6 within the time domains we consider.Nevertheless, we examine the dependence of our results on the MPS bond dimension χ.For this purpose we compute the deviation of the value of an observable O for two different values of χ, At any fixed value of α ini and dt we find this to be the most reliable estimator for the accuracy of our combined MPS+MF approach (assuming the latter parameter is chosen to be sufficiently small) and focus in the following on this quantity.In Figs. 6 and 7 we present results for the observables |α(t)| and E(t) obtained with two different bond dimensions χ 1 = 500 and χ 2 = 1000 for the 12-site system, and for χ = 250 and χ = 500 for the 30-site system, respectively, and also the difference of the respective results.For the larger system it was necessary to substantially reduce the values of χ, since otherwise the numerical expenses would exceed the available resources.We find that the deviation of the results is ∼ 10 −6 for the values of |α(t)| and ∼ 10 −4 for the total energy E(t), in the case of the 12-site system.For both observables, this is small compared to the order of magnitude of the observables themselves, so that we conclude these values of χ suffice to provide quantitatively accurately results, within the dynamical MPS+MF framework.
For the 30-site system, however, the deviation is ∼ 10 −3 for |α(t)| and ∼ 10 −2 for E(t).This is rather large in comparison to the order of magnitude of the observables themselves.The data obtained from these calculations is hence only trustworthy in regards to the qualitative physics, but for the larger chain lengths one needs a larger bond dimension to obtain a better quantitative convergence of the results.

Conclusion
This work presents a self-consistent real-time MPS+MF approach for investigating the time evolution of a 3D extended Hubbard model after a fast ramp.By combining perturbation theory with a MF ansatz, we construct an effective 1D Hamiltonian Eq. ( 13) capable of capturing the dynamical build-up of SC correlations for this 3D model system, when quenching or rapidly ramping into a Hamiltonian parameter regime corresponding to SC order in equilibrium.This approach is generic to any 3D system composed out of gapped 1D systems of fermions, as long as coupling between 1D systems is weak for single-fermion tunneling in-between 1D systems to be suppressed.For concrete demonstration of the performance of this approach, we chose systems of 1D extended Hubbard chains, arranged in parallel in a 2D square array, forming a 3D system with weak interchain tunneling t ⊥ , negative onsite repulsion U, and nearest-neighbor interaction V along each chain.We benchmark the self-consistent algorithm introduced on the simplest possible version Eq. ( 16) of the resulting effective MF Hamiltonian, only taking onsite pairing into account and neglecting the particle-hole terms Eq. (15).We test our approach on systems where each chain is up to L = 30 sites long.Using this algorithm we compute the time evolution of the BCS order parameter for SC order α(t), as a direct indicator of dynamically induced SC.The results show that SC order sets in after a fast ramp from V = 0.25 to V = −0.25,where the initial V -value realizes an insulating CDW state, and the final value would correspond to SC order at equilibrium.These results are broadly comparable to previous 1D results [19] and represent a best-case scenario, in which double occupancies already present in the CDW help to form the non-equilibrium SC state after the ramp.
Performing infinite-size extrapolations and studying the effect of the microscopic initial kernel of SC order α ini shows that dynamically induced SC is not merely a trivial size effect, but actually present in the thermodynamic limit, and even the smallest yet finite magnitude for α ini will result in establishing order within a finite window of time.At the same time, we find that resource requirements increase substantially with chain length L, but several tens of sites and time frames between one and two orders of magnitude in units of inverse fermion tunneling t −1 are accessible already with the modest resources employed for the present proofof-principle work.
The present work presents multiple avenues for interesting and potentially valuable follow-up work.One of these would be to move towards a regime that is physically more realistic as far as solid state systems are concerned, in which the pair-binding energies ∆E p would be significantly smaller than in the present work.This would entail either lowering U, or working directly with a 1D model offering repulsively mediated pairing, such as a doped two-leg Hubbard ladder [37,38].This would require retaining more particle-particle terms Eq. ( 14) than we have done for the present proof-of-principle, as well as incorporating the particle-hole terms Eq. ( 15) into the self-consistent time-evolution step, see Fig. 2.This would be straightforward, as a generic ansatz for the first iteration of these terms is practically imposed by the physics of these 1D systems.As detailed in, e.g., [39], both classes of terms decay with an exponential envelope function characterized by the spin-correlation length, which in turn is easy to obtain from static correlators via density-matrix renormalization group (DMRG) [40,41] calculations for the isolated systems.Finally, we observe that while in the present work the initial state was a pure zero-temperature one, extensions of our new approach to finite temperature would be entirely straightforward.Being ultimately MPS-based, it can tap into practically all previous refinements of MPS-numerics, such as the simulation of non-equilibrium dynamics of finite-temperature states that minimizes the growth of entanglement with time [42].
With the above refinements, our approach opens new and potentially exciting avenues for advancing the understanding of dynamically induced SC in solids jointly with the domain of analog quantum simulation.Existing experiments with ultracold atomic lattice gases of fermions in 3D already offer the capability of realizing arrays of weakly coupled negative-U Hubbard chains initially in the analogue of the PG regime.Results from ref. [25] (c.f.Fig. 11a there) illustrate that with today's experimental control over the parameters of these latticegas experiments such a regime can practically always be found.That is, inside each chain upand down-pseudospin atoms would be paired, with an associated pairing energy ∆E p , from the effective attractive on-site interaction, but the coupling t ⊥ between these chains would be deliberately chosen to be too weak for the 3D array to enter the SC regime at the working temperature of the experiment.In this manner, these analog quantum simulators would mimick the basic presumed initial conditions of the layered high-T c compounds in pump-probe experiments on dynamically induced SC.However, there would be two key differences to the solid state: the active sub-units would be 1D chains instead of 2D planes, and as a result of that these analog quantum simulations would be amenable to be modeled by quantitative many body numerics, i.e., the approach developed in the present work and incorporating some of the refinements laid out in the preceding paragraphs, in a way that the solid state experiments are not.Then t ⊥ would be quenched or rapidly increased to a value that at equilibrium would correspond to the SC regime, analogous to other types of experimentally performed rapid parameter changes (c.f.[43,44] and references therein).Any SC emerging dynamically from this change in t ⊥ could be directly tracked in real time and unambiguously via the coherence peak of paired atoms, which would constitute direct proof of the onset of SC off-diagonal long-range order in these 3D ultracold atomic lattice gases [45].The dynamical MPS+MF framework developed here could then be used to directly model these experiments, ascertaining whether any predicted onset of dynamically induced SC qualitatively and quantitatively matches the observations.It would further allow to check whether the natural choice for initializing α ini , the residual pairing correlations measured experimentally for the initial non-SC state, would yield the correct time scale for SC to emerge dynamically.This modeling would thus serve the dual purpose of validating the numerical framework itself, which would be necessary as there is currently no alternative classical numerical approach that could approach the system sizes, time scales and temperatures accessible to our dynamical MPS+MF approach.In this way, comparison between experiment and our theory framework would also advance the field of analog quantum simulation of demanding non-equilibrium many-body problems.
Once these comparisons work sufficiently well, a natural extension would be to advance to experiments in which 2D planes of negative-U Hubbard models are initially weakly coupled with a t ⊥ -term, s.t.one would again start in a PG-regime, but this time even closer to the basic physics at the root of the experiments of the layered high-T c materials.Anchored by the previous experiments on the quasi 1D arrays and the quantitative comparison made possible by the theory framework developed here, the outcomes of such a second stage of experiments could provide a strong input to the ongoing discussion as to how to interpret experiments on dynamically induced SC in layered high-T c materials.showing the difference between the two bond dimensions.We gain an accuracy of the order of 10 2 via the fast ramp, as opposed to the instant quench.dimensions χ, as it was done in Sec.4.2 as a check of accuracy.We find that difference is two orders of magnitude smaller for the ramp compared to the case of the instantaneous quench.This is why we chose to use ramps for all our calculations presented in this paper.

Figure 3 :
Figure3: Evolution of the considered parameters in time during and after a ramp on a 30-site system.The plots at the bottom ((c) and (f)) show the nearest-neighbor interaction, which decreases from V = 0.25 to V = −0.25 during a time window of ∆t ramp = 3.0.Graphs (a) and (b) show the evolution of the MF parameter α split up into magnitude and phase.Graphs (d) and (e) show the evolution of the total energy per site of the system and of the total density.The inset in (d) shows the evolution of the energy per site after V was decreased.The dip towards the end of the evolution coincides with the rapid increase in |α(t)|, and signals the system overall lowering it's energy by entering a SC state.The legend is valid for all plots.All the data shown here were obtained with a bond dimension of χ = 250, an initial guess of the MF parameter of α ini = 10 −4 /dt, and the chemical potential was taken from table 1.We compare the ramp scenario (solid violet and dashed green) with an evolution during which we keep the nearest neighbor interaction at V = 0.25 constant (dotted blue).Even though this state is given the same initial α-value as a seed as in the ramp scenario, α(t) rapidly decays and remaining at near-constant and near-zero values.This shows that no SC state develops without a quench to parameters corresponding to SC in equilibrium.For this calculation we chose a time step of dt = 0.01.

Figure 4 :
Figure 4: Evolution of SC order-parameter α, energy per site E/L, and density ρ in time during and after a ramp for sizes L = 12, 16, 20, 30.Magnitude and phase of α are shown in (a) and (b), respectively.The inset in (a) shows the time at which the first local maximum in |α| occurs plotted against the inverse of the system size 1/L.Both the linear and the quadratic fit suggest a finite and comparable value of t SC in the limit L → ∞.Evolution of the total energy of the system and the total density are shown in (c) and (d).All data were obtained with a bond dimension of χ = 500, an initial guess of the MF parameter of α ini = 10 −4 /dt, a ramp time window ∆t ramp = 3, and the chemical potential was taken from table 1.

Figure 6 :
Figure 6: Evolution of (a) the magnitude of the MF parameter |α| and (c) the total energy E for two different bond dimensions χ 1 = 500 and χ 2 = 1000 in a 12site system.(b) and (d) show the difference δ O,χ 1 ,χ 2 between the observables we measure for these two different bond dimensions.All calculations were done with α ini = 10 −4 /dt, a ramp time window ∆t ramp = 3, and dt = 0.01.

Figure 9 :
Figure 9: Comparing an instantaneous quench (left hand side) with a fast ramp in finite time (right hand side).Top row showing the time evolution of |α| for two different bond dimensions χ = 500 and χ = 1000 in a 12-site system, bottom rowshowing the difference between the two bond dimensions.We gain an accuracy of the order of 10 2 via the fast ramp, as opposed to the instant quench.