Kibble-Zurek exponent and chiral transition of the period-4 phase of Rydberg chains

Chains of Rydberg atoms have emerged as an amazing playground to study quantum physics in 1D. Playing with inter-atomic distances and laser detuning, one can in particular explore the commensurate-incommensurate transition out of density waves through the Kibble-Zurek mechanism, and the possible presence of a chiral transition with dynamical exponent z > 1. Here, we address this problem theoretically with effective blockade models where the short-distance repulsions are replaced by a constraint of no double occupancy. For the period-4 phase, we show that there is an Ashkin-Teller transition point with exponent ν = 0.78 surrounded by a direct chiral transition with a dynamical exponent z = 1.11 and a Kibble-Zurek exponent μ = 0.41. For Rydberg atoms with a van der Waals potential, we suggest that the experimental value μ = 0.25 is due to a chiral transition with z ≃ 1.9 and ν ≃ 0.47 surrounding an Ashkin-Teller transition close to the 4-state Potts universality.

U nderstanding the nature of quantum phase transitions in low-dimensional systems is one of the central topics in condensed matter physics 1,2 . Over the last decades, the combination of conformal field theory in 1+1D 3,4 and advanced numerical techniques such as the density matrix renormalization group (DMRG) algorithm [5][6][7][8] has proven to be extremely powerful in coming up with theoretical predictions for numerous fascinating critical phenomenas. In that respect, modern quantum simulators based on Rydberg atoms trapped with optical tweezers offer a remarkably rich experimental playground to further investigate quantum physics in 1D. In particular, in a recent experiment 9 , the quantum critical dynamics of a chain of Rubidium atoms 87 Rb with programmable interactions has been probed. The atoms are excited to a Rydberg state by a homogeneously applied laser with Rabi frequency Ω, while the laser detuning Δ controls the population of excited atoms. The quantum many-body Hamiltonian of the system can be written in terms of hard-core bosons (i.e. bosons with no more than one particle per site) as where the van der Waals potential between Rydberg states decays as The competition between the detuning term Δ that favors a high density of Rydberg states and the blockade leads to a sequence of lobes of density-wave phases with fixed periodicities. In general, these periodicities can be any rational number, and in the classical limit Ω → 0 they form a Devil's staircase 10 , but at finite values of Ω, the phase diagram is dominated by lobes of integer periodicities p = 2, 3, 4, … 9,11,12 , surrounded at least partially by a critical floating phase for p ≥ 3 (ref. 11 ). However, according to recent experiments in which the detuning frequency has been swept for various interatomic distances a, this floating phase cannot be present along the whole boundary for p = 3 and 4 since a direct transition with a non-integer dynamical exponent z larger than 1 has been detected in the vicinity of the tip of the lobe 9 .
The transition out of a period-p phase is an example of commensurate-incommensurate transition, a problem with a long history that goes back to the investigation of adsorbed monolayers on surfaces [13][14][15] . In these systems the role of Rydberg atoms is played by domain walls between periodic phases, and quite remarkably the melting of these periodic phases is a very subtle problem that has not yet received a full solution. For p = 2, the transition is known to be generically Ising, while for p ≥ 5 it is a two-step process through a Luttinger liquid phase (called a floating phase in the context of adsorbed monolayers) if it is not first order. The difficult cases are precisely p = 3 and p = 4. In these cases, along the commensurate line, i.e. the line along which, in the disordered phase, the wave vector keeps the value q = 2π/p of the ordered phase, the transition is expected to be continuous in the three-state Potts universality class for p = 3, and in the Ashkin-Teller universality class (see below) for p = 4. Away from this line, the disordered phase is incommensurate. As pointed out by Ostlund 16 and Huse 17 , this introduces a chiral perturbation, and the open problem is to understand the effect of this chiral perturbation on the transition.
For p = 3, the chiral perturbation is always relevant, and the question is whether it immediately opens a floating phase away from the Potts point, or whether the transition remains direct and continuous for a while, but in a new chiral universality class, as suggested by Huse and Fisher 18 , with a dynamical exponent z > 1. Numerical [19][20][21][22][23] and experimental evidence 14,15 in favor of this possibility has been obtained in the 1980s and early 1990s in the context of adsorbed layers, and very recently in the context of Rydberg atoms 9,12,24-26 . For p = 4, the situation is even richer because the chiral perturbation is not always relevant 13 . With four degrees of freedom, there is in fact a family of universal classes described by the Ashkin-Teller model in which the local degrees of freedom are described by two Ising spins σ i ⊗ τ i coupled by an interaction σ i σ j + τ i τ j + λσ i τ i σ j τ j . The asymmetry parameter λ controls the relevance of the chiral perturbation. Indeed, according to Schulz 27 , the crossover exponent ϕ of the chiral perturbation for the Ashkin-Teller model is given by where ν is the exponent of the correlation length The chiral per- ' 0:683, irrelevant otherwise. Now, the exponent ν is known exactly as a function of λ 28,29 : At λ = 0, the model is known as the four-state clock model and corresponds to two decoupled transverse-field Ising chains. In that case, ν = 1: The chiral perturbation is relevant, and it is known to drive the system immediately into a critical phase. By contrast, at the four-state Potts model (λ = 1), ν = 2/3 < ν c : The chiral perturbation is irrelevant. The critical value of λ below which the chiral perturbation becomes relevant is given by As long as the chiral perturbation is irrelevant, a line of continuous transition in the Ashkin-Teller universality class can be expected around the commensurate line until the chiral perturbation becomes relevant. Then the situation is similar to the p = 3 case, with again the possibility of a chiral transition before a floating phase appears, as emphasized by Huse and Fisher 30 . These two possibilities are summarized in the generic phase diagrams of Fig. 1. In this figure, λ denotes the parameter of the Ashkin-Teller model that describes the transition along the commensurate line, and δ stands for the amplitude of the chiral perturbation. For Rydberg atoms, δ should be understood as the distance to the Ashkin-Teller point along the transition into the disordered phase. For the Ashkin-Teller model itself, the chiral perturbation is given by δ(σ i τ j − τ i σ j ), the form used in ref. 27 to derive the crossover exponent ϕ. Now, for p = 4, the experimental results on Rydberg atoms 9,12 are compatible with a continuous transition, with a Kibble-Zurek exponent μ ≃ 0.25. This exponent is related to ν by the relation μ = ν/(1 + νz), where z is the dynamical exponent. Since between the clock model (λ = 0) and the Potts model (λ = 1) the exponent ν decreases from 1 to 2/3, the Kibble-Zurek exponent should be between 1/2 and 2/5 if the dynamical exponent was equal to 1. So, according to these experiments, the dynamical exponent has to be larger than 1. This implies that the transition should be a chiral Huse-Fisher transition, i.e., that the scenario of Fig. 1a is realized.
In this paper, we investigate this problem in the context of an effective model for the period-4 phase, the quantum hard-boson model with two-site blockade (see below). We show that: (i) the transition along the commensurate line is sufficiently far from the four-state Potts point to ensure that the chiral perturbation is relevant; (ii) there is an intermediate floating phase far enough from this point; (iii) there is evidence in favor of a small region of chiral transition in between for which we estimate the dynamical exponent and the Kibble-Zurek exponent. Implications for the original model of Eq. (1) and for the experiments on Rydberg atoms are also discussed.

Results
The blockade models. Because of the very steep increase of the van der Waals potential at short distance, the simultaneous excitation of atoms at a distance smaller than the so-called Rydberg blockade radius R b ðV 1 =ΩÞ 1=6 is essentially excluded, a phenomenon known as Rydberg blockade. This means that, on a chain with lattice parameter a, the interaction between sites i and j can be considered to be infinite if i − j ≤ r, where r is the largest integer satisfying r < R b /a. Keeping only the dominant next-to-blockade repulsion leads to the following effective Hamiltonian: where d i (d y i ) is an annihilation (creation) operator that acts in a constrained Hilbert space: n i ðn i À 1Þ ¼ n i n iþ1 ¼ ::: ¼ n i n iþr ¼ 0: We will refer to this model as the r-site blockade model. When r = 1, it reduces to the original hard-boson model introduced by Fendley et al. 31 . Note also that a constrained Hilbert space equivalent to r = 2 has been introduced by Huijse et al. 32 in the context of a supersymmetric model on a zig-zag ladder. Quite generally, the r-site blockade model allows one to discuss period p = r + 1 and p = r + 2 phases and their surrounding (see Supplementary Note 1). The main advantage of these constrained models is that their Hilbert space grows much more slowly than that of the original model of Eq. (1), and simulations can be performed on systems large enough to keep track of small changes in the incommensurability and to identify the critical behavior at the transition. Interestingly, the r-site blockade model can be seen as the limit of Rydberg atoms with infinitely fast decaying interactions. Indeed, if we consider the model of Eq. (1) with the interaction the r-site blockade model corresponds to m → + ∞ while the 1/R 6 model of Eqs. (1) and (2) is recovered for m = 6.
Overview of the phase diagram for p = 4. As a first step towards the period-4 phase of Rydberg atoms, let us now turn to the properties of the two-site blockade model. Our numerical results have been obtained with a state-of-the-art DMRG algorithm 5-8 that explicitly implements the constraints (see "Methods" for details about the algorithm). They are summarized in the phase diagram of Fig. 2. There are three main phases: a disordered phase with incommensurate short-range correlations, and two ordered commensurate phases with period 3 and 4, respectively. Note that these three main phases have been accessed in recent Rydberg atom experiments 9,12 . There are also small floating phases close to the ordered phases. In particular, for large values of Δ, there are two floating phases at the boundaries of the period-three and period-four phases that come closer and create an area of extremely high correlation length. It is therefore probable that the disordered phase eventually disappears and that, for some parameter range, the two ordered phases are connected through a single floating phase, as suggested in refs. 11,12 for the model of Eq. (1). Due to the exponential growth of the correlation length at the Kosterlitz-Thouless 33 phase transition, an accurate investigation of this scenario would require simulations beyond our current limitations.
Commensurate line. The transition out of the period-four phase is the main focus of the rest of this section. Our first task is to locate the point on the phase boundary where the chiral perturbation vanishes, hence where the transition can be expected to be described by a conformal field theory. Note that for the original hard-boson model of Fendley et al. 31 this was not necessary because the three-state Potts belongs to an integrable line, and its location is known exactly. Here this is not the case, but we can expect this point to be located at the intersection of the phase boundary and of the line with wave vector q = π/2 since along this line the correlations remain commensurate in the disordered phase so that there is no chiral perturbation. To achieve this goal, we have extracted the wave vector q (see "Methods") and have mapped the results over the disordered phase to determine the lines of constant incommensurate wave vector q. These lines are depicted in Fig. 3a. The line q = π/2 enters the period-four phase at Δ/Ω ≃ 1.593. An accurate estimate of the second coordinate has been obtained by a finite-size scaling of the order parameter. It turns out that open boundary conditions favor a boson on the first and last sites. This effectively acts as a conformally invariant fixed boundary condition at the critical point and induces Friedel oscillations in the local boson density. According to boundary conformal field theory, the profile of these oscillations on a finitesize chain is proportional to ½N sinðπj=NÞ Àd , where the scaling dimension d = 1/8 for the Ashkin-Teller model 34 . By scanning V 3 /Ω for Δ/Ω = 1.593, we identify a separatrix in the log-log scaling at V 3 /Ω = 1.2839 as shown in Fig. 3b. The slope corresponds to d ≃ 0.124, in excellent agreement with the scaling dimension d = 1/8. As a further check that this is a critical point, we have extracted the central charge by fitting the profile of the reduced entanglement entropy to the Calabrese-Cardy formula (see "Methods"), leading to a central charge c ≃ 0.96, within 4% of the conformal field theory prediction c = 1. The correlation length of the hard-boson model can be simply obtained by fitting correlations, a straightforward task along the commensurate line (see "Methods"). The resulting correlation diverges at the critical point with an exponent ν ≃ 0.78. This is the first indication that λ must be significantly smaller than 1. This is actually quite natural. Indeed, when λ = 1, the model corresponds to the four-state Potts model with the same amplitude for all flipping processes, while for λ < 1 two processes are favored over the third one by the transverse field term. Such an asymmetry naturally appears in the hard-boson model due to the two-site blockade. In the p = 4 phase every fourth site is occupied by a boson. So each of the ground states, let us call them A, B, C, and D, corresponds to the location of the occupied sites mod 4. From Fig. 4 one can see that domains B and D shifted by one site with respect to the bulk A cost less energy than the domain C shifted by two sites.
One can also estimate λ directly by comparing the excitation spectrum of the two-site blockade model with that of the quantum 1D version of the Ashkin-Teller model (see "Methods"). This leads to λ ≃ 0.57, in excellent agreement with ν = 0.78. At that point, the chiral transition is relevant, with a crossover exponent ϕ ≃ 0.33. This means in particular that, away from that point, the transition cannot be a standard continuous transition in the Ashkin-Teller universality class. Either a floating phase opens or the transition becomes chiral.
Chiral transition versus floating phase. Quite generally, the incommensurate wave vector q is expected to approach the commensurate value π/2 with a critical exponent called β. To the best of our knowledge the exact value of this critical exponent is not known for the Ashkin-Teller model, but Huse and Fisher 30 argue that β > ν. This implies that the product ξ × |π/2 − q| decays to zero upon approaching the Ashkin-Teller transition. By contrast, if the transition is chiral, the equality β ¼ ν should hold, and ξ × |π/2 − q| is expected to go to some finite value 30 . When the transition is Ashkin-Teller or chiral, the exponents of the correlation length ν in the disordered phase and ν 0 in the ordered phase should satisfy ν ¼ ν 0 . By contrast, in the presence of an intermediate floating phase, the correlation length in the disordered phase diverges exponentially at a Kosterlitz-Thouless 33 transition, while the wave vector q remains incommensurate, so that the product ξ × |π/2 − q| diverges. The commensurate-incommensurate transition between the floating and the ordered phases is then expected to be in the Pokrovsky-Talapov 35 universality class with critical exponent β ¼ ν 0 ¼ 1=2.
In Fig. 5 we take a closer look at three cuts across the transition. Let us start with the vertical cut through the Ashkin-Teller point identified above at Δ/Ω = 1.593. The critical exponents ν and ν 0 are in good agreement with each other, and they are also in reasonable agreement with the value obtained for ν along the commensurate line and with the value of λ. An accurate estimate of β is very difficult due to the proximity of the commensurate value of q in the disordered phase. Nevertheless it is clear qualitatively, just looking at the curvature, that β is significantly larger than ν, in agreement with Huse and Fisher 30 . As a consequence, the product ξ × |π/2 − q| goes to zero at the critical point as shown in Fig. 5c.  For p = 4, domains with B or D inside A cost an energy V 3 while domains with C inside cost an energy Δ > V 3 since there is one particle less, leading to an asymmetry in the effective transverse field term.
The next cut at V 3 /Ω = 1.35, slightly away but very close to this Ashkin-Teller point, is presented in Fig. 5d-f. The correlation length diverges as a power law with similar exponents on both sides of the transition, but, by contrast to the Ashkin-Teller point, the critical exponent β is much smaller than 1, a clear indication that the chiral perturbation changes the physics immediately away from the Ashkin-Teller point. Its value is comparable to ν and ν 0 , and accordingly, even if it increases slightly towards the transition, the product ξ × |π/2 − q| seems to remain finite. The absence of divergence of the product ξ × |π/2 − q| is a clear indication in favor of the Huse-Fisher universality class. However, as in the case p = 3, an extremely narrow floating phase cannot be excluded.
Further away from the commensurate point, the inverse of the correlation length decays in a very asymmetric way, as we show for the horizontal cut at V 3 /Ω = 3.5 in Fig. 5g-i. The numerically extracted critical exponent ν 0 is in reasonable agreement with the Pokrovsky-Talapov value 1/2, while the product ξ × |π/2 − q| clearly diverges towards the transition. By fitting the divergence of the correlation length ξ with the predictions for Kosterlitz-Thouless and Pokrovsky-Talapov transitions, we estimate the width of the floating phase to be dΔ/Ω ≈ 3 × 10 −3 . The physics is very similar on the other side of the Ashkin-Teller line (see Supplementary Note 3 for more data).
As a further check, we have investigated the behavior of the second derivate of the ground-state energy, the equivalent of the specific heat for quantum systems. If the transition is continuous, it is expected to diverge with the same exponent α on both sides of the transition, while if there is an intermediate floating phase it is expected to diverge with exponent 1/2 at the Pokrovsky-Talapov transition when coming from the incommensurate phase, and to saturate with a logarithmic singularity on the other side 30 . As can be seen in Fig. 6, the results are fully consistent with a single transition at and close to the commensurate line, and with an asymmetric behavior far enough from it. According to hyperscaling, α should be related to ν at the Ashkin-Teller point by α = 2 (1 − ν) ≃ 0.44, in good agreement with the numerical results of Fig. 6a. Interestingly, α barely changes as long as the transition is continuous, a fact already noticed and rigorously established for integrable and selfdual versions of the three-state chiral Potts model [36][37][38] .
Kibble-Zurek mechanism and dynamical exponent. To estimate the Kibble-Zurek exponent μ = ν/(1 + νz), we need both the dynamical exponent z and the correlation length exponent ν. Along the commensurate line, the transition is in the Ashkin-Teller universality class and has conformal invariance, so z = 1. The estimate ν = 0.78 then leads to a Kibble-Zurek exponent μ = 0.44. Away from the Ashkin-Teller point, the dynamical exponent z can Inverse correlation length 1/ξ, wave vector q/π, and product ξ × |π/2 − q| along three different cuts across the transition. a-c Vertical cut through the Ashkin-Teller point at Δ/Ω = 1.593; d-f, g-i Horizontal cuts at V 3 /Ω = 1.35 and V 3 /Ω = 3.5, respectively. Inside the p = 4 phase, the correlation length is fitted with a power law with critical exponent ν 0 . In the disordered phase, the correlation length is fitted either with a power law with critical exponent ν (a, d) or with the KT form ξ / expðC = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi g KT À g p Þ (g), where g is the coordinate along the cut. The wave vector q is fitted with a power law with exponent β (dotted lines). Wave vectors q are defined within the error bars ±πξ/N 2 ; and ξ × |π/2 − q| is defined up to ±πξ 2 /N 2 . For points without error bars, the error bar is smaller than the size of the symbol. In the lower panels, the red lines indicate the boundary of the ordered phase. NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-20641-y ARTICLE NATURE COMMUNICATIONS | (2021) 12:414 | https://doi.org/10.1038/s41467-020-20641-y | www.nature.com/naturecommunications be extracted from ν and α according to the hyperscaling relation for anisotropic systems: ν + zν = 2 − α. For the cut of Fig. 5d, ν = 0.74. Assuming that α keeps the value α = 0.44 of the Ashkin-Teller point, in agreement with the discussion above, leads to z = 1.11. Across this cut, the Kibble-Zurek exponent is thus given by μ = 0.41, smaller than across the Ashkin-Teller transition. This conclusion remains true even if we assume that α slightly increases away from the Ashkin-Teller point, as suggested by Fig. 6.

Discussion
To make contact with experiments, let us briefly discuss the implications of the present results for the model with 1/R 6 longrange interactions. As explained above, both models belong to the same family of models defined by Eq. (8). Let us estimate the effect of reducing m from +∞ to 6 for r = 2 in the vicinity of the Ashkin-Teller transition. The critical values of the two-site blockade are given by Δ/Ω = 1.593 and V 3 /Ω = 1.2839. This value of V 3 corresponds to a Rydberg blockade radius R b =a ¼ 3ðV 3 =ΩÞ 1=6 ¼ 3:1276, near the tip of the p = 4 lobe where the experiments have been carried out 9 , and where there is no evidence of a floating phase 11 . The critical value of Δ/Ω = 1.593 is different from that of the Rydberg model at the tip of this lobe (around 2.39), but this is not surprising since Δ is the chemical potential in the bosonic language, and its critical value must be strongly affected by the details of the interactions. Now let us turn to the nature of the transition, assuming that there is a portion of boundary without floating phase. The physical reason behind λ < 1 is the difference in energy cost of domains shifted by one or two sites with respect to the bulk (see Fig. 4). In the very simple classical approximation δE B,D ≃ V 3 while δE C ≃ Δ. At the Ashkin-Teller critical point, these expressions lead to δE B,D ≃ 1.2839 and δE C ≃ 1.593 for the blockade model. Taking into account longer-range interactions, the energy of domain walls for the Rydberg model can be estimated as δE B,D ≃ V 3 − 2V 4 + V 5 and δE C ≃ Δ − 3V 4 + 2V 7 . Assuming V 3 = V ≃ 1.2839 and a 1/R 6 decay, we get δE B,D ≃ 0.885 and δE C ≃ 0.925. The asymmetry is still present, but it is smaller, implying that the point where the chiral perturbation vanishes gets closer to the four-state Potts point. Therefore, there are two possibilities: (i) λ is still smaller than λ c = 0.9779. Then the chiral perturbation remains relevant, and the transition immediately becomes chiral until a floating phase emerges; (ii) the long-range interactions bring the Ashkin-Teller point close enough to the four-state Potts point so that the chiral perturbation is irrelevant; then there will be an extended region of direct Ashkin-Teller transition, followed on both sides by a chiral transition, and ultimately by a floating phase. Since λ c ≃ 0.9779 is very close to 1, the first possibility (i) is more likely. More importantly, the fact that the asymmetry can be expected to be reduced by long-range interactions and not increased implies that, if anything, the Rydberg model is further away from the clock limit λ = 0 where there would be an intermediate floating phase all along the boundary. So our conclusion that there is a portion of the boundary to the period-4 phase where the transition is direct and continuous in the chiral universality class before a floating phase opens can be considered as a prediction for the Rydberg model with 1/R 6 interactions.
Note that the finite-size effects associated with the restricted number of Rydberg atoms in experiments 9,12 will, if anything, enlarge the portion without the floating phase. Indeed, if, coming from the disordered phase, the floating phase starts at an incommensurate wave vector q, its detection requires the size of the chain to be significantly larger than the period necessary to form at least one helix N > 2π/(q − π/2). So, for a finite-size system, the floating phase can only be detected further away from the commensurate line than in the thermodynamic limit, and the transition will look continuous in a larger parameter range, making the observation of this direct transition easier. Finally, let us discuss briefly the consequences for the Kibble-Zurek experiment. If the transition is chiral, the scaling becomes anisotropic, but if hyperscaling applies, the correlation exponent along the chains ν and the dynamical exponent z are related by ν(1 + z) = 2 − α. Let us further assume that, as for the two-site blockade model and the self-dual three-state chiral Potts model, the specific heat exponent keeps the value it has at the Ashkin-Teller critical point α = 2(1 − ν λ ). Then we get ν(1 + z) = 2ν λ , where the Ashkin-Teller critical exponent ν λ is given by Eq. (4). This implies that ν and z can be deduced from the Kibble-Zurek exponent μ and the asymmetry parameter λ according to ν ¼ μð1 þ 2ν λ Þ ½ =ð1 þ μÞ and z ¼ ð2ν λ À μÞ= μð1 þ 2ν λ Þ ½ . Taking the experimental value μ ≃ 0.25 and assuming that λ is close to 1, as suggested by the small asymmetry of domain walls for Rydberg atoms, we get z ≃ 1.9 and ν ≃ 0.47. It will be interesting to see if these values can be confirmed by a direct numerical investigation of the model of Eqs. (1) and (2).

Methods
Details about the algorithm. The size of the Hilbert space for a model with twosite Rydberg blockade can be calculated using a recursive relation HðNÞ ¼ HðN À 1Þ þ HðN À 3Þ, with the first three elements of the sequence Hð1Þ ¼ 2, Hð2Þ ¼ 3, and Hð3Þ ¼ 4. So the growth of the Hilbert space with the system size HðNÞ / 1:466 N is much slower than HðNÞ / 2 N for an unconstrained model. In order to fully profit from the restricted Hilbert space we implement the blockade explicitly into the DMRG. Recently it has been shown that the hardboson model with r = 1 can be rigorously mapped onto a quantum dimer model on a two-leg ladder 39 that provides a simple and intuitive way to encode the constraint into DMRG. Although this mapping is not valid for r > 1, we can rely on the idea of auxiliary quantum numbers that would preserve the block-diagonal structure of the local tensors. This is achieved by a rigorous mapping onto an effective model that spans the local Hilbert space over three consecutive sites on the original lattice as shown in Fig. 7b. The new local Hilbert space contains four states listed in Fig. 7c. Because of the overlap, the three possible states of two shared sites can be used as a quantum label for the auxiliary bond between two consecutive sites of the new model. By adding a site, for example, by increasing the left environment, one can change the quantum labels according to the fusion graph shown in Fig. 7d. The fusion graph for the right environment can be obtained by inverting the arrows. An example of the label assignment is provided in Fig. 7e.
At the next step, one has to rewrite the hard-boson model given by Eq. (1) of the main text in terms of new local variables h i j i. For example, the boson occupation number operator n i , which is also equal to (1 − n i−1 )n i (1 − n i+1 ), can be written in the new local Hilbert space as a 4 × 4 matrixñ i with the only non-zero elementñ i ð3; 3Þ ¼ 1. The term V 3 n i−1 n i+2 can be written in the new Hilbert space as a nearest-neighbor interaction V 3piqiþ1 where the only non-zero matrix elements of the operatorsp andq are given bypð4; 4Þ ¼ 1 andqð2; 2Þ ¼ 1. Finally the constrained flip term À Ω 2 ð1 À n iÀ2 Þð1 À n iÀ1 Þðd y i þ d i Þð1 À n iþ1 Þð1 À n iþ2 Þ can be rewritten as a three-site operator where the only non-zero matrix elements of the operatorsã,b, andc are given bỹ að1; 2Þ ¼ 1,bð1; 3Þ ¼ 1, andcð1; 4Þ ¼ 1.
With these definitions, the matrix product operator in the bulk takes the following simple form:Ĩ : : : : : : q : : : : : : c : : : : : : c y : : : : : : : :b : : : : : : : : where dots mark zero entries of the tensor. Close to the edges one has to carefully modify the MPO to properly encode the boundary terms. This requires the definition of local operators slightly different from those used in the bulk.
There is yet another crucial point that we want to mention. The labels that we have introduced split the Hilbert space into blocks or sectors and therefore correspond to some conserved quantity. For the hard-boson model with a singlesite blockade, the quantum labels correspond to the parity of the domain walls. In the present case, the physical meaning of the conserved quantity is not as obvious. However, the only relevant information for us is that the conservation of this abstract quantity requires at least three sites. In other words, by acting with any term (read flip term) on a two-site MPS, one necessary changes one of the outgoing labels, while the flip term applied on three consecutive MPS keeps all external labels fixed. As a consequence, neither single-nor two-site DMRG routines are compatible with the presented constraint implementation, and one has to go for at least three-site updates. At a glance this might look costly with a local Hilbert space of dimension 4 since it leads in principle to an MPO operator of size 7 × 7 × 64 × 64. However, taking into account all the constraints on three sites, the projected three-site MPO is only of size 7 × 7 × 9 × 9.
The explicit implementation of two-sites blockade allows us to reach systems with up to N = 3001 sites systematically (and N = 4801 sites occasionally), keeping up to 2000 states. where dðlÞ ¼ 2L π sin πl L À Á is the conformal distance; s 1 and log g are non-universal constants. The presence of Friedel oscillations caused by the fixed boundary conditions is also reflected in the entanglement entropy profile. In order to remove the oscillations we follow ref. 41 and construct the reduced entanglement entropy: where ζ is a non-universal constant in front of the leading local correlations between nearest allowed neighbors adjusted to best remove the oscillations. The fits are performed using sites sufficiently far from the edges (l, L − l ≫ 1).
Comparison with the Ashkin-Teller model and estimate of λ. To estimate λ directly, one can compare the excitation spectrum of the two-site blockade model with that of the quantum 1D version of the Ashkin-Teller model defined in terms of Pauli matrices σ x,z and τ x,z by the Hamiltonian: at its critical point β = 1. The spectra have been obtained by targeting several states (up to 11) at every DMRG iteration 42 . In Fig. 8a we show the energy spectrum of the Ashkin-Teller model for N = 60 with fixed A-A boundary conditions. We compare these results with the spectrum of the constrained model with N = 201 sites. Since the velocity is a non-universal constant, one cannot compare the absolute values of the gap. However we find that the structure of the spectrum in the hard-boson model corresponds to the structure of the Ashkin-Teller spectrum at λ ≃ 0.57 (red line in Fig. 8a). In Fig. 8c we further compare the finite-size scaling for hard-boson (red) and Ashkin-Teller model at λ = 0.57 (green) and at λ = 1 (four-state Potts, blue) and the agreement with λ = 0.57 is quite good. The comparison can be made even more systematic by re-scaling both spectra with respect to the lowest excitation energy as explained in the Supplementary Note 2.
We compute the energy spectrum in a chain with open and fixed boundary conditions. There are two reasons for that. First, DMRG is well known to be more efficient for open boundary conditions than for periodic ones. Second, the number of conformal towers of states that appears in the spectra of periodic or anti-periodic chains are usually larger than the number of towers selected by fixed boundary conditions. However, we have to establish the correspondence between the different boundary conditions in the hard-boson model and in the original Ashkin-Teller model. In the hard-boson model, the simplest way to fix the boundary is to force the first and the last sites to be occupied. If the total number of sites is 4k + 1 the same state is favored at each edge, corresponding to the A-A boundary condition in the Ashkin-Teller model. If the total number of sites is 4k or 4k + 2, we expect A-B and A-D boundary conditions. They are expected to give the same spectrum (assuming that states B and D have equal weight in the transverse field applied on A, while C has a factor λ). Finally, if the total number of sites is 4k + 3, we expect to observe the spectrum of the A-C boundary condition. Numerical results for A-C and A-B/A-D boundary conditions are provided in the Supplementary Note 2.
We extract the correlation length critical exponent along the commensurate line which, close to the transition, is given by V 3 /Ω = 0.3645Δ/Ω + 2.825. Since we expect a direct transition the critical exponent has to be the same on both sides of the critical point. However, the pre-factor is non-universal. We therefore fit our numerical data with: jx À Δ c j ν ½aθðx À Δ c Þ þ bθðΔ c À xÞ; where a, b, Δ c , and ν are fitting parameters; and θ(x) is the Heaviside function: θ(x) = 1 if x > 0 and zero otherwise. The results are presented in Fig. 9a.
We compare the values of λ and ν obtained to fit the hard-boson model with the conformal field theory result of Kohmoto et al. 28,29 in Fig. 9b. The agreement is very good.
Extraction of the critical exponent α. In order to extract the specific heat critical exponent α we look at the divergence of the second derivative of the ground-state energy density d 2 e/d(Δ/Ω) 2 . In order to check the finite-size effects we take the energy per site e = E/N extracted from the total ground-state energy E for finite chains with two values of the number of sites: N = 1201 and N = 2101. We also consider the difference between the two ground-state energies E 2101 − E 1201 to suppress the edge effects and get a better estimate for the bulk energy per site as (E 2101 − E 1201 )/900. Approaching the Ashkin-Teller point the specific heat should diverge as |Δ − Δ c | α . The effective exponent α eff close to the transition can thus be obtained as the slope of log d 2 e=dðΔ=ΩÞ 2 with respect to log jΔ À Δ c j. The results are presented in Fig. 6a, where the pink line shows the location of the critical point. According to the hyperscaling relations α = 2 − 2ν and to our estimate of the correlation critical exponent ν ≈ 0.78 ± 0.02, the specific heat critical exponent is expected to be α ≈ 0.44 ± 0.04. This corresponds to the gray area in Fig. 6a, showing that our results for α are in reasonable agreement with this estimate at the Ashkin-Teller point.
Far enough from the Ashkin-Teller point, the transition is expected to take place through an intermediate floating phase. At the Pokrovsky-Talapov point, the second derivative of the energy is expected to be very asymmetric, with a divergence with exponent 1/2 on the incommensurate side and no divergence on the commensurate side 30 . The results of Fig. 6d obtained for V 3 /Ω = 2 are in good agreement with these predictions.
Extraction of the correlation length and of the wave vector. In order to extract the correlation length and the wave vector q, we fit the boson-boson correlation function to the Ornstein-Zernicke form 43 where the correlation length ξ, the wave vector q, and the initial phase φ 0 are fitting parameters. In order to extract the correlation length and the wave vector with a sufficiently high precision, we fit the correlation function in two steps. First, we discard the oscillations and fit the main slope of the decay as shown in Fig. 10. This  allows us to perform a fit in a semi-log scale logCðx ¼ ji À jjÞ % c À x=ξ À log ðxÞ=2, which in general provides more accurate estimates of the correlation length on a long scale. Second we define a reduced correlation functioñ and fit it with a cosineC i;j % a cosðqji À jj þ φ 0 Þ as shown in Fig. 10b. The agreement is almost perfect: The DMRG data (blue dots) are almost completely behind the fit (red dots). Fitting the correlations over different windows shows that the error on the correlation length does not exceed 3%, and that the wave vector q is determined with a precision O(10 −6 ). However, the main source of error in the case of q is not the fit itself, but finite-size effects. If the correlation length was infinite, q would exhibit finite-size steps of width 2π/N, leading to an error bar of π/N. But if the correlation length is smaller than the number of sites, this is a clear overestimate. Indeed the q vector adapts close to the boundary, and the steps in the q vector in the bulk are significantly rounded and disappear in the limit of small correlation length. To take this effect into account, we include a factor ξ/N into the error, leading to an error bar of the order πξ/N 2 .

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
The code that support the findings of this study is available from the corresponding author upon reasonable request.
Received: 9 July 2020; Accepted: 11 December 2020; Fig. 10 Example of fit of the correlation function to the Ornstein-Zernicke form. In the first step a, we extract the correlation length discarding the oscillations. In the second step b, we fit the reduced correlation function to extract the wave vector q.