A Path Integral Ground State Monte Carlo Algorithm for Entanglement of Lattice Bosons

Based on the comments of the referee, we acknowledge that the resulting rather long paper, combined with its overly technical title and abstract, may have "buried the lead." Motivated by previous and on-going experimental searches able to measure entanglement entropy in ultracold lattice bosons, we aimed to devise an exact ground state numerical method able to make direct predictions of quantum information quantities in these systems. Numerical methods previously reported in the literature for itinerant lattice soft-core bosons to date are based on exact diagonalization or the density matrix renormalization group (DMRG) and are thus limited to small systems or one spatial dimension (or quasi-1D, highly anisotropic 2D systems). DMRG studies also require imposing a local restriction on the bosonic Hilbert space that can lead to large errors in entanglement for weak interactions. Previous quantum Monte Carlo methods for entanglement work only for distinguishable local degrees of freedom (e.g. spin systems) or within a continuous space path integral framework. Moreover, the lack of any open source quantum Monte Carlo code applicable to extended Bose-Hubbard systems has prevented ab initio studies of experimentally accessible quantum simulators. Our intent was to write a pedagogical manuscript on the general algorithim and method, while at the same time, highlighting our novel results on both spatial mode and symmetry resolved entanglement in the Bose-Hubbard model. The latter has generated considerable interest in recent years due to its connection to the accessible entanglement in a quantum many-body phase that could be utilized as a resource for quantum information processing applications. The resulting length of our submitted manuscript (>50 pages), combined with an admitted lack of emphasis on the new results, led the referee to question the suitability of our paper for SciPost Physics.

We disagree with this sentiment. The new contributions of this work include:

1. A ground state QMC algorithm that, while potentially understood by experts, had never appeared in the literature. We derived new Monte Carlo estimators able to measure Rényi entanglement entropies within the lattice path integral framework.
2. New results for the spatial entanglement entropy in the 1D and 2D Bose-Hubbard model both at the quantum critical point, and across the superfluid-insulator phase diagram for considerably larger system sizes than had been previously measured (up to 1024 bosons in 2D).
3. An analysis of subleading corrections to the area law at the 2D quantum critical point. In the current version of the manuscript, we now substantially improve our analysis here (motivated by the referee) and are able to extract the number of Goldstone modes in the resulting gapless phase.
4. Particle number distributions and the symmetry resolved entanglement in the superfluid, critical, and Mott insulating regimes for a small system size (N = 8 particles, Figure 13). We interpreted the observed strong n dependence in the Mott phase in terms of the interaction between holes and doublons. To further explore this physics, we have performed additional simulations for a much larger system (N = 64) far outside the reach of exact diagonalization, and included a new Figure 14.

We now better highlight the physical focus and new results in the substantially modified introduction.

Focusing on the SciPost Physics Acceptance Criteria, we believe our paper "Open(s) a new pathway in an existing or a new research direction, with clear potential for multipronged follow-up work."

While the SWAP estimator had been previously implemented in other quantum Monte Carlo flavors, this has only included those with localized and distinguishable degrees of freedom, e.g. spin models or in the spatial continuum. The extension to lattice models of itinerant bosons that are clearly the most relevant for experimental studies of entanglement (see e.g. Refs. [17-18]) had yet to be accomplished. As evidenced by the lengthy derivation of the SWAP in terms of weights in the continuous imaginary time quantum Monte Carlo algorithm presented here, the generalization to this case is non-trivial. Unlike a spin system, the number of degrees of freedom (worldlines) in a region defined by some subset of spatial modes is a dynamical quantity. As mentioned above, while entanglement in the Bose-Hubbard model has been previously studied in DMRG, these simulations were limited to 1D systems and always implemented a local occupation number restriction (sometimes of only a few bosons per site), resulting in potentially extremely larger errors (of order 100%) in the weakly interacting regime (see e.g. https://journals.aps.org/prb/abstract/10.1103/PhysRevB.105.L121116).

Another benefit of the ground state QMC method presented here is its ability to scale to systems in general dimensions, and based on the comment of referee 1, we have more closely analyzed the sub-leading logarithmic correction to the area law in the 2D Bose-Hubbard model in a 32x32 system composed of 1024 bosons at unit filling. It is predicted that logarithmic corrections should arise due to the existence of a spontaneously broken continuous symmetry in the thermodynamic limit, contributing a universal coefficient due to the presence of a low-energy "tower of states" spectrum and a Goldstone boson. The universal prefactor of this term is known to be related to the number of Goldstone modes, and we report that the value we extract from the QMC is consistent with the value of one for the Bose-Hubbard model. This represents an important first step in studying these, and higher order corrections to the area law which can encode information on the central charges that characterize higher-dimensional conformal field theory, believed to be the fundamental constants that quantify how entropy monotonically decreases under renormalization group flows.

In addition, the growing recent interest in the symmetry resolved entanglement (the number of papers with this in the title/abstract has doubled every year for the last few years) is a strongly motivating factor for our work. The vast majority of these studies have been limited to non-interacting models, as there was simply no technology or method available beyond exact diagonalization for its measurement in interacting systems until now. Computing this quantity, which can have a very large signal to noise ratio in quantum Monte Carlo is no small task, and utilizing the presented algorithm we believe that there is clear potential for follow up work. This may be especially important given the tendency of this quantity to signal quantum phase transitions and its experimental measurement through post-selection. It can also display highly non-trivial dependence on the number of particles in the subsystem (n) and we accentuate this with a new figure 14 in a system of N = 64 bosons.

For these reasons, we believe our substantially modified and refocused manuscript now better highlights our novel physics contribution and thus goes beyond the scope of SciPost Physics Codebases. With the specific reply to the referees we have provided, in combination with the updates to the manuscript, we plan to resubmit our manuscript to SciPost Physics and hope you will continue its consideration for that venue.

- Modified title and abstract to better reflect physics goals of paper
- Rewrote introduction to include a more physics-based introduction and included a list of specific physics-based novel contributions
- When introducing the fits of Eqs. (41) and (54), which depend exponentially on the projection length, we now clarify the origin of the exponential behavior.
- Added 3d diagram illustrating the two-replica configuration space in which we measure second Rényi entropies.
- For consistency, we now report a rescaled projection length β → βt when reporting results in terms of the dimensionless interaction strength U/t.
- Numerical values of the fit parameters a,b,c are now reported. The universality of the logarithmic correction is now mentioned and its direct proportionality to the number of Goldstone modes. Relevant literature has been cited.
- We have added citations to large number of references of these one-dimensional critical point estimations. We also now comment on how a finite size scaling analysis of our S_2 vs. U/t results for large systems could be used to identify the critical point.
- We now mention that error bars are computed using the jackknife method in the main text.
- The terminal output is now thoroughly explained in Appendix A and the code's documentation.
- We have now included an MIT open source license to the code repositories.
- The code now compiles with g++-12. We have fixed all other issues that were causing compiler warnings. After testing with clang++ and g++-12, no warnings should be given now.
- All indicated typos have been fixed.