The attractive Hubbard model as an $SO(3)$ system of competing phases: supersolid order and its thermal melting

Competition between superconductivity and charge order is a recurring theme in contemporary condensed matter physics. This is quintessentially captured in the attractive Hubbard model, a simple theoretical model where the competition can be directly tuned. In previous studies by the current authors, it has been suggested that the Hubbard model maps to an $SO(3)$ non-linear sigma model, where the phase competition becomes manifest. In this article, we rigorously demonstrate this mapping and use it to study thermal disordering of a supersolid. Starting with the attractive Hubbard model in the presence of an orbital field, we take the limit of strong coupling where a pseudospin description emerges. The in-plane pseudospin components represent superconducting pairing while the out-of-plane component encodes charge density wave order. We obtain an effective spin-$1/2$ Hamiltonian with ferromagnetic in-plane couplings and antiferromagnetic z-z couplings. In addition, the orbital field gives rise to a textured Dzyaloshinskii-Moriya interaction that has the same periodicity as the magnetic unit cell. In order to examine the nature of ordering in this spin model, we consider it in the classical limit. We assume slowly varying fields, leading to the $SO(3)$ non-linear sigma model description. As an application of these ideas, we study the nature of ordering using simulated annealing and classical Monte Carlo simulations. The ground state represents a supersolid with coexisting superconductivity and charge order. It can be viewed as a `meron crystal', a regular arrangement of superconducting vortices with charge-ordered cores. The coherent overlap of core regions gives rise to coherent long-ranged charge order. As the temperature is raised, this charge order is lost via a sharp phase transition in the Ising universality class.


I. INTRODUCTION
Experiments on the underdoped cuprates have fuelled renewed interest in phase competition, with superconductivity competing with charge density wave (CDW) order 1- 16 .
Studies have highlighted several manifestations of phase competition including ordered vortex cores 11 , coexistence 5 , strong impurity response 16 , non-monotonic evolution of critical fields 9 , etc. Manifestations of phase competition have also been seen in other material families such as transition metal dichalcogenides [17][18][19][20][21][22] , pnictides 23 and Ba 1−x K x BiO 3 24-26 . However, the physics in these materials is obscured by complications such as disorder and incommensurate ordering vectors. We require a simple model system where the effects of phase competition can be first understood in a clean setting. The attractive Hubbard model fits this requirement as a simple system that is amenable to various theoretical approaches.
A particularly interesting consequence of phase competition occurs in the presence of an orbital field. The superconducting order parameter forms vortices with superconductivity suppressed within the core region of each vortex. This allows for the competing CDW to arise locally. When the vortex density is large, the overlap between neighbouring vortex cores leads to coherent longranged CDW order 7,27 . This leads to a 'supersolid' phase that has coexisting superconductivity and CDW orders. Supersolidity has remained elusive in experiments despite intense studies in several contexts 28 . Overlap of ordered vortex cores provides a new mechanism that could allow for robust and verifiable supersolidity. In previous studies by the current authors, this mechanism has been demonstrated in the attractive Hubbard model using a meanfield approach 27 . Remarkably, the results were consistent with an SO(3) field theory for competing phases. Based on this observation, it was conjectured that these two models were equivalent. This conjecture was supported by further studies 29,30 . In this article, we establish this equivalence by way of a rigorous mapping. In the process, we find interesting results concerning a 'meron crystal', its stability to thermal fluctuations and melting.
One of our key results is the derivation of an SO (3) non-linear sigma model. This has strong similarities to the well-known SO(5) theory proposed in the context of the cuprates. The SO(5) theory is written in terms of a five component vector field: two corresponding to superconductivity and three to antiferromagnetism 31 . Although several consequences were worked out 32,33 , the model has not received support from experiments on the cuprates. In addition, it remains a phenomenological construct as no microscopic origin has been demonstrated. Here, we present a similar, but simpler, theory with a three-component order parameter. We start with a precise microscopic model and derive an effective SO (3) theory. This is potentially directly testable in experiments, with several proposals for realizing the attractive Hubbard model in ultracold atomic gases 34 .

II. THE ATTRACTIVE HUBBARD MODEL AT STRONG COUPLING
We consider particles on a square lattice with nearestand next-nearest-neighbour hopping, t and t . When two particles are on the same site, they lower the energy of the system due to an attractive interaction. Due to Pauli exclusion, this is only possible if they carry opposite spin. This leads to the Hamiltonian We have introduced Peierls' phases in the hopping amplitudes, θ ij , that originate from a uniform orbital magnetic field. This can be viewed as an Aharonov-Bohm phase accrued by the particle as it hops from one site to the next. We have introduced a chemical potential µ to fix particle density. In the rest of this article, we will restrict our attention to half-filling and use t as a handle to tune phase competition. The competition between superconductivity and CDW orders arises from a special SU (2) symmetry, first pointed out by C. N. Yang [35][36][37] . It requires three conditions: (a) a bipartite lattice with hopping between opposite sublattices, (b) absence of the orbital field, with all θ ij = 0, and (c) the density being fixed at half-filling. The square lattice with purely nearest neighbour hopping (t = 0) meets these requirements. At half-filling, it has perfect degeneracy between superconductivity and CDW orders. Upon introducing a next-nearest neighbour t hopping, this degeneracy is lost with a superconducting ground state. Nevertheless, CDW order remains as low-lying competitor with the energy cost scaling as ∼ t 238 .
The SU (2) symmetry is best seen in the strong coupling limit of the Hubbard model (U t). Previous studies have shown that the Hubbard model maps to a S = 1/2 pseudospin XXZ model 39 . This can be further mapped to a Heisenberg model using a sublatticedependent spin rotation. With a non-zero t , the CDW state manifests as a low-lying 'roton' excitation in the spin wave spectrum 40,41 . In this article, we derive the pseudospin Hamiltonian in the presence of an orbital field. Going further, we derive a coarse-grained field theory from the spin model. We study the role of thermal fluctuations by investigating the pseudospin model using Monte Carlo simulations.

III. STRONG COUPLING PSEUDOSPIN MODEL
We consider the strong coupling limit of the model with U t, t , following the superexchange scheme that has been presented in Refs. 39,40. If we only keep this dominant U -term in the Hamiltonian, the sites decouple from one another, leaving a purely on-site problem. The spectrum for the single site problem is shown in Fig. 1. The energy of the singly occupied states is −µ + U/4 as can be seen from the Hamiltonian above. The energy of the empty state is −U/4. Likewise, the energy of the doubly occupied state is −2µ − U/4.
At half-filling, the empty and doubly occupied sites must have the same energy so that they are occupied with the same probability. To ensure this, we set µ = 0. The spectrum splits into two pairs of states as shown in Fig. 1. The empty and doubly occupied states have lower energy, while the singly occupied states have higher energy. The energy difference between the pairs of states is U/2. The hopping terms in the Hamiltonian act as small perturbations on these states. Their effect is seen at second order in perturbation theory where they couple two sites at a time. To see this explicitly, we consider a two-site problem next.

A. Two site problem
We consider two sites labelled A and B. They may represent nearest neighbours or next-nearest neighbours on the square lattice. Apart from the dominant on-site terms, the Hamiltonian contains inter-site hopping terms, The hopping amplitude t AB can be complex with its phase given by the Peierl's substitution scheme. We reexpress it as t AB = τ AB e iθ AB . We now consider the low energy Hilbert space of the two-site problem. We introduce a pseudospin notation for the low energy states on a given site. We denote the empty state as pseudospin-down (| ⇓ ) and the doubly occupied state as pseudospin-up (| ⇑ ). In the two-site Hilbert space, we have four low energy states with each site having pseudospin-up or -down. As the hopping term takes us out of this subspace, we treat it within perturbation theory. Indeed, there are second order processes that connect low energy states, as shown in Fig. 2. In each path in the figure, the intermediate state has two singly occupied states. As a result, it has an energy cost given by 2 × U/2 = U .
The two-site states with parallel pseudospins (both empty or both doubly occupied) are unaffected within second order. In states with antiparallel spins, we find two processes: one that preserves pseudospins and one that exchanges them. We obtain the following Hamiltonian where FIG. 2: Superexchange pathways: The two panels depict processes that involve two sites -one doubly occupied and one empty. In the top panel, the initial and final states are the same, i.e., the pseudospin at each site is preserved. In the bottom panel, the initially empty state becomes doubly occupied and vice versa. This represents an exchange of pseudospins.

The 4 × 4 Hamiltonian matrix is given by
AB e −2iθ AB /U . The diagonal terms have a contribution from second order perturbation theory, in the form of C. In this term, the two hopping processes contribute with opposite phases that cancel out. In contrast, the phases add in the off-diagonal D term, imbuing it with a phase of 2θ AB .
We now add a constant shift of U/2 + τ 2 AB /U along the diagonals. The resulting Hamiltonian can be expressed in terms of an effective exchange coupling, J = 4τ 2 AB /U , This matrix has a simple interpretation in terms of spin operators. It can be written as whereŜ are pseudospin-1/2 operators. This can be rewritten as follows, The term proportional to sin(2θ AB ) can be expressed 43 . The term proportional to cos(2θ AB ) represents an XY-like exchange coupling between in-plane components. Note that the coupling constant, J, depends on the hopping strength on the bond. For example, it will have different strengths along nearest and next-nearest bonds.

B. Pseudospin model on the lattice
We have defined a pseudospin operator on each site. Its z-component represents the local CDW order parameter. To see this, we note that a site with pseudospin-up is doubly occupied with positive deviation from half-filling, whereas a site with pseudospin-down is empty with negative deviation. A state with maximal CDW order corresponds to an alternating arrangement of empty and doubly-occupied sites. This corresponds to an antiferromagnetic pseudospin arrangement with moments pointing alternately along ±ẑ. On the other hand, the inplane pseudospin components represent superconductivity. More precisely, the x and y components represent the real and imaginary parts of the pairing order parameter. This can be seen from the SU (2) pseudospin operators,

Superconductivity is signalled by non-zero expectation values for these operators.
Extending the two-particle effective Hamiltonian to the lattice, we arrive at a square lattice spin problem with the Hamiltonian The coupling strengths are given by J = 4t 2 /U and J = 4t 2 /U on nearest and next-nearest neighbours respectively. The bond-dependent exchange and Dzyalonshinskii-Moriya coefficients are given γ ij = cos(2θ ij ) and δ ij = sin(2θ ij ). The latter two depend on θ ij , the Peierls' phase associated with the bond (ij).
In the initial Hubbard model, the Peierls' phases encode a uniform orbital magnetic field. They are given by θ ij = (e/ ) j i A · dl, where A is the magnetic vector potential. Several studies have explored ways to realize this physical setup in ultracold atomic gases [44][45][46][47][48][49] . In a superconductor, strictly speaking, the orbital field must be self-consistently determined using Maxwell's equations. For the sake of simplicity, we assume a uniform orbital magnetic field below. This is a reasonable assumption in strongly type-II superconductors. The results discussed in Sec. IV below hold regardless of this assumption.
As the vector potential is not unique, neither is the assignment of Peierls' phases. If the vector potential is altered by a gauge transformation, this can be absorbed into the in-plane spin components by a suitable redefinition. This can be seen from Eq. 8, where in-plane pseudospin components couple to the θ ij 's while the z components do not. This is consistent with the identification of the in-plane components with the superconducting order parameter. In this sense, the effective model of Eq. 8 should not be thought of as a true spin model, as the inplane spin components are gauge-dependent quantities.
Traditionally, spin models are studied on finite lattices using periodic boundary conditions. Taking such an approach to Eq. 8 leads to some fundamental issues. We first note that the Peierls' phases necessarily contain singularities. To see this, we note that the square lattice forms a closed surface (a torus) due to periodic boundary conditions. A net flux through the lattice corresponds to having a magnetic monopole charge inside the torus. As argued by Dirac 50 , the vector potential cannot be smoothly defined on a surface enclosing a magnetic monopole. It necessarily includes flux tubes, called Dirac strings, that impart an Aharanov-Bohm phase of 2π. The number of Dirac strings is equal to the number of flux quanta that pierce the lattice. It follows that the Peierls' phases (θ ij 's) cannot have the same periodicity of the underlying lattice. They must necessarily form a large unit cell. The smallest possible unit cell corresponds to the area that contains a single Dirac string, i.e., the area carrying a single flux quantum. In other words, it is the 'magnetic unit cell'. One such phase assignment is shown in Fig. 3. This leads to the Hamiltonian in Eq. 8 with translational symmetry such that the unit cell is the same as the magnetic unit cell.
Using periodic boundaries has a second important consequence. Considering a charged particle on a surface enclosing a magnetic monopole, Dirac showed that its wavefunction cannot be defined in a smooth manner 50,51 .
In the system at hand, the superconducting order parameter cannot be smoothly defined on the torus. This can be seen as a consequence of having a non-zero number of vortices and no compensating anti-vortices. In the spin model, the in-plane spin components will not vary smoothly on the square lattice. They will invariably contain singularities or jumps. This serves as an additional caveat in viewing Eq. 8 as a spin problem.

IV. THE SO(3) EFFECTIVE FIELD THEORY
In the previous section, we arrived at an effective pseudospin description, assuming half-filling and U t, t . We now show that this pseudo-spin problem gives rise to a non-linear sigma model in the low energy limit.
We begin with the pseudospin Hamiltonian of Eq. 8 on an infinite square lattice. Promoting the spins to the classical limit, we have a lattice problem with threedimensional vector moments. We make two further assumptions: (a) at low energies (low temperatures), the spin configurations are 'smooth' with small gradients, and (b) with a weak orbital field, the Peierls' phase on each bond is small. We now note that Eq. 8 has antiferromagnetic z − z couplings between nearest neighbours. In contrast, the in-plane couplings are ferromagnetic (for small θ ij 's). This indicates that, in low-energy configurations, the spins are of the form, S × (∆ x , ∆ y , (−1) r ρ), where ∆ x , ∆ y and ρ are slowly varying quantities satisfying ∆ 2 x + ∆ 2 y + ρ 2 = 1. Here, S denotes the spin length. We henceforth set S = 1 for simplicity. The z-component carries a rapid oscillation given by (−1) r , which varies in a checkerboard fashion on the square lattice. As we expect ∆ x , ∆ y and ρ to vary smoothly on the scale of the lattice constant, we elevate them to slowly varying fields, ∆ x (r), ∆ y (r) and ρ(r) respectively. Note that the spatially-averaged z-moment vanishes. This corresponds to the assumption of half-filling, as the z-moment represents the local deviation from half-filling. We now calculate the contribution from each term in Eq. 8 within the language of coarse-grained fields.

A. CDW terms
We first consider the z − z couplings in Eq. 8 that resemble those of an Ising model on the square lattice, where (m, n) represents a site on the square lattice. The contribution from each site is given by The factors of 1/2 have been added to avoid double counting. We reinterpret this energy density in terms of the coarse-grained ρ(r) field. We useŜ z m,n ≈ (−1) m+n ρ(r m,n ), where (m, n) denotes a site of the square lattice. As with the standard Ising model, we elevate the summation over (m, n) to an integral and reexpress the integrand using ρ(r) and its derivatives. We obtain where a ρ = 2{J − J } and χ ρ = (J−2J ) 2

2
. Here, denotes the lattice constant of the square lattice. In these two coefficients, the J and J appear with opposite sign. This stems from the rapidly oscillating (−1) m+n factor that takes the opposite (same) sign on (next-) nearest neighbours. In addition, their relative amplitudes are different in a ρ and χ ρ , i.e., we have a ρ ∼ (J − J ) while χ ρ ∼ (J − 2J ). This difference arises from the differing bond lengths for nearest ( ) and next-nearest ( √ 2 ) neighbours.

B. Superconducting terms
We now consider the in-plane terms in the pseudospin Hamiltonian. In order to get a better understanding, we first take the vector potential to be zero, i.e., we ignore the Peierls' phases. This leads to a two-component spin model on the square lattice with ferromagnetic XY couplings. We have where S i, = (S x i , S y i ). Taking the in-plane components to be described by the slowly-varying fields ∆ x (r) and ∆ y (r), we obtain the field theory of an XY ferromagnet, where ∆ ≡ (∆ x , ∆ y ), a ∆ = 2{J + J } and χ ∆ = (J+2J ) 2 2 . This can be seen in direct analogy with the CDW term above, by replacing ρ with ∆ x/y . Unlike the CDW terms, the J and J contributions have the same sign here.
We now draw an analogy to the problem of a free particle in two-dimensional space. We take its wavefunction to be ∆(r) = ∆ x (r) + i∆ y (r). Taking its mass to be 1/2χ ∆ and assuming a constant potential (−a ∆ ), its Hamiltonian is given byĤ = χ ∆p 2 + a ∆ . The expectation value of the Hamiltonian is then precisely given by Eq. 13. A discrete form of this Hamiltonian can be constructed using a tight-binding-like approach. Discretizing the space as a square mesh with sites denoted by (m, n), we take ∆(r) → S x m,n + iS y m,n . This leads to the Hamiltonian in Eq. 12. This analogy provides a simple interpretation for in-plane terms in the Hamiltonian: the superconducting order parameter represents the wavefunction of a free particle (the Cooper pair).
We now introduce an orbital magnetic field. By comparing the Eqs. 12, 8 and 6, we see that the orbital field enters as Peierls' phases in a tight binding Hamiltonian. The superconducting wavefunction couples to the vector potential as a charged particle with charge 2e. It can immediately be deduced that the vector potential enters Eq. 12 via the well known minimal coupling prescription, where D ≡ ∇ + i 2e A. Note that the charge here is 2e, that of a Cooper pair. Indeed, we find the same result by a systematic analysis of the in-plane terms. The orbital field modifies Eq. 12 to give where δ and η sum over the nearest and next-nearest neighbour vectors respectively. Here, A's denote Peierls' phases, e.g., A m,n,δ = 2e (m+δx,n+δy) (m,n) A · dl. Assuming slow variations in the ∆'s and small values of the Peierls' angles, we precisely recover Eq. 14. This follows the usual derivation of the long-wavelength minimalcoupling Hamiltonian from a tight binding model with Peierls' phases.

2
. Rescaling allows us to write a simpler form, (17) where = 2J /(J + J ), χ = a 2 4 {1 + /2} and ξ = 4J /(J + 2J ). If t is small in the microscopic Hubbard problem, we have , ξ ∼ (t /t) 2 with ξ ≈ 2 . Here, and ξ reflect the anisotropy between superconductivity and CDW order. The SO(3) character of this model can be seen by setting t and the orbital field to zero. In this limit, Eq. 17 reduces to the Hamiltonian density of a symmetric Heisenberg ferromagnet. When weak anisotropies are introduced, the physics retains signatures of the proximate SO(3) point.
This form is closely related to the previously conjectured model in Ref. 27, where the anisotropy in the gradient term was ignored (i.e., ξ was set to zero). Nevertheless, this does not lead to any qualitative change in the physics of phase competition. We see this below in the nature of the ground state.

V. SIMULATING THE NON-LINEAR SIGMA MODEL
We have shown that the attractive Hubbard model reduces to an SO(3) non-linear sigma model. Using this equivalence, we seek to study its physics in the presence of an orbital field. The energy of the system is given by the Hamiltonian density of Eq. 17. The ground state can be found by minimizing the energy, subject to the uniform length constraint (ρ 2 + |∆| 2 = 1). However, minimizing Eq. 17 on the infinite two-dimensional plane is a non-trivial task. Likewise, thermal properties of the non-linear sigma model can be found by averaging over configurations with a suitable Boltzmann weight. Once again, this is a difficult task on the infinite plane.
We approach this problem by reversing the arguments put forward in the previous sections. We now view Eq. 8, the pseudospin model on the square lattice, as a regularization of the non-linear sigma model in Eq. 17. We will study the pseudospin model on finite lattices with periodic boundary conditions and look for results that remain consistent upon increasing system size. This opens the door to well established techniques from the field of magnetism. In particular, we use simulated annealing to find the ground state of Eq. 8. We will interpret the result in terms of the smooth fields of the non-linear sigma model. We will then study the role of thermal fluctuations using classical Monte Carlo simulations.
The pseudospin model of Eq. 8 is defined on the square lattice. As explained in Sec. III B above, the Hamiltonian depends on the choice of the Peierls' phases. We present results using the scheme depicted in Fig. 3. We assume a 12 × 12 magnetic unit cell so that the Peierls' phases do not vary too rapidly from one bond to the next. We consider a lattice composed of an n × n array of magnetic unit cells, giving rise to a 12n × 12n lattice with periodic boundaries. We approach the thermodynamic limit by increasing n.

A. Supersolidity in the ground state
To find the lowest energy state, we perform simulated annealing of the pseudospin model. We use two types of single-site moves: Metropolis and microcanonical (overrelaxation). At each site, we find the effective field that arises from the neighbouring moments. The Metropolis move corresponds to changing the inclination with respect to the effective field. The microcanonical move rotates the spin about the effective field so as to preserve the energy.
The lowest energy state found from simulated annealing is shown in Fig. 4. We have used a 24 × 24 lattice containing four 12×12 magnetic unit cells. The net magnetic flux through the lattice thus corresponds to four flux quanta. At each site, we interpret the in-plane components of the pseudospin as the superconducting order parameter. From Fig. 4(left), we see that the superconducting amplitude vanishes at regularly spaced points, indicating a vortex lattice. The number of vortices is eight, with two vortices for each flux quantum. We have defined a flux quantum with respect to the charge e of the particle hopping on the lattice. As a Cooper pair has charge 2e, we find two vortices for each flux quantum.
The competition with CDW order is clearly seen in Fig. 4(centre) which shows the z-component of the spins in the ground state. We see strong CDW order appearing in each vortex core. The CDW order percolates through the inter-vortex space and covers the entire lattice. This leads to a 'meron crystal' as shown in Fig. 4(right). Here, we plot (S x m,n , S y m,n , (−1) m+n S z m,n ) vs. (m, n), i.e., position on the lattice. This conveys the variation of the pseudospin orientation in space. We have removed a rapidly oscillating phase in the z-component of the pseudospin (see discussion in Sec. IV A above). Each superconducting vortex takes the form of a 'meron' in the pseudospin. The in-plane components wind by 2π as we move around the vortex. Within the core region, an out-of-plane component develops to preserve the spin length. Due to overlap between adjacent merons, the out-of-plane component is non-zero everywhere. It has the same sign at all sites, indicating coherent CDW order.
This picture is consistent with the results of Ref. 27 where the Hubbard model was directly studied using Bogoliubov-deGennes mean field simulations. In particular, the low energy state here represents a 'supersolid'. It has well-defined superconducting order that is reflected in the formation of a vortex lattice. At the same time, it has long-ranged CDW order.

B. Classical Monte Carlo simulations
We have established that ground state of the non-linear sigma model in Eq. 17 is a supersolid with coexisting superconductivity and CDW order. The superconductivity sector encapsulates an additional layer of ordering in the form of a vortex lattice with discrete translational symmetry. Upon increasing the temperature, we may see multiple phase transitions where these orders melt independently. To study thermal fluctuations, we study the pseudospin model of Eq. 8 using classical Monte Carlo simulations. We use single-site Metropolis and microcanonical (overrelaxation) moves. We start from a random initial configuration on an L × L lattice at high temperature and progressively decrease the temperature. At each temperature, we perform 8 ×10 6 sweeps, each with L 2 single-site moves, with the ratio of Metropolis to microcanonical fixed at 4:3. The first 2 × 10 5 moves are discarded to allow for equilibration.
We first discuss the thermal evolution of the CDW order. In the non-linear sigma model of Eq. 17, the CDW order parameter shows an Ising-like character with the energy being invariant under ρ(r) → −ρ(r). This originates from the Hubbard model where the CDW order represents a checkerboard-like modulation in density. The Ising degree of freedom corresponds to choosing one of the two sublattices as that with higher density. In the pseudospin model of Eq. 8, the CDW order parameter is the staggered z-magnetization, given by M = 1 The temperature dependence of the order parameter is shown in Fig. 5(a). We find a profile that is typical of an Ising magnet. Starting from zero at high temperatures, it approaches a non-zero value at low temperatures. Unlike the standard Ising magnet, the magnetization in the zero-temperature-limit is not unity. This can be understood from the ground state configuration in Fig. 4. The CDW order is not uniform; rather, it has maximal intensity at vortex cores and weak order at inter-vortex positions. Nevertheless, we see a clear indication of an Ising-like phase transition. Fig. 5(a) shows the order parameter for various system sizes with L = 12n, where n = 2, 3, 4, 5. We choose L to be multiples of 12 so that we can construct the pseudospin Hamiltonian using a 12 × 12 magnetic unit cell. The flux density is the same for all system sizes. We find further evidence for a phase transition in the form of a peak in susceptibility as shown in Fig. 5(b). The peak height grows with system size as expected.
To determine the precise location of the CDW phase transition, we examine the Binder cumulant for various system sizes, shown in Fig. 5(c). We find a crossing at T c ≈ 0.145 ± 0.001. We surmise that this transition belongs to the universality class of the 2D Ising model. To verify this, we perform a scaling analysis of the data. In Fig. 5(d-f), we plot the rescaled order parameter, susceptibility and Binder cumulant vs. reduced temperature (using T c as obtained from the Binder cumulant cross-ing). We find good scaling collapse using the well known critical exponents of the 2D Ising model 52 , viz., ν = 1, β = 1/8 and γ = 7/4. Based on this finding, we assert that CDW order vanishes via a continuous phase transition in the 2D Ising universality class.
We next discuss thermal evolution of the superconducting order. We do not find a distinct phase transition within our Monte Carlo scheme. We believe this is due to technical limitations, as discussed below. Nevertheless, a qualitative understanding can be gained by examining typical configurations extracted from the Monte Carlo simulations, shown in Fig. 6. At low temperatures, we see a vortex lattice, albeit with small distortions. The distortions increase with increasing temperature. Beyond T ∼ 0.03, the vortex lattice is lost as some vortices come close to one another and essentially fuse. At this point, we may view the system as being deep inside a vortex liquid phase.
We believe that a vortex melting transition occurs at T 0.005. However, this is not discernible in our simulations as the spins do not relax adequately at low temperatures. Indeed, we do not find a perfect vortex lattice even at the lowest temperatures as some distortions persist (see Fig. 4). This could be a consequence of our single-site update scheme. At a more subtle level, this could be a consequence of the gauge structure. Our system with periodic boundaries cannot support a smoothly varying superconducting field. As discussed in Sec. III B above, it must necessarily contain singularities or jumps. On account of these discontinuities, a single-site update scheme may not be able to explore the space of all low energy configurations.

VI. DISCUSSION
We have presented a study of phase competition in the attractive Hubbard model at strong coupling. We demonstrate a mapping to a pseudospin problem and further onto an SO(3) field theory. This brings out phase competition as an inherent feature of this model. It also reveals an interesting role for an orbital magnetic field as it induces vortices in the superconducting order, but with CDW-ordered cores. Indeed, we find a supersolid ground state with phase coexistence arising from vortexcore-overlap. In the language of spins, we find a meron crystal -an emergent crystalline phase with a mesoscopic lattice scale, analogous to the well known skyrmion crystal phase. With increasing temperature, superconductivity and CDW orders melt independently with a sharp Ising phase transition in the CDW sector.
Our results bear similarities with disordering transitions in other systems with coexisting orders. We mention two examples from the field of magnetism: (a) The square J 1 − J 2 antiferromagnet with J 2 > J 1 /2 breaks O(3) ⊗ Z 2 symmetry, where the Z 2 character corresponds to a choice between vertical and horizontal stripes 53 . While the O(3) rotational symmetry is restored at an infinitesimal temperature, the Z 2 order persists up until a critical temperature where it is lost via an Ising transition. (b) The triangular lattice XY antiferromagnet breaks Z 2 ⊗ U (1) symmetry in the ground state, where the Z 2 character corresponds to a local chirality degree of freedom. The Z 2 order is lost via an Ising transition 54 . In the context of the attractive Hubbard model, we have presented an effective field theory for competing orders. This could be used to potentially develop a renormalization group scheme to understand the physics of disordering. For example, the Ising transition temperature can be lowered by increasing t , i.e., the energy cost of the CDW phase. At a critical value of t , the Ising transition will compete with the vortex lattice melting transition. This can potentially give rise to an interesting combined melting transition.
The pseudospin model derived in Sec. III B above is essentially a quantum model with S = 1/2 moments. We have studied this model in the classical limit, taking into account thermal fluctuations. An interesting future direction is to investigate the role of quantum fluctuations. In analogy with the thermal state immediately below the Ising transition, quantum fluctuations may disrupt superconductivity while preserving CDW order. Such a state would represent a 'pairing liquid' in analogy with a spin liquid. The pairing liquid offers two advantages over typical spin liquid models: (i) it has an additional tuning handle in the form of an orbital field, (ii) fluctuations of the pairing liquid are intrinsically coupled to the CDW order parameter due to the non-linear uniform length constraint. This offers a new route to probe fluctuations in the liquid phase. These issues may be explored within a quantum treatment of the pseudospin model.