Further Computations of Quantum Fluid Triplet Structures at Equilibrium in the Diffraction Regime

1. Introduction

The equilibrium description of the quantum monatomic fluid structures at the pair level $\left(n=2\right),$ in the diffraction and Bose-Einstein regimes, can be achieved using Feynman’s path integrals (PI) [1,2] combined with computer simulation methods, i.e., path integral Monte Carlo (PIMC) and path integral molecular dynamics (PIMD) [3,4,5,6,7,8,9,10,11,12,13,14]. (The Fermi-Dirac regime is out of the conventional PI practical applications and requires a special PI formulation [15,16]; see below). The usual PI framework is directly expressed in the coordinate representation and its success at the pair level prompts the interest in undertaking the study of the quantum triplet-structural level $(n=3)$ [17,18,19]. Although fluid triplet structures in general cannot be obtained via radiation scattering experiments today [20,21,22,23,24], the equilibrium structures in statistical mechanics behave in a hierarchical manner [25,26,27,28] and the triplet step forward in the quantum domain is needed. This task implies the numerical determination of an involved variety of structural functions, i.e., $g_3(r_{12},r_{13},r_{23})$ correlation functions in real space (r-space) and $S^{\left(3\right)}(k_1,k_2,\varphi )$ structure factors in the reciprocal Fourier space (k-space) [19,29,30,31,32,33,34,35]. The different nature of these functions depends on the external field $\Psi$ applied that makes them show up.

As a matter of fact, the triplet-structural task for monatomic fluids may be regarded as already accomplished in the classical domain [36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59], where computer simulation methods (Monte Carlo (MC) and molecular dynamics (MD) [41,45,46,48,51,52,53,54,55,57]) and theories based on integral equations and closures [17,20,25,36,37,38,39,40,42,43,44,45,46,47,49,50,52,53,55,57,58] have been utilized. (Closures are cost-effective theoretical approaches that try, in general, to infer n-level structures from the knowledge of the lower-level $n^{\text{'}}<n$ structures; $n\text{'}$ and $n$ denote the elemental number of particles involved). Given that PI computer simulations can be regarded as appropriate “translations” of their classical counterparts [4,5,9], the parallel experience accumulated in the classical domain is a precious asset to tackling the quantum triplet-structural challenge. The same can be said of closures, which were used to deal with classical and quantum structures alike, regardless of their original motivations and derivations [17,20,25,36,60,61]. Nevertheless, the complexities of the quantum domain have led to consider some special features that escape the classical analogies [19].

Therefore, the quantum fluid triplet program to be followed not only shares the same general reasons that guided the corresponding classical developments, but also must include the new aspects arising from the distinct quantum behaviors. Among the general reasons, one may mention the following: statistical thermodynamics questions beyond the usual pairwise approach, the characterization of the freezing transition, the understanding of the selection between crystal lattice sites, the discussion of glass-forming liquid properties, multiple scattering phenomena, and the calculations of transport coefficients and time-dependent properties (e.g., the dynamic structure factor). Among the new aspects, which in a sense extend the scope of the latter reasons, one may mention the following quantum problems [6,9,14,15,16,17,18,19,62,63,64]: the variety of the distinct fluid responses to external fields, the effect of quantum fluctuations on the formations of crystals and glasses, the role of phonon-phonon interactions in superfluid systems, and the fixing of spin-resolved fluid structures. In connection with these problems, from the scarce initial results on quantum fluid triplets obtained so far, one observes intriguing behaviors of order parameters on the crystallization lines of liquid para-H₂ [33] and the hard-sphere fluid [34]. In addition, triplet closures have been proven to capture more quantum traits than expected and their usefulness deserves further investigation [19,34,35].

Now, some comments on the difficulties that one encounters when tackling quantum fluid triplets are in order. In the classical and the quantum domains, the dimensionality of the triplet structural functions for a homogeneous and isotropic monatomic fluid is 4-D, with conceptual and computational reasons increasing the complications in the quantum case. The simplest thermal quantum behavior is that of the diffraction effects (i.e., interference phenomena among delocalized atoms, whose magnitude cannot be disregarded); this behavior ignores any possible spin feature of the indistinguishable atoms (or the model one-site particles) composing the fluid. Such diffraction behavior can be observed in every system subjected to low-temperature conditions, but the applied temperatures are to be sufficiently high as to make the spin features negligible [1,2,4,9,12,13,14]. The spin features become fundamental at very low temperatures and add further intricacies to the fluid descriptions: (a) for integer spin atoms, one faces Bose-Einstein exchange statistics (BE) [2,9,10,11,65]; and (b) for half-odd-integer spin atoms, one faces Fermi-Dirac exchange statistics (FD) [2,15,16,65]. Typical examples of monatomic systems that can show these exchange regimes are liquid ⁴He (zero-spin atoms, BE) and liquid ³He (one-half spin atoms, FD), for which diffraction effects dominate their behaviors so long as the temperatures are $T\gtrsim 2.17 \mathrm{K}$ and $T>1\mathrm{K},$ respectively [65]. Below these temperature limits the corresponding BE and FD behaviors cannot be neglected (³He even enters the BE regime for $T\lesssim 2.7 \mathrm{m}\mathrm{K})$[65].

Interestingly, and focusing on the structural questions, the diffraction and zero-spin BE regimes admit a common general framework in that, by paying attention to their distinct peculiarities, they can be dealt with by using PI in its conventional original form [9,10,11,12,13,14,19]. However, nonzero-spin cases require special PI developments; for BE statistics see Reference [66], whereas for FD statistics see the recent works in References [15,16] that use Wigner’s formulation of quantum mechanics [67,68,69] combined with PIMC simulations (WPIMC). (Note that, when studying FD conditions with PI, any proposed method must cope with the so-called “sign problem” [15,16,70,71,72,73]). Therefore, with the addition of WPIMC, the basic methods for computing monatomic fluid structures via PI have been put forward and are, in principle, applicable to the foregoing quantum regimes. However, although the problem of pair-level structures involves affordable PI computations, the current situation for triplet structures is not so nice, even in the thermal quantum diffraction regime. This work deals with this latter regime and the presentation concentrates on it hereafter.

The thermal quantum diffraction regime may be visualized through the well-known PI image of necklaces composed of beads for representing the actual quantum atoms/particles. Despite the many PI mathematical intricacies, this allows one to study the quantum fluid and its interactions with external fields in a very intuitive way [4,9,12,13,14,74,75,76]. Thus, for example, an actual fluid composed of N atoms is represented, in the end, by a PI model consisting of N necklaces with P beads apiece ($P>1$ is an integer that becomes greater with increasing quantum effects and is to be optimized); thus, there is a one-to-one correspondence between thermally quantum delocalized atoms and necklaces. Through various applications of linear response theory [27,28,77], adapted to the treatment of the quantum domain [14,19,32,78,79], three general classes of fluid structures can be identified and defined in terms of inter-bead/inter-necklace distances in real space. One finds the following classes: (a) centroid ${\left\{g_{CMn},S_{CM}^{\left(n\right)}\right\}}_{n=\mathrm{2,3},\dots },$abbreviated to CMn (a centroid is the “center of mass” of its necklace); (b) instantaneous ${\left\{g_{ETn},S_{ET}^{\left(n\right)}\right\}}_{n=\mathrm{2,3},\dots },$ abbreviated to ETn; and (c) total thermalized-continuous linear response ${\left\{G_{TLRn},S_{TLR}^{\left(n\right)}\right\}}_{n=\mathrm{2,3},\dots },$ abbreviated to TLRn. As seen, the basic pattern per n-level is the same for every fluid class: one correlation structure $g_n/G_n$ associated with real space plus one structure factor $S^{\left(n\right)}$ associated with its corresponding reciprocal Fourier space. Each of the foregoing classes is related to the fluid response to a type of external weak field $\Psi$; within a class, the connections among its functions are far from straightforward, since their formal complexity increases with the order n [19,32]. A warning is in order here: for structural purposes, depending on the PI-scheme selected to carry out the fluid study, there may be beads that are not significant and must be disregarded; for generality, the number of structurally significant beads per necklace will be denoted by X hereafter (in normal practice, one finds either $X=P$ or $X=P/2)$[9,12,14,76]. The centroid and instantaneous quantities are formulated in ways similar to that of the familiar classical domain [26,27,28]; however, in the total thermalized-continuous linear response case, one faces more entangled formulations [14,19,32].

Furthermore, it is worth remarking that radiation scattering experiments allow one to obtain the pair-level ET2 and TLR2 structures [21,22,23,24]; these structures are defined by the actual atoms composing the sample, because atoms interact effectively with the external field. However, given that the (small) intensity of the triplet contribution is hidden in the whole intensity of the outgoing signal [20], no structural function can be experimentally determined beyond $n=2.$ As for the centroid class CMn, although not even the pair CM2 structures can be directly obtained in any experiments (centroids cannot actually couple with any fields), its importance as an intermediate theoretical object cannot be overemphasized, as its usefulness is astonishing. In this regard, centroid structures do follow the same rules as those of a classical fluid [19,30,78], and their applications cover from the exact connections with the fluid equation of state to order parameters for charactering quantum samples to appealing approximate formulations of quantum dynamics, and even to certain BE exchange questions [6,7,12,14,19,30,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95].

The foregoing comments on the distinct quantum structural classes point to the greater complexity of the quantum studies as compared to the more reduced scope of their classical counterpart, where only one class of structures ${\left\{g_{Cn},S_C^{\left(n\right)}\right\}}_{n=\mathrm{2,3},\dots },$ abbreviated to Cn, exists. (Note that Cn can be extended further to incorporate, in a formally exact way, the n-body direct correlation functions $c_{Cn}$ [37,38,40,45]). Nevertheless, these theoretical reasons are just half the issue. The other half includes the practical causes behind the current scarcity of quantum triplet structure applications using “exact” path integral methodologies, and it is considered below.

PIMC computational studies on triplet structures in homogeneous and isotropic quantum monatomic fluids [19,29,30,31,32,33,34,35] have focused mostly on the equilateral and isosceles features of the instantaneous ET3 and the centroid CM3 structures. PI simulation work needs powerful computational resources and time-consuming runs; the results so obtained are “exact”, that is, with controllable statistical errors. Therefore, the desired accuracy of the numerical answers can be increased, but at the expense of more calculations. Such a fact has its ramifications. In this connection, one must consider the PI sample size, which leads to the analysis of the expanded configurational space resulting from the introduction of the necklaces and their associated beads: for a system of N actual quantum atoms, the dimensionality of the configuration space is $3N,$ whereas for its PI model the corresponding dimensionality is $3NP.$ Moreover, the number of different quantum structural classes mentioned above magnify the task. Roughly speaking, one finds in PI work that: (a) centroid structure calculations scale with N, as in their classical counterparts [41,48]; and (b) the instantaneous and the total thermalized-continuous linear response calculations, apart from the usual N-dependence, also scale with X and $X^n,$ respectively [19,29,32]. Furthermore, aside from the 4-D character of the triplet functions, structure factors in simulation work are fixed via the scanning of appropriate sets of wave vectors, because the latter must: be commensurate with the basic box in which the sample is contained [46], and target independently the three different classes of structures [33]. Therefore, the quantum diffraction situation contrasts sharply with that of the classical monatomic fluid composed of N atoms, where the configuration space has a dimensionality $3N$ and only one class of structures is to be targeted. Owing to these difficulties, in the quantum treatments of the classes CM3 and ET3, the introduction of closures and direct correlation functions $c_3$ [17,25,38,45,47] can be expected to play a significant role [19,29,30,31,32,33,34,35]. However, although these mathematical tools augment the reach of the conceptual framework and are intended to diminish the computational overload, their quantum applications are not free from uncertainties [19,45], as happens (to a lower extent) in the classical domain.

The aim of this work is to continue the study of fluid triplet structures in the thermal quantum diffraction regime. Two systems that were analyzed by this author in pilot investigations are selected: helium-3 under supercritical conditions [19,35,96], and the quantum hard-sphere fluid (QHS) on its crystallization line [19,34,97]. PIMC simulations are performed for both systems in the canonical ensemble. The PI propagators selected are: (a) for helium-3, the fourth-order propagator in Voth et al’s form [12] (SCVJ, based on the Suzuki-Chin’s developments [98,99]), its application involving the use of Janzen-Aziz’s SAPT2 interatomic potential [100]; and (b) for QHS, Cao-Berne’s pair action [76]. Also, complementary closure calculations are reported; the following closures are utilized: Jackson-Feenberg convolution [17], Kirkwood superposition [25], the intermediate closure AV3 [34], and Denton-Ashcroft approximation in symmetrized form [47]. Obviously, exhaustiveness in this study is not attempted since, for the time being, the computational load involved would be immense. Rather, the current work adds up to the knowledge of these two systems by exploring quantum triplet properties under very different fluid conditions: far and near the first-order fluid-solid transition. The results focus on salient triplet features (equilateral and isosceles) of the centroid and instantaneous structures, by covering spatial correlations in helium-3, and structure factors in helium-3 and QHS. Path integral and closure results are compared, and a number of further procedural and structural conclusions are drawn.

The outline of this article is as follows. Section 2 contains the basic PI theory. Section 3 describes the PI triplet centroid and instantaneous structural concepts. Section 4 focuses on the triplet closures utilized in this work. Section 5 gives the computational details, and Section 6 the results and their discussion. Finally, Section 7 collates a number of closing remarks and perspectives for future work.

2. Path Integral Background

In the study of quantum monatomic fluids, one utilizes the usual statistical ensembles [1,2,26,27,28]. For constraints defining a closed system there are the canonical $(N,\mathcal{V},T)$ and the isothermal-isobaric $(N,p,T)$ ensembles, while for constraints defining an open system there is the grand canonical ensemble $(\mu ,\mathcal{V},T)$. The variables, whose values are held fixed in each case, are the number of particles N, the volume $\mathcal{V}$, the temperature T, the pressure $p$, and the chemical potential $\mu$. Within the conventional statistical mechanical treatments, the atoms composing the system under study are considered structureless particles, their spatial positions in the system being defined by those of their nuclei ${\left\{{\mathbf{r}}_j\right\}}_{j=\mathrm{1,2},\dots ,N,\dots }.$ The latter reduction is consistent with the experimental techniques that reveal the actual existence of fluid structures [20,21,22,23,24,27], and it implies that the electronic degrees of freedom have undergone a quantum mechanical averaging process for setting an adequate potential energy interaction function $V^{\left(N\right)}\left({\mathbf{r}}_1,{\mathbf{r}}_2,\dots ,{\mathbf{r}}_N\right)$ (non-collapsing and tempered [28]). Such an operation is rooted in the Born-Oppenheimer approximation [101]. In the absence of any external field $\Psi ,$ the whole situation for a quantum system is contained essentially in the canonical density matrix given by the operator $\exp \left(-\beta H_0^{\left(N\right)}\right)$ [2,4,9,28]. $H_0^{\left(N\right)}$ stands for the isolated system Hamiltonian, which is built as $H_0^{\left(N\right)}=T^{\left(N\right)}+V^{\left(N\right)},$ where one includes the operators for the kinetic energy, $T^{\left(N\right)}=-\left({\hslash }^2/2m\right){\sum }_j{\nabla }_j^2,$ and the potential energy, $V^{\left(N\right)}\left({\mathbf{r}}_1,{\mathbf{r}}_2,\dots ,{\mathbf{r}}_N\right).$ Once the form of $H_0^{\left(N\right)}$ is fixed, the canonical partition function $Z^Q$ arises as the trace of the density matrix, $Z^Q(N,\mathcal{V},T)=Tr\left[\exp \left(-\beta H_0^{\left(N\right)}\right)\right]$ [2,26,28], which allows one to set the basic thermodynamic connection $A=-k_{B}Tln{Z}^Q,$ where A stands for Helmholtz free energy. The definitions of the partition functions for the isothermal-isobaric or the grand canonical ensembles are straightforward using $Z^Q$ [26,27]. Within this context, recall that there is a Wick rotation contained in the formal equivalence $\mathrm{e}\mathrm{x}\mathrm{p}(-\beta H)\equiv \mathrm{e}\mathrm{x}\mathrm{p}(-iH\vartheta /\hslash$), where $H$ stands for a time-independent Hamiltonian and $\vartheta$ for real time. The definition of an imaginary time as $\tau =\beta \hslash$ follows from that equivalence [1,2], and $\tau$ emerges as a thermal quantum variable for the characterization of equilibrium states. Furthermore, to give operational forms to the partition functions, one uses the PI approach and obtains their proper adaptations [4,8,9,11,12,76,102,103]. To introduce the basic concepts, attention will be given to the canonical ensemble questions in this Section.

2.1. PI Canonical Partition Function

When using the PI formalism [2], every atom/nucleus j is represented by an elastic path ${\mathbf{r}}_j\left(\tau \right)$ in imaginary time $0\le \tau \le \beta \hslash$. The general form of the path-integral canonical partition function of a monatomic quantum fluid, e.g., composed of helium atoms, with all the atoms in the same spin state, can be written in the canonical ensemble as [2]:

$$\begin{gathered} Z^Q\left(N,\mathcal{V},T\right)=Tr\left[\mathrm{e}\mathrm{x}\mathrm{p}(-\beta H_0^{\left(N\right)})\right]={\displaystyle \frac{1}{N!}}\sum _{\mathcal{P}}{\pi }_{\mathcal{P}}\int {\oint }_{\begin{array}{c}{\mathbf{r}}_j\left(0\right)={\mathbf{r}}_j\\ {\mathbf{r}}_j(\beta \hslash )={\mathcal{P}\mathbf{r}}_j\end{array}}\prod _{j=1}^N\mathcal{D}{\mathbf{r}}_j\left(\tau \right)d{\mathbf{r}}_j\left(0\right)\times \\ \mathrm{e}\mathrm{x}\mathrm{p}\left[-{\displaystyle \frac{1}{\hslash }}{\int }_0^{\beta \hslash }\left({\sum }_{j=1}^N{\displaystyle \frac{m{\left(\dot{{\mathbf{r}}_j}\right)}^2}{2}}+V^{\left(N\right)}\left({\mathbf{r}}_1\left(\tau \right),{\mathbf{r}}_2\left(\tau \right),\dots ,{\mathbf{r}}_N\left(\tau \right)\right)\right)d\tau \right], \end{gathered}$$

(1)

The symbol $\mathcal{P}$ runs over the $N!$ permutations among the $N$ indistinguishable atoms, ${\pi }_{\mathcal{P}}$ is the sign of the permutation $\mathcal{P}$ (i.e., +1 for every boson permutation; $\pm 1$ for fermion permutations that, depending on their parity, is $+1$ for even $\mathcal{P}$ or $-1$ for odd $\mathcal{P}$). In Equation (1a), the integrations cover all the configuration space associated with the particle paths ${\left\{{\mathbf{r}}_j\left(\tau \right)\right\}}_{j=\mathrm{1,2},\dots ,N}$, the symmetry constraints are denoted by the conditions at $\tau =0$ and $\tau =\beta \hslash ,$ which imply that the particle paths may interlink for $\mathcal{P}\ne$Identity in a large variety of ways. The exp-factor contains the action in imaginary time, where: $\dot{{\mathbf{r}}_j }$stands for the derivative with respect to $\tau ,$ and $V^{\left(N\right)}$ is the potential energy function acting at equal-$\tau$ “instants”. Equivalently, $Z^Q$ can also be cast in the coordinate representation as [2,4,9]:

$$Z^Q\left(N,\mathcal{V},T\right)=Tr\left\{\mathrm{e}\mathrm{x}\mathrm{p}(-\beta H_0^{\left(N\right)})\right\}={\displaystyle \frac{1}{N!}}\sum _{\mathcal{P}}{\pi }_{\mathcal{P}}\int d{\mathbf{r}}^N⟨{\mathbf{r}}^N|\mathrm{e}\mathrm{x}\mathrm{p}(-\beta H_0^{\left(N\right)})|{\mathcal{P}\mathbf{r}}^N⟩,$$

(1b)

where $d{\mathbf{r}}^N=d{\mathbf{r}}_{1}d{\mathbf{r}}_2\dots d{\mathbf{r}}_N$ is the $3N-$dimensional volume element of the configurational space of the $N$actual atoms, $\left|{\mathbf{r}}^N⟩\right.$ is the ket of position states $\left|{\mathbf{r}}_1{,\mathbf{r}}_2,\dots ,{\mathbf{r}}_N⟩\right.$, and the permutations act on the particle position states as ${|\mathcal{P}\mathbf{r}}^N⟩=$ $\left|\mathcal{P}{\mathbf{r}}_1{,\mathcal{P}\mathbf{r}}_2,\dots ,{\mathcal{P}\mathbf{r}}_N⟩\right..$ At sufficiently high temperature, only the identity permutation contributes to effectively shape Equations (1), the BE or FD quantum statistics features can be neglected, and both exchange situations lead to the same general result in the form of “distinguishable” delocalized atoms (this is nicely illustrated by helium systems [9,96]). Such a situation is that of thermal quantum diffraction/dispersion effects, for which the PI canonical partition function reads as [2,4]:

$$\begin{gathered} Z^D\left(N,\mathcal{V},T\right)=Tr{\left\{\mathrm{e}\mathrm{x}\mathrm{p}(-\beta H_0^{\left(N\right)})\right\}}_D={\displaystyle \frac{1}{N!}}\int {\oint }_{{\mathbf{r}}_j\left(0\right)={\mathbf{r}}_j(\beta \hslash )}\prod _{j=1}^N\mathcal{D}{\mathbf{r}}_j\left(\tau \right)d{\mathbf{r}}_j\left(0\right)\times \\ \mathrm{e}\mathrm{x}\mathrm{p}\left[-{\displaystyle \frac{1}{\hslash }}{\int }_0^{\beta \hslash }\left({\sum }_{j=1}^N{\displaystyle \frac{m{\left(\dot{{\mathbf{r}}_j}\right)}^2}{2}}+V^{\left(N\right)}\left({\mathbf{r}}_1\left(\tau \right),{\mathbf{r}}_2\left(\tau \right),\dots ,{\mathbf{r}}_N\left(\tau \right)\right)\right)d\tau \right], \end{gathered}$$

(2a)

where every atom j is represented by an elastic closed path ${\mathbf{r}}_j\left(\tau \right)$ in imaginary time, such that ${\mathbf{r}}_j\left(0\right)={\mathbf{r}}_j\left(\beta \hslash \right).$ In the coordinate representation, Equation (2a) reads as follows [4,9]:

$$Z^D\left(N,\mathcal{V},T\right)=Tr{\left\{\mathrm{e}\mathrm{x}\mathrm{p}(-\beta H_0^{\left(N\right)})\right\}}_D={\displaystyle \frac{1}{N!}}\int d{\mathbf{r}}^N⟨{\mathbf{r}}^N|\mathrm{e}\mathrm{x}\mathrm{p}(-\beta H_0^{\left(N\right)})|{\mathbf{r}}^N⟩.$$

(2b)

In Equations (2), although only the identity permutation is retained and the atoms can be “numbered”, the presence of $N!$ as the remaining factor of the true indistinguishability among the atoms is to be noticed; the presence of $\hslash ,$ as a factor necessary to guarantee the correct dimensionality, is also to be noticed. (Both $N!$ and $\hslash$ are indispensable for the formulation of the classical partition function [2]). In general, for many-body system studies, Equations (1b) and (2b) are the practical forms of the partition function that serve as starting points to develop the PI-discretized numerical approaches [4,9,12,76].

In this work, attention is focused on Equation (2b) and its connections with the incorporation of weak external fields ${\Psi }^{\left(N\right)}$ acting on the atoms/nuclei/particles of the quantum fluid. The fields define the three classes of equilibrium structures $\left\{g_n,S^{\left(n\right)}\right\}$ in a fluid under quantum diffraction conditions, with linear response theory being instrumental in the corresponding derivations. In this context, the specific nature of the field and its interactions with the atoms/nuclei/particles are determinant, which contrasts with its classical counterpart [19,27,32,77,78,79].

In this connection, some cautionary remarks are to be made. The discussion in this Section considers the path-integral case involving weak continuous fields, which is a situation amenable to classical-like treatments and allows one to introduce most of the basic concepts. It so happens that the total thermalized-continuous linear response (TLRn) and centroid (CMn) classes are directly linked to the following developments [14,19,32]. However, the instantaneous class (ETn) escapes such a direct treatment based on continuous fields, owing to the localizing fields involved in its definition. Despite the previous remark, Equation (2b) serves equally well the purposes of formulating the instantaneous structural functions at any n-level, because the underlying arguments involved are also associated with linear response theory [9,19,27]. The structural facts relevant to this work (CMn and ETn) are considered in the next Section.

2.2. The Action of Weak External Continuous Fields

In the presence of a weak external continuous field $\Psi ,$ acting as ${\Psi }^{\left(N\right)}={\sum }_{j=1}^N\Psi \left({\mathbf{r}}_j\right),$ the canonical partition function for diffraction effects can be cast in the coordinate representation as [14]:

$$\begin{gathered} Z^D\left(N,\mathcal{V},T;\Psi \right)=Tr{\left\{\mathrm{e}\mathrm{x}\mathrm{p}(-\beta {(H}_0^{\left(N\right)}+{\Psi }^{\left(N\right)}))\right\}}_D= \\ {\displaystyle \frac{1}{N!}}\int d{\mathbf{r}}^N⟨{\mathbf{r}}^N|\mathrm{e}\mathrm{x}\mathrm{p}(-\beta {(H}_0^{\left(N\right)}+{\Psi }^{\left(N\right)}))|{\mathbf{r}}^N⟩. \end{gathered}$$

(3)

The operators $H_0^{\left(N\right)}$ and ${H^{\left(N\right)}=H}_0^{\left(N\right)}+{\Psi }^{\left(N\right)}$ are assumed to be self-adjoint, having complete sets of eigenvectors and eigenvalues. (To generalize the right-hand side of Equation (2a) by including $\Psi ,$ the action integral should contain the term ${\Psi }^{\left(N\right)}\left({\mathbf{r}}_1\left(\tau \right),{\mathbf{r}}_2\left(\tau \right),\dots ,{\mathbf{r}}_N\left(\tau \right)\right))$. To obtain general operative forms, one applies to the density matrix operator the identity [9]:

$$\exp \left(-\beta H^{\left(N\right)}\right)=\stackrel{X factors}{\overbrace{ \exp \left(-\beta H^{\left(N\right)}/X\right)\exp \left(-\beta H^{\left(N\right)}/X\right)\dots \exp \left(-\beta H^{\left(N\right)}/X\right) }}.$$

(4)

Next, one inserts in Equation (4) $(X-1)$ spectral resolutions of the identity in the coordinate representation, i.e., $\int d{\mathbf{r}}^{N,t^{*}}|{\mathbf{r}}^{N,t^{*}}⟩⟨{\mathbf{r}}^{N,t^{*}}|=1$, where $t^{*}=\mathrm{1,2},\dots ,X,$ $|{\mathbf{r}}^{N,t^{*}}⟩=|{\mathbf{r}}_1^{t^{*}},{\mathbf{r}}_2^{t^{*}},...,{\mathbf{r}}_N^{t^{*}}⟩$ and $d{\mathbf{r}}^{N,t^{*}}={d\mathbf{r}}_1^{t^{*}}d{\mathbf{r}}_2^{t^{*}}...d{\mathbf{r}}_N^{t^{*}}.$ Then, by denoting ${\mathbf{r}}^{N,1}\equiv {\mathbf{r}}^N,$ Equation (4) reads as:

$$\begin{gathered} Z^D\left(N,\mathcal{V},T;\Psi \right)={\displaystyle \frac{1}{N!}}\int d{\mathbf{r}}^{N,1}d{\mathbf{r}}^{N,2}\dots d{\mathbf{r}}^{N,X}⟨{\mathbf{r}}^{N,1}|\mathrm{e}\mathrm{x}\mathrm{p}(-\beta H^{\left(N\right)}/X)|{\mathbf{r}}^{N,2}⟩\times \\ ⟨{\mathbf{r}}^{N,2}|\mathrm{e}\mathrm{x}\mathrm{p}(-\beta H^{\left(N\right)}/X)|{\mathbf{r}}^{N,3}⟩⟨{\mathbf{r}}^{N,3}|\mathrm{e}\mathrm{x}\mathrm{p}(-\beta H^{\left(N\right)}/X)|{\mathbf{r}}^{N,4}⟩\dots ⟨{\mathbf{r}}^{N,X}|\exp \left(-\beta H^{\left(N\right)}/X\right)|{\mathbf{r}}^{N,1}⟩. \end{gathered}$$

(5)

Now, Equation (5), which is an exact expansion, needs explicit forms for the propagator $⟨{\mathbf{r}}^{N,t^{*}}|\exp \left(-{\displaystyle \frac{\beta H^{\left(N\right)}}{X}}\right)|{\mathbf{r}}^{N,t^{*}+1}⟩$, where the cyclic property $t^{*}+1=X+1\equiv 1$ is to be applied. In addition, for every atom j, the latter property brings forth discretized closed trajectories: ${\mathbf{r}}_j^{t^{*}=1}\to {\mathbf{r}}_j^{t^{*}=2}\to \dots {\mathbf{r}}_j^{t^{*}=X}\to {\mathbf{r}}_j^{t^{*}=1}$. In the limit $X\to \infty ,$ any ${\mathbf{r}}_j(t^{*}$) trajectory transforms into a closed continuous path ${\mathbf{r}}_j\left(\tau \right)$ in imaginary time [2,9]. Therefore, any ${\mathbf{r}}_j(t^{*}$) trajectory takes place in real space and evolves in imaginary time according to the “time”-step $\beta \hslash /X$ (see further details below). The foregoing general propagator, defined in the coordinate representation, takes nonnegative values [9], a decisive fact that guarantees the existence of probability densities in the related PI formulations. Note that the selection of the size of the imaginary time-step will influence the propagator form. Also, one must bear in mind that the kinetic energy operator $T^{\left(N\right)}$ does not commute with the potential energy operators $V^{\left(N\right)}$ and ${\Psi }^{\left(N\right)};$ hence, $H_0^{\left(N\right)}$ does not commute with ${\Psi }^{\left(N\right)}$ either.

Given that the PI partition function in the associated configurational space is consistent with a proper statistical analysis, the set of PI elastic closed paths is equivalent to a statistical distribution of closed paths. Factors that affect such statistical distribution are: (a) the symmetry among the N descriptions ${\left\{{\mathbf{r}}_j\left(t^{*}\right)\right\}}_{j=\mathrm{1,2},\dots ,N}$ of the particles, and (b) the packing of the particles and/or other parameters relevant for the system description [4,5,6,7,8,9,12,13,14]. In relation to this, the actual bulk density of a system is taken as the number density ${\rho }_N$ (e.g., for N particles in a fixed volume $\mathcal{V},$ the density is ${\rho }_N={\displaystyle \frac{N}{\mathcal{V}}}$).

The formulation of structures based on developments of Equation (5) is a delicate matter, and one should deal first with the questions involving $H_0^{\left(N\right)}$ and ${\Psi }^{\left(N\right)}$. One can proceed by applying the straightforward operator relationship [14]:

$$\exp \left(-{\displaystyle \frac{\beta H^{\left(N\right)}}{X}}\right)=\exp \left(-{\displaystyle \frac{\beta }{X}}\left[H_0^{\left(N\right)}+{\Psi }^{\left(N\right)}\right]\right)=\exp \left(-{\displaystyle \frac{\beta {\Psi }^{\left(N\right)}}{2X}}\right)\exp \left(-{\displaystyle \frac{\beta H_0^{\left(N\right)}}{X}}\right)\exp \left(-{\displaystyle \frac{\beta {\Psi }^{\left(N\right)}}{2X}}\right)+O\left(X^{-3}\right),$$

(6)

which arises from the Baker-Campbell-Haussdorf formula for the exponential operator built as the sum of two non-commuting operators [104]. Potential energy operators are diagonal in the coordinate representation, i.e., ${\Psi }^{\left(N\right)}|{\mathbf{r}}^{N,t^{*}+1}⟩=|{\mathbf{r}}^{N,t^{*}+1}⟩{\Psi }^{\left(N\right)}\left({\mathbf{r}}^{N,t^{*}+1}\right)$, a fact that leads to the PI-discretized approximation $Z_{PI}^D$ for the canonical partition function [14]:

$$\begin{gathered} Z^D\left(N,\mathcal{V},T;\Psi \right)\approx Z_{PI}^D\left(N,\mathcal{V},T;\Psi \right)= \\ {\displaystyle \frac{1}{N!}}\int {\prod }_{t^{*}=1}^{X^{\text{'}}}d{\mathbf{r}}^{N,t*}⟨{\mathbf{r}}^{N,t^{*}}|\mathrm{e}\mathrm{x}\mathrm{p}(-\beta H_0^{\left(N\right)}/X)|{\mathbf{r}}^{N,t^{*}+1}⟩\times \exp \left(-{\displaystyle \frac{\beta }{X}}{\sum }_{j=1}^N{\sum }_{t^{*}=1}^X\Psi \left({\mathbf{r}}_j^{t^{*}}\right)\right), \end{gathered}$$

(7)

where $X^{\text{'}}$ implies the abovementioned cyclic property $t^{*}+1=X+1\equiv 1$ for the $t^{*}-$product. The foregoing partition function is accurate up to terms $O\left(X^{-2}\right)$. In the limit $X\to \infty$, one retrieves the actual $Z^D(\Psi ),$ since the operators $H_0^{\left(N\right)}$ and ${\Psi }^{\left(N\right)}$ are self-adjoint and make sense separately, as required by Trotter’s analysis [3,9] (the same applies to $T^{\left(N\right)}$ and $V^{\left(N\right)}$ regarding questions that involve the isolated $H_0^{\left(N\right)}$ operator [9]). Given that the propagators contained in Equation (7) must be nonnegative everywhere in the configuration space corresponding to the set of variables $\left\{{\mathbf{r}}^{N,t^{*}}\right\}$, it is easy to identify the underlying probability density of this formulation. Therefore, classical-like statistical developments and calculations involving $Z_{PI}^D$ can be performed, and $X$ represents a compromise between statistical convergence and theoretical accuracy (i.e., $X\to \infty )$ [3], which means that $X$ can be properly optimized for the problem under study. Remarkably, linear response theory utilizes the structural functions of the isolated-from-$\Psi$ fluid, that is, those functions derived from the consideration of $H^{\left(N\right)}=H_0^{\left(N\right)}$ alone. Hence, questions on the accuracy of the structures obtained through the use of Equations (6) and (7), related to the number of “instants” $X$, are merely formal, because one can optimize $X$ by analyzing just the isolated-from-$\Psi$ fluid structures.

2.3. Propagators for the PI Discretized Canonical Partition Function

At this point, it is convenient to discuss the options for the propagator $⟨{\mathbf{r}}^{N,t^{*}}|\mathrm{e}\mathrm{x}\mathrm{p}(-\beta H_0^{\left(N\right)}/X)|{\mathbf{r}}^{N,t^{*}+1}⟩$ in the absence of $\Psi ,$ which will give final forms to $Z_{PI}^D$. For structural studies, there are three main types of propagators: primitive [4,5,6,7,8,9], pair actions [9,76], and the fourth-order SCVJ [12,98,99]. (See also another fourth-order propagator version in Reference [105]). Two main interrelated facts take part in their derivations. First, there is the form in which the imaginary-time step is chosen. Thus, if $\beta$ is divided into $X=P$ equally-spaced intervals, with the step-unit set to $\beta /P$ $(P=$ arbitrary integer $>1),$ one will deal with the primitive or a pair action propagator. However, if a double step $2\beta /P$ is selected, making $P$ be an even positive integer, i.e., $X=P/2,$ one can deal with the SCVJ propagator. In the end, for practical and consistency reasons, the final number of intervals will be P, regardless of the type of propagator selected, although the roles of the interval ends may not be equivalent in the final description achieved (e.g., the SCVJ case). Second, there are different choices to deal with the noncommutativity between $T^{\left(N\right)}$ and $V^{\left(N\right)}.$ These choices are in the roots of the picture of $N$ necklaces with $P$ beads apiece for the discretized closed paths in imaginary time, ${\mathbf{r}}_j(t^{*}$), of the particles: such $N\times P$ picture is common to the three types of propagators considered. In this connection, for a sufficiently large $P,$ the three propagators yield the same final description; thus, a critical issue in PI calculations is the distinct statistical convergence rates that, depending on the optimal P, can be achieved with each type of propagator [9,12,14,76].

Although for specific details the reader is referred to the References quoted above, for the current purposes, it is worth stressing that the general form of $Z_{PI}^D$ can be cast as [9,12,14]:

$$\begin{gathered} Z_{PI}^D\left(N,\mathcal{V},T;\Psi \right)= \\ {\displaystyle \frac{1}{N!}}{\left({\displaystyle \frac{mP}{2\pi \beta {\hslash }^2}}\right)}^{{\displaystyle \frac{3NP}{2}}}\int \prod _{t=1}^{P}d{\mathbf{r}}^{N,t}\times \mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta W_{NP}({\mathbf{r}}_1^1,{\mathbf{r}}_1^2,\dots {\mathbf{r}}_1^P,\dots ,{\mathbf{r}}_N^P;\beta ,\hslash ,m\right))\times \\ \exp \left(-{\displaystyle \frac{\beta }{X}}{\sum }_{j=1}^N{\sum }_{t^{*}=1}^X\Psi ({\mathbf{r}}_j^{t^{*}})\right), \end{gathered}$$

(8)

where $X=P$ or $P/2$ and, for notational simplicity, two symbols to denote the beads, t and $t^{*},$ are employed. The relations between the bead sets $\left\{t=\mathrm{1,2},\dots ,P\right\}$ and $\left\{t^{*}=\mathrm{1,2},\dots ,X\right\}$ are explained below (the $t^{*}$ elements coincide with choices extracted from the $t$ set) [14].

In Equation (8), the integrations cover the whole $3NP$-dimensional configurational space associated with the $N$ closed necklaces, $j=\mathrm{1,2},\dots ,N;$ each necklace is composed of $P$ beads, $t=\mathrm{1,2},\dots ,P.$ $W_{NP}$ contains the specific connections/interactions among all the beads in the sample, thereby being an “effective potential” for the resulting classical-like set of $N\times P$ beads. Thus, one finds the so-called (semi) classical isomorphism [4], in which the $W_{NP}$ form containing the bead interactions depends on the propagator selected and, also, on $\beta ,m,$ and $\hslash .$ While for the primitive and pair action propagators the situation is straightforward $\left(X=P\right),$ the SCVJ derivation doubles necessarily the initial number $X$ of instants, $t^{*}=\mathrm{1,2},\dots ,X=P/2,$ by adding $P/2$ intermediate instants to the initially closed X-trajectory [12]. Therefore, the complete SCVJ bead set is renumbered according to the total number of beads $P,$ $t=\mathrm{1,2},3,\dots ,P-1,P,$and the odd-numbered ones are associated with the field-bead interactions, as stated in Equations (7) and (8). Linear response theory [14,77] uses such interactions in its derivations and, therefore, one can anticipate that the significant beads for the definitions of structures will be as follows: (a) for the primitive and pair action propagators, $X=P$and $t^{*}=t=\mathrm{1,2},,3\dots ,P,$ with all the beads being significant and equivalent; and (b) for SCVJ, $X=P/2$ and $t=\mathrm{1,2},,3\dots ,P,$ with the odd-numbered beads playing the role of the $t^{*}$ initial ones, $t^{*}\equiv t=\mathrm{1,3},5,\dots ,P-1$, these being the significant ones and also equivalent among themselves. (Equivalent means that the mathematical treatment is the same for every bead in the sample).

It may be worthwhile to highlight some facts contained in the foregoing discussion. Each closed path contains $X$ marked imaginary-time instants $(X=P$ or $X=P/2),$ such that: conventionally, $t^{*}=1$ corresponds to ${\tau }_1=0\equiv \beta \hslash ,t^{*}=2$ to ${\tau }_2={\displaystyle \frac{\beta \hslash }{X}},$ $t^{*}=3$ to ${\tau }_3={\displaystyle \frac{2\beta \hslash }{X}},$…, and$t^{*}=X$ to ${\tau }_X={\displaystyle \frac{(X-1)\beta \hslash }{X}}$. The continuous description is retrieved in the limit $X\to \infty ,$ which may seem obvious but runs deeper than one might think at first sight, as no errors creep in the process (i.e., Trotter’s theoretical accuracy) [3,9]. Moreover, note that the taking of the imaginary-time origin at ${\tau }_1=0$ has nothing special; it could have been taken at any other ${\tau }_n$ instant of the X-sequence, since translational $\tau -$invariance must hold. In this regard, the only requirement is that once the origin is chosen, it must be applied to every atom/nucleus/particle path in the system for consistency reasons [9]. The imaginary-time evolution is periodic in that $0\equiv \beta \hslash ,$ as is customary in Fourier analysis (this does not mean that $\beta =0$ whatsoever!). Finally, the recipe for structural studies: all the beads are significant in the primitive and pair action cases, although only the odd-numbered beads are significant in the SCVJ case.

For the reader’s convenience, the explicit forms of the effective potential $W_{NP}$ dealt with in this article are worth giving. The optimal P discretization is assumed in each case. For helium-3, the SCVJ propagator leads to [12]:

$$\begin{gathered} W_{NP}^{SCVJ}={\displaystyle \frac{{(m}_{{}_3\mathrm{H}\mathrm{e}})P}{2{\beta }^2{\hslash }^2}}{\sum }_{j=1}^N{\sum }_{t=1}^{P^{\text{'}}}{\left({\mathbf{r}}_j^t{-\mathbf{r}}_j^{t+1}\right)}^2+{\displaystyle \frac{2}{3P}}\sum _{j<l}\left\{\sum _{\begin{array}{c}t-odd\\ \mathrm{1,3},5,\dots ,P-1\end{array}}u\left(r_{jl}^t\right)+2\sum _{\begin{array}{c}t-even\\ \mathrm{2,4},6,\dots ,P\end{array}}u\left(r_{jl}^t\right)\right\}+ \\ {\displaystyle \frac{{\beta }^2{\hslash }^2}{9\left(m_{{}_{3}He}\right)P^3}}{\sum }_{j=1}^N\left\{\alpha \sum _{\begin{array}{c}t-odd\\ \mathrm{1,3},5,\dots ,P-1\end{array}}{\left(\sum _{l\ne j}{\displaystyle \frac{\partial u\left(r_{jl}^t\right)}{\partial r_{jl}^t}}{\mathsf{\eta }}_{jl}^t\right)}^2+(1-\alpha )\sum _{\begin{array}{c}t-even\\ \mathrm{2,4},6,\dots ,P\end{array}}{\left(\sum _{l\ne j}{\displaystyle \frac{\partial u\left(r_{jl}^t\right)}{\partial r_{jl}^t}}{\mathsf{\eta }}_{jl}^t\right)}^2\right\}, \end{gathered}$$

(9)

where $m=m_{{}_3\mathrm{H}\mathrm{e}}$ is the atom mass, use is made of the pairwise interaction approach $V^{\left(N\right)}={\sum }_{j<l}u\left(r_{jl}\right)$ involving the pair potentials $u\left(r_{jl}\right)$, $r_{jl}^t=\left|{\mathbf{r}}_j^t-{\mathbf{r}}_l^t\right|,$ the primed $P\text{'}$ denotes the cyclic property $t+1=P+1\equiv 1$, the parameter $\alpha$ is a real number in $\left[\mathrm{0,1}\right],$ and ${\mathsf{\eta }}_{jl}^t$ is the unit vector ${\displaystyle \frac{{\mathbf{r}}_{jl}^t}{r_{jl}^t}}.$ Note that the interactions between beads occur if and only if the beads share the same t-label, and that a recommended value for $\alpha$ is $\alpha =1/3$ [12]$.$ Furthermore, the nonequivalence between even- and odd-numbered beads is apparent; it influences not only the structural calculations but also the thermodynamic ones, although in quite different practical forms (e.g., in thermodynamic evaluations using $\alpha =1/3,$ all the beads play a role) [12,14,96,97]. The SCVJ propagator is accurate up to $O\left(P^{-5}\right)$ and its associated partition function is up to $O\left(P^{-4}\right)$. For the quantum hard-sphere fluid, the effective potential $W_{NP}$ is built with the Cao-Berne propagator (CBHSP) [76] that, for hard spheres with mass $m=m_{HS}$ and diameter $\sigma ,$ can be cast as [81,82,97]:

$$\begin{gathered} W_{NP}^{CBHSP}={\displaystyle \frac{{(m}_{HS})P}{2{\beta }^2{\hslash }^2}}{\sum }_{j=1}^N{\sum }_{t=1}^{P^{\text{'}}}{\left({\mathbf{r}}_j^t{-\mathbf{r}}_j^{t+1}\right)}^2+{\displaystyle \frac{1}{P}}{\sum }_{j<l}{\sum }_{t=1}^P{u}_{HS}\left(r_{jl}^t\right)- \\ {\displaystyle \frac{1}{\beta }}ln{\prod }_{j<l}{\prod }_{t=1}^{P^{\text{'}}}\left\{1-{\displaystyle \frac{\sigma (r_{jl}^t+r_{jl}^{t+1}-\sigma )}{r_{jl}^t r_{jl}^{t+1}}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{{(m}_{HS})P}{2\beta {\hslash }^2}}{(r}_{jl}^t-\sigma \left)\right(r_{jl}^{t+1}-\sigma \left)\right(1+\mathrm{c}\mathrm{o}\mathrm{s}{ \phi }_{jl}^{t,t+1})\right)\right\}, \end{gathered}$$

(10)

where, once again, the cyclic property implied by $P\text{'}$ applies to the corresponding sum and product that run over t. The pairwise approach for the interactions between equal-t beads is used; $u_{HS}$ is the usual hard-sphere potential that vanishes for $r_{jl}^t>\sigma$ and becomes infinity for $r_{jl}^t<\sigma ,$ and ${ \phi }_{jl}^{t,t+1}$ is the angle defined by ${\mathbf{r}}_{jl}^t$ and ${\mathbf{r}}_{jl}^{t+1}.$ The equivalence among all the beads in the sample is clear in this case. Also, there are no specific rules for the accuracy of pair-action propagators; their effectiveness is to be checked by increasing P, but they are extremely efficient in reducing the optimal number of beads [9,70,71,76,106,107,108]. Equations (9) and (10) illustrate the analogies and differences between the effective potentials $W_{NP}$ for both types of systems. Both possess a formally identical first term on the right-hand side that contains the image of P-membered closed necklaces (i.e., free-particle contributions [2,4]). However, the rest of the terms clearly differ from one another, owing to: (a) the characteristics of the interparticle potentials involved in their constructions, which bring about the appearance of forces in SCVJ [12] and of kinetic correlation effects in CBHSP (the contributions from adjacent beads in different necklaces) [97]; and (b) the bead symmetry and asymmetry present in CBHSP and SCVJ, respectively.

3. Quantum Fluid Structures

The whole set of equilibrium structures of a monatomic fluid in the thermal quantum diffraction regime can be deduced from the basic form given in Equation (8). In achieving this task, the related developments are to be complemented with the linear response considerations of external weak fields $\Psi$ acting on the fluid [14,19,32,78,79]. Thus, one finds the following classes and associations: (a) total thermalized-continuous linear response TLRn, associated with continuous external fields ${\Psi }^{\left(N\right)}={\sum }_j{\Psi (\mathbf{r}}_j);$ (b) centroids CMn, as a particular case of TLRn, since its corresponding continuous field acts as ${\Psi }_F^{\left(N\right)}={\sum }_j{\mathbf{f}·\mathbf{r}}_j,$where f is a constant strength; and (c) instantaneous ETn, associated with singular fields that cause the localization of the thermally delocalized quantum atoms (nuclei) in the fluid sample, i.e., the fields used in elastic X-ray diffraction or in the inelastic scattering of neutrons [20,21,22,24], for which collision processes are involved [27,79]. (As stressed earlier, zero-spin-atom BE fluids also follow the same systematics [19,78,79]).

Remarkably, the classes CMn and ETn keep a classical-like pattern in that the thermal quantum delocalization of the atoms (nuclei), although accounted for in their developments at every structural n-level, is not obviously patent in the resulting analytic equations for their structure factors [9,19]. In sharp contrast, the class TLRn includes nontrivial particle self-correlations in the formulations at every n-level [14,19,32],

The analysis involving fields completes the standard probabilistic method for the definition of structures [9,12,13,26]. For conceptual reasons [26,27,28,37,45], these quantum developments can be extended fully if they are carried out in the grand canonical ensemble [14,19,30,32], although the use of the canonical ensemble also serves the primary purpose of establishing basic formal definitions of the structures [26,27,28]. In this connection, it is worthwhile to stress that, albeit the parallel reasoning employed in the classical domain, based on functional calculus operating with $\Psi$, is an excellent guide [26,37,45], the quantum complexity does not allow one to follow such a procedure in its entirety. The specific nature of $\Psi$ may make certain classical-like manipulations either meaningful or void in the quantum domain. Moreover, some special mathematical objects that are perfectly defined in the classical domain, as associated with the grand canonical ensemble (i.e., the direct correlation functions $c_{Cn})$[37,45], may not be defined in an exact manner for every quantum structural case [19,32]. Given that, for computational reasons, this work will be focused on CM3 and ET3, attention to their main theoretical features is given in this Section.

3.1. The Centroid Structures

3.1.1. Opening Centroid Facts

The centroid concept in quantum statistical mechanics arose within Feynman’s path integral formulation (PI) of the equilibrium behavior of many-body monatomic systems at nonzero temperatures [1,2]. For the centroid definition to be made, thermal quantum diffraction effects were required to completely dominate the system behavior. The centroid associated with atom j is defined as the “center of mass” ${\mathbf{R}}_{CM,j}$ of its representative closed path ${\mathbf{r}}_j\left(\mathsf{\tau }\right)$ that, for this purpose, is always considered to have a uniform “mass-density” along its contour. Thus, the centroid concept can be formulated per atom j as follows [1,2]:

$${\mathbf{R}}_{CM,j}={\displaystyle \frac{1}{\beta \hslash }}{\int }_0^{\beta \hslash }{\mathbf{r}}_j\left(\tau \right)d\tau ; j=\mathrm{1,2},\dots ,N,\dots .$$

(11)

The same definition may be applied to every delocalized particle composing neat many-body systems, so long as they may be described by one-site models (e.g., liquid deuterium [109], liquid para-H₂ [97,110,111,112], liquid N₂ [113], the hard-sphere fluid or solid [29,81,82,97], etc.).

Application of a variational principle leads to approximate (semiclassical) forms of the partition function in terms of the foregoing centroid variables [1,2,114,115]. The main result of these approximations is contained in the centroid-effective interparticle potentials $w_F\left(R_{jl};\beta ,\hslash ,m\right)$ that shape the variational semiclassical partition functions (intercentroid distance $R_{jl}=\left|{\mathbf{R}}_{CM,j}-{\mathbf{R}}_{CM,l}\right|).$ The $w_F$ potentials correct the (usual) interatomic potential energies $u\left(r\right)$ but differ dramatically from them in that they depend on the parameters $\beta ,\hslash ,$and m. There is, obviously, formal equivalence among the centroids that can be defined in this simplified way to model a given monatomic fluid. However, it is worth realizing that these objects are not representations of the true particles/atoms of the fluid [14]; in fact, centroids mimic a fluid at a higher density than the actual one, and the conclusions drawn from their only use do not necessarily apply to the properties of the quantum system under study [12,13,14,80,85,90,94,114,115,116,117,118].

The best illustration of the foregoing fact is, perhaps, given by the centroid structures, which can be directly calculated from the centroid-effective partition functions in classical-like ways (MC or MD) [14,84,85,117]. The determination of the actual particle correlation structures does need special convolutions [84,85,117]; the latter involve particle thermal quantum packets (i.e., representations of a delocalized particle around its centroid, which may be nontrivial Gaussians) and the centroid structures, producing approximations to the true path integral TLRn [14,85]. These PI centroid-based approaches are far more useful than the well-known Wigner-Kirkwood expansion [67,119,120], even when the latter is extended to highly sophisticated forms [121,122,123]. Nonetheless, the thermodynamic and structural predictive powers of the $w_F$ effective potential pictures are limited: they cannot yield complete results, nor lead to accurate pictures of the real fluid under increasing quantum effects [14,85].

Despite the incompleteness of the centroid-effective approximations, the usefulness of the centroid concept as seen from the full PI perspective exceeds expectations [12,19,34,78,90,94,112] (see below). The starting point is the discretized version of the PI formulation, Equation (8), which involves a faithful probability density (i.e., nonnegative everywhere). Note that no variational approximation is involved here and, therefore, this “exact” realization of the centroid concept is not the same as that arising from the $w_F$ schemes [85]. Through PIMC and PIMD simulations, one can gather statistics related to the n-body general structural quantities (correlation functions $g_n/G_n$ and structure factors $S^{\left(n\right)})$, in particular the centroid structures [19,32,33,34,35,78,97]. The grand canonical ensemble derivations for the centroid CMn class need to consider a field ${\Psi }_F$ of constant strength f$,$ and its main features are given below [19,32,78,94].

3.1.2. PI Centroid Formulations

The PI grand canonical partition function in the presence of ${\Psi }_F$ can be written as:

$${\mathsf{\Xi }}_{PI}^D\left(\mu ,\mathcal{V},T;{\Psi }_F\right)={\sum }_{N\ge 0}\mathrm{e}\mathrm{x}\mathrm{p}\left(\beta \mu N\right) Z_{PI}^D\left(N,\mathcal{V},T;{\Psi }_F^{\left(N\right)}\right),$$

(12)

For the SCVJ and CBHSP frameworks, the centroid variables are given by [12,14]:

$${\mathbf{R}}_{CM,j}={\displaystyle \frac{2}{P}}{\sum }_{t=\mathrm{1,3},\dots ,P-1}{\mathbf{r}}_j^t, \mathrm{f}\mathrm{o}\mathrm{r}\mathrm{S}\mathrm{C}\mathrm{V}\mathrm{J},$$

(13a)

$${\mathbf{R}}_{CM,j}={\displaystyle \frac{1}{P}}{\sum }_{t=\mathrm{1,2},\dots ,P}{\mathbf{r}}_j^t, \mathrm{f}\mathrm{o}\mathrm{r}\mathrm{C}\mathrm{B}\mathrm{H}\mathrm{S}\mathrm{P},$$

(13b)

where each definition can be viewed as the proper discretization of Equation (11). Therefore, the general form of the partition function can be reformulated as [19,30,78]:

$$\begin{gathered} {\mathsf{\Xi }}_{PI}^D\left({\Psi }_F\right)=\sum _{N\ge 0}{\displaystyle \frac{{\mathfrak{z}}^N}{N!}}\int \left\{\prod _{j=1}^{N}d{\mathbf{R}}_j\right\}\int \prod _{t=1}^{P}d{\mathbf{r}}^{N,t}\times \left\{\prod _{j=1}^N\mathsf{\delta }\left({\mathbf{R}}_j-{\mathbf{R}}_{CM,j}\right)\right\}\times \\ \mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta W_{NP}\right)\times \exp \left(-\beta {\sum }_{j=1}^N{\Psi }_F({\mathbf{R}}_j)\right). \end{gathered}$$

(14)

In Equation (14), $\mathfrak{z}=\exp \left(\beta \mu \right){\left(\sqrt{P}{\lambda }_B^{-1}\right)}^{3P},$ the de Broglie wavelength is ${\lambda }_B=h/\sqrt{2\pi m{k}_{B}T}$, and $W_{NP}$ stands for the corresponding PI effective potential. The foregoing form of the partition function is akin to a (semi-)classical partition function, and it is ready to carry out the functional derivatives with respect to the field, i.e., ${\displaystyle \frac{{\delta }^{n}ln{\mathsf{\Xi }}_{PI}^D\left({\Psi }_F\right)}{\delta {\Psi }_F\left({\mathbf{R}}_1\right)\delta {\Psi }_F\left({\mathbf{R}}_2\right){\dots \mathsf{\delta }\Psi }_F\left({\mathbf{R}}_n\right)}},$ by applying the standard classical procedures [27,28,45,77]. Accordingly, the structural results for centroids are of a clear classical-like nature [19,30,78]. Hereafter, the grand canonical ensemble averages at zero field $({\Psi }_F=0)$ will be denoted by $〈\dots 〉$, and those in the canonical ensemble by ${〈\dots 〉}_{Z_{PI}^D}.$The first three functional derivatives of Equation (14) with respect to ${\Psi }_F\left({\mathbf{R}}_j\right)$ lead to the centroid structural functions up to the triplet level. They can be cast as [19]:

$$\left(-k_{B}T\right){\displaystyle \frac{\delta ln{\mathsf{\Xi }}_{PI}^D\left({\Psi }_F\right)}{\delta {\Psi }_F\left({\mathbf{R}}_1\right)}}={\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\Psi }_F\right)={\mathsf{\Gamma }}_{CM1}\left({\mathbf{R}}_1;{\Psi }_F\right),$$

(15a)

$${\left(-k_{B}T\right)}^2{\displaystyle \frac{{\delta }^{2}ln{\mathsf{\Xi }}_{PI}^D\left({\Psi }_F\right)}{\delta {\Psi }_F\left({\mathbf{R}}_1\right){\mathsf{\delta }\Psi }_F\left({\mathbf{R}}_2\right)}}=\left(-k_{B}T\right){\displaystyle \frac{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\Psi }_F\right)}{\delta {\Psi }_F\left({\mathbf{R}}_2\right)}}={\mathsf{\Gamma }}_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_2;{\Psi }_F\right),$$

(15b)

$${\left(-k_{B}T\right)}^3{\displaystyle \frac{{\delta }^{3}ln{\mathsf{\Xi }}_{PI}^D\left({\Psi }_F\right)}{\delta {\Psi }_F\left({\mathbf{R}}_1\right){\mathsf{\delta }\Psi }_F\left({\mathbf{R}}_2\right){\mathsf{\delta }\Psi }_F\left({\mathbf{R}}_3\right)}}=\left(-k_{B}T\right){\displaystyle \frac{\delta {\mathsf{\Gamma }}_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_2;{\Psi }_F\right)}{\delta {\Psi }_F\left({\mathbf{R}}_3\right)}}={\mathsf{\Gamma }}_{CM3}\left({\mathbf{R}}_1,{\mathbf{R}}_2,{\mathbf{R}}_3;{\Psi }_F\right).$$

(15c)

Equations (15) define the generalized quantities ${\mathsf{\Gamma }}_{CMn}$ $\left(n=\mathrm{1,2},3\right),$ which contain the in-the-field correlation functions $g_{CMn\text{'}}\left({\mathbf{R}}_1,..,{\mathbf{R}}_{n\text{'}};{\Psi }_F\right)$ up to the very same upper level n $\left(n^{\text{'}}\le n\right).$ Application of Yvon’s linear response developments [27,77] implies to take ${\Psi }_F=0$ on the right-hand sides of Equations (15). Thus, the trivial result at $n=1$ can be cast as:

$${\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\Psi }_F=0\right)={\rho }_N={\displaystyle \frac{〈N〉}{\mathcal{V}}},$$

(16a)

where ${\rho }_N$ is the mean number density for the isolated-from-$\Psi$ quantum fluid. At the levels $n=2$ and $n=3,$ one obtains the generalized zero-field centroid structures that, with the use of the total correlation function $h_{CM2}\left(R\right)=g_{CM2}\left(R\right)-1,$ can be written as [19]:

$${\mathsf{\Gamma }}_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_2;{\Psi }_F=0\right)={\rho }_N^2{h}_{CM2}\left(R_{12}\right)+{\rho }_N\delta ({\mathbf{R}}_1-{\mathbf{R}}_2),$$

(16b)

$$\begin{gathered} {\mathsf{\Gamma }}_{CM3}\left({\mathbf{R}}_1,{\mathbf{R}}_2,{\mathbf{R}}_3;{\Psi }_F=0\right)= \\ {\rho }_N^3\left[g_{CM3}\left(R_{12},R_{13},R_{23}\right)-h_{CM2}\left(R_{12}\right)-h_{CM2}\left(R_{13}\right)-h_{CM2}\left(R_{23}\right)-1\right]+ \\ {\rho }_N^2\left[h_{CM2}\left(R_{12}\right)\delta ({\mathbf{R}}_2-{\mathbf{R}}_3)+h_{CM2}\left(R_{13}\right)\delta ({\mathbf{R}}_1-{\mathbf{R}}_2)+h_{CM2}\left(R_{23}\right)\delta ({\mathbf{R}}_1-{\mathbf{R}}_3)\right]+ \\ {\rho }_N\delta ({\mathbf{R}}_1-{\mathbf{R}}_3)\delta \left({\mathbf{R}}_2-{\mathbf{R}}_3\right), \end{gathered}$$

(16c)

where the $g_{CM2}$ and $g_{CM3}$ correlation functions are analogous to the usual two- and three-body correlation functions in the theory of classical homogeneous and isotropic fluids [26,27,28]. The latter are defined by the following ensemble averages involving Equation (12) at ${\Psi }_F=0:$

$${{\rho }_N^{2}g}_{CM2}\left(R_{12}\right)=〈\sum _{j\ne l}\delta \left({\mathbf{R}}_{CM,j}-{\mathbf{q}}_1\right)\delta \left({\mathbf{R}}_{CM,l}-{\mathbf{q}}_2\right)〉,$$

(17)

$${{\rho }_N^{3}g}_{CM3}\left(R_{12},R_{13},R_{23}\right)=〈\sum _{j\ne l\ne m\ne j}\delta \left({\mathbf{R}}_{CM,j}-{\mathbf{q}}_1\right)\delta \left({\mathbf{R}}_{CM,l}-{\mathbf{q}}_2\right)\delta \left({\mathbf{R}}_{CM,m}-{\mathbf{q}}_3\right)〉,$$

(18)

where the generic intercentroid distances are defined in terms of auxiliary q variables, i.e., $R_{ab}=\left|{\mathbf{q}}_a-{\mathbf{q}}_b\right|.$ These spatial averages can be determined via PI simulation by following closely procedures parallel to those applied in the study of classical monatomic fluids [29,30,31,32,33,34,41,46,48].

The formulations of the static structure factors, the pair $S_{CM}^{\left(2\right)}\left(\mathbf{k}\right)$ and the triplet $S_{CM}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right),$ are obtained through the Fourier transforms of ${\mathsf{\Gamma }}_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_2;{\Psi }_F=0\right)$ and ${\mathsf{\Gamma }}_{CM3}\left({\mathbf{R}}_1,{\mathbf{R}}_2,{\mathbf{R}}_3;{\Psi }_F=0\right),$respectively. Thus, one finds [19]:

$$S_{CM}^{\left(2\right)}\left(\mathbf{k}\right)=1+{\rho }_N\int d\mathbf{R}\exp \left(i\mathbf{k}·\mathbf{R}\right)h_{CM2}\left(R\right),$$

(19)

$$\begin{gathered} S_{CM}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)=1+{\rho }_N\left[h_{CM2}\left({\mathbf{k}}_1\right)+h_{CM2}\left({\mathbf{k}}_2\right)+h_{CM2}\left({\mathbf{k}}_1+{\mathbf{k}}_2\right)\right]+ \\ {\rho }_N^2\left[g_{CM}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)-h_{CM2}\left({\mathbf{k}}_1\right){\left(2\pi \right)}^3\delta \left({\mathbf{k}}_1+{\mathbf{k}}_2\right)-h_{CM2}\left({\mathbf{k}}_2\right){\left(2\pi \right)}^3\delta \left({\mathbf{k}}_1\right)-\right. \\ -h_{CM2}\left({\mathbf{k}}_1\right){\left(2\pi \right)}^3\delta \left({\mathbf{k}}_2\right)-\left.{\left(2\pi \right)}^3\delta \left({\mathbf{k}}_1\right){\left(2\pi \right)}^3\delta \left({\mathbf{k}}_2\right)\right]. \end{gathered}$$

(20)

In the foregoing formulas, the k wave vectors define the momentum transfers from the field to the fluid (i.e., $\hslash \mathbf{k}).$ The operative expressions take the forms $S_{CM}^{\left(2\right)}\left(k\right)$ and ${ S}_{CM}^{\left(3\right)}\left(k_1,k_2,\mathrm{c}\mathrm{o}\mathrm{s}({\mathbf{k}}_1,{\mathbf{k}}_2\right))\equiv S_{CM}^{\left(3\right)}\left(k_1,k_2,\varphi \right),$ in which a modulus-variable $k=\left|\mathbf{k}\right|$ denotes the corresponding wavenumber and $\varphi$ is the angle between the wave vectors ${\mathbf{k}}_1$ and ${\mathbf{k}}_2.$ From the PI computational standpoint, the centroid pair and triplet structure factors can be fixed through the ensemble averages [18,19,32,96]:

$$S_{CM}^{\left(2\right)}\left(\mathbf{k}\right)\to {\displaystyle \frac{1}{〈N〉}}〈\sum _{j=1}^N\sum _{l=1}^N\mathrm{e}\mathrm{x}\mathrm{p}\left(i\mathbf{k}·\left[{\mathbf{R}}_{CM,j}-{\mathbf{R}}_{CM,l}\right]\right)〉,$$

(21)

$$S_{CM}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)\to {\displaystyle \frac{1}{〈N〉}}〈\sum _{j=1}^N\sum _{l=1}^N\sum _{m=1}^N\mathrm{e}\mathrm{x}\mathrm{p}\left[i\left({\mathbf{k}}_1·\left({\mathbf{R}}_{CM,j}-{\mathbf{R}}_{CM,m}\right)+{\mathbf{k}}_2·\left({\mathbf{R}}_{CM,l}-{\mathbf{R}}_{CM,m}\right)\right)\right]〉,$$

(22)

where the reader’s attention should be drawn to the fact that the $\delta \left(\mathbf{k}\right)$-evaluations at $\mathbf{k}=0,{\mathbf{k}}_1=0,{\mathbf{k}}_2=0,$ and ${\mathbf{k}}_1+{\mathbf{k}}_2=0$ are avoided, as they are not directly accessible in simulation work [27,46].

At higher orders ($n\ge 3),$ the process can be iterated giving the generalized correlation functions ${\mathsf{\Gamma }}_{CMn}$ that describe, from the centroid standpoint, the loss of homogeneity of the quantum fluid brought about by ${\Psi }_F$. In general, one finds the linear response relationships $(n\ge 1)$:

$${\left(-k_{B}T\right)}^{n+1}{\displaystyle \frac{{\delta }^{n+1}ln{\mathsf{\Xi }}_{PI}^D\left({\Psi }_F\right)}{\delta {\Psi }_F\left({\mathbf{R}}_1\right)\dots {\mathsf{\delta }\Psi }_F\left({\mathbf{R}}_n\right){\mathsf{\delta }\Psi }_F\left({\mathbf{R}}_{n+1}\right)}}=\left(-k_{B}T\right){\displaystyle \frac{\delta {\mathsf{\Gamma }}_{CMn}\left({\mathbf{R}}_1,\dots ,{\mathbf{R}}_n;{\Psi }_F\right)}{\delta {\Psi }_F\left({\mathbf{R}}_{n+1}\right)}}\approx {\mathsf{\Gamma }}_{CM,n+1}\left({\mathbf{R}}_1,\dots ,{\mathbf{R}}_{n+1};{\Psi }_F=0\right),$$

(23)

$$\delta {\mathsf{\Gamma }}_{CM}^{\left(n\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2,\dots ,{\mathbf{k}}_n;{\Psi }_F\right)\approx -\beta \delta {\Psi }_F\left({\mathbf{k}}_1+{\mathbf{k}}_2+\dots +{\mathbf{k}}_n\right) {\mathsf{\Gamma }}_{CM}^{\left(n+1\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2,\dots ,{\mathbf{k}}_n;{\Psi }_F=0\right),$$

(24)

$$S_{CM}^{(n+1)}\left({\mathbf{k}}_1,{\mathbf{k}}_2,\dots ,{\mathbf{k}}_n;{\Psi }_F=0\right)={\displaystyle \frac{1}{{\rho }_N}}{\mathsf{\Gamma }}_{CM}^{\left(n+1\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2,\dots ,{\mathbf{k}}_n;{\Psi }_F=0\right).$$

(25)

Equation (24) results from Fourier transforming the center-right-hand side of Equation (23), which takes advantage of the translational invariance [28] of the fluid structural quantity ${\mathsf{\Gamma }}_{CM,n+1}({\Psi }_F=0$). As seen, in Fourier space, the variation $\delta {\mathsf{\Gamma }}_{CM}^{\left(n\right)}$ is a linear functional of the field applied. Also, note that the centroid structures are intermediate quantities, since the application of a weak field of constant strength would lead to the TLR2 response function as the, in principle, measurable object [19,30,32]; centroids as such do not couple physically with external fields, actual particles do.

From the practical standpoint, the correlation function and structure factor calculations are more affordable using the canonical ensemble, involving a sample size $N_S$ in a box of volume ${\mathcal{V}}_S$ at temperature T. In relation to this, the corresponding averages for the centroid structures given above, Equations (17) to (22), only need to reflect the changes associated with the reduction $〈\dots 〉\to {〈\dots 〉}_{Z_{PI}^D}$; for example, the triplet Equation (22) in the canonical ensemble reads as:

$$S_{CM}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)\to {\displaystyle \frac{1}{N}}{〈\sum _{j=1}^N\sum _{l=1}^N\sum _{m=1}^N\mathrm{e}\mathrm{x}\mathrm{p}\left[i\left({\mathbf{k}}_1·\left({\mathbf{R}}_{CM,j}-{\mathbf{R}}_{CM,m}\right)+{\mathbf{k}}_2·\left({\mathbf{R}}_{CM,l}-{\mathbf{R}}_{CM,m}\right)\right)\right]〉}_{Z_{PI}^D},$$

(26)

As for the zero-wavenumber situations $(k\to 0)$, costly PI extrapolation procedures must be employed to obtain approximate answers: by increasing consistently $N_S$ and ${\mathcal{V}}_S,$ so as to keep the number density ${\rho }_N=N_S/{\mathcal{V}}_S$ constant [46], lower k can be reached and extrapolations to zero wavenumbers may be obtained.

3.1.3. The Exact Centroid OZn Framework

In the quantum diffraction regime, the important point regarding PI centroids, Equations (13), is that they behave according to the same rules for structures applicable to classical monatomic fluids at equilibrium [27,40,45]. This is of great practical importance when addressing k-space questions. In fact, it has been demonstrated [19] that the (hierarchies of) centroid direct correlation functions and structure factors $\left\{c_{CMn},S_{CM}^{\left(n\right)}\right\}$ can be defined in the same fashion as in the classical domain. Therefore, the underlying classical Ornstein-Zernike framework (OZn) [37,38,40,45,124,125,126,127,128], which is consistently defined in the grand canonical ensemble, can be safely applied to the PI centroid correlations and their structure factors. Hence, one may deal with the centroid zero-wavevector problems considered above in an essentially “exact” manner. In particular, at the pair OZ2 and triplet OZ3 levels, in the absence of the external field (${\Psi }_F=0),$ one can write rigorously the following basic equations for PI centroids [19,30]:

$$h_{CM2}\left(R_{12}\right)=c_{CM2}\left(R_{12}\right)+{\rho }_N\int d{\mathbf{R}}_3{c}_{CM2}\left(R_{13}\right)h_{CM2}\left(R_{23}\right),$$

(27)

$${\left({\displaystyle \frac{\partial c_{CM2}(R_{12};{\rho }_N)}{\partial {\rho }_N}}\right)}_T=\int d{\mathbf{R}}_3{c}_{CM3}\left({\mathbf{R}}_1,{\mathbf{R}}_2,{\mathbf{R}}_3;{\rho }_N\right),$$

(28a)

$${\left({\displaystyle \frac{\partial c_{CM}^{\left(2\right)}({\mathbf{k}}_1;{\rho }_N)}{\partial {\rho }_N}}\right)}_T=c_{CM}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2=0;{\rho }_N\right),$$

(28b)

$$S_{CM}^{\left(2\right)}\left(\mathbf{k}\right)={\left(1-{\rho }_N{c}_{CM}^{\left(2\right)}\left(\mathbf{k}\right)\right)}^{-1},$$

(29)

$$S_{CM}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)=S_{CM}^{\left(2\right)}\left({\mathbf{k}}_1\right)S_{CM}^{\left(2\right)}\left({\mathbf{k}}_2\right)S_{CM}^{\left(2\right)}\left(\left|{\mathbf{k}}_1{+\mathbf{k}}_2\right|\right)\left\{1+{\rho }_N^2{c}_{CM}^{\left(3\right)}({\mathbf{k}}_1,{\mathbf{k}}_2)\right\},$$

(30)

where $c_{CM}^{\left(2\right)}\left(\mathbf{k}\right)=c_{CM}^{\left(2\right)}\left(k\right)$ and $c_{CM}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)=c_{CM}^{\left(3\right)}\left(k_1,k_2,\mathrm{c}\mathrm{o}\mathrm{s}({\mathbf{k}}_1,{\mathbf{k}}_2)\right)$ are the Fourier transforms of the functions $c_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_2\right)=c_{CM2}\left(R_{12}\right)$ and $c_{CM3}\left({\mathbf{R}}_1,{\mathbf{R}}_2,{\mathbf{R}}_3\right)=c_{CM3}\left(R_{12},R_{13},R_{23}\right),$respectively. Equation (27) is the pair-level Ornstein-Zernike equation that defines the pair direct correlation function $c_{CM2}\left(R_{12}\right),$ which is a short-ranged function that shows a quick decay to zero with increasing distances [46]. Equations (28) belong to Baxter’s exact hierarchy, under uniform changes of density [38,45], for the direct correlation functions in the quantum $c_{CMn}$ extension [19]. Equations (29) and (30) are the analytically closed reformulations of the pair and triplet structure factors defined in Equations (19) and (20). The advantage of Equations (29) and (30) is clear: none of them explicitly contains the simulation-intractable $\delta \left(\mathbf{k}\right)$ terms [46] that appear in the “raw” analytic formulations based on the Fourier transforms of the $g_2$and $g_3$ functions.

Note that all the structural functions depend on the density and the temperature, although such dependence is only written when necessary [20,27,60,61]. Equations (27)-(30) are mathematically exact and yield the formal solutions to the problematic zero-wavevector situations: (a) at the pair level, via $c_{CM}^{\left(2\right)}\left(k=0\right)$ [78]; and (b) at the triplet level, via Equation (28b) plus the symmetry property $c_{CM}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)=c_{CM}^{\left(3\right)}\left({\mathbf{k}}_1,-{\mathbf{k}}_1-{\mathbf{k}}_2\right)$ for the cases ${\mathbf{k}}_1=-{\mathbf{k}}_2$ [45].

However, once the function $h_{CM2}\left(R_{12}\right)$ in real space is known (via PI simulation), the question arises as to how to determine the pair and triplet structure factors using the direct correlation functions $c_{CM2}\left(R_{12}\right)$ and $c_{CM3}\left(R_{12},R_{13},R_{23}\right).$ These latter functions are defined by integral equations, and a hierarchical structure is in the way. In this regard, there are theoretical methods derived in the classical domain that, being based on closures, provide answers to these key problems [14,40,45,46,125,126,127,128]. These methods are highly accurate and reliable at the level $n=2$, and they can also be successfully applied to quantum calculations (regardless of the structural class!) [14,78,85,96,97]. However, the same favorable situation is not generally found at the level $n=3$, neither in the classical case [44,45,47] nor in the quantum case [19,33,35]. In any event, one has two complementary routes to compute the foregoing structure factors: (a) the computationally “exact”, though incomplete (and expensive), based on the PI simulation of Equations (21)-(22) -or better, of their more affordable translations into the canonical ensemble-; and (b) the theoretically exact, though subjected to the closure uncertainties in the triplet case, based on the general OZn framework summarized in Equations (27)-(30).

An observation of practical value is in order here. The standard PI simulation of the structural functions in r-space using the canonical ensemble, i.e., $Z_{PI}^D\left(N,V,T;\Psi =0\right)$, also poses well-known theoretical drawbacks (e.g., inadequate asymptotic behaviors and finite sample size) [19,81]. To cope with this general problem [26,46,126,127], some corrections [127,128] allow one to improve the pair-level canonical results, thus yielding good approximations to the grand canonical pair structures [34,35,81,96,97]. These corrections are based on the use of the Ornstein-Zernike OZ2 framework and go both ways: r-space $\leftrightarrows$ k-space, that is, involving jointly the pair of structures $\left\{g_{CM2},S_{CM}^{\left(2\right)}\right\},$ thereby benefiting the whole centroid pair structural determination issue.

3.1.4. The Centroid Usefulness

Although centroid structures are not directly accessible via experimental techniques, it is important to make some further comments on the global usefulness of the centroid quantities. In addition to the abovementioned $w_F$ effective potential approximations, one finds the following applications of the set ${\left\{{\mathbf{R}}_{CM,j}\right\}}_{j=\mathrm{1,2},\dots ,N,\dots }$:

(i) Keeping track of the closed paths (or necklaces) associated with the atoms j (or the one-site particles) throughout general PI simulations. Without any loss of generality, the closed paths can be expressed as ${\mathbf{r}}_j\left(t\right)={\mathbf{R}}_{CM,j}+{\mathsf{\xi }}_j\left(t\right),$ where use of the auxiliary general path-vector ${\mathsf{\xi }}_j\left(t\right)$ is made $(t$ runs over the bead labels $\mathrm{1,2},3,\dots ,P).$ Thus, the displacements of the particles $∆{\mathbf{r}}_j\left(t\right)$ can be referred to those of their moving centroids $∆{\mathbf{R}}_{CM,j}$ plus those of the closed paths ${∆\mathsf{\xi }}_j\left(t\right)$ around the centroids. This also yields the useful PI image of the centroid-constrained paths for every particle j [6,81,82,85,90,91,92,93,94,95,129].

(ii) Fixing, via OZ2-$c_{CM2}\left(R\right),$ the equations of state for fluids with thermal quantum diffraction behavior [19,78,81], since in the grand canonical ensemble one finds:

$$S_{CM}^{\left(2\right)}\left(k=0\right)={\displaystyle \frac{〈N^2〉-{〈N〉}^2}{〈N〉}}{=\rho }_N{k}_{B}T{\chi }_T,$$

(31a)

where ${\chi }_T=$ $-{\displaystyle \frac{1}{\mathcal{V}}}{\left({\displaystyle \frac{\partial \mathcal{V}}{\partial p}}\right)}_T$ is the isothermal compressibility.

(iii) Formulating higher-order number fluctuations, as for example [19]:

$$S_{CM}^{\left(3\right)}\left({\mathbf{k}}_1=0,{\mathbf{k}}_2=0\right)={\displaystyle \frac{〈{\left(N-〈N〉\right)}^3〉}{〈N〉}}={\displaystyle \frac{〈N^3〉}{〈N〉}}-3〈N^2〉+2{〈N〉}^2,$$

(31b)

where the exact centroid relationships with OZ2 and OZ3 are to be noticed [32].

(iii) Characterizing global order in quantum samples through the use of centroid adapted quantities: Steinhardt et al’s parameters [83], configurational pair structure factor information, etc. [6,34,82,97].

(iv) Designing approaches to deal with quantum dynamics (e.g., centroid molecular dynamics [90,91,92,112] and other approximations [80,94]).

(v) Dealing with coarse-graining approaches in quantum statistical mechanics [95].

(vi) Studying the number fluctuations under zero-spin BE conditions [19,78,79]. As a result of the action of an external field of constant strength ${\Psi }_F,$ given the algebraic group character of the $N!$ permutations, one can derive a partition function closely similar to (14) involving the conventional PI centroids given by Equations (13):

$$\begin{gathered} {\mathsf{\Xi }}_{PI}^{BE}\left({\Psi }_F\right)=\sum _{N\ge 0}{\displaystyle \frac{\mathrm{e}\mathrm{x}\mathrm{p}\left(\beta \mu N\right)}{N!}}\int \left\{\prod _{j=1}^{N}d{\mathbf{R}}_j\right\}\int \prod _{t=1}^{P}d{\mathbf{r}}^{N,t}\times \left\{\prod _{j=1}^N\mathsf{\delta }\left({\mathbf{R}}_j-{\mathbf{R}}_{CM,j}\right)\right\}\times \\ \underset{{\mathsf{\Omega }}_N^{BE}}{\underbrace{{{\sum }_{{\mathcal{P}}_N}C}_{N,{\mathcal{P}}_N}\mathrm{e}\mathrm{x}\mathrm{p}\left\{-\beta W_{NP}\left({\mathcal{P}}_N\right)\right\}}}\times \exp \left(-\beta {\sum }_{j=1}^N{\Psi }_F({\mathbf{R}}_j)\right), \end{gathered}$$

(32)

where the coefficients within the permutation sum are $C_{N,{\mathcal{P}}_N}>0$, and the probability density function at zero field behaves as ${\mathsf{\Omega }}_N^{BE}\ge 0$ [19,78]. Consequently, with the proviso that ${\mathsf{\Omega }}_N^{BE}$ is used, number fluctuations can be counted as in Equations (31).

3.2. The Instantaneous Structures

Conceptually, the instantaneous ETn class is linked to the context of elastic scattering of radiation. The pair level ET2 can be understood by following the standard quantum treatments [21,27] of X-ray diffraction and inelastic neutron scattering. In X-ray scattering [22,27], the collisions of the X-ray photons with the electrons of the atoms define the nuclei positions, thereby yielding the actual pair quantum structures associated with this phenomenon: the radial correlation function $g_{Q2}\left(r\right),$ and the structure factor $S_Q^{\left(2\right)}\left(k\right)$ [9]. In neutron scattering [21,24,27], the neutron-nucleus interactions are described through Fermi’s pseudopotential (a $\delta -$collision process); this time-dependent quantum phenomenon is characterized by the dynamic structure factor that, through its “elastic reduction” (frequency sum rule), also yields the same pair structures $g_{Q2}\left(r\right)$ and $S_Q^{\left(2\right)}\left(k\right)$ [27]. For both experimental techniques, the corresponding arguments belong to the general linear response theory [27], hence the statistical mechanics treatments give, in the end, the same response function from the fluid: $S_Q^{\left(2\right)}\left(k\right)$. Such a response is formulated in terms of $g_{Q2}\left(r\right)$ in the absence of the external field [79], the $S_Q^{\left(2\right)}\left(k\right)$ mathematical form being essentially the Fourier transform of $(g_{Q2}\left(r\right)-1)$ [21,22,27] (i.e., a “classical-like” formulation). Thus, one finds the instantaneous ensemble average recipes:

$$S_Q^{\left(2\right)}\left(\mathbf{k}\right)=1+{\rho }_N\int d\mathbf{r}\exp \left(i\mathbf{k}·\mathbf{r}\right)\left(g_{Q2}\left(r\right)-1\right),$$

(33a)

$${{\rho }_N^{2}g}_{Q2}\left(r_{12}\right)=〈\sum _{j\ne l}\delta \left({\mathbf{r}}_j-{\mathbf{q}}_1\right)\delta \left({\mathbf{r}}_l-{\mathbf{q}}_2\right)〉,$$

(33b)

where $r=r_{12}=\left|{\mathbf{q}}_1-{\mathbf{q}}_2\right|.$ If the PI formalism is used, $g_{Q2}\left(r_{12}\right)$ directly translates into the following approximation [9,19,30,32]:

$$g_{Q2}\left(r_{12}\right)\approx g_{ET2}\left(r_{12}\right)={\displaystyle \frac{1}{X}}〈{\sum }_{t^{*}=1}^X\sum _{j\ne l}\delta ({\mathbf{r}}_j^{t^{*}}-{\mathbf{q}}_1)\delta ({\mathbf{r}}_l^{t^{*}}-{\mathbf{q}}_2)〉,$$

(34a)

where the $t^{*}$ notation for beads is retained, and the associated structure factor reads as:

$$S_{Q2}^{\left(2\right)}\left(\mathbf{k}\right)\approx S_{ET}^{\left(2\right)}\left(\mathbf{k}\right)=1+{\rho }_N\int d\mathbf{r}\exp \left(i\mathbf{k}·\mathbf{r}\right)\left(g_{ET2}\left(r\right)-1\right),$$

(34b)

where the grand canonical partition function ${\mathsf{\Xi }}_{PI}^D$ at zero field is involved [19]:

$${\mathsf{\Xi }}_{PI}^D\left(\mu ,\mathcal{V},T;\Psi =0\right)=\sum _{N\ge 0}{\mathrm{e}\mathrm{x}\mathrm{p}\left(\beta \mu N\right) Z}_{PI}^D\left(N,\mathcal{V},T;\Psi =0\right).$$

(35)

(The forms of Equations (33) and (34) can be applied to both the grand canonical and the canonical ensembles [27], provided that the proper partition function is employed).

Although the basic PI canonical partition function $Z_{PI}^D$, Equation (8), already defines the beads that are going to be significant in the structure calculations, centroid-like functional calculus manipulations are inapplicable to the instantaneous class; this fact forces the $\Psi =0$ condition in Equation (35). If one wished to apply functional derivatives, setting $\Psi \ne 0,$ to obtain the generalized structure functions ${\mathsf{\Gamma }}_{ETn}$ $\left(n=\mathrm{1,2},3,\dots \right)$, a fundamental difficulty would arise: the collision phenomena localize the actual quantum particles, which is not compatible with the structure of $Z_{PI}^D$ [19,79]. Within the current $\left\{{\mathbf{r}}_j\right\}$ framework of nucleus positions, one faces the “collapse” of the atom quantum thermal packets, and the functional derivatives turn out to be meaningless; in fact, if this sort of manipulations were “formally” carried out, one would obtain the total thermalized-continuous linear response structures TLRn [19,32]. The functions $g_{ETn},$ ${\mathsf{\Gamma }}_{ETn},$ and $S_{ET}^{\left(n\right)}$ at levels $n>2$ can be defined (using $\Psi =0)$ by PI-adapting their classical counterparts (e.g., at $n=3,$ Equations (16c), (18), (20), etc.) [19]. The reasons behind this procedure are: (a) it is a straightforward consequence of the averaging of the dynamical n-body number density functions/operators $\left(\prod \delta \left(\mathbf{r}-\mathbf{q}\right)\right)$ that define the elemental $g_{Qn}$ structures in real space [28] (e.g., see Equation (33b)); and (b) consistency with the transition to the classical limit at every n-level for ${\mathsf{\Gamma }}_{ETn}$ [19]. In this work, attention will be focused on the triplet-level instantaneous structure functions $g_{ET3}(r_{12},r_{13},r_{23})$ and $S_{ET}^{\left(3\right)}\left(k_1,k_2,\cos \left({\mathbf{k}}_1,{\mathbf{k}}_2\right)\right).$As remarked earlier, grand canonical PI simulations for obtaining structures pose a hard problem because of the computational effort involved. In general, the ETn class computations scale with $X$ (i.e., P/2 for SCVJ, or P for CBHSP), and the practical alternative is based on the canonical ensemble applications [9,19,33,34,96]. Therefore, depending on the propagator selected, the instantaneous triplet structure averages read as [19,30,32]:

- SCVJ (fourth-order propagator)

$${{\rho }_N^{3}g}_{ET3}\left(r_{12},r_{13},r_{23}\right)={{\displaystyle \frac{2}{P}}〈\sum _{odd-t}\sum _{j\ne l\ne m\ne j}\delta \left({\mathbf{r}}_j^t-{\mathbf{q}}_1\right)\delta \left({\mathbf{r}}_l^t-{\mathbf{q}}_2\right)\delta \left({\mathbf{r}}_m^t-{\mathbf{q}}_3\right)〉}_{{\mathrm{Z}}_{PI}^D},$$

(36)

$$S_{ET}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)={\displaystyle \frac{2}{NP}}{〈\sum _{odd-t}\sum _{j=1}^N\sum _{l=1}^N\sum _{m=1}^N\mathrm{e}\mathrm{x}\mathrm{p}\left[i\left({\mathbf{k}}_1·\left({\mathbf{r}}_j^t-{\mathbf{r}}_m^t\right)+{\mathbf{k}}_2·\left({\mathbf{r}}_l^t-{\mathbf{r}}_m^t\right)\right)\right]〉}_{Z_{PI}^D}.$$

(37)

- CBHSP (pair action)

$${{\rho }_N^{3}g}_{ET3}\left(r_{12},r_{13},r_{23}\right)={{\displaystyle \frac{1}{P}}〈{\sum }_{t=1}^P\sum _{j\ne l\ne m\ne j}\delta \left({\mathbf{r}}_j^t-{\mathbf{q}}_1\right)\delta \left({\mathbf{r}}_l^t-{\mathbf{q}}_2\right)\delta \left({\mathbf{r}}_m^t-{\mathbf{q}}_3\right)〉}_{{\mathrm{Z}}_{PI}^D},$$

(38)

$$S_{ET}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)={\displaystyle \frac{1}{NP}}{〈\sum _{t=1}^P\sum _{j=1}^N\sum _{l=1}^N\sum _{m=1}^N\mathrm{e}\mathrm{x}\mathrm{p}\left[i\left({\mathbf{k}}_1·\left({\mathbf{r}}_j^t-{\mathbf{r}}_m^t\right)+{\mathbf{k}}_2·\left({\mathbf{r}}_l^t-{\mathbf{r}}_m^t\right)\right)\right]〉}_{{\mathrm{Z}}_{PI}^D.}$$

(39)

The classical-like OZn framework for the instantaneous structures is not exact but an approximation [19,32,78]. At the pair level, however, OZ2 and its corrections to the canonical ensemble results have proven to be highly accurate in the thermal quantum diffraction regime [81,85,96,97,117]. Therefore, applications of OZ3 to the instantaneous triplet structures are well-worth exploring in depth [19,33,35]. In this context, the related basic equations can be translated from the centroid class and written as $\left(r_{jl}=\left|{\mathbf{q}}_j-{\mathbf{q}}_l\right|\right)$:

$$h_{ET2}\left(r_{12}\right)\approx c_{ET2}\left(r_{12}\right)+{\rho }_N\int d{\mathbf{q}}_3{c}_{ET2}\left(r_{13}\right)h_{ET2}\left(r_{23}\right),$$

(40)

$${\left({\displaystyle \frac{\partial c_{ET2}(r_{12};{\rho }_N)}{\partial {\rho }_N}}\right)}_T\approx \int d{\mathbf{q}}_3{c}_{ET3}\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3;{\rho }_N\right),$$

(41a)

$${\left({\displaystyle \frac{\partial c_{ET}^{\left(2\right)}({\mathbf{k}}_1;{\rho }_N)}{\partial {\rho }_N}}\right)}_T\approx c_{ET}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2=0;{\rho }_N\right),$$

(41b)

$$S_{ET}^{\left(2\right)}\left(\mathbf{k}\right)\approx {\left(1-{\rho }_N{c}_{ET}^{\left(2\right)}\left(\mathbf{k}\right)\right)}^{-1},$$

(42)

$$S_{ET}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)\approx S_{ET}^{\left(2\right)}\left({\mathbf{k}}_1\right)S_{ET}^{\left(2\right)}\left({\mathbf{k}}_2\right)S_{ET}^{\left(2\right)}\left(\left|{\mathbf{k}}_1{+\mathbf{k}}_2\right|\right)\left\{1+{\rho }_N^2{c}_{ET}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)\right\}.$$

(43)

(Recall that the use of closures is compulsory at every OZn structural level [37,38,39,40,44,45,46,47,125,126,127,128]).

4. Triplet Closures and Associated Features

Given the classical-like forms of $\left\{g_{CM3},S_{CM}^{\left(3\right)}\right\}$ and $\left\{g_{ET3},S_{ET}^{\left(3\right)}\right\}$ , the closures may be applied indistinctly to both of them. For the calculations carried out in this work, the employed closures are: (a) for triplets in real space, Jackson-Feenberg convolution (JF3) [17] that contains Kirkwood superposition (KS3) [25], and the intermediate average AV3 [34]; and (b) for triplets in Fourier space, Jackson-Feenberg convolution [17,40,45], and the symmetrized version of Denton-Ashcroft approximation (DAS3) [47]. Obviously, the pair information needed arises from PIMC simulations and OZ2 treatments.

The mathematical expressions for the real-space closures are given below.

-

JF3

$$\begin{gathered} g_{JF3}\left(d_{12},d_{13},d_{23}\right)=g_{KS3}\left(d_{12},d_{13},d_{23}\right)-h_2\left(d_{12}\right)h_2\left(d_{13}\right)h_2\left(d_{23}\right)+ \\ {\rho }_N\int d{\mathbf{q}}_4{h}_2\left(d_{14}\right)h_2\left(d_{24}\right)h_2\left(d_{34}\right), \end{gathered}$$

(44a)

where the distances $d_{ab}$ stand for $\left|{\mathbf{R}}_a-{\mathbf{R}}_b\right|$ for centroids or $\left|{\mathbf{q}}_a-{\mathbf{q}}_b\right|$ for actual quantum particles (the atoms of helium-3, or the one-site hard spheres), and $g_{KS3}$ is Kirkwood superposition [25], which reads as:

-

KS3

$$g_{KS3}=g_2\left(d_{12}\right)g_2\left(d_{13}\right)g_2\left(d_{23}\right).$$

(44b)

-

AV3

$$g_{AV3}\left(d_{12},d_{13},d_{23}\right)={\displaystyle \frac{1}{2}}\left(g_{KS3}\left(d_{12},d_{13},d_{23}\right)+g_{JF3}\left(d_{12},d_{13},d_{23}\right)\right).$$

(45)

The Fourier-space closures read as:

-

JF3

$$c_{JF3}\left(d_{12},d_{13},d_{23}\right)=0\to c_{JF}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)=0,$$

(46)

-

DAS3

$$c_{DAS}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)={\displaystyle \frac{1}{3}}\left[c_{DA}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)+c_{DA}^{\left(3\right)}\left({\mathbf{k}}_1,\left|{\mathbf{k}}_1+{\mathbf{k}}_2\right|\right)+c_{DA}^{\left(3\right)}\left({\mathbf{k}}_2,\left|{\mathbf{k}}_1+{\mathbf{k}}_2\right|\right)\right],$$

(47a)

$$c_{DA}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)=c_{DA}^{\left(3\right)}\left(k_1,k_2\right)=$$

(47a)

$${\displaystyle \frac{1}{c^{\left(1\right)\text{'}}}}\left\{c^{\left(2\right)}{(k}_1\left)c^{\left(2\right)\text{'}}{(k}_2\right)+c^{\left(2\right)}{(k}_2\left)c^{\left(2\right)\text{'}}{(k}_1\right)\right\}-{\displaystyle \frac{c^{{\left(1\right)}^{\text{'}\text{'}}}}{{\left(c^{{\left(1\right)}^{\text{'}}}\right)}^2}}{c}^{\left(2\right)}{(k}_1\left)c^{\left(2\right)}{(k}_2\right),$$

(47b)

$$c^{\left(1\right)\text{'}}=c^{\left(2\right)}\left(k=0\right), c^{{\left(1\right)}^{\text{'}\text{'}}}={\left({\displaystyle \frac{\partial c_2\left(k=0\right)}{\partial {\rho }_N}}\right)}_T, { c}^{{\left(2\right)}^{\text{'}}}\left(k\ne 0\right)={\left({\displaystyle \frac{\partial c_2\left(k\right)}{\partial {\rho }_N}}\right)}_T.$$

(47c)

This closure satisfies Baxter’s exact behavior Equation (28b) only at ${\mathbf{k}}_1={\mathbf{k}}_2=0$.

The use of other closures to deal with triplet quantum calculations can be seen in References [19,29,30,31,32,33,34,35]. In particular, when considering Fourier-space applications, AV3 and the elaborate Barrat-Hansen-Pastore factorization ansatz (BHP3) [45] appear to need very detailed numerical treatments [19,45]; this clarifying task is left for future work.

The structure factors $S_{CM}^{\left(3\right)}$ and $S_{ET}^{\left(3\right)}$ given by Equations (30) and (43), respectively, need their corresponding pair and triplet direct correlation functions, and some remarks are necessary here. First, the pair direct correlation functions $c_{CM2}\left(R\right)$ and $c_{ET2}\left(r\right)$ can be determined in a highly accurate way [81,96,97] with the use of the OZ2 numerical procedure put forward by Baxter, Dixon and Hutchinson (BDH) [125,126]. As stated earlier, the BDH application to $c_{CM2}\left(R\right)$ comes from an exact OZ2 framework, whereas to $c_{ET2}\left(r\right)$ is an approximation. Second, the final closure results for $S^{\left(3\right)}$ are, in general, subjected to inaccuracies owing to the hierarchical structure of the direct correlation functions [38], which is illustrated for triplets by Equations (28) and (41)$.$ Third, the formal use of direct correlation functions yields answers to questions that cannot be exactly solved using PI simulations: the behaviors at zero-momentum transfers $\mathbf{k}=0$ and $({\mathbf{k}}_1=0,{\mathbf{k}}_2=0),$ which are contained in the theoretical relationships [19,32,78]:

$$S_{CM}^{\left(2\right)}\left(k=0\right)=S_{ET}^{\left(2\right)}\left(k=0\right)={\rho }_N{k}_{B}T{\chi }_T={\left(1-{\rho }_N{c}_{CM}^{\left(2\right)}\left(k=0\right)\right)}^{-1},$$

(48a)

$$S_{CM}^{\left(3\right)}\left({\mathbf{k}}_1=0,{\mathbf{k}}_2=0\right)=S_{ET}^{\left(3\right)}\left({\mathbf{k}}_1=0,{\mathbf{k}}_2=0\right)=$$

(48b)

$$S_{CM}^{\left(2\right)}(k=0)\left[S_{CM}^{\left(2\right)}\left(k=0\right)+{\rho }_N{\left({\displaystyle \frac{\partial S_{CM}^{\left(2\right)}\left(k=0\right)}{\partial {\rho }_N}}\right)}_T\right].$$

(48b)

It is worthwhile to make some useful observations here. (a) Equations (48) can be extended to include the TLR case [19,32,78], with $S_{TLR}^{\left(2\right)}\left(k=0\right)$ in Equation (48a) and $S_{TLR}^{\left(3\right)}\left({\mathbf{k}}_1=0,{\mathbf{k}}_2=0\right)$ in Equation (48b). (b) Given the identities in Equation (48a), the ET and TLR extensions of Equation (48b) in terms of their pair $S^{\left(2\right)}(k=0$) components are equally valid. (c) Highly accurate estimates of the exact values can be obtained from the OZ2 centroid quantities, as Equations (48) arise from the consideration of the exact role played by $c_{CM}^{\left(2\right)}\left(k=0\right).$ (d) Although the parallel extensions of the whole Equation (48b) to ET and TLR are exact, their practical applications via OZ2 schemes will produce approximations to the exact behavior of the $S^{\left(3\right)}\left(\mathrm{0,0}\right)$ quantity. And (e) by applying the underlying k-space symmetry properties [45], one can also obtain from Equation (28b) the exact values of the triplet CM3 components $\left(\mathbf{k},-\mathbf{k}\right)$:

$$c_{CM}^{\left(3\right)}\left(\mathbf{k},-\mathbf{k}\right)=c_{CM}^{\left(3\right)}\left(-\mathbf{k},\mathbf{k}\right)=c_{CM}^{\left(3\right)}\left(\mathbf{k},0\right)=c_{CM}^{\left(3\right)}\left(0,\mathbf{k}\right)={\left({\displaystyle \frac{\partial c_{CM}^{\left(2\right)}\left(k\right)}{\partial {\rho }_N}}\right)}_T,$$

(48c)

where the presence of single zero-momentum transfers is to be noticed. The general type of symmetry expressed in Equation (48c) allows one to set sum rules for the $S_{CM}^{\left(3\right)}\left(\mathbf{k},0\right)$ and $S_{ET}^{\left(3\right)}\left(\mathbf{k},0\right)$ components [32], which are exact for the centroid case but approximate for the instantaneous case.

5. Computational Details

5.1. State Points

For helium-3, the atom mass is $m_{{}_3\mathrm{H}\mathrm{e}}=3.01603 \mathrm{a}\mathrm{m}\mathrm{u}$ and the critical point is located at ${(T}_C=3.3157 \mathrm{K}$; ${\rho }_{N,C}=0.0082246 {Å}^{-3})$ [130]. The following supercritical state points are studied [131]: SP1$\left(T=4.21 \mathrm{K}; {\rho }_N=0.0228717687 {Å}^{-3}\right),$ SP2$\left(T=4.21 \mathrm{K}; {\rho }_N=0.0272988971 {Å}^{-3}\right),$ SP3$\left(T=8.99 \mathrm{K}; {\rho }_N=0.0228717687 {Å}^{-3}\right),$ and SP4$\left(T=4.2 \mathrm{K}; {\rho }_N=0.02286713 {Å}^{-3}\right).$ At state points SP1, SP2, and SP3, the pair and triplet centroid correlation structures in r-space are investigated. At state point SP4, attention is focused on the evaluation of the instantaneous and centroid triplet response functions in k-space. In this connection, to carry out closure calculations, advantage is taken of previous work along the isotherm $T=4.2 \mathrm{K}$[96], where the pair direct correlation functions were obtained at four state points about SP4 defined by the density variations $\pm ∆{\rho }_N=\pm 0.002{ Å}^{-3}, \pm 0.004{ Å}^{-3}.$For quantum hard spheres, the characterization of the state points only needs two parameters, namely the reduced density ${\rho }_N^{*}={\rho }_N{\sigma }^3$ and the reduced de Broglie wavelength ${\lambda }_B^{*}={\displaystyle \frac{{\lambda }_B}{\sigma }},$ a fact that unifies the description of this special system [34,132,133,134]. The state points on the crystallization line studied are [82]: QHS1$({\lambda }_B^{*}=0.2$; ${\rho }_N^{*}=0.789)$, QHS2$({\lambda }_B^{*}=0.4$; ${\rho }_N^{*}=0.672)$, and QHS3$({\lambda }_B^{*}=0.6$; ${\rho }_N^{*}=0.589)$. For comparison purposes, the actual parameter values used in these calculations are $m_{HS}=28.0134 \mathrm{a}\mathrm{m}\mathrm{u}$ and $\sigma =3.5 Å.$ In this regard, the temperatures and the number densities are: QHS1($T=22.2024 \mathrm{K};{\rho }_N=0.018402332$ ${Å}^{-3}),$ QHS2($T=5.5506 \mathrm{K};{\rho }_N=0.015673469$ ${Å}^{-3}),$ and QHS3($T=2.4669 \mathrm{K};{\rho }_N=0.013737609$ ${Å}^{-3}).$ These state points are ascribed to the quantum diffraction regime. For example, note that at QHS3 one finds a low value for the degeneracy parameter $\gamma ={\left({\lambda }_B^{*}\right)}^3{\rho }_N^{*}\approx 0.13.$ If fluid helium is modelled with QHS, even at ${(\lambda }_B^{*}=1.9832{; \rho }_N^{*}=0.348)$ that is on the crystallization line of this model [82,97,134], one would deal with temperatures far from FD or BE exchange interactions: (a) $T\approx 5.3 \mathrm{K}$in helium-3, or (b) $T\approx 4 \mathrm{K}$ in helium-4 [65,134]. To carry out closure calculations, two additional QHS state points about each of the listed above are studied; ${\lambda }_B^{*}$ is kept constant and density variations $\pm {\mathsf{\Delta }\rho }_N^{*}=\pm 1.715\times {10}^{-2}$ are taken, which for comparison with the density variations in the helium-3 calculations amount to $\pm ∆{\rho }_N=\pm 0.0004 {Å}^{-3}$.

5.2. PIMC Simulations

The PIMC simulations of both systems follow the same patterns explained elsewhere [29,34,35,81,96,97]. However, for the reader’s convenience, the following operating procedures are worth mentioning. The canonical ensemble is selected, using a central cubic box of side L and sample size $N_S\times P$ (see below for the different choices made). The necklace normal-mode algorithm is utilized [75,85,135]. Therefore, any necklace in the simulation sample can move in P different modes: one of them is the translation of the necklace as a whole, and the rest (breathing modes or “collective excitations” of the beads [135]) describe independent ways for the P beads to displace by keeping their centroid fixed (recall the role of the relative path-vector ${\mathsf{\xi }}_j\left(t\right)$). In a uniform sampling of modes, this yields a reduced ratio $\left({\displaystyle \frac{1}{P}}\right)$ for the centroid displacements, which indicates that the sampling of centroid quantities is less effective and will need more PIMC moves to attain reduced error bars. When evaluating structures, this fact becomes more critical in k-space than in r-space, because of the special features of the k-space sampling (see below) [46]. In building the Markov chains, the acceptance ratio for each of the different normal modes is set to 50%. The propagators are: SCVJ $(\alpha ={\displaystyle \frac{1}{3}})$ for helium-3, involving the SAPT2 interatomic potential [100], and CBHSP for quantum hard spheres. Given a simulation sample size $N_S\times P,$ the run lengths are quantified with the use of kpass and Mpass units: 1 kpass = ${10}^3{N}_S\times P$ attempted bead moves, 1 Mpass = ${10}^6{N}_S\times P$ attempted bead moves. The selected P values are known to be optimal [35,81,82,96,97], and the computations are focused on equilateral and isosceles features of the targeted centroid and instantaneous triplet $g_3$ and $S^{\left(3\right)}$ structures. The specific details are given in what follows.

5.2.1. Real r-Space

The PIMC simulations of helium-3 utilize the following $N_S\times P$ values for the calculations of $g_{CM2}$ and $g_{CM3}:$ $1372\times 66$ at SP1, 1372$\times 80$ at SP2, $1372\times 22$ at SP3 (also, $1024\times 66$ at SP4 [96]). (Results for the instantaneous $g_{ET2}$ and $g_{ET3}$ can be found in References [19,35,96]. The sampling of the pair and triplet structures in r-space uses a basic distance spacing $∆=0.1 Å$ for defining the histogram-bins of the functions [29,46]. The run lengths for the pair structures, are: 2000 kpasses at SP1, 1800 kpasses at SP2, and 1300 kpasses at SP3. The run lengths for the equilateral and isosceles triplet structures are: 3660 kpasses at SP1, 3560 kpasses at SP2, and 2560 kpasses at SP3. (The simulation sampling at SP3 is enhanced with respect to Reference [35]: from 500 to 1300 kpasses for $g_{ET2 },$and from 750 to 2560 kpasses for $g_{ET3 }).$ To fix the error bars (one-standard deviation), the run lengths are distributed in superblocks that give partial subaverages, which in turn serve to estimate the corresponding standard deviations. The resulting error bars are small. For example, in the close vicinities of the main peaks (max) obtained in the simulations, one finds: (a) at the pair-level $g_{CM2}(R_{max}$) and $g_{ET2}\left(r_{max}\right)$, the error bars are lower than $0.6\%$ of the corresponding mean heights; and (b) at the triplet-level $g_{CM3}(R_{max},R_{max},R_{max}$) and $g_{ET3}\left(r_{max},r_{max},r_{max}\right),$the error bars remain below $1\%$ of the corresponding mean heights. The pair radial functions $g_{CM2}\left(R\right)$ (also $g_{ET2}\left(r\right)$ at SP3) are further processed with the OZ2 methodology [81,96] for fixing the pair structure factors $S_{CM}^{\left(2\right)}\left(k\right)$ (and $S_{ET}^{\left(2\right)}\left(k\right)$ at SP3); BDH variational calculations [125,126] are carried out, and their accuracy is increased with grand canonical corrections by using five Baumketner-Hiwatari iterations (BHw) [128]).

For the additional PIMC simulations, on the isothermal fluid branches, of the auxiliary QHS state points, i.e., those about the targeted on the crystallization line, the fixing of the pair radial structures $g_{CM2}$ and $g_{ET2}$ uses: $N_S\times P=$ $864\times 12,$ and a spacing $∆=0.1 Å$ as the width of the histogram-bins. The run lengths are circa 3700 kpasses, and the error bars are very small (e.g., below 1% for the main peak heights).

5.2.2. Fourier k-Space

These PIMC simulations utilize the following sample sizes $N_S\times P$: (a) for helium-3, $128\times 66$ (SP4) and (b) for QHS, $250\times 12$ (QHS1, QHS2, and QHS3). Note that for QHS the number of actual particles $N_S=250$ is greater than that of helium-3. The studied QHS state points are on the fluid side of the melting-freezing transition of this system, and the spurious fluid-solid flipping that is known to happen for small $N_S$ values [19,46] is to be avoided. Gaseous helium-3 state point SP4 does not present this problem, and $N_S=128$ serves well the purposes of this investigation.

The PIMC sampling of the triplet structures in k-space, for both helium-3 and QHS, involves the scanning of wave vectors commensurate with the simulation box, i.e., $\mathbf{k}=\left({\displaystyle \frac{2\mathsf{\pi }}{L}}\right)\left(k_x{,k}_y,k_z\right),$ where $k_x{,k}_y,$ and $k_z$ are integers. For the equilateral components $S^{\left(3\right)}(k_q,k_q,{\displaystyle \frac{\pi }{3}})$, twelve sets ${\left\{\left({\mathbf{k}}_{q1},{\mathbf{k}}_{q2}\right);\dots \right\}}_{q=1,\dots ,12}$ are scanned, where each set is composed of eight pairs of wave vectors; these wave vectors have the same modulus, $\left|{\mathbf{k}}_{q1}\right|=\left|{\mathbf{k}}_{q2}\right|=\dots =k_q,$ and every pair defines an angle $\varphi =\pi /3.$ The wave vectors are such that $2\le k_x^2+k_y^2+k_z^2\le 200,$ and it is easy to understand how the equilateral sets are built. For example, taking $k_x^2+k_y^2+k_z^2=32$, a pair of vectors is $\left[\left(\mathrm{0,4},-4),(\mathrm{4,4},0\right)\right]$, another is $\left[\left(\mathrm{0,4},4),(-\mathrm{4,0},4\right)\right]$, and so on. The corresponding ranges of simulated wavenumbers $k_q,$ given in ${Å}^{-1}$ for comparison purposes, turn out to be: (a) for helium-3, $0.5\lesssim k_q/{Å}^{-1}\lesssim 5$ (SP4); and (b) for QHS, $0.74\lesssim k_q/{Å}^{-1}\lesssim 3.72$ (QHS1), $0.70\lesssim k_q/{Å}^{-1}\lesssim 3.53$ (QHS2), and $0.67\lesssim k_q/{Å}^{-1}\lesssim 3.38$ (QHS3). For the isosceles components, only helium-3 at SP4 is investigated. In this regard, the performed analysis covers the components $S_{CM}^{\left(3\right)}(k_m,k_m,\varphi )$ and $S_{ET}^{\left(3\right)}\left(k_m,k_m,\varphi \right)$ defined at wavenumber $k_m=2.123236 {Å}^{-1}.$ (Such value $k_m$ is near the maxima of the pair structure factors $S_{CM}^{\left(2\right)}\left(k_M\right)$ and $S_{ET}^{\left(2\right)}\left(k_M^{\text{'}}\right)$, where $k_M\approx k_M^{\text{'}}\approx 2.06 {Å}^{-1}$[96]). Ten angles ${\left\{{\varphi }_i\right\}}_{i=1,\dots ,10},$ defined by groups of pairs $\left\{\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)\right\}$ of modulus $k_m,$ are scanned in the interval $0\le \varphi <\pi .$ For completeness, an extra eleventh component defined by $k_m^{\text{'}}=2.093539 {Å}^{-1}$ and ${\varphi }_{11}=2.077976\approx 2\pi /3$ is also added. Each angle ${\varphi }_i$ contains in its sampling group a number of independent pairs of wave vectors, which are in between 4 and 16 in this work. It is also easy to see how the groups that define the angles ${\varphi }_i$are formed. For example, taking $k_x^2+k_y^2+k_z^2=36$, isosceles geometries at different angles ${\varphi }_i$ can be obtained by extracting pairs of vectors from $\left[\left(\mathrm{6,0},0);\left(\mathrm{0,0},-6\right);(\mathrm{4,4},2\right);\left(-\mathrm{4,4},-2\right);\dots \right].$ (Recall that Equations (28b) and (41b) serve to determine the values of the $S^{\left(3\right)}\left(\mathrm{0,0}\right)$ components and, also, of the isosceles components at $\varphi =\pi$).

It is important to stress that the PIMC simulations of triplet structure factors are subjected to large fluctuations. Therefore, in the present applications to helium-3 and QHS, the gathering of statistics for the k-space triplet structures is not equally distributed among the different components, as some of them receive an increased scanning attention (e.g., the components of maximal amplitude arising from the simulations).

The helium-3 run lengths are as follows. (a) Centroid case: in between 30 Mpasses and 409 Mpasses for the equilateral components, and in between 53 Mpasses and 278 Mpasses for the isosceles components. (b) Instantaneous case: in between 30 Mpasses and 120 Mpasses for the equilateral components, and in between 30 Mpasses and 48 Mpasses for the isosceles components. Helium-3 statistics are gathered from the Markov chain every 5000 configurations. As for QHS, the run lengths for the equilateral components are as follows. (c) Centroid cases: in between 6 Mpasses and 148 Mpasses. (d) Instantaneous cases: in between 6 Mpasses and 123 Mpasses. QHS statistics are gathered from the Markov chains every 8000 configurations. The statistical errors (one-standard deviation), fixed from block-subaverages, are found to remain controlled. For example, for helium-3 at the most sensitive of the equilateral components analyzed, that is, the obtained maximal peaks of the centroid and instantaneous structure factors (located at $k_{qm}\approx 2 {Å}^{-1}),$ the error bars turn out to be $\approx 0.5\%$ and $1.8\%$ of the mean amplitudes $S_{CM}^{\left(3\right)}(k_{qm},k_{qm},\varphi ={\displaystyle \frac{\pi }{3}})$ and $S_{ET}^{\left(3\right)}\left(k_{qm},k_{qm},\varphi ={\displaystyle \frac{\pi }{3}}\right),$ respectively. As regards QHS, the situation is reasonably good but not so satisfactory: the error bars found for the equilateral maximal peak amplitudes, in both the centroid and the instantaneous cases, are circa 3.4% at QHS1, 5.4% at QHS2, and 5.6% at QHS3. (Complete results can be seen in the next Section).

5.3. Closure Calculations

The closure calculations with JF3 and AV3 in r-space (KS3 is trivial), and JF3 and DAS3 in k- space, are carried out in ways that have been described in detail elsewhere [29,30,31,32,33,34,35]. The main facts of interest to this work are given here.

5.3.1. Helium-3

For helium-3, the following points are to be remarked:

(i) The $g_{CM3}(R,S,S)$ correlations, as approximated by the JF3 and AV3 r-space closures, involve the convolution integral included in Equation (44a). Such a convolution can be expanded as a Legendre polynomial series [29,43,45]; truncation, keeping the first thirty-one terms/polynomials $P_{\nu }\left(x\right),$ is utilized $(-1\le x\le 1;$ $0\le \nu \le 30).$ The pair radial function $g_{CM2}\left(R\right)$ at each of the state points (SP1, SP2, and SP3) is used as data input, and the basic distance spacings for tabulating the triplet $g_3$ correlations are taken equal to 0.1 $Å.$(ii) The pair structure factors $S_{CM}^{\left(2\right)}\left(k\right)$ and $S_{ET}^{\left(2\right)}\left(k\right)$ related to the helium-3 calculations are needed. In relation to this, OZ2 treatments apply BDH+BHw (5-iterations) to the pair radial structures obtained with PIMC [96] (see Reference [35] for calculations of $S_{ET}^{\left(2\right)}\left(k\right)$at SP1, SP2, and SP3). Although, for the triplet structure factors at SP4, the closure JF3 only requires the pair structure factors at this state point, DAS3 closure calculations do need their corresponding OZ2 direct correlation functions at the neighboring four state points to fix the corresponding isothermal density derivatives. The basic information on SP4 at the pair level is taken from the helium-3 study reported in Reference [96]. With regard to all these questions, the following OZ2-related information is worth giving here.

At a given state point, it is known that the combined method BDH+BHw, applied to a pair of functions $\left\{g_2,S^{\left(2\right)}\right\}$ produces generally a sequence of solutions ${\left\{g_2\left(d\right);c_2\left(d\right);S^{\left(2\right)}\left(k\right);R_{Zi}\right\}}_{i=\mathrm{1,2},\dots ,}$ [81]. These solutions are characterized by the values of the trial distances $R_{Zi}$ (the so-called “zeros”) that make $c_2\left(d\ge R_{Zi}\right)\equiv 0$. Practical convergence of the obtained functions is observed as the $R_{Zi}$ increase, and methods to deal with this feature have been discussed [81,96,97]. To carry out the related calculations in this article, the CM2 and ET2 solutions associated with their respective longest $R_{Zi}$ are selected. In this connection, it is worthwhile to compare the results at the most sensitive component, $k=0$ (i.e., the isothermal compressibility ${\chi }_T$, also including the experimental data [131]. These values turn out to be:

(a)

At SP1, ${\chi }_T\left({bar}^{-1}\right)$is 0.005288 for CM2 $(R_Z=18.507 Å)$, 0.005202 for ET2 $(R_Z=16.974 Å)$, with 0.005715 being the experimental value.

(b)

At SP2, ${\chi }_T\left({bar}^{-1}\right)$is 0.002401 for CM2 $(R_Z=17.108 Å)$, 0.002844 for ET2 $(R_Z=15.187 Å)$, with 0.002383 being the experimental value.

(c)

At SP3, ${\chi }_T\left({bar}^{-1}\right)$is 0.005526 for CM2 $(R_Z=16.083 Å)$, 0.005552 for ET2 $(R_Z=16.738 Å)$, with 0.005003 being the experimental value.

The structural improvement achieved herein at SP3 explains the difference from its ET2 value with respect to Reference [35], where ${\chi }_T\left({bar}^{-1}\right)=$ 0.005071 was reported; as seen, the order of magnitude of ${\chi }_T$ does not change and a better agreement between the theoretical CM2 and ET2 estimates is obtained (Equation (48a)). Also, for completeness, the main amplitudes of the centroid and instantaneous pair structure factors at SP4 are [96]: $S_{CM}^{\left(2\right)}\left(k_M=2.054\right)=1.947$, and $S_{ET}^{\left(2\right)}\left(k_M=2.067\right)=1.438.$(iii) The isothermal density derivatives at SP4 needed for the k-space DAS3 calculations (CM3 and ET3) can be calculated through Richardson extrapolation [136], using four selected $c^{\left(2\right)}\left(k;{\rho }_N;R_Z\right)$ functions, one per each of the auxiliary state points about SP4 $(T=4.2 K$ and the uniform density spacing ${∆}_b{\rho }_N=0.002{ Å}^{-3}$) [96]. By taking advantage of such uniform spacing, this algorithm reads as:

$${\left({\displaystyle \frac{\partial c^{\left(2\right)}(k;{\rho }_N)}{\partial {\rho }_N}}\right)}_T\approx {\displaystyle \frac{4D_1\left({∆}_b{\rho }_N\right)-D_1\left(2{∆}_b{\rho }_N\right)}{3}},$$

(49a)

where, for a given density spacing ${\mathbb{s}}_{\rho },$ the simple Stirling numerical estimates $D_1$ are given by:

$$D_1\left({\mathbb{s}}_{\rho }\right)={\displaystyle \frac{c^{\left(2\right)}\left(k;{\rho }_N+{\mathbb{s}}_{\rho }\right)-c^{\left(2\right)}\left(k;{\rho }_N-{\mathbb{s}}_{\rho }\right)}{2{\mathbb{s}}_{\rho }}}.$$

(49b)

In these DAS3 calculations at SP4, the selection of the $c_2$(and its associated $c^{\left(2\right)}$ and $S^{\left(2\right)})$ is made by utilizing these quantities at the longest $R_{Zi}$ zeros obtained via BDH+BHw(5-iterations). Such $R_Z$ zeros fall into $14.3<{\displaystyle \frac{R_Z}{Å}}<17.3$ for the centroid functions, and into $10.2<{\displaystyle \frac{R_Z}{Å}}<14.1$ for the instantaneous functions. Also, it is interesting to note that the values of the (four) real-space pair direct correlation functions at the origin $(d=0)$ are such that:

About SP4:

$$-29.8\lesssim c_{CM2}\left(R=d=0\right)\lesssim -10.9,\mathrm{a}\mathrm{n}\mathrm{d}-9.2\lesssim c_{ET2}\left(r=d=0\right)\lesssim -5,$$

which indicate the depths of the respective $c_2-$“bowls”. To be highlighted are the significantly deeper depth and the larger variation range in the centroid case.

(iv) The basic wavenumber spacings for tabulating the $S^{\left(3\right)}$ components calculated with the closures are taken equal to 0.1 ${Å}^{-1}.$ The evaluations of the double-zero momentum transfer quantities $S^{\left(3\right)}\left(\mathrm{0,0}\right)$ are carried out through: Equation (48b) for the centroid cases using $S_{CM}^{\left(2\right)}\left(k=0\right),$ and its parallel (approximate) relationship for the instantaneous cases using $S_{ET}^{\left(2\right)}\left(k=0\right).$ Richardson extrapolations akin to Equation (49a) are employed for the density derivatives of the pair structure factor $k=0$ values. Alternatively, these quantities can be evaluated using the direct method based on the $c^{\left(2\right)}\left(k=0;{\rho }_N\right)$ values, Equations (28b), (41b), and (49a). These two different methods, (i.e., based on $S^{\left(2\right)}\left(k=0\right)$ or $c^{\left(2\right)}(k=0)),$ are theoretically equivalent, although their numerical estimates may show slight differences because of the underlying density dependence of the input functions.

5.3.2. QHS

Analogously, for the QHS k-space applications, the following items are to be remarked:

(i) For DAS3 calculations, the auxiliary QHS pair structural quantities needed are obtained via application of OZ2 to the PIMC $(N_S\times P=864\times 12)$ pair radial structures fixed in this work; applications of BDH+BHw(5- iterations) to the centroid structures, and of BDH alone to the instantaneous structures, are made. As mentioned earlier, OZ2 for the instantaneous correlations is an approximation, and the related application of the simplified BHw procedure in the region of the fluid-solid change of phase leads to unclear numerical behaviors of the iterations [97]. That is why BHw iterations are not included in the instantaneous analysis at any QHS state point. Whether more advanced grand canonical corrections [127] might cope with this drawback remains to be investigated. (Recall that for the JF3 application at a given state point only its pair structure factor is needed).

(ii) The DAS3 isothermal density derivatives for the cases CM3 and ET3 are calculated through the two-point differentiation Equation (49b). This involves the $c^{\left(2\right)}(k;{\rho }_N;R_Z)$ functions obtained via OZ2, i.e., those arising from the sequences ${\left\{g_2\left(d\right);c_2\left(d\right);S^{\left(2\right)}\left(k\right);R_{Zi}\right\}}_{i=\mathrm{1,2},\dots },$ about each QHS state point investigated (${\lambda }_B=$ constant; basic density spacing $=0.0004{ Å}^{-3},$ or $1.715\times {10}^{-2}$ in reduced units). For each state point on the crystallization line, one representative solution of the foregoing CM2 or the ET2 sequences, is selected from the respective regions of significant convergence. The pair structure factors so fixed at the state points QHS1, QHS2, and QHS3 are employed in the final calculation of their triplet structure factors.

For the sake of comparison, it is worthwhile to quote the following QHS data associated with these OZ2 applications. The solutions selected correspond to $R_Z$ values that fall into the intervals: CM3 applications, $14<{\displaystyle \frac{R_Z}{Å}}<15;$ ET3 applications, $8<{\displaystyle \frac{R_Z}{Å}}<12.5$. These intervals are global in that they contain the $R_Z$ of all the state points: QHS1 and its two neighbors, QHS2 and its two neighbors, etc. Also, the depths of the real-space $c_2(d=0)-$“bowls” of the centroid and instantaneous cases are such that:

(a) About QHS1, $-60.6\lesssim c_{CM2}\left(R=0\right)\lesssim -49.7,$ and $-52.8\lesssim c_{ET2}\left(r=0\right)\lesssim -44.7.$

(b) About QHS2, $-58.7\lesssim c_{CM2}\left(R=0\right)\lesssim -47.3,$ and $-45.1\lesssim c_{ET2}\left(r=0\right)\lesssim -37.3.$

(c) About QHS3, $-58.7\lesssim c_{CM2}\left(R=0\right)\lesssim -46.6,$ and $-36.8\lesssim c_{ET2}\left(r=0\right)\lesssim -30.9.$The depths of these QHS “bowls” are significantly deeper than those at helium-4 state point SP4 (see Reference [19] for a discussion of the apparent negative repercussions of this feature on some closure-scheme evaluations). The wavenumber spacings for tabulating the $S^{\left(3\right)}$ components calculated with the closures are taken equal to $∆k^{*}=0.35$ (or $∆k=0.1 {Å}^{-1}).$ The evaluations of the double-zero momentum transfer quantities $S^{\left(3\right)}\left(\mathrm{0,0}\right)$ are carried out via the density derivatives of $c^{\left(2\right)}\left(k=0;{\rho }_N\right),$ by employing Stirling density derivatives Equation (49b)

(iii) As a complementary test, extended DAS3 calculations are carried out at QHS1 using five densities. On the isothermal fluid branch ${\lambda }_B^{*}=0.2,$ the state points about QHS1 at densities ${\rho }_N^{*}=0.789\pm 4.2875\times {10}^{-3}, {\rho }_N^{*}=0.789\pm 8.575\times {10}^{-3}$ are also studied at the pair level with PIMC $(864\times 12)$ and OZ2. The procedures for fixing the necessary CM2 and ET2 properties are the same as those described earlier, and the impact of this extension on the DAS3 results is discussed in the next Section.

6. Results and Discussion

6.1. Helium-3 Triplet Results in Real r-Space

Figure 1 focuses on relevant equilateral PIMC results in r-space at state points SP1, SP2 and SP3. In Figure 1a one observes that more pronounced centroid features are obtained as the quantum effects get stronger (i.e., by increasing $\gamma ={\lambda }_B^3{\rho }_N),$ that is, SP2 > SP1 > SP3; at constant T (SP1 and SP2), the density marks the relative intensity of the diffraction effects, whilst at constant ${\rho }_N$ (SP1 and SP3) the inverse temperature does it. In this regard, see Figure 1b for a comparison between the centroid and instantaneous equilateral structures at SP1, where for instantaneous triplets (ET3) one notices a smaller distance of approach and a far smoother structure [19], as expected from delocalization arguments.

Figure 2 displays, at state point SP1, typical features of the present centroid triplet correlations. Results obtained with PIMC, AV3, KS3, and JF3, are shown. Figure 2a,b focus on the equilateral $g_{CM3}(R,R,R)$ correlations: PIMC-AV3 comparison in Figure 2a, and PIMC-KS3-JF3 comparison in Figure 2b. As seen, AV3 yields excellent results for the range of significant-correlation distances (i.e., PIMC nonzero $g_{CM3}$ values), whereas neither KS3 nor JF3 does the work within the region of the main maximum (JF3 also fails within the range of short distances $R\lesssim 3 Å;$ this failure influences the AV3 behavior) [34,35,60]). Figure 2c,d focus on some salient isosceles features of $g_{CM3}\left(R,S,S\right)$ by selecting the R-slices $3.55 Å$ and $6.75Å,$ which are located in the vicinities of the equilateral main and secondary maxima, respectively. By discarding the negative influence of JF3 at short range, one observes that the overall performance of AV3 is excellent, albeit, as R increases, AV3 needs improvement in the short-range region of the available S distances (Figure 2d). For completeness, Figure 3 shows results for the PIMC centroid triplet correlations at state points SP2 and SP3, where trends analogous to those discussed above are observed. Besides, the updated results at SP3 obtained for the instantaneous $g_{ET3}(r,r,r)$ do not differ significantly on the scale of the graph from those previously reported in Reference [35], and they are not included in this work. More information on these triplet correlations is contained in several graphs within the Supplementary Material.

6.2. Helium-3 Triplet Results in Fourier k-Space

The computations carried out at state point SP4 are summarized in Figure 4 and Figure 5 and Table 1 and Table 2, where, as expected, the centroid features are more pronounced than the instantaneous ones. In more detail, from the intercomparison of the results obtained with PIMC, JF3, and DAS3, the following points can be remarked.

For the centroid and instantaneous equilateral components, $S_{CM}^{\left(3\right)}(k,k,{\displaystyle \frac{\pi }{3}})$ and $S_{ET}^{\left(3\right)}\left(k,k,{\displaystyle \frac{\pi }{3}}\right),$ the two closures produce results remarkably close to PIMC data, as shown in Figure 4 and Table 1. On the scale of the graphs, the corresponding amplitudes follow patterns only slightly different from each other, and the maxima and minima for $k\gtrsim 2 {Å}^{-1}$ are reasonably well located: see the comparisons PIMC-JF3 in Figure 4a and PIMC-DAS3 in Figure 4b (the resulting spacings $∆k$ used in the PIMC calculations preclude further analysis). A more detailed comparison can be seen in Table 1, which reveals some distinctive features. PIMC yields for the centroid and instantaneous cases negative values up to $k\lesssim 1 {Å}^{-1}$, a trait that is obviously absent from JF3 (owing to its null $c_3)$ but, most interestingly, not from DAS3. In particular, the sign of the double-zero momentum transfer is $S^{\left(3\right)}\left(\mathrm{0,0}\right)=S_{CM}^{\left(3\right)}\left(\mathrm{0,0}\right)<0$, as evaluated via Equation (48b) and the corresponding Richardson derivative of the type (49a). In this connection, it is interesting to compare the results obtained for $S_{CM}^{\left(3\right)}\left(\mathrm{0,0}\right)$ and $S_{ET}^{\left(3\right)}\left(\mathrm{0,0}\right)$ when using the two methods available: the based on the pair single zero-components $S_{CM}^{\left(2\right)}(k=0;{\rho }_N)$ and $S_{ET}^{\left(2\right)}\left(k=0;{\rho }_N\right),$ and the based on the pair single zero-components $c_{CM}^{\left(2\right)}(k=0;{\rho }_N)$ and $c_{ET}^{\left(2\right)}\left(k=0;{\rho }_N\right).$ Thus, one finds consistent results that indicate a negative value of the actual $S^{\left(3\right)}\left(\mathrm{0,0}\right):$

(a)

Method $S^{\left(2\right)}\left(k=0;{\rho }_N\right),$, Equation (48b) and its instantaneous parallel,

$$S_{CM}^{\left(3\right)}\left(\mathrm{0,0}\right)=-0.03211, \mathrm{a}\mathrm{n}\mathrm{d}{S}_{ET}^{\left(3\right)}\left(\mathrm{0,0}\right)=-0.02153.$$

(b)

Method $c^{\left(2\right)}\left(k=0;{\rho }_N\right),$, Equations (28b) and (41b),

$$S_{CM}^{\left(3\right)}\left(\mathrm{0,0}\right)=-0.02755, \mathrm{a}\mathrm{n}\mathrm{d}{S}_{ET}^{\left(3\right)}\left(\mathrm{0,0}\right)=-0.02111.$$

In addition, the differences between JF3 and DAS3 are noticeable within the regions of low wavenumbers. DAS3 definitely takes negative values within: $k\lesssim 1.1 {Å}^{-1}$ in the centroid case, and $k\lesssim 0.3 {Å}^{-1}$ in the instantaneous case. The two closures JF3 and DAS3 stay very close to each other for wavenumbers: $k\gtrsim 2.3 {Å}^{-1}$ in the centroid case, and $k\gtrsim 1.5 {Å}^{-1}$ in the instantaneous case.

Figure 5 and Table 2 show the isosceles components $S_{CM}^{\left(3\right)}(k_m,k_m,\cos \varphi )$ and $S_{ET}^{\left(3\right)}\left(k_m,k_m,\cos \varphi \right)$, where $k_m\approx 2.12{ Å}^{-1}$ is near the maximum amplitude positions of the respective pair structure factors $k_M\approx 2.06 {Å}^{-1}$ [96]. PIMC results are compared to those of JF3 and DAS3, and one observes in Figure 5 that both closures yield similar representations of the PIMC centroid and instantaneous behaviors; the closures appear to work better for smoother structures, as displayed in Figure 5b that contains the instantaneous case. By separating the $\varphi$-angular dominion into two regions, ${\mathsf{\Phi }}_1\left(-1\le \cos \varphi \lesssim 0.4\right)$ and ${\mathsf{\Phi }}_2\left(0.4\lesssim \cos \varphi \le 1\right)$, the following features are to be noticed. Within ${\mathsf{\Phi }}_1$, both closures capture reasonably well the shapes of the triplet CM3 and ET3 functions; recall that the main amplitudes located near $\cos \varphi =-0.5$ are fixed at wavenumber $k_m^{\text{'}}\approx 2.09 {Å}^{-1},$ which is slightly different from $k_M$ and $k_m.$ Nevertheless, within ${\mathsf{\Phi }}_2$, both closures depart from PIMC, which appears to show a levelled behavior (see also Table 2): the closures display larger amplitude values and a remarkably undulating behavior in the centroid case. (For isosceles components at low wavenumbers, DAS3 yields negative values within certain angular ranges; not having PIMC-calculated components to be compared with, these results are not reported here).

The current closure results for the instantaneous equilateral components are somewhat similar to those obtained with Barrat-Hansen-Pastore (BHP3) closure in Reference [35]. The BHP3 isosceles components in ${\mathsf{\Phi }}_2$ also showed a levelled behavior that significantly stayed below most of the JF3 results. However, initial BHP3 calculations of the centroid triplet components at SP4 lead to not so accurate results, which may be linked to the deeper bowl of the centroid pair direct correlation function quoted earlier (see Reference [19] for a related discussion). Therefore, no BHP3 results are reported in this article.

6.3. QHS Triplet Results in Fourier k-Space

Selected PIMC numerical results for the three state points, QHS1, QHS2 and QHS3, are given in Table 3. (More detailed numerical results are contained in the Supplementary Material). For visualization purposes, Figure 6 collects the computed equilateral components of the centroid and instantaneous triplet structure factors $S^{\left(3\right)}\left(k^{*},k^{*},{\displaystyle \frac{\pi }{3}}\right)$; PIMC, DAS3 and JF3 results are plotted with the use of reduced units ${(k}^{*}=k\sigma ).$ Interestingly, at QHS3, the current PIMC mean maximal values are consistent with those previously reported in Reference [19]. However, the extended sampling lengths involved in the current computations, which double those in Reference [19], are not sufficient for reducing the corresponding error bars. Although the convergence of the Monte Carlo sampling is expected to behave with the inverse square root of the number of trials, the insensitiveness obtained is to be ascribed to the adverse role of the strong triplet fluctuations in k-space. This points to the expected magnitude of the task ahead when tackling this sort of triplet study on quantum crystallization lines. Also, recall that DAS3 calculations are based on the the use of Equations (28b), (41b) and (49b), involving the groups of three state points studied with PIMC+OZ2 along each of the isotherms. For the three targeted fluid state points, one obtains values $S_{CM}^{\left(3\right)}\left(\mathrm{0,0}\right)<0,$ which agrees with the independent similar behavior seen at the fluid helium-3 state point SP4 earlier; this behavior will be considered ed in more detail later on.

In these QHS applications, the following equilateral component $S^{\left(3\right)}\left(k^{*},k^{*},{\displaystyle \frac{\pi }{3}}\right)$ features are to be remarked:

(i)

As compared to PIMC, the graphs indicate that DAS3 and JF3 appear to be slightly shifted towards larger wavenumbers (see also the Supplementary Material). The differences between DAS3 and JF3 are most noticeable within the low $k^{*}$-wavenumber region, since DAS3 gives negative values, whereas JF3 cannot yield that behavior. Despite the inaccuracies in the closure results, the current representations of the triplet structure factors turn out to be very valuable approximations to PIMC results. In the behavior of the equilateral components calculated with closures, JF3 takes over as $k^{*}$increases (JF3 appears to yield a limit behavior).

(ii)

From the graphs, one can observe at each state point that the triplet structure factor components show more pronounced features in the centroid $S_{CM}^{\left(3\right)}$cases than in the instantaneous $S_{ET}^{\left(3\right)}$cases. Furthermore, as the quantum effects increase (i.e., larger $\gamma ={\lambda }_B^{*3}{\rho }_N^{*},$QHS1<QHS2<QHS3), one observes that: the CM3-ET3 differences increase, and the structures are shifted towards lower wavenumbers. All of these traits are consistent with their analogous at the pair level when comparison of $S_{CM}^{\left(2\right)}\left(k\right)$and $S_{ET}^{\left(2\right)}\left(k\right)$is made [97].

(iii)

As regards the centroid absolute maxima that can be derived from the equilateral PIMC results, if one uses a simple quadratic fitting of the mean amplitudes in the related regions, the obtained values hint (once again [33]) at the existence of a “constancy” in the maximal centroid equilateral amplitudes. The related centroid estimates turn out to be: (a) 11.153 (QHS1), (b) 11.259 (QHS2), and (c) 11.559 (QHS3). Roughly speaking, this would amount to having an interval $\approx 11.3\pm 0.2$for the absolute maxima of the equilateral centroid components of the QHS model, which could be applicable within the present range of conditions investigated.

(iv)

The foregoing centroid pilot result might be related to a hypothetical maximal $S_{CM}^{\left(3\right)}(k_M,k_M^{\text{'}},{\varphi }_M)$-constancy for real systems that, because of the existence of triple point, show a liquid-solid coexistence line T-bounded from below (e.g., para-hydrogen, neon) [32,33,137,138]. However, to ascertain this question, one should carefully analyze the following issues:

(a) The combined effects of a more exhaustive PI sampling, seeking to reduce the present error bars $(\approx 3.5-5.5\%)$, and a closer spacing $∆k^{*}$ for the fixing of the interpolating points; the sets of commensurate wave vectors should also include general triangular geometries.

(b) The possible validity of such constancy hypothesis under stronger quantum effects; this implies to extend the study to a wider range of QHS fluid-solid coexistence conditions [82,134].

(c) The situations presented by systems that, described by realistic interatomic potentials, present a noteworthy discordance between their associated crystallization lines, the modelled and the genuine [33,97].

This search seems appealing, in that if proven successful, it would extend to the quantum triplet level the classical Hansen-Verlet’s rule for crystallization at the pair level, which gives the interval 2.85$\pm 0.1$ for the main maximum amplitudes of the fluid pair structure factors $S_C^{\left(2\right)}(k$) [27,139]. This sort of result could be very useful for condensed matter studies (e.g., complex real systems, in which interparticle repulsions dominate their behaviors), even if the quantum intervals found were dependent on the specific fluid and no general rule could be established (e.g., the observed trend at the quantum pair level [97]). In addition, other centroid triplet signatures of a quantum fluid-solid change of phase may emerge from this search [34], and this possibility should not be disregarded.

(v)

In relation to the instantaneous maximal equilateral components, by quadratically fitting the PIMC data, as done in the centroid case, one obtains the maximal amplitude estimates: 10.724 (QHS1), 10.060 (QHS2), and 9.238 (QHS3). Therefore, the maxima of the instantaneous components $S_{ET}^{\left(3\right)}(k^{*},k^{*},{\displaystyle \frac{\pi }{3}})$appear to decrease monotonically with increasing quantum effects. Regardless of the magnitude of the error bars and of the $∆k^{*}$spacing available for defining these maximum regions, this behavior agrees with the general pattern of less-structured instantaneous functions arising when going from higher to lower temperatures and densities (QHS1 to QHS3) [33,34,97].

(vi)

The existence of equilateral negative values for low wavenumbers is an observable fact for both classes of components, $S_{CM}^{\left(3\right)}$and $S_{ET}^{\left(3\right)}.$From the current calculations, PIMC results show just a glimpse of this behavior for $k^{*}\lesssim 2.6$(Table 3), whilst DAS3 gives such a behavior within larger ranges, which go from $k^{*}\lesssim 4.9$at QHS1-CM3 to $k^{*}\lesssim 3.5$at QHS3-ET3 ${(∆k}^{*}=0.35).$Moreover, the double-zero momentum transfer components turn out to be negative, with the exception of $S_{ET}^{\left(3\right)}\left(\mathrm{0,0}\right)=+0.00005$at QHS1. Also, see the Supplementary Material for illustrative results at QHS3.

(vii)

Given that density-differentiation algorithms more accurate than Equation (49b) can influence greatly the final estimates, particularly in situations involving zero-momentum transfers (i.e., values $k^{*}=0),$this point deserves further consideration. Thus, it is worth noting that the additional DAS3 five-point closure calculations at QHS1, based on Equation (49a), do not show alterations in the magnitude of the foregoing CM3 and ET3 ranges of negative values resulting from Equation (49b). However, use of Equation (49a) yields at QHS1: $S_{CM}^{\left(3\right)}\left(\mathrm{0,0}\right)=-1.81\times {10}^{-4}$and $S_{ET}^{\left(3\right)}\left(\mathrm{0,0}\right)=2.5\times {10}^{-5},$which are to be compared with the $k^{*}=0$entries in Table 3 at QHS1: $-0.0011$and 0.00005, respectively. (Recall that the calculation of the exact $S^{\left(3\right)}\left(\mathrm{0,0}\right)$-component is obtained through the CM3 scheme). Finally, for completeness, the effect of the differentiation algorithm on DAS3 results becomes less important as $k^{*}$increases, which can be observed within the maximal amplitude regions $k^{*}\approx k_M^{*}.$For example, at QHS1 for $k_M^{*}=6.65,$one sees that:

(a) Use of Equation (49b) yields $S_{CM}^{\left(3\right)}\left(k_M^{\mathrm{*}},k_M^{\mathrm{*}},{\displaystyle \frac{\pi }{3}}\right)=11.2462$, whilst use of Equation (49a) yields $S_{CM}^{\left(3\right)}\left(k_M^{\mathrm{*}},k_M^{\mathrm{*}},{\displaystyle \frac{\pi }{3}}\right)=11.0842.$

(b) Use of Equation (49b) yields $S_{ET}^{\left(3\right)}(k_M^{\mathrm{*}},k_M^{\mathrm{*}},{\displaystyle \frac{\pi }{3}})=10.4487$, whilst use of Equation (49a) yields $S_{ET}^{\left(3\right)}(k_M^{\mathrm{*}},k_M^{\mathrm{*}},{\displaystyle \frac{\pi }{3}})=10.4435.$

7. Final Remarks and Perspectives

The present investigation on quantum fluid structures adds new triplet results to the knowledge of fluids composed of helium-3 atoms and of quantum hard spheres, for which initial information can be found in previous works by this author. Both systems are studied under thermal quantum diffraction conditions, although the respective fluid state points differ qualitatively: supercritical for helium-3, and on the crystallization line for quantum hard spheres. Path integral Monte Carlo simulations and closure computations, at the pair and triplet levels of the centroid and instantaneous structures, are performed in the real and the Fourier spaces, and the following remarks can be made.

In real space, the path-integral triplet structures (equilateral and isosceles configurations) obtained for helium-3 show that the centroid structure features are more pronounced than their instantaneous counterparts, in agreement with the general expected behavior since centroids are known to mimic a fluid at a higher density than the actual one. The related whole path-integral computational load depends critically on the simulation sample size $N_S\times P$ -canonical ensemble-, which includes the number of necklaces (particles) $N_S$ and the number of beads $P;$ even under strong quantum diffraction effects, this task is known to be perfectly affordable at the pair level, and, by focusing on the abovementioned specific types of triplet configurations, it is also at the triplet level (run lengths $~3000$kpasses). Most interestingly, the application of the real-space intermediate closure AV3 to the determination of centroid and instantaneous triplet structures allows one, via comparison with exact path integral data, to confirm the great usefulness of this approximation: by analyzing structural conditions more extreme (the centroid correlations in helium-3) than those previously studied elsewhere, the AV3 closure performance is, once again, excellent and far less expensive than full path integral simulations.

AV3 is the straightforward average of two basic closures: Kirkwood superposition and Jackson-Feenberg convolution; both closures show clear deficiencies in the treatment of the triplet features in real space. Therefore, although AV3 partly inherits the unphysical short-range-distance failure of the Jackson-Feenberg convolution, the compensating effect of the abovementioned average makes AV3 capture most of the salient features of the path-integral centroid correlations (equilateral and isosceles). Such a surprising triplet fact in real space was observed in recent works by this author (e.g., centroid and instantaneous triplets in the quantum hard-sphere fluid, and instantaneous triplets in supercritical helium-3). Within the significant correlation ranges of distances (i.e., those showing nonzero path-integral structure values), the AV3 applicability to fluids with quantum behavior appears to cover diffraction effects characterized by substantial values of the degeneracy parameter $\gamma$, including singular model systems (e.g., quantum hard spheres for $\gamma \lesssim 2.7)$ and real systems (e.g., helium-3 for $\gamma \lesssim 3.2).$ However, note that AV3 worsens as the density increases along isotherms, and it does not capture relevant maximum region behaviors in the isosceles configurations either. All in all, closures may give good and cost-effective “pictures” of complex interrelations among the actual atoms/particles in fluid phases. In this regard, the further centroid triplet results obtained in this work indicate that, in the study of quantum fluids via closures, the improvement of AV3 is well-worth exploring, which should also incorporate general three-body configurations in the study. When trying to fix the abovementioned density effects that AV3 suffers from, the time-honored Abe’s developments may be a good guidance.

In the study reported in this work, the Fourier space applications to the fluid centroid and instantaneous triplet structures (equilateral and isosceles wave-vector geometries) illustrate the nature of these highly demanding calculations. In general, path integral Monte Carlo simulations present strong fluctuations in the structure factor component values along the Markov chain; this fact becomes critical in the vicinities of the main peak amplitudes investigated for both the centroid and the instantaneous structure factors. This imposes very long run lengths ($~{10}^2-4\times {10}^2$ Mpasses) to achieve reasonably reduced statistical errors in such regions. Such a hard computational situation is aggravated by some factors: (a) the number and sizes of the sets of commensurate wave vectors, which are needed to obtain significant descriptions of the components analyzed; and (b) the conditions under study, which put lower limits to the simulation sample size $N_S\times P$ to be utilized (i.e., in general, changes of phase influence $N_S$ and diffraction effects influence $P$). In this connection, selective sampling of the path-integral structure factor components (e.g., equilateral main peak amplitudes) can be used for increasing the accuracy of the results in specific regions. In any event, currently, the whole quantum fluid triplet study in Fourier space via path integrals, even in the diffraction regime, with its diversity of classes and 4-D intricate numerical behavior, poses a formidable challenge in every respect. Accordingly, the use of closures in Fourier space is a useful way to obtain insights into the quantum triplet structure factor issues.

The current path integral results clearly indicate that, for low and not large wavenumbers, the Fourier equilateral components take on (small) negative values in both the centroid and the instantaneous cases. This exact equilateral feature is observed at every fluid state point investigated in this work (helium-3 and quantum hard spheres). However, as regards the present closures, such a feature is only captured by the symmetrized Denton-Ashcroft and not in its full details (Jackson-Feenberg convolution cannot give this behavior). Thus, with respect to these below-zero regions, the centroid triplets are in general better represented than the instantaneous triplets by such Denton-Ashcroft closure, a feature to be ascribed to the classical-like structural nature of centroids (i.e., the exact theoretical applicability of Ornstein-Zernike schemes, in particular the pair OZ2, to centroids). In addition, as compared to the data obtained via path integrals, both closures give reasonably good equilateral results beyond the below-zero regions, with Jackson-Feenberg convolution closure becoming dominant for large wavenumbers.

Continuing with the equilateral components, the path-integral centroid calculations in Fourier space for the quantum hard-sphere fluid lead, via interpolation in the computed data tables, to equilateral main-peak amplitude estimates that suggest a “constancy” in these centroid quantities along the section of the crystallization line studied. This surmise needs more elaboration before stating its validity and range of applicability. Therefore, extended simulation work is needed to: (a) reduce the current statistical errors, and (b) analyze more extreme quantum conditions. (As an aside, note that the analogous maxima of the instantaneous components cannot be candidates for “constancy”, because they show a definite decay pattern with increasing quantum effects). Moreover, other triplet centroid features in Fourier space (e.g., the maxima of salient isosceles components, as those near angles $\varphi \approx 2\pi /3,$ etc.) might be subjected to this sort of “constancy” search. Even though the analysis of these possibilities entails costly simulations, this issue seems well-worth studying; its potential applicability to complex real systems, as a triplet way to provide more insights into quantum freezing phenomena, should be sufficient motivation to justify the effort. Furthermore, in the search for triplet-structure signatures of distinctive quantum fluid behaviors, the role of the centroid component $S_{CM}^{\left(3\right)}\left(\mathrm{0,0}\right)$ should not be overlooked (recall its relationship with the third-order number fluctuations). The present results for the fluid phases analyzed give $S_{CM}^{\left(3\right)}\left(\mathrm{0,0}\right)<0,$ a fact that indicates a negative skewness for all of the particle number (marginal) distributions involved in this work. Therefore, the ranges of values that this quantity may take under different thermodynamic conditions deserve further investigation; this is essentially a pair-level task, which could be undertaken by paying careful attention to the accuracy in the density derivatives involving just the centroid pair direct correlation functions (Equation (28b)).

Although, for low wavenumbers, the isosceles components may be expected to take negative values (as suggested by closure calculations), the current path-integral applications to helium-3 for the centroid and instantaneous structures, recorded at a wavenumber close to those of the maximum amplitudes of the pair structure factors, show that the exact path integral values stay above zero within the complete angular $\varphi -$range $[0,\pi ]$. The symmetrized Denton-Ashcroft and Jackson-Feenberg convolution closures produce results remarkably close to the path integral treatment within the $\varphi -$subinterval $\left[1.16,\pi \right],$ but they fail in its complement $\left[\mathrm{0,1.16}\right],$ where Jackson-Feenberg’s behaves in a wrong way that propagates into Denton-Ashcroft’s. The foregoing facts apply to both classes of isosceles results, centroid and instantaneous, although the discrepancies from the path-integral reference values are especially pronounced in the centroid case (i.e., the closure significant undulations against the path-integral leveled behavior). These results suggest that the study with other more numerically-involved triplet closures (e.g., Barrat-Hansen-Pastore) may be highly significant and also help to clarify how to improve closure approaches.

The quantum fluid triplet topic presents one with many challenges, theoretical and computational. The coming and extended use of increasingly efficient computational resources (i.e., exascale and, hopefully, quantum computers) can thrust the topic of fluid n-body structures forward, especially in the quantum domain. In any event, the research on triplet closures in both the real and the Fourier spaces will remain a valuable complement to the thorough understanding of the quantum fluid structures. It is expected that not only Statistical Mechanics but also disciplines involved in materials design can benefit from all these developments.