Fragmented superconductivity in the Hubbard model as solitons in Ginzburg–Landau theory

Intertwined orders come in many guises in strongly correlated electron systems^1,2,3. Continuous and discrete symmetries are broken in a manner in which order parameters transform genuinely non-trivially in the entire symmetry group. A prominent example is constituted by stripe order, where both the translational and spin rotation symmetry are broken to a state which intertwines a charge density wave (CDW) and a spin density wave. Initially proposed in Hartree-Fock studies^4,5,6,7, modern numerical techniques have in the last years succeeded in firmly establishing it as the ground state in certain parameter regimes of the strongly correlated repulsive Hubbard model in two spatial dimensions^8,9,10,11,12. As such, these orders are relevant to the physics of high-temperature superconductivity where electron-electron interactions in the form of the Hubbard model are understood to play a crucial and likely decisive role^13,14,15.

Historically, the competition between superconductivity and stripe order has been discussed thoroughly^16,17. Both stripe order and superconductivity have by now been found to be realized in different parameter regimes of the paradigmatic Hubbard and t–J models^{8,9,10,11,12,18,19,20,21,22,23,24,25}, and have also been found when further neighbor hopping processes are included^26,27,28. Moreover, numerical studies in the last years have discovered that superconductivity can also be intertwined with CDW order^{26,27,28,29,30,31,32,33,34}. A particular form of intertwined CDW order and superconductivity is the so-called pair-density wave state^{30,35,36,37,38,39,40}, for which recently evidence has been reported in several unconventional superconductors^41,42. The coexistence of CDW order, which breaks translational symmetry, and superconductivity is sometimes also referred to as supersolid order. Previously, one of the authors reported that the CDW can “fragment” the superconducting condensate⁴³. Fragmentation refers to the phenomenon when fermion pairs (or bosons in Bose-Einstein condensates) condense into not just one, but multiple macroscopic wave functions. The ground state in a certain regime of the t–({t}^{{prime} })–J model exhibits exactly one fragment per charge stripe of the superconducting condensate, interpreted as an emergent array of coupled Josephson junctions.

In this manuscript, we strengthen our understanding of this intriguing form of order. Employing density matrix renormalization group (DMRG)^44,45,46,47 simulations, we show the fragmented superconducting stripe order is also realized in the paradigmatic Hubbard model on the square lattice for systems up to width W = 6, both with open and cylindrical boundary conditions. Additionally, we investigate the effect of an orbital magnetic field on the fragments of the Cooper condensate. Moreover, we propose a minimalistic macroscopic model to describe the superconducting condensates which is given by an intertwined Ginzburg–Landau theory⁴⁸. Importantly, the mass term is chosen to be position dependent and exactly proportional to the hole-density modulation of the CDW. A one-to-one comparison between the solution of the Ginzburg–Landau equation and the numerical data from DMRG is performed and yields detailed agreement. The fragmentation of the superconducting condensate is then understood as a quantum tunneling between several soliton solutions of the Ginzburg–Landau equation.

Model and observables

We study the (t-{t}^{{prime} }-U) Hubbard model on a square lattice defined by the Hamiltonian,

where t_ij = t for nearest neighbor lattice sites i and j and ({t}_{ij}={t}^{{prime} }) for next-nearest neighbor sites. Here, ({c}_{isigma }^{dagger },{c}_{isigma }) denote the fermionic creation and annihilation operators of spin σ = ↑, ↓ and ({n}_{isigma }={c}_{isigma }^{dagger }{c}_{isigma }) denotes the number operator. We additionally study the effect of a uniform orbital magnetic field B = ∇ × A, orthogonal to the plane of the lattice. We employ the Landau gauge,

where r = (x, y) denotes a real space coordinate. The total magnetic flux through the simulation cluster is denoted by Φ = ∫_S dS B. The Landau gauge is well suited for open and cylindrical boundary conditions, since in the latter case it remains translationally invariant in the y-direction. The coupling of the model to the background magnetic field is implemented via a Peierls substitution,

where e denotes the elementary charge. In the following, we set e = 1, ℏ = 1 and focus on the model parameters U/t = 10 and ({t}^{{prime} }/t=0.2). For our geometry, the Peierls substitution reduces to t → te^i2πϕx in the y-direction, and ({t}^{{prime} }to {t}^{{prime} }{e}^{i2pi phi (x+1/2)}) as explained in ref. ⁴⁹.

In large orbital magnetic fields, this model has been studied in the context of cold-atoms experiments^50,51,52. As worked out extensively in ref. ²⁹ for U/t = 12, on six-leg cylinders this set of parameters stabilizes a superconducting state coexisting with a CDW at small hole-doping. Analogously, the coexistence of CDW and superconducting order has been reported at larger hole-doping by a combination of auxiliary-field quantum Monte Carlo and DMRG³¹. Moreover, the phase diagram at small hole-doping has been found to be in close agreement with the t–({t}^{{prime} })–J model, where similarly a coexisting CDW and superconducting state have been reported^26,27,43,53. We would like to emphasize that our findings are in full agreement with these previous studies.

Our aim is to elucidate the relation between charge and pairing degrees of freedom. We consider the hole density,

where r denotes the lattice position and n(r) the electron density. As the pairing mechanism in the lightly-doped Hubbard model has been diagnosed in previous studies to be of singlet-pairing type, we consider,

where the singlet-pairing operators ({Delta }_{{{boldsymbol{r}}}_{i}{{boldsymbol{r}}}_{j}}^{dagger }) are defined as,

Moreover, as the pairing in this superconducting phase was established as local in real space^{26,27,28,29,31,43,53}, we restrict ourselves to investigating the nearest-neighbor singlet-pairing density matrix,

where (mu =hat{{boldsymbol{x}}},hat{{boldsymbol{y}}}) (resp. ν) denote the vectors connecting nearest-neighbors on the square lattice. We exclude site-local contributions from density and spin correlations by setting ρ_S(r_i, μ∣r_j, ν) = 0 whenever the two links (r_i, r_i + μ) and (r_j, r_j + ν) are overlapping, see Eq. (8) in ref. ⁴³. Because of this, ρ_S(r_i, μ∣r_j, ν) ceases to be positive definite and therefore can have negative eigenvalues. Off-diagonal long-range order occurs whenever,

However, the information contained in ρ_S(r_i, μ∣r_j, ν) in Eq. (7) is richer and allows for additional insights beyond diagnosing long-range order. As a Hermitian matrix, ρ_S(r_i, μ∣r_j, ν) has an eigendecomposition with real eigenvalues ε_n,

The onset of Cooper pair condensation can then be inferred by the presence of either one or multiple dominant eigenvalues^43,54,55. This criterion is also known as the Penrose-Onsager criterion⁵⁵ of superconductivity. In case a single dominant eigenvalue ε₀ is observed the condensate is called simple. Otherwise, if multiple ε_k are dominant, the condensate is said to be fragmented. The dominant eigenvalues ε_k are referred to as the condensate fractions, while the corresponding eigenvectors χ_k(r, μ) are referred to as macroscopic wave functions. Notice, that the eigenvectors depend on the position of the lattice r as well as the orientation of the bond μ and are, hence, defined on the bonds of the lattice.

Ground state properties from DMRG

To investigate the onset of intertwined order in the Hubbard model we perform DMRG calculations to obtain the ground state properties of the system. In this scheme, the 2D system is mapped to a 1D chain suitable for DMRG, with the drawback of the x-direction hopping term in the Hamiltonian becoming long-range. This leads to an exponential increase in the computational cost when increasing the width W of the system, limiting the simulations to narrow cylinders. For this reason, in this work, we choose W = 6.

We obtain results for systems of length L = 24 with open boundary conditions and L = 32 with periodic boundary conditions in the y-direction. On a cylindrical geometry the Hamiltonian in Eq. (1) is not gauge invariant as a gauge transformation introduces a magnetic flux piercing the cylinder, corresponding to a ground state with a finite momentum k_y along the y-direction. As this is unphysical when looking for ground state properties we fix the gauge in order to stay in a sector with zero net k_y momentum. Analogously to previous works^26,31 we restrict to the ({S}_{tot}^{z}equiv langle {sum }_{i}{n}_{iuparrow }-{n}_{idownarrow }rangle =0) sector of total spin projection in the z direction. Exploiting this symmetry, together with the conservation of the total charge N ≡ 〈∑_in_i↑ + n_i↓〉 drastically reduces the computational load, allowing us to access larger bond dimensions. For the initial state of DMRG, we choose a physically motivated Néel product state with the holes evenly distributed in pairs along the x direction. Special care had to be taken to make sure that this choice of the initial state does not lead to local minima of the DMRG optimization. As the convergence properties are highly dependent on the initial ansatz, we also ran the calculations with different starting states to check that we are converging to a global minimum.

For the L = 24 system with open boundary conditions we insert n_h = 4 holes in the system, corresponding to a hole-doping p = 1/36. We start by stabilizing the ground state without any magnetic flux, Φ = 0, at a bond dimension D = 3500, and then scale up the magnetic flux to Φ = 4π at D = 4000, using the previously obtained ground state as the initial state for DMRG. For the L = 32 system with cylindrical boundary conditions, we insert n_h = 12 holes, corresponding to a hole-doping of p = 1/16, and we focus on the Φ = 0 case only, to better showcase how the different superconducting condensates arrange themselves over the different stripes and demonstrate consistency with the intertwined Ginzburg–Landau results. For this larger system we show results at a bond dimension D = 8000.

We plot the spectrum ε_k of ρ_s as defined in Eq. (9) in Fig. 1 for both parameter sets. In analogy with our previous results on the t–({t}^{{prime} }-J) model⁴³, we observe multiple dominant eigenvalues above a continuum. Again, the number of dominant eigenvalues matches exactly the number of CDW maxima of the ground state. Thus, the system exhibits the presence of n_h/2 fragmented condensates, irrespective of the boundary conditions and the presence or absence of a background magnetic field at these parameter values.

a 24 × 6 square lattice (open boundary conditions) doped with n_h = 4 holes away from half-filling and total magnetic flux Φ = 4π through the full sample.Two dominant eigenvalues are observed. b 32 × 6 cylinder doped with n_h = 12 holes without magnetic flux. We observe six dominant eigenvalues. In both cases, the number of dominant eigenvalues exactly matches the number of CDW maxima.

Full size image

To show the interplay between the CDW and the superconducting order we plot the hole density n_h(r) superimposed on the macroscopic wave functions. For the L = 24 system with open boundary conditions, we show the macroscopic eigenvectors χ₀(r, μ), χ₁(r, μ) in Fig. 2(a), (b), where the angle and the color of the arrows represent the phase of χ_k(r, μ) and the length of its amplitude respectively. We observe the macroscopic wave functions to be supported on the hole-rich regions of the ground state. Compared to an expected d-wave pattern for Φ = 0, a finite magnetic field induces rotation in the complex phase that increases by increasing magnetic flux Φ.

a, b Macroscopic wave functions of ground state of the (t-{t}^{{prime} }-U) Hubbard model from DMRG with total magnetic flux Φ = 4π on a 24 × 6 lattice with open boundary conditions when doped with n_h = 4 holes away from half-filling. Here we take ({t}^{{prime} }/t=0.2), U/t = 10. a (resp. b) shows the values of χ_k(r, μ) for the k = 0 (resp. k = 1) fragment of the condensate. The length of the arrows is proportional to ∣χ_k(r, μ)∣ and the color indicates the complex angle of χ_k(r, μ), where for (mu =hat{{boldsymbol{y}}}) we shift the phase by π, to reflect the d-wave nature of the order parameter. The grayscale of the background squares shows the value of the local hole density n_h(r). c, d Both soliton solutions ψ_k(r, μ) (k = 0, 1) of the intertwined Ginzburg–Landau theory for α(r) = αn_h(r), where n_h(r) is taken from the values of the DMRG ground state and the free parameters α = − 2.331, β = 14.465 are found by matching the DMRG data to Ginzburg–Landau theory. We observe close agreement between both χ_k(r, μ) and ψ_k(r, μ).

Full size image

For the L = 32 system with cylindrical boundary conditions the three macroscopic eigenvectors χ₀(r, μ), χ₂(r, μ), χ₅(r, μ) are shown in Fig. 3, both in the x − y plane and averaged over the y-direction via,

Since χ_k(r, μ) is real for Φ = 0, we omit the arrows and only employ the color to define the sign. In all cases, the d-wave pattern of the pairing and the disposition of the superconducting pairing along the density stripes are clearly visible in alternating red and blue colors for (mu =hat{{boldsymbol{x}}},hat{{boldsymbol{y}}}). In particular, for every different condensate, the macroscopic wave function has a different spatial modulation and correspondingly a different number of nodes.

a, c, e Macroscopic wave functions χ_k(r, μ) as in Eq. (9) with index k = 0, 2, 5. The color and width of the bars denote the value of χ_k(r) on the links of the lattice. b, d, f y-averaged hole-density ({bar{n}}_{h}=({sum }_{y}{n}_{h}(x,y))/W) of the ground state (gray, dotted) and comparison between the y-average of the macroscopic wave functions ({bar{chi }}_{k}({boldsymbol{r}})) from DMRG and model wave functions ({bar{psi }}_{k}({boldsymbol{r}})) from the intertwined Ginzburg–Landau theory with α = − 1.3 and β = 13. We report close agreement between the DMRG data and the solutions from Ginzburg–Landau theory.

Full size image

Intertwined Ginzburg–Landau theory

In the following, we demonstrate that the macroscopic wave functions χ_k(r, μ) of the fragmented Cooper condensates are well described by solutions of an intertwined Ginzburg–Landau functional of the form,

The key aspect in Eq. (11) is the position dependence of the mass term α(r), which is chosen to be,

so the charge and phase degrees of freedom are intertwined through directly coupling to the hole density n_h(r). This action is invariant under a gauge transformation,

reflecting the transformation rule of Cooper pairs. Minimizing the functional ({mathcal{F}}[psi ]) by solving,

yields the time-independent Ginzburg–Landau equations,

Importantly, the criterion Eq. (14) from which Eq. (15) is derived is valid for every local minimum of the functional Eq. (11). Hence, whenever there are multiple solutions to Eq. (15) there are also multiple local minima of Eq. (11). In the absence of both the nonlinear term (β = 0) and the vector potential (A = 0), the Ginzburg–Landau equation (Eq. (15) reduces to the linear time-independent Schrödinger equation. For periodic boundary conditions, in the presence of a periodic modulation of the mass term,

Bloch’s theorem⁵⁶ states that this equation has multiple solutions labeled by the wave number k,

where u(r + λ) = u(r). Naturally, the question arises whether an analog of Bloch’s theorem still holds in the presence of a nonlinearity (β ≠ 0). This has already been addressed in the literature, where generalizations of Bloch’s theorem with additional nonlinearities have been proven^57,58,59. Note that nonlinear Schrödinger equations in periodic potentials naturally occur when describing Bose-Einstein condensates in optical lattices, see e.g.^59,60,61,62.

In the presence of both a nonlinearity and the magnetic vector potential in two dimensions, we solve Eq. (15) numerically. Our approach is to solve the equivalent problem of minimizing ({mathcal{F}}[psi ]) via a local minimization algorithm and verifying that the functional derivative in Eq. (14) vanishes. Different solutions corresponding to distinct Bloch waves can be obtained by starting the optimization procedure with distinct initial configurations within the basin of attraction of the particular solution, cf. Sec. Details on solving the Ginzburg-Landau equations. We employ the BFGS algorithm^63,64,65,66 for local optimization and verify the minimization condition Eq. (14) by monitoring the norm of the Jacobian matrix, which we require to attain a numerical value < 10⁻⁷ in all studied geometries.

To compare the solutions of the intertwined Ginzburg–Landau equation (Eq. (15) to the macroscopic wave functions from DMRG, we solve for ψ_k(r) on the vertices of the square lattice. We then define the d-wave order parameter ansatz on the links of the lattice as,

With this definition, the d-wave nature of the order parameter ψ_k(r, μ) is predetermined. We compare the intertwined Ginzburg–Landau wave functions with the numerical results in Fig. 2b for a finite magnetic flux ϕ. At this set of parameters, DMRG yields two fragments of the condensate while the Ginzburg–Landau equation exhibits two distinct solutions, which we find using different initial states when minimizing the Ginzburg–Landau functional, cf. Fig. 6. Thus, there is a one-to-one correspondence between the superconducting fragments and the distinct solutions of the Ginzburg–Landau theory. Moreover, we observe that both the periodicity in χ₀(r, μ) and χ₁(r, μ) and the rotation due to the magnetic field are correctly captured by the Ginzburg–Landau solutions ψ₀(r, μ) and ψ₁(r, μ) and observe accurate agreement in Fig. 2. To obtain such an agreement, the two model parameters α and β have to be optimized. We performed a full parameter scan to determine optimal values of α = − 2.331, β = 14.465, cf. Supplementary Fig. 4.

Next, we consider the case of multiple stripes without a magnetic field (Φ = 0) on a 32 × 6 cylinder. As we only obtain DMRG results for open boundary conditions in the x-direction, the corresponding Ginzburg–Landau equation coupling to the hole density from DMRG is not translationally invariant. Thus, we do not immediately expect Bloch waves as solutions. Instead, we find solutions that are localized on the individual stripes and compute their Fourier transform at a momentum k of the superlattice set by the charge density wave which serves as ansatz wave functions ψ_k(r, μ) to compare to DMRG. A detailed discussion of this construction is given in Supplementary Note 4. Figure 3 again shows a close agreement between the ansatz wave functions ψ_k(r, μ) from the intertwined Ginzburg–Landau theory and the condensate fragments χ_k(r, μ). Moreover, this comparison reveals that the different fragments of the condensate can be labeled by the quasi-momentum of the wave function on the superlattice given by the charge density wave. This is the reason we used the label ‘k’ for both enumerating the fragments χ_k(r, μ) and model wave functions ψ_k(r, μ). In fact, we observe a ‘dispersion’, in the sense that the smaller values of k have a larger condensate fraction, i.e. ({k}^{{prime} }, < ,k) implies ({varepsilon }_{{k}^{{prime} }} ,> ,{varepsilon }_{k}), reflecting the fact that the uniform condensate is the most dominant condensate.

The intertwined Ginzburg–Landau theory allows us to make predictions for larger systems not accessible by DMRG simulations. As an example, we report in Fig. 4 the solutions of Eq. (11) on a 16 × 16 grid with a CDW of the form ({n}_{h}({boldsymbol{r}})={sin }^{2}(2pi x/L)). For this larger size we can see how, on the k = 0 condensate, vortices are effectively pinned between the stripes for the chosen set of parameters α = − 2.331 and β = 14.465. This set of parameters is the optimal choice when fitting to the DMRG results in Fig. 2. However, we expect the behavior of the intertwined Ginzburg–Landau equation to be more complex in general. Besides the conventional coherence length and penetrations depth, another natural length scale that can be defined is the period of the charge modulation. Moreover, we also expect the depth of the potential to play an important role. A more detailed study of the behavior when varying these scales will be the subject of future investigations.

Both solutions ψ_k(r, μ) of the intertwined Ginzburg–Landau equation (Eq. (11) with k = 0 (a) and k = 1 (b). We choose α = − 2.331 and β = 14.465 on an L × L lattice with L = 16 with a total magnetic flux of Φ = 1.64 ⋅ 2π through the sample and open boundary conditions. The CDW is assumed to be of the form ({n}_{h}({boldsymbol{r}})={sin }^{2}(2pi x/L)). The grayscale of the background squares is set by the hole density n_h(r). The color and length denote the phase and amplitude of ψ_k(r, μ), where for (mu =hat{{boldsymbol{y}}}) we shift the phase by π to reflect the d-wave nature of the order parameter. We observe pinned vortices between the stripes in the k = 0 condensate.

Full size image

The coexistence of CDWs and superconductivity has by now been reported by numerous numerical studies of doped Hubbard- and t–J-like models^{26,27,28,29,30,31,32,33,34}. By employing the Penrose-Onsager criterion for superconductivity, one of the authors previously reported that this coexistence yields a fragmentation of the superconducting condensate in the case of the t–({t}^{{prime} })–J model on the square lattice⁴³. This is seen from investigating the eigenvalues and eigenvectors of a suitably chosen two-body density matrix, which correspond to the condensate fractions and the macroscopic wave functions of the Cooper pairs, where more than one dominant eigenvalue is observed. In the present manuscript, we corroborated these findings by demonstrating that fragmentation is not just a peculiarity of the t–J model but is observed in the phase diagram (t-{t}^{{prime} }-U) Hubbard model. More importantly, we propose a macroscopic field theory that describes the wave functions of the fragmented condensate with remarkable precision. The key proposition is the intertwined Ginzburg–Landau theory Eq. (11), where the mass term of the superconducting order parameter is coupled to a static field proportional to the hole density of the system. The CDW is seen as a periodic background potential that leads to the existence of multiple solutions of the Ginzburg–Landau theory, which correspond to Bloch waves at different momenta. These distinct solutions of the intertwined Ginzburg–Landau theory have been shown to precisely describe the individual fragments of the superconducting condensates, cf. Figs. 2, 3. We would like to emphasize that the model we propose only features two parameters α and β (up to an overall constant factor m^*) and can therefore be considered a minimal model coupling the charge density to the superconducting order.

Our findings give rise to an intuitive picture of the behavior of pairs in the stripe-fragmented superconductor. Here, a Cooper pair is in an unequal-weight superposition of the condensate wave functions at different momenta. The individual macroscopic wave functions constitute local minima of the Ginzburg–Landau free energy functional Eq. (11) and do not necessarily yield the same value of the Ginzburg–Landau free energy. Nevertheless, the Cooper pair can “tunnel” between these different minima. An energy difference between these minima is then reflected by a different occupation number in these modes, just as we observe in our DMRG data. We numerically verified that the uniform condensate solution is a global minimum of the Ginzburg–Landau functional, which explains why this mode is populated with the highest condensate fraction, whereas finite momentum solutions have higher values of the Ginzburg–Landau free energy and thus a smaller condensate fraction.

A coupling between a CDW and the superconducting order parameter has been suggested previously^{67,68,69,70,71,72}. Importantly, novel pump-probe experiments allow for a precise measurement of this coupling employing Fano interference^73,74, and a coupling between the CDW and the superconducting order has recently been measured in NbSe₂⁷⁵ and the cuprate superconductors^73,74. Thus, there is strong evidence that coupling between CDWs and superconducting order parameters is realized in actual strongly correlated superconductors.

The final question we would like to address is whether a fragmented Cooper condensate could be observed in materials. While our particular choice of parameters in the (t-{t}^{{prime} }-U) might not be the set of couplings expected in e.g. the cuprate superconductors, it is very likely that such supersolid phases with fragmented condensates can occur in more complex scenarios, which then capture the physics of actual materials more accurately. Several studies employing scanning Josephson tunneling spectroscopy have recently reported the emergence of a pair-density wave state³⁷ in Bi₂Sr₂CaCu₂O_8+x^41,76, in transition metal dichalcogenides⁷⁷, or the heavy-fermion superconductor UTe₂⁴². A peak in the superconducting order parameter at finite momentum has been detected, indicating the emergence of a pair density wave. Besides this finite momentum peak, however, a zero momentum contribution was also measured. We would like to point out that this scenario bears a strong resemblance to our findings, as we propose that the Cooper pairs are in a superposition of condensates with different momentum quantum numbers. Hence, if a fragmented superconductor were to be realized, further peaks at different momenta could indicate this state. Moreover, a prediction of our intertwined Ginzburg–Landau theory would be that the superconducting vortices in a magnetic field are pinned between the CDWs, as shown in Fig. 4, which we propose as a potential hallmark experimental signature of this phase.

Details on DMRG simulations

To ensure consistent convergence of DMRG for the systems we have studied, some care must be taken in choosing the initial states of the variational optimization. Here we describe our strategy for obtaining initial states. For systems of width W = 4, we found that DMRG is not sensitive to the initial starting states and we use random initial states. Independent of the starting state, we find that the ground state is antiferromagnetic with a checkerboard spin pattern, and the holes bind in pairs and distribute evenly along the width of the system.

For width W = 6, we use starting states that are inspired by the ground states we find with DMRG at W = 4, which we find provide better convergence than random initial states. We have tried different starting states to test that our choice of starting state doesn’t lead to a biased result. As an example of this, in Fig. 5(a), we show results from DMRG calculations of an L = 24 system with open boundary conditions and an equal number of total holes which were biased by the starting states to have different hole configurations (two holes per stripe with more stripes or four holes per stripe fewer stripes). By extrapolating the energy as a function of DMRG truncation error⁷⁸, we can conclude that the state with 2 holes per stripe is lower in energy.

a Scaling of the energy with the truncation error for two different ansatzes in a 24 × 6 system with four holes. b Energy behavior for different ansatzes at fixed bond dimension χ = 4000 for a 32 × 6 system with 12 holes.

Full size image

In Fig. 5b, for an L = 32 system with cylindrical boundary conditions and n_h = 12 holes, we compare the variational minimum energy for three different hole configurations (two, four, and six holes per stripe) at fixed bond dimension and as before we see that the energy is lowest for the state with two holes per stripe. Throughout this work, we therefore show results using states that have two holes per stripe.

A finite magnetic field ϕ in the Hamiltonian Eq. (1) introduces complex coefficients. In DMRG, this leads to tensors that have complex elements, and therefore the computation time increases by a factor of approximately four compared to the case of zero magnetic field (ϕ = 0) where tensors with real elements can be used. To improve the convergence time of calculations with finite magnetic field, we first compute the ground state at some fixed bond dimension without a magnetic field and use the state found by DMRG as a starting state for calculations at finite magnetic field.

Care must be taken to efficiently compute the two-body singlet-pairing density matrix ρ_S. Naively, without the use of caching and sparsity, computing every element of a general two-body density matrix on a system of size (N={mathcal{O}}(L)) (where for simplicity we assume a fixed width W) would scale as ({mathcal{O}}({N}^{5})) since there are ({mathcal{O}}({N}^{4})) elements and computing a single element requires ({mathcal{O}}(N)) tensor contractions, which would quickly become impractical. Luckily, there are multiple ways we can reduce this scaling. First of all, we only need to compute the pairing on neighboring sites, which reduces the scaling to ({mathcal{O}}({N}^{3})) if no caching is used. If caching of intermediate tensor contractions is used across the computation of different elements, this scaling can be further reduced to ({mathcal{O}}({N}^{2})). The code we use, which can compute general n-body correlators and automatically cache intermediate tensors involved in computing different elements, can be found at⁷⁹.

Details on solving the Ginzburg–Landau equations

We describe how to numerically attain the distinct solutions corresponding to the distinct Bloch wave instanton solutions of Eq. (15) in the main text. These are computed by numerically minimizing the parent Ginzburg–Landau functional Eq. (11) with different initial starting vectors. The procedure is illustrated in Fig. 6. To obtain a uniform solution we start with a constant initial state ({psi }_{0}^{0}), e.g.

as shown in Fig. 6a. ψ(r) is a function with complex values on the vertices of the lattice. The values of hole-density n_h(r) on the lattice are taken from measurements of the ground state from DMRG. We obtain a local minimum of the Ginzburg–Landau functional Eq. (11) by employing the BFGS algorithm^63,64,65,66 as implemented in SciPy⁸⁰. Convergence is tested by postulating the norm of the Jacobian matrix to be of size < 10⁻⁷. This allows us to obtain the wave function ψ₀(r), as shown in Fig. 6b. To compare this to the macroscopic wave function from DMRG, we compute the corresponding d-wave order parameter on the links of the lattice by Eq. (18). The final model wave function ψ₀(r, μ) is then shown in Fig. 6(c).

a, b Initial ansatz wave functions ({psi }_{k}^{0}({boldsymbol{r}})) used as a starting point in the minimization of the Ginzburg–Landau functional Eq. (11). c, d Obtained local minima ψ_k(r) from these initial states after optimization. e, f Resulting d-wave order parameters ψ_k(r, μ) as in Eq. (18) defined on the links of the lattice.

Full size image

To compute different solutions we initialize the minimization procedure with a starting state of the form,

The solution shown in Fig. 2d of the main text is obtained when setting k = 1. The corresponding initial configuration is shown in Fig. 6d. After finding the local minimum associated with this initial solution, we obtain ψ₁(r) as shown in Fig. 6e, and finally the d-wave order parameter ({psi }_{1}^{d}({boldsymbol{r}})) shown in Fig. 6f.