Stochastic Series Expansion Quantum Monte Carlo for Rydberg Arrays

- Detailed description of the algorithm
- Finite-temperature calculations are interesting and could compliment previous studies
- The paper provides a MonteCalro alternative to the existing ED- and tensor network-based numerical tools usually used for Rydberg atoms

- The reference list in the introduction is a bit outdated (see the report)
- The manuscript does not report new physics
- The limitations of the algorithm have not been analyzed
- No comparison to alternative numerical methods previously used for Rydberg atoms (MPS, iDMRG, constrained DMRG, ED)

The paper is purely methodological: it does not report new phenomena but introduces the numerical tool capable to reproduce the results obtained in the experiments. As for the methodological paper it however lacks a few essential ingredients: First, the manuscript does not provide a careful analysis of the limitations in terms of computational costs, systems sizes, accuracy etc. The authors limited themselves to system sizes available in experiments already a few years ago. Since the experiments on Rydberg atoms are progressing fast, in particular, in the number of trapped atoms, it would be interesting to have more systematic analysis of the limitations of the method in terms of reachable system sizes and other parameters that would put the method in a perspective. Second, I was lacking a comparison in terms of an accuracy and computational costs with respect to other available methods reporting results also on much larger system sizes (e.g. MPS results from the original Nature papers; iDMRG from Rader&Lauchli arXiv:1908.02068; and from Guidici et al. Phys. Rev. B 99, 094434; constrained DMRG from SciPost Phys. 6, 033 (2019)).

In addition, I would like to point out, that references in the introduction list some early papers and refer to the phase diagrams that have been recently corrected with more advanced numerical simulations in:
Phys. Rev. Lett. 122, 017205 (2019)
Phys. Rev. B 99, 094434 (2019)
Phys. Rev. Research 3, 023049 (2021)
And there are also relevant unpublished works:
Rader&Lauchli: arXiv:1908.02068
Chan et al: https://www.youtube.com/watch?v=bptOdSHo2dI&t=1664s

Finally, note, that even for the specific set of parameters chosen in the paper (\Omega=1; R_b=1.2) in 1D there is still an open question: the transition to the ordered Z_2 phase could either be 1st order or Ising and to the best of my knowledge the two possibilities have not been resolved yet for the van der Waals potential. By the way, Ref.14 at the bottom of p.14 is quite misleading since it refers to another transition out of Z_3 phase that takes place at larger values of R_b>2.

To summarize, the manuscript describes numerical tool complimentary to the exiting numerical methods used for Rydberg atoms and undoubtedly deserves to be published. However, the paper in its current form misses a few essential methodological ingredients that has to be addressed before the manuscript can be accepted for publication in SciPost.

See the report

1. Provides a detailed explanation of the formulation of the stochastic series expansion quantum Monte Carlo (SSEQMC) method for currently relevant Rydberg atom arrays.
2. Step by step implementation protocol for the algorithm is presented for both finite temperature and the zero temperature projection method, along with the energy estimator, which can be easily generalized to other similar systems.
3. The line cluster algorithm proposed in this manuscript is shown to perform significantly better than a naive cluster algorithm, especially for the phase space of interest in the Rydberg atom arrays.
4. The appendices describe in detail the calculations of the probabilities involved in the Monte Carlo simulation, and further optimizations to improve performance. As this can also be easily extended to more general systems, practitioners of this method will benefit directly from these explanations.

1. As the main focus of this manuscript is the line cluster update, it is essential to note that the expected range of parameters in which the update is expected to lead to efficient sampling is relatively small, and depends crucially on the values of W^{(1)-(4)}. The region of parameter space which corresponds to long line clusters is expected to suffer from small acceptance probabilities (this effect does not show up in a significant way at the system sizes considered here).
2. The performance of the line cluster is strongly tied to the region of parameter space being sampled, and thus is not naively expected to generalize to other Ising-like systems with interactions which break Z_2.

The authors have presented the implementation of a quantum Monte Carlo method for Hamiltonians describing Rydberg atom arrays. The central contribution of this manuscript is the development of a line cluster method which performs significantly better than a naive cluster method for the parameter range of interest.
All the general acceptance criteria are met by this manuscript. From the list of expectations provided by the journal, this manuscript satisfies "Open a new pathway in an existing or a new research direction, with clear potential for multipronged follow-up work". One must note however that such line cluster methods have been used (arXiv:1812.05326,2008.11206), without a study of the improvement in ergodicity. The current manuscript fills this gap and holds promise for further improvements and developments in sophisticated cluster methods.
Typo: On page 5, above equation 7, I think you mean "operator sequence is less than 80%"

1. As far as cluster methods are concerned, one of the first checks for efficiency is usually performed through an examination of the success rate of proposed moves. This quantity can be easily extracted from simulations with small error bars and usually proves to be much less noisy than auto-correlation times. The manuscript will benefit from a short discussion of the behavior of the success rate as a function of Hamiltonian paramters.
2. (Optional (up to the authors' discretion); may prove to be too computationally involved to add to the current manuscript) Long line clusters are invariably bound to fail due to the large number of two-site operators which contribute to them (especially away from the special point in parameter space where the values of W^i are favorable). An alternative to this strategy would be to allow a Swendsen-Wang style cluster building algorithm, where operators are added to the cluster using a probability calculated from the different matrix elements of the operator. Something along these lines is suggested in arXiv:2009.03249 for a quantum clock model, and was found to work better than naive cluster of line updates. It would be intriguing to try to extend the same to the Hamiltonian studied in the current manuscript, and may lead to improvements in ergodicity.

Arrays of Rydberg atoms are a powerful platform to realize strongly-interacting quantum many-body systems. It can be translated to an effective transverse field Ising model (TFIM) with long range interactions and vertical field, without sign problem. In this work, the authors develop a cluster scheme for this model. This work is an improvement of a previous SSE algorithm of normal TFIM.

I think this work sounds interesting and worthy to be published. However, I still have a question and hope authors can answer it. A similar algorithm for TFIM with a vertical field is also used in the supplementary materials of Phys. Rev. B 103, 104416 (2021). In these methods, the weights before/after flipping cluster are not equal, so I think you have to create the cluster first and calculate its weight to get the accept-probability. If you accept this update, then you can replace the configuration memory, else you have to keep the old arrays of operator string and spins. I think this step will cost more time than normal TFIM with accept-probability 0.5. Because in normal TFIM case, people can judge accept or not first (probability is independent on cluster), then update related arrays. I doubt it makes the computation complexity of this algorithm is no longer \beta*size as in normal SSE. I am just interested in whether you have good way to deal with this problem?

1- clear description of the algorithm
2- quality of the numerics
3- connection to experimental results

1- presentation of the limitations of the algorithm

This paper extends quantum Monte Carlo methods previously used for the transverse-field Ising model to a Hamiltonian describing an array of Rydberg atoms. In contrast to the pure transverse-field Ising model, this model lacks a spin-flip symmetry so that the conventional cluster updates change the configuration weight. The authors remedy this by introducing an a-posteriori acceptance probability for each cluster flip and find that the construction of space-local clusters improves ergodicity. Further, an efficient estimator for the ground state energy in the projector algorithm is derived.

As a proof of concept, one and two-dimensional systems are simulated close to the checkerboard quantum phase transition.

These results are of interest and should inspire further investigations given the recent interest in Rydberg systems. The paper is also otherwise clear and well written, up to some concerns detailed below. After these concerns are addressed, I think that the paper satisfies all criteria for publication in SciPost Physics.

1- The system sizes considered are small compared to what is possible e.g. in Ref. [30]. In the current manuscript, the reason for this is not clear. Refs. [8,9] suggest the choices $N=51$ and $N=256$ were inspired by experimental setups. This should be clearly stated in the main text. Is $N=256$ already the limit of what is possible? Are the autocorrelation times the bottleneck for this or something else?

2- The authors present data for a one-dimensional model. Due to the method focus of this paper, they should comment on the applicability of DMRG in this case. In particular, in Refs. [14-16], DMRG is already used to study Rydberg systems. (Why) should we use QMC in one dimension at $R_b = 1.2$?

3- On p.11 there is the statement “the operators in the sequence $S_M$ will most likely be dominated by one of $\{W_{i,j}^{(1)}, W_{i,j}^{(4)}\}$”. Is this really the case? In the checkerboard phase, I would naively expect that operators with differing states on their legs ($W^{(2)}$ and $W^{(3)}$) proliferate. In the shown example (Fig. 3), $W^{(3)}$ also has the highest weight out of all four. Please correct me if I am wrong.

4- Before introducing the ground state energy estimator in Sec. 3.3, it would be helpful to reference the already available ground state estimators such as the overlap method detailed in [A. W. Sandvik, PRL 95, 207203]. Does the new estimator have advantages over existing methods?

5- In the discussion of the autocorrelation times in the two-dimensional model (top of p. 17), there is an interesting point about normalizing autocorrelation times by the number of clusters formed. What exactly is meant by the “number clusters formed”? Does this number count the number of attempted cluster flips or the number of actually flipped clusters? Please clarify.
I would expect the actual size of the cluster in terms of flipped operators to play a role as well. Is there something one can say about the scaling of the acceptance probability for the multibranch and the line updates? My guess would be exponential decay in the number of operators per cluster. Then the line update would work better because its line clusters contain fewer operators.

6- Have the authors considered adding a constant $\epsilon$ to $C_{ij} = |\min(\dots)| + \epsilon$? This is routinely done in the context of the SSE directed loop update (e.g. Ref. [26]) and sometimes improves ergodicity by avoiding such W=0 situations. I doubt that this will fundamentally cure the autocorrelation issues, so there is no need to redo simulations. Mentioning the possibility of such an offset may nevertheless be helpful to the reader.

7- In the discussion of the results, a quantum phase transition and the checkerboard phase is mentioned. This should be shortly introduced to give the reader some overview over what is being studied.

8- currently, in the Figures, $\delta/\Omega$ is given as a ratio of $\Omega=1$ but $E$ and $T$ in Fig. 4 are not. The authors may consider to make it consistent by writing $E/Ω$ and $T/Ω$.

9- on page 15, the inline equation for $|M_s|$ enters into the margin. This should be fixed.