Dynamical properties of quantum many-body systems with long-range interactions

,

Introduction.-Quantummany-body systems with longrange (LR) interactions exhibit different and exotic properties compared with their short-range cousins, as the LR nature of the interaction differentiates them from many universally accepted long-wavelength and lowenergy customs governing the short-range ones over the years.For example, the well-known Hohenberg-Mermin-Wagner theorem [1,2] that forbids spontaneous symmetry-breaking of continuous symmetry at finite temperature in low dimensions can be easily circumvented and LR interactions can generate interesting finite temperature transitions [3][4][5][6][7][8][9] and new critical phenonema [10][11][12][13][14][15].The bedrock in the research of highly entangled quantum matter -the area law scaling of the entanglement entropy -can also be bypassed in LR systems, and the consequent new scaling behavior points towards new guiding principle of quantum entanglement that awaits to be worked out [8,[16][17][18][19][20][21].
Moreover, recently the field of LR quantum-many systems becomes even more active due to their fast experimental realizations, such as the Rydberg atom arrays with long-range van der Waals or dipolar couplings where topological ordered phases, emergent glassy behavior and quantum criticality have been suggested and realized [22][23][24][25][26], magic angle twisted bilayer Graphene (TBG) and 2D quantum moiré materials in which flat-band topology and long-range Coulomb interaction give rise to a plethora of correlated phases beyond semi-classical or band-theory description , as well as the quantum gases coupled to optical cavities [70] and programmable quantum simulators [26,[71][72][73][74][75].
Despite such fast developments, theoretical and numerical investigations on the dynamical properties of the LR quantum many-body systems are however still lacking.This is mainly due to the fact that dynamical properties, such as spectral functions [63,[76][77][78][79][80][81][82][83][84], are usually difficult to compute without approximation in analytic the-ory and numerical simulations, even for the short-ranged systems.And therefore by now there only exist few perturbative works such as Refs.[85][86][87][88][89], which are mainly valid at various mean-field limits where the fluctuations are suppressed, and previous algorithmic developments in non-perturbative numerical approaches for LR system are mainly focused on classical and 1D systems [90,91].However, in the aforementioned experiments of 2D quantum LR systems, it is actually the dynamical information that can be easily detected by means of neutron scattering, nuclear magnetic resonance, scanning tunneling spectroscopy, nonlinear and non-equilibrium transport and optical probes, etc.
To overcome the dilemma between the fast experimental developments and the slow theoretical reality in LR quantum many-body systems, the need to develop and carry out unbiased approaches such as large-scale quantum Monte Carlo (QMC) simulations, to systematically investigate the dynamical properties therein is obvious.And only in this way, can one fully reveal the interplay between the LR interaction and quantum topology and fluctuations to explain the fascinating experimental outcomes and predict new ones.This is the focus of our paper.Here we develop and employ the stochastic series expansion (SSE) QMC [26,83,92,93] simulation for the LR quantum many-body systems, to compute the energy spectra of the 2D spin-1/2 Heisenberg model with 1/r α interaction where α is the decay exponent, as shown in Fig. 1.With the interaction types of ferromagnetic (see Fig. 1 (a)) and antiferromagnetic (staggered without introducing frustration, Fig. 1 (b)), we find the explicit range in α for the short-range Goldstone-type (gapless), anomalous Goldstone-type (gapless but with varying dynamical exponent) and Higgs-type (where the spectra are gapped) spectra.As shown in Fig. 1 (d) and (e), accompanied by spin-wave theory (SWT) analysis [12,[94][95][96][97][98], our re- sults reveal how the long-range interactions induce a mass to the Goldstone mode via the generalized Higgs mechanism [85].Moreover, different from the conventional wisdom for systems with long-range interactions [88,[99][100][101], we find explicit case -the staggered antiferromagnetic model -where the gapped excitation exists even when the Hamiltonian is extensive.Therefore our work provides the first set of unbiased dynamical data of LR quantum many-body systems where universally accepted low-energy physics are substantially modified.Implications of ongoing experiments in quantum simulators and 2D quantum moiré materials are discussed.Model and Method.-Weconsider the 2D spin-1/2 LR Heisenberg model with power-law decaying couplings on the square lattice.The Hamiltonians for the ferromagnetic and the antiferromagnetic cases (with staggered interaction to avoid the sign problem [102,103]) are given by The schematic spin configurations and the decaying LR interactions of J(r) = 1/r α for both cases are shown in Fig. 1(a-c).Here we set J = 1 and simulate the system sizes upto L = 64, the inverse temperature β = L/2 and the decay exponent α from 1.5 to 100, with the focus on α ≤ 6.We note that to probe the ground-state properties in finite space-time size QMC simulations, one usually scale β ∼ L z where z is the dynamical exponents.In our problem, z is the largest at the short-ranged cases (for example, z = 1 for short-range antiferromagnetic case) and our results are well converged.Since the long-range cases will only have smaller z, β = L/2 is more than sufficient to extract the ground state properties therein.Detailed implementation and finite size analysis of the obtained dispersions in QMC are shown in Sec.II of Supplementary Materials (SM) [104].
As discussed in Refs.[85,105], for H F M , the SWT analysis accompanied by a continuum approximation concludes that for α > d + 2 (denoted as standard Goldstone regime) where d is the spatial dimension, the lowmomentum dispersion of the LR model reduces to the short-range case with ω ∼ |q| 2 , and for d < α < d + 2 (denoted as anomalous Goldstone regime) the dispersion is ω ∼ |q| α−d .For α ≤ d (denoted as Higgs regime) the system becomes gapped because of the generalized Higgs mechanism [85].As will be shown below, our QMC results are consistent with this picture as we reveal three different regimes via fitting ω ∼ |q| s(α) and finite size analysis.When α is large, the system is in the standard Goldstone region with dispersion ω ∼ |q| 2 .As α decreases, LR interaction brings the system into the anomalous Goldstone regime, and we find exactly the same Higgs regime as in Ref. [85] which is α ≤ 2 (see Fig. 1 (d)) As for the antiferromagnetic case, it is worth noting that Ref. [85]   gapped spectra, and the anomalous and standard Goldstone regimes are d − 2 < α < d and α > d respectively.However, we consider a sign-problem-free Hamiltonian of Eq. (2) which does not host frustrations and we get different boundaries of the three regimes.The system returns to the standard Goldstone mode when α is large enough, and our QMC results show the anomalous Goldstone regime is α > 2.2 and the Higgs regime is α ≤ 2.2 (see Fig. 1 (e)).
It is interesting to see both H F M and H AF M are superextensive when α ≤ 2 and it can be seen from Fig. 2 (a) that the gap diverges with system sizes at α = 1.5 in the ferromagnetic case.Between 2 < α < 2.2 for H AF M , the Hamiltonian is extensive and yet the spectra are gapped (as shown in Fig. 3 (a) and (b)).Therefore, our results go beyond the conventional wisdom for long-range systems [88,[99][100][101], and suggesed the disentanglement of the gapped spectrum and the superextensiveness of the Hamiltonian.We further considered the Kac construction which couples a normalization factor (N − 1)/ i =j 1 |ri−rj | α to Eqs. (1) and perform SWT analysis on the normalized Hamiltonian and confirm the gapped spectra with the same boundaries.SWT data are presented in Sec.I in SM [104].
Results.-Fig. 2 shows the obtained QMC spectra along the high symmetry path Γ(0, 0) → X(π, 0) → M (π, π) → Γ(0, 0) for the ferromagnetic case.At α = 100 (panel (f)), the system reduces to the short-ranged case with only nearest-neighbor couplings [10][11][12], and our QMCobtained spectra matches well with the SWT spectra.Both of them show a ω q ∼ |q| 2 dispersion close to Γ, and they match well along the whole path.As α gets smaller, as shown in panels (c), (d) and (e), we find the dispersion enters the anomalous Goldstone region [85], i.e., the dispersion close to Γ deviates from a quadratic one.We use ω q ∼ |q| s(α) to fit the dispersion close to Γ and find the power s(α) gradually decreases as α gets smaller.Insets of these three panels demonstrate the power law fitting of s(α) using L = 64 QMC data.We find, at α = 3 (panel (d)), s = 1.076 which agrees well with the relation of s(α) = α − 2 suggested in Ref. [85].However, for α = 4 (panel (e)) and α = 2.5 (panel c)) our results show apparent derivations from s(α) = α − 2. Fig. 1(d) collects the fitted power s(α) by QMC (red dots) at various α and we observe a satisfactory match with our SWT results.At α = 2, we find ω q near Γ, i.e., q = ( 2π L , 0) for different system sizes converge to a large and finite value of ω ≈ 3.42 as indicated in the inset of Fig. 2 (b).This phenomenon is fundamentally different from a gapless excitation in which the finite size gap ω (2π/L,0) converges to zero as L → ∞ and results in a continuous spectra.Our result reveals that at α = 2 the system enters the Higgs regime where the Goldstone mode acquires mass due to the LR interation and the excitation spectrum becomes gapped.For α < 2 (α = 1.5 in panel (a)) we find the gaps begin to diverge with the system size L due to the aforementioned super-extensive Hamiltonian.Therefore, we conclude that α = 2 is the separation power between the Higgs-type and Goldstonetype spectra in H F M from our QMC results.
Fig. 3 illustrates the QMC dispersion relation for H AF M along the high symmetry path.Similarly, in panel (f), we benchmark the spectrum at α = 100 with SWT result for the short-range antiferromagnetic Hamiltonian (with an extra coefficient ∼1.158 multiplied to approximate the second order spin wave effects [76,98,106,107]) and find QMC results agree well with SWT dispersion close to M with ω q ∼ |q|.As α decreases, the system also enters the anomalous Goldstone region with ω q ∼ |q| s(α) and 0 < s(α) < 1 close to M .Fitted powers via QMC at various α are displayed in Fig. 1(e) and agree well with the SWT results.In Fig. 3 (b) at α = 2.2, ω q close to M converges to a large and finite value of ω = 4.57.This means H AF M is in the Higgs-regime with gapped spectra when α ≤ 2.2.And in fact, panel (a) with α = 2.1 clearly forecasts a non-vanishing gap close to M at thermaldynamic limit.Interestingly, performing similar analysis as done in Ref. [85], one would obtain the same Higgs boundary as in the ferromagnetic case, i.e. α = 2 which deviates from the boundary at α = 2.2 obtained from unbiased QMC simulations.The deviation seems to suggest that antiferromagnetic quantum fluctuation plays a non-negligible role in this gap generating process and makes it different from the ferromagnetic case.Therefore, it is of interest to conduct further theoretical analysis for Eq. ( 2) to understand this discrepancy and reveal the subtle working of the disentanglement of the gapped spectrum and the super-extensiveness of the Hamiltonian.
Discussions.-With the unbiased large-scale QMC simulations and SWT analysis, we systematically investigate the dynamical properties of 2D spin-1/2 Heisenberg model with LR interactions.We find in contrast to the well accepted low-energy customs such as Hohenberg-Mermin-Wagner theorem and gapless Goldstone mode, the LR quantum many-body systems offer richer tun-ability and exhibit new phenomena.As the interaction exponent α varies, the Goldstone modes can be strongly modified, in that they can be either distorted (in the anomalous Goldstone regime), or even be gapped via a generalized Higgs mechanism.We also find explicit case where the gapped excitation exists even when the Hamiltonian is extensive.
These dynamical properties have immediate relevance to the ongoing experiments with ultracold atom arrays and quantum moiré materials.For example, the longrange Coulomb interaction in quantum moiré systems can be easily tuned by varying dielectric environment, electrostatic gating and twisting angles, and in this way observed thermodynamical and dynamical properties (such as switching between gapped and gapless spectra) [34,59,60,65,66,68] can be identified with different LR interaction types and regimes, when compared unbiased results such as ours.Similar tunability can also be realised in dressed Rydberg atom arrays whose interaction can be modified [108], one can then compare different responses from experiments with our results to identify the LR interaction and the novel phases.

Supplementary Materials
In this supplementary material, we present the linear spin wave analysis for the LR FM and staggered AFM Heisenberg model, in which the dispersion relation of the low-lying magnetic excitations at different decaying power α are extracted.From here, we make comparison with the dispersion obtained from the QMC simulations in the main text.Spin wave analysis with Kac construction is presented.Moreover, we outline the QMC procedure and provide detailed data on the fitting of the excitation gaps from the dynamical correlation functions obtained in QMC simulations.

Linear spin wave analysis
We applied the linear spin wave theory (SWT) to analyze the dispersion of the low energy excitation in the LR spin-1/2 Heisenberg model with power law decaying couplings in both ferromagnetic and staggered antiferromagnetic cases [12,[94][95][96][97][98].Taking the staggered antiferromagnetic cases as an example, it calls for the definition of two sublattices, A and B. The spin on each sublattice is pointing in the same direction.Then, we rewrite the spin operators by S + = S x + iS y and S − = S x − iS y and apply the Holstein-Primakoff transformation up to order S that for sublattice and for sublattice B Here we take S = 1/2 and the Hamiltonian in the momentum space is given by in which γ † (q) = (a † q , b q ) and a † q is the Fourier transformed that a † q = N 1/2 q a † i e −iqr .And, J s q = r s ∈same e −iqr s J s r refers to the coupling between the spins belong to the same sublattice, and r to that of the different sublattices.J d(s) r = 1/|∆r| α is the coupling strength.Finally, the single magnon dispersion relation of the LR Heisenberg model with staggered antiferromagnetic power-law decaying couplings is given by Similarity, in the ferromagnetic case, the dispersion relation of single magnon can be read as ω FM q = |J 0 − J q | with J q = r e −iqr J r .
To capture the dependence of the dispersion relation to α in Eqs. ( 1) and ( 2) in the main text, we numerically calculate the linear SWT results by applying a cut-off of longest range coupling as r max , meaning that we only consider the coupling between the sites (r x , r y ) and (r x + ∆r x , r y + ∆r y ) with ∆r x and ∆r y ranging from −r max /2 to r max /2, and we have computed r max up to 1000.
Fig. 4 (a) describes the dispersion relation of the linear SWT along the momentum path Γ → X → M → Γ with α changing from α = 2.0 to 4.0.For large α, the LR coupling rapidly decays which makes its disperion relation of single magnon similar to that of the typical antiferromagnetic square lattice only with nearest-neighbor coupling.Here, the single magnon dispersion relation is gapless only at Γ and M .As α decreases, the LR couplings strongly distort the dispersion relation, which makes the magnon excitation cost more energy and the dispersion relation goes higher as α decreasing.But this dispersion relation obtained from the linear SWT theory still remains gapless at Γ and M .Actually, the single magnon dispersion relation would become discrete at Γ and M in the limit r max → ∞.
With ∆q the relative momentum away from the M point, we plot the dispersion relation along the M → Γ direction in Fig. 4 (b) with a double logarithmic scale at α = 3.5 with varying r max , which shows the power-law dependence between ∆q and ω.We fit the linear SWT results with ω q = A|∆q| s near the M point.Note that for small r max (r max ≤ 100), the single magnon dispersion around the M point still depends on r max , which is shown in Fig. 4 (b) with r max changing from 16 to 1000.Such a dependence would disappear and s would finally converge in the limit r max → ∞.In order to demonstrate this convergence process as r max → ∞, we plot two black dashed lines that ω q ∝ |∆q| 0.48 in (b), where s = 0.48 comes from the fitting of the linear SWT result with r max = 1000.In Fig. 4 (b), as r max increases, s converges to 0.48.Finally, with r max = 1000, Fig. 4 (c) presents our fitting about the relation between the power s(α) and the decay exponent α, which suggests the limit r max → ∞ and is also plotted as the blue dots in Fig. 1 (d) in the main text.
Similarly, we also plot our linear SWT results of the ferromagnetic case in Fig. 4 (d), (e), and (f).Fig. 4(f) also shows the relation between the power s(α) and α taking r max = 1000, which is also given as the blue dots in Fig. 1 (e) in the main text.where H|n = E n |n and E 0 is the ground state energy of the system.When β∆E 1 1 where ∆E n = E n − E 0 , we can estimate G q (τ ) ≈ n=1 | 0|S z q |n | 2 e −∆En(q)τ + e −∆En(q)(β−τ ) .
(12) When the imaginary time is sufficiently large, we assume that the system will gradually evolve to the ground state, so that the correlation function can be further approximated by If | 0|S z q |1 | 2 is finite (which is usually the case), we can then extract the energy gap for each q point by fitting G q (τ ) with an exponentially decaying function.
We fit the QMC data of G q (τ ) by the relation G q (τ ) ∝ e −∆qτ and the fitting process is shown in Fig. 6.We first choose the data points for fitting according to their relative errors.If the relative error of one data point is less than 0.2, then the data point is chosen to be used for fitting.In the fitting process, we gradually omit the first N τ data points and then do the curve fitting to find the most probable gap.As shown in the inset of Fig. 6, the fitting error becomes intolerant when N τ = 10 and the fitted gap converges around ∆ = 2.35 when N τ gradually decreases to 0. In this case, we choose ∆ = 2.35 to be the fitted gap for the data.Note that we find that for all the q points at different α, the fitted gap does not change evidently with N τ , which means that higher excited states have much bigger energy gaps then the first excited states(∆E 2 ∆E 1 ) so that e −∆E1τ term in G q (τ ) contributes much more than other terms for the range of τ we consider. .Plotting under double logarithm scale, it is demonstrated that the power of low-momentum dispersion s(α) depends on the system size L.However, as the system size L increases, the dispersion gradually converge and s(α) will finally remain unchanged as L → ∞.Such process can already be seen as the vast majority of L = 56 and L = 64 data collapse onto the same curve in both FM and AFM cases, we thus obtain s(α) by fitting L = 64 data and y ∝ x s(α) is plotted as red lines in Fig.

FIG. 1 .
FIG. 1. 2D LR Heisenberg models with Higgs, anomalous and short-range Goldstone spectra.The schematic plots of 2D Heisenberg model with LR ferromagnetic interaction (a) and staggered antiferromagnetic interaction (b).The power-law decay of J(r) ∼ 1/r α for three different α is shown in (c).(d) and (e) show the power s(α) of low-energy spectra ω ∼ |q| s(α) obtained from QMC and SWT versus α for both the ferromagnetic and antiferromagnetic cases.The greenshaded area represents the Higgs regime where the spectra are gapped, the yellow shaded area represents the anomalous Goldstone regime where the dispersion powers change with α, and the white area is the standard short-range Goldstone regime where s = 1 for antiferrmagnetic and s = 2 for ferromagnetic cases.The red dots are fitting results from QMC (L = 64) and the red stars denote QMC boundaries at α = 2 (for ferromagnetic) and α = 2.2 (for antiferromagnetic) which speparate the Higgs and Goldstone regimes.The blue dashed lines are fitting results from SWT, with a cut-off of longest coupling distance rmax = 1000.
predicts the Higgs regime occurs at α ≤ d−2 for the Hamiltonian H = J i =j 1 |ri−rj | α S i •S j .Therefore for d = 2 there will be no finite α values with

FIG. 2 .
FIG. 2. Dynamical properties of 2D LR ferromagnetic Heisenberg model.Dispersion relations along the path (Γ → X → M → Γ) with panels (a)-(f) for different decay exponents α. Results of various sizes L are plotted together in each panel and share the same legend on the top.Inset of panel (b) indicates that at α = 2 the first excitation gaps near Γ for various sizes converge to a finite value and the system has a gaped spectrum, i.e., inside the Higgs-regime.Insets of (c), (d) and (e) show the fitting of power-law dispersions ωq ∼ |q| s(α) near Γ (with |∆q| denotes the relative momentum away from Γ) in range 2 < α ≤ 4. Red dashed line in (f) is the SWT dispersion for 2D nearest neighbor FM Heisenberg model with ωq = |J| zS(1 − γq) where S = 1/2, the coordination number z = 4, and γq = 1 z δ e iqδ .

FIG. 3 .
FIG. 3. Dynamical properties of 2D LR (staggered) antiferromagnetic Heisenberg model.Dispersion relations along the path (Γ → X → M → Γ) with panels (a)-(f) for different decay exponents α. Results of various sizes L are plotted together in each panel and share the same legend on the top.Inset of panel (b) indicates that the first excitation gaps near M = (π, π) for various sizes converge to a finite value and thus the system is inside the Higgs regime at α ≤ 2.2.Insets of (c), (d) and (e) show the fitting of power-law dispersion as ω ∼ |q| s(α) near M (with |∆q| denotes the relative momentum away from M ) in range 2.2 < α < 4.5.Dashed red line in (f) is the nearest-neighbor SWT dispersion ωq = |J| zS (1 − γq)2 with an additional coefficient ∼1.158 to approximate the second order spin wave effects[76,98,106,107].

7 FMFIG. 5 .FIG. 6 .
FIG. 5.The linear SWT results considering the Kac normalization.Panel (a) is for the the staggered antiferromagnetic coupling and describes the spin wave dispersion relation near the M point while (b) stands for the ferromagnetic case and near the Γ point.|∆q| = 2π/rmax with rmax ranging from 16 to 1000.The blue, green and red line refer to α = 1.5, 1.6, 1.7.The dashed line is the energy gap ∆ observed by fitting ω(∆q) with a|q| b + ∆.

Fig. 7 (
Fig.7(a) shows the dispersion of H F M near Γ and Fig.7(b)  shows the dispersion of H AF M near M for various system sizes L at α = 3. |∆q| denotes the relative momentum away from the Γ in (a) (M in (b)).Plotting under double logarithm scale, it is demonstrated that the power of low-momentum dispersion s(α) depends on the system size L.However, as the system size L increases, the dispersion gradually converge and s(α) will finally remain unchanged as L → ∞.Such process can already be seen as the vast majority of L = 56 and L = 64 data collapse onto the same curve in both FM and AFM cases, we thus obtain s(α) by fitting L = 64 data and y ∝ x s(α) is plotted as red lines in Fig.7.The same red lines are shown in the insets (c), (d) and (e) in Figs.2 and 3in the main text.

7
Fig.7(a) shows the dispersion of H F M near Γ and Fig.7(b)  shows the dispersion of H AF M near M for various system sizes L at α = 3. |∆q| denotes the relative momentum away from the Γ in (a) (M in (b)).Plotting under double logarithm scale, it is demonstrated that the power of low-momentum dispersion s(α) depends on the system size L.However, as the system size L increases, the dispersion gradually converge and s(α) will finally remain unchanged as L → ∞.Such process can already be seen as the vast majority of L = 56 and L = 64 data collapse onto the same curve in both FM and AFM cases, we thus obtain s(α) by fitting L = 64 data and y ∝ x s(α) is plotted as red lines in Fig.7.The same red lines are shown in the insets (c), (d) and (e) in Figs.2 and 3in the main text.
Fig.7(a) shows the dispersion of H F M near Γ and Fig.7(b)  shows the dispersion of H AF M near M for various system sizes L at α = 3. |∆q| denotes the relative momentum away from the Γ in (a) (M in (b)).Plotting under double logarithm scale, it is demonstrated that the power of low-momentum dispersion s(α) depends on the system size L.However, as the system size L increases, the dispersion gradually converge and s(α) will finally remain unchanged as L → ∞.Such process can already be seen as the vast majority of L = 56 and L = 64 data collapse onto the same curve in both FM and AFM cases, we thus obtain s(α) by fitting L = 64 data and y ∝ x s(α) is plotted as red lines in Fig.7.The same red lines are shown in the insets (c), (d) and (e) in Figs.2 and 3in the main text.

FIG. 7 .
FIG. 7. Dispersion relation at α = 3 for various system sizes L. (a) Dispersion of HF M near Γ with ∆q denotes the relative momentum away from Γ. Red line y ∝ x 1.076 shows the fitted power s(α = 3) = 1.076 using L = 64 QMC data.(b) Dispersion of HAF M near M with ∆q denotes the relative momentum away from M .Red line y ∝ x 0.469 shows the fitted power s(α = 3) = 0.469 using L = 64 QMC data.