# Emergent U(1) symmetry in non-particle-conserving one-dimensional models

Zakaria Jouini, Natalia Chepiga, Loïc Herviou, Frédéric Mila

► To cite this version:

Zakaria Jouini, Natalia Chepiga, Loïc Herviou, Frédéric Mila. Emergent U(1) symmetry in non-particle-conserving one-dimensional models. Physical Review B, 2023, 108 (205145), 10.1103/phys-revb.108.205145. hal-04727006

HAL Id: hal-04727006

https://hal.science/hal-04727006v1

Submitted on 9 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Emergent U(1) symmetry in non-particle-conserving one-dimensional models

Zakaria Jouini,¹ Natalia Chepiga,² Loïc Herviou,¹ and Frédéric Mila¹

¹Institute of Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

²Kavli Institute of Nanoscience, Delft University of Technology, Lorentzweg 1, 2628 CJ Delft, The Netherlands

(Dated: December 8, 2023)

The properties of stable Luttinger liquid phases in models with a non-conserved number of particles are investigated. We study the Luttinger liquid phases in one-dimensional models of hard-core boson and spinless fermion chains where particles can be created and annihilated three by three on adjacent sites. We provide an intuitive and systematic method based on the flow equation approach, which accounts for additional terms in the correlations generated by the $Z_3$-symmetric interactions. We find that despite the emergence of U(1) symmetry under renormalization, the observables are still affected by its breaking in the bare Hamiltonian. In particular, the standard bosonization mapping becomes insufficient to capture the full behavior of correlation functions.

# I. INTRODUCTION

The observation of density-wave order in recent experiments on one-dimensional systems of Rydberg atoms¹,² has brought back unsolved questions about the nature of the commensurate melting of period-$p$ phases³⁻⁸. For $p > 2$, a floating phase, characterized by incommensurate and algebraic correlations, separates the $Z_p$-ordered and the disordered phases⁹,⁸⁻¹⁰. For $p = 3$ and $p = 4$, the extension of the floating phase is still debated⁸,¹¹⁻²⁰. One open question is the existence of a direct and continuous transition in the chiral universality class that would occur before the floating phase develops⁴,⁸. The quantum version of this problem can be formulated in terms of hard-core bosons associated with domain walls in the commensurate structure of the density wave²¹. The resulting Hamiltonian exhibits a $Z_p$ symmetry as it contains terms that create and annihilate $p$ adjacent particles. When these perturbations to the free-fermion fixed point are irrelevant, the U(1) charge conservation is restored at the scaling limit in an extended incommensurate Luttinger liquid phase¹⁷,²², equivalent to the floating phase of the 2D classical problem⁸.

The properties of the Luttinger liquid phase are described by a bosonic conformal field theory²³,²⁴. The correspondence between the operators at the lattice scale and the bosonic fields in the continuum limit can be in principle constructed from selection rules dictated by the symmetry of the lattice Hamiltonian. When the latter exhibits U(1) symmetry, this correspondence is generally given by the standard bosonization mapping²⁴. For models that do not conserve the number of particles, the latter becomes insufficient to capture the long-distance behavior of the correlations. A similar problem was discussed recently in the case of nonsymmetric 1D models²⁵,²⁶. We provide in this paper a way to derive a bosonic representation of lattice operators in $Z_p$-symmetric 1D models. Using the flow equation approach²⁷, introduced by Wegner in 1994, a continuous unitary transformation is designed to restore the U(1) symmetry perturbatively in an effective Hamiltonian. The action of these transformations on the lattice opera-

tor generates an expansion in terms of the bosonic fields of the Luttinger liquid theory. The bosonization mappings obtained from this procedure are sufficient to capture the long-distance behavior of the correlations. We illustrate this method on two chains of spinless fermions and hard-core bosons where particles are created and annihilated three by three on adjacent sites. Our findings are assessed by density matrix renormalization group (DMRG)³⁰⁻³³ simulations.

The paper is organized as follows. In Section II, we discuss the phase diagrams of the hard-core boson and spinless fermion models. Section III provides a brief introduction to the flow equation approach. In Section IV, the flow equation approach is applied to the fermionic model to derive perturbatively a U(1)-symmetric effective Hamiltonian. A modified bosonic representation of the single-fermion operator is then derived and used to calculate correlation functions inside the Luttinger liquid phase. In Section V, the same procedure is applied to the continuum limit of the hard-core boson model, using the generator that diagonalizes the dual sine-Gordon model. The results are summarized in Section VI.

# II. MODELS AND PHASE DIAGRAMS

We consider one-dimensional models of hard-core bosons and spinless fermions where particles are created and annihilated three by three on adjacent sites. The two models share the feature of a Luttinger liquid phase that remains stable when the $Z_3$-symmetric interaction is turned on. In this phase, the low-energy properties of the system are described by the Luttinger liquid Hamiltonian

$$H_{LL} = \frac{1}{2\pi} \int dx v K[\partial_x \theta(x)]^2 + \frac{v}{K}[\partial_x \phi(x)]^2, \quad (1)$$

where $K$ is the Luttinger parameter and $v$ is the velocity. The fields $\phi$ and $\theta$ are bosonic in nature and satisfy the commutation relation $[\phi(x), \theta(y)] = i\pi \text{sign}(y-x)/2$. The correlations decay algebraically with an exponent controlled by $K$ and oscillate with an incommensurate wave vector proportional to the Fermi wave vector $k_F$.
We present in this section the arguments for the stability of the Luttinger liquid phase and discuss the nature of the transitions out of it.

## A. Hard-core bosonic model

The Hamiltonian of the hard-core boson model is given by

$$H = \sum_i -t(b_{i+1}^\dagger b_i + \text{h.c.}) - \mu n_i + \lambda(b_i^\dagger b_{i+1}^\dagger b_{i+2}^\dagger + \text{h.c.}), \quad (2)$$

where $b_i^\dagger$ and $b_i$ are respectively the creation and annihilation operators of hard-core bosons at site $i$, and $n_i = b_i^\dagger b_i$ is the density operator. The hard-core constraint amounts to a restriction of the occupation number to $n_i = 0$ and $n_i = 1$. Accordingly, the operators satisfy the commutation relation $[b_i, b_j^\dagger] = (1 - 2n_i)\delta_{i,j}$. The hard-core boson model was recently introduced as a dual description of the transition to phases with a density-wave order of period 3 in chains of Rydberg atoms$^{15,22}$. Its phase diagram is studied extensively in Ref. 17. We recall here the main results. Without loss of generality, the hopping amplitude is set to $t = 1$. At $\lambda = 0$, Eq. (2) is a free fermion Hamiltonian that can be mapped in the continuum limit to the Luttinger liquid Hamiltonian (1) with a Luttinger exponent $K = 1$ and a velocity $v = 2 \sin(k_F)$. When $\lambda \neq 0$, the U(1) symmetry is reduced to a $\mathbb{Z}_3$ symmetry. The stability of the Luttinger liquid phase follows from the scaling analysis of the Hamiltonian in the continuum limit, obtained by applying the bosonization mapping$^{24}$. Since the bosonic representation of the hard-core boson operators takes the form $b \sim e^{i\theta}$, the Hamiltonian reduces, at half filling ($\mu = 0$) and up to the most relevant term, to

$$H \sim H_{\text{LL}} + \frac{g}{\pi\alpha^2} \int dx \cos(3\theta(x)), \quad (3)$$

where $g$ is a dimensionless coupling and $\alpha$ is a real-space cutoff. A standard renormalization group (RG) analysis$^{24}$ yields the RG equations

$$\begin{aligned} \frac{dK(l)}{dl} &= \frac{9}{4}g^2(l), \\ \frac{dg(l)}{dl} &= \left(2 - \frac{9}{4K(l)}\right)g(l). \end{aligned} \quad (4)$$

When $K < 9/8$, the coupling constant decays exponentially under renormalization such that the Luttinger liquid Hamiltonian is recovered with effective parameters $K^*$ and $v^*$ that depend on the bare coupling constants in Eq. (3). The U(1) symmetry is thus restored at the scaling limit. It manifests itself as a symmetry under translation of the dual field $\theta$. At $K = 9/8$, the system undergoes a Kosterlitz-Thouless (KT) transition$^{28}$ into a $\mathbb{Z}_3$-ordered phase by pinning the dual field in the minima of the cosine, i.e., at $\theta_n = 2\pi n/3$. Upon varying the

chemical potential $\mu$ inside the Luttinger liquid phase, a commensurate-incommensurate transition into a disordered phase occurs. This transition is in the Pokrovsky-Talapov (PT) universality class$^{29}$, characterized by a dynamical exponent $z = 2$ and an incommensurate correlations wave vector that approaches its commensurate value with a singularity proportional to $|\mu - \mu_c|^{1/2}$. The KT and PT lines are expected to meet at a Lifshitz point that would appear before the three-state Potts point$^8$. Between the two points, the commensurate melting is direct and takes place in the chiral universality class. These predictions are confirmed numerically in Ref. 17.

## B. Fermionic model

The Hamiltonian of the fermionic model is given by

$$H = \sum_i -t(c_{i+1}^\dagger c_i + \text{h.c.}) - \mu n_i + \lambda(c_i^\dagger c_{i+1}^\dagger c_{i+2}^\dagger + \text{h.c.}), \quad (5)$$

where $c_i^\dagger$ and $c_i$ are the creation and annihilation operator of spinless fermions. It differs from the hard-core boson model (2) by a string operator when the Jordan-Wigner transformation is applied.

<figure><img src="image_2.png" /><figcaption>FIG. 1. (a)-(b) Scaling of the incommensurate wave vector as a function of the chemical potential close to the PT transition for (a) $\lambda = 0.2$ and (b) $\lambda = 2$. (c) Luttinger parameter $K$ as a function of $\lambda$ at $\mu = 0$, extracted numerically from the decay of the two-fermion correlations. This procedure is only valid when $K$ is larger than $\sqrt{3}$ (see Section IV C 2). The system is simulated with $N = 600$ sites.</figcaption></figure>

Similar to the bosonic case, the stability of the Luttinger liquid phase in Eq. (5) can be investigated using
the bosonization mapping<sup>24</sup>. In terms of the fields $\phi$ and $\theta$, the Hamiltonian reads

$$H \sim H_{LL} + \frac{g}{\pi\alpha^2} \int dx \cos(3\phi(x) - 3k_F x) \cos(3\theta(x)), \quad (6)$$

where terms that oscillate at $k_F$ are neglected as they vanish under the integration for all fillings of the band inside the Luttinger liquid phase. Eq. (6) becomes non-oscillating at $k_F = 2\pi/3$. The perturbation to the Luttinger liquid Hamiltonian has a scaling dimension $\Delta = 9(K + 1/K)/4$ and becomes relevant when $\Delta < 2$. This inequality has however no solution and the coupling constant $g$ decays exponentially under renormalization for all values of $K$. It should be noted that terms such as $\cos(6\theta)$ and $\cos(6\phi)$ are expected to arise along the RG flow and can lead to a KT transition when $K > 9/2$ or $K < 2/9$. DMRG computations of $K$ suggest however that it increases slowly enough from $K = 1$ for the generated terms to remain irrelevant even at large values of $\lambda$. For instance, the numerical results at $\mu = 0$ indicate an extension of the Luttinger liquid phase at least up to $\lambda \sim 50$ (Fig.1). Nevertheless, a KT transition at larger values of $\lambda$ cannot be excluded. As the chemical potential is varied in the Luttinger liquid phase, a PT transition line can be identified numerically using the scaling behavior of the wave vector (Fig.1).

### III. FLOW EQUATION APPROACH

The general idea behind the flow equation approach introduced by Wegner<sup>27</sup> is to apply continuous unitary transformations to the Hamiltonian in order to bring it into a more band-diagonal form. The approach consists in a renormalization scheme where states with large energy differences are first decoupled while smaller energy differences are later suppressed along the flow.

The formalism of the method is based on the parametrization of a set of unitarily equivalent Hamiltonians $H(l) = U(l)H(0)U^{\dagger}(l)$. By taking the derivative with respect to $l$, the problem is recast into the differential equation

$$\frac{dH(l)}{dl} = [\eta(l), H(l)], \quad (7)$$

where

$$\eta(l) = \frac{dU(l)}{dl} U^{\dagger}(l) \quad (8)$$

is the anti-hermitian generator of the flow. The latter can be chosen appropriately to diagonalize the bare Hamiltonian $H(0)$. The canonical choice proposed by Wegner<sup>27</sup> is given by $\eta(l) = [H_d(l), H_{od}(l)]$, where $H_d(l)$ and $H_{od}(l)$ are the diagonal and the off-diagonal parts of the flowing Hamiltonian $H(l)$. From this definition of the generator, it can be shown<sup>34</sup> that

$$\frac{d\mathrm{Tr}(H_{\mathrm{od}}^2(l))}{dl} = -2\mathrm{Tr}(\eta^\dagger(l)\eta(l)) \le 0, \quad (9)$$

which indicates that the flow gradually brings the Hamiltonian into a more band-diagonal form. A fixed point of the flow is reached when $\eta(l)$ vanishes. The diagonal and off-diagonal parts of the flowing Hamiltonian then commute and the Hamiltonian becomes block-diagonal with respect to the symmetry of the non-interacting part. Thus, the flow equation approach provides a systematic way to design a unitary transformation that recovers the U(1) symmetry in models that do not conserve the number of particles. Its drawback is the proliferation of terms along the flow. Truncations schemes are hence needed to keep the calculations tractable.

Once the diagonal Hamiltonian is obtained from the flow equation procedure, the change of basis associated with $\eta(l)$ can be applied to the operators. Given an operator $O$, its transformation along the flow is dictated by the flow equation

$$\frac{dO(l)}{dl} = [\eta(l), O(l)]. \quad (10)$$

The expectation value in the ground state of the bare Hamiltonian can be evaluated using the relation

$$\langle \psi_{\mathrm{gs}} | O | \psi_{\mathrm{gs}} \rangle = \langle \psi_{\mathrm{gs}}(\infty) | O(\infty) | \psi_{\mathrm{gs}}(\infty) \rangle, \quad (11)$$

where $|\psi_{\mathrm{gs}}(\infty)\rangle = U^{\dagger}(\infty)|\psi_{\mathrm{gs}}\rangle$ is the ground state of the diagonal Hamiltonian $H(\infty)$.

In the following sections, we apply the flow equation procedure to the models in Eqs. (2),(5). The transformed bosonic representations of hard-core bosonic and fermionic operators is derived from the generator of the flow that restores the U(1) symmetry in the Hamiltonians.

### IV. FLOW EQUATION APPROACH TO THE FERMIONIC MODEL

We derive in this section a U(1)-symmetric effective Hamiltonian that describes the low-energy properties of the fermionic Hamiltonian (5).

#### A. Flow of the Hamiltonian

We proceed by writing the Hamiltonian in the Fourier basis in order to separate it into a diagonal part $H_0$ and an interaction part $\lambda H_3$. We have

$$H = \sum_k \xi_k c_k^\dagger c_k + \frac{\lambda}{3!\sqrt{N}} \sum_{k,q} B_{k,q} (c_k^\dagger c_q^\dagger c_{-k-q}^\dagger - \mathrm{h.c.}), \quad (12)$$

where

$$\xi_k = -2 \cos(k) - \mu \\ B_{k,q} = 2i[\sin(2k+q) - \sin(2q+k) - \sin(k-q)]. \quad (13)$$

Along the flow, other interaction terms that are not initially present are generated. They are incorporated in
the following ansatz for the flowing Hamiltonian $H(l)$:

$$H(l) = H_0(l) + \lambda H_3(l) + \lambda^2 H_U(l), \quad (14)$$

where

$$H_U(l) = \frac{1}{N} \sum_{k,q,p} U_{k,q,p}(l) c_{k+p}^{\dagger} c_{q-p}^{\dagger} c_{q} c_{k}. \quad (15)$$

The normal-ordering of terms in Eq. (14) with respect to the Fermi sea of the diagonal Hamiltonian $H_0(l)$ is implicit. It is systematically carried out during the calculation to truncate the generated terms in a controlled manner. We note that three-body interaction and six-fermion terms also arise along the flow. They are neglected in the ansatz (14) due to their large scaling dimension. By taking the generator as $\eta(l) := [H_0(l), \lambda H_3(l)]$, which reads

$$\eta(l) = \frac{\lambda}{3!\sqrt{N}} \sum_{k,q} B_{k,q}(l) \alpha_{k,q}(l) (c_k^\dagger c_q^\dagger c_{-k-q}^\dagger + \text{h.c.}), \quad (16)$$

with $\alpha_{k,q} = \xi_k + \xi_q + \xi_{-k-q}$, we obtain a set of flow equations:

$$\begin{aligned} \frac{dB_{k,q}(l)}{dl} &= -\alpha_{k,q}^2(l)B_{k,q}(l), \\ \frac{dU_{k,q,p}(l)}{dl} &= -\frac{1}{2}[\alpha_{k+p,q-p}(l) + \alpha_{k,q}(l)]B_{k,q}(l)B_{k+p,q-p}(l) \\ & \qquad \times [2\Theta(k_F - |k+q|) - 1] \\ \frac{d\xi_k(l)}{dl} &= -\frac{\lambda^2}{N} \left[ \sum_{|q|<k_F} \alpha_{k,q}(l)B_{k,q}^2(l)\Theta(k_F - |k+q|) \right. \\ & \qquad \left. + \sum_{|q|>k_F} \alpha_{k,q}(l)B_{k,q}^2(l)[1 - \Theta(k_F - |k+q|)] \right] \end{aligned} \quad (17)$$

with the initial conditions $B_{k,q}(0) = B_{k,q}$, $\xi_k(0) = \xi_k$, and $U_{k,q,p}(0) = 0$. We have omitted for brevity the $l$-dependence in the bare Hamiltonian, i.e. $H \equiv H(0)$. The derivation of the flow equations is detailed in Appendix A. To the leading order in $\lambda$, it is sufficient to take the bare value $\alpha_{k,q}$ in the flow equation of $B_{k,q}(l)$. The solution shows that the three-site term decays to zero along the flow as $B_{k,q}(l) = B_{k,q}e^{-\alpha_{k,q}l}$. On the other hand, the generated terms take finite values at $l = \infty$ given by

$$\begin{aligned} U_{k,q,p}(\infty) &= -\frac{1}{2} \frac{\alpha_{k+p,q-p} + \alpha_{k,q}}{\alpha_{k,q}^2 + \alpha_{k+p,q-p}^2} B_{k,q} B_{k+p,q-p} \\ & \qquad \times [2\Theta(k_F - |k+q|) - 1] \\ \xi_k(\infty) &= \xi_k - \frac{\lambda^2}{2N} \left[ \sum_{|q|<k_F} \frac{B_{k,q}^2}{\alpha_{k,q}} \Theta(k_F - |k+q|) \right. \\ & \qquad \left. + \sum_{|q|>k_F} \frac{B_{k,q}^2}{\alpha_{k,q}} [1 - \Theta(k_F - |k+q|)] \right] \end{aligned} \quad (18)$$

By construction, the two-body interaction matrix $U_{k,q,p}$ is symmetrized with respect to the permutations of

fermion operators in Eq. (15). The divergence as $\alpha_{k,q}$ goes to zero in the renormalized quantities (18) is an artifact of the limit $l \to \infty$. In fact, the contribution of the flow to the dispersion vanishes when $\alpha_{k,q} = 0$. Finally, we note that Wegner's prescription is relaxed by taking in the definition (16) of $\eta(l)$ only the three-site term and not all the interacting part of $H(l)$. As a result, the approach becomes perturbative in $\lambda$. Terms that are generated along the flow have a second-order dependence in $\lambda$ and can be eliminated by including them in a redefinition of the generator. This induces higher-order corrections to the flow equations (17). Thus, the flow equation procedure enables to push the $\lambda$-dependence of terms that break the U(1) symmetry to higher orders, thereby restoring the symmetry perturbatively.

## B. Bosonization of the effective Hamiltonian

We investigate in this section the low-energy properties of the Hamiltonian (5) starting from the effective Hamiltonian $H_{\text{eff}} = H_0(\infty) + \lambda^2 H_U(\infty)$. After taking the long-range limit in the two-body interaction term, only particle-hole excitations within a momentum range $\Lambda$ around the Fermi points $\pm k_F$ are retained. The renormalized dispersion relation can then be linearized around the Fermi point and the fermion operators separated into right and left modes:

$$c_k = \Theta(\Lambda - |k - k_F|)c_{R,k} + \Theta(\Lambda - |k + k_F|)c_{L,k}, \quad (19)$$

where $\Theta$ is the Heaviside function. After constructing fermionic fields from these modes, i.e.,

$$c(x) = \frac{1}{\sqrt{N}} \sum_k e^{ikx} c_k := c_R(x) + c_L(x), \quad (20)$$

we obtain in the continuum limit a Hamiltonian that describes the low-lying states of $H_{\text{eff}}$. It is given by

$$H_{\text{eff}} = -i\tilde{v}_F \int dx [c_R^\dagger(x) \partial_x c_R(x) - c_L^\dagger(x) \partial_x c_L(x)] \\ + 4\lambda^2 g_2 \int dx \rho_R(x) \rho_L(x), \quad (21)$$

where $\tilde{v}_F = \partial\xi_k(\infty)/\partial k|_{k=k_F}$ is the renormalized Fermi velocity, $\rho_R(x)$ and $\rho_L(x)$ are respectively the density operators of the right and left branches, and $g_2 = U_{k_F, -k_F, 0}$ is the forward scattering matrix element. Due to the symmetry of the two-body interaction matrix, $g_4$ scattering processes that couple fermions at the same branch vanish. Finally, the factor 4 accounts for the two possible $g_2$ processes and the two backscattering processes, i.e., $g_1 = -U_{k_F, -k_F, 2k_F}$, which for spinless fermions coincide with $g_2$ processes. We note that since $g_2 > 0$, the interaction is attractive, which indicates that the three-site term in the bare Hamiltonian favors the occupation of three adjacent sites.
We now apply the bosonization mapping between the fermionic fields $c_R, c_L$ and the bosonic fields $\phi, \theta$. It is given by$^{24}$

$$c_r(x) = \frac{F_r}{\sqrt{2\pi\alpha}} e^{i\pi k_F x} e^{-i[r\phi(x) - \theta(x)]}, \quad (22)$$

where $r = 1$ for $r = R$ and $r = -1$ for $r = L$. Here, $\alpha \sim 1/\Lambda$ is a short-distance cut-off and $F_r$ are unitary operators called Klein factors. They follow the commutation relations $\{F_r, F_{r'}^\dagger\} = 2\delta_{r,r'}$ and ensure the anti-commutation of fermions from different species. Using Eq. (22), the Hamiltonian (21) can be reduced to the Luttinger liquid Hamiltonian (1) with a renormalized velocity $u$ and a Luttinger parameter, given to the second order of $\lambda$ by

$$K = 1 + \frac{4}{\pi} \sin(k_F) \sin^2\left(\frac{k_F}{2}\right) \lambda^2. \quad (23)$$

The derivation of Eq. (23) is detailed in Appendix C. Here, $k_F$ is defined by the filling of the renormalized band $\xi_k(\infty)$ and differs from its value at the non-interacting limit $\lambda = 0$. The relation between $k_F$ in Eq. (23) and the bare chemical potential can be obtained by setting $\xi_{k_F}(\infty) = 0$. This leads to the auto-coherent equation

$$\mu = -2 \cos(k_F) - \lambda^2 I(k_F, \mu), \quad (24)$$

where $I$ is the sum over the $q$ modes in Eq. (18).

<figure><img src="image_3.png" /><figcaption>FIG. 2. (a) Example of the local density profile at $\lambda = 0.1$ and $\mu = 1.2$. (b) Analytical calculation and DMRG simulations of the Luttinger parameter $K$ at $\lambda = 0.1$ as a function of the bare chemical potential $\mu$, obtained by solving Eq. (24).</figcaption></figure>

The behavior of the Luttinger parameter is compared to DMRG results, along a small $\lambda$ cut, where the perturbative calculations still hold (Fig.2). $K$ is obtained numerically by fitting the profile of the local density, which for open boundary conditions exhibits Friedel oscillations. According to conformal field theory, we have$^{35,36}$

$$\langle n_j \rangle \propto \frac{\cos(2k_F j + \beta)}{[(N/\pi) \sin(\pi j/N)]^K}, \quad (25)$$

with a phase shift $\beta$. The DMRG results for $K$ follow the form described by Eq. (23). The absence of a reflection symmetry with respect to $\mu$ is due to the breaking of particle-hole symmetry by the three-site interaction in Eq. (5). The maximal deviation between the analytical and the numerical results is of order $10^{-3}$. Its origin is twofold. First, it should be noted that by considering only scattering processes with momentum transfer $p \sim 0, 2k_F$ in the effective Hamiltonian (21), we have neglected terms that oscillate with a wave vector $2k_F$. These terms contribute to the renormalization of $K$ at the 4th order of $\lambda$, mainly at half-filling. Secondly, three-body and six-fermion interaction terms are discarded in the early stages of the calculations. Including them in the generator of the flow also provides corrections of order 4 in $\lambda$ to the two-body interaction matrix.

When the edges of the band $\xi_k(\infty)$ are crossed, the system undergoes a PT transition as the density of fermions respectively vanishes or saturates. For small $\lambda$, the shape of the transition line can be extracted from the solution of Eq. (24) at $k_F = 0, \pi$ (Fig.3). As the transition is approached from the Luttinger liquid phase, the two-body interaction matrix vanishes, i.e., $U_{k_F, -k_F, p} = 0$ for all $p$ and the Luttinger parameter tends, up to the second order in $\lambda$, towards its non-interacting value $K = 1$. This result indicates that although the particles are not conserved in the bare Hamiltonian, they still behave as free fermions in the low density limit.

<figure><img src="image_4.png" /><figcaption>FIG. 3. (a)-(b) Theoretical prediction for the Pokrovsky-Talapov transition line at (a) $k_F = 0$ and (b) $k_F = \pi$, as a function of $\lambda$.</figcaption></figure>
### C. Transformed Fermionic operator

In this section, a modified bosonic representation of the single fermion operator is derived from the flow to the effective Hamiltonian (21). Consequences of these modifications on the correlation functions of bare operators are discussed.

#### 1. Flow of single fermion operators

In order to evaluate observables in the ground state of $H_{\text{eff}}$, the operators need also to be transformed. In particular, the form of the single fermion operator under the transformation $U(\infty)$ that block-diagonalizes the Hamiltonian can be obtained by solving the flow equation

$$ \frac{dc_k(l)}{dl} = [\eta(l), c_k(l)]. \quad (26) $$

Similar to the flow of the Hamiltonian, a closed form of the solution to Eq. (26) is not tractable and truncations should be carried out. We take here the ansatz

$$ c_k(l) = c_k + \frac{\lambda}{\sqrt{N}} \sum_q \gamma_{k,q}(l) c_q^\dagger c_{-k-q}^\dagger \\ - \frac{2\lambda}{3\sqrt{N}} \sum_{q,p} \gamma_{q,p}(l) [c_q^\dagger c_p^\dagger c_{-p-q}^\dagger - c_q c_p c_{-p-q}] c_k. \quad (27) $$

It will later be argued that higher-order terms in $\lambda$ do not bring any qualitative change to the behavior of the correlation functions. The flow equation reads

$$ \frac{d\gamma_{k,q}(l)}{dl} = -\frac{1}{2} B_{k,q}(l) \alpha_{k,q}(l), \quad (28) $$

with the initial condition $\gamma_{k,q}(0) = 0$. The solution is given by $\gamma_{k,q}(\infty) = -B_{k,q}/2\alpha_{k,q}$. Since the U(1) symmetry is recovered in the rotated basis, the standard bosonization mapping (22) can be applied. We start by constructing a fermionic field $\tilde{c}(x)$ from the transformed modes $c_k(\infty)$ in the same fashion as in Eq. (20). Given that the mode $k$ decouples from the other modes in the last term of Eq. (27), the Fourier transform yields a local operator that will at most lead to a renormalization of the numerical prefactors in correlation function. Hence, we only consider

$$ \tilde{c}(x) \sim c(x) + \frac{\lambda}{N} \int dy dz \Gamma(z, z-y) c^{\dagger}(x+y) c^{\dagger}(x+z), \quad (29) $$

where

$$ \Gamma(x, y) = \frac{1}{N} \sum_{k,q} e^{-ikx} e^{-iqy} \gamma_{k,q}(\infty). \quad (30) $$

The most relevant operators in the bosonic representation of Eq. (29) are extracted by carrying out operator

product expansions (OPE) in the product of vertex operators. We have

$$ c^{\dagger}(x+y)c^{\dagger}(x+z) \sim \frac{1}{\pi\alpha} \sin(k_F(z-y))e^{-2i\theta(x)}. \quad (31) $$

By inserting this expression back into Eq. (29) and carrying out the double integration, we obtain

$$ \tilde{c}(x) \sim \frac{e^{ik_F x}}{\sqrt{2\pi\alpha}} e^{-i[\phi(x)-\theta(x)]} + \frac{e^{-ik_F x}}{\sqrt{2\pi\alpha}} e^{i[\phi(x)+\theta(x)]} \\ + \frac{4\lambda}{\pi\alpha} \sin(k_F) e^{-2i\theta(x)}. \quad (32) $$

From the structure of the generator $\eta$, the form of higher-order terms in the transformed operator can be guessed. We will only consider operators of order $\lambda^2$ as these can contribute to the second-order expansion of the correlations by combining with the zeroth-order part of Eq. (27). Given an operator $O(c, c^\dagger)$ generated at the second order in $\lambda$, it transforms under U(1) rotations as $O(e^{i\alpha c}, e^{-i\alpha c^\dagger}) = e^{i\alpha c}O(c, c^\dagger)$, where $n$ can only take the values 1, 7, -5. Since the U(1) rotations act on the bosonic fields as $\phi \to \phi$ and $\theta \to \theta+\alpha$, this indicates that operators with a scaling dimension smaller than those in Eq. (32) cannot be generated. The terms that transform as a single-fermion operator will merely add a $\lambda^2$-dependence to the prefactor of the corresponding vertex operators. Thus, the transformed fermion operator truncated to the second order in $\lambda$ takes the form

$$ \tilde{c}(x) \sim C_1 e^{ik_F x} e^{-i[\phi(x)-\theta(x)]} \\ + C_1 e^{-ik_F x} e^{i[\phi(x)+\theta(x)]} + C_2 e^{-2i\theta(x)}, \quad (33) $$

with $C_1 \sim 1 + O(\lambda^2)$ and $C_2 \sim O(\lambda)$.

#### 2. Fermionic correlation functions

Using the effective bosonic representation of the transformed fermion operator (33), we can now compute correlation functions in the Luttinger liquid phase. Consider the point-split product of bare-fermion operator

$$ F_p(x) = \lim_{\Delta \to 0} \prod_{n=0}^{p-1} c(x + n\Delta). \quad (34) $$

The correlations in the ground state of the bare Hamiltonian can be evaluated in the Luttinger liquid ground state of $H_{\text{eff}}$ using the relation

$$ \langle F_p(x)^{\dagger} F_p(y) \rangle_{\text{GS}} = \langle \tilde{F}_p^{\dagger}(x) \tilde{F}_p(y) \rangle_{\text{LL}}, \quad (35) $$

where $\tilde{F}_p(x) = U(\infty)F_p(x)U^{\dagger}(\infty)$ is obtained by replacing $c$ with the transformed field $\tilde{c}$ in Eq. (34). The correlations of $p$ fermions are deduced from the well-known result$^{24}$ for the correlation function of vertex operators in the Luttinger liquid Hamiltonian (1):

$$ \langle e^{i[n\phi(x)+m\theta(x)]} e^{-i[n\phi(y)+m\theta(y)]} \rangle_{\text{LL}} \sim \frac{1}{|x-y|^{\frac{n^2K}{2} + \frac{m^2}{2K}}}. \quad (36) $$
Accordingly, Eqs. (33) and (35) yield the one-fermion correlations

$$ \langle F_1^\dagger(x) F_1(y) \rangle_{\text{GS}} \sim 2(C_1)^2 \frac{\cos(k_F r)}{r^{\frac{1}{2K} + \frac{K}{2}}} + (C_2)^2 \frac{1}{r^{\frac{1}{2K}}}, \quad (37) $$

where $r = |x - y|$. Similarly, the two-fermion correlations can be obtained from the bosonic representation of $F_2$. We have

$$ \tilde{F}_2(x) \sim (C_1)^2 e^{2i\theta(x)} + C_1 C_2 e^{ik_F x} e^{-i[\phi(x)+\theta(x)]} + C_1 C_2 e^{-ik_F x} e^{i[\phi(x)-\theta(x)]}, \quad (38) $$

which leads to

$$ \langle F_2^\dagger(x) F_2(y) \rangle_{\text{GS}} \sim (C_1)^4 \frac{1}{r^{\frac{1}{2K}}} + 2(C_1 C_2)^2 \frac{\cos(k_F r)}{r^{\frac{1}{2K} + \frac{K}{2}}}. \quad (39) $$

<figure><img src="image_5.png" /><figcaption>FIG. 4. One fermion and two-fermion correlation functions in the ground state of the bare Hamiltonian (5) at $\mu = 1.86$ and $\lambda = 3$. The dashed line are provided as a guide to the eye. The DMRG calculations are performed on a system of size $N = 2100$. The Luttinger parameter $K$ is extracted from Friedel oscillations as discussed in section IV B and the coefficients in Eqs. (37) and (39) are obtained from a least square fit.</figcaption></figure>

We note that a single-fermion operator appears in the bosonic representation of $F_2$. It is an example of a term that is not initially present in the standard bosonization mapping but is generated by the $Z_3$-symmetric interaction along the flow. Its consequence is a crossover between two power-laws in the two-fermion correlations. At short distance, the correlations decay with the standard exponent $2/K$ of the two-fermion operator. At large distance, a single-fermion part with an exponent $\frac{K}{2} + \frac{1}{2K}$, oscillating with a wave vector $k_F$, takes over the correlations when $K < \sqrt{3}$. The crossover takes place at length scale $l \sim 1/\lambda^{4K/(3-K^2)}$. Numerically, it can be observed near the PT transition line, where the Luttinger parameter tends to $K = 1$ while $\lambda$ remains small enough for $l$ to remain smaller than the system size. The numerical results for the correlations (Fig.4), fitted with Eqs. (37), (39) by the least-square method, are in concordance with the analytical calculations. Finally, it should be noted that a two-fermion operator is also generated in the bosonic representation of the single fermion operator (see Eq. (33)). Its effect on the correlation remains however small compared to the single-fermion part.

## V. FLOW EQUATION APPROACH TO THE HARD-CORE BOSONIC MODEL

We turn now to the hard-core boson model (2). The commutation relation of the hard-core boson operator renders the calculation of the flow on the lattice difficult. It is more convenient to consider the Hamiltonian in the continuum limit, where it reduces to the dual sine-Gordon model:

$$ H = \frac{v_F}{2\pi} \int dx [\partial_x \theta(x)]^2 + [\partial_x \phi(x)]^2 + \frac{g}{\pi\alpha^2} \int dx \cos(\beta\theta(x)), \quad (40) $$

with $\beta = 3/\sqrt{K}$. Eq. (40) is obtained from the bosonization mapping (3) by absorbing the Luttinger parameter $K$ into a redefinition of the bosonic fields, i.e., $\phi \to \phi\sqrt{K}$ and $\theta \to \theta/\sqrt{K}$. The duality transformation $\phi \leftrightarrow \theta$ and $K \to 1/K$ recovers the sine-Gordon Hamiltonian, which is studied with the flow equation approach in Refs. 37 and 38. We review here the main result of this work and establish the duality correspondence with Eq. (40) to obtain a low-energy effective Hamiltonian.

### A. Flow of the sine-Gordon Hamiltonian

We introduce the vertex operators $V_r(\beta, x) =: e^{i\beta[r\phi(x)-\theta(x)]}$, where $r = R, L$ denotes the left and right species. The relation between our notations and those used in the reference can be found in Appendix B. The columns refer to the normal ordering with respect to the ground state of the non-interacting part in Eq. (40). After performing the duality transformation, the inter-
action part of the Hamiltonian (40) reads

$$H_3 = \frac{g}{2\pi\alpha^2} \left(\frac{2\pi\alpha}{L}\right)^{\frac{\beta^2}{4}} \int dx [V_R(\tilde{\beta}/2, x)V_L(-\tilde{\beta}/2, x) + \text{h.c.}], \quad (41)$$

with L the total length of the chain and $\tilde{\beta}(K) = \beta(1/K)$. Combined with the non-interacting part, Eq. (41) is the starting point of calculations of the flow equations carried out in Ref. 37 and 38. The Hamiltonian is diagonalized by a generator $\eta(l) = \eta^{(1)}(l) + \eta^{(2)}(l)$, where

$$\eta^{(1)}(l) = -2iv_F \int dx dy \frac{\partial u(y, l)}{\partial y} \times [V_R(\tilde{\beta}/2, x)V_L(-\tilde{\beta}/2, x - y) + \text{h.c.}], \quad (42)$$

and $u(x, l)$ is obtained from the Fourier transform of $u(k, l) = \frac{g(l)}{4\pi^2\alpha^2} \left(\frac{2\pi\alpha}{L}\right)^{\frac{\beta^2}{4}} e^{-4v_F^2 k^2 l}$. Here, $g(l)$ is a running coupling that flows to zero in the weak-coupling regime, i.e., inside the Luttinger liquid phase. It is initially given by the bare coupling constant in Eq. (40). $\eta^{(2)}(l)$ generates the flow of the parameter $\beta$. Its expression can be found in Ref. 38. At the end of the flow, an effective Hamiltonian $H(\infty) = H_0 + H_d(\infty)$ is obtained, where $H_0$ denotes the non-interacting part in the bare Hamiltonian and

$$H_d(\infty) = \sum_{k>0} \omega_k(\infty) \left[ \tilde{P}_R(-k) \tilde{P}_R^\dagger(-k) + \tilde{P}_R^\dagger(k) \tilde{P}_R(k) + \tilde{P}_L(k) \tilde{P}_L^\dagger(k) + \tilde{P}_L^\dagger(-k) \tilde{P}_L(-k) \right], \quad (43)$$

with $\omega_k(\infty) = -v_F g^2 \frac{\cos(\pi\tilde{\beta}^2/4)}{2\Gamma^2(\tilde{\beta}^2/4)} k |\alpha k|^{(\tilde{\beta}^2-8)/2}$. $\tilde{P}_R(k)$ and $\tilde{P}_L(k)$ are soliton and antisoliton creation and annihilation operator defined as the Fourier transform of vertex operators (see Appendix B). The effective Hamiltonian of the dual sine-Gordon Hamiltonian (40) is deduced from the dual transformation of Eq. (43). The latter acts on the soliton and antisoliton operators as $\tilde{P}_R(k) \rightarrow P_R^\dagger(-k)$ and $\tilde{P}_L(k) \rightarrow P_L(k)$, where $P_r(k)$ is obtained from $\tilde{P}_r(k)$ by replacing $\tilde{\beta}$ with $\beta$. Moreover, the flow equations for $g$ and $K$ in Eq. (40) can be deduced (see Appendix A):

$$\frac{dK(l_{RG})}{dl_{RG}} = \frac{9}{4\Gamma(\frac{9}{4K(l_{RG})} - 1)} g^2(l_{RG}), \quad (44)$$

$$\frac{dg(l_{RG})}{dl_{RG}} = \left(2 - \frac{9}{4K(l_{RG})}\right) g(l_{RG}),$$

where $l_{RG}$ is the parameter of the RG flow. It is related to the parameter of the flow equation approach $l$ by $l_{RG} = \frac{1}{2} \ln(\frac{32l}{\alpha^2})$. The RG equations, derived in Eq. (4), are recovered by an expansion of the Gamma function around the critical value $K_c = 9/8$. Finally, we note that since the soliton and antisoliton operators transform under U(1) rotations as $P_r(k) \rightarrow e^{i\alpha} P_r(k)$, the U(1) symmetry is restored in the effective Hamiltonian.

## B. Transformed hard-core boson operator

From the Jordan-Wigner transformation and the bosonization mapping (22), a bosonic representation of the hard-core boson can be derived$^{24}$. It is given by

$$b(x) = \frac{e^{i\theta(x)}}{\sqrt{2\pi\alpha}} [1 + \cos(2\phi(x) - 2k_F x)]. \quad (45)$$

We calculate in this section the transformation of Eq. (45) along the flow that diagonalizes the dual sine-Gordon model.

### 1. Flow of the hard-core boson operator

The flow of the hard-core boson operator is evaluated using the generator of the dual sine-Gordon Hamiltonian. The latter is given by the dual of Eq. (42), i.e.,

$$\eta_{\text{dual}}^{(1)}(l) = -2iv_F \int dx dy \frac{\partial u(y, l)}{\partial y} \times [V_R(\beta/2, x)V_L(\beta/2, x - y) + \text{h.c.}]. \quad (46)$$

We take the following ansatz for the flowing hard-core boson operator, truncated to the most relevant terms:

$$b(x, l) = [C_1(l) + C_2(l) \cos(2\phi(x) - 2k_F x)]e^{i\theta(x)} + [C_3(l) + C_4(l) \cos(2\phi(x) - 2k_F x)]e^{-2i\theta(x)}, \quad (47)$$

with $C_1(0) = C_2(0) = 1/\sqrt{2\pi\alpha}$, and $C_3(0) = C_4(0) = 0$. Terms with larger scaling dimensions can be neglected in describing the behavior of the correlation functions. Since the bosonic fields in Eq. (47) are those of the original Hamiltonian (3), the flow of the rescaled operator need to be considered. We calculate here the flow of $e^{i\theta(x)}$, from which the flow equation for $C_2(l)$ can be deduced. We have

$$e^{i\theta(x)/\sqrt{K}} = \left(\frac{2\pi\alpha}{L}\right)^{\beta^2/36} V_R(-\beta/6, x) V_L(-\beta/6, x). \quad (48)$$

The details of the calculation can be found in Appendix D. We summarize here the main steps of the derivation. First, the commutator of Eq. (48) with the hermitian conjugate part of the generator (46) leads to less relevant terms that are truncated in the ansatz (47). The most relevant operator in the remaining part of $[\eta_{\text{dual}}^{(1)}(l), e^{i\theta(x)/\sqrt{K}}]$ is extracted by an OPE of the vertex operators. The parameters of the vertex operators in Eqs. (46) and (48) combine to produce the operator $e^{-2i\theta(x)}$ in the ansatz (47). Namely,

$$e^{-2i\theta(x)/\sqrt{K}} = \left(\frac{2\pi\alpha}{L}\right)^{\beta^2/9} V_R(\beta/3, x) V_L(\beta/3, x). \quad (49)$$
The flow equation for $C_3(l)$ is then given by

$$ \frac{dC_3(l)}{dl} = -C_1(l) \frac{4v_F g(l)}{\Gamma(\beta^2/12)^2} \frac{2\pi}{L} \sum_{k>0} k |\alpha k|^{\beta^2/6-2} e^{-4v_F^2 k^2 l}. \quad (50) $$

To the leading order in the bare coupling constant $g$, we can replace $C_1(l)$ by its bare value. Moreover, the running coupling constant can be replaced by its approximate solution in the weak-coupling regime$^{39}$. From the RG equations (4), we have

$$ g(l_{\text{RG}}) \sim g e^{(2-\beta^2/4)l_{\text{RG}}} = g \left(\frac{32l}{\alpha^2}\right)^{1-\beta^2/8} . \quad (51) $$

Finally, the solution of Eq. (50) at the end of the flow reads

$$ C_3(\infty) = -\frac{4v_F g}{\sqrt{2\pi\alpha}} D_\beta \int_0^\infty dk \, k^{\beta^2/6-1} f(2v_F k), \quad (52) $$

with $D_\beta = \frac{(32)^{1-\beta^2/8}}{\Gamma(\beta^2/12)}$ and $f(k) = k^{\beta^2/4-4}\Gamma(2-\beta^2/8, k^2)$.

## 2. Bosonic correlation functions

We consider the correlation functions of the point-split product of p hard-core bosonic operators

$$ B_p(x) = \lim_{\Delta \to 0} \prod_{n=0}^{p-1} b(x + n\Delta). \quad (53) $$

As in section IV C, the correlation functions in the ground state of the bare Hamiltonian (2) can be evaluated from the bosonic representation of the transformed operator $\tilde{B}_p(x) = U(\infty)B_p(x)U^\dagger(\infty)$. It is derived from $B_p(x)$ by substituting the hard-core boson operators with their transformed counterpart $b(x, \infty)$. The one-boson and two-boson operators are given by

$$ \begin{aligned} \tilde{B}_1(x) &\sim C_1 e^{i\theta(x)} + C_3 e^{-2i\theta(x)}, \\ \tilde{B}_2(x) &\sim (C_1)^2 e^{2i\theta(x)} + 2C_1 C_3 e^{-i\theta(x)}, \end{aligned} \quad (54) $$

where the coefficients are taken at the end of the flow, i.e. at $l = \infty$. The oscillating terms in Eq. (47) are neglected since they have larger scaling dimensions. Using the result (36) for the correlation functions in the Luttinger liquid Hamiltonian, we obtain

$$ \langle B_1^\dagger(x) B_1(y) \rangle_{\text{GS}} \sim (C_1)^2 \frac{1}{r^{\frac{1}{2K}}} + (C_3)^2 \frac{1}{r^{\frac{2}{K}}}, \quad (55) $$

$$ \langle B_2^\dagger(x) B_2(y) \rangle_{\text{GS}} \sim (C_1)^4 \frac{1}{r^{\frac{2}{K}}} + 4(C_1 C_3)^2 \frac{1}{r^{\frac{1}{2K}}}, \quad (56) $$

with $r = |x - y|$. For $p = 3$, a combination of $e^{i\theta(x)}$ in the product (53) can lead to a vanishing exponent. In this case, a higher-order expansion in the point-splitting parameter $\Delta$ need to be carried out. Moreover, the Klein factors become necessary as their anti-commutation prevents artificial cancellations. They are reintroduced in the ansatz by replacing $\cos(2\phi)$ with $F_R e^{2i\phi} + F_L e^{-2i\phi}$. The details of the derivation are presented in Appendix D. Finally, the correlation function takes the form

$$ \langle B_3^\dagger(x) B_3(y) \rangle_{\text{GS}} \sim D_0 \frac{1}{r^{\frac{9}{2K}}} + D_1 \frac{1}{r^2} + D_2 \frac{\cos(2k_F r)}{r^{2K}} \\ + D_3 \frac{\cos(2k_F r)}{r^{2K+2}} + D_4 \frac{1}{r^4}, \quad (57) $$

where $D_1, D_2, D_3$ and $D_4$ depend on the coefficient of the transformed operator $b(x, \infty)$.

<figure><img src="image_6.png" /><figcaption>FIG. 5. (a) One-boson and two-boson correlation functions at $\mu = 1$, $\lambda = 0.06$ and (b) Connected three-boson correlation functions at $\mu = 1.9$, $\lambda = 0.2$, in the ground state of the bare Hamiltonian. The dashed line are provided as a guide to the eye. The DMRG calculations are performed on a system of size $N = 1200$. The Luttinger parameter $K$ and the Fermi wave vector $k_F$ are extracted from Friedel oscillations as discussed in section IV B and the coefficients in the correlations are obtained from a least square fit.</figcaption></figure>

As shown in Fig. 5, the correlations decay algebraically at short distance with the exponent $p^2/2K$, associated
with the operator $e^{ip\theta(x)}$ in the standard bosonic representation of the product of $p$ hard-core bosons. At large distances, the two-boson and three-boson correlations undergo a crossover to power-laws with smaller exponents, induced by the terms generated along the flow. In particular, the two-boson correlations exhibit a crossover to the one-boson correlations. Similarly, the three-boson correlations acquire, among other terms, an oscillating part with a wave vector $2k_F$ that decays with an exponent $2K$. We also note that the presence of gradients of bosonic fields in $\tilde{B}_3(x)$ leads to a non-vanishing expectation of the three-boson operator inside the Luttinger liquid phase. Since these terms are generated to the first order of the coupling $g \propto \lambda$, the expectation decays linearly with the interaction strength and vanish at the non-interacting point $\lambda = 0$. This result is confirmed by the numerical simulations (Fig. 6).

<figure><img src="image_7.png" /><figcaption>FIG. 6. DMRG calculations of the expectation value of the three-boson operator $B_3$, in the ground state of the bare Hamiltonian (2), as a function of the three-site interaction coupling $\lambda$ and at $\mu = 1$.</figcaption></figure>

## VI. CONCLUSION

Using the flow equation approach, we have provided an analytical derivation of modified bosonic representations of hard-core boson and spinless fermion operators that take account of the $Z_3$ symmetry of the Hamiltonian. As a result, the correlation functions of $p$-particles are dominated in the long distances by the power-law decay of operators that are not initially present in the bosonic representation of the operators. These calculations can be straightforwardly generalized to $Z_n$-symmetric models. Since the generator of the flow exhibits the symmetry of the bare Hamiltonian by construction, terms that transform covariantly under $Z_n$-rotations will be generated in the bosonic representation of single-particle operators. For instance, in the hard-core boson model with creation/annihilation operators on four adjacent sites, which is of interest in the open problem of commensurate melting of density-waves, the bosonic representation of a product of $p$ bosons is expected to contain terms that are associated with $p-4m$ particles, where $m \in Z$.

Some of these terms have smaller scaling dimension than the bare operator and will dominate the correlations at long distances.

More generally, this paper demonstrates that even if an emergent $U(1)$ symmetry is present, the long distance correlation functions can have a behaviour that differs from the naive expectation and that the flow equation approach is a very useful tool to identify the correct bosonic representation of lattice operators. It would be interesting to investigate the extension of this method to models with an emergent non-abelian symmetry.

## ACKNOWLEDGMENTS

The work has been supported by the Swiss National Science Foundation (FM, ZJ) grant 212082 and by the Delft Technology Fellowship (NC). Numerical simulations have been performed on the Dutch national e-infrastructure with the support of the SURF Cooperative and the support of the Scientific IT and the Application Support Center of EPFL. We thank Melvyn Nabavi for proof reading the manuscript.
# APPENDIX A: DERIVATION OF THE FLOW EQUATIONS

## 1. Flow equations for the fermionic Hamiltonian

An efficient way of computing the flow equations is to use Wick's theorem for operator products to collect the contributions to the different generated terms in $[\eta, H_3]$. For convenience, we write the interaction term as

$$H_3 = \frac{1}{3!\sqrt{N}} \sum_{\vec{k}} \tilde{B}_{k_1,k_2,k_3} (c_{k_1}^{\dagger} c_{k_2}^{\dagger} c_{k_3}^{\dagger} + c_{k_1} c_{k_2} c_{k_3}), \qquad (A1)$$

where $\vec{k} = (k_1, k_2, k_3)$ and

$$\tilde{B}_{k_1,k_2,k_3} = \delta_{k_1+k_2+k_3,0} \sum_{\sigma \in S_3} \epsilon(\sigma) e^{-ik_{\sigma(2)}} e^{-2ik_{\sigma(3)}}, \qquad (A2)$$

is purely imaginary. $S_3$ is the symmetric group of order 3, and $\epsilon(\sigma)$ is the sign of the permutation. The contribution from $[\eta, H_3]$ to the flow of $U_{k,q,p}$ and $\xi_k$ stems from the following term:

$$[\eta, H_3] \rightarrow \frac{1}{(3!)^2 N} \sum_{\vec{k}} \tilde{B}_{k_1,k_2,k_3} \tilde{B}_{q_1,q_2,q_3} \tilde{\alpha}_{k_1,k_2,k_3} \left( [c_{k_1}^{\dagger} c_{k_2}^{\dagger} c_{k_3}^{\dagger}, c_{q_1} c_{q_2} c_{q_3}] + (\vec{k} \leftrightarrow \vec{q}) \right), \qquad (A3)$$

where $\tilde{\alpha}_{k_1,k_2,k_3} = \sum_{i=1}^3 \xi_{k_i}$. To proceed with the calculation, it is useful to apply Wick's theorem in order to normal order the generated terms with respect to the Fermi sea in the ground state of the free Hamiltonian $H_0$. The contractions of two-fermion operators are given by

$$\overline{c_k^\dagger c_q} = \Theta(k_F - |k|) \delta_{k,q}, \quad \overline{c_k c_q^\dagger} = [1 - \Theta(k_F - |k|)] \delta_{k,q}. \qquad (A4)$$

To compute the flow equation of the two-body interactions, we need to collect all the normal-ordered terms that result from single contractions in Eq. (A3). These are given by

$$\begin{aligned} & : c_{k_1}^{\dagger} c_{k_2}^{\dagger} c_{k_3}^{\dagger} :: c_{q_1} c_{q_2} c_{q_3} : \rightarrow \Theta(k_F - |k_1|) \delta_{k_1, q_1} : c_{k_2}^{\dagger} c_{k_3}^{\dagger} c_{q_2} c_{q_3} : + \text{other single contractions} \\ & = \left(\frac{1}{2!}\right)^2 \sum_{\sigma, \sigma' \in S_3} \epsilon(\sigma) \epsilon(\sigma') \Theta(k_F - |k_{\sigma(1)}|) \delta_{k_{\sigma(1)}, q_{\sigma'(1)}} : c_{k_{\sigma(2)}}^{\dagger} c_{k_{\sigma(3)}}^{\dagger} c_{q_{\sigma'(2)}} c_{q_{\sigma'(3)}} : \end{aligned} \qquad (A5)$$

where the columns denote the normal-ordering with respect to the Fermi sea and the factor $(1/2!)^2$ avoids over-counting contractions. Similarly,

$$: c_{q_1} c_{q_2} c_{q_3} :: c_{k_1}^{\dagger} c_{k_2}^{\dagger} c_{k_3}^{\dagger} : \rightarrow \left(\frac{1}{2!}\right)^2 \sum_{\sigma, \sigma' \in S_3} \epsilon(\sigma) \epsilon(\sigma') [1 - \Theta(k_F - |k_{\sigma(1)}|)] \delta_{k_{\sigma(1)}, q_{\sigma'(1)}} : c_{k_{\sigma(2)}}^{\dagger} c_{k_{\sigma(3)}}^{\dagger} c_{q_{\sigma'(2)}} c_{q_{\sigma'(3)}} : . \qquad (A6)$$

Therefore, the contribution to the two-body interaction term reads

$$\begin{aligned} [\eta, \lambda H_3] & \rightarrow \left(\frac{1}{2!}\right)^2 \frac{\lambda^2}{(3!)^2 N} \sum_{\vec{k}, \vec{q}} \tilde{B}_{k_1, k_2, k_3} \tilde{B}_{q_1, q_2, q_3} (\tilde{\alpha}_{k_1, k_2, k_3} + \tilde{\alpha}_{q_1, q_2, q_3}) \\ & \times \sum_{\sigma, \sigma' \in S_3} \epsilon(\sigma) \epsilon(\sigma') [2\Theta(k_F - |k_{\sigma(1)}|) - 1] \delta_{k_{\sigma(1)}, q_{\sigma'(1)}} : c_{k_{\sigma(2)}}^{\dagger} c_{k_{\sigma(3)}}^{\dagger} c_{q_{\sigma'(2)}} c_{q_{\sigma'(3)}} : \\ & = \left(\frac{1}{2!}\right)^2 \frac{1}{(3!)^2 N} \sum_{\vec{k}, \vec{q}} \tilde{B}_{k_{\sigma'-1}(1), k_{\sigma'-1}(2), k_{\sigma'-1}(3)} \tilde{B}_{q_{\sigma'-1}(1), q_{\sigma'-1}(2), q_{\sigma'-1}(3)} (\tilde{\alpha}_{k_{\sigma'-1}(1), k_{\sigma'-1}(2), k_{\sigma'-1}(3)} + \tilde{\alpha}_{q_{\sigma'-1}(1), q_{\sigma'-1}(2), q_{\sigma'-1}(3)}) \\ & \times \sum_{\sigma, \sigma' \in S_3} \epsilon(\sigma) \epsilon(\sigma') [2\Theta(k_F - |k_1|) - 1] \delta_{k_1, q_1} : c_{k_2}^{\dagger} c_{k_3}^{\dagger} c_{q_2} c_{q_3} : \\ & = \left(\frac{3!}{2!}\right)^2 \frac{1}{(3!)^2 N} \sum_{\vec{k}, \vec{q}} \tilde{B}_{k_1, k_2, k_3} \tilde{B}_{q_1, q_2, q_3} (\tilde{\alpha}_{k_1, k_2, k_3} + \tilde{\alpha}_{k_1, q_2, q_3}) [2\Theta(k_F - |k_1|) - 1] : c_{k_2}^{\dagger} c_{k_3}^{\dagger} c_{q_2} c_{q_3} : . \end{aligned} \qquad (A7)$$
In the second line, the change of indices $k_{\sigma'(i)} \to k_i$ for $i = 1, 2, 3$ is made in order to shift the permutations dependence to the indices of the interaction matrix. Finally, the last line is obtained from the anti-symmetry of the interaction matrix under permutation, i.e., $B_{k_{\sigma(1)},k_{\sigma(2)},k_{\sigma(3)}} = \epsilon(\sigma)B_{k_1,k_2,k_3}$. A similar result can be obtained for terms generated from double contractions, which contribute to the dispersion $\xi_k$. After writing Eq. (A7) in terms of $B_{k,q} = \tilde{B}_{k,q,-k-q}$ and $\alpha_{k,q} = \tilde{\alpha}_{k,q,-k-q}$, the contribution to the two-body term reads

$$[\eta, \lambda H_3] \to \frac{\lambda^2}{4N} \sum_{k,q,p} B_{k+p,q-p} B_{q,k} (\alpha_{k+p,q-p} + \alpha_{k,q}) [2\Theta(k_F - |k+q|) - 1] : c_{k+p}^{\dagger} c_{q-p}^{\dagger} c_q c_k : . \quad (A8)$$

Thus, the flow equation of the two-body interaction is given by

$$\frac{dU_{k,q,p}(l)}{dl} = -\frac{1}{4} B_{k+p,q-p}(l) B_{k,q}(l) (\alpha_{k+p,q-p}(l) + \alpha_{k,q}(l)) [2\Theta(k_F - |k+q|) - 1]. \quad (A9)$$

## 2. Flow equations for the hard-core bosonic Hamiltonian

The flow equations of the sine-Gordon model are derived in Ref. 38. We have

$$\begin{aligned} \frac{d\tilde{\beta}^2(l_{\text{RG}})}{dl_{\text{RG}}} &= -\frac{\tilde{\beta}^4(l_{\text{RG}})}{4\Gamma(\frac{\tilde{\beta}^2(l_{\text{RG}})}{4} - 1)} g^2(l_{\text{RG}}), \\ \frac{dg(l_{\text{RG}})}{dl_{\text{RG}}} &= \left(2 - \frac{\tilde{\beta}^2(l_{\text{RG}})}{4}\right) g(l_{\text{RG}}). \end{aligned} \quad (A10)$$

Eqs. (44) are deduced from the relation $\tilde{\beta} = 3\sqrt{K}$ and the duality transformation $K \to 1/K$.

# APPENDIX B: BOSONIZATION DICTIONARY

## 1. Bosonic fields

The bosonic field $\phi, \theta$ are constructed from the density modes $\rho_r^{\dagger}(p)$, where $r = R, L$ denotes their species. We have

$$\begin{aligned} \phi(x) &= -\frac{i\pi}{L} \sum_{p \neq 0} \frac{e^{-\frac{\alpha|p|}{2} - ipx}}{p} [\rho_R^\dagger(p) + \rho_L^\dagger(p)], \\ \theta(x) &= \frac{i\pi}{L} \sum_{p \neq 0} \frac{e^{-\frac{\alpha|p|}{2} - ipx}}{p} [\rho_R^\dagger(p) - \rho_L^\dagger(p)], \end{aligned} \quad (B1)$$

where

$$[\rho_r^\dagger(p), \rho_{r'}^\dagger(-q)] = -\delta_{r,r'} \delta_{p,q} \frac{rpL}{2\pi}. \quad (B2)$$

The relation between the density modes and the normal modes of the sine-Gordon model is given by

$$\begin{aligned} \rho_R^\dagger(p) &= \sqrt{|p|} \sigma_1(p), \\ \rho_L^\dagger(p) &= \sqrt{|p|} \sigma_2(p). \end{aligned} \quad (B3)$$

The relation between the fields in Eq. (B1) and the fields $\tilde{\phi}$ and $\tilde{\theta}$ introduced in the reference is given by

$$\begin{aligned} \phi(x) &= \sqrt{\pi} \tilde{\phi}(x), \\ \theta(x) &= -\sqrt{\pi} \tilde{\theta}(x). \end{aligned} \quad (B4)$$
## 2. Vertex operators

We define the vertex operators as

$$V_r(\beta, x) := e^{i\beta[r\phi(x) - \theta(x)]} : . \qquad (B5)$$

Note that this definition coincides with the vertex operators $\tilde{V}_r(\beta, x) := e^{i\sqrt{\pi}\beta[r\tilde{\phi}(x) - \tilde{\theta}(x)]}$: introduced in the reference. In terms of the density modes, Eq. (B5) reads

$$V_r(\beta, x) := \exp \left( r\beta \frac{2\pi}{L} \sum_{p \neq 0} \frac{e^{-\alpha|p|/2 - ipx}}{p} \rho_r^\dagger(p) \right) : . \qquad (B6)$$

The normal ordering with respect to the vacuum defined by

$$\begin{aligned} \rho_R^\dagger(p < 0) |0\rangle &= 0, \\ \rho_L^\dagger(p > 0) |0\rangle &= 0, \end{aligned} \qquad (B7)$$

yields

$$V_r(\beta, x) = \left(\frac{L}{2\pi\alpha}\right)^{\beta^2/2} e^{i\beta[r\phi(x) - \theta(x)]}. \qquad (B8)$$

The bosonization mapping is then given by

$$\psi_r(x) = \frac{e^{irk_Fx}}{\sqrt{L}} V_r(-1, x). \qquad (B9)$$

### a. Operator product expansion

The operator product expansion of a product of vertex operators is given by

$$\begin{aligned} V_R(\beta, x)V_R(-\gamma, y) &\sim \left(\frac{L/2\pi}{i(y-x)+\alpha}\right)^{\beta\gamma} V_R(\beta-\gamma, x), \\ V_L(\beta, x)V_L(-\gamma, y) &\sim \left(\frac{L/2\pi}{i(x-y)+\alpha}\right)^{\beta\gamma} V_L(\beta-\gamma, x). \end{aligned} \qquad (B10)$$

### b. Exchange relations

The order of vertex operators can be exchanged using the following relations:

$$\begin{aligned} V_R(-\gamma, y)V_R(\beta, x) &\sim V_R(\beta, x)V_R(-\gamma, y) \times \frac{[i(y-x)+\alpha]^{\beta\gamma}}{[i(x-y)+\alpha]^{\beta\gamma}}, \\ V_L(-\gamma, y)V_L(\beta, x) &\sim V_L(\beta, x)V_L(-\gamma, y) \times \frac{[i(x-y)+\alpha]^{\beta\gamma}}{[i(y-x)+\alpha]^{\beta\gamma}}, \end{aligned} \qquad (B11)$$

and $[V_R(\gamma, x), V_L(\delta, y)] = 0$ for all $\gamma, \delta$.

### c. Soliton and antisoliton operators

The soliton and antisoliton operators in the effective Hamiltonian (43) are defined as the Fourier transform of the vertex operators:

$$\tilde{P}_r(k) = \left[ \frac{\Gamma(\tilde{\beta}^2/4)}{2\pi L} \left(\frac{L|k|}{2\pi}\right)^{1-\frac{\tilde{\beta}^2}{4}} \right]^{1/2} \int dx e^{-ikx} V_r(-\tilde{\beta}/2, x), \quad \text{for } r = R, L. \qquad (B12)$$
# APPENDIX C: DERIVATION OF THE LUTTINGER LIQUID HAMILTONIAN

We derive here the Luttinger parameter and the velocity inside the Luttinger liquid phase of the model (5). We consider the particle-hole excitations close to the Fermi point $\pm k_F$ in the two-body interactions of the effective Hamiltonian $H_{\text{eff}} = H_0(\infty) + \lambda^2 H_U(\infty)$. They consist in two $g_2$ processes and two $g_1$ processes. Since the fermions are spinless, these two processes are indistinguishable and the two-body interaction term reduces to

$$H_U(\infty) = \frac{4g_2}{N} \sum_p \rho_R(p) \rho_L(-p), \qquad (C1)$$

where

$$g_2 = U_{k_F, -k_F, 0} = -\frac{16 \sin^2(k_F) \sin^4(k_F/2)}{2 + 4 \cos(k_F) + 3\mu}, \qquad (C2)$$

and

$$\rho_r(p) = \sum_k c_{r,k+p}^{\dagger} c_{r,k} \qquad (C3)$$

are the Fourier components of the density operators at the right ($r = R$) and left branches ($r = L$). Using the expressions of the density operators in terms of the fields $\phi$ and $\theta$:

$$\begin{aligned} \rho_R(x) &= -\frac{1}{2\pi} [\partial_x\phi(x) - \partial_x\theta(x)], \\ \rho_L(x) &= -\frac{1}{2\pi} [\partial_x\phi(x) + \partial_x\theta(x)], \end{aligned} \qquad (C4)$$

the two-body interaction becomes

$$H_U(\infty) = \frac{4g_2}{(2\pi)^2} \int dx [\partial_x\phi(x)]^2 - [\partial_x\theta(x)]^2. \qquad (C5)$$

Eq. (C5) is then combined with the non-interacting Luttinger liquid Hamiltonian such that the Hamiltonian remains quadratic. We obtain

$$H_U(\infty) = \frac{u}{2\pi} \int dx K[\partial_x\phi(x)]^2 + \frac{1}{K}[\partial_x\theta(x)]^2, \qquad (C6)$$

where

$$\begin{aligned} K &= \left[ \frac{1 - 2g_2\lambda^2/\pi\tilde{v}_F}{1 + 2g_2\lambda^2/\pi\tilde{v}_F} \right]^{1/2}, \\ u &= \tilde{v}_F \left[ 1 - \left( \frac{2g_2\lambda^2}{\pi\tilde{v}_F} \right)^2 \right]^{1/2}, \end{aligned} \qquad (C7)$$

and $\tilde{v}_F = \partial\xi_k/\partial k|_{k=k_F}$ is the renormalized velocity. Eq. (23) is obtained from an expansion of Eq. (C7) to the second order of $\lambda$.

# APPENDIX D: SOME DETAILS OF THE CALCULATIONS

## 1. Derivation of the flow of the hard-core boson operator

We give here the detailed calculation of the flow equation for $C_3$ in Eq. (47). We have

$$[\eta_{\text{dual}}^{(1)}(l), e^{i\theta(x)/\sqrt{K}}] = -2iv_F \left(\frac{2\pi\alpha}{L}\right)^{\beta^2/36} \int dy dz \frac{\partial u(z, l)}{\partial y} [V_R(\beta/2, y)V_L(\beta/2, y-z) + \text{h.c.}, V_R(-\beta/6, x)V_L(-\beta/6, x)]. \qquad (D1)$$
Since the terms obtained from the hermitian conjugation in Eq. (D1) are less relevant terms, we only consider the commutator

$$
\begin{align}
[V_R(\beta/2, y)V_L(\beta/2, y-z), V_R(-\beta/6, x)V_L(-\beta/6, x)] &= V_R(\beta/2, y)V_R(-\beta/6, x)V_L(\beta/2, y-z)V_L(-\beta/6, x) \nonumber \\
&\quad \times \left\{ 1 - \left[ \frac{i(x-y) + \alpha}{i(y-x) + \alpha} \right]^{\beta^2/12} \left[ \frac{i(y-z-x) + \alpha}{i(x-y+z) + \alpha} \right]^{\beta^2/12} \right\} \tag{D2} \\
&\sim \left(\frac{L}{2\pi}\right)^{\beta^2/6} V_R(\beta/3, x)V_L(\beta/3, x) \nonumber \\
&\quad \times \left\{ [i(x-y) + \alpha]^{-\beta^2/12} [i(y-z-x) + \alpha]^{-\beta^2/12} - [i(y-x) + \alpha]^{-\beta^2/12} [i(x-y+z) + \alpha]^{-\beta^2/12} \right\}, \nonumber
\end{align}
$$

where an OPE of the vertex operator is carried out in last line. By inserting this expression in Eq. (D1), we obtain,
after a few steps of calculations,

$$
\begin{align}
[\eta_{\text{dual}}^{(1)}(l), e^{i\theta(x)/\sqrt{K}}] &\rightarrow -2iv_F \alpha^{\beta^2/6} \frac{g(l)}{4\pi^2\alpha^2} \frac{2\pi}{L} \sum_k (-ik)e^{-4v_F^2 k^2 l} e^{-2i\theta(x)/\sqrt{K}} \nonumber \\
&\quad \times \int dy dz e^{-ikz} \{[i(x-y) + \alpha]^{-\beta^2/12} [i(y-z-x) + \alpha]^{-\beta^2/12} - [i(y-x) + \alpha]^{-\beta^2/12} [i(x-y+z) + \alpha]^{-\beta^2/12}\} \nonumber \\
&= -\frac{v_F g(l)}{\pi^2 L} \sum_k k e^{-4v_F^2 k^2 l} f_\beta(k) e^{-2i\theta(x)/\sqrt{K}}, \tag{D3}
\end{align}
$$

with

$$
\begin{align}
f_{\beta}(k) &= \int dy dz e^{-ik\alpha z} [1 - iy]^{-\beta^2/12} [1 + i(y-z)]^{-\beta^2/12} \tag{D4} \\
&= |\alpha k|^{\beta^2/6-2} \frac{4\pi^2}{\Gamma(\beta^2/12)^2} \Theta(k). \nonumber
\end{align}
$$

Therefore, the flow equation for $C_3(l)$ is given by

$$
\frac{dC_3(l)}{dl} = -C_1(l) \frac{4v_F g(l)}{\Gamma(\beta^2/12)^2} \frac{2\pi}{L} \sum_{k>0} k |\alpha k|^{\beta^2/6-2} e^{-4v_F^2 k^2 l}. \quad (D5)
$$

2. Derivation of $\tilde{B}_3(x)$

We give here the detailed calculation of the three-boson correlations in Eq. (57). We start by reordering the vertex
operators in the product $\tilde{B}_3(x)$. This leads to

$$
\tilde{B}_3(x) = \lim_{\Delta \to 0} \sum_{n,m,l=0,1} \left[ A_n + B_n \cos(2\tilde{\phi}(x)) \right] \left[ A_m + B_m (-1)^{p_n} \cos(2\tilde{\phi}(x+\Delta)) \right] \tag{D6} \\
\times \left[ A_l + B_l (-1)^{p_n+p_m} \cos(2\tilde{\phi}(x+2\Delta)) \right] e^{i[p_n\theta(x)+p_m\theta(x+\Delta)+p_l\theta(x+2\Delta)]}
$$

where $\tilde{\phi}(x) = \phi(x) - k_F x$, $p_n = 1 - 3n$, $A_0 = C_0$, $B_0 = C_1$, $A_1 = C_3$ and $B_1 = C_4$. To avoid artificial cancellations, it
is necessary to reintroduce the Klein factors in the definition of the ansatz. This amounts to the following substitution:

$$
\begin{align}
B_n \cos(2\phi(x)) &\rightarrow \frac{1}{2} B_n [F_R e^{2i\phi(x)} + F_L e^{-2i\phi(x)}], \tag{D7} \\
A_n &\rightarrow \frac{1}{2} A_n (F_R + F_L). \nonumber
\end{align}
$$

By carrying out an expansion in the splitting parameter $\Delta$, we obtain for a term at position $x + (n - 1)\Delta$ in the
product of Eq. (D6):

$$
\begin{align}
& \frac{1}{2} [F_R e^{2i\phi(x+(n-1)\Delta)} + F_L e^{-2i\phi(x+(n-1)\Delta)}] \nonumber \\
& \sim \frac{1}{2} F_R \left\{ 1 + 2i \left[ (n-1)\Delta \partial_x \phi + \frac{(n-1)^2 \Delta^2}{2} (\partial_x \phi)^2 \right] - 2(n-1)^2 \Delta^2 \partial_x^2 \phi \right\} e^{2i\phi} \tag{D8} \\
& \qquad + \frac{1}{2} F_L \left\{ 1 - 2i \left[ (n-1)\Delta \partial_x \phi + \frac{(n-1)^2 \Delta^2}{2} (\partial_x \phi)^2 \right] - 2(n-1)^2 \Delta^2 \partial_x^2 \phi \right\} e^{-2i\phi}, \nonumber
\end{align}
$$
where the position of the fields is at $x$ and is omitted for brevity. Similarly, we have

$$ e^{i[p_n\theta(x)+p_m\theta(x+\Delta)+p_l\theta(x+2\Delta)]} \sim e^{i[p_n+p_m+p_l]\theta} e^{i[(p_m+2p_l)\Delta\partial_x\theta+(\frac{1}{2}p_m+2p_l)\Delta^2\partial_x^2\theta]} \\ \sim e^{i[p_n+p_m+p_l]\theta} [1 + i[(p_m+2p_l)\Delta\partial_x\theta + (\frac{1}{2}p_m+2p_l)\Delta^2\partial_x^2\theta] - \frac{1}{2}(p_m+2p_l)^2\Delta^2(\partial_x\theta)^2]. \quad (D9) $$

After collecting all the terms that have a scaling dimension smaller than the bare operator $e^{3i\theta(x)}$, we obtain

$$ \tilde{B}_3(x) \sim (C_0)^3 e^{3i\theta(x)} + [8(C_1)^2 C_4 - (C_2)^2 C_4] \cos(2\tilde{\phi}) - 12\Delta C_1 C_2 C_4 \partial_x \theta \\ + \Delta[32C_1C_2C_3 - 8(C_1)^2C_4 + (C_2)^2C_4]\partial_x\phi \sin(2\tilde{\phi}) - 72\Delta^2(C_1)^2C_3(\partial_x\theta)^2 \\ + i\Delta^2[156(C_1)^2C_3 - 2(C_2)^2C_3 + 12C_1C_2C_4]\partial_x^2\theta + 16\Delta^2(C_2)^2C_3(\partial_x\phi)^2. \quad (D10) $$

It should be noted that for simplicity of the calculation, the normal-ordering of the operators in Eq. (D6) is not taken before carrying out the Taylor expansion. The latter modifies the prefactors of the generated terms in Eq. (D10).

1 A. Keesling, A. Omran, H. Levine, H. Bernien, H. Pichler, S. Choi, R. Samajdar, S. Schwartz, P. Silvi, S. Sachdev, P. Zoller, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Nature 568, 207 (2019).

2 H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Nature 551, 579 (2017).

3 S. Ostlund, Phys. Rev. B 24, 398 (1981).

4 D. A. Huse and M. E. Fisher, Phys. Rev. Lett. 49, 793 (1982).

5 D. A. Huse, A. M. Szpilka, and M. E. Fisher, Physica A: Statistical Mechanics and its Applications 121, 363 (1983).

6 S. F. Howes, Phys. Rev. B 27, 1762 (1983).

7 M. den Nijs, in Phase transitions and critical phenomena, edited by J. L. Lebowitz and C. Domb (Academic Press, 1988).

8 D. A. Huse and M. E. Fisher, Phys. Rev. B 29, 239 (1984).

9 F. D. M. Haldane, P. Bak, and T. Bohr, Phys. Rev. B 28, 2743 (1983).

10 H. Katsura, D. Schuricht, and M. Takahashi, Phys. Rev. B 92, 115137 (2015).

11 W. Selke and J. M. Yeomans, Zeitschrift für Physik B Condensed Matter 46, 311 (1982).

12 P. M. Duxbury, J. Yeomans, and P. D. Beale, Journal of Physics A: Mathematical and General 17, L179 (1984).

13 G. Giudici, A. Angelone, G. Magnifico, Z. Zeng, G. Giudice, T. Mendes-Santos, and M. Dalmonte, Phys. Rev. B 99, 094434 (2019).

14 N. Chepiga and F. Mila, Phys. Rev. Lett. 122, 017205 (2019).

15 S. Whitsitt, R. Samajdar, and S. Sachdev, Phys. Rev. B 98, 205118 (2018).

16 N. Chepiga and F. Mila, Nature Communications 12, 414 (2021).

17 N. Chepiga and F. Mila, Phys. Rev. Research 3, 023049 (2021).

18 M. Rader and A. M. Läuchli, “Floating phases in one-dimensional Rydberg Ising chains,” ArXiv:1908.02068, 2019, arXiv:1908.02068 [cond-mat.quant-gas].

19 S. Nyckees, J. Colbois, and F. Mila, Nuclear Physics B 965, 115365 (2021).

20 S. Nyckees and F. Mila, Phys. Rev. Res. 4, 013093 (2022).

21 H. J. Schulz, B. I. Halperin, and C. L. Henley, Phys. Rev. B 26, 3797 (1982).

22 P. Fendley, K. Sengupta, and S. Sachdev, Phys. Rev. B 69, 075106 (2004).

23 F. D. M. Haldane, Journal of Physics C: Solid State Physics 14, 2585 (1981).

24 T. Giamarchi, Quantum Physics in One Dimension (Oxford University Press, 2003).

25 W. Yang, C. Xu, S. Xu, A. Nocera, and I. Affleck, “Non-symmetric luttinger liquids in the generalized antiferromagnetic kitaev spin-1/2 chain,” (2022).

26 W. Yang, C. Xu, A. Nocera, and I. Affleck, Phys. Rev. B 106, 064425 (2022).

27 F. Wegner, Annalen der Physik 506, 77 (1994).

28 J. M. Kosterlitz and D. J. Thouless, Journal of Physics C: Solid State Physics 6, 1181 (1973).

29 V. L. Pokrovsky and A. L. Talapov, Phys. Rev. Lett. 42, 65 (1979).

30 S. R. White, Phys. Rev. Lett. 69, 2863 (1992).

31 U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).

32 S. Östlund and S. Rommer, Phys. Rev. Lett. 75, 3537 (1995).

33 U. Schollwöck, Annals of Physics 326, 96 (2011).

34 S. Kehrein, The Flow Equation Approach to Many-Particle Systems (Springer Berlin, Heidelberg, 2006).

35 S. R. White, I. Affleck, and D. J. Scalapino, Phys. Rev. B 65, 165122 (2002).

36 M. Fabrizio and A. O. Gogolin, Phys. Rev. B 51, 17827 (1995).

37 S. Kehrein, Phys. Rev. Lett. 83, 4914 (1999).

38 S. Kehrein, Nuclear Physics B 592, 512 (2001).

39 J. Sabio and S. Kehrein, New Journal of Physics 12, 055008 (2010).