The model
We realize the topologically ordered Laughlin state on a quantum processor through constructing an HVA for its parent Hamiltonian defined by the following effective one-dimensional fermion chain model25,26 on a cylinder geometry (see Methods)
$$H={\sum }_{j}{\sum }_{k > m}{V}_{km}{c}_{j+m}^{{\dagger} }{c}_{j+k}^{{\dagger} }{c}_{j+k+m}{c}_{j},$$
(1)
where \({c}_{j}^{{\dagger} }\) and cj are the fermionic creation and annihilation operators corresponding to the single-particle orbitals under the Landau gauge. Physically, the index j specifies the x-coordinate of Gaussian-localized electron wave functions (Fig. 1a). The interaction matrix elements Vkm implement the Haldane-Trugman-Kivelson pseudopotential27,28, under which the ν = 1/3 Laughlin state (referred to as exact state throughout this work) is an exact ground state. This repulsive interaction decays at different rates for different interaction ranges (k + m) as the cylinder’s circumference Ly increases.
Fig. 1: Cylinder geometry and interaction truncation effect on Laughlin state.
The alternative text for this image may have been generated using AI.
a Schematic of cylinder geometries in Tao-Thouless (thin-cylinder) limit Ly → 0 and the isotropic geometry limit Lx ≈ Ly corresponding to Ly ≈ 10 in (b). The Gaussian peaks illustrate the localized orbitals of the lowest Landau level along the axial direction, with spacing \(2\pi {l}_{B}^{2}/{L}_{y}\) where lB is the magnetic length. Opacity of the Gaussian peaks represent local electron density. b Fidelity between the exact state and the ground state of the effective Hamiltonian for various truncation ranges of interactions (k + m ≤ 3, 4, and 5) in Eq. (1) for system with number of electrons Ne = 6, 7, and 8. The cylinder’s height Lx is determined through the constraint NΦ = LxLy/(2π) where NΦ is the number of flux quanta in the system and satisfies NΦ = 3Ne − 2 (see “Methods”). Lines are guide to the eye.
It is important to recognize that the exact state’s defining behaviors, such as incompressible quantum liquid correlations and long-range entanglement, are not universally captured by the ground state of Eq. (1) for arbitrary Ly. Its characteristics are hosted by the ground state of Eq. (1) only near the isotropic geometry limit when the cylinder’s circumference (Ly) matches its height (Lx)25. Strong deviations from it, such as the Tao-Thouless (TT) limit (Ly → 0), where the ground state becomes a charge-density-wave (CDW) state \(\left|{\Psi }_{{{\rm{CDW}}}}\right\rangle=\left|100100100…\right\rangle \) (Fig. 1(a)), and the squeezed cylinder limit (Ly → ∞), where the system is collapsed into a one-dimensional Luttinger liquid, lead to unfaithful description of exact state’s physical behavior.
Due to the two-body interactions in Eq. (1), the full Hamiltonian H contains \({{\mathcal{O}}}({N}^{3})\) terms for N orbitals, making variational ansatz based on the full Hamiltonian impractical for large system sizes. To address this, we develop an efficient and scalable protocol that constructs a HVA with an effective Hamiltonian Heff which retains only the dominant terms for correlated topological electronic systems (see Methods).
In this protocol, the terms in Heff are selected and validated following two criteria: (i) quantitative fidelity of wavefunction, and (ii) qualitative preservation of topology, entanglement, and symmetry. The first criteria is universal for quantum simulations of molecules and solids. The terms in Heff may be identified heuristically by their large ∣Vkm∣, which determines the term’s energy scale. Their validity can be further verified via ED within computationally viable regimes by comparing the wavefunction overlap and low-energy spectra of Heff and H. The second criteria is specific for the topologically ordered states. Qualitatively, we ensure the target state retains its defining properties—such as symmetry and topological order by verifying that Heff belongs to the same topological class as H, using topological invariants, entanglement entropy, or symmetry classifications.
Since FQH states are governed by short-range correlations, we expand Eq. (1) by interaction range (k + m) and evaluate the fidelity \({{\mathcal{F}}}\), defined as the wavefunction overlap between the exact state and the ground state of Heff consisting of truncated interactions as a comparative diagnostic across truncation ranges, rather than as an absolute threshold. This quantifies how well Heff captures the exact state’s key features. Figure 1b shows that from the TT limit to Ly < 7, all truncations regardless of the interaction range yield high fidelity. But as we approach the isotropic geometry regime Ly ≈ 10, the exact state’s strong correlation and long-range entanglement kicks in. As a result, \({{\mathcal{F}}}\) drops at significantly different rate depending on the truncations range. With only the lowest-order scattering (k + m≤3), \({{\mathcal{F}}}\) drops to 0.8 at Ly = 10 for system with number of electrons Ne = 6, whereas including longer-range interactions (k + m ≤ 4, 5) increases \({{\mathcal{F}}}\) to 0.95 and essentially 1.0, respectively.
Following the second criterion, we study how the interaction truncation range affects topology and entanglement. With only the lowest-order scattering (k + m≤3) included, the action of the effective Hamiltonian HTT on the CDW state \(\left|{\Psi }_{{{\rm{CDW}}}}\right\rangle \) forms a Krylov subspace \({{\mathcal{K}}}({H}_{{{\rm{TT}}}},\left|{\Psi }_{{{\rm{CDW}}}}\right\rangle )\). As an example of Hilbert space fragmentation29, this can be used to map FQH model, such as the Laughlin state’s parent Hamiltonian, under TT limit onto exactly solvable spin models30,31. This Krylov subspace \({{\mathcal{K}}}\) is significantly smaller than the full Hilbert space of a generic exact state. As a result, the second Rényi entanglement entropy \({S}_{A}^{(2)}=-ln{{\rm{Tr}}}{\rho }_{A}^{2}\) of the HTT ground state, computed for a subsystem A of the cylinder, rapidly saturates to a finite value as the subsystem boundary Ly increases toward the isotropic limit. This behavior signals a breakdown of area law scaling and the loss of the exact state’s correlation structure. In contrast, extending the truncation range to (k + m≤4) or higher restores the linear scaling of \({S}_{A}^{(2)}\) with Ly, recovering the expected area law behavior of a topological quantum liquid (See Supplementary Information).
Based on the quantitative criteria of fidelity and qualitative criteria of topology and entanglement, we choose k + m ≤ 4 as the truncation range of interactions in Heff. While incorporating longer-range interactions (k + m ≥ 5) can marginally improve fidelity, it does not qualitatively affect the topology or entanglement properties of the ground state. On the other hand, it significantly increases the complexity of the HVA circuit, pushing it beyond the capabilities of current NISQ devices. Thus, we conclude the minimal Heff for constructing the HVA for the ν = 1/3 Laughlin state includes the following interaction terms
$${H}_{{{\rm{eff}}}}= {\sum }_{j}\left[{V}_{10}{\widehat{n}}_{j}{\widehat{n}}_{j+1}+{V}_{20}{\widehat{n}}_{j}{\widehat{n}}_{j+2}+{V}_{30}{\widehat{n}}_{j}{\widehat{n}}_{j+3}\right.\\ + \left.({V}_{21}{c}_{j+1}^{{\dagger} }{c}_{j+2}^{{\dagger} }{c}_{j+3}{c}_{j}+{V}_{31}{c}_{j+1}^{{\dagger} }{c}_{j+3}^{{\dagger} }{c}_{j+4}{c}_{j}+\,{{\rm{H.c.}}})\right],$$
(2)
where \({\widehat{n}}_{j}={c}_{j}^{{\dagger} }{c}_{j}\) is the density operator. We note that at the interaction range k + m = 4, we retain only the off-diagonal scattering term V31 in Heff, which plays a crucial role in shaping the wavefunction structure and avoiding Hilbert space fragmentation. In contrast, V40, despite falling within the same interaction range, is a diagonal electrostatic term that primarily results in energy shifts without significantly influencing the wavefunction. To further reduce circuit depth, we exclude V40 from Heff (see Supplementary Information).
Quantum circuit for state preparation
With Heff identified based on our selection criteria, we construct the corresponding state preparation circuit in HVA fashion to simulate the Laughlin state on a quantum processor, with the expected HVA repetition p scaling linearly with the system size; the number of variational parameters per repetition being constant, so the total parameter count scales as \({{\mathcal{O}}}(p)\).
We interpret the HVA as a digitized adiabatic protocol generated by a local effective Hamiltonian32. Lieb-Robinson bounds on the spread of correlations under local dynamics imply an effective light cone with finite velocity33. We therefore expect that, for our Laughlin state HVA, the number of repetitions p must grow at least linearly with the system size in order to faithfully reproduce the long-range entanglement structure of the topological phase. As we show below, our ansatz also achieves a linear scaling of the total number of variational parameters by generalizing parameters in an HVA layer across the lattice. This avoids the quadratic or worse parameter growth that would result from assigning independent parameters to every microscopic term and aligns with previous work showing that constrained HVA remains expressive while improving trainability32,34,35.
The state preparation circuit \({|\psi (\{{\beta }_{j}\})\rangle }_{{{\rm{eff}}}}\), shown in Fig. 2, is given by the following unitaries
$${\widehat{U}}_{km}={\prod }_{j}\exp [-i{\beta }_{km}({c}_{j+m}^{{\dagger} }{c}_{j+k}^{{\dagger} }{c}_{j+k+m}{c}_{j}+\,{{\rm{H.c.}}})],$$
(3)
where βkm are variational parameters. The sum of indices are implicitly bound by the system size. The construction and optimization of \({|\psi (\{{\beta }_{j}\})\rangle }_{{{\rm{eff}}}}\) is guided by two fundamental principles. Firstly, we generalize the variational parameters βkm throughout the lattice, due to the similarity in mathematical structures at different j [Eq. (1)]. In practice, this means that all gates within the same unitary \({\widehat{U}}_{km}\) share the same parameter βkm, yielding a constrained HVA32,34 with one variational parameter per physical generator \({\widehat{U}}_{km}\) rather than one per microscopic term. As a result, each HVA repetition uses five independent parameters, independent of the system size. This dimensionality reduction of parameter space not only simplifies the variational optimization but also ensures the total number of parameters grows only through the number of HVA repetitions p, i.e., \({N}_{{{\rm{param}}}} \sim {{\mathcal{O}}}(p)\propto {N}_{e}\) (see Supplementary Information for explicit circuit gate and depth). Secondly, the squeezing rule in FQH36 requires \({\widehat{U}}_{21}\) as the first layer of the circuit which only contains terms with \(j=3n,n\in {\mathbb{Z}}\).
Fig. 2: Schematic N-qubit Hamiltonian variational ansatz circuit for preparing the ν = 1/3 Laughlin state.
The alternative text for this image may have been generated using AI.
The initial state is taken as the charge-density wave state \(\left|{\Psi }_{0}\right\rangle=\left|100100….1001\right\rangle \), where we use the Jordan-Wigner transformation in this work68. Commuting operators in \({\widehat{U}}_{km}\) are executed in parallel. We show the structure of \({\widehat{U}}_{20}\) layer as an example (see Supplementary Information for a full state preparation circuit at Ne = 6).
Using classical simulator (noiseless), we optimize βkm for the exact state in the isotropic geometry regime (see “Methods”), and demonstrated that the optimized parameters obtained with Ne = 6 can be transferred to larger systems as warm starts, assuming a fixed HVA repetition p. The optimized parameters βkm achieves \({{\mathcal{F}}}=0.93\) compared with the exact state, the ground state of the full Hamiltonian (1) obtained by ED at Ne = 6. Since the fidelity between the ground state of Heff and the exact state decays naturally with system size Ne (Fig. 1), we expect the fidelity between \({|\psi (\{{\beta }_{j}\})\rangle }_{{{\rm{eff}}}}\) and the exact state to follow the same trend when we transfer the optimized parameters to larger systems. Figure 3a shows the fidelity scales as expected for larger systems up to Ne = 10. Optimizing \({|\psi (\{{\beta }_{j}\})\rangle }_{{{\rm{eff}}}}\) with larger system size did not achieve higher fidelity (see Supplementary Information), further supporting parameter transferability and our construction’s resilience to barren plateau35. This smooth transferability suggests that parameters optimized on smaller systems provide high-quality warm starts for larger systems, reducing classical optimization costs and mitigating trainability issues when one subsequently increases the HVA repetition p with system size.
Fig. 3: Finite-depth scaling of fidelity and intensive quantities for the optimized protocol in the isotropic geometry regime.
The alternative text for this image may have been generated using AI.
a Fidelity between the state preparation circuit and ground state obtained by ED for system with number of particle Ne = 6–10. (Blue triangle) Fidelity between \({|\psi (\{{\beta }_{j}\})\rangle }_{{{\rm{eff}}}}\) and \(\left|{\Psi }_{{{\rm{eff}}}}\right\rangle \), ground state of Heff. (Red circle) Fidelity between \({|\psi (\{{\beta }_{j}\})\rangle }_{{{\rm{eff}}}}\) and the exact state \(|{\Psi }_{{{\rm{exact}}}}\rangle \). b Average deviation of local density δ〈nj〉. (c) Average deviation of two-point correlation function δ〈Cij〉. In (b, c), deviation of the quantity 〈x〉 is defined as \(\delta \langle x\rangle=| {\langle x\rangle }^{{\prime} }-{\langle x\rangle }_{{{\rm{exact}}}}| \), where 〈x〉exact is the exact state’s value and \({\langle x\rangle }^{{\prime} }\) corresponds to \({\left|\Psi \right\rangle }_{{{\rm{eff}}}}\) or \({|\psi (\{{\beta }_{j}\})\rangle }_{{{\rm{eff}}}}\). All error bars indicate the 16th and 84th percentiles. Lines are guide to the eye.
Notably, the average deviation of intensive quantities, such as the local density and two-point correlation between \({|\psi (\{{\beta }_{j}\})\rangle }_{{{\rm{eff}}}}\) and the exact state, remain constant with increasing system size (Fig. 3b, c). This observation strengthens the smooth transferability and suggests that for Heff considered here, reproducing local physics with high accuracy, does not require large prefactors in the linear depth scaling of the HVA. As such, our protocol can be extended sensibly to near-term quantum simulations of strongly correlated topological systems at scale.
Lastly, the Hamiltonian in Eq. (1) exhibits both particle number conservation \(\widehat{N}={\sum }_{j}{\widehat{n}}_{j}\) and center-of-mass coordinate conservation \(\widehat{K}={\sum }_{j}j{\widehat{n}}_{j}\,(\,{{\rm{mod}}}\,\,N)\). The unitaries \({\widehat{U}}_{km}\) composing our state preparation circuit naturally respect these symmetries, constraining the subspace of the variational search. Similarly, the final state \({|\psi (\{{\beta }_{j}\})\rangle }_{{{\rm{eff}}}}\) must transform identically under these symmetries as the initial state \(\left|{\Psi }_{0}\right\rangle \), enabling symmetry-verification protocols for robust error-mitigation37,38.
Edge and bulk density structure
We next proceed to prepare and probe the Laughlin state on quantum processors. A key question we sought to address was whether a deep quantum circuit, involving hundreds of two-qubit gates but only a few variational parameters, could successfully capture the physics of strongly correlated topological states on NISQ devices. While the cost of storing and manipulating many-body wavefunctions grows exponentially on classical hardware, this experiment, if successful, would be an important step toward scalable quantum simulations for materials-intrinsic topological order on near-term quantum processors. Given the depth of the circuit, i.e., 369 two-qubit gates for Ne = 6, we selected a trapped-ion quantum processor (IonQ’s 25-qubit Aria-1) for its relatively high two-qubit gate fidelity and low readout error rates, both of which are critical for mitigating noise and enabling effective post-selection strategies (see Methods).
One of the defining features of the quantum Hall states is the existence of chiral edge modes. On the cylinder geometry, the bulk-boundary correspondence39,40 guarantees the presence of chiral edge modes, which emerge from the bulk’s nontrivial topological order and appear as oscillatory deviations in the local density structure near the physical boundary41. We can directly probe this edge structure in the prepared state by measuring the local density operator \(\langle {n}_{j}\rangle=\langle {c}_{j}^{{\dagger} }{c}_{j}\rangle \) where \({n}_{j}=\frac{1}{2}(1-{Z}_{j})\) under Jordan-Wigner transformation.
In Fig. 4, we present the measured 〈nj〉 obtained by executing our state preparation circuit for Ne = 6 on Aria-1. Despite the limitation of current NISQ devices, the edge density structure is distinctly identified with an overdensity near the system boundaries (j = 0, 15) and subsequent oscillatory deviations of 〈nj〉 from the bulk filling fraction ν = 1/3. Away from the boundaries, the bulk region exhibits a relatively uniform density plateau, signaling the incompressibility and homogeneity nature of the topologically ordered Laughlin state. This spatial structure – a compressible, gapless edge surrounding an incompressible bulk – is an emblematic signature of FQH liquids.
Fig. 4: Probing edge and bulk density structure.
The alternative text for this image may have been generated using AI.
〈nj〉 is the observed electron occupation at site j, obtained by sampling 5000 shots on IonQ’s Aria-1 quantum computer with symmetry-verification postselection (PS) and debiasing error-mitigation (red triangle), which leads to a 10% selection rate. Error bars indicate 68% confidence intervals obtained by means of percentile bootstrap. These results are compared with noiseless simulation of state preparation circuit (orange square) and exact values obtained by ED (blue circle). Lines are guide to the eye.
The ability to resolve these edge structures relies critically on the symmetry-verification error mitigation that is naturally supported by our state preparation circuit. On the day of execution, Aria-1 reports a mean two-qubit gate fidelity of 98.5%. With approximately 300 two-qubit gates per qubit’s light-cone, a naive estimate implies a circuit fidelity of 1%, making error mitigation crucial to retrieve meaningful information from experiments on NISQ device. To address this challenge, we employ a combined error mitigation strategy: a custom symmetry-verification postselection protocol alongside IonQ’s debiasing mitigation scheme42. The postselection depends on the conservation of particle number and center-of-mass coordinate that are both respected by our state preparation circuit. Any measured bitstrings violating either of these two symmetries are deemed unphysical and thus discarded during postselection.
With IonQ’s debiasing mitigation alone, the result displays a systematic drift towards 〈nj〉 = 0.5, corresponding to the expectation value from a maximally mixed state, though the overall trend aligns qualitatively with the exact value obtained by ED. The application of symmetry-verification postselection significantly improves the fidelity of the results, eliminating the drift and confirming the observation of Laughlin state’s edge density structure (see Supplementary Information for debiasing only data and details on postselection).
Spatial correlation and topological entanglement entropy
After establishing the presence of edge modes, we turn to investigate the incompressible bulk region of the prepared Laughlin state. In the bulk region, the Laughlin state behaves as an interacting incompressible quantum liquid. This results in a uniform featureless bulk density but leaves nontrivial spatial fingerprints in the wavefunction. To investigate such spatial characteristics, we measure the two-point correlation function Cij = 〈ninj〉 − 〈ni〉〈nj〉 between site i and j. By construction, Cij is inversion-symmetric, that is, Cij = Cji and approaches 1 (−1) when the electron densities are correlated (anticorrelated).
With debiasing mitigation alone, we observe clear spatial signatures of anticorrelation in the first two off-diagonal elements of Cij, consistent with repulsive interactions (see Supplementary Information). After applying symmetry-verification postselection (Fig. 5a), we fully resolve the spatial correlation contrast of the correlated electron liquid. Additionally, long-wavelength density fluctuations are strongly suppressed as Cij converges rapidly to zero as ∣i − j∣ increases. The long-range correlation remains negligible in the bulk, except near the system’s boundaries where edge effects dominate.
Fig. 5: Spatial correlations and incompressibility of the prepared Laughlin state.
The alternative text for this image may have been generated using AI.
a Two-point correlation function Cij between site i and j obtained from results after debiasing and postselection (PS) closely align with ED benchmark. We set Cij = 0 for i≤j. b Site-averaged correlation C(d) over sites separated by d = ∣i − j∣. We include only site index i, j ∈ [2, 13] when calculating C(d) to avoid boundary effect. Error bars indicate 68% confidence intervals obtained by means of percentile bootstrap. Lines are guide to the eye.
We further compute the site-averaged correlation function \(C(d)=\overline{{C}_{j,j+d}}\) as a function of the separation distance d = ∣i − j∣ and observe characteristic fluctuations in the short-range correlation of the prepared Laughlin state. The first two sites near each boundary are excluded to minimize edge effects. The results, shown in Fig. 5(b), reveal a strong correlation hole C(d) < 0 at short distances (d < 4), signifying the underlying repulsive nature of Laughlin state. The medium-range oscillations in C(d) reflect a short-range solid-like order, characteristic of a strongly coupled plasma. Such oscillations are a hallmark of the strongly correlated FQH liquid43. Beyond d ≥ 7, C(d) decays rapidly to zero, representing a featureless and homogeneous liquid at long range. Not only does C(d) from our prepared Laughlin state exhibit qualitative agreement across all distance ranges, but it also quantitatively captures the precise maxima and minima, as well as the spatial extent of the correlation hole.
To demonstrate entanglement behavior beyond pairwise correlation, we directly measured the topological entanglement entropy γtopo44,45 of our prepared state via geometric deformation of the cylinder circumference Ly on the quantum processor. This quantity, which reflects the quantum dimension of anyonic excitations, serves as a robust diagnostic of topological order. We optimized the HVA ansatz \({\left|\psi (\{{\beta }_{j}\})\right\rangle }_{{{\rm{eff}}}}\) for a range of Ly ∈ [6, 10] near the isotropic geometry limit, and applied a randomized measurement protocol46 to estimate the second-order Rényi entropy \({S}_{A}^{(2)}=-ln{{\rm{Tr}}}\,{\rho }_{A}^{2}\) for three different subsystem partition A in the bulk region. (see Supplementary Information for details)
In Fig. 6, the experimentally measured \({S}_{A}^{(2)}\) shows the expected area-law scaling \({S}_{A}^{(2)}=\alpha {L}_{y}-{\gamma }_{{{\rm{topo}}}}\) with a systematic drift to higher entropy due to hardware noise when compared to noiseless simulator benchmark. Fitting the measured second-order Rényi entropy to the area-law scaling, we extracted \(-{\gamma }_{\exp }=-0.92\pm 0.17\) (68% confidence interval by bootstrap resampling of finite-shot randomized measurement estimator, see Methods). For the ideal ν = 1/3 Laughlin state, \(-{\gamma }_{{{\rm{topo}}}}=-ln\sqrt{3}\)47,48,49 and because our system partition introduces two entanglement boundaries, the expected value is \(-{\gamma }_{{{\rm{topo}}}}=-2ln\sqrt{3}\approx -1.10\). The consistent behavior in \({S}_{A}^{(2)}\) and γtopo between our experiments and the theory provides compelling evidence of the topological order of the prepared ν = 1/3 Laughlin state. Our pairwise correlation and entanglement entropy measurements demonstrate the ability to access microscopic structures that underlies topologically ordered states on a quantum processor.
Fig. 6: Topological entanglement entropy.
The alternative text for this image may have been generated using AI.
The second-order Rényi entropy \({S}_{A}^{(2)}\) of the six-qubit subsystem as a function of cylinder circumference Ly. Red triangles represent experimental data obtained on IonQ’s Forte-1 quantum computer using randomized measurements with an ensemble size of NU = 200 unitaries and NM = 300 shots per unitary. The result is compared with noiseless simulation of the variationally optimized HVA (orange square). Dashed lines indicate linear fits to the area law form S2(Ly) = αLy − γ. The noiseless simulation yields αHV A = 0.249 and − γHV A = − 1.09. Experimental fit yields \({\alpha }_{\exp }=0.245\pm 0.021\) and \(-{\gamma }_{\exp }=-0.92\pm 0.17\). Error bars indicate 68% confidence intervals obtained by means of percentile bootstrap. Inset: Schematic of the orbital partition. The system is partitioned into a bulk subsystem A and the environment B, illustrating the two spatial cuts contributing to the entanglement entropy.