text:
# Ring exchange, the Bose metal, and bosonization in two dimensions

## I Introduction

Despite intense experimental effort exploring the phase diagram of the cuprates, the nature of the pseudo-gap regime remains mysterious. In the very underdoped normal state there are strong experimental hints of local pairing and superconducting correlationspseudo1 ; pseudo2 ; pseudo3 , despite the absence of phase coherent superconductivity. Within this picturephasefluc , the pseudo-gap regime supports a pair field with non-zero amplitude but with strong phase fluctuations which disrupt the superconductivity. A theoretical approach then requires disordering the superconductivity by unbinding and proliferating vortices. Unfortunately, such an approach will likely face a most worrisome dilemma. Proliferation and condensation of single $`hc/2e`$ vortices necessarily leads to a confined insulatorBFN2 , which should show sharp features in the electron spectral function in apparent conflict with ARPES experimentspseudo3 . On the other hand, if pairs of $`hc/2e`$ vortices condenseBFN2 ; Z2A , the pseudo-gap phase must necessarily support gapped unpaired vortex excitationsvison1 ; vison2 ; vison3  (“visons”)Z2B ; Z2C , which have yet to be observedvisexp1 ; visexp2 . Is there an alternative possibility? In this paper we find and explore a truly novel quantum fluid phase of 2d bosons, which we refer to as a “Bose metal”, which answers this question in the affirmative. The Bose metal phase supports gapless charge excitations, but is not superconducting \- the Cooper pairs are not condensed. The $`hc/2e`$ vortices are likewise gapless and uncondensed. Being a stable quantum fluid with no broken spatial or internal symmetries, the Bose metal phase has many properties reminiscent of a 2d Fermi liquid - a 1d locus of gapless excitations (a “Bose surface”), a finite compressibility, an almost $`T`$linear specific heat with logarithmic corrections.

Neutron scattering measurements in the undoped cuprates have revealed a zone boundary magnon dispersion which can best be accounted for by presuming the presence of appreciable four-spin exchange processesAeppli . Recent theoretical work points to the importance of such ring exchange processes driving fractionalization and spin-charge separationZ2A ; ringnum ; BGF . Motivated by these considerations, we focus on a class of model Hamiltonians describing a square lattice of 2d Bosons (Cooper pairs) with appreciable ring-exchange. It is hoped that such microscopic models will be appropriate in describing the charge sector in the under-doped cuprates. When the ring exchange processes are sufficiently strong, we find that the Bose metal phase is the stable quantum ground state.

Lattice models of interacting bosons in two dimensions, such as the boson Hubbard model, have been studied quite exhaustively during the past several decadessupins ; subirbook , primarily as models for Josephson junction arrays and superconducting films but also in the context of quantum magnets with an easy plane $`U(1)`$ symmetry and most recently in the context of trapped Bose condensates moving in an optical latticeBEC . Typically, the ground state phase diagram consists of a superconducting phase and one or more insulating statessupins . At fractional boson densities commensurate with the lattice, say at half-filling, the insulating behavior is driven by a spontaneous breaking of translational symmetry. When accessed from the superconductor, such insulating states can be fruitfully viewed as a condensation of elementary vorticesMPAFvort . Very recent work on a Kagome lattice boson model with ring exchangeBGF  (which arises in the context of frustrated quantum magnets) has revealed the existence of a more exotic fully gapped insulating state with no broken symmetries whatsoever. In this state the charged excitations are “fractionalized”, carrying one-half of the bosons chargeSubirZ2 , and there are also gapped vortex excitations called visons. Here, we study the simpler system of bosons with ring exchange moving on a 2d square lattice. Our central finding is the existence of a critical gapless quantum phase of $`2d`$ bosons which is stable over a particular but generic parameter range. This phase is in some respects very similar to the familiar gapless Luttinger liquid phase of interacting bosons and Fermions moving in one spatial dimensionLutt1 ; Lutt2 ; Lutt3 . Indeed, to describe this phase we introduce a 2d generalization of bosonization. Specifically, we introduce a new duality mapping which transforms a square lattice model of 2d bosons (in a rotor representation) with ring exchange, into a 2d theory of “vortices” hopping on the sites of the 2d dual lattice. Then, based on the symmetries we construct a low energy effective description, in terms the two (dual) phases of the boson and vortex creation operators. The Bose metal phase has a Gaussian fixed point description in terms of these fields. In addition, there are various non-quadratic interactions, which scale to zero in the Bose metal. Instabilities towards superconducting and insulating states are triggered when one or more of these interactions become relevant, closely analogous to $`2k_F`$ instabilities in a 1d Luttinger liquid.

Since this paper is quite long, the remainder of the introduction is devoted to a brief summary of the ring model Hamiltonian and the bosonized description of the Bose metal. The Hamiltonian we focus on describes 2d bosons in a rotor representation:

$$_{\mathrm{}}=\frac{U}{2}\underset{𝐫}{}(n_𝐫\overline{n})^2K\underset{𝐫}{}\mathrm{cos}(\mathrm{\Delta }_{xy}\varphi _𝐫),$$
(1)
with

$$\mathrm{\Delta }_{xy}\varphi _𝐫\varphi _𝐫\varphi _{𝐫+\widehat{𝐱}}\varphi _{𝐫+\widehat{𝐲}}+\varphi _{𝐫+\widehat{𝐱}+\widehat{𝐲}}.$$
(2)
Here $`𝐫`$ labels sites of the 2d square lattice and $`\varphi _𝐫`$ and $`n_𝐫`$ are canonically conjugate variables:

$$[\varphi _𝐫,n_𝐫^{}]=i\delta _{𝐫,𝐫^{}}.$$
(3)
Representing the phase of the boson wavefunction, $`\varphi _𝐫`$ is taken to be $`2\pi `$ periodic $`\varphi _𝐫=\varphi _𝐫+2\pi `$, so that the eigenvalues of the conjugate boson number operator $`n_𝐫`$ are integers. The mean boson density is set by $`\overline{n}`$, and we will primarily focus on the case with half-filling: $`\overline{n}=\frac{1}{2}`$. The argument of the cosine in the second term is a lattice second derivative, involving the four sites around each elementary square plaquette. This term hops two bosons on opposite corners of a plaquette, clockwise (or counterclockwise) around the plaquette, and is an XY analog of the more familiar $`SU(2)`$ invariant 4-site ring exchange term for spin one-half operators (which arises in the context of solid $`3He`$oldring1 ; oldring2 ).

In addition to the conventional spatial, particle-hole, and total number conservation symmetries (see Sec. III), this ring Hamiltonian has an infinite set of other symmetries. Specifically, the dynamics of the Hamiltonian conserves the number of bosons on each row and on each column of the 2d square lattice - a total of $`2L`$ associated symmetries for an $`L`$ by $`L`$ system. This is fewer than a gauge theory, which has an extensive number of local symmetries, but these symmetries will nevertheless play a crucial role in constraining the dynamics of the model and stabilizing a new phase.

We will also consider adding a near neighbor Boson hopping term,

$$_t=t\underset{𝐫,𝐫^{}}{}\mathrm{cos}(\varphi _𝐫\varphi _𝐫^{}),$$
(4)
where the summation is taken over near-neighbor sites on the square lattice. This term breaks the $`2L`$ U(1) symmetries corresponding to conserved boson numbers on each row and column, leaving only the globally conserved total boson number.

The remarkable result of this paper is that, over a particular but generic parameter range, ring models of this type sustain the novel Bose metal phase. The effective description of the Bose metal is in terms of low energy “coarse grained” variables $`\phi \varphi `$ and an additional field $`\vartheta `$ (see Sec. III), which is related via $`n\overline{n}\pi ^1\mathrm{\Delta }_{xy}\vartheta `$ to the long-wavelength density, and further can be used to construct a “vortex creation operator” $`ve^{i\vartheta }`$. To fix notation, we define the $`n_𝐫,\varphi _𝐫,\phi _𝐫`$ operators on the sites of the original square lattice, with $`x,y`$ coordinates taking half-integer values, and $`\vartheta _𝐫`$ (with $`N_𝐫,\theta _𝐫`$ operators to be defined later) on the dual square lattice with integer $`x,y`$ coordinates. The $`\phi ,\vartheta `$ fields are governed by an approximately Gaussian (Euclidean) effective action $`𝒮=_0^\beta 𝑑\tau `$, with $`=_0+_1`$, and $`\tau `$ is imaginary time, $`\beta =1/k_BT`$. The Gaussian part of the Lagrangian is

$$_0=_0[\phi ,\vartheta ]+\underset{𝐫}{}\frac{i}{\pi }_\tau \phi _𝐫^{}\mathrm{\Delta }_{xy}\vartheta _𝐫,$$
(5)
where the sum is over sites $`𝐫`$, $`𝐫^{}=𝐫+\widehat{x}/2+\widehat{y}/2`$ on the dual and original lattices, respectively, and the effective Hamiltonian is

$$_0=_𝐤\left[\frac{𝒦(𝐤)}{2}|(\mathrm{\Delta }_{xy}\phi )_𝐤|^2+\frac{𝒰(𝐤)}{2\pi ^2}|(\mathrm{\Delta }_{xy}\vartheta )_𝐤|^2\right].$$
(6)
Here the momentum ($`𝐤`$) integral is taken over the Brillouin zone $`|k_x|,|k_y|<\pi `$, and $`(\mathrm{\Delta }_{xy}\phi )_𝐤`$ denotes the Fourier transform of $`\mathrm{\Delta }_{xy}\phi _𝐫`$ (and similarly for $`\mathrm{\Delta }_{xy}\vartheta `$). The functions $`𝒦(𝐤),𝒰(𝐤)`$ are non-vanishing finite periodic, and analytic. Their values along the $`k_x`$ and $`k_y`$ axes parameterize the Bose liquid, much like the effective mass and Fermi liquid parameters do a Fermi liquid. Experts will note a strong similarity to the bosonized effective action for a one dimensional Luttinger liquidLutt2 ; Lutt3 , which is explored in Sec. III. Like in a Luttinger liquid, there are additional non-quadratic terms in the action. It is sufficient to keep only those of the form $`_1=_{1\phi }+_{1\vartheta }`$, with

$`_{1\vartheta }`$ $`=`$ $`{\displaystyle \underset{𝐫}{}}{\displaystyle \underset{q,𝐬}{}}t_q^{m,n}\mathrm{cos}[q(\phi _𝐫\phi _{𝐫+𝐬})],`$ (7)
with $`𝐫`$ on the original lattice, and

$`_{1\phi }`$ $`=`$ $`{\displaystyle \underset{𝐫}{}}\{{\displaystyle \underset{q=1}{\overset{\mathrm{}}{}}}\upsilon _{2q}\mathrm{cos}(2q\vartheta _𝐫+qQxy)`$ (8)
$`+{\displaystyle \underset{q𝐬}{}}w_{2q}^{m,n}\mathrm{cos}[2q(\vartheta _𝐫\vartheta _{𝐫+𝐬})qQ(nx+my+mn)]\},`$
with $`𝐫`$ on the dual lattice. In both Eqs. (7-8), $`Q=2\pi \overline{n}`$, $`𝐫=(x,y)`$, $`𝐬=(m,n)`$. Here, $`w_{2q}^{m,n}=w_{2q}^{n,m}`$ gives the amplitude to hop $`2q`$ vortices by the translation vector $`(m,n)`$ (or $`(n,m)`$), and $`t_q^{m,n}=t_q^{n,m}`$ gives the hopping amplitude for a $`q`$ bosons along the same vector. In the Bose metal phase, all these terms are “irrelevant”, i.e. give only perturbative corrections to physical properties, though these corrections can be significant. Also like in a Luttinger liquid, there is a non-trivial relation between microscopic quantities such as the energy and particle densities and the coarse grained variables of the low energy theory. In particular, for the energy and particle densities, one finds

$`\delta n_{x+\frac{1}{2},y+\frac{1}{2}}`$ $``$ $`c_0^\rho \mathrm{\Delta }_{xy}\vartheta _{xy}`$ (9)
$`+{\displaystyle \underset{q=1}{\overset{\mathrm{}}{}}}c_{2q}^\rho \mathrm{\Delta }_{xy}\mathrm{sin}(2q\vartheta _{xy}+qQxy),`$
$`\epsilon _{xy}`$ $``$ $`c_0^\epsilon \dot{\vartheta }_{xy}^2+{\displaystyle \underset{q=1}{\overset{\mathrm{}}{}}}c_{2q}^\epsilon \mathrm{cos}(2q\vartheta _{xy}+qQxy),`$ (10)
where $`\delta n=n\overline{n}`$, $`\epsilon `$ is the energy density, and $`c_{2q}^{\rho /\epsilon }`$ are non-universal constants.

The paper is organized as follows. In Section 2 we treat the ring exchange term within a “spin-wave” approximation. This reduces the Hamiltonian to a quadratic form which can be readily diagonalized. Within the two-dimensional Brillouin zone there are gapless excitations along the lines $`k_x=0`$ and $`k_y=0`$, which are associated with an infinite set of conservation laws possessed by the ring exchange model. Dual “vortex” variables are then introduced via a new plaquette duality transformation. This dual representation is well suited to numerical simulations since it is free of any “sign problems”. We implement a quantum Monte Carlo simulation, and show that the Bose metal phase is present over appreciable regions of the phase diagram of the ring-exchange rotor model.

In Section III we construct the low energy effective model in terms of the dual vortex fields, and extract the universal properties of the Bose metal phase in Section IV. These properties are confirmed by the quantum Monte Carlo results. Section V is devoted to an analysis of the stability of the Bose metal phase. Remarkably, we find that for generic (incommensurate) boson densities, there are regions of parameters where the Bose metal phase is stable towards all perturbations, even those that break the row/column symmetries. The situation is reminiscent of the Fermi liquid phase, whose “fixed point” description possesses an enormous set of “emergent symmetries” - the number of Fermions on each patch of the Fermi surface are independently conserved. At commensurate fillings we explore in some detail the instabilities of the Bose metal towards various insulating states.

In Section 6 we consider a “hard-core” version of the boson Hamiltonian which allows us to obtain an exact zero energy wavefunction when the couplings are carefully tuned. We then perturb away from the soluble point, and compute properties of the adjacent quantum phase which we thereby identify as the Bose metal. Finally, in Section 7 we conclude with a discussion of how various ring exchange processes emerge from a $`Z_2`$ gauge theory formulation of interacting electrons and then make some speculative remarks about the possible relevance of the Bose metal phase in the context of the under-doped cuprates.

## II  Spin Waves and Plaquette Duality

### II.1 Spin Waves, massless modes and symmetry

Here we consider a spin wave approximation to the microscopic ring Hamiltonian which leads to a harmonic and soluble theory. The resulting energy spectrum vanishes along two lines in momentum space. This remarkable feature is then shown to follow directly from the existence of the infinite set of conservation laws of the ring Hamiltonian. As we detail in later sections, such a harmonic description underlies an effective theory of the Bose-metal phase, in close analogy to bosonization of the 1d Luttinger liquid.

To this end, it is useful to employ a combination of path integral and Hamiltonian methods. A standard path integral representation is constructed in the usual way using $`\varphi `$ eigenstates. In the time continuum limit, the partition function for the pure ring Hamiltonian takes the form,

$$Z=\mathrm{Tr}e^{\beta H_{\mathrm{}}}=[d\varphi _{𝐫\tau }]e^{_0^\beta 𝑑\tau L^\varphi },$$
(11)
with the Lagrangian

$$L^\varphi =\underset{𝐫}{}[\frac{1}{2U}(_\tau \varphi _𝐫)^2+i\overline{n}_\tau \varphi _𝐫K\mathrm{cos}(\mathrm{\Delta }_{xy}\varphi _𝐫)].$$
(12)
The resemblance of the above Lagrangian with that for the standard Boson-Hubbard modelsupins  (which has a single, rather than double, lattice derivative inside the cosine) suggests that one might try expanding the cosine potential to quadratic order. Doing so gives a soluble Harmonic theory for the action, which can be readily diagonalized as,

$$𝒮_{\mathrm{spinwave}}=\frac{1}{2U}\frac{d^2k}{(2\pi )^2}_{\mathrm{}}^{\mathrm{}}\frac{d\omega }{2\pi }[\omega ^2+E_𝐤^2]|\varphi (𝐤,\omega )|^2,$$
(13)
with $`𝐤=(k_x,k_y)`$ living in the first Brillouin zone, $`k_x,k_y[\pi ,\pi ]`$. Remarkably, the energy dispersion,

$$E_𝐤=4\sqrt{UK}|\mathrm{sin}(k_x/2)\mathrm{sin}(k_y/2)|,$$
(14)
vanishes on both the $`k_x`$ and $`k_y`$ axis. The presence of these gapless excitations can be traced directly to the existence of the infinite set of symmetries of the ring Hamiltonian, which conserves the number of bosons on each row and on each column of the 2d square lattice. Specifically, these symmetries imply an invariance of the energy (and action) under,

$$\varphi _𝐫\varphi _𝐫+\mathrm{\Phi }_x(x)+\mathrm{\Phi }_y(y),$$
(15)
for arbitrary functions $`\mathrm{\Phi }_x(x)`$ and $`\mathrm{\Phi }_y(y)`$. This invariance dictates that in the Harmonic spin wave form the energy must vanish whenever $`k_x=0`$ or $`k_y=0`$.

Since the plaquette term involves a lattice second derivative rather than an ordinary lattice gradient, however, it is clear that phase ($`\varphi `$) fluctuations will be large in the ring model, and the spin wave expansion is suspect. Indeed, even in the classical limit $`U0`$, one can readily see that “vortex” configurations in which $`\varphi _𝐫`$ winds by $`2\pi `$ around some plaquette (e.g. for a vortex with center at $`x=y=0`$, $`\varphi _{xy}=\pi /2\mathrm{\Theta }(x)\mathrm{\Theta }(y)+\pi \mathrm{\Theta }(x)\mathrm{\Theta }(y)+3\pi /2\mathrm{\Theta }(x)\mathrm{\Theta }(y)`$, where $`\mathrm{\Theta }(z)`$ is the Heavyside step function) are finite in energy rather than logarithmic as in an ordinary XY model. Further, “double vortex” configurations in which this winding is $`4\pi `$ can be smoothly deformed into zero energy conformations (e.g. $`\varphi _{xy}=\pi \mathrm{\Theta }(x)+\pi \mathrm{\Theta }(y)`$). This suggests that for non-zero $`U`$ there will be substantial “vorticity” in the low-lying states.

To address the legitimacy of the spin wave expansion, it is necessary to account for the periodicity of the cosine potential. This is most readily accessed by transforming to a dual form, just as one transforms to a dual bosonized representation for a system of 1d particlesLutt3 . Here we introduce a new 2d quantum duality transformationoldduality  which is specifically tailored for the ring model.

### II.2 Plaquette Duality

We consider dual fields living on the sites of the dual lattice, denoted $`\theta _𝐫`$ and $`N_𝐫`$, which are canonically conjugate variables:

$$[N_𝐫,\theta _𝐫^{}]=i\delta _{𝐫,𝐫^{}}.$$
(16)
Here we take $`\pi N_𝐫`$ to be $`2\pi `$ periodic so that $`\theta _𝐫/\pi `$ has integer eigenvalues. The reason for this unusual choice of normalization should become apparent below. For notational ease we are denoting the sites of both the original and dual square lattices as $`𝐫`$. In analogy to Eq. (2), it is convenient to introduce an operator defined on the plaquettes of the dual lattice (i.e. sites of the original lattice) as:

$$\mathrm{\Delta }_{xy}\theta _𝐫\theta _𝐫\theta _{𝐫+\widehat{𝐱}}\theta _{𝐫+\widehat{𝐲}}+\theta _{𝐫+\widehat{𝐱}+\widehat{𝐲}}.$$
(17)
These two new dual fields are related to the original fields by the relations,

$$\pi N_𝐫^{}=\mathrm{\Delta }_{xy}\varphi _𝐫;\mathrm{\Delta }_{xy}\theta _𝐫=\pi n_𝐫^{},$$
(18)
where $`𝐫^{}=𝐫+\widehat{x}/2+\widehat{y}/2`$. One can check that with $`N_𝐫`$ and $`\theta _𝐫`$ as conjugate fields, the original variables satisfy the required commutation relations.

To interpret the new variables, consider now a vortex centered on some site $`(x,y)`$ of the dual lattice (the “core”). Classically, for the four sites on plaquette of the original lattice surrounding this site,

$$e^{i\varphi _{x+a/2,y+b/2}}=\left(\frac{a+ib}{\sqrt{2}}\right)^{N_v},$$
(19)
where $`a,b=\pm 1`$, and $`N_v`$ denotes the number of vortices (vorticity) on this plaquette. Comparing to Eq. (18), one finds $`N=N_v(\mathrm{mod}\mathrm{\hspace{0.17em}2})`$, so $`N`$ can be interpreted as a vortex number operator, modulo two. Since $`\theta `$ and $`N`$ are canonically conjugate, the operators $`e^{i\theta }`$ perform canonical “translations” of $`N`$, and hence can be regarded as vortex creation and annihilation operators, respectively. Note that since $`\theta \pi 𝒵`$, $`e^{2i\theta }=1`$, consistent with the ambiguity in $`N=N_v`$ under even integer shifts. Clearly periodicity of $`\varphi `$ in the original variables is encoded in the discreteness of $`\theta `$ in the dual description.

The plaquette Hamiltonian when re-expressed in the dual variables reads,

$$H_{\mathrm{}}=K\underset{𝐫}{}\mathrm{cos}(\pi N_𝐫)+\frac{U}{2\pi ^2}\underset{𝐫}{}[\mathrm{\Delta }_{xy}\theta _𝐫\pi \overline{n}]^2.$$
(20)
We will for the time being neglect the tunneling term $`H_t`$, which will be returned to later. Note the strong similarity between the dual and original forms. To bring this out more clearly, and for the numerical simulations to be considered shortly, it is useful to go to a path integral formulation. In particular, consider the partition function

$$Z_{\mathrm{}}=\mathrm{Tr}e^{\beta H_{\mathrm{}}}.$$
(21)
Expanding $`Z_{\mathrm{}}`$ in the usual Trotter fashion with a time-slice $`\mathrm{\Delta }\tau =ϵ0^+`$, using the discrete basis of eigenstates of $`\theta _𝐫`$, one finds (Appendix A)

$$Z=\underset{\{\theta _𝐫(\tau )\}}{}e^{\left(ϵ_\tau \right)L_{\mathrm{}}},$$
(22)
with the “Lagrangian”

$$L_{\mathrm{}}[\theta ]=\underset{𝐫}{}\left\{\frac{ϵ}{\pi ^2}\mathrm{ln}(\frac{2}{ϵK})(_\tau \theta _𝐫)^2+\frac{U}{2\pi ^2}[\mathrm{\Delta }_{xy}\theta _𝐫\pi \overline{n}]^2\right\},$$
(23)
where $`_\tau \theta _𝐫=[\theta _𝐫(\tau +ϵ)\theta _𝐫(\tau )]/ϵ`$. The $`ϵ`$-dependence of the time derivative term in Eq. (23) is familiar from the “time-continuum limit” relating e.g. the $`d+1`$-dimensional classical and $`d`$-dimensional quantum transverse-field Ising models. Note the strong similarity of Eq. (23) to Eq. (12), which emphasizes their nearly self-dual nature.

The formulation in Eqs. (22-23) is quite convenient for numerical simulations. For the simulations, we define an integer-valued “height” field $`\stackrel{~}{\theta }_𝐫(\tau )`$ through

$$\mathrm{\Delta }_{xy}\theta _𝐫(\tau )=\pi n_𝐫^{(B)}+\pi \mathrm{\Delta }_{xy}\stackrel{~}{\theta }_𝐫(\tau )$$
(24)
where $`n_𝐫^{(B)}=[1+(1)^{x+y}]/2`$ is the background staggered density, chosen such that $`\stackrel{~}{\theta }_𝐫(\tau )`$ obeys periodic boundary conditions. In this language, we obtain a generalized (anisotropic) solid-on-solid model in $`2+1`$ dimensions. The correlation functions for the $`\theta `$-fields are easily re-expressed (and numerically evaluated) in terms of these height variables, using Monte Carlo methods. An indication of the sort of results obtained in shown in Fig. 1 which presents the phase diagram of this model as a function of $`U/K`$ and $`ϵ`$, based on an evaluation of the density correlations. For the simulations, we worked on a $`L_x\times L_x\times L_\tau `$ lattice, with various system sizes indicated in the paper. We used a Metropolis algorithm with a single-site update, $`\stackrel{~}{\theta }_𝐫(\tau )\stackrel{~}{\theta }_𝐫(\tau )\pm 1`$, which corresponds to a ring-exchange move for the boson density. We checked for equilibration of various quantities, and averaged the data over $`10^410^6`$ sweeps of the lattice depending on the correlation function and location in the phase diagram. Notice that in the time continuum limit ($`ϵ0`$), the simulations reveal two phases as the dimensionless ratio $`K/U`$ is varied, separating an ordered charge density wave state at large $`U`$ from the Bose metal phase when the ring exchange term is large.

Before obtaining an effective low-energy description, it is convenient to rewrite the dual partition function, Eq. (22) in terms of continuous $`\theta `$ variables using the Poisson summation formula,

$$Z=[d\theta _𝐫(\tau )]\underset{s_𝐫(\tau )}{}e^{_𝐫(\tau )(2is_𝐫(\tau )\theta _𝐫(\tau )\delta |s_𝐫(\tau )|)}e^{ϵ_\tau L_{\mathrm{}}[\theta ]}.$$
(25)
In Eq. (25), we have included a parameter $`\delta 1`$, which “softens” the discrete-$`\theta `$ constraint (exactly imposed for $`\delta =0`$). Carrying out the sum over the integer-valued $`s_𝐫(\tau )`$ variables, one finds

$$Z=Z_0[d\theta _𝐫(\tau )]e^{S_0S_1},$$
(26)
where $`Z_0`$ is a constant (divergent for $`\delta 0`$),

$`S_0`$ $`=`$ $`{\displaystyle \underset{\tau }{}}ϵL_{\mathrm{}}[\theta ],`$ (27)
$`S_1`$ $`=`$ $`{\displaystyle \underset{𝐫\tau }{}}{\displaystyle \underset{q=1}{\overset{\mathrm{}}{}}}v_{2q}\mathrm{cos}2q\theta _𝐫(\tau )`$ (28)
and the $`v_{2q}`$ are O(1) coefficients whose precise values are not important. It is convenient to drop the constant and introduce infinitesimal source fields $`h^\rho `$, $`h^\epsilon `$, for the density and energy density, respectively,

$$Z[h^\rho ,h^\epsilon ]=[d\theta _𝐫(\tau )]e^{S_0S_1S_h},$$
(29)
with

$$S_h=\underset{𝐫,\tau }{}\left[h_𝐫^{}^\rho (\tau )\left(\mathrm{\Delta }_{xy}\theta _𝐫\pi \overline{n}\right)+h_𝐫^\epsilon (\tau )(_\tau \theta _𝐫)^2\right],$$
(30)
where $`𝐫^{}=𝐫+\widehat{x}/2+\widehat{y}/2`$.

## III Effective Model

We now make an educated but perhaps bold guess as to the nature (and existence!) of a low energy effective description. In the spirit of the renormalization group and effective field theory, we imagine defining a temporally coarse-grained variable $`\vartheta `$ in which the high-frequency modes of $`\theta `$ are averaged over

$$\vartheta _𝐫(\tau )=\left[\theta _𝐫(\tau )\right]_f\pi \overline{n}xy,$$
(31)
where the square brackets indicate an average over “fast” high-frequency modes, and we have for convenience removed the mean “curvature” in $`\vartheta `$ by a non-fluctuating spatially dependent shift. Provided that we average only over high-frequency modes, we expect that the resulting effective description of $`\vartheta `$ will remain local both in space and (imaginary) time. Inspection of Eq. (26) shows that the dual microscopic action for $`\theta `$ can be written as a sum of a quadratic part ($`_{\tau ,𝐫}ϵ_{\mathrm{}}`$) and non-quadratic corrections. We postulate that the low-energy effective action for $`\vartheta `$ is a sum of a renormalized quadratic term similar to $`_{\mathrm{}}`$ and small non-quadratic corrections. Although we have not implemented a detailed renormalization group treatment, the latter statement is tantamount to declaring the system is controlled by a stable Gaussian fixed point. Obviously such a fixed point can have at best a large finite basin of attraction, and we will return to the problem of determining the region of stability of the Gaussian theory in Sec. V.

### III.1 Symmetries

Formally, the above postulate of the effective field theory description can be formulated as the statement that the generating functional can be approximated for low-frequency source fields by $`Z[h^\rho ,h^\epsilon ]𝒵_0𝒵[h^\rho ,h^\epsilon ]`$, where $`𝒵_0`$ is an irrelevant constant, and

$$𝒵=[d\vartheta _𝐫(\tau )]e^{𝒮_0[\vartheta ]𝒮_1[\vartheta ]𝒮_h[\vartheta ]},$$
(32)
where $`𝒮_0[\vartheta ]`$ is a renormalized quadratic form, $`𝒮_1`$ contains small renormalized non-quadratic perturbations, and $`𝒮_h`$ is linear in $`h^\rho `$, $`h^\epsilon `$ (we drop higher order terms in the infinitesimal sources). In general, we expect all three parts of the action to take the most general form possible consistent with locality, analyticity (since only high-frequency modes are averaged over), and the symmetries of the problem, which should be found from the microscopic partition function, Eq. (29). These are

* Row/Column symmetry: $`\theta _{x,y}(\tau )\theta _{x,y}(\tau )+\pi [m_x(x)+m_y(y)]`$, where $`m_x(x)`$ and $`m_y(y)`$ are arbitrary integer-valued functions. This corresponds to conservation of the number of vortices (modulo 2) on each row and column, a dual version of the conserved boson number on each row/column of the original lattice.
* Space and time translations: $`\theta _𝐫(\tau )\theta _{𝐫+𝐑}(\tau +\tau _0)`$, where $`𝐑`$ is an arbitrary lattice vector and $`\tau _0`$ is an arbitrary real number. Under this translation, the sources must also be translated, $`h_𝐫^{\rho /\epsilon }(\tau )h_{𝐫+𝐑}^{\rho /\epsilon }(\tau +\tau _0)`$.
* Reflections across a row or column containing a site of the dual lattice (or bonds of the original lattice): $`\theta _{x,y}ss^{}\theta _{sx,s^{}y}`$, where $`s,s^{}=\pm 1`$. The sources transform as $`h_{x,y}^\rho ss^{}h_{sx,s^{}y}^\rho `$, $`h_{x,y}^\epsilon h_{sx,s^{}y}^\epsilon `$.
* Reflections across a row or column containing bonds of the dual lattice (or sites of the original lattice): This is not independent, and can be obtained from a composition of a translation and a site-reflection. But for completeness, $`\theta _{x,y}ss^{}\theta _{1/2+s(x1/2),1/2+s^{}(y1/2)}`$, $`h_{x,y}^\rho ss^{}h_{1/2+s(x1/2),1/2+s^{}(y1/2)}^\rho `$, $`h_{x,y}^\epsilon h_{1/2+s(x1/2),1/2+s^{}(y1/2)}^\epsilon `$.
* Time-reversal: $`\theta _𝐫(\tau )\theta _𝐫(\tau )`$, $`h_𝐫^{\rho /\epsilon }(\tau )h_𝐫^{\rho /\epsilon }(\tau )`$.
* Four-fold rotations: $`\theta _{x,y}\theta _{y,x}`$, $`h_{x+1/2,y+1/2}^\rho h_{y+1/2,x1/2}^\rho `$, $`h_{x,y}^\epsilon h_{y,x}^\epsilon `$.
* Particle-hole symmetry: This is an invariance only for integer or half-integer densities ($`2\overline{n}𝒵`$). For such values, the symmetry operation is $`\theta _{xy}2\pi \overline{n}xy\theta _{xy}`$, $`h_𝐫^\rho h_𝐫^\rho `$, $`h_𝐫^\epsilon h_𝐫^\epsilon `$.

### III.2 Gaussian action

The form of $`𝒮`$ is dictated by these symmetries and the relation between the microscopic and coarse-grained fields, Eq. (31). Note that the shift by $`\pi \overline{n}xy`$ implies that $`\vartheta `$ is not a scalar. Consider first $`𝒮_0`$. By space and time translational symmetry, it is diagonal in momentum and frequency space,

$$𝒮_0=\frac{1}{2}_𝐤_\omega (𝐤,\omega )|\vartheta (𝐤,\omega )|^2,$$
(33)
where $`_𝐤d^2𝐤/(2\pi )^2`$, $`_\omega _{\mathrm{}}^{\mathrm{}}𝑑\omega /(2\pi )`$, and the momentum integral is taken over the Brillouin zone $`|k_x|,|k_y|<\pi `$. By row/column translational symmetry,

$$(k_x=0,k_y,\omega =0)=0\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{all}k_y,$$
(34)
and similarly for $`k_xk_y`$. The latter condition implies the existence of gapless excitations along the $`k_x`$ and $`k_y`$ axes in momentum space. These are the only such gapless states mandated by symmetry, and we expect the low energy physics to be dominated therefore by momenta near these axes and by low frequencies. Quite generally, the kernel $``$ can be expanded near the zero modes at $`k_x=\omega =0`$, and at lowest order takes the form,

$$(𝐤,\omega )=A(k_y)\omega ^2+B(k_y)k_x^2,$$
(35)
where the expansion functions $`A(k_y)`$ and $`B(k_y)`$ are even and $`2\pi `$-periodic in $`k_y`$. A similar expansion is of course possible around $`k_y=\omega =0`$, with the identical expansion coefficients due to rotational symmetry. Analyticity at $`𝐤=0`$ further implies $`B(0)=0`$. A convenient representation of $``$ which satisfies these requirements is

$$(𝐤,\omega )=\frac{\omega ^2+E_𝐤^2}{\pi ^2𝒦(𝐤)},$$
(36)
with mode energy,

$$E_𝐤4\sqrt{𝒦(𝐤)𝒰(𝐤)}|\mathrm{sin}(k_x/2)\mathrm{sin}(k_y/2)|.$$
(37)
Here, $`𝒰(𝐤)`$ and $`𝒦(𝐤)`$ are positive $`2\pi `$ periodic functions which characterize the Bose metal phase, but it is only their behavior for $`k_x1`$ or $`k_y1`$ that determine the universal low-energy properties of the theory.

### III.3 Interaction terms

Consider next the interaction terms, $`𝒮_1=𝑑\tau _1`$. Locality requires that they couple combinations of $`\vartheta _𝐫`$ fields only with nearby points $`𝐫`$. Hence we consider a successive sequence of terms, $`_1=_m_{1;m}`$, coupling $`\vartheta `$ fields at a total of $`m`$ distinct points. We expect that the $`_{1;m}`$ becomes (exponentially) increasingly small with increasing $`m`$. First consider the single-site terms, which are highly constrained, particularly by translational invariance, under which $`\vartheta _{xy}+\pi \overline{n}xy`$ transforms as a scalar. Generally,

$$_{1;1}=\underset{xy}{}\left[\underset{q=1}{\overset{\mathrm{}}{}}\upsilon _{2q}\mathrm{cos}(2q\vartheta +qQxy)+\underset{q=2}{\overset{\mathrm{}}{}}𝒦_q^1(_\tau \vartheta )^{2q}\right],$$
(38)
where $`Q=2\pi \overline{n}`$, and has the physical meaning of the smallest reciprocal lattice vector of a one dimensional lattice with density $`\overline{n}`$ (like $`2k_F`$ for a charge density wave). The second set of terms involve more time-derivatives than the analogous quadratic interactions in $`𝒮_0`$, and hence are negligible at low frequency, so we will take $`𝒦_q^1=0`$ in the following. The $`\upsilon _{2q}`$ terms will play an important role in Sec. V. Note that, except when the density takes special commensurate values (such as the interesting case $`\overline{n}=1/2`$), they are strongly oscillatory in space.

Next consider two-site terms. Similar to the $`K_q^1`$ terms in Eq. (38), a variety of spatial (lattice) and time derivative invariants are possible, but are negligible relative to $`𝒮_0`$, so we do not include them here. More interesting are sinusoidal terms, which must take the form

$`_{1;2}`$ $`=`$ $`{\displaystyle \underset{xy;m,n}{}}{\displaystyle \underset{qq^{}}{}}\{w_{2q,2q^{},mn}^{}\mathrm{cos}[2q\vartheta _{xy}2q^{}\vartheta _{x+m,y+n}+Q(qxyq^{}(x+m)(y+n))]`$ (39)
$`+w_{2q,2q^{},mn}^+\mathrm{cos}[2q\vartheta _{xy}+2q^{}\vartheta _{x+m,y+n}+Q(qxy+q^{}(x+m)(y+n))]\},`$
where rotational and reflection invariance require $`w_{2q,2q^{},nm}^\pm =w_{2q,2q^{},mn}^\pm `$, and without loss of generality we may also impose $`w_{2q^{},2q,nm}^\pm =w_{2q,2q^{},mn}^\pm `$. Like the $`\upsilon _{2q}`$ terms in Eq. (38), most of the operators in Eq. (39) are highly oscillatory for generic $`\overline{n}`$. An important exception arises in $`w_{2q,2q,mn}^{}`$ for rational densities $`\overline{n}=z^{}/z`$, where $`z,z^{}`$ are integers. Then these terms are non-oscillatory if $`z^{}mq`$ and $`z^{}nq`$ are $`z`$ times integers. We consider in particular the case $`z^{}=1`$, for which $`m,q`$ and $`n,q`$ can be chosen as all possible factors of $`z\mathrm{}`$ into two integers, with arbitrary integer $`\mathrm{}`$. Keeping only these terms, one has

$`_{1;2}`$ $`=`$ $`{\displaystyle \underset{xy}{}}{\displaystyle \underset{mq=nq=z\mathrm{}}{}}w_{2q}^{m,n}\mathrm{cos}[2q(\vartheta _{xy}\vartheta _{x+m,y+n})],`$ (40)
where, $`w_{2q}^{m,n}=w_{2q}^{n,m}`$ by rotational invariance. Since $`e^{i\vartheta }`$ acts as a vortex creation operator, the terms in Eq. (40) can be thought of as hopping $`2q`$ vortices a distance $`𝐬=(m,n)`$ on the 2d dual lattice.

Finally, consider $`𝒮_h=𝑑\tau _h`$. As for $`_1`$, we expect that the largest contributions will occur in terms involving $`\vartheta _𝐫`$ fields at a small number of distinct points $`𝐫`$. For simplicity, we keep only the most local of these. Employing the symmetries above, one finds

$`_h`$ $`=`$ $`{\displaystyle \underset{xy}{}}\{h_{x+1/2,y+1/2}^\rho [c_0^\rho \mathrm{\Delta }_{xy}\vartheta _{xy}`$ (41)
$`+{\displaystyle \underset{q=1}{\overset{\mathrm{}}{}}}c_{2q}^\rho \mathrm{\Delta }_{xy}\mathrm{sin}(2q\vartheta _{xy}+qQxy)]`$
$`+h_{xy}^\epsilon [c_0^\epsilon \dot{\vartheta }_{xy}^2+{\displaystyle \underset{q=1}{\overset{\mathrm{}}{}}}c_{2q}^\epsilon \mathrm{cos}(2q\vartheta _{xy}+qQxy)]\},`$
where $`c_{2q}^{\rho /ϵ}`$ are constants determined by the high-energy physics. Here we have kept only the single-plaquette terms in the charge density, and the single-site terms in the energy density. Note that the operators multiplying $`h^\rho `$ involve sines rather than cosines, which is a consequence of particle-hole and reflection symmetries. The form in Eq. (41) can be derived perturbatively in $`v_{2q}`$ using an explicit one-step coarse-graining procedureunpub . The quantitative reliability of this constructive procedure is, however, limited to small bare $`v_{2q}`$, so for the problem of interest the constants $`c_{2q}^{\rho /\epsilon }`$ should be viewed as phenomenological parameters. Eq. (41) should be interpreted analogously to the non-trivial bosonization rules for e.g. the charge density and other operators in one dimension. Very generally, a low energy effective action must be accompanied by rules specifying the operator content of various physical fields. In analogy to the notations used there, one may rewrite Eq. (41) as a scaling equality between the microscopic and coarse-grained fields, as in Eq. (9).

### III.4 Self-dual form - Bosonization analogy

Having established the form of the effective action for the $`\vartheta `$ variables, we are also interested in being able to compute quantities involving the boson phase $`\phi `$. Since the effective theory is perturbative in $`𝒮_1`$, it is straightforward to re-introduce $`\phi `$ using Gaussian integration. In particular, the quadratic action $`𝒮_0`$ in Eq. (33) with the kernel in Eq. (36) can be transformed into the form shown in Eq. (5) in the introduction using a Hubbard-Stratonovich transformation:

$`e^{\frac{1}{2}_{𝐤\tau }\frac{1}{\pi ^2𝒦(𝐤)}|_\tau \vartheta (𝐤,\tau )|^2}`$ $`=`$ (42)
$`{\displaystyle [d\phi ]e^{\frac{1}{2}_{𝐤\tau }𝒦(𝐤)|(\mathrm{\Delta }_{xy}\phi )_𝐤|^2}e^{_\tau _𝐫\frac{i}{\pi }_\tau \vartheta _𝐫^{}\mathrm{\Delta }_{xy}\phi _𝐫}},`$
where $`𝐫^{}=𝐫+\widehat{x}/2+\widehat{y}/2`$ as usual, and $`(\mathrm{\Delta }_{xy}\phi )_𝐤`$ denotes the Fourier transform of $`\mathrm{\Delta }_{xy}\phi _𝐫`$. After an integration by parts in the last term on the right hand side above, the full Gaussian part of the action takes a particularly transparent form, $`S_0=_\tau _0`$, with Lagrangian $`_0`$ defined in Eqs. (5,6).

The partition function is now represented as a path integral over both sets of “low energy” fields, $`\phi _𝐫`$ and $`\vartheta _𝐫`$, in an appealing self-dual form. The Gaussian theory above gives a fixed point description of the Bose metal phase. When augmented by the operator content of the fields, as in Eq. (9), together with the irrelevant non-linear interaction terms in Eqs. (38,39), it gives a complete description of the universal properties of the Bose metal phase, as detailed in the next Section.

It will sometimes be convenient to integrate out the dual field, $`\vartheta `$, leaving a Gaussian action just in terms of the phase of the boson wavefunction:

$$𝒮_0=\frac{1}{2}_𝐤_\omega _\phi (𝐤,\omega )|\phi (𝐤,\omega )|^2,$$
(43)
with Kernel,

$$_\phi (𝐤,\omega )=\frac{\omega ^2+E_𝐤^2}{𝒰(𝐤)},$$
(44)
and with $`E_𝐤`$ as given in Eq. (37). Notice that this form is identical to the Gaussian fixed point theory in terms of $`\vartheta `$ in Eq. (33), except that $`𝒰(𝐤)`$ has replaced $`\pi ^2𝒦(𝐤)`$ in the denominator of Eq. (44).

The above self-dual representation makes particularly apparent the close analogy between this theory and Bosonization theory in one-dimensionLutt3 . In particular, underlying both theories is a pair of dual scaler fields, one the phase of the wavefunction and the other related to the density. The Luttinger liquid fixed point is a Gaussian theory in the two fields, and has a form which is nearly identical to the Bose metal action above:

$$_{1d}=_{1d}+\frac{i}{\pi }_x_\tau \varphi _x\theta ,$$
(45)
with Gaussian Hamiltonian:

$$_{1d}=_x[\frac{K}{2}(_x\varphi )^2+\frac{U}{2\pi ^2}(_x\theta )^2].$$
(46)
Notice that the long-wavelength particle density in Bosonization theory, $`n_x\theta /\pi `$, has simply been replaced by a (lattice) second derivative in our $`2+1`$-dimensional theory: $`n\mathrm{\Delta }_{xy}\vartheta /\pi `$. The commutation relations between the two dual fields has also been modified, as is apparent by inspecting the Berry’s phase term involving both fields in the above Lagrangians. The expressions in Eq. (9) relating the bare boson density and energy density to the low-energy fields are a generalization of the more familiar Bosonization expressions, where for example the Boson density near $`2k_F`$ is proportional to $`\mathrm{cos}(2\theta )`$.

A key strength of Bosonization in 1d is that it allows one to study the instabilities of the Luttinger liquid towards various types of ordered phases, such as a charge density wave state. Similarly, the above Gaussian representation of the Bose metal fixed point is particularly suitable for studying instabilities both towards insulating states with broken translational symmetry and towards a superfluid. But before undertaking this analysis, we explore the Bose metal phase in some detail.

## IV The Bose Metal phase

We now turn to the physical properties of the Bose metal phase, using the effective low energy theory developed in the previous section. For the present, we will assume stability of the Bose metal, so that all physical quantities are accessible via perturbation theory (in $`𝒮_1`$) around the Gaussian theory for $`\phi ,\vartheta `$. The validity of this assumption is discussed in Sec. V.

### IV.1 Thermodynamic properties

Consider first the thermodynamic properties of the Bose metal. The simplest is the compressibility $`\kappa =n/\mu `$. This is trivially given from the Gaussian theory, and one finds $`\kappa =1/𝒰(𝐤=\mathrm{𝟎})`$. Numerically, the compressibility may be obtained from an extrapolation of the density-density correlation function,

$$\chi _{nn}(𝐤,\omega _n)=\underset{𝐫}{}𝑑\tau \frac{1}{\pi ^2}\mathrm{\Delta }_{xy}\theta _𝐫(\tau )\mathrm{\Delta }_{xy}\theta _\mathrm{𝟎}(0)e^{i\omega _n\tau i𝐤𝐱},$$
(47)
from which $`\kappa =\chi _{nn}(𝐤0,\omega =0)`$. The compressibility $`\kappa `$ determined in this way from the simulations were used to define the phase diagram in Fig. (1)

Also interesting is the specific heat. Since the Bose metal is essentially a free boson theory, this is determined completely from the density of states for the collective free boson modes

$$\rho (E)=_𝐤\delta (EE_𝐤),$$
(48)
where the boson energy is

$$E_𝐤=\mathrm{\Omega }(𝐤)4|\mathrm{sin}\frac{k_x}{2}\mathrm{sin}\frac{k_y}{2}|,$$
(49)
and $`\mathrm{\Omega }(𝐤)=\sqrt{𝒰(𝐤)𝒦(𝐤)}`$. Like an ordinary (Fermi liquid) metal, the density of states gets a large finite contribution from the low energy modes near the “Bose surface” along the coordinate axes in momentum space. Unlike a Fermi liquid, however, there is a weak logarithmic divergence associated with the crossing point at $`k_x=k_y=0`$. In particular,

$$\rho (E)\frac{1}{\pi ^2\mathrm{\Omega }_0}\left[\mathrm{ln}(1/ϵ)+𝒞_0+O(ϵ)\right],$$
(50)
where $`\mathrm{\Omega }_0=\mathrm{\Omega }(0,0)`$, $`ϵ=E/\mathrm{\Omega }_0`$ and

$$𝒞_0=4\mathrm{ln}2+_0^\pi \frac{dk}{\mathrm{sin}k/2}\left(\frac{\mathrm{\Omega }_0}{\mathrm{\Omega }(k,0)}1\right).$$
(51)
Note that the constant term in the density of states depends upon the full form of the dispersion $`\mathrm{\Omega }(k,0)=\mathrm{\Omega }(0,k)`$ all along the Bose surface, while the logarithmic term depends only upon the behavior at $`𝐤=\mathrm{𝟎}`$. The specific heat is then

$$C_v=\frac{T}{4}_0^{\mathrm{}}𝑑x\frac{x^2\rho (Tx)}{\mathrm{sinh}^2x/2},$$
(52)
which to the same accuracy gives

$$C_v\frac{T}{4\pi ^2\mathrm{\Omega }_0}\left[I_0\mathrm{ln}\frac{\mathrm{\Omega }_0}{T}+𝒞_0I_0+I_1\right],$$
(53)
and $`I_n=_0^{\mathrm{}}𝑑xx^2(\mathrm{ln}(1/x))^n/\mathrm{sinh}^2(x/2)`$, or $`I_0=4\pi ^2/3`$, $`I_1=2(32\gamma _E)\pi ^2/3+8\zeta ^{}(2)4.643`$, where $`\zeta (z)`$ is the Riemann zeta function, and $`\gamma _E0.5772`$ is Euler’s constant.

### IV.2 Boson correlation functions

It is interesting to contrast the large density of states for collective excitations with the tunneling density of states for the original bosons, which we will see is strongly suppressed. In particular, consider

$$G_\varphi (𝐫,\tau )=e^{i\varphi _𝐫(\tau )}e^{i\varphi _\mathrm{𝟎}(0)}.$$
(54)
By symmetry, since the boson number is conserved on each row and column, $`G_\varphi (𝐫,\tau )`$ can be non-zero only for $`𝐫=\mathrm{𝟎}`$. To further determine the behavior of $`G_\varphi `$ in the Bose metal, we must relate the microscopic Bose creation/annihilation operators to the low energy modes. By symmetry,

$$e^{i\varphi _𝐫}e^{i\phi _𝐫}\left[1+A\mathrm{cos}(2\vartheta _𝐫+2\pi \overline{n}xy)+\mathrm{}\right].$$
(55)
We expect the low energy properties to be dominated by the first term, i.e. simply replacing $`\varphi \phi `$. Since $`\phi `$ is a Gaussian variable,

$$G_\varphi (𝐫,\tau )e^{F^\varphi (\tau )}\delta _{𝐫,\mathrm{𝟎}},$$
(56)
with

$$F^\varphi (\tau )=\frac{1}{2}(\phi _\mathrm{𝟎}(\tau )\phi _\mathrm{𝟎}(0))^2=_𝐤\frac{𝒰(𝐤)}{2E_𝐤}\left(1e^{E_𝐤\tau }\right).$$
(57)
Since $`E_𝐤`$ vanishes linearly both as $`k_x0`$ and as $`k_y0`$, inspection of Eq. (57) shows that $`F^\varphi (\tau )`$ has a double logarithmic divergence as $`\tau \mathrm{}`$. Indeed, a careful calculation shows that

$$F^\varphi (\tau )\frac{F_2^\varphi }{2}|\mathrm{ln}\mathrm{\Omega }_0\tau |^2+F_1^\varphi \mathrm{ln}(\mathrm{\Omega }_0\tau )+\mathrm{}$$
(58)
for $`\mathrm{\Omega }_0\tau 1`$, with

$`F_2^\varphi `$ $`=`$ $`{\displaystyle \frac{1}{2\pi ^2}}\sqrt{{\displaystyle \frac{𝒰_0}{𝒦_0}}},`$ (59)
$`F_1^\varphi `$ $`=`$ $`{\displaystyle \frac{1}{\pi ^2}}{\displaystyle _0^\pi }𝑑k\left[\sqrt{{\displaystyle \frac{𝒰(k,0)}{𝒦(k,0)}}}{\displaystyle \frac{1}{2\mathrm{sin}k/2}}\sqrt{{\displaystyle \frac{𝒰_0}{𝒦_0}}}{\displaystyle \frac{1}{k}}\right]`$ (60)
$`{\displaystyle \frac{1}{2\pi ^2}}\sqrt{{\displaystyle \frac{𝒰_0}{𝒦_0}}}\left[2\mathrm{ln}\pi +\gamma _E\right],`$
where $`𝒰_0=𝒰(0,0)`$, $`𝒦_0=𝒦(0,0)`$. The $`\mathrm{ln}^2\mathrm{\Omega }_0\tau `$ behavior of $`F^\varphi (\tau )`$ at large $`\tau `$ implies that $`G_\varphi (\tau )`$ decays faster than any power law. This translates into a similar singular behavior for the tunneling density of states. Writing the boson Green’s function in a spectral representation,

$$G_\varphi (\mathrm{𝟎},\tau )=_0^{\mathrm{}}𝑑E\rho _{\mathrm{tun}}(E)e^{E|\tau |},$$
(61)
one finds using a simple saddle-point analysis that the above behavior at large $`\tau `$ requires

$$\rho _{\mathrm{tun}}(E)\mathrm{exp}\left\{\frac{\alpha }{2}\mathrm{ln}^2\frac{\mathrm{\Omega }_0}{E}\mathrm{ln}\frac{\mathrm{\Omega }_0}{E}(\beta \mathrm{ln}\mathrm{ln}\frac{\mathrm{\Omega }_0}{E}+\lambda )\right\},$$
(62)
where $`\alpha =\beta =F_2^\varphi `$, $`\lambda =F_1^\varphi 1F_2^\varphi +\mathrm{ln}F_2^\varphi `$.

Thus, although the Bose metal possesses a large set of gapless modes, the density of states for adding a boson into the system vanishes at low energy. Again, this is analogous to a Luttinger liquid in 1DLutt3 ; the conservation of particles per row or column means that the added particle affects all the particles in a particular row or column, leading to a suppressed amplitude for such tunneling events. Remarkably, the tunneling density of states actually vanishes faster than in a Luttinger liquid, indeed faster than any power law in energy. This behavior can be roughly understood as arising from two orthogonality catastrophes occurring simultaneously in the row and column in which the boson is added or removed.

It is also instructive to consider the boson four-point correlation function, for simplicity at equal times. Due to the row/column symmetries this is non-vanishing only when the four points sit at the corners of a rectangle:

$`G_\varphi ^{(4)}(x,y)`$ $`=`$ $`e^{i\varphi _{oo}(\tau )}e^{i\varphi _{xo}(\tau )}e^{i\varphi _{oy}(\tau )}e^{i\varphi _{xy}(\tau )}.`$ (63)
Upon replacing $`\varphi `$ with the low energy field $`\phi `$, this can be readily evaluated using the Gaussian Bose metal action. For any fixed $`y`$, one finds power-law behavior in $`x`$ (and vice versa since $`G_\varphi ^{(4)}(y,x)=G_\varphi ^{(4)}(x,y)`$)

$$G_\varphi ^{(4)}(x,y)\frac{1}{|x|^{\eta (y)}}\mathrm{for}|x|1,$$
(64)
with

$$\eta (y)=\frac{1}{\pi ^2}_0^\pi 𝑑k\sqrt{\frac{𝒰(0,k)}{𝒦(0,k)}}\frac{\mathrm{sin}^2(ky/2)}{\mathrm{sin}(k/2)}.$$
(65)
Note that Eq. (64) requires only $`|x|1`$ and places no restriction on $`y`$. Hence it obtains also when both arguments are large, and thus

$$G_\varphi ^{(4)}(x,y)e^{\frac{1}{\pi ^2}\sqrt{\frac{𝒰_0}{𝒦_0}}(\mathrm{ln}x)(\mathrm{ln}y)C(\mathrm{ln}x+\mathrm{ln}y)},$$
(66)
as $`x,y\mathrm{}`$, with

$`C`$ $`=`$ $`{\displaystyle \frac{1}{\pi ^2}}\sqrt{{\displaystyle \frac{𝒰_0}{𝒦_0}}}(\gamma _E+\mathrm{ln}\pi )`$ (67)
$`+{\displaystyle \frac{1}{\pi ^2}}{\displaystyle _0^\pi }\left(\sqrt{{\displaystyle \frac{𝒰(0,k)}{𝒦(0,k)}}}{\displaystyle \frac{1}{2\mathrm{sin}k/2}}\sqrt{{\displaystyle \frac{𝒰_0}{𝒦_0}}}{\displaystyle \frac{1}{k}}\right).`$
### IV.3 Vortex correlation functions

The exact plaquette duality makes it possible to define a number of characteristic “vortex” correlators in terms of the field $`e^{i\theta }`$. As for the boson correlator above, the behavior in the Bose metal depends upon the expression for the vortex operator in the low energy variables. By symmetry, we expect the leading terms (involving the smallest exponentials of $`\vartheta `$) to be

$$e^{i\theta _𝐫}\mathrm{Re}\left[Ae^{i\pi \overline{n}xy}e^{i\vartheta _𝐫}\right],$$
(68)
where in general $`A`$ is complex. Note that, since $`\theta /\pi `$ is an integer, the microscopic vortex field satisfies $`e^{i\theta }=e^{i\theta }`$, but this is not true for $`\vartheta `$. In the special case of $`\overline{n}=1/2`$, particle-hole symmetry ($`\vartheta \vartheta `$) further implies $`A`$ is real, whence

$$e^{i\theta _𝐫}_{\overline{n}=1/2}2A\mathrm{cos}(\vartheta _𝐫\frac{\pi }{2}xy).$$
(69)
With this relation in hand, let us consider first the vortex two-point function,

$$G_\theta (𝐫,\tau )=e^{i\theta _𝐫(\tau )}e^{i\theta _\mathrm{𝟎}(0)}.$$
(70)
As for the boson correlator $`G_\varphi `$, the vortex two-point functions vanishes unless $`𝐫=\mathrm{𝟎}`$ due to the dual row and column symmetries. Evaluating it using Eq. (68) and the Gaussian action, one finds similar results to the boson tunneling density of states,

$$G_\theta (\mathrm{𝟎},\tau )\mathrm{exp}\left\{\frac{F_2^\theta }{2}|\mathrm{ln}\mathrm{\Omega }_0\tau |^2F_1^\theta |\mathrm{ln}\mathrm{\Omega }_0\tau |\right\},$$
(71)
where

$`F_2^\theta `$ $`=`$ $`{\displaystyle \frac{1}{2}}\sqrt{{\displaystyle \frac{𝒦_0}{𝒰_0}}},`$ (72)
$`F_1^\theta `$ $`=`$ $`{\displaystyle _0^\pi }𝑑k\left[\sqrt{{\displaystyle \frac{𝒦(k,0)}{𝒰(k,0)}}}{\displaystyle \frac{1}{2\mathrm{sin}k/2}}\sqrt{{\displaystyle \frac{𝒦_0}{𝒰_0}}}{\displaystyle \frac{1}{k}}\right]`$ (73)
$`{\displaystyle \frac{1}{2}}\sqrt{{\displaystyle \frac{𝒦_0}{𝒰_0}}}\left[2\mathrm{ln}\pi +\gamma _E\right].`$
In fact, these coefficients are obtained directly from Eq. (59) using the duality transformation $`𝒦(𝐤)𝒰(𝐤)/\pi ^2`$. In Fig. 2 we show numerical results for the two-point vortex correlation function obtained from the quantum Monte Carlo simulation in the parameter regime corresponding to the Bose metal phase. The downward curvature of the data is consistent with a decay more rapid than a power law, and as shown in the inset can be fit to the form in Eq. (71) with $`\frac{𝒰_0}{𝒦_0}1`$.

Next consider the vortex four-point function, which due to the dual row/column symmetry is non-vanishing only when the four points sit at the corners of a rectangle:

$`G_\theta ^{(4)}(x,y)`$ $`=`$ $`e^{i\theta _{oo}(\tau )}e^{i\theta _{xo}(\tau )}e^{i\theta _{oy}(\tau )}e^{i\theta _{xy}(\tau )}.`$ (74)
This can be re-expressed in terms of the low-energy field, $`\vartheta `$, using Eq. (68), and then evaluated with the Gaussian Bose metal action. For any fixed $`y`$, one finds power-law behavior in $`x`$,

$$G_\theta ^{(4)}(x,y)\frac{\mathrm{cos}(\pi \overline{n}xy)}{|x|^{\eta _v(y)}}\mathrm{for}|x|1,$$
(75)
with

$$\eta _v(y)=_0^\pi 𝑑k\sqrt{\frac{𝒦(0,k)}{𝒰(0,k)}}\frac{\mathrm{sin}^2(ky/2)}{\mathrm{sin}(k/2)}.$$
(76)
Interestingly, one can see from Eq. (75) that for the case $`\overline{n}=1/2`$, the four boson correlator vanishes exactly whenever $`xy`$ is odd as a consequence of particle-hole symmetry. This behavior, and the associated power-law correlations, are shown in Fig. 3. In the limit when both $`x,y\mathrm{}`$, the vortex four-point correlator vanishes faster than any power law:

$$G_\theta ^{(4)}(x,y)\mathrm{cos}(\pi \overline{n}xy)e^{\sqrt{\frac{𝒦_0}{𝒰_0}}(\mathrm{ln}x)(\mathrm{ln}y)C_v(\mathrm{ln}x+\mathrm{ln}y)},$$
(77)
as $`x,y\mathrm{}`$, with

$`C_v`$ $`=`$ $`\sqrt{{\displaystyle \frac{𝒦_0}{𝒰_0}}}(\gamma _E+\mathrm{ln}\pi )`$ (78)
$`+{\displaystyle _0^\pi }\left(\sqrt{{\displaystyle \frac{𝒦(0,k)}{𝒰(0,k)}}}{\displaystyle \frac{1}{2\mathrm{sin}k/2}}\sqrt{{\displaystyle \frac{𝒦_0}{𝒰_0}}}{\displaystyle \frac{1}{k}}\right).`$
With the exception of the $`\mathrm{cos}(\pi \overline{n}xy)`$ prefactor in Eq. (77), all the results in this subsection can be obtain from those in the previous one by the duality transformation $`\phi \vartheta `$ and $`𝒰(𝐤)\pi ^2𝒦(𝐤)`$.

### IV.4 Collective correlation functions

Next, let us consider the correlations of “two-particle” operators such as the boson and energy densities. As above, we require the operator correspondences between the microscopic and effective variables. These were already worked out in Sec. III, and summarized in Eqs. (9-10) in the introduction. For simplicity, we focus on the most interesting case of $`\overline{n}=1/2`$, for which $`Q=\pi `$, and furthermore keep only the lowest non-trivial harmonics with $`q=1`$. In this limit, Eqs. (9-10) reduce to

$`\delta n_{x+\frac{1}{2},y+\frac{1}{2}}`$ $``$ $`c_0^\rho \mathrm{\Delta }_{xy}\vartheta _{xy}`$ (79)
$`+c_2^\rho {\displaystyle \underset{a,b=0,1}{}}(1)^{(x+a)(y+b)+a+b}\mathrm{sin}(2\vartheta _{x+a,y+b}),`$
$`\epsilon _{xy}`$ $``$ $`c_0^\epsilon \dot{\vartheta }_{xy}^2+c_2^\epsilon (1)^{xy}\mathrm{cos}(2\vartheta _{xy}).`$ (80)
Consider first the density-density correlation function,

$$\chi _{nn}(𝐫,\tau )=\delta n_{x+\frac{1}{2},y+\frac{1}{2}}(\tau )\delta n_{\frac{1}{2},\frac{1}{2}}(0).$$
(81)
Substituting for $`\delta n`$ using Eq. (79), one obtains three contributions to $`\chi _{nn}`$:

$$\chi _{nn}(𝐫,\tau )(c_0^\rho )^2\chi _{nn}^{00}(𝐫,\tau )+(c_2^\rho )^2\chi _{nn}^{22}(𝐫,\tau )+c_0^\rho c_2^\rho \chi _{nn}^{20}(𝐫,\tau ).$$
(82)
The cross term $`\chi _{nn}^{20}`$ is negligible. The first contribution, $`\chi _{nn}^{00}(𝐫,\tau )`$ is just the correlator between “coarse-grained” densities $`\mathrm{\Delta }_{xy}\vartheta _{xy}`$. This term is non-zero in the Gaussian theory, and gives a smooth function of $`𝐫`$, with a power-law behavior at large arguments. For instance, at equal times and large $`|x|,|y|1`$,

$$\chi _{nn}^{00}(𝐫,0)\sqrt{\frac{𝒦_0}{𝒰_0}}\frac{1}{x^2y^2}.$$
(83)
More generally, $`\chi _{nn}^{00}`$ has a smooth Fourier transform,

$$\chi _{nn}^{00}(𝐪,\omega _n)=\frac{\pi ^2}{𝒰(𝐤)}\frac{E_𝐤^2}{\omega _n^2+E_𝐤^2}.$$
(84)
Perturbative corrections from $`𝒮_1`$ to $`\chi _{nn}^{00}`$ do not modify this qualitative behavior.

The remaining contribution to the density-density correlator comes from the $`\mathrm{sin}2\vartheta `$ terms in Eq. (79). Naively, this contribution is ultra-local (and hence uninteresting), i.e. vanishes unless $`|x|,|y|1`$, as a consequence of the fact that the discrete row/column symmetries are promoted to continuous ones at the Gaussian level. One may interpret this as meaning that the “vorticity” on each row or column is conserved exactly in the Gaussian model. This conclusion, however, is incorrect once the non-quadratic corrections in $`𝒮_1`$ are taken into account, since expanding factors of $`\upsilon _2(1)^{xy}\mathrm{cos}2\vartheta `$ can “supply” vorticity in units of $`2`$ to a particular site. Hence, $`\chi _{nn}^{22}`$ becomes non-trivial at second order in $`\upsilon _2`$. Provided $`|x|,|y|>1`$, the appropriate second order perturbative expression is

$`\chi _{nn}^{22}(𝐫,\tau )`$ $`={\displaystyle \frac{\upsilon _2^2}{8}}{\displaystyle \underset{abcd=0,1}{}}`$ $`(1)^{(x+ca)(y+db)+a+b+c+d}`$ (85)
$`\times C_{x+ca,y+db}^{(4)}(\tau ),`$
where

$$C_{x,y}^{(4)}(\tau )=𝑑\tau _1𝑑\tau _2e^{2i[\vartheta _{oo}(0)+\vartheta _{xy}(\tau )\vartheta _{x0}(\tau _1)\vartheta _{0y}(\tau _2)]}_0,$$
(86)
where the expectation value indicates an average calculated in the Gaussian theory. At large distances, $`|x|,|y|1`$, one expects $`C_{x+ca,y+db}^{(4)}(\tau )C_{x,y}^{(4)}(\tau )`$, independent of $`a,b,c,d`$. In this limit, only the prefactor in Eq. (85) depends upon $`a,b,c,d`$, and the sums can be carried out explicitly. To do so, it is convenient to employ the identity

$$(1)^{xy}=\frac{1}{2}[1+(1)^x+(1)^y(1)^{x+y}],$$
(87)
valid for integer $`x,y`$. Applying this identity to Eq. (85), only the last term survives the sum, and gives

$$\chi _{nn}^{22}(𝐫,\tau )\upsilon _2^2(1)^{x+y}C_{x,y}^{(4)}(\tau ),$$
(88)
for $`|x|,|y|1`$. This indicates the presence of $`(\pi ,\pi )`$ correlations in the boson density. More generally, if only $`|x|1`$ but $`|y|>1`$ but still of $`O(1)`$, one finds

$$\chi _{nn}^{22}(𝐫,\tau )\frac{\upsilon _2^2}{4}(1)^{x+y}\underset{bd=0,1}{}(1(1)^{y+db})C_{x,y+db}^{(4)}(\tau ),$$
(89)
so at any fixed $`y`$ the density-density correlator oscillates at wavevector $`\pi `$ as a function of $`x`$ (and vice-versa by rotational symmetry). To establish the range of these correlations, we must consider $`C_{x,y}^{(4)}(\tau )`$ in some detail. Using properties of Gaussian fields, we have

$$C_{x,y}^{(4)}(\tau )=𝑑\tau _1𝑑\tau _2\mathrm{exp}\left[c_{x,y}(\tau ,\tau _1,\tau _2)\right],$$
(90)
where

$`c_{x,y}(\tau ,\tau _1,\tau _2)`$ $`=`$
$`2\left(\vartheta _{xy}(\tau )+\vartheta _{00}(0)\vartheta _{x0}(\tau _1)\vartheta _{0y}(\tau _2)\right)^2_0.`$
The calculation and analysis of $`c_{xy}`$ is somewhat involved, and is described in detail in Appendix B. There we derive a useful approximation for $`c_{xy}`$ which captures the behavior in all the relevant limits ($`|x|1,\mathrm{\Omega }_0|\tau \tau _1|,\mathrm{\Omega }_0|\tau _2|`$, and $`|y|<|x|`$):

$`c_{xy}`$ $``$ $`2\eta _v(y)\mathrm{ln}(x^2+\mathrm{\Omega }_y^2\tau ^2)`$ (92)
$`+\sqrt{{\displaystyle \frac{𝒦_0}{𝒰_0}}}\left([{\displaystyle \frac{\mathrm{\Omega }_0|\tau \tau _1|}{|y|}}]+[{\displaystyle \frac{\mathrm{\Omega }_0|\tau _2|}{|y|}}]\right),`$
where $`\mathrm{\Omega }_y\mathrm{\Omega }_0`$ (but in general depends weakly upon the full form of $`𝒦(0,k_y),𝒰(0,k_y),`$ and $`y`$) and the crossover function $`(𝒳)`$ satisfies

$$(𝒳)\{\begin{array}{cc}0\hfill & 𝒳1\hfill \\ \mathrm{ln}^2𝒳\hfill & 𝒳1\hfill \end{array}.$$
(93)
Using Eqs. (92,93,90), we find

$$C_{xy}^{(4)}(\tau )\stackrel{}{\left|x\right|1}\left(\frac{y}{\mathrm{\Omega }_0}\right)^2e^{2\left(\frac{𝒦_0}{𝒰_0}\right)^{1/4}}\left(x^2+\mathrm{\Omega }_y^2\tau ^2\right)^{2\eta _v(y)},$$
(94)
and of course the same behavior with $`xy`$.

This power-law behavior of $`C_{xy}^{(4)}`$, and hence $`\chi _{nn}(𝐫,\tau )`$ translates into singularities in the static structure factor,

$`S_{nn}(𝐤)`$ $``$ $`{\displaystyle \underset{xy}{}}e^{i(k_xx+k_yy)}\chi _{nn}(𝐫,\tau =0),`$ (95)
for wavevectors near $`𝝅(\pi ,\pi )`$. In particular, let $`𝐤=𝝅+𝐪`$, and consider the limit $`q_x1`$ with $`q_y`$ fixed. In this limit the singular behavior of the structure factor is dominated by large $`|x|`$ but $`y`$ of $`O(1)`$. Hence we may apply Eq. (89) to write

$`S_{nn}(𝝅+𝐪)`$ $`\stackrel{}{\left|q_x\right|1}`$
$`\upsilon _2^2{\displaystyle \underset{\mathrm{odd}y}{}}\mathrm{cos}q_yy(1+\mathrm{cos}q_y){\displaystyle \underset{x}{}}e^{iq_xx}C_{xy}^{(4)}(0).`$
From the power-law behavior in Eq. (94), one can readily see that the Fourier transform in $`x`$ leads to singularities for small $`q_x`$. Indeed, if any $`\eta _v(y)<1/4`$ for any odd $`y`$, $`\chi _{nn}`$ diverges as $`q_x0`$ at fixed $`q_y`$, while for $`1/4<\eta _v<3/4`$, the structure factor remains finite but has a divergent second derivative at $`q_x=0`$:

$`S_{nn}(𝝅+𝐪)`$ $``$ $`S_{nn}^0(𝐪)+{\displaystyle \underset{\mathrm{odd}y}{}}A_y^\rho (q_y)\mathrm{sgn}(\alpha _y)|q_x|^{\alpha _y},`$ (97)
where $`S_{nn}^0(𝐪)`$ is a smooth function, $`\alpha _y=14\eta _v(y)`$, and $`A_y^\rho (q_y)\upsilon _2^2\mathrm{cos}q_yy(1+\mathrm{cos}q_y)`$ is a positive amplitude peaked at $`q_y=0`$. Hence the behavior for small $`q_x`$ is dominated by the minimum (over odd $`y`$) value of $`\eta _v(y)`$ (maximum of $`\alpha _y`$).

The zero frequency susceptibility, $`\chi _{nn}(𝐤)\chi _{nn}(𝐤,\omega _n=0)`$ has a similar but stronger divergence due to the extra time integration:

$`\chi _{nn}(𝝅+𝐪)`$ $``$ $`\chi _{nn}^0(𝐪)+{\displaystyle \underset{\mathrm{odd}y}{}}\stackrel{~}{A}_y^\rho (q_y)\mathrm{sgn}(\stackrel{~}{\alpha }_y)|q_x|^{\stackrel{~}{\alpha }_y},`$ (98)
where $`\chi _{nn}^0(𝐪)`$ is another smooth function, $`\stackrel{~}{A}_y^\rho (q_y)A_y^\rho (q_y)`$, and $`\stackrel{~}{\alpha }_y=24\eta _v(y)=1+\alpha _y`$ signals the stronger divergence. The difference in exponents implies that if $`S_{nn}(𝝅+𝐪)`$ has a divergent slope at $`q_x=0`$, $`\chi _{nn}(𝝅+𝐪)`$ itself divergences there. Conversely, if $`\stackrel{~}{\chi }_{nn}(𝝅+𝐪)`$ has a slope divergence at $`q_x=0`$ (occurs for $`1/2<\eta _v<3/4`$), the static structure factor does not. Numerical results for the density susceptibility in the Bose metal phase are shown in Fig. 4, and reveal a singular cusp-like behavior at wavevector $`𝝅`$. This form is consistent with that predicted by Eq. (98) with a maximum value $`1<\stackrel{~}{\alpha }_y^{max}<0`$.

Now consider the energy-energy correlation function,

$$\chi _{\epsilon \epsilon }(𝐫,\tau )=\epsilon _{xy}(\tau )\epsilon _{xy}(0)_c,$$
(99)
where the $`c`$ subscript indicates the cumulant expectation value. As for the density-density correlator, we can employ Eq. (80) to expand $`\chi _{\epsilon \epsilon }`$ as

$$\chi _{\epsilon \epsilon }(𝐫,\tau )(c_0^\epsilon )^2\chi _{\epsilon \epsilon }^{00}(𝐫,\tau )+(c_2^\epsilon )^2\chi _{\epsilon \epsilon }^{22}(𝐫,\tau )+c_0^\epsilon c_2^\epsilon \chi _{\epsilon \epsilon }^{20}(𝐫,\tau ).$$
(100)
As for the density case, $`\chi _{\epsilon \epsilon }^{00}`$ has a smooth Fourier transform, and $`\chi _{\epsilon \epsilon }^{20}`$ is negligible. We focus on $`\chi _{\epsilon \epsilon }^{22}`$, which (for $`|x|,|y|>0`$) to second order in $`\upsilon _2`$ is

$$\chi _{\epsilon \epsilon }^{22}=\frac{\upsilon _2^2}{8}(1)^{xy}C_{xy}^{(4)}(\tau ).$$
(101)
Using Eq. (87,94), one straightforwardly sees that there are singularities in $`\chi _{\epsilon \epsilon }^{22}(𝐤)`$ as $`k_x(k_y)0`$ and $`k_x(k_y)\pi `$. In particular,

$`S_{\epsilon \epsilon }(𝐤)`$ $`\stackrel{}{\left|k_x\right|1}`$ $`S_{\epsilon \epsilon }^0(𝐤)+{\displaystyle \underset{\mathrm{even}y}{}}A_y^\epsilon (k_y)\mathrm{sgn}(\alpha _y)|k_x|^{\alpha _y},`$ (102)
$`\chi _{\epsilon \epsilon }(𝐤)`$ $`\stackrel{}{\left|k_x\right|1}`$ $`\chi _{\epsilon \epsilon }^0(𝐤)+{\displaystyle \underset{\mathrm{even}y}{}}\stackrel{~}{A}_y^\epsilon (k_y)\mathrm{sgn}(\stackrel{~}{\alpha }_y)|k_x|^{\stackrel{~}{\alpha }_y},`$ (103)
and

$`S_{\epsilon \epsilon }(\pi +q_x,k_y)`$ $`\stackrel{}{\left|q_x\right|1}`$
$`S_{\epsilon \epsilon }^0(\pi +q_x,k_y)+{\displaystyle \underset{\mathrm{odd}y}{}}A_y^\epsilon (k_y)\mathrm{sgn}(\alpha _y)|q_x|^{\alpha _y},`$
$`\chi _{\epsilon \epsilon }(\pi +q_x,k_y)`$ $`\stackrel{}{\left|q_x\right|1}`$
$`\chi _{\epsilon \epsilon }^0(\pi +q_x,k_y)+{\displaystyle \underset{\mathrm{odd}y}{}}\stackrel{~}{A}_y^\epsilon (k_y)\mathrm{sgn}(\stackrel{~}{\alpha }_y)|q_x|^{\stackrel{~}{\alpha }_y},`$
where $`A_y^\epsilon (k_y)\stackrel{~}{A}_y^\epsilon (k_y)\mathrm{cos}k_yy`$. Similar formulae for $`k_y0,\pi `$ are obtained by rotation. Note that because $`A_y^\epsilon (k_y)\mathrm{cos}k_yy`$ is negative for $`k_y\pi `$ and $`y`$ odd, the energy-energy correlator has a singular dip near $`𝐤=𝝅`$ and singular peaks near $`𝐤=(0,0),(\pi ,0),`$ and $`(0,\pi )`$. Numerical results for the energy susceptibility in the Bose metal phase are shown in Fig. 5. Notice the singular cusp-like peaks at wavevectors $`(0,\pi )`$ and $`(\pi ,0)`$, as predicted in Eq. (IV.4).

### IV.5 Electrical Conductivity

We finally consider the electrical conductivity in the Bose metal. To obtain an expression for the current operator requires coupling in a vector potential. The usual minimal coupling prescription for a lattice model of bosons would have one replace a discrete lattice derivative as,

$$\varphi _{𝐫+𝝁}\varphi _𝐫\varphi _{𝐫+𝝁}\varphi _𝐫+A_𝐫^\mu ,$$
(106)
with $`𝝁=\widehat{𝐱},\widehat{𝐲}`$, and define the current operator by differentiating with respect to the vector potential. For the boson ring model this prescription is ambiguous due to the second derivative form, and a family of gauge inequivalent forms are possible:

$$\mathrm{\Delta }_{xy}\varphi _𝐫\mathrm{\Delta }_{xy}\varphi _𝐫+\alpha \mathrm{\Delta }_xA_𝐫^y+\beta \mathrm{\Delta }_yA_𝐫^x,$$
(107)
where $`\mathrm{\Delta }_xf_𝐫f_{𝐫+\widehat{𝐱}}f_𝐫`$ denotes a discrete derivative. Gauge invariance requires $`\alpha +\beta =1`$, whence this can be re-expressed as,

$$\mathrm{\Delta }_{xy}\varphi _𝐫\mathrm{\Delta }_{xy}\varphi _𝐫+\mathrm{\Delta }_xA_𝐫^y+(\alpha 1)\mathrm{\Phi },$$
(108)
where $`\mathrm{\Phi }=\mathrm{\Delta }_xA^y\mathrm{\Delta }_yA^x`$ is the gauge invariant flux through the plaquette. If one derives the boson ring model by starting with a model of electrons reformulated in terms of a $`Z_2`$ gauge field coupled to spinons and chargons as detailed in Appendix D, one arrives at an appealing symmetric form for the ring term with $`\alpha =\beta =1/2`$.

Henceforth we focus on the zero wavevector conductivity, $`\sigma (\omega )`$. In this case, the above ambiguity is irrelevant, since a spatially uniform vector potential does not enter for any value of $`\alpha `$. But the vector potential will still of course enter into the boson hopping term in Eq. (4), and upon differentiation generates the usual current operator,

$$I_𝐫^\mu =t\mathrm{sin}(\varphi _{𝐫+𝝁}\varphi _𝐫).$$
(109)
When we coarse grain the theory, we should write down the current operator in terms of the slow field, $`\phi `$:

$$I_𝐫^\mu =c_It\mathrm{sin}(\phi _{𝐫+𝝁}\phi _𝐫)+\mathrm{},$$
(110)
with $`c_I`$ a dimensionless constant. The other contributions will include terms which hop a single boson several lattice spacings and terms which hop several bosons - generally all local terms allowed by the symmetries. Such terms will generically be subdominant at low frequencies, and so we retain only the leading contribution.

With the current operator in hand it is straightforward to obtain the conductivity from the usual Kubo expression:

$$\sigma _1(\omega )=\mathrm{Re}\sigma (i\omega _n\omega ),$$
(111)
with

$$\sigma (i\omega _n)=\frac{1}{\omega _n}𝑑\tau \underset{𝐫}{}e^{i\omega _n\tau }I_𝐫^y(\tau )I_\mathrm{𝟎}^y(0).$$
(112)
We have dropped the diamagnetic contribution since it will not contribute to the real part of the conductivity. Evaluating the correlator using the Gaussian Bose metal action gives at zero temperature,

$$I_𝐫^y(\tau )I_\mathrm{𝟎}^y(0)_0t^2\delta _{y,0}(x^2+\tau ^2)^\mathrm{\Delta },$$
(113)
with scaling dimension $`\mathrm{\Delta }=\eta (1)/2`$ where $`\eta (y)`$ is given explicitly in terms of the Bose liquid parameters in Eq. 65. Performing the time integration and spatial summation gives a power law singular contribution,

$$\omega _n\sigma (i\omega _n)=At^2(1)^{\mathrm{int}(\mathrm{\Delta })}\omega _n^{2\mathrm{\Delta }2}+\mathrm{\Sigma }_{reg}(i\omega _n),$$
(114)
where $`A`$ is a positive constant and $`\mathrm{\Sigma }_{reg}`$ is analytic in it’s argument and thus does not contribute to the real part of the conductivity. Analytic continuation to real frequencies gives a singular contribution to the complex conductivity of the form,

$$\sigma (\omega )=At^2(1)^{\mathrm{int}(\mathrm{\Delta })}(i\omega )^{2\mathrm{\Delta }3}+i\sigma _2(\omega ),$$
(115)
with real $`\sigma _2(\omega )`$. We note that causality places strong constraints on the phase angledorsey  from the singular contribution, consistent with the above form. Finally, taking the real part gives the optical conductivity,

$$\sigma _1(\omega )=At^2|\mathrm{sin}(\pi \mathrm{\Delta })||\omega |^{2\mathrm{\Delta }3}.$$
(116)
Notice that the amplitude of this contribution vanishes for integer $`\mathrm{\Delta }`$. For these special cases, the higher order contributions to the current operator should be kept, and will contribute a similar form but generally with larger scaling dimensions.

As we discuss in the next section, stability of the Bose metal phase to boson hopping requires that $`\mathrm{\Delta }>2`$, implying an optical conductivity vanishing rapidly at low frequencies, $`\sigma _1(\omega )\omega ^\alpha `$ with $`\alpha >1`$. For $`\mathrm{\Delta }<2`$, the Bose metal will be unstable to superconductivity, but for small boson hopping amplitude the transition temperature would be low. In that case, the optical conductivity above $`T_c`$ might still be well described by the above power law form.

One may also directly calculate the non-linear dc conductivity, $`\sigma (E_y,T)=I_y/E_y`$ as a function of $`E_y`$ and $`T`$. This is mathematically complicated, but formally quite similar to the perturbative calculations of tunneling conductance between parallel Luttinger liquids, as carried out e.g. in Ref. Carpentier, . We forgo this calculation here for the sake of brevity, quoting only the resulting scaling form

$$\sigma (E_y,T)=At^2T^{2\mathrm{\Delta }3}𝒢(E_y/T),$$
(117)
where $`𝒢()1`$ as $`0`$, and $`()^{2\mathrm{\Delta }3}`$ for $`1`$. This implies in particular that the DC (linear response) conductivity at non-zero temperature behaves as $`\sigma (T)At^2T^{2\mathrm{\Delta }3}`$. These results obtain for the pure ring model, but will doubtless be altered somewhat in the presence of impurity scattering. This we leave for future study.

## V Instabilities of the Bose metal

There are two classes of perturbations that one can add to the Bose metal fixed point (described by the quadratic action $`𝒮_0`$) that can potentially destabilize the phase. The first are terms involving the hopping of bosons. When relevant, such boson hopping terms “stiffen” the phase fluctuations of the bosons, leading to off-diagonal long range order and a superconducting state. As we shall see, the perturbative relevance of the boson hopping terms is determined by the Bose liquid parameters which enter into the fixed point action which characterizes the Bose metal phase. There are regions of parameters where all such boson hopping terms are irrelevant and the Bose metal phase is stable to such superconducting perturbations, as we detail in Subsection A below.

The second class of perturbations involve hopping or motion of vortices, conveniently expressed in the dual representation. When relevant, these perturbations signal a condensation of vortices which typically leads to a breaking of translational symmetry and drives the system into an incompressible insulating state. The presence of these instabilities and the precise form of the translational symmetry breaking depend sensitively on the boson density, generally requiring boson densities commensurate with the underlying lattice. In subsection B below we focus on half-filling ($`\overline{n}=1/2`$), and study the nature of the resulting commensurate insulating states.

If the boson density is incommensurate with the lattice, on the other hand, small vortex hopping terms are unimportant. Provided the Bose liquid parameters are in the regime where Boson hopping is likewise irrelevant, the Bose metal exists as a completely stable critical phase. The gapless Bose metal is then a 2d analog of the stable 1d Luttinger liquid.

### V.1 Boson hopping and superconductivity

Consider then the stability of the Bose metal phase in the presence of boson hopping operators. To be general, we consider processes where $`q`$ bosons hop along a displacement vector $`𝐬`$:

$$_t=\underset{𝐫}{}\underset{q,𝐬}{}t_q^𝐬\mathrm{cos}[q(\phi _𝐫\phi _{𝐫+𝐬})].$$
(118)
To assess the perturbative relevance of such processes, we compute the two-point function of the tunneling operators,

$$𝒯_q^𝐬(𝐫)=\mathrm{cos}[q(\phi _𝐫\phi _{𝐫+𝐬})],$$
(119)
using the Gaussian action for the Bose metal. The row/column symmetries in the Bose metal phase greatly constrain the spatial dependence of these correlators. We consider first the special class of tunneling operators which hop bosons along the $`x`$ or $`y`$ axes with $`𝐬=(m,0)`$ or $`𝐬=(0,m)`$. For this class of operators one finds a power law decay both in time and in one spatial dimension:

$$𝒯_q^{m\widehat{𝐲}}(𝐫,\tau )𝒯_q^{m\widehat{𝐲}}(\mathrm{𝟎},0)_0=\delta _{y,0}[x^2+v^2\tau ^2]^{\mathrm{\Delta }_{qm}},$$
(120)
with scaling dimension,

$$\mathrm{\Delta }_{qm}=\frac{q^2}{2\pi }_\pi ^\pi \frac{dk_y}{2\pi }\sqrt{\frac{𝒰(0,k_y)}{𝒦(0,k_y)}}\frac{\mathrm{sin}^2(mk_y/2)}{|\mathrm{sin}(k_y/2)|}=\frac{q^2}{2}\eta (m).$$
(121)
Notice that these scaling dimensions vary as $`q^2`$, and will generally increase (slowly – $`\mathrm{\Delta }_{qm}[q^2/(2\pi ^2)]\sqrt{𝒰_0/𝒦_0}\mathrm{ln}|m|`$) with hopping distance $`m`$. As expected these scaling dimensions increase with $`U/K`$.

For the remaining boson tunneling operators with vector $`𝐬`$ connecting two different rows and two different columns, the correlator is further constrained by the row/column symmetry, being spatially local,

$$𝒯_q^𝐬(𝐫,\tau )𝒯_q^𝐬(\mathrm{𝟎},0)_0=G(\tau )\delta _{𝐫,\mathrm{𝟎}}\mathrm{for}s_xs_y0,$$
(122)
with $`G(\tau )|\tau |^{2\mathrm{\Delta }_q^𝐬}`$.

General renormalization group reasoning implies that operators are irrelevant about a given fixed point when the associated scaling dimension exceeds the space-time dimension, $`\mathrm{\Delta }>D=d+1`$, with dynamical exponent $`z=1`$. But due to the constrained form of the above correlators which only exhibit a power law decay in a reduced set of space-time dimensions, $`D_{red}`$, one expects that the condition for irrelevance should be modified to be $`\mathrm{\Delta }>D_{red}`$. Thus, when $`\mathrm{\Delta }_{q,m}>2`$ and $`\mathrm{\Delta }_q^𝐬>1`$, one expects that all of the boson hopping operators should be unimportant as one scales down in energy, and the Bose metal phase will be stable. This can always be achieved (in principle) by increasing the ratio of $`U/K`$.

It is tempting to strengthen this expectation by constructing an explicit renormalization group transformation, but this is somewhat problematic due to the peculiar existence of zero energy states at both $`k_x=0`$ and $`k_y=0`$ in the Bose metal phase. One could try to integrate out gapped modes away from the zero energy “cross” in the Brillouin zone, and then successively integrate out “shells” of modes pinching down onto the cross. A difficulty arises, however, in the rescaling transformation of the momenta, because the range of both $`k_x`$ and $`k_y`$ – the interval $`[\pi ,\pi ]`$ – is not invariant. Similar but less severe difficulties are encountered when one tries to implement a momentum shell RG procedure for a 2d Fermi surface, due to the possible modifications of the shape/size of the Fermi surface. It might be possible to circumvent this difficulty along the lines of Shankar et. al.Shankar , or by an RG procedure in frequency space.unpub  Some insight can be gleaned by ignoring the zero modes along the $`k_x=0`$ axis, and considering a 1d RG transformation where the integration is over a shell of momenta in $`k_x`$ for all $`k_y[\pi ,\pi ]`$ and frequency, and then rescaling both $`k_x`$ and $`\omega `$ (ie. with $`z=1`$) but not $`k_y`$. The resulting perturbative (linearized) RG equation for the Boson hopping amplitude in the $`y`$direction, $`t_{q,m}t_q^{m\widehat{𝐲}}`$ is then simply,

$$_{\mathrm{}}t_{q,m}=(1+z\mathrm{\Delta }_{q,m})t_{q,m},$$
(123)
with $`z=1`$ and $`\mathrm{\Delta }_{q,m}`$ the scaling dimension given explicitly above. This argument indeed supports the expectation that weak Boson hopping will be irrelevant provided $`\mathrm{\Delta }_{q,m}>D_{red}=2`$. In the absence of a fully controlled 2d RG procedure, we can verify that this conclusion in fact correct by resorting to perturbation theory. In Appendix C, we do just this by formally computing corrections to the two-point function, $`e^{i\phi _𝐫(\tau )}e^{i\phi _𝐫(0)}`$ perturbatively in powers of the Boson hopping $`t_{q=1,m=1}`$. Stability of the Bose metal phase requires that the long-time behavior be unmodified from that calculated within the Gaussian theory. Carrying out this expansion to second order in $`t_{1,1}`$, we show explicitly that this is in fact the case provided $`\mathrm{\Delta }_{1,1}>2`$, thereby confirming the expected result.

When the scaling dimension of a single boson hopping term is sufficiently small, $`\mathrm{\Delta }_{1m}<2`$, one expects the Bose metal phase to be unstable to a superfluid state in which the bosons condense, with $`e^{i\phi }0`$.bosecondensenote  If the original lattice Bosons are supposed to represent the Cooper pairs (say in a model of the cuprates), then this will of course be the superconducting phase.

### V.2 Vortex hopping and insulating states

#### V.2.1 Stability

We next consider the effects of the various non-linear interactions involving the vortex operators in Eq. (38) and (39), which can potentially destabilize the Bose metal phase. At a generic density, for which $`\overline{n}`$ is irrational, all these cosine terms are oscillatory, and, if weak, cannot lead to any long-wavelength divergences. Hence we expect that away from commensurate densities, and provided $`\mathrm{\Delta }_{qm}>2`$ for all $`q,m`$ so that boson tunneling is irrelevant, the Bose metal is a stable zero temperature phase of matter.

For rational densities, some of the vortex operators will be non-oscillatory, and must be considered more carefully. In particular, for very commensurate boson densities (i.e. $`\overline{n}`$ equal to a small-denominator rational fraction), it is not apriori obvious whether the Bose metal can even in principle be stable to both boson hopping and vortex operators. To demonstrate the issues, we present a stability analysis for the special set of rational fillings, $`\overline{n}=1/z`$, where $`z`$ is a positive integer.

At $`\overline{n}=1/z`$, the on-site terms take the form,

$$_v=\underset{𝐫}{}\underset{q=1}{\overset{\mathrm{}}{}}\upsilon _{2q}\mathrm{cos}(2q\vartheta +\frac{2\pi q}{z}xy).$$
(124)
Note that although for $`qjz`$, with integer $`j1`$, these terms are spatially-varying, they are not wholly oscillatory, in the sense that they all have a non-zero spatial average for constant $`\vartheta `$, e.g. for prime $`z`$, and $`qjz`$, $`[\mathrm{cos}(\frac{2\pi q}{z}xy)]_{xy}=1/z`$. Hence even for $`qjz`$, they cannot be argued away simply on the grounds that they are oscillatory. Instead, to assess the importance of these perturbations, we again consider the two-point function of these operators evaluated with the Gaussian Bose metal action,

$$G_\vartheta ^{2q}(𝐫,\tau )=e^{i2q\vartheta _𝐫(\tau )}e^{i2q\vartheta _\mathrm{𝟎}(0)}_0=\delta _{𝐫,\mathrm{𝟎}}F(\tau ),$$
(125)
with $`F(\tau )e^{q^2\sqrt{𝒦_0/𝒰_0}\mathrm{ln}^2(\tau )}`$ at large times. Due to the dual row/column symmetry this correlator is spatially local, and is also “short-ranged” in time, vanishing faster than any power law. The associated scaling dimension is thus infinite, and these operators are strongly irrelevant and will not destabilize the Bose metal. However, as we shall see, they will play an important role in determining which of the various insulating phases is selected when the Bose metal is driven unstable by the two-site terms.

When the two-site terms in Eq. (39) are weak, we can drop those that are spatially oscillating and focus on the rest which take the form of vortex hopping and creation terms:

$$_w=\underset{𝐫,a,\pm }{}w_{a,\pm }𝒪_{a,\pm }(𝐫),$$
(126)
with operators, $`𝒪_{a,\pm }(𝐫)=\mathrm{cos}(2q\vartheta _𝐫\pm 2q^{}\vartheta _{𝐫+𝐬})`$, where “a” labels the various values of integers $`q,q^{}`$ and the hopping vector $`𝐬`$. Again, to establish the perturbative relevance of such terms in the Bose metal phase, we evaluate the two-point correlators with the Gaussian fixed point action. For the operators $`𝒪_+`$, the dual row/column symmetry is especially restrictive and we find that the associated correlator is once again spatially local:

$$𝒪_{a,+}(𝐫,\tau )𝒪_{a,+}(\mathrm{𝟎},0)_0=F(\tau )\delta _{𝐫,\mathrm{𝟎}},$$
(127)
and “short-ranged” in time with $`F(\tau )e^{c_{a,q}ln^2(\tau )}`$ (with constant $`c_{a,q}`$), so that like the on-site terms, these operators cannot destabilize the Bose metal. Similar behavior is found for $`𝒪_{}`$, except for the special class of operators with $`q=q^{}`$, which correspond physically to $`2q`$ vortices hopping along a vector $`𝐬`$. Of these, except for the special case with $`𝐬=(m,0)`$ and $`𝐬=(0,n)`$, the two-point function is again spatially local, although it is now a power law in time, $`F(\tau )|\tau |^{2\mathrm{\Delta }}`$. While such operators could potentially destabilize the Bose metal if the power $`\mathrm{\Delta }<1`$, they will generally be less singular than the remaining class of operators with vortices hopping along the $`x`$ or $`y`$ axes:

$`_w=`$ $``$ $`w_{2q,m}{\displaystyle \underset{xy}{}}{\displaystyle \underset{qm=jz}{}}\{\mathrm{cos}[2q(\vartheta _{xy}\vartheta _{x+m,y})]`$ (128)
$`+`$ $`\mathrm{cos}[2q(\vartheta _{xy}\vartheta _{x,y+m})]\}.`$
The operators with $`qmjz`$ (with integer $`j`$) are truly oscillatory (i.e. have zero spatial average for constant $`\vartheta `$) at $`\overline{n}=1/z`$ and have been dropped. We expect the coefficients $`w_{2q,m}`$ to be positive (with the minus sign out front), since these operators will be generated at second order from the irrelevant on-site terms: $`w_{2q,m}\upsilon _{2q}^2`$.

For this last class of operators the two-point function is delta-correlated in one-spatial direction but a power law in the other and in time:

$$𝒪_{2q,m}^y(𝐫,\tau )𝒪_{2q,m}^y(\mathrm{𝟎},0)_0\delta _{y,0}(x^2+v^2\tau ^2)^{\mathrm{\Delta }_{2q,m}^v},$$
(129)
where $`𝒪_{2q,m}^y(x,y)=\mathrm{cos}[2q(\vartheta _{xy}\vartheta _{x,y+m})]`$. The associated scaling dimension is given in terms of the Bose liquid parameters:

$$\mathrm{\Delta }_{2q,m}^v=2q^2_0^\pi 𝑑k_y\sqrt{\frac{𝒦(0,k_y)}{𝒰(0,k_y)}}\frac{\mathrm{sin}^2(mk_y/2)}{\mathrm{sin}(k_y/2)}=2q^2\eta _v(m).$$
(130)
As for the boson hopping operators, the vortex hopping scaling dimensions vary as $`q^2`$, and also increase slowly with increasing hopping distance, $`m`$. In general, we expect that for large $`z`$ the most relevant of these will be $`w_{2,z}`$, with scaling dimension $`\mathrm{\Delta }_{2,z}^v=2\eta _v(z)2\sqrt{𝒦_0/𝒰_0}\mathrm{ln}|z|`$, so for sufficiently large $`z`$ this can be rendered arbitrarily large, and hence all the vortex hopping terms can be made irrelevant for any $`𝒰(𝐤),𝒦(𝐤)`$. Thus for sufficiently large $`z`$, such that $`\mathrm{\Delta }_{2q,m}>D_{red}=2`$ for all $`qm=jz`$, there is certainly a domain of stability for the Bose metal phase (i.e. where the Bose liquid parameters are further tuned to make all hopping terms irrelevant).

For small $`z`$, however, this is not clear. Indeed, it is straightforward to show that for $`z=2`$ (i.e. half-filling, $`\overline{n}=1/2`$), a choice of uniform functions $`𝒰(𝐤)=𝒰_0,𝒦(𝐤)=𝒦_0`$ does not lead to a stable regime. In particular, in this case, the leading vortex hopping operator has scaling dimension $`\mathrm{\Delta }_{2,2}^v=2\eta _v(2)=4\sqrt{𝒦_0/𝒰_0}`$, while the leading boson hopping operator has scaling dimension $`\mathrm{\Delta }_{11}=\frac{1}{\pi ^2}\sqrt{𝒰_0/𝒦_0}`$. Hence $`\mathrm{\Delta }_{22}^v\mathrm{\Delta }_{11}=4/\pi ^2`$, which implies $`\mathrm{Min}(\mathrm{\Delta }_{11},\mathrm{\Delta }_{22}^v)2/\pi <2`$, so that at least one or the other operator is relevant. We have not determined whether a stable Bose metal phase might be possible at half filling when $`𝒰(𝐤),𝒦(𝐤)`$ are momentum-dependent.

#### V.2.2 Instabilities

As indicated above, for small $`z`$, either a vortex or boson hopping instability may be inevitable. In any case, it is interesting to study the nature of the state resulting from relevant vortex hopping terms. Here, we briefly study the nature of the resulting phase for half-filling, taking into account the presence of the two most potentially relevant operators, with $`q=1,m=2`$ and $`q=2,m=1`$. In general, at half-filling, there are Bose liquid parameters for which both $`w_{22}`$ and $`t_{1,1}`$ are relevant, and more complex behavior may well occur in this regime. We will however neglect boson hopping completely, as appropriate for large $`𝒰/𝒦`$.

To this end, let us assume that $`w_{2,2}`$ is the most relevant operator, with $`\mathrm{\Delta }_{2,2}^v<2`$. We also assume that $`w_{4,1}`$ is the next most relevant (it could be irrelevant, but still the least irrelevant remaining operator). Then provided all bare couplings are small, we imagine integrating down in energy until $`w_{2,2}`$ becomes comparable to the Gaussian terms in $`𝒮_0`$. This requires $`w_{2,2}\mathrm{\Lambda }^{\mathrm{\Delta }_{2,2}^v}𝒰(\pi ,0)\mathrm{\Lambda }^2`$, which always occurs for sufficiently small $`\mathrm{\Lambda }`$ with $`\mathrm{\Delta }_{2,2}^v<2`$. At this point it seems appropriate to minimize the potential $`w_{2,2}[\mathrm{cos}(2\vartheta _{xy}2\vartheta _{x+2,y})+\mathrm{cos}(2\vartheta _{xy}2\vartheta _{x,y+2})]`$, simultaneously with $`𝒮_0`$. This minimization somewhat underestimates fluctuation effects, which will be commented upon later. The most general form for $`\vartheta _{xy}`$ which minimizes the $`w_{2,2}`$ term and keeps the Gaussian action small can be written

$$\vartheta _𝐫=\mathrm{\Theta }_0(𝐫)+(1)^x\mathrm{\Theta }_1(𝐫)+(1)^y\mathrm{\Theta }_2(𝐫)+\frac{\pi ax}{2}+\frac{\pi by}{2},$$
(131)
where $`\mathrm{\Theta }_{0/1/2}(𝐫)`$ are slowly-varying functions of space. Inserting this expression for $`\vartheta `$ into the action including $`𝒮_0`$ and the $`w_{2,2}`$ term only gives

$`𝒮_0^{\mathrm{eff}}`$ $`=`$ $`{\displaystyle _{𝐫\tau }}\{{\displaystyle \frac{1}{\pi ^2𝒦_0}}(_\tau \mathrm{\Theta }_0)^2+{\displaystyle \frac{1}{\pi ^2𝒦_\pi }}[(_\tau \mathrm{\Theta }_1)^2+(_\tau \mathrm{\Theta }_2)^2]`$ (132)
$`+`$ $`{\displaystyle \frac{A_{}}{2}}(\mathrm{\Theta }_0)^2+{\displaystyle \frac{A_{}}{2}}(_x\mathrm{\Theta }_1)^2+{\displaystyle \frac{A_+}{2}}(_y\mathrm{\Theta }_1)^2`$
$`+`$ $`{\displaystyle \frac{A_+}{2}}(_x\mathrm{\Theta }_2)^2+{\displaystyle \frac{A_{}}{2}}(_y\mathrm{\Theta }_2)^2\},`$
where $`𝒦_\pi =𝒦(\pi ,0)`$, $`A_+=8𝒰(\pi ,0)/\pi ^2`$, $`A_{}w_{22}\mathrm{\Lambda }^{\mathrm{\Delta }_{2,2}^v}w_{2,2}^{2/(2\mathrm{\Delta }_{2,2}^v)}[𝒰(\pi ,0)]^{\mathrm{\Delta }_{2,2}^v/(2\mathrm{\Delta }_{2,2}^v)}`$, where $`\mathrm{\Lambda }`$ is the reduced low energy momentum cutoff. Note that we are not rescaling any fields or coordinates in this schematic RG treatment, just taking into account the fluctuation corrections to $`w_{2,2}`$. Thus the low-energy continuous fluctuations around the minima are described by three massless fields $`\mathrm{\Theta }_j`$ ($`j=0,1,2`$), which unlike $`\vartheta `$ in the absence of $`w_{2,2}`$ have an ordinary “relativistic” dispersion. In addition to these continuous variables, the discrete degeneracy of distinct minima allowed by periodicity of the cosine is indexed by the integers $`a,b`$.

Having taken into account $`w_{2,2}`$ already in Eq. (132), we next include the effects of a weak, renormalized $`w_{4,1}`$ coupling, which after renormalization is of order $`w_{4,1}^Rw_{4,1}\mathrm{\Lambda }^{\mathrm{\Delta }_{4,1}^v}w_{4,1}(w_{2,2}/𝒰(\pi ,0))^{\mathrm{\Delta }_{4,1}^v/(2\mathrm{\Delta }_{2,2}^v)}`$. If all the bare non-linear couplings are small, then after renormalization, this will be the largest remaining term, which is why we treat it next. Inserting the decomposition of Eq. (131) into the four vortex hopping term, one obtains

$$𝒮_1^{\mathrm{eff}}=w_{4,1}^R_{𝐫\tau }\left[\mathrm{cos}8\mathrm{\Theta }_1+\mathrm{cos}8\mathrm{\Theta }_2\right],$$
(133)
assuming small gradients of $`\mathrm{\Theta }_{1/2}`$, as mandated by Eq. (132). Eqs. (132,133) describe a 3d “height” model for $`\mathrm{\Theta }_1,\mathrm{\Theta }_2`$. As is well known in such models, the fluctuations of the free scalar fields $`\mathrm{\Theta }_1,\mathrm{\Theta }_2`$ are bounded, so that the cosine terms in Eq. 133 are always relevant and pin these fields to integer multiples of $`\pi /4`$, which we can take, without loss of generality, to be zero. Thus the net effect of $`w_{4,1}`$ is to leave only $`\mathrm{\Theta }_0`$ as a low-energy mode.

Setting therefore $`\mathrm{\Theta }_1=\mathrm{\Theta }_2=0`$ in the pinned (“flat”) phase, we have

$$e^{2i\vartheta _{x,y}}=e^{2i\vartheta _0}(1)^{ax+by}.$$
(134)
Here, since the fluctuations of $`\mathrm{\Theta }_0`$ are likewise bounded, we have replaced it by its average (zero mode) value $`\mathrm{\Theta }_0\vartheta _\mathrm{𝟎}`$. This average can be determined by minimizing the on-site Lagrangian in Eq. (124). In principle, the parameters in Eq. (124) should also be renormalized, e.g. $`\upsilon _{2q}\stackrel{~}{\upsilon }_{2q}\upsilon _{2q}\mathrm{exp}[q^2c_q\mathrm{ln}^2(w_{2,2}/𝒰(\pi ,0))]`$. This gives,

$$e^{2i\vartheta _\mathrm{𝟎}}=(1)^{ab},$$
(135)
for $`\stackrel{~}{\upsilon }_4>\stackrel{~}{\upsilon }_2/8`$ (with $`\stackrel{~}{\upsilon }_2>0`$ assumed) and

$$e^{2i\vartheta _\mathrm{𝟎}}=(1)^{ab}[(\stackrel{~}{\upsilon }_2/8|\stackrel{~}{\upsilon }_4|)\pm i\sqrt{1(\stackrel{~}{\upsilon }_2/8\stackrel{~}{\upsilon }_4)^2}],$$
(136)
for $`\stackrel{~}{\upsilon }_4<\stackrel{~}{\upsilon }_2/8`$.

The spatial orderings that are implied by these mean field solutions follow readily from the expressions in Eq. (9,10) relating the bare boson density and plaquette energy density to the “low-energy” field $`\vartheta `$. In the former case above, the boson density is uniform since $`\mathrm{sin}(2q\vartheta _{x,y})=0`$, whereas there is a plaquette energy density wave with,

$$ϵ_{x,y}c_2^ϵ(1)^{(x+a)(y+b)}.$$
(137)
The four states with $`a,b=0,1`$ correspond to the four plaquette density wave states of the original lattice Boson model in which one out of every four plaquettes is resonating more strongly.

When $`\stackrel{~}{\upsilon }_4<\stackrel{~}{\upsilon }_2/8`$, on the other hand, since $`\mathrm{sin}(2\vartheta )0`$ the plaquette energy density wave states (with reduced order $`c_2^ϵ(\stackrel{~}{\upsilon }_2/8|\stackrel{~}{\upsilon }_4|)c_2^ϵ`$) are co-existing with a charge-density-wave (CDW) state ordered at $`(\pi ,\pi )`$:

$$\delta n_{x+\frac{1}{2},y+\frac{1}{2}}\pm c_2^\rho \sqrt{1(\stackrel{~}{\upsilon }_2/8\stackrel{~}{\upsilon }_4)^2}(1)^{x+y}.$$
(138)
In the limit in which $`\stackrel{~}{\upsilon }_4\mathrm{}`$, the amplitude of the plaquette density wave vanishes, leaving only CDW order. We believe that an absence of a pure CDW state more generally is an artifact of the mean-field treatment being employed, which ignores all fluctuations in the $`\vartheta `$ field. In particular, one can imagine domain walls forming between domains of different plaquette-density-wave order, which cost an energy $`E_{wall}w_{2,2}(\stackrel{~}{\upsilon }_2/\stackrel{~}{\upsilon }_4)^2`$ for small $`w_{2,2}`$ and large $`|\stackrel{~}{\upsilon }_4|`$. When this energy is small, fluctuations will undoubtedly disorder the plaquette-density-wave order, leaving the pure CDW.

Our quantum Monte Carlo simulations on the boson ring model in the dual representation reveal that in the large $`U`$ limit with $`U/K>(U/K)_c2.5`$, the Bose metal phase is unstable to the formation of a CDW state. This is apparent in Fig. 6 which shows that the density structure factor peak at $`(\pi ,\pi )`$ grows as $`L^2`$, indicating the presence of long-ranged CDW order. This order is, however, weak, with a staggered magnetization $`m0.1`$ (see Fig. 6)), much smaller than the classical value $`m=0.5`$. It is also instructive to examine the vortex two-point function in the CDW state, which is shown in Fig. 7. Notice that the vortex correlation function is long-ranged indicating a “vortex condensation” ($`e^{i\theta }0`$), as expected in such a conventional insulating state.

Both the plaquette and charge density wave states will be insulators, with a charge gap. This follows since in both cases the field $`\vartheta `$ is “pinned” by the cosine potentials and is not fluctuating. Since the bare boson density is $`n=\mathrm{\Delta }_{xy}\theta /\pi `$, adding a particle at the origin can be achieved by shifting $`\theta _{xy}\theta _{xy}+\pi `$ for all $`x,y>0`$. This will cost a finite amount of energy (coming from the plaquette at the origin) when the field $`\vartheta `$ is pinned.

## VI Exact Wavefunction

In this Section we consider a variant of the boson ring Hamiltonian which allows us to obtain an exact zero energy wavefunction when the couplings are carefully tuned. We then perturb away from the soluble point, using the exact wavefunction to compute properties of the adjacent quantum phase. Specifically, we find a translationally invariant fluid phase with a finite compressibility, behavior consistent with both a superfluid and the Bose metal. But a calculation of the boson tunneling gap in a finite system shows a $`1/L`$ scaling as expected for the Bose metal, and inconsistent with the $`1/L^2`$ dependence of a 2d superfluid. Thus, we confidently conclude that the Bose metal phase exists over a finite portion of the phase diagram adjacent to the soluble point.

The model we consider is a “hard-core” version of the boson ring model, in which only zero or one boson is allowed per site, $`b^{}b=0,1`$. In the usual way we represent the hard-core bosons as Pauli matrices,

$`b^{}`$ $`=`$ $`\sigma ^+,`$ (139)
$`b`$ $`=`$ $`\sigma ^{},`$ (140)
$`n=b^{}b`$ $`=`$ $`{\displaystyle \frac{1}{2}}(1+\sigma ^z),`$ (141)
where $`\stackrel{}{\sigma }`$ is the standard vector of Pauli matrices, and $`\sigma ^\pm =\frac{1}{2}(\sigma ^x\pm i\sigma ^y)`$. The Hamiltonian we consider consists of two terms:

$`H_{1/2}`$ $`=`$ $`{\displaystyle \underset{xy}{}}\{J_4(\sigma _{x,y}^+\sigma _{x+1,y}^{}\sigma _{x+1,y+1}^+\sigma _{x,y+1}^{}+\mathrm{h}.\mathrm{c}.)`$ (142)
$`+u_4\widehat{P}_{\mathrm{flip}}(x,y)\},`$
where

$$\widehat{P}_{\mathrm{flip}}(x,y)=\frac{1}{16}\underset{\alpha =\pm 1}{}\underset{\beta ,\gamma =\pm 1}{}\left(1+\alpha \beta \gamma \sigma _{x+\frac{1}{2}+\frac{\beta }{2},y+\frac{1}{2}+\frac{\gamma }{2}}^z\right).$$
(143)
The first term with $`J_4>0`$ in Eq. (142) is the hard-core analog of the ring hopping term proportional to $`K`$ in the rotor model, Eq. (1). The operator $`\widehat{P}_{\mathrm{flip}}(x,y)`$ is a projection operator onto the two flippable configurations of the square plaquette whose lower-left corner is at the site $`(x,y)`$. For $`u_4>0`$, this term competes with the ring term by disfavoring configurations with flippable plaquettes.

Remarkably, the ground state of $`H_{1/2}`$ can be found exactly for the special Rokhsar-Kivelson (RK) point $`u_4=J_4`$, following a general construction in the spirit of the Rokhsar-Kivelson point of the square lattice quantum dimer model,RK  and employed more recently in Ref. BGF,  for a similar spin model on the Kagome lattice. The solution can be seen by rewriting $`H_{1/2}`$ as

$$H_{1/2}=\underset{xy}{}\widehat{P}_{\mathrm{flip}}(x,y)\{J_4\left(1\underset{u,v=0,1}{}\sigma _{x+u,y+v}^x\right)v\},$$
(144)
where $`v=J_4u_4`$ measures the deviation from RK point. For $`v=0`$, an obvious (zero energy) ground state of $`H_{1/2}`$ is the fully-polarized state

$$|0=\underset{x,y}{}|\sigma _{xy}^x=1.$$
(145)
A useful alternate representation follows by rewriting $`|\sigma ^x=1=\frac{1}{\sqrt{2}}(|\sigma ^z=1+|\sigma ^z=1)`$, and expanding out the direct product,

$$|0=\frac{1}{2^{N/2}}\underset{\{\sigma _{xy}=\pm 1\}}{}\underset{xy}{}|\sigma _{xy}^z=\sigma _{xy},$$
(146)
which demonstrates that $`|0`$ is a uniform real superposition of all configurations in the $`S^z`$ (boson number) basis.

This state, however, has uncertain $`\sigma ^z`$ (boson number), and so can be decomposed into many distinct ground states by projection. In particular, at a given average $`\sigma _{xy}^z=2n_{xy}1=m`$, we may project onto spatially uniform states according to

$$|0;m=\frac{1}{\sqrt{\stackrel{~}{Z}}}\underset{x,y=1}{\overset{L}{}}\widehat{P}_x(m)\widehat{P}_y(m)|0,$$
(147)
where $`\stackrel{~}{Z}`$ is a normalization constant, since the row/column projection operators,

$`\widehat{P}_x(m)`$ $`=`$ $`{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\theta }{2\pi }}\mathrm{exp}\left[i\theta _x{\displaystyle \underset{y}{}}\left(\sigma _{xy}^zm\right)\right],`$
$`\widehat{P}_y(m)`$ $`=`$ $`{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\varphi }{2\pi }}\mathrm{exp}\left[i\varphi _y{\displaystyle \underset{x}{}}\left(\sigma _{xy}^zm\right)\right],`$ (148)
commute with $`H_{1/2}`$ and amongst themselves, reflecting the conservation of boson number on each row and column by the ring dynamics.

Of course, other non-uniform ground states may be obtained by choosing $`m`$ differently on different rows and columns. This vast degeneracy signals a pathology of the RK point, which is in fact at the boundary of a first-order transition to a phase-separated “frozen” regime, for $`v<0`$. We therefore focus instead on the behavior for infinitesimal $`v>0`$, which splits this degeneracy. For boson densities near half-filling, i.e. $`|m|1/2`$, we expect that the uniform states will be favored energetically as $`v`$ is increased to slightly positive values, since these states have in this density regime more flippable plaquettes (see below). For larger $`|m|`$, this assumption is certainly violated, however, since e.g. at very low boson density the ring moves clearly do not connect all possible configurations. At $`m=0`$, however, the set of uniform configurations does form a single ergodic component under the ring move (as can be straightforwardly shown numericallyunpublished  and probably argued analytically). So we expect the set of uniform states to be an adequate description for small $`|m|`$.

To make contact with the Bose metal fixed point description, we calculate some energetic properties for infinitesimal positive $`v`$ using first-order perturbation theory. From the splitting of the different projected states at $`O(v)`$, we calculate two quantities: (1) the ground state energy density $`ϵ(m)=E(m)/L^2`$ as a function of boson density $`n=(m+1)/2`$, from which we obtain the compressibility $`\kappa (m)`$, and (2) the “single-particle” gap $`\mathrm{\Delta }_1(m,L)`$ for a finite-size $`L\times L`$ system, essentially the addition energy for a boson onto a particular row and column in the grand canonical ensemble, defined more precisely below.

Consider first the ground state energy density. For the projected state at magnetization $`m`$, the first order energy shift is

$`E(m)`$ $`=`$ $`v{\displaystyle \underset{xy}{}}0;m|\widehat{P}_{\mathrm{flip}}(x,y)|0;m`$ (149)
$`={\displaystyle \frac{v}{\stackrel{~}{Z}}}{\displaystyle \underset{xy}{}}0|\widehat{P}_{\mathrm{flip}}(x,y){\displaystyle \underset{x,y=1}{\overset{L}{}}}\widehat{P}_x(m)\widehat{P}_y(m)|0,`$
where the latter equality follows from the fact that $`\widehat{P}_{\mathrm{flip}}`$ is diagonal in the $`\sigma ^z`$ basis and hence commutes with the projection operators. Using Eq. (145), one readily sees that the energy shift can be rewritten in a form reminiscent of classical statistical mechanics,

$$E(m)=v\underset{xy}{}P_{\mathrm{flip}}(x,y;\{\sigma \})_{\{\sigma _{xy}\}}$$
(150)
where

$$𝒪_{\{\sigma _{xy}\}}=\frac{1}{Z}\underset{\{\sigma _{xy}=\pm 1\}}{}𝒪\underset{x,y=1}{\overset{L}{}}P_x(m;\{\sigma \})P_y(m;\{\sigma \}).$$
(151)
Here $`P_{\mathrm{flip}}`$, $`P_x`$ and $`P_y`$ are the classical functions obtained from $`\widehat{P}_{\mathrm{flip}}`$, $`\widehat{P}_x`$, $`\widehat{P}_y`$, respectively, by replacing the operator $`\sigma _{xy}^z`$ by $`\sigma _{xy}`$, and $`Z`$ is chosen such that $`1_{\{\sigma _{xy}\}}=1`$.

We see that Eq. (151) simply defines an expectation value for an “infinite temperature” Ising model with a constrained magnetization on each row and column. By proceeding along similar lines, one can calculate the energy shift for a state with an extra boson on a single particular row and column, and hence obtain $`\mathrm{\Delta }_1`$. Amusingly, both problems can be solved exactly using a saddle-point technique (see Appendix E). The results demonstrate that the constraints are “nearly irrelevant”, and the lattice gas behaves nearly as its unconstrained counterpart in the large system limit. In particular, the energy density, as $`L\mathrm{}`$, becomes

$$ϵ(m)=\frac{v}{8}(1m^2)^2.$$
(152)
This can be easily understood by assuming that each site is completely independent, and that the only effect of the constraint is to determine the relative probabilities of the two spin states, to wit $`p_{}=(1+m)/2,p_{}=(1m)/2`$. On a given four-site plaquette, only the two (of 16 total) configurations in which the spin alternates around the plaquette are flippable. Hence the average flippability per plaquette is $`2p_{}^2p_{}^2=(1m^2)^2/8`$, in agreement with the result above.

The “single-particle” gap is somewhat less intuitive. To be precise and avoid ambiguities due to an unspecified chemical potential, we define

$$\mathrm{\Delta }_1(m)=\frac{1}{2}\left[E_1(m;+1)2E_1(m;0)+E_1(m,1)\right],$$
(153)
where $`E_1(m;\lambda )`$ is the ground-state energy of a state in which the number of bosons on row and column $`x=y=1`$ is increased by $`\lambda `$ relative to the the uniform state at magnetization $`m`$:

$`|0;m,\lambda `$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{Z^{}}}}\widehat{P}_{x=1}(m+2\lambda /L)\widehat{P}_{y=1}(m+2\lambda /L)`$ (154)
$`\times {\displaystyle \underset{x,y=2}{\overset{L}{}}}\widehat{P}_x(m)\widehat{P}_y(m)|0,`$
and $`E_1(m;\lambda )=v_{xy}0;m,\lambda |\widehat{P}_{\mathrm{flip}}(x,y)|0;m,\lambda `$. One finds

$$\mathrm{\Delta }_1(m)=\frac{2v}{L}(1m^2)+O(1/L^2).$$
(155)
These two quantities can be compared to the general predictions expected from our theory of the Bose metal phase. Consider first the compressibility. From Eq. (152), one has $`\mu =dϵ/dn=2dϵ/dm=vm(1m^2)`$, and hence

$$\kappa =\frac{dn}{d\mu }=\frac{1}{2v}\frac{1}{13m^2}.$$
(156)
Note that the compressibility diverges and becomes negative for $`1/\sqrt{3}<|m|<1`$, indicative of an instability of the uniform state well away from half-filling. For $`|m|<1/\sqrt{3}`$, however, $`\kappa `$ is finite, consistent with the Bose metal phase. Indeed, from our general Harmonic description of the Bose metal phase as in Eq. ??, we have $`\kappa ^1=𝒰(𝐤=0)`$.

In order to rule out a superfluid phase, we now compute the single-particle gap, $`\mathrm{\Delta }_1(L)`$, using the Bose metal Gaussian fixed point theory to show that it varies as $`1/L`$ as in Eq. (155). This gap can be extracted from the spatially local correlator,

$$G_\phi (\tau )=e^{i\phi _𝐫(\tau )}e^{i\phi _𝐫(0)},$$
(157)
evaluated in a finite $`L\times L`$ system. From the spectral representation, one has, at zero temperature,

$$\mathrm{\Delta }_1(L)=\underset{\tau \mathrm{}}{lim}\frac{\mathrm{ln}G_\phi (\tau )}{\tau }.$$
(158)
Performing the Gaussian integral using the effective action in Eq. (43), one has

$$\mathrm{ln}G_\phi (\tau )=_{\mathrm{}}^{\mathrm{}}\frac{d\omega }{2\pi }\frac{1}{L^2}\underset{𝐤}{}\frac{𝒰(𝐤)}{\omega ^2+E_𝐤^2}\left(1e^{i\omega \tau }\right),$$
(159)
where the wavevector sum is over inequivalent values with $`k_x,k_y(2\pi /L)𝒵`$ in the first Brillouin zone, and with the mode energy $`E_𝐤`$ given explicitly in Eq. (37). As $`\tau \mathrm{}`$, the $`\omega `$ integral in Eq. (159) is dominated by those terms for which $`E_𝐤=0`$. This occurs along the Bose surface, i.e. for $`k_x=0`$ or $`k_y=0`$. In the finite size system there are $`2L12L`$ such points (for $`L1`$), hence

$`\mathrm{ln}G_\varphi (\mathrm{𝟎},\tau ){\displaystyle \frac{2\overline{𝒰}}{L}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi }}{\displaystyle \frac{1}{\omega ^2}}\left(1e^{i\omega \tau }\right),`$ (160)
where

$$\overline{𝒰}=\frac{1}{2}\left[_\pi ^\pi \frac{dk_x}{2\pi }𝒰(k_x,0)+_\pi ^\pi \frac{dk_y}{2\pi }𝒰(0,k_y)\right],$$
(161)
is the Bose-surface average of $`𝒰(𝐤)`$. Thus, we finally have,

$$\mathrm{ln}G_\varphi (\mathrm{𝟎},\tau )\frac{\overline{𝒰}|\tau |}{L},$$
(162)
giving the $`1/L`$ behavior for the single-particle gap:

$$\mathrm{\Delta }_1=\frac{\overline{𝒰}}{L}.$$
(163)
In a 2d superfluid, the single-particle gap is much smaller, vanishing as $`1/L^2`$ in the large $`L`$ limit.

The agreement of the finite compressibility and $`1/L`$ scaling of $`\mathrm{\Delta }_1`$ between the soluble model and our fixed point theory of the Bose metal strongly argues that the ground state of $`H_{1/2}`$ is the Bose liquid for $`0<v1`$ and $`|m|<1/\sqrt{3}`$. If this postulate is correct, certain combinations of the Bose liquid function $`𝒰(𝐤)`$ can be obtained explicitly for the soluble model. In particular, we find

$`𝒰(𝐤=0)`$ $`=`$ $`2v(13m^2),`$ (164)
$`\overline{𝒰}`$ $`=`$ $`2v(1m^2).`$ (165)
## VII Discussion

### VII.1 Fractionalization, the $`Z_2`$ gauge theory and the high-$`T_c`$ cuprates

In this paper we have described a novel phase of quantum matter, the Bose metal, which we argue occurs in a class of square lattice boson ring models with XY symmetry. Having done so, it is reasonable to reflect upon the context in which such models are physically appropriate and the consequences perhaps observable. As discussed briefly in the introduction, the primary motivation for these models comes from the high temperature cuprate superconductors. Both the remarkably high critical temperatures and the strange behaviors in the “normal” state of these materials motivated a number of theorists early on to the radical suggestion that spin-charge separation might underly these peculiarities.Anderson ; KRS ; LN ; SenMPAF  In this picture, the electron charge, liberated both from its spin and its Fermi statistics as a “chargon” (or “antiholon”), is relatively free to Bose condense to form a superconducting state. At higher temperatures, or when the material is under-doped, the chargons and spin-carrying fermionic “spinons” form a novel fluid, which would no doubt behave very differently from a normal metal.

Subsequent theoretical and consequent experimental work has both advanced the status of these ideas and posed a severe challenge to their applicability to the cuprates. A recent formulationZ2A ; Z2B  of interacting electrons in terms of spinons and chargons minimally coupled to a $`Z_2`$ gauge field has helped to support the concept of electron fractionalization with a very concrete formal framework. The essential elements of the $`Z_2`$ gauge theory are reviewed in Appendix D. As with earlier $`U(1)`$ and $`SU(2)`$ gauge theory formulationsU1A ; U1B ; U1C , the “spinons” are taken as Fermions carrying the spin of the electron but are electrically neutral. The bosonic “chargons” carry the electrons charge. The $`Z_2`$ gauge theory provides a convenient phenomenology for describing a fractionalized phase in which the spinons and chargons are deconfined, and live as well-defined particle excitations. Because of its concrete formulation, and the virtue that the $`Z_2`$ gauge theory has a well-understood confinement-deconfinement phase transition, it is possible to calculate from it in a simple way the qualitative (universal) properties of the fractionalized state. One of the most fundamental of these properties is the existence of gapped topological excitationsvison1 ; vison2 ; vison3  called “visons”. The visons act as sources of $`Z_2`$ flux, and may be thought of as remnants of unpaired superconducting vortices in the fractionalized insulatorZ2A , which in turn is viewed as a paired-vortex condensateBFN2 . Understanding these excitations led to a rather direct proposal for a “vison trapping” experiment,Z2B  designed to trap and detect visons by cycling through the superconducting–normal transition. Current experimentsvisexp1 ; visexp2  apparently imply that the gap to such vison excitations in the insulator, if it is non-zero at all, is less than 190$`K`$. Since the presence of a vison gap is a necessary condition for the existence of a true spin-charge separated ground state, this unnaturally low energy scale presents a difficult obstacle to theories of fractionalization as applied to at least these particular cuprate materials (samples of underdoped BSSCO and YBCO).

The spin-charge separation scenario is nonetheless extremely appealing theoretically, and it is interesting to consider the possibility of retaining some degree of this physics locally. One may imagine for simplicity an undoped model which interpolates between at one extreme a conventional “Hubbard-like” insulator with gapped charge degrees of freedom and Heisenberg-interacting spins, and at the other extreme a fully fractionalized insulator. The $`Z_2`$ gauge theory is well-suited for this purpose. The interpolation between the above two limits is accomplished in this model by varying a coefficient $`K`$, which controls the strength of fluctuations in the gauge field. The deconfined phase is obtained in the large $`K`$ limit, while we consider here the first deviations away from $`K=0`$, which is deep within the confined, Hubbard-like phase. As we show in Appendix D, when the gauge theory is deep within it’s confined phase (with $`K`$ small), the gauge fields can be formally integrated out and one recovers a Hamiltonian expressible in terms of electron operators and composites built from the electron such as the spin operator, together with a Cooper pair field. In the spin sector, we find that the leading terms obtained in this limit are an antiferromagnetic nearest-neighbor Heisenberg exchange (for $`K=0`$) and a plaquette ring term,

$$H_{\mathrm{}}^s=J_{\mathrm{}}^s\underset{\mathrm{}}{}_{\mathrm{}}^s,$$
(166)
with $`J_{\mathrm{}}^sO(K)`$ ($`J_{\mathrm{}}^s=K(t_s^4+\mathrm{\Delta }^4)/(h+U)^4`$, in terms of the parameters in Appendix D). Here, $`_{\mathrm{}}^s`$ is a plaquette ring operator defined in terms of the four spins on a given plaquetteoldring1 ; oldring2 :

$`_{\mathrm{}}^s`$ $`=`$ $`{\displaystyle \underset{i<j=1}{\overset{4}{}}}𝐒_i𝐒_j+4(𝐒_1𝐒_2)(𝐒_3𝐒_4)`$ (167)
$`+`$ $`4(𝐒_1𝐒_4)(𝐒_2𝐒_3)4(𝐒_1𝐒_3)(𝐒_2𝐒_4).`$
In the charge sector, to $`O(K)`$ one finds the leading term

$$H_{\mathrm{}}^c=J_{\mathrm{}}^c\underset{\mathrm{}}{}_{\mathrm{}}^c,$$
(168)
with $`J_{\mathrm{}}^cO(K)`$ ($`J_{\mathrm{}}^c=Kt_c^4/(h+U)^4`$, derived in Appendix D). Here, $`_{\mathrm{}}^c`$ is a ring operator in the charge sector expressed simply in terms of Cooper pair operators, $`e^{i\varphi _𝐫}b_𝐫^2`$ as,

$$_{\mathrm{}}^c=\underset{i<j=1}{\overset{4}{}}\mathrm{cos}(\varphi _i\varphi _j)+\mathrm{cos}(\mathrm{\Delta }_{xy}\varphi _𝐫).$$
(169)
Notice that the second term above is precisely of the form taken in our starting Hamiltonian. As $`K`$ is increased from zero, such ring-exchange processes play more and more important roles. We thus view the ring Hamiltonian studied in this paper as a suitable model that offers an intermediate ground between the fully spin-charge separated scenario and conventional phases of matter.

Should this bosonic ring model for the charge sector have some relevance for the cuprates, what might be the consequences and interpretation? As we have seen in Sec. V, the ring model sustains the Bose metal phase for sufficiently repulsive interactions (large $`U`$) and incommensurate densities. Indeed, as we saw in Sec. V, a simplistic estimate indicates that there is no stable regime for the Bose metal at half-filling. Hence, if we assume the strong interaction condition obtains, one would expect an evolution from a conventional insulator with a charge gap at half-filling to the Bose metal upon doping sufficiently away. If one further presumes (as seems natural) that the generalized “stiffness” $`𝒦`$ increases with doping, then one would expect further doping to lead to a superconducting (boson hopping) instability. This scenario thus naturally associates the Bose metal with the pseudo-gap regime of the high-$`T_c`$ cuprates.

### VII.2 Extensions

To determine if there is any truth to the above scenario requires considerable extensions of the present work. Most significantly, spin and quasiparticle degrees of freedom are important components of the high-$`T_c`$ materials, and should be incorporated into the description. It will be interesting to consider interactions of the Fermionic quasiparticles with the strongly fluctuating collective modes of the Bose metal. It seems to us quite likely that the Bose metal can remain stable in the presence of these interactions, while perhaps at the same time producing rather strong modifications of the fermionic degrees of freedom. In any case, it should be possible to consider photoemission spectra and local electron tunneling density of states in a model with the quasiparticles coupled to the Bose metallic modes.

We have also left a number of issues within the purely bosonic description unanswered. To understand transport measurements generally will require an understanding of the effects of disorder. Several experiments attempting to access the “normal” state of the cuprates employ large local or uniform magnetic fields to suppress superconductivityTaillefer . Such magnetic fields indeed tend to suppress the superconducting instabilities of the Bose metal, and one should explore the ramifications of a Bose metallic description for e.g. transport measurements in high fields.

There are a number of other possible applications of the present formalism to frustrated magnets with appreciable ring exchangeringnum . Magnetic ring exchange processes are believed important for instance in the Wigner crystal phase of the two-dimensional electron gasUCLA ; Ceperley  at $`r_s40`$. The effect of such processes are generally difficult to analyze, but a perhaps illuminating approach may be to consider the easy-plane limit of such spin ring models. A straightforward calculation shows that such a limit recovers XY ring models of the sort studied here (representing the spins by hard-core bosons). The plaquette duality constructed here for the square lattice can be straightforwardly generalized to other lattices, such as the triangular and Kagome cases. For the bosonic Kagome case, one can directly in this way construct a dual formulation from which it is straightforward to construct the visons and spinons of Ref. BGF, . Whether the spinons are confined or deconfined in the pure XY ring model on this lattice can be easily established numerically using the dual $`\theta `$ variables, and this investigation is underway.hermele  On the triangular lattice with four-site (parallelogram) ring exchange, the plaquette duality shows that there is no stable metallic state in this case, but a numerical approach will be necessary to determine the ultimate nature of the ground state.

On the purely theoretical side, we have introduced the mathematics necessary for a $`2+1`$ dimensional “bosonization” scheme, applied here to “bosonize bosons” in terms of collective $`\phi ,\vartheta `$ modes. We expect it should be possible to bosonize fermionic ring models in a similar way. Consider spinless fermions $`\{f_i,f_j^{}\}=\delta _{ij}`$ (not at this point to be associated with cuprate quasiparticles) living on the sites of the original square lattice. One may formally define Jordan-Wigner hard-core bosons by introducing a “string” $`\widehat{S}`$ winding around the lattice,

$$f_𝐫=b_𝐫\widehat{S}_𝐫,$$
(170)
where

$$\widehat{S}_{x+\frac{1}{2},y+\frac{1}{2}}=\underset{x^{}<x}{}(1)^{n_{x^{}+\frac{1}{2},y+\frac{1}{2}}}\underset{y^{}>y}{}\underset{x^{}}{}(1)^{n_{x^{}+\frac{1}{2},y^{}+\frac{1}{2}}}.$$
(171)
With this definition, the boson operators obey proper canonical commutation relations at different sites, and have a hard-core interaction on the same site. For an ordinary dynamics, this extremely non-local string presents unsolvable difficulties for analytic treatment. However, using the exact plaquette duality appropriate for ring dynamics, one finds that one can write

$$\widehat{S}_{x+\frac{1}{2},y+\frac{1}{2}}=e^{i(\theta _{x,y+1}\theta _{x,y})}.$$
(172)
Here we have ignored any possible boundary terms. Up to this proviso, the string becomes local in the $`\theta `$ variables! Thus, in a rotor representation with $`b=e^{i\varphi }`$ (this requires a large $`Un(n1)`$ term to maintain the hard-core constraint), one has

$$f_{x+\frac{1}{2},y+\frac{1}{2}}=e^{i\varphi _{x+\frac{1}{2},y+\frac{1}{2}}}e^{i(\theta _{x,y+1}\theta _{x,y})},$$
(173)
which is directly analogous to the bosonization formula for spinless fermions in one dimension! This approach may lead to an understanding of novel non-Fermi liquid states in fermionic ring models.

Clearly there are a range of applications to be explored based on the present work. We hope that some of these may ultimately enhance our understanding of experimentally accessible strongly correlated materials.

###### Acknowledgements.

We are grateful to Anders Sandvik, Doug Scalapino, and T. Senthil for sparkling conversations. L.B. and A.P. were supported by NSF grant DMR-9985255, and the Sloan and Packard foundations, M.P.A.F. by NSF Grants DMR-9704005 and PHY-9907949.

## Appendix A Trotter decomposition

The usual Trotter decomposition of $`Z_{\mathrm{}}`$ gives

$$Z=\underset{\{\theta _{𝐫\tau }\}}{}\underset{\tau =0}{\overset{\tau =\beta }{}}\theta _{𝐫\tau +ϵ}|e^{ϵK_𝐫\mathrm{cos}(\pi N_𝐫)}e^{\frac{ϵU}{2\pi ^2}_𝐫[\mathrm{\Delta }_{xy}\theta _𝐫\pi \overline{n}]^2}|\theta _{𝐫\tau },$$
(174)
where we have used $`ϵ1`$ to separate the imaginary-time evolution operator into to factors. The second can be directly evaluated to give

$$Z=\underset{\{\theta _{𝐫\tau }\}}{}\underset{\tau =0}{\overset{\tau =\beta }{}}e^{\frac{ϵU}{2\pi ^2}_𝐫[\mathrm{\Delta }_{xy}\theta _{𝐫\tau }\pi \overline{n}]^2}\underset{𝐫}{}\theta _{𝐫\tau +ϵ}|e^{ϵK\mathrm{cos}(\pi N_𝐫)}|\theta _{𝐫\tau }.$$
(175)
The latter matrix element can be written

$`\theta ^{}|e^{ϵK\mathrm{cos}(\pi N_𝐫)}|\theta `$ $`=`$ $`\theta ^{}|1+ϵK\mathrm{cos}(\pi N_𝐫)|\theta `$ (176)
$`=`$ $`\delta _{\theta ^{},\theta }+{\displaystyle \frac{ϵK}{2}}\left(\delta _{\theta ^{},\theta +\pi }+\delta _{\theta ^{},\theta \pi }\right)`$
$`=`$ $`\mathrm{exp}\left[{\displaystyle \frac{1}{\pi ^2}}\mathrm{ln}\left({\displaystyle \frac{ϵK}{2}}\right)\left(\theta ^{}\theta \right)^2\right],`$
correct to $`O(ϵK)`$. This leads directly to Eq. (22)

## Appendix B Asymptotics

In this appendix we calculate the asymptotics of the correlator $`c_{xy}(\tau ,\tau _1,\tau _2)`$, introduced in Sec. IV. Using the general rules of Gaussian theories, one has

$`c_{x,y}(\tau ,\tau _1,\tau _2)`$ $`=`$ $`2\left(\vartheta _{xy}(\tau )+\vartheta _{00}(0)\vartheta _{x0}(\tau _1)\vartheta _{0y}(\tau _2)\right)^2_0`$ (177)
$`=2\pi ^2\mathrm{Re}{\displaystyle _𝐤}𝒦(𝐤){\displaystyle \frac{2+\left(e^{E_𝐤|\tau |}+e^{E_𝐤|\tau _1\tau _2|}\right)e^{i(k_xx+k_yy)}\left(e^{E_𝐤|\tau _1|}+e^{E_𝐤|\tau \tau _2|}\right)e^{ik_xx}\left(e^{E_𝐤|\tau _2|}+e^{E_𝐤|\tau \tau _1|}\right)e^{ik_yy}}{E_𝐤}}.`$
The above integral was considered earlier for the special case $`\tau =\tau _1=\tau _2=0`$, for which it was shown to grow logarithmically with $`x`$ at fixed $`y`$ (and vice versa). This growth with $`x`$ and $`y`$ persists when $`\mathrm{\Omega }_0|\tau |,\mathrm{\Omega }_0|\tau _1|,\mathrm{\Omega }_0|\tau _2||x|`$ (or $`|y|`$). To see this, we differentiate

$`_xc_{xy}`$ $`=`$ $`2\pi ^2\mathrm{Re}{\displaystyle _𝐤}{\displaystyle \frac{ik_xe^{ik_xx}𝒦(𝐤)}{E_𝐤}}\left[e^{E_𝐤|\tau _1|}+e^{E_𝐤|\tau \tau _2|}e^{ik_yy}\left(e^{E_𝐤|\tau |}+e^{E_𝐤|\tau _1\tau _2|}\right)\right]`$ (178)
$`2{\displaystyle _0^\pi }𝑑k_y{\displaystyle \frac{𝒦(0,k_y)}{\omega (k_y)}}\left[{\displaystyle \frac{x}{x^2+\omega ^2(k_y)\tau _1^2}}+{\displaystyle \frac{x}{x^2+\omega ^2(k_y)(\tau \tau _2)^2}}\mathrm{cos}k_yy\left({\displaystyle \frac{x}{x^2+\omega ^2(k_y)\tau ^2}}+{\displaystyle \frac{x}{x^2+\omega ^2(k_y)(\tau _1\tau _2)^2}}\right)\right],`$
where $`\omega (k_y)=\mathrm{\Omega }(0,k_y)|2\mathrm{sin}(k_y/2)|`$, and in the second line we assumed (as valid for large $`|x|`$) the $`k_x`$ integral is dominated by small $`k_x`$ and integrated it in that regime. Hence for $`|x|\mathrm{\Omega }_0|\tau |_1,\mathrm{\Omega }_0|\tau |`$, one has

$$c_{xy}4\eta _v(y)\mathrm{ln}|x|+f(\tau ,\tau _1,\tau _2,y),$$
(179)
where $`f`$ is an unknown (but $`x`$-independent) function of $`y`$ and the imaginary times. We expect that $`c_{xy}`$ should cross over to a $`\mathrm{ln}^2|\mathrm{\Omega }_0\tau _i|`$ behavior if any of the $`\mathrm{\Omega }_0|\tau _i|`$ become large compared to the spatial distances involved. However, it is unclear without further calculation whether this crossover occurs when the imaginary time argument becomes large compared to the smaller (e.g. $`y`$ above) or the larger (e.g. $`x`$ above) of the two spatial coordinates. To ascertain this information, we differentiate with respect to the largest time (which we choose positive without loss of generality), e.g. for $`\tau >|\tau _1|,|\tau _2|`$,

$`_\tau c_{xy}`$ $`=`$ $`2\pi ^2{\displaystyle _𝐤}𝒦(𝐤)[e^{E_𝐤(\tau \tau _1)ik_yy}`$ (180)
$`+e^{E_𝐤(\tau \tau _2)ik_xx}e^{E_𝐤\tau i(k_xx+k_yy)}]`$
$``$ $`2\pi ^2{\displaystyle _𝐤}𝒦(𝐤)e^{E_𝐤(\tau \tau _1)}e^{ik_yy},`$
where in last line we have kept the dominant term for $`|x||y|,\mathrm{\Omega }_0\tau `$. If we assume $`|y|O(1)`$, then for $`\mathrm{\Omega }_0(\tau \tau _1)1`$, the integral is dominated by $`𝓀0`$ (due to the logarithmic singularity in the collective density of states), and one may approximate $`e^{ik_yy}1`$, leading to $`_\tau c_{xy}2\sqrt{𝒦_0/𝒰_0}\mathrm{ln}|\mathrm{\Omega }_0(\tau \tau _1)|/(\tau \tau _1)`$, whence $`f\sqrt{𝒦_0/𝒰_0}\mathrm{ln}^2|\mathrm{\Omega }_0(\tau \tau _1)|`$. More generally, one can show that this behavior obtains provided $`\mathrm{\Omega }_0(\tau \tau _1)|y|`$. By identical arguments, $`f`$ also grows like $`\sqrt{𝒦_0/𝒰_0}\mathrm{ln}^2|\mathrm{\Omega }_0\tau _2|`$ for $`\mathrm{\Omega }_0|\tau _2||y|`$. Thus, the integral defining $`C_{xy}^{(4)}(\tau )`$ is dominated by $`\tau _1\tau `$ and $`\tau _20`$. Hence, we can focus on this case, and consider the behavior in general as a function of $`|x|1`$ and $`\tau 1`$, but with no particular relation between them. One finds

$`c_{xy}(\tau ,\tau ,0)`$ $``$ $`4\pi ^2\mathrm{Re}{\displaystyle _𝐤}{\displaystyle \frac{𝒦(𝐤)}{E_𝐤}}\left(1e^{ik_yy}\right)(1e^{E_𝐤|\tau |ik_xx}).`$
It is straightforward to see that this integral recovers the result, Eq. (179) for $`|x||\tau |`$, and conversely, $`c_{xy}(\tau ,\tau ,0)4\eta _v(y)\mathrm{ln}\tau `$ for $`|\tau ||x|`$. Putting all the above analysis together, one arrives at the general approximate result of Eq. (92) in Sec. IV.

## Appendix C Relevance of hopping

In this Appendix we compute the corrections to the Bose metal imaginary-time equal-space correlator,

$$G(\tau )=e^{i\phi _𝐫(\tau )}e^{i\phi _𝐫(0)},$$
(181)
to leading (quadratic) order in the Boson hopping term:

$$𝒮_t=t_y_\tau \underset{𝐫}{}\mathrm{cos}(\phi _{𝐫+\widehat{y}}\phi _𝐫).$$
(182)
Perturbing around the Gaussian Bose metal theory gives, $`G(\tau )=G_0(\tau )+G_2(\tau )+O(t_x^4)`$, with

$$G_2(\tau )=\frac{1}{2}e^{i\varphi _𝐫(\tau )i\varphi _𝐫(0)}𝒮_t^2_{0,c},$$
(183)
where the connected correlator is taken with respect to the Gaussian action. Using Wick’s theorem this can be re-expressed as, $`G_2(\tau )=t_x^2G_0(\tau )I(\tau )`$, with

$$I(\tau )=\underset{x_1,x_2,y_1}{}_{\tau _1,\tau _2}𝒟_0(x_{12},y_{12})[e^A1],$$
(184)
where the function $`A=A(x_1,x_2,y_1,\tau _1,\tau _2,\tau )`$ is given by,

$$A=_{𝐤,\omega }e^{i𝐤𝐫_1i\omega \tau _1}B(𝐤,\omega ,\tau )(1e^{ik_xx_{12}+i\omega \tau _{12}})+h.c.,$$
(185)
$$B(𝐤,\omega ,\tau )=𝒰(𝐤)\frac{(e^{i\omega \tau }1)(1e^{ik_y})}{\omega ^2+E_𝐤^2}.$$
(186)
Here, $`x_{12}=x_1x_2`$, $`\tau _{12}=\tau _1\tau _2`$ and $`𝒟_0(x,\tau )(x^2+\tau ^2)^\mathrm{\Delta }`$ is the two-point correlator for the tunneling operator evaluated in the Bose metal. To isolate the $`\tau `$ dependence of $`I(\tau )`$, we re-scale $`x_1,x_2,\tau _1,\tau _2`$ by $`\tau `$, and $`k_x`$ and $`\omega `$ by $`1/\tau `$ (leaving $`y`$ and $`k_y`$ alone). For large $`\tau `$ this gives, $`I(\tau )\tau ^{2(2\mathrm{\Delta })}I(1)`$, or equivalently,

$$\frac{G_2(\tau )}{G_0(\tau )}t_y^2I(1)\tau ^{2(2\mathrm{\Delta })}.$$
(187)
We thus conclude that provided $`I(1)`$ is finite, the large $`\tau `$ behavior of the Bose metal correlator, $`G_0(\tau )`$, is unmodified if the scaling dimension $`\mathrm{\Delta }>2`$. To show that $`I(1)`$ is finite, it suffices to expand the exponential to second order in $`A`$, which gives:

$$I(1)=8_{𝐤,\omega }𝒰(𝐤)𝒟(k_x,\omega )\frac{(1\mathrm{cos}(\omega ))(1\mathrm{cos}(k_y))}{(\omega ^2+E_𝐤^2)^2},$$
(188)
with

$$𝒟(k_x,\omega )=_{x,\tau }𝒟_0(x,\tau )(1e^{i(k_xx\omega \tau )}).$$
(189)
Provided $`\mathrm{\Delta }>1`$ the integrals over $`x`$ and $`\tau `$ converge, giving $`𝒟(k_x,\omega )(k_x^2+\omega ^2)^{\mathrm{\Delta }1}`$. Inserting this into Eqn. 188, one readily sees that the $`𝐤`$ and $`\omega `$ integrals are likewise convergent, confirming that $`I(1)`$ is finite whenever $`\mathrm{\Delta }>1`$.

## Appendix D Z2 Gauge Theory and Ring Exchange

A recent formulationZ2A  of interacting electrons in two dimensions has been developed which re-expresses the electron operator in terms of “spinons” and “chargons” which are minimally coupled to a $`Z_2`$ gauge field. As with earlier $`U(1)`$ and $`SU(2)`$ gauge theory formulations, the “spinons” are taken as Fermions carrying the spin of the electron but are electrically neutral. The bosonic “chargons” carry the electrons charge. The $`Z_2`$ gauge theory provides a convenient phenomenology for describing a fractionalized phase in which the spinons and chargons are deconfined, and live as well-defined particle excitations. The deconfined phase is most readily accessed by increasing the strength of a term in the gauge theory Hamiltonian (with coefficient $`K`$, below) which suppresses the fluctuations in the gauge field. When the gauge theory is deep within it’s confined phase (with $`K=0`$), on the other hand, the gauge fields can be formally integrated out and one recovers a Hamiltonian expressible in terms of electron operators and composites built from the electron such as the spin operator, together with a Cooper pair field. Here, we show that upon integrating out the gauge field with small but non-zero $`K`$, one generates ring-exchange terms in the effective Hamiltonian. In the spin-sector of the theory, the leading order ring term involves four-spin operators around an elementary (square) plaquette, and is the familiar quartic form shown explicitly below. But there are also ring terms generated in the charge sector, and the dominant one takes precisely the form in Eq. (1), where two Cooper pairs are destroyed on opposite corners of an elementary square plaquette, hopping to the other two corners.

Consider then a Hamiltonian version of the $`Z_2`$ gauge theory:

$`H`$ $`=`$ $`H_c+H_\sigma +H_s,`$ (190)
$`H_c`$ $`=`$ $`t_c{\displaystyle \underset{\mathrm{𝐫𝐫}^{}}{}}\sigma _{\mathrm{𝐫𝐫}^{}}^z(b_𝐫^{}b_𝐫^{}+\mathrm{h}.\mathrm{c}.)`$ (191)
$`+`$ $`U{\displaystyle \underset{𝐫}{}}\left(b_𝐫^{}b_𝐫2\overline{n}\right)^2,`$
$`H_\sigma `$ $`=`$ $`K{\displaystyle \underset{\mathrm{}}{}}{\displaystyle \underset{\mathrm{}}{}}\sigma _{\mathrm{𝐫𝐫}^{}}^zh{\displaystyle \underset{<𝐫j>}{}}\sigma _{\mathrm{𝐫𝐫}^{}}^x,`$ (192)
$`H_s`$ $`=`$ $`{\displaystyle \underset{\mathrm{𝐫𝐫}^{}}{}}\sigma _{\mathrm{𝐫𝐫}^{}}^z[t_s(f_𝐫^{}f_𝐫^{}+\mathrm{h}.\mathrm{c}.)`$ (193)
$`+`$ $`\mathrm{\Delta }_{\mathrm{𝐫𝐫}^{}}(f_𝐫f_𝐫^{}f_𝐫^{}f_𝐫+\mathrm{h}.\mathrm{c}.)],`$
where $`𝐫`$ and $`𝐫^{}`$ denote sights of a 2d square lattice. Here $`b_𝐫^{}`$ creates a chargon at site $`𝐫`$ while $`f_{𝐫\alpha }^{}`$ creates a spinon with spin $`\alpha =,`$ at site $`𝐫`$. The operator $`b_𝐫^{}b_𝐫`$ measures the number of charge $`e`$ bosons at site $`𝐫`$, and the Hubbard-like $`U`$ term sets the mean Cooper pair density, $`n_𝐫b_𝐫^{}b_𝐫/2`$, to be $`\overline{n}`$. With this convention, the half-filled electron model with charge $`e`$ per site corresponds to $`\overline{n}=1/2`$. The constant $`\mathrm{\Delta }_{\mathrm{𝐫𝐫}^{}}`$ contains the information about the pairing symmetry of the spinons. The $`\sigma _{\mathrm{𝐫𝐫}^{}}^z,\sigma _{\mathrm{𝐫𝐫}^{}}^x`$ are Pauli spin matrices that are defined on the links of the lattice. The full Hamiltonian is invariant under the $`Z_2`$ gauge transformation $`b_𝐫b_𝐫,f_𝐫f_𝐫`$ at any site $`𝐫`$ of the lattice accompanied by letting $`\sigma _{\mathrm{𝐫𝐫}^{}}^z\sigma _{\mathrm{𝐫𝐫}^{}}^z`$ on all the links connected to that site. This Hamiltonian must be supplemented with the constraint equation

$$𝒢_𝐫=\underset{𝐫^{}𝐫}{}\sigma _{\mathrm{𝐫𝐫}^{}}^xe^{i\pi \left(f_𝐫^{}f_𝐫^{}+b_𝐫^{}b_𝐫^{}\right)}=1.$$
(194)
Here the product over $`\sigma _{\mathrm{𝐫𝐫}^{}}^x`$ is over all links that emanate from site $`𝐫`$. The operator $`𝒢_𝐫`$, which commutes with the full Hamiltonian, is the generator of the local $`Z_2`$ gauge symmetry. Thus the constraint $`𝒢_𝐫=1`$ simply expresses the condition that the physical states in the Hilbert space are those that are gauge invariant.

Generally, the gauge field dynamics is very complicated, but it simplifies considerably when either $`h`$ or $`K`$ greatly exceeds the strength of the couplings to the matter fields, $`t_c,t_s`$ and $`\mathrm{\Delta }_{\mathrm{𝐫𝐫}^{}}`$. For example, when $`K\mathrm{}`$, fluctuations in $`\sigma _{\mathrm{𝐫𝐫}^{}}^z`$ are suppressed completely, and one can choose a gauge with the $`\sigma _{\mathrm{𝐫𝐫}^{}}^z=1`$ on every link. This corresponds to a deconfined phase in which the spinons and chargons exist as bona fide particle excitations. On the other hand, when $`h\mathrm{}`$ one has $`\sigma _{\mathrm{𝐫𝐫}^{}}^x1`$, and it is convenient to integrate out the gauge field perturbatively. Focusing for simplicity on the Mott insulator at half-filling, and in the latter large $`h`$ limit, with $`\overline{n}=1/2`$ with $`U\mathrm{}`$, at lowest order in $`t_s`$ and $`\mathrm{\Delta }`$ one recovers the Heisenberg spin-model,

$$H_2=J_2\underset{<\mathrm{𝐫𝐫}^{}>}{}𝐒_𝐫𝐒_𝐫^{},$$
(195)
with $`J_2=(t_s^2+\mathrm{\Delta }^2)/h`$. Here, $`𝐒_𝐫\frac{1}{2}f_{𝐫\alpha }^{}𝐬_{\alpha \beta }f_{𝐫\beta }`$, with $`𝐬`$ a vector of Pauli spin matrices. Note that in the $`h,U\mathrm{}`$ limit the gauge constraint becomes a single-occupancy constraint: $`e^{i\pi f_𝐫^{}f_𝐫^{}}=1`$. Upon inclusion of a small non-zero coupling $`K`$, one generates the spin ring exchange term in Eq. (166), upon integrating out the gauge field, with $`J_{\mathrm{}}^s=K(t_s^4+\mathrm{\Delta }^4)/(h+U)^4`$. Notice that the strength of the 4-spin ring exchange interaction, $`J_{\mathrm{}}^s`$, is proportional to the gauge theory coupling $`K`$. This suggests that spin models with appreciable ring exchange interactions are good candidate models to exhibit fractionalized phases.

It is very instructive to examine the charge sector of the theory one generates upon integrating out the gauge field by similarly expanding perturbatively in small $`t_c`$. Generally, one will generate various Cooper pair hopping processes. To leading order in $`K`$ the Cooper pair plaquette ring term of Eq. (168) is generated, with $`J_{\mathrm{}}^c=Kt_c^4/(h+U)^4`$.

## Appendix E RK Manipulations

We begin by inserting the integral representations of $`P_{x/y}`$ and the explicit representation of $`P_{\mathrm{flip}}`$, inherited from Eqs. (148) and Eq. (143), respectively. The sums over $`\{\sigma _{xy}\}`$ can be explicitly performed, and we find

$$P_{\mathrm{flip}}(x,y;\{\sigma \})_{\{\sigma _{xy}\}}=F_{xy}[\{\theta ,\varphi \}]_{\{\theta _x,\varphi _y\}},$$
(196)
where

$$𝒪_{\{\theta _x,\varphi _y\}}=\frac{1}{Z_{\theta \varphi }}\underset{a=1}{\overset{L}{}}_0^{2\pi }\left(\frac{d\theta _a}{2\pi }\frac{d\varphi _a}{2\pi }\right)𝒪e^{s[\{\theta _x,\varphi _y\}]}.$$
(197)
In these expressions, the transformed “flippability” function is

$$F_{xy}[\{\theta ,\varphi \}]=\frac{1}{8}\underset{u,v=0,1}{}\mathrm{sec}(\theta _{x+u}\varphi _{y+v}),$$
(198)
and the classical “action” is given by

$`s[\{\theta _x,\varphi _y\}]`$ $`=`$ $`{\displaystyle \underset{x,y=1}{\overset{L}{}}}\mathrm{ln}\mathrm{cos}(\theta _x\varphi _y)+imL{\displaystyle \underset{a=1}{\overset{L}{}}}(\theta _a\varphi _a)`$ (199)
$`+{\displaystyle \frac{M^2}{2}}\left({\displaystyle \underset{a}{}}\theta _a\varphi _a\right)^2.`$
In Eq. (199), we had added the final $`M^2`$ term to fix a redundancy in the description. In it’s absence, the action $`s`$ possesses a continuous translational symmetry inherited from the definition of the projection operators, $`\theta _x\theta _x+\lambda `$, $`\varphi _y\varphi _y+\lambda `$, for all $`x,y`$. This represents a redundancy in the constraints due to the fact that the total particle number of the system is given both by the sum of the row particle numbers and by the sum of the column particle numbers. The addition of the $`M^2`$ term has no effect on physical quantities, and has the benefit of lifting the unnecessary zero-mode of $`s`$. We will in fact take $`M^2\mathrm{}`$ ultimately for simplicity at an appropriate stage of the calculation.

Note the $`L`$ prefactor in front of the second term in Eq. (199), which suggests a saddle point approximation for large $`L`$ (in the first term, note also that there are $`L^2`$ components in the sum but only $`2L`$ variables). The saddle point conditions, $`\frac{s}{\theta _x}=\frac{s}{\varphi _y}=0`$, are

$`{\displaystyle \underset{y=1}{\overset{L}{}}}\mathrm{tan}(\theta _x\varphi _y)`$ $`=`$ $`imLM^2{\displaystyle \underset{a}{}}(\theta _a+\varphi _a),`$ (200)
$`{\displaystyle \underset{x=1}{\overset{L}{}}}\mathrm{tan}(\theta _x\varphi _y)`$ $`=`$ $`imLM^2{\displaystyle \underset{a}{}}(\theta _a+\varphi _a).`$ (201)
These are solved by the uniform (imaginary) solution $`\theta _x=\varphi _y=i\eta /2`$, with

$$\mathrm{tanh}\eta =m.$$
(202)
Next, expanding around this saddle, we let $`\theta _x=i\eta /2+\frac{1}{2}(\psi _{+x}+\psi _x)`$, $`\varphi _y=i\eta /2+\frac{1}{2}(\psi _{+x}\psi _x)`$. Expanding the action to cubic order in $`\psi _\pm `$, as justified by the following analysis, one finds $`s=s_0+s_2+s_3`$, with

$`s_0`$ $`=`$ $`{\displaystyle \frac{1}{2}}L^2\left(m\mathrm{ln}\left[{\displaystyle \frac{1+m}{1m}}\right]+\mathrm{ln}(1m^2)\right),`$ (203)
$`s_2`$ $`=`$ $`{\displaystyle \frac{1}{4}}(1m^2)[L{\displaystyle \underset{a}{}}(\psi _{+a}^2+\psi _a^2)`$ (204)
$`+{\displaystyle \underset{a,b}{}}((2M^21)\psi _{+a}\psi _{+b}+\psi _a\psi _b)],`$
$`s_3`$ $`=`$ $`{\displaystyle \frac{i}{12}}m(1m^2)[L{\displaystyle \underset{a}{}}(\psi _a^3+3\psi _{+a}^2\psi _a)`$ (205)
$`+3{\displaystyle \underset{a,b}{}}(\psi _a^2\psi _b+\psi _{+a}^2\psi _b2\psi _{+a}\psi _{+b}\psi _b)].`$
The constant term $`s_0`$ drops out of the observables in which we are interested. The quadratic action, $`s_2`$, governs small fluctuations of $`\psi _\pm `$, treating treating cubic and higher order terms as perturbations. The resulting Gaussian averages obey Wick’s theorem, governed by the propagators that are obtained in the usual way by inverting the quadratic form in $`s_2`$.

$`G_{ab}^+=\psi _{+a}\psi _{+b}_{s_2}`$ $`=`$ $`{\displaystyle \frac{2}{1m^2}}\left({\displaystyle \frac{\delta _{ab}}{L}}{\displaystyle \frac{1}{L^2}}\right),`$ (206)
$`G_{ab}^{}=\psi _a\psi _b_{s_2}`$ $`=`$ $`{\displaystyle \frac{2}{1m^2}}\left({\displaystyle \frac{\delta _{ab}}{L}}{\displaystyle \frac{1}{2L^2}}\right),`$ (207)
where in Eq. (206) we have taken for simplicity the $`M^2\mathrm{}`$ limit. Eqs. (206-207) imply that the fluctuations (variance) of individual $`\psi _{a\pm }`$ fields are small ($`O(1/L)`$), and moreover the correlations between fields at different $`ab`$ are even smaller ($`O(1/L^2)`$).

For this reason, the energy density $`ϵ(m)`$ is determined by the saddle-point value alone, i.e.

$$P_{\mathrm{flip}}(x,y;\{\sigma \})_{\{\sigma _{xy}\}}=F_{xy}[\theta =\varphi =i\eta /2],$$
(208)
leading directly to Eq. (152)

The “single-particle” gap $`\mathrm{\Delta }_1`$ is slightly more involved, since we need the energy up to terms of order $`1/L`$. To calculate it, we require $`E_1(m;\lambda )`$, which is obtained from the expectation value of the flippability in a state with $`\lambda `$ additional bosons on row and column 1. This is obtained, according to Eq. (154), by slightly modifying the projection operators on this row and column. This amounts to adding an additional source term to the action, $`s[\{\theta ,\varphi \}]s[\{\theta ,\varphi \}]+2i\lambda (\theta _1\varphi _1)`$, but otherwise calculating the same expectation value as in Eq. (196). Transforming to the $`\psi _\pm `$ variables as before, this additional linear term can be removed from the action by shifting $`\psi _a\psi _a2i\lambda G_{a1}^{}`$, which “completes the square” in the action using $`s_2`$ from Eq. (204). Dependence upon $`\lambda `$ thereby moves into $`F`$ and $`s_3`$:

$`E_1(m;\lambda )`$ $`=`$ $`v{\displaystyle \underset{xy}{}}F_{xy}[\theta _a={\displaystyle \frac{i}{2}}\eta +{\displaystyle \frac{1}{2}}(\psi _{+a}+\psi _a2i\lambda G_{a1}^{}),\varphi _a={\displaystyle \frac{i}{2}}\eta +{\displaystyle \frac{1}{2}}(\psi _{+a}\psi _a+2i\lambda G_{a1}^{})]^{},`$ (209)
where the prime indicates the expectation value is as defined in Eq. (197), except that $`s_3[\psi _{+a},\psi _a]`$ is replaced by $`s_3^{}[\psi _{+a},\psi _a]=s_3[\psi _{+a},\psi _a2i\lambda G_{a1}^{}]`$. This can be evaluated by expanding $`F_{xy}`$ above in $`\psi _\pm `$ and $`\lambda `$ (the latter since $`G_{a1}^{}1`$), and further evaluating the expectation value perturbatively in $`s_3`$. A careful examination of these expansions shows that there are only two contributions to $`\mathrm{\Delta }_1`$. The first comes at $`O(\psi _\pm ^0)`$ from the pure saddle-point contribution to $`F_{xy}`$, expanded to second order in $`\lambda `$. The second comes from expanding $`F_{xy}`$ to $`O(\psi _\pm ^2)`$ at $`\lambda =0`$ and evaluating the expectation value in Eq. (209) to first order in $`s_3`$ (itself expanded to $`O(\lambda ^2\psi _\pm )`$). Adding the two gives the result in Eq. (155).