Vortex core order and field-driven phase coexistence in the attractive Hubbard model

TL;DR

This paper investigates how magnetic fields induce coexistence of charge density wave order and superconductivity in the attractive Hubbard model, revealing a transition to a supersolid phase and mapping the phase diagram.

Contribution

It demonstrates the formation and evolution of vortex-induced charge order and identifies a field-driven supersolid phase in the attractive Hubbard model.

Findings

Vortices seed charge density wave regions.

Charge-ordered vortex cores overlap at high flux.

Superconductivity transitions to a charge-ordered state.

Abstract

Superconductivity occurs in the proximity of other competing orders in a wide variety of materials. Such competing phases may reveal themselves when superconductivity is locally suppressed by a magnetic field in the core of a vortex. We explore the competition between superconductivity and charge density wave order in the attractive Hubbard model on a square lattice. Using Bogoliubov-deGennes mean field theory, we study how vortex structures form and evolve as the magnetic flux is tuned. Each vortex seeds a CDW region whose size is determined by the energy cost of the competing phase. The vortices form a lattice whose lattice parameter shrinks with increasing flux. Eventually, their charge-ordered vortex cores overlap, leading to a field-driven coexistence phase exhibiting both macroscopic charge order and superconductivity -- a `supersolid'. Ultimately, superconductivity disappears via a first-order phase transition into a purely charge ordered state. We construct a phase diagram containing these multiple ordered states, using $t'$, the next-nearest neighbour hopping, to tune the competition between phases.

Equations

$$H=\sum _{⟨ij⟩,\sigma }\left\{-t_{ij}{c}_{i,\sigma }^{†}{c}_{j,\sigma }+h.c.\right\}-U\sum _i{\hat{n}}_{i,\uparrow }{\hat{n}}_{i,\downarrow }-\mu \sum _{i,\sigma }{\hat{n}}_{i,\sigma }.$$

$$\mathcal{L}=\frac{\rho }{2}{\left|\left(\nabla -\frac{ie}{\hslash c}\mathbf{A}\right)\Delta (\mathbf{r})\right|}^2+\frac{1}{8\pi }{\left(\nabla \times \mathbf{A}\right)}^2$$

$$+\frac{\rho }{2}|\nabla \tilde{\varphi }(\mathbf{r}){|}^2-|\Delta (\mathbf{r}){|}^2-(1-g{t}^{\prime 2})|\tilde{\varphi }(\mathbf{r}){|}^2.$$

$$H_{MFT}$$

$${\Delta }_i$$

$$n_{i\uparrow }$$

$$n_{i\downarrow }$$

Full text

Introduction:

Superconductivity is often obtained in proximity to other ordered ground states. The most prominent example being the high Tc cuprates, where superconductivity competes with antiferromagnetism and with charge orderLake et al. (2001); Chang et al. (2012a). A particularly interesting way to stabilize underlying competing phases is to apply a magnetic field, locally suppressing superconductivity to create vortices. The core region of the vortex can then host competing correlationsHoffman et al. (2002); Curran et al. (2011); Machida et al. (2016). Indeed, experiments with scanning tunnelling microscopy have revealed charge-orderedHoffman et al. (2002); Machida et al. (2016) vortex cores in the cuprates. NMR studies of YBa2Cu3Oy indicate that as the magnetic field increases, the inter-vortex distance decreases; at a critical field strength, vortex cores overlap leading to charge order throughout the systemWu et al. (2013). These and related experiments motivate the study of vortex core order and field-driven coexistence in the attractive Hubbard model, the simplest model to show competition between superconductivity (SC) and charge density wave (CDW) order.

Hubbard model and $SO(3)$ symmetry:

We consider fermions on a square lattice, described by

[TABLE]

where $\mu$ is the chemical potential and $U$ is the strength of the on-site attractive interaction ($U>0$). The hopping parameter $t_{ij}$ takes the value $t$ (henceforth set to unity) for nearest neighbours and is zero otherwise. When $\mu$ is tuned to half-filling, this model possesses a remarkable symmetry with SC and CDW order becoming degenerate; the order parameters form an enlarged space having $SO(3)$ symmetry as shown in Fig. 1(left)Yang (1989); Zhang (1990); Yang and Zhang (1990); Burkov and Paramekanti (2008); Ganesh et al. (2009). This is a delicate symmetry arising from the bipartite nature of the square lattice with hoppings connecting sites of different sublattices. We can tune away from this $SO(3)$ degenerate point by introducing a next-nearest neighbour hopping, $t^{\prime}$. The $t^{\prime}$ term lowers the energy of the SC phase relative to CDW phase.

The $SO(3)$ degeneracy leads to a local pseudospin order parameter whose components are $\left\{\mathrm{Real}(\Delta_{i}),\mathrm{Imag}(\Delta_{i}),\tilde{\phi}_{i}\right\}$ as shown in Fig. 1. Here, $\Delta_{i}$ and $\tilde{\phi}_{i}$ (defined below) are the local superconducting and CDW order parameters. This $SO(3)$ symmetry is directly analogous to the hypothesized $SO(5)$ symmetryDemler et al. (2004) in the cuprates which groups SC and antiferromagnetism into an enlarged order parameter space. As a testable consequence of $SO(5)$ theory, it was proposed that vortex cores would have antiferromagnetic orderArovas et al. (1997); Hu and Zhang (2002). Analogously, the Hubbard model in Eq. 1 will possess CDW order in the vortex core. In the language of $SO(3)$ pseudospins, a vortex corresponds to a ‘meron’, as shown in Fig. 1(right) – in the core region, the moments cant out of the plane to locally give rise to CDW order. Here, unlike in the cuprates, we have direct control over the $SO(3)$ symmetry breaking in the form of the $t^{\prime}$ hopping. We study the Hubbard model in an applied field demonstrating CDW order in the vortex core. Our key result is a field-driven SC-CDW coexistence regime which arises from the overlap of vortex cores. This state simultaneously breaks translational symmetry and $U(1)$-gauge symmetry, demonstrating a new route to ‘supersolidity’Boninsegni and Prokof’ev (2012).

Bogoliubov deGennes mean-field theory:

We perform simulations on an $L\times L$ lattice with periodic boundary conditions, with $L$ up to $30$. To introduce an orbital magnetic field, we add a complex phase to the hopping amplitudes $t_{ij}$ given by $\theta_{ij}=e\int_{\bf r_{i}}^{\bf r_{j}}{\bf A.dr}$, where ${\bf A}({\bf r})$ is the vector potential, see Supplementary Materials for details. The net magnetic flux through a closed surface must be quantized in units of $h/e$Dirac (1931). We take the net flux through our system to be $\alpha h/e$ where $\alpha$ is an integer. As each vortex carries a flux $\Phi_{0}=h/2e$, we will always have an even number of vortices in the system. In particular, the lowest magnetic flux we can have is $2\Phi_{0}$, corresponding to two vortices.

We decompose the on-site interaction term in pairing and density channels. The SC order parameter is complex-valued, defined as $\Delta_{i}=U\langle c_{i\downarrow}c_{i\uparrow}\rangle$. The density order parameter is defined as $\phi_{i}=\frac{U}{2}\langle\hat{n}_{i\uparrow}+\hat{n}_{i\downarrow}\rangle=\frac{U}{2}(\langle c_{i\uparrow}^{\dagger}c_{i\uparrow}\rangle+\langle c_{i\downarrow}^{\dagger}c_{i\downarrow}\rangle)$. The local CDW order parameter can be defined as $\tilde{\phi}_{i}=(-1)^{\mathbf{r}_{i}}\{\phi_{i}-U/2\}$, which measures the local deviation from half-filling. With these mean-field parameters, the Hamiltonian takes the form of a $2L^{2}\times 2L^{2}$ matrix, which can be diagonalized using the Bogoliubov-Valatin transformation Chen and Ting (2002, 2005); Han (2010). We obtain self-consistent values of $\Delta_{i}$ and $\phi_{i}$ on every site.

We find several self-consistent mean-field configurations, of which the one with lowest energy is to be chosen. One solution is a pure CDW state in which $\Delta_{i}=0$ for all $i$ and $\phi_{i}=\left\{\phi_{0}+(-1)^{\mathbf{r}_{i}}\tilde{\phi}\right\}$, corresponding to uniform CDW order. In the absence of a magnetic field and in the presence of a non-zero $t^{\prime}$, this state has higher energy than the uniform SC phase. When a field is imposed, this state is not affected as it is insulating – its energy remains constant, independent of the flux (see Supplementary Materials). In contrast, the SC phase necessarily develops vortices when a field is imposed. As the number of vortices increases with flux, so does the energy of the SC. As seen from these energetic arguments, an applied magnetic field induces competition between SC and CDW orders.

The $SO(3)$ symmetry of the attractive Hubbard model only exists precisely at half-filling. As we are interested in phase competition, all results presented here are at half-filling. We present results for $U=10t$ for the following reason. At large $U$, the Hubbard model can be mapped to a spin problem with antiferromagnetic superexchange interactionsBurkov and Paramekanti (2008); Ganesh et al. (2009). The local order parameter is, in fact, the $SO(3)$ spin whose components are $\left\{\mathrm{Real}(\Delta_{i}),\mathrm{Imag}(\Delta_{i}),\tilde{\phi}_{i}\right\}$ as shown in Fig. 1. At low temperatures, we expect the system to have uniform spin length, ie., $|\Delta_{i}|^{2}+\tilde{\phi}_{i}^{2}=c$, a constant independent of position. SC and CDW order parameters are not independent, having to satisfy this uniform-length constraintArovas et al. (1997). With these considerations, the appropriate Landau Ginzburg free energy density is given byArovas et al. (1997)

[TABLE]

The order parameters are coupled by the uniform length constraint: $|\Delta(\mathbf{r})|^{2}+\tilde{\phi}^{2}(\mathbf{r})=c$. SC and CDW orders become degenerate when $t^{\prime}=0$ and the magnetic field is turned off, revealing the underlying $SO(3)$ symmetry. At $U=10t$, we find that the mean-field results always satisfy the uniform spin length constraint and the above Landau Ginzburg theory applies. We find the same qualitative results extending to small $U$ values as well.

Motivated by recent experiments revealing charge order in the cuprates, several authors have studied field theories similar to Eq. 2Efetov et al. (2013); Hayward et al. (2014); Wachtel and Orgad (2014) with Ref. Meier et al., 2013 also incorporating an orbital magnetic field. Our study of the attractive Hubbard model at strong coupling can be viewed as an ultraviolet regularization of such a field theory.

Vortex profile:

Setting $\alpha=1$, we obtain the lowest flux configuration with two well-separated vortices. As $t^{\prime}$ is increased, we find CDW order in the vortex core until $t^{\prime}\lesssim 0.5t$. For larger $t^{\prime}$ values, we find a normal core with no CDW correlations. Figs. 2(a) and 2(b) show the profiles of superconducting and CDW order at selected values of $t^{\prime}$. The lower panels of Fig. 2 show the spatial maps of SC and CDW order parameters around a single vortex for $t^{\prime}=0.3t$: CDW correlations can be clearly seen in the vortex core region. The same information is presented in spin language in Fig. 1(right).

The SC and CDW profiles are, in fact, set by the same length scale, $\xi$, as the order parameters satisfy the uniform spin-length constraint. We obtain $\xi$ by fitting the SC profile to $\Delta(x)\sim\Delta_{0}\tanh(x/\xi)$; this functional form is consistent with the free energy in Eq. 2Arovas et al. (1997). The resulting $\xi$ is plotted in the inset to Fig. 2(a). Separately, we define two length scales, $L_{\Delta}$ and $L_{\phi}$, as the full-widths at half-maximum of SC and CDW profiles respectively. We find that $L_{\phi}$ is always larger than $L_{\Delta}$ as shown in the inset to Fig. 2(b), although there is only one underlying length scale, $\xi$. Superconductivity is suppressed in the vortex core out to a radius set by $L_{\Delta}$. However, each vortex also hosts CDW correlations which extend to a larger distance, $L_{\phi}$. This sets the stage for a coexistence phase, which can be understood as follows. As the field is increased, the vortices are packed more and more tightly. Naïvely, superconductivity persists until the inter-vortex distance approaches $L_{\Delta}$. Much before this, when the inter-vortex distance reaches $L_{\phi}$, the CDW regions around each vortex overlap. CDW order percolates throughout the system on top of the SC background, leading to a field-driven coexistence phase as we demonstrate below. Similar arguments have recently been put forward in the cuprates based on NMR resultsWu et al. (2013).

Vortex lattice evolution:

Fig. 3 shows our results for $t^{\prime}=0.2t$ with varying $\alpha$ (the total magnetic flux being $\alpha h/e$) on a $24\times 24$ lattice. The panels show real space maps of $\Delta_{i}$ and $\phi_{i}$, showing the evolution of a vortex lattice with increasing flux. In addition, we plot the Fourier transform of the SC order parameter, defined as $S_{\Delta}({\bf q})=(1/N)\sum_{i}|\Delta_{i}|^{2}e^{i\bf q\cdot\bf{r}_{i}}$. The distribution of peaks in $S_{\Delta}({\bf q})$ reveals the geometry of the vortex lattice. The figure also plots the electronic density of states in the mean field ground state.

For small fields, with $\alpha=2,4$, we find well separated vortices forming an anisotropic triangular lattice. At $\alpha=6$, the lattice becomes near-isotropic. Upon increasing the field to $\alpha=8$, the vortices form a square lattice. This suggests an underlying phase transition driven by tuning vortex density. A similar transition into a second vortex lattice phase has been suggested in YBCO from torque magnetometry resultsYu et al. (2016). For $\alpha>8$, we find phase separation into square and triangular vortex lattices. Finally, at $\alpha=12$, we find a first order phase transition to a pure CDW phase, which has lower energy than solutions with SC order. Thus, at mean-field level, $H_{c2}$ is set by the competing CDW phase, unlike in conventional superconductors.

Phase coexistence:

As seen in Fig. 3, every vortex core nucleates CDW correlations which begin to overlap when $\alpha=4$. The CDW order becomes progressively stronger with increasing field as vortex cores overlap more and more. For $\alpha\gtrsim 8$, we have near-uniform CDW order. With increasing field, the SC order weakens while the CDW order parameter grows. As a result, the electronic gap never closes, as shown in Fig. 3.

Fig. 4 shows the in-field phase diagram for $t^{\prime}=0.2t$. As we force a magnetic flux through Peierl’s substitution, there is no $H_{c1}$ in our simulations. To quantify the strength of CDW order, we define the parameter $\rho=\rho_{(\pi,\pi)}/\rho_{(0,0)}$ where $\rho_{\bf q}$ is the Fourier component of the density order parameter: $\rho_{\bf q}=\sum_{i}\phi_{i}e^{i{\bf q}\cdot{\bf r}_{i}}$. In a pure CDW state with site occupation oscillating between 0 and 2, $\rho$ takes the value unity. For small $\alpha$ values, we find $\rho$ to be small, indicating weak CDW order arising from well separated vortex cores. With increasing $\alpha$, $\rho$ increases monotonically. Within our mean-field theory, we find that CDW correlations begin to span the system when $\rho\sim 0.3$. Based on this observation, we use $\rho\gtrsim 0.3$ as a heuristic criterion to signal macroscopic CDW order.

While CDW and SC compete spatially, they both serve to open an electronic gap. At zero magnetic field, we have a uniform SC state with a large gap of the order of $U$. With an applied field, we induce vortices with CDW correlations, leading to a spatially textured $SO(3)$ order parameter field. Initially, for small magnetic fields, the order parameter gradients reduce the electronic gap. Once CDW order percolates throughout the system, the CDW order parameter no longer suffers sharp gradients and strengthens the gap once again.

Discussion:

Working in the context of the attractive Hubbard model, we have shown that competing CDW order emerges in vortex cores. At large fields, vortex cores overlap leading to a coexistence phase with SC and CDW correlations spanning the system – a lattice version of a supersolid. The existence of a supersolid has been heavily debated in the context of liquid HeBoninsegni and Prokof’ev (2012). Our study suggests that superconductors with competing phases are strong candidates for supersolidity.

Our mean field theory in the strong coupling limit can also be seen as a spin problem with chiral interactions introduced by the orbital magnetic field. The vortex lattice ground state can be interpreted as a ‘meron crystal’ – a pseudospin state with a chiral textureVolovik and Krusius (2012). When CDW order spans the system, this state spontaneously breaks a $\mathbb{Z}_{2}$ symmetry corresponding to two possible checkerboard density patterns, or equivalently to the choice between $z$ and $-z$ ordering of the spins. Our estimate for the superfluid stiffness (see Supplementary Material) indicates that this state is stable to fluctuations at intermediate fields. At large magnetic fields, the CDW order becomes much stronger than SC. In this regime, fluctuations may destabilize SC while leaving the CDW order intact. This suggests a ‘vortex liquid’ phase with remnant CDW order and vanishing superfluid stiffness. This is an exciting direction for future study. In fact, field-induced coexistence and a charge-ordered vortex-liquid state have been reported in YBCOLeBoeuf et al. (2013); Gerber et al. (2015); Yu et al. (2016).

We present a schematic phase diagram for the Hubbard model in the $t^{\prime}$-$\alpha$ plane in Fig. 5. The plot shows the $t^{\prime}<t_{c}^{{}^{\prime}}$ region, where $t_{c}^{\prime}\sim 0.5t$ is the critical value beyond which CDW order is energetically unfavourable at all magnetic fields. This phase diagram provides a reference point for understanding phase competition in the cuprates. For example, we find that SC is lost at $H_{c2}$ via a first order transition to the competing CDW phase; this may have been seen in YBCO using thermal conductivity measurementsGrissonnanche et al. (2014). We see that vortices themselves show interesting ordering phenomena with the vortex lattice changing from triangular to square configuration. This resembles the suggestion of two vortex solid states in YBCO from torque magnetometry measurementsYu et al. (2016).

While our study focusses on a simple model Hamiltonian, our results broadly apply to several material families which host competing orders. We have used $t^{\prime}$ as a convenient handle to tune phase competition, this role could be played by experimentally tunable parameters such as doping in the cuprates Chang et al. (2012b), pressure in TiSe2Kusmartseva et al. (2009), etc. Our results provide a theoretical paradigm to understand phase competition in these systems.

Acknowledgments: We thank Arun Paramekanti, David Hawthorn, David Broun and Anton Burkov for discussions. MK acknowledges the use of High performance computing facility at Harish Chandra Research Institute, Allahabad, India.

I SUPPLEMENTARY MATERIAL

I.1 Peierl’s substitution

We have two types of hopping terms in our Hamiltonian: nearest-neighour hopping with amplitude $t$ and next-nearest hopping with amplitude $t^{\prime}$. The phase of each hopping element represents a line integral of the vector potential, in accordance with the principle of Peierl’s substitution. As the vector potential is not uniquely defined, there are several possible ways to assign the complex phases. The gauge invariant quantity is the magnetic flux: the sum of the hopping-phases along closed loops on the lattice.

In our periodic $L\times L$ lattice, we assign hopping-phases so as to obtain a uniform magnetic flux. We use the scheme shown in Fig. S1. We have introduced a parameter $\phi$ which encodes the phase picked up by an electron when hopping around any square plaquette, i.e., the magnetic flux through each square plaquette is $\hbar\phi/e$. Our square lattice system with periodic boundary conditions is equivalent to a torus, a closed surface. As shown by Dirac, the magnetic flux through any closed surface must be quantized in units of $h/e$ so that $\phi=2\alpha\pi/L^{2}$, with $\alpha$ being an integer. The parameter $\alpha$ determines the total flux through the system, given by $\alpha h/e$.

Dirac further argued that in the presence of a non-zero flux, we cannot define electronic wavefunctions on the surface smoothly. The phase of the wavefunctions must wind around a singularity, which is called the ‘Dirac string’. In our scheme, the Dirac string passes through the red region within the top-right square plaquette in Fig. S1. The sum of the hopping-phases around a contour which encloses this region has an additional contribution of $L^{2}\phi=2\pi\alpha$. This does not indicate an increased magnetic flux through this region; rather, it reflects a singularity in the definition of the electronic wavefunctions. Indeed, as this phase is a multiple of $2\pi$, it does not lead to any observable consequences.

I.2 Bogoliubov-de-Gennes formalism

The Hamiltonian for the Hubbard model is given in the main text. We decompose the on-site interaction term in pairing and density channels via a mean field decomposition. The complex superconducting order parameter is defined as, $\Delta_{i}=U\langle c_{i\downarrow}c_{i\uparrow}\rangle$, while the charge order parameter is defined as, $\phi_{i}=\frac{U}{2}(n_{i\uparrow}+n_{i\downarrow})=\frac{U}{2}(\langle c_{i\uparrow}^{\dagger}c_{i\uparrow}\rangle+\langle c_{i\downarrow}^{\dagger}c_{i\downarrow}\rangle)$. The resulting effective Hamiltonian is given by

[TABLE]

The hopping phases $\theta_{ij}$ and $\chi_{ij}$ are assigned according to the Peierl’s substitution scheme described above. We diagonalize this Hamiltonian using a Bogoliubov-Valatin transformation given by $c_{i\sigma}=\sum_{m}(u_{mi\sigma}\gamma_{m\sigma}-s_{\sigma}v_{mi\sigma}^{*}\gamma_{m,-\sigma}^{\dagger})$, where $\gamma_{m\sigma}^{\dagger}$ ($\gamma_{m\sigma}$) creates (annihilates) a quasiparticle with spin $\sigma$ with energy $\epsilon_{m}^{\sigma}$ and wavefunctions $u_{mi\sigma}$ and $v_{mi\sigma}$. We have introduced a spin index $s_{\uparrow}=1$ and $s_{\downarrow}=-1$. The resulting gap and number equations are

[TABLE]

where $f(\epsilon_{m})=1(0)$ if $\epsilon_{m}<0(>0)$ is the Fermi function at zero temperature. Starting from initial guess values, we iterate these equations to obtain self-consistent values of $\Delta_{i}$ and $\phi_{i}$. The chemical potential is tuned to fix the density at half-filling.

I.3 Phase competition

For a given set of parameters $t^{\prime}$, $U$ and $\alpha$, we find several self-consistent solutions. In particular, we find a pure CDW state with $\Delta_{i}=0$. To illustrate the phase competition in the Hubbard model, we compare the energy of this CDW state with that of the SC state in Fig. S2. In the absence of a magnetic field ($\alpha=0$), the two states are degenerate when $t^{\prime}=0$, while a non-zero $t^{\prime}$ lowers the energy of the SC state.

Competition between orders can also be tuned by increasing the magnetic field. This is depicted in Fig. S2(right). When $\alpha$ is increased at fixed $t^{\prime}$, the energy of the CDW state does not change as the CDW state is an insulator. However, in the SC phase, increasing $\alpha$ introduces more vortices and increases the energy. The energy of the SC state steadily rises and eventually crosses the CDW energy. This signals $H_{c2}$ at the mean-field level, with CDW order becoming energetically favourable over a superconducting vortex state.

I.4 Higher $t^{\prime}$ regime

With increasing $t^{\prime}$, the radius of the CDW region in each vortex core shrinks. Consequently, the threshold field for coexistence increases. Fig. S3 shows spatial maps of the SC and CDW order parameters for $t^{\prime}=0.4t$ which can be compared with the $t^{\prime}=0.2t$ data in Fig. 3 of the main text. The percolation of CDW correlations is slow to occur with a coexistence state only setting at $\alpha\sim 11$. Finally, $H_{c2}$ is encountered at $\alpha\sim 48$, when the CDW state becomes energetically favourable. Unlike the case of $t^{\prime}=0.2$ discussed in the main text, it is difficult to discern changes in the geometry of the vortex lattice here due to the high density of vortices.

I.5 Superfluid stiffness

Our mean-field results indicate coexistence of SC and CDW orders, forming a supersolid state. To check if this phase is stable to fluctuations, the standard diagnostic is superfluid stiffness – which measures the energy cost of imposing a smooth gradient in the SC order parameter. To estimate the stiffness, we take the following route. We introduce an additional component of the vector potential $\mathbf{A}_{tangential}=2\pi\hat{x}/L$. If we were to turn off the orbital magnetic field, this vector potential leads to a ‘flowing’ superfluid solution with $\Delta_{i}\sim\Delta_{0}e^{i4\pi x_{i}/L}$. The resulting energy cost is a measure of superfluid stiffness. For a simple superfluid with no competing order, this energy cost (the increase in energy per site) is proportional to $\rho/L^{2}$, where $\rho$ is the superfluid stiffness and $L$ is the system size. We define $\eta(L)=E_{2\pi}-E_{0}$, where $E_{2\pi}$ is the energy (per site) of the flowing state. This is calculated by using the $\Delta_{i}$ values obtained after inclusion of the tangential vector potential to evaluate the expectation value of the mean-field Hamiltonian.

The obtained $\eta(L)$ values are plotted as a function of $1/L^{2}$ in Fig. S4. In the $\alpha=0$ case (no orbital magnetic field), we see that $\eta(L)$ indeed scales as $1/L^{2}$, with a positive slope. This slope is proportional to the superfluid stiffness.

In the presence of the orbital field, we seek to plot $\eta(L)$ for configurations with the same magnetic flux density across different system sizes. In our calculations, the magnetic flux density is $\alpha h/eL^{2}$, with $\alpha$ being an integer and $L$ ranging from $18-30$ (for smaller sizes, we see strong finite size effects). Generically, it is not possible to find multiple $\{\alpha,L\}$ values for which $\alpha/L^{2}$ is a constant. In Fig. S4, we plot $\eta(L)$ for $\{\alpha,L\}=(\{2,18\},\{4,24\},\{6,30\})$ and $(\{4,18\},\{8,24\},\{12,30\})$ which correspond to approximately constant $\alpha/L^{2}$ values. The resulting $\eta(L)$ values also scale linearly with $1/L^{2}$ with a positive slope. We conclude that these superfluid stiffness is positive for these flux densities. We also note that the stiffness decreases with increasing flux density. In particular, we note that the stiffness is positive for the flux densities corresponding to $\alpha=4$ and $\alpha=8$ on a 24$\times$24 lattice. As discussed in the main text, these parameters have macroscopic CDW order with CDW correlations spanning the entire system. This suggests that the coexistence phase is stable to fluctuations.

It is possible that that stiffness may vanish at higher flux densities, perhaps close to $H_{c2}$. In this regime, the CDW order parameter becomes approximately constant while the SC order parameter suffers large gradients due to the presence of vortices. It is then conceivable that fluctuations can wash out the in-plane order while preserving order in the $z$-direction. This will lead to a ‘pairing liquid’ state (analogous to a spin liquid) with remnant CDW ordering. Equivalently, this can be understood as melting of the vortex lattice. Below this threshold, the inter-vortex distances are fixed by strong interactions between vortices. A small amount of disorder will then suffice to pin the entire vortex lattice so as to generate a robust SC state. However, when fluctuations wash out the coherence in the SC, the vortices become mobile giving rise to a ‘vortex liquid’. This is an interesting direction for future study.