REVIEW 3 major objections 6 minor 33 references
Efficient classical simulation of two-dimensional long-range systems: Rydberg arrays and beyond
T0 review · 3 major / 6 minor · reviewed 2026-07-08 · glm-5.2
Pith's one-line read O(N³) to O(N): Classical simulation of 2D long-range quantum systems
desk verdict Adaptive-basis sampling reduces TNS-tVMC local energy cost from O(N³) to O(N); the 10×10 Rydberg application is promising but under-benchmarked at D=4 with no error bars. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
Adaptive-basis sampling for tensor network VMC: decomposing a long-range two-body Hamiltonian into terms diagonal in distinct product bases, sampling each term in its own basis via separate Markov chains, and combining the parameter updates using the linearity of the time-dependent variational principle.
What would settle it
Attempt to apply the adaptive-basis framework to a long-range Hamiltonian whose off-diagonal terms cannot be simultaneously diagonalized in any finite set of product bases. If no such decomposition exists, the O(N) reduction does not apply.
Extended reading notes
Core claim
The central mechanism is the observation that the local energy bottleneck in long-range VMC — the O(N²) off-diagonal terms each costing O(N) to evaluate — can be eliminated by sampling in bases adapted to each component of the Hamiltonian. For a tensor network state, changing the sampling basis requires only applying a product of one-site rotation operators, which is efficient. By splitting the XY Hamiltonian into H_X (diagonal in the x-basis) and H_Y + H_Z (diagonal in the y-basis), running two Markov chains, and combining the resulting parameter derivatives, the overall local energy evaluation becomes O(N) per sweep. This enables scalable real-time dynamics and ground-state optimization in
Load-bearing premise
The method requires that the long-range Hamiltonian can be decomposed into terms that are each diagonal in some product basis, and that the corresponding basis-change operators are products of one-site rotations applicable to the tensor network. This holds for the XY model but the breadth of Hamiltonians satisfying this condition is not characterized.
Editorial extensions
If this is right
- Classical simulation of 2D Rydberg and trapped-ion quantum simulators with long-range 1/r³ or 1/r⁶ interactions becomes feasible at system sizes relevant to current experiments, providing benchmarks that were previously unavailable.
- The 10×10 dipolar XY result clarifies that the ideal adiabatic protocol produces strong antiferromagnetic correlations, directing attention to experimental imperfections (such as reduced detuning amplitude) as the source of suppressed correlations in the actual experiment.
- The O(N) scaling opens the door to studying real-time dynamics, spectral functions, and sign-problem-afflicted ground states in two-dimensional long-range systems that are currently inaccessible to deterministic tensor network methods.
- The method establishes a workflow where classical simulation first explores parameter space to identify interesting regimes, and quantum simulators take over when entanglement growth exceeds tensor network representability.
Reading between the lines
- The adaptive-basis technique is specific to tensor network states because the basis-change operator is a product of one-site rotations; neural quantum states, which lack this structure, would not benefit from the same speedup. This creates a performance asymmetry between TNS-VMC and NQS-VMC for long-range systems.
- The applicability of the method depends on whether a given Hamiltonian can be decomposed into terms each diagonal in some product basis. Hamiltonians with off-diagonal couplings along multiple non-commuting directions simultaneously would require a different or more general approach.
- The trade-offs identified — loss of zero-variance principle and loss of symmetry-sector sampling — suggest that for ground-state calculations in fixed-basis-diagonal Hamiltonians, the adaptive approach may not always be superior, and a hybrid strategy switching between adaptive and fixed-basis sampling depending on the task could be optimal.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. This manuscript introduces an adaptive-basis framework for tensor-network variational Monte Carlo (TNS-tVMC) that reduces the local energy evaluation cost for arbitrary long-range two-body Hamiltonians from O(N³) to O(N) per Monte Carlo sweep. The key idea is to decompose the Hamiltonian into components diagonal in different product bases (e.g., x-basis and y-basis for the XY model), sample each basis independently, and combine the parameter updates linearly. The method is applied to simulate the 10×10 dipolar XY Rydberg dynamics from the experiment of Ref. [20] (Chen et al., Nature 2023), which was previously beyond classical reach. The authors also address the TDVP pathology at product states using projected tVMC (p-tVMC) for the initial step. The 4×4 benchmark against exact diagonalization shows systematic convergence with bond dimension, and the 10×10 results show sample-size convergence. The central physical conclusion is that the ideal adiabatic protocol does produce long-range AFM correlations, attributing the experimental suppression to imperfections in the realization.
Significance. The algorithmic contribution — reducing the local energy cost from O(N³) to O(N) for long-range two-body Hamiltonians in TNS-tVMC — is significant and addresses a genuine computational bottleneck for 2D long-range systems. The derivation is grounded in the linearity of the TDVP equation and the efficiency of basis changes for tensor networks, making the core claim internally consistent and parameter-free in its complexity analysis. The application to the 10×10 dipolar XY Rydberg protocol provides a concrete classical benchmark for a state-of-the-art quantum experiment. The p-tVMC initialization for product states is a useful technical contribution. The authors provide open data (Ref. [30]) and use reproducible code (JAX-based implementation), which strengthens the work. The physical conclusion about the origin of suppressed correlations in the experiment is of direct interest to the Rydberg array community.
major comments (3)
- The 10×10 results use only D=4 (Table I, Fig. 3), with no bond-dimension convergence study at this system size. The 4×4 benchmark (Fig. 2) shows convergence from D=4 to D=5, but for 2D real-time dynamics evolved to T=9.68, entanglement growth may require larger D. The conclusion that 'the suppression observed experimentally is unlikely to result solely from an intrinsic breakdown of the ideal adiabatic protocol' (end of Results section) rests on the D=4 correlations being trustworthy. Without at least one additional bond dimension at 10×10, or a quantitative argument for why D=4 suffices at this system size and evolution time, this load-bearing claim is not robustly established.
- No statistical error bars are shown on any of the 10×10 observables (Fig. 3b, Fig. 3c, Fig. 4). Fig. 4 demonstrates sample-size convergence between N_s=81920 and 163840 for m²_x, but this is a consistency check, not a statistical uncertainty estimate. The residual energy gap above the ground state (Fig. 3a, 'slightly above the ground-state energy' for t≳6) could reflect either genuine non-adiabaticity or insufficient variational expressivity; error bars on the observables and the energy would help distinguish these explanations.
- The title and conclusion state the method applies to 'arbitrary long-range two-body Hamiltonians,' but the adaptive-basis approach requires that the Hamiltonian can be decomposed into terms each diagonal in some product basis. The Conclusion acknowledges this is limited to 'Hamiltonians whose terms can be organized in locally accessible product bases,' but the breadth of this class is not characterized. For Hamiltonians with off-diagonal terms in multiple non-commuting directions (e.g., Heisenberg with all three σ^α σ^α terms), it is unclear whether the method applies or how many Markov chains would be needed. The scope claim in the title and abstract should be qualified to match what is demonstrated.
minor comments (6)
- The abstract contains a grammatical error: 'the pathology of evolving from product states within of tensor network VMC' — 'within of' should be 'within'.
- In the 'Adaptive basis' section, the notation for the two Markov chains (sample_x and sample_y) could be stated more precisely. It would help to explicitly state that U is the tensor product of Hadamard-like rotations mapping the x-basis to the y-basis, and to clarify that the two chains are run independently.
- Fig. 2 axis labels contain Unicode rendering issues (e.g., /u1D461, /u1D43B). These should be fixed for the final version.
- The wall-time data in Table I (367–793 s/step for 10×10) provides useful scaling information, but it would strengthen the O(N) claim to include a brief comparison with the O(N³) scaling of a naive implementation at the same system size, even if only estimated.
- The trade-offs section (bottom of p. 5) notes that adaptive sampling loses the zero-variance principle and Sz conservation via sampling. It would be useful to briefly quantify the overhead: how many more samples does adaptive sampling require compared to fixed-basis sampling in the ground-state VMC?
- Ref. [24] (Wu and Nys, arXiv:2512.06768) appears to be a closely related concurrent work on TNS-tVMC in 2D. The relationship between the present work and Ref. [24] should be clarified, particularly regarding the novelty of the adaptive-basis technique versus the tVMC framework itself.
Simulated Author's Rebuttal
We thank the referee for a careful and constructive report. The referee raises three major comments: (1) the absence of bond-dimension convergence at 10×10, (2) the lack of statistical error bars on 10×10 observables, and (3) the scope of the 'arbitrary long-range two-body Hamiltonians' claim. We agree that all three points warrant revision. For (1), we will add a D=6 calculation at 10×10 to demonstrate convergence of the central physical conclusion. For (2), we will add statistical error bars throughout. For (3), we will qualify the title, abstract, and conclusion to accurately characterize the applicable Hamiltonian class. We believe these revisions substantially strengthen the manuscript without altering the core algorithmic contribution.
read point-by-point responses
-
Referee: The 10×10 results use only D=4, with no bond-dimension convergence study at this system size. The conclusion that the suppression observed experimentally is unlikely to result solely from an intrinsic breakdown of the ideal adiabatic protocol rests on the D=4 correlations being trustworthy. Without at least one additional bond dimension at 10×10, or a quantitative argument for why D=4 suffices, this load-bearing claim is not robustly established.
Authors: The referee is correct that bond-dimension convergence at 10×10 is needed to robustly establish the central physical claim. We will address this by adding a D=6 calculation at 10×10 for the full AFM protocol. We note several factors that support the adequacy of D=4 even before the new data: (i) the 4×4 benchmark shows convergence between D=4 and D=5, with the differences already small; (ii) the system is evolved from a product state with a relatively short total time T=9.68 (corresponding to 2 μs in the experiment), limiting entanglement growth; (iii) the energy at late times plateaus only slightly above the variational ground-state energy, suggesting the ansatz remains reasonably expressive; and (iv) the spatial symmetry constraint reduces the effective parameter space. Nevertheless, we agree that these arguments are not a substitute for an explicit convergence check at the target system size. The D=6 data will be included in the revised manuscript, and we will explicitly state whether the physical conclusion is preserved. revision: yes
-
Referee: No statistical error bars are shown on any of the 10×10 observables. Fig. 4 demonstrates sample-size convergence between N_s=81920 and 163840 for m²_x, but this is a consistency check, not a statistical uncertainty estimate. The residual energy gap above the ground state could reflect either genuine non-adiabaticity or insufficient variational expressivity; error bars on the observables and the energy would help distinguish these explanations.
Authors: We agree. The sample-size comparison in Fig. 4 serves as a convergence check but does not constitute a proper statistical uncertainty estimate. We will add error bars to all 10×10 observables (Fig. 3b, 3c, and Fig. 4) computed from independent Monte Carlo runs or binning analysis. We will also add error bars to the energy evolution in Fig. 3a. Regarding the residual energy gap above the ground state for t≳6: we agree that this could reflect either genuine non-adiabaticity or variational expressivity limitations, and we will discuss this distinction explicitly in the revised text. The error bars will help quantify the statistical component, though we note that disentangling non-adiabaticity from variational bias requires the bond-dimension convergence study addressed in the previous point. revision: yes
-
Referee: The title and conclusion state the method applies to 'arbitrary long-range two-body Hamiltonians,' but the adaptive-basis approach requires that the Hamiltonian can be decomposed into terms each diagonal in some product basis. The Conclusion acknowledges this is limited to 'Hamiltonians whose terms can be organized in locally accessible product bases,' but the breadth of this class is not characterized. For Hamiltonians with off-diagonal terms in multiple non-commuting directions (e.g., Heisenberg with all three σ^α σ^α terms), it is unclear whether the method applies or how many Markov chains would be needed. The scope claim in the title and abstract should be qualified to match what is demonstrated.
Authors: The referee correctly identifies an overstatement in the title and abstract. The method applies to Hamiltonians whose two-body terms can be decomposed into groups, each diagonal in a single product basis — not to literally arbitrary two-body Hamiltonians. For the XY model, two bases (x and y) suffice. For a Heisenberg model with σ^xσ^x + σ^yσ^y + σ^zσ^z, three Markov chains (one per Pauli basis) would be needed, and the method still applies with O(N) cost per chain. More generally, any Hamiltonian whose interaction terms are sums of commuting two-body operators in a finite number of product bases is covered. However, Hamiltonians with two-body terms that are simultaneously off-diagonal in multiple non-commuting directions at each site would not be directly amenable. We will revise the title, abstract, and conclusion to replace 'arbitrary long-range two-body Hamiltonians' with more precise language such as 'long-range two-body Hamiltonians decomposable into product-basis-diagonal components,' and we will add a brief discussion characterizing the applicable class, including the Heisenberg example and the scaling of the number of Markov chains with the number of required bases. revision: yes
Circularity Check
No circularity found: the O(N) adaptive-basis derivation is self-contained from TDVP linearity, and self-citations are for supporting tools, not load-bearing premises.
full rationale
The paper's central algorithmic claim — reducing local energy evaluation from O(N³) to O(N) — is derived from first principles: the linearity of the TDVP equation (Eq. 6: Sθ̇ = −iF) allows decomposing the force into components, each evaluated in a basis where the corresponding Hamiltonian component is diagonal (Eq. 9: θ̇ = θ̇_X + U†θ̇_Y). No step in this chain reduces to its own inputs by construction. The basis-change operator U is a tensor product of one-site rotations, which is independently verifiable as efficient for TNS. The Rydberg simulation parameters come from the external benchmark Ref. [20] (Chen et al., Nature 2023), not from the authors' own prior work. Self-citations [24] (tVMC for TNS), [25] (column direct sampling), and [19] (p-tVMC) provide supporting algorithmic tools but are not load-bearing for the central complexity-reduction claim, which stands on the mathematical argument alone. The 4×4 benchmark against exact diagonalization (Fig. 2) provides independent validation. The lack of bond-dimension convergence for the 10×10 system is a correctness risk, not a circularity issue. No fitted parameter is renamed as a prediction, no uniqueness theorem is invoked to forbid alternatives, and no ansatz is smuggled through self-citation.
Assumptions & free parameters
free parameters (4)
- TNS bond dimension D =
4 (10×10), 4-5 (4×4)
- MC sample size Ns =
81920-163840
- Time step dt =
0.02
- δ₀ =
19.48
assumptions (4)
- standard math The TDVP equation (Eq. 6) Sθ̇ = -iF governs the variational parameter evolution.
- domain assumption Tensor network states admit efficient basis changes via tensor products of one-site operators.
- domain assumption The Hamiltonian can be decomposed into terms diagonal in locally accessible product bases.
- standard math The p-tVMC fidelity maximization (Eq. 10) provides a valid initial step from product states.
Cite this review
Pith. "Pith review of Efficient classical simulation of two-dimensional long-range systems: Rydberg arrays and beyond." pith.science (2026). https://pith.science/paper/NMSXWAFN
@misc{pith2026260705178,
author = {Pith},
title = {Pith review of: Efficient classical simulation of two-dimensional long-range systems: Rydberg arrays and beyond},
year = {2026},
howpublished = {\url{https://pith.science/paper/NMSXWAFN}},
note = {Machine review of arXiv:2607.05178}
}
abstract
In variational Monte Carlo (VMC) calculations of $N$-site quantum systems with arbitrary all-to-all two-body interactions, evaluating the local energy generally costs $O(N^3)$. We introduce a new framework that reduces this cost to $O(N)$ for tensor network states, capable of scalable and accurate computation of real-time dynamics and ground states. As a result, we obtain accurate simulations of the adiabatic real-time protocol of a $10\times10$ dipolar XY model realized in a Rydberg simulator [C. Chen et al., Nature 616, 691 (2023)], which was previously beyond the reach of classical simulation. Going beyond quantum experiments, we also directly perform ground state VMC to compare with the adiabatic state preparation. Our work demonstrates tensor network VMC as a powerful classical simulator for long-range quantum platforms such as Rydberg and ion-trap simulators, which are currently in urgent need of scalable classical benchmarking tools. As a separate technical contribution, we resolve the pathology of evolving from product states within of tensor network VMC.
Figures
Reference graph
Works this paper leans on
-
[20]
C. Chen, G. Bornet, M. Bintz, G. Emperauger, L. Leclerc, V. S. Liu, P. Scholl, D. Barredo, J. Hauschild, S. Chatterjee, M. Schuler, A. M. L¨ auchli, M. P. Zaletel, T. Lahaye, N. Y. Yao, and A. Browaeys, Continuous symmetry breaking in a two-dimensional Rydberg array, Nature616, 691 (2023)
work page 2023
-
[30]
efficient classical sim- ulation of two-dimensional long-range systems: Ry- dberg arrays and beyond
Y. Wu, Data repository for “efficient classical sim- ulation of two-dimensional long-range systems: Ry- dberg arrays and beyond”,https://github.com/ yantaow/open_data/tree/main/wu2026efficient (2026). 7 Supplemental Material S-1. UNITS CONVERSION IN THE R YDBERG DYNAMICS Using notation in Ref. [20], the Rydberg dynamics is prescribed by HXY =− J 2 X i<j (...
work page 2026
-
[1]
C. Gross and I. Bloch, Quantum simulations with ultracold atoms in optical lattices, Science357, 995 (2017)
work page 2017
-
[2]
A. Browaeys and T. Lahaye, Many-body physics with individually controlled Rydberg atoms, Nature Physics16, 132 (2020)
work page 2020
- [3]
- [4]
-
[5]
M. Saffman, T. G. Walker, and K. Mølmer, Quan- 6 tum information with Rydberg atoms, Reviews of Modern Physics82, 2313 (2010)
work page 2010
-
[6]
G. Vidal, Efficient simulation of one-dimensional quantum many-body systems, Physical Review Let- ters93, 040502 (2004)
work page 2004
Show all 33 references
-
[7]
Verstraete and J
F. Verstraete and J. I. Cirac, Renormalization algo- rithms for quantum-many body systems in two and higher dimensions (2004), arXiv:cond-mat/0407066
2004 arXiv
-
[8]
Carleo, L
G. Carleo, L. Cevolani, L. Sanchez-Palencia, and M. Holzmann, Unitary dynamics of strongly inter- acting Bose gases with the time-dependent varia- tional Monte Carlo method in continuous space, Physical Review X7, 031026 (2017)
2017
-
[9]
Schmitt and M
M. Schmitt and M. Heyl, Quantum many-body dy- namics in two dimensions with artificial neural net- works, Physical Review Letters125, 100503 (2020)
2020
-
[10]
Liu, Y.-Z
W.-Y. Liu, Y.-Z. Huang, S.-S. Gong, and Z.-C. Gu, Accurate simulation for finite projected entangled pair states in two dimensions, Phys. Rev. B103, 235155 (2021)
2021
-
[11]
Vieijra, J
T. Vieijra, J. Haegeman, F. Verstraete, and L. Van- derstraeten, Direct sampling of projected entangled- pair states, Physical Review B104, 235141 (2021)
2021
-
[12]
Wu and Z
Y. Wu and Z. Dai, Algorithms for variational Monte Carlo calculations of fermion projected entangled pair states in the swap gates formulation and the detailed balance of tensor network sequential sam- pling, Chinese Physics B35, 020502 (2026)
2026
-
[13]
Carleo and M
G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science355, 602 (2017)
2017
-
[14]
Sharir, Y
O. Sharir, Y. Levine, N. Wies, G. Carleo, and A. Shashua, Deep autoregressive models for the ef- ficient variational simulation of many-body quan- tum systems, Physical Review Letters124, 020503 (2020)
2020
-
[15]
Mendes-Santos, M
T. Mendes-Santos, M. Schmitt, and M. Heyl, Highly resolved spectral functions of two-dimensional sys- tems with neural quantum states, Phys. Rev. Lett. 131, 046501 (2023)
2023
-
[16]
Mauron, Z
L. Mauron, Z. Denis, J. Nys, and G. Carleo, Predict- ing topological entanglement entropy in a Rydberg analogue simulator, Nature Physics21, 1332 (2025)
2025
-
[17]
Sorella, Generalized Lanczos algorithm for vari- ational quantum Monte Carlo, Physical Review B 64, 024512 (2001)
S. Sorella, Generalized Lanczos algorithm for vari- ational quantum Monte Carlo, Physical Review B 64, 024512 (2001)
2001
-
[18]
Sinibaldi, C
A. Sinibaldi, C. Giuliani, G. Carleo, and F. Vicen- tini, Unbiasing time-dependent variational Monte Carlo by projected quantum evolution, Quantum7, 1131 (2023)
2023
-
[19]
Gravina, V
L. Gravina, V. Savona, and F. Vicentini, Neural projected quantum dynamics: a systematic study, Quantum9, 1803 (2025)
2025
-
[21]
J.-L. Chen, T. Xiang, and Y. Wu, Supplemental ma- terial (2026)
2026
-
[22]
Sbierski, M
B. Sbierski, M. Bintz, S. Chatterjee, M. Schuler, N. Y. Yao, and L. Pollet, Magnetism in the two- dimensional dipolar XY model, Physical Review B 109, 144411 (2024)
2024
-
[23]
Haegeman, J
J. Haegeman, J. I. Cirac, T. J. Osborne, I. Piˇ zorn, H. Verschelde, and F. Verstraete, Time-dependent variational principle for quantum lattices, Physical Review Letters107, 070601 (2011)
2011
-
[24]
Wu and J
Y. Wu and J. Nys, Real-time dynamics in two dimensions with tensor network states via time-dependent variational monte carlo (2026), arXiv:2512.06768 [cond-mat.str-el]
2026 arXiv
-
[25]
T. Chen, J. Liu, Y. Wu, P. Zhang, and Y. Deng, Variational monte carlo (vmc) with row- update projected entangled-pair states (peps) and its applications in quantum spin glasses (2026), arXiv:2601.20608 [cond-mat.dis-nn]
2026
-
[26]
Assaraf and M
R. Assaraf and M. Caffarel, Zero-variance principle for Monte Carlo algorithms, Physical Review Let- ters83, 4682 (1999)
1999
-
[27]
Hauschild and F
J. Hauschild and F. Pollmann, Efficient numerical simulations with tensor networks: Tensor Network Python (TeNPy), SciPost Physics Lecture Notes5, 10.21468/SciPostPhysLectNotes.5 (2018)
2018 doi
-
[28]
Hauschild, J
J. Hauschild, J. Unfried, S. Anand, B. Andrews, M. Bintz, U. Borla, S. Divic, M. Drescher, J. Geiger, M. Hefel, K. H´ emery, W. Kadow, J. Kemp, N. Kirchner, V. S. Liu, G. M¨ oller, D. Parker, M. Rader, A. Romen, S. Scalet, L. Schoonderwoerd, M. Schulz, T. Soejima, P. Thoma, Y....
2024
-
[29]
Bradbury, R
J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, JAX: composable transformations of Python+NumPy programs (2018)
2018
-
[31]
˜τ= τ J ℏ = 0.3µs×2π×0.77MHz = 1.45
-
[32]
˜δ0 = 2π×15MHz ℏ J = 15M Hz 0.77M Hz = 19.48
-
[33]
Dropping tilde gives the dimensionless variables
˜T= JT ℏ = 9.68 The initial state is prepared as the ground in the limit ˜δ0 =∞. Dropping tilde gives the dimensionless variables. S-2. P-TVMC FOR THE FIRST STEP OF THE R YDBERG DYNAMICS We use p-tVMC only for the first time step, where the initial N´ eel product state |ϕ⟩ ≡ |...
Reviewed July 8, 2026 · model on record in the stance chip above.
Discussion (0). Sign in to comment.