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_{\mathrm{b}}\equiv {\left(V_1/\mathrm{\Omega }\right)}^{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:

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³³ phase transition, an accurate investigation of this scenario would require simulations beyond our current limitations.

Fig. 2

Phase diagram of the hard-boson model with two-site blockade.

At the commensurate transition point located at Δ/Ω ≃ 1.593 and V₃/Ω ≃ 1.2839 the transition is in the Ashkin–Teller universality class with λ ≃ 0.57 (open green circle). Away from it but not too far (for V₃/Ω ≲ 1.8 and for Δ/Ω ≲ 1.7), our results are consistent with a chiral transition in the Huse–Fisher³⁰ universality class. Further away from the Ashkin–Teller point we detect intermediate floating phases bounded by Pokrovsky–Talapov (PT) and Kosterlitz–Thouless (KT) transitions. Above the Ashkin–Teller point the width of the floating phase is always smaller than the width of the gray line. The transition into the p = 3 phase is a direct chiral transition from the disordered phase at small Δ/Ω and a PT transition from the floating phase at large Δ/Ω²⁶. For large Δ/Ω, the two floating phases adjacent to the p = 3 and p = 4 phases eventually merge into a single floating phase connecting the two ordered phases. Dotted lines are constant correlation length lines with ξ = 50 (yellow), 100 (purple), and 200 (green).

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.³¹ 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 finite-size chain is proportional to $$ {\left[N\sin \left(\pi j/N\right)\right]}^{-d}$$, where the scaling dimension d = 1/8 for the Ashkin–Teller model³⁴. By scanning V₃/Ω for Δ/Ω = 1.593, we identify a separatrix in the log–log scaling at V₃/Ω = 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.

Fig. 3

Identification of the conformal point.

a Phase diagram with equal-q lines in the disordered phase extracted for N = 601 (systematically) and for N = 1201 (in the vicinity of the critical lines). Colors used for equal-q lines are guides to the eye. The location of the Ashkin–Teller point has been determined as the crossing point of the critical (green) line and the q = π/2 line (black, open circles). b Finite-size scaling of the amplitude of the oscillations in on-site boson density in the middle of the finite-size chain. The separatrix corresponds to the critical point. Inset: Scaling of the entanglement entropy with the conformal distance d(n) after removing the Friedel oscillations, leading to a central charge c ≃ 0.96.

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 $$ \mathrm{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.

Fig. 4

Asymmetry of domain walls in the model with two-site blockade.

For p = 4, domains with B or D inside A cost an energy V₃ while domains with C inside cost an energy Δ > V₃ since there is one particle less, leading to an asymmetry in the effective transverse field term.

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 $$ \overline{\beta }$$. 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³⁰ argue that $$ \overline{\beta }>\nu $$. 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 $$ \overline{\beta }=\nu $$ should hold, and ξ × ∣π/2 − q∣ is expected to go to some finite value³⁰. When the transition is Ashkin–Teller or chiral, the exponents of the correlation length ν in the disordered phase and $$ {\nu }^{\prime }$$ in the ordered phase should satisfy $$ \nu ={\nu }^{\prime }$$. By contrast, in the presence of an intermediate floating phase, the correlation length in the disordered phase diverges exponentially at a Kosterlitz–Thouless³³ 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³⁵ universality class with critical exponent $$ \overline{\beta }={\nu }^{\prime }=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 $$ {\nu }^{\prime }$$ 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 $$ \overline{\beta }$$ 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 $$ \overline{\beta }$$ is significantly larger than ν, in agreement with Huse and Fisher³⁰. As a consequence, the product ξ × ∣π/2 − q∣ goes to zero at the critical point as shown in Fig. 5c.

Fig. 5

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₃/Ω = 1.35 and V₃/Ω = 3.5, respectively. Inside the p = 4 phase, the correlation length is fitted with a power law with critical exponent $$ {\nu }^{\prime }$$. In the disordered phase, the correlation length is fitted either with a power law with critical exponent ν (a, d) or with the KT form $$ \xi \propto \exp \left(C/\sqrt{g_{KT}-g}\right)$$ (g), where g is the coordinate along the cut. The wave vector q is fitted with a power law with exponent $$ \overline{\beta }$$ (dotted lines). Wave vectors q are defined within the error bars ±πξ/N²; and ξ × ∣π/2 − q∣ is defined up to ±πξ²/N². 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.

The next cut at V₃/Ω = 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 $$ \overline{\beta }$$ 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 $$ {\nu }^{\prime }$$, 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.5 in Fig. 5g–i. The numerically extracted critical exponent $$ {\nu }^{\prime }$$ 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⁻³. 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³⁰. 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 self-dual versions of the three-state chiral Potts model^36–38.

Fig. 6

Behavior of the second derivative of the energy close to the transition line.

Effective critical exponent α across a the Ashkin–Teller point and b, c the chiral transition across oblique cuts perpendicular to the critical line. d Second derivative of the energy per site with respect to Δ/Ω for V₃/Ω = 2 around the Pokrovsky–Talapov transition. The results are extracted from the ground-state energy of a chain with N = 1201 sites (blue circles) and N = 2101 sites (red squares), and from the difference between the two (black diamonds). The gray area indicates the expected value of α for a critical exponent ν ≈ 0.78 ± 0.02.

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 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.

Details about the algorithm

The size of the Hilbert space for a model with two-site Rydberg blockade can be calculated using a recursive relation $$ \mathcal{H}\left(N\right)=\mathcal{H}\left(N-1\right)+\mathcal{H}\left(N-3\right)$$, with the first three elements of the sequence $$ \mathcal{H}\left(1\right)=2$$, $$ \mathcal{H}\left(2\right)=3$$, and $$ \mathcal{H}\left(3\right)=4$$. So the growth of the Hilbert space with the system size $$ \mathcal{H}\left(N\right)\propto 1.466^N$$ is much slower than $$ \mathcal{H}\left(N\right)\propto 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 hard-boson model with r = 1 can be rigorously mapped onto a quantum dimer model on a two-leg ladder³⁹ 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.

Fig. 7

Rigorous mapping onto a model that preserves the block-diagonal structure of tensors.

a Local Hilbert space of the original model $$ ∣l_i⟩$$. The open (filled) circle stands for an empty (occupied) site. b Rigorous mapping onto a model with a local Hilbert space spanned over three consecutive hard bosons that consist of four states sketched in c. The index of the new site corresponds to the index of the middle site. d Fusion graph for the recursive construction of the left environment; for the right environment the direction of the arrows should be inverted. e Example of the label assignment in MPS representation on two consecutive tensors (green circles) written for the selected state.

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⟩$$. 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 $$ {\tilde{n}}_i$$ with the only non-zero element $$ {\tilde{n}}_i\left(3,3\right)=1$$. The term V₃n_i−1n_i+2 can be written in the new Hilbert space as a nearest-neighbor interaction $$ V_3{\tilde{p}}_i{\tilde{q}}_{i+1}$$ where the only non-zero matrix elements of the operators $$ \tilde{p}$$ and $$ \tilde{q}$$ are given by $$ \tilde{p}\left(4,4\right)=1$$ and $$ \tilde{q}\left(2,2\right)=1$$. Finally the constrained flip term

With these definitions, the matrix product operator in the bulk takes the following simple form:

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 single-site 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 out-going 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.

Calabrese–Cardy formula

According to Calabrese and Cardy⁴⁰ the entanglement entropy in finite-size chain with open boundary conditions scales with the block size l as

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:

Fig. 8

Identification of the Ashkin–Teller parameter λ from the conformal tower of states.

a Excitation spectrum of the Ashkin–Teller model as a function of λ for N = 60 and A−A boundary conditions (black). Reg line indicate the point where the Ashkin–Teller spectrum resembles the structure presented in b. b Excitation spectrum of the hard-boson model with N = 201 sites (red). c Conformal towers of states for hard-boson (red), Ashkin–Teller at λ = 0.57 (green), and four-state Potts as Ashkin–Teller with λ = 1 (blue). The tower is plotted with respect to the lowest excitation energy.

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₃/Ω = 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:

Fig. 9

Extraction of the correlation length exponent ν and comparison with Ashkin–Teller model.

a Inverse of the correlation length along the commensurate line. b Critical exponent ν as a function of the Ashkin–Teller asymmetry parameter λ. The dark blue line shows the exact result of refs. ^28,29. The open red circle is the numerical result of this work ν ≃ 0.777 (result for N = 3001 shown in panel a, finite-size error do not exceed ±0.02) and λ ≃ 0.57 ± 0.01 as shown in Fig. 8.

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²e/d(Δ/Ω)². 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₂₁₀₁ − E₁₂₀₁ to suppress the edge effects and get a better estimate for the bulk energy per site as (E₂₁₀₁ − E₁₂₀₁)/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{\left(\mathrm{\Delta }/\mathrm{\Omega }\right)}^2$$ with respect to $$ \log \mid \mathrm{\Delta }-{\mathrm{\Delta }}_c\mid $$. 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³⁰. The results of Fig. 6d obtained for V₃/Ω = 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⁴³

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.