Solving 2D and 3D Lattice Models of Correlated Fermions — Combining Matrix Product States with Mean-Field Theory

Correlated electron states are at the root of many important phenomena including unconventional superconductivity (USC), where electron pairing arises from repulsive interactions. Computing the properties of correlated electrons, such as the critical temperature T c for the onset of USC, efficiently and reliably from the microscopic physics with quantitative methods remains a major challenge for almost all models and materials. In this theoretical work, we combine matrix product states (MPS) with static mean field (MF) to provide a solution to this challenge for quasi-one-dimensional (Q1D) systems: two-and three-dimensional materials comprised of weakly coupled correlated 1D fermions. This MPS þ MF framework for the ground state and thermal equilibrium properties of Q1D fermions is developed and validated for attractive Hubbard systems first, and further enhanced via analytical field theory. We then deploy it to compute T c for superconductivity in 3D arrays of weakly coupled, doped, and repulsive Hubbard ladders. The MPS þ MF framework thus enables the quantitative study of USC and high-T c superconductivity — and potentially many more correlated phases — in fermionic Q1D systems based directly on their microscopic parameters, in ways inaccessible to previous methods. This approach further allows one to treat competing macroscopic orders, such as superconducting and insulating ones, on an equal footing. Benchmarks of the framework using auxiliary-field quantum Monte Carlo techniques show that the overestimation of, e.g., T c due to its mean-field component, is near constant in microscopic parameters. These features of the MPS þ MF approach to correlated fermions openupthe possibility ofdesigningdeliberately optimized Q1D superconductors, from experiments in ultracold gases to synthesizing new materials.


I. INTRODUCTION
Obtaining quantitative and reliable solutions to strongly correlated fermionic models from first principles is one of the greatest challenges to the theory of quantum matter, arising in many different areas, from solid state physics to ultra cold atomic gases.Its prominence derives, to a large degree, from the fundamental and technological importance of many of the collective quantum states emerging from electronic correlations.
Unconventional superconductivity (USC) epitomizes this challenge: Repulsive interactions lead to electrons forming correlated pairs which attain macroscopic phase coherence and enter a superconducting state at the critical temperature T c .Crucially, the high-T c superconducting models and materials belong to this group [1][2][3].These can be divided into two classes: the quasitwo-dimensional (Q2D) models and materials, and the quasi-one-dimensional (Q1D) ones.In both cases, the relevant three-dimensional (3D) system is comprised of weakly coupled lower-dimensional sub-units, 2D and 1D ones for Q2D and Q1D respectively.
The present work introduces a quantitative, efficient and reliable theoretical framework for Q1D systems, capable of calculating properties of unconventional and high-T c superconductivity as well as other, competing correlated phases of fermions from microscopic parameters.It thus allows to deliberately engineer high-T c su-perconductivity from the bottom up in such systems.This is achieved by leveraging two unique advantages with respect to USC that the 1D sub-units within a Q1D system can have compared to the known 2D sub-units in the Q2D ones: (1) The microscopic mechanism of pair formation from repulsive interactions can typically be worked out [4].(2) The pairing energy resulting from it can be computed accurately from microscopic parameters [5,6].These pairing energies can moreover be very high.For doped two-leg Hubbard ladders, pairing energies up to about 15% of electron tunneling have already been demonstrated [5], as the 1D nature of these ladders naturally enhances the repulsively mediated effective attraction.As fluctuations preclude superconductivity in isolated 1D sub-units, the coupling of these into a threedimensional Q1D array is essential to enable macroscopic ordering of pair-phases.
Analysis of the resulting arrays has so far been limited to Tomonaga-Luttinger liquid (TLL) field theory [4] combined with static mean-field (MF) techniques [5,7].The TLL approach has many uses, such as explaining repulsively-mediated pairing, or mapping the lowenergy theory of Hubbard-ladders to that of the simpler negative-U Hubbard chains.However, it allows neither to compute T c of a Q1D array quantitatively, nor to determine from microscopic parameters whether it realizes USC or another correlated phase.Furthermore, this approach also disregards exchange processes between 1D sub-units.Thus, it cannot inform the search for USC and high-T c systems in the Q1D-space, be it for synthesizing candidate materials or guiding quantum simulations towards analogue states of USC within current experimental capabilities.Both are intensely sought goals that remain highly challenging for the theory of Q2D superconductors with repulsively-mediated pairing [1,2,[8][9][10][11][12][13][14].The latter issue is acute for finding USC-like states in highly controlled ultracold gas experiments with fermions, i.e., analog quantum simulators, for the 2D Hubbard models [15][16][17][18][19][20].As current theory cannot efficiently determine where and whether these realize a USC state within the achievable entropies, the impasse on classical computational theory for 2D fermions still impedes progress on the quantum simulation front.
Shifting focus to Q1D systems instead, this work establishes a numerical theory framework capable of meeting these aims.We further enhance the framework via novel TLL+MF theory that allows to compute T c from easierto-obtain zero-temperature numerics.Briefly summarized, the framework exploits that in the weak coupling regime between the 1D sub-units we can apply perturbation theory in the ratio of the single-fermion tunneling in-between 1D sub-units to the pairing energy.We show how this ratio naturally controls the possible superconducting T c .In this perturbative regime, which is different from USC and high-T c superconductivity in Q2D system such as the cuprates, the full 3D array can be decoupled via static MF theory.This maps the problem to that of a single 1D sub-unit with multiple MF-amplitudes, which is then solved self-consistently at the microscopic level.As summarized in Fig. 1, these represent the various possible injections/ejections of fermion pairs into/out of the 1D sub-unit stemming from the 3D array, as well as the previously neglected exchange processes.Singlefermion tunneling in-between 1D sub-units is naturally suppressed in this perturbative regime.
This framework is not limited to superconductivity.It can be applied to any fermionic Q1D setup in which the 1D sub-units have some type of gap sufficient for inter-system tunneling to be relevant first at second order.Such a regime is realized for instance in the various insulating phases of the Bechgaard and Fabre salts [21,22].Alternatively, such regimes may be found, with minor modifications, in the recently proposed and partially realized mixed-dimensional large-J systems in ultracold fermionic lattice gases [23,24].At the same time, the quantitative study of USC regimes of the same compounds, as well as of other Q1D materials such as Chromium pnictide [25,26] and the telephonenumber compounds [27,28], would require including single-electron tunneling in-between 1D sub-units.This would be beyond the MPS+MF framework described here, necessitating an approach akin to chain-dynamical mean field theory [29,30].
The present work hinges on an efficient, reliable numerical method at the microscopic level for correlatedfermion 1D sub-units including MF-amplitudes.Algo-rithms using matrix product states (MPS) [31], such as the density matrix renormalization group (DMRG) [32,33] are uniquely suited here.
Frameworks using MPS+MF for Q1D spin [34,35] and bosonic systems [36] have been successfully employed, e.g., on experiments in spin-ladder materials.They compare well to Quantum Monte Carlo (QMC) algorithms, despite the MF approximation in the weak-coupling directions, at greatly increased efficiency.
Applying the MPS+MF approach to fermions is more demanding than for spins or bosons, as the MFamplitudes are more numerous than for those systems, which typically have one such amplitude.Moreover, these amplitudes range over multiple sites, raising the bond dimension of the matrix product operator (MPO) representing the total Hamiltonian.When the 1D subunits have internal structure, as the two-leg Hubbard ladder does, even higher performance is required.Modern MPS implementations can treat even such complex models with many long-range coupling terms, despite the potentially large bipartite entanglement, which sets computational complexity in MPS-methods.Work on DMRG for ground states of 2D Hubbard models demonstrate this fact [9][10][11][12][13]37].A recent DMRG-implementation for distributed-memory architectures, parallel DMRG (pDMRG), can treat very large clusters by MPS standards, of e.g., the 2D U-V Hubbard models of the Bechgaard and Fabre salts [14].Yet, calculations for correlated fermions in Q1D-systems cannot be handled in this way, the bipartite entanglement in 3D would be far too large.Auxiliary-field QMC (AFQMC) methods [38][39][40] are applied in this work for the cases with attractive interactions, but they are typically computationally intensive, and in general the fermionic sign problem prevents exact solutions.
The MF-approximation within the MPS+MF scheme for fermions thus allows the present work to study the ground state and thermal properties of much larger correlated Q1D systems, including those with repulsive interactions, than either brute-force MPS-methods or AFQMC.Being primarily MPS-based, the framework can also be extended to non-equilibrium real-time problems, such as the study of dynamically induced superconductivity in a Q1D system [41].
The present work is structured as follows: Sec.II describes a basic Q1D model of superconductivity, an array of negative-U Hubbard-chains, and how to obtain a 1D problem with self-consistent MF-amplitudes from it.Sec.III defines the MPS+MF algorithm used to solve this effective 1D model, and motivates various optimizations.The MF-order treated here is a superconducting one, but we briefly discuss how the framework can be extended to address multiple competing orders simultaneously.Then, Sec.IV features an analytical treatment of the same model.In Sec.V, results of the MPS+MF framework are shown.These are compared to analytical methods, and we show how these can be used to obtain T c from ground-state calculations.Furthermore,

II. MODEL
The two related models that are studied in this paper are illustrated in Fig. 1.Both are Q1D arrays in 2D and 3D, comprised of many identical 1D sub-units in parallel (e.g., chains and ladders).They are described by a Hamiltonian H 0 , and connected by a perturbing Hamiltonian H ⊥ with the general form The first model we consider is shown in Fig. 1a and permits extensive testing and validation and is treated in sections Secs.II to V. It features a 3D Hamiltonian with anisotropic tunneling and interactions where H 0 and H ⊥ take the following forms: and Here, {R i } indicates a 2D subspace spanned by ŷ, ẑ.The operator of c † n,Ri,σ creates a fermion on site n in a chain at R i with spin σ.The parameters t, µ and U denote hopping, chemical potential and on-site interaction respectively.The only higher-dimensional effect in this model is given by H ⊥ .
The second model we consider is a generalization of the first.It is made up of weakly coupled doped Hubbard ladder systems with repulsive interactions shown in Fig. 1b.The definition of H 0 and H ⊥ differ somewhat in appearance from the first and is specified further in Sec.VI and App.E. The methods developed to treat the first model in Secs.II to V will be used to treat this latter, more complicated case in Sec.VI.
For both models, there are two fundamental energy scales: the spin gap and the pairing energy.The former measures the cost of energy to flip a spin in a single 1D system at any position R i .It is defined by where E(N, S) is the ground state energy of H 0 (R i ) (i.e., a single 1D sub-unit) at charge and spin quantum numbers N and S respectively.Conversely, the pairing energy is the cost to move one spin from a S = 0 chain to a neighbouring 1D sub-unit, also at S = 0, creating two chains at S = 1/2 As will be seen in the following section these energy scales and in particular ∆E p , for which generally ∆E s ≤ ∆E p (see for instance [5]), will govern the strength of higher dimensional effects.

A. Perturbation theory
The models presented in the previous section are typically a challenge to solve numerically and analytically.In particular, the doped Hubbard ladder with repulsive interactions would be impossible to treat in a 3D array of any meaningful size.However, restricting ourselves to the case of the gaps in Eqs. ( 4) and ( 5) being large in comparison to the strength of the perturbing Hamiltonian H ⊥ we can construct a perturbation theory in the ratio t ⊥ /∆E p , that is, in order to solve Hamiltonians of the form Eq. ( 1) we specialize to the case where t ⊥ ∆E p .Specifically, when U < 0 and t ⊥ = 0 this model describes a set of disconnected 1D Hamiltonians H 0 (R i ).Each such Hamiltonian exhibits a spectrum of the form [4] where α, β indicate changes to the state that induce a large energy shift and i, j small shifts.The labeling is used to distinguish energy states i, j which lie in the same manifold and α, β which denote what manifold the state is in.In the case of our model this is related to spin singlet pair-forming by the parameter U < 0. The subsequent break-up of such pairs by spin-flip costs ∆E s in the isolated 1D system, giving rise to (6) which is illustrated in Fig. 2. Conversely, when a pair is broken by a process like Eq. ( 3) two well-separated and unpaired spins are left which changes this cost to ∆E p .When t ⊥ is much smaller than the pairing energy it is possible to produce an effective Hamiltonian with perturbation theory acting in manifold α with matrix elements [42] The action of H ⊥ on |j, α changes the manifold of the state.We assume that the energy differences between states inside manifolds are much smaller than that of states in different manifolds: where we let α = 0 indicate the lowest manifold and α = 2 the one reached by H ⊥ .In other words, we approximate the energy difference between states of different manifolds to be independent of the exact state in that manifold.This yields the effective Hamiltonian where P α = i |i, α i, α| is a projector onto manifold α.
By design, H ⊥ puts a state in manifold 0 into manifold 2. This means that Thus, the first order and odd order contribution of the effective Hamiltonian disappears and we are left with In order to understand how H 2 ⊥ acts on a state we expand: obtaining a Hamiltonian describing two-particle tunneling events.What characterizes the α = 0 manifold is pair formation of opposite spins due to attractive interaction.Some terms within H 2 ⊥ will put the state in a manifold β > α and will subsequently be projected out.Most importantly, any terms flipping spins in two chains simultaneously, such as c † n,Ri,↑ c n,Ri+ŷ,↑ c † n,Ri+ŷ,↓ c n,Ri,↓ , which would move two chains initially in their S = 0 manifolds into their S = ±1 manifolds, will be projected out due to describing a state with at least energy 2∆E s above the lowest-energy manifold.In this manner, it is clear that each chain in our model has to conserve spin and in particular pair all spins such that total spin S = 0.
The parts of H 2 ⊥ which remain after projecting to α = 0 must either move particles as pairs (simultaneous tunneling of an up-spin and down-spin particle) or exchange them between chains (an up/down-spin particle is moved out of the chain and another of the same spin is moved in).Any of these processes involve at most two chains in order to conserve spin.While H 2 ⊥ allows for processes involving chains at arbitrary distance, these would concern four separate chains, with a final state at least 2∆E p above the low-energy manifold, and thus removed by projection.
With these restrictions H 2 ⊥ becomes heavily reduced: The first set of operators, H pair , correspond to pairs of fermions jumping from one chain to another.This conserves spin but not number of particles within a chain.
The second set of operators, H exc , correspond to fermions of like spin switching chains (not necessarily at the same site within a chain).Within this Hamiltonian there is still the degree of freedom to pick how large |n − m| we include.This is something that cannot be restricted by symmetry arguments or the smallness of t ⊥ /∆E p .

B. Mean-field theory
Notably, the perturbation theory has produced an effective quartic interaction from the single-particle tunneling Hamiltonian.This can be further simplified by casting Eq. ( 13) in a form where there is no explicit mention of other chains.This allows us to solve an effective 1D model instead of the full 3D system.We use an ansatz of quasi-free states: which produces a mean-field Hamiltonian from the quartic operators in Eq. (13).We further approximate the expectation value of any operator creating/annihilating particles on two different chains to zero, i.e., pair constituents cannot live on different chains.This is motivated by choosing small enough t ⊥ /∆E p , amounting to the standard mean-field approximation.

Pairing terms
For each R i of H pair we obtain where z c = 4 is the coordination number for three dimensions.Notably the mean-field is performed on two dimensions instead of the full three.We have assumed that (16) i.e., that all chains are identical and the mean-field is real-valued.

Exchange terms
For each R i of H exc we obtain where we have done similarly as in Eq. ( 16).

C. Effective 1D Hamiltonian
To summarize: within the MPS+MF approach, the total Hamiltonian is described by its one-dimensional subsets in a mean-field sense.It is sufficient to consider the effectively 1D Hamiltonian where are the pair-MF parameters, describing pair-tunneling into/out of the 1D system mediated by H pair , and are the exchange processes mediated by H exc (c.f.Eq. ( 13)).Notably, α ik is proportional to the bound state pair-wavefunction.In practice, keeping all ranges will make the problem intractable and limits must be introduced.Fortunately, the existence of a spin gap causes α and β to decay exponentially with distance controlled by the spin correlation length [4].This allows us to include a cut-off for which amplitudes to keep.The choice of this value ultimately depends on microscopic parameters (in particular interaction strength) and will be exemplified in Sec.V.In addition, it is important to note that the zero-range terms β i0σ ∝ n iσ are absent in H MF as we assume constant density in the superconducting order and exclude potential insulating orders.We have already extended the MPS+MF framework to study the competition of insulating charge-orders with superconducting ones, based on adding β i0σ -parameters, and are applying this to the Hubbard-ladder arrays in forthcoming work.
Finally, as discussed in Sec.II A, MF-terms such as S + i S − i cannot occur in H MF by construction, as such magnetic exchange-terms are suppressed through the condition that t ⊥ ∆E s .

III. MPS+MF -NUMERICAL SOLUTIONS TO Q1D SYSTEMS USING MEAN-FIELD THEORY
MPS+MF has been developed to solve Q1D models by relying on mean-field approximations of the full Hamiltonian which converts it to an effectively 1D system.The produced Hamiltonian in Eq. ( 18) may be solved iteratively until self-consistency is reached for the mean-field amplitudes as shown schematically in Fig. 3.This framework has been developed as an approach to anisotropic systems in two or more dimensions for which the 1D correlations are most important and where single-particle tunneling and spin-flipping processes between chains can be neglected on account of t ⊥ ∆E p .The primary cost of the routine comes from the repeated solutions of effective 1D mean-field Hamiltonian such as (18) using MPS-based methods.Notably, any MPS-method may be used to iterate mean-field amplitudes, for example original DMRG, MPS-based DMRG or imaginary time evolution on purified states to obtain thermal states.In practice, the only requirement is that the mean-field amplitudes converge.Thus, the framework scales as the utilized DMRG method does with bond dimension and system size [31].In this paper we utilize MPS-based DMRG to solve for ground states [31].For thermal states we utilize both trotterized imaginary time evolution and the time-dependent variational principle (TDVP) [43].
One common feature of states obtained with DMRG ground state algorithms is that the error of local quantities scales linearly with the discarded weight [44].Given that the system under study is 1D, albeit representing a 3D system, we might expect that such linear scaling holds for MPS+MF as well.In Fig. 4 we show an example of our general finding that this holds true for quantities such as energy and order parameter.

A. Achieving self-consistent convergence
When iterating a MF-theory to self-consistency, the required number of iterations becomes crucial as each one  may have significant cost.The number of solutions required to achieve a self-consistent state is not constant over parameter space as shown in Fig. 5a.However, it only varies modestly with Hamiltonian parameters, such as t ⊥ .The notable increase in loop number with larger t ⊥ is primarily due to increased difficulty fixing the density (see Sec. IIID).Furthermore, we find that the number of required iterations peaks strongly around the phase transition from superconducting to normal phase as shown in Fig. 5b where the transition occurs around T = 0.105 [36,45,46].When the superconductor is close to its transition to a metal the MF-coupling between chains is weak but nonzero.We find this slows down convergence rate explaining the increase in the number of required loops for convergence.
In both this work and the primary uses to which the MPS+MF framework would be applied, the interest in phase transitions is central and we shall have to resolve the points which are difficult to obtain.This has prompted the development of several heuristics in this method to produce a faster convergence.

B. Extrapolation
Convergence implies that the mean-field amplitudes approach a set which no longer changes with further iteration, i.e., the change of these amplitudes with each iteration decreases.As can be seen from Fig. 6 we find such exponential behaviour after an initial fluctuation (related to the initial guess).Given the clear trend an extrapolation can be performed to attempt a prediction of the converged amplitudes.
In practice, such fits seldom give precise results but typically lead to a more converged amplitude compared to simply iterating one more loop.As such, we employ a strategy of repeated extrapolation in an attempt to speed up the self-consistent iteration as shown in Fig. 6.The result of 5 consecutive iterations is checked for sufficient exponential behaviour and is fitted to an exponential function where the fitting parameters A, B, C are determined by a least-squares fit.In this manner, C is the fitted result when approaching infinite iterations and is chosen as the amplitude to use for the next iteration.Typically, we find that this scheme reduces the number of necessary iterations to varying degree.In particular, close to transitions the slow convergence can be greatly aided by skipping ahead using extrapolations.In effect, the scheme is a method of breaking off the iterative procedure in order to generate an improved initial guess using exponential extrapolation.The amount of speedup depends on how evenly and slowly the mean-field amplitudes converge.If the convergence is quick, fitting to an exponential form is less faithful.Conversely, if conver-gence is too slow the extrapolated result makes predictions far outside any region of reliability (e.g., using 10 iterations to predict what would happen after 1000 iterations).In the best case scenarios we find the number of iterations reduced by up to a factor of 5 and in the worst case it only executes once or not at all leading to small or no speedup.

C. Initial guesses
As shown in Fig. 3 the routine has to be started with some input mean-field amplitudes and chemical potential.This indicates that the number of DMRG solutions is dependent on the quality of the initial guess.In certain cases bosonization can be used to make estimations of what the converged amplitudes would be, in particular where analytical control is good.However, in general such estimations will differ from the DMRG result and are often harder to obtain than simple heuristic guesses.
What is always available for initial values is that of previous solutions using the framework.Indeed, the best guess is the set of converged amplitudes leading to an immediate solution.This becomes useful when considering that mean-field amplitudes should often only change marginally when subject to a marginal shift in parameter space (notable exception is that of first-order transitions).Hence, if we seek a point in parameter space that is close to one which is already computed the converged amplitudes of the computed point serve as a guess which should be close to the converged result.
This observation is of particular use when considering the bond dimension.DMRG scales cubically in this quantity [31] (which may be improved using conserved quantum numbers).It is imperative to keep this parameter large enough in order for the obtained state to reasonably approximate the targeted state.With the MPS+MF framework it is possible to compute the selfconsistent amplitudes at a lower bond dimension starting with unguided guesses at low cost.When high precision is desired we use the cheap results to compute amplitudes at higher bond dimension.We find that the number of required iterations drops significantly with this strategy, as exemplified in Fig. 7.
Further, this strategy may be applied to the issue of phase transitions.For such cases a dense grid of data is often necessary with amplitudes varying only modestly between points.Once a grid of some sparsity has been generated a tighter one is cheaper due to the possibility of interpolating between existing points, providing good initial guesses.

D. Density-fixing
While the original Hamiltonian Eq. ( 1) is explicitly particle number conserving the derived effective 1D Hamiltonian Eq. ( 18) loses this property.Since carrier density has been shown to affect superconductivity it is important to consider what density the converged solution obtains.In practice, being able to fix the density is important, since there is no guarantee that converged solutions at different parameters have the same density thereby making comparison difficult.
In Fig. 3 two classes of amplitudes are considered: {α, β}, which represent the set of mean-field amplitudes, and µ, the chemical potential, which may be used to control the density of the system.However, there is no way to determine which chemical potential to use for a given density.In previous applications this issue was resolved by sampling a grid of exact data from exact diagonalization on small systems in order to obtain a range of µ in which interpolation to high precision was possible [36].
For models considered in this work the number of mean-field amplitudes is too great to attempt such a solution.Instead, a heuristic algorithm has been designed to obtain the appropriate density as shown in Fig. 8.In essence, the algorithm attempts to find a chemical potential µ target such that where computing the density n(µ) for different µ necessitates a full DMRG solution leading to a significant time cost.To alleviate this issue we assume that n(µ) is a monotonic function of µ and look for a range in which n(µ) = n target .Once obtained we use the secant method to narrow the range until precision is achieved.
It is notable that the density-fixing routine is most important when the mean-field amplitudes are strongly varying with each iteration.When approaching the self-consistently converged solution the density typically changes modestly, allowing for loops without running the density-fixing subroutine.

E. Self-consistent excited states
In a superconducting system, much information can be obtained from the energy gap between the ground and the first excited state.This is especially true in Q1D systems, as is shown in Sec.IV analytically and confirmed numerically in Sec.V C. Strikingly, finite-temperature properties such as superconducting T c can be obtained from this gap with minimal computational effort, which is readily exploited in Sec.VI.Thus, we employ the ability of DMRG to compute the lowest-lying eigenstates orthogonal to the ground state.This allows the computation of the first excited state energy where |ψ 1 is the first excited state which minimizes the energy E exc with the constraint that it is orthogonal to |ψ 0 , the ground state.The excitation gap may then be defined as The excited state will not generally obey the selfconsistency constraint, i.e., mean-field amplitudes measured in the excited state where α is the self-consistent amplitude obtained for the ground state.This means our numerics indicate a depletion of the ground state of condensed pairs, in line with earlier analytical theory [46].Depletion of pairing can be a concern as the Hamiltonian Eq. ( 18) does not conserve density.Using the above procedure will yield whatever state of lowest energy that is orthogonal to the ground state and could well be at another density, as is exemplified in Fig. 9.
At the same time, the data also shows the deviation decreasing with system size.An extrapolation may thus be performed using a general form of finite size behavior Thus, we assume that energy gaps extrapolated to infinite size systems are comparisons of energy at the same density despite the fact that finite size system densities generally differ.This strategy is verified in Sec.V C and Fig. 11, as two methods to obtain superconducting T c , one purely numeric, the other using ∆ and the field theory of Sec.IV, are shown to coincide.The upshot is that once a self-consistent Hamiltonian has been produced from MPS+MF for ground states, excited states will cost no more than one additional DMRGrun, making the computation of ∆ relatively cheap.

IV. FIELD THEORY DESCRIPTION
Let us now turn to a field theory analysis of the Hamiltonian Eq. ( 2).In 1D systems the effects of interactions are dramatically amplified.Additionally, in 1D the quantum and thermal fluctuations prevent the breaking of continuous symmetries [47].The combination of these effects leads to a unique universality class for interacting 1D quantum systems, known as Tomonaga-Luttinger liquids (TLL) [48].
The low-energy physics of TLLs can be described in terms of two bosonic fields φ and θ related to collective excitations of density and currents.These field are re-lated by the canonical relation which expresses the duality in 1D between density and phase fluctuations.In this bosonized representation, the single-particle fermionic operator of fermions with spin ν reads where ψ R,L,ν (x) are slowly varying field describing excitations close to the Fermi points ±k F (right and left movers) and ν =↑, ↓ denotes the spin.These fields are expressed as [4] where U r,ν are Klein factors, r = R, L and α is a cutoff proportional to the lattice spacing, which simulates a finite bandwidth.The fields φ ρ,σ are given by and similarly for the field θ.In the basis {θ ρ,σ , φ ρ,σ }, the 1D Hubbard model has the peculiarity of decoupling charge and spin sectors: Away from half-filling (one particle per site) and for negative interactions U < 0, the spin sector is gapped while the charge sector is gapless.In the bosonic language, the gapless Hamiltonian is the universal TLL Hamiltonian and the spectrum is linear ω = u ρ |k|.Such a Hamiltonian is fully described by two non-universal parameters which depend on the microscopic model: u ρ the charge velocity of the collective excitation and K ρ the Luttinger parameter that controls the algebraic decay of the correlations [4].On the contrary, the gapped sector is described by a sine-Gordon Hamiltonian which with respect to Eq. ( 32) has an additional term proportional to cos (2 √ 2φ σ ) that wants to lock the field φ σ to one of its minima: Physically, it means that the system tends to form Cooper pairs with opposite spins and as a direct consequence suppresses spin excitations.The energy of the bound state is the gap ∆ σ in the spin sector.Notably, ∆ σ is kept distinct in notation from ∆E s in Eq. ( 4) as their definitions differ.The field theory spin gap can be computed from the sine-Gordon model by various methods [4] and is exactly known for the microscopic attractive Hubbard model with U < 0 by Bethe ansatz [49].
As described in the previous sections we assume here that the spin gap ∆ σ is larger than the interchain coupling t ⊥ .Our system has thus a low and high-energy sector.The former relative to hopping in the transverse direction (t ⊥ t) and the latter to break the pairs.Let us first consider the case when ∆ σ t (or |U | t).In this case it is necessary to eliminate the spin sector before bosonizing the Hamiltonian.This can be achieved by a Schrieffer-Wolff transformation [50].The resulting Hamiltonian is rewritten as where H 0 is the 1D quadratic Hamiltonian Eq. (32).The effective coupling is now proportional to the gained energy over the cost of breaking the pair t 2 ⊥ /∆ σ and indeed it expresses local pairs hopping in the transverse direction.The spin excitations are exponentially suppressed, while the charge sector is massless.In Eq. ( 34), we neglect the terms corresponding to the formation of chargedensity-wave (CDW), because the superconductive correlation (SS) decays slower.Indeed, the corresponding correlations are [4] where the corresponding operators are defined as For negative U , the TLL charge parameter is larger than one K ρ > 1 , meaning that CDW formation is a sub-dominant instability.By using Eq. ( 34) and considering the leading terms n = n , since the spin gap is larger than the bandwidth (approximation of local pairs) we now use the mean-field approximation whose order parameter reads with C a constant that depends on the spin gap ∆ σ but is of order one in this regime (see Appendix G).In this limit, since only the charge sector survives, the problem maps onto a system of hard-core bosons (Cooper pairs) with a transverse hopping described by the field √ 2θ(x) instead of θ(x).Indeed, the physics is similar: the Pauli principle forbids to have two pairs in the same site and two hard-core bosons never occupy the same site.
All chains are now identical and, in the bosonized version, the 1D effective Hamiltonian reads ) where z c is the number of nearest neighbours in the transverse direction, L is the system size and ρ 0 is the unperturbed density.We find a sine-Gordon-like Hamiltonian and therefore, at T < T c , t ⊥ opens a gap in the spectrum because the cosine wants to lock the field θ ρ to one of the minima.In the thermodynamic limit, the zero temperature gap ∆ ρ is known analytically [51].This gap in the charge sector should equal the gap to the first excited state in Eq. ( 24) yet notation is kept distinct due to the differing definitions of the gaps.
The dimensional crossover [36] that occurs in such systems is represented by the mean-field critical temperature T c above which ψ † ↑ (x)ψ † ↓ (x) = 0 and the system is made of incoherent and decoupled 1D chains [46].This means that for T > T c , the thermal fluctuations wash out the transverse coherence due to the presence of t ⊥ .The system behaves essentially as if it was an isolated chain.Notably, the critical temperature scales like the charge gap at zero temperature, see appendix F for more details.Even though the prefactors in both T c and the charge gap are partially unknown, because of the constant C, the ratio is completely controlled by the Luttinger parameter K ρ only which can be computed from numerical calculations (DMRG, Bethe Ansatz, . . . ) with B(x, y) the Beta function and κ(K) a combination of gamma functions Γ(K), defined in Eq. (F6).
A more challenging case arises when |U | t.Indeed in that case since ∆ σ t the pairs are non local and the full bosonized form of the Hamiltonian with both charge and spin sectors (Eq.(32) and Eq. ( 33)) must be used.To deal with this situation we use a renormalization group (RG) procedure (see Appendix G), in which we eliminate all degrees of freedom from the initial bandwidth of the system, down to the spin gap.At that scale, since the running ultra-violet cutoff is now identical to the spin gap we are back to the situation where "local" pairs (hard-core bosons) hop in the transverse direction.The single particle hopping thus disappears at that scale and simply leads to an effective Josephson coupling between the various 1D units.This coupling can be computed from the RG as detailed in Appendix G.The limit ∆ σ W , where W ∼ 2t is the bandwidth of the 1D system, corresponds naturally to the case of weak interactions |U | W .In this regime, the spin gap either obtained from Bethe Ansatz for the Hubbard model, or more generally from the RG [52], is naturally exponentially small in the interactions ∆ σ ∼ |U |e −1/|U | .Thus, in order to be in the relevant regime for the whole procedure to work a very weak interchain hopping t ⊥ is needed.Specifically, the RG flow needs to be cut by the spin gap and it may not be cut by the interchain hopping.If this condition is satisfied we recover a model analogous to the one of the strong coupling albeit with a different Josephson parameter.As discussed above the ratio is insensitive to the precise value of the coupling and thus the ratio Eq. ( 38) is expected to hold for all values of U .

V. RESULTS
The primary focus of this section is to make use of and test the developed MPS+MF framework on the study case of negative-U Hubbard chain arrays.The methods described in Sec.II are tested thoroughly to produce a robust routine for determining critical temperature for the onset of superconductivity in Q1D systems of fermions with a gapped spin sector.These results are subsequently leveraged in Sec.VI to obtain T c for USC in weakly coupled doped Hubbard ladders with repulsive interactions.The negative-U Hubbard chain array results are split into 4 sub-sections which are ordered as follows: (A) The Hamiltonian Eq. ( 18) depends on a range parameter of the mean-fields which is studied in this sub-section, (B ) a numerical study of the ground state superconducting energy gap and critical temperature from thermal states, as well as (C ) comparing those numerical results with a more efficient mixture of analytics and numerics, and finally (D) benchmarking MPS+MF against AFQMC in a 2D system where the latter approach yields quasi-exact results.
As there are enough parameters to consider already, in the following results for chain arrays we are targeting a fixed density i.e., a quarter-filled system.We expect other close-lying densities will not yield markedly different results due to the nature of the isotropic case having weak dependence on density around this filling in the 2D case [53].

A. Hamiltonian range dependence
Considering the effective 1D Hamiltonian of Eq. ( 18) the first question to be answered is how long-range the mean-field terms can be.As discussed at the end of Sec.II, the MF-amplitudes will always decay with an exponential envelope as a function of distance between the operators appearing in them.Since longer range terms are more difficult to simulate, we want to find a minimal range for each parameter set with which longer-ranged Hamiltonians agree.As a metric for determining agreement we use the order parameter and energy gap ∆: The former defined by with r = 0 and where i f and i l chosen to avoid boundary effects of the open boundary conditions typical in DMRG.As can be seen from Fig. 10 the minimal range required varies with the strength of interaction.This is natural since weaker attraction makes electron pairs more dispersive and thus less localized.We note that longer-ranged terms can still be finite and ignoring them should yield at the very least a difference in ground state energies.However, the primary question is rather the degree at which these terms affect the ground state wave function and the critical temperature.In Fig. 10 we show that the order parameter and excited state gap ∆ do not change appreciably beyond the ranges that are displayed.While we find α-terms at longer ranges to be finite, ∆ and the onsite order parameter are largely unaffected.FIG. 10.Range dependence of the order parameter defined in Eq. ( 40) and first excitation energy defined in Eq. ( 24) for different values of interaction U/t = −2, −4, −10 and density n = 0.5.Notably the difference between the two largest ranges is not visible which indicates sufficient range in the Hamiltonian.

B. Numerical results
Using the minimal ranges we may compute ground and finite temperature states for the Hamiltonian Eq. ( 18) using the MPS+MF framework.As shown in Fig. 11a we find the critical temperature decreasing with transverse tunneling t ⊥ , vanishing as expected for t ⊥ → 0.
The excited state gap ∆ disappears in the same manner.Notably, the zero-temperature gap in field theory ∆ ρ should have the same meaning as the excited state gap ∆, which has been verified.Thus, from field theory it is expected that both T c and ∆ scale with t ⊥ with the same exponent, as can be seen from Eq. (F4) and Eq.(F5).
In determining T c numerically via state-purification within the MPS-approach, we have to contend with the increase in inverse temperature β = 1/T required as t ⊥ is decreased.This results not just in longer imaginarytime evolutions, but also in increased finite size effects (see App. A).
Reaching large system sizes for these thermal state calculations can be challenging.We used two different approaches.In the first, used for U = −10t, the infinitetemperature purified starting state of the imaginary time evolution is constructed to be in the S = 0 subspace, which allows to keep exploiting the conserved spin quantum number during the evolution.The drawback is that this β = 0 initial state has entanglement growing strongly with system size, limiting the practically attainable system lengths.
The second approach, employed for U = −2t, −4t, is to sacrifice the spin-conservation.This makes the purified initial state into a trivial-to-construct product state, and thus arbitrary system lengths are accessible.However, this makes the imaginary time evolution more costly, as there are no conserved quantum numbers anymore.In order to alleviate this issue we use the PP-DMRG framework [54] expounded upon in App.B. We moreover rely on the ability of H MF to filter out the S = 0 subspace at temperatures at and below T c .We ascertain this to be correct, with violations of S z = 0 reaching at most 10 −2 , and typically much less, across all calculations.

C. Numerical-analytical hybrid results
For the field theory of Sec.IV it is difficult to quantify the superconducting T c or the excited state gap ∆, due to the unknown pre-factors arising from the massive spinsector.However, forming the ratio ∆/T c = R(K ρ ) is free from these unknown pre-factors, depending just on the TLL parameter K ρ and is, strikingly, constant in t ⊥ .
Generating R from the numerically determined data of Fig. 11a and the analytical expression Eq. ( 38) yields Fig. 11b.Notably, while the analytical ratio is not agreeing with numerical estimates exactly, the constant nature of the analytics is likely approximative.Achieving a ratio which lies close enough to the data is sufficient in order to obtain critical temperatures.With this knowledge the  new T c estimate becomes Since the primary issues of T c computation came from the imaginary time evolution we may now obtain T c estimates from ground-state DMRG by computing ∆(T = 0).As shown in Fig. 11a the estimation scheme Eq. ( 41) is agreeing very well with the numerical T c values obtained from thermal state calculations.
Since the more efficient ground-state DMRG is the main tool of this alternative way to obtain T c , it can be brought to even lower values of t ⊥ and greater system sizes, allowing for greater precision at larger parameter ranges.

D. Comparison with AFQMC
In Eq. ( 18) the coordination number z c tracks the dimensionality of the underlying Q1D array.For the calculations on 3D Q1D systems, as performed so far, we have z c = 4. Lattices of other dimension can be simulated just as well, just by changing z c to the appropriate value.
We exploit this for benchmarking the MPS+MF approach against quasi-exact results, which AFQMC is able to obtain in the absence of a sign problem.These benchmarks are necessarily done for 2D models, as the finite temperature algorithm used scales cubically in the number of lattice sites [40].We stress that the MPS+MF method describes ordering in the Ginzburg-Landau sense, via the pairing mean-fields α ik , while AFQMC detects the actual 2D BKT transition.We interpret the mean-field T c of MPS+MF as an approximation of the T BKT which occurs for the 2D model.
We compare the two algorithms both for ground state and finite temperature calculations.For ground states, which do achieve superconducting LRO even in 2D, the onsite order parameter is compared between the two algorithms as shown in Fig. 12a.Since 2D is the lower critical dimension, quantum fluctuations around any mean will be especially strong.Considering this fact we note that the smallest simulated t ⊥ values have a modest overestimation as MPS+MF neglects transverse quantum fluctuations by design.
For finite temperature states we compute the BKT transition termperature, T BKT , using AFQMC and com- pare it to T c from MPS+MF.The strategy for obtaining the T BKT temperature from AFQMC is standard and elaborated in App.D. Both T BKT and T c , as well as their ratio, are shown in Fig. 12b.Over the range of t ⊥ for which we simulate, the near-constant ratio between T c and T BKT is striking, being approximately This is in line with previous work where we also found such a ratio to be robust to changes in parameters [36].Similarly to the comparison to zerotemperature AFQMC, MPS+MF will overestimate the transition due to the neglect of both quantum and thermal fluctuations, and these will again be especially pronounced, given that these are 2D systems.Previous work on bosonic systems, as well as the generally known dependence of phase transitions on spacial dimensionality of a system, indicates that these will be strongly reduced for a Q1D 3D system [36].Specifically, there we found that T QMC c /T MPS+MF c ≈ 0.7.For 3D fermionic Q1D systems we thus expect that a correction factor for T c computed from the MPS+MF framework to obtain the true T c will lie somewhere between these two extremes, and probably closer to 0.7 than to 0.25.

VI. 3D ARRAY OF WEAKLY DOPED REPULSIVE-U HUBBARD LADDER
With the MPS+MF framework for fermions developed on 3D Q1D systems of negative-U Hubbard chains this section applies it to a much more demanding system: 3D arrays of weakly coupled, doped, repulsive-U Hubbardladders.These systems have been investigated via field theory as an alternative to Q2D systems in the study of USC and high-T c superconductivity.The microscopic mechanism of repulsively mediated pairing is understood from field theory, and the strength of this pairing can be quantified reliably via MPS-based methods.Despite these critical advantages over the Q2D models, there was no quantitative method to study these 3D arrays that incorporates the pairing physics at the microscopic level, as the TLL field theory in practice yields largely qualitative results.
The MPS+MF framework supplies that ability.Just like the negative-U Hubbard chains, the isolated, doped repulsive-U Hubbard ladders have finite ∆E s and ∆E p , manifesting the repulsively mediated pairing.Further, when analyzed via renormalization group theory within the TLL approach, the low-energy physics of both these 1D sub-units is structurally analogous.Both exhibit an ungapped charge sector, characterized by a TLL coefficient K ρ for the chain, and K ρ,+ for the ladder.Furthermore, both can be computed from the microscopic Hamiltonians via MPS-methods.For the Hubbard ladders, there are three additional gapped modes, a charge one and two spin ones, the smaller having minimal energy ∆E s , as for the single spin mode in the negative-U chains.
The MPS+MF framework thus can be applied to the Hubbard-ladder arrays.For a proof-of-principle treatment, we will focus on plain Hubbard ladders depicted in Fig. 1b, with U = 8t and average density fixed at n = 0.9375.The characteristic energy scales and TLL parameters can be computed with DMRG which results in ∆E s ≈ 0.078t, ∆E p ≈ 0.134t and K ρ,+ = 0.77 (the latter extracted from [6]).The study of optimized Hubbardladder arrays, engineered for high T c 's via deliberate optimization of ∆E s , ∆E p and K ρ , as well as examining the possibility of charge-density order competing with USC within two-channel MPS+MF, is the subject of forthcoming future work.
As ladder geometries require much larger MPS resources than chains, direct calculation of T c , while feasible, is challenging.We thus use MPS+MF combined with the analytics developed and tested in Sec.II-Sec.V.In this manner we obtain T c for USC in these arrays by computing excited state energy gaps ∆ using DMRG for ground states, then applying Eq. ( 41).Compared to Sec.IV and Sec.V, this procedure requires only marginal adjustments for the ladder as negative-U chain and Hubbard ladder look largely identical in the lowenergy parts of their respective field theories.For the former, the order-parameter scales as e −i √ 2θρ(x) in the phase-operator of the ungapped charge mode, while for the latter it is e −iθρ,+(x) .From that, it follows that the ratio-function for the ladder is where R is given by Eq. ( 38).Thus, R ladder will retain its dependence on a single TLL parameter: K ρ,+ .This is due to the gapped spin-sectors entering the orderparameter in the same manner both for T c and ∆.These non-universal contributions thus cancel when forming R, analogous to the derivation in Sec.IV and App.F. The same analysis yields that for the ladder we have The effective MF-Hamiltonian for the Q1D array of the ladders is given by where i, j are leg and rung indices respectively and H P air and H P H are derived in App.E and defined by and the pairing amplitudes are given by whereas the exchange terms are given by In this work we exclude the possibility for a CDW phase such that density is independent of which rung we measure and is constant throughout the system.This is obtained by the restriction β i1i1,j1j1,σ = β i2i2,j1j1,σ .We note that inclusion of the exchange terms β ii ,ll ,σ has previously not been possible by analytical methods.
In isolation, Hubbard ladders typically require large bond dimensions for accurate simulations (see, e.g., [6]).With the included superconducting MF ordering channel we have found this requirement to be relaxed.This is expected, as any long-range order, be it in real space or momentum space, requires fewer retained Schmidt components than that same system without such order.However, at these lowered bond dimensions, converged MF-amplitudes exhibit non-negligible dependence on the bond dimension despite modest truncation errors.Regardless of this, the linear scaling of energy with truncation error, typically found in DMRG, remains intact even here, as shown in Fig. 13a-c.Thus, it is possible to compute the mean-field amplitudes at modest bond dimension and obtain energy measurements, and thus ∆, extrapolated to zero truncation error at different system sizes.
Finally, we obtain ∆, and thus T c after rescaling with Eq. ( 43), in the thermodynamic regime via infinitesize extrapolation as shown in Fig. 13d.We note that t ⊥ is chosen quite close to ∆E p and ∆E s .The primary reason are the small energy gaps we need to resolve.At smaller values of t ⊥ the discrete energy gaps of the finitelength systems mask ∆ even for ladders with L = 80 or even larger.However, using Eq. ( 44) it is possible to extrapolate thermodynamic T c to smaller, physically more reasonable values of t ⊥ as shown in Fig. 13e.

VII. RESOURCE REQUIREMENTS
In the present section we summarize time and resource consumption for the different algorithms used in this work.The MPS+MF routines feature repeated DMRG solutions.Generally, the total number of CPU-core hours τ scales with the same parameters as DMRG does, i.e., where L is the number of lattice sites (and thus MPStensors), d the size of the local Hilbert space, and χ is R ladder (K ρ,+ ) extrapolated to infinite size where R ladder (Kρ,+) is given by Eq. ( 43) and Kρ,+ = 0.77 [6], while (e) shows the critical temperature vs. t ⊥ exploiting the scaling known from Eq. (F4), using the Tc value obtained in sub-figure (d) (red cross).
bond dimension.The number of loops required to reach self-consistency, N tot is shown for representative examples in Fig. 5 and subject to the optimizations mentioned in Sec.III C. The total resources consumed by the algorithm are thus given by For the ground states of the negative-U Hubbard chain with MF-amplitudes we have used MPS-based DMRG from the Matrix Product Toolkit package [55].The data has been computed using Intel Xeon E5 2630 v4 at 2.20 GHz CPU-cores.Most results are obtained at a bond dimension of χ = 300.A certain speed-up was obtained using 2 CPU-cores and threads used in the LAPACK and BLAS routines on which the algorithm rests.For a typical run a single MF-loop takes about 8000s of wallclock time for χ = 300 and L = 100.Notably d = 4 for the negative-U Hubbard chain.With 2 CPU-cores in use for this calculation, the total required resources for solution are τ M P S+M F ∼ 5N tot CPU core-hours.Additionally, in order to obtain critical temperatures directly we have calculated thermal states using imaginary time evolution of purified states.For simpler Hamiltonians where a shorter range (maximum of 1) for both α ik and β irσ is possible, trotterized time evolution suffices [31,43].For this case we have used a time step of δτ = 0.1 and a fourth order Trotter discretization.With additional linear scaling in the length of the imaginary time simulated (equal to half of inverse temperature β) a typical solution requires about 12000s of wallclock time for L = 60, χ = 200 and β = 9.5.Having used 2 CPU-cores, the required resources in total are τ M P S+M F ∼ 7N tot CPU-core hours.Notably, the range of inverse temperature β to be simulated to determine T c changes markedly with t ⊥ , leading to commensurate changes in τ .
In the case of Hamiltonians with long-range terms it is necessary to apply more advanced time evolution schemes.With that the resource requirements increase significantly, i.e., a typical solution requires about 32000s on 6 CPU cores for L = 60, β = 9.5, and χ = 100, leading to an overall τ M P S+M F = 53N tot CPU-core hours.In part, this increase is due to the longer-range couplings.Another cause is the increased effective system size.This increase in turn is down to the specific MPS-implementation [56] used, which always requires full quantum number conservation.In our usecase, where both charge-and spin-conservation are discarded (c.f.subsec.V B), quantum-numbers are restored artificially via the use of projected purified density-matrix renormalization group (PP-DMRG) [54], at the price of a larger effective system.Further details on obtaining finite temperature results are provided in the Apps.A and B. Exact ground-state order parameters in finite lattices were computed with an AFQMC method using generalized Metropolis with force bias [57].This algorithm scales quadratically with the number of electrons, N e , and linearly with the number of lattice sites N l : For a N l = 40 • 4 = 160 system with 80 electrons, we performed calculations of 2000 sweeps for measurements after 10 sweeps of thermalization, with 100 independent repeats.Such a calculation, with an imaginary propagation time β set to 64, took 49 hours of wallclock time on 100 Intel Xeon E5-2640 v4 2.4 GHz cores yielding a requirement of τ AF QM C,GS 4900 CPU-core hours.
In order to obtain exact values for T BKT in twodimensional Q1D for benching MPS+MF against (c.f.subsec.V D), we have used finite temperature AFQMC.The package utilized for this purpose is called Algorithms for Lattice Fermions (ALF) [40].Data was obtained with an algorithm which scales cubically in the number of lattice sites, N , and linearly with inverse temperature β: An improvement of this scaling is available [58] but was not implemented for our calculations.We find that a single instance of sampling requires 200s for N = 48 • 8 = 384 and β = 20 using a Intel Xeon E5-2698 v3 2.30GHz CPU-core.Statistical error bars are sufficiently small for a sample size of ∼ 100000 for the parameter set we study, requiring τ AF QM C 5000 CPU-core hours.The ALF-package, like most QMCimplementations, can of course parallelize this workload near-perfectly.
When performing MPS+MF on the weakly coupled repulsive-U Hubbard ladders, MPS+MF becomes more resource-intensive: as earlier, one conserved quantum number is lost and the Matrix Product Operator representing the Hamiltonian is significantly larger than any for the negative-U Hubbard chains.We use a DMRGimplementation offering distributed-memory parallelism (pDMRG) to obtain faster solutions [14].Different from QMC-type algorithms, any such parallelisation will inevitably show non-trivial communication overheads, and thus not scale linearly in the number of MPI-processes.For a single converged ground state from pDMRG we thus require 27000s using 32 • 8 = 256 Intel Xeon E5-2698 v3 2.30GHz CPU-cores, at χ = 1000 and L = 96.The total cost thus becomes τ ∼ 2000N tot CPU corehours.The optimizations of Sec.III C particularly apply to the Hubbard ladder systems.Without these optimizations, N tot can be as high as N tot = 25, but by employing them, where possible, can drop as low as N tot = 6.

VIII. CONCLUSION
In this work we have developed a numerical framework combining MPS-based numerics with MF and perturbation theory to solve correlated quasi-one-dimenional fermionic systems, constructed out of weakly coupled 1D sub-units, in two and three spacial dimensions.This method relies on MF-approximating tunneling processes occurring transverse to the 1D sub-units with amplitude t ⊥ .The requirement for this approximation being reasonable is that t ⊥ be weaker than any gap on the 1D sub-units that suppresses first-order tunneling between 1D sub-units.Using the example of superconductivity in such Q1D arrays, we show how this framework allows to map otherwise difficult or even intractable correlatedfermion models in 2D and 3D onto a self-consistent 1D problem.We then demonstrate how these can be effectively solved both for ground states and thermal states.
We test the framework on a model of attractive fermions on a 1D chain extensively comparing to both analytical methods and AFQMC.We obtain analytical expressions for superconducting T c of the model and the gap, ∆, to the first excited state.Utilizing that a ratio R = ∆/T c of these two quantities remains constant over t ⊥ , we obtain R analytically and greatly speed up calculation of T c via the use of ∆ and R. Comparing this value with T c obtained directly from thermal state calculations shows that obtaining T c from ∆ and R yields excellent agreement, especially at small t ⊥ .This allows MPS+MF to obtain T c without using imaginary time evolution, which is numerically more costly than obtaining ∆ via ground state calculations.
Subsequently, we use the gap and ratio method to obtain T c from MPS+MF and compare with T BKT from AFQMC.We find a semi-constant ratio of these over a range of t ⊥ , in line with previous comparisons to QMC [36].As expected, the MF-approximation yields greater over-estimation of the ordering temperature for lower-dimensional systems.With this in mind the method seems able to provide reliable estimates of T c in fermionic systems for an appropriate choice of parameters.
Utilizing the developed MPS+MF framework we treat the case of a 3D array of weakly coupled, doped Hub-bard ladders with strong repulsive interaction.With the tools developed in this work we are able to calculate T c quantitatively for the first time for these systems.The MPS+MF framework thus allows the simulation of a subgroup of 3D systems of strongly correlated fermions, namely the Q1D models, which were previously out of reach for any quantitative method.
Notably, the 3D arrays of weakly coupled Hubbard ladders studied in this work have not been optimized to yield large critical temperatures.Previous work has indicated that by modifying the ladder-parameters larger T c 's may be achieved [5].With the MPS+MF framework it is possible to systematically search for improved critical temperatures starting from the microscopic models.This allows not just to deliberately search for optimal high-T c prototype-materials in the Q1D space, something that remains elusive for Q2D materials.It likewise permits to design ultra cold gas experiments capable of observing analogue states of high-T c superconductivity within current or near-future experimental constraints.We are pursuing both possibilities in current follow-up work.
In this work we have focused on the physics of superconductivity using the MPS+MF routine.However, MPS+MF can be used for any Q1D system in which tunneling in-between 1D sub-units is suppressed at first order by a gap .gap can be of any physical nature, such as, e.g., the charge gap present in the insulating phases of the Bechgaard and Fabre salts.The MPS+MF framework can thus also be deployed to understand, e.g., the antiferromagnet to spin density wave transition in these materials.This potential application of the framework highlights its power to incorporate multiple ordering channels simultaneously at the mean-field level, and thus its ability to resolve the competition between competing orders.
Finally, the capacity of MPS-numerics to address realtime dynamics of many-body systems both near and far from equilibrium opens the possibility to use the MPS+MF algorithm to study real-time dynamics of correlated fermions in high-dimensional Q1D systems.Such forthcoming work is currently in preparation on dynamically induced superconductivity in such systems [41].The Q1D models we consider in this paper suffer from finite size effects like all numerics on finite systems.In particular, as the connection between 1D systems weaken so does the strength of the resulting superconductor resulting in increased healing lengths.Additionally, we have found that finite size effects persist even at large t ⊥ albeit reduced in size.In order to accurately simulate these systems extrapolation has to be performed on the finite size observables to the limit of infinite size.
Using MPS-numerics we compute the observables to be measured for several system sizes.We then utilize various strategies to obtain infinte size values depending on the type of finite size effect and observable.The strategies are outlined in this appendix.

Local observables
For local quantities such as the order parameter and energy we use a common heuristic where O(L) is a measurement of any local observable for a system size L and O ∞ the thermodynamic limit of that observable.Thus, we fit measurements at finite size to a quadratic polynomial in inverse size.We find data fits the heuristic pretty well as shown in Fig. 14.

Finite temperature at criticality
Several results in this work obtain the critical temperature of a system by evolving in imaginary time.The phase transition point is dependent on system size.For the case of significant finite size effects we consider the critical behaviour of the superconducting order parameter in particular as that determines when the system enters superconductivity.

Significant finite size effects
When finite size effects must be considered we follow a common strategy used, e.g., on QMC results [59].We study a second-order phase transition, where the critical behaviour of the order parameter (here named ψ for simpler notation) is given by The reduced temperature, t, is given by and critical exponents β, ν are given by where ξ is the correlation length of the ordering field.Notably, on the unordered side the order parameter is zero.In order to extract T c we assume that our systems critical behaviour belongs to the mean-field universality class which is consistent with fits to Eq. (A4) close to transition such that β = ν = 0.5.Using these exponents the finite size order parameter is rescaled by L β/ν and plotted over L 1/ν t The function ψ(x) is the scaling function of the order parameter and is system size independent.Thus, for a correct choice of T c all curves overlap close to transition as shown in Fig. 15.The quality of the collapse has been determined using where N is the number of system sizes, s 2 ψ is the variance integrated over a range [x min , x max ] and ψL (x) is the scaling function for size L given a T c .The range is chosen to lie around the proposed value of T c .The critical temperature is then obtained for different widths of the integration interval and extrapolated to zero.The final error bar shown in values for T c is taken to be the fitting error for decreasing integration range.This ultimately gives a small error bar.Added data for larger sizes might move the result outside this error bar.
For parameter sets where t ⊥ is particularly small generating data close to transition is tricky due to slow convergence of the mean-field amplitudes.One consequence of this is that the data used for collapse can end up too far from transition for finite size scaling to apply.On such occasion the analysis fails to produce reliable collapse of data and another strategy is needed.

Critical temperature interval
When the previous analysis fails to produce a reasonable collapse we generate a grid around the estimated phase transition.Since the convergence close to the phase transition is especially demanding, we focus on surrounding temperatures.In order to reduce the number of necessary self-consistency iterations we extrapolate the order parameter via O(n) = O n→∞ + a • exp(−bn) in the number of self-consistency iterations n.Those extrapolations are shown as an example in the inset of Fig. 16 for L = 50, but for convenience plotted over the inverse number of iterations.The results of these extrapolations are then used within to fit the data with as can be seen in Fig. 16.The final result interval for the critical temperature is then given by the estimated values Example obtaining Tc from, fits without finite size scaling for = −4t, t ⊥ = 0.4t, n = 0.5.Symbols denote the onsite order parameter extrapolated to infinite number of iterations, for L = 20 (plus), L = 50 (cross) and L = 60 (star) computed via MPS+MF for thermal states.Solid lines are guide to the eye.Dashed lines are fits with ansatz Eq. (A8).Inset shows an example of extrapolating the order parameter to infinite numbers of iterations for L = 50 Stucture of MPS thermal state calculations using state purification and imaginary time evolution via TDVP, which requires introducing ancilla sites.The use of the SymMPS-package [56] for TDVP also requires using PP-DMRG [54], meaning adding auxiliary sites for both physical and ancilla sites to recover conserved charge and spin quantum numbers.
for different system sizes.This procedure causes significant errors compared to the previous method and fails to account for finite size effects that may occur.Nevertheless, the errors are small enough to permit analysis as can be seen in Fig. 11.served U (1) quantum numbers are necessary.Since the MPS+MF Hamiltonian does violate those symmetries, we applied PP-DMRG [54] to restore them.That means we not only double the system size in order to represent density matrices instead of states, as it is usually the case in state purification [60,61], but also double it again in order to have an auxiliary bath to restore the conservation-laws via these.Hence, as can be seen in Fig. 17, we have the physical and the ancilla system, which represent density matrices and the auxiliary and the ancilla-auxiliary system, which restore the quantum number conservation.
Since the long-range interactions are thus increased in their range by a factor of four, another obstacle needs to be circumvented: this is the loss of particles from the physical (and its auxiliary) system into the ancilla (and its auxiliary) system [43].This leakage occurs due to the accummulation of numerical errors, which we prevent by having the U (1) symmetries conserved on the physical and the ancilla system separately.To achieve that we increased the local basis from four to 16 states, by adding a separate fermionic particle species, which is supposed to only occur on the ancilla systems.
Evolving a product state, such as an infinite temperature state, with single-site time-dependent variational principle (TDVP) lead to significant projection errors [43].Hence one needs to increase the bond dimension first via a different time-evolution method [62,63].We chose to use two-site TDVP for the first ten time steps, which are chosen to be very small δt = 10 −7 .This way the bond dimension grows, but the projection error, which is scaled within the exponential by the time step, stays small.Afterwards one time step is performed to go to the usual time step grid δt = 0.05.Then all following time steps are executed by single-site TDVP, since that is faster.
apply Trotter-Suzuki breakup for each small imaginary step, where K is the kinetic part of the Hamiltonian, containing only one-body terms, while V is the potential part, which consists of two-body terms.
To rewrite all two-body terms into one-body terms, we apply the Hubbard-Stratonovich (HS) transformation in a charge decomposition form: where cosh(γ) = exp(−∆τ U/2).With the above transformation, the short-time propagator can be written as where B(x) = e −∆τ K/2 i bi (x i ) e −∆τ K/2 is now a onebody propagator, and p(x) is a probability density function, which in the form above is uniform in the auxiliaryfield (AF) configurations x = {x 1 , x 2 , ..., x N l }.
Ground state observables are given by where In our calculations with the generalized Metropolis algorithm, we choose a total projection time β, which defines a fixed length of the imaginarytime path integral.The location along the path where Ô is measured moves with our sampling; for example, as we sweep from left to right, we start measuring when β L > β eq , where β eq < β/2 is a parameter which ensures that the asymptotic limit in Eq. (D1) is reached (in a numerical sense), and β R ≡ β − β L .Conversely, when we sweep from right to left, the measurement starts when β R > β eq and stops when β L ≡ β − β R reaches β eq .The expectation O is expressed as path integrals in AF space: ), with M R ≡ β R /∆τ and M L ≡ β L /∆τ .A heat-bath like algorithm and a cluster update scheme are incorporated in our generalized Metropolis algorithm, which is described in detail in Appendix A of Ref. [57].
The pair-correlator (the pair-pair correlation) can be computed by the path-integral above in Eq. (D7).For each path, if we denote the expectation value be decomposed by Wick's theorem into pair products of one-body Green's functions: where Φ L and Φ R are the matrix representation of the kets |φ L and |φ R , respectively.We choose a reference site and then compute the paircorrelator between the reference site and different other lattice sites, as is shown in Fig. 18.We then average the pair-correlation of the sites that have the longest distance.To further reduce statistical error, the reference point is averaged over the whole lattice, since each lattice site is equivalent under periodic boundary condition.

Kosterlitz-Thouless transition temperature
For finite temperature results we use ALF.The structure of the algorithm is similar to that of ground state AFQMC previously outlined and we refer the reader to the ALF documentation [40].
Using ALF we obtain results for any value of imaginary time β.From linear response theory it is possible to relate the superfluid weight to current-current correlators and kinetic energy [66]: For the KT-transition we expect that [68] lim which is then extrapolated to the thermodynamic limit using the form [69] T An example of this fit is shown in Fig. 19b.

Density
The AFQMC finite temperature algorithm from the ALF collaboration uses the Blancenbecler-Scalapino-Sugar (BSS) algorithm [40,70].This algorithm works in the grand canonical ensemble as is necessary and will not have fixed density.In general, we are interested in specifying the density to work at as T BKT will have some dependence on this quantity.
Precisely fixing the density requires simulation of a large number of chemical potential values.In order to alleviate this problem we run simulations for a small number of lattice points and determine the correct chemical potential for a given temperature.We then use this value of chemical potential for all lattice sizes and temperatures of that parameter set.This will yield a notable error in density as shown in Fig. 20.At the same time we find that T BKT is only modestly affected by density.

Appendix E: Array of Hubbard ladders
We derive the effective 1D model for the weakly coupled Hubbard ladders starting from the 3D array given by where The vector, R kl , denotes the position of the ladder in a 2D grid.We have used the definition (R kl suppressed) and c ijσ follow the usual anti-commutation relations.
So far, we have only described a set of Fermi-Hubbard ladders.The added Hamiltonian H ⊥ is defined by Note that movement to neighbouring ladders may change what leg you are on.This is due to half of neighbouring ladders being side-by-side neighbours and the other half are front-and-back neighbours.

Effective Hubbard ladder Hamiltonian
When U is strongly repulsive and density is close to unit filling H HL and thus also the full set of ladders have a spectrum which contains clusters of energy eigenstates separated by large gaps.Thus, analogous to Eq. ( 9) it is possible to derive an effective Hamiltonian where ∆E p = 2E(N + 1, 1/2) − E(N, 0) − E(N + 2, 0), (E5) is the pairing energy for a single ladder at particle number N and E(N, S) is the energy of a ladder at particle number N and spin S. The operator P 0 is a projector to the lowest energy manifold of the total system.As in Sec.II A, this removes certain terms within H 2 ⊥ , such as two particles moving to two separate ladders.
Expanding P 0 H 2 ⊥ P 0 yields a new operator which is quartic and acts like an effective interaction.Each interaction involves particles on two different ladders, e.g., moving two particles from one ladder to an adjacent one.

Mean-field Hamiltonian
With a quartic interaction where half the operators involve one ladder and the other half involve the other we can make an ansatz of quasi-free states: (E6) This allows us to create a mean-field Hamiltonian which produces expectation values of this form.In the process we assume that expectation values involving operators on different ladders are of negligible size and ignore them.This leads to the mean-field Hamiltonian (for one ladder) where ii ,jj α ii ,jj c † ij↓ c † i j ↑ + c i j ↑ c ij↓ , (E8) ii ,jj ,σ β ii ,jj ,σ c † ijσ c i j σ , (E9) and the pairing amplitudes are given by whereas the exchange terms are given by Note that the hermiticity of Eq. (E9) is not apparent from the expression but is hidden in the sum.
1. First step RG: renomalization of charge and spins In order to find the RG equations, we compute perturbatively (in the couplings) the correlation ψ † (x)ψ(x) [4,71].Let us start from the Luttinger parameter K σ and the dimensionless interaction g = U πv F , with v F the Fermi velocity dK σ (l) dl = − K 2 σ (l)g 2 (l) 2 dg(l) dl = (2 − 2K σ (l))g(l) (G1) Because the system is spin-isotropic, the equations are equivalent to the ones from the XY problem [72].For U < 0, we flow towards larger g and smaller K σ and we need to stop the flow when g ∝ O(1), say at the RG length l 1 .This fictitious length is defined from α(l) = α(l = 0)e l with α(l = 0) the original cut-off (lattice spacing).
In presence of interchain tunnelling we need to complete the above equations by the ones generated by the interchain tunnelling: where Kν = K ν + K ν −1 .The dimensionless couplings are defined as J = πα 2 4u ρ (ρ 0 A F ) 2 J and Js = α 2 2u 2 ρ t 2 ⊥ (G3) and the subscript s stands for source term.Note that the transverse hopping is also, in principle, contributing to the renormalization of the other parameters Eq. (G1) and of K ρ .In practice, because we consider that the interchain tunnelling is the smallest quantity in the problem and in particular that t ⊥ ∆ σ , we will neglect such renormalization in the first step of the RG.In particular K ρ can be considered as essentially constant in the first step of the RG.Note also that the combination of interchain hopping and interactions lead to a modification of the naive scaling of the interchain tunnelling.In addition to its own renormalization, the single particle interchain tunnelling also generates via RG the pair tunnelling.This is due to the fact that pairs of electrons that hop within a distance |r 1 − r 2 | < α(l) are to be considered local.Moreover, we need to enforce the condition J(l = 0) = 0 because the original Hamiltonian Eq. ( 1) has only single-particle hopping, not pair hopping.It is clear that, as we renormalize, the pair-hopping term is generated and eventually will be the relevant coupling.It easy to see that if from Eq. (G1) we flow towards smaller K σ , then the 1/K σ term in Eq. (G2) makes t ⊥ an irrelevant coupling.
We have to stop this first step of the flow when the coupling constant in Eq. (G1) is of order one.At that scale (l = l 1 ) the microscopic cutoff is of the order of the pair size.Another estimation of l 1 , is the pair size α(l 1 ) ∼ u σ /∆ σ where ∆ σ is the spin gap and u σ is the velocity of the spin sector, before pairs become local.At that scale the single particle tunnelling is suppressed because of the gap in the spin sector which can be formally seen in the RG equations Eq. (G2) by the fact that for small K σ , t ⊥ is formally irrelevant.

Second step RG: pair hopping and dimensional crossover
In the second RG step l > l 1 , we have only pairhopping and spin excitations and single particle hopping are suppressed.The spin sector is out of the picture, and the RG equation expressing the pair-hopping coupling becomes We thus see that at the scale l 1 we are now left with only the charge sector as a massless sector and an effective Josephson coupling between the chains.The situation is thus similar to the one we had in the large spin gap limit but with a different Josephson coupling than the strong coupling limit J ∼ t 2 ⊥ ∆σ .This has consequences for the absolute values of the T c and charge gap at zero temperature but the ratio is unchanged compared to Eq. (38).
The absolute values of the T c or the charge gap can simply be computed by continuing the flow of Eq. (G4) until the Josephson coupling itself become of order one.This defines a second RG length l * .The condition is set by J(l * ) ∝ O(1), which means that for l > l * the coupling is so large that we are back to a 2D/3D system.From this length, we can define either the critical temperature or the charge gap, because their ratio is fixed from Eq. (38).By integrating Eq. (G4),, we find that the RG length l * at which the DC occurs is α(l * ) = α(l 1 ) 1 J(l 1 ) with α(l 1 ) the pair size.
Finally let us note that if we look directly at the pair operator we have Averaging over the massive spin sector would lead to where the prefactor C could be related to Because of the gap in the spin sector this average is nonzero.

FIG. 1 .
FIG. 1. Overview of MPS+MF for fermions.A 3D array composed of weakly coupled 1D sub-units is mapped onto a selfconsistent 1D problem with mean-field amplitudes α ik and βirσ.(a) Array of negative-U Hubbard chains, with onsite attraction U , tunneling t along the chain and inter-chain-tunneling t ⊥ .(b) Array of repulsive-U doped Hubbard ladders.Doping levels are δ, n.n.-tunneling inside the ladders is t and t ⊥ in-between them.

FIG. 2 .
FIG.2.Schematic representation of the spectrum for Hamiltonians of form Eq. (1) which are considered.The gaps highlighted are the spin gap ∆Es and pairing energy ∆Ep.Notably, the spin gap and pairing energy are approximated as independent of small variations in energy.

FIG. 4 .
FIG.4.Example of truncation error extrapolation of order parameter and ground state energy for a chain of length L = 100 at attractive interaction U = −4t and density n = 0.5 and t ⊥ = 0.05t.The y-axes represent difference of order parameter δα = α( ψ ) − α(0) and similarly for energy.The range of MF terms taken into account is r = 4.

FIG. 8 .FIG. 9 .
FIG. 8. Schematic representation of the density fixing routine.Chemical potential at different stages of the loop are named µi and density ni.

FIG. 13 .
FIG. 13.Energy gap to first excited state and critical temperature of the Hubbard ladder for U = 8t, n = 0.9375 and t ⊥ = 0.0489t.The sub-figures show ground state energy E0 (blue circles) and excited state energy E1 (red circles) extrapolated to zero truncation error (blue/red lines) for a ladder with length (a) L = 80, (b) L = 96, (c) L = 112.Sub-figure (d) shows the critical temperature of the 3D ladder array Tc = ∆/R ladder (Kρ,+) = E 1 −E 0
FIG.16.Example obtaining Tc from, fits without finite size scaling for = −4t, t ⊥ = 0.4t, n = 0.5.Symbols denote the onsite order parameter extrapolated to infinite number of iterations, for L = 20 (plus), L = 50 (cross) and L = 60 (star) computed via MPS+MF for thermal states.Solid lines are guide to the eye.Dashed lines are fits with ansatz Eq. (A8).Inset shows an example of extrapolating the order parameter to infinite numbers of iterations for L = 50

1 )
}, which is an M N l -dimensional vector, denotes the AF configuration of the entire path, with M ≡ β/∆τ being the number of time slices in the path, and the probability function P (X) = M m=1 p(x (m) ).The wave functions |φ R and |φ L are single Slater determinants (if we choose |Ψ I to be a single Slater determinant), and have the form
Schematic representation of the MPS+MF framework.Green boxes highlight the density fixation part of the routine.Red boxes denote the mean-field amplitude fixation part.The mean-field amplitude set is denoted by {α, β} where {α} represent pairing amplitudes and {β} the exchange amplitudes.
(38)les show data obtained via direct calculation of thermal states, while crosses show Tc computed from zero-temperature calculations, using the excitation gap ∆ and Eq.(38).Sub-figure (b) shows the ratio R = ∆/Tc of ground state excitation gap to critical temperature, both computed separately, compared with the analytical ratio Eq. (38) (dashed lines) vs t 2 ⊥ /∆Ep.
FIG. 11.A comparison of analytical methods and combination with numerical methods for several values of interaction.Subfigure (a) shows critical temperature Tc vs transverse tunneling t 2 ⊥ /∆Ep.TABLE I. Results from Bethe Ansatz solution of 1D Hubbard model with attractive interaction at quarter-filling.Here, we compute the TLL parameter of charge sector Kρ and the field theory spin gap ∆σ as a function of interactions.