Continuous easy-plane deconfined phase transition on the kagome lattice

TL;DR

This paper investigates a quantum phase transition in a kagome lattice model, providing evidence for a continuous transition with unconventional scaling, indicative of deconfined quantum criticality, using large-scale quantum Monte Carlo simulations.

Contribution

The study introduces a theory based on an easy-plane NCCP^1 gauge model and demonstrates, through extensive simulations, that the transition at 1/3 filling is continuous, supporting deconfined quantum criticality.

Findings

Evidence for a continuous phase transition at 1/3 filling.

Unconventional finite size scaling behavior observed.

Supports the deconfined quantum criticality scenario.

Abstract

We use large scale quantum Monte-Carlo simulations to study an extended Hubbard model of hardcore bosons on the kagome lattice. In the limit of strong nearest-neighbor interactions at 1/3 filling, the interplay between frustration and quantum fluctuations leads to a valence bond solid ground state. The system undergoes a quantum phase transition to a superfluid phase as the interaction strength is decreased. It is still under debate whether the transition is weakly first order or represents an unconventional continuous phase transition. We present a theory in terms of an easy-plane NCCP$^1$ gauge theory describing the phase transition at 1/3 filling. Utilizing large scale quantum Monte-Carlo simulations with parallel tempering in the canonical ensemble up to 15552 spins, we provide evidence that the phase transition is continuous at exactly 1/3 filling. A careful finite size scaling analysis reveals an unconventional scaling behavior hinting at deconfined quantum criticality.

Equations

$$H$$

$$\mathcal{L}=\sum _{\alpha =1}^2\left[|({\partial }_{\mu }-i{a}_{\mu })z_{\alpha }{|}^2+s|z_{\alpha }{|}^2+u|z_{\alpha }{|}^4\right]+v|z_1{|}^2|z_2{|}^2+\cdots ,$$

$$b_i^{†}$$

$$\sum _{i\in △}{n}_i$$

$$A(L,\delta )$$

$${\int }_1^L{C}_L(r)rdr/L^2=\frac{a(L)}{1-\eta (L)}(L^{-1-\eta (L)}-L^{-2}).$$

$${\eta }_f(L)={\eta }_0+{\eta }_1/L+{\eta }_2/L^2,$$

$$f_s\equiv \frac{a_f(L)}{1-{\eta }_f(L)}(L^{-1-{\eta }_f(L)}-L^{-2}),$$

$$f_a$$

$$f_2$$

$$f_{\infty }$$

Full text

Introduction. Understanding universal and non-universal properties of quantum phase transitions in strongly correlated systems is a key topic in modern physics Sachdev (2007). In many cases, quantum phase transitions can be described by Landau’s theory of spontaneous symmetry breaking just like classical ones. On the other hand, there appear to exist exotic quantum phase transitions beyond the Landau’s paradigm such as continuous deconfined phase transitions (DCPs) between phases with different, incompatible symmetry breakings. A well known example is the transition between a Néel state and a valence bond solid (VBS) Senthil et al. (2004a, b).

Contrary to conventional phase transitions, a deconfined phase transition exhibits fractionalized quasiparticles that couple to emergent gauge fields Senthil et al. (2004a, b). Deconfined phase transition are generically described by strongly interacting gauge theories. One example is the non-compact $\mathbb{CP}^{1}$ (NCCP1) model with a bosonic $\mathbb{CP}^{1}$ field $z_{\alpha}$ (describing spinons with ${\mathrm{SU}(2)}$ flavors $\alpha=1,2$), which couples to a non-compact $U(1)$ gauge field $a_{\mu}$. Depending on the symmetries of the field $z_{\alpha}$, the NCCP1 models are divided into $\mathrm{SU}(2)$-NCCP1 and the easy-plane-NCCP1, which describe the Néel to VBS transition in ${\mathrm{SU}(2)}$ or XY magnets Senthil et al. (2004a, b); Motrunich and Vishwanath (2004), respectively.

The concept of DCPs leads to several interesting questions: First, to which extent do these emergent gauge fields and fractionalized excitations appear at critical points in concrete model systems? Second, what is the fate of the NCCP1 model in the infrared (IR) limit? Recent progress in the understanding of dualities of gauge theories has brought new perspectives to deconfined phase transitions Wang and Senthil (2015); Metlitski and Vishwanath (2016); Mross et al. (2016); Xu and You (2015); Motrunich and Vishwanath (2004); Senthil and Fisher (2006); Seiberg et al. (2016); Karch and Tong (2016); Wang et al. (2017); Mross et al. (2017). It has been conjectured that the bosonic easy-plane NCCP1 theory is dual to a widely studied fermionic $N_{f}=2$ QED3 theory Senthil and Fisher (2006); Xu and You (2015); Seiberg et al. (2016); Karch and Tong (2016); Wang et al. (2017); Mross et al. (2017). Significant effort has been put into the investigation of the IR fate of QED3, but it remains an open issue after several decades of study Appelquist et al. (1986, 1988); Hands et al. (2004); Grover (2014); Braun et al. (2014); Hands et al. (2002, 2004); Karthik and Narayanan (2016). Studying concrete realizations of DCPs helps to deepen the understanding of this long-standing problem.

Numerical work Sandvik (2007); Melko and Kaul (2008); Lou et al. (2009); Banerjee et al. (2010); Sandvik (2010); Shao et al. (2016); Pujari et al. (2013); Harada et al. (2013); Jiang et al. (2008); Chen et al. (2013a); Nahum et al. (2015a, b); Motrunich and Vishwanath (2008); Kuklov et al. (2008); Charrier et al. (2008); Chen et al. (2009); Charrier and Alet (2010); Sreejith and Powell (2015); Kragset et al. (2006); Sen et al. (2007); D’Emidio and Kaul (2016, 2017); Qin et al. (2017) has studied both the $\mathrm{SU}(2)$ and easy-plane DCPs—most of them focused on the J-Q model  Sandvik (2007); Melko and Kaul (2008); Lou et al. (2009); Banerjee et al. (2010); Sandvik (2010); Shao et al. (2016); Pujari et al. (2013); Harada et al. (2013); Jiang et al. (2008); Chen et al. (2013a); D’Emidio and Kaul (2016, 2017) and classical loop models Sreejith and Powell (2015); Nahum et al. (2015b, a). It is still controversially discussed if $\mathrm{SU}(2)$ DCPs are continuous Kuklov et al. (2008), and an emergent $\mathrm{SO}(5)$ symmetry is observed Nahum et al. (2015b) between Néel and VBS phases. The easy-plane case on the other hand appears to be a first order transition in all previous numerical studies on various candidate model systems Kragset et al. (2006); Sen et al. (2007); D’Emidio and Kaul (2016, 2017); Kuklov et al. (2006). The question arises whether the easy-plane-NCCP1 is intrinsically first order or if it is specific to the models that have been studied so far.

In this paper, we provide numerical evidence for the existence of a continuous easy-plane-DCP using large scale quantum Monte-Carlo simulations, our results are in agreement with a parallel work Qin et al. (2017). Specifically, we study an extended Hubbard model of hardcore bosons on the kagome lattice,

[TABLE]

at $1/3$ filling with $t,V>0$. The system is known to form a VBS ground state in the limit $V\gg t$ and a superfluid for $V\ll t$; where both phases are separated by a quantum phase transition Isakov et al. (2006); Damle and Senthil (2006); Zhang and Eggert (2013); Carrasquilla et al. (2017). We first discuss the easy-plane NCCP1 theory Senthil et al. (2004a, b) that describes the superfluid-VBS transition. In particular, we highlight the difference between our system and other systems hosting DCPs (e.g. the J-Q model). By using large scale quantum Monte Carlo methods with parallel tempering (QMC-PT) in the canonical ensemble, we find that the phase transition between VBS to the superfluid is anomalously continuous at exactly $1/3$ filling. Several hallmarks of DCP are found: (i) At the critical point, the superfluid density decays slower than at regular continuous phase transitions. Comparing with different scenarios Shao et al. (2016); Pujari et al. (2013); Nahum et al. (2015a), we adopt logarithmic corrections to fit this drift. (ii) A direct analysis of two point correlations reveals that the anomalous critical exponent $\eta\approx 0.3$ is relatively large. (iii) We identify a lattice operator for a conserved charge (i.e. the spinon density) of NCCP1, and numerically show that its scaling dimension is close to two, as expected for a 2+1D conformal field theory (CFT) Francesco et al. (2012). (iv) An emergent $U(1)$ symmetry is identified at the critical point.

Effective theory and phases. Similar to the much studied Néel-VBS transition in antiferromagnets Senthil et al. (2004a, b), the superfluid-VBS transition in our system is also described by the NCCP1 theory,

[TABLE]

where the $z_{\alpha=1,2}$ are bosonic ($\mathbb{CP}^{1}$) fields (or spinon operators) carrying half the charge of the physical bosons, and they are coupled to an emergent dynamical $U(1)$ gauge field $a_{\mu}$. The mass term $s|z_{\alpha}|^{2}$ with $s\approx V/t-7$ controls the phases: (i) if $z_{\alpha}$ condenses, a superfluid phase is formed; (ii) if $z_{\alpha}$ is gapped, a VBS state forms due to the proliferation of monopoles of the gauge field $a_{\mu}$ Read and Sachdev (1990) (iii) the case of $z_{\alpha}$ being gapless corresponds to the critical point. The quartic terms with $u\approx O(t^{3})$ and $v\approx O(t^{3})$ control the putative IR fixed point to which the theory flows under renormalization. When $v=2u$, there is a $\mathrm{SU}(2)$ symmetry between $z_{1}$ and $z_{2}$, and the theory is called $\mathrm{SU}(2)$ NCCP1. Usually the $\mathrm{SU}(2)$ symmetry will be manifest as a global $\mathrm{SU}(2)$ symmetry of the spin system. On the other hand if $v<2u$, the theory flows to the easy-plane NCCP1 fixed point where the $\mathrm{SU}(2)$ symmetry is broken. Our hardcore boson model Eq. (1) naturally falls into the easy-plane NCCP1 class. The same field theory also describes the Néel-VBS transition in other related spin models (e.g. J-Q model), and the hard-core boson model we are studying can be exactly mapped to a spin-1/2 model by $b^{\dagger}\rightarrow S^{+}$, $n=b^{\dagger}b\rightarrow S^{z}+1/2$.

In our system, however, the relation between the continuous field operator and the lattice operators is very different from the usual DCP in spin models. In usual spin models (e.g. the J-Q model), one would have $S^{+}\sim z_{1}^{*}z_{2}$, $S^{z}\sim|z_{1}|^{2}-|z_{2}|^{2}$. In our case, the relations are

[TABLE]

Here $E_{i}$ represents the electric fields of the dynamical gauge field $a_{\mu}$, $\mathcal{M}_{a}$ is the monopole operator while $\theta_{i}$ is a phase factor ($=0,\pm 2\pi/3$) depending on the sublattice index. $\sum_{i\in\bigtriangleup,\bigtriangledown}n_{i}$ refers to the summation of the density of three sites in the up or down triangles of the kagome lattice. The difference originates from the different fractionalization schemes of the spin operator $\vec{S}$ into the $\mathbb{CP}^{1}$ (spinon) field $\mathbf{z}=(z_{1},z_{2})$. Usually at DCPs, the spin operator is fractionalized via the $\mathbb{CP}^{1}$ representation $\vec{S}_{i}=\mathbf{z}_{i}^{*}\vec{\sigma}\mathbf{z}_{i}^{T}$ Senthil et al. (2004a, b), and such a spinon operator is argued to capture the low energy physics. In contrast, our kagome model can be faithfully mapped onto a lattice gauge model defined on the medial honeycomb lattice Nikolić and Senthil (2005); He et al. (2015), in which spinons ($z_{1,2}$) live on honeycomb sites (i.e. center of kagome triangles) and $\mathrm{U}(1)$ gauge fields live on the honeycomb links. Then we can straightforwardly take the continuum limit of the lattice gauge model, which precisely gives the easy-plane NCCP1 theory.

The relations in Eq. (3)-(4) call for a slightly different way of extracting critical exponents. Specifically, the anomalous dimension $\eta_{\textrm{VBS}}$ of the VBS order parameter should be extracted from the density operator $n_{i}$, instead of the dimer operator in the J-Q model. The operators $s_{\bigtriangleup,\bigtriangledown}=\sum_{i\in\bigtriangleup,\bigtriangledown}n_{i}$, correspond to conserved charges of the gauge theory, $z^{*}_{\alpha}\partial_{t}z_{\alpha}$. For any $2+1$D CFT, such a conserved charge will always have scaling dimension two Francesco et al. (2012), providing an additional numerical check.

Numerical results. We use a stochastic cluster series expansion with parallel tempering Sandvik (1999); Syljuåsen and Sandvik (2002); Louis and Gros (2004); Sengupta et al. (2002) and adopt periodic boundary conditions with $L_{x}=L_{y}$. To reach the ground state, we use half million steps of thermalization before producing two million samples for measuring and consider temperatures down to $\beta V/L=25/3$ ($\beta=1/T$). We identify the diagonal order in the VBS phase using the structure factor $S(\mathbf{Q})=\langle|n(\mathbf{Q,\tau})|^{2}\rangle=\langle|\sum_{k=1}^{N}n_{k,\tau}e^{\mathbf{iQ}\cdot\mathbf{r}_{k}}|^{2}\rangle/N^{2}$ at $\mathbf{Q}=(4\pi/3,0)$ where $N$ is the number of sites. For the superfluid phase, we consider the superfluid density $\rho_{s}=\langle W^{2}\rangle/\beta t$ where $W$ is the winding number Pollock and Ceperley (1987) and also the condensate fraction $\rho_{0}=\langle\sum_{i,j}b_{i}^{\dagger}b_{j}\rangle/N^{2}$ to characterize long range off-diagonal correlations.

It turns out that a continuous phase transition only occurs at exactly $1/3$ filling where the system has particle-hole symmetry. In previous studies Isakov et al. (2006); Damle and Senthil (2006), a grand canonical ensemble was used for the QMC simulations which made it difficult to fine tune to exactly $1/3$ filling. Here, we restrict our simulations to the canonical ensemble by tuning the chemical potential to minimize the deviation from $1/3$ filling during the loop-update and then only accept samples with exactly $1/3$ filling. From Fig. 2, we find (i) the structure factor $S(\mathbf{Q})$ does not show any discontinuity for sizes up to a linear dimension of $L=72$ ($N=15552$ spins); (ii) its Binder cumulant $S(\mathbf{Q})_{b.c.}=1-\frac{\langle S(\mathbf{Q})^{2}\rangle}{3\langle S(\mathbf{Q})\rangle^{2}}$ is always positive and crosses at approximately same point $t_{c}/V\approx 0.1303$; (iii) at variance from Ref.Isakov et al. (2006), at larger size $L=72$ near the critical point, we do not find any double peak structure in the probability distribution of kinetic energy. Since the parameter $t/V=0.1283$ in Ref. Isakov et al. (2006) is actually far from $t_{c}/V$, it reflects the weakly first order phase transition at the upper/lower boundary of the lobe, but not at the tip. These three findings strongly support a continuous phase transition up to system size $L=72$.

Next we perform finite size scaling (FSS) for different variables to extract the critical behavior. For a continuous phase transition, the scaling function takes the form:

[TABLE]

where $\nu$ and $\kappa$ are related to the universality class of the phase transition, and $\delta=t/V-t_{c}/V$. Because the form of the scaling function $f$ is not known, we choose the method of Kawashima and Ito proposed to do the data collapse Kawashima and Ito (1993); Houdayer and Hartmann (2004).

An anomalous behavior of FSS of physical quantities has been observed in all previous numerical works on DCP (see examples Sandvik (2007); Melko and Kaul (2008); Banerjee et al. (2010); Sandvik (2010); Shao et al. (2016); Pujari et al. (2013); Nahum et al. (2015a, b); Motrunich and Vishwanath (2008)). It has been suggested that these anomalous scaling behaviors arise due to finite size effects of dangerously irrelevant operators Nahum et al. (2015a). For example, the superfluid density $\rho_{s}$ shows a drift Shao et al. (2016); Pujari et al. (2013); Nahum et al. (2015a) compared to the scaling of the conventional phase transition, $\rho_{s}(L,\delta)=L^{-1}f(\delta L^{1/\nu})$. To resolve the drift, two schemes have been proposed: (i) logarithmic corrections (LCs) $\rho_{s}(L,\delta)=L^{-1}\log(L/L_{0})f(\delta L^{1/\nu})$ Sandvik (2010); Nahum et al. (2015a) and (ii) two-length scales $\rho_{s}(L,\delta)=L^{-\nu/\nu^{\prime}}f(\delta L^{1/\nu})$ Shao et al. (2016). In our work we use the LCs and find a good data collapse with $1/\nu=2.37(0.04)$, as shown in Fig. 3a. Using two-length scales also gives a reasonably good collapse sup .

Scaling violations are also observed in the diagonal structure factor $S(\mathbf{Q})$ and condensate fraction $\rho_{0}$, whose FSS has previously been used to extract the anomalous dimension $\eta_{\text{VBS}}$ and $\eta_{\text{SF}}$. The anomalous dimensions $\eta_{\mathrm{VBS}}=0.015(0.014)$ and $\eta_{\mathrm{SF}}=0.200(0.006)$ extracted from $S(\mathbf{Q})$ and $\rho_{0}$, respectively, strongly deviate from each other (shown in Fig. 3b and Fig. 3c). This is not expected for the easy-plane NCCP1 theory as both anomalous exponents are the same due to self-duality Motrunich and Vishwanath (2004). Previous studies Block et al. (2013); Nahum et al. (2015a) find a large drift of the critical exponent and anomalous dimension due to the $1/L^{2}$-correction foo . Therefore, a simple data collapse does not give good results, but by using the two-point correlator a size-dependent anomalous dimension can be extracted, which shows a systematic convergence to the thermodynamic limit as discussed in the following and in the Supplemental Material sup .

We find two different scaling behaviors as we approach the continuous quantum phase transition from the two neighboring phases: In the disordered phase $C_{s}(\delta,r)=ar^{-1-\eta}\exp(-r/\xi(\delta))$ with the correlation length $\xi(\delta)\propto\delta^{-\nu}$, while approaching the critical point from the ordered phase, $C_{l}(\delta,r)=ar^{-1-\eta}+b(\delta)$ ($b(0)=0$) Sachdev (2007). At the critical point we then expect a power law decay $C_{l}=ar^{-1-\eta}$. As shown in Fig. 4a, the off-diagonal correlation function $\langle b_{i}^{{\dagger}}b_{j}\rangle$ decays very fast in the VBS phase ($t/V=0.1$) which hints at an exponential behavior, while it decays slowly to a constant in the superfluid phase ($t/V=0.137$). Near the critical point ($t/V=0.126$), it shows a clear power law behavior. To approach the thermodynamic limit (TDL), we calculate the correlation function near the critical point ($t/V=0.1303$), and perform a FSS analysis on the exponent. As shown in Fig. 4b, we identify a power law decay with increasing system size. The inset of Fig. 4b shows strong finite size effects of the anomalous exponent, and these size effects can extremely depress the exponent obtained from the data collapse of the condensate fraction sup . With second order polynomial fitting, we get $\eta_{\text{SF}}=0.305(0.020)$ in the TDL. Fig. 4c shows the density correlation function for different parameters. Contrary to the off-diagonal correlations, the density correlation functions have a density modulation due to translational symmetry breaking in the VBS. We thus subtract its mean value and divide by the density modulation $\cos(\textbf{Q}r_{i,j})$. While the correlations show strong fluctuations at short distance, a smooth power law decay emerges at long distances and we thus neglect the first ten points for the fitting. Comparing to the off-diagonal correlations, the error is larger and the anomalous exponent in the TDL is $\eta_{\text{VBS}}=0.313(0.057)$.

From Eq. (4) we identify the spinon density $s_{\bigtriangleup,\bigtriangledown}=\sum_{i\in\bigtriangleup,\bigtriangledown}n_{i}$ as a conserved charge. Such conserved charge should have scaling dimension $\Delta=2$ for any $2+1$D CFT. As shown in Fig. 4d, the corresponding correlation function shows a fast power law decay and rather small finite size effects. The extracted exponent is $-4.02(0.18)$, which is relatively close to $-2\Delta$, strongly supporting the scenario that the easy-plane-NCCP1 is a CFT.

A hallmark of DCPs are emergent symmetries Senthil et al. (2004b); Nahum et al. (2015b); Wang et al. (2017). For example, at the critical point, the lattice $Z_{3}$ rotation symmetry will be enlarged to a continuous $U(1)$ rotation symmetry. To check this, we consider the resonant valence bond order parameter $\Xi(\mathbf{Q}_{h})=\sum_{k=1}^{N_{h}}\Xi_{k}e^{\mathbf{iQ_{h}}\cdot\mathbf{r}_{k}}/\sqrt{N_{h}}$ with $\mathbf{Q_{h}}=(2\pi/3,0)$ where $\Xi_{k}$ is the density operator of a hexagon of a resonant configuration, and $N_{h}$ is the number of hexagons. In the VBS phase, the $Z_{3}$ degeneracy implies the phases of $\Xi(\mathbf{Q}_{h})$ are [math], $2\pi/3$ and $4\pi/3$. We define $\Xi(\mathbf{Q}_{h})_{x}$ ($\Xi(\mathbf{Q}_{h})_{y}$) as the real (imaginary) part $\text{Re}[\Xi(\mathbf{Q}_{h})]$ ($\text{Im}[\Xi(\mathbf{Q}_{h})]$) of $\Xi(\mathbf{Q}_{h})$. From the histogram shown in Fig. 5 we find that it has three peaks in the VBS phase, indicating $Z_{3}$ symmetry, which shrinks to one point in the superfluid phase, reflecting no solid order. Near the critical point, the distribution approaches a uniform circle which reveals emerging $U(1)$ symmetry. In order to quantitatively check this, we introduce a $Z_{3}$ anisotropy parameter $W_{3}=\langle\cos(3\arg(\Xi(\mathbf{Q}_{h})))\rangle$. Fig. 5d shows this quantity to increase in the VBS phase and to vanish in the superfluid phase with increasing system size.

Conclusions and discussions. We have studied the easy-plane deconfined phase transition of a hard-core Bose-Hubbard model using QMC. Finite size simulations of clusters up to $L=72$ indicate an anomalous critical point separating the VBS and superfluid phase. We estimate the critical point is at $t_{c}/V\approx 0.1303$. Following the approach in Ref.Nahum et al. (2015a), we extract the anomalous exponents $\eta_{\text{SF}}\approx 0.305$ and $\eta_{\text{VBS}}\approx 0.313$ from the two-point correlation functions. In addition, we identify a lattice operator for the conserved charge of NCCP1, and we numerically show its scaling dimension is \Delta$$\approx$$2. At last the emergent U(1) rotation symmetry is found at the critical point.

Comparing with another easy-plane NCCP1 model Kuklov et al. (2006), our model can be viewed as a different way to regularize an easy-plane NCCP1 continuous field theory on a discrete lattice. For example, in our system there is only one $U(1)$ global symmetry, while in the paper by Kuklov et. al. there are two $U(1)$ global symmetries. This difference leads to the winding number of two type spinons in our model are equal $W_{-}=W_{\bigtriangleup}-W_{\bigtriangledown}=0$ which means no super-counter fluid phase exists. Such difference may also strongly change the type of phase transition. Altogether, our results strongly support the presence of easy-plane deconfined criticality. Similar to previous works, our data shows some scaling violation that require further studies.

We are thankful for useful discussions with Adam Nahum, Arnab Sen, Yuan Wan, G. J. Sreejith, Wenan Guo, Stefan Wessel, and Chong Wang. This work was supported in parts by the German Research Foundation (DFG) via the Collaborative Research Centers SFB/TR49, SFB/TR173, SFB/TR185, SFB 1143 and Research Unit FOR 1807 through grants no. PO 1370/2-1. The authors gratefully acknowledge the computing time granted by the John von Neumann Institute for Computing (NIC) on the supercomputer JURECA at Jülich Supercomputing Centre (JSC), by the Allianz für Hochleistungsrechnen Rheinland-Pfalz (AHRP) and by the Max-Planck Computing and Data Facility (MPCDF). YCH is supported by a postdoctoral fellowship from the Gordon and Betty Moore Foundation, under the EPiQS initiative, GBMF4306, at Harvard University.

Appendix A Finite size scaling

At the critical point of a second order phase transition, the two-point correlator follows a power lay decay $C_{L}(r)=a(L)r^{-1-\eta(L)}$. In large systems, the structure factor or condensate fraction is proportional to its integration per site:

[TABLE]

As mentioned in main text, we can use $C_{L}(r)$ to fit the two-point correlator function. For the off-diagonal correlator $b_{i}^{\dagger}b_{j}$ in Fig.6, we can find the prefactor $a(L)$ changes less than $\eta(L)$. Then, we use a second order polynomial function $\eta_{f}(L)$ to fit $\eta(L)$:

[TABLE]

which gives $\eta_{0}=0.305(0.020)$. We also considered higher order polynomial fitting, but the third order $\eta_{0}=0.299(0.113)$ has much less accuracy and higher orders are even overfitted.

In order to check how finite size effects change the anomalous exponent, we substitute $\eta(L)$ with fitting function $\eta_{f}$ in Eqn.(6) and get

[TABLE]

where $a_{f}(L)$ is a second order polynomial fitting function for $a(L)$, and condensate fraction should be proportional to $f_{s}$. Then, we analyze the finite size effect of $a_{f}(L)$, $L^{-2}$ and $\eta_{f}(L)$ separately by defining:

[TABLE]

As shown in Fig.7, $f_{s}$ matches well with $\rho_{0}$ after rescaling its magnitude ($\rho_{0}^{RN}=C\rho_{0}$). Both $a_{f}$ and $L^{-2}$ terms can change the shape of the curve. However, the finite size effect of $\eta_{f}$ bends the curve which explains why the anomalous critical exponent obtained from a data collapse of the condensate fraction deviates strongly. A similar phenomenon also happens for the density correlator and structure factor show in Fig.8 and Fig.9. Therefore the size independent critical exponents got directly from structure factor and condensate fraction are less convincing.

In addition, we also consider the data collapse of superfluid density with two-length scale scenario. As shown in Fig.10, it is reasonable as good as LC scenario, so we can not conclude which one is better.

Appendix B numerical flowgram

The numerical flowgram method was introduced by A.B. Kuklov, et. al. Kuklov et al. (2006, 2008); Chen et al. (2013b) to study the DCPs. They obtain the finite size critical point $t_{c}(L)/V$ from the condition that the ratio of probabilities of having zero and non-zero winding numbers is some fixed number of the order of unity. If the transition is continuous, they claim the winding number at this critical point will approach a universal value when enlarging the system, otherwise, it will linearly scale with system size $L$.

We also implement numerical flowgram method for our case. We fix the probability of zero winding numbers $P_{0}$ equal to $0.25$ which marked in Fig. 11. We determine the size-dependent hopping values $t_{P_{0}}(L)$ for this probability. The corresponding winding numbers for those parameters can then be analyzed as a function of length $L$. As shown in inset of Fig.11, the winding numbers don’t flow to a universal value, but seem to linearly depend on the system size. Similar behaviour is also found in the J-Q model Chen et al. (2013b), and it may be directly related to the drift of the superfluid density or weakly first order phase transition.