Dynamical structure factor of one-dimensional hard rods
DDynamical structure factor of one-dimensional hard rods
M. Motta, E. Vitali, M. Rossi,
2, 3
D. E. Galli, and G. Bertaina Department of Physics, College of William and Mary, Williamsburg, Virginia 23187-8795, USA Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy International Center for Theoretical Physics (ICTP), Strada Costiera 11, I-34154 Trieste, Italy Dipartimento di Fisica, Universit`a degli Studi di Milano, via Celoria 16, I-20133 Milano, Italy
The zero-temperature dynamical structure factor S ( q, ω ) of one-dimensional hard rods is com-puted using state-of-the-art quantum Monte Carlo and analytic continuation techniques, comple-mented by a Bethe Ansatz analysis. As the density increases, S ( q, ω ) reveals a crossover from theTonks-Girardeau gas to a quasi-solid regime, along which the low-energy properties are found inagreement with the nonlinear Luttinger liquid theory. Our quantitative estimate of S ( q, ω ) extendsbeyond the low-energy limit and confirms a theoretical prediction regarding the behavior of S ( q, ω )at specific wavevectors Q n = n π/a , where a is the core radius, resulting from the interplay ofthe particle-hole boundaries of suitably rescaled ideal Fermi gases. We observe significant similar-ities between hard rods and one-dimensional He at high density, suggesting that the hard-rodsmodel may provide an accurate description of dense one-dimensional liquids of quantum particlesinteracting through a strongly repulsive, finite-range potential.
I. INTRODUCTION
One-dimensional (1D) quantum systems are subjectto intense research, due to their theoretical and exper-imental peculiarities [1–3]. On the theoretical side, thereduced dimensionality enhances quantum fluctuationsand interaction, giving rise to unique phenomena like thenonexistence of Bose-Einstein condensation [4–8] and thebreakdown of Fermi-liquid behavior [9]. On the experi-mental side, a 1D system is realized when a 3D systemis loaded into an elongated optical trap or confined intoa narrow channel, and the transverse motion is frozen tozero-point fluctuations.Remarkably, several 1D many-body models of consid-erable conceptual and experimental relevance are exactlysolvable [10–15], and provide a precious support for un-derstanding static and dynamic properties of interacting1D systems in suitable regimes [16–28].In particular, the behavior of 1D systems with a hard-core repulsive interaction, like He or other gases ad-sorbed in carbon nanotubes [28–32], can be understoodmaking the assumption that particles behave like a gasof impenetrable segments or hard rods (HRs). Indeed, athigh density, the principal effect of a short-range hard-core repulsive interaction is volume exclusion. Therefore,a reasonable approximation of the actual microscale be-havior of the system can be obtained by taking into ac-count the volume exclusion phenomenon only, neglectingall other details of the interaction: within this approach,the system is described as an assembly of HRs of a suit-able length a .The recognition that volume exclusion is the most im-portant factor in analyzing short-range hard-core repul-sive interactions in high-density classical systems datesback to the seminal work by van der Waals [33] and Jeans[34]. It was later recognized [35, 36] that the statisticalmechanics of a system of classical HRs is exactly solv-able. In 1940 Nagamiya proved [10] that also a system ofquantum HRs is exactly solvable using the Bethe Ansatz technique, and imposing a special system of boundaryconditions. Nagamiya’s treatment was later adapted bySutherland [13] to the more familiar periodic boundaryconditions.It is remarkable that local properties of the HR modelare independent of the particles being bosons or fermions[7], since in 1D the hard core interaction creates nodes inbosonic wavefunctions which can be completely mappedto the nodes of fermionic wavefunctions. Only non-localproperties differ, such as the momentum distribution [31].Even if the eigenfunctions and eigenvalues of the HRmodel can be determined exactly, so far the only system-atic way to obtain a complete description of the ground-state correlation functions of the model has been theVariational Monte Carlo (VMC) method [31, 32]. Dy-namical properties have also been addressed by using thevariational Jastrow-Feenberg theory in [37].In the present work, we resort to state-of-the-art pro-jective quantum Monte Carlo (QMC) [38–40] and an-alytic continuation [41] techniques to compute the dy-namical structure factor, S ( q, ω ), of a single-componentsystem of HRs. This analysis is supported by the BetheAnsatz solution of the elementary excitations of themodel, following [12].The dynamical structure factor characterizes the linearresponse of the system to an external field which weaklycouples to the density. In the context of quantum liquids,it can be probed via inelastic neutron scattering [42, 43],while in the ultracold gases field it can be probed withBragg scattering [27, 44], also implemented via digitalmicromirror devices [45], or cavity-enhanced spontaneousemission [46].While the low-energy properties of S ( q, ω ) are univer-sal and can be described by the Tomonaga-Luttinger liq-uid (TLL) theory [47–51] and its recent and remarkablegeneralization, called the nonlinear TLL theory[3, 52],high-energy properties depend explicitly on the shape ofthe interaction potential, and lie in a regime beyond thereach of those theoretical approaches. Due to such limita- a r X i v : . [ c ond - m a t . o t h e r] A ug tion, we rely on QMC to estimate S ( q, ω ) for all momentaand energies.The HR model, and its solution by Bethe Ansatz, isdescribed in Section II. The methods used to computethe dynamical structure factor are reviewed in SectionIII. Results are presented and discussed in Section IV,and conclusions are drawn in the last Section V. II. THE HARD-RODS MODEL
Hard rods are the 1D counterpart of 3D hard spheres[31, 32]. The interparticle hard-rod potential is V HR ( r ) = (cid:26) ∞ | r | ≤ a | r | > a , (1)where a is the rod size. The Hamiltonian of a systemof N particles inside an interval [0 , L ] of length L withinterparticle HR potential is H = − (cid:126) m N (cid:88) i =1 ∂ ∂r i + N (cid:88) i The solution of the HR Hamiltonian (2) was first ad-dressed by Nagamiya [10], relying on the Bethe Ansatzmethod [53]. The author substituted PBC (3) withslightly different boundary conditions, motivated by thestudy of particles arranged on a circle [10]. The solu-tion of the HR Hamiltonian by Bethe Ansatz was subse-quently addressed by Sutherland in [13], applying PBC.In the present Section, we provide a detailed reviewof the solution of the HR model following the methodof Refs.[11, 12], and a detailed description of its elemen-tary excitations. This is a key ingredient that permits tocharacterize the singularities of S ( q, ω ) predicted by thenonlinear Luttinger liquid theory (Sec. II D). In order to solve the HR Hamiltonian (2), followingRef.[10], let us concentrate on the sector S of the config-uration space C where0 < r < r − ar i − + a < r i < r i +1 − a i = 2 . . . N − r N − + a < r N < L − a , (4)which is related to all other sectors of the configurationspace by a combination of permutations and translationsof the particles, and eliminate the rod size a by the trans-formation x i = r i − ( i − a . (5)The rod coordinates x i lie in the set0 < x < x < · · · < x N < L (cid:48) , (6)where L (cid:48) = L − N a is called the unexcluded volume. TheHR Hamiltonian (2) then takes the form H = − (cid:126) m N (cid:88) i =1 ∂ ∂x i , (7)and the third condition (3), imposing that particles col-lide with each other as impenetrable elastic rods, can becorrespondingly expressed as˜Ψ( x . . . x i . . . x j . . . x N ) = 0 if x i = x j , (8)where we introduce the notationΨ( r r . . . r N ) = ˜Ψ( x , . . . x N ) , ( r r . . . r N ) ∈ S (9)to express Ψ in terms of the rod coordinates. Eigenfunc-tions of (7) satisfying the condition (8) have the form[10] ˜Ψ( x . . . x N ) = 1 √ N ! det (cid:18) e ik i x j √ L (cid:48) (cid:19) , (10)where k . . . k N are a set of quantum numberscalled quasi-wavevectors, that will be identified later.The energy eigenvalue corresponding to (10) is E { k } = (cid:126) m (cid:80) Ni =1 k i . Moreover, (10) is identically zeroif and only if any two quasi-wavevectors coincide. Thevalues of the quasi-wavevectors k . . . k N are fixed impos-ing PBC to the wavefunctions (10). Practically, imposingPBC means requiring thatΨ(0 r . . . r N ) = Ψ( L r . . . r N ) (11)for all r . . . r N . Merging (10) and (11) one finds thatPBC are satisfied if, for all quasi-wavevectors k i , the fol-lowing condition holds [13]( k i − K ) a = k i ( L − ( N − a ) − πn i + ξ B,F ( N ) , (12)where i = 1 . . . N , K = (cid:80) Ni =1 k i , n i ∈ Z is an integernumber, ξ F ( N ) = 0 and ξ B ( N ) = (cid:26) N odd π for N even . (13)Equation (12) leads easily to k i = 2 πL (cid:48) n i − L (cid:48) ξ B,F ( N ) − aKL (cid:48) . (14)Remarkably, even if the quasi-wavevectors k i are con-structed with both L and L (cid:48) , the total momentum K isan integer multiple K = 2 πL N (cid:88) i =1 n i − N ξ B,F ( N ) L (15)of πL . To summarize, the eigenfunctions of the HRHamiltonian are in one-to-one correspondence with com-binations of N integer numbers without repetition. B. Ground-state properties For a system of N Bose hard rods, the ground-statewavefunction is characterized by quasi-wavevectors k i,GS = 2 πL (cid:48) n i,GS n i,GS = − n F + ( i − 1) (16)symmetrically distributed around 0, with n F = ( N − / 2. The ground-state wavevector isnaturally K = 0. The ground-state energy reads [10, 31]: E GS = (cid:126) m (cid:18) πL (cid:48) (cid:19) n F ( n F + 1)(2 n F + 1)3 (17)In the thermodynamic limit of large system size N atconstant linear density ρ = N/L , the ground-state energyper particle converges to E ∞ = lim N →∞ E GS N = (cid:126) k F m (1 − ρa ) , (18)where k F = πρ is defined in analogy with the fermioniccase. The reduced dimensionality is responsible for the fermionization of impenetrable Bose particles: the strongrepulsion between particles mimics the Pauli exclusionprinciple [7]. In particular, the limit ρa = 0 corre-sponds to the well-known Tonks-Girardeau gas, namelythe hard-core limit of the Lieb-Liniger model [11], whereall local properties are the same as for the ideal Fermigas. At finite ρa , we can think of HRs as evolving fromthe Tonks-Girardeau gas, in that the infinitely strong re-pulsive interaction is accompanied by an increasing vol-ume exclusion. HRs are therefore a model for the super Tonks-Girardeau gas, which has been predicted and ob-served [54–58] as a highly excited and little compressible k − k F k F FIG. 1. Pictorial representation of the Lieb- I excitation with p = 6 (top) and of the Lieb- II excitation (middle) with h = 5for N = 7 HRs. The Lieb- II excitation with h = 1 is calledumklapp excitation (bottom). state of the attractive Lieb-Liniger Bose gas, in which nobound states are present.In the case of hard rods, the eigenfunctions of bothBose and Fermi systems have the same functional formin the sector S of the configuration space; away from S ,they differ from each other only by a sign associated toa permutation of the particles [7]. Therefore, the matrixelements of local operators like the density fluctuationoperator ρ q = N (cid:88) i =1 e − iqr i (19)are identical for Bose and Fermi particles. This, in par-ticular, implies that the dynamical structure factor S ( q, ω ) = (cid:90) ∞−∞ dt e iωt πN (cid:104) Ψ | e itH/ (cid:126) ρ q e − itH/ (cid:126) ρ − q | Ψ (cid:105) (20)and the static structure factor S ( q ) = (cid:90) ∞ dωS ( q, ω ) = 1 N (cid:104) Ψ | ρ q ρ − q | Ψ (cid:105) (21)are independent of the statistics. Quite usefully for thepurpose of QMC simulations in configuration space, theunnormalized bosonic ground-state wavefunction can bewritten in a Jastrow form for any a [7, 32, 37]:Ψ GS ( r · · · r N ) = (cid:89) i In the previous Subsection II B, we have recalled thatthe ground-state wavefunction of the HR system, onceexpressed in terms of the rod coordinates, has the func-tional form of a Fermi sea with renormalized coordinatesand wavevectors k i,GS , specified in terms of integer num-bers n i,GS . The excited states of the system are obtainedcreating single or multiple particle-hole pairs on top ofthis pseudo Fermi sea [59]. The simplest excitations, il-lustrated in Figure 1, consist in the creation of a singleparticle-hole pair n i,P H = n i,GS + ( p − n h,GS ) δ i,h , (23)where h ∈ { . . . N } is the index of the original quantumnumber to be modified (“hole”) and p > n F is the newinteger quantum number of index h (“particle”). Amongsingle particle-hole excitations, a role of great importancein the interpretation of S ( q, ω ) is played by the following Lieb-I n i,I = n i,GS + ( p − n F ) δ i,N (24)and Lieb-II n i,II = n i,GS + ( n F + 1 − n h,GS ) δ i,h (25)modes, which are reminiscent of the corresponding ex-citations of the Lieb-Liniger model [12]. As shown inFigure 1, in the Lieb-I excitation, a rod is taken from theFermi level h = N to some high-energy state associatedto an integer p > n F , while in the Lieb-II excitation arod is taken from a low-energy state h II to just abovethe Fermi level p = n F + 1. In both cases, in view of thecollective nature of the quasi-wavevectors k i , the excita-tion of that rod provokes a recoil of all the other rods,according to Equations (14) and (15).A simple calculation shows that the dispersion relationof the Lieb-I excitation is given by E ( N ) I ( q ) = E I ( q ) + ∆ E ( N ) I ( q ) (26)where q ≥ E I ( q ) in the thermodynamic limit reads E I ( q ) E F = 4 K L (cid:0) x + x (cid:1) , (27)In the previous equation E F = (cid:126) k F m , x = q k F andthe physical meaning of the Luttinger parameter K L =(1 − ρa ) will be elucidated in Section II D. Size effects∆ E ( N ) I ( q ) have the form∆ E ( N ) I ( q ) E F = − N K L x (1 + (1 − K L ) x ) (28)whenever ξ B,F ( N ) = 0. A similar calculation shows thatthe dispersion relation of the Lieb-II excitation is givenby E ( N ) II ( q ) = E II ( q ) + ∆ E ( N ) II ( q ) with 0 ≤ q ≤ k F and E II ( q ) E F = 4 K L (cid:0) x − x (cid:1) , ∆ E II ( q ) E F = E II ( q ) N + x N (29)Other relevant excitations are those producing super-current states [9, 12, 60] n i,SC = n i,GS + s s ∈ N (30) with momenta q = 2 sk F and excitation energies E sc ( q ) = (cid:126) m q N independent of the rod length a , and vanishingin the thermodynamic limit. Supercurrent states corre-spond to Galilean transformations of the ground statewith velocities v SC = (cid:126) k F m sN . The first supercurrentstate, in particular, is also termed umklapp excitation[9, 12, 60].In the thermodynamic limit, the particle-hole excita-tions (23) span the region ω ∗− ( q ) ≤ ω ≤ ω ∗ + ( q ) of the( q, ω ) plane, where (cid:126) ω ∗± ( q ) E F = 4 K L (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) q k F ± (cid:18) q k F (cid:19) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (31)The curves (cid:126) ω ∗± ( q ) have the same functional form ofthe ideal Fermi gas particle-hole boundaries (cid:126) ω ± ( q ) = | (cid:126) k F q/m ± (cid:126) q / m | , except for the substitution of thebare mass m with m ∗ = mK L < m , as we argued inRef. [28].The upper branch of this renormalized particle-holecontinuum coincides with the Lieb-I mode. For q ≤ k F ,its lower branch coincides with the Lieb-II mode and, for q ≥ k F , with the particle-hole excitations n i,P H = n i,GS + 1 + ( p − n N,GS − δ i,N , (32)resulting from the combination of the Lieb-I and theumklapp modes.It is worth noticing that the Lieb-II dispersion rela-tion constitutes the energy threshold for excitations for0 < q < k F . Away from this basic region, the energythreshold in the thermodynamic limit can be obtainedby a combination of inversions and shifts [52], and cor-responds to a combination of a Lieb-II mode and multi-ple umklapp excitations. To summarize, the low-energythreshold is given by (cid:126) ω th ( q ) E F = 4 K L (cid:32) q ∗ n k F − (cid:18) q ∗ n k F (cid:19) (cid:33) , (33)where 2 nk F ≤ q ≤ n + 1) k F and q ∗ n = q − nk F . Finitesize corrections are the same as in Eq. (29) [28].Remarkably, (33) corresponds also to the dispersionrelation of dark solitons of composite bosons in Yang-Gaudin gases of attractively interacting fermions in thedeep molecular regime, even though in that case themolecular scattering length is negative (corresponding toa repulsive Lieb-Liniger molecular gas) [61]. D. Comparison with Luttinger liquid theories The low-energy excitations of a broad class of inter-acting 1D systems are captured by the phenomenologicalTLL field theory [3, 47–51]. The TLL provides a uni-versal description of interacting Fermi and Bose particlesby introducing two fields, φ ( x ) and θ ( x ) representing thedensity and phase oscillations of the destruction operatorΨ( x ) (cid:39) (cid:112) ρ + ∂ x φ ( x ) e iθ ( x ) , and a quadratic low-energyHamiltonian describing the dynamics of those fields H LL = (cid:126) π (cid:90) dx (cid:18) cK L ∂ x θ ( x ) + cK L ∂ x φ ( x ) (cid:19) . (34)For Galilean-invariant systems, the sound velocity c isrelated to the positive Luttinger parameter K L through c = v F K L . The quadratic nature of (34) allows for thecalculation of correlation functions and thermodynamicproperties in terms of c and K L . Within the TLL the-ory, in the low-momentum and low-energy regime S ( q, ω )features collective phonon-like excitations ω LL ( q ) = c | q | with sound velocity c .The TLL theory has been recently extended [3, 52] be-yond the low-energy limit, where the assumption of lin-ear excitation spectrum ω LL ( q ) is not sufficient for accu-rately predicting dynamic response functions. Assumingthat, for any momentum q , S ( q, ω ) has support above alow-energy threshold ω th ( q ) and interpreting excitationswith momentum q between 2 nk F and 2 nk F + 2 k F as thecreation of mobile holes of momentum q ∗ n = q − nk F coupled with the TLL [3, 52], it is possible to show thatfor a broad class of Galilean-invariant systems S ( q, ω )features a power-law singularity close to the low-energythreshold ω th ( q ) with the following functional form: S ( q, ω ) = θ ( ω − ω th ( q ∗ n )) | ω − ω th ( q ∗ n ) | − µ n ( q ) , (35)where the exponent µ n ( q ) = 1 − (cid:18) (2 n + 1) (cid:112) K L + δ + ( q ∗ n ) + δ − ( q ∗ n )2 π (cid:19) − (cid:18) √ K L + δ + ( q ∗ n ) − δ − ( q ∗ n )2 π (cid:19) (36)is specified in terms of the phase shifts δ ± ( q )2 π = √ K L (cid:16) (cid:126) qm + ∂ω th ( q ) ∂q (cid:17) ± √ K L (cid:16) v s K L − π ∂ω th ( q ) ∂ρ (cid:17) (cid:16) ∓ ∂ω th ( q ) ∂q − v s (cid:17) (37)The only phenomenological inputs required by the non-linear TLL theory are the Luttinger parameter K L andthe low-energy threshold ω th ( q ), which in the case of hardrods are exactly known. Namely, we recall that the Lut-tinger parameter K L [50, 51] can be computed from thecompressibility κ − S = ρ ∂∂ρ (cid:18) ρ ∂E ∞ ∂ρ (cid:19) (38)through the formula mK L = (cid:126) k F ρ κ S . The resultingexact expression, K L = (1 − ρa ) [32], provides a Lut-tinger parameter always smaller than 1, and convergingtowards 0 as the excluded volume N a converges towards L . Notice that both Lieb-I and Lieb-II dispersions ap-proach q = 0 with slope equal to the sound velocity -4-3-2-1 0 1 0 1 2 µ ( q ) q (2k F ) FIG. 2. (color online) Nonlinear Luttinger theory exponents µ ( q ) for hard rods at the studied densities from ρa = 0 . ρa = 0 . 900 (red solid line). Momenta aremeasured in units of 2 k F = 2 πρ . The ρa → ρa ≥ . Q between 2 k F and 4 k F . c = v F K L > v F . The low-energy threshold, Eq. (33),has been described in the previous Section.Knowledge of these two quantities permits to computethe exponent µ n ( q ) exactly from (36) and (37). We find[28] µ n ( q ) = − q − n ) (˜ q − ( n + 1)) , (39)with ˜ q = qa π . In Fig. 2 we show the power-law exponentsfor momenta 0 < q < k F and different densities. Noticethat the functional form of (39) is that of a sequence ofparabola arcs, intersecting null values at the special mo-menta Q n = n π/a , with integer n < ρa/ (1 − ρa ). Suchmomenta, even for larger n , have already been recognizedto be special [32], in that they admit the exact calcula-tion of S ( Q n ) and S ( Q n , ω ). Namely, for those specialmomenta, the HRs at density ρ behave as an ideal Fermigas with increased density ρ (cid:48) = ρ/ (1 − ρa ).It is worth pointing out that TLL theories have limitsof applicability, and thus do not exhaust our understand-ing of 1D substances [27, 62, 63]. The investigation ofdynamical properties like S ( q, ω ) beyond the limits of ap-plicability of Luttinger liquid theories, where the Physicsis non-universal, typically requires numerical calculationsor QMC simulations [28, 62, 64]. III. METHODS In the present work, the zero-temperature dynamicalstructure factor of a system of Bose HRs is calculatedusing the exact Path Integral Ground State (PIGS) QMCmethod to compute imaginary-time correlation functionsof the density fluctuation operator, and the state-of-the-art Genetic Inversion via Falsification of Theories (GIFT)analytic continuation method to extract the dynamicalstructure factor.This approach, which we briefly review in this Section,has provided robust calculations of dynamical structurefactors for several non-integrable systems like 1D, 2D and3D He atoms [28, 41, 65, 66] and hard spheres [67].The PIGS method is a projection technique in imagi-nary time that, starting from a trial wavefunction Ψ T ( R ),where R = ( r . . . r N ) ∈ C denotes a set of spatial coordi-nates of the N particles, projects it onto the ground-statewavefunction Ψ GS ( R ) after evolution over a sufficientlylong imaginary-time interval τ [38–40]. In typical situa-tions, the functional form of Ψ T ( R ) is guessed combiningphysical intuition and mathematical arguments based onthe theory of stochastic processes [68]. Ψ T ( R ) is thenspecified by one or more free parameters, that are chosenusing suitable optimization algorithms [69–71].In the case of HRs, knowledge of the exact ground-state wavefunction (22) makes the projection of a trialwavefunction Ψ T ( R ) approximating the ground state ofthe system unnecessary. However, the PIGS method canbe used to give unbiased estimates of the density-densitycorrelator F ( q, τ ) = (cid:104) Ψ GS | e τH ρ − q e − τH ρ q | Ψ GS (cid:105) == (cid:82) dR M dR p ( R M , R ) ρ − q ( R M ) ρ q ( R ) (cid:82) dR M dR p ( R M , R ) , (40)with p ( R M , R ) = Ψ GS ( R M ) G ( R M , R ; τ )Ψ GS ( R ) and G ( R M , R ; τ ) = (cid:104) R M | e − τH | R (cid:105) .The propagator G ( R (cid:48) , R ; τ ) is in general not known,but suitable approximate expressions are available forsmall δτ = τ /M , where M is a large integer number.Using one of these expressions in place of the exact propa-gator is the only approximation characterizing the calcu-lations of the present work. The method is exact though,since this approximation affects the computed expecta-tion values to an extent which is below their statisticaluncertainty and such regime is always attainable by tak-ing δτ sufficiently small. Then, the convolution formulapermits to express G ( R M , R ; τ ) as G ( R M , R ; τ ) = (cid:90) dR M − . . . dR M − (cid:89) i =0 G ( R i +1 , R i ; δτ ) , (41)whence the PIGS estimator of F ( q, τ ) takes the form F ( q, τ ) = (cid:82) dX p ( X ) ρ − q ( R M ) ρ q ( R ) (cid:82) dX p ( X ) . (42)In (42), X = ( R . . . R M ) denotes a path in the configu-ration space C of the system, and p ( X ) = Ψ GS ( R M ) M − (cid:89) i =0 G ( R i +1 , R i ; δτ )Ψ GS ( R ) (43) can be efficiently sampled using the Metropolis algorithm[72]. In the present work, we have employed the pair-product approximation [73] to express the propagatorrelative to a small time step δτ as G ( R, R (cid:48) ; δτ ) = N (cid:89) i =1 G ( r i , r (cid:48) i ; δτ ) N (cid:89) i 0) = S ( q ). For finite valuesof τ , F ( q, τ ) is instead related to S ( q, ω ) by the Laplacetransform F ( q, τ ) = (cid:90) ∞ dω e − τω S ( q, ω ) . (47)Equation (47) should be inverted in order to determine S ( q, ω ) from F ( q, τ ). However, it is well-known that suchinverse problem is ill-posed , in the sense that many differ-ent trial dynamical structure factors, ranging from fea-tureless to rich-in-structure distributions, have a forwardLaplace transform which is compatible with the QMC re-sults for F ( q, τ ): there is not enough information to finda unique solution of (47) [41, 76–81]. Different method-ologies have been used to extract real-frequency responsefunctions from imaginary-time correlators; in the presentwork, we rely on the Genetic Inversion via Falsificationof Theories (GIFT) method [41]. The aim of GIFT is tocollect a large collection of such dynamical structure fac-tors in order to discern the presence of common features(e.g. support, peak positions, intensities and widths).The GIFT method has been applied to the study of liq-uid He [28, 82, 83], 3D hard spheres [67], 2D YukawaBosons [84], liquid He [65], 2D soft disks [85, 86], the2D Hubbard model [87] and 1D soft rods [88], in all casesproviding very accurate reconstructions of S ( q, ω ) or thesingle particle spectral function. Recently [28], we haveshown that in 1D, when ω th ( q ) is known, also the shapeclose to the frequency threshold can be approximatelyinferred. This is the reason why in subsection II C weinsisted on the calculation of ω th ( q ) for a finite system,which is a most useful quantity in our approach. Detailsof the GIFT method can be found in [28, 41]. As in [28], -0.0050.0000.0050.0100.0150.0200.025 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 | - S N ( q ) / S ∞ ( q ) | q (2k F ) N=100N=50N=25 FIG. 3. (color online) We quantify the finite-size effectson S ( q ) at the representative density ρa = 0 . 700 computingthe relative error | S N ( q ) − S ∞ ( q ) | /S ∞ ( q ), where S ∞ ( q ) =lim N →∞ S N ( q ) is extrapolated. Away from the points q =2 k F , k F , the relative error is below 1% for N = 50 particles. we have used genetic operators which are able to betterdescribe broad features typical of 1D systems. Moreover,the set of discrete frequencies of the model spectral func-tions used in the algorithm has been extended to non-equispaced frequencies, in order to better describe theregions where most of the weight accumulates. IV. RESULTS We computed F ( q, τ ) for systems of N = 50 hard rodsat densities ρa listed in Table I, and wavevectors q ≤ k F .Before describing our results on the dynamical structurefactor, we demonstrate the accuracy of our calculationsby analyzing in detail finite-size effects on static proper-ties which have already been studied in Ref. [32]. A. Assessment of accuracy Results are affected by very weak finite-size effects, andthus are well representative of the thermodynamic limit.For example, (17) and (18) yield the following finite-size ρa K L K L studied in the presentwork. corrections to the ground-state energy E GS N = E ∞ (cid:18) − N (cid:19) (48)whence E GS = 0 . E ∞ . To further assess the finite-size effects on our results, in Figure 3 we compute thestatic structure factor of N = 50, 100 hard rods at ρa =0 . 700 using the VMC method, which is an exact methodwhen the exact (ground state) wave function is known,as in this case.At q = 2 k F , k F , the static structure factor displayspeaks of diverging weight as predicted by the TLL theory[32, 89]: S (2 mk F ) = S smooth (2 mk F ) + S peak (2 mk F ) == S smooth (2 mk F ) + C m N − m K L (49)Away from those points, the VMC estimates of the staticstructure factor are compatible with each other and withthe extrapolation of S ( q ) to the thermodynamic limit,within the error bars of the simulations, reflecting theweakness of finite-size effects. Moreover, in Figure 4 weshow that static structure factors of N = 50 rods permitto compute the Luttinger parameter K L without appre-ciable finite-size effects.The same favorable behavior is exhibited by F ( q, τ ). InFigure 5, we show F ( q, τ ) for N = 50 , 100 rods at the rep-resentative density ρa = 0 . q = 2 k F thetwo systems have statistically compatible F ( q, τ ), con-firming the modest entity of finite-size effects. B. Dynamical structure factors In Figure 6 we show the dynamical structure factor atdensities ranging from ρa = 0 . 005 to 0 . S ( q ) q (2k F ) 0.0050.0770.3210.4710.6420.7000.900 -0.02-0.01 0 0.01 0.02 0 0.5 1 - K L , N / K L , ∞ ρ a FIG. 4. (color online) Static structure factor of N = 50 rodsat ρa = 0 . . . . . . . S ( q ) (cid:39) K L q/ (2 k F ). In all cases, the relative error is below 1%. F ( q , τ ) τ [1/E F ] (a)(b)(c)(d)(e)(f)(g) FIG. 5. (color online) F ( q, τ ) of N = 50 , 100 (blue circles, redopen diamonds) rods at ρa = 0 . 700 for q k F = 0 . 02, 0 . 22, 0 . . 82, 0 . 96, 1 . 00, 1 . 08 (a to g). At q = 2 k F we observe strongfinite-size effects, originating from the quasi-Bragg peaks in S ( q ). Away from q = 2 k F , in the relevant imaginary-timeinterval τ E F ≤ . F ( q, τ ) are in satisfactory agree-ment. studied densities, for momentum q < k F , S ( q, ω ) hasmost of the spectral weight inside the particle-hole band ω ∗− ( q ) ≤ ω ≤ ω ∗ + ( q ) spanned by single particle-hole exci-tations, Eq. (31). Contributions from multiple particle-hole excitations become relevant at high densities andmomenta.At the lowest density, panel (a), the spectral weight isbroadly distributed inside the particle-hole band, show-ing a behavior reminiscent of the Tonks-Girardeau modelof impenetrable point-like bosons, to which the HR modelreduces in the ρa → q < k F it predictsa power-law behavior for S ( q, ω ), with an exponent (39)slightly larger than zero. This prediction is consistentwith the spectrum in panel (a), showing a weak concen-tration of spectral weight close to the low-energy thresh-old for q < k F . For q > k F , the support of S ( q, ω )departs from the low-energy threshold as it can be seenin panel (a), where the spectral weight remains concen-trated inside the particle-hole band for all q . Correspond-ingly, for 2 k F < q < k F , the nonlinear TLL predicts alarge negative exponent, suggesting the absence of spec-tral weight in the proximity of the low-energy threshold.A similar behavior is shown at density ρa = 0 . q < k F .This is again in agreement with the nonlinear TLL the-ory, predicting a larger negative exponent µ ( q ).The dynamical structure factors in panels (a), (b)are also in qualitative agreement with numeric calcu-lations for the super Tonks-Girardeau gas, for which0 . < ∼ K L < ρa (cid:39) . (cid:126) ω F A ( q ) = (cid:126) q mS ( q ) (50)breaks down beyond q k F (cid:39) . 1. Interestingly, around q k F (cid:39) . ω th ( q ) with a linearfunction of q ceases to be adequate. The simultaneous ap-pearance of nonlinear terms in ω th ( q ) and corrections tothe Feynman approximation in S ( q, ω ) are in fact deeplyrelated phenomena, as explained by the nonlinear TLLtheory.When K L < / 2, Eq. (49) indicates that a peak man-ifests in the static structure factor at q = 2 k F . Thischange in behavior is also reflected in S ( q, ω ), as shownin panels (c), (d). The spectral weight concentrates closeto the lower branch ω ∗− ( q ) of the particle-hole band, in aregion of dense spectral weight that we call lower modefollowing [28], where a similar behavior was observed in1D He at high density. In both HRs and He, abovethe lower mode stretches a high-energy structure gath-ering a smaller fraction of spectral weight. Such high-energy structure has a minimum at q = 2 k F close to thefree-particle energy E = 4 E F , and is symmetric around q = 2 k F , (see panel (d) in Figure 6).Panel (d) also shows that, as K L decreases below 1 / S ( q, ω ) extends below ω ∗− ( q ) for q > k F ,still remaining above ω th ( q ). In the high-density regime ρa ≥ . S ( q, ω )as it can be read from panel (e), where the shape of S ( q, ω ) changes considerably for 2 k F < q < k F .For q < . k F , the spectral weight concentrates ina narrow region of the momentum-energy plane thatgradually departs from the low-energy threshold. For2 . k F < q < . k F , S ( q, ω ) suddenly and considerablybroadens and flattens. Finally, for 3 . k F < q < k F ,the spectral weight again concentrates close to the low-energy threshold. This highly non-trivial behavior is inqualitative agreement with the nonlinear Luttinger liq-uid theory, predicting a negative exponent for q < Q ,where Q k F = aρ = 1 . q = Q the non-linearLuttinger liquid theory predicts a flat dynamical struc-ture factor close to the low-energy threshold, in agree-ment with an exact prediction by F. Mazzanti et al. [32]and with our observations for the 1D HR system, butalso for a 1D system of He atoms[28]. Beyond Q thenon-linear Luttinger liquid theory predicts a positive ex-ponent and we observe the spectral weight concentratingclose to the low-energy threshold, as for q < k F . No-tice that, although the condition for having quasi-Braggpeaks is K L < / 2, the first special momentum with flat 0 0.5 1 1.5 2 0 2 4 6 8 10 12 14 ω ( E F / − h ) q (2k F ) 0 0.5 1 1.5 2 0 2 4 6 8 10 12 14 0 0.5 1 1.5 2 0 2 4 6 8 10 12 14 0 0.5 1 1.5 2 2.5 3 0 2 4 6 8 10 12 14 0 0.5 1 1.5 2 2.5 3 0 5 10 15 20 25 30 35 0 0.5 1 1.5 2 2.5 3 3.5 4 0 50 100 150 200 250 300 (a) ρ a = 0.005K L = 0.990 ω th ω ∗ ± ω FA (b) ρ a = 0.077K L = 0.852(c) ρ a = 0.321K L = 0.461 (d) ρ a = 0.471K L = 0.280(e) ρ a = 0.642K L = 0.128 (f) ρ a = 0.900K L = 0.010 FIG. 6. (color online) Color map of the dynamic structure factor, in units of (cid:126) /E F , at ρa = 0 . . . . . . 900 (left to right, top to bottom). Momenta are measured in units of 2 k F = 2 πρ and energies in units of E F . The low-energythreshold (blue dotted line), the branches ω ∗± ( q ) (purple solid lines) and the Feynman approximation for the excitation spectrum(green dashed line) are drawn for comparison. Panels (e), (f) also show the special wavevectors Q n (red arrows). 0 2 4 6 8 10 12 14 16 18q (2k F ) 0 200 400 600 800 1000 ω ( E F / − h ) ρ a = 0.900K L = 0.010 ω R ± ω ∗ ± ω th ω FA FIG. 7. (color online) Color map of the dynamic structurefactor at ρa = 0 . Q n (redarrows) are shown by long (red) arrows. The role of boththe HR renormalized particle-hole frequencies (solid curves)and of the reduced-volume ideal Fermi gas band (dot-dashedcurves) are shown. spectrum appears only for ρa > / 2, namely K L < / Q = 2 π/a < πρ .It is well-known [31, 32] that, in the high-densityregime K L (cid:28) 1, HRs show up a packing order leadingto a quasi-solid phase, crystallization being prohibitedby the reduced dimensionality and by the range of theinteraction. The emergence of the quasi-solid phase issignaled by the peaks of the static structure factor, thatapproach the linear growth with the system size, peculiarprerogative of the Bragg peaks, only in the ρa → S ( q, ω ) is concerned, the increase of the den-sity transfers spectral weight to the low-energy thresholdat wave vectors q = 2 k F , k F . We interpret this phe-nomenon as revealing that S ( q, ω ) is gradually approach-ing translational invariance q → q +2 k F in the variable q .This conjecture is corroborated by the observation that,as the density is further increased, the exponent µ n ( q )pointwise converges to the periodic function: µ n ( q ) = − (cid:18) q k F − n (cid:19) (cid:18) q k F − ( n + 1) (cid:19) (51)with 2 nk F ≤ q ≤ n + 1) k F (illustrated in Figure 2).In this respect, it is interesting to observe panel (f) ofFigure 6, showing the dynamical structure factor of HRsat ρa = 0 . q < k F , the spectral weight al-most always concentrates around ω th ( q ), except in thesmall ranges of wavevectors 2 nk F < q < Q n . This makes S ( q, ω ) resemble the dispersion relation of longitudinalphonons of a monoatomic chain, in this range of mo-menta. However, the non-commensurability of Q with2 k F renders the spectrum only quasi-periodic, a behav-ior which is more and more manifest at higher momenta.To analyze this intriguing regime in more detail, we have reconstructed the spectra at ρa = 0 . 900 up to q = 38 k F .The results are shown in Fig. 7 and indicate a crucialrole of the reduced-size ideal Fermi gas in drawing large-scale momentum and frequency boundaries for the HRspectrum. If we define the density ρ (cid:48) = ρ/ (1 − ρa ), weobserve that above q = 20 k F = 2 πρ (cid:48) , namely twice theFermi momentum of the reduced-size IFG, the spectrumis never peaked along the low-energy HR threshold, how-ever it presents a stripe structure, repeating the Lieb-Imode plus multiple umklapp excitations, and becomingagain flat at momenta Q n . At the level of accuracy ofour GIFT reconstructions, the stripes are bounded bythe particle-hole band of the reduced-size IFG (cid:126) ω R ± ( q ) = (cid:126) m (cid:12)(cid:12) πρ (cid:48) q ± q (cid:12)(cid:12) . (52)To corroborate this observation, we find that the specialmomenta Q n also analytically correspond to the cross-ings of the lower reduced-size IFG boundary and the HRthreshold (for q < πρ (cid:48) ) or the HR repeated Lieb-I modes( q > πρ (cid:48) ). It is tempting to conclude that the HR spec-trum can be almost completely described by the synergyof two rescaled ideal Fermi gases: one with the same den-sity, but renormalized mass m ∗ = mK L , the other withthe same mass, but increased density ρ (cid:48) = ρ/K / L .Only in the unphysical ρa → S ( q, ω )would attain the translational invariance observed for in-stance in half-filled Hubbard chains with strong on-siterepulsion [90, 91].Finally, in discussing Figure 6, we stressed in multi-ple occasions that low-energy properties of S ( q, ω ) arecaptured by the nonlinear TLL theory. Following the ap-proach of [28], in Figures 8 and 9 we show that the agree-ment is quantitative, focusing on density ρa = 0 . 700 andtwo momenta close to Q , representative of negative andpositive power-law exponents, and fitting multiple spec-tral reconstructions with the functional form given byEq. (35). The agreement with the analytical expressionfor the power-law exponent (39) is good, even though theaccuracy is strongly dependent on the quality of the origi-nal F ( q, τ ) and it is particularly delicate to fit the spectrawhen the spectrum goes to zero close to the threshold. V. CONCLUSIONS AND OUTLOOKS We have computed the zero-temperature dynamicalstructure factor of one-dimensional hard rods by meansof state-of-the-art QMC and analytic continuation tech-niques. By increasing the rod length, the dynamicalstructure factor reveals a transition from the Tonks-Girardeau gas, to a super Tonks-Girardeau regime andfinally to a quasi-solid regime.The low-energy properties of the dynamical structurefactor are in qualitative agreement with the nonlinearLL theory. However, the methodology provides a quan-titative estimation of the dynamical structure factor alsoin the high-energy regime, lying beyond the reach of LL1 -6 -5 -4 -3 -2 -1 F S ( q , ω ) ( − h / E F ) ω - ω th (E F / − h )power-law fit FIG. 8. (color online) Power-law fit of reconstructed spectrafor ρa = 0 . 700 and q = 2 . k F . The fitted exponent is µ = − . µ = − . -4 -3 -2 -1 F S ( q , ω ) ( − h / E F ) ω - ω th (E F / − h )power-law fit FIG. 9. (color online) Power-law fit of reconstructed spectrafor ρa = 0 . 700 and q = 3 . k F . The fitted exponent is µ =0 . µ = 0 . theories. Our study reveals strong similarities betweenthe dynamical structure factor of HRs and 1D He atlinear densities ρ ≥ . 150 ˚A − [28], extending to thehigh-energy regime (well above the low-energy thresh-old). In particular, both systems show a flat dynami-cal structure factor in correspondence of the wavevectors Q n , in agreement with a previous theoretical prediction[32], and feature a high-energy structure overhanging thelower mode around the umklapp point q = 2 k F . We havealso unveiled a peculiar structure of the spectrum in thehigh-density and high-momentum regime, which can bedescribed in terms of the particle-hole boundaries of tworenormalized ideal Fermi gases.At this point we want to remark that an intrigu-ing feature of 1D He (and arguably of all 1D quan-tum liquids which admit a two-body bound state) isthat its Luttinger parameter K L spans all positive val-ues 0 < K L < ∞ as the linear density is increased[28], not only the K L ≤ ≤ K L ≤ ∞ by tuning the interaction strength. Fi-nally also the Calogero-Sutherland model reproduces allpossible K L , however by tuning interaction only [93].We hope the present work will encourage further ex-perimental research in 1D systems with volume exclusioneffects, and theoretical investigation of the HR model.Our results may be relevant also for the linear responsedynamics of resonant Rydberg gases in 1D configurations[94]. On the theoretical side, possible directions for fur-ther developments might be the alternative calculationsof S ( q, ω ), based e.g. on the VMC evaluation of thematrix elements of ρ q , and/or the calculation of finite-temperature equilibrium and dynamical properties. ACKNOWLEDGMENTS We acknowledge useful discussions with G. As-trakharchik. We thank M. Panfil and co-authors forproviding us with their data on the super Tonks-Girardeau gas [54]. We acknowledge the CINECA awardsIsC29-SOFTDYN, IsC32-EleMEnt and CINECA andRegione Lombardia LISA award LI05p-PUMAS, for theavailability of high-performance computing resources andsupport. We also acknowledge computing support fromthe computational facilities at the College of William andMary and at the Physics Department of the Universityof Milan. M.M. and E.V. acknowledge support from theSimons Foundation and NSF (Grant no. DMR-1409510).M.R. acknowledges the EU QUIC project for fundings. [1] T. Giamarchi, Quantum Physics in One Dimension , Ox-ford University Press (2004)[2] M.A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac andM. Rigol, Rev. Mod. Phys. , 1405 (2011)[3] A. Imambekov, T. L. Schmidt and L. I. Glazman, Rev.Mod. Phys. , 1253 (2012)[4] P.C. Hohenberg, Phys. Rev. , 383 (1967) [5] N. D. Mermin and H. Wagner, Phys. Rev. Lett. , 1133(1966)[6] S. Coleman, Commun. Math. Phys. , 259 (1973)[7] M. Girardeau, J. Math. Phys. , 516 (1960)[8] L. Pitaevskii and S. Stringari, Phys. Rev. B , 10915(1993)[9] G. F. Giuliani and G. Vignale, Quantum Theory of the Electron Liquid , Cambridge University Press (2005)[10] T. Nagamiya, Proc. Phys. Math. Soc. Jpn. , 705 (1940)[11] E. H. Lieb and W. Liniger, Phys. Rev. , 1605 (1963)[12] E.H. Lieb, Phys. Rev. , 1616 (1963).[13] B. Sutherland, J. Math. Phys. , 251 (1971)[14] B. Sutherland Phys. Rev. A , 2019 (1971)[15] B. Sutherland, Phys. Rev. B , 6689 (1988)[16] P . Paredes, Nature , 277 (2004)[17] T. Kinoshita, T. Wenger and D. S. Weiss, Science ,1125 (2004)[18] E. Haller, M. Gustavss, M. J. Mark, J. G. Danzl, R. Hart,G. Pupillo and H. C. Negerl, Science , 1224 (2009)[19] H. P. Stimming, N. J. Mauser, J. Schmiedmayer and I.E. Mazets, Phys. Rev. Lett. , 015301 (2010)[20] M. Kuhnert, R. Geiger, T. Langen, M. Gring, B. Rauer,T. Kitagawa, E. Demler, D. Adu Smith and J. Schmied-mayer, Phys. Rev. Lett. , 090405 (2013)[21] S. Hofferberth, I. Lesanovsky, B. Fischer, T. Schumm andJ. Schmiedmayer, Nature , 324 (2007)[22] T. Langen, R. Geiger, M. Kuhnert, B. Rauer and J.Schmiedmayer, Nature Physics , 640 (2013)[23] E. Witkowska, P. Deuar, M. Gajda and K. Rzazewski, Phys. Rev. Lett. , 135301 (2011)[24] T. Karpiuk, P. Deuar, P. Bienias, E. Witkowska, K.Pawlowski, M. Gajda, K. Rzazewski and M. Brewczyk, Phys. Rev. Lett. , 205302 (2012)[25] V. Guarrera, D. Muth, R. Labouvie, A. Vogler, G. Baron-tini, M. Fleischhauer and H. Ott, Phys. Rev. A ,021601 (2012)[26] J.P. Ronzheimer, M. Schreiber,S. Braun,S.S. Hodg-man,S. Langer,I.P. McCulloch,F. Heidrich-Meisner,I.Bloch,U. Schneider, Phys. Rev. Lett. , 205301 (2013)[27] N. Fabbri, M. Panfil, D. Cl´ement, L. Fallani, M. Inguscio,C. Fort, and J. S. Caux Phys. Rev. A , 043617 (2015)[28] G. Bertaina, M. Motta, M. Rossi, E. Vitali and D. E.Galli, Phys. Rev. Lett. , 135302 (2016)[29] J. V. Pearce, M. A. Adams, O. E. Vilches, M. R. Johnsonand H. R. Glyde, Phys. Rev. Lett. , 185302 (2005)[30] M. Mercedes Calbi, M. W. Cole, S. M. Gatica, M. J.Bojan and G. Stan, Rev. Mod. Phys. , 857 (2001)[31] F. Mazzanti, G. E. Astrakharchik, J. Boronat and J. Ca-sulleras, Phys. Rev. A , 043632 (2008)[32] F. Mazzanti, G. E. Astrakharchik, J. Boronat and J. Ca-sulleras, Phys. Rev. Lett. , 020401 (2008)[33] J. D. Van der Waals, The equation of state for gases andliquids , Nobel Lectures in Physics 254 (1910)[34] J. H. Jeans, The Dynamical Theory of Gases , CambridgeUniversity Press (1916)[35] L. Tonks, Phys. Rev. , 955 (1936)[36] B. R. A. Nijboer and L. Van Hove, Phys. Rev. , 777(1952)[37] E. Krotscheck, M.D. Miller, and J. Wojdylo, Phys. Rev.B , 13028 (1999).[38] A. Sarsa, K. E. Schmidt, and W. R. Magro, J. Chem.Phys. , 1366 (2000)[39] D. E. Galli and L. Reatto, Mol. Phys. , 1697 (2003)[40] M. Rossi, M. Nava, L. Reatto, and D.E. Galli, J. Chem.Phys. , 154108 (2009)[41] E. Vitali, M. Rossi, L. Reatto and D. E. Galli, Phys. Rev.B , 174510 (2010)[42] R.A. Cowley and A.D.B. Woods, Can. J. Phys. , 177(1971).[43] K. Beauvois, C.E. Campbell, J. Dawidowski, B. F˚ak, H.Godfrin, E. Krotscheck, H.-J. Lauter, T. Lichtenegger, J. Ollivier, and A. Sultan, Phys. Rev. B , 024504 (2016).[44] R. Ozeri, N. Katz, J. Steinhauer, and N. Davidson, Rev.Mod. Phys. , 187 (2005).[45] L.-C. Ha, L.W. Clark, C.V. Parker, B.M. Anderson, andC. Chin, Phys. Rev. Lett. , 055301 (2015).[46] R. Landig, F. Brennecke, R. Mottl, T. Donner, and T.Esslinger, Nat. Commun. , (2015).[47] S. Tomonaga, Prog. Theor. Phys. , 544 (1950)[48] J. M. Luttinger, J. Math. Phys. , 1154 (1963)[49] D. C. Mattis and E. H. Lieb, J. Math. Phys. , 304 (1965)[50] F. D. M. Haldane, Phys. Rev. Lett. , 1840 (1981)[51] F. D. M. Haldane, Phys. Rev. Lett. , 569 (1982)[52] A. Imambekov and L. I. Glazman, Phys. Rev. Lett. ,126405 (2009)[53] H. Bethe, Zeit. Phys. , 205 (1931)[54] M. Panfil, J. De Nardis and J.-S. Caux, Phys. Rev. Lett. , 125302 (2013)[55] G. E. Astrakharchik, J. Boronat, J. Casulleras and S.Giorgini, Phys. Rev. Lett. , 190407 (2005)[56] M. T. Batchelor, M. Bortz, X. W. Guan and N. Oelkers, J. Stat. Mech. L10001 (2005)[57] E. Tempfli, S. Z¨ollner and P. Schmelcher, New J. Phys. , 103021 (2008)[58] E. Haller, M. Gustavsson, M.J. Mark, J.G. Danzl, R.Hart, G. Pupillo, and H.-C. Ngerl, Science , 1224(2009).[59] A.H. Castro Neto, H.Q. Lin, Y.-H. Chen, and J.M.P.Carmelo, Phys. Rev. B 50 , 14032 (1994).[60] A. Y. Cherny, J. S. Caux and J. Brand Frontiers ofPhysics , 54 (2012)[61] S.S. Shamailov and J. Brand, New J. Phys. , 75004(2016)[62] J. S. Caux and P. Calabrese, Phys. Rev. A , 031605(2006)[63] F. Meinert, M. Panfil, M.J. Mark, K. Lauber, J.-S. Cauxand H.-C. N¨agerl, Phys. Rev. Lett. , 085301 (2015)[64] S. De Palo, E. Orignac, R. Citro, and M.L. Chiofalo, Phys. Rev. B , 212101 (2008)[65] M. Nava, D.E. Galli, S. Moroni and E. Vitali, Phys. Rev.B , 144506 (2013)[66] F. Arrigoni, E. Vitali, D. E. Galli and L. Reatto, LowTemp. Phys. , 793 (2013)[67] R. Rota, F. Tramonto, D.E. Galli and S. Giorgini, Phys.Rev. B , 214505 (2013)[68] M. Holzmann, D. M. Ceperley, C. Pierleoni and K. Esler, Phys. Rev. E , 046707 (2003)[69] M. H. Kalos and P. A. Whitlock Quantum Monte Carlo ,in Monte Carlo Methods , Wiley (1986)[70] J. Toulouse and C. J. Umrigar, J. Chem. Phys. ,084102 (2007)[71] M. Motta, G. Bertaina, D. E. Galli and E. Vitali, Comp.Phys. Comm. , 62 (2015)[72] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A.H. Teller and E. Teller, J. Chem. Phys. Rev. Mod. Phys. , 279 (1995)[74] J. Cao and B.J. Berne, J. Chem. Phys. , 2382 (1992)[75] F. De Soto and M.C. Gordillo, Phys. Rev. A , 013607(2012)[76] A. Tarantola, Nature Physics , 492 (2006)[77] R. N. Silver, J.E. Gubernatis, D. S. Sivia and M. Jarrell, Phys. Rev. Lett. , 496 (1990)[78] M. Jarrell and J.E. Gubernatis, Phys. Rep. , 133(1996)[79] M. Boninsegni, and D.M. Ceperley, J. Low Temp. Phys. , 339 (1996)[80] A. Roggero, F. Pederiva and G. Orlandini Phys. Rev. B , 094302 (2013)[81] A.S. Mishchenko, N.V. Prokofev, A. Sakamoto, and B.V.Svistunov, Phys. Rev. B , 6317 (2000).[82] M. Rossi, E. Vitali, L. Reatto and D.E. Galli, Phys. Rev.B , 014525 (2012)[83] M. Nava, D.E. Galli, M.W. Cole, L. Reatto, J. LowTemp. Phys. , 699 (2013)[84] S. Molinelli, D. E. Galli , L. Reatto and M. Motta, J. LowTemp. Phys. DOI: 10.1007/s10909-016-1628-3 (2016)[85] S. Saccani, S. Moroni, E. Vitali, and M. Boninsegni, Mol.Phys. , 2807 (2011).[86] S. Saccani, S. Moroni, and M. Boninsegni, Phys. Rev.Lett. , 175301 (2012).[87] E. Vitali, H. Shi, M. Qin, and S. Zhang, Phys. Rev. B , 085140 (2016).[88] M. Teruzzi, D.E. Galli, and G. Bertaina,arXiv:1607.05308 (2016).[89] G. E. Astrakharchik and J. Boronat, Phys. Rev. B ,235439 (2014)[90] G. Roux, A. Minguzzi and T. Roscilde, New J. Phys. ,055003 (2013)[91] R. G. Pereira, K. Penc, S. R. White, P. D. Sacramentoand J. M. P. Carmelo, Phys. Rev. B , 165132 (2012)[92] R. Citro, E. Orignac, S. De Palo, and M.L. Chiofalo, Phys. Rev. A , 051602 (2007)[93] G.E. Astrakharchik, D.M. Gangardt, Y.E. Lozovik, andI.A. Sorokin, Phys. Rev. E , 021105 (2006)[94] D. Petrosyan, M. H¨oning, and M. Fleischhauer, Phys.Rev. A87