Supersymmetry and multicriticality in a ladder of constrained fermions

Supersymmetric lattice models of constrained fermions are known to feature exotic phenomena such as superfrustration, with an extensive degeneracy of ground states, the nature of which is however generally unknown. Here we address this issue by considering a superfrustrated model, which we deform from the supersymetric point. By numerically studying its two-parameter phase diagram, we reveal a rich phenomenology. The vicinity of the supersymmetric point features period-4 and period-5 density waves which are connected by a floating phase (incommensurate Luttinger liquid) with smoothly varying density. The supersymmetric point emerges as a multicritical point between these three phases. Inside the period-4 phase we report a valence-bond solid type ground state that persists up to the supersymmetric point. Our numerical data for transitions out of density-wave phases are consistent with the Pokrovsky-Talapov universality class. Furthermore, our analysis unveiled a period-3 phase with a boundary determined by a competition between single and two-particle instabilities accompanied by a doubling of the wavevector of the density profiles along a line in the phase diagram.


Abstract
Supersymmetric lattice models of constrained fermions are known to feature exotic phenomena such as superfrustration, with an extensive degeneracy of ground states, the nature of which is however generally unknown. Here we address this issue by considering a superfrustrated model, which we deform from the supersymetric point. By numerically studying its two-parameter phase diagram, we reveal a rich phenomenology. The vicinity of the supersymmetric point features period-4 and period-5 density waves which are connected by a floating phase (incommensurate Luttinger liquid) with smoothly varying density. The supersymmetric point emerges as a multicritical point between these three phases. Inside the period-4 phase we report a valence-bond solid type ground state that persists up to the supersymmetric point. Our numerical data for transitions out of density-wave phases are consistent with the Pokrovsky-Talapov universality class. Furthermore, our analysis unveiled a period-3 phase with a boundary determined by a competition between single and two-particle instabilities accompanied by a doubling of the wavevector of the density profiles along a line in the phase diagram. 1 Introduction: superfrustration and multicriticality

Contents
Generally speaking, the physics of a quantum many-body system is determined by a competition between various terms, in particular the kinetic and interaction terms, of its Hamiltonian. For lattice models, endowing the kinetic term with additional density dependent constraints leads to so-called kinetically constrained Hamiltonians, for instance, the East [1][2][3] and PXP [4][5][6][7] models. Focusing on bosonic or fermionic models, such constraints turn out to have dramatic consequences for their physics, resulting for instance in topological ordering [8] or non-thermalizing behaviour following a quantum quench [9][10][11][12][13]. Among lattice models with kinetic constraints, supersymmetric models of spin-less fermions play a special role due to an enhanced mathematical control offered by the supersymmetry. Despite this fact, these models come with many outstanding questions. In particular, when going beyond one spatial dimension, it is not only the dynamics, but even the nature of the ground states which remains essentially unexplored. In this work, we consider such a model with a simplest geometry beyond a strictly 1D chain -a zig-zag ladder -and study its ground state phase diagram as we now describe.
Supersymmetric lattice models It has long been known that (space-time) supersymmetry in a quantum field theory (QFT) leads to special features in the physics described by such QFT, and to an enhanced mathematical control. An example of the latter is the Witten index [14] for theories with a complex (N = 2) supersymmetry, which guarantees the existence of zero-energy ground states, without the need of diagonalizing the Hamiltonian. The papers [15][16][17] traced the N = 2 supersymmetry in critical (conformal) or massive QFT's in 1+1 dimensions to an exact supersymmetry in associated 1D lattice models of spin-less fermions. These supersymmetric lattice models, dubbed the M k models, have a characteristic constraint, allowing at the most k nearest neighbour sites to be occupied.
Superfrustration Among the M k models, the M 1 model is particularly interesting -it features spin-less fermions with nearest neighbour exclusion principle, forbidding the simultaneous occupation of nearest neighbour sites. Consequently, this allows the M 1 model to be formulated on any graph G and there exist proposals how to engineer it in 1D systems using Rydberg atoms exploiting the blockade mechanism to implement the kinetic constraint [18][19][20]. Evaluating the Witten index for such models beyond 1D led to a surprise: many choices of G, such as ladders and 2D grids with underlying triangular, hexagonal or kagome lattice geometries, have a Witten index that is extensive in the system size N (number of vertices of G) [21][22][23]. This implies that the number of ground states grows exponentially with N , and that the ground state entropy per site is finite. This phenomenon was dubbed superfrustration in [22]. The ground state counting problem for a variety of choices of G has been addressed, with varying degree of success. For the 2D square lattice with toroidal boundary conditions (BC) precise results have been obtained, establishing that, depending on the choice of toroidal BC, the ground state entropy is at the most sub-extensive, scaling with the linear dimension of the system [24]. For the triangular lattice many results were obtained (bounds on ground state fermion densities [25,26], ground states on ladders [27] and finite patches [28]) but the ground state counting problem remains, to the best of our knowledge, still unsolved.
Multicriticality The massive degeneracy of ground states (which typically come with a range of fermion densities) indicates that supersymmetric points are highly singular points in the ground state phase diagrams of these lattice models. To study this phenomenon, we focus on a relatively simple case, which is the M 1 model on a zig-zag ladder. This model, while simple enough to allow powerful numerics, does display superfrustration, with ground states coming in a range of fermion densities 1/5 ≤ N f /N ≤ 1/4. Here and in what follows, N f is the fermion number. In this paper we study how this model is embedded in a phase diagram set by a Hamiltonian H, given in eq. (4), with parameters t, U , V 3 and V 4 , and reducing to the M 1 model for U = −t, V 3 = t, V 4 = t. We explore the vicinity of the supersymmetric point, which turns out to be a multicritical point connecting both gapped and floating phases.
Methods Our main tool has been numerical simulations performed with the state-of-the-art density matrix renormalization group (DMRG) algorithm [29][30][31][32] operating directly within the constrained Hilbert space [33,34]. The explicit implementation of the nearest-neighbor blockade on a zig-zag ladder allows us to reach convergence for critical systems with up to N = 1201 sites, keeping up to 1500 states (bond dimension of local tensors) and truncating singular values below 10 −9 . Further details of the algorithm will be discussed in Appendix A. Additionally, we supplement the numerics by analytic arguments that are possible in special regions of the phase diagram.
Specifically, considering a parameter space specified by the two interaction strenghts V 3 and V 4 , cf. Sec. 2, we probe a phase diagram described in Sec. 3, where we identify a floating phase (Sec. 3.2) and period-3,4 and 5 density wave phases and analyze the transitions out of them (Sec. 3.3-3.5). We also discuss in Sec. 3.3 how the boundaries of the period-3 phase can be estimated from single-and two-particle instabilities. Furthermore, we provide analytical explanation of the observed valence-bond type state in period-4 and of the density profiles in period-5 phases in Sec. 3.4 and Sec. 3.5 respectively. Sec. 3.6 explores the effect of changing the chemical potential U . We conclude in Sec. 4. Figure 1: Sketch of the system governed by the Hamiltonian Eq. (4). The hopping t and the interaction terms V 3 , V 4 are highlighted as solid blue and dashed red lines respectively. The orange region indicate the radius of the blockade where double occupancy is excluded, n i (1 − n i ) = n i n i+1 = n i n i+2 = 0, ∀i.

The model
We start by recalling the construction of a supersymmetric model for constrained fermions on a lattice. N = 2 supersymmetry is explicitly realized in the following Hamiltonian built out of two fermionic generators Q + and Q − Supercharges Q + and Q − are constructed from constrained fermions, i.e. the fermions with nearest-neighbor blockade on a selected lattice graph: where P i = <i,j> (1 − n j ) is a projector onto a constrained Hilbert space, with n j = c † j c j the local occupation operator and < i, j > denoting nearest neighbours on the lattice. On a generic graph the N = 2 supersymmetric Hamiltonian can therefore be rewritten as: In this paper we investigate the emergence of the supersymmetric point on a zig-zag ladder. Our many-body Hamiltonian acts on a restricted Hilbert space with n i (1 − n i ) = n i n i+1 = n i n i+2 = 0 and is given by The first two terms describe constrained nearest-neighbor hopping, the third term plays the role of a chemical potential and controls the density, while the last two terms describe the repulsion beyond the blockade. A sketch of the lattice geometry and the Hamiltonian terms can be found in Fig.1. Supersymmetry emerges when −U/t = V 3 /t = V 4 /t = 1, and (4) reduces to (3). Without loss of generality we fix t = 1. In addition we fix U = −1 and explore the vicinity of the supersymmetric point in the remaining two-dimensional parameter space of coupling constants V 3 and V 4 . In quasi-1D systems (chains or ladders) with open boundary conditions (OBC) the fermionic nature of the particles does not manifest itself due to local constraints. So fermionic operators in the Hamiltonian eq. (4) can be replaced by hard-boson operators with nearest and next-nearest neighbor blockade. We also note that the Hamiltonian (4) preserves the total number of particles N f = i n i .
The supersymmetric point U = −1, V 3 = 1, V 4 = 1 is accompanied with massive degeneracy of zero-energy ground states. For periodic boundary conditions (PBC), the supersymmetric Hamiltonian is H SUSY,PBC = H + N , with N the number of lattice sites. For OBC, supersymmetry requires the following boundary terms As a simple example, consider N = 5. For PBC, the Witten index is For N f = 1 and in the single particle basis {|1 , . . . , |5 }, where |i denotes a fermion at site i, the Hamiltonian is given by The translationally invariant state 1 √ 5 (1, 1, 1, 1, 1), with energy E = 5, is the superpartner of the state with zero particles. The supersymmetric groundstates can be chosen to be the eigenstates of the translation operator T with eigenvalues t l = e 2πil/5 with l = 1, . . . , 4.
For N = 5, OBC, the Witten index is We thus expect at least one supersymmetric ground state with an odd N f . As for (7), for N f = 1 and in the single particle basis {|1 , . . . , |5 }, the Hamiltonian is given by The coefficients of the unique supersymmetric ground state (1, 1, −4, 1, 1). In [35,36] an exact expression was found for the generating function P N (z) = Tr GS z N f of the ground state multiplicity for the M 1 model (with OBC) on the N -site zig-zag ladder. Unlike in the Witten index Eqs. (6), (8), where the trace is evaluated over the whole Hilbert space, here the trace is taken only over the subspace spanned by the ground states. It turns out that the generating function satisfies the recursion P N (z) = zP N −4 + zP N −5 and a general formula can be found and reads We note that P 5 (z) = z, i.e. there is a single ground state for N = 5 and OBC, in agreement with the explicit example discussed above, cf. Eqs. (8) and (9). The formula (10) reveals a massive (extensive) degeneracy of E = 0 ground states at the supersymmetric point, at densities in the interval 1/5 ≤ N f /N ≤ 1/4. Their total number grows like 1.167 N [35,36]. This clearly raises the question of the phase diagram in the vicinity of the supersymmetric point, which should be such that ground states at densities in the given interval all meet at a highly singular point in phase space.

Mapping to spins
Before analyzing the phase diagram we present a reformulation, first explored in [35], of the model (with OBC) in terms of unconstrained spin-1 2 degrees of freedom. Allowed configurations of the zig-zag ladder model are sequences of 0 and 1 with each 1 accompanied by at least two 0 on both sides. Now add '0' at ficticious sites i = 0 and i = N + 1 and then substitute 010 →↑ and then 0 →↓ for the remaining 0. This gives a sequence of N ↑ = N f up spins and We remark that spin Hamiltonians with tri-linear coupling similar to those in H SUSY,OBC spin have been proposed in the context of superconducting circuits [37] and recently realized experimentally [38]. It would be thus interesting to investigate whether one could engineer and quantum simulate the specific Hamiltonian (11) in these systems.

Phase diagram 3.1 Overview
We explore the vicinity of the supersymmetric point by tuning the coupling constants V 3 and V 4 of the Hamiltonian Eq. (4) (with t = 1, U = −1). Our main results are summarized in the phase diagram in Fig. 2. It contains three gapped density-wave phases with periodicity three, four and five; and a critical Luttinger liquid phase with the density not frozen to a specific value but continuously tuned by two coupling constants. The latter results into incommensurate quasi-long-range order, known in the literature as a floating phase [39]. The supersymmetric point emerges as a multicritical point between the gapped period-4 and period-5 phases and the critical floating phase with the density range 1/5 < N f /N < 1/4. In the period-3 phase every third site is occupied by a fermion while the hopping is completely suppressed due to local constraints. Using a notation (n i , . . . , n i+p ) for a repeated pattern of period p of particle densities, the resulting three ground states are classical of the form (1,0,0) (and (0,1,0) and (0,0,1) respectively) and are decoupled in the sense they are invariant under the action not only of the full Hamiltonian Eq. (4) but also of its individual terms. They also span (in the thermodynamic limit) the N f /N = 1/3 sector and are thus decoupled from each other and the rest of the Hilbert space. An example of the local density profile in the period-3 phase is presented in Fig. 3(a). In the period-5 phase the density is fixed to 1/5. Here, every fifth site is occupied, however, by contrast to the period-3 phase, quantum fluctuations are not suppressed. As a result the density profile is different from (1, 0, 0, 0, 0), as shown in Fig. 3(c). The density of the period-4 phase is fixed to 1/4. The density profile is very different from the previous two cases: single particles resonate between two nearestneighbor sites followed by two empty sites (see Fig. 3(b)). Based on the numerical results, we have phenomenologically established that the valence-bond solid type pattern (0.5, 0.5, 0, 0)

Floating phase
In the thermodynamic limit the density changes continuously inside the floating phase. We extract the density by averaging the local density over a number of sites as The interval [i min , i max ] over which we average always lies between two local maxima, as indicated by the red arrows in Fig. 3(d)-(f). This way, even if the wave-vector q is close to a commensurate value (which is in particular relevant in the vicinity of the periodic densitywave phases, cf. also Figs. 4,7,8,11), we obtain meaningful results. To reduce the edge effects we start with the maxima located at a distance of 20-80 sites from the edges. In Fig. 4 we demonstrate how the density changes across the floating phase for two selected cuts across the phase diagram. In Fig. 4(a) the cut at V 4 = 0.8 goes through the period-4 phase reflected in a pronounced plateau at 1/4 filling. Fig. 4(b) shows how the density changes between the period-3 and period-5 phases in the upper part of the phase diagram along V 4 = 2. By mapping the density profiles throughout the phase diagram we extract equal-density lines shown in the phase diagram in Fig.2. Surprisingly, part of the floating phase with the density range 1/5 <n < 1/4 collapses into the supersymmetric point in the "wedge" delimited by the period-4 and period-5 phases, cf. Fig. 2, and re-emerges on the other side of the wedge, within the same density range. The observed collapse leads to an extensive degeneracy of the ground states within the specified density range, resulting in the superfrustration at the supersymmetric point.
The floating phase is a critical phase in the Luttinger liquid universality class. We check this by extracting the central charge from the scaling of the entanglement entropy with the block size in open systems. According to conformal field theory (CFT), the entanglement entropy scales as [40] S N (l) = c 6 ln d(l) + s 1 + log g, where d(l) is the conformal distance d(l) = 2N π sin πl N , l is a size of the sub-block (1 l N ) and s 1 and g are non-universal. Friedel oscillations of the local density cause significant oscillations in the entanglement entropy profile. Nevertheless, when the system size is sufficiently large to accommodate many oscillations they can be safely averaged out, for example, by fitting the local maxima as shown in Fig. 5. The obtained values for the central charge at various points in the floating phase agree within 5% with the CFT prediction c = 1 for a Luttinger liquid.
We close this section by commenting on the criticality of ground states at the supersymmetric point. The issue was analysed in [23,35], by studying the response of the model with PBC to a twist in the boundary conditions. The conclusion was that, while many of the supersymmetric groundstates are gapped, some show critical behavior compatible with the number k unitary minimal model of N = 2 superconformal field theory, with central charge c k = 3k k+2 , with k even. Our present analysis is complementary to this approach, but it does suggest that, at the supersymmetric point, critical states exist at all (rational) fillings between 1/5 and 1/4. Such states can be obtained by following the floating phase state at given filling into the supersymmetric point.

Boundary of the period-3 phase
An estimate of the location of the boundary of the period-3 phase can be obtained via the following argument. Given the exclusion rule, the only possible density pattern at filling N f /N = 1/3 is of the form ..100100... Let us now imagine taking one particle out N f → N f −1. Furthermore, each of the defects d has a unit-strength hopping amplitude, giving a kinetic energy of ∆K = −2 per defect in the large-N limit. Assuming that the defects are far apart and independent gives as energy difference between the defect and the fully packed states where we used that in leading order each of the holes can hop with unit amplitude. This gives as phase boundary V 2d,d 3 = V 4 6 . The relative potential energy for a size-5 hole (triple defect 3d) is We note that the triple defect does not involve a direct hopping term, i.e. hopping between two 3d defects. More specifically, it is connected by the kinetic term to adjacent pairs of defects provide a first approximation to the boundary of the period-3 phase. Interestingly, we find that for V 4 sufficiently positive, the leading instability of the fully packed period-3 phase is to a configuration with N f = N 3 − 2 particles. Assuming that the 6 resulting defects group as 2d, 2d, 2d, we find a potential energy difference ∆V = 8−10V 3 . The nearest neighbour hopping terms in H of the bare particles, cf. Fig. 1, lead to unit-strength hopping of the double defects 2d, so that the estimated energy difference is, with the numerically found phase boundaries is excellent for sufficiently large system sizes.
In order to get further insight into the nature of the transition out of the period-3 phase we look at the energy level crossings in the immediate vicinity of the phase boundary. Because the total number of fermions is a conserved quantum number, on a finite size system we expect to see a set of explicit level crossings between the fully-packed period-3 state, and the states with one, two etc. particles less. We perform the simulation on systems with OBC and N = 31 and N = 49. OBC favors first and last sites to be occupied by the fermion, thus the fully packed period-3 state has in total N f = N −1 3 + 1 fermions. We look at the level crossing between this state and the states with one and two particles less. Based on the results presented in Fig. 6 we conclude that the single particle instability is a relevant excitation out of period-3 phase for V 4 1.2, while for V 4 1.2 the system is driven out of the period-3 phase by the double-fermion instability. In section 3.5 we find that this change of behaviour is connected to a doubling of the dominant density profile wave-vector q in the floating phase. In the vicinity of the period-3 phase, it changes from q = 2πN f /N for V 4 1.2 The commensurate-incommensurate transition between the ordered and the floating phases is expected to be in the Pokrovsky-Talapov universality class [41]. In Fig. 7 we show how density changes in the vicinity of the period-3 phase boundary. We fit our DMRG data with A|V 3 − V c 3 |β, whereβ = 1/2 is a Pokrovsky-Talapov critical exponent, and A and V c 3 are the two fitting parameters. Agreement with the field-theory prediction is remarkable for all cuts, below and above V 4 = 1.2, indicating that the double-fermion instability discussed above does not change the nature of the transition in the thermodynamic limit.

Period-4 phase
By contrast to the period-3 and period-5 phases, where (up to quantum fluctuations) a fermion occupies a single site followed then by 2 or 4 empty sites, the density profile of the period-4 3 + 1 (green), one fermion less (blue) and two fermions less (red). The interval of the single-particle instability δV 3 is indicated on panels (a)-(e) and is decreasing with N , indicating its finite-size nature. Upon increasing V 3 for V 4 = 1 (a)-(b) the ground state changes from the fully packed state to a state with one fermion less, then with two fermions less etc. However, starting from V 4 ≈ 1.2 and above the width of the single-fermion instability vanishes rapidly with the system size and for sufficiently large system the period-3 state changes directly into a state with two fermions less. For V 4 = 1.4 we do not observe a single-fermion instability as a ground state for N = 49 and larger; for V 4 = 2 a single-fermion instability does not show up for systems as small as N = 31. phase is fundamentally different: the fermion resonates between two nearest-neighbor sites followed by two nearly empty sites. This state is reminiscent of a valence-bond solid (VBS) state in quantum magnets.
Phenomenologically we have established that the pattern (0.5, 0.5, 0, 0) corresponds to an (approximate) ground state along a line V 3 = V 4 starting from V 3 , V 4 ≈ 0.8 and up to the supersymmetric point V 3 , V 4 = 1, cf. the dash-dotted line in Fig. 2. At the supersymmetric point, it is easily checked that state (assuming N = 4N f − 2 and OBC) is annihilated by both Q and Q † and is thus a E = 0 supersymmetric ground state of H SUSY,OBC [35]. We have observed that within the period-4 phase, essentially the same VBS state survives as ground state of H along the line V 3 = V 4 = a, 0.8 a < 1. This is easiest established for periodic boundary conditions (PBC) and N = 4N f , where the VBS state (together with three similar states obtained by translating (19) over 1, 2 or 3 sites) is an exact For OBC the situation is more subtle. Away from the edges the VBS pattern (0.5, 0.5, 0, 0) is again favored by the Hamiltonian with V 3 = V 4 = a, 0.8 a < 1, but there are corrections near the edges. In addition, for N odd (N = 4N f ± 1), the VBS pattern requires a breaking of the left-right reflection symmetry i ↔ N − i of the ground state (we indeed observe such symmetry breaking produced in our numerical DMRG procedure). For N even (N = 4N f −σ, σ = 0, 2), the finite size corrections to the PBC ground state can be understood in two ways: either one has to deform the Hamiltonian H at the boundaries to make the VBS state [such as (18)] an eigenstate or one has to modify the state itself to be an eigenstate of H. Assuming N = 4N f − 2 for specificity, in the former case if we modify the Hamiltonian into H a , defined as H a = H + a(2n 1 + n 2 + n N −1 + 2n N ), (20) the VBS state eq. (18) is once again an exact eigenstate, with E = −N + (N f + 2)(a − 1), and it is a ground state for a < 1 close enough to a = 1. In the latter case, considering H instead, the VBS state is modified near the edges. For example, for N = 10, N f = 3, and V 3 = V 4 = 0.9, the ground state has densities recovering thus the pattern ≈ (0.5, 0.5, 0, 0) in the bulk. As shown in Fig. 8(b) the exact line terminates at V 3 = V 4 ≈ 0.8 with the Pokrovsky-Talapov transition into the floating phase with density N f /N > 1/4. The density scales towards the transition with the Pokrovsky-Talapov critical exponentβ = 1/2. The location of the critical point is affected by noticeable finite-size effects. The transition across the remaining two sides of the period-4 phase is into the floating phase with lower density N f /N < 1/4. In Fig.8(a) we show how the density scales approaching the period-4 phase from above (larger V 4 ). The results suggest that the transition remains in the Pokrovsky-Talapov universality class.

Period-5 phase
The density profile in the period-5 phase has a characteristic pattern with peaks every 5 sites and small but non-zero values in between. Below we argue that the change in the local  density curve of the four quasi-unoccupied sites from concave to convex (see Fig. 9) is due to the dominant density wave-vector q changing its value from 2π/5 to 4π/5. The observed doubling of the wave-vector can be understood from the following analysis of a finite size system with PBC, where we add a symmetry breaking term that induces the density profile.
As an example, consider the ladder with PBC, 15 sites, 3 particles. Starting at the supersymmetric point V 3 = 1, V 4 = 1, we observe 7 degenerate E = 0 ground states, with eigenvalues t l of the translation operator given by In each of these states the density profile is flat at value 1/5. We then break the symmetry (both the supersymmetry and the translational invariance) by adding a term H = − (n 1 + n 6 + n 11 ).
Due to the strict ground state degeneracy, any 0 < 1 immediately breaks the symmetry and leads to a ground state with characteristic period-5 density profile. This same pattern arises in the setting with OBC for large enough system sizes. The qualitative form of the fluctuations in the density profile after breaking the symmetry can be understood as follows. The perturbation away from the supersymmetric point breaks the ground state degeneracy. For example, it turns out that V 3 = 1.2, V 4 = 2 gives a pair of ground states with l = ±3, while V 3 = 1.6, V 4 = 1.6 has two ground states with l = ±6. Turning on leads to a ground state that is a superposition of states with all five momenta satisfying t 5 l = 1, but we observe from the numerical solution that the convexity of the density profile follows the pattern at = 0.
For V 3 = 1.2, V 4 = 2, among the states with l = ±3, the symmetry breaking term favors the following combination of these states which also provides the leading contribution to the density profile with q-vector of 4π/5. For V 3 = 1.6, V 4 = 1.6, the analogous analysis with l = ±6 leads to a density pattern with q-vector 8π/5, which is equivalent to 2π/5. This explains the change of the behaviour of the local profiles between the dominant density peaks observed in Fig. 9 from concave to convex. Scanning the density profiles in the period-5 phase, we find that the q-vector is always 2π/5 or 4π/5 depending on V 3 and V 4 . This pattern extends to the rest of the phase diagram, cf. Fig. 10: for low enough values of V 4 the q-vector equals 2πn with n = N f /N the fermion density, but for higher values of V 4 we find double that value, q = 4πn. Interestingly, the line separating the two behaviours connects to the period-3 phase at the cusp point (V 3 = 1/5, V 4 = 6/5) where the leading instability out of the period-3 phase changes from 1-particle to a 2-particle instability, as discussed in section 3.3.
In analogy with the commensurate-incommensurate transitions out of the period-3 and period-4 phases, the critical line between the floating phase and the period-5 phase is also expected to be in the Pokrovsky-Talapov universality class. Fig. 11 depicts how the density in the floating phase scales towards the transition. For V 4 = 1.6 we show in Fig.11(a) that the density approaches its commensurate value 1/5 with critical exponent consistent with the field-theory predictionβ = 1/2. We extract an "effective" critical exponent by includinḡ β in the set of fitting parameters. The obtained valueβ ≈ 0.62 agrees within 25% with the Pokrovsky-Talapov value. For V 3 = 1.6 presented in Fig.11(b) the computed critical exponentβ ≈ 0.93 differs significantly from the theoretical expectation. Likely, this is because the presented data points are still too far from the transition to resolve the correct critical behavior. This is supported by the fact that the density in the immediate vicinity of the transition shows fast decrease down to 1/5 with an increasing system size, consistent with the shift of the critical value V c 4 towards the smaller values and followed by the decrease of an effective critical exponentβ. Unfortunately, accurate resolution of the Pokrovsky-Talapov transition in these two cases would require an access to much larger system sizes inside the floating phase which is currently beyond the limitation of our computing capacity.

Changing the chemical potential
In order to make the problem tractable, so far we have been focusing at a particular cut of the phase space of the Hamiltonian (4) at a chemical potential U = −1. It is thus interesting to investigate the ground state properties where also the chemical potential is varied. To this end Fig. 12 shows the situation along two cuts in the V 3 − V 4 plane, namely for V 4 = 1 (Fig. 12a) and V 4 = 2V 3 − 1 (Fig. 12b). Both cuts are also indicated as dash-dotted lines in Fig. 10, where the former crosses the supersymmetric point from and to a floating phase, while the latter from a crystalline period-4 to a period-5 phase. It is apparent from the results in Fig. 12 that the multicritical nature of the supersymmetric point remains preserved in other planes in the parameter space, namely that the period-4 and period-5 phases remain always separated by a floating phase with varying particle density.

Conclusion
The model studied in the paper is a prototypical example of a model displaying superfrustration: an exponentially large degeneracy of supersymmetric ground states or, equivalently, a ground state entropy that is extensive in the system size [21][22][23]. By studying a twoparameter phase diagram in the vicinity of the supersymmetric point we find that it emerges as a multicritical point connecting period-4 and period-5 many-body groundstates. Superfrustration arises through the collapse of a floating phase, with intermediate densities, into the supersymmetric point.
It would be extremely interesting to generalize our conclusion beyond the models that are supersymmetric by construction. This can be achieved, for example, by weakening the blockade. Both gapped phases as well as the floating phase are expected to survive in the soft-blockade regime and there is no immediate argument that would prevent the multicritical point to appear. It would be thus interesting to investigate how the properties such as the number of the ground states or the supersymmetric nature of this multicritical point get modified when softening the blockade. In light of our results, it might be wise to look further into the surroundings of the superconformal points in other systems, for instance, in chains of SU(2) k anyons [42]. In particular, it would be interesting to study what will happen to the superconformal critical point separating Z 2 phase from the Haldane phase in the presence of perturbations that breaks translation of topological symmetry. It would also be interesting to investigate whether various Z n phases can be connected via a floating phase.

A Constrained DMRG
Explicit next-nearest-neighbor blockade implies that the total dimension of the Hilbert space grows with the length of the chain as H(N ) ∝ 1.466 N [34] which is much slower than H(N ) ∝ 2 N for an unconstrained model of spin-less fermions. In order to fully profit from a restricted Hilbert space of constrained fermions, the next-nearest-neighbor blockade has been explicitly encoded into DMRG. Here we briefly recap the main set of implementation details.
First, we perform 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. 13(b). Empty circles indicate original local degrees of freedom -empty l 1 or occupied l 2 . Green and blue circles corresponds to the left and right normalized MPS tensors. By spanning local degrees of freedom over three sites (dotted lines) the local Hilbert space increases from two states |l i sketched in Fig. 13(a) to four states |h i sketched in Fig. 13(c). States with two and more occupied sites are forbidden by the constraint. The new local Hilbert space |h i corresponds to the physical leg of the local tensor shown as a vertical solid line in Fig. 13(b). From the Fig. 13(b) it is obvious that any of the two consecutive tensors have a pair of common (or shared) sites. Three possible states of these two sites can be used as quantum labels (00), (10) and (01) for auxiliary legs that naturally create a block-diagonal structure of local tensors. The latter drastically reduces the computational cost of simulations. In the bulk, quantum labels of the left environment are changing according to the fusion graph shown in Fig. 13(d). Fusion graph can be interpreted in the following way: to the chain that ends with two empty sites on the right and therefore labeled by (00) on the right side of the chain one can add either an occupied site that would corresponds to the local state |h 2 and change the label of the chain to (01), or one can add another empty site that corresponds to the local Hilbert space |h 1 and does not change the label (loop on the top of the fusion graph); if the chain ends with an occupied site it is labeled by (01) and only an empty site with |h 3 can be added on its right side, in this case the label changes to (10); finally, if the chain ends with an occupied site followed by an empty site it is labeled with (10), due to two-site blockade one can add only empty site next to it that corresponds to |h 4 and the new label will be (00). For the right environment, i.e. when we add sites on the left side of the chain, the direction of the arrows in the fusion graph has to be reversed. An example of the label assignment on a pair of consecutive tensors is provided in Fig. 13(e). At the next step, we have to write down the Hamiltonian Eq. (4) in terms of local matrix product operators (see e.g. Ref. [32] for the review on the MPO approach). For this one has to rewrite the spin-less fermion model given by Eq.4 in terms of the new local variables |h i . For example, the local occupation number operator n i 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 interaction terms 2V 3 n i n i+3 and V 4 n i n i+4 are transformed into two-and three-body terms 2V 3f1f2 and V 4f1f3f2 correspondingly, where each of thef -matrix has only one non-zero entry:f 1 (4, 4) = 1,f 2 (2, 2) = 1, andf 3 (1, 1) = 1.
With these definitions, the MPO in the bulk takes the following form: where the dots mark zero entries of the tensor. Note that all entries are 4 × 4 matrices, so resulting MPO is a rank-4 tensor with dimensions 4 × 4 × 14 × 14. 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. Further details on constrained DMRG can be found in Refs. [33,34].