Observation of Hilbert space fragmentation and fractonic excitations in 2D

Experimental details

Here we briefly describe the initial-state preparation common to all measurements. Experiments were performed in a single plane of a vertical one-dimensional optical lattice. For the in-plane lattice, we used the folded lattice described in ref. ³⁹. As the in-plane lattice is subject to disorder and harmonic confinement, we used a digital micromirror device to shape the horizontal on-site potential, allowing us to achieve approximately homogeneous trapping depths and tunnelling energies throughout the system. Using a second digital micromirror device, we additionally projected a tapered, rectangular box in the centre of this corrected system, to achieve reliable loading and high filling in a central area of about 15 × 15 lattice sites.

Starting from these Mott insulators, to prepare the initial states of interest, we then performed local addressing over the entire area^41,42, whereas the data analysis was performed in a smaller region of interest (ROI) of either 8 × 8 or 10 × 10 lattice sites at the centre of the system. In addition, working with larger systems than the size of the ROI minimizes the influence of finite-size and boundary effects. With this preparation sequence, we achieved a filling of 0.88(2) per site on the addressed sites and a filling of 0.04(2) on the non-addressed sites in the ROI. These values were averaged over all datasets and initial configurations.

Magnetic-field gradient calibration

The potential tilt in our experiments was realized by global magnetic fields, which allowed us to induce the most homogeneous gradients. We calibrated the magnitude and the orientation of the magnetic gradient using spatially resolved microwave spectroscopy on the magnetic-field-sensitive transition between the \(| F=1,{m}_{F}=-\,1\rangle \) and the \(| F=2,{m}_{F}=-\,2\rangle \) hyperfine ground states.

To this end, we prepared a large Mott insulator with all atoms in the \(| F=1,{m}_{F}=-\,1| \rangle \) state. We then adiabatically ramped up the magnetic field to its target configuration and performed narrow microwave sweeps at variable centre frequencies. As a consequence, atoms were addressed resonantly within a narrow stripe subjected to the same magnetic-field strength and flipped into the \(| F=2,{m}_{F}=-\,2\rangle \) state, which were then removed before imaging. We fitted a two-dimensional Gaussian to these stripes of missing atoms, which allowed us to map the field strength and gradient orientation versus their position (Extended Data Fig. 1a).

To be able to continuously vary the applied gradient strength, we used a combination of coils: a single vertical gradient coil and a pair of vertical offset coils in Helmholtz configuration with reversed field polarity to realize a quadrupole field near the plane of the atoms. For the initial calibration, we worked with a fixed gradient coil setpoint and tuned the vertical offset and additional in-plane offset fields such that the magnetic zero point was at the location of the atoms; subsequently, we shifted the zero point by a fixed amount using the in-plane offset fields, resulting in an in-plane gradient at the correct angle. We then proceeded to calibrate the gradient strength for various gradient coil setpoints as described above; for technical reasons, we tuned the gradient coil instead of the offset coils. We interpolated between the calibrated values by fitting them with the function

where r is the displacement of the magnetic field zero to the atoms, g₀ is the maximal gradient strength, B₀ describes background fields and ΔB is the change of the setpoint of the gradient coil. As we changed ΔB by only a few per cent, we can assume g₀(B) = constant, which is also supported by the fact that the fit function describes the data well, as shown in Extended Data Fig. 1b.

On the basis of this curve, we can then rescale the x axis in Extended Data Fig. 2 and obtain an absolute value for the gradient strength.

Hubbard parameters

To extract the Hubbard parameters of our folded optical lattice³⁹, we made use of two methods. First, we performed amplitude modulation spectroscopy to calibrate the lattice depth. The results were then compared with a band-structure calculation to obtain the values for the on-site interaction U and the tunnelling energy J. Here we found U/J = 21(2) with U = h × 275(5) Hz and J = h × 13(1) Hz. The error bars arise from the uncertainty of the lattice-depth calibration itself as well as the slightly anisotropic hopping along the two lattice axes³⁹. Second, we can independently calibrate the Hubbard parameters using the quench dynamics of isolated dimers (see Extended Data Fig. 3 and below). As a result, we extracted τ = ħ/J = 10.0(3) ms, equivalent to J = h × 16.0(5) Hz. Comparing again with our band-structure calculation, this corresponds to U/J = 17(1) with U = h × 260(5) Hz. We attribute the deviations between these two calibrations to day-to-day drifts of the lattice beam alignment over the entire data-taking period.

For data evaluation, we used τ = 11 ms for all datasets, motivated by the long data-taking period of several days for a given dataset. Theory calculations (see below) were performed for U/J = 18, which was chosen as an intermediate value between the two calibrations.

Tuning the gradient to resonance

For the presented studies, it is important that the applied gradient matches the on-site interaction, that is, Δ = U. We benchmarked the resonance location by measuring the dimer imbalance as a function of the gradient strength for various tunnelling times, as illustrated in Extended Data Fig. 2. Here we expected the strongest decay of the dimer imbalance, as defined in the main text, when the resonance condition is fulfilled. For smaller gradients, we expected a slower drop in imbalance, whereas for much stronger gradients, we expected all processes to be off-resonant and no dynamics to occur at all, leading to high imbalance even at later times.

Our experimental results match the described expectation qualitatively. To confirm that we were not accidentally probing at a time where the imbalance shows any Δ-dependent oscillations, we probed for multiple fixed evolution times (up to t/τ = 40), observing consistent behaviour for all of the chosen evolution times. The resonance width is inherently limited by the finite tunnelling bandwidth and residual potential disorder. Our chosen operation point was located at the centre of the resonance and showed the strongest decay, as marked by the vertical dashed line in Extended Data Fig. 2. On the basis of our gradient calibration presented above and in Extended Data Fig. 1b, this point corresponds to a value of Δ = h × 238(3) Hz.

Comparing with our independent band-structure calculation, we found a qualitative agreement within 15% to the value of U for both calibration methods of the Hubbard parameters described above. In particular, U changes only very slowly with the lattice depth and varies by less than J for our calibrations. As such, this gradient setpoint remains valid throughout all measurements.

Data analysis

All data, unless specified differently, were analysed as explained in the following: we calculated the quantity of interest (imbalance, Fourier components, diagonal sums) on the individual experimental shots, then averaged over these results to obtain the data shown in the figures. For the reference-subtracted defect occupations (Fig. 4b, middle, and Extended Data Fig. 4b), we subtracted the densities averaged over all shots. To calculate the imbalances in Fig. 5c,d, we chose the boundary of the respective ROIs such that atoms close to the interface boundary that could be part of either the chequerboard or the dimer were counted as belonging to the dimer part of the system. As such, we obtained the same number of atoms for both halves of the system and the imbalance can, in principle, reach its typical limits of ±1; this explains the perhaps unintuitive shape of the ROIs shown in Fig. 5b.

Fourier analysis

To analyse the Fourier components of the average densities, we calculated the discrete Fourier transform according to

with a_n,m the average densities at site n, m for an ROI of size N × M and \({k}_{x}=\frac{2{\rm{\pi }}j}{N},{k}_{y}=\frac{2{\rm{\pi }}l}{M}\). Here the index j runs from \(-\left(\frac{N-1}{2}\right),\ldots ,0,\ldots ,\frac{N-1}{2}\) for odd N and \(\lceil \frac{N-1}{2}\rceil ,\ldots ,0,\ldots ,\frac{N}{2}\) for even N and analogously for l. The value at (k_x, k_y) = (0, 0) is just the sum of the signal; it contains no additional relevant information and is thus neglected (white rectangles in the insets of Fig. 3).

It is noted that the discrete Fourier transform obeys point reflection symmetry, that is, F(k) = F(−k). Therefore, in the main text, we plot only the parts of the momentum space (k_x, k_y) that contain non-redundant information.

Isolated dimer dynamics

To further understand and investigate the decay of the dimer pattern on a microscopic level, we prepared isolated dimers and tracked their evolution after a sudden quench. We isolated the dimers by adding empty columns between the atom pairs, as illustrated in the inset of Extended Data Fig. 3. For this configuration, the dimers were, including only first-order processes, completely decoupled from one another, allowing us to study the formation of the horizontally oriented dimers described in Fig. 3a. The change in orientation can be understood intuitively. Starting from a dimer, the upper atom can tunnel onto the neighbouring site by forming a doublon, as illustrated in the middle inset of Extended Data Fig. 3. From there, the atoms can either rearrange into the original dimer or into the flipped dimer, which is energetically degenerate to the original dimer configuration.

Extended Data Fig. 3 shows the time evolution of the isolated dimers. Here we plot the populations of the three possible states: the vertical dimer, the doublon and the horizontal dimer. Although the dimer states can be detected unambiguously, we assigned the doublon if all three sites were empty. To correct, on average, for cases where no atoms were initially prepared, we subtracted the value obtained analogously from a reference measurement tracking the initial-state preparation. We observed a clear oscillation between the two cases of vertically and horizontally oriented dimers, which quickly dephases owing to residual potential disorder. We compared the measured data with a numerical simulation of a three-state model given by

where δ_i, δ_j, δ_k describe the disorder strength between adjacent sites. For the calculation, we sampled δ_i, δ_j, δ_k from a normal distribution around zero and averaged over N = 100 such realizations. The additional factor of \(\sqrt{2}\) for the hopping has to be taken into account owing to the bosonic enhancement characteristic for indistinguishable bosons. We then fitted the calculations to the measured occupation of the vertical and horizontal dimers to generate the solid lines in Fig. 3b. Here we allowed for the disorder strength, the difference U − Δ, an overall amplitude (which respects normalization) as well as the timescale as free fit parameters. The initial time offset was kept fixed at zero. It is noted that the doublon occupation was not included in the fits, instead the solid line in Extended Data Fig. 3 is given by the model expectation using the fit values obtained from fitting the two other curves. We observed good agreement between the doublon occupation as obtained from our measured data and the numerical model using the fit parameters for the two other curves, validating our method of extracting the doublon occupation.

From the fit, we extracted the standard deviation of the disorder distribution σ = 1.2(1) × J, a deviation from resonance of U − Δ = 0.0(3) × J and a timescale of τ = 10.0(3) ms. The latter can serve as a secondary way to calibrate the Hubbard parameters of our system (see above).

Negative defect and additional analysis

Here we present our measurements on the negative defect and describe the data presented in Fig. 4 and Extended Data Fig. 4 in more detail. We also present an alternative way of evaluating the data for the positive defect and directly compare the spreading of the defect holes for both the negative and the positive defects.

The spreading of the defects can be observed directly in the average occupations (Fig. 4b and Extended Data Fig. 4b, leftmost column), through a reduced contrast of the (background) chequerboard on sites accessible to the defect atoms. This is owing to the following processes. First, the defect atoms can move to initially empty sites of the chequerboard, thereby increasing the average density on these sites. The defect atoms can also move onto initially occupied sites of the chequerboard, where we then observe a reduced average density owing to parity projection. Finally, nearby atoms from the background chequerboard can become mobile owing to the presence of the defect and move onto the site occupied by the defect atom, thus reducing the average density on their original sites as well as on the site of the defect atom owing to parity projection. For the negative defect in particular, the motion of the hole can be observed by an increase of the average occupation on its initial site (see also Extended Data Fig. 5b in the following), and a simultaneous decrease of the average occupation on the neighbouring sites on its equipotential line. By contrast, the hole site for the positive defect remains unoccupied. These effects are highlighted by subtracting the occupation of a chequerboard state without deterministically created defects. Initially empty sites accessible to the defect atoms will feature a positive, reference-subtracted value, whereas initially occupied sites accessible to the defect atoms will show negative values. The latter, as explained above, is due to either parity projection, atoms becoming mobile owing to the defect or, for the case of the negative defect, also the spreading of the defect hole. When comparing with theory, we observed good agreement, especially for the negative defect (Extended Data Fig. 4b,c). For the simulations, we did not include any experimental imperfections such as disorder and initial-state preparation fidelities. We further quantified the directional spreading of the defects by summing along the diagonals of the reference-subtracted occupations. When summing parallel to the equipotential lines, the occupation is only different from zero on the diagonals on which the defect atom and hole were initially placed. The growth by one additional diagonal for times t/τ > 0 can be explained by the above-mentioned processes, that is, the defect’s influence on the neighbouring atoms. As an additional characterization of the positive defect, we also studied the spreading on the zigzag-shaped equipotential line (Extended Data Fig. 5a, inset) instead of summing the reference-subtracted signal along the ROI diagonals. The result of this analysis is shown in Extended Data Fig. 5a. Here we again observe that the spread occurs along only one direction, as the immobile hole prevents the spread in the opposite direction. The latter is expected, as the hole can only move by second-order processes⁸. This is also evidenced by the density on the site of the hole remaining nearly unchanged.

Looking at this further, by comparing the increase of the densities on the sites initially occupied by the defect holes, we can also clearly observe the difference between the positive and the negative defects (Extended Data Fig. 5b). For the positive defect, the density increased only slightly, whereas for the negative defect we observed an immediate, fast increase, as here the hole is mobile in first order. Specifically, the hole of the negative defect can move in processes where neighbouring particles located on the equipotential line above the defect atom hop onto the defect atom and then to the site of the hole (Extended Data Fig. 4a, bottom right). The hole’s motion is restricted to its initial equipotential line. As for all other measurements on the spreading of defects on top of the chequerboard background, we attribute deviations from the theoretical expectations to disorder in the system and imperfect initial-state preparation, that is, the presence of additional, non-deterministic defects.

Numerical methods for defects

The underlying physics in the Bose–Hubbard model is described by an effective Hamiltonian derived in ref. ⁸, which features HSF. In Extended Data Fig. 6, we compare the time evolution of the positive defect under this effective Hamiltonian with the time evolution of the original Bose–Hubbard model. We show the parity-projected on-site occupation and have subtracted a perfect chequerboard state (at t/τ = 0, that is, without time evolution) to better highlight the differences. It is noted that in Fig. 4 and Extended Data Figs. 4 and 5a, we instead subtract the theory calculations with a time-evolved version of the chequerboard for better comparison with the experimental data. In contrast to the effective model, the background chequerboard state is not completely frozen under time evolution with the Bose–Hubbard model. Nevertheless, this additional dynamics of the background does not strongly influence the dynamics of the mobile defect compared with the effective model. Therefore, we conclude that the underlying physics of the Bose–Hubbard model in the chosen limits are well captured by the effective model featuring HSF.

Numerical methods for convergence

The data were obtained using tensor-network methods and exact diagonalization. All data were calculated by matrix-product-operator time evolution using the TeNPy package^46,47, except for the time evolution with the effective model in Extended Data Fig. 6, which was performed with exact diagonalization. In Extended Data Fig. 7, convergence in the bond dimension and in the Trotter step is studied. In Extended Data Fig. 7a, the evolution of the imbalance for the chequerboard and dimer states shows perfect overlap for bond dimensions χ = 256 and χ = 300. For both curves, a Trotter step size of dt = 0.001 was used. In Extended Data Fig. 7b, the imbalance is compared for Trotter step sizes of dt = 0.001 and dt = 0.0005 at a fixed bond dimension of χ = 256.

Numerical methods for imbalance

In Extended Data Fig. 7c, we compare the imbalance of the perfect case to the time evolution under imperfect conditions, similar to those of the experiment. For the latter, we have included deviations of all relevant quantities away from optimum, fidelities for state preparation and an additional random on-site potential (see the caption of Extended Data Fig. 7 for details). Each time step is averaged over N_av,dimer ∈ [29, 100], N_av,squares ∈ [17, 100], N_av,CHB ∈ [10, 100] different preparations. We find that the effect of state-dependent dynamics is still clearly visible also for experimental conditions. For the dimer state, the impact of experimental conditions is the strongest, which we attribute to the highest sensitivity to imperfect state preparation. In the case of the dimer state, all atoms have only one nearest neighbour. Removing this neighbour directly leads to a decrease in mobility and can induce frozen particles. By contrast, the squares state does not suffer from this effect on the same level. Each atom has three nearest neighbours, and therefore one missing neighbour does not lead to frozen sites.