Understanding repulsively mediated superconductivity of correlated electrons via massively parallel DMRG

TL;DR

This paper demonstrates that the 2D U-V Hubbard model, a minimal model for unconventional superconductivity, exhibits signs of superconducting order through large-scale parallelized DMRG calculations, ruling out other phases.

Contribution

It provides the first nonperturbative evidence supporting superconductivity in the 2D U-V Hubbard model using massively parallel DMRG, with systematic size and error control.

Findings

Spin gap compatible with fully gapped or point-zero spectrum

Exclusion of Fermi liquid and magnetic order as ground states

Enhanced short-range correlations in the d_{xy} pairing channel

Abstract

The so-called minimal models of unconventional superconductivity are lattice models of interacting electrons derived from materials in which electron pairing arises from purely repulsive interactions. Showing unambiguously that a minimal model actually can have a superconducting ground state remains a challenge at nonperturbative interactions. We make a significant step in this direction by computing ground states of the 2D \mbox{U-V} Hubbard model - the minimal model of the quasi-1D superconductors - by parallelized DMRG, which allows for systematic control of any bias and that is sign-problem-free. Using distributed-memory supercomputers and leveraging the advantages of the \mbox{U-V} model, we can treat unprecedented sizes of 2D strips and extrapolate their spin gap both to zero approximation error and the thermodynamic limit. Our results for the spin gap are shown to be compatible with a spin excitation spectrum that is either fully gapped or has zeros only in discrete points, and conversely that a Fermi liquid or magnetically ordered ground state is incompatible with them. Coupled with the enhancement to short-range correlations that we find exclusively in the $d_{xy}$ pairing-channel, this allows us to build an indirect case for the ground state of this model having superconducting order in the full 2D limit, and ruling out the other main possible phases, magnetic orders and Fermi liquids.

Tables

Table 1. Table 1: (a) Summary of spin gaps at L → ∞ → 𝐿 L\rightarrow\infty for different N ch subscript 𝑁 ch N_{\rm ch} , assuming the 1 / L 1 𝐿 1/L -scaling of Δ ¯ s ​ [ L , N ] subscript ¯ Δ 𝑠 𝐿 𝑁 \overline{\Delta}_{s}[L,N] we appear to observe continues to hold beyond the strip lengths we computed. “Extr. + estimates” is based on estimating ground state energies at zero truncation error ( ϵ = 0 italic-ϵ 0 \epsilon=0 ) for cases where ground states cannot be systematically extrapolated to ϵ = 0 italic-ϵ 0 \epsilon=0 , because converged simulations were only available for a single χ 𝜒 \chi -value (see text for details). (a) : only one ground state energy of the finite- L 𝐿 L ensemble used to obtain this value via linear fit was estimated. (b) : multiple finite- L 𝐿 L ground state energies were estimated. (c) : All finite- L 𝐿 L ground state energies were estimated.

	${\overline{\Delta }}_s\left[N_{\mathrm{ch}}\right]$ [$t\times {10}^{-3}$]
	Straight extrapolation		Extr. + estimates
$N_{\mathrm{ch}}$	$V/t=0.5$	$V/t=1$	$V/t=0.5$	$V/t=1$
$2$	$5.59$	$7.59$	$5.36$	$7.31$
$4$	$9.68$	$10.2$	${}_{}{}^{\left(a\right)}10.77\left(10\right)$	$10.268\left(11\right)$
$6$	$9.37$	$11.4$	${}_{}{}^{\left(b\right)}4.3\left(1.4\right)$	$11.20\left(97\right)$
$8$	$9.45$		${}_{}{}^{\left(c\right)}7.6\left(1.3\right)$

Table 2. Table 2: Decay of correlations at short range, r ≈ r s = 4 ​ a 𝑟 subscript 𝑟 𝑠 4 𝑎 r\approx r_{s}=4a , in the effective 2D regime (c.f. main text and example shown in Fig,  8 a). The strong enhancement in the d x ​ y subscript 𝑑 𝑥 𝑦 d_{xy} -channel compared to all other channels for N ch ≥ 4 subscript 𝑁 ch 4 N_{\rm ch}\geq 4 is highlighted in red.

	$N_{\mathrm{ch}}=2$		$N_{\mathrm{ch}}=4$		$N_{\mathrm{ch}}=6$		$N_{\mathrm{ch}}=8$
short-range decay	$V/t=0.5$	$V/t=1$	$V/t=0.5$	$V/t=1$	$V/t=0.5$	$V/t=1$	$V/t=0.5$
$\frac{{\max }_{\left|r-r_s\right|\le 1}\left[d_{xy}\left(r\right)\right]}{d_{xy}\left(1\right)}$	$0.17$	$0.15$	$1.08$	$0.55$	$0.67$	$0.79$	$1.05$
$\frac{{\max }_{\left|r-r_s\right|\le 1}\left[d_{x^2-y^2}\left(r\right)\right]}{d_{x^2-y^2}\left(1\right)}$	$0.09$	$0.08$	$0.08$	$0.09$	$0.09$	$0.07$	$0.09$
$\frac{{\max }_{\left|r-r_s\right|\le 1}\left[s\left(r\right)\right]}{s\left(1\right)}$	$0.97$	$0.87$	$0.22$	$0.09$	$0.10$	$0.16$	$0.14$
$\frac{{\max }_{\left|r-r_s\right|\le 1}\left[C\left(r\right)\right]}{C\left(1\right)}$	$0.06$	$0.07$	$0.06$	$0.06$	$0.06$	$0.06$	$0.06$

Table 3. Table 3: Spin gaps Δ ¯ s ​ [ L , N ch ] subscript ¯ Δ 𝑠 𝐿 subscript 𝑁 ch \overline{\Delta}_{s}[L,N_{\rm ch}] (in units of t × 10 − 3 𝑡 superscript 10 3 t\times 10^{-3} ) upon which extrapolations L → ∞ → 𝐿 L\rightarrow\infty were performed in Figs.  5 and for Fig.  7 .

	$V/t=0.5$				$V/t=1$
$L$	$N_{\mathrm{ch}}=2$	$N_{\mathrm{ch}}=4$	$N_{\mathrm{ch}}=6$	$N_{\mathrm{ch}}=8$	$N_{\mathrm{ch}}=2$	$N_{\mathrm{ch}}=4$	$N_{\mathrm{ch}}=6$
$20$	$95.953$	$96.945$	$96.956$	$97.023$	$87.977$	$89.311$	$89.516$
$32$	$66.272$	$65.597$	$64.519$	$64.799$	$61.726$	$60.503$	$61.011$
$40$	$53.981$	$53.943$	$54.949$	$52.908$	$51.298$	$50.607$	$51.774$
$48$	$45.559$	$46.421$	$45.634$		$43.520$	$45.186$	$42.838$
$56$	$39.519$	$41.736$	$39.809$		$37.810$	$37.049$
$60$	$37.113$				$35.515$
$64$	$35.024$	$35.634$			$33.508$	$34.508$
$72$	$31.646$				$30.193$
$80$	$29.213$				$27.619$
$96$	$24.326$				$24.061$
$120$	$19.618$				$19.231$
$160$	$14.991$

Equations

$$\hat{H}$$

$$-t_{\perp }\sum _{i=1}^L\sum _{n=1}^{N_{\mathrm{ch}}-1}\sum _{\sigma =\uparrow ,\downarrow }\left({\hat{c}}_{i,n,\sigma }^{†}{\hat{c}}_{i,n+1,\sigma }+\text{mbox}h.c\right)$$

$$+U\sum _{i=1}^L\sum _{n=1}^{N_{\mathrm{ch}}}{\hat{n}}_{i,n,\uparrow }{\hat{n}}_{i,n,\downarrow }+V\sum _{i=1}^{L-1}\sum _{n=1}^{N_{\mathrm{ch}}}{\hat{n}}_{i,n}{\hat{n}}_{i+1,n}$$

$${\chi }_s^{-1}(L,N_{\mathrm{ch}})=L\frac{E_{\mathrm{GS}}^{L,N_{\mathrm{ch}}}(S_z=s)+E_{\mathrm{GS}}^{L,N_{\mathrm{ch}}}(S_z=-s)-2E_{\mathrm{GS}}^{L,N_{\mathrm{ch}}}(S_z=0)}{(2s{)}^2}$$

$$A(L,N_{\mathrm{ch}},s)$$

$${\overline{\Delta }}_s[L,N_{\mathrm{ch}}]:=\frac{E_{\mathrm{GS}}^{L,N_{\mathrm{ch}}}(S_z=N_{\mathrm{ch}})-E_{\mathrm{GS}}^{L,N_{\mathrm{ch}}}(S_z=0)}{N_{\mathrm{ch}}},$$

$${\Delta }_s[L,N_{\mathrm{ch}}]:=E_{\mathrm{GS}}^{L,N_{\mathrm{ch}}}(S_z=1)-E_{\mathrm{GS}}^{L,N_{\mathrm{ch}}}(S_z=0),$$

$$d_{x^2-y^2}(r)$$

$$d_{xy}(r),s(r)$$

$$C(r)$$

$$\hat{H}=\sum _{b_l,b_c,b_r=1}^D{\hat{H}}_{\mathrm{left}}[b_l]\otimes {\hat{h}}_i[b_l,b_c]\otimes {\hat{h}}_{i+1}[b_c,b_r]\otimes {\hat{H}}_{\mathrm{right}}[b_r]$$

Full text

I Introduction

Understanding the impact of inter-particle repulsions in fermionic quantum systems remains one of the most challenging problems in physics. As was shown by Landau, electron-electron repulsions in a solid may often be treated by redefining excitations as quasiparticles. The resulting physics is thus essentially of free particles (the so-called Landau quasiparticles), with interactions mostly leading to a renormalization of some physical parameters such as the mass. This Fermi liquid (FL) approach has seen very broad success and allowed to understand the effects caused by perturbations on top of the Fermi liquid, such as instabilities towards e.g. a Bardeen-Cooper-Schrieffer (BCS) superconducting (SC) state, a magnetically ordered (MO) state or a charge density wave (CDW) state Ziman (1972).

There are however situations where the FL approach breaks down, leading to so-called “non-Fermi liquid” behavior. At high spatial dimensionality, this can happen if the interactions are particularly strong, or if the filling of the systems is commensurate with the lattice, with e.g. one particle per site leading to a Mott (MI) state. Reduced dimensionality of the system can further enhance the effect of interactions. In one dimension in particular, the interactions always have drastic effects and lead to physical properties very different from the ones of a Fermi liquid. The FL is replaced by another set of universal features, called the Tomonaga-Luttinger liquid (TLL), where the elementary excitations are the collective excitations of charge and spin in the system. Such systems are critical and at zero temperature $T=0$ possess various competing quasi-long range orders, ranging from antiferromagnetism to superconductivity Giamarchi (2003).

Non-FLs have thus been prime candidates to search for novel physical phases, in particular unconventionally superconducting (USC) phases mediated by excitations other than the electron-phonon one, of which the high-temperature (high-$T_{c}$) superconductors Anderson (1997); Orenstein (2000). are just one example. The observed proximity between USC and MO phases has made fluctuations around the MO state prominent candidates for such pairing mechanisms. On the experimental side, many systems ranging from heavy fermions Varma (1985); Si and Steglich (2010), organic superconductors Bourbonnais and Jérome (2007); Jérome (2012), high-$T_{c}$ superconductors and pnictides Damascelli and Shen (2003); Fischer et al. (2007); Scalapino (2012) show USC phases. Yet, obtaining any widely accepted theory on the nature of such exotic pairing for any of these material-groups has proven to be difficult; even the simplest physical models, the so-called minimal models, which are abstracted from the materials’ full physical structure are fundamentally hard to solve: the 2D U-V model at half filling (organics), and the weakly doped 2D Hubbard model (cuprates). So far, any solution of these minimal models has almost always necessitated additional technical approximations that introduce errors of unknown magnitude. It is thus difficult to determine whether phases predicted for such models in the framework of a given approximation scheme are truly present or an artefact of the approximation, or whether to explain the experiments, different, more extended models are required, such as e.g. in the three-band models of the high-$T_{c}$ cuprates Varma C. M. and Giamarchi (1991), or in those approaches that investigate the role of phonons in USC materials Salje et al. (1996).

Even when moving to the more tractable regime of weak repulsion, SC and MO instabilities interfere at all orders of the perturbative renormalization group (RG) treatments applied so far Schulz (1987); Furukawa et al. (1998); Läuchli et al. (2004), necessitating complex numerical methods Deng et al. (2015) or sophisticated, and still approximate, functional RG procedures Nickel et al. (2005, 2006); Bourbonnais and Sedeki (2011); Sedeki et al. (2012); Metzner et al. (2012) . More uncertain yet is the situation in the most important regime, actually relevant to the models’ underlying materials: when repulsion becomes large and competes with kinetic energy. For this regime, numerical methods are a route of choice. However, in two or more spatial dimensions these too are faced with steep challenges. Quantum Monte Carlo (QMC) treatments suffer from the so-called sign problem for fermions leading to an exponential degradation of the signal to noise ratio as e.g. the temperature is lowered LeBlanc et al. (2015). And while innovative techniques like diagrammatic QMC have delivered intriguing insights into USC pairing for the 2D Hubbard model, these are still constrained to intermediate interactions at most and doping very far from unity Deng et al. (2015).

Our numerical approach here is different, and rests on two observations: (1) The fact that for the Density Matrix Renormalization Group (DMRG) technique, and other tensor-network methods such as PEPS, any possible bias can be controlled against via extrapolation, and that these methods are free of the sign problem Schollwöck (2011); Stoudenmire and White (2012). (2) That the one minimal 2D model in which a USC ground state from strong electron repulsions can unambiguously be shown to exist, is comprised of quasi-1D electrons (doped 2-leg Hubbard ladders) weakly coupled in parallel Carlson et al. (2000); Karakonstantakis et al. (2011). This minimal model however, as yet, does not correspond to any extant material.

These two observations have led us to conclude that the minimal 2D U-V model of the organic superconductors - which has an analogous structure, of many 1D systems coupled weakly in parallel to each other - offers an unique opportunity for insight on the strongly interacting ground state of a candidate model for USC pairing. This would be done by swapping out the perturbative theory used on the Hubbard-ladders in the low-energy field-theory limit, which would not be applicable to the U-V model setting, with quantitatively reliable DMRG.

Yet, using DMRG in this particular setting comes with issues to resolve first. While DMRG delivers highly accurate results both for the statics and the dynamical correlations in 1D systems while using only modest computational resources, the 2D setting of the U-V model is different. On such lattices, DMRG is known tho require resources that increase quasi-exponentially with lattice-width when trying to maintain any pre-set accuracy Stoudenmire and White (2012). Tractable sizes for systems of itinerant electrons so far usually have been on the order of about 100 sites, depending on the specific problem White and Scalapino (2009); Scalapino and White (2012); Liu et al. (2012); LeBlanc et al. (2015); Ehlers et al. (2017); Zheng et al. (2017) - pushing towards around 200 sites, while done, is generally challenging and goes along with a decrease in accuracy. As the sum of the evidence indicates that states with SC order are in close energetic competition with MO ones, even changes in geometry (i.e. boundary effects) could easily result in different orders winning out for lattices of such size. As a result, no unbiased approach has so far been able to show unambiguously that a ground state of either of the two minimal models may exhibit USC order when repulsion is non-perturbatively large.

In the present work, we take a significant step in that direction. We build on our development of the numerical parallelized density matrix renormalization group (pDMRG), a distributed-memory version of one of the standard DMRG formulations Schollwöck (2011) capable of exploiting modern supercomputer architectures. With this method we investigate the properties of the 2D U-V Hamiltonian at half filling of the bands. We consider an array of parallel chains of strongly correlated 1D electrons, each described by such a U-V Hamiltonian and weakly coupled to each other by transverse tunnelling. In addition to being the model Hamiltonian for organic superconductors Giamarchi (2004), and thus potentially containing the unconventionally superconducting phase observed in this system, this Hamiltonian has two critical advantages for numerical study compared to the doped 2D Hubbard model LeBlanc et al. (2015): (1) it is trivial to maintain fixed ratios of electrons to lattice sites and thus to extrapolate to infinite system size here, since no doping is required. The pDMRG then allows doing this with the accuracy required to extrapolate energies of large systems, such as to enable the reliable calculation of energy differences, as appearing in e.g. the most important observable we use, the spin gap. Underlying this accuracy is that pDMRG can exploit the strong anisotropy of electron tunneling in the 2D U-V model by spreading out a single ground state calculation across many supercomputer nodes. This anisotropy ties directly into the U-V models’ second advantage: (2) the strip or cylinder geometry, that DMRG is especially good at handling intrinsically, suits the U-V model naturally if aligned with the strong-tunneling direction, as correlations across the strips’ width will be naturally weak due to the small perpendicular tunneling, which is not the case for the doped 2D Hubbard model.

We show that in the thermodynamic limit this systems’ averaged spin gap $\overline{\Delta}_{s}$ is incompatible with either a FL or an MO ground state, for strips of finite width. The behaviour of $\overline{\Delta}_{s}$ as $L\rightarrow\infty$ appears most compatible with $\overline{\Delta}_{s}>0$ at infinite length, and being either non-decaying or even increasing with the strips’ width, i.e. pointing towards a fully gapped spin excitation spectrum. But we also consider the alternate possibility that $\overline{\Delta}_{s}$ might ultimately scale to zero beyond the system sizes we can access. However, we show this alternative is still at most compatible with a spin excitation spectrum with isolated zero-energy points, i.e. in line with weak-coupling theories of USC systems. But either scenario is incompatible with either antiferromagnetic order or Fermi liquid-like ground state. Together with the strong enhancement exclusive to the $d_{xy}$-channel of correlations that we find with increasing strip width, our results make the case for SC singlet pairing with $d_{xy}$-order being the dominant component in the ground state of the U-V model at half filling when repulsion is competitive with kinetic energy.

The structure of the paper is as follows: in section II we outline the U-V model and our approach to analysing the system and choosing the model parameters; in section III we describe the key features of pDMRG relevant to this work; in section IV we detail our scaling procedure; in section V we analyse the behaviour of the spin gap and the correlation functions as the width of the 2D U-V strips is increased, making the case for a SC ground state with $d_{xy}$-order and against any MO or FL competitor-order; in section VI we discuss these results and how they connect to experiment; in Appendix A we provide details for the used scaling procedure; finally, in Appendix B we provide a broader technical overview on pDMRG.

II Model and physical observables

As depicted in Fig. 1a, the U-V model in 2D is formed by a parallel array of $N_{\rm ch}$ 1D chains of lattice electrons, each of length $L$, with inter-chain tunneling. Its lattice-Hamiltonian is

[TABLE]

Here, $\hat{c}^{\dagger}_{i,n,\sigma}$ ($\hat{c}_{i,n,\sigma}$) denotes the electron creation (annihilation) operator on site $i$ of chain $n$ of spin $\sigma=\uparrow,\downarrow$, and $\hat{n}_{i,n,\sigma}:=\hat{c}^{\dagger}_{i,n,\sigma}\hat{c}_{i,n,\sigma}$, $\hat{n}_{i,n}:=\hat{n}_{i,n,\uparrow}+\hat{n}_{i,n,\downarrow}$ denote local electron density for spin $\sigma$ and total electron density operators respectively. The tunneling amplitude along the chains is given by $t$, while in-between adjacent chains it is given by $t_{\perp}$, and generally $t\gg t_{\perp}$. We point out that the $t_{\perp}$-term above obeys open boundary conditions, just as the $t$-term does. We have made this choice because obtaining ground state energies accurate enough to compute reliable differences between them is an indispensable requirement to extrapolate this systems spin gaps in the limit $L\rightarrow\infty$. Cylindrical 2D-lattices however are well known to immediately suffer from lower accuracy, and are thus avoided here. Finally, the onsite Coulomb repulsion is given by $U$, while intrachain nearest-neighbour repulsion is $V$, with $V<U/2$. The entire system is at half-filling, i.e. $N_{\uparrow}=N_{\downarrow}=LN_{\rm ch}/4$, and thus $k_{F}^{\uparrow}=k_{F}^{\downarrow}=k_{F}=\pi/4a$, where $a$ is the lattice spacing.

This minimal model originally arose from the study of the organic Bechgaard and Fabre salts (‘the organics’), the first materials discovered to support an USC phase Jérome et al. (1989); Bourbonnais and Jérome (2007). These compounds show a very rich phase diagram, and especially a change in effective dimensionality with temperature - 1D TLL-like at high temperature, then 2D-like and finally 3D-like at near zero - as quantum coherence can increasingly be established along the three orthogonal and successively weaker directions for electron-tunneling. Their most striking feature is the phase transition from a magnetically ordered to a unconventionally superconducting phase at low temperature, as the proximity, and thus the inter-chain tunneling $t_{\perp}$, between the strong-coupling chains of the U-V model is increased 111Experimentally, this distance is controlled either by applying external pressure to the sample, or chemically (choice of anions between the cationic stacks of molecules that make up the 1D chains). To understand this competition, and whether the minimal U-V model (II) of the organics can actually capture it at least in a weak-interaction scenario (which does not correspond to the actual materials), perturbative RG has been repeatedly applied Nickel et al. (2005, 2006); Bourbonnais and Sedeki (2011); Sedeki et al. (2012). While beset with the same technical limitations this approach has for the doped 2D Hubbard-model of the cuprates (c.f. Sec. I), it does suggest that the minimal 2D U-V model, might indeed support a USC phase transition from a MO phase when extended with e.g. a next-to-nearest neighbour tunneling perpendicular to the chains, and that pairing would be in the $d_{x^{2}-y^{2}}$-channel (gapless superconductivity would be compatible with experiments on the organics Shinagawa et al. (2007); Yonezawa et al. (2012)). However, so far there has been no complementary method to either validate these results, and especially no numerics that could deal with the U-V model (II) for the experimentally relevant case of strong interactions.

In order to be able to address the case $t_{\perp}\neq 0$, $N_{\rm ch}>2$ for this Hamiltonian with reliable numerics, we have developed pDMRG to exploit the distributed-memory architectures of modern supercomputers. Based on the tensor-network formulation widely used in modern DMRG, the code has major advantages over other quantitative numerical approaches for this problem: it inherits the lack of sign-problem and the ability to use extrapolation in the known error to filter out possible implicit bias towards any particular order (which is inherent e.g. to any wavefunction ansatz of variational quantum Monte Carlo) from conventional DMRG. And further, as described in Sec. III and Appendix B, the parallelized nature of pDMRG makes the necessary ground state calculations tractable in the first place, by distributing the many non-local tunneling terms of the model across different nodes of the supercomputer. As laid out in Sec. III, the large number of such terms becomes unavoidable when keeping the amount of bipartite entanglement at manageable levels, which is performance-critical for any 2D model.

As detailed in Sec. II.1, the spin gap is the most important observable for our study of this system, as its scaling as $L\rightarrow\infty$ allows the crucial elimination of possibilities for ground state ordering in the thermodynamic limit that has so far proven elusive for any minimal model of USC systems. This is because ground states with either MO order or in a conventional FL state will show a linear vanishing of any spin-gap. In contrast, standard analytical theories of USC systems combine weak-coupling RG, for identifying the leading instability and its symmetry, with mean-field descriptions of the ordered phase. For these theories, the unavoidable zero-nodes of the order parameter result in sublinear scaling to zero for a spin gap Sigrist and Ueda (1991); Hasegawa (1996); Moriya and Ueda (2003); Shinagawa et al. (2007); Chubukov et al. (2008); Vladimirov et al. (2011), while experiments on the actually strongly coupled regime in USC materials have found finite spin gaps (c.f. Sec. VI). It is therefore vital that the known maximal value of the approximation error of the DMRG technique allows us to extrapolate the computed spin gaps to zero truncation error for any given lattice size, as described in IV.2.

To outline our strategy of analysis, we now discuss in Sec. II.1 the concrete ground state observables we considered, providing the roadmap for ruling out the main competitor states to any superconducting ground states, the FL and the MO states. In Sec. II.2 we then discuss our choice of parameters for Hamiltonian (II), and why these should reduce the competition with charge-ordered ground states.

II.1 Observables

One of our key priorities is to establish whether any set of parameters for Hamiltonian (II) could or could not result in one of the main competitor orders to a USC ground state in the 2D limit, a MO or FL state, being realized. One quantity that would allow making these distinctions is the spin susceptibility at zero temperature

[TABLE]

where $E_{\rm GS}^{L,N_{\rm ch}}$ denotes the ground state energy of (II) at half-filling, and $S_{z}=N_{\uparrow}-N_{\downarrow}$ is the difference between spin-populations. Here, $s$ can be any spin number, as long as it is negligible compared to a macroscopic number of spins (proportional to $L$). Typically, $s$ is chosen so that other quantum numbers in the system can be kept identical between $s$ and $s=0$. The usual choice is $s=1$, corresponding to the spin flip of a single spin. However, as discussed in subsection IV.1, the choice $s=N_{\rm ch}$ has a crucial technical advantage for extrapolating to infinite length $L\rightarrow\infty$.

When $L\to\infty$, $\chi_{s}(L,N_{\rm ch})$ becomes the thermodynamic spin susceptibility. In order to get a finite susceptibility, as expected for a FL or a MO state, one would thus need the quantity

[TABLE]

to scale as $1/L$. Other scalings of this quantity thus directly signal a non-FL, non-MO behavior for the spin susceptibility. The simplest case is that the averaged spin gap

[TABLE]

remains finite when $L\to\infty$, where $\overline{\Delta}_{s}[L,N_{\rm ch}]$ is the generalization of the standard spin gap

[TABLE]

Then the spin susceptibility would scale as ${(L\overline{\Delta}_{s}[L\rightarrow\infty,N_{\rm ch}])^{-1}}$, leading to zero susceptibility in the thermodynamic limit, as expected for fully gapped systems. For SC systems with unconventional electron pairing such fully gapped spin-excitations have been actually been found experimentally (see Section VI).

More complex behavior can be expected if (3) has a different scaling with the size of the system such as e.g. $1/L^{\alpha}$ with $0<\alpha<1$. The thermodynamic spin susceptibility remains zero in this case, yet the systems’ spin excitations are not fully gapped. In mean-field theories of USC systems all this stems from the line- or point-nodes of the SC gap Sigrist and Ueda (1991); Hasegawa (1996); Moriya and Ueda (2003); Shinagawa et al. (2007); Chubukov et al. (2008); Vladimirov et al. (2011). At finite temperatures, the scaling with $L$ would be converted into a non-trivial temperature dependence of the spin susceptibility corresponding to a pseudogap behavior James et al. (2013); Chen et al. (2017).

Establishing the $L$-dependence of $\overline{\Delta}_{s}[L\rightarrow\infty,N_{\rm ch}]$ through our numerical results therefore allows to differentiate between three possibilities for the ground state of the two dimensional system in the $T=0$ limit: FL/MO, or fully spin gapped USC, or mostly spin-gapped USC - to the extent that we can extrapolate trends as $N_{\rm ch}$ is increased.

The use of pDMRG allows meaningful extrapolation of $\overline{\Delta}_{s}[L\rightarrow\infty,N_{\rm ch}]$ in $1/L$ for the first time, enabling the calculation of ground states of multiple coupled, long U-V chains ($N_{\rm ch}\leq 8$, $L\leq 64$) with quantitatively accurate numerics (c.f. Sec. III) 222 $N_{\rm ch}$ larger than presented in this work can be handled by pDMRG. As a single calculation can be spread out across dozens of nodes i.e. thousands of CPU cores, a given problem is less limited by the capacity of the computer as such, but by the available time on the supercomputer.. As discussed in Sec. V.1, this allows us to make a strong case that for our choice of Hamiltonian parameters any magnetic order or Fermi liquid behaviour is suppressed in the ground state of the 2D U-V model, due to the way $\overline{\Delta}_{s}[L,N_{\rm ch}]$ behaves with $L\rightarrow\infty$, and showing this beviour to be robust as $N_{\rm ch}$ grows.

Our point of reference for this are 2D arrays of coupled spin-chains, which demonstrate the opposite scenario, of the spin gap decaying with $N_{\rm ch}$. If infinitely many 1D Heisenberg-chains were to be coupled in parallel, the result is a gapless 2D Heisenberg model - but when only two such chains are coupled this system exhibits a finite spin gap Haldane (1983); Dagotto and Rice (1996). Multiple techniques have since shown how these two extremes are bridged: as the number of coupled Heisenberg chains increases, the spin gap remains finite (for an even numbers of chains) but continuously decreases as $N_{\rm ch}$ grows, scaling to zero in the limit of infinitely many chains Schulz (1986); Affleck et al. (1989).

The information gained from spin gaps can be supplemented by studying the change in the short-range behaviour of certain correlation functions as $N_{\rm ch}$ is increased. When magnetic orders can be ruled out, the functions of interest here are those of the two other main competitor groups: SC-ordering on the one hand - with the possibility of pairing in the $d_{x^{2}-y^{2}}$-, $d_{xy}$- and extended $s$-wave channels - and CDW-ordering on the other hand. The corresponding correlation functions considered for this are

[TABLE]

where $\hat{D}$ denotes the nearest-neighbour spin-singlet operators $\hat{D}_{i,n,j,m}:=\hat{c}_{i,n,\uparrow}\hat{c}_{i+j,n+m,\downarrow}-\hat{c}_{i,n,\downarrow}\hat{c}_{i+j,n+m,\uparrow}$. As illustrated in Fig. 8b, we consider the SC correlation functions in the two central chains $N_{\rm ch}/2$, $N_{\rm ch}/2+1$, and $C(r)$ on chain $N_{\rm ch}/2$, in each case at distance $r$ from the middle of the chain(s). Based on these correlation functions, we present further evidence in Sec. V.3 beyond the behaviour of the spin gap that for the U-V model we study here, electrons do in fact pair, with a dominant component in the $d_{xy}$-channel.

II.2 Choosing model parameters - $U=4t$, $V/t=0.5,1$

The basic results serving as background to any choice of Hamiltonian parameters are summarized in Fig. 1b. In the limit $t_{\perp}=0$, the chains decouple and are described by TLL-theory for sufficiently small $V$. As long as $U$ is fixed to any value $\geq 2t$, there exists a value $V_{c}$ above which a charge gap opens in the system, turning it into a Mott insulator (MI) Mila and Zotos (1993); Ejima et al. (2005). The phase boundary at $t_{\perp}=0$ is shown as a green line in Fig. 2. Now, in order to maximise the chance of finding SC order, we attempt to target a parameter regime in which the decoupled chains would be in the TLL phase at the same time as the coupled chains ($t_{\perp}>0$) would remain in a delocalized regime. Ultimately, only the computational results can indicate whether one succeeds with the latter condition. But given the limited supercomputing time at hand, we try to maximise our chances for this to happen from the outset, as described within the rest of this section.

A look at the minimal U-V model of the organic superconductors (Bechgaard and Fabre salts) is instructive here. Experimental probes of these show that here the 1D chains of the U-V model have a Luttinger-liquid parameter $K_{\rho}$ (which characterizes interactions and correlations of the charge mode of an isolated chain) in the vicinity of $0.25$ Giamarchi (2004). When $K_{\rho}=0.25$, theory predicts the single chain at half-filling to develop a charge gap and become a MI. In terms of the microscopic parameters, this vicinity to $K_{\rho}=0.25$ for the single chain would correspond to $U/t=10$, $V/t\approx 2-3$, c.f. yellow area in Fig. 2. When the ratio $t_{\perp}/t$ is small enough, experiments show the physics of the single chain extending to the whole material, i.e. the system is an MI as well, and magnetically ordered on top of it. Only as $t_{\perp}/t$ grows and $t_{\perp}$ becomes competitive with the Mott gap of the single chain is superconductivity found in experiments, typically at $t_{\perp}/t\approx 0.1$.

It may be tempting to try to straightforwardly replicating this setting for our calculations. However, there has recently been new insight on the two-chain limit by two co-authors of the present work Kantian and Giamarchi , which cautions against such an approach. Combining analytical RG treatment of the bosonized chains with DMRG, we find that finite $t_{\perp}$ for two chains increases the U-V parameter space in which the system develops a Mott gap (red line in Fig. 2). While early exploratory work indicates that this parameter space might shrink again for $N_{\rm ch}=4$, it is beyond the scope of this work or of Ref. Kantian and Giamarchi, for definite statements how exactly the few-chain regime evolves in this regard.

In the following, we therefore focus on a regime in which these complications of few-chain physics are unlikely to apply, at $U/t=4$ and $V/t=0.5,1$, and with $t_{\perp}/t=0.1$ (blue shaded region in Fig. 2). This regime realizes the essential physics present in the superconductivity of the organics, i.e. a regime where inter-chain tunneling is not effectively suppressed by a single-chain gap. A priori, this regime thus appears to offer the highest chance of finding a superconducting ground state due to repulsively mediated pairing. We have targeted two values of $V/t$ to be certain our results will be representative for a finite area of parameter space.

III The Method: parallel DMRG

As stated in Sec. II, exploiting the physical structure of the 2D U-V-model is essential for being able to handle the bipartite entanglement of this system when mapping the physical lattice to the DMRG-chain. Being comprised of many strongly correlated 1D chains in parallel with weak coupling in-between, it would seem intuitive that performing the mapping along the direction of strong tunneling (c.f. Fig. 3a) delivers lower overall bipartite entanglement in the 1D DMRG chain. This argument is complemented by theoretically estimating the performance of the alternate mapping along the $t_{\perp}$-direction that is used conventionally when applying DMRG to what is effectively a long and narrow 2D-strip. The bond dimension, denoted as $\chi$ in the following, is the central quantity that controls the accuracy of DMRG. It denotes the number of optimally chosen basis states in which any bipartition of the system is represented. If $\chi_{\rm chain}$ was required to represent a single U-V-chain of a given length with truncation error $\epsilon$, achieving the same accuracy for $N_{\rm ch}$ chains would require $\chi=\chi_{\rm chain}^{N_{\rm ch}}$ - even in case $\chi_{\rm chain}=10$, this scaling would be hopeless beyond $N_{\rm ch}=2$. A simple test comparing the two mappings for a representative system of free fermions, summarized in Fig. 3b, confirms the total disparity between the performance of these mappings. Shown is the energy vs. the number of DMRG iterations for $L=20$, $N_{\rm ch}=5$, $N_{\uparrow}=N_{\downarrow}=50$, $t_{\perp}/t=0.2$, $U=V=0$, $\chi=5000$, once with mapping along $t$ (blue line), once with mapping along $t_{\perp}$(green line). The $t$-mapping converges quickly to the exact ED solution, while the $t_{\perp}$-mapping flounders even for this small system.

The price paid for the $t$-direction mapping, necessary to have manageable $\chi$ in the first place, is the very large number of long-range tunneling terms. In the two-site DMRG we use, for every adjacent pair of sites ${(i,i+1)}$ on the DMRG-chain for which the tensors are jointly optimized, the Hamiltonian is represented as a matrix product operator (MPO) Schollwöck (2011):

[TABLE]

Here, each sum runs from $1$ to $D$, the MPO bond-dimension, which for the $t$-direction mapping scales as $D\approx 4L$. For the range of system lengths we treat, $L=20$ to $64$, this amounts to $D\approx 80$ to $256$, implying an equal number of Hamiltonian contributions $\hat{H}_{\rm left}[b_{l}]$ ($\hat{H}_{\rm right}[b_{r}]$). These describe the action of the Hamiltonian to the left (right) of the sites $i$, $i+1$, while for every $b_{l}$, $b_{c}$ ($b_{c}$, $b_{r}$) $\hat{h}_{i}[b_{l},b_{c}]$ ($\hat{h}_{i+1}[b_{c},b_{r}]$) is a purely local operator on site $i$ ($i+1$). Combined with the fact that we utilize $\chi=10000$ to $18000$, this large number of Hamiltonian contributions result in memory requirements on the order of several TB for any single ground state calculation and a commensurate amount of computational effort.

Parallelized DMRG was developed to make calculations of such magnitude tractable Bauer et al. . We provide a more comprehensive overview of its technical features in Appendix B, and highlight here its first parallelization layer, which is the most relevant to this work, shown schematically in Fig. 3c: the largest objects of DMRG are the boundary terms. These are the expressions of $\hat{H}_{\rm left}[b_{l}]$, $\hat{H}_{\rm right}[b_{r}]$ in the $\chi$ optimal basis states of the system to the left or right of sites $i$, $i+1$ respectively. Denoting these basis states as $|\alpha_{l}\rangle$, $|\alpha_{r}\rangle$ respectively, the left boundary is $L[b_{l}]:=\sum_{\alpha_{l},\alpha_{l}^{\prime}}\left(\langle\alpha_{l}^{\prime}|\hat{H}_{\rm left}[b_{l}]|\alpha_{l}\rangle\right)|\alpha_{l}^{\prime}\rangle\langle\alpha_{l}|$, and the right boundary $R$ is defined analogously using $\hat{H}_{\rm right}[b_{r}]$ and $|\alpha_{r}\rangle$. The parallelization layer distributes each of the $D$ elements of both these vectors (each element being a block-sparse matrix due to the use of conserved quantum numbers to enhance performance), to nodes of a parallel supercomputer. For this project, the computer cluster was the ’Piz Daint’ system of the Swiss National Supercomputing Center (CSCS). In this manner, we have spread the calculations out to up to many dozens of nodes (for the largest lattices).

IV Scaling analysis

As discussed in Sec. II, our indirect strategy for arguing that a ground state of the 2D U-V model has SC order, by eliminating competitor phases, relies on obtaining the correct behaviour of the systems spin gap as $L\rightarrow\infty$. In the following we discuss the two necessary steps required in service of this goal: sharply reducing the oscillations of the spin gap as $1/L$ decreases, and extrapolating ground states energy to zero error in the DMRG-approximation.

IV.1 Removing oscillations of $\Delta_{s}[L]$

The standard definition of the spin gap, eq. (5), is based on the difference between two ground state energies. However, due to the small value of $t_{\perp}/t$ in the regime of interest, one encounters oscillations of $\Delta_{s}[L]$ with $L$, which presents a complication for obtaining clean extrapolation of $\Delta_{s}[L]$ in $L$, as shown in Fig. 4a. That weak-interchain coupling is evidently the cause is demonstrated by studying the U-V model in the non-interacting limit $U=V=0$. While clearly $\Delta_{s}=0$ for $L\rightarrow\infty$ (at any $N_{\rm ch}$) in this case, the $N_{\rm ch}$ bands are only weakly split. Then, the shifting positions of both the highest occupied as well as of the lowest unoccupied states in each band with changing $L$ (c.f. Fig. 4b) will invariably cause such oscillatory behaviour in the definition (5).

Based on this analysis, clearly one should smooth this finite-size effect by using a more general definition of the spin gap. The natural generalization, given by eq. (4), is the one that averages over all the $N_{\rm ch}$ bands close to the Fermi-level in the non-interacting limit, and that applies unchanged to the interacting regime. As depicted in Fig. 4a and b, we find this definition removes the oscillatory behaviour to a large degree for all our data. As discussed in Appendix A, this definition of the averaged spin gap results in a linear polynomial in $1/L$ as the only consistent method for extrapolating $\overline{\Delta}_{s}[L,N_{\rm ch}]$ to $L^{-1}=0$; linear extrapolation can recover physically expected behaviour correctly, while the oscillatory remnant overwhelms any vestiges of quadratic dependency in $1/L$ that might remain after the averaging.

We further note that this procedure is similar to adding two particles instead of one when computing the compressibility of spinless particles, to avoid the unwanted oscillations provoked by the breaking of $k\to-k$ symmetry that adding a single particle would entail. Analogously, one usually adds four particles to stay in the sector of total spin zero for spinful systems.

IV.2 Extrapolating in the DMRG truncation error

In DMRG, one attempts to compute the ground state wavefunction $|\psi_{\rm GS}\rangle$ of a lattice-Hamiltonian, by mapping the physical lattice to a 1D chain. At every bipartition of this DMRG chain, the algorithm achieves the required reduction of the exponentially-scaling full Hilbert space by approximating $|\psi_{\rm GS}\rangle$ through two sets of $\chi$ optimally chosen basis states, yielding a wavefunction $|\psi^{\chi}_{\rm GS}\rangle$. A given $\chi$ corresponds to a particular truncation error $\epsilon:=\langle\delta\psi|\delta\psi\rangle$, ${|\delta\psi\rangle:=|\psi_{\rm GS}\rangle-|\psi^{\chi}_{\rm GS}\rangle}$, where $\epsilon$ is known in principle in DMRG.

The controlled accuracy and known error $\epsilon$ enables local quantities, and thus energies, to be extrapolated to zero error White and Chernyshev (2007). Specifically for energies, it is known that the approximation error in the energy, $\delta E$, is $\propto\epsilon$. As we obtain $\overline{\Delta}_{s}[L,N_{\rm ch}]$ from the difference between the $S_{z}=N_{\rm ch}$ and $S_{z}=0$ ground state energies, wherever we have at least two ground states with different $\chi$ we use this technique to extrapolate the energy to zero $\epsilon$ before computing $\overline{\Delta}_{s}$. In this, we always employ $\epsilon_{\rm max}$, the maximal truncation error committed in the final sweep as a proxy for $\epsilon$. For $V/t=1$, $N_{\rm ch}=4$ we have obtained multiple energies, by computing ground states for $\chi=10000$ and $18000$ independently, then obtained $\chi=12000,14000$ and $16000$ using the $\chi=18000$ ground state as an initial state and applying two complete sweeps at the lower $\chi$. An example is of this is shown in Fig. 4c. As illustrated in Fig. 4d, the change in $\overline{\Delta}_{s}[L\rightarrow\infty,N_{\rm ch}]$ relative to values at finite-$\chi$ is typically small anyway, on the orders of a few percent at most.

Outside these parameters, we have aimed to produce both $\chi=10000$ and $\chi=18000$ ground states in independent simulations wherever possible, but particularly for $V/t=0.5$ available computing ressources proved ultimately insufficient. Thus, many results for this parameter are based partly or fully on a single $\chi=10000$ ground state.

In order to obtain any particular value of $\overline{\Delta}_{s}[L,N_{\rm ch}]$ as the basis for an extrapolation to $L\rightarrow\infty$ at fixed $N_{\rm ch}$, we therefore pursue two separate protocols:

Straight extrapolation - Here, we extrapolate both $S_{z}=0$ and $S_{z}=N_{\rm ch}$ ground state energies to zero $\epsilon_{\rm max}$ when possible. If scaling is only possible in one of the spin sectors, we form $\overline{\Delta}_{s}$ from the two lowest-energy ground states with comparable $\epsilon$.

Extrapolation plus estimates - Whenever a ground state is available at $\chi=10000$ exclusively, we estimate the correction heuristically. We have done this based on the following observations, made on states where extrapolations were possible: (i) the fractions by which energies further decrease upon extrapolation seem to roughly double going from one $L$ to the next, (ii) they also seem to, roughly double with every increase of $N_{\rm ch}$. (iii) for the same parameters, fractions seem to be about $1.3$ larger in the $S_{z}=N_{\rm ch}$ sector compared to the $S_{z}=0$ sector. We then assume an uncertainty of $25\%$ in these estimated correction fractions.

V Extrapolated spin gaps and inferring SC-order from correlation functions

The exact behaviour of $\overline{\Delta}_{s}[L,N_{\rm ch}]$ as $1/L$ scales to zero is critical to our approach. In V.1, we take the data as they appear to be, linearly increasing in $1/L$, with only weak remnants of the oscillatory behaviour that $\overline{\Delta}_{s}$ was defined to eliminate (c.f. Appendix A for details). It then follows that the spin gap is finite for the infinitely long strips, and non-decaying or even increasing in $N_{\rm ch}$. In V.2, we consider the possibility that what appear to be remnants of oscillations is actually the onset of a scaling of $\overline{\Delta}_{s}[L,N_{\rm ch}]$ to zero as $1/L$ decreases. In V.3, we then analyse the corrollary information that correlation functions offer.

V.1 The spin gap as a function of $N_{\rm ch}$

Applying the two distinct protocols for computing $\overline{\Delta}_{s}[L,N_{\rm ch}]$ outlined in Sec. IV.2, the behaviour of the spin gap under increasing $N_{\rm ch}$ emerges. In Figs. 5a and b, we show how we obtain spin gaps in the thermodynamic limit $L\rightarrow\infty$ through linear fits, for both $V/t=0.5$ and $V/t=1$, using the straight extrapolation protocol. Not shown are the equivalent plots for the extrapolation plus estimates protocol. There, we also employ linear fits, but now we obtain the maximal range of possible outcomes of $\overline{\Delta}_{s}[L\rightarrow\infty,N_{\rm ch}]$ by performing separate linear fits for every possible combination of extremal values of the $\overline{\Delta}_{s}[L,N_{\rm ch}]$ at different $L$. The resulting $\overline{\Delta}_{s}[L\rightarrow\infty,N_{\rm ch}]$, respectively the mean values and maximal/minimal values, of both protocols are summarized in Tab. 1 and shown in Figs. 6a and b.

Using the straight extrapolation protocol, there seems to be little room for doubt that $\overline{\Delta}_{s}[N_{\rm ch}]$ is at least a non-decaying function beyond $N_{\rm ch}=2$ for $V/t=0.5$, while it is unambiguously a monotonically increasing function at $V/t=1$. At this latter value of the nearest-neighbour repulsion, the extrapolation plus estimates protocol is benefitting from relatively small uncertainty, as here we have a greater collection of states that we could extrapolate to $\epsilon_{\rm max}=0$, and thus heuristic extrapolations were only necessary for a few states. With these smaller uncertainties, this second protocol supports the straight extrapolation protocol. While the non-decaying nature of $\overline{\Delta}_{s}[N_{\rm ch}]$ beyond $N_{\rm ch}=6$ cannot be guaranteed, a scenario in which the the spin gap suddenly reverses at larger $N_{\rm ch}$ and starts decaying towards zero seems highly unlikely and we are not aware of any mechanism that would support such a reversal. We find further evidence for this view with the significant enhancement of short-range $d_{xy}$-correlations as $N_{\rm ch}$ increases, which we discuss further in Sec. V.3.

For $V/t=0.5$ the situation is possibly more uncertain when comparing the $\overline{\Delta}_{s}[L\rightarrow\infty,N_{\rm ch}]$ resulting from the two different extrapolation protocols. The deviations are noteworthy for $N_{\rm ch}=6$, but it is very easily possible for our heuristic model of the corrections (described in Sec. IV.2) to be just off in this case. But we will observe that neither protocol is compatible with a steady decrease of the spin gap with increasing $N_{\rm ch}$. And, as for $V/t=1$, a significant and consistent enhancement of short-range $d_{xy}$ correlations at $N_{\rm ch}>2$ (see below) further supports this reading of the spin gap behaviour We thus can summarize that for $V/t=1$ the extended U-V model is highly likely to have a finite spin gap for $N_{\rm ch}\rightarrow\infty$, while for $V/t=0.5$ this likelihood, while reduced, remains high.

V.2 Testing the alternative: could $\overline{\Delta}_{s}[L\rightarrow\infty,N_{\rm ch}]$ scale to zero with $N_{\rm ch}$?

As discussed, the results for $\overline{\Delta}_{s}[L,N_{\rm ch}]$ summarized in Figs. 5a,b appear to be most congruent with a linear scaling in $1/L$ to a finite intersect at $1/L=0$, with weak oscillations surviving the averaging inherent to definition (4) (c.f. also Appendix A).

Returning to the discussion of Sec II.1, the main follow-up question then must be: could there be alternate fits that make physical sense, associated with either nodal superconductivity, or with a Fermi liquid, or a magnetically ordered state? The first possibility would be marked by $\overline{\Delta}_{s}[L,N_{\rm ch}]$ scaling to zero as $A/L^{\alpha}$, $\alpha<1$ in the large-$N_{\rm ch}$ regime, while for the latter two possibilities $\overline{\Delta}_{s}[L,N_{\rm ch}]$ should also scale to zero in the large-$N_{\rm ch}$ regime, but as $B/L+C/L^{2}$. We allow for a quadratic correction term to the ideal linear vanishing of $\overline{\Delta}_{s}[L,N_{\rm ch}]$ expected for a FL or MO state (c.f. Sec II.1) because our results from Sec. V.1 show that a purely linear scaling to zero is impossible given our data.

Analyzing our data under these alternative scenarios requires recognizing that at $N_{\rm ch}=2$ we start in a regime where we are guaranteed a finite spin gap Giamarchi (2003) (c.f. also Appendix A). That makes it extremely likely that some finite spin gap could survive for at least some larger $N_{\rm ch}$ - our linear-polynomial ansatz of the previous sections yields just that.

For this reason one cannot just blindly fit the data in Tab. 3 with these two alternate scaling ansatzes; these ansatzes will always “fit” any data like in those tables, appearing linear within the fitting region via parameter-tuning, then bending down to zero outside of it. They easily “fit” even the $N_{\rm ch}=2$ data, a nonsensical result considering these systems are known to have a finite spin gap.

Thus, the only meaningful question is whether or not these fits converge on a physically consistent scenario as $1/N_{\rm ch}$ decreases. For this purpose, we plot the fitting parameters of both ansatzes to the raw data of Tab. 3 against $1/N_{\rm ch}$ in Fig. 7. For the FL/MO-state ansatz $B[N_{\rm ch}]/L+C[N_{\rm ch}]/L^{2}$ in Figs. 7a and b, neither coefficient shows discernable convergence: $B[N_{\rm ch}]$ drops monotonically, albeit slowly, while $C[N_{\rm ch}]$ increases monotonically, at a much faster rate. Thus, within the range of accessible $N_{\rm ch}$-values, the quadratic correction term appears unstable. Rather, the rapid growth of $C[N_{\rm ch}]$ with growing $N_{\rm ch}$, coupled with the slow flattening of the linear slope controlled by $B[N_{\rm ch}]$, is that of an ansatz that increasingly struggles to approximate a behaviour of $\overline{\Delta}_{s}[L,N_{\rm ch}]$ that is actually linear or near-linear in $1/L$.

The linear-polynomial ansatz was already treated in Sec. V.1, yielding a finite spin gap. There is a motivated alternative however. A near-linear ansatz $A[N_{\rm ch}]/L^{\alpha[N_{\rm ch}]}$ with $\alpha<1$, that scales to zero gap would be suggested by earlier theory in contradiction to our result. These prior works were based on analytical perturbative RG and mean-field theory; for our 2D U-V model here, a $d_{x^{2}-y^{2}}$-order has been predicted Nickel et al. (2005, 2006); Bourbonnais and Sedeki (2011); Sedeki et al. (2012). Any order other than isotropic $s$-wave results in nodal lines intersecting with the Fermi surface (c.f. 9). The intersects result in gap-nodes, isolated points in the Brillouin zone where the SC gap closes. At these nodes zero-energy quasiparticles may thus be produced, according to the mean-field treatment of these systems, and these theories thus predict an ungapped spin excitation spectrum. This nodal structure results in the spin susceptibility, eq. (2), scaling as $1/L^{\alpha}$, with $\alpha<1$ Sigrist and Ueda (1991); Hasegawa (1996); Moriya and Ueda (2003); Shinagawa et al. (2007); Chubukov et al. (2008); Vladimirov et al. (2011).

The results of fitting $A[N_{\rm ch}]/L^{\alpha[N_{\rm ch}]}$ to our spin-gap data, depicted in Figs. 7c and d, then show a difference between $V/t=0.5$ and $V/t=1$ as $N_{\rm ch}$ increases. In the latter case both coefficients saturate quickly to nearly stable values. For $\alpha[N_{\rm ch}]$, this is crucially well below $1$. A ground state with nodal USC in the 2D regime could not be completely ruled out on this basis. To what extent this is an artefact of our data being confined to smaller $L$-values at larger $N_{\rm ch}$ must remain as a key question to be answered in future work.

For $V/t=0.5$ on the other hand, at $N_{\rm ch}\geq 4$ we cannot find saturation of the growing $\alpha[N_{\rm ch}]$ within the accessible $N_{\rm ch}$. Taken together with the linear polynomial fitting done in Sec. V.1 and the behaviour of $C[N_{\rm ch}]$ in the FL/MO-scaling ansatz above (Fig. 7b), it may not be possible to rule out that the $V/t=0.5$ data could move towards a fully linear scaling to zero in the large $N_{\rm ch}$-limit. Three conjoined caveats apply however: firstly, $N_{\rm ch}$, is limited to $\leq 8$ here, it is clearly possible that $\alpha[N_{\rm ch}]$ stabilizes below $1$ at larger $N_{\rm ch}$. Secondly, the behaviour of $\alpha[N_{\rm ch}]$ with $N_{\rm ch}$ could easily be an artefact of the oscillatory remnant in $1/L$ that $\overline{\Delta}_{s}[L,N_{\rm ch}]$ retains despite the averaging, analogous to the discussion in Appendix A for the case of a full second-order polynomial fit to the raw data. Thirdly, the results of the straight extrapolation protocol for the spin gap at $L\rightarrow\infty$ (Tab. 1, Fig. 6a) would point in the opposite direction.

We thus conclude that our data most strongly supports a fully gapped spin spectrum. As laid out in Appendix A, the averaged spin gap $\overline{\Delta}_{s}[L,N_{\rm ch}]$ appears to result in a linear structure in $1/L$ with a superimposed weak oscillatory remnant, of unknown analytical form, with its amplitude slowly decreasing with $1/L$. It is those weak downwards oscillations that the two alternate scaling ansatzes we now studied are, erroneously, susceptible to. A conservative reading of our results would still mandate that at $V/t=1$ we cannot completely rule out a spin excitation spectrum with nodal zeros, which would be compatible with the nodal USC systems described by mean-field theories. At the same time, within the confines of the $N_{\rm ch}$-values that we have treated in this work, there is no consistent scenario for a FL or MO state at $V/t=1$ due to the apparent instability of the quadratic correction as shown in Fig. 7b.

For $V/t=0.5$ the situation is less clear. The upward growth of $\alpha[N_{\rm ch}]$ parametrizing the fit to a nodal USC-regime, combined with the instability of the quadratic correction term $C[N_{\rm ch}]$ in the fit to the FL/MO-regime opens up the possibility that there may not be enough data to reliably exclude any of the three regimes for this parameter, when viewed in conjunction with the outcome of the extrapolation plus estimates protocol for the linear-polynomial extrapolation (Fig. 6a and Tab. 1). However, viewed together with the enhancement exclusive to the $d_{xy}$-correlation channel discussed in the next section, it turns out that either the nodal USC scenario or the fully spin-gapped one should be favored for $V/t=0.5$ as well.

V.3 Correlation functions: pairing in the $d_{xy}$-channel as likely candidate-order in the 2D limit

From the ground states, any desired correlation function between site ${\bf r}_{1}$ and site ${\bf r}_{2}$ in the lattice can in principle be computed. With the form-factor of the lattices we consider, strips of length $L$ and width $N_{\rm ch}$, where $L\gg N_{\rm ch}$, generally a cross-over of behaviours is to be expected with growing $r:=|{\bf r}_{1}-{\bf r}_{2}|$. As $N_{\rm ch}$ increases, at short distances the physics will become increasingly dominated by that of the 2D limit ($N_{\rm ch}\rightarrow\infty$). Conversely, at sufficiently large distances, any correlation function will exhibit the typical behaviour of a 1D-system, i.e. it will decay algebraically as a function of $r$ - as long as no discrete symmetry breaks spontaneously 333a transition to a CDW would happen in our systems eventually if $V/t$ were to be increased further. This behaviour is generic to the ground state of any 1D quantum system, even one of finite width Giamarchi (2003). The length-scale at which the cross-over between the two regimes takes place will increase as $N_{\rm ch}$ is raised, with the 1D-regime starting at larger and larger distances, until it disappears completely when $N_{\rm ch}\rightarrow\infty$.

Even with pDMRG we cannot access the full 2D regime for the time being, but we can get new and valuable insight into how the various potential channels for ordering of the U-V model in the 2D limit respond to the systematic increase of $N_{\rm ch}$. As we already evaluate the possibility of magnetic order or a Fermi liquid based on the scaling of $\overline{\Delta}_{s}[L,N_{\rm ch}]$, we limit ourselves to the main competitor-orders, the SC $d_{x^{2}-y^{2}}(r)$-, $d_{xy}(r)$-, $s(r)$- as well as the charge-density-wave channel. We compute the correlation functions (6) - (8) connected to each of these, as defined in Sec. II.1.

We then specifically study the dependency of the short-range behaviour, $r\lesssim r_{s}:=\pi/k_{F}=4a$ of the above correlation functions on $N_{\rm ch}$. An example that is highly representative of the general behaviour for $N_{\rm ch}\geq 4$ is shown in Fig. 8: when each correlation function is normalized to its value at $r=1$, the first maximum within the short-distance quasi-2D regime decays strongly in all channels except in the $d_{xy}$ channel. Table 2, which summarizes the short-range decay for all studied $N_{\rm ch}$-values, confirms this short-range enhancement as exclusive to the $d_{xy}$ channel. As is to be expected, once $r$ increases beyond the short-range regime, $d_{xy}(r)$, like all the other correlation functions, behaves according to 1D physics, i.e. it decays algebraically.

We do not take the above short-range behaviour of these systems at $N_{\rm ch}\leq 8$ to constitute definite proof of $d_{xy}$-ordering in the full 2D-limit. At the same time however, we observe that it would align with a simple order-of-magnitude estimate, namely that the diagonal singlet pairing of $d_{xy}$-order is energetically advantageous at strong coupling - naive perturbative analysis suggests that the gain in energy should be $\mathcal{O}(t_{\perp}^{2}/V)$ across a diagonal, while it is $\mathcal{O}(t_{\perp}^{2}/U)$ for singlet pairing across a rung in $d_{x^{2}-y^{2}}$ order. And of course, this would be in line with the high likelihood for a finite spin gap and the attendant elimination of competing MO or FL phases in these systems that we discussed in Secs. V.1 and V.2; the open questions raised by a system having a finite spin gap jointly with SC $d_{xy}$-order we will discuss now.

VI Discussion and connections to experiment

Our results show a fully gapped spin excitation spectrum for finite-width strips of the 2D U-V model at the given parameters, as evidenced by the finite ${\overline{\Delta}_{s}[L\rightarrow\infty,N_{\rm ch}]}$ we find.

We then show that the behaviour of the gap with the strip-width is most consistent with a fully-gapped spin spectrum also in the 2D limit when taking the apparent structure of $\overline{\Delta}_{s}[L,N_{\rm ch}]$ as consisting of a linear part with a superimposed weak oscillatory remnant into account (c.f. Sec. V.2 and Appendix A). At the same time, the behaviour of the various correlation functions at short range with increasing width, indicative of the behaviour of the 2D system, lends support to pairing happening predominantly in the $d_{xy}$-channel, rather than in the extended $s$- or $d_{x^{2}-y^{2}}$-channel.

Such a conclusion however is at variance with existing analytical descriptions of USC systems, such as from combining the results of perturbative functional RG with mean-field descriptions. Besides these methods predicting $d_{x^{2}-y^{2}}$ symmetry for the order parameter instead Nickel et al. (2005, 2006); Bourbonnais and Sedeki (2011); Sedeki et al. (2012), the more far-reaching difference to our findings lies in the effects these theories predict from the nodal lines, which any plausible order parameter will exhibit (these originate from the systems attempt to minimize the Coulomb repulsion). Specifically, they forecast that node lines’ intersections with the Fermi surface will always create point nodes for gapless quasiparticle excitations (c.f. Fig. 9). In the mean-field theory, these point nodes in turn necessitate an ungapped spin excitation spectrum.

There are two ways to look at the apparent discrepancy between this mean-field treatment and our findings of a finite spin gap in Sec. V.1. We considered the first way in Sec. V.2, observing that disregarding the apparent structure of $\overline{\Delta}_{s}[L,N_{\rm ch}]$ (linear with oscillatory remnant, c.f. Appendix A) allows for sublinear scaling of $\overline{\Delta}_{s}[L,N_{\rm ch}]$ to be a physically consistent fit within the accessible $N_{\rm ch}$-range, at least for $V/t=1$. As explained there, this scaling of $\overline{\Delta}_{s}[L,N_{\rm ch}]$ taken together with the boost to the $d_{xy}$-correlations from Sec. V.3 would be consistent with the prediction of nodal USC derived from mean-field theory just the exact symmetry of the order parameter would be different.

The second way is to consider the structure of the $\overline{\Delta}_{s}[L,N_{\rm ch}]$-data as linear + weak oscillations. Then one has to reconcile the finite spin gap in the 2D limit that ultimately results from our numerics at strong coupling with the weak-coupling analytical theory. Taking stock of the experimental situation allows just such a reconciliation. Firstly, we notice that our result of dominant $d_{xy}$ correlations in the 2D regime of the strips’ ground state would be in line with the general argument that it is energetically advantageous for the system to minimize the number of nodes, which $d_{xy}$ would achieve better than a $d_{x^{2}-y^{2}}$-order given the Fermi-surface topology (c.f. Fig. 9). Secondly, measurements on actual USC materials, which are certainly in the strongly interacting regime, provide evidence for USC order with zero-gap point nodes, yet which exhibit fully gapped spin excitations at the same time. This has been demonstrated for LSCO Lake et al. (1999) close to and away from optimal doping, as well as YBCO Dai et al. (2001) over a range of doping. To our knowledge, this fundamental discrepancy between the standard mean-field descriptions of USC models with zero-gap nodes and actual measurements so far has neither been theoretically explained nor even replicated. Thirdly, strongly underdoped yet still superconducting LSCO has been shown to have not just gapped spin excitations but to be fully gapped overall Razzoli et al. (2013). The possibilities discussed for this observation are either another phase coexisting with USC order or topological superconductivity of the chiral $d_{x^{2}-y^{2}}\pm id_{xy}$-type, which the authors show to fit the measured gap. Fourthly, the loss of $D_{4}$ symmetry in the anisotropic 2D lattice results in a mixing of orbitals, such as $d_{xy}$ and $g_{xy(x^{2}-y^{2})}$ (c.f. Fig. 9) Cho et al. (2013). When viewed in light of the numerous proposals to explain USC phases at low or zero temperature via mixing of different order symmetries, it cannot be ruled out that the remaining nodes resulting from a pure $d_{xy}+\alpha g_{xy(x^{2}-y^{2})}$ orbital, sketched in Fig. 9, are eliminated when admixed with yet other orders, even without invoking the possibility of topological superconductivity.

Given the available data we cannot currently distinguish among these possibilities, but provide them to make it plain that the most consistent interpretation of our results, a fully spin-gapped USC order (with dominant $d_{xy}$ symmetry), has clear precedents, the predictions of weak-coupling analytics notwithstanding.

VII Conclusion

We have performed large-scale calculations with parallel DMRG to obtain ground states of the 2D U-V Hubbard model, the presumed minimal model of unconventional superconductivity in the organic Bechgaard and Fabre salts. The use of parallel DMRG here allows to study large systems at sufficient accuracies to reliably compute ground state energy differences between different quantum number sectors, i.e. the spin gap.

Considering strips of increasing width up to $8$ chains, we have computed this models’ spin gap as a function of strip-length and -width, averaged so as to strongly reduce oscillations that would preclude extrapolations to infinite strip-length. After extrapolating the ground state energies in the discarded weight, we arrived at data consistent with a finite spin gap at the thermodynamic limit (i.e. at infinite length). Taking the apparent structure of the computed spin gaps as a function of inverse length into account, we conclude that the spin gap of the infinite strips is non-decreasing or even increasing with the strips width. A finite spin gap in the full 2D limit is therefore the most likely possibility, thus making the case for singlet pairing. Disregarding key features of the computed spin gaps (i.e. its decomposition into linear and weak oscillatory remnant) still yields scaling results at most consistent with nodal unconventional superconductivity, where the spin excitation spectrum would touch zero in isolated points. Even then however our results do not admit a physically consistent scaling scaling solution in line with a Fermi liquid or magnetically ordered states.

As we steer away from model parameters that could result in charge density wave order of the ground state, the remaining candidate by exclusion is unconventional (i.e. repulsively mediated) superconducting order of singlet type. This reading is further supported by the enhancement of short range correlations (i.e. where the strips’ resemble a 2D system) with the strips’ width to take place only in the $d_{xy}$-channel.

These results suggest several lines for further inquiry. One would be to test for the possibility of topological $d_{x^{2}-y^{2}}\pm id_{xy}$ superconductivity that is raised by the presence of the spin gap at strong coupling. Another would be to move more closely to the parameter regime of the Bechgaard and Fabre salts, i.e. to increase both $U/t$ and $V/t$ further yet. This will have to entail a careful study of the enlarged parameter space for the onset of the Mott-insulating state that we discussed in Sec. II.2 for the case $N_{\rm ch}=2$, and the behaviour of which at $N_{\rm ch}>2$ is currently unclear.

VIII Acknowledgements

A.K. would like to thank Ronny Thomale for useful discussions. This work was supported in part by the Swiss NSF under MaNEP and Division II, and has further received funding through an ERC Starting Grant from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 758935. The development of pDMRG was made possible by support of the Swiss government under its HP2C initiative. Carrying out the large-scale ground state calculations was made possible through two grants of compute time on the “Piz Daint” cluster of the Swiss National Supercomputing Center (CSCS), proposal IDs 20131011-18256183-84 and 20141010-173654696-445. The authors would like to thank CSCS for supporting this work. Calculations for $N_{\rm ch}=1$ and for some values at $N_{\rm ch}=2$ were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the PDC, KTH Stockholm.

Appendix A Spin gap data at finite length and scaling procedure to $L\rightarrow\infty$.

In this section we provide background on the linear scaling procedure we have used in Fig. 5 to obtain the infinite-sized gaps summarized in Tab. 1 and shown in Fig. 6. We also list all $\overline{\Delta}_{s}[L,N_{\rm ch}]$-values that entered into performing those scalings, in Tab. 3.

There are three related reasons why we perform the fitting of $\overline{\Delta}_{s}[L,N_{\rm ch}]$ in Sec. V.1 with a linear polynomial in $1/L$: (i) we find that in the noninteracting regime, $U=V=0$, the averaged spin gap $\overline{\Delta}_{s}[L,N_{\rm ch}]$ shows behaviour congruent with purely linear scaling with weak superimposed oscillations, while a fitting with a quadratic polynomial yields physically nonsensical answers. (ii) mirroring this effect for the interacting systems, we find that the amplitude of the oscillatory remnant that survives in the averaged spin gap seems to decrease so slowly with $1/L$ that the oscillations swamp any weak quadratic dependency, should there be any, For this reason, fitting the limited number of available data points with a quadratic polynomial yields unstable results, while a linear fit works stably. (iii) we know by comparing to the existing theory of the $N_{\rm ch}=2$ ladders that a linear polynomial fit yields physically consistent results that will, at minimum, provide correct trends for the behaviour of $\overline{\Delta}_{s}[L\rightarrow\infty,N_{\rm ch}]$ with $N_{\rm ch}$.

To take each of the points (i) - (iii) in turn:

(i): Consider the noninteracting version of our model Hamiltonian, $U=V=0$, at half-filling and $t_{\perp}/t=0.1$. To take $N_{\rm ch}=8$ as an arbitrary example, analysing the averaged spin gap $\overline{\Delta}_{s}[L,8]$ as a function of $1/L$ reveals an universal pattern in these tunneling-anisotropic systems: while much reduced over the non-averaged gap $\Delta_{s}[L,8]$, weak oscillations persist in $\overline{\Delta}_{s}$, as shown for this example in Fig. 10a. We find this data to be naturally decomposing into a subset of points obeying perfect linear scaling (blue crosses in Fig. 10a, linear fit superimposed as dash-dotted line), and the complementary subset of points that is affected by the downwards oscillatory deviations from this scaling (red crosses in Fig. 10a). That the seemingly linear subset of points is actually a linear one becomes apparent if one tries to fit it with a second order polynomial $C_{0}+C_{1}/L+C_{2}/L^{2}$. In Fig. 10b we plot the resulting $C_{2}$ as a function of the maximal upper fitting range, $L_{\rm max}$. It is apparent that $C_{2}$ cannot be any underlying physical property of the system like e.g. a remnant of a band-curvature that survives the averaging, in which case it should converge towards a constant value with increasing $L_{\rm max}$. Rather, $C_{2}$ continues to drop even at the large $L_{\rm max}$ shown here, because that is the only way a second-order polynomial fit can maintain optimal approximation to a purely linear distribution of data points as its fitting range grows. If any dependence on $1/L^{2}$ remains within $\overline{\Delta}_{s}[L,N_{\rm ch}]$ in the noninteracting regime, the oscillatory remnant clearly dominates it completely.

(ii): Just as for the noninteracting system, remnants of oscillatory behaviour remain visible as we compute $\overline{\Delta}_{s}[L,N_{\rm ch}]$ at $U/t=4$, $V/t=0.5,1$. But unlike the noninteracting regime, each data point now comes with substantial computational cost. The limited number of available data points compared to the non-interacting regime means that they cannot be reliably divided into a purely linear subset and an oscillatory remnant. As a consequence, the impact of the weak oscillatory remnant on any scaling $L\rightarrow\infty$ must be considered, especially as the amplitude of oscillations appears to decay only weakly with $1/L$, as in the non-interacting case (i). In Figs. 10c and d we contrast the instability of any second-order polynomial fit depending on the fitting range with the respective stability of a linear fit over the same fitting ranges. This illustrates again, in a way different from (i), how the oscillatory remnant overwhelms any quadratic component in $1/L$, if such were to exist in the interacting regime (as opposed to the non-interacting one, where we can rule out its existence within the linear subset).

(iii): Complementing the argument of (i) and (ii) so far, that $\overline{\Delta}_{s}[L,N_{\rm ch}]$ is dominated by linear behaviour in $1/L$ with weak oscillatory remnants that overshadow a potential quadratic contribution, is the physical consistency of linear extrapolation. For this, we can turn to the $N_{\rm ch}=2$ case, the 2-leg U-V ladder. The bosonized field theory that describes this ladder at low energy predicts that both its spin sectors are gapped Giamarchi (2003). It also shows that both spin gaps must grow with both $U/t$ and $V/t$ (eqs. (8.23), (8.24) and (8.26) in ref. Giamarchi (2003)). This effect is clearly recovered by our linear extrapolation, as shown in Tab. 1. This further supports that a linear extrapolation, even one that necessarily averages over the weak oscillatory remnant, correctly captures the spin gaps behaviour in the thermodynamic limit and in particular its change with the systems’ parameters.

Appendix B Handling large $\chi$ and large MPO bond dimension simultaneously through pDMRG

To make DMRG capable of handling simulations that would exhaust the RAM of any single compute node, and in order to bring the many CPU-cores of modern parallel supercomputers to bear on such large-scale problems, the groups of T. Giamarchi (U. Geneva) and M. Troyer (ETHZ) have developed pDMRG Bauer et al. , with financial support from the Swiss government under its’ HP2C initiative, with a focus on easy extendability and optimizability. Key features are (i) a clean tensor network framework for development of new algorithms; (ii) use of a matrix-like data storage type for direct use of BLAS/LAPACK routines; (iii) implementation of an arbitrary number of Abelian symmetries, reducing the matrix complexity by use of block-sparse matrices in models with conserved quantum numbers; (iv) bindings to the ALPS library Albuquerque et al. (2007); Bauer et al. (2011) for a generic model description. Parallelization uses the Message Passing Interface (MPI) standard for distributed memory and Intel® Cilk™ Plus tasks for multi-threading. The shared memory version of the code has been already published Dolfi et al. (2014), and is used extensively in 1D and 2D condensed matter physics  Dolfi et al. (2012); Bauer et al. (2012, 2013, 2014); Dolfi et al. (2015). The resulting parallelism of pDMRG is fully scalable and can be run on anything from a single CPU core to hundreds of compute nodes.

The stages of pDMRG have been summarized in Fig. 11. As in regular DMRG, the physical problem is exactly mapped onto a 1D chain. In any realistic problem requiring the use of pDMRG, there will be many long-range interaction and/or tunneling terms along the chain. As a result, the left-/right-boundary pairs $L_{n}$, $R_{n}$ at the $n=0,\dots,M$ bipartitions of the DMRG-chain will become objects requiring very large amounts of memory for any problem with a substantial bond dimension $\chi$ (Fig. 11a). For every bipartition of the DMRG-chain the DMRG-Hamiltonian will be a set of one left- and one right-boundary, together with the (usually small) local MPOs $W_{n}$, $W_{n+1}$ for the sites $n$ and $n+1$ to the left/right of the $n$-th bipartition (c.f. Fig. 11b). The technical challenge for any DMRG-implementation aiming to parallelize the repeated application of any of these Hamiltonians to a tensor-decomposed wave-function, which is the linear-algebra operation at the heart of any ground-state DMRG, is to distribute the boundaries across the nodes of a parallel supercomputer. Our pDMRG implementation offers three layers for this, summarized in Fig. 11c, making extensive use of the Message Passing Interface (MPI), the de-facto standard for parallel supercomputing:

(i) Layer 1 exploits that the boundaries are effectively vectors of block-sparse matrices (c.f. eq. (9) and following, and Fig.  3c). It distributes the $D$ elements of both these vectors, each element being a block-sparse matrix, in a size-ordered, round-robin fashion within the global group of MPI-groups. Each of the MPI-groups within the global group may encompass several nodes of the compute cluster.

(ii) Inside each MPI-group, Layer 2 then may further distribute the dense matrices comprising each boundary element among MPI-subgroups into which the main group can be split.

(iii) As dense matrices can themselves become large, they are always stored in tiled format. For Layer 3, these tiles may be distributed to different compute nodes within a subgroup.

In the Jacobi-Davidson eigensolver of pDMRG, the tensor contractions become dot-products of vectors of matrices. For this, MPI groups communicate via MPI_Allreduce. When shifting to the next pair of sites, the number of boundaries can change, which results in a redistribution of elements among the MPI-processes, using asynchronous point-to-point transfers. On each node, operations are executed in two steps. A dry-run collects all linear algebra operations (mostly GEMM and AXPY). In the second step the directed acyclic operation-graph is executed in parallel on several threads. This abstract approach optimizes performance, autonomously exploiting hidden operational parallelism. The basic linear algebra operations are then explicitly unrolled between tiles of the dense matrix and dispatched to the underlying BLAS library.