Introduction

Ultracold helium constitutes one of the most actively investigated quantum systems and has been of central relevance for our understanding of important physical concepts such as superfluidity1 and Bose–Einstein condensation2. Due to its nature as a strongly correlated quantum liquid, helium exhibits an intricate interplay of non-ideality effects, quantum statistics, and thermal excitations. Naturally, an accurate description of physical effects such as the lambda phase transition of \(^4{\text {He}}\) must capture all of these effects simultaneously—a challenging task beyond simple mean-field models and perturbative approaches.

This challenge was met by Feynman3 in terms of the path integral formalism that exactly maps the interacting quantum system of interest onto an effective classical system of ring polymers4. Specifically, this quest for an accurate description of helium5 has given rise to the widely used path integral Monte Carlo (PIMC) simulation method6,7,8, one of the most successful tools in statistical physics, quantum chemistry, and related disciplines (We note that Monte Carlo methods in general are applied in a gamut of different contexts, like solid state physics9,10 or the investigation of magnetic properties11,12,13,14,15). Recently the worm algorithm16,17 elegantly addressed the challenge of sampling the permutation space ergodically. This is required to take into account the effect of quantum statistics, which is nontrivial because of the strong short-range repulsion between two He atoms.

The PIMC method gives straightforward access to important physical observables like the superfluid fraction18, the momentum distribution1, and the static structure factor, which has resulted in excellent agreement between theory and experiments for \(^4{\text {He}}\); see the review by Ceperley1 for details. In addition, PIMC simulations can be used as the starting point for an analytic continuation19 that provides access to the dynamic structure factor \(S({\mathbf {q}},\omega )\)20,21,22—a key quantity in neutron scattering experiments23,24,25,26. In particular, PIMC-based data for \(S({\mathbf {q}},\omega )\) have given important insights into the connection between superfluidity and roton-like quasi-particle excitations.

In stark contrast, PIMC simulations of \(^3{\text {He}}\) are substantially hampered by the notorious fermion sign problem27,28,29, which leads to an exponential increase of the computation time with increasing the system size N or decreasing the temperature T. Therefore, Ceperley has introduced the uncontrolled fixed-node approximation into fermionic PIMC simulations27 and derived the first results for \(^3{\text {He}}\). This investigation was restricted to the total energy, and the agreement to experimental data30 was inconclusive. Few PIMC studies of \(^3{\text {He}}\) have been published since31. To our knowledge, no data have been presented for either the structural properties or the spectrum of collective excitations.

This is unfortunate, as ultracold \(^3{\text {He}}\) offers a potential wealth of interesting physical effects such as a superfluid phase transition at \(T\lesssim 2.5\;{\text {mK}}\)32 that is triggered by the formation of Cooper pairs. Furthermore, it has been recently demonstrated33,34 that two-dimensional \(^3\text {He}\) exhibits a rich phonon–maxon–roton dispersion relation that phenomenologically resembles the well-known dispersion of \(^4{\text {He}}\). We also note that bulk \(^3{\text {He}}\) is notoriously difficult to study with laboratory experiments35. A thorough theoretical approach is thus indispensable to understand the underlying physical mechanisms. Furthermore, progress has been made in characterizing \(^3{\text {He}}\)–\(^4{\text {He}}\) mixtures with PIMC simulations36,37,38 but some disagreements regarding the kinetic energy remain between theoretical and experimental results.

In this work, we remedy this unsatisfactory situation by carrying out extensive direct PIMC simulations of normal liquid \(^3{\text {He}}\) without any nodal restrictions. Therefore, our simulations are free of uncontrolled approximations, but computationally extremely costly when the temperature is decreased, cf. the discussion of Fig. 1 below. This allows us to present highly accurate results for the temperature dependence of important properties such as the static structure factor \(S({\mathbf {q}})\), the momentum distribution function \(n({\mathbf {q}})\), and the static density response function \(\chi ({\mathbf {q}})\). Furthermore, we are able to unambiguously characterize the impact of Fermi statistics, which is comparably small for \(S({\mathbf {q}})\) and \(\chi ({\mathbf {q}})\), but very pronounced on \(n({\mathbf {q}})\) in the small-momentum range.

Figure 1
figure 1

Average sign, S, (left) and increase in computational cost, \(C=1/S^2\), (right) for \(N=14\) and \(N=38\) \(^3{\text {He}}\) atoms as a function of the inverse temperature, \(T^{-1}\).

Full size image

We compute the imaginary time density–density correlation function \(F({\mathbf {q}},\tau )\) for the same parameters, which gives us access to the dynamic structure factor \(S({\mathbf {q}},\omega )\). First and foremost, we indeed find the familiar phonon–maxon–roton dispersion relation33 in these spectra, which is qualitatively similar to normal liquid \(^4{\text {He}}\)22 at comparable conditions. Our new PIMC data for the spectrum of collective excitations are in excellent agreement with results from neutron scattering experiments24. If they are available, experimental and theoretical results are found to be in excellent agreement.

To our knowledge, this constitutes the first comprehensive PIMC study of an ultracold atomic bulk system of fermions at finite temperature without the fixed-node approximation27, thereby opening up new avenues for the investigation of other applications such as quantum-dipole systems39,40, bilayer structures41,42, or isotopic mixtures of helium36,37,38.

Results

We consider spin-unpolarized \(^3{\text {He}}\) at a number density \(n=N/V=0.016355\) Å\(^{-3}\). Convergence with the number of imaginary-time steps has been carefully checked; see the “Methods” section for additional details. We note that the spin-polarization has negligible effects on the properties of \(^3{\text {He}}\) at these conditions; see also the discussion of Fig. 3 below.

We begin our investigation by analyzing the fermion sign problem, which constitutes the main computational bottleneck. In Fig. 1, we show our PIMC results for the average sign S (see e.g. Ref.29) for \(N=14\) and \(N=38\) atoms interacting via the Aziz-2 potential43. S constitutes a measure for the amount of cancellation of positive and negative terms in the simulations. It monotonically decreases with T. In the high-temperature limit, S approaches 1 because the effects of quantum statistics vanish, whereas \(S\rightarrow 0\) as the system approaches the ground state44. Furthermore, the error bar of an observable \({{\hat{A}}}\) scales as \(\Delta A/A\sim 1/S\), resulting in an increase in computational cost of \(C=1/S^2\)29. This trend is shown in the right panel of Fig. 1 and can be interpreted as follows: For \(T=5\;{\text {K}}\), which is close to the Fermi temperature of \(^3{\text {He}}\), the effect of quantum statistics is negligible and there is no increase in computational cost, i.e., \(C\sim 1\). In contrast, we find \(C\sim 10^3\) for \(T=2\;{\text {K}}\) and \(N=38\), which means we need 1000 times as much computer time for the fermionic calculations as we do for the bosonic case without a sign problem. While this is still feasible with \({\mathcal {O}}(10^5)\) CPU hours, it constitutes the limit of the present investigation.

In Fig. 2, we study the temperature dependence of various structural properties of normal liquid \(^3{\text {He}}\), for \(T=20, 5, 4, 3,\) and 2 K. Panel (a) shows \(S({\mathbf {q}})\) that exhibits a weak dependence on T for \(T\in [2,5]\). Here, the most pronounced temperature effect manifests in the long wavelength limit, which is determined by the isothermal compressibility22. We note that the exact \(q\rightarrow 0\) limit cannot be accessed in our simulations due to the finite simulation cell45,46. Apart from this momentum quantization effect45, we find no finite-size effects in our PIMC results; see the “Methods” section for a corresponding analysis. Both the position and the shape of the peak are hardly affected by T in this regime, which is consistent with earlier findings for \(^4{\text {He}}\)22. The solid yellow line in panel (a) shows experimental data for \(S({\mathbf {q}})\) at \(T=0.41\;{\text {K}}\) by Hallock47, which is in excellent agreement with our PIMC results for the lowest temperature. Finally, the grey triangles depict results for a substantially larger temperature of 20 K. These results are clearly distinct from the other curves, and we find temperature dependence in the small q regime and also around the peak, which decreases in magnitude and shifts to larger q.

Figure 2
figure 2

Static properties of \(^3{\text {He}}\) for different temperatures T: Shown are PIMC results for the static structure factor S(q) (a), momentum distribution function n(q) (b), and static density response function \(\chi (q)\) (c), see Eq. (1). Solid yellow: experimental data for S(q) by Hallock47.

Full size image

Panel (b) shows the momentum distribution function \(n({\mathbf {q}})\) that we have derived with the procedure described in Ref.48. For this property, the temperature plays an important role because the \(^3{\text {He}}\) atoms are propelled to larger momenta by thermal excitations. This can be seen particularly well for \(T=20\;{\text {K}}\), where n(0) is reduced by more than a factor of four compared to 2 and 3 K results. Furthermore, \(n({\mathbf {q}})\) does not resemble a step function even for the lowest depicted temperature (\(T=2\;{\text {K}}\)) for these reasons. The reduced temperature, \(\Theta =k_{{\text {B}}}T/E_{{\text {F}}}\approx 0.4\), is not quite low enough, which is exacerbated by the strong pair interaction.

Lastly, panel (c) shows the static density response function49, which we derived from the relation50,51,

$$\begin{aligned} \chi ({\mathbf {q}}) = -n\int _0^\beta {{\text {d}}}\tau \ F({\mathbf {q}},\tau ), \end{aligned}$$
(1)

with the definition of the imaginary-time correlation function

$$\begin{aligned} F({\mathbf {q}},\tau ) = \left\langle {\hat{n}}({\mathbf {q}},0) {\hat{n}}(-{\mathbf {q}},\tau )\right\rangle , \end{aligned}$$
(2)

where \({\hat{n}}({\mathbf {q}},\tau )\) is the density operator in Fourier space evaluated at \(\tau \in [0,\beta ]\); see also Ref.50 for a generalization. We find that \(\chi ({\mathbf {q}})\) exhibits an interesting, non-monotonous structure: (1) In the limits of large and small q, the response function only weakly depends on temperature for \(T\in [2,5]\text {K}\). Further, \(\chi ({\mathbf {q}})\) does not approach zero in the long-wavelength limit, as there is no perfect screening52,53 for helium because of the strong short-range interactions43. (2) The density response function exhibits a pronounced peak around \(q\approx 1.8\) Å\(^{-1}\), which corresponds to \(q\approx 2.25q_{{\text {F}}}\) (where \(q_{{\text {F}}}\) is the Fermi wave number54). This feature closely resembles recent results55,56 for the uniform electron gas at warm dense matter conditions57 at similar values of the reduced temperature \(\Theta\) and reduced wavenumber \(x=q/q_{{\text {F}}}\). Finally, the peak of the density response substantially depends on T. Specifically, the peak location is directly connected to the attractive minimum of the inter-atomic potential. Increasing the temperature reduces the correlation among the atoms, there is less of a collective behaviour, which manifests itself in a weaker density response. This trend is confirmed by \(T=20\) K results. \(\chi ({\mathbf {q}})\) is substantially reduced compared to the other temperatures, especially for \(q\lesssim 3q_{{\text {F}}}\).

In Fig. 3, we study the impact of quantum statistics on the momentum distribution of \(^3{\text {He}}\), \(n({\mathbf {q}})\). For \(T=2\;{\text {K}}\), we compare the results from our fermionic PIMC simulations (red) with the corresponding results for distinguishable particles, so-called boltzmannons. Evidently, the most pronounced differences appear for the zero-momentum state. Its occupation differs by approximately 50% because it is the momentum state with the highest occupation and therefore Pauli exclusion effects matter the most. Conversely, quantum statistical effect only change \(S({\mathbf {q}})\) and \(\chi ({\mathbf {q}})\) by little more than 1% at these conditions, which is consistent with previous findings for \(^4{\text {He}}\)1.

Figure 3
figure 3

Impact of Fermi statistics on the momentum distribution n(q). Results from simulations of fermions and boltzmannons are compared for \(T=2\;{\text {K}}\) in the upper panel and their relative deviations are shown below.

Full size image

Finally we investigate is the dynamic structure factor \(S({\mathbf {q}},\omega )\) of \(^3{\text {He}}\), which we obtain by numerically inverting the equation19,

$$\begin{aligned} F({\mathbf {q}},\tau ) = \int _{-\infty }^\infty {{\text {d}}}\omega \ S({\mathbf {q}},\omega ) e^{-\tau \omega }. \end{aligned}$$
(3)

We employ a genetic algorithm similar to the scheme presented in Ref.58, which simultaneously minimizes the \(\chi ^2\) measure of Eq. (3) as well as the first and inverse frequency moments; see the “Methods” section for more details. To our knowledge, no experimental measurements of \(S({\mathbf {q}},\omega )\) exist in the temperature range \(T\ge 2\text {K}\) that is accessible to our direct PIMC simulations. However, Sköld et al.24 presented neutron scattering data for \(T=1.2\) K. Since our investigation of \(S({\mathbf {q}})\) and \(\chi ({\mathbf {q}})\), both of which are closely related to \(S({\mathbf {q}},\omega )\), revealed no significant impact of quantum statistics, we have carried out PIMC simulations of \(N=100\) \(^3{\text {He}}\) atoms using Boltzmann statistics at this temperature.

The results of the numerical inversion of Eq. (3) are shown a heatmap in Fig. 4. All depicted datasets exhibit the phonon–maxon–roton dispersion relation that is well known for \(^4\text {He}\). The blue circles are the experimental peak positions of \(S({\mathbf {q}},\omega )\) from Ref.24. The experimental data are in excellent agreement to our results. Since we used Boltzmann statistics in our simulations, this is a strong indication that \(S({\mathbf {q}},\omega )\) is predominantly shaped by the particle interactions rather than quantum statistical effects. Furthermore our results fully corroborate previous findings for \(^3{\text {He}}\) in two dimensions33,34. The black diamonds show \(S({\mathbf {q}},\omega )\) from a theoretical investigation of \(^4{\text {He}}\) at the same T22 that differ substantially from our results for \(q\gtrsim 1.3\) Å\(^{-1}\). The differences between the two helium isotopes are mainly caused by the different particle masses because the heavier \(^4{\text {He}}\) atoms are more strongly coupled than the \(^3{\text {He}}\) species. Furthermore, our new PIMC results for \(S({\mathbf {q}},\omega )\) clearly indicate the wave-number range for future neutron scattering experiments to resolve the precise location of the roton minimum in the spectrum.

Figure 4
figure 4

(a) Heatmap of reconstructed dispersion \(S({\mathbf {q}},\omega )\) of \(^3{\text {He}}\) at \(T=1.2\;{\text {K}}\). Good agreement is found with the experimental data from Ref.24 (blue circles). Because of the higher mass, the roton minimum is shifted to higher q for \(^4{\text {He}}\) as the theoretical results from Ref.22 show (black diamonds). (b, c) Dynamic structure factor \(S({\mathbf {q}},\omega )\) for \(q=1.62\) Å−1 and \(q=2.01\) Å−1, solid purple: reconstructed PIMC data; yellow triangles: experimental data by Sköld et al.24.

Full size image

Finally we investigate the \(\omega\) dependence of \(S({\mathbf {q}},\omega )\), which is shown in Fig. 4b,c for \(q=1.62\) and 2.01 Å\(^{-1}\). A comparison of the full \(\omega\) dependence is of high importance to assess the quality of a theoretical method because the peak positions alone provide only limited information59. The yellow triangles depict the experimental measurements while the solid purple line shows the theoretical prediction from the analytic continuation. The experimental and theoretical dataset are in excellent agreement over the entire frequency range. This implies that the PIMC simulation and subsequent reconstruction are not only capable to reproduce the correct dispersion \(\omega _{{\text {max}}}(q)\) but also gives access to the actual shape of \(S({\mathbf {q}},\omega )\). From a physical perspective, we note that both considered wave numbers are located around the roton minimum (cf. Fig. 4), but do not exhibit the sharp quasi-particle excitation peak that is characteristic for the superfluid phase of \(^4{\text {He}}\)22.

To elucidate the temperature dependence of the dynamic structure factor, we include results for \(T=3\) K (dotted green) and \(T=20\) K (dashed grey) into panels (b) and (c) of Fig. 4. We find the same qualitative trend for both wave numbers, i.e., a broadening of the spectra with increasing temperature, and a shift of the maximum towards higher frequencies \(\omega\) for \(T=20\;{\text {K}}\). In addition, there appears a small positive feature around \(\omega =0\) for \(T=3~{\text {K}}\), which is particularly pronounced for \(q=1.62\) Å\(^{-1}\). We note that such a diffusive feature can be a consequence of the attractive minimum in the pair interaction60,61,62,63,64. On the other hand, the \(\omega =0\) peak could easily be an artifact due to the ill-posed nature of the reconstruction of \(S({\mathbf {q}},\omega )\) based on \(F({\mathbf {q}},\tau )\). A more detailed investigation of this effect must include the combination of different methods and will be pursued in a future publication.

Discussion

We have presented an extensive set of new PIMC results for normal liquid \(^3{\text {He}}\). Specifically, we have carried out direct PIMC calculations over the temperature range from 2 to 20 K, which are computationally challenging because of the fermion sign problem, but are numerically exact within the given error bars. We computed rigorous theoretical results for different static properties of bulk \(^3{\text {He}}\) like the static structure factor \(S({\mathbf {q}})\), momentum distribution function \(n({\mathbf {q}})\), and static density response \(\chi ({\mathbf {q}})\). We found the correlation-induced peak of \(\chi ({\mathbf {q}})\) to strongly depend on T, whereas \(S({\mathbf {q}})\) remains almost unchanged over the 2–5 K temperature interval. The absence of nodal restrictions in our simulations allowed us unambiguously quantify the impact of Fermi statistics, which is pronounced for \(n({\mathbf {q}})\) for small momenta, but practically negligible for \(S({\mathbf {q}})\) and \(\chi ({\mathbf {q}})\) at these temperature conditions.

Since quantum statistical effects change the latter quantities, we carried out PIMC simulations of \(^3{\text {He}}\) using Boltzmann statistics at a lower temperature of 1.2 K because direct fermionic PIMC simulations would have been unfeasible due to the sign problem. This has allowed us to compare the dynamic structure factor \(S({\mathbf {q}},\omega )\) to experimental measurements. The two sets of results were found to be in excellent agreement when the peak position \(\omega _{{\text {max}}}(q)\) and the actual shape of the entire spectra are compared. We have found the familiar phonon–maxon–roton dispersion relation that is well known for \(^4{\text {He}}\). This substantiates the previous findings for \(^3{\text {He}}\) in two dimensions33,34, where it has been reported that the shape of the dispersion is predominantly shaped by the interaction and not by quantum statistical effects. In addition, we analyzed the temperature dependence of \(S({\mathbf {q}},\omega )\). Primarily the spectra broaden with increasing temperature and we have found a possible diffusive feature at \(T=3\;{\text {K}}\) that will be investigated in more detail in future works.

Overall, our new results considerably extend the current understanding of one of the most important and widely studied quantum systems in the literature, which is important in its own right. Our highly accurate results for different properties constitute a benchmark for other methods such as the fixed-node approximation, and may guide the development of new, computationally less expensive approximations. In addition, we have predicted the behaviour of \(^3{\text {He}}\) at a number of temperatures, which can be checked in future experimental investigations. Possible extensions of our work might include simulations of quantum dipole systems and other types of fermionic ultracold atoms. This might provide additional insights into important physical effects like the formation of composite bosons in bilayer systems and indirect excitons42, and help to explicitly resolve the impact of quantum statistics on collective excitations39,65.

In principle, it is possible to study \(^3{\text {He}}\) based on direct PIMC simulations in the grand-canonical ensemble66, which would give access to additional physical properties such as the compressibility and the single-particle spectrum \(A({\mathbf {q}},\omega )\)65.

Finally, we note that the direct PIMC simulation of 3He along the superfluid phase transition32,67 is presently not feasible because the fermion sign problem is too severe at the required temperatures but simulations of this transition remain an import open question for future investigations.

Methods

Verification of our PIMC implementation of helium

Accurate data for bulk \(^3{\text {He}}\) at the conditions that we consider in the present work are sparse in the literature, which makes the verification of our numerical implementation with experimental data difficult. On the other hand, a gamut of theoretical works have been devoted to the bosonic simulation of \(^4{\text {He}}\), which is typically simulated with the same pair potential. In addition, direct PIMC simulations of fermions may actually be derived from a reference systems of bosons. The (exact) information about the fermionic system is then extracted by taking into account the cancellation of positive and negative contributions due to the fermionic antisymmetry; see Refs.29,68 for details. Therefore, if we are able to obtain the correct permutation structure for \(^4{\text {He}}\), we can verify the predictions of our direct PIMC simulations of fermionic \(^3{\text {He}}\).

Therefore we have simulated \(N=64\) spin-polarized \(^4{\text {He}}\) atoms at \(T=2\) K. The probability to find an individual particle involved in a permutation cycle of length l, P(l)l, is plotted in Fig. 5. The red circles show results from Boninsegni69 and the blue diamonds display the present work. First and foremost, we find excellent agreement between the two independent data sets over several orders of magnitude in the probability. In addition, we note that the curve from the present work is smoother compared to Ref.69 and includes results for longer permutation cycles, probably because we invested more computer time.

Figure 5
figure 5

Comparing permutation-cycle probabilities for \(^4\text {He}\) at \(T=2\;{\text {K}}\) to previous results by Boninsegni69.

Full size image

Finite-size effects

The direct PIMC simulations that we employ throughout this work are capable to give exact expectation values for a given combination of temperature T, volume V, and particle number N. On the other hand, we are interested in the properties of bulk \(^3{\text {He}}\) in the thermodynamic limit (\(N\rightarrow \infty\) where the number density \(n=N/V\) is being kept constant.) While the comparison of PIMC results to experimental measurements for both \(S({\mathbf {q}})\) and \(S({\mathbf {q}},\omega )\) presented in the main text is convincing, we feel that an additional analysis of finite-size effects in our PIMC data is pertinent. To this end, we show PIMC results for different properties and different numbers of particles in Fig. 6. The top left panel shows results for the static density response function \(\chi ({\mathbf {q}})\). The red circles, black stars, and green crosses have been obtained for \(N=38\), \(N=20\), and \(N=14\) \(^3{\text {He}}\) atoms, respectively. For completeness, we note that we take into account the pair interaction between atoms both within the original simulation cell, and the \(N_\mathbf{I }=3^3-1\) nearest images, which has a noticeable impact for the simulations with \(N=14\) particles despite the relatively short-range nature of the employed Aziz potential43. Still all three datasets are in remarkably good agreement. No dependence on the system size can be resolved both for small and large q, and the small deviations between the black stars and green crosses around the peak in \(\chi ({\mathbf {q}})\) do not exceed a few per cent.

Figure 6
figure 6

Finite-size effects in PIMC simulations of \(^3{\text {He}}\) at \(T=2\) K for the static density response function \(\chi ({\mathbf {q}})\) (top left), the static structure factor \(S({\mathbf {q}})\) (top right), and the momentum distribution function \(n({\mathbf {q}})\) (bottom row).

Full size image

The top right panel of Fig. 6 shows the corresponding results for the static structure factor \(S({\mathbf {q}})\). The finite-size effects are found to be even smaller. Finally, the bottom row shows results for the momentum distribution function \(n({\mathbf {q}})\) both on a linear (left) and logarithmic (right) scale. The linear scale allows us to analyse the behaviour at small momenta, whereas the logarithmic plot is well suited to resolve effects at large q. Yet, no significant deviations can be resolved between the three datasets in either panel.

Convergence with the number of imaginary-time slices

An additional practical issue in PIMC simulations of either fermions, bosons, or boltzmannons is the convergence with the number of steps in imaginary time, P. In real space, the matrix elements of the density operator, \({\hat{\rho }}=e^{-\beta {\hat{H}}}\), are given by,

$$\begin{aligned} \rho ({\mathbf {R}},{\mathbf {R}}',\beta ) = \left\langle {\mathbf {R}}|e^{-\beta {\hat{H}}}|{\mathbf {R}}'\right\rangle , \end{aligned}$$
(4)

where \({\mathbf {R}}=({\mathbf {r}}_1,\ldots ,{\mathbf {r}}_N)^T\) contains the coordinates of all N \(^3{\text {He}}\) atoms. These matrix elements are not exactly known because the kinetic (\({\hat{K}}\)) and potential (\({\hat{V}}\)) contributions to the full Hamiltonian \({\hat{H}}\) do not commute, \(e^{-\beta {\hat{H}}}\ne e^{-\beta {\hat{K}}}e^{-\beta {\hat{V}}}\). As a practical workaround, we employ the exact semi-group property of the density operator

$$\begin{aligned} e^{-\beta {\hat{H}}} = \prod _{\alpha =0}^{P-1} e^{-\epsilon {\hat{H}}}, \end{aligned}$$
(5)

where \(\epsilon =\beta /P\) is the time step of the path integral. Ultimately, Eq. (5) implies that we have to evaluate P density matrices of a temperature that is P times higher than the original one. This allows one to introduce a suitable high-temperature approximation, which becomes exact in the limit of large P.

In practice, we employ the primitive approximation

$$\begin{aligned} e^{-\epsilon {\hat{H}}} \approx e^{-\epsilon {\hat{K}}} e^{-\epsilon {\hat{V}}}, \end{aligned}$$
(6)

and the associated factorization error decays as \(P^{-2}\); see Refs.70,71 for more detailed information on this point.

As a practical example for the convergence with time step number P, we show results for the static structure factor \(S({\mathbf {q}})\) for \(N=14\) and \(T=2.5\;{\text {K}}\) in Fig. 7. This observable constitutes a representative example, as it is directly connected to the dynamic structure factor \(S({\mathbf {q}},\omega )\) [cf. Eq. (10)], to the imaginary-time correlation function \(F({\mathbf {q}},\tau )\) by the relation \(S({\mathbf {q}})=F({\mathbf {q}},0)\), and indirectly also to the static density response function \(\chi ({\mathbf {q}})\), cf. Eq. (1). Evidently, a pronounced factorization error can only be resolved for \(P=50\), while the results that were obtained with the higher P values are in very good agreement. \(P=50\) is clearly not sufficient for the description of a strongly correlated quantum liquid. We used \(P=1000\) throughout this work. Therefore, any residual factorization error are much smaller than the Monte Carlo error bars, and our results are thus well converged.

Figure 7
figure 7

Convergence with the number of imaginary-time slices P. Shown are PIMC data for the static structure factor \(S({\mathbf {q}})\) for \(N=14\) \(^3{\text {He}}\) atoms at \(T=2.5\;{\text {K}}\).

Full size image

Analytic continuation

The goal for the analytic continuation19 is to find a suitable trial solution for the dynamic structure factor \(S_{{\text {trial}}}({\mathbf {q}},\omega )\) that, when it is being inserted into Eq. (3) of the main text, reproduces the PIMC data for the imaginary-time correlation function \(F({\mathbf {q}},\tau )\) for all \(\tau \in [0,\beta ]\). Yet, these constraints are typically not sufficient to fully constrain the space of possible solutions of \(S({\mathbf {q}},\omega )\), and additional information are needed39,72. For this purpose, we consider the frequency moments of the dynamic structure factor, which are defined as

$$\begin{aligned} \left\langle \omega ^\nu \right\rangle = \int _{-\infty }^\infty {{\text {d}}}\omega \ S({\mathbf {q}},\omega ) \omega ^\nu . \end{aligned}$$
(7)

Two moments are known from the corresponding sum rules: (1) the first moment is determined by the well-known f-sum rule54

$$\begin{aligned} \left\langle \omega ^1\right\rangle = \frac{\hbar {\mathbf {q}}^2}{2m}, \end{aligned}$$
(8)

and the inverse moment is given by73

$$\begin{aligned} \left\langle \omega ^{-1}\right\rangle = -\frac{\chi ({\mathbf {q}})}{2n}. \end{aligned}$$
(9)

In addition, the zero-moment is automatically satisfied, since it holds

$$\begin{aligned} F({\mathbf {q}},0) = S({\mathbf {q}}) = \int _{-\infty }^\infty {{\text {d}}}\omega \ S({\mathbf {q}},\omega ) = \left\langle \omega ^0\right\rangle . \end{aligned}$$
(10)

The final exact property of \(S({\mathbf {q}},\omega )\) that we consider in this work is the detailed balance relation between positive and negative frequencies74,

$$\begin{aligned} S({\mathbf {q}},-\omega ) = e^{-\beta \omega }S({\mathbf {q}},\omega ), \end{aligned}$$
(11)

which is automatically fulfilled by our trial solutions \(S_{{\text {trial}}}({\mathbf {q}},\omega )\).

Many practical approaches have been suggested to accomplish goals of the analytic continuation. One family of methods is based on Bayes’ theorem, and is capable of producing smooth solutions without any unphysical sawtooth instabilities19. Yet, such maximum entropy methods21 might introduce an artificial a-priori bias into the solutions, although notable progress is continually being made in this area20,75. A second paradigm for the analytic continuation is based on the averaging over \(N_{{\text {S}}}\sim 10^3\)–\(10^4\) noisy, independent trial solutions to compose a smooth result22,58,65. While computationally more demanding, this method has the advantage that unexpected physical features of \(S({\mathbf {q}},\omega )\) might still be recovered because no explicit bias is introduced into the solution.

In the present work, we pursue the latter strategy and employ a genetic algorithm that maximizes the fitness function,

$$\begin{aligned} \Phi= & {} \left\{ \frac{1}{P/2+1}\sum _{\alpha =0}^{P/2}\frac{|F({\mathbf {q}}, \alpha \epsilon )-F_{{\text {trial}}}({\mathbf {q}},\alpha \epsilon )|}{\Delta F({\mathbf {q}},\alpha \epsilon )} + a_\chi \frac{|\chi ({\mathbf {q}})-\chi _{{\text {trial}}}({\mathbf {q}})|}{\Delta \chi ({\mathbf {q}})} + a_1 \frac{|\left\langle \omega ^1\right\rangle -\left\langle \omega ^1\right\rangle _{{\text {trial}}}|}{\Delta \left\langle \omega ^1\right\rangle } \right\} ^{-1}. \end{aligned}$$
(12)

We note that \(F({\mathbf {q}},\tau )\) is symmetric around \(\beta /2\) so that we only have to consider the interval \(\tau \in [0,\beta /2]\). Furthermore, the weights \(a_\chi\) and \(a_1\) control the respective influence of the individual constraints and are chosen empirically. We find that a reasonable choice is given by \(a_\chi = a_1 = 1/2\), which means that the frequency moments are of the same importance as \(F({\mathbf {q}},\tau )\). The f-sum rule, Eq. (8), is actually known exactly, and we empirically set \(\Delta \left\langle \omega ^1\right\rangle =\left\langle \omega ^1\right\rangle \times 10^{-3}\).

A practical example for the analytical continuation is shown in the left panel of Fig. 8 for the case considered in Fig. 4b of the main text. The red (blue) curve has been composed by averaging over \(N_{{\text {s}}}\sim 10^2\) (\(N_{{\text {s}}}\sim 10^4\)) individual noisy trial solutions \(S_{{{\text {trial}}},i}({\mathbf {q}},\omega )\). A frequency grid with \(\delta \omega \approx 0.6\;{\text {K}}\) and \(N_\omega =400\) points was employed.

Figure 8
figure 8

Left: Averaging of the dynamic structure factor \(S({\mathbf {q}},\omega )\) for \(T=1.2\;{\text {K}}\) and \(q=1.62\) Å\(^{-1}\). Red (blue) curve: average of \(N_{{\text {s}}}\sim 10^2\) (\(N_{{\text {s}}}\sim 10^4\)) noisy individual solutions. Right: Comparing PIMC results for \(F({\mathbf {q}},\tau )\) (black) and the reconstructed solution (red) for the same case as in Fig. 8. Note that \(F({\mathbf {q}},\tau )\) is symmetric with respect to \(\beta /2\), i.e., \(F({\mathbf {q}},\tau )=F({\mathbf {q}},\beta -\tau )\).

Full size image

To conclude this discussion, we consider the main input of the analytic continuation method: the \(\tau\)-dependence of the imaginary-time correlation function \(F({\mathbf {q}},\tau )\) for a specific wave number. This is shown for the present example in the right panel of Fig. 8, with the black stars and red line corresponding to the original PIMC results and the reconstructed solution, respectively. Evidently, the latter curve perfectly agrees with the PIMC data within the given Monte Carlo error bars that are on the order of \(\Delta F\sim 10^{-3}\).

All \(S({\mathbf {q}},\omega )\) results in this work have been obtained with the procedure that we just described but a note of caution is warranted. Despite being defined as an average value, it is not straightforward to estimate the true uncertainty in the final result for the reconstructed solution for \(S({\mathbf {q}},\omega )\). Rather, the space of possible trial solutions remains a property of a specific optimization method (genetic algorithm, simulated annealing, etc.) and might, or might not, include the true, physical curve with finite probability. Therefore, we abstain from using the variance of the trial solutions as an error bar, as they might underestimate the true uncertainties.