Quantum coarsening and collective dynamics on a programmable simulator

Experimental platform

A detailed description of our experimental platform is given in refs. ^11,60. All measurements are realized using a two-dimensional programmable quantum simulator based on Rydberg atom arrays. Single ⁸⁷Rb atoms are stochastically loaded into optical tweezers shaped by a spatial light modulator (SLM), and then rearranged into defect-free patterns using a pair of crossed acousto-optic deflectors. Both sets of tweezers use 852-nm-wavelength light. The atoms are then laser-cooled and optically pumped to the \(| 5\rmS_1/2,F=2,m_F=-\,2\rangle \) state, which we denote as \(| g\rangle \) in the main text. A pair of counterpropagating lasers at 420-nm and 1,013-nm wavelengths couple \(| g\rangle \) to the highly excited Rydberg state \(| r\rangle =| 70\rmS_1/2,J=\frac12,m_J=-\frac12\rangle \) through a two-photon transition through the \(6\rmP_3/2\) orbital, blue-detuned by approximately 2π × 2.4 GHz with respect to the \(5\rmS_1/2\to 6\rmP_3/2\) transition.

We use 16 × 16 lattices of atoms, and maintain an R_b/a ratio of 1.12–1.15, such that only nearest neighbours lie within the blockade radius. Owing to the rapid fall-off of the van der Waals term, only nearest- and next-nearest-neighbour interactions meaningfully contribute to our observations. The data shown in the main text are taken with two-photon Rabi frequencies of either Ω/2π = 3.8 MHz, Ω/2π = 6.0 MHz or Ω/2π = 3.1 MHz, with corresponding lattice spacings a = 6.8 μm, a = 6.45 μm or a = 7.15 μm and R_b/a = 1.15, R_b/a = 1.12 or R_b/a = 1.13. The 1,013-nm and 420-nm single-photon Rabi frequencies are approximately balanced (Ω_1,013 ≈ Ω₄₂₀). The experiment time T, sweep parameters and lattice spacings are all rescaled for different Ω, such that Ω/Δ, R_b/a and the total phase accumulated ΩT are constant when comparing different experimental configurations.

To observe coarsening, the sweep rate through the critical point must fall within a certain range. A sweep rate that is too slow would create a low-energy state, and consequently, the coarsening dynamics may be too slow to measure. In contrast, a very fast sweep or an instantaneous quench could inject too much energy and bring the system out of the ordered phase (which persists up to a finite energy density). Here, all measurements use linear sweeps with sweep rates of \(\frac\varDelta \,/\,\varOmega \varOmega T\approx \frac12\rm\pi \times 3.0\), which, we find, falls within the desired range.

For the global sweeps, we apply post-selection based on the success of the rearrangement protocol, selecting shots with ≤4 defects. This threshold retains an average of 93% of shots over all end-detuning values presented in Figs. 1 and 2a,b. In addition, we post-select on measurement results to discard runs in which we suspect large-scale errors have occurred. Owing to the high energy of the Rydberg blockade, a large number of blockade violations are extremely unlikely to be naturally generated by the coherent dynamics of equation (1). Furthermore, as our projective measurement cannot differentiate \(| r\rangle \) occupation from loss induced by other mechanisms, we attribute the presence of a large number of apparent blockade violations to manifestations of unwanted noise processes, such as blackbody-induced avalanche decays⁶¹. We therefore discard runs where the longest chain of consecutive atoms in a single row or column detected in state \(| r\rangle \) has a length of more than four sites. Small-scale exact diagonalization (ED) numerical simulations support that the probability of reaching such states from purely unitary dynamics is negligibly small. Imposing this post-selection threshold retains 92% of the data. Imposing both rearrangement and avalanche post-selection, we retain 86% of the global-sweep shots across all end-detunings presented. We note that the post-selection is most significant for the longest times and largest detunings presented (t = 1.16 μs, Δ/Ω = 4 in Extended Data Fig. 2d), as larger number of shots are corrupted through avalanche decays and we only retain 63% of the shots at the longest time step.

Local control

To enable individual single-site addressing of atoms with a local light shift, we use an SLM (Hamamatsu LCOS-SLM X15213-02) to generate optical tweezers in arbitrary spatial patterns with a beam waist of 1 μm, ensuring robustness to atomic position fluctuations (Extended Data Fig. 1). The wavelength we choose to operate at, 784 nm, achieves a measured differential a.c. Stark shift between the 5S_1/2 and 5P_3/2 states of 12.2(3) MHz with about  160 μW per spot, but a negligible scattering rate (about 35 Hz) (the scaling of the light shift with laser amplitude is shown in Extended Data Fig. 1c). The light is linearly polarized to minimize vector light shifts on the ground-state hyperfine manifold. We further measure the shift on the \(| g\rangle \to | r\rangle \) transition and find that the light shift is well approximated by the differential 5S_1/2 → 5P_3/2 light shift as δ₀ = −2π × 12(2) MHz at the same power per spot.

The phase holograms for the SLM are generated using the phase-fixed weighted Gerchberg–Saxton algorithm, taking into account the desired position and relative intensity of the local light-shift pattern⁶². We first generate a local addressing pattern that closely matches the positioning of the atomic tweezer array; however, perfect matching of the two arrays is computationally expensive as it requires an extremely high sampling rate of the image plane of the local addressing pattern. To overcome this computational barrier, after creating an initial local addressing pattern, we align it to the atom positions by transforming the phase hologram. By stretching, rotating and applying tilts and defocus, we can match the two patterns with feedback on the atom signal. The latter three can be easily controlled using Zernike polynomials, whereas the stretching and rotation require more care to preserve the intensity homogeneity of the desired pattern. We find that naive rescaling or rotating of the hologram results in unwanted distortion of the intensity pattern, attributed to software interpolation when working with a pixelated hologram. This is mitigated by applying the computational corrections in the image plane. We take the Fourier transform of the hologram, convolve the intensity profile with a two-dimensional Gaussian to broaden each spot over several pixels (to minimize effects of interpolation), and then apply the rotation and stretching. Lastly, we apply an inverse Fourier transform back to the Fourier plane and use the resultant phase hologram for the SLM. Using this procedure, we first coarsely align the individual addressing pattern to the tweezers on a camera, and then precisely align the two using a spin-echo measurement of the light shift (Extended Data Fig. 1b) to optimize the alignment parameters such that the intensity is maximized at the atom sites. Good alignment is also crucial to prevent atom loss from turning on a misaligned potential. Finally, we correct the tweezer intensities as required using the fitted light shifts to feedback on the target weights in the hologram generation.

Examples of states prepared using such tweezer profiles, where the detuning is used to strongly pin atoms to the ground state, are shown in Extended Data Fig. 1d. At the boundary between different antiferromagnetic orders and the edges of the array, the mean-field repulsive interaction strength decreases for sites with fewer Rydberg neighbours; we therefore weight the local detuning strength inverse-proportionally to the number of neighbours. It is noted that when arbitrary weighting is used, the total power remains constant (number of addressed sites × 2π × 12 MHz), but the power is redistributed in the tweezers accordingly. Nevertheless, particularly at large Δ/Ω, neighbouring Rydberg excitations start to be energetically favoured (antiblockaded) along the domain boundaries. Excluding such edge effects, the preparation probability of preparing the single-atom ground state on the pinned sites is 93–95% (Extended Data Fig. 1e).

For other realizations of single-site addressing using light shifts in atom arrays, see, for example, refs. ^63,64.

Theoretical background of coarsening dynamics

In this section, we summarize the theoretical details of the different kinds of coarsening processes that govern the dynamics of the system as long-range order is formed. Although our focus will be on the Rydberg atom array, to begin, let us consider the generic situation of a system driven through a continuous QPT by tuning some parameter of the Hamiltonian, g, linearly with time. Without loss of generality, we assume that the quantum critical point is located at g = 0 and the zero of time is set such that g(t) = t/τ; hence, the system crosses the quantum critical point at t = 0. For the specific case of the neutral atom array considered in this work, the time-dependent parameter g can be defined as g(t) = (Δ(t) − Δ_c)/Ω.

As the system approaches the quantum critical point, its relaxation time diverges and it necessarily falls out of equilibrium. However, when it does so depends on the velocity of the linear ramp, \(\mathopg\limits^.(t)=1/\tau \). The quantum KZM posits that the time at which the system’s evolution ceases to be adiabatic is t = −t_KZ with \(t_\rmKZ\approx t_(\tau /t_)^\nu z/(\nu z+1)\), where ν is the correlation length exponent, z is the dynamical critical exponent and t₀ is some microscopic timescale. Thereafter, as the system cannot dynamically respond fast enough to the changing parameter of the Hamiltonian, it remains ‘frozen’ through a so-called impulse regime until a later time t = +t_KZ, when it unfreezes on the other side of the QPT. During this impulse regime, the KZM presumes that the system’s correlation length remains the same as when it initially froze: \(\xi _\rmKZ\approx l_(\tau /t_)^\nu /(\nu z+1)\), where l₀ is some microscopic length scale. As a consequence, in this picture, the correlation length in the ordered phase is also set by ξ_KZ with no subsequent dynamics.

However, the non-equilibrium correlation length of the system, ξ(t), can and does grow in both the impulse regime and the ordered phase as the long-range correlations take time to develop. In the experiments described in the main text, this occurs via a two-step process. First, as the system passes through the quantum critical regime, it undergoes quantum critical coarsening, which is governed by the dynamical critical exponent z of the particular quantum critical point; for the (2+1)-dimensional Ising transition, z = 1. Then, as time progresses and the ramp continues, the system eventually enters the ordered phase. Here, once the growing non-equilibrium correlation length ξ(t) exceeds the equilibrium correlation length of the quantum ground state (which, recall, scales as ξ_q ≈ ∣g∣^−ν), the dynamics cross over to a regime of non-critical coarsening, for which

where ε is the many-body gap between the ground and the first excited states. The dynamical exponent z_d is dependent on the dimensionality and conservation laws of the system. For curvature-driven coarsening dynamics with a non-conserved scalar order parameter—as is indeed the case experimentally—z_d = 2 > z. A particular feature of non-critical coarsening worth emphasizing is the dependence of the dynamics on the distance to the quantum critical point encoded in equation (3). Specifically, the ground-state equilibrium correlation length scales as ξ_q ≈ ∣g∣^−ν and the gap ε ≈ ∣g∣^νz. Plugging in the exponents of the (2+1)-dimensional Ising QPT, ν = 0.629 and z = 1, along with z_d = 2, in equation (3), we find the growth law

This relation can be observed in Fig. 3, which studies the rate at which a locally introduced domain in the centre of the array shrinks. The area of such a domain decreases at a rate dr²/dt ≈ −ξdξ/dt, which scales as \((\varDelta -\varDelta _\rmc)^-0.629\), consistent with the behaviour observed in Fig. 3d.

For a ramp that continues indefinitely without stopping, the entire dynamical evolution of the correlation length can be described by a single universal scaling function encompassing the adiabatic, quantum critical coarsening and non-critical coarsening regimes^6,17,32:

where f(x) is some universal function, and ξ_KZ and t_KZ depend on the ramp rate τ as specified earlier. The scaling variable x delineates the three regimes discussed above as

More generally, if the ramp is stopped at a time t_s, the dynamical scaling form is altered to

for x ≤ x_s, one recovers the earlier scaling as \(\mathcalF(x,x_\rms)=f(x)\). If the ramp is stopped at x_s = t_s/t_KZ ≫ 1 in the non-critical coarsening regime, the behaviour of the universal scaling function for x > x_s describes the physics during the hold time and is given by

for some \(\mathcalO(1)\) constants \(\mathcalC > \mathcalC_\rms\). It is noted that because z < z_d, the coarsening speeds up as we stop earlier in the ordered phase, closer to the quantum critical point (which results in a lower x_s). This is indeed what we observe experimentally during the hold time following global sweeps across the phase transition, as shown in the inset of Fig. 1c. Intuitively, this is because a smaller Δ/Ω corresponds to a greater relative influence of critical coarsening, which is faster than non-critical coarsening.

In contrast, near the thermal phase boundary, the system can undergo an interval of classical critical coarsening, which is described by a growth law \(\xi (t)\approx t^1/\barz\) with a distinct dynamical exponent. For the two-dimensional classical Ising phase transition, \(\barz\approx 2.16 > z_\rmd\) (ref. ⁴³), so the growth of correlations through classical critical coarsening is slower than for non-critical coarsening. Correspondingly, the dynamics should decelerate as one approaches the classical critical point, in sharp contrast to the speed-up outlined above in the vicinity of a quantum critical point.

Structure factor and correlation length

To extract the structure factor and the correlation length⁶⁵, we first calculate the two-point connected correlation function \(G(\bfr_1,\bfr_2)\,=\)\(\langle \widetildeZ_\bfr_1\widetildeZ_\bfr_2\rangle -\langle \widetildeZ_\bfr_1\rangle \langle \widetildeZ_\bfr_2\rangle \) and then average over all pairs of points with identical displacements r:

We first derive the standard structure factor by computing the Fourier transform of G(r)

and then calculate the radially averaged structure factor

To extract a correlation length, we fit S(k) to:

and we factorize S₀ as S₀ = bξ²/π. This form of equation (12) is equivalent to assuming that the position–space correlations follow an exponential decay, \(G(r)\approx A\exp \left(-\fracr\xi \right)\) (ref. ¹¹), up to finite-size corrections. Although equilibrium considerations for an infinite-size system suggest that for small k, S(k) should obey the Ornstein–Zernike form \(S(k)\,\approx \)\(\fracS_1+\xi ^2k^2\) (ref. ⁶⁶), we empirically find that equation (12) better captures our observed non-equilibrium distributions.

For universal coarsening dynamics, we theoretically expect b to be constant. Although our data indeed show scaling collapse as in equation (12), we find that b varies during the dynamics, indicating the presence of an additional length scale(s). We observe that b is correlated with ξ and depends on Δ/Ω, as shown in Extended Data Fig. 3c. Additional length scales that may affect the dynamics include the finite system size, the finite width of the domain walls (which depends on Δ/Ω), spatial inhomogeneity in Δ and/or Ω, and length scales introduced by decoherence effects such as decay owing to the finite lifetime of the Rydberg state.

Specifically, the expected universal scaling regime is expected to hold for distances r and correlation lengths ξ such that l ≪ r and ξ ≪ L, where l is the width of a domain wall and L is the system size²⁴, indicating that finite-size effects probably have an important role in the present experiments. Hence, observing the theoretically expected universal coarsening behaviour in global quenches would probably require access to larger system sizes and correspondingly longer experiment times. Although recent experimental advances in neutral atom array platforms suggest that lattices more than an order of magnitude larger than the one presented in this paper are within reach⁶⁷, elsewhere we describe the use of local control to deterministically nucleate and study domain dynamics away from the system’s boundaries, allowing us to study universal properties of coarsening under present experimental conditions.

Analysis of domains in global sweeps

Using single-site-resolved detection, we can map out the domains in each snapshot. First, we calculate the local staggered magnetization. Each domain is then identified as a region of the array where the same ordering AF₁ or AF₂ is connected by nearest neighbours. We do not consider single spins of opposite order as a separate domain. For Fig. 2a,b, we therefore first identify and correct individual spin flips. These are identified as single atoms that are of the opposite order compared with all of their nearest and next-nearest neighbours. Only after we have identified single spin flips and corrected them to match their surrounding bulk order do we identify the domain boundaries. A domain’s area is defined as the total number of atoms comprising the domain. For the probability distribution of domain occupations presented in Fig. 2a, the frequency of each domain area is weighted by the area of that domain. We normalize the distribution by the sum of all area-weighted frequencies at each time step.

Classical energy analysis

To calculate the classical energy per single shot of the experiment, we first perform the single-spin-flip correction as described above. We then identify regions of the array that do not belong to either antiferromagnetic ordering by calculating a coarse-grained local staggered magnetization with a similar approach to previous studies¹¹. In this work, specifically, we calculate the convolution C_x,y of the Rydberg occupation n_x,y with the kernel \(W=\left(\beginarrayccc0 & 1 & 0\\ 1 & 0 & 1\\ 0 & 1 & 0\endarray\right)\) for each snapshot. The output values of C_x,y range from 0 to 4, where the extremal values correspond to atoms surrounded by nearest and next-nearest neighbours that all belong to the same AF orderings, as shown in Extended Data Fig. 4a. We consider an atom to be at a boundary if n(x, y) = 1(0) and C_x,y ≠ 0(4) (Extended Data Fig. 4b). In the raw array (not single-spin-flip corrected), we then compute the classical energy using equation (2) for each snapshot. The value of the interaction energy of nearest-neighbour sites (V_nn) is calculated from the lattice spacing (a) and V₀ as V_nn = V₀/a⁶. For the dataset presented in Fig. 2c, V_ij = V_nn/2π = 11.69 MHz for Ω/2π = 6 MHz. We also account for next-nearest-neighbour contributions. For each next-nearest-neighbour Rydberg atom per snapshot, an additional V_nnn/2π = 1.46 MHz is considered in the classical energy. The effects of longer-range interactions are negligible. By using the spin-flip-corrected lattice for domain identification and the uncorrected one for the subsequent energy calculation, single spin flips that occur contribute to the classical energy in identified domain walls and bulk orderings but not as separate domains. We exclude the layer of atoms closest to the edge of the array for all contributions to the classical energy. It is noted that for the classical energy calculation in Fig. 2c, we post-select such that the maximum number of directly adjacent Rydberg atoms can be no more than three (compared with four used throughout this work). Owing to the sensitivity of the boundary identification procedure used here to correlated decays, this post-selection is slightly stricter. In the data presented in Fig. 2c, as the data are at high end-detuning, we retain 69% of data owing to the avalanche post-selection and 84% owing to rearrangement post-selection. Overall, we retain 58% of data.

Analysis of locally prepared domains

We estimate the local radius of curvature R given in the equation ∂_tR² = −v_a as the radius, r, of the central domain in Fig. 3. The radius r is defined as the Manhattan distance, d_m, at which the radially averaged local staggered magnetization, m_s, crosses zero (Extended Data Fig. 4c). We consider Manhattan distances instead of Euclidian distances from the centre of the injected domains as the former is more representative of the nearest-neighbour interactions dominating the dynamics for very short times (on the order of one Rabi cycle). For long times, both measurements of distances in our lattice reveal a collapse to the linear form shown in Fig. 3c. For this analysis, we consider only the unpinned sublattice. The local state-preparation protocol prepares states where atoms that are locally detuned are prepared in \(| g\rangle \) with high probability for all values of Δ/Ω (Extended Data Fig. 1d). We observe that for Δ/Ω close to the QPT, there is a larger discrepancy of the local order parameter m_s between the pinned and unpinned sublattices. Therefore, when considering both sublattices, the radius is less clearly defined by a single point at which the order parameter crosses zero. We fit a linear relationship to extract dr²/dt at each Δ/Ω.

A similar procedure is followed for the analysis of the coarsening dynamics in the zigzag domain wall in Fig. 4. Here too, the mean value of the local order parameter m_x,y is calculated at each lattice site per time step. For Fig. 4c, the domain wall’s horizontal position is calculated as the point at which the linearly interpolated line between points crosses m_s = 0. In Extended Data Fig. 4d, we show the variation of the domain-wall position with hold time for two additional values of Δ/Ω. These data points reinforce the strong Δ/Ω dependence of the domain-wall velocity already seen in Fig. 3c,d. It is noted that for this analysis, we include atoms that were initially locally detuned as we are considering each row separately. We therefore see larger uncertainty in the domain walls’ positions for low Δ/Ω (Extended Data Fig. 4d).

Errors in Figs. 3c and 4c, and Extended Data Fig. 4d are calculated using bootstrapping. From the full set of experimental single snapshots of size N ≈ 600 shots, we sample N times with replacement and calculate the value of interest (r² or the horizontal domain-wall position x) on each sample. The plotted error bar is the standard deviation of the value of interest, calculated from 1,000 repetitions of the above procedure.

Numerical simulations of local domains

We simulate the dynamics of locally prepared domains using the time-dependent variational principle (TDVP)^46,68. We use a two-site variant of this algorithm, which allows the bond dimension to grow with the evolution time at the expense of forgoing strict energy conservation owing to the truncation step involved. In our calculations, we find that the energy is conserved to within 0.004% of that of the initial state up to the longest times simulated.

The initial states for the numerics are chosen to be a mean-field approximation of the experimentally prepared state. Specifically, we pin certain lattice sites to \(| g\rangle \), as specified by the target configuration, and the remaining ones are set to the vector on the Bloch sphere that minimizes the system’s mean-field energy for a given Δ/Ω (instead of the fully polarized \(| r\rangle \) state). The many-body evolution is simulated using a maximum bond dimension of χ = 1,200 with a time step Δt of 0.2 Ω⁻¹ (the dynamics are also consistent with those for a smaller Δt = 0.1 Ω⁻¹). The results thus obtained are showcased in Extended Data Fig. 5 and are found to be in good agreement with the experiments.

Amplitude/‘Higgs’ mode

Background and theory

To describe the observed amplitude mode, we consider the low-energy effective action that describes the transition between the disordered and antiferromagnetic phases. Its Lagrangian is the ϕ⁴ theory

where ϕ corresponds to the coarse-grained order parameter of the antiferromagnetic phase. Although the phase transition is described via the Wilson–Fisher fixed point, here we perform a simple mean-field treatment of equation (13) to capture the physics away from the immediate vicinity of the transition. In particular, the classical mean-field equation of motion for the order parameter’s expectation value is given by

which corresponds to a classical anharmonic oscillator. The stationary value of the order parameter is given by ϕ₀ = 0 in the disordered phase (q > 0) and \(\phi _=\pm \sqrt-q/\lambda \) on the ordered side of the transition (q < 0). Expanding equation (14) for small amplitudes, ϕ = ϕ₀ + δϕ, around the potential minima leads to harmonic oscillations of the order parameter with frequencies \(\omega (q > 0)=\sqrtq\) and \(\omega (q < 0)=\sqrt \).

Numerical simulations

We first investigate the low-energy spectrum of the two-dimensional Rydberg Hamiltonian at different values of Δ/Ω using the excited-state density-matrix renormalization group method, which iteratively finds the lowest-energy eigenstate that is orthogonal to previous lower-energy eigenstates. Our simulations take into account van der Waals interactions up to third-nearest neighbours on the square lattice. We identify the gapped paramagnetic and spontaneous-symmetry-breaking phases, and obtain the ground-state energy as well as the energy gaps ΔE₁ and ΔE₂ of the first two excited states above the ground state (in the spontaneous-symmetry-breaking phase, the first excited state that we identify is the symmetry-related ground state). In Extended Data Fig. 6a, we perform bond-dimension scaling for ΔE₁ and ΔE₂ on a 10 × 10 lattice with open and periodic boundary conditions up to bond dimensions χ = 512 and χ = 256, respectively. We note that although the energy gap in the paramagnetic (disordered) phase is robust to boundary conditions, in the spontaneous-symmetry-breaking (ordered) phase, we identify a distinct boundary energy gap (Extended Data Fig. 6a(i)) that is smaller than the bulk energy gap (Extended Data Fig. 6a(ii)). In the dynamics, the coupling of the order parameter to the boundary mode vanishes with increasing system size, and we thus consider the bulk gap extracted from periodic boundary conditions as the relevant frequency of the amplitude (Higgs) mode.

Furthermore, we perform density-matrix renormalization group calculations to obtain the initial state of the quench dynamics shown in Fig. 5, which is the ground state of the Hamiltonian with additional local detunings ∣δ_l∣/Ω = 0.7 that pin one sublattice of the chequerboard to the ground state. We evaluate the energy expectation value of this initial state with respect to the unpinned two-dimensional Rydberg Hamiltonian, and compare this with the energy of a \(\mathbbZ_2\) product state, confirming that the pinned state is indeed a low-energy state (Fig. 5b).

Dynamics from the low-energy pinned initial state lead to amplitude oscillations at frequencies matching the ground-state gap, which we simulate using TDVP^46,68 on a 10 × 10 lattice at a bond dimension χ = 256. We use time steps of 0.25 Ω⁻¹ and have verified that the resulting dynamics agrees with smaller time steps of 0.1 Ω⁻¹. As seen in Extended Data Fig. 6b, for small systems with open boundary conditions, the dynamics in the ordered phase (here, Δ/Ω = 2) are dominated by a slow mode with frequency matching the boundary gap shown in Extended Data Fig. 6a(i). On top of this slow mode, we see the presence of a mode with higher frequency matching the bulk gap shown in Extended Data Fig. 6a(ii). Therefore, in the following, we consider dynamics with periodic boundary conditions, which allow us to isolate the bulk mode.

In Extended Data Fig. 6c, we observe clear oscillations of the order parameter deep in either phase, which we fit to a damped harmonic oscillator \(\phi (t)\approx \phi _+A\cos (\omega t+\theta _)\rme^-\gamma t\) with frequency ω, damping γ, offset ϕ₀ and amplitude A. We show the dependence of these parameters on the ratio Δ/Ω in Extended Data Fig. 6d. Away from the phase transition, the oscillation frequencies overlap with the previously obtained ground-state energy gaps and are robust to system size. We further see that the damping and amplitude become larger towards the transition, where the offset acquires a non-zero value. Close to the transition, the TDVP dynamics fails to converge in bond dimension and fit to the damped harmonic oscillator’s functional form, as apparent in Extended Data Fig. 6c(ii). Moreover, even away from the critical point, the limited bond dimension does not capture a finite damping rate of the oscillations.

Dependence on local detuning strength

We experimentally investigate the dependence of the amplitude mode on the strength of the applied local detunings at two points, Δ/Ω = 0 and Δ/Ω = 1.1. These data are obtained by performing the state-preparation sequence illustrated in Fig. 3a for varying strengths of the applied local Stark shift δ. For all other uses of the local detunings in this work, the pinning is applied at a constant magnitude, δ₀ = −2π × 12(2) MHz, and the corresponding state preparation for various Δ/Ω is documented in Extended Data Fig. 1c,d. Here, this strength is varied to much lower values than for the saturated pinning of the ground state (0.04–0.1 δ₀) in the rest of the work (as indicated in Extended Data Fig. 9a,b) before significant differences in the resultant oscillations are observed.

We find that at both values of Δ/Ω considered, the amplitude of the oscillations is progressively reduced as ∣δ∣ is decreased. The frequency of the oscillations also decreases with ∣δ∣. The change in the oscillation frequency is more pronounced for oscillations near the critical point than in the disordered phase. The δ dependence of the oscillation frequency, as well as the increased sensitivity near the phase transition, are qualitatively consistent with the behaviour of an anharmonic oscillator, as predicted by the mean-field equation (14).

Amplitude mode in global sweeps

The ‘Higgs’-mode dynamics are also apparent in parallel to coarsening, as manifested in the oscillations of the magnetization n_i (Extended Data Fig. 8a). In addition, they can be clearly discerned by observing the dynamics of the total magnetization of the two-point correlation function in position space C(r), as shown in Extended Data Fig. 8b–d. The frequency of the oscillations of the correlation length closely support those extracted from the quench protocol described below (Extended Data Fig. 9) and the calculated ground-state energy gap (the global-sweep data are plotted in purple in Fig. 5). As mentioned earlier, the interplay of the amplitude mode with coarsening dynamics is generally unexplored theoretically. Therefore, we note further exploration of the two processes occurring in parallel, as a possible future extension of this work.

Quenches in the ordered phase

Although ‘Higgs’ oscillations in the ordered phase can be extracted by the state-preparation sequence through deterministic preparation with local detunings, as described in Fig. 3a, the amplitude of the oscillations is substantially reduced when compared with that inside the disordered phase (see Δ/Ω = 1.5 in Fig. 5). We therefore perform an alternative state-preparation sequence to extract the amplitude-mode frequencies deep in the ordered phase, as shown in Extended Data Fig. 9a. First, local detunings are applied in a chequerboard pattern as the global detuning Δ is swept from negative values to Δ/Ω = 3.3, a point far inside the ordered phase. We hold the global detuning constant while quenching off the site-dependent δ. At this point, we quench Δ to its final detuning value in the ordered phase, but closer to the phase transition. By way of this protocol, we observe, as with all other sweep protocols presented thus far, long-lived oscillations of the order parameter and correlation lengths as shown in Extended Data Fig. 9b,c. The extracted oscillation frequencies ω are in close agreement with the ground-state energy gap in the ordered phase. It is noted that in Fig. 5d, the points in red at Δ/Ω = 2.0 and Δ/Ω = 2.5 are measured through this quench protocol.

Frequency doubling

From the above-mentioned quenches to the ordered phase as well as in data from the local protocol, we extract the oscillations in both the order parameter and the correlation lengths. We find, as shown in Extended Data Fig. 9d, that in the ordered phase, these two frequencies are approximately equal, whereas in the disordered phase they vary by \(\omega _\xi /\omega _m_\rms\approx 2\). The changing relationship between the two observables can be understood by the following symmetry argument.

We begin by taking into account the dynamics of order-parameter fluctuations within a Gaussian approximation. Neglecting corrections to the effective mass due to fluctuations, the relevant equations of motion are^33,69:

where π_k(t) ≡ ∂_tϕ_k(t) and \(\mathcalD_\bfk,\bft^\phi \phi \equiv \langle \phi (-\bfk,\bft)\phi (\bfk,\bft)\rangle _c\). From the correlation function \(\mathcalD_\bfk,\bft^\phi \phi \), which corresponds to the structure factor discussed in the main text, one can extract the evolution of the correlation length.

Expanding ϕ = ϕ₀ + δϕ and \(\mathcalD_\bfk,\bft^\phi \phi =\mathcalD_\bfk^\phi \phi +\delta \mathcalD_\bfk,\bft^\phi \phi \), in the disordered phase, an eigenmode analysis of equations (15)–(17) yields a frequency spectrum \(2\sqrtk^2+q\). The smallest frequency, which is expected to set the correlation-length oscillations⁶⁹, is thus \(2\sqrtq\), that is, twice that of the order parameter. In contrast, in the ordered phase, equation (16) contains a term \(6\lambda \phi _\mathcalD_\bfk^\phi \phi \delta \phi (t)\), and the oscillation of the order parameter thus acts as a linear drive on the dynamics of two-point correlation functions. As such, the correlation length will oscillate at the corresponding frequency \(\sqrt \) of the order parameter.

Frequency ratio of oscillations

In Extended Data Fig. 10, we present the full dataset of amplitude-mode oscillations (also shown also in Fig. 5) as a function of the distance from the phase transition. As described above, Landau mean-field theory predicts the relationship between oscillation frequencies to be \(\omega (-| q| )/\omega (| q| )=\sqrt2\) (refs. ^13,47). However, beyond mean-field theory, the phase transition is described by the Wilson–Fisher fixed point, and the universal frequency ratio is shifted accordingly. Theoretical estimates based on both analytical and numerical methods yield a frequency ratio around ω(−∣q∣)/ω(∣q∣) ≈ 1.9 (refs. ^48,70,71,72). We find that both the experimental data and matrix-product-state simulations deviate from the mean-field prediction and are suggestive of a similarly higher ratio. However, very close to the critical point, ∣(Δ − Δ_c)/Ω∣ ≲ 0.3, we observe deviations from this theoretical ratio. Possible explanations (besides the limited bond dimension for the matrix-product-state data) include finite-size effects (resulting, for example, in a non-vanishing gap at the QPT point), possible errors in the QPT location, and the possibility that sufficiently close to the transition, the overdamped oscillations may no longer track the ground-state excitation gap. In future work, a detailed exploration of the region near the critical point through the amplitude mode could allow for higher precision tests of this universal ratio.