{'fraction_non_alphanumeric': 0.04714713191662018, 'fraction_numerical': 0.018623522848303493, 'mean_word_length': 4.2584809046298275, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Collective and single-particle excitations in 2D dipolar Bose gases

The Berezinskii-Kosterlitz-Thouless transition in 2D dipolar systems has been studied recently by path integral Monte Carlo (PIMC) simulations [A. Filinov et al., PRL 105, 070401 (2010)]. Here, we complement this analysis and study temperature-coupling strength dependence of the density (particle-hole) and single-particle (SP) excitation spectra both in superfluid and normal phases. The dynamic structure factor, S(q,omega), of the longitudinal excitations is rigorously reconstructed with full information on damping. The SP spectral function, A(q,omega), is worked out from the one-particle Matsubara Green's function. A stochastic optimization method is applied for reconstruction from imaginary times. In the superfluid regime sharp energy resonances are observed both in the density and SP excitations. The involved hybridization of both spectra is discussed. In contrast, in the normal phase, when there is no coupling, the density modes, beyond acoustic phonons, are significantly damped. Our results generalize previous zero temperature analyses based on variational many-body wavefunctions [F. Mazzanti et al., PRL 102, 110405 (2009), D. Hufnagl et al., PRL 107, 065303 (2011)], where the underlying physics of the excitation spectrum and the role of the condensate has not been addressed.


I. INTRODUCTION
Dipolar bosonic systems are of increasing interest for recent experiments studying the onset of superfluidity in nonideal Bose systems and its connection with correlation and quantum degeneracy effects. Examples include dipolar gases, as in recent studies of ultra-cold dipolar gases, 1-6 indirect excitons in coupled quantum wells 7-10 in external electric fields or exciton-polaritons in quantum wells embedded in optical micro-cavities. 11,12 Due to significant experimental achievements, a superfluid regime of polar molecules is expected to be observed in the near future. [14][15][16][17][18][19][20] There is also an increasing activity on alkali atom gases where the long-long range dipolar interactions can be generated by excitations to high energy atomic Rydberg states (one electron is excited to a very high principal quantum number). The applied moderate electric field can result in a large polarizability and dipole moment. 21,22 All these achievements have a direct impact on understanding the properties of Bose-condensates dominated by long-range correlation effects.
From the theoretical side, the static properties of dipolar bosons have been quite extensively analyzed in recent years. Strongly correlated phases of dipolar bosons confined in a harmonic potential have been studied from first principle simulations by Lozovik et al., 23 Nho et al., 24 Pupillo et al., 21 Golomedov et al. 25 and Jain et al. 26 2D homogeneous systems have been analyzed by Astrakharchik et al. 27 and Büchler et al. 28 with the prediction that above a critical density a superfluid dipolar gas undergoes a crystallization transition. Properties of a dipole solid have been addressed by Kurbakov et al. 29 The presence of a finite fraction of vacancies and interstitials, in incommensurate crystal, leads to a superfluid response and quasi-equilibrium supersolid phase. In PIMC simulations by Filinov et al. 30 the density-dependence of the transition temperature from a superfluid to a normal gas phase has been analyzed. The composite dipolar bosons, such as indirect excitons (pairs of electrons and holes spatially sep-arated in semiconductor bilayers) have been recently studied by Böning et al. 31 with focus on the effective exciton-exciton interaction and its consequence for the complete phase diagram of excitonic system. Relevant experimental parameters to observe exciton crystallization, in a semiconductor quantum well, have been predicted. All these results provides a reliable estimate of the critical temperature and degeneracy parameters required to observe a combined effect of interparticle correlations and Bose statistics in experiment.
On the other side, theoretical predictions for dynamic properties, in particular excitation spectra of density fluctuations, are more challenging if performed on a microscopic level. For 3D dipolar superfluids at low densities the rotonization of the excitation spectrum has been predicted from mean-field analyses. [32][33][34][35] A cigar-shaped trap is a promising candidate for observation of a roton minimum, where the dipole-dipole interaction is anisotropic, and partially attractive. In a quasi 2D trap (a pancake geometry), a roton has a different physical nature and is due to a strong in-plane repulsion of similar oriented dipoles. This regime assumes a high density and goes beyond the mean-field predictions. This calls for a more accurate treatment. The correlated basis functions theory (CBF), originally developed for strongly correlated systems like 4 He, has been renewed recently by Mazzanti et al. 36 to analyze in detail the phonon-roton dispersion of an infinite 2D dipolar gas at T = 0. The density range has been identified, where the CBF dispersion (which includes three-phonon interactions) deviates from the Bijl-Feynman 27,37 and Bogolubov predictions. The upper bound for the phonon-maxon-roton dispersion, derived from the frequency-sum rules, [39][40][41] is mainly in a good agreement. 30 Important extension of the CBF theory to 3D geometry and full anisotropy of the dipole-dipole interaction has been performed by Hufnagl et al. 42 This study builds an important bridge between the nature of two kinds of rotons caused either by attractive or repulsive tail of the dipole interaction. Still the current analyses of the excitation spectrum lack the important regime of finite temperatures (in-cluding T ∼ T c ) and treat excitation damping in an approximate way. [43][44][45] With the present analysis we aim to fill this gap and, in addition, discuss the nature of excitations in Bosecondensed systems and the role played by a condensate.
We start the discussion about the nature of excitations from the model developed for a one-component neutral BEC.
The properties of the density fluctuation spectrum have been discussed extensively in relation with superfluid liquid 4 He. In particular, the temperature (T ) dependence of S(q, ω) from inelastic-neutron scattering experiments was compared vs. Glide-Griffin (GG) model based on the dielectric function formalism. 46,47 The relation with the dynamic structure factor is written as χ(q, ω) = Λ(q, ω)G 1 (q, ω)Λ(q, ω) +χ ρ (q, ω) whereχ is the irreducible part of dynamic susceptibility χ, separated into a singular condensate part (related to singleparticle excitations) and a regular thermal (multiparticle) component, present also above T c . The coupling to the SP spectrum is described by the vertex function Λ ∝ n 0 (T ) and the single-particle Green function G 1 . The first term in (3) vanishes in the normal phase (T > T c ), with the resultχ =χ ρ , withχ ρ being the density response of weakly interacting particle-hole excitations, e.g. treated via a Lindhard function. The excitation spectrum described by χ(q, ω) is expected to include for small q a zero sound and acoustic phonons, for intermediate q a maxon-roton branch, and for large q a combination of sharp SP resonances and a broad distribution close to the recoil energy q = q 2 /2m due to the multiple scattering of weakly interacting particle-hole excitations.
To fit the experimental observations for 4 He by Eqs.(1)-(3), i.e. a single phonon-maxon dispersion, one suggests that the imaginary parts of χ and G 1 share a common denominator due to the vertex Λ(q, ω). 48 As a result the SP and collective excitations can not be distinguished in S(q, ω) in the superfluid phase. For T > T c the SP spectrum decouples and only the thermal part remains (χ =χ ρ ). The dispersion gets broadened and remains well defined only in the phonon (low-q) part of the spectrum.
In its simple form the dielectric function model failed to predict the double roton feature observed in 4 He. Therefore, it was necessary to include the mechanism of the quasiparticledecay processes. This has lead to development of theories beyond the roton minimum. [51][52][53][54][55][56][57] In his pioneer work, Pitaevskii 51 considered a semi-emperical ansatz for a singleparticle Green's function which includes the unhybridized Feynman-Cohen phononmaxon-roton branch, 58 ω q = ω FC q , and the self energy, Σ 12 = n 0 Jχ 2 J, expressed in terms of the two-(quasi)particle response function χ 2 and the three-point interaction vertex J.
The self energy accounts for the typical repulsion by anticrossing of two modes, i.e. the interaction between a single and a pair of quasiparticle excitations. As a result the SP dispersion, beyond the roton, is renormalized and bends to the energy of two roton minimum. This concept has been further successfully developed [52][53][54][55][56]59 by improving both ω FC q and Σ 12 .
In particular, Zawadowski et al. 53 developed the theory of a bound two-roton state. Due to hybridization with this state the dispersion splits into two branches. The upper branch consists of heavily damped density excitations with the energy close to the recoil energy q . The low-energy branch is due to the SP excitations and saturates slightly below the energy of twice the roton minimum. At low T with increasing q its weight is transfered to the upper branch. While Zawadowski's model provides a good description for low-temperature 4 He data by using several fit parameters, it becomes unsatisfactory at high temperatures. The expected decrease of the intensity of the lower (SP) branch, by vanishing of the condensate n 0 (T ), is not compensated by contribution of atoms above a condensate. This is in a contradiction with the sum rule, S L (q) + S H (q) = S(q), as the static structure factor S(q) for 4 He was found to be nearly temperature independent for T ∼ T c . 49 Further, in some range of q-wavevectors a fit to the T -dependent experimental S(q, ω)-data resulted in the repulsive interaction between two roton excitations, 60 opposite to the negative coupling constant suggested in the original model. This point, has been resolved by a more elaborated theory 55 using a T-matrix which allowed the coupling constant to oscillate with q.
The dielectric function models interpret the roton and the two-roton state as a specific feature of the SP spectrum coupled in S(q, ω) in a superfluid phase. The standard procedure to follow is to fit the involved key parameters (ω q , J, χ 2 ) of the Green's function (4) to be consistent with the low temperature S(q, ω)-data. The T -dependence is neglected and enters only via n 0 (T ) and the Bose function in Eq. (1). In a similar way, the density response functionχ ρ (3) is treated also as Tindependent, and, therefore, is fitted by high temperature data, where the singular condensate part does not contribute.
The weak point of this model is the discrepancy with experiments near T c . The intensity of the sharp peak (two-roton state) assigned with the SP excitations (poles of G 1 ) becomes too low as its spectral weight vanishes with n 0 (T ). In contrast, a typical experiment on 4 He demonstrates a continuous broadening of this peak and nearly conserved integrated spectral weight. Also a broad multi-excitation background, related withχ ρ , extends significantly to low frequencies, as T is increased, being in contradiction with the model assumption on its temperature independence. These discrepancies have been reviewed in Refs. [59][60][61][62][63] The experiments 61,62 provided no indication of a welldefined mode (SP excitations) suddenly appearing in S(q, ω) for T < T c due to the first term in Eq. (10). In contrast, they demonstrated that temperature and density dependence of the spectral density in the phonon and roton part of the spectrum can be fitted with just a single-mode susceptibility. In particular, the roton mode demonstrates a rapid attenuation when approaching T c from below and continues as an overdamped diffusive mode of zero frequency above T c . The authors concluded "our results cannot be explained using the GG model, unless of course the two components in the GG model hybridize into one (having one lifetime and excitation energy) at all temperatures and pressures, independent of the value of n 0 ". 64 This suggests, that the phonon-roton spectrum must be considered as unified branch as was done in the original papers by Landau, However the origin of the excitations and their structure can be different. By rigorous considerations, Nepomnyashchy 68 has demonstrated that a simultaneous presence in the spectrum of both zero sound (ZS) and SP branch is not possible. At T < T c the ZS is suppressed and replaced by the SP branch defined by poles of the SP Green function. The phonons become identical with the elementary excitations. The SP branch, however, accumulates some properties of the ZS branch, such as linear spectrum and independence on n 0 . This is in a principle contradiction with the ansatz (3) when one uses the substitution Λ ∝ n 0 (T ). Finally, we note that the Bogolubov type spectrum can be also derived without involving the gauge symmetry breaking, i.e. also for a normal phase with the conserved gauge symmetry. 69 This implies that the phonon-maxon-roton dispersion, observed in a superfluid, should be also preserved in the normal phase, once the damping is small.
For a 2D dipolar gas, we found that its properties near T c qualitatively follow that of 3D 4 He discussed above. From the comparison of S(q, ω) at T ∼ 0.7T c , we can confirm that the intensity of the two-roton "plateau" is unexpectedly too large and does not scale with the zero-temperature condensate fraction once the density is varied. The condensate demonstrates a broad variation 27 (n 0 = 0.70 . . . 0.025) and this implies, that for Λ ∝ n 0 (T ) in (3), the two-roton peak should be strongly suppressed at high densities, if we assume a roton being solely a SP excitation. In our simulations we do not find any evidence of this behavior. The intensity of the peak is not much reduced when approaching T c either. Formally, in 2D geometry, as considered here, at T = 0 the condensate fraction should vanish in the thermodynamic limit and no coupling to the SP spectrum is expected. Still the two-roton state is observed for T ≤ T c and even for T > T c at high densities (but strongly damped). The finite size effects are negligible for the considered system sizes (N ∼ 165, 576) and, therefore, do not influence the imaginary time correction functions used for the reconstruction of the spectral densities (Sec. V). All observed characteristic spectrum features should remain also in the thermodynamic limit (except the limit q → 0). This brings us to a preliminary conclusion that if the two-roton state is a feature of SP spectrum the sufficient condition for the hybridization in S(q, ω) is the presence of at least off-diagonal quasi-long range order and a non-zero superfluid density.
It is important to mention that, there exists a second viewpoint on roton and different combinations (roton+roton, ro-ton+maxon, maxon+maxon) as the intrinsic density excitations related with the short-range particle correlations. [70][71][72] To clarify these important issues one needs to independently access S(q, ω) and the SP spectral density A(q, ω) for temperatures below and above T c to indentify how the exci-tation branches of both spectra are coupled in the superfluid phase and lead to rotonization. This gives motivation for the present analysis. The path integral Monte Carlo simulations in combination with the stochastic optimization method are used to reconstruct collective and single-particle excitations at finite temperatures from first principles (in the linear response regime). Both should be coupled or form a unified branch in the density fluctuation spectrum S(q, ω) once a system has a condensate or the Bose broken symmetry. 73 Whether this prediction remains also valid for 2D Bose systems, with the off-diagonal quasi-long range order and only local Bose condensate, remains an open issue and calls for further analysis on the level of microscopic simulations.
Such analyses in application to 2D dipolar gas will be presented in Sec. VIII. In particular, we found that at weak and intermediate coupling (Sec. VIII A,VIII B) the single-particle and the collective (density fluctuations) modes have similar dispersion relations in the normal phase (when they are independent). As a result, in the superfluid phase due to the coupling (Eq. 3) both branches become hybridized. In general, the dispersion relation ω q observed in S(q, ω) splits into two branches, but the high energy (multiparticle) part gets significant broadening in the normal phase. At intermediate q (the maxon-roton region) the lower branch is shifted to lower energies and its damping is enhanced for T > T c .
For strongly correlated dipolar gas (D = 12.5) our results support the viewpoint on the roton as a collective density mode. While, the SP spectrum A(q, ω) shows a well pronounced roton minimum (denoted asẼ R ), the corresponding energiesẼ R ∼ 10 and 2Ẽ R ∼ 20 are larger than those of the roton-minimum (E R ∼ 5) and two-roton state (Pitaevskii "plateau" 2E R ∼ 10) observed in S(q, ω). According to the Zawadowski's model, this suggests a strong binding mechanism which currently lacks experimental confirmation in other strongly correlated Bose systems ( 4 He).
Independent on the presence of the two-roton state, the frequency sum rules (see Appendix B) predict a branch near the recoil energy q . Indeed, the simulations reproduce this branch, however, it is strongly damped, as expected from weakly interacting density quasiparticles. This multiparticle feature is present both below and above T c . In contrast to 4 He, in a dipolar system at high density (e.g. D = 12.5) the Pitaevskii "plateau" appears to be well separated from this multiparticle continuum. Similar free-particle like branch is observed in the SP spectrum A(q, ω).
Based on our results we interpret the two-roton state as a combination of two density rotons (the poles of the dynamic susceptibility) as we can observe it also for T > T c . The detailed discussion will be presented in Sec. VIII B,VIII C.
In Sec. II we introduce the model of a 2D dipolar gas. In Sec. III physical realizations are listed and experimental parameters are compared with our model. In Sec. IV the dynamic structure factor S(q, ω) is introduced for a system with Bose condensate. In Sec. V we analyze the structure of the imaginary time correlation functions and their relation with the SP spectral density A(q, ω) and S(q, ω). In Sec. VI reconstruction of A(q, ω) and S(q, ω) using the stochastic optimization procedure is discussed. In Sec. VII we analyze temperature/density dependence of the momentum distribution and the spectral densities. In Sec. VIII we give our interpretation of the different excitation branches observed in A(q, ω), S(q, ω). In Sec. IX we list our conclusions.

II. MODEL
Dipolar BEC has been first realized in atomic systems of magnetic dipoles. 1-6 New exciting results have been recently reported for dipolar molecules [14][15][16][17][18] with the major advantage of significantly larger electric dipole moments. The mean field description of such systems predicts that the dipolar condensates in 3D geometry become dynamically unstable when the dipole interaction dominates the s-wave scattering contact interaction. [32][33][34][35] This effect is a consequence of the rotonization of the Bogolubov-type spectrum 32,33,74 and leads to an instability of the time-dependent Gross-Pitaevskii equation. 34 To avoid this problem and study the behavior of dipolar systems well beyond the mean-field regime with a stable condensate and/or a superfluid phase, we consider a quasi-2D geometry. The dipole moments are oriented along the direction of the strong trap confinement and produce only repulsive interaction p 2 /r 3 . The dipolar forces are assumed to completely dominate the contact interaction, therefore, the latter is neglected in the present model. We assume 2D spatial homogeneity by considering an infinite pancake trap geometry. The evaluated spectrum for the in-plane momentum q will be valid for the excitation wavelengths, λ q = 2π/q, bounded by the perpendicular L z and the in-plane L ρ system size, i.e L z < λ q < L ρ . Otherwise, the excitations acquire 3D character. For a quasi-2D model of dipolar condensate, 33 the perpendicular extension can be estimated by L z ∼ l z = /mω z . The dipole moment p = ed and particle density n = 1/a 2 are free parameters controlled by external fields and in-plane trap frequencies. The system Hamiltonian can be made dimensionless by using the length and energy units: a = 1/ √ n (the mean interparticle distance) and In the canonical ensemble, the thermodynamic properties of (5) are completely defined by the dipole coupling D = p 2 / b a 3 E 0 , temperature T = k B T /E 0 and particle number N . The Bose broken-symmetry condition, i.e. Ψ (r) , Ψ + (r) = 0, and evaluation of imaginary-time Green's functions suggest to use the grand-canonical ensemble. In this case, the worm algorithm 75 has been employed by fixing a value of the chemical potential µ and system volume V instead of N (see note 76 ). Indeed, existence of the order parameter, Φ 0 (r) = Ψ (r) , is related with the fluctuations of a particle number in a condensate state. For a weakly interacting dipolar gas (D = 0.1), we observed significant particle number fluctuations, ∆N 2 ∼ 0.17 N , with a condensate fraction, n 0 ∼ 0.7. This can lead to some differences between the canonical ensemble (experimental BEC realized in traps) and the grand-canonical results.
The parameters of PIMC simulations are from Ref. 30 The dipole coupling (or density), D ∼ n 1/2 , was varied in the range 0.1 ≤ D ≤ 12.5 and temperature 0.7 ≤ T ≤ 3.3. With the known D-dependence of the critical temperature 30 (T c ≈ 0.70 . . . 1.4) we scan both the superfluid and the normal phase.

III. PHYSICAL REALIZATIONS
Next we discuss what D values are accessible in current experimental realizations. Consider the effective radius of dipole-dipole interaction, a d = mp 2 / b 2 , corresponding to the distance when the dipolar energy reaches the value of the zero-point kinetic energy, i.e. p 2 / b a 3 d = 2 /ma 2 d . The introduced dipolar parameter equals the ratio of two characteristic length scales, D = a d /a. Given a value of a magnetic/electric dipole moment (which enters in a d ) we can estimate the density required for a specific D. The effective dipole interaction radius can be expressed as where m is the mass in the unified atomic mass unit, p is the dipole moment in the debye. An external electric field aligns dipolar particles leading to a repulsive 1/r 3 -interaction. Several examples are considered below.
• Cold bosonic atoms with a permanent magnetic moment 1-3 in tight pancake traps. Magnetic dipoles are aligned by a magnetic field. A dipolar gas of 52 Cr with p = 6µ B and n ∼ 10 11 cm −2 (a ∼ 316 Å) has a d ∼ 24 Å. The dipolar interaction will dominate, when the s-wave scattering length (a s ∼ 100a B ) is suppressed to a s < a d using a high magnetic field Feshbach resonance. The resulted coupling, D ∼ 7.6 · 10 −2 , is close to the analyzed value D = 0.1. The superfluid transition temperature is estimated 30 to be T c ∼ 12µK (in a 2D geometry).
Very recently, a Bose-Einstein condensate of erbium ( 168 Er) has been reported, 6 being another promising candidate for experimental realization of strongly correlated dipolar quantum gases. Similar to dysprosium, erbium has a large magnetic moment p = 7µ B and atomic mass which enhance the dipole interaction radius (6). A large dipole moment also allows to reach a high efficiency both in a magneto-optical trapping and subsequent evaporative laser cooling to temperatures T ∼ 200mK.
• Polar molecules [16][17][18][19][20] 16 aligned by an electric field. For 41 (40) K 87 Rb gas with n ∼ 10 9 . . . 10 11 cm −2 and p = 0.57 Debye, we estimate a d ∼ 6118 Å, D ∼ 1.9 . . . 19, T c (D = 1.9) ∼ 0.053µK and T c (D = 15) ∼ 3.4µK. High densities and D-values can be achieved by localizing Feshbach molecules in an optical lattice to suppress their collisions and then quenching to a ground state via stimulated Raman adiabatic passage. 19 Besides molecule collisions, an upper bound for quantum degeneracy and the dipole coupling D is limited by chemical reactions and associated molecule loss rates. 13 • Rydberg-dressed atoms excited to high principal quantum numbers k and confined in two dimensions. 21,22 A large dipole moment can be reached in the linear Stark regime with the result p = (3/2)ea B k(k − 1) ∼ 1450 Debye for k = 20. We consider an ensemble of ground state atoms excited to high Rydberg states with the single-atom excitation probability P = (Ω/∆) 2 ∼ 0.01 (where Ω and ∆ are a Rabi-frequency and a detuning laser frequency). The effect of a local dipole blockade, 78 i.e. a zero probability to excite more than one atom within a distance R s , can be included via the reduced Hamiltonian of "superatoms", 78,79 with the av- s is a typical number of atoms (in 2D) within a superatom radius R s . The latter depends on the properties of the excitation laser (P, Ω, ∆), the gas density n and the interaction potential V (r). For estimation, we treat an excitation probability P (n, R s ) of a many-body system (modified by correlation effects and defined by P (n, R s )N s ∼ 1) as that of an uncorrelated gas, i.e. P (n, R s ) ∼ P ∼ 0.01. For a gas of rubidium Rydberg atoms ( 87 Rb) with the number density n ∼ 10 9 cm −2 we estimate N s ∼ 100 (R s ∼ 3.1µm), but only a half of the superatoms are interacting as they perform Rabi oscillations between their ground and excited states. This results in the superatom spacing a ∼ √ 2R s , and, correspondingly, a d ∼ 272µm and D = a d /a ∼ 61. For low enough temperature a Rydberg gas will undergo a crystallization transition, as the coupling is above the critical value D c ∼ 17. 27,28,30 Quantum effects between Rydberg-dressed atoms will be relevant for T < 1 nK, due to low densities.
These estimations show that the systems of ultracold dysprosium (erbium), polar molecules and indirect excitons are good candidates for experimental validations of the present analyses for 0.1 ≤ D ≤ 12.5. We assume that in an experiment a long-lived quantum gas is created and, hence, a thermodynamic description remains valid. Further, there should be no noticeable heating, e.g. due to breaking of molecular bonds or recombination of excitonic states. By this assumption, experimental parameters which correspond to different D-values are listed in Tab. I. Table I: Relevant parameters for 2D systems of the dysprosium 164 Dy (Dy), 41 K 87 Rb polar molecules (M) and indirect excitons (Ex). The assumed dipole moments are, correspondingly, p Dy = 10µB, p M = 0.6 Debye and p Ex = 960 Debye. The dipole coupling D is expressed in terms of the in-plane density n [cm −2 ]. For T ≤ Tc a system is in a superfluid phase. 30 The frequency ω0 ( ω0 ≡ E0) denotes a characteristic energy scale of collective and single-particle excitations.

IV. SINGLE PARTICLE AND DENSITY RESPONSE EXCITATIONS IN A BOSE LIQUID
Existence of a Bose-Einstein condensate is introduced by the Bose broken symmetry condition, i.e. by a finite value of the order parameter, â (+) 0 (r) = √ N 0 e iφ , where the phase φ can be set to zero for a homogeneous condensate. As introduced by Beliaev 80 the number density operator in momentum space can be decomposed aŝ The first two terms describe the scattering process which creates or destroys an atom with zero momentum. The last term includes all thermal atoms outside the condensate (k = 0). In the presence of a BEC the density operator is factorized which directly demonstrates the coupling of the singleparticle operator,Â q , and the operator of particle-hole excitations of thermal atoms,ρ q . The result of this coupling in the density fluctuation spectrum was written down by Hugenholtz and Pines. 73 The dynamic structure factor S(q, ω) describing propagation of the particle-hole excitations (9) decouples into three parts S(q, ω) = S A (q, ω) + S I (q, ω) +S(q, ω), where Correspondingly, S A describes the density fluctuations involving condensate atoms, S I is a result of the condensate-induced intermixing of single-particle and density excitations, andS defines the dynamic response of thermal atoms.
In the long wavelength limit and T = 0 Gavoret and Nozières 46 showed explicitly that the acoustic phonons are the poles of both S A (q, ω) and S(q, ω). At finite temperatures these analyses have been extended by Griffin and Glide 47 based on the dielectric function formalism. As a Bose liquid is cooled down to a critical temperature and a Bose condensate is formed (N 0 = 0), the single-particle excitations, i.e. the poles of S A (q, ω), get a finite weight in the density response spectrum. Independent on the representation (11)- (13), as discussed in Sec. I, the microscopic calculations 68 predict that S A (q, ω) and S(q, ω) have the same poles in the presence of the condensate and that the weight of the SP excitations should be independent on the absolute value of the condensate fraction.
Hence, the theory allows to interpret the sharp component observed beyond the maxon region in superfluid helium as a single-particle excitation branch, which vanishes in the normal phase with the condensate N 0 . This in turn allows to indirectly probe the single-particle spectrum of 4 He atoms by S(q, ω), as the latter can be measured experimentally by neutron-scattering. The development of this concept was stimulated by a set of high-resolution measurements. 81 Our goal is to verify this interpretation and identify the SP excitations in S(q, ω). The reconstructed spectra contain possible errors coming from the statistical noise in the evaluated imaginary time correlation functions and the stochastic reconstruction procedure discussed in Sec. VI.

V. QUANTUM CORRELATION FUNCTIONS
In linear response theory an external field produces a weak perturbation of a system in thermodynamic equilibrium. In this regime dynamical properties can be evaluated via timecorrelation functions (TCF) of the corresponding dynamical operators. For classical many-body systems, TCF can be efficiently computed in the canonical ensemble by molecular dynamics (MD) simulations. A brute-force solution of classical equations of motion nowadays is possible for large ensembles (N ∼ 10 6 ) of particles. This allows to make accurate extrapolation to the thermodynamic limit, study in detail topological phase transitions 82 (e.g. BKT in 2D) involving large-scale density fluctuations and divergence of the correlation length near a critical point.
In contrast, for quantum systems first principle timedependent simulations can be done only for few-particle atomic and molecular systems, e.g. by time-dependent multi-configurational Hartree Fock and CI methods. 83 For condensed matter systems the continuous-time Monte Carlo method proved to be an efficient and powerful technique. 84 While predictions from time-dependent methods for equilibrium can be accurately checked by quantum Monte Carlo (QMC) methods for ground state and finite temperatures, the resolved dynamic correlations are more difficult to control for self-consistency. In particular for Bose-condensed systems, the importance to accurately treat dynamic correlations between thermal and condensate atoms has been discussed recently by Griffin, Nikuni and Zaremba. 85 Here, we follow the approch based on quantum correlation functions evaluated via path integral Monte Carlo.
The correlation function of single-particle operators in Eq. (10) (for real time t > 0) contains the contributions of two normal and two anomalous Green's functions The poles of G nm are given by zeros of the common denominator as follows from the formal solution of the Dyson-Beliaev equation in terms of the self-energies. 73,86 Hence, characteristic single-particle excitations can be extracted from one (normal) Green's function G 11 as other spectral densities A nm share the common poles. Therefore, we can concentrate on comparison between A(q, ω) and S(q, ω), and their possible hybridization in a Bose condensed phase. The general problem with the time-dependent ensemble averages using stochastic methods, like QMC, is the evaluation Figure 2: The same as in Fig. 1 but for the density-density correlation function ρ(r2, t2)ρ(r1, t1) (in units a 4 ).
of the high-dimensional integrals from rapidly oscillating exponentials in real-time propagators. This well known "dynamical sign problem" 84 is the main obstacle for first principle evaluation of TCF in Eq. (15) in the framework of the particle orbital-free continuous phase-space methods. Another possibility to compute the spectral function is the analytical continuation of an imaginary time correlation function to real frequencies. 86 The imaginary-time formulation of QMC (e.g. DMC, PIMC, GFPIMC) is the best suited tool. In particular, the Matsubara Green's function can be efficiently evaluated via PIMC in the grand canonical ensemble on a set of imaginary time points τ ∈ (0, β = 1/k B T ]. Below we present several examples. Fig. 1 (Fig. 2) shows the changes in the structure of the Matsubara Green's function (density correlation function) with the increase of the dipole coupling D on the spatial and imaginary time axis. The simulation parameters correspond to a superfluid (the superfluid fraction is ρ s /ρ > 0.80). In both cases, we observe that an enhancement of the inter-particle interaction leads to a more complicated structure with oscillations. These oscillations decay both in space and time, which is a general feature related to damping of one-particle or density (collective) excitations of a specific wave-length and energy. As we will see later, the interaction also leads to additional high-energy excitation branches which are absent in the weakly coupled regime, e.g. D = 0.1.
The behavior of the spatially isotropic function G 1 (r, τ ), with r = |r 2 − r 1 |, can be easily understood in two limiting cases. As τ → 0 and r → 0, the many-body effects have only a small influence on the one-particle propagator, which is then close to the free-particle density matrix, e −(r2−r1) 2 /2λ 2 τ with λ 2 τ = 2 τ /m. Indeed, in Fig. 1 near the origin we observe a similar Gaussian-shaped peak irrespective on the coupling D. After few oscillations in real space it Table II: The D and T -dependence of the average particle number N and the zero momentum occupation Nq=0 in the volume V1 = 165. Simulation parameters: T1 = 1.0, T2 = 0.714 and chemical potential µ1. The critical temperature Tc is from Ref. 30 The statistical error is given in the brakets or if not specified is in the last significant digit.  evolves into a flat distribution. Another recognizable behavior is recovered as τ → β. The Green's function G 1 (r, β) becomes the one-particle density with two well-known limits, i.e. G 1 (r = 0, β) = n, where n = N /V is the average particle density, and G 1 (r = L/2, β) being an upper bound for the condensate density n 0 = N q=0 /V which decays with a power law with system size (in 2D and T = 0). Both values depend on µ, V, T and can be read out from Tab. II.
The listed values explain the observed off-set in the region where G 1 (r, τ ) is almost flat. For temperatures above the critical, T c , the off-diagonal quasi-long range order is lost and G 1 (r, τ ) decays exponentially with r. With our choice 76 of the chemical potential µ(D) (see Tab. II) the particle density n and the Green's function G 1 (0, β) = n take a similar value, i.e. n ≈ 1, independent on D.
The behavior of the density correlator ρ(r 2 , t 2 )ρ(r 1 , t 1 ) (Fig. 2) has one simple limit. At τ = 0 it is a superposition of a δ-function (at r = 0) and of the pair distribution function. At finite but small τ the δ-function turns into the free-particle like density matrix, which at strong coupling (large D) gets more localized due to the inter-particle interactions.
Next we consider the reconstruction problem of spectral density from the imaginary-time.

VI. STOCHASTIC OPTIMIZATION METHOD
In this section we present the stochastic optimization method for reconstruction of the spectral densities. The method is free of difficulties involved in the analytic continuation of imaginary time correlation functions. As application, the dispersion relations for a 2D dipolar Bose system will be presented in Sec VII.
We start from the general relation between the spectral density and the single-particle/density-density propagator Fourier transformed to the momentum-space The spectral densities satisfy two normalization conditions The inversion of equations similar to (17)-(18) is known to be an ill-conditioned problem and results in a non-uniqueness of solution. By Monte Carlo simulations the reconstructed spectral densities get affected by the statistical errors present in G n (q, τ ). The standard tool used to partially overcome this problem is the Maximum Entropy (ME) method. 87 However, the reconstructed spectral densities for Bose liquids appear to be too smooth and the important information on a sharp δlike quasi-particle feature present in the excitation spectra is typically lost. 88 Recently, the Stochastic Optimization (SO) method has been introduced by Mishchenko et al., 89 which allows to overcome this difficulty. Its main advantages compared to the standard regularization methods, like the ME, are the continuous parametrization in frequency space, instead of a predefined mesh, and non-suppression of high derivatives of the spectral function, performed by the regularization methods. As a result, sharp peaks and edges are not lost during the reconstruction. In its core, with the SO one solves by stochastic sampling the minimization problem of the least deviation where δG(q, τ i ) is the statistical error of G(q, τ i ) at the imag- andG n is generated from withS n (q, ω) being a trial spectral density parameterized into some basis set. The deviation D n [G n ] is optimized by a random sequence of updates which can change both the parameters of the basis set and its size. By increasing n (the number of independent solutions) we evaluate the corresponding variance σ 2 Dn = D 2 n (with the zero mean). In the end, we select only "good" solutions from the whole sample which satisfy the condition D n ≤ D min with D min = 1.5 σ Dn . The final solution is constructed from their linear combination (100 ≤ M ≤ 400) to take advantage of a self-averaging of the noise. For the parametrization of trialG n we use a set of rectangles 89 with height h, width w and center c being the optimization parameters for a fixed value of the q-vector. The basis size N was also varied during the optimization in the range 80 N 400.
For the reconstruction of the dynamic structure factor S(q, ω) we have used only positive frequencies, 0 < ω < ω CO , with the cut-off frequency ω CO ≈ 400, by taking into account explicitly the relation between negative and positive energy transfers, i.e. S(−q, −ω) = e −βω S(q, ω). This results iñ In the basis of the rectangular functions the trial imaginary time density correlation function takes the form which is symmetric with respect to the mid-point t = β/2. Therefore, it is sufficient to evaluate G 2 (q, τ ) (see Eq. (18) for the imaginary times τ ∈ [0, β/2]. The normalization condition (19) results in the constraint where the factor S corresponds to the occupation of the negative frequencies with the meaning of the energy transfer ω from a system to a scattered probe particle, i.e. from the excitations existing in the system. Typically, this contribution is important in the q-region of the acoustic phonons and rotons, and one can start the optimization first by neglecting S and find an optimized solutionS 1 (q, ω). Then the normalization of a second solution is corrected by i.e. the correction S is evaluated based on the spectral densities obtained in the previous iterations (M ≥ 1). With n increasing the on fly estimation of both S est (q, ω, M ) and S is improved. This results in a fast convergence of S within few (n ∼ 5) iterations. The reconstruction of the spectral density A(q, ω) of the one-particle Green's function (17) gets a bit more involved. First, there are no simple relations between the densities at They should be worked out independently in the frequency range, −ω CO < ω < ω CO . We found that using ω CO ∼ 600 is sufficient to fit a fast drop of G 1 (q, τ ) near τ = 0, see Fig. 3. Second, G 1 (q, τ ) is not symmetric relative to τ = β/2 and should be evaluated on the whole interval, 0 < τ ≤ β. Third, the Laplace transform in Eq. (17) contains the additional Bose factor which leads to the rigorous result, 90 A simple solution which allows to use the same SO procedure as for S(q, ω) is to introduce two spectral densities which are both positive and defined for ω ≥ 0 This results in the following decomposition Using Eqs. (17), (19) we end up with two normalization constrains At low temperature and/or high excitation energies, βω 1, the integral terms in Eqs. (38), (39) can be neglected. Then both spectral densities become independent with the normalization given by the momentum distribution The latter can be directly evaluated via Fourier transform of the one-particle density matrix In contrast, at high temperatures, the normalizations (38)- (39) are mutually dependent and can be treated iteratively, similarly to Eq.
The procedure converges within few iterations.
For the detailed description of the stochastic optimization method and types of the update we refer to Ref. 89 We found of particular importance to implement the annealing allowing to escape from local minimum and minimization of the deviation (20) using parabolic interpolation. During the reconstruction, the updates involving changes in two rectangles (26) (i.e. change a weight of two rectangles, add a new rectangle, remove a rectangle) correctly redistribute the total spectral density A(q, ω) between A > (q, ω) and A < (q, ω), even if initially they are chosen to be equal. The decay of spectral weight A < (q, ω) at large q-vectors should follow that of the momentum distribution and, therefore, practically vanishes beyond the roton region (qa ∼ 6).
To speed up the optimization process and convergence of the finite-temperature corrections in the normalizations (34), (44), (45) the centers of the basis set rectangles {c i }, representingS(q, ω), A > (q, ω), A < (q, ω), have been initially normally distributed around the frequency ω χ (q), i.e. the upper bound of the lower excitation branch, see Appendix B. The optimization process was started with N, N > , N < ≈ 30 basis functions with the initial width, w = w min , i.e. equal to the frequency resolution w min /E 0 = 0.5. The optimization was performed in several iterations each of 120000 steps. In each step one of the following update types has been randomly chosen: (1) shift of a rectangle, (2) change a rectanglewidth, (3) change the heights of two rectangles, (4) add a new rectangle, (5) remove a rectangle, (6) split a rectangle into two, (7) glue two rectangles. For each update involving a change in one or couple of parameters {h i }, {w i }, {c i } we try to converge to a local minimum by parabolic interpolation. In the first 60000 steps we perform several annealing sequences with a duration 1000 − 5000 steps and temporary accept the updates increasing the deviation D n (20). Once an iteration is finished, the actual deviation D n and the optimized solution are saved and checked for the acceptance criterion, D n ≤ D min . If not accepted, we proceed to the next iteration simultaneously increasing the minimal deviation by some factor, D min → γD min , e.g. γ = (1.1 + ν), with the random  For the density correlation function we get a typical error, δG 2 ∼ 10 −4 , which is one order of magnitude smaller compared to the single-particle Green's function, δG 1 ∼ 10 −3 . The solid line is the result of the optimization, i.e. the correlation function (24) evaluated from one of the solutionsS n (q, ω) orÃ n (q, ω) which enters in the final estimation (25). The right panel shows the relative error for a set of solutionsG n (n = 1, 100). The dif-ficulty to distinguish individual curves demonstrates the convergence of the optimization process. All solutions are obtained from independent (randomly chosen) initial spectral densities. The relative error depends on the q-vector and its value is mainly determined not by the optimization result but by the statistical fluctuations near τ = β/2, when the value of the correlation function significantly drops. Still the optimization process converges. A well-behaved decay at short times (for τ or β − τ ) fixes some of excitation energies ω in the spectral density and, hence, partially the decay at large times when approaching the mid-point τ = β/2. As follows from Eq. (20) the fitted points contribute to the least deviation D n not equally but with the weight factor determined by the statistical error. Hence, the middle points with large statistical fluctuations are less important for the fit.
The influence of the statistical noise on the accuracy of the reconstruction procedure is discussed in detail in Appendix A.

A. Momentum distribution
Before analyzing of the spectral densities, we now discuss the momentum distribution which enters in the normalization (44)- (45). These data will characterize the normal and superfluid phase of 2D dipolar gas.
In Fig. 5a,b,c the low (T < T c ) and high (T > T c ) temperature momentum distribution is shown for D = 0.1, 1.75, 12.5. These are referenced in the following as weak, intermediate and strong coupling, correspondingly. The T -dependence of the one-particle density matrix (42) is shown in Fig. 5d,e,f. For low temperatures (T = 0.714, 1.0) in the superfluid phase it demonstrates only minimal changes. Our system has a finite volume and satisfies the periodic boundary conditions (PBC), therefore, the momentum distribution is evaluated at a discrete set of wavevectors, q = 2πn/L (n ∈ Z), shown by symbols in Fig. 5. For comparison, the solid line is the result (shown for T = 1.0, 3.3) obtained by extension of Eq. (42) to the limit V → ∞ with the assumption that G 1 (r, β) decays as a power-law (exponent) below (above) T c beyond the simulation box size L = √ V . The fitting parameters for T < T c are given in Tab. IV. This method agrees with the finite-size results (shown by symbols). As the system size L gets larger, more discrete values of q will get occupied dwelling on this interpolation curve.
We get a full agreement in the normal phase, when G 1 (r, β) decays fast to a small value at r = L/2 (e.g. to ∼ 10 −6 for T = 3.3 for D = 12.5), and finite-size effects are of minor importance for systems with V ≥ 165.
Some discrepancies appear in the superfluid phase. In particular, for T = 1, D = 0.1 (Fig. 5a) and qa > 5 we observe some statistical noise in n q related with the Fourier transform of G 1 (r, β). In the long wave-length limit (qa < 0.5) the discrepancies are induced by the interpolation. The density matrix is influenced by the PBC and deviates from the expected power-law decay r −ν(T ) near r = L/2. Therefore, the fit with the critical exponent ν(T ) was applied in the range 4 ≤ r ≤ 6   for V 1 = 165 and 4 ≤ r ≤ 10 for V 2 = 576. This partially allows to exclude the effects of the short-range correlations and of the finite-size errors on the expected asymptotic decay of G 1 . As predicted by the interpolation curve, in a system with off-diagonal quasi-long range order, the momentum distribution diverges when approaching zero momentum. For a finite system and discrete q we can only see the onset of this regime. For V = 165, the smallest wavevector (qa = 0.486)  is too large to capture the divergence. The low-momentum behavior will be discussed in more detail later. The information about condensate or occupation of the zero-momentum state is given by the value of the Matsubara Green's function (41), G 1 (q = 0, β). The momentum distribution in Eq. (42) is normalized by the average particle number N in volume V . Therefore, the number of particles at zero momentum N q=0 also depends on V and N , see Tab. V. For a macroscopic system this value will diverge as V → ∞, in agreement with the interpolation curves in Fig. 5a,b,c. On the other hand, at a finite temperature the condensate fraction should vanish in the thermodynamic limit, 40 lim V →∞ N q=0 /V = 0, as the zero momentum state will be depleted by thermal fluctuations. A number of occupied phonon modes with q = 2πn/L ≤ q c gets larger with increasing L.
How the thermal depletion of the low-momentum states proceeds can be analyzed for T = 1.0, 2.0 and 3.3. For D = 0.1 broadening of n q is observed in a wide range of momenta. The formation of the condensate feature at small q is accompanied by suppression of the high-momentum states. For D = 1.75 this is noticeable for qa < 1, and for D = 12.5 for qa < 0.5. Here, formation of the condensate is due to suppression of the states with qa > 6. We conclude that the increase of the coupling/density narrows the interval of momenta where a fast divergence can be observed in the superfluid phase. This, certainly, will complicate an experimental detection of a superfluid transition in a strongly correlated Bose system based only on the specific features of the momentum distribution. In particular, for D = 12.5, in a broad range of momenta (1 < qa < 6) the distribution is practically temperature independent (0 ≤ T 3.3). Here the main depletion mechanism (also at T = 0) is a strong interparticle interaction. Only a small fraction of particles occupies the q = 0 state. The depletion out of the condensate is enhanced by the presence of low-energy excitations, e.g. rotons. More insight should be given by the spectral densities of single-particle excitations and their dependence on the interaction strength.
The average particle number in the zero-momentum state N q=0 and the condensate fraction The T -dependence of the condensate is more pronounced for D = 0.1 and 12.5. In the first case, it is due to occupation of the higher-momentum states (see Fig. 5a). In the second case, T 1 and T 2 are close to the critical temperature of an infinite system (T c = 1.01, Tab. II) and thermal fluctuations play an important role.
The condensate fraction is also estimated for a larger system (V 2 = 576) to demonstrate the finite size effect. For D = 1.75, 7.5 and 12.5 we findG L/2 1 (T 1 ) = 0.220(5), 0.0438 (8) and 0.0136(5), correspondingly (G 1 ≡ G 1 /n). Compared to n 0 (T ) in Tab. V, the relation, n 0 (T ) G L/2 1 (T ), always holds and connects the reduction of the condensate with the boundary value of G 1 . The difference of both decreases with ν(T ). For D = 0.1, 0.5 and 1.75 (Tab. IV) the difference is within few percents. For D = 7.5 and 12.5 for a better agreement T should be lowered to decrease the critical exponent ν(T ).
TakingG L/2 1 as the estimation of the condensate fraction at T = 0 (when ν is small), some predictions can be made for experimental systems with a number of bosons beyond direct numerical treatment (see Appendix C).
Next, we consider the divergence of the momentum distribution as q → 0. Small system size limits the resolution in this important region to q ≥ 2π/L. To overcome this limitation we employ the "1/p 2 sum rule" [93][94][95] which holds independently on the coupling strength. In the collisionless (T T c , ρ s = n) and low-frequency (hydrodynamic) regime we can apply the ansatz by Gavoret and Nozières 46 where c is the compressional sound speed. Substitution of this ansatz in (46) yields both the spectral weight Z and the T -dependence of the momentum distribution 46 The sound speed c(T ) can be evaluated from the particle number fluctuations where κ T is the isothermal compressibility. The sound speed c and the superfluid fraction are given in Tab. VI.   The comparison with the PIMC results for T 1 = 1 and T 2 = 0.714 is shown in Fig 6. Two system sizes (V 1 , V 2 ) are considered to demonstrate the finite size effects. They are found to be negligible and are within the statistical errors of n q . The agreement with the asymptotic (48) can be confirmed for the weak (D = 0.1, 0.5) and intermediate coupling (D = 1.75). For D ≥ 7.5 we observe the onset of the slope predicted by (48). Larger system sizes (V > 10 3 ) are required to access smaller q-values. Independently, one can check the reconstructed spectral density A(q, ω) in Sec. VIII A-VIII C. In the limit q, ω → 0, for all coupling strengths A(q, ω) has a pole at ω = cq, where c coincides with the compressional sound speed. This explains a good agreement in Fig. 6a,b. However, for D ≥ 1.75, a second excitation branch appears in A(q, ω) (see Fig. 8). Its spectral weight increases with q and coupling D. This can be the reason for the observed systematic deviation in Fig. 6c,d and its onset at smaller q as D is increased.

B. Collective excitations
As discussed in Sec. IV, in the superfluid phase, the singleparticle (SP) spectrum is coupled to the spectrum of density fluctuations. The first term S A (q, ω) in Eq. (10) describes the scattering of quasiparticles in and out of a condensate, with a sharp δ-like peak quasiparticle dispersion expected at T T c . These sharp features should be present in S(q, ω) in the superfluid phase and vanish for T > T c . Hence, it is instructive to analyze the T -dependence of S(q, ω) to see any difference in both phases. We first discuss some general features observed for different couplings (0.1 ≤ D ≤ 12.5) or densities (D ∼ √ n), and then go into detail in Sec. VIII A-VIII C.
The dynamic structure factor reconstructed by the SO method is represented in Fig. 7. In the superfluid (T = 1) we observe at least two excitation branches with well defined dispersions (with the sharp energy resonances in the phonon and roton part of the spectrum). With D increasing the onset of the second high energy (H) branch systematically shifts to smaller q-vectors: for D = 0.1 (D = 12.5) it gets a significant spectral weight at qa 6 (qa 3). The H-branch does not vanish in the normal phase (T = 2.0, 3.3), but is sufficiently damped or merge with the lower branch. This implies that it is closely related to the multiparticle excitations.
The dynamic structure factor of the upper branch S H (q, ω) (see Fig. 7) is strongly influenced by the statistical noise in the imaginary correlation function G 2 (q, τ ) (18) as discussed in Appendix A. The contribution of the high energy features is exponentially damped by the factor e −τ ω . To accurately resolve the form of S(q, ω) at large frequencies requires higher accuracy. Moreover, the reconstruction procedure (Sec. VI) has a tendency to underestimate the half-width of the highfrequency structure (Sec. A 2). See also a note in Appendix D.
In contrast, the low-energy (L) features can be reproduced more accurately (Sec. A 2). The L-branch remains in the spectrum both at low and high temperatures and shows the temperature induced broadening of the half-width of the spectral peaks characterizing the inverse excitation-lifetime. The dispersion remains well defined in the acoustic range of the spectrum and gets significant broadening at large wavevectors (qa 8). Near the origin (q = 0) the half-width is increased again due to the off-resonant excitations by thermal fluctuations, ω(q) k B T . For simulated temperatures T = 1.0, 2.0 and 3.3 in Fig. 7 this effect is well observed for D ≤ 1.75.
Now we discuss the features which build up due to the correlation effects. For D ≥ 1.75 and qa ≥ 3 (see Fig. 7g,k,n) the L-branch dispersion bends down forming a local maximum -a maxon. At larger wavevectors qa ∈ (6, 8) the rotonminimum is observed. The critical coupling for the roton formation, D = 1.0 − 1.75, is in agreement with the previous analyses. 27,30,36 As the present reconstructed spectrum is free of any approximations the roton-depth is found to be lower then reported before. In Ref. 30 it was found that the upper bound estimate ω χ (q) (Eq. B20) is at least as good as the correlated basis function result (CBF) from Ref. 36 in the rotonregion and predicts a deeper roton-minimum for strong coupling D ≥ 5. The reconstructed dispersion (see Figs. 10a and 11a) shows that the correct maxon-roton dispersion is even lower than ω χ (q). According to the sum rules (B8)-( refcomp) presence of an upper H-branch and the increase of its spectral weight in S(q, ω) should push the L-branch to lower energies and, correspondingly, deepen the roton feature. We also do not exclude that the discrepancy with the CBF result (T = 0) is a temperature effect. As was shown in Ref. 61 for superfluid 4 He one observes a softening of the roton mode once approaching T c from below, while the peak in S(q, ω) shifts to lower frequencies.
In 2D dipolar systems the roton-feature is a pure correlation effect which cannot be reproduced by the Bogolubov dispersion, ω 2 (q) = 2 q + 2n 0 (T )V (q) q . The Fourier amplitude V (q) of the dipole potential 96 is positive at all q-vectors and cannot lead to the rotonization of the spectrum at weak coupling (low densities), in contrast to 3D geometry. 32,42,74 Even at the lowest coupling considered (D = 0.1) the dispersion deviates from the Bogolubov result. The basic assumption -a small depletion of the condensate, is not satisfied for D ≥ 0.1 (see Tab. V). The zero temperature analysis 36 ends up with the same conclusion.
Next, we analyze the splitting of the dispersion curve into the L and H-branch. The onset of splitting can be predicted based on the f-sum rules (B7)-(B8). In the superfluid phase both branches are well defined (Fig. 7a,d,g,k,n) and we can consider the ansatz where ω L(H) and S L(H) define the dispersion and the spectral weight for two branches. Substituted in Eqs. (B7)-(B8) the system of coupled equations can be solved with the result The solutions depend on the dispersion of the L-branch ω L (q) assumed to be known from the SO reconstruction. Two upper bounds ω F (q) = q 2 /2mS(q, T ) and ω χ (q) = 2nS(q, T )/|χ(q, T )| are defined in Appendix B and can be evaluated via the static structure factor S(q) and the static density response function χ(q). The energies of the lower branch ω L , due to a slower decay in the imaginary time, can be resolved by the SO-reconstruction more accurately than ω H and, therefore, are considered as an input. In addition, ω L are less damped at large q-vectors compared to the frequencies of the collective modes. The solution for the upper branch ω H is compared in Figs. 9a,10a,11a (solid gray line) with the full SO-reconstruction (solid symbols show position of the maxima including the half-width). The SO data are in good agreement except for the q-vectors when the frequencies ω L and ω χ overlap or the H-branch is significantly damped. Therefore, such analysis is not applied at T = 2.0 and 3.3. In summary, the SO-spectrum provides more information (and more complicated structure) than suggested by the simple ansatz (51). New parameters (additional excitation branches and their spectral weights) can be added in Eq. (51), however, to resolve them one needs to know the additional frequency moments ω n .
The solutions (52)- (53) can also predict a q-dependence of the spectral weights. For D = 0.1, S H (q) is enhanced around qa ≈ 6.36 with ω H being slightly above the recoil energy q . Next, S H (q) is slightly increasing while the intensity of the L-branch is decreasing. The SO data confirm this variation, see Fig. 7a. Such a behavior is imposed by almost a constant value of S(q) in Eq. (54) for large momenta.
The theoretical interpretation of the splitting into two branches will be further discussed in Sec. VIII A. Here we mention two possible scenarios, considering as an example the dispersion for D = 0.1 (Fig. 7a): I) After hybridization at qa ≈ 6.4 the observed lower branch ω L (q) is the continuation of the a single-particle (SP) dispersion or a dispersion of collective density modes. The high energy resonances (H-branch) correspond to combinations of two or more quasiparticle excitations, II) both dispersions of the SP and collective modes for qa > 6.4 continue above (not lower) the freeparticle dispersion q . The L-branch appears due to decay processes, then the quasi-particle energy becomes larger than the energy of a quasiparticle pair, i.e. ω L (q) > ω L (q 1 ) + ω L (q 2 ) with q = q 1 + q 2 .

C. Single-particle excitations
The reconstructed spectral density |A(q, ω)| is demonstrated in Fig. 8. The density at positive A > (q, ω) ≥ 0 and negative A < (q, ω) ≤ 0 frequencies satisfies the normalization (38)- (39) and characterizes the excitations that are either internally excited by thermal and quantum fluctuations (ω < 0) or excited by the energy transfer (ω > 0). The decay of the negative amplitude |A < | with the q-vector is due to the normalization (39) and follows the decay of the momentum distribution n q = G 1 (q, β) depending both on the dipole strength D and temperature T (see Fig. 5). The weak T -dependence of n q observed for D ≥ 1.75 and qa > 1 can be directly linked with A > and A < . Rewritten in the form 90 the T -dependence enters through the first antisymmetric component. If the spectral densities A > and A < are antisymmetric, the second term will dominate being proportional to n q . In this case the T -dependence does not enter explicitly and appears only indirectly via multi-excitation and damping effects contained in A < (q, ω). Therefore, the observed symmetry of the lowest excitation branches ω >,L , ω <,L for D = 1.75, 12.5 in Fig. 8 should result in a weak T -dependence of the momentum distribution, which is in agreement with n q (T ) in Fig. 5c,e. The SP excitation spectrum has a complicated structure and many high-frequency harmonics compared to the collective density modes. For the strong coupling D = 12.5, to fit the imaginary time decay of the Matsubara Green function (17) the reconstruction was performed on the enlarged frequency interval. Some high-frequency harmonics are out of the scale shown in Fig. 8. The detailed interpretation of the SP spectrum and check of the corresponding f -sum rules is out of the scope of the present paper and requires further analyses.
Below we discuss the differences observed in the reconstructed spectral densities S(q, ω) and A(q, ω) in the superfluid and normal phases in three coupling regimes. Possible experimental realizations were listed in Sec. III.

VIII. DISCUSSION OF THE SPECTRA
A. Weak coupling First, we analyze S(q, ω) at D = 0.1. In Fig. 9a-c we plot the position of the peaks and their half-width. The black-solid symbols denote the main spectral features in S(q, ω). Three temperatures are compared. In the superfluid phase (T = 1) in the range qa < 8 we observe a sharply peaked single excitation branch. The dispersion is accurately reproduced by the upper bound ω χ (q) derived from the f -sum rules, see Appendix B. The half-width of the peaks characterize the damping effects. The smallest half-width in Fig. 9 is limited by the frequency-resolution used in the SO-reconstruction ( ∆ω/E 0 = 1, Sec. VI). Except of the region, ω(q) k B T , the half-width stays close to this lower bound for qa ≤ 5. For larger q the damping starts to increase systematically (see T = 2, 3.3). In the region (qa ≈ 5) the dispersion curve crosses the free-particle branch, q = q 2 /2m, and broadens. This can be interpreted as hybridization and level repulsion which occurs whenever two branches cross. In the superfluid (T = 1) this broadening is only slightly increased with the q-vector and simultaneously we find the splitting into the Land H-branches. The L-branch goes well below the recoil energy and follows the dispersion of the zero sound (ZS). Interestingly, this feature can be only observed in the superfluid phase. Based on the theory of the hybridization of G 1 (q) and χ(q) in the presence of a condensate (their poles are reproduced in each function), this branch should be due to the coupling with the SP spectrum. However, at large q this branch is not observed in the spectral function A(q, ω). For a further test, one needs to evaluate the two-particle Green function and check the two-particle spectrum where this branch should get a significant spectral weight.
This behavior is not reproduced in the normal phase, see Fig. 9b-c. The damping is systematically increasing with the q-vector which is typical for a normal gas/liquid. In addition, we observe that the L-branch is now damped and saturates near a "plateau" well below the ZS-dispersion. This shift to lower frequencies is accompanied by the shift of the H-dispersion to higher excitation energies. This result is expected from the sum rule (B8). We remark that the half-width of the H-branch can be underestimated by the reconstruction (see Appendix A).
What is common for both superfluid and normal phase is the linear acoustic phonon dispersion, ω(q, T ) = c(T )q, which extends to qa ≈ 5 at T = 1.0, qa ≈ 4 at T = 2.0 and qa ≈ 3 at T = 3.3. The results for the isothermal sound speed c(T ) are given in Tab. VI. Its non-monotonic Tdependence (included in Fig. 9) can be explained by suppression of the density fluctuations in the pre-superfluid regime (T 2) and the corresponding non-monotonic behavior of the compressibility κ T . The variation of the zero sound 97 δc = |1 − c(2)/c(1)| evaluated from Tab. VI is most pronounced (δc ∼ 7%) for weak coupling D = 0.1. At larger coupling D the density fluctuations due to formation of a local condensate and superfluid density are strongly suppressed by correlation effects (the relative change of the compressibility is also reduced).
To test the hybridization with the SP spectrum in the superfluid, Fig. 9d-f shows the spectral function A(q, ω). The statistical error in evaluation of G 1 (q, τ ) is larger then in G 2 (q, τ ) (compare Fig. 3 and Fig. 4), as a result the reconstructed SP dispersion curve at low q-vectors (Fig. 9d) shows  some statistical fluctuations around the ZS dispersion. At high temperatures the convergence to the linear dispersion (in the range qa 5) is obscured by the increased damping (Fig. 9e,f). However, at larger q the accuracy of the reconstruction should be improved as the Matsubara Green function is evaluated more accurately with less statistical noise. In comparison with the collective excitations, the overall slope of the SP dispersion is less influenced by temperature.
The linear dispersion observed as q → 0 is in agreement with the Hugenholtz-Pines sum rule for the self energies 73 with the result that the spectrum A(q, ω) (at T = 0) is gapless in the long-wavelength limit. Similar result for Bose liquids at T = 0 has been worked out by Cheung and Griffin. 98 Due to the hybridization the density-density response and oneparticle Green functions in the superfluid share the same singularities. At q → 0 these are the phonon poles as was shown by the field-theoretical calculations of Gavoret and Nozières 46 with the velocity precisely equal to the ZS speed. This result is confirmed for superfluid helium both experimentally 81,99 and theoretically. 46,49,50,73 This justifies the Landau-Feynman interpretation of the sound waves as elementary excitations in Bose systems. Excited out of a Bose condensate a quasiparticle excites the collective oscillations in a superfluid.
In the dipole system, the linear dispersion terminates at qa ≈ 3−5.4 (depending on T ) with the splitting and broadening. Then the dispersion curve goes slightly above the recoil energy q . Fig. 9a shows that the same branch (H-branch) is observed in the dynamic structure factor. This comparison confirms, that, indeed, in the superfluid phase both spectra are coupled. In the superfluid (Fig. 9a,d) at D = 0.1 both dispersions are linear and coincide with the compressional sound. The broad high-energy branch at T ≥ 2 and large wavevectors (qa 5.4) is identified with the multiparticle excitations. The halfwidth of this dispersion is reduced in the superfluid (Fig. 9a) but not drastically. Still the damping is much larger compared with the lifetime of the SP-mode in Fig. 9d.
Finally, we interpret the L-branch in Fig. 9a as a continuation of the ZS mode which terminates (loses its intensity) for qa 9.

B. Intermediate coupling
In Sec. VII B we discussed formation of the maxon-roton branch for D = 1 . . . 1.75 and failure of the Bogolubov spectrum. At T = 1 and qa ≤ 3 the dispersion relation is accurately reproduced by ω χ , see Fig. 10a. Similar to D = 0.1 the deviations appear together with the Hbranch. This branch has a weak q-dependence and spans the frequency range ω H ≈ 20 . . . 30. The resonant frequencies can be explained by a combination of low laying excitations, ω(q 1 )+ω(q 2 ) ≈ 10+10 ≈ 20. The q-dependence is in agreement with the solution ω H (q) from Eq. (52). At qa > 8 there is an abrupt broadening of the linewidths. For qa > 9 the Land H-branches merge into a single dispersion which goes slightly below q , see Fig. 10a. The peak position is accurately predicted by ω χ . A similar shift below q was experimentally observed for liquid 4 He 60 .
In the presence of the H-branch (Fig. 10a) the low-energy dispersion has a well defined maxon-roton feature at q ∈ (6, 8). Its depth is lower than found in the previous analyses. 30,36 The present approach is not limited by three-particle decay processes of CBF theory 36,42 and from first principles includes the spectral weights and energies of all excitation branches, here the Land H-branch. They are coupled via the f -sum rules (B7)-(B9) as in the simple example (52)- (54). Therefore, the under(over)estimation of one of the branches has a strong influence on the full spectrum.
In agreement with Gavoret and Nozières,46 in the q, ω → 0 limit the density response and SP spectra converge to the ZS dispersion. The sound speed c is given in Tab. VI. At temperature T = 1 in the SP spectrum, we observe, in addition, a second excitation branch with a finite gap ω g ≈ 40 at q = 0, see Figs. 8g, 10d. The gap value can be explained by the excitation of two quasiparticle with the opposite momentum, ω(q) + ω(−q) ≈ 40 with ω(q) ≈ 20 for qa ∈ (2, 6).
We can confirm that the lower dispersions in S(q, ω) and A(q, ω) coincide very well up to a maxon region (Fig. 10a,d).
In the maxon-roton region, both spectra are not much similar. Starting from qa ∼ 4 the H-branch in Fig. 10a bends up, while the SP dispersion in Fig. 10d goes down approaching a local minimum near the roton wavevector qa ≈ 7. Both spectra coincide again only for large momenta and follow the ω χ -dispersion. From these results we can not confirm existence of the unified excitation branch in the whole range of q-vectors, as was discussed in Sec. I and argued in Ref. 68 Here, we note that the SP spectra in Fig. 10d-f are more difficult to reconstruct as accurately as S(q, ω). The statistical noise of the Matsubara Green function G 1 (q, τ ) is by factor 10 larger than in G 1 (q, τ ) (see Fig. 3), while positions of the energy resonances are influenced by the noise level as is shown in Appendix A. For A(q, ω) one also needs to accurately reconstruct the additional resonances at ω < 0 (Fig. 8). This can be a reason why for some q-vectors the spectrum in Fig. 10d-f does not behave like a continuous dispersion relation.
Next, we continue to discuss the collective excitations in the normal phase, Fig. 10b,c. The splitting of the dispersion curve, previously observed for D = 0.1, here is also reproduced. At T = 2 the roton peak shifts to low frequencies (Fig. 10b). Similar behavior was observed in the normal phase of 4 He, 61,62 where above T c the softening of the roton mode was found with the roton energy approaching the zero frequency. At high temperature (T = 3) the roton-feature gets broader and eventually becomes nearly dispersionless and saturates as a "plateau" (Fig. 10c). Simultaneously the upper branch demonstrates a large damping and gradually transforms into the free-particle dispersion. In conclusion, our analyses show that the roton mode, which is very sharp in the superfluid phase (T = 1), remains also in the normal phase of the dipolar gas. Its lifetime drops significantly when the temperature is increased from T = 1 to T = 2. Then the damping stays practically constant by further increase to T = 3.
The final form of the dispersion (lower branch) is typical for a normal gas/fluid. After the phonon part the dispersion slightly bends down. The self-consistent field produced by the density fluctuations cannot support stable shortwavelength collective modes. In the normal phase the dispersion relation can be compared with the results of MD simulations for 2D classical dipoles 72 and predictions of the quasilocalized charge approximation, 100 after the correct mapping of the quantum dipolar coupling D = p 2 m/ b 2 a on the classical coupling parameter Γ D = p 2 /a 3 k B T .
Interestingly, that in the SP spectrum in Fig. 10d also resembles a roton mode but with larger excitation energies, compared to the roton in the density fluctuation spectrum. The slope of the L-branch is close to ω χ (see Fig. 8). Besides the L-branch, there is a well distinguished upper branch. This is a common spectral feature for all coupling strengths (Fig. 8).
The H-branch is very pronounced for D ≥ 1.75 and accumulates most of the spectral density for ω > 0.  Fig. 9 for further details.

C. Strong coupling and rotonization
As a third example we consider strong coupling. First, the dispersion in the superfluid phase will be discussed. For simulated temperatures, at D = 7.5 and 12.5, there is a local crystalline ordering with fluctuating orientation. This regime is followed by the freezing transition at D = 17 (1). 27,28 At D = 12.5 the static structure is peaked at qa = 6.65 and the phonon-maxon-roton feature is well pronounced, see Fig. 7k,n. Due the 0−sum rule (B7) the rotons are the dominant excitations (the integrated spectral density is proportional to S(q)) and the depth of the roton minimum affects the critical temperature of the superfluid transition. 30 The full reconstruction clearly demonstrates the temperature effect, see Fig. 7n. Near the roton-minimum at E R = ω L (6.36) ≈ 4.9 the dynamic structure factor S(q, ω) shows an asymmetric broadening to lower energies. This also applies to D = 7.5 ( Fig. 7k) with the roton gap E R = ω L (6.36) ≈ 6.95. For D = 12.5, the simulated temperature T = 1 is close to the BKT-temperature of an infinite system (T c = 1.01) and the superfluid fraction is reduced to ρ s = 0.81 (see Tab. VI). Therefore, we observe some additional damping compared to D = 7.5 with ρ s = 0.95 and T c = 1.22.
Except the maxon region, the upper bound ω χ (q) reproduces well the L-branch for qa ≤ 8 including the roton minimum, see Fig. 11a,b,c. The deviations at large q are due to decay processes. Once the quasiparticle energy exceeds the energy of a quasiparticle pair, it becomes unstable and would decay into this pair. For a strongly correlated system the lowest quasiparticle energy corresponds to a roton, E R , and, hence, the dispersion curve should not exceed 2E R .
Surprisingly, at D = 12.5 the roton minimum is so deep, that the energy of the two-roton state 2E R ≈ 10 is much lower than the maxon energy. As a result we observe this state symmetrically on both sides of the roton minimum at qa ≈ 6.36 (see horizontal dashed lines in Fig. 11a). A similar effect has been discussed for 4 He at high pressures (above 18 bars). 101 For D = 7.5 (see Fig. 7k) the double roton-feature at 2E R ≈ 14 is also reproduced. The two-roton state, first predicted by Pitaevskii, 90 is well known experimentally for superfluid 4 He and was partially discussed in Sec. I. Its experimental verification was complicated, since the sharp energy resonance merges with a broad multi-excitation background and the measurements required a high instrumental resolution. 102 In contrast, for 2D dipolar gases with D ≥ 7.5 there is no such complication, since both components are well separated (see Figs. 11a and 7k,n).
Next, we discuss the maxon region, qa ∈ (2, 4) (Fig. 11a). Here, we observe a splitting of the dispersion. The lowest branch is the two-roton state. Its energy is well below ω χ (q) (solid black curve) which corresponds to a prediction of a continuous dispersion which neglects decay processes. It provides an accurate prediction only up to the maxon. Here, the upper branch starts which continues with few oscillations (see Fig. 7n). In the roton region, q ∈ (5, 7), its energy is in the range ω ≈ 20−25 and there is a large gap to the low frequency roton state. At large momenta the H-branch merges with the strongly damped dispersion of multiparticle excitations near the recoil energy. Below we give a possible interpretation of the observed behavior.
In the interpretation proposed by Glyde and Griffin 47 the hybridization of the dispersion near a maxon is due to the crossing of two different branches, i.e the acoustic phonons which dominate at low q and the maxon-roton mode. This cross-over behavior should be characterized by a double peak structure. Indeed, in the simulations we observe a double peak structure, but its origin is different. It comes from the decay of quasiparticles into pairs of lower energy (here the two-roton state). If this decay is not allowed by kinematics then we always observe a continuous dispersion relation (see Fig. 9a-c and 10a-c), in agreement with Ref. 68 The dispersion does not lose its spectral intensity or increases its halfwidth as it would be for the case for hybridization of two distinct branches. Further analyses for the case when 2E R > E maxon will be helpfull to finally clarify wether the maxon for strongly coupled systems is a transition region between phonons and rotons or a part of a continuous dispersion relation.
Another important issue concerns the origin of the roton minimum. For D = 1.75 we have observed a weak rotonization of the spectrum both in the superfluid and normal phases. Its existence is independent on the presence of a condensate. For strong coupling, D = 12.5, the maxon-roton feature is more pronounced including the normal phase (Fig. 11b,c). Temperature, practically, does not influence the slope of the low-frequency roton-two roton state, but has a noticeable effect on the damping characterized by the half-width. We can conclude, that the roton minimum is not related with the coupling to the SP excitation spectrum, being an intrinsic col- lective density response mode. Its formation is solely due to the short-range correlations which supports earlier suggestions 70,71 and is in agreement with the analysis of 2D dipolar system in the classical limit. 72 The next question to be addressed is whether the roton feature can be also found in the SP spectrum A(q, ω). The peak positions of A(q, ω) and their half-widths are shown in Fig. 11d-f. We do not discuss the uppermost branch. In contrast to S(q, ω), the high frequency branches carry the main part of the spectral density both in the superfluid and normal phases. They also dominate over the L-branch in the rotonregion, see Figs. 8n-p. Two lower branches are well separated at T = 1, see Fig. 11d. Temperature systematically reduces the energy of the upper branch (H) and flattens the lower branch (L).
In the superfluid phase (Fig. 11d) the L-branch has a roton-like dispersion, resembling the density response roton in S(q, ω). Similar to the density-roton, the left roton-branch does not reach the maxon maximum (shown by ω χ ), but for qa ∈ (2, 4) decays into quasiparticles with lower energy.
The T -dependence of the roton in A(q, ω) is different from the density-roton in S(q, ω). Both in the superfluid and normal phases it stays very near or above ω χ . There is no contradiction as ω χ is the upper bound for the collective density modes. The dispersion becomes almost flat at T = 3.3 (Fig. 11f) and the roton mode completely vanishes. In contrast, the density-roton mode and two-roton states do not change much except for the broadening. In the long wave-length the dispersion converges to the ZS. The lowest branch is gapless (ω, q → 0) in agreement with the Hugenholtz-Pines sum rule.
On the right side (qa 7) the roton-branch follows the dispersion predicted by ω χ for the density modes. Simultaneously, we observe that there are no excitations at the freeparticle dispersion q . This is in contrast, with the behavior in the normal phase (Fig. 11e,f). Now, the density of the SP excitations is centered around the free-particle branch q . The dispersion is broadened significantly compared with the one in the superfluid phase. From this result we can argue that the hybridization scenario, discussed in Sec. I, is realized in our 2D system of bosonic dipoles. The poles of the dynamic response function χ(q, ω) predicted based on the sum rules result ω χ (q) (see Sec. B) can be directly observed in the elementary excitation spectrum A(q, ω) but only in the superfluid phase, when this is allowed by the hybridization. On the other side, we do not see exactly the same poles in S(q, ω), as the collective spectrum is strongly renormalized (perturbed) by the quasiparticle decay processes. They shift the positions of the energy resonances and lead to a slightly different dispersion relation. This can be a possible explanation for our numerical results.
Finally, we note that many additional modes are observed in A(q, ω) which, probably, are intrinsic to the single-particle excitations. Their interpretation requires further analysis which goes beyond the present work.

IX. SUMMARY
Using the path integral Monte Carlo technique in the grand canonical ensemble we evaluated the imaginary time density response and the one-particle Matsubara Green's function. The stochastic optimization method was used to accurately reconstruct the dynamic structure factor S(q, ω) and the SP spectral density A(q, ω) in a wide range of wavevectors.
As a result, a first-principle treatment of the dynamic structure factor was presented both for weakly and strongly interacting Bose-condensed systems. The temperature-density dependence of the phonon-maxon-roton dispersion, its damping and hybridization was resolved by the simulations. The dispersion relations were compared vs. the upper bound derived from the frequency sum rules. The introduced approach allows to discuss the physical origin of the roton minimum and is a useful benchmark for theoretical studies and interpretations of experiments on quantum liquids/gases in different physical realizations. As an example, a 2D dipolar Bose gas at weak, intermediate and strong coupling has been analyzed in detail.
Our goal was to demonstrate that the single-particle and density-response spectra are coupled and share common poles due to hybridization in a superfluid, as was predicted theoretically. 46,49,50,73 Once dipolar bosons become superfluid (T ≤ T c ), we observe a similar dispersion relation in both spectral densities A(q, ω) and S(q, ω). However, we cannot confirm that this dispersion can be observed at all wavevectors as a unified excitation branch 68 when the coupling is strong.
In this case, the collective spectrum S(q, ω) is strongly renormalized by the quasiparticle decay processes. This, according to the frequency sum rules, should necessarily lead to a shift of the S(q, ω)-maxima from the position of the original energy resonances. One prominent example is the decay into the two-roton state observed in the strongly correlated Bose gas.
Further, we found that the two-roton state exists independently of the presence of a condensate and, therefore, conclude that the roton mode has a classical origin due to shortrange particle correlations. Its observation in classical liquids is, probably, prevented by the softening and overdamping. 61,62 For quantum dipolar liquids at high coupling/density (here D ≥ 7.5) the reconstructed S(q, ω) predicts a pronounced roton minimum. The sharp resonant roton feature is transfered from a superfluid to a normal phase but with a significant damping. Our reconstruction of the SP spectrum A(q, ω) in a superfluid also predicts the roton-like feature which, however, vanishes in the normal phase. Here it is substituted by the excitation mode peaked around the recoil energy.
The common features of the density fluctuation spectrum S(q, ω) are the following. For large momenta (qa > 8) we observe a strongly damped mode near the recoil energy both in the superfluid and normal phases. As expected, the phonon dispersion at small wavevectors exists independently of the presence of spatial coherence in the one-particle density matrix.
These interesting aspects of the excitation spectrum near T c may be addressed in the near future in the ongoing experiments on dipolar systems. 3,4 The analyzed dipolar coupling strengths can be realized in 2D systems of magnetic atoms, polar molecules and indirect excitons as discussed in Sec. III.
Appendix A: Model spectral densities: accuracy of the reconstruction To demonstrate accuracy of the reconstruction procedure, below we consider a test for several spectral densities S(ω). The goal is to analyze how a level of statistical noise, present in the imaginary correlation function, influences positions of energy resonances and their half-width. To characterize deviations, in addition, we evaluate the frequency moments ω n to check fulfillment of the sum rules (B7)-(B9). In our model, they are substituted by the exact integral properties of S(ω).
Here, we analyze a dynamic structure model and assume that the detailed balance holds, S(−ω) = e −βω S(ω).
Further, for convenience, S(ω) is substituted by a similar distribution S 0 (ω) discretized into rectangles of width ∆ω By discretization evaluation of the correlation function can be easily performed as As a second benefit, the integral properties (frequency moments) in the basis of rectangles {h i = S(ω i ), w i = ∆ω, c i = ω i } can be written down explicitly and evaluated via with Ei(x) = x −∞ dt e t /t. To avoid the divergence in Eq. (A5) when (c 1 − w 1 /2) = 0, we set S 0 to be non-zero for ω ∈ [δ, ω CO + δ] (with δ = 10 −3 and the cut-off frequency ω CO = 100) and choose, correspondingly, ω i = (i− 1 2 )∆ω+δ. Eqs. (A3)-(A5) can be used to evaluate the frequency moments (f -moments) of the ensemble average (25) by applying to each term in the ensemble. This result was used to check the sum rules of the reconstructed S(q, ω) for 2D dipolar bosons. However, as we discuss below, an accurate fulfillment some of the sum rules [relative error 10 −5 . . . 10 −3 ] does not necessarily guarantee a correct reconstruction of the shape of S(ω).
In our tests S 0 (ω) is specified by a linear combination of M Gaussians which vary in the peak position Ω j , the variance σ j and the weight factors p j . Similar functions has been used in Ref. 38 to test a GIFT reconstruction procedure (the genetic inversion via falsification of theories). We refer to Ref. 38 for the complementary tests and discussion of limitations by reconstruction from imaginary times.
To simulate the statistical noise inherent to QMC simulations, the imaginary time correlation function is perturbed by a random noise ν sampled from the normal Gaussian distribution and multiplied by the statistical error δG(τ i ) at imaginary time τ i . The time dependence δG(τ i ) was taken from a typical PIMC run (e.g. a dipolar gas at D = 1.75 and β = 1) and rescaled by factor γ (1/5 ≤ γ ≤ 5) to model the situation when the reconstruction is performed for a larger or smaller MC sample size. Following typical PIMC data presented in Fig. 2, as setup parameters we used: the inverse temperature β = 1, the time step δτ = 1/100 (τ i = iδτ ∈ [0, β/2] and i = 0 . . . 50) and the statistical error δG(τ i ) of the order 10 −5 . . . 10 −3 [with G(τ = 0) ∼ 1 due to normalization in Eq. (A6)].
How particular spectral features can be downgraded by statistical errors is clarified by few examples.

Spectral density with double peak
This example is common for experimental and theoretical studies of superfluid helium when two distinct excitation branches can be resolved in the low temperature (T < T c ) spectrum. The reference spectral density S 0 (ω) is specified by the following parameters: discretization width in Eq. (A1) ∆ω = 0.5, two Gaussians with Ω 1 = 10, Ω 2 = 25, σ 1 = 0.5, σ 2 = 2.0, p 1 = 1 and p 2 = 3. Parameters used in the reconstruction [see Eq. (29) in Sec. VI]: the cut-off frequency ω CO = 100, the allowed width of rectangles w i ≥ 0.5 (i = 1, N ) in the basis of size N = 80, and M = 200 − 400 solutions in the ensemble average (25). Fig. 12a-f demonstrates the results. In the error-free case (δG(τ i ) = 0, first column) the reconstructed spectrum (solid red line) quite accurately reproduces the position and halfwidth of the first peak in S 0 (ω) (step-like curve) even on the logarithmic scale, see Fig. 12d. The center of the second peak is accurately captured but the shape is perturbed. The reconstruction overestimates the spectral density at the peak center and the wings. The corresponding relative deviation from the reference correlation function G 0 (τ ) is shown in Fig. 12k and is of the order of 10 −6 . . . 10 −5 . Here we show the deviations of the correlation functions G n (τ i ) corresponding to difference spectral densities S n (ω) used in the ensemble average [only 25 samples are shown from typical n = 200 . . . 400]. This is the best result we can achieve with the stochastic optimization method (Sec. VI) using as an input the error-free correlation function G 0 (τ ).
An increase of the basis size in Eq. (A2) to N = 160 does not bring any substantial improvement. The reconstruction result is mainly determined by the level of statistical noise. We also tried to fix the f -moments (A3)-(A5) to their reference values ω n S 0 (n = −1, 0, 1). They were included in (20) via a set of Lagrange multipliers and subjected to minimization. However, the spectrum shape was downgraded by appearance of systematically shifted peaks which lead to an exact fulfillment of the sum rules at a price of a larger deviation from the reference function G 0 (τ ) along the imaginary time. Surprisingly, even if the moments are not fixed, their relative error δ ω n = |1 − ω n S(ω) / ω n S 0 (ω) | is small, see lower panels of Fig. 12.
By turning on the noise in the input data (A7) the peak po-sitions are slightly perturbed, see Fig. 12b,c,e,f. The accuracy in the peak position is acceptable even in the worst case when the error δG(τ i )/G(τ i ) is in the range of 2% for τ → 0.5, Fig. 12i,m. However, the half-width of the second peak is underestimated and should not be taken seriously in a reconstruction from QMC data with similar error level. The advantage to use the ensemble average (25) instead of individual solutions S n (ω) becomes evident from the panels k-m. While one solution G n (τ ) has a large error at a particular imaginary time, it is reduced by the averaging with other solutions.
The errors in δ ω n are due to a slightly different shape of S(ω) and the discretization parameter ∆ω in Eq. (A1). The step-like shape of the reference distribution S 0 can not be reproduced in the smoothed ensemble average result. This results in different f -moments even if the shape of both densities is close.
Second and third columns in Figs. 13,14 show the reconstruction with the statistical errors in G(τ ) (cf. the error bars δG(τ i ) in Fig. 13h,i and the solid black line δG(τ i )/G(τ i ) in panels l,m).
In Fig. 13b,c we see a trend similar to the double peak structure, cf. Fig. 12b,c. There is no appreciable difference in the position and the half-width of the first peak, but the noise shifts the second peak and results in a systematic underestimation of its half-width. We can conclude, that with similar relative errors 10 −3 . . . 10 −2 in the imaginary correlation function from QMC, a spectral shape of a high-frequency branch should not be taken too seriously as it is obscured by statistical noise. Nevertheless, the relative error to the reference f -moments [see Tab. VII] is within the range δ ω n = 10 −5 . . . 5·10 −4 (n = −1, 0, 1) and is surprisingly small. This demonstrates a difficulty to qualify a reconstructed spectrum based only on fulfillment of few sum rules. We return to this point again in the summary. Fig. 14 shows an example when a sharp resonance is placed on a broad distribution which extends to low frequencies. Again in the error-free case the reconstructed density is reliable including the log-scale.
It becomes difficult to reconstruct this sample when a tiny noise is introduced. The relative error (cf. Fig. 14i) is by two orders of magnitude smaller (3 · 10 −4 ) then in Fig. 13m, though it sufficient to obscure the sharp resonance, cf. Fig. 14c, and two distributions merge into one. However, the reconstruction does reproduce the half-width of the background density and its asymptotic decay with frequency (cf. panel f). The errors in the f -moments for the reconstruction in Fig. 14a (n = −1, 0, 1).
In second and third columns of Fig. 15 we study effect of the noise. In both cases, the relative error for short (τ < 0.1) and long (τ → 0.5) imaginary times differs by orders of magnitude and reach 7% and 30% at τ = 0.5 (cf. Fig. 15l,m). Still the properties of the reference spectral density are not completely lost. It is possible qualitatively to reproduce the number of peaks and their half-width in Fig. 15b,e and the reminiscence of the first peak, the overall extension over a broad frequency interval and the asymptotic decay in Fig. 15c,f. Similar to the previous tests, the high-frequency peaks can be shifted from their correct position. The shift does not exceed 17%.

Summary
We conclude, from the performed tests, that the reconstruction reliability crucially depends on the relative noise in the input correlation function, δG(τ i )/G(τ i ). This quantity can be used to qualify reconstructed spectra, in particular, presented in Sec. VII B, based on the downgrade of the spectral features analyzed in this section. In general, the relative noise in the density-density correlation function increases with τ and the q-vector due to a faster decay of G(q, τ ) (see Fig. 4a). As a result the reconstructed S(q, ω) in Fig. 7 at large wavevectors becomes less accurate compared to the low momentum part of the spectrum.
When only two well separated energy resonances are reconstructed, the peak positions remain stable against the noise (cf. Fig. 12). A noticeable effect is observed in the half-widths of the peaks.
The reconstruction results are more sensitive to the noise level for the overlapping bimodal (Figs. 13 and 14) and the multipeak distribution (Fig. 15). As expected, details on how one individual feature (a peak position and a half-width) is reconstructed have a strong influence on the rest of the spectral density. Even if the low frequency moments are fulfilled with high accuracy 10 −6 − 10 −4 , this is not sufficient to guarantee a reliable reconstruction. A large class of solutions G n (τ ) can fit within the error-bars the correlation function (A7). The broader is the reference distribution S 0 in the frequency space, more "degrees of freedom" are allowed to be used in the reconstruction. They can be distributed in many different ways to satisfy the sum rules. As a result the ensemble average can smooth specific features of the reference spectrum.
Next, consider a Taylor expansion of the correlation function G(τ ) = Several first frequency moments can uniquely specify only a short time asymptotic of G, and, in general, all moments are required to reconstruct G at the full time scale. Vise versa, an exact fit to G(τ ) at all τ necessarily means that all frequency moments are uniquely defined. However, if there is some uncertainty, i.e. G(τ )+δG(τ ), it can be fulfilled by different expansions in ω n , where the successive terms mutually compensate, due to the sign-alternating prefactor in (A10). Now the shape of S(ω) is not unique.

Appendix B: Sum rules upper bound for the lower excitation branch
The upper bound for the lower excitation branch with a negligible damping can be derived using the decomposition of the density-density response function (imaginary part) [χ(q, ω)] = χ 1 (q, ω) + χ 2 (q, ω). (B1) It is assumed that the δ-peak quasiparticle contribution with the energy ω q > 0 is not overlapping with the multi-particle excitations at higher energies We note that χ 1 (2) is odd in ω. The fluctuation dissipation theorem and the detailed balance are written as [χ(q, ω)] = −π 1 − e −βω S(q, ω), S(q, −ω) = e −βω S(q, ω).
In the following, the 0-sum rule, the longitudinal f −sum rule 39,40 and the compressibility sum rule (all valid at any T ) will be employed These constraints can be used to characterize the quality of experimental/theoretical spectra, i.e. check how well they satisfy the sum rules.
For a dispersion consisting of a single sharp excitation branch the upper bound ω χ becomes almost exact, e.g. D = 0.1 and Fig. 9a-c. For other couplings, the phonon-maxon and roton features are also accurately reproduced (Fig. 9a, 10a, 11a). The reason is the damping of high-energy contributions by the factor 1/ω. This significantly improves the Feynman upper bound (B16), where the spectral weight of high frequencies scales as ω. This leads to the general relation, ω q (T ) ≤ ω χ (q, T ) ≤ ω F (q, T ), valid in a superfluid phase for the lower excitation branch. In the long-wavelength limit, we expect that both upper bounds converge to the isothermal sound dispersion, i.e. ω χ (q) ≈ ω F (q) ≈ cq. Then the product of (B16) and (B20) results in the compressibility sum rule lim q→0 |χ(q, T )| = n mc 2 (T ) .
The offset from zero increases with temperature and decreases with the coupling strength/density (the sound speed increases with D, see Tab. VI).