Contribution to the Statistical Mechanics of Static Triplet Correlations and Structures in Quantum Fluids

1. Introduction

The study of the correlations between the particles of a system is key to the understanding of its behavior, classical or quantum. Without trying to be exhaustive, one may mention the following main lines in which are involved: (a) equilibrium thermodynamics [1,2,3,4,5,6,7,8,9,10,11,12,13]; (b) time-dependent phenomena [12,13,14,15,16,17,18,19,20,21,22]; (c.1) the experimental [23,24,25,26,27,28,29,30,31] and/or (c.2) computational [7,8,9,10,11,12,13,17,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46] techniques for their practical determinations; and (d) the fundamental questions related to entanglement and separability [47,48,49,50,51,52,53]. Each of these lines forms a whole body of knowledge with direct/indirect practical applications, ranging from phase equilibria and stability to quantum networks and information. Furthermore, all of them taken together display useful and revealing intersection areas. In this connection, the present article will focus on statistical mechanics issues at the intersection between lines (a) and (c.2), by selecting the somewhat forgotten topic of the triplet correlations in equilibrium quantum monatomic fluids (in 3-D space) at nonzero temperatures. This situation contrasts sharply with that of the classical fluid, for which insightful theoretical and computational developments have been available in the literature for a long time (see for instance References [37,38,41,42,43,44,45]).

As a convention, the concept “monatomic fluids” can imply either actual fluids composed of atoms or model fluids composed of structureless particles (one-site). In this work, both terms atom and particle will be used indistinctly. Also, it is convenient to stress that in a broad sense the general concept of equilibrium “correlations” encompasses the usual fluid structures $g_n$ in the real r-space and the response functions $S^{(n)}$ in the reciprocal Fourier k-space (n refers to the number of particles involved in their formulations). Both types of static structures are connected essentially through Fourier transforms, which are far from trivial for $n\ge 3$ [3,7,9,11,41]. One also notes that an extended use has assigned the names correlation functions and structure factors to the $g_n$ and the $S^{(n)},$ respectively. This is the convention followed in this article, although depending on the context “correlations” or “structures” will be used occasionally in the global sense.

For homogeneous and isotropic monatomic fluids, at the pair level $(n=2)$ one deals with 2-D structural functions, whilst at the triplet level$(n=3)$ one deals with 4-D structural functions. Contrarily to the pair functions, the triplet functions cannot be extracted from radiation scattering experiments, because the triplet contribution to the scattering intensity is negligible as compared to the rest of the contributions that shape the differential cross section [24]. Therefore, a comprehensive quantum fluid triplet study, which goes beyond the time-honored closure-theory attempts [3,39,40], entails an appropriately developed theoretical framework, which must be complemented with very powerful computational means. In this way, one can perform the extensive numerical calculations needed (the experiments here). From recent path-integral Monte Carlo plus approximate closure works by the present author [54,55,56,57,58,59,60], one may feel that the time is becoming ripe to start undertaking the triplet topic. As in the classical domain the interest in this quantum task does justify by itself, if only because it makes statistical mechanics go beyond the usual pair level. But for those more application-orientated, one may point out that there are many questions which can benefit from such development (e.g., further fluid structural characterizations of the fluid-solid transition [6,41,61,62], phonon-phonon interactions in superfluids [40], time-dependent phenomena [16,17,18,19,20,21,22,23,33,34] (see Reference [24] to infer the possible connections), colloid suspension properties [63,64], etc.). In relation to this, it may be motivating to mention the empirical relationships found for the maxima of the equilateral correlations in r-space of the quantum hard-sphere (QHS) fluid along the crystallization line [59], and the surprising success of closures in capturing significant quantum fluid triplet features [59,60]. The reader should be aware however that the quantum fluid complexity is far greater than that of its classical fluid counterpart, as will be discussed in detail in this article.

The outline of this article is as follows. Section 2 is devoted to giving an introductory global description of the context in which the quantum triplet topic is inserted, together with its current state of development. Section 3 contains a condensed presentation of the pair and triplet structures in classical statistical mechanics. Why? Because, although there are radical differences between the classical and the quantum domains, both share a good deal of the general notations involved, plus a significant number of basic concepts and tools used (e.g., functional differentiation and closures). This excursion into the classical background is expected to ease the way to the quantum fluid discussion. Section 4 reviews the basic theory of Feynman’s path integrals PI in (thermal) imaginary time [4,7,9]. When complemented with computer simulations (resembling somehow the classical ones), PI can be utilized for the complete study of quantum fluid static triplets (r-space and k-space) in the diffraction and the bosonic exchange regime (more on the exchange issues later). Section 5 focuses upon the higher complexity of the quantum fluid structures as revealed by PI, concentrating on the quantum diffraction regime by reason of its fundamental role. Thus, one finds that the single structural class present in the classical domain, C, composed of pairs ${\left\{g_n,S^{(n)}\right\}}_C,$ splits into three different structural classes, namely instantaneous ETn, ${\left\{g_n,S^{(n)}\right\}}_{ET},$ total continuous linear response TLRn, ${\left\{G_n,S^{(n)}\right\}}_{TLR},$ and centroids CMn, ${\left\{g_n,S^{(n)}\right\}}_{CM}$ [7,9,11]. The discussion focuses on $n=3$, although paying due attention to $n=2$, by reason of its importance when studying triplets. The parallel extension to zero-spin particles and bosonic fluids is deferred to an Appendix. To illustrate some of the main points discussed, a set of various results are included in this Section 5. The studied systems are helium-3 under supercritical conditions, and the quantum hard-sphere fluid on its crystallization line. For helium-3 the results are related to the r- and k-spaces, while for quantum hard spheres are on k-space. Section 6 gives a summary description of the systems, and the methods used to facilitate the understanding of the reported results. The methods are path integral Monte Carlo simulations (PIMC) and a significant number of closures, such as Kirkwood superposition (KS3) [37], Jackson-Feenberg (JF3) [40], the intermediate AV3 $=$(KS3+JF3)/2 [59], Barrat-Hansen-Pastore (BHP3) [41], and Denton-Ashcroft [65]. Section 7 contains a series of final remarks together with directions for future research. Finally, two Appendices, I the boson fluid and II a list of acronyms and basic associated references for the reader’s convenience, close this work.

2. Setting Quantum Fluid Triplets in Their Context

Before entering the mathematical formulations, it  seems useful to give an introductory guide for situating the triplet topic in  its statistical mechanics context. This will also serve the purpose of  introducing some basic concepts and notations. Unrestricted reference to the  classical domain is made in what follows. By doing so, the fundamental quantum  departures from the common ground are expected to be fully grasped by readers  not very familiar with the structural issues.

2.1. The Tackling of the Classical and the Quantum Domains

The use of computers has given a definitive impulse to every aspect of modern research activities. As regards condensed matter physics, the related computer applications have been guided by the numerical implementations of two central lines of thought: approximate theories about the fluid behavior (e.g., see References [2,3,4,5,6]), and the powerful simulation techniques which started from the works by Metropolis et al. (Monte Carlo MC) [66] and Alder and Wainwright (molecular dynamics MD) [67]. In fact, all these methodologies have made it possible to circumvent both theoretical and experimental difficulties that impeded further developments in statistical mechanics. In this connection, suffice it to recall two examples directly relevant to the triplet question: (a) the practical solution of the problem posed by the hierarchical nature of the n-particle correlation functions $g_n$ in r-space (i.e., the exact Bogolyubov-Born-Green-Kirkwood-Yvon BBGKY hierarchy [2,3,4,5,6]), enabling one to compute via simulations the $g_n$ independently of one another [13,42]; and (b) the possibilities of devising computational alternatives for remedying the lack of experimental techniques [24] to determine structure factors $S^{(n)}$in k-space (and thus $g_n$ correlation functions in r-space) beyond the pair level $n=2$.

The approximate theories in their most useful forms utilize closures, which allow one to tackle the calculations of the $g_n$ and the $S^{(n)}$ functions. These two types of objects are in general connected through nontrivial Fourier transforms [3,7,8,9,11,13,24,41]. Closures in r-space hypothesize relationships between a correlation function ${\mathrm{o}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{r}n,g}_n,$ and the set of lower-order correlation functions (e.g., KS3 [37], JF3 [40], or other choices [26,27,38]). To compute the $S^{(n)},$ which are far more complicated objects, additional correlation functions in r-space are introduced. The latter are the so-called direct correlation functions $c_n,$ which are based on Ornstein-Zernike OZn concepts [1,68] (e.g., for $c_2:$ Percus-Yevick [5,69,70,71], hypernetted chain [6,69,70,71], or Baxter’s partition [72]; and for $c_3:$ Barrat-Hansen-Pastore [41], or Denton-Ashcroft [65]). Nowdays the relative importance of closures has obviously decreased in favor of simulation, but one should not adopt a cavalier attitude towards them. Indeed, these closure-theories may yield, at a much lower cost, excellent results as compared to simulation or experiments for key pair quantities (e.g., classical and quantum k-space pair structure factors $S^{(2)})$ [13,46,72] and, also, insightful pictures of the complex structural problems at the triplet level [37,38,39,40,41,58,59,60,65].

The two basic simulation techniques, MC and MD, exploit the properties of the statistical mechanics ensembles using a finite sample to model the real fluid behavior (e.g., in the canonical ensemble $\left(N_S,{\mathcal{V}}_S,T\right),N_S=$number of particles, ${\mathcal{V}}_S=$volume of the simulation box, $T=$ temperature) [2,5,6,12,13,17,71]. These applications are consistent with the concept of thermodynamic limit, $T-lim$: for $\left\{N_S,{\mathcal{V}}_S\right\}\to \infty$, such that ${\displaystyle \frac{N_S}{{\mathcal{V}}_S}}=finite\ne 0$, the intensive quantities are independent of the system size [5] (surface effects from a container are not considered). MC is an elaborate stochastic Markov process that generates particle configurations compatible with the equilibrium state of the fluid being modelled. MD is based on the solving of the equations of motion (e.g., Newton’s) for the particles in the modelled fluid, either in equilibrium or in nonequilibrium states. MC provides “exact” results, since it is essentially a multi-dimensional numerical integration method, and it can be extended to deal with quantum many-body systems at non-zero temperatures [9,13,16]. However, the MD reliability is different depending on the classical or the quantum nature of the application. Thus, although for the time-evolution in classical fluids MD provides a perfectly defined calculational framework [13,17], the same does not occur for quantum fluids because of Heisenberg’s uncertainty principle. Hence, the related quantum dynamics applications are of limited validity, this being a general issue that remains an open problem. (See the appealing formulations and results, based on the PI formalism [4], reported in References [18,19,21,22,33,34]). It must be stated that within the quantum domain PI is the theoretical framework that has led to the implementation of the most powerful simulation techniques to date for investigating equilibrium states. The two basic techniques are path-integral Monte Carlo (PIMC) and path-integral molecular dynamics (PIMD), the latter being an equilibrium construction that bears no connection with the actual time evolution of the system [7,9,10,11,35,36,46,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97]. All these developments of Feynman’s PI owe much to the pioneering works by Barker [73], Chandler and Wolynes [7], Berne et al. [16,74,81], and Ceperley and Pollock [9,77,78]. It is worthwhile to mention that, in addition to the diffraction and the bosonic exchange regimes [9,35,36,78], a further PI extension involving Wigner’s formulation (WPIMC) has been used recently by Filinov et al. to study in a practical way fermionic exchange [98,99], an issue that remained far from accurate PIMC computations owing to the so-called “sign problem” [98,100,101].

2.2. Hamiltonians and Interparticle Potentials

Any of the foregoing computer applications must incorporate as a key ingredient a model for describing the interactions between the actual particles composing the fluid. These models are built by employing any (or all) of the following methods [13,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118]: geometrical arguments (e.g., hard spheres), experimental data, and quantum mechanical reasoning/computations. In this regard, it is well-known that the Born-Oppenheimer approximation (BOA) [119], aided again by the availability of increasing computational resources, has been instrumental in developing the quantum study of real systems. Among its many applications, those of interest in the structural study of fluids may be listed as follows: interatomic/intermolecular interactions [102,103,104,105,106,107,108,109], molecular properties in condensed phases [109,110,111,112], properties of atomic/molecular clusters or far larger assemblies, etc. [113,114,115,116,117,118]. Note that for an isolated polyatomic system BOA defines its electronic energy (including obviously the nuclear repulsions) as the potential energy for the motions of the nuclei within the overall atomic structure. Consequently, for a set of N atoms, from the consideration of a number of geometries of the nuclei at fixed positions ${\mathbf{q}}_j$ in 3-D space, operational forms of the global potential energy $V^{(N)}$ can be determined quantum mechanically as sums of n-body contributions (two-body, three-body, etc.) that depend only on the positions ${\mathbf{q}}_j.$ The point is that a quantum average over the electronic coordinates is to be performed, and this leads to a model of the polyatomic system in which the effective centers for the interactions are placed at the positions of the nuclei, i.e., the non-collapsing and tempered function $V^{\left(N\right)}({\mathbf{q}}_1,{\mathbf{q}}_2,\dots ,{\mathbf{q}}_N)$ [5].

Therefore, for an isolated monatomic fluid one can write the model Hamiltonian in the conveniently simplified form:

$$H_0^{(N)}=K_{kin}^{(N)}+V^{(N)}={\sum }_{j=1}^N{\displaystyle \frac{{\mathit{p}}_j^2}{2m}}+V^{\left(N\right)}({\mathbf{q}}_1,{\mathbf{q}}_2,\dots ,{\mathbf{q}}_N),$$

(1)

where ${\mathbf{p}}_j$ is the momentum of atom j and m is the atom mass. (The ${\mathbf{p}}_j$ and ${\mathbf{q}}_j$ are to be interpreted as classical dynamical functions or as quantum operators, depending on the application; explicit spin interactions are disregarded). Among the $V^{(N)}$ constructions, the success of the pairwise approach, $V^{\left(N\right)}={\sum }_{j<l}v\left(r_{jl}\right),$ has been fundamental. The latter form usually contains “averaged” many-body effects in the “pair” potential $v\left(r_{jl}\right),$ $r_{jl}=\left|{\mathbf{q}}_j-{\mathbf{q}}_l\right|,$ [13] and has yielded a thorough understanding of fluid systems. Particularly, this is true of monatomic fluids at equilibrium, regardless of their classical or quantum behavior [9,13,17]. One also notes that the use of explicit many-body sums for $V^{\left(N\right)}$ becomes necessary when the relative orientations of three [102,103,108] or more atoms (or particles in general) play an important role, as for example at very high densities in the solid state [115,116,117], or when the fixing of the n-body terms is strictly individualized with the use of quantum chemical computations [108,114].

The final test of an interaction model is the comparison with experiment of the calculated properties (e.g., for fluids: the pair structures $S^{\left(2\right)}(k)$ and $g_2(r),$ the internal energies E, the pressures p, etc.). In this regard, one must realize that the form Equation (1) for the Hamiltonian is perfectly suited for comparisons of the related fluid pair structures, fixed via computations, with those determined through radiation scattering experiments (X-ray, neutron diffraction) [5,6,23,28,29,30]. In turn, from these comparisons one can improve the predictive power of the interaction model via further fittings of its parameters [13].

Also, it is important to point out that although the fluid particle correlations are induced by the particle interactions (including those coming from quantum symmetries), both concepts are generally not the same thing. Thus, the many-body correlations may not be specifically associated with an explicit inclusion in the calculations of three-body or higher-order potential energy contributions. Note, for example, that in a monatomic gas described with just a pair potential $v\left(r_{jl}\right)$, one can identify correlations of any order between the atoms. Within such a context, a key fact that may guarantee the high quality of the computed triplet structures is associated with the fluid description achieved via $v\left(r_{jl}\right),$ albeit this might not be sufficient and evaluations with alternative interaction models, or calculations of special fluid properties involving these structures [3,40,44], could be necessary to settle the question.

2.3. Classical Calculations

As mentioned earlier, the classical and PI-quantum applications share several basic features and problems (thermodynamic evaluations can be carried out alongside those of the equilibrium structures and/or utilize the latter [9,13]). Therefore, consideration is given first to the classical case, and the specific quantum characteristics will be discussed afterwards.

The equilibrium classical developments have focused mainly on the pair and the triplet levels, the two basic sets being: ${\left\{g_2\left(r\right);S^{\left(2\right)}(k)\right\}}_C,$ ${\left\{g_3\left(r,s,u\right);S^{\left(3\right)}(k_1,k_2,\cos \varphi )\right\}}_C$ [13,24,41,42,43,44,45,120,121,122,123,124]. The distances r, s, and u define a triplet configuration of three generic atoms: $r=\left|{\mathbf{q}}_1-{\mathbf{q}}_2\right|,$ $s=\left|{\mathbf{q}}_1-{\mathbf{q}}_3\right|,$ $u=\left|{\mathbf{q}}_2-{\mathbf{q}}_3\right|.$The k variables of the response Fourier functions are the wavenumbers $k=\left|\mathbf{k}\right|$ of the related (elastic) momentum transfers $\mathbf{\hslash }\mathbf{k}$ from an external weak field to the fluid $(\varphi$is the angle between ${\mathbf{k}}_1$ and ${\mathbf{k}}_2$ and can be used in lieu of $\cos \varphi$ in$S^{\left(3\right)})$ [5,6,41]. Recall that in going from $n=2$ to $n=3$ the expected computational load must increase greatly, since the dimensionality of the problem goes from 2-D to 4-D.

The classical pair level can be solved via simulation to a high degree of accuracy. However, the $g_2$ calculations are more affordable and straightforward than those of $S^{\left(2\right)},$ which encounter some difficulties. The latter are brought about by the necessity of: (a) scanning wavevectors k commensurate with the simulation box, i.e., for the usual cubic box $\mathbf{k}=\left({\displaystyle \frac{2\mathsf{\pi }}{\mathrm{L}}}\right)\left(k_x,k_y,k_z\right),$ with $k_x,k_y,k_z$ integers and L the box length; and (b) the use of equal-modulus sets of vectors $\left\{{\mathbf{k}}_m\right\}$ for a given value $k=\left|{\mathbf{k}}_m\right|$. Accordingly, the $S^{\left(2\right)}-$ simulation task is far more expensive than that of $g_2$ and, more critically, does not permit the direct calculation of the $\mathbf{k}=0$ component, which is a central thermodynamic quantity: the dimensionless isothermal compressibility, $S^{\left(2\right)}\left(k=0\right)={\rho }_N{k}_{B}T{\chi }_T,$ where ${\rho }_N$ is the bulk (number) density and ${\chi }_T=-{\displaystyle \frac{1}{\mathcal{V}}}{\left({\displaystyle \frac{\partial \mathcal{V}}{\partial p}}\right)}_T$ [13]. These features are well known by statistical mechanics practitioners and, instead of employing extended simulation schemes looking for reliable (and unnecessarily expensive) $k\to 0$ extrapolations$,$ one can utilize methods based on the direct MC/MD simulation of $g_2$ complemented with well-grounded closure procedures [13,71,72,125,126,127]. The latter involve the OZ2 framework, which yields the pair direct correlation function $c_2\left(r\right),$ a function that shows a very rapid decay to zero with increasing distances. This function is decisive in circumventing the difficulties mentioned above and provides, through its Fourier transform, a highly accurate and cost-effective answer to the whole pair structure factor problem: that is $S^{\left(2\right)}(k)$ for $k\ge 0.$

As regards triplets, the whole situation is not so well developed. The simulation techniques are certainly possible in this context [41,42,44,45,120,121,122,123,124], albeit the abovementioned drawbacks at the pair level are clearly augmented. To diminish the cost of full 4-D computations in both r-space and k-space, the triplet studies have concentrated mainly on the considered as the most characteristic features, namely the statistical properties of the equilateral and isosceles features. The descriptions so obtained are not complete, but one can expect them to be sufficiently informative for many practical purposes [41,120,121,122,123,124]. However, if more detailed triplet results are needed, they can be fixed for $g_3$ in r-space through affordable simulations [44], albeit for $S^{(3)}$ in k-space the computational load increases largely. At this juncture one might think of using a combination of simulation and closures [41,120,121,122,123,124], as achieved at the pair level, that is, using the simulated $g_3$ (and $g_2$!) to fix $S^{(3)}.$ Nevertheless, the computations of $g_3(r,s,u)$ would have to be quite complete over the three spatial dimensions for the Fourier transforms to be obtained accurately [120,121]. In addition, the alternative based on OZ3 direct correlation functions, i.e., Baxter’s $c_n$ hierarchy [68], remains in general uncertain due to the lack of a univocal formulation that cuts the hierarchy in a generally significant way [41,65,68,128,129]. Despite these difficulties, the use of closures is a very useful source of information (e.g., for $g_3:$ Kirkwood superposition KS3 [37], Jackson-Feenberg JF3 [40], the intermediate AV3 $=$(KS3+JF3)/2 [59], and/or for $c_3:$ Barrat et al. BHP3 [41], Denton-Ashcroft DAS3 [65]).

For the current purposes, it seems natural to characterize the structural complexity associated with a classical monatomic fluid by the class Cn = ${\left\{g_n;S^{\left(n\right)}\right\}}_C.$ Such class might be extended to include $c_n,$ thus separating the two ways to obtain $S^{\left(n\right)}$ (simulation and direct correlation functions), but given the direct relationships between the $c_n$ and the $S^{\left(n\right)}$ that measure does not seem necessary. By doing so, the comparison with the quantum case is expected to be clearer by avoiding unnecessary complications to the essence of the related discussion. Moreover, direct correlation functions in the quantum case do not have the same general status as in the classical domain, which makes their absence from Cn a reasonable choice.

2.4. Quantum Calculations: Path Integrals

The consideration of the structural study of quantum monatomic fluids at equilibrium is in order now. The classical approach is an idealization in which thermal quantum effects can be neglected at sufficiently high temperature T. However, as T is lowered, the quantum diffraction regime becomes non-negligible, its importance being greater with diminishing temperatures and increasing densities. Furthermore, for sufficiently low T quantum exchange regimes (Bose-Einstein or Fermi-Dirac spin statistics) show up and must be taken into account [9,36,98,99]. All the classical simulation difficulties are also present in this context, with the added intricacies arising from the variety of quantum effects. In this regard, diffraction effects and zero-spin bosonic statistics can be tackled, broadly speaking, along lines of thought that are somewhat parallel. These are the subjects of the present article, whose interest is clear if one considers that every fluid can undergo diffraction effects, and that helium-4 as being composed of zero-spin particles is a paradigm for boson superfluids. Nevertheless, the treatment of quantum fluids with non-zero spin statistics is far more complex, owing to the roles that can be played by: (a) the multiplicity of spin states, (b) the intrinsic half-integer spin statistics difficulties (fermionic fluids), or (c) the magnetic interactions with external fields. All these more general cases will not be addressed in this work, which is focused on the current state of equilibrium triplet structure studies.

To begin with, the quantum monatomic fluid under diffraction effects displays three structural classes, which can be seen as arising from the particle delocalizations that produce a splitting of the classical Cn. These new classes are associated with three different forms of linear response from the quantum fluid to an external weak field. This means that if the structural complexity of the classical monatomic fluid is taken as a reference, the complexity of this quantum counterpart is (at least) three times higher. These quantum classes are revealed by the PI formalism [7,9,18,19,55,84,130], and are known as: instantaneous (ETn), total continuous linear response (TLRn), and centroids (CMn). The general scheme for any of them is just the same as before, correlation functions and structure factors, but with TLRn involving self-particle correlations; such a fact makes this latter class a distinct one. For notational convenience [11] these quantum classes will be denoted as: $\mathrm{E}\mathrm{T}\mathrm{n}=$ ${\left\{g_{ETn};S_{ET}^{(n)}\right\}}_Q,$ $\mathrm{T}\mathrm{L}\mathrm{R}\mathrm{n}={\left\{G_{TLRn};S_{TLR}^{\left(n\right)}\right\}}_Q,$and $\mathrm{C}\mathrm{M}\mathrm{n}={\left\{g_{CMn};S_{CM}^{\left(n\right)}\right\}}_Q$. Their connections with external weak fields are as follows: (a) as ETn is a close “analogous” to Cn, then the corresponding fields are those of elastic neutron scattering or X-ray diffraction [23,28]; (b) TLRn arises from the action of continuous fields [11]; and (c) CMn is a particular case of the latter and corresponds to a field of constant strength [19,130], also showing a truly deep classical-like character [55,130]. On the one hand, the first two classes can be measured experimentally at the pair level $n=2$, that is, ET2 and TLR2; from a wider perspective, these measurements also show connections with dynamic quantities via sum rules [23,29,30]. On the other hand, CM2 cannot be fixed directly through experiments, although it is very useful as an auxiliary object for theoretical developments [18,19,33,34,130,131,132,133,134,135,136,137]. (Beyond $n=2$ recall that no structural measurement is experimentally obtainable today).

To get a first feeling of the PI-origin of these classes, a simple image will suffice [4]. Under diffraction conditions at temperature T every actual particle with mass m delocalizes in the form of a thermal packet, its size being related to the thermal de Broglie wavelength ${\lambda }_B={\displaystyle \frac{h}{\sqrt{2\pi m{k}_{B}T}}}$. Accordingly, there are many possibilities for defining interparticle distances, but within the current response context three sorts of physically significant PI interparticle distances arise [7,9,11,84,130], each associated with one of the three PI classes. Although their exact nature will be discussed in Section 5 of this article, just for visualization purposes notice that the CMn class is built on the distances between the thermal packet centers, which are termed centroids or “centers of mass”. The other two classes involve more elaborate forms of distances between the particle thermal packets. Although the higher complexity of the structural quantum case should be clear at this point, there is still more to this issue that deserves consideration, as the following PI calculational facts reveal.

The quantum PI simulations in the diffraction-effects regime (PIMC, PIMD) describe the delocalized atoms in real 3-D space through a statistical distribution of (closed) elastic “necklaces” with P “beads” apiece (i.e., the thermal packets: one P-necklace per actual atom) [4,7,9]. $P$ is a positive integer, greater than 1, which is to be optimized and represents a compromise between statistical convergence and theoretical accuracy. The latter is given by the Trotter’s limit $P\to \infty$ [138], in which necklaces transform into continuous closed paths in Euclidean (imaginary) time $\beta \hslash =\hslash /k_{B}T,$ and an exact representation of the actual atoms and the fluid is retrieved [9]. (There is an alternate and equivalent PI formulation based on the Fourier expansions of the paths [139,140], but it will not be analyzed in this work). An actual fluid composed of N atoms is thus modelled by a set of $N\times P$ beads. Therefore, P plays a decisive role in the simulation studies of the quantum structures (of course, using a reduced sample $N_S\times P):$ it makes them far more expensive than those of the classical case where $P=1.$ Furthermore, apart from such an increase in the computational effort just to generate the equilibrium configurations, and of the usual N-dependences of any of the n-th order structure calculations [2,13], to obtain the structures it is crucial to point out that: the instantaneous (ETn) calculations scale as P, the total continuous linear response (TLRn) do as $P^n,$ and the centroid (CMn) do as in the classical case.

As stressed, when diffraction effects dominate the fluid behavior, and thus exchange interactions can be neglected, the PI necklaces are closed. However, if quantum exchange cannot be neglected, one faces complicated mixtures of closed and open necklaces, the latter interlinking with one another, in response to the corresponding permutations among atoms [7,9,35,36,78]. Within the PI-description, the physically significant interparticle distances are defined by specific inter-bead distances (ETn and TLRn), or by constructions that use the bead positions (CMn), as will be shown in detail in Section 5 (diffraction effects) and Appendix I (zero-spin bosonic exchange) [9,35,130,132,133,134,135].

Therefore, the quantum structural complexity extremely exceeds in practice that of the classical case. Even if one only wishes to perform a full structural study including up to the triplet level, the PI computational effort appears today as a daunting task. As an aside, one may mention the approximate frameworks arising from the exact PI formalism (i.e., Feynman-Hibbs potentials [4,140] and their generalizations [141,142]), which can be utilized to carry out thermodynamic and structural studies of fluids with weak quantum behavior [11,131,143,144,145]. Although their applicability is certainly limited, they are far from being computationally expensive, and turn out to be very useful due to the intuitive structural pictures they provide at both the pair and the triplet levels [11].

2.5. Quantum Calculations: Closures

As regards closure-theories, their applications to quantum fluids may be carried out by following/adapting the same methodologies as those employed for classical fluids in both the r-space and the k-space. Actually, these adaptations have gone both ways (classical ↔ quantum), which has been a normal practice throughout the years [3,37,38,39,40,41,65,68,69,70,71,72,125,126,127,128]. However, the related theoretical quantum grounds deserve especial consideration, which becomes clear when applying OZn schemes.

Thus, by focusing on diffraction effects, it is very interesting that OZn turns out to be a rigorous framework for the centroid class CMn [55], albeit OZn is only approximate for ETn and TLRn [11]. Surprisingly, despite such mathematical observation, excellent results, as compared to experiment and simulation, have been obtained when fixing the three pair structure factors, $S_{ET}^{(2)}(k)$, $S_{TLR}^{\left(2\right)}\left(k\right),$ and $S_{CM}^{\left(2\right)}\left(k\right),$ by OZ2-treating the corresponding PI-simulated correlation functions [11,46,55,60,62,131,136,144,145]. The OZ2 calculations are much less demanding than those carried out with PIMC in k-space, and the results for $k>0$ obtained with PIMC and OZ2 are almost indistinguishable from one another [46]. Moreover, the OZ2 evaluations of the $\mathbf{k}=0$ components are obtained with great accuracy as part of the whole $S^{\left(2\right)}(k)$ determinations. Furthermore, according to the extended compressibility theorem [130], there is the exact formulation of the isothermal compressibility ${\chi }_T$ of the fluid that can be achieved with OZ2 through the centroid component $S_{CM}^{\left(2\right)}\left(k=0\right)={\rho }_N{k}_{B}T{\chi }_T$. This implies that the OZ2 applications to the CM2 correlations $g_{CM2}(R)$ produce, in principle, exact values for ${\chi }_T$ (within the accuracy and precision of the calculations performed) [130,136]. For the cases ET2 and TLR2, which in rigor must share that same centroid $k=0$ value [130], one can obtain approximate, though very valuable, ${\chi }_T$ estimates via OZ2 [46,60,62]. As regards triplet closures, the question is open and must be settled through comparison with exact PIMC/PIMD results. Although uncertainties larger than in the classical domain might be expected here, some promising results have been obtained lately [58,59,60] as summarized below. In relation to bosonic exchange, even though closures were employed in the past to deal with zero-spin/charged fluids [3,39,40,69], their use is not so straightforward, and some consideration will be given in Appendix I.

2.6. Some Recent Quantum Triplet Facts

To close this overview, three motivating pieces of information on quantum fluid triplets in the diffraction regime, obtained recently with PIMC and closures, may be worth giving here.

(i) For the quantum hard-sphere fluid on the crystallization line [59], the PIMC absolute maximum amplitudes of the instantaneous (ET) and centroid (CM) pair and equilateral structures, i.e., $g_2(r_M)$ and $g_3\left(r_M,r_M,r_M\right),$ do follow empirical behaviors of the type $a{\gamma }^{-b},$ where $\gamma ={\rho }_N{\lambda }_B^3$ is the degeneracy parameter, and $a$ and b are positive constants.

(ii) A preliminary PIMC triplet study of liquid para-hydrogen in k-space [58] suggests that on the quantum crystallization line there might exist an almost “constancy” (i.e., values in a relatively narrow interval) of a salient trait in the centroid triplet structure factor (e.g., the maximum of the equilateral components $S_{CM}^{\left(3\right)}(k_M,k_M,{\displaystyle \frac{\pi }{3}})$ or an isosceles amplitude located nearby). If confirmed, this could be a parallel to the quantum freezing constancy observed at the pair level of the centroid maximum amplitude $S_{CM}^{\left(2\right)}(k_M)$ for a given fluid [62]. (At the pair level there is no general quantum extension of the classical Hansen-Verlet rule [6] for the centroid maximum amplitudes $S_{CM}^{\left(2\right)}(k_M)$).

(iii) The applicability of closures to quantum triplets could have a broader scope than suspected and be very helpful to provide physically significant pictures of the different quantum triplet behaviors in cost-effective ways. Even though the available information on closure applications to quantum triplets is scarce [54,55,56,57,58,59,60], for a wide range of quantum diffraction conditions (quantum hard spheres [59] and supercritical helium-3 [60], the intermediate triplet closure AV3 = (KS3+JF3)/2 in r-space captures significant features of the correlation functions $g_{ET3}\left(r,s,s\right)$ and $g_{CM3}\left(r,s,s\right)$. Consequently, improvements on the AV3 closure, e.g., via adaptations of Abe’s expansion [38], should be well-worth exploring.

3. Equilibrium Classical Fluid Triplets

A homogeneous and isotropic classical monatomic  fluid is considered, which will serve as a reference for basic notations,  concepts, and tools, when presenting the quantum fluid discussion. A summary  description of both the pair and the triplet levels is given because a complete  study of the triplet level $(n=3)$ already involves the consideration of the pair  level $\left(n=2\right)$ [37,38,39,40,41]. The  grand canonical ensemble $\left(\mu ,\mathcal{V},T\right)$is employed in all the derivations, with the mean  number density being ${\rho }_N=〈N〉/\mathcal{V}$, and the interatomic distances denoted by $r_{jl}=\left|{\mathbf{q}}_j-{\mathbf{q}}_l\right|.$ The partition function reads as [2]:

$${\mathsf{\Xi }}_C={\sum }_{N\ge 0}\mathrm{e}\mathrm{x}\mathrm{p}\left(\beta \mu N\right)Z_{N,C}(N,\mathcal{V},T)={\sum }_{N\ge 0}{\displaystyle \frac{z^N}{N!}}\int d{\mathbf{r}}^N\times \mathrm{e}\mathrm{x}\mathrm{p}\left[-\beta V^{(N)}({\mathbf{r}}^N)\right],$$

(2)

where $\mu$ is the chemical potential, $\beta$ the inverse temperature ${\displaystyle \frac{1}{k_{B}T}},Z_{N,C}$ the partition function of the canonical ensemble $\left(N,\mathcal{V},T\right),$ z the activity ${\lambda }_B^{-3}\mathrm{e}\mathrm{x}\mathrm{p}\left(\beta \mu \right)$, and $d{\mathbf{r}}^N=d{\mathbf{r}}_{1}d{\mathbf{r}}_2\dots d{\mathbf{r}}_N$ is the 3N-dimensional element of volume for  N atoms. Throughout this article, the ensemble averages are denoted by $〈\dots 〉,$ which in their definitions contain the division by  the corresponding partition function.

3.1. r-Space Correlation Functions

The one-point density dynamical function is simply  given by ${\sum }_j\delta \left({\mathbf{r}}_j-{\mathbf{q}}_1\right),$its ensemble average being ${\rho }_N^{\left(1\right)}\left({\mathbf{q}}_1\right)={\rho }_N.$ The pair $g_2$ and triplet $g_3$ correlation functions are probability density  functions (not normalized to unity) derived from the partition function ${\mathsf{\Xi }}_C$ [2]. These  functions can be defined with the help of Dirac’s $\delta$ as [5]:

$${{\rho }_N^{\left(2\right)}\left({\mathbf{q}}_1,{\mathbf{q}}_2\right)={\rho }_N^{2}g}_2\left({\mathbf{q}}_1,{\mathbf{q}}_2\right)={\rho }_N^2{g}_2\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)〉=$$

(3a)

$${\displaystyle \frac{1}{{\mathsf{\Xi }}_C}}{\sum }_{N\ge 0}{\displaystyle \frac{z^N}{N!}}\int d{\mathbf{r}}^N\times \mathrm{e}\mathrm{x}\mathrm{p}\left[-\beta V^{(N)}({\mathbf{r}}^N)\right]\times {\sum }_{j\ne l}\delta \left({\mathbf{r}}_j-{\mathbf{q}}_1\right)\delta \left({\mathbf{r}}_l-{\mathbf{q}}_2\right),$$

(3b)

$${{\rho }_N^{\left(3\right)}\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3\right)={\rho }_N^{3}g}_3\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3\right)={\rho }_N^3{g}_3\left(r_{12},r_{13},r_{23}\right)=$$

(4a)

$$〈{\sum }_{j\ne l\ne m\ne j}\delta \left({\mathbf{r}}_j-{\mathbf{q}}_1\right)\delta \left({\mathbf{r}}_l-{\mathbf{q}}_2\right)\delta \left({\mathbf{r}}_m-{\mathbf{q}}_3\right)〉=$$

(4b)

$${\displaystyle \frac{1}{{\mathsf{\Xi }}_C}}{\sum }_{N\ge 0}{\displaystyle \frac{z^N}{N!}}\int d{\mathbf{r}}^N\times \mathrm{e}\mathrm{x}\mathrm{p}\left[-\beta V^{(N)}({\mathbf{r}}^N)\right]\times {\sum }_{j\ne l\ne m\ne j}\delta \left({\mathbf{r}}_j-{\mathbf{q}}_1\right)\delta \left({\mathbf{r}}_l-{\mathbf{q}}_2\right)\delta \left({\mathbf{r}}_m-{\mathbf{q}}_3\right).$$

(4c)

Note that $g_2$ and $g_3$ are intensive quantities which depend on ${\rho }_N$ and T and, essentially, they are given by  the ensemble averages of the two- and three- point density dynamical functions  (notice the avoidance of self-contributions). Their normalizations read as [2]:

$${\rho }_N^2\int d{\mathbf{q}}_{1}d{\mathbf{q}}_2{g}_2\left(r_{12}\right)=〈N(N-1)〉,$$

(5)

$${\rho }_N^3\int d{\mathbf{q}}_{1}d{\mathbf{q}}_{2}d{\mathbf{q}}_3{g}_3\left(r_{12},r_{13,}{r}_{23}\right)=〈N(N-1)(N-2)〉.$$

(6)

A powerful alternate way to formulate the foregoing  correlation functions is based on the use of the functional differentiation of $\ln {\mathsf{\Xi }}_C$ with respect to the variations $\delta \mathsf{\Psi }$ in an external weak field $\mathsf{\Psi },$ which acts on every particle as $\mathsf{\Psi }\left({\mathbf{r}}^N\right)={\sum }_j\mathsf{\Psi }\left({\mathbf{r}}_j\right)$ [5,6,146,147,148]. The resulting equations are somewhat involved, but they lead to the physical interpretation of these functions as directly related to the response from the fluid to $\mathsf{\Psi }.$ By starting from${\Xi }_C(\mathsf{\Psi })$ given by

$${\mathsf{\Xi }}_C(\mathsf{\Psi })={\sum }_{N\ge 0}\mathrm{e}\mathrm{x}\mathrm{p}\left(\beta \mu N\right)Z_{N,C}(\mathsf{\Psi })={\sum }_{N\ge 0}{\displaystyle \frac{z^N}{N!}}\int d{\mathbf{r}}^N\times \mathrm{e}\mathrm{x}\mathrm{p}\left[-\beta \left(V^{\left(N\right)}({\mathbf{r}}^N)+\mathsf{\Psi }\left({\mathbf{r}}^N\right)\right)\right],$$

(7)

one obtains the following set of equations:

$$-k_{B}T{\displaystyle \frac{\delta \ln {\Xi }_C(\mathrm{}\mathsf{\Psi })}{\delta \mathsf{\Psi }\left({\mathbf{q}}_1\right)}}={\rho }_N^{\left(1\right)}\left({\mathbf{q}}_1;\mathsf{\Psi }\right)={\mathsf{\Gamma }}_1\left({\mathbf{q}}_1;\mathsf{\Psi }\right),$$

(8a)

$${\left(-k_{B}T\right)}^2{\displaystyle \frac{{\delta }^2\ln {\Xi }_C(\mathrm{}\mathsf{\Psi })}{\delta \mathsf{\Psi }\left({\mathbf{q}}_1\right)\delta \mathsf{\Psi }\left({\mathbf{q}}_2\right)}}=-k_{B}T{\displaystyle \frac{\delta {\mathsf{\Gamma }}_1\left({\mathbf{q}}_1;\mathsf{\Psi }\right)}{\delta \mathsf{\Psi }\left({\mathbf{q}}_2\right)}}={\mathsf{\Gamma }}_2\left({\mathbf{q}}_1,{\mathbf{q}}_2;\mathsf{\Psi }\right),$$

(8b)

$${\left(-k_{B}T\right)}^3{\displaystyle \frac{{\delta }^3\ln {\Xi }_C(\mathrm{}\mathsf{\Psi })}{\delta \mathsf{\Psi }\left({\mathbf{q}}_1\right)\delta \mathsf{\Psi }\left({\mathbf{q}}_2\right)\delta \mathsf{\Psi }\left({\mathbf{q}}_3\right)}}=-k_{B}T{\displaystyle \frac{\delta {\mathsf{\Gamma }}_2\left({\mathbf{q}}_1,{\mathbf{q}}_2;\mathsf{\Psi }\right)}{\delta \mathsf{\Psi }\left({\mathbf{q}}_3\right)}}={\mathsf{\Gamma }}_3\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3;\mathsf{\Psi }\right).$$

(8c)

The foregoing development can be easily generalized  to obtain the hierarchy of functionals $\left\{{\mathsf{\Gamma }}_n(\mathsf{\Psi })\right\},$which can be seen as the recursive definition of  the properties${\mathrm{}\mathsf{\Gamma }}_n(\mathsf{\Psi }$) of the inhomogeneous fluid in the presence of the  external field $\mathsf{\Psi }.$ By letting $\mathsf{\Psi }\to 0$, the ${\mathsf{\Gamma }}_n$ functionals/properties on the right-hand sides of  Equations (8a)-(8c) transform into properties of the isolated fluid

$${\mathsf{\Gamma }}_1\left({\mathbf{q}}_1;\mathsf{\Psi }=0\right)={{\rho }_N^{\left(1\right)}\left({\mathbf{q}}_1;\mathsf{\Psi }=0\right)=〈{\sum }_{j=1}^N\delta \left({\mathbf{r}}_j-{\mathbf{q}}_1\right)〉=\rho }_N,$$

(9a)

$${\mathsf{\Gamma }}_2\left({\mathbf{q}}_1,{\mathbf{q}}_2;\mathsf{\Psi }=0\right)=〈\left({\sum }_{j=1}^N\delta \left({\mathbf{r}}_j-{\mathbf{q}}_1\right)-{\rho }_N\right)\left({\sum }_{l=1}^N\delta \left({\mathbf{r}}_l-{\mathbf{q}}_2\right)-{\rho }_N\right)〉,$$

(9b)

$${\mathsf{\Gamma }}_3\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3;\mathsf{\Psi }=0\right)=$$

$$〈\left({\sum }_{j=1}^N\delta \left({\mathbf{r}}_j-{\mathbf{q}}_1\right)-{\rho }_N\right)\left({\sum }_{l=1}^N\delta \left({\mathbf{r}}_l-{\mathbf{q}}_2\right)-{\rho }_N\right)\left({\sum }_{m=1}^N\delta \left({\mathbf{r}}_m-{\mathbf{q}}_3\right)-{\rho }_N\right)〉.$$

(9c)

which, as shown in Equations (9b)-(9c), give the  two- and three- body spontaneous fluctuations of the fluid density. The  explicit formulas for ${\mathsf{\Gamma }}_2$ and ${\mathsf{\Gamma }}_3$ are standard knowledge in classical statistical  mechanics [41,146,147,148], and for brevity are  not written here. It is easy to see, however, that ${\mathsf{\Gamma }}_3$ contains contributions of the ${\mathsf{\Gamma }}_2$ type. Note that the pair correlation function $g_2$ is contained in the formula of ${\mathsf{\Gamma }}_2\left({\mathbf{q}}_1,{\mathbf{q}}_2;\mathsf{\Psi }=0\right),$ and that both correlation functions $g_2$ and $g_3$ are in that of ${\mathsf{\Gamma }}_3\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3;\mathsf{\Psi }=0\right)$. As an aside, note that for the canonical  component $Z_{N,C}(\mathsf{\Psi })$ the following complementary identities hold:

$$\int d{\mathbf{q}}_{1}d{\mathbf{q}}_2{\displaystyle \frac{{\delta }^2\ln Z_{N,C}(\mathrm{}\mathsf{\Psi })}{\delta \mathsf{\Psi }\left({\mathbf{q}}_1\right)\delta \mathsf{\Psi }\left({\mathbf{q}}_2\right)}}\equiv 0,$$

(9d)

$$\int d{\mathbf{q}}_{1}d{\mathbf{q}}_{2}d{\mathbf{q}}_3{\displaystyle \frac{{\delta }^3\ln Z_{N,C}(\mathrm{}\mathsf{\Psi })}{\delta \mathsf{\Psi }\left({\mathbf{q}}_1\right)\delta \mathsf{\Psi }\left({\mathbf{q}}_2\right)\delta \mathsf{\Psi }\left({\mathbf{q}}_3\right)}}\equiv 0,$$

(9e)

which is the expected behavior of this type of  integrals, since there are no fluctuations in the number of particles in a  canonical ensemble. The analogous integrals involving $\ln {\mathsf{\Xi }}_C(\mathrm{}\mathsf{\Psi })$ are nonvanishing and are related to the true  number fluctuations, which are consistently defined in a grand canonical  ensemble. In this connection, given the formal equivalence between the  structural classes classical Cn and PI-centroid CMn [55,130],  the explicit ${\mathsf{\Gamma }}_2$ and ${\mathsf{\Gamma }}_3$ formulas and many other features are deferred to  the treatment of such PI quantum correlations. Hence, only some important facts  and consequences are to be considered in what follows.

3.2. k-Space Structure Factors

By approximating the right-hand sides of Equations  (8) by the zero-field Equations (9) and Fourier transforming, one obtains the  response functions in k-space from the fluid to the action of $\mathsf{\Psi }$ at the pair and the triplet levels [6,41]. These functions are the static structure  factors that read as:

$$S^{(2)}\left(\mathbf{k}\right)=S^{(2)}\left(k\right)={\displaystyle \frac{1}{{\rho }_N}}{\mathsf{\Gamma }}^{(2)}\left(\mathbf{k};\mathsf{\Psi }=0\right)={\displaystyle \frac{1}{〈N〉}}〈{\left|{\sum }_{\mathit{j}=1}^N\mathrm{e}\mathrm{x}\mathrm{p}\left(i\mathbf{k}·{\mathbf{r}}_j\right)\right|}^2〉-{\left(2\pi \right)}^3{\rho }_N\delta \left(\mathbf{k}\right)=$$

(10a)

$$1+{\rho }_N\int d\mathbf{r}\mathbf{}\mathrm{e}\mathrm{x}\mathrm{p}\left(i\mathbf{k}·\mathbf{r}\right)\mathrm{}\left(g_2\left(r\right)-1\right),$$

(10b)

$$S^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)=S^{(3)}\left(k_1,k_2,\cos \left({\mathbf{k}}_1,{\mathbf{k}}_2\right)\right)={\displaystyle \frac{1}{{\rho }_N}}{\mathsf{\Gamma }}^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2;\mathsf{\Psi }=0\right)=$$

$${\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·{\mathbf{r}}_j+{\mathbf{k}}_2·{\mathbf{r}}_l-\left({\mathbf{k}}_1+{\mathbf{k}}_2\right)·{\mathbf{r}}_m\right)\right]〉-\left\{\delta \left(\mathbf{k}\right)-\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{s}\right\}.$$

(11)

In writing Equations (10) and (11) use of the  homogeneity and isotropy of the fluid is made. The response functions $S^{(2)}$ and $S^{(3)}$ are thus related to the density-density and the  density-density-density correlation functions of the isolated fluid,  respectively.

Note that in a broad sense one enters a general  linear response framework when takes $\mathsf{\Psi }=0$ on the right-hand sides of the hierarchy and  carries out the Fourier transforms on both sides of the equations. In this  regime, the fluid “does not know” if what is being undergone is due to an  external field or is just the result of its own spontaneous fluctuations.  Therefore, each of the (variations in the) fluid properties${\mathsf{\Gamma }}^{(n)}\left(\mathsf{\Psi }\right)$ in k-space can be formulated in terms of  its corresponding isolated fluid property ${\mathsf{\Gamma }}_n\left(\mathsf{\Psi }=0\right)$ in r-space, and in such a way that the  functional dependence on the external field is linear [5,6]. In this connection, for clarity, it is worth writing here the well-known relationship at the pair level [6,148]: ${\delta \rho }_N^{\left(1\right)}\left(\mathbf{k};\mathsf{\Psi }\right)=\delta {\mathsf{\Gamma }}^{(1)}\left(\mathbf{k};\mathsf{\Psi }\right)\approx -\beta {\rho }_N{S}^{\left(2\right)}\left(\mathbf{k}\right)\delta \mathsf{\Psi }\left(\mathbf{k}\right),$ and also that at the triplet level: $\delta {\mathsf{\Gamma }}^{(2)}\left({\mathbf{k}}_1,{\mathbf{k}}_2;\mathsf{\Psi }\right)\approx -\beta {\rho }_N{S}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)\delta \mathsf{\Psi }\left({\mathbf{k}}_1+{\mathbf{k}}_2\right).$ The static structure factors, or response  functions, are thus the proportionality factors between the variations in the  properties and the variations in the external field. Accordingly, the  associated wavenumbers k, $k_1$,and $k_2$ are the moduli of the wavevectors that  define the (elastic) momentum transfers ${\mathbf{p}}_{\mathsf{\Psi }}$ from the field to the fluid i.e., ${\mathbf{p}}_{\mathsf{\Psi }}=\mathbf{\hslash }\mathbf{k}$. These characteristics linked to the use of the  linear response regime will be maintained when studying the quantum domain.

Some important points to be noticed are the  following:

(i) $S^{(2)}\left(\mathbf{k}\right)$ can be determined experimentally, while $S^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)$ cannot [24].

(ii) $S^{(2)}\left(\mathbf{k}\right)$ and $S^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)$ are real-valued quantities [41].

(iii) The general form of $S^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)$ expanded in terms of the Fourier transforms of the  different contributions built with $g_2\left(r\right)$ and $g_3\left(r_{12},r_{13},r_{23}\right)$ is quite involved and will be given in connection  with the quantum CMn discussion.

(iv) For its direct theoretical and experimental  interest, note that the $\delta \left(\mathbf{k}\right)-$term in $S^{(2)}\left(\mathbf{k}\right)$ Equation (10a) represents the subtraction  of the forward scattering of radiation (elastic scattering). This operation is  an integral part of the complete definition of this pair structure function  (note that $S^{(2)}\left(k\right)\to 1$ for large wavenumbers). Obviously, a  generalization of this operation also applies to the $S^{(3)}$ case, as will be shown later on.

Some Concomitant Problems and Their OZ Solutions

In this connection a deeper consideration is  necessary. Note that the $\delta \left(\mathbf{k}\right)-$term subtraction in Equation (10a) is also linked  to the theoretical asymptotic value of $g_2\left(r\right),$ i.e., $g_2\left(r\right)\to 1$ for $r\to \infty$ in the grand canonical ensemble [2], which gives full meaning to the Fourier  transform in Equation (10b). Moreover, it is customary in simulation work to  omit this $\delta \left(\mathbf{k}\right)-$term, as the $\mathbf{k}=0$ component is strictly unattainable due to the  finite sample size (i.e., finite $N_S$ and/or ${\mathcal{V}}_S),$ which imposes commensurability with the simulation  box on the wavevectors to be scanned (clearly, this also negatively affects the  calculation of the low-k region). Such omission is acceptable when  dealing with the simulation techniques but neglects the key fact that the $\mathbf{k}=0$ component makes the connection with the  thermodynamics of the fluid [3,6,13], since $S^{(2)}\left(\mathbf{k}=0\right)={\rho }_N{k}_{B}T{\chi }_T,$ as stressed earlier. The previous connection  yields an appealing way to determine the equation of state of the fluid via ${\chi }_T$ integration, an operation that is formally  independent of the many-body form selected to represent the potential energy $V^{(N)}$ (obviously, $g_2(r)$ does need$V^{(N)}$ to be calculated). However, in calculating $S^{(2)}\left(\mathbf{k}\right)\mathit{}$via Equation (10b) one faces the problem of  the long-ranged oscillations of $g_2(r)$ about unity [2],  which forces one to have a knowledge of the pair correlations over a  sufficiently long range of interparticle distances Otherwise, accuracy in the  Fourier transform is doomed to fail [13]  (e.g., unphysical estimates for ${\chi }_T$ can result from that naïve use of the transform).  In addition, as the optimal determination of $g_2(r)$ is accomplished via simulation, this already poses  a practical problem derived from the need for more expensive computations by  increasing the simulation sample sizes [13].  Furthermore, although the canonical ensemble is a very convenient  calculational tool, the necessarily finite $N_S$ combined with the intrinsically bad asymptotic  behavior of $g_2(r)$ in this ensemble, i.e., $g_2\left(r\to \infty \right)\to 1+O\left({\displaystyle \frac{1}{N}}\right)$ [2,3], present  one with another important drawback to be dealt with.

The same discussion, but with augmented  difficulties, applies to the computation of $S^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right),$ in which the $\left\{\delta \left(\mathbf{k}\right)-\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{s}\right\}$ in Equation (11) involve ${\mathbf{k}}_1$ and/or ${\mathbf{k}}_2.$ Under these circumstances, irrespective of  the mighty computational facilities at one’s disposal, targets that appear as  most desirable are: the saving of electrical power, the shortening of the time  to get answers through more rapid calculations, and the possibilities of  obtaining physical pictures of the underlying processes.

The obvious advantage of using functional  techniques is that they give the physically complete answers of the classical  structural problems at one stroke, although, inevitably, some problems stand  out and others arise. Nevertheless, there is more to using functional  techniques in this context, since they can be employed to build insightful  solutions to these drawbacks. In relation to this, the additional role played  by the inverse functional derivatives (e.g., $\delta \mathsf{\Psi }({\mathbf{q}}_2)/\delta {\rho }_N^{\left(1\right)}({\mathbf{q}}_1))$ in the derivation of the Ornstein-Zernike (OZn)  equations is to be highlighted [41,65,68,69,70,71,72,146,147].  The OZn equations define a further type of structural functions: the direct  correlation functions $c_2\left(r\right),$ $c_3\left(r_{12}{,r}_{13},r_{23}\right),$ etc. By restricting up to the triplet level for  simplicity, $c_2$ and$c_3$ yield compact and efficient formulations of the  structure factors $S^{(2)}$ and $S^{(3)}.$ The key point is that the direct correlation  functions are short-ranged (i.e., they decay to zero rapidly with increasing  distances), which allows one to fix very accurately the Fourier transforms  involved.

Given that the OZn framework at the pair and the  triplet levels will be considered in detail when dealing with the isomorphic  PI-centroid CMn class, for brevity, the classical situation may be summarized  here in the following basic equations [1,41,68]:

-

OZ2  equation

$$h_2\left(r_{12}\right)=c_2\left(r_{12}\right)+{\rho }_N\int d{\mathbf{r}}_3{h}_2\left(r_{13}\right)c_2\left(r_{23}\right);\mathrm{}{h}_2\left(r_{12}\right)=g_2\left(r_{12}\right)-1,$$

(12)

where the total correlation function $h_2\left(r_{12}\right)$ is introduced and a convolution integral over the  positions of a generic third particle is included. (Although $c_2(r)$ may show oscillatory tails [149], the observed behavior is that the decay $c_2(r)\to 0$ for increasing distances occurs effectively).

-

OZ2  pair structure factor

$$S^{(2)}\left(\mathbf{k}\right)={\left(1-{\rho }_N{c}_2\left(\mathbf{k}\right)\right)}^{-1},$$

(13a)

$$S^{(2)}\left(\mathrm{k}=0\right)={\left(1-{\rho }_N{c}_2\left(k=0\right)\right)}^{-1}={\rho }_N{k}_{B}T{\chi }_T={\displaystyle \frac{〈N^2〉-{〈N〉}^2}{〈N〉}},$$

(13b)

where $c_2\left(\mathbf{k}\right)$ is the Fourier transform of (the short-ranged) $c_2\left(r_{12}\right).$(Note that even at the critical point $c_2\left(\mathbf{k}=0\right)$ remains finite, its value being ${\rho }_N^{-1}$ [5,6]).

-

Baxter’s  hierarchy related results

$${\displaystyle \frac{\partial c_2\left({\mathbf{q}}_{1,}{\mathbf{q}}_2\right)}{\partial {\rho }_N}}=\int d{\mathbf{q}}_3{c}_3\left({\mathbf{q}}_{1,}{\mathbf{q}}_2{,\mathbf{q}}_3\right);T=\mathrm{constant},$$

(14a)

$${\displaystyle \frac{\partial c_2\left({\mathbf{k}}_1\right)}{\partial {\rho }_N}}=c^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2=0\right);T=\mathrm{constant},$$

(14b)

where $c_2\left({\mathbf{q}}_{1,}{\mathbf{q}}_2\right)=c_2\left(r_{12}\right),$ $c_3\left({\mathbf{q}}_{1,}{\mathbf{q}}_2{,\mathbf{q}}_3\right)=c_3\left(r_{12}{,r}_{13},r_{23}\right),$ and $c^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)=c^{(3)}\left(k_1,k_2,\mathrm{c}\mathrm{o}\mathrm{s}({\mathbf{k}}_1,{\mathbf{k}}_2)\right)$ is the Fourier transform of $c_3\left(r_{12}{,r}_{13},r_{23}\right)$.

-

OZ3  triplet structure factor

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

(15)

Some further comments are in order now. First, the  OZ2 equation takes a single form, whereas one finds four different, though  equivalent, versions of the OZ3 equation [70].  Second, the $\delta \left(\mathbf{k}\right)-$terms appearing in Equations (10a) and (11) are  fully included in their reformulations given by Equations (13) and (15).  Therefore, the $\mathbf{k}=0$ incompleteness associated with simulation  schemes appears formally OZn-circumvented. Third, note the important role of  the pair level in k-space, as shown in the triplet Equations (14)-(15)  where the pair quantities give a solution to the simulation-intractable  situations at ${\mathbf{k}}_1=0$ and/or ${\mathbf{k}}_2=0$ (also ${\mathbf{k}}_1=-{\mathbf{k}}_2)$, and are indispensable building blocks. Fourth,  however, the OZn equations are void in themselves, and they need further  approaches (closures) to be solved. This is easy to grasp by observing  the structure of the Equations (12) and (14): both are integral equations for  the unknown $c_2$ and $c_3$, which add new difficulties to the whole  structural question [24,26,27,37,38,39,40,41,45,65,68,69,70,71,72,121,122,123,124,125,126,127,128,129].  Further formal considerations (triplet symmetries and asymptotics, consistency  relationships, etc.) and analyses to fix the pending issues (e.g., closures,  the vexing canonical $g_2\left(r_{12}\right)$ long-range behavior, etc.) are given in the next  subsections.

3.3. Grand Canonical Triplet Symmetries and Asymptotic Behaviors

(i) In r-space, by renumbering the atom  labels, the triplet $g_3$ and $c_3$ functions satisfy the following symmetries:

$$g_3\left(r_{12},r_{13},r_{23}\right)=g_3\left(r_{12},r_{23},r_{13}\right)=g_3\left(r_{23},r_{12},r_{13}\right)=\dots ,$$

(16a)

$$c_3\left(r_{12},r_{13},r_{23}\right)=c_3\left(r_{12},r_{23},r_{13}\right)=c_3\left(r_{23},r_{12},r_{13}\right)=\dots ,$$

(16b)

or using vector notation (e.g., ${\mathbf{r}}_{12}=\mathbf{u},$ ${\mathbf{r}}_{13}=\mathbf{s})$ [41]:

$$g_3\left(\mathbf{u},\mathbf{s}\right)=g_3\left(\mathbf{s},\mathbf{u}\right)=g_3\left(-\mathbf{u},\mathbf{s}-\mathbf{u}\right),$$

(16c)

$$c_3\left(\mathbf{u},\mathbf{s}\right)=c_3\left(\mathbf{s},\mathbf{u}\right)=c_3\left(-\mathbf{u},\mathbf{s}-\mathbf{u}\right).$$

(16d)

Also, when two particles get very close together  the triplet function $g_3$ obviously vanishes [150,151].  Technically, one writes

$${\lim }_{\left|{\mathbf{q}}_1-{\mathbf{q}}_2\right|\to 0}{g}_3\left(r_{12},r_{13},r_{23}\right)=0,$$

(17a)

although, as at the pair level, the zero value is  reached for interparticle distances greater than zero. The $g_3$ asymptotic behaviors are [2,151]:

$${\lim }_{\left|{\mathbf{q}}_3\right|\to \mathbf{\infty }}{g}_3\left(r_{12},r_{13},r_{23}\right)=g_2\left(r_{12}\right),$$

(17b)

$${\lim }_{\begin{array}{c}\left|{\mathbf{q}}_2\right|,\left|{\mathbf{q}}_3\right|\to \infty \\ \left|{{\mathbf{q}}_2-\mathbf{q}}_3\right|\to \infty \end{array}}{g}_3\left(r_{12},r_{13},r_{23}\right)=1,$$

(17c)

(conventionally, one may take particle 1 as the  origin of coordinates for the three generic particles). For the equilateral and  isosceles correlations these conditions are usually written as

$${\lim }_{s\to \mathbf{\infty }}{g}_3\left(r,s,s\right)=g_2\left(r\right),$$

(17d)

$${\lim }_{r\to \mathbf{\infty }}{g}_3\left(r,r,r\right)=1,$$

(17e)

In the canonical ensemble Equation (17e) and its  analogous at higher orders are not so neat because of terms $O({\displaystyle \frac{1}{N}})$ [2,3], but this  theoretical feature is not central at this stage. (More on asymptotic  properties later).

(ii) In k-space the symmetries of interest  are [41]

$$c^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)=c^{(3)}\left({\mathbf{k}}_2,{\mathbf{k}}_1\right)=c^{(3)}\left({\mathbf{k}}_1,{-{\mathbf{k}}_1-\mathbf{k}}_2\right),$$

(18a)

$$S^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)=S^{(3)}\left({\mathbf{k}}_2,{\mathbf{k}}_1\right)=S^{(3)}\left({\mathbf{k}}_1,{-{\mathbf{k}}_1-\mathbf{k}}_2\right).$$

(18b)

(iii) Note that for $c_3$ in r-space one expects a decay to zero when  any of the particles increasingly separates from the other two, because this is  a reasonable physical behavior that can be drawn from Equations (14). However,  for $c^{(3)}$and $S^{(3)}$ in k-space their explicit behaviors may  become involved, since they depend on the angles between ${\mathbf{k}}_1$ and ${\mathbf{k}}_2$ (e.g., for $k\to \infty ,$ $S^{(3)}\left(k,k,\pi /3\right)\to 1$ is the observed behavior).

3.4. The Interplay between Simulation Techniques and Closures

There are well known MC and MD simulation  techniques to obtain $g_2,$ $g_3,S^{(2)}$ and $S^{(3)}$ [13,17,41,42,43,44,45,120,121,122,123,124],  and they will not be described here. Rather, the interest is to be focused on  some significant facts that show the interplay and complementarity between the  simulation calculations and the approximations that can be obtained via  closures.

(i) In the OZ2 case there are highly accurate  methods that yield $c_2(r)$ and $S^{\left(2\right)}(k)$ [5,6,13,71,72,125,126,127,152].  It is useful to point out that the knowledge of the latter functions allows one  to extend the range of distances obtained for $g_2(r)$ via simulation which, for instance, using a cubic  box is limited to half the box length $L/2.$ In addition to Equations (13), the basic formulas  are:

$$c^{(2)}\left(k\right)={\displaystyle \frac{4\pi }{k}}{\int }_0^{\mathrm{\infty }}dr\mathrm{}r\mathrm{}{c}_2\left(r\right)\sin kr,$$

(19)

$$g_2\left(r\right)=1+{\displaystyle \frac{1}{2{\pi }^2{r\rho }_N}}{\int }_0^{\mathrm{\infty }}dk\mathrm{}k\left(S^{\left(2\right)}\left(k\right)-1\right)\sin kr;r>0.$$

(20)

As stressed earlier, the access to the fluid  equation of state via $c_2\left(k=0\right)$ and ${\chi }_T$ is a key OZ2 feature. Moreover, the knowledge of $c^{(2)}\left(k\right),$ and hence of an extended $g_2\left(r\right),$ is useful for further studies, such as especial asymptotic decay properties of the pair correlations (pure exponential and exponentially damped oscillatory decays) [153,154,155] or some triplet correlation evaluations [3,41,156].

Among the OZ2 methods, Baxter’s partition of Equation (12) deserves especial attention [72,125], since it per se does not need any explicit knowledge of the underlying $u\left(r\right)$ interactions between the particles: the only input  parameters to carry out calculations are $\left\{g_2\left(r\right),{\rho }_N\mathrm{}\right\}$ (and T to extract the actual value of ${\chi }_T$) [46,130,136]. This contrasts sharply with typical OZ2 approaches, such as Percus-Yevick or the hypernetted chain [5,6,71,152], and makes Baxter’s applicability more general and powerful. Although the latter is a theoretical development in classical statistical mechanics, a good deal of its applications have been performed in the context of fluids with quantum behavior. Therefore, for the quantum purposes of this work, it is convenient at this point to give some specific details of Baxter’s partition [72] and its practical implementation, which involving a minimization procedure [125] (BDH $=$ Baxter-Dixon-Hutchinson) also needs a further  analysis of its results [46,60,62,127,131,136].

Baxter’s partition for a disordered fluid and its  practical implementation (BDH) are based on: (a) $c_2\left(r\right)$ vanishes for distances greater than a cut-off  distance $r>R_{CO};$ (b) $h_2\left(r\right)$ obtained via simulation is known for $r\le R_{CO}$; (c) $R_{CO}$ is limited by the range of distances set by half  the simulation box length $L/2$, and is fixed through a minimization procedure  that preserves the continuity of $h_2\left(r\right)$ at such cut off; (d) the minimization yields in  general a set of possible cut-off distances ${\left\{R_{Zv}\right\}}_{v=\mathrm{1,2},\dots ,v(M)},$ which are the “zeros” of an auxiliary function; a  practical convergence of the observable property, i.e., $S^{\left(2\right)}\left(k\right),$ occurs within an upper subset of these zeros; (e)  there are ways to treat the situation in which no zeros arise [125] and, also, to extract the physically  significant information from the upper $R_{Zv}$ zeros [46,62,136].  Clearly, the longer the L value is, the better the OZ2 results should  be. In this connection, to minimize the defects of the $g_2\left(r\right)$ fixed with a finite $N_S$ size in a simulation (e.g., canonical ensemble  calculations), an additional combination of OZ2 and simulation is highly  recommended. As a matter of fact, the iterative improvement on $g_2\left(r\right)$ with grand-canonical corrections [126,127] leads to a better-defined $S^{\left(2\right)}\left(k\right)$ convergence along the $R_{Zv}-$sequence$.$ An example of the latter iterations can be  summarized in the useful algorithm put forward by Baumketner and Hiwatari (BHw)  [127]:

$$g_2^{\left(GC\right)}\left(r;n+1\right)=g_2\left(r\right)\left(1+{\displaystyle \frac{S^{\left(2\right)}\left(k=0;n\right)}{N_S}}\right),n=\mathrm{0,1},2,\dots .,$$

(21)

where the initial conditions $(n=0)$ are the $N_S-$simulated $g_2\left(r\right)$ and its associated OZ2-$S^{\left(2\right)}\left(k=0;n=0\right)$ that yield the first correction $g_2^{\left(GC\right)}\left(r;n+1=1\right)$, which in turn is OZ2 analyzed to obtain its $S^{\left(2\right)}\left(k=0;n=1\right)$ giving the second correction $g_2^{\left(GC\right)}\left(r;n+1=2\right)$ and so on, until convergence in $S^{\left(2\right)}\left(k=0\right)$ is reached. When using OZ2-(BDH+BHw), the $S^{\left(2\right)}\left(k=0;n\right)$ value can be fixed at every step as the mean value  of the $S^{\left(2\right)}\left(k=0\right)$ arising from the upper set associated with the  resulting $R_{Zv}$ zeros [136]. For  further applications (e.g., triplet calculations [41,55,60])  involving a specific single function as an intermediate quantity, i.e., $c_2\left(r\right)$ or $S^{\left(2\right)}\left(k\right),$ the selection within the final upper $R_{Zv}$-range of convergence of one of them as a  representative is sufficient to give perfectly consistent results.

(ii) For OZ3 the situation involves the design of  accurate closure methods for obtaining $g_3\left(r_{12}{,r}_{13},r_{23}\right),$ or$c_3\left(r_{12}{,r}_{13},r_{23}\right)$ and $S^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right),$ which still remains a non-closed problem. For the $g_3$ functions key closures are:

-

Kirkwood superposition (KS3) [37]

$${\mathrm{}g}_{KS3}\left(r_{12}{,r}_{13},r_{23}\right)=g_2(r_{12})g_2\left(r_{13}\right)g_2\left(r_{23}\right).$$

(22)

-

Jackson-Feenberg convolution (JF3) [3,40]

$${\mathrm{}g}_{JF3}\left(r_{12}{,r}_{13},r_{23}\right)={\mathrm{}g}_{KS3}\left(r_{12}{,r}_{13},r_{23}\right)-h_2\left(r_{12}\right)h_2\left(r_{13}\right)h_2\left(r_{23}\right)+{\rho }_N\int d{\mathbf{q}}_4{h}_2\left(r_{14}\right)h_2\left(r_{24}\right)h_2\left(r_{34}\right).$$

(23)

-

The intermediate average (AV3) [59,60]

$${\mathrm{}g}_{AV3}\left(r_{12}{,r}_{13},r_{23}\right)={\displaystyle \frac{1}{2}}\left({\mathrm{}g}_{KS3}\left(r_{12}{,r}_{13},r_{23}\right)+{\mathrm{}g}_{JF3}\left(r_{12}{,r}_{13},r_{23}\right)\right).$$

(24)

The interesting point is that KS3 and JF3 are  pivotal references in general triplet studies. In this regard, one may  distinguish between triplets in r-space and  triplets in k-space: KS3 plays a central role in r-space [38,39,40,41,121,150,151], whereas JF3 does it in k-space  because $c_{JF3}\left(r_{12}{,r}_{13},r_{23}\right)\equiv 0.$ [41,70]. The  forms adopted by the KS3 and JF3 triplet structure factors can be found in  Reference [41], from which the case AV3 is a  trivial matter (note that for JF3 the $S^{(3)}$ formula reduces to the product of the three $S^{(2)}$ factors in Equation (15)). KS3 was derived in the context of classical statistical mechanics, while JF3 was in the context of quantum statistical mechanics. Despite these different conceptual frameworks, these two closures have been applied to both classical and quantum fluids indistinctly [26,27,32,37,38,39,40,41,54,55,56,57,58,59,60,121,150,151]. (In relation to this, recent results for fluids with quantum behavior show that AV3 performs very well as compared to KS3 or JF3 [59,60]). As regards their behaviors under limiting conditions, note that KS3 and JF3 satisfy the long range Equations (17b)-(17e); KS3 also satisfies the short-range Equation (17a), whilst JF3 fails to describe this latter situation [150,151].

Other triplet closures are orientated towards k-space [41,65,122,123,124,128]. Among them the insightful Barrat-Hansen-Pastore (BHP3) closure deals with $c_3\left(r_{12}{,r}_{13},r_{23}\right)$ [41]. By starting  from Equation (14a), the BHP3 closure introduces an auxiliary pair function $t(r)$ such that:

$${\mathrm{}c}_{BHP3}\left(r_{12}{,r}_{13},r_{23}\right)=t\left(r_{12}){t(r}_{13})t(r_{23}\right),$$

(25a)

which is fixed through the minimization of the  following nonnegative functional:

$$\mathsf{\Theta }\left[t\left(r\right)\right]=\int d\mathbf{r}\mathbf{}{\left\{{\displaystyle \frac{\partial c_2(r)}{\partial {\rho }_N}}-t(r)\int d\mathbf{s}\mathbf{}t\left(\mathbf{s}\right)t\left(\left|\mathbf{r}-\mathbf{s}\right|\right)\right\}}^2.$$

(25b)

To initiate the iterations $t^{(0)}\left(r_{12}\right)=h_2(r_{12})$ is the recommended choice [41]. The convergence is assessed via the quotient ${\displaystyle \frac{\mathsf{\Theta }\left[t\left(r\right)\right]}{{‖\frac{\partial c_2(r)}{\partial {\rho }_N}‖}^2}},$which involves the value of the functional being  minimized and the (squared) norm of the density derivative of the pair direct  correlation function, which serves as a reference. For the pertinent formal  details of this type of minimization, the reader is referred to Reference [41] (gradients), or to References [45,55,56] (conjugate gradients). After fixing $t\left(r\right),$ the application of Fourier transforms to obtain $c_{BHP}^{\left(3\right)}({\mathbf{k}}_1,{\mathbf{k}}_2)$ yields the approximation $S_{BHP}^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)$ [41]. BHP3  satisfies Equations (14) and has been applied to classical and quantum fluids [41,45,55,56,57,58,60,121,122,123,124].

Along the same k-line of thought, another  example is Denton-Ashcroft symmetrized approach (DAS3) [65], which can be cast as:

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

(26a)

in which the following definitions are to be used:

$$c_{DA}^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)={\displaystyle \frac{1}{c^{(1)\mathrm{\text{'}}}}}\left(c^{\left(2\right)}\left(k_1\right)c^{\left(2\right)\mathrm{\text{'}}}\left(k_2\right)+c^{\left(2\right)}\left(k_2\right)c^{\left(2\right)\mathrm{\text{'}}}\left(k_1\right)\right)-{\displaystyle \frac{c^{(1)\mathrm{\text{'}}\mathrm{\text{'}}}}{{{(c}^{(1)\mathrm{\text{'}}})}^2}}{c}^{\left(2\right)}\left(k_1\right)c^{\left(2\right)}\left(k_2\right),$$

(26b)

$$c^{(1)\mathrm{\text{'}}}=c^{\left(2\right)}\left(k=0\right),$$

(26c)

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

(26d)

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

(26e)

DAS3 satisfies the key Equation (14b) only at ${\mathbf{k}}_1={\mathbf{k}}_2=0$, and has been applied to classical and  quantum fluids [55,65].

(iii) The validity of the closure-based approaches  is to be established by comparison with exact simulation results. On the one  hand, closures at the pair level have proven their value as was mentioned  earlier. On the other hand, there are missing traits when triplets are  described in terms of pair properties, and one would not expect much success in  utilizing closures for triplets. Nonetheless, the triplet results reveal that  closures may capture salient features in both the r-space and the k-space.  Therefore, their usefulness as interpretive tools may be greater than imagined  (see References [41] and [122] for the classical domain and [59] and [60] for the quantum domain).

3.5. Other Theoretical Features

(i) There is another functional calculus way to  obtain the formulation of the pair $g_2$ Equation (3). Using the pairwise approach for the  interactions between particles in the fluid, this alternative is based on the  functional derivative of $\mathrm{l}\mathrm{n}{\mathsf{\Xi }}_C$ with respect to the variations in the  interparticle potential $\delta u\left(r_{jl}\right)$: such first functional derivative $\delta \mathrm{l}\mathrm{n}{\mathsf{\Xi }}_C/\delta u\left(r_{jl}\right)$ leads to $g_2$. If one iterates the procedure, the second  functional derivative ${\delta }^2\mathrm{l}\mathrm{n}{\mathsf{\Xi }}_C/\left(\delta u\left(r_{jl}\right)\delta u\left(r_{mn}\right)\right)$ produces a highly intricate expression. Given that  no external field is considered, the usefulness of these formal manipulations  seems rather reduced. The reader is referred to Reference [6] for the derivation of $g_2$ in the canonical ensemble. This issue will be  reconsidered when analyzing the quantum fluid instantaneous structures (Section 5).

(ii) An interesting relationship is obtained from  the grand canonical average of a function that depends on the particle  coordinates, $\mathsf{\Phi }\left({\mathbf{q}}_1,{\mathbf{q}}_2,\dots ,{\mathbf{q}}_N,\dots \right),$ (e.g., the n-point density dynamical  functions), and reads as [150,151,157]:

$$z{\left({\displaystyle \frac{\partial 〈\mathsf{\Phi }〉}{\partial z}}\right)}_{V,T}=〈N\mathsf{\Phi }〉-〈N〉〈\mathsf{\Phi }〉,$$

(27)

For the two-point density $\mathsf{\Phi }={\sum }_{j\ne l}\delta \left({\mathbf{r}}_j-{\mathbf{q}}_1\right)\delta \left({\mathbf{r}}_l-{\mathbf{q}}_2\right)$ one arrives at an exact connection between $g_2$

 and $g_3$ that can be cast as [26,32,150,151]:

$$S^{\left(2\right)}\left(k=0\right)\mathrm{}{\left({\displaystyle \frac{\partial \mathrm{l}\mathrm{n}{g}_2(r_{12})}{\partial {\rho }_N}}\right)}_T=\int d{\mathbf{q}}_3\left\{{\displaystyle \frac{g_3(r_{12},r_{13},r_{23})}{g_2(r_{12})}}-1\right\}+{\displaystyle \frac{2}{{\rho }_N}}\left(1-S^{\left(2\right)}\left(k=0\right)\right),$$

(28)

where $g_2\left(r_{12}\right)>0.$This relationship may be useful for checking the  consistency of the results computed for $g_2$ and $g_3,$ or for testing approximate triplet closures. This  latter issue will also be discussed in connection with fluids with quantum  behavior.

4. Basic Path Integral Concepts

4.1. PI Partition Functions

The starting point is the quantum grand-canonical  partition function ${\mathsf{\Xi }}_Q,$which for an isolated quantum system composed of N  identical particles can be cast as [4,5,7]

$${\mathsf{\Xi }}_Q\left(\mu ,\mathcal{V},T\right)={\sum }_{N\ge 0}\mathrm{e}\mathrm{x}\mathrm{p}\left(\beta \mu N\right)Z_{N,Q}(N,\mathcal{V},T)={\sum }_{N\ge 0}\mathrm{e}\mathrm{x}\mathrm{p}\left(\beta \mu N\right)Tr\left\{\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta H_0^{(N)}\right)\right\},$$

(29)

where Tr denotes the trace operation, $H_0^{(N)}$ is given by Equation (1), the mean number  density is given again by ${\rho }_N=〈N〉/\mathcal{V}$, and the other symbols retain the meanings seen  earlier.

If attention is focused on monatomic fluids at low  temperatures, such that exchange interactions can be neglected, the partition  function ${\mathsf{\Xi }}_Q$ for diffraction effects can be written in the  coordinate representation $\left\{{\mathbf{r}}^N\right\}$ as [5,7]:

$${\mathsf{\Xi }}_{QD}\left(\mu ,\mathcal{V},T\right)={\sum }_{N\ge 0}\mathrm{e}\mathrm{x}\mathrm{p}\left(\beta \mu N\right)Z_{N,QD}={\sum }_{N\ge 0}{\displaystyle \frac{\mathrm{e}\mathrm{x}\mathrm{p}\left(\beta \mu N\right)}{N!}}\int d{\mathbf{r}}^N⟨{\mathbf{r}}^N|\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta H_0^{(N)}\right)|{\mathbf{r}}^N⟩,$$

(30)

where $d{\mathbf{r}}^N=d{\mathbf{r}}_{1}d{\mathbf{r}}_2\dots d{\mathbf{r}}_N$ stands for the volume element of the N  distinguishable particles, and $|{\mathbf{r}}^N⟩=\mathrm{}|{\mathbf{r}}_1,{\mathbf{r}}_2,\dots ,{\mathbf{r}}_N⟩$ is the configurational state-ket. Now, with the  use of the path-integral approach (PI), ${\mathsf{\Xi }}_{QD}$ can be reduced to [7,9,11,74,77,78,79,81,84,90]

$${\mathsf{\Xi }}_{QD}\approx {\mathsf{\Xi }}_{PI,D}={\sum }_{N\ge 0}\mathrm{e}\mathrm{x}\mathrm{p}\left(\beta \mu N\right)Z_{N,QD}^{PI}={\sum }_{N\ge 0}{\displaystyle \frac{{\mathfrak{z}}^N}{N!}}\int {\prod }_{j=1}^N{\prod }_{t=1}^{P}d{\mathbf{r}}_j^t\times \mathrm{e}\mathrm{x}\mathrm{p}\left[-\beta W_{NP}({\mathbf{r}}^{NP})\right],$$

(31a)

in which every actual particle j (=1,2,…,N,…)  is transformed into a closed elastic necklace j with P beads  labelled $t=\mathrm{1,2},\dots ,P.$Thus, one finds a PI-model system composed of $N\times P$-sized canonical ensembles$,$ each ruled by an  “effective potential” $W_{NP}({\mathbf{r}}^{NP})$ depending upon all the corresponding bead  coordinates $\left\{{\mathbf{r}}_j^t\right\}$, such that

$$W_{NP}\left({\mathbf{r}}^{NP}\right)=W_{NP}\left({\mathbf{r}}_1^1,{\mathbf{r}}_1^2,\dots ,{\mathbf{r}}_1^P,{\mathbf{r}}_2^1,{\mathbf{r}}_2^2,\dots ,{\mathbf{r}}_2^P,{\dots ,{\mathbf{r}}_N^1,{\mathbf{r}}_N^2,\dots ,\mathbf{r}}_N^P;\beta ,\mathrm{\hslash },m\right).$$

(31b)

The parameter $\mathfrak{z}$ is a generalized activity that reads as follows:

$$\mathfrak{z}=\mathrm{e}\mathrm{x}\mathrm{p}\left(\beta \mu \right){\left({\displaystyle \frac{\sqrt{P}}{{\lambda }_B}}\right)}^{3P},$$

(31c)

which helps to identify the expression of the  canonical $Z_{N,QD}^{PI}.$ A crucial step in PI work is to give suitably  operational forms to $W_{NP}({\mathbf{r}}^{NP})$ by resorting to propagator calculus. In addition,  the formal generalizations to tackle quantum exchange regimes rely heavily on  the basic development for distinguishable particles [4,7,9].

4.2. Sum over Histories and Propagators

The PI-P picture (31a) becomes exact in the  limit $P\to \infty$ (Trotter’s theoretical accuracy [9,138]), although one can find optimal finite P values  that yield sufficiently precise and accurate representations of the actual  system. The optimal P is thus a compromise between theoretical accuracy  and statistical convergence for properties. All of this is related to the  design of (thermal) propagators that serve to build $W_{NP}$ [7,9,75,81,82,83,84]. A  brief explanation of how the practical image of beads-and-necklaces arises is  worth giving here.

The basic idea is to split the density matrix  diagonal elements $⟨{\mathbf{r}}^N|\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta H_0^{(N)}\right)|{\mathbf{r}}^N⟩$ into, in general, a number X of factors  involving operators $\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta H_0^{\left(N\right)}/X\right),$ and to employ the completeness relationship in the  coordinate representation $\int d{\mathbf{r}}^N|{\mathbf{r}}^N⟩⟨{\mathbf{r}}^N|=1$. By doing so, one obtains the exact convolution  property [9]:

$$⟨{\mathbf{r}}^N|\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta H_0^{(N)}\right)|{\mathbf{r}}^N⟩=\int d{\mathbf{r}}^{N,2}\dots d{\mathbf{r}}^{N,X}⟨{\mathbf{r}}^{N,1}|\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta H_0^{\left(N\right)}/X\right)|{\mathbf{r}}^{N,2}⟩\dots ⟨{\mathbf{r}}^{N,X}|\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta H_0^{\left(N\right)}/X\right)|{\mathbf{r}}^{N,1}⟩,$$

(32)

which leads to the canonical partition function:

$$Z_{N,QD}={\displaystyle \frac{1}{N!}}\int d{\mathbf{r}}^{N,1}\dots d{\mathbf{r}}^{N,X}{\prod }_{\tau =1}^{X(\mathrm{*})}⟨{\mathbf{r}}^{N,\tau }|\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta H_0^{\left(N\right)}/X\right)|{\mathbf{r}}^{N,\tau +1}⟩,$$

(33)

where $\tau =\mathrm{1,2},\dots ,X,$ $d{\mathbf{r}}^{N,\tau }={d\mathbf{r}}_1^{\tau }d{\mathbf{r}}_2^{\tau }\dots d{\mathbf{r}}_N^{\tau }$, $|{\mathbf{r}}^{N,\tau }⟩=|{\mathbf{r}}_1^{\tau },{\mathbf{r}}_2^{\tau },\dots ,{\mathbf{r}}_N^{\tau }⟩,$ $|{\mathbf{r}}^{N,1}⟩=|{\mathbf{r}}^N⟩,$ and the (*) in the product implies the cyclic property $\tau +1=X+1\equiv 1$. The thermal propagator is the generic nonnegative matrix element $⟨{\mathbf{r}}^{N,\tau }|\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta H_0^{\left(N\right)}/X\right)|{\mathbf{r}}^{N,\tau +1}⟩$ [9], which clearly corresponds to dealing with a density matrix at a higher temperature than the actual.

The whole process contained in Equations (32) and (33) can be visualized in an appealing way by resorting to the image of the periodic “motion” of the system in imaginary time (period $\beta \hslash )$ [4,140]. Thus, in Equation (32) the system may be viewed as if it “travelled” from ${\mathbf{r}}^N$ to ${\mathbf{r}}^N$ by following every possible broken path (in configuration space) defined by the intermediate stages that can be associated with the equally spaced $\tau$-“instants” in imaginary time. Every particle j follows then a closed path (necklace), ${\mathbf{r}}_j^1\to {\mathbf{r}}_j^2\to \dots \to {\mathbf{r}}_j^X\to {\mathbf{r}}_j^1,$ marked out by the $\tau$-“instants” (beads). In the language of continuous paths $\left(P\to \infty \right)$ one speaks of the basic interval $\left[0,\beta \hslash \right]$ for the motion in imaginary time to take place, and there is the obvious equivalence $0\equiv \beta \hslash ,$ as is customary in every Fourier series analysis. (Warning: this does not imply by any means that $\beta =0!).$The partition function as given in Equation (33) closes the situation by integrating over every possible starting point ${\mathbf{r}}^{N,1}$ and, in this way, the “sum over histories” in imaginary time arises. It is easy to see that there is equivalence among all the X imaginary-time instants $\tau$ (translational invariance in imaginary time).

To give operational forms to the propagator $⟨{\mathbf{r}}^{N,\tau }|\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta H_0^{\left(N\right)}/X\right)|{\mathbf{r}}^{N,\tau +1}⟩$ one needs approximations to tackle the problem posed by the nonzero commutator $\left[K_{kin}^{(N)},V^{(N)}\right]\ne 0$ [7,9,81,82,83]. In relation to this, several possibilities are available in the PI literature. The accuracy attainable with the approximation $Z_{N,QD}\approx Z_{N,QD}^{PI}$ depends strongly on the propagator selected and the X value employed ${(Z}_{N,QD}^{PI}$ is defined in Equations (31)). For the purposes of this work, only a few guiding lines are to be commented, the reader being referred to the pertinent references cited below for complete details.

In the study of quantum many-body systems three main types of constructions for the propagator stand out, which lead to: (i) the primitive propagator (PP), for which conventionally the number of beads is $P=X$ and yields an accuracy $O(P^{-2})$ for $Z_{N,QD}^{PI}$ [4,7,9,61,80]; (ii) the fourth-order propagator (SCVJ), arising from the works by Suzuki [82], Chin [83], and Voth et al. [84], for which the final number of beads P must be even $(P=2X)$ and yields an accuracy $O(P^{-4})$ for $Z_{N,QD}^{PI};$ and (iii) the pair actions (PA) based on the developments put forward by Ceperley [9] and by Cao and Berne (a specific PA for quantum hard spheres, CBHSP) [81], for which the number of beads is $P=X$ again and there is no systematic recipe for their $Z_{N,QD}^{PI}$ accuracies. It is useful to recall at this point that another version of a fourth-order propagator, the Takahasi-Imada propagator [75,76], is not very useful for computing correlations and structures [10,11,75]. Consequently, the PI computational effort depends on the size of the optimal P value, which depends on the propagator effectiveness, a property that varies from type to type. One may summarize this issue by saying that: pair actions are more efficient since they appear to require lower optimal P values, whereas high-order propagator expressions can be improved, seeking to reduce P, in a systematic way. (For early hard-sphere PA’s constructions and applications based on Barker’s work [73] see References [79,158,159], and also [11] for a general discussion of this issue).

The $W_{NP}$ expressions of the PP, CBHSP and SCVJ propagators are given in what follows. Hereafter, the P value is assumed to be the optimal, and the (final) P beads are numbered correlatively $t=\mathrm{1,2},\dots ,P.$ The three cases can be written as

$$W_{NP}=W_{1F}+W_2+W_3,$$

(34)

where $W_{1F}$ is common to the three and can be cast as [4,7]

$$W_{1F}\left(\mathrm{P}\mathrm{P}\right)=W_{1F}\left(\mathrm{C}\mathrm{B}\mathrm{H}\mathrm{S}\mathrm{P}\right)=W_{1F}\left(\mathrm{S}\mathrm{C}\mathrm{V}\mathrm{J}\right)={\displaystyle \frac{mP}{2{\beta }^2{\mathrm{\hslash }}^2}}{\sum }_{j=1}^N{\sum }_{t=1}^{P(\mathrm{*})}{\left({\mathbf{r}}_j^t-{\mathbf{r}}_j^{t+1}\right)}^2,$$

(35)

with the (*) implying the cyclic property $t+1=P+1\equiv 1.$ This contribution arises from the free particle behaviors, which within PI take the form of harmonic couplings between adjacent beads in every necklace.

By assuming for simplicity pairwise interactions between the actual particles, and denoting the inter-necklace bead-bead distances by $r_{jl}^t=\left|{\mathbf{r}}_j^t-{\mathbf{r}}_l^t\right|,$ $t=\mathrm{1,2},\dots ,P,$ (the equal t label is to be noticed) the contributions $W_2$ and $W_3$ read as [7,9,62,81,84]:

(a) Primitive propagator

$$W_2\left(PP\right)={\displaystyle \frac{1}{P}}{\sum }_{j<l}{\sum }_{t=1}^{P}u(r_{jl}^t),$$

(36a)

$$W_3\left(PP\right)=0.$$

(36b)

(b) Cao-Berne propagator for quantum hard spheres (diameter ${\sigma }_{HS})$

$$W_2\left(\mathrm{C}\mathrm{B}\mathrm{H}\mathrm{S}\mathrm{P}\right)={\displaystyle \frac{1}{P}}{\sum }_{j<l}{\sum }_{t=1}^P{u}_{HS}\left(r_{jl}^t\right)=\left\{\begin{array}{c}0if{r}_{jl}^t>{\sigma }_{HS}\\ \infty if{r}_{jl}^t<{\sigma }_{HS}\end{array}\right.,$$

(37a)

$$W_3\left(\mathrm{C}\mathrm{B}\mathrm{H}\mathrm{S}\mathrm{P}\right)=-{\displaystyle \frac{1}{\beta }}\mathrm{l}\mathrm{n}{\prod }_{j<l}{\prod }_{t=1}^{P(\mathrm{*})}\left[1-{\displaystyle \frac{{\sigma }_{HS}\left(r_{jl}^t+r_{jl}^{t+1}-{\sigma }_{HS}\right)}{r_{jl}^t{r}_{jl}^{t+1}}}\times \right.$$

$$\left.\mathrm{e}\mathrm{x}\mathrm{p}\left(\left.-{\displaystyle \frac{mP}{2\beta {\mathrm{\hslash }}^2}}\left(r_{jl}^t-{\sigma }_{HS}\right)\left(r_{jl}^{t+1}-{\sigma }_{HS}\right)\left(1+\cos \left({\mathbf{r}}_{jl}^t,{\mathbf{r}}_{jl}^{t+1}\right)\right)\right)\right.\right].$$

(37b)

(c) SCVJ propagator

$$W_2\left(\mathrm{S}\mathrm{C}\mathrm{V}\mathrm{J}\right)={\displaystyle \frac{2}{3P}}{\sum }_{j<l}\left\{{\sum }_{odd-t}u\left(r_{jl}^t\right)+2{\sum }_{even-t}u\left(r_{jl}^t\right)\right\},$$

(38)

$$W_3\left(\mathrm{S}\mathrm{C}\mathrm{V}\mathrm{J}\right)={\displaystyle \frac{{\beta }^2{\mathrm{\hslash }}^2}{9m{P}^3}}{\sum }_{j=1}^N\left\{\alpha {\sum }_{odd-t}{\left({\sum }_{j\ne l}{\displaystyle \frac{du(r_{jl}^t)}{d{r}_{jl}^t}}{\mathsf{\eta }}_{jl}^t\right)}^2+(1-\alpha ){\sum }_{even-t}{\left({\sum }_{j\ne l}{\displaystyle \frac{du(r_{jl}^t)}{d{r}_{jl}^t}}{\mathsf{\eta }}_{jl}^t\right)}^2\right\},$$

(38b)

where ${\mathsf{\eta }}_{jl}^t={\displaystyle \frac{{\mathbf{r}}_{jl}^t}{r_{jl}^t}}$ is a unit vector, and the parameter $\alpha \in \left[\mathrm{0,1}\right],$ although for a good deal of applications $\alpha =1/3$ is a recommended value [84]. (In Equations (38) note that necessarily $u\left(r_{jl}^t\right)\ne u_{HS}\left(r_{jl}^t\right)$).

4.3. The Classical Isomorphism, Bead Roles, and Notational Conventions

(i) The comparison of the quantum ${\mathsf{\Xi }}_{PI,D}$ given in Equations (31) with the classical ${\mathsf{\Xi }}_C$ given in Equation (2) leads to the so-called classical isomorphism [7,9,16], or more properly semi-classical isomorphism, due to the dependences shown in Equation (31b) (see References [2,11,13] for the semiclassical approach forms). This fact makes the PI computational study of quantum fluids amenable to using well-known basic classical simulation techniques (MC and MD), such as path integral Monte Carlo (PIMC) and path integral molecular dynamics (PIMD). The latter basic techniques (and its variants) have become standard tools in modern quantum condensed matter research [9,10,11,18,19,21,22,33,34,35,36,54,55,56,57,58,59,60,61,62,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101]. Naturally, there are unavoidable differences between the classical and the PI developments and computations, since the language of classical dynamical functions is radically different from that of quantum operators [4,5]. However, as regards the studies of the structures of fluids with quantum behavior, the relevance of such isomorphism runs deeper than one might think at first sight [9,11,18,19,55,93].

(ii) The inter-necklace interactions are uniquely defined as averages involving equal imaginary time t-bead interactions, which is connected to the diagonal character of the potential energy operator $V^{(N)}.$(iii) $W_3(\mathrm{C}\mathrm{B}\mathrm{H}\mathrm{S}\mathrm{P})$ contains kinetic effects between consecutive imaginary times (beads) in different necklaces.

(iv) The singularity of the $u_{HS}$ potential (37a) in CBHSP prevents the numerical evaluation of its derivative $d{u}_{HS}/dr$ for a given (“frozen”) configuration of necklaces [136,160].

(v) An inspection of the PP, SCVJ, and CBHSP formulas indicates that, when using these propagators, careful attention must be paid to the equivalence among the beads. The point is that while there is full equivalence between all the $N\times P$ beads in PP and CBHSP applications, such symmetry is broken when use of SCVJ is made. The SCVJ beads are classified into two groups, conventionally the odd-numbered $t=\mathrm{1,3},\dots ,P-1,$ and the even-numbered $t=\mathrm{2,4},\dots ,P.$ In thermodynamic calculations the general roles played by the beads are as follows: the same in PP [7,9,74,80], the same in CBHSP [9,81], but different in SCVJ depending on the odd/even group they belong to [62,84]. Moreover, for the current purposes of structural evaluations, it is worth stressing that: (a) in PP and CBHSP all the beads play the same role in the corresponding ensemble averages [9,11,79,80,81], whereas (b) in SCVJ only the odd-numbered beads are physically significant, thereby playing the same role in such ensemble averages [11,84] (the even-numbered beads must be discarded). Technically, these different traits arise from the use of a double-time step$\left(\mathrm{i}.\mathrm{e}.,\beta \hslash /({\displaystyle \frac{P}{2}})\right)$ when constructing SCVJ [82,83,84] instead of the usual single-time step $\left(\mathrm{i}.\mathrm{e}.,\beta \hslash /P\right)$ applied in the PP and CBHSP derivations (the reader is referred to the quoted references for full details). These structural bead features will be formalized in Section 5.

In order not to burden the notations, the conventions for structures followed in this article are worth summarizing now. Given the optimal number of beads $P$ per necklace:

- The three symbols for denoting beads are: t, $\tau ,$ and ${\tau }^{*},$ depending on the context.

- Label t runs over $t=\mathrm{1,2},3,\dots ,P,$ and is associated with every necklace j.

- Label $\tau$ runs over $\tau =\mathrm{1,2},3,\dots ,X,$ and is associated with every necklace j.

- Label ${\tau }^{*}$ runs over ${\tau }^{*}=\mathrm{1,2},3,\dots ,NX$, and serves the purpose of a whole renumbering of all the structurally significant beads in a canonical sample, such that they are associated with the necklaces as ${\tau }^{*}={{\mathrm{1,2},\dots ,X}_{⏟}}_{j=1},\dots ,{{X\left(N-1\right)+1,\dots ,NX}_{⏟}}_{j=N}$.

Depending on the propagator the above conventions mean:

- For the PP and CBHSP propagators, in which $X=P\to \tau =t=\mathrm{1,2},3,\dots ,P,$ and ${\tau }^{*}=\mathrm{1,2},3,\dots ,NP.$

- For the SCVJ propagator, in which $X={\displaystyle \frac{P}{2}}\to \tau =\mathrm{1,2},3\dots ,X$ $\equiv t=\mathrm{1,3},5,\dots ,P-1$, and ${\tau }^{*}=\mathrm{1,2},3,\dots ,NX\equiv$ the global picking of the odd-numbered t-beads in every necklace. (Note that the even-numbered beads represent intermediate stages that fit mathematically into the whole picture of P-membered necklaces).

4.4. Quantum Exchange Interactions

The situation gets much more complex when one deals with monatomic fluids in which quantum exchange interactions cannot be neglected. For now, it is worthwhile to remark that bosonic exchange can be successfully tackled with PI computations, as shown about forty years ago by Ceperley and Pollock [78] and by others later [9,35,36,91,96,97]. However, the known as the “sign problem” [100,101] made fermionic exchange remain out of the PI reach. However, this situation has changed for the better very recently, thanks to the appealing approach put forward by Filinov et al. using the Wigner extension of the usual PI formalism (WPIMC) [98,99].

The intricacies of the exchange issue can be grasped by observing the form of the canonical partition function for the general case of N particles not only identical, but also indistinguishable, all of them in the same spin state [4,9,100]:

$$Z_{N,Qexch.}(N,\mathcal{V},T)={\displaystyle \frac{1}{N!}}{\sum }_{\mathrm{\wp }}{\mathfrak{F}}_{\mathrm{\wp }}\int d{\mathbf{r}}^N⟨{\mathbf{r}}^N|\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta H_0^{(N)}\right)|\mathrm{\wp }{\mathbf{r}}^N⟩,$$

(39)

where $\wp$ runs over the $N!$ permutations of the particles, ${\mathfrak{F}}_{\wp }=+1$ for every $\wp$ if there is bosonic exchange, ${\mathfrak{F}}_{\wp }=\pm 1$ if there is fermionic exchange and such that $+1$ holds for even $\wp$ and $-1$ holds for odd $\wp ,$ and $\left|\wp {\mathbf{r}}^N⟩\right.$ is defined as $\left|\wp {\mathbf{r}}_1,\wp {\mathbf{r}}_2,\dots ,\wp {\mathbf{r}}_N⟩\right.,$ which is meant to denote the action of the permutation $\wp$ on the ordered ket $\left|{\mathbf{r}}_1,{\mathbf{r}}_2,\dots ,{\mathbf{r}}_N⟩\right.$ denoting the position states of N distinguishable particles. As seen, a pure summation of nonnegative terms arises from the bosonic case, whereas an alternating summation does from the fermionic case. Accordingly, feasible calculations can be carried out for bosons since a well-defined probability density is always involved. However, given the enormous size that $N!$ can reach, the same does not occur for fermions, for which large uncertainties plague in practice PI direct calculations (it seems out of the question that one could have access today to unlimited computational resources that would yield the exact answer to the whole fermionic sum). In their WPIMC works, Filinov et al. overcome the “sign problem” by treating the fermionic exchange determinant via a positive semidefinite Gram determinant [161].

Some observations regarding Equation (39) are in order. (a) By keeping the identity permutation alone, one retrieves the pure diffraction effects $Z_{N,QD}$ as given in Equations (30)-(31); this results in the non-classical Boltzmann statistics associated with $Z_{N,QD}^{PI}.$ (b) Note that the presence of ${\left(N!\right)}^{-1}$ guarantees the correct transition to the classical partition function and also the access to thermal properties (free energies, entropy) [162]; this is a factor which must always be written when dealing with $Z_{N,QD}.$ (c) If more than one spin state were involved, the expression (39) would have to be generalized by including every spin possibility [100,163].

The PI triplet structural issues of a zero-spin bosonic fluid will be addressed in Appendix I. However, this does not yet seem the time to enter a detailed discussion of the PI fermionic fluid case, and the interested reader is referred to the following references for basic information [98,99,164,165,166,167,168,169].

4.5. The PI-Centroid Variable

A particularly interesting PI quantity is the centroid variable ${\mathbf{R}}_{CM,j}$ which is the necklace “center of mass” (all the necklaces are equivalent, so are their centroids) [4,18,19,140,141,142]. The corresponding definitions for the propagators discussed above are the following [61,80,84,159]:

$${\mathbf{R}}_{CM,j}=\left\{\begin{array}{c}{P}^{-1}{\sum }_{t=1}^P{\mathbf{r}}_j^t;(\mathrm{P}\mathrm{P},\mathrm{C}\mathrm{B}\mathrm{H}\mathrm{S}\mathrm{P})\\ {\left(P/2\right)}^{-1}{\sum }_{odd-t}{\mathbf{r}}_j^t;(\mathrm{S}\mathrm{C}\mathrm{V}\mathrm{J})\end{array}.\right.$$

(40)

This variable is of paramount importance in both the thermodynamic and the structural studies of fluids with quantum behavior. As regards equilibrium situations, one may mention the structural relationships (OZ2!) with the equation of state [130,136], the structural extension of the semi-classical isomorphism [11,55,130], the connections with observable quantities and special situations [62,135,137,143,144,145], or the applications to deal with quantum exchange interactions [130,132,133,134], as will be discussed later for the zero-spin bosonic case in Appendix I. Furthermore, centroids are key to building interesting approaches to quantum dynamics in condensed phases [18,19,33,34,93].

5. Theory of Equilibrium Quantum Fluid Triplet Structures under Diffraction Effects

A homogeneous and isotropic quantum monatomic fluid is considered. Exchange interactions are neglected and the Hamiltonian $H_0^{(N)}$ Equation (1) is selected. The correlation functions and their associated structure factors can be grouped into three classes [7,9,11,84]: centroids (CMn), total continuous linear response (TLRn), and instantaneous (ETn). Each class is associated with the linear response from the fluid to an external weak field $\mathsf{\Psi }$. When use of functional calculus techniques is made, a good deal of the developments for the CMn class are parallel to that of the classical monatomic fluid, and the same may be said of the TLRn class, though to a significantly lesser extent [7,11,55,130]. The ETn case has a different nature, the functional techniques that can be applied to CMn and TLRn are inapplicable to ETn, and this latter class is to be treated separately [3,6,9,84,135].

The grand-canonical partition function for the study of the CMn and TLRn classes incorporates the action of $\mathsf{\Psi }$ and reads as follows [11,62]:

$${\mathsf{\Xi }}_{QD}\left(\mu ,\mathcal{V},T;\mathsf{\Psi }\right)={\sum }_{N\ge 0}{\displaystyle \frac{\mathrm{e}\mathrm{x}\mathrm{p}\left(\beta \mu N\right)}{N!}}\int d{\mathbf{r}}^N⟨{\mathbf{r}}^N|\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta H_0^{\left(N\right)}-\beta {\mathsf{\Psi }}^{(N)}\right)|{\mathbf{r}}^N⟩,$$

(41)

where the operator ${\mathsf{\Psi }}^{(N)}$ is defined as ${\mathsf{\Psi }}^{(N)}={\sum }_{j=1}^N\mathsf{\Psi }({\mathbf{r}}_j).$ The key point here is how to deal with the exponential of the sum of two non-commutative operators, i.e., $\left[H_0^{(N)},{\mathsf{\Psi }}^{(N)}\right]\ne 0$, which can be accomplished in the standard way via the Baker-Campbell-Hausdorff formula [163]. For the current purposes, the density matrix elements in Equation (41) are split into X steps and one applies the basic approximation:

$$\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{\beta }{X}}(H_0^{\left(N\right)}+{\mathsf{\Psi }}^{\left(N\right)})\right)\approx \mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{\beta }{2X}}{\mathsf{\Psi }}^{\left(N\right)}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{\beta }{X}}{H}_0^{\left(N\right)}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{\beta }{2X}}{\mathsf{\Psi }}^{\left(N\right)}\right),$$

(42)

which is accurate up to $O(X^{-3})$ terms (as occurs in the PP propagator derivation). This approximation leads to the partition function:

$${\mathsf{\Xi }}_{QD}\left(\mathsf{\Psi }\right)={\sum }_{N\ge 0}{\displaystyle \frac{\mathrm{e}\mathrm{x}\mathrm{p}(\beta \mu N)}{N!}}\int {\prod }_{\tau =1}^{X(\mathrm{*})}\left\{d{\mathbf{r}}^{N,\tau }⟨{\mathbf{r}}^{N,\tau }|\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{\beta }{X}}{H}_0^{\left(N\right)}\right)|{\mathbf{r}}^{N,\tau +1}⟩\right\}\times \mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{\beta }{X}}{\sum }_{j=1}^N{\sum }_{\tau =1}^X\mathsf{\Psi }({\mathbf{r}}_j^{\tau })\right).$$

(43)

$({\mathsf{\Psi }}^{(N)}$ is diagonal in the coordinate representation). By observing Equation (43) one notes that the nonnegative density matrix elements may be given the final PI propagator form desired or needed for PI convergence reasons. In relation to this, notice that: (a) functional derivations involving $\mathsf{\Psi }$ are independent of such form, (b) the linear response from the fluid is characterized by functions in the limit $\mathsf{\Psi }\to 0,$ and (c) the optimal value for X can therefore be conveniently adapted. By recalling that, regardless of the propagator employed, the final optimal number of beads in PI calculations is conventionally taken as P, the latter remark (c) is unimportant for PP and CBHSP (or PA’s in general), since $X=P$. However, it turns out to be crucial for SCVJ for which $X=P/2,$ and the physically significant beads for structural purposes are just those contained in Equation (43): $\tau =\mathrm{1,2},3,\dots ,X$. (Note that ${\mathsf{\Xi }}_{QD}\left(\mathsf{\Psi }\right)$ is the primary object that serves to define the structural ensemble averages). As stressed earlier, such significant $\tau$ beads are the odd-numbered t ones when SCVJ is fully developed and the intermediate beads completing the whole P-set do arise, i.e., $t=1,2,3,\dots ,\mathit{P}-1,P$ [11,84]. Consequently, ${\mathsf{\Xi }}_{QD}\left(\mathsf{\Psi }\right)$ can be approximated by the general PI form [11,62]:

$${\mathsf{\Xi }}_{PI,D}\left(\mathsf{\Psi }\right)={\sum }_{N\ge 0}{\displaystyle \frac{{\mathfrak{z}}^N}{N!}}\int {\prod }_{j=1}^N{\prod }_{t=1}^{P}d{\mathbf{r}}_j^t\times \mathrm{e}\mathrm{x}\mathrm{p}\left[-\beta W_{NP}({\mathbf{r}}^{NP})\right]\times \mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{\beta }{X}}{\sum }_{j=1}^N{\sum }_{\tau =1}^X\mathsf{\Psi }({\mathbf{r}}_j^{\tau })\right).$$

(44)

Now, it is worthwhile to point out that the action of the field expressed in Equation (44) is consistent with a thermalized $\mathsf{\Psi }$-interaction, but not with a “sudden” interaction such as that of the radiation-atom elastic collisions. On the one hand, examples of the fields in 3-D space involved in Equation (44) are: static continuous fields in r-space and continuous fields in imaginary time, such as a gravitational field or the special case of a set of neutrons thermalized with a fluid sample composed of zero-spin particles, as is the case of fluid helium-4 [30]. (For the latter, such continuous interaction takes place in Trotter’s limit $X\to \infty$, which is PI treated in (44) via the finite number of beads, P or $P/2)$. The classes CMn and TLRn are associated with these types of continuous $\mathsf{\Psi }$-fields. On the other hand, examples of related fields not involved in Equation (44) correspond to the elastic scattering of radiation experiments, e.g., neutrons for fluid helium-4, X-rays for fluids helium-4 and helium-3 [6,28] (from a conceptual statistical mechanical point of view there is no difference between the static response functions to both radiation probes [6]). Thus, within the framework set in subsection 2.2, the inadequacy of the static picture related to Equation (44) can be illustrated, for simplicity, by considering the interaction between a neutron and a spinless particle, which is essentially given by the neutron-nucleus Fermi’s potential (a Dirac-$\delta$) [6,23]: the (elastic) collision between an incoming neutron and a thermally delocalized atom localizes the atom at just a position in the sample (the “collapse” by the position measurement); afterwards the atom must delocalize again according to the fluid equilibrium state. “Viewed” from the foregoing PI approach, the neutron-atom collision should take place at the position of a given bead, which would imply the immediate disappearance of the associated necklace distribution and its ensuing reappearance [11,170]. This quantum phenomenon cannot therefore be analyzed via the functional derivatives of (44), which cannot cope with the disappearance/ reappearance of the atom quantum thermal packet, and a full quantum time-dependent treatment is the way to understand it [6]. Therefore, for the class ETn, which is associated with localizing $\mathsf{\Psi }$-fields, Equation (44) is not appropriate (more about this in subsections 5.3 and 5.4). In sharp contrast, under classical conditions the radiation-atom elastic collision phenomena pose no problems to the functional derivations and their static linear response interpretation [5,6,148]. Now, a detailed discussion of the three classes of equilibrium quantum structures follows.

5.1. The PI-Centroid CMn Class

This class behaves from the structural standpoint in the classical way. This is a remarkable result that arises from the explicit consideration of a continuous weak field ${\mathsf{\Psi }}_F$ of constant strength, $\mathbf{f}=\left(f_x,f_y,f_z\right),$ acting on the fluid [19,33,55,84,130]:

$${\mathsf{\Psi }}_F^{(N)}={\sum }_{j=1}^N{\mathsf{\Psi }}_F\left({\mathbf{r}}_{\mathit{j}}\right)={\sum }_{j=1}^N\mathbf{f}·{\mathbf{r}}_{\mathit{j}}.$$

(45)

Inclusion of this field leads to the general PI partition function [130]:

$${\mathsf{\Xi }}_{PI,D}({\mathsf{\Psi }}_F)={\sum }_{N\ge 0}{\displaystyle \frac{{\mathfrak{z}}^N}{N!}}\int {\prod }_{j=1}^N{\prod }_{t=1}^{P}d{\mathbf{r}}_j^t\times \mathrm{e}\mathrm{x}\mathrm{p}\left[-\beta W_{NP}({\mathbf{r}}^{NP})\right]\times \mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{\beta }{X}}{\sum }_{j=1}^N{\sum }_{\tau =1}^X{\mathbf{f}·\mathbf{r}}_j^{\tau }\right)=$$

(46a)

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

(46b)

where use of Equation (40) is made. For clarity reasons, it is convenient to rewrite the previous formula as:

$${\mathsf{\Xi }}_{PI,D}\left({\mathsf{\Psi }}_F\right)={\sum }_{N\ge 0}{\displaystyle \frac{{\mathfrak{z}}^N}{N!}}\int {\prod }_{t=1}^{P}d{\mathbf{r}}^{N,t}\times \mathrm{e}\mathrm{x}\mathrm{p}\left[-\beta W_{NP}\left({\mathbf{r}}^{NP}\right)\right]\times$$

$$\prod _{j=1}^{N}d{\mathbf{R}}_j\delta ({\mathbf{R}}_j-{\mathbf{R}}_{CM,j})\times \mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta \sum _{j=1}^N{\mathsf{\Psi }}_F\left({\mathbf{R}}_j\right)\right).$$

(47)

It is straightforward to show that the functional derivatives of $\ln {\mathsf{\Xi }}_{PI,D}\left({\mathsf{\Psi }}_F\right)$ are formally identical to those of the classical case contained in Equations (8)-(9). Thus, one finds up to the third order [11,57,130]:

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

(48a)

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

(48b)

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

(48c)

where the inhomogeneities caused by ${\mathsf{\Psi }}_F$ are reflected in the fluid centroid functionals/properties ${\mathsf{\Gamma }}_{CMn}.$The CMn hierarchy can be obtained by iterating the procedure to higher orders. As explained earlier in the classical case, if one makes ${\mathsf{\Psi }}_F=0$ on the right-hand sides, the isolated fluid properties arise:

$${\mathsf{\Gamma }}_{CM1}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)={\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)\to {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F=0\right)={\rho }_N,$$

(49a)

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

(49b)

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

(49c)

where $R_{jl}=\left|{\mathbf{R}}_j-{\mathbf{R}}_l\right|,$ and use is made of the total pair correlation function of centroids $h_{CM2}=g_{CM2}-1.$ Incidentally, note that the factor accompanying ${\rho }_N^3$ may be abbreviated to $h_{CM3}\left(R_{12},R_{13},R_{23}\right)$, consistently with the subtraction from $g_{CM3}$ of the asymptotic behaviors$.$ Note that Equations (49b) and (49c) can be formulated in a concise manner with the use of the spontaneous fluctuations of the centroid density (analogous with Equations (9b) and (9c)).

The formal expressions for the grand-canonical averages giving $g_{CM2}$ and $g_{CM3}$ are identical to those of the classical case Equations (3) and (4) (Figure 1):

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

(50)

$${{\rho }_{N,CM}^{\left(3\right)}\left({\mathbf{R}}_1,{\mathbf{R}}_2,{\mathbf{R}}_3\right)={\rho }_N^{3}g}_{CM3}\left({\mathbf{R}}_1,{\mathbf{R}}_2,{\mathbf{R}}_3\right)={\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{R}}_1\right)\delta \left({\mathbf{R}}_{CM,l}-{\mathbf{R}}_2\right)\delta \left({\mathbf{R}}_{CM,m}-{\mathbf{R}}_3\right)〉.$$

(51)

5.1.1. PI-Centroid Linear Response

By integrating (48b) and (48c) with the substitutions (49b) and (49c) on the right-hand sides, and taking the Fourier transforms, the two- and three-body centroid response functions arise as [11,55]:

$${\delta \rho }_{N,CM}^{\left(1\right)}\left(\mathbf{k};{\mathsf{\Psi }}_F\right)\approx -\beta {\rho }_N\delta {\mathsf{\Psi }}_F\left(\mathbf{k}\right)S_{CM}^{\left(2\right)}\left(\mathbf{k}\right),$$

(52a)

$$S_{CM}^{(2)}\left(\mathbf{k}\right)={\displaystyle \frac{1}{{\rho }_N}}{\mathsf{\Gamma }}_{CM}^{(2)}\left(\mathbf{k};\mathsf{\Psi }=0\right)={\displaystyle \frac{1}{〈N〉}}〈{\left|{\sum }_{j=1}^N\mathrm{e}\mathrm{x}\mathrm{p}\left(i\mathbf{k}·{\mathbf{R}}_{CM,j}\right)\right|}^2〉-{\left(2\pi \right)}^3{\rho }_N\delta \left(\mathbf{k}\right)=$$

(52b)

$$1+{\rho }_N\int d\mathbf{R}\mathrm{e}\mathrm{x}\mathrm{p}\left(i\mathbf{k}·\mathbf{R}\right)\left(g_{CM2}\left(R\right)-1\right),$$

(52c)

$${\delta \mathsf{\Gamma }}_{CM}^{\left(2\right)}\left({\mathbf{k}}_1{,\mathbf{k}}_2;{\mathsf{\Psi }}_F\right)\approx -\beta {\rho }_N\delta {\mathsf{\Psi }}_F\left({\mathbf{k}}_1{+\mathbf{k}}_2\right)S_{CM}^{\left(3\right)}\left({\mathbf{k}}_1{,\mathbf{k}}_2\right),$$

(53a)

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

(53b)

$${\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·{(\mathbf{R}}_{CM,j}-{\mathbf{R}}_{CM,m})+{\mathbf{k}}_2·({\mathbf{R}}_{CM,l}-{\mathbf{R}}_{CM,m})\right)\right]〉-\left\{\delta \left(\mathbf{k}\right)-\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{s}\right\}=$$

(53c)

$${\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·{\mathbf{R}}_{CM,j}+{\mathbf{k}}_2·{\mathbf{R}}_{CM,l}-\left({\mathbf{k}}_1+{\mathbf{k}}_2\right)·{\mathbf{R}}_{CM,m}\right)\right]〉-\left\{\delta \left(\mathbf{k}\right)-\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{s}\right\}.$$

(53d)

The same dependences on the moduli of the wavevectors (e.g., $k=\left|\mathbf{k}\right|),$ and similar comments as those made in connection to the classical Equations (10)-(11) do apply here. In particular, the $\delta \left(\mathbf{k}\right)-$terms are key to formulating completely $S_{CM}^{\left(2\right)}\left(\mathbf{k}\right)$ and $S_{CM}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)$, and correspond to the especial situations: (a) $\mathbf{k}=0$; (b) ${\mathbf{k}}_1=0$ and/or ${\mathbf{k}}_2=$ $0$, and ${\mathbf{k}}_1+{\mathbf{k}}_2=0.$ Explicitly, the CM3 triplet structure factor can be cast as:

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

$${\rho }_N^2\left({g_{CM}^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)-h}_{CM}^{\left(2\right)}\left(k_1\right)(2{\pi )}^3\delta \left({\mathbf{k}}_1+{\mathbf{k}}_2\right)-h_{CM}^{\left(2\right)}\left(k_2\right)(2{\pi )}^3\delta \left({\mathbf{k}}_1\right)-\right.$$

$$h_{CM}^{\left(2\right)}\left(k_1\right)(2{\pi )}^3\delta \left({\mathbf{k}}_2\right)-(2{\pi )}^3\delta \left({\mathbf{k}}_1\right)(2{\pi )}^3\left.\delta \left({\mathbf{k}}_2\right)\right),$$

(54)

where $g_{CM}^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)$ stands for the Fourier transform of $g_{CM3}\left(R_{12},R_{13},R_{23}\right),$ $h_{CM}^{\left(2\right)}\left(k_1\right)$ does for that of $h_{CM2}\left(R\right)$, etc. Given the analogies between the structural formulations for quantum centroids and classical monatomic particles, it would be a significant step forward if such parallelism could be extended to also include the OZn direct correlation function schemes for centroids. This is indeed the case as discussed in detail below [55,130].

5.1.2. PI-Centroid Direct Correlation Functions

There is the direct axiomatic way, consisting in checking that the OZn relationships satisfied in the classical domain by the set $\left\{c_n\right\}$ (see the insightful work by Lee [70]) are also satisfied by a set of PI-centroid functions $\left\{c_{CMn}\right\},$ which is defined accordingly in the same manner (see Reference [130] for OZ2). This procedure finds support in: (a) the formal equivalence between the PI-centroid and the classical hierarchies arising from their partition functions, and (b) the fact that the external field is switched off $({\mathsf{\Psi }}_F=0)$ when deriving the usual OZn equations [70], which serve to analyze the structures of the isolated fluid equilibrium states. Nonetheless, it seems useful to obtain the centroid OZn equations via a free-energy argument closely related to the standard procedures [171,172,173,174,175]. The reasoning below expands in a more comprehensive manner the early theoretical sketch reported in [55], and hinges heavily on the external field of constant force ${\mathsf{\Psi }}_F$ as the cause for the PI-centroid structures to show up. Recall that the PI-centroid (a “center-of-mass” variable) is an auxiliary mathematical object which, as such, cannot interact with external fields in an experimentally measurable direct way. Therefore, no claim can be made regarding any general validity of the centroid OZn($\mathsf{\Psi }$) derivations under other weak external potentials $\mathsf{\Psi }\ne {\mathsf{\Psi }}_F$, as is usual in the general functional OZn-context of classical fluids. The following development is an adaption to PI-centroids of the procedure that can be found in Haymet et al.’s insightful works [173,174,175].

(i) The starting point is the differential of the fluid energy:

$$dE=TdS-pd\mathcal{V}+\mu \mathsf{\delta }N+\int d\mathbf{q}{\rho }_N^{\left(1\right)}\left(\mathbf{q};{\mathsf{\Psi }}_F\right)\mathsf{\delta }{\mathsf{\Psi }}_F\left(\mathbf{q}\right),$$

(55)

where S is the entropy, p the pressure, ${\rho }_N^{\left(1\right)}\left(\mathbf{q};{\mathsf{\Psi }}_F\right)$ the one-body number density distribution of actual particles, and $\mathsf{\delta }{\mathsf{\Psi }}_F$ a (constant) variation in the field acting on such particles. The key point is that under these conditions Equation (55) can be rewritten in terms of the centroids:

$$dE=TdS-pd\mathcal{V}+\mu \mathsf{\delta }N+\int d\mathbf{R}{\rho }_{N,CM}^{\left(1\right)}\left(\mathbf{R};{\mathsf{\Psi }}_F\right)\mathsf{\delta }{\mathsf{\Psi }}_F\left(\mathbf{R}\right),$$

(56)

because it so happens that

$$-k_B\mathit{T\delta }\ln {\Xi }_{PI,D}({\mathsf{\Psi }}_F)=\int d\mathbf{q}{\rho }_N^{\left(1\right)}\left(\mathbf{q};{\mathsf{\Psi }}_F\right)\mathsf{\delta }{\mathsf{\Psi }}_F\left(\mathbf{q}\right)=\int d\mathbf{R}{\rho }_{N,CM}^{\left(1\right)}\left(\mathbf{R};{\mathsf{\Psi }}_F\right)\mathsf{\delta }{\mathsf{\Psi }}_F\left(\mathbf{R}\right).$$

(57)

Equation (57) can be obtained directly by integrating the functional derivative ${\displaystyle \frac{\delta \ln {\Xi }_{PI,D}({\mathsf{\Psi }}_F)}{\delta {\mathsf{\Psi }}_F\left(\mathbf{q}\right)}}$ (see Equation (48a) and subsection 5.2). However, at this point it seems better to consider the role of the normalized thermal packet $f_{inh.}$ representing a delocalized atom in the inhomogeneous fluid. The function $f_{inh.}$ is a probability density necessarily normalized to unity:

$$\int d\mathbf{q}{f}_{inh.}(\mathbf{q}-\mathbf{R};\mathbf{R};{\mathsf{\Psi }}_F)=1,$$

(58)

where $\mathbf{q}\left(x,y,z\right)$ is the position vector of the atom, and $\mathbf{R}(R_{\mathit{x}},R_y,R_z)$ the position vector of the centroid. To grasp this concept, consider the simpler example of the homogeneous Gaussian thermal packet provided by the semiclassical Feynman-Hibbs picture (GFH) [4,11,143,144]), in which a particle of mass m in a fluid is described as a thermal packet of width $\hslash /\sqrt{12m{k}_{B}T}$ [4]. The corresponding $f_{inh.}$ may be visualized as a deformation of the Gaussian packet adopting the symmetry of the external field. The derivation of Equation (57) uses the following facts: (a) the centroid is the “center of mass” of the probability density $f_{inh.}$; (b) the actual atom density ${\rho }_N^{\left(1\right)}\left(\mathbf{q};{\mathsf{\Psi }}_F\right)$ is obviously related to the centroid density ${\rho }_{N,CM}^{\left(1\right)}\left(\mathbf{R};{\mathsf{\Psi }}_F\right)$ through a smearing out operation involving the inhomogeneous $f_{inh.}$ packet; and (c) the external force undergoes a constant variation $\left(\mathsf{\delta }{f}_x,\mathsf{\delta }{f}_y,\mathsf{\delta }{f}_z\right).$ The proof goes as follows $(\mathsf{\xi }=\mathbf{q}-\mathbf{R})$:

$$\int d\mathbf{q}{\rho }_N^{\left(1\right)}\left(\mathbf{q};{\mathsf{\Psi }}_F\right)\mathsf{\delta }{\mathsf{\Psi }}_F\left(\mathbf{q}\right)=\int d\mathbf{q}\int \mathit{d}\mathbf{R}d\mathsf{\xi }{f_{inh.}(\mathsf{\xi };\mathbf{R};{\mathsf{\Psi }}_F)\rho }_{N,CM}^{\left(1\right)}\left(\mathbf{R};{\mathsf{\Psi }}_F\right)\mathsf{\delta }{\mathsf{\Psi }}_F\left(\mathbf{q}\right)\delta \left(\mathbf{q}-\mathbf{R}-\mathsf{\xi }\right)=$$

$$\int d\mathbf{R}{\rho }_{N,CM}^{\left(1\right)}\left(\mathbf{R};{\mathsf{\Psi }}_F\right)\int d\mathbf{q}{f}_{inh.}(\mathbf{q}-\mathbf{R};\mathbf{R};{\mathsf{\Psi }}_F)(x\mathsf{\delta }{f}_x+y\mathsf{\delta }{f}_y+z\mathsf{\delta }{f}_z)=$$

$$\int d\mathbf{R}{\rho }_{N,CM}^{\left(1\right)}\left(\mathbf{R};{\mathsf{\Psi }}_F\right)\left[R_x\mathsf{\delta }{f}_x\right.+R_y\mathsf{\delta }{f}_y+R_z\left.\mathsf{\delta }{f}_z\right]=\int d\mathbf{R}{\rho }_{N,CM}^{\left(1\right)}\left(\mathbf{R};{\mathsf{\Psi }}_F\right)\mathsf{\delta }{\mathsf{\Psi }}_F\left(\mathbf{R}\right).$$

(59)

(ii) Following the steps in References [173,174,175] one defines an auxiliary energy via a Legendre transformation:

$$\mathcal{U}=E-\int d\mathbf{R}{\rho }_{N,CM}^{\left(1\right)}\left(\mathbf{R};{\mathsf{\Psi }}_F\right){\mathsf{\Psi }}_F\left(\mathbf{R}\right),$$

(60)

and a Helmholtz free energy for the fluid in the form:

$$\mathcal{F}=\mathcal{U}-TS=-p\mathcal{V}+\int d\mathbf{R}\left(\mu -{\mathsf{\Psi }}_F\left(\mathbf{R}\right)\right){\rho }_{N,CM}^{\left(1\right)}\left(\mathbf{R};{\mathsf{\Psi }}_F\right).$$

(61)

(iii) A discussion of the selection of reference ideal systems in classical and quantum studies can be found in References [172,173,174,175]. This is a very important matter, for the application of functional differentiation to the corresponding excess free energy will yield a hierarchy of direct correlation functions. The direct correlation functions so obtained may or may not coincide with the usual for classical simple fluids [70,146,147,171,172,173,174,175]. In the current case of PI-centroids, which are quantum intermediate objects [4,11,18,19,140], given the classical-like form of the derivatives of the grand canonical potential, $-p\mathcal{V}={-k}_{B}T\ln {\mathsf{\Xi }}_{PI,D}\left({\mathsf{\Psi }}_F\right)$, with respect to the field ${\mathsf{\Psi }}_F$ given in Equations (48), the selection of the ideal Boltzmann system in the field ${\mathsf{\Psi }}_F$ seems the appropriate choice. By doing so, a classical-like hierarchy of PI-centroid direct correlation functions will be obtained. For the present structural purposes, focused on the usual Ornstein-Zernike framework at ${\mathsf{\Psi }}_F=0,$ this choice is going to be rewarding. Therefore, the Boltzmann-system reference [171] yields [55]:

$${\mathcal{F}}_{id}\left({\rho }_{N,CM}^{\left(1\right)}\left(\mathbf{R};{\mathsf{\Psi }}_{\mathrm{F}}\right)\right)=\int d\mathbf{R}{\rho }_{N,CM}^{\left(1\right)}\left(\mathbf{R};{\mathsf{\Psi }}_F\right){\beta }^{-1}\left\{\ln \left({\lambda }_B^3{\rho }_{N,CM}^{\left(1\right)}\left(\mathbf{R};{\mathsf{\Psi }}_F\right)\right)-1\right\},$$

(62)

and the excess free energy ${\mathcal{F}}_{exs}$ reads as:

$$-\beta (\mathcal{F}-{\mathcal{F}}_{id})=-\beta {\mathcal{F}}_{exs}\left({\rho }_{N,CM}^{\left(1\right)}\left(\mathbf{R};{\mathsf{\Psi }}_{\mathrm{F}}\right)\right)=\ln {\mathsf{\Xi }}_{PI,D}\left({\mathsf{\Psi }}_F\right)+\int d\mathbf{R}{\rho }_{N,CM}^{\left(1\right)}\left(\mathbf{R};{\mathsf{\Psi }}_F\right)\times$$

$$\left\{-\beta \mu +\beta {\mathsf{\Psi }}_F\left(\mathbf{R}\right)+\ln \left({\lambda }_B^3{\rho }_{N,CM}^{\left(1\right)}\left(\mathbf{R};{\mathsf{\Psi }}_F\right)\right)-1\right\},$$

(63)

(iv) The functional differentiation of ${\mathcal{F}}_{exs}$ leads to the hierarchy of the PI-centroid direct correlation functions $\left\{c_{CMn}\right\}$ which, up to the third order, reads as follows:

$$-{\displaystyle \frac{\delta \left(\beta {\mathcal{F}}_{exs}\right)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)}}=-\beta \mu +\beta {\mathsf{\Psi }}_F\left({\mathbf{R}}_1\right)+\ln \left({\lambda }_B^3{\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)\right)=c_{CM1}\left({\mathbf{R}}_1;{\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1);{\mathsf{\Psi }}_F({\mathbf{R}}_1\right)\right),$$

(64a)

$$-{\displaystyle \frac{{\delta }^2\left(\beta {\mathcal{F}}_{exs}\right)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)}}={\displaystyle \frac{\delta c_{CM1}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)}}=c_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_2;{\rho }_{N,CM}^{\left(1\right)}{,\mathsf{\Psi }}_F\right),$$

(64b)

$$-{\displaystyle \frac{{\delta }^3\left(\beta {\mathcal{F}}_{exs}\right)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_3;{\mathsf{\Psi }}_F\right)}}={\displaystyle \frac{\delta c_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_3;{\mathsf{\Psi }}_F\right)}}=c_{CM3}\left({\mathbf{R}}_1,{\mathbf{R}}_2,{\mathbf{R}}_3;{\rho }_{N,CM}^{\left(1\right)}{,\mathsf{\Psi }}_F\right).$$

(64c)

Of particular interest is the pair function $c_{CM2},$ which can be developed in the convenient form:

$$c_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_2;{\rho }_{N,CM}^{\left(1\right)}{,\mathsf{\Psi }}_F\right)={\displaystyle \frac{\delta c_{CM1}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)}}={\displaystyle \frac{\delta \left[-\beta \mu +\beta {\mathsf{\Psi }}_F\left({\mathbf{R}}_1\right)+\ln \left({\lambda }_B^3{\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)\right)\right]}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)}}=$$

$${\displaystyle \frac{\delta \left({\beta \mathsf{\Psi }}_F\left({\mathbf{R}}_1\right)\right)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)}}+{\displaystyle \frac{1}{{\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)}}{\displaystyle \frac{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)}}={\displaystyle \frac{\delta \left({\beta \mathsf{\Psi }}_F\left({\mathbf{R}}_1\right)\right)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)}}+{\displaystyle \frac{\delta {(\mathbf{R}}_1-{\mathbf{R}}_2)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)}},$$

(65)

where attention should be drawn to the inverse derivative ${\displaystyle \frac{\delta \left({\beta \mathsf{\Psi }}_F\left({\mathbf{R}}_1\right)\right)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)}},$ which is perfectly consistent with the use of the grand canonical ensemble, for density variations independent of the external field can be made [147].

5.1.3. OZ2 and OZ3 Frameworks

The OZ2 framework for PI-centroids arises in the usual way by applying the matrix identity [5,147]

$$\int d{\mathbf{R}}_3{\displaystyle \frac{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)}{\delta \left({-\beta \mathsf{\Psi }}_F\left({\mathbf{R}}_3\right)\right)}}{\displaystyle \frac{\delta \left(-{\beta \mathsf{\Psi }}_F\left({\mathbf{R}}_3\right)\right)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)}}=\delta {(\mathbf{R}}_1-{\mathbf{R}}_2).$$

(66)

Inserting Equations (48b) and (65) into Equation (66), one arrives at the inhomogeneous OZ2 equation in the field ${\mathsf{\Psi }}_F$ ${(h}_{CM2}\left({\mathsf{\Psi }}_F\right)=g_{CM2}\left({\mathsf{\Psi }}_F\right)-1):$

$$h_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)=c_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)+\int d{\mathbf{R}}_3{\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_3;{\mathsf{\Psi }}_F\right)h_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_3;{\mathsf{\Psi }}_F\right)c_{CM2}\left({\mathbf{R}}_3,{\mathbf{R}}_2;{\mathsf{\Psi }}_F\right),$$

(67)

which transforms into the conventional homogeneous version if the field is switched off $\left({\mathsf{\Psi }}_F=0\right):$

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

(68)

By Fourier transforming Equation (68) and combining the result with Equation (52c), the static pair structure factor for PI-centroids is:

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

(69)

formally identical to the classical form Equation (13a). It can be demonstrated that for the quantum fluid [130]

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

(70)

a result that should not be surprising, since the fluctuations in the number of particles can be counted with the centroids. Therefore, the PI-centroid $\mathbf{k}=0$ component fixed via $c_{CM}^{\left(2\right)}(\mathbf{k}=0)$ formally provides an exact way to obtain the equation of state of a fluid with quantum diffraction behavior. As in the classical case, $c_{CM2}\left(R_{12}\right)$ is expected to decay effectively to zero within a short range of distances, which is indeed the case and leads to highly accurate results in k-space [46,62]. In this connection, recall that direct correlation functions show generally a tail in their decays [149], but its importance is negligible beyond a certain cut off (see References [130,136] and Figure 2).

At this point, and following Reference [41], it is worth introducing a new hierarchy $\left\{{\mathsf{{\rm K}}}_{CMn}\right\}$ arising from the functional derivatives of the total free energy $\mathcal{F}$. As before, the following formulas are limited to the third order, the iteration to higher orders being straightforward. For brevity, detailed reference to the various dependences of the functionals are omitted in what follows. At first order one finds:

$${\displaystyle \frac{\delta (\beta \mathcal{F})}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)}}=\beta \mu -\beta {\mathsf{\Psi }}_F\left({\mathbf{R}}_1\right)={\mathsf{{\rm K}}}_{CM1}({\mathbf{R}}_1;{\mathsf{\Psi }}_F)$$

(71)

$${\displaystyle \frac{\delta {\mathsf{{\rm K}}}_{CM1}({\mathbf{R}}_1;{\mathsf{\Psi }}_F)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)}}={\displaystyle \frac{\delta {(\mathbf{R}}_1-{\mathbf{R}}_2)}{{\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)}}-c_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_2;{\rho }_{N,CM}^{\left(1\right)},{\mathsf{\Psi }}_F\right)={\mathsf{{\rm K}}}_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_2;{\mathsf{\Psi }}_F\right),$$

(72)

$${\displaystyle \frac{\mathit{\delta }{\mathsf{{\rm K}}}_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_3;{\mathsf{\Psi }}_F\right)}}=-{{\displaystyle \frac{\delta {(\mathbf{R}}_1-{\mathbf{R}}_2)\delta {(\mathbf{R}}_1-{\mathbf{R}}_3)}{{\left[{\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)\right]}^2}}-c}_{CM3}\left({\mathbf{R}}_1,{\mathbf{R}}_2,{\mathbf{R}}_3;{\rho }_{N,CM}^{\left(1\right)}{,\mathsf{\Psi }}_F\right)={\mathsf{{\rm K}}}_{CM3}\left({\mathbf{R}}_1,{\mathbf{R}}_2,{\mathbf{R}}_3;{\mathsf{\Psi }}_F\right).$$

(73)

where the terms containing $\delta {(\mathbf{R}}_j-{\mathbf{R}}_l)$ come from the ideal system contributions. In Equation (71) one identifies the actual chemical potential $\mu$ as the algebraic sum of the external field and an intrinsic contribution, in agreement with the expected constancy of the chemical potential in the presence of the external field [171,172].

The interest of the $\left\{{\mathsf{{\rm K}}}_{CMn}\right\}$ hierarchy lies in the fact that, together with the hierarchy $\left\{{\mathsf{\Gamma }}_{CMn}\right\},$ they allow the OZn frameworks to be formulated in compact ways [41]. Thus, OZ2 as given in Equation (67) can be cast as:

$$\int d{\mathbf{R}}_3{\mathsf{{\rm K}}}_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_3;{\mathsf{\Psi }}_F\right){\mathsf{\Gamma }}_{CM2}\left({\mathbf{R}}_2,{\mathbf{R}}_3;{\mathsf{\Psi }}_F\right)=\delta {(\mathbf{R}}_1-{\mathbf{R}}_2),$$

(74)

which reduces to Equation (68) in the limit $\mathsf{\Psi }\to 0.$ Furthermore, via functional differentiation of Equation (74) one obtains the OZ3 framework. One of the four equivalent OZ3 equations can be derived by following closely References [41] and [70]. First, by taking the functional derivative with respect to the one-particle density a new identity arises:

$${\displaystyle \frac{\delta }{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_4;{\mathsf{\Psi }}_F\right)}}\int d{\mathbf{R}}_3{\mathsf{{\rm K}}}_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_3;{\mathsf{\Psi }}_F\right){\mathsf{\Gamma }}_{CM2}\left({\mathbf{R}}_2,{\mathbf{R}}_3;{\mathsf{\Psi }}_F\right)=$$

$$\int d{\mathbf{R}}_3{\displaystyle \frac{\delta {\mathsf{{\rm K}}}_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_3;{\mathsf{\Psi }}_F\right)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_4;{\mathsf{\Psi }}_F\right)}}{\mathsf{\Gamma }}_{CM2}\left({\mathbf{R}}_2,{\mathbf{R}}_3;{\mathsf{\Psi }}_F\right)+\int d{\mathbf{R}}_3{\mathsf{{\rm K}}}_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_3;{\mathsf{\Psi }}_F\right){\displaystyle \frac{\delta {\mathsf{\Gamma }}_{CM2}\left({\mathbf{R}}_2,{\mathbf{R}}_3;{\mathsf{\Psi }}_F\right)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_4;{\mathsf{\Psi }}_F\right)}}=0,$$

(75)

Second, with application to ${\displaystyle \frac{\delta {\mathsf{\Gamma }}_{CM2}}{\delta {\rho }_{N,CM}^{\left(1\right)}}}$ of intermediate differentiation with respect to $\delta {\mathsf{\Psi }}_F$, and renumbering the particle labels, one arrives at the final form:

$$\int d{\mathbf{R}}_3{{\mathsf{{\rm K}}}_{CM3}\left({\mathbf{R}}_1,{{\mathbf{R}}_2,\mathbf{R}}_3;{\mathsf{\Psi }}_F\right)\mathsf{\Gamma }}_{CM2}\left({\mathbf{R}}_3,{\mathbf{R}}_4;{\mathsf{\Psi }}_F\right)+$$

$$\int d{\mathbf{R}}_{3}d{\mathbf{R}}_5{\mathsf{{\rm K}}}_{CM2}\left({\mathbf{R}}_2,{\mathbf{R}}_3;{\mathsf{\Psi }}_F\right){{\mathsf{\Gamma }}_{CM3}\left({\mathbf{R}}_3,{\mathbf{R}}_4,{\mathbf{R}}_5;{\mathsf{\Psi }}_F\right)\mathsf{{\rm K}}}_{CM2}\left({\mathbf{R}}_5,{\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)=0,$$

(76)

which, again, reduces to the homogeneous OZ3 equation in the limit ${\mathsf{\Psi }}_F\to 0.$The triplet linear response function for PI-centroids of the homogeneous quantum fluid $({\mathsf{\Psi }}_F=0)$ arises from the Fourier transform of Equation (76) combined with the definition given in Equation (53b). Thus, by noting that:

$${\mathsf{{\rm K}}}_{CM}^{\left(2\right)}\left(\mathbf{k}\right){\mathsf{\Gamma }}_{CM}^{\left(2\right)}\left(\mathbf{k}\right)=1,$$

(77a)

one finds the triplet structure factor:

$${\rho }_N{S}_{CM}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)={\mathsf{\Gamma }}_{CM}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)=-{\displaystyle \frac{{\mathsf{{\rm K}}}_{CM}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right){\mathsf{\Gamma }}_{CM}^{\left(2\right)}\left(\left|\left({\mathbf{k}}_1+{\mathbf{k}}_2\right)\right|\right)}{{\mathsf{{\rm K}}}_{CM}^{\left(2\right)}\left({\mathbf{k}}_1\right){\mathsf{{\rm K}}}_{CM}^{\left(2\right)}\left({\mathbf{k}}_2\right)}},$$

(77b)

$$S_{CM}^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)=S_{CM}^{\left(2\right)}\left({\mathbf{k}}_1\right)S_{CM}^{\left(2\right)}({\mathbf{k}}_2)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)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)\right\},$$

(77c)

which is equivalent to Equation (54) but formulated in terms of direct correlation functions [55], with all the advantages and problems associated with these objects commented earlier. In this connection, if the variations in the density are uniform, the general hierarchy $\left\{c_{CMn}\right\}$ sketched in Equations (64) takes the simpler and operative Baxter’s form [68]. In particular, the analogous to Equations (14) read as:

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

(78a)

$${\displaystyle \frac{\partial c_{CM}^{(2)}\left(k=0\right)}{\partial {\rho }_N}}{=c}_{CM}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2=0\right);T=\mathrm{constant}.$$

(78b)

5.1.4. Some Additional Relationships and Further Remarks

For the hierarchy $\left\{c_{CMn}\right\}$ the partition function ${\mathsf{\Xi }}_{PI,D}({\mathsf{\Psi }}_F)$ leads to all the relationships that are applicable in the classical context [70]. A straightforward calculation yields:

$${\displaystyle \frac{{\delta }^s\ln {\mathsf{\Xi }}_{PI,D}({\mathsf{\Psi }}_F)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)\dots \delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_s;{\mathsf{\Psi }}_F\right)}}=-\left(s-1\right)c_{CMs}\left({\mathbf{R}}_{1,}{\mathbf{R}}_2{,\dots ,\mathbf{R}}_s;{\mathsf{\Psi }}_F\right)-$$

$$\int d{\mathbf{R}}_{s+1}{\rho }_N^{\left(1\right)}\left({\mathbf{R}}_{s+1};{\mathsf{\Psi }}_F\right)c_{CMs+1}\left({\mathbf{R}}_{1,}{\mathbf{R}}_2{,\dots ,\mathbf{R}}_s,{\mathbf{R}}_{s+1};{\mathsf{\Psi }}_F\right);s\ge 2,$$

(79a)

which at zero field gives the corresponding centroid equations for the isolated quantum fluid. Completing with the first order derivative $(s=1)$, the first three of the related equations read as:

$${\left.{\displaystyle \frac{\delta \ln {\mathsf{\Xi }}_{PI,D}({\mathsf{\Psi }}_F)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)}}\right|}_{{\mathsf{\Psi }}_F=0}=1-{\rho }_N\int d{\mathbf{R}}_2{c}_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_2\right),$$

(79b)

$${\left.{\displaystyle \frac{{\delta }^2\ln {\mathsf{\Xi }}_{PI,D}({\mathsf{\Psi }}_F)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)}}\right|}_{{\mathsf{\Psi }}_F=0}=-c_{CM2}\left({\mathbf{R}}_1,{\mathbf{R}}_2\right)-{\rho }_N\int d{\mathbf{R}}_3{c}_{CM3}\left({\mathbf{R}}_1,{\mathbf{R}}_2,{\mathbf{R}}_3\right),$$

(79c)

$${\left.{\displaystyle \frac{{\delta }^3\ln {\mathsf{\Xi }}_{PI,D}({\mathsf{\Psi }}_F)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_2;{\mathsf{\Psi }}_F\right)\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_3;{\mathsf{\Psi }}_F\right)}}\right|}_{{\mathsf{\Psi }}_F=0}=-2c_{CM3}\left({\mathbf{R}}_1,{\mathbf{R}}_2,{\mathbf{R}}_3\right)-$$

$${\rho }_N\int d{\mathbf{R}}_4{c}_{CM4}\left({\mathbf{R}}_1,{\mathbf{R}}_2,{\mathbf{R}}_3,{\mathbf{R}}_4\right).$$

(79d)

Besides, note the classical-like additional identities that hold in the PI-centroid context:

$${\displaystyle \frac{\delta \ln {\mathsf{\Xi }}_{PI,D}({\mathsf{\Psi }}_F)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)}}+\int d{\mathbf{R}}_2{\rho }_{N,CM}^{\left(1\right)}({\mathbf{R}}_2;{\mathsf{\Psi }}_F){\displaystyle \frac{\delta \left(\beta {\mathsf{\Psi }}_F\left({\mathbf{R}}_2\right)\right)}{\delta {\rho }_{N,CM}^{\left(1\right)}\left({\mathbf{R}}_1;{\mathsf{\Psi }}_F\right)}}\equiv 0,$$

(80a)

$$\int d{\mathbf{R}}_{1}d{\mathbf{R}}_2{\displaystyle \frac{{\delta }^2\ln Z_{N,QD}^{PI}({\mathsf{\Psi }}_F)}{\delta {\mathsf{\Psi }}_F({\mathbf{R}}_1)\delta {\mathsf{\Psi }}_F({\mathbf{R}}_2)}}\equiv 0;{\mathsf{\Psi }}_F\ne 0,$$

(80b)

$$\int d{\mathbf{R}}_{1}d{\mathbf{R}}_{2}d{\mathbf{R}}_3{\displaystyle \frac{{\delta }^3\ln Z_{N,QD}^{PI}({\mathsf{\Psi }}_F)}{\delta {\mathsf{\Psi }}_F({\mathbf{R}}_1)\delta {\mathsf{\Psi }}_F({\mathbf{R}}_2)\delta {\mathsf{\Psi }}_F({\mathbf{R}}_3)}}\equiv 0;{\mathsf{\Psi }}_F\ne 0,$$

(80c)

$${\left(-k_{B}T\right)}^2\int d{\mathbf{R}}_{1}d{\mathbf{R}}_2{\displaystyle \frac{{\delta }^2\ln {\mathsf{\Xi }}_{PI,D}({\mathsf{\Psi }}_F)}{\delta {\mathsf{\Psi }}_F({\mathbf{R}}_1)\delta {\mathsf{\Psi }}_F({\mathbf{R}}_2)}}=〈N^2〉-{〈N〉}^2;{\mathsf{\Psi }}_F\ne 0,$$

(80d)

$${\left(-k_{B}T\right)}^3\int d{\mathbf{R}}_{1}d{\mathbf{R}}_{2}d{\mathbf{R}}_3{\displaystyle \frac{{\delta }^3\ln {\mathsf{\Xi }}_{PI,D}({\mathsf{\Psi }}_F)}{\delta {\mathsf{\Psi }}_F({\mathbf{R}}_1)\delta {\mathsf{\Psi }}_F({\mathbf{R}}_2)\delta {\mathsf{\Psi }}_F({\mathbf{R}}_3)}}=〈N^3〉-{3〈N^2〉〈N〉+2〈N〉}^3;{\mathsf{\Psi }}_F\ne 0,$$

(80e)

where $Z_{N,QD}^{PI}({\mathsf{\Psi }}_F)$ is the in-the-field PI canonical partition function, and the relations to the number fluctuations are clearly displayed. These number fluctuations refer to the in-the-field fluid; the formulas on the right-hand sides remain the same when ${\mathsf{\Psi }}_F=0,$but obviously the meaning is not the same.

Finally, several remarks are in order. It is worthwhile to point out that the formal OZn derivations applied in the classical domain and related to the bijective connection external field$\leftrightarrow$equilibrium density [171,172,175,176,177], hold here for ${\mathsf{\Psi }}_F\leftrightarrow {\rho }_{N,CM}^{(1)}$. This is guaranteed by the form of the quantum partition function Equation (44) which is isomorphic to its classical counterpart Equation (2). Another point is that there is no defined interparticle potential $v(R_{jl})$ between PI centroids (compare to the variational approximations arising from the PI formalism [4,140,141,142], which yield versions of such a potential), but this is not relevant to this discussion [171]. In this regard, nor is the fact that under the same conditions two different densities are included in Equation (57), namely the true particle density and the centroid density, because: (a) one is the actual and the other is an auxiliary quantity, (b) each corresponds to a distinct object, and (c) both are intimately related. (The situation is analogous to a distribution of mass and the consideration of its center of gravity). To bring this subsection to a close, note that the applicability of the classical OZn numerical schemes to the quantum CMn class is perfectly possible and rewarding: (i) at the pair level, OZ2(BDH+BHw) is a powerful way to fix the equations of state of quantum fluids [46,62,130,136,137]; and (ii) at the triplet level, first OZ3 applications were reported in References [55,56,58].

5.2. The PI Total Continuous Linear Response TLRn Class

For this class the fluid undergoes a general external weak field $\mathsf{\Psi }$ which is continuous either spatially or in imaginary time (no magnetic interactions are considered). The PI partition function can be cast as [7,11,57]:

$${\mathsf{\Xi }}_{PI,D}\left(\mathsf{\Psi }\right)={\sum }_{N\ge 0}{\displaystyle \frac{{\mathfrak{z}}^N}{N!}}\int {\prod }_{j=1}^N{\prod }_{t=1}^{P}d{\mathbf{r}}_j^t\times \mathrm{e}\mathrm{x}\mathrm{p}\left[-\beta W_{NP}\left({\mathbf{r}}^{NP}\right)\right]\times \mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{\beta }{X}}{\sum }_{j=1}^N{\sum }_{\tau =1}^X{\mathsf{\Psi }(\mathbf{r}}_j^{\tau })\right).$$

(81)

The functional derivatives for this model now can be cast as:

$$-k_{B}TX{\displaystyle \frac{\delta \mathrm{l}\mathrm{n}{\mathsf{\Xi }}_{PI,D}\left(\mathsf{\Psi }\right)}{\delta \mathsf{\Psi }\left({\mathbf{q}}_1\right)}}={\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_1;\mathsf{\Psi }\right),$$

(82a)

$${\left(-k_{B}TX\right)}^2{\displaystyle \frac{{\delta }^2\mathrm{l}\mathrm{n}{\mathsf{\Xi }}_{PI,D}\left(\mathsf{\Psi }\right)}{\delta \mathsf{\Psi }\left({\mathbf{q}}_1\right)\delta \mathsf{\Psi }\left({\mathbf{q}}_2\right)}}={\mathsf{\Gamma }}_{TLR2}\left({\mathbf{q}}_1,{\mathbf{q}}_2;\mathsf{\Psi }\right)={\rho }_{NX}^{\left(2\right)}\left({\mathbf{q}}_1,{\mathbf{q}}_2;\mathsf{\Psi }\right)-{\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_1;\mathsf{\Psi }\right){\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_2;\mathsf{\Psi }\right)+$$

$${\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_1;\mathsf{\Psi }\right)\delta \left({\mathbf{q}}_1-{\mathbf{q}}_2\right),$$

(82b)

$${\left(-k_{B}TX\right)}^3{\displaystyle \frac{{\delta }^3\mathrm{l}\mathrm{n}{\mathsf{\Xi }}_{PI,D}\left(\mathsf{\Psi }\right)}{\delta \mathsf{\Psi }\left({\mathbf{q}}_1\right)\delta \mathsf{\Psi }\left({\mathbf{q}}_2\right)\delta \mathsf{\Psi }\left({\mathbf{q}}_3\right)}}={\mathsf{\Gamma }}_{TLR3}\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3;\mathsf{\Psi }\right)=$$

$${\rho }_{NX}^{\left(3\right)}\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3;\mathsf{\Psi }\right)-{\rho }_{NX}^{\left(2\right)}\left({\mathbf{q}}_1,{\mathbf{q}}_2;\mathsf{\Psi }\right){\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_3;\mathsf{\Psi }\right)-$$

$${\rho }_{NX}^{\left(2\right)}\left({\mathbf{q}}_1,{\mathbf{q}}_3;\mathsf{\Psi }\right){\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_2;\mathsf{\Psi }\right)-{\rho }_{NX}^{\left(2\right)}\left({\mathbf{q}}_2,{\mathbf{q}}_3;\mathsf{\Psi }\right){\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_1;\mathsf{\Psi }\right)+$$

$$2{\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_1;\mathsf{\Psi }\right){\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_2;\mathsf{\Psi }\right){\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_3;\mathsf{\Psi }\right)+$$

$${\rho }_{NX}^{\left(2\right)}\left({\mathbf{q}}_1,{\mathbf{q}}_2;\mathsf{\Psi }\right)\delta \left({\mathbf{q}}_1-{\mathbf{q}}_3\right)-{\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_1;\mathsf{\Psi }\right){\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_2;\mathsf{\Psi }\right)\delta \left({\mathbf{q}}_1-{\mathbf{q}}_3\right)+$$

$${\rho }_{NX}^{\left(2\right)}\left({\mathbf{q}}_1,{\mathbf{q}}_3;\mathsf{\Psi }\right)\delta \left({\mathbf{q}}_3-{\mathbf{q}}_2\right)-{\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_1;\mathsf{\Psi }\right){\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_3;\mathsf{\Psi }\right)\delta \left({\mathbf{q}}_3-{\mathbf{q}}_2\right)+$$

$${\rho }_{NX}^{\left(2\right)}\left({\mathbf{q}}_2,{\mathbf{q}}_3;\mathsf{\Psi }\right)\delta \left({\mathbf{q}}_2-{\mathbf{q}}_1\right)-{\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_2;\mathsf{\Psi }\right){\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_3;\mathsf{\Psi }\right)\delta \left({\mathbf{q}}_2-{\mathbf{q}}_1\right)+$$

$${\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_1;\mathsf{\Psi }\right)\delta \left({\mathbf{q}}_1-{\mathbf{q}}_3\right)\delta \left({\mathbf{q}}_2-{\mathbf{q}}_3\right).$$

(82c)

In the foregoing formulas a vector ${\mathbf{q}}_i$ stands for the position vector of a structurally significant bead in the grand-canonical sample (i.e., ${\mathbf{q}}_i$ is any of the ${\mathbf{r}}_j^{\tau }).$Recall that for structural evaluations each canonical ensemble has $NX$ significant particles, and that X is held fixed throughout the grand canonical ensemble. Therefore, by following the standard procedures [2,11,57] one sets the definitions:

$${\rho }_{NX}^{\left(1\right)}\left({\mathbf{q}}_1;\mathsf{\Psi }\right)=(X{\rho }_N)G_{TLR1}\left({\mathbf{q}}_1;\mathsf{\Psi }\right),$$

(83a)

$${\rho }_{NX}^{\left(2\right)}\left({\mathbf{q}}_1,{\mathbf{q}}_2;\mathsf{\Psi }\right)={\left(X{\rho }_N\right)}^2{G}_{TLR2}\left({\mathbf{q}}_1,{\mathbf{q}}_2,;\mathsf{\Psi }\right),$$

(83b)

$${\rho }_{NX}^{\left(3\right)}\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3;\mathsf{\Psi }\right)={\left(X{\rho }_N\right)}^3{G}_{TLR3}\left({\mathbf{q}}_1{\mathbf{q}}_2,{\mathbf{q}}_3;\mathsf{\Psi }\right).$$

(83c)

and obtains the relationships:

$$\int d{\mathbf{q}}_{1}d{\mathbf{q}}_2{\displaystyle \frac{{\delta }^2\ln Z_{N,\mathit{QD}}^{\mathit{PI}}(\mathsf{\Psi })}{\delta \mathsf{\Psi }({\mathbf{q}}_1)\delta \mathsf{\Psi }({\mathbf{q}}_2)}}\equiv 0;\mathsf{\Psi }\ne 0,$$

(83d)

$$\int d{\mathbf{q}}_{1}d{\mathbf{q}}_{2}d{\mathbf{q}}_3{\displaystyle \frac{{\delta }^3\ln Z_{N,\mathit{QD}}^{\mathit{PI}}(\mathsf{\Psi })}{\delta \mathsf{\Psi }({\mathbf{q}}_1)\delta \mathsf{\Psi }({\mathbf{q}}_2)\delta \mathsf{\Psi }({\mathbf{q}}_3)}}\equiv 0;\mathsf{\Psi }\ne 0,$$

(83e)

$${\left(-k_{B}T\right)}^2\int d{\mathbf{q}}_{1}d{\mathbf{q}}_2{\displaystyle \frac{{\delta }^2\ln {\mathsf{\Xi }}_{PI,D}(\mathsf{\Psi })}{\delta \mathsf{\Psi }({\mathbf{q}}_1)\delta \mathsf{\Psi }({\mathbf{q}}_2)}}=〈N^2〉-{〈N〉}^2;\mathsf{\Psi }\ne 0,$$

(83f)

$${\left(-k_{B}T\right)}^3\int d{\mathbf{q}}_{1}d{\mathbf{q}}_{2}d{\mathbf{q}}_3{\displaystyle \frac{{\delta }^3\ln {\mathsf{\Xi }}_{PI,D}(\mathsf{\Psi })}{\delta \mathsf{\Psi }({\mathbf{q}}_1)\delta \mathsf{\Psi }({\mathbf{q}}_2)\delta \mathsf{\Psi }({\mathbf{q}}_3)}}=〈N^3〉-{3〈N^2〉〈N〉+2〈N〉}^3;\mathsf{\Psi }\ne 0,$$

(83g)

where the number fluctuations refer to the in-the-field fluid. The formulas on the right-hand sides remain the same when $\mathsf{\Psi }=0,$although the meaning is therefore not the same.

The linear response cases in this case arise by making $\mathsf{\Psi }=0$ on the right-hand sides of Equations (82). In r-space this can be summarized in the following equations [57]:

$$-k_{B}T{\displaystyle \frac{\delta \mathrm{l}\mathrm{n}{\mathsf{\Xi }}_{PI,D}\left(\mathsf{\Psi }\right)}{\delta \mathsf{\Psi }\left({\mathbf{q}}_1\right)}}\approx {\rho }_N^{\left(1\right)}\left({\mathbf{q}}_1;\mathsf{\Psi }=0\right)={\rho }_N,$$

(84a)

$$-k_{B}T{\displaystyle \frac{\delta {\rho }_N^{\left(1\right)}\left({\mathbf{q}}_1;\mathsf{\Psi }\right)}{\delta \mathsf{\Psi }\left({\mathbf{q}}_2\right)}}\approx {\displaystyle \frac{1}{X^2}}{\mathsf{\Gamma }}_{\mathit{TLR}2}\left({\mathbf{q}}_1,{\mathbf{q}}_2;\mathsf{\Psi }=0\right)={\rho }_N^2\left(G_{TLR2}\left({\mathbf{q}}_1,{\mathbf{q}}_2\right)-1\right)+{\rho }_N{\displaystyle \frac{\delta \left({\mathbf{q}}_2-{\mathbf{q}}_1\right)}{X}},$$

(84b)

$$-k_{B}T{\displaystyle \frac{\delta {\mathsf{\Gamma }}_{TLR2}\left({\mathbf{q}}_1,{\mathbf{q}}_2;\mathsf{\Psi }\right)}{\delta \mathsf{\Psi }\left({\mathbf{q}}_3\right)}}\approx {\displaystyle \frac{1}{X^3}}{\mathsf{\Gamma }}_{\mathit{TLR}3}\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3;\mathsf{\Psi }=0\right)=$$

$${\rho }_N^3\left[G_{TLR3}\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3\right)-G_{TLR2}\left({\mathbf{q}}_1,{\mathbf{q}}_2\right)-G_{TLR2}\left({\mathbf{q}}_1,{\mathbf{q}}_3\right)-G_{TLR2}\left({\mathbf{q}}_2,{\mathbf{q}}_3\right)+2\right]+$$

$${\displaystyle \frac{1}{X}}{\rho }_N^2\left[\left(G_{TLR2}\left({\mathbf{q}}_1,{\mathbf{q}}_2\right)-1\right)\delta \left({\mathbf{q}}_1-{\mathbf{q}}_3\right)+\left(G_{TLR2}\left({\mathbf{q}}_1,{\mathbf{q}}_3\right)-1\right)\delta \left({\mathbf{q}}_3-{\mathbf{q}}_2\right)\right.+$$

$$\left.\left(G_{TLR2}\left({\mathbf{q}}_2,{\mathbf{q}}_3\right)-1\right)\delta \left({\mathbf{q}}_2-{\mathbf{q}}_1\right)\right]+{\displaystyle \frac{1}{X^2}}{\rho }_N\delta \left({\mathbf{q}}_1-{\mathbf{q}}_3\right)\delta \left({\mathbf{q}}_2-{\mathbf{q}}_3\right),$$

(84c)

where it is important to note the role played by X in Equations (84). Equations (84b) and (84c) can be formulated in concise forms with Equations (82) and the spontaneous fluctuations of the bead density at the two- and three-body levels respectively.

For notational simplicity, by ${\tau }^{*}-$arranging correlatively the NX beads present in each canonical ensemble$,$the pair $G_{TLR2}$and triplet $G_{TLR3}$ correlation functions are defined by the ensemble averages [57]:

$${\rho }_N^2{G}_{TLR2}\left({\mathbf{q}}_1,{\mathbf{q}}_2\right)={\rho }_N^2{G}_{TLR2}\left(r_{12}\right)={\displaystyle \frac{1}{X^2}}〈{\sum }_{{\tau }_1^{\mathrm{*}}\ne {\tau }_2^{\mathrm{*}}}\delta ({\mathbf{r}}^{{\tau }_1^{\mathrm{*}}}-{\mathbf{q}}_1)\delta ({\mathbf{r}}^{{\tau }_2^{\mathrm{*}}}-{\mathbf{q}}_2)〉,$$

(85a)

$${\rho }_N^3{G}_{TLR3}\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3\right)={\rho }_N^3{G}_{TLR3}\left(r_{12},r_{13},r_{23}\right)=$$

$${\displaystyle \frac{1}{X^3}}〈{\sum }_{{\tau }_1^{\mathrm{*}}\ne {\tau }_2^{\mathrm{*}}\ne {\tau }_3^{\mathrm{*}}\ne {\tau }_1^{\mathrm{*}}}\delta ({\mathbf{r}}^{{\tau }_1^{\mathrm{*}}}-{\mathbf{q}}_1)\delta ({\mathbf{r}}^{{\tau }_2^{\mathrm{*}}}-{\mathbf{q}}_2)\delta ({\mathbf{r}}^{{\tau }_3^{\mathrm{*}}}-{\mathbf{q}}_3)〉.$$

(85b)

Note that the correlations between beads include both: equal and different imaginary times. The foregoing averages (85a)-(85b) can be split into actual atom components: one-particle, two-particle, and three-particle. In this regard, (a) for $G_{TLR2}$: the bead-bead pair correlations can correspond to beads in the same necklace (one-atom self-correlations $s_{SC\mathrm{1,2}}$) or to beads in different necklaces (two-atom correlations $g_{LR2}$) [11}:

$$G_{TLR2}\left({\mathbf{q}}_1,{\mathbf{q}}_2\right)=s_{SC\mathrm{1,2}}\left(r_{12}\right)+g_{LR2}\left(r_{12}\right);$$

(86a)

and (b) for $G_{TLR3}:$ the bead-bead-bead triplet correlations can correspond to three beads in the same necklace (one-atom self correlations $s_{SC\mathrm{1,3}}$), two beads in the same necklace and the third in another necklace (mixture of one-atom self-correlations and two-atom correlations ${\mathsf{\Sigma }}_{LR\mathrm{2,3}}$), and the three beads in different necklaces (three-atom correlations $g_{LR3})$. At the triplet level this can be cast as [11,57]:

$$G_{TLR3}\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3\right)=s_{SC\mathrm{1,3}}\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3\right)+{\mathsf{\Sigma }}_{LR\mathrm{2,3}}({\mathbf{q}}_1,{\mathbf{q}}_2|{\mathbf{q}}_3)+{\mathsf{\Sigma }}_{LR\mathrm{2,3}}({\mathbf{q}}_1,{\mathbf{q}}_3|{\mathbf{q}}_2)+$$

$${\mathsf{\Sigma }}_{LR\mathrm{2,3}}({\mathbf{q}}_2,{\mathbf{q}}_3|{\mathbf{q}}_1)+g_{LR3}\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3\right),$$

(86b)

where for clarity reasons the vector notation and vertical bars to separate necklaces in the $\mathsf{\Sigma }$ cases are used.

The pair and triplet static structure factors arise from Equations (84b)-(84c) via the Fourier transform process explained earlier. They can be cast as [11,55,56,57]:

$$S_{TLR}^{\left(2\right)}\left(\mathbf{k}\right)={\displaystyle \frac{1}{{\rho }_N}}{\mathsf{\Gamma }}_{TLR}^{\left(2\right)}\left(\mathbf{k};\mathsf{\Psi }=0\right)={\displaystyle \frac{1}{〈N〉X^2}}〈{\left|\sum _{{\mathit{\tau }}^{\mathbf{*}}=1}^{NX}\mathrm{e}\mathrm{x}\mathrm{p}\left(i\mathbf{k}·{\mathbf{r}}^{{\mathit{\tau }}^{\mathbf{*}}}\right)\right|}^2〉-{\left(2\pi \right)}^3{\rho }_N\delta \left(\mathbf{k}\right)=$$

$${\displaystyle \frac{1}{X}}+{\rho }_N\int d\mathbf{r}\mathrm{e}\mathrm{x}\mathrm{p}\left(i\mathbf{k}·\mathbf{r}\right)\left(G_{TLR2}\left(r\right)-1\right),$$

(87a)

$$S_{TLR}^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)={\displaystyle \frac{1}{{\rho }_N}}{\mathsf{\Gamma }}_{TLR}^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2;\mathsf{\Psi }=0\right)=$$

$${\displaystyle \frac{1}{〈N〉X^3}}〈{\sum }_{{\tau }_1^{\mathrm{*}}=1}^{NX}{\sum }_{{\tau }_2^{\mathrm{*}}=1}^{NX}{\sum }_{{\tau }_3^{\mathrm{*}}=1}^{NX}\mathrm{e}\mathrm{x}\mathrm{p}\left[i\left({\mathbf{k}}_1·{\mathbf{r}}^{{\mathit{\tau }}_1^{\mathbf{*}}}+{\mathbf{k}}_2·{\mathbf{r}}^{{\mathit{\tau }}_2^{\mathbf{*}}}-\left({\mathbf{k}}_1+{\mathbf{k}}_2\right)·{\mathbf{r}}^{{\mathit{\tau }}_3^{\mathbf{*}}}\right)\right]〉-\left\{\delta \left(\mathbf{k}\right)-\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{s}\right\},$$

(87b)

which in expanded form can be written as:

$$S_{TLR}^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)={\displaystyle \frac{1}{X^2}}+{\displaystyle \frac{{\rho }_N}{X}}\int d\mathbf{r}\left\{\exp \left(i{\mathbf{k}}_1·\mathbf{r}\right)+\exp \left(i{\mathbf{k}}_2·\mathbf{r}\right)+\exp \left(i{(\mathbf{k}}_1+{\mathbf{k}}_2)·\mathbf{r}\right)\right\}\left(G_{TLR2}\left(r\right)-1\right)$$

$$+{\rho }_N^2\left(s_{SC1}^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)+g_{LR}^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)+{\mathsf{\Sigma }}_{LR2}^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)+{\mathsf{\Sigma }}_{LR2}^{(3)}\left({\mathbf{k}}_1+{\mathbf{k}}_2,{\mathbf{k}}_1\right)+{\mathsf{\Sigma }}_{LR2}^{(3)}\left({\mathbf{k}}_2,{\mathbf{k}}_1+{\mathbf{k}}_2\right)\right)$$

$$-{\rho }_N^2\left({2\pi )}^3\delta ({\mathbf{k}}_2\right)\int d\mathbf{r}\left\{\exp \left(i{\mathbf{k}}_1·\mathbf{r}\right)\right\}\left(G_{TLR2}\left(r\right)-1\right)$$

$$-{\rho }_N^2\left({2\pi )}^3\delta ({\mathbf{k}}_1\right)\int d\mathbf{r}\left\{\exp \left(i{\mathbf{k}}_2·\mathbf{r}\right)\right\}\left(G_{TLR2}\left(r\right)-1\right)$$

$$-{\rho }_N^2\left({2\pi )}^3\delta ({\mathbf{k}}_1+{\mathbf{k}}_2\right)\int d\mathbf{r}\left\{\exp \left(i{\mathbf{k}}_1·\mathbf{r}\right)\right\}\left(G_{TLR2}\left(r\right)-1\right)$$

$$-{\rho }_N^2\left({2\pi )}^3\delta ({\mathbf{k}}_1\right)\left({2\pi )}^3\delta ({\mathbf{k}}_2\right).$$

(87c)

More information about these structural functions in the r- and the k- spaces can be found in works by this author [11,46,54,55,56,57,62].

The fixing of the $G_{TLR2}$ and $G_{TLR3}$ spatial structures can be accomplished with PIMC or PIMD simulations. Note that their computations scale as $X^2$ and $X^3,$respectively, which makes $G_{TLR3}$a demanding task. To date these calculations have focused on full evaluations of $G_{TLR2}$ [11,46,62], but only on restricted evaluations of $G_{TLR3}$ covering equilateral and isosceles features of the three-particle part $g_{LR3}$ of quantum hard spheres [54]. As for the structure factors $S_{TLR}^{\left(2\right)}\left(\mathbf{k}\right)$ and $S_{TLR}^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right),$ the PIMC/PIMD situation regarding the X-scaling is similar but further aggravated by the necessity to scan sets of commensurate wavevectors and, also, by the drawbacks associated with the low-k regions. Again, a significant part of $S_{TLR}^{\left(2\right)}$, not very close to $k=0,$ can be determined with reasonable PI simulation effort [46]. Nonetheless, the simulations of $S_{TLR}^{(3)}$ over a significant range of wavevectors, excluding obviously the unattainable ${\mathbf{k}}_1=0$ and/or${\mathbf{k}}_2=0$, and ${\mathbf{k}}_1=-{\mathbf{k}}_2$, remain today daunting (at least for most researchers).

OZ2 and OZ3 Frameworks

In view of the present-day PI-simulation difficulties in obtaining even restricted descriptions of the structure factors TLR2 and (above all) TLR3, one might think of resorting to direct correlation functions. However, the delocalization of quantum particles brings about several complications if one tries to construct OZn frameworks by following classical-like ways. A couple of examples will suffice to illustrate this knotty issue. Firstly, just at the pair level, one can consider Haymet et al.’s approach [173,174,175]. Here, the pair part of the TLR2 case Equation (87a), which depends only on $g_{LR2}$ and is normalized to unity for large wavevectors (see below), was the object of interest as applied to freezing in fluid helium-4 [175]. Broadly speaking, some of the main features are: (a) the ideal system was defined as the Feynman ideal system; (b) the direct use of inverse derivatives ${\displaystyle \frac{\delta \mathsf{\Psi }({\mathbf{r}}_1)}{\delta {\rho }_N^{\left(1\right)}({\mathbf{r}}_2)}}$ led to a need for effective particle masses in the ideal system to deal with the self correlations, though this measure could not be implemented under diminishing temperatures; and (c) in the end, a generalized OZ2 framework was obtained, but “effective” structure factors had to be defined, and unnatural amplifications in that pair part of TLR2 were detected for large wavevectors k. Secondly, another generalized OZ2 framework is that put forward by Shinoda et al. [178,179], in which the time-honored molecular approach RISM [7,71] was adapted to tackling the PI-model of a fluid with quantum behavior, as if it were a classical molecular fluid (RISM = reference interaction site model). This approach also focused on the same pair part of TLR2 mentioned above, and self interactions did appear when working the related OZ2 scheme arising from RISM (i.e., a given particle “interacts” with itself at different imaginary times). Accordingly, further approaches to cope with this fact of the model were to be taken.

At this juncture, the following two observations may be useful to help grasp the difficulties in this context. (i) First, a simple exercise. By assuming that for the TLR2 unrestricted bead-bead correlations one might define a direct correlation function in the conventional manner [5]:

$$c_{TLR2}\left({\mathbf{q}}_1{,\mathbf{q}}_2;\mathsf{\Psi }\right)={\displaystyle \frac{\delta (\beta X^{-1}\mathsf{\Psi }\left({\mathbf{q}}_1\right))}{\delta {\rho }_{NX}^{\left(1\right)}({\mathbf{q}}_2;\mathsf{\Psi })}}+{\displaystyle \frac{\delta ({\mathbf{q}}_1-{\mathbf{q}}_2)}{{\rho }_{NX}^{\left(1\right)}({\mathbf{q}}_1;\mathsf{\Psi })}},$$

(88)

it is straightforward to check the strange $S_{TLR}^{(2)}$ behavior that application of an identity analogous to (66) would yield with increasing X$.$ As in the two foregoing approaches, one can rederive the classical OZ2 equation for $X=1,$ but the expression for the Trotter limit $X\to \infty$ turns out to be meaningless. (ii) Second, a mathematical fact. Note that most of the OZ2 frameworks (e.g., Percus-Yevick, hypernetted chain, and their variants [69,71]) are formulated with reference to interatomic potentials $u\left(r_{jl}\right).$ The point is that the potential energy operators $V^{\left(N\right)}({\mathbf{r}}_1,{\mathbf{r}}_2,\dots ,{\mathbf{r}}_N)$ are diagonal in the coordinate representation, which means that beads do interact via $u(r_{jl})$ if and only if they belong to different necklaces and are at the same imaginary time step (e.g., as shown in Equations (36a), (37a), and (38)). Hence, there cannot be any particle self interactions in a PI context.

The approximation proposed by the present author to deal with $S_{TLR}^{(2)}$ consists in separating the intra-necklace bead-bead correlations from the inter-necklace bead-bead correlations and applying directly the classical OZ2 framework to the latter (the pair part mentioned above) [180]. This (lucky) choice found some support by the applicability of this procedure to the GFH picture of fluids with (weak) quantum behavior [11,143,144,180]. This approach can be summarized as follows:

$$S_{TLR}^{\left(2\right)}\left(k\right)={\displaystyle \frac{1}{X}}+{\rho }_N\int d\mathbf{r}{s}_{SC\mathrm{1,2}}\left(r\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(i\mathbf{k}·\mathbf{r}\right)+{\rho }_N\int d\mathbf{r}\left(g_{LR2}\left(r\right)-1\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(i\mathbf{k}·\mathbf{r}\right)=$$

(89a)

$$F_{LR}^{\left(1\right)}\left(k\right)+\left(S_{LR}^{\left(2\right)}\left(k\right)-1\right),$$

(89b)

where $F_{LR}^{\left(1\right)}(k)$ is defined by the first two terms in Equation (89a), being a sort of “form factor” for the one-particle thermal packet, and $S_{LR}^{\left(2\right)}\left(k\right)$ adopts the well-known classical form associated with the overall inter-necklace correlations $g_{LR2}\left(r\right).$ The self correlations do not contribute to the isothermal compressibility of the fluid (i.e., the $\mathbf{k}=0$ component), which can also be formulated in the usual way only with $S_{LR}^{\left(2\right)}\left(k\right)$ [130]. Application to $g_{LR2}\left(r\right)$ of the standard OZ2 Equation (12), which defines its associated $c_{LR2}\left(r\right),$ yields after Fourier transforming:

$$S_{TLR}^{\left(2\right)}\left(k\right)\approx F_{LR}^{\left(1\right)}\left(k\right)+{\displaystyle \frac{{\rho }_N{c}_{LR}^{\left(2\right)}\left(k\right)}{1-{\rho }_N{c}_{LR}^{\left(2\right)}\left(k\right)}}.$$

(90a)

where one notes that $F_{LR}^{\left(1\right)}\left(k\right)$ admits approximate analytical representations [131,180], a very useful one being that derived within the GFH picture [131]

$$F_{LR}^{\left(1\right)}\left(k\right)\approx F_{GFH}^{\left(1\right)}\left(k\right)={\displaystyle \frac{6m}{\beta {\mathrm{\hslash }}^2{k}^2}}\left[1-\mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \frac{\beta {\mathrm{\hslash }}^2{k}^2}{6m}}\right)\right].$$

(90b)

One way or another, quantum approximations that involve classical $c_2(r)$ schemes provide this context with reasonable alternatives to handle k-space difficulties [46,173,178,180], although their corresponding ranges of validity are to be analyzed a posteriori. Surprisingly, for the approximation given by Equation (90a), such range turns out to exceed expectations, since even under very strong diffraction effects the comparisons with experiment and PI-simulations are excellent. In this respect, at $T=4.2\mathrm{K},$ for supercritical helium-3 and for normal liquid helium-4 at and near SVP (saturated vapor pressure) conditions, the TLR2 results based on OZ2(BDH+BHw) were found to be almost indistinguishable from those obtained via PIMC, and/or very close to experiment [30,46,170,181]. Under diffraction effects and at the pair level, the conclusion when conducting OZ2 treatments is clear: the intra- and inter-necklace bead-bead correlations can be separated safely for most practical purposes. Two more observations: (a) with increasing quantum effects, the very low-k region may be affected by the OZ2 approximations used (see subsection 5.4); and (b) when exchange interactions are not negligible, this sort of bead-bead separation cannot be applied.

So far, to this author’s knowledge, no related OZ3 schemes have been reported for TLR3. The intricacy of these triplet features is certainly greater than at the pair level. In relation to this, one notes that the straightforward procedure followed in the pair linear response structures, Equations (86a)-(89)-(90), cannot be used in the functions leading to $S_{TLR}^{\left(3\right)}$, as shown in Equations (86b) and (87c). Accordingly, whether this type of OZ ideas might be fruitful when studying the triplet $S_{TLR}^{\left(3\right)}$ remains to be investigated. The use of approximations, extending those in Equations (90), for treating the terms displayed in Equation (86b) will be necessary.

5.3. The PI Instantaneous ETn Class

This last class is not amenable to being treated with functional derivatives involving an external weak field, because of the collapse of the thermal packet under an instantaneous localization process. To obtain the definitions of the related pair correlation function $g_{Q,ET2}\left(r\right)$ and its associated structure factor $S_{Q,ET}^{\left(2\right)}\left(k\right),$ which are in many ways analogues to their classical counterparts$,$ the standard quantum reasoning is based on the study of: (a) neutron scattering, with the time-dependent treatment using Born’s approximation and Fermi’s potential, plus the elastic sum rule applied to the dynamic structure factor $S_Q^{dyn}(\mathbf{k},\omega )$[6,23,29]; or (b) the elastic scattering of X-rays, plus the consideration only of the pair coherent part associated with the nuclei (i.e., the single-atom actual form factor due to the electrons is dealt with separately) [6,28]. For the diffraction quantum effects the final expressions are [3,6,9,11]:

$$S_{Q,ET}^{\left(2\right)}\left(k\right)={\displaystyle \frac{1}{〈N〉}}〈{\sum }_{j=1}^N{\sum }_{l=1}^N\mathrm{e}\mathrm{x}\mathrm{p}\left(i\mathbf{k}·({\mathbf{r}}_j-{\mathbf{r}}_l)\right)〉-{\left(2\pi \right)}^3{\rho }_N\delta \left(\mathbf{k}\right)=$$

$$1+{\rho }_N\int d\mathbf{r}\mathrm{e}\mathrm{x}\mathrm{p}\left(i\mathbf{k}·\mathbf{r}\right)\left(g_{Q,ET2}\left(r\right)-1\right),$$

(91)

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

(92)

where the ensemble averages refer to ${\mathsf{\Xi }}_Q$ Equation (29), both being formally equivalent to the classical definitions $\left(r=\left|{\mathbf{q}}_1-{\mathbf{q}}_2\right|\right)$. The latter remark implies that both quantum averages reduce to the usual formulas in the classical limit. Explicitly, the pair correlations are defined via the ensemble average:

$${{\rho }_N^{2}g}_{Q,ET2}\left(r\right)={\displaystyle \frac{1}{{\mathsf{\Xi }}_Q}}{\sum }_{N\ge 0}{e}^{\beta \mu N}Tr\left\{\mathrm{e}\mathrm{x}\mathrm{p}(-\beta H_0^{(N)})\times {\sum }_{j\ne l}\delta ({\mathbf{r}}_j-{\mathbf{q}}_1)\delta ({\mathbf{r}}_l-{\mathbf{q}}_2)\right\},$$

(93)

and the PI approach gives:

$${{\rho }_N^{2}g}_{Q,ET2}\left(r\right)\approx {{\rho }_N^{2}g}_{ET2}\left(r\right)={\displaystyle \frac{1}{{\mathsf{\Xi }}_{PI,D}}}{\sum }_{N\ge 0}{\displaystyle \frac{{\mathfrak{z}}^N}{N!}}\int {\prod }_{j=1}^N{\prod }_{t=1}^{P}d{\mathbf{r}}_j^t\times \mathrm{e}\mathrm{x}\mathrm{p}(-\beta W_{NP})\times$$

$$\left({\displaystyle \frac{1}{X}}{\sum }_{j\ne l}{\sum }_{\tau =1}^X\delta ({\mathbf{r}}_j^{\tau }-{\mathbf{q}}_1)\left.\delta ({\mathbf{r}}_l^{\tau }-{\mathbf{q}}_2)\right)\right.,$$

(94)

or in compact form

$${{\rho }_N^{2}g}_{ET2}\left(r\right)=〈{\displaystyle \frac{1}{X}}{\sum }_{j\ne l}{\sum }_{\tau =1}^X\delta ({\mathbf{r}}_j^{\tau }-{\mathbf{q}}_1)\delta ({\mathbf{r}}_l^{\tau }-{\mathbf{q}}_2)〉,$$

(95a)

where the selection of the structurally significant (and equivalent) equal $\tau$-time beads is to be highlighted; this is the distinctive trait of the ETn class. Likewise, the instantaneous structure factor can be calculated within PI as [9,11]:

$$S_{Q,ET}^{\left(2\right)}\left(k\right)\approx S_{ET}^{\left(2\right)}\left(k\right)={\displaystyle \frac{1}{〈N〉X}}〈\sum _{j=1}^N\sum _{l=1}^N\sum _{\tau =1}^X\mathrm{e}\mathrm{x}\mathrm{p}\left(i\mathbf{k}·\left({\mathbf{r}}_j^{\tau }-{\mathbf{r}}_l^{\tau }\right)\right)〉-{\left(2\pi \right)}^3{\rho }_N\delta \left(\mathbf{k}\right)=$$

$$1+{\rho }_N\int d\mathbf{r}\mathrm{e}\mathrm{x}\mathrm{p}\left(i\mathbf{k}·\mathbf{r}\right)\left(g_{ET2}\left(r\right)-1\right).$$

(95b)

It is a simple matter to extend the formulation for the instantaneous triplets in r-space, the PI average being (Figure 3):

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

(96)

The spontaneous fluctuations in the density at the pair and the triplet levels can be expressed by formulas, involving ${\mathsf{\Xi }}_{QD}$ Equation (30), analogous to the classical Equations (9b)-(9c):

$${\mathsf{\Gamma }}_{Q,ET2}\left({\mathbf{q}}_1,{\mathbf{q}}_2;\mathsf{\Psi }=0\right)=〈\left({\sum }_{j=1}^N\delta \left({\mathbf{r}}_j-{\mathbf{q}}_1\right)-{\rho }_N\right)\left({\sum }_{l=1}^N\delta \left({\mathbf{r}}_l-{\mathbf{q}}_2\right)-{\rho }_N\right)〉,$$

(97a)

$${\mathsf{\Gamma }}_{Q,ET3}\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3;\mathsf{\Psi }=0\right)=$$

$$〈\left({\sum }_{j=1}^N\delta \left({\mathbf{r}}_j-{\mathbf{q}}_1\right)-{\rho }_N\right)\left({\sum }_{l=1}^N\delta \left({\mathbf{r}}_l-{\mathbf{q}}_2\right)-{\rho }_N\right)\left({\sum }_{m=1}^N\delta \left({\mathbf{r}}_m-{\mathbf{q}}_3\right)-{\rho }_N\right)〉.$$

(97b)

and, therefore, the actual and the PI definitions of the pair and triplet instantaneous structure factors can be written as:

$$S_{Q,ET}^{\left(2\right)}\left(k\right)={\displaystyle \frac{1}{{\rho }_N}}{\mathsf{\Gamma }}_{Q,ET}^{(2)}\left(\mathbf{k};\mathsf{\Psi }=0\right){\to }_{PI}{\approx S}_{ET}^{\left(2\right)}\left(k\right)={\displaystyle \frac{1}{{\rho }_N}}{\mathsf{\Gamma }}_{ET}^{(2)}\left(\mathbf{k};\mathsf{\Psi }=0\right),$$

(98a)

$$S_{Q,ET}^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)={\displaystyle \frac{1}{{\rho }_N}}{\mathsf{\Gamma }}_{Q,ET}^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2;\mathsf{\Psi }=0\right){\to }_{PI}\approx S_{ET}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)={\displaystyle \frac{1}{{\rho }_N}}{\mathsf{\Gamma }}_{ET}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2;\mathsf{\Psi }=0\right).$$

(98b)

Once again, this definition of the triplet structure factor is very convenient, since it reduces to the classical formulation when neglecting quantum effects [3]. For completeness, note that within PI one finds the explicit triplet formulas [55]:

$${\mathsf{\Gamma }}_{ET3}\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3;\mathsf{\Psi }=0\right)=$$

$$〈{\displaystyle \frac{1}{X}}{\sum }_{\tau =1}^X\left\{\left({\sum }_{j=1}^N\delta \left({\mathbf{r}}_j^{\tau }-{\mathbf{q}}_1\right)-{\rho }_N\right)\left({\sum }_{l=1}^N\delta \left({\mathbf{r}}_l^{\tau }-{\mathbf{q}}_2\right)-{\rho }_N\right)\left({\sum }_{m=1}^N\delta \left({\mathbf{r}}_m^{\tau }-{\mathbf{q}}_3\right)-{\rho }_N\right)\right\}〉,$$

(99)

$$S_{ET}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)={\displaystyle \frac{1}{〈N〉X}}〈{\sum }_{j=1}^N{\sum }_{l=1}^N{\sum }_{m=1}^N{\sum }_{\tau =1}^X\mathrm{e}\mathrm{x}\mathrm{p}\left(i\left({\mathbf{k}}_1·{\mathbf{r}}_j^{\tau }+{\mathbf{k}}_2·{\mathbf{r}}_l^{\tau }-\left({{\mathbf{k}}_1+\mathbf{k}}_2\right)·{\mathbf{r}}_m^{\tau }\right)\right)〉-\left\{\delta \left(\mathbf{k}\right)-\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{s}\right\}$$

(100)

Some related remarks are the following. (a) From the simulation standpoint: note that the ETn structural computations scale as X, which make them more expensive than those of the classical case ($X=1),$but far less expensive than those of TLRn which scale as $X^3,$ as shown in Equation (87b). (b) The development of Equation (100) is formally the same as that of the classical case or, equivalently, that of the PI centroids given in Equation (54), where one only has to substitute the CMn quantities for the ETn ones. (c) As stressed earlier, there is a full formal equivalence among the Cn, CMn, and ETn basic formulations of the isolated fluid quantities in r-space and k-space $\left\{g_n,S^n\right\}.$

5.3.1. OZ2 and OZ3 Frameworks

At this stage of the presentation, it is clear that no exact OZn frameworks are available for analyzing the structural functions of the ETn class. However, as occurred in the TLRn class, the adoption of classical schemes/ideas, either as such [3] or modified [69], seems unavoidable if one tries to use the device of direct correlation functions. In this regard, and again based on the applicability of the classical-OZn to the PI-based variational pictures [11,143,144,145,180], there was nothing to lose by trying the direct use of the classical OZ2 and OZ3 schemes. It goes without saying that the OZ2(BDH+BHw) results were beyond expectations. Thus, under very strong quantum diffraction effects (e.g., helium-4 and helium-3 fluids at $T\approx 4.2\mathrm{K}$), the overall agreement with experiment and PI simulation was excellent [28,29,46,55,170,181], the very low-k region being the most sensitive to this approximation (see subsection 5.4). As regards OZ3, applications of the classical scheme have been made to helium-3 [55,60], the quantum hard sphere fluid (bare of with Yukawa attractions) [56], and liquid parahydrogen [57,58], but no full tests of validity have been conducted.

The two basic equations here are the adaptations of the classical Equations (13) and (15), i.e., formally the same as in CMn, although now they are approximations:

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

(101a)

$$S_{ET}^{(3)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)\approx S_{ET}^{\left(2\right)}\left({\mathbf{k}}_1\right)S_{ET}^{\left(2\right)}({\mathbf{k}}_2)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\},$$

(101b)

where the direct correlation functions $c_{ET}^{\left(2\right)}$ and $c_{ET}^{\left(3\right)}$ are introduced as arising from the standard OZ2 and OZ3 classical procedures (Figure 4 and Figure 5).

5.3.2. An ETn Functional Digression

Before going any further, it is worthwhile to remark that there is a functional approach to the ET2 correlation function [6,135]. The point is that, by focusing on the variations of the partition function $\delta {\mathsf{\Xi }}_{PI,D}$ with respect to the variations in the interparticle potential $\delta u\left(r_{jl}^t\right),$ one obtains an expression in which the correlation function $g_{ET2}$ shows up. The general usefulness of this development seems scarce because there is no response from the fluid to an external field. However, for completeness, it is noticeable that one can obtain useful relations at the pair-level involving both the ET2 and the CM2 structures [135]. To simplify the following discussion, only the PP propagator is used in what follows.

The PP grand canonical partition function for the homogeneous and isotropic fluid can be cast as [7,90]:

$${\mathsf{\Xi }}_{PI,D}\left(\mathrm{P}\mathrm{P}\right)={\sum }_{N\ge 0}{\displaystyle \frac{{\mathfrak{z}}^N}{N!}}\int {\prod }_{j=1}^N{\prod }_{t=1}^{P}d{\mathbf{r}}_j^t\times \mathrm{e}\mathrm{x}\mathrm{p}\left[-\beta W_{1F}-{\displaystyle \frac{\beta }{P}}{\sum }_{j<l}{\sum }_{t=1}^{P}u(r_{jl}^t)\right],$$

(102)

A straightforward calculation yields the related first functional derivative:

$$-k_{B}T{\displaystyle \frac{\delta \ln {\mathsf{\Xi }}_{PI,D}(PP)}{\delta u\left(r_{12}\right)}}={\displaystyle \frac{\mathcal{V}}{2}}{\rho }_{N,ET}^{\left(2\right)}\left({\mathbf{q}}_1,{\mathbf{q}}_2\right)={\displaystyle \frac{\mathcal{V}}{2}}{{\rho }_N^{2}g}_{ET2}\left(r_{12}\right),(r_{12}\equiv r_{12}^t\text{ equal-time distance),}$$

(103)

which is formally identical to the classical or to the PI-ET2 counterpart in the canonical ensemble that can be found in [6] and [135], respectively.

As regards the extensions to higher orders, the situation is rather inconclusive. To get a feeling of the drawbacks posed by this type of formulation, it seems instructive to consider the development of ${\delta }^2{\mathsf{\Xi }}_{PI,D}\left(\mathrm{P}\mathrm{P}\right)$:

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

$${\left(-{\displaystyle \frac{\beta }{P}}\right)}^2\left\{{\sum }_{j<l}{\sum }_{m<n}{\sum }_{t=1}^P{\sum }_{t^{\mathrm{\text{'}}}=1}^P\delta {u(r}_{jl}^t)\delta u(r_{mn}^{t^{\mathrm{\text{'}}}})\right\},$$

(104)

the main problem being that different imaginary times are mixed. One might even think of extracting artificially equal-imaginary time quantities by “tinkering” with the t-sums, but these manipulations are of no interest to ET correlations. Furthermore, beyond the pair level expressed in Equation (103), these derivatives become increasingly entangled. Although the latter might be related to spontaneous effects within the quantum fluid, they seem to be of no use for the current study of static quantum fluid structures.

5.4. A Joint Consideration of CMn, TLRn, and ETn

(i) Apart from the general consideration that the three quantum classes of linear response functions are directly related to different forms of density-correlation functions in the isolated fluid (as in the classical case), there are some properties exactly shared by the three classes. These properties obviously belong to the isolated quantum fluid and can be identified at zero-field, $\mathsf{\Psi }={\mathsf{\Psi }}_F=0$. A most interesting case arises for the structure factors at zero momentum transfer(s) $(\mathbf{k}=0)$, because the distinct response functions of the (isolated) fluid must take the same value regardless of the pair, or the triplet, structure considered.

At the pair level there is the extended compressibility theorem [130] that states:

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

(105)

It is worth remarking that the foregoing exact result can be distorted when classical-like OZ2 calculations with ET2 and LR2 are carried out. On the one hand, the discrepancies from the OZ2-CM2 exact ${\chi }_T$ estimate remain in general controlled and do not alter the overall representation of the k-space structures determined. On the other hand, the importance of the disagreement may be a matter to analyze in certain cases (e.g., near changes of phase), as both the density and the temperature play roles in this issue [46,62,130,136,137,180].

The generalization of Equation (105) to the triplet level, involving the double-zero momentum transfer, reads as [57]:

$$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)=S_{TLR}^{\left(3\right)}\left({\mathbf{k}}_1=0,{\mathbf{k}}_2=0\right)=S^{\left(3\right)}\left(\mathrm{0,0}\right),$$

(106)

By taking advantage of the OZn formulation for centroids, Equation (106) can be written in terms of $S^{\left(2\right)}\left(0\right).$Using Equations (69), (77c) and (78b), one finds

$$c_{CM}^{\left(3\right)}\left({\mathbf{k}}_1=0,{\mathbf{k}}_2=0\right)={\displaystyle \frac{\partial c_{CM}^{\left(2\right)}(k=0)}{\partial {\rho }_N}}=-{\displaystyle \frac{1}{{\rho }_N^2}}+{\displaystyle \frac{1}{{\rho }_N^2{S}_{CM}^{\left(2\right)}\left(k=0\right)}}+{\displaystyle \frac{1}{{\rho }_N{\left[S_{CM}^{\left(2\right)}\left(k=0\right)\right]}^2}}{\displaystyle \frac{\partial S_{CM}^{\left(2\right)}\left(k=0\right)}{\partial {\rho }_N}};\text{T=constant},$$

(107)

which in view of Equations (105) and (106) leads to the following result valid irrespective of the response class [57]:

$$S^{\left(3\right)}\left(\mathrm{0,0}\right)=S^{\left(2\right)}\left(0\right)\left\{S^{\left(2\right)}\left(0\right)+{\rho }_N{\displaystyle \frac{\partial S^{\left(2\right)}\left(0\right)}{\partial {\rho }_N}}\right\};\text{T=constant; CM, ET, TLR.}$$

(108a)

The quantity $S^{\left(3\right)}\left(\mathrm{0,0}\right)$ is a quantum fluid property related to higher-order fluctuations in the number of the atoms (particles), and it is a simple matter to prove that:

$$S^{\left(3\right)}\left(\mathrm{0,0}\right)={\displaystyle \frac{〈N^3〉}{〈N〉}}-3〈N^2〉+2{〈N〉}^2;\text{T=constant; CM, ET, TLR,}$$

(108b)

The trivial connections of the zero-field Equations (105) and (108b) with the functional Equations (80d)-(80e) and (83f)-(83g), at zero field, are to be noticed (ET and CM are completely parallel in this regard).

Also, for the single-zero momentum transfer one finds the exact CM relationship [57]

$$S_{CM}^{(3)}\left(k,0\right)=S_{CM}^{(2)}\left(0\right)\left\{S_{CM}^{(2)}\left(k\right)+{\rho }_N{\displaystyle \frac{\partial S_{CM}^{(2)}\left(k\right)}{\partial {\rho }_N}}\right\};k\ge 0,T=\mathrm{constant},$$

(109)

whilst for ET an analogous equation can also be written, but owing to the OZn properties utilized it is an approximation:

$$S_{ET}^{(3)}\left(k,0\right)\approx S_{ET}^{(2)}\left(0\right)\left\{S_{ET}^{(2)}\left(k\right)+{\rho }_N{\displaystyle \frac{\partial S_{ET}^{(2)}\left(k\right)}{\partial {\rho }_N}}\right\};k>0,T=\mathrm{constant}.$$

(110)

It seems worthwhile to explore in future work the possibility of obtaining an equation for $S_{ET}^{(3)}\left(k,0\right)$ based on the connection between the instantaneous $g_{ET}^{(2)}$ and centroid $g_{CM}^{(2)}$ pair correlation functions [135]. If successful, this could help to solve the ET problems for low-k wavenumbers, generalizing in this way the exact result connecting the properties of the centroid and the delocalized particle within the GFH picture [131,180]. So far, no further formulas based on direct correlation functions for TLR3 have been derived.

(ii) The symmetry properties in r-space and in k-space of the structures CMn, TLRn, and ETn, are translations of the classical ones seen in subsection 3.3. CMn and ETn behave in the classical way regarding the limits of increasing/decreasing distances between atoms and, also, in the limit of increasing wavenumbers (see also Figures 1, 2b, and 3). However, the class TLRn (unrestricted bead-bead-… correlations) presents some striking differences. In relation to this, suffice it to mention that for the TLR2 functions one finds in r-space: (a) $G_{TLR2}(r)$ at $r=0$ increases with the P discretization, because of its one-atom component $s_{SC\mathrm{1,2}}$, and tends to unity for long distances due to the two-atom component $g_{LR2}$; (b) $s_{SC\mathrm{1,2}}\to 0$ for long distances, and $g_{LR2}$ does not necessarily vanish when two atoms come close together (it may even be nonzero at $r=0!$); (c) $S_{TLR}^{\left(2\right)}(k)\to 0$ for large wavenumbers. The magnitudes of these discrepancies from the standard classical behaviors depend on the intensity of the quantum effects, and the reader is referred to [46,54,55,62,170,180] for related results. Consequently, one may expect an interesting variety of patterns when studying TLR3 functions.

(iii) Another set of exact PI relationships can be obtained [54] by following the classical developments [26,150,151,157] given in Equations (27) and (28). It is easy to establish that for the two-point densities (referring to actual two-particle correlations) given by:

$${\mathsf{\Phi }}_{CM2}={\sum }_{j\ne l}\delta \left({\mathbf{R}}_{CM,j}-{\mathbf{R}}_1\right)\delta \left({\mathbf{R}}_{CM,l}-{\mathbf{R}}_2\right),$$

(111a)

$${\mathsf{\Phi }}_{ET2}=X^{-1}{\sum }_{\tau =1}^X{\sum }_{j\ne l}\delta \left({\mathbf{r}}_j^{\tau }-{\mathbf{q}}_1\right)\delta \left({\mathbf{r}}_l^{\tau }-{\mathbf{q}}_2\right),$$

(111b)

$${\mathsf{\Phi }}_{LR2}=X^{-2}{\sum }_{\tau =1}^X{\sum }_{\tau \mathrm{\text{'}}=1}^X{\sum }_{j\ne l}\delta \left({\mathbf{r}}_j^{\tau }-{\mathbf{q}}_1\right)\delta \left({\mathbf{r}}_l^{\tau \mathrm{\text{'}}}-{\mathbf{q}}_2\right),$$

(111c)

one obtains the classical-like relationships:

$${\mathfrak{z}\left({\displaystyle \frac{\partial 〈\mathsf{\Phi }〉}{\partial \mathfrak{z}}}\right)}_T=〈N\mathsf{\Phi }〉-〈N〉〈\mathsf{\Phi }〉;\text{ CM2, ET2, LR2.}$$

(112)

The latter formula leads to the classical-like equations:

$${{\rho }_N{S}^{\left(2\right)}\left(0\right){\left({\displaystyle \frac{\partial 〈\mathsf{\Phi }〉}{\partial {\rho }_N}}\right)}_T=\rho }_N^3\int d{\mathbf{q}}_3\left\{{\phi }_3\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3\right)-{\phi }_2\left({\mathbf{q}}_1,{\mathbf{q}}_2\right)\right\}+2{\rho }_N^2{\phi }_2\left({\mathbf{q}}_1,{\mathbf{q}}_2\right),$$

(113a)

$$S^{\left(2\right)}\left(0\right){\left({\displaystyle \frac{\partial \mathrm{l}\mathrm{n}{\phi }_2(r_{12})}{\partial {\rho }_N}}\right)}_T=\int d{\mathbf{q}}_3\left\{{\displaystyle \frac{{\phi }_3\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3\right)}{{\phi }_2\left({\mathbf{q}}_1,{\mathbf{q}}_2\right)}}-1\right\}+{\displaystyle \frac{2}{{\rho }_N}}\left(1-S^{\left(2\right)}\left(0\right)\right);{\phi }_2\left(r_{12}\right)>0,$$

(113b)

where the functions $\left\{{\phi }_2,{\phi }_3\right\}$ stand for:

$$\left\{{\phi }_2,{\phi }_3\right\}=\left\{g_{CM2},g_{CM3}\right\},\left\{g_{ET2},g_{ET3}\right\},\left\{g_{LR2},g_{LR3}\right\},$$

(113c)

and, for consistency, $S^{\left(2\right)}\left(0\right)$ corresponds to its associated pair structure ${\phi }_2$ (recall that in practice the $S^{\left(2\right)}\left(0\right)$ values may show deviations from the exact Equation (105)).

Note that the total continuous linear response TLR is not contained in the LR relationships, since the self-correlations do not fit into this scheme. Equations (113) arise only from the pure actual two- and three-atom correlations.

For the purposes of this work, it is also worthwhile to remark that Equation (113b) can be transformed and split into two parts$\left[\mathrm{150,151}\right]$ as $\left({\phi }_2\left(r\right)>0\right)$:

$${\tilde{\Lambda }}_{exact}\left(r_{12}\right)\equiv \mathsf{\Lambda }\left(r_{12}|{\phi }_3\right);\mathrm{if}{\phi }_3=exacttripletfunction,$$

(113d)

where

$${\tilde{\Lambda }}_{exact}\left(r_{12}\right)=S^{\left(2\right)}\left(0\right){\left({\displaystyle \frac{\partial \mathrm{l}\mathrm{n}{\phi }_2(r_{12})}{\partial {\rho }_N}}\right)}_T-\int d{\mathbf{q}}_3\left({\phi }_2\left(r_{13}\right)-1\right)\left({\phi }_2\left(r_{23}\right)-1\right),$$

(113e)

$$\mathsf{\Lambda }\left(r_{12}|{\phi }_3\right)=\int d{\mathbf{q}}_3\left\{{\displaystyle \frac{{\phi }_3\left(r_{12},r_{13},r_{23}\right)}{{\phi }_2\left(r_{12}\right)}}-{\phi }_2\left(r_{13}\right){\phi }_2\left(r_{23}\right)\right\}.$$

(113f)

Given that ${\tilde{\Lambda }}_{exact}\left(r_{12}\right)$ is just a pair quantity that can be fixed with high accuracy, whilst $\mathsf{\Lambda }\left(r_{12}|{\phi }_3\right)$ is a quantity that depends explicitly on ${\phi }_3,$ the use of these two quantities yields consistency tests of possible closures for ${\phi }_3$. In this connection, one uses KS3 as a reference by expressing the triplet closure in the form:

$${\left[{\phi }_3\left(r_{12},r_{13},r_{23}\right)\right]}_{clos.}={\phi }_{KS3}\times {\stackrel{̿}{\phi }}_3\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3\right)={\phi }_2\left(r_{12}\right){\phi }_2\left(r_{13}\right){\phi }_2\left(r_{23}\right)\times {\stackrel{̿}{\phi }}_3\left({\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3\right).$$

(113g)

By substituting a closure to be tested in Equation (113f), two comparison tests are between: (a) ${\tilde{\Lambda }}_{exact}\left(r_{12}\right)$ and $\mathsf{\Lambda }\left(r_{12}|{\phi }_3\right)$, and (b) the values of the density derivative ${\left(\partial {\phi }_2/\partial {\rho }_N\right)}_T$ determined with the computed pair results and the corresponding estimates arising from the closure.

In relation to this, for completeness, the three basic closures KS3 Equation (22), JF3 Equation (23), and AV3 Equation (24), lead to the following expressions for their quantities in Equation (113d) (Figure 6):

-

KS3:

$$\mathsf{\Lambda }\left(r_{12}|{\phi }_{KS3}\right)\equiv 0,$$

(113h)

$${\left(\partial {\phi }_2(r_{12})/\partial {\rho }_N\right)}_{T,KS3}=\left({\phi }_2/S^{\left(2\right)}(0)\right)\int d{\mathbf{q}}_3\left({\phi }_2\left(r_{13}\right)-1\right)\left({\phi }_2\left(r_{23}\right)-1\right).$$

(113i)

-

JF3:

$$\mathsf{\Lambda }\left(r_{12}|{\phi }_{JF3}\right)=\left({\displaystyle \frac{S^{\left(2\right)}\left(0\right)}{{\phi }_2}}-1\right)\int d{\mathbf{q}}_3\left({\phi }_2\left(r_{13}\right)-1\right)\left({\phi }_2\left(r_{23}\right)-1\right),$$

(113j)

$${\left(\partial {\phi }_2(r_{12})/\partial {\rho }_N\right)}_{T,JF3}=\int d{\mathbf{q}}_3\left({\phi }_2\left(r_{13}\right)-1\right)\left({\phi }_2\left(r_{23}\right)-1\right).$$

(113k)

-

AV3:

$$\mathsf{\Lambda }\left(r_{12}|{\phi }_{AV3}\right)={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{S^{\left(2\right)}\left(0\right)}{{\phi }_2}}-1\right)\int d{\mathbf{q}}_3\left({\phi }_2\left(r_{13}\right)-1\right)\left({\phi }_2\left(r_{23}\right)-1\right),$$

(113l)

$${\left(\partial {\phi }_2(r_{12})/\partial {\rho }_N\right)}_{T,AV3}={\displaystyle \frac{1}{2}}\left(1+{\phi }_2/S^{\left(2\right)}(0)\right)\int d{\mathbf{q}}_3\left({\phi }_2\left(r_{13}\right)-1\right)\left({\phi }_2\left(r_{23}\right)-1\right).$$

(113m)

(iv) Furthermore, as regards OZn schemes, application of the classical and quantum reasonings [38,39,40,41,59,60,72,125,126,127] is a possible and rewarding way to obtain information on the quantum triplet CM3 and ET3 structures. Recall that the CM2 calculations have, in the end, an exact OZ2 framework, but that under any other circumstances approximations are involved. Despite the latter remark, the usual PIMC/PIMD (“raw”) pair radial structures CM2, ET2, LR2, arising from the canonical simulations, can be: (a) utilized to carry out OZ2 studies of the asymptotic decay properties [183,184], and (b) improved using OZ2-based treatments and corrections (e.g., the BDH+BHw analysis considered in Section 3.4). By doing so, grand-canonical approximations for these pair functions can be determined, which may be viewed as better intermediate quantities for an intensive use of direct correlation functions in triplet calculations via closures. In this connection, although for r-space triplets the PI simulations in the canonical ensemble are sufficient to give very useful information, note that for k-space triplets there are no complete OZ3 tests regarding the (density-temperature) ranges of validity. The huge PI computational load involved, covering significant ranges of the two wavevectors and physical conditions, and the variety of closures for triplets are reasons for explaining this lack of knowledge (Figure 4, Figure 5 and Figure 7).

(v) Finally, it is worthwhile to write the connections (“sum rules”) of the dynamic structure factor $S_Q^{dyn}\left(\mathbf{k},\omega \right),$ which can be obtained through inelastic neutron diffraction experiments [23], with the PI quantities $S_{ET}^{\left(2\right)}\left(\mathbf{k}\right)$ and $S_{TLR}^{\left(2\right)}\left(\mathbf{k}\right):$

$$S_{Q,ET}^{\left(2\right)}\left(\mathbf{k}\right)={\int }_{-\mathbf{\infty }}^{\mathbf{\infty }}d\omega S_Q^{dyn}(\mathbf{k},\omega ){\to }_{PI}{S}_{ET}^{\left(2\right)}\left(\mathbf{k}\right)$$

(114a)

$$S_{Q,TLR}^{\left(2\right)}\left(\mathbf{k}\right)=2{\int }_0^{\mathbf{\infty }}d\omega {\displaystyle \frac{1-\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta \mathrm{\hslash }\omega \right)}{\beta \mathrm{\hslash }\omega }}{S}_Q^{dyn}(\mathbf{k},\omega ){\to }_{PI}{S}_{TLR}^{\left(2\right)}\left(\mathbf{k}\right),$$

(114b)

where $\omega$ stands for the angular frequencies of the neutron radiation. $S_{Q,TLR}^{\left(2\right)}\left(\mathbf{k}\right)$ arises as a density-density relaxation function [23] (Equation (87a)). The latter expressions have been utilized in fluid helium-4 studies (e.g., see References [29,30,170,178]).

6. Systems Studied in This Work, Computational Details and Related Observations

The main structural concepts reviewed in the previous Sections are illustrated in this work with further numerical applications to supercritical helium-3, and the QHS fluid on the crystallization line, thus expanding the scope of previous works on these key systems [54,55,59,60]. For helium-3 the focus is on pair and triplet structures in r-space and triplet structures in k-space, whereas for quantum hard spheres is on triplets in k-space. The computational techniques employed are PIMC simulations in the canonical ensemble $\left(N,\mathcal{V},T\right),$ (or better $\left(N\times P,\mathcal{V},T\right)),$ complemented with OZ2 and OZ3 treatments, by using the closures and relationships addressed in this work. The simulations utilize sample sizes $N_S\times P$ inside a central cubic box of side L, with the usual periodic boundary conditions, and the Metropolis sampling scheme involving the normal-mode algorithm (acceptance ratio 50%) [74,145]. One PIMC kpass is defined as $\left({{10}^{3}N}_S\times P\right)$ attempted bead moves, and one PIMC Mpass is ${10}^3$kpasses. The computed structural classes are the centroid (CMn) and the instantaneous (ETn) at the pair and triplet levels, the latter being focused on some equilateral and isosceles features. Representative results are contained in the Supplementary Material.

The state points studied are listed below:

(i) For helium-3: ($T=4.21\mathrm{K}$; ${\rho }_N=0.0025{Å}^{-3})$, ($T=4.2\mathrm{K}$; ${\rho }_N=0.02286713{Å}^{-3}),$ ($T=(8.99\mathrm{K}$; ${\rho }_N=0.0228717687{Å}^{-3}),$ The pair interatomic potential selected is Janzen-Aziz’s SAPT2 [106], and PIMC calculations involve the SCVJ propagator with $\alpha =1/3$ [84]. Information obtained in Reference [46] on auxiliary state points about the state point at $4.2\mathrm{K}$: $\mathsf{\Delta }{\rho }_N=\pm 0.002{Å}^{-3},\pm 0.004{Å}^{-3},$ is also employed.

(ii) For the QHS fluid [59,62,137] only k-space results for one state point are reported, $\left({{\lambda }_B^{*};\rho }_N^{*}\right)=$ $\left(0.6;0.589\right).$ Complementary structures for two state points about the latter, by varying ${\rho }_N^{*}$ at constant ${\lambda }_B^{*}$, are also obtained in the present work; these variations are: ${\mathsf{\Delta }\rho }_N^{*}=\pm 0.01715$, which for a system with ${\sigma }_{HS}=3.5Å$ amount to $\pm 0.0004{Å}^{-3}$. (Other pilot applications along the crystallization line, at ${\lambda }_B^{*}=0.2,0.8,1.2543,1.9832$ [62,79], have also been carried out and will be commented later). PIMC calculations involve the CBHSP propagator [81]. Recall that there is a unique way to characterize a state point in the QHS system by using reduced units [59,62]: i.e., the length unit is the hard-sphere diameter, ${\rho }_N^{*}={\rho }_N{\sigma }_{HS}^3,$ ${\lambda }_B^{*}={\displaystyle \frac{{\lambda }_B}{{\sigma }_{HS}}}$.

PIMC canonical simulations focus on the instantaneous and centroid pair $g_2(r)$ and/or triplet $g_3(r,s,s)$ correlation functions. The sample sizes $\left(N_S\times P\right)$ employed are as follows: (a) for $g_2$ and $g_3$ in helium-3, $\left(1372\times 66\right)$ at $T=4.21\mathrm{K}$, $\left(1372\times 22\right)$ at $T=8.99\mathrm{K}$, and $\left(1024\times 66\right)$ at $T=4.2\mathrm{K}$ [46]; (b) for $g_2$ in the QHS fluid, $\left(864\times 12\right)$ at the auxiliary state points analyzed, and the same at the target state point $\left({{\lambda }_B^{*};\rho }_N^{*}\right)=$ $\left(0.6;0.589\right)$ taking its pair information from calculations in Reference [59].

The pair structures in r-space are computed in the usual way utilizing histograms [13], the width of the bins being $\mathsf{\Delta }=0.1Å.$ Most of the “raw” pair canonical structures are subsequently OZ2-analyzed and improved in the form already described [46,136]: application of the OZ2 Baxter-Dixon-Hutchinson’s (BDH) procedure [72,125] plus Baumketner-Hiwatari (BHw) corrections [127]. By doing so, access to every k-space related property at the pair level is obtained with high accuracy (Figure 2).

As regards the MC/MD calculations of $g_3$ triplets, they present a number of nontrivial subtleties, as discussed by Tanaka and Fukui [42], Baranyai and Evans [43,44], and Barrat et al. [41]. In the current PIMC applications, the method followed in r-space is based on References [42,43,44], which for the study of quantum fluids was firstly used in Reference [54]. The interested reader is referred to the latter reference for specific details. To facilitate the understanding of the reported results, however, suffice to say that triplet correlations are determined via:

$$g_3\left(R,S,U\right)={\displaystyle \frac{\left(\mathsf{\Delta }{n}_T\right)}{N{\rho }_N^2{\left(\mathsf{\Delta }\mathcal{V}\right)}_{RSU}}};\text{ CM3, ET3,}$$

(115)

where $\left(\mathsf{\Delta }{n}_T\right)$ is the number of times that mutual distance triplets lie within the ranges $R-\mathsf{\Delta }\le r\le R,$ $S-\mathsf{\Delta }\le s\le S,$ and $U-\mathsf{\Delta }\le u\le U$ (again: $\mathsf{\Delta }=0.1Å);$ the previous notation should be obvious: $R=r_{12},$ $S=r_{13},$ $U=r_{23}.$ In addition, to avoid redundancies in the counting, the conservative (and less time-consuming) condition $r,s,u<L/4$ is imposed [42]. The run lengths of the PIMC computations for helium-3 in r-space are as follows: in between 1233-1300 kpasses for fixing the pair structures $g_{CM2}(r)$and $g_{ET2}\left(r\right),$ gathering statistics every 8000-12500 accepted configurations; and in between 2470-2560 kpasses for fixing the triplet structures $g_{CM3}(r,s,s)$and $g_{ET3}\left(r,s,s\right)$, where $r=s$ is included, gathering statistics every 8000-10000 accepted configurations. Statistical errors (one-standard deviation) are fixed with block subaverages and remain very small: at the main peaks they are below 0.8% (see main details in Figure 1 and Figure 3).

The calculations of the triplet structures in k-space turn out to be much more involved, and a few remarks are to be made. The PIMC simulations involve the sample sizes $\left(N_S\times P\right)$: (a) $\left(128\times 66\right)$ for helium-3 at $T=4.2\mathrm{K}$(Figure 5); and (b) $\left(250\times 12\right)$ at the QHS fluid state point $\left({{\lambda }_B^{*};\rho }_N^{*}\right)=$ $\left(0.6;0.589\right)$(Figure 7). These calculations are focused only on the equilateral components $S^{\left(3\right)}\left(k,k,{\displaystyle \frac{\pi }{3}}\right):$ instantaneous for helium-3, and instantaneous and centroid for the QHS state point. Twelve sets of ${\left\{{\mathbf{k}}_1,{\mathbf{k}}_2\right\}}_i$ vectors commensurate with the simulation box, are scanned; each set is composed of 8 pairs such that: $\left|{\mathbf{k}}_1\right|=\left|{\mathbf{k}}_2\right|=k_i,$ the vectors being at an angle $\pi /3,$ and with the constraint $2\le k_x^2+k_y^2+k_z^2\le 200$ [58]. This implies that for helium-3 at $T=4.2\mathrm{K}$the resulting moduli (wavenumbers) are in $0.5\lesssim k_i/{Å}^{-1}\lesssim 5$, whilst for the QHS fluid state point, the moduli are in ${2.3\lesssim k}_i^{*}\le 11.8$ $\left(k_i^{*}=k_i{\sigma }_{HS}\right).$ Also, the sampling is more focused on the low-k wavenumbers and the regions of the main peaks. PIMC run lengths are as follows: for helium-3 in between 30 and 120 Mpasses, gathering statistics every 5000 accepted configurations; and for the QHS fluid state point in between 6.5 and 60 Mpasses, gathering statistics every 8000 accepted configurations. Statistical errors (one-standard deviation) are fixed with block subaverages and remain controlled: at the main amplitudes they are: 1.8% for the instantaneous case in helium-3; $4.8\%$ for the centroid case and $4\%$ for the instantaneous case in the QHS fluid (see main details in Figure 5 and Figure 7)). As a complementary result, the (mean) imaginary components, which must be identically zero [41], at the main amplitudes turn out to be: $0.0045\pm 0.0114$ for the instantaneous case at the helium-3 state point; $-0.0218\pm 0.4219$ for the centroid case and $0.0485\pm 0.2518$ for the instantaneous case at the QHS fluid state point. Clearly, although the QHS sampling is highly informative, it must be improved. In addition, all these results give an idea of the slowness of this sort of calculations when dealing with fluids with quantum behavior. This drawback becomes more acute if the fluid state point investigated belongs to the actual crystallization line of the potential model selected to describe the fluid. In this case, the $N_S$ sample size cannot be small (i.e., about one hundred particles, as the used in helium-3): the “flipping” back and forth in the simulation between the fluid and solid branches [13] must be avoided. This is the case in the current QHS fluid application. After monitoring such an effect, it has not been observed with $N_S=250,$but it is rather conspicuous with $N_S=128$. This monitoring has involved the configurational centroid CM2 structure factor $(\mathbf{k}\ne 0)$ defined as normalized to unity at its crystal-structure maxima [137]. For the current QHS study with $N_S=250,$ the latter configurational quantity remains below the maximal fluid value estimate $0.2$ [61].

Closure calculations reported here for triplets (instantaneous and centroid) in r-space are carried out utilizing KS3, JF3 and AV3, whereas in k-space are with BHP3 [41,55,56] and DAS3 [65]. These calculations necessitate the corresponding pair structures, $g_2(r)$ and $S^{\left(2\right)}(k)$. Detailed descriptions of the related procedures can be found in works by this author [54,55,56,57,58,59,60], although for the reader to better grasp the reported results it may be worthwhile to mention the following facts.

(i) JF3 and BHP3 share a common type of integral [39,41,156], which in each case is utilized in this work as follows:

- JF3 (r-space) [54]

$$I\left(r_{12},r_{13},r_{23}\right)=\int d{\mathbf{q}}_4{h}_2(r_{14})h_2(r_{24})h_2(r_{34})\approx$$

$$\pi {\sum }_{\nu =0}^{{\nu }_{max}}(2\nu +1)P_{\nu }\left(\cos \theta \left[{\mathbf{r}}_{13},{\mathbf{r}}_{23}\right]\right){\int }_0^{\mathrm{\infty }}dy{y}^2{h}_2(y)f_{\nu }(y,r_{13})f_{\nu }(y,r_{23}),$$

(116a)

$$f_{\nu }\left(y_1,y_2\right)={\int }_{-1}^{+1}dx{P}_{\nu }\left(x\right)h_2\left(\sqrt{y_1^2+y_2^2-2x{y}_1{y}_2}\right),$$

(116b)

- BHP3 (k-space) [55]

$$c_{BHP}^{\left(3\right)}\left(k_1,k_2,\theta \left[{\mathbf{k}}_1,{\mathbf{k}}_2\right]\right)={\displaystyle \frac{1}{{\left(2\pi \right)}^3}}\int d{\mathbf{k}}_{3}t\left({\mathbf{k}}_3\right)t\left({\mathbf{k}}_3-{\mathbf{k}}_1\right)t({\mathbf{k}}_3+{\mathbf{k}}_2)\approx$$

$${\displaystyle \frac{1}{{\left(2\pi \right)}^2}}{\sum }_{\nu =0}^{{\nu }_{max}}\left(\nu +{\displaystyle \frac{1}{2}}\right)P_{\nu }\left(\cos \theta \left[{\mathbf{k}}_1,{\mathbf{k}}_2\right]\right){\int }_0^{\mathrm{\infty }}dk{k}^{2}t\left(k\right)f_{\nu }\left(k_1,k\right)f_{\nu }\left(k_2,k\right),$$

(117a)

$$f_{\nu }\left(k_a,k\right)={\int }_{-1}^{+1}dw{P}_{\nu }\left(w\right)t\left(\sqrt{k_a^2+k^2-2w{k}_{a}k}\right),$$

(117b)

where the $P_{\nu }$ stand for the standard Legendre polynomials, and the upper limit ${\nu }_{max}=30$ is utilized. (In the limit${\nu }_{max}\to \infty$ the $\nu$-expansions are exact [41,156]).

(ii) For the closure analyses at constant temperature, the numerical density-derivatives of the pair $g_2(r)$ and the direct correlation functions, $c_2\left(r\right)$ and $c_2\left(k\right),$ are evaluated with the centered algorithms: Stirling’s (two-point) and Richardson’s (four-point) [182] (Figures 5 to 7). By denoting a generic radial pair function by ${\mathfrak{f}}_2\left(r;{\rho }_N\right)$, the previous algorithms can be written as:

$${{\left({\displaystyle \frac{\partial {\mathfrak{f}}_2\left(r;{\rho }_N\right)}{\partial {\rho }_N}}\right)}_T^{Stirl.}\approx D}_1\left(\mathsf{\Delta }{\rho }_N\right)={\displaystyle \frac{{\mathfrak{f}}_2\left(r;{\rho }_N+\mathsf{\Delta }{\rho }_N\right)-{\mathfrak{f}}_2\left(r;{\rho }_N-\mathsf{\Delta }{\rho }_N\right)}{2\mathsf{\Delta }{\rho }_N}},$$

(118a)

$${{\left({\displaystyle \frac{\partial {\mathfrak{f}}_2\left(r;{\rho }_N\right)}{\partial {\rho }_N}}\right)}_T^{Rich.}\approx D}_2\left(\mathsf{\Delta }{\rho }_N;2\mathsf{\Delta }{\rho }_N\right)={\displaystyle \frac{{4D}_1\left(\mathsf{\Delta }{\rho }_N\right)-D_1\left(2\mathsf{\Delta }{\rho }_N\right)}{3}}.$$

(118b)

(iii) As regards JF3 and AV3, the following observations are worth making. Comparison with PIMC results for the equilateral components of the helium-3 triplet structure factors $S_{ET}^{\left(3\right)}(k,k,{\displaystyle \frac{\pi }{3}})$ indicates that (Figure 5 and Table 1): (a) JF3 cannot capture the negative behavior within the low-k region and its usefulness reduces to the large-k wavenumbers, since it gives the asymptotic behavior of the structure factor; and (b) AV3 as reported in Reference [60] yields lower values for both the negative depth in the low-k region and the height of the main amplitude. In view of the results in Figure 6, it is highly likely that AV3 can be improved using Abe’s developments [38] and gives better triplet results in both r-space and k-space.

(iv) BHP3 minimizations for ET3 helium-3 at the state point at $T=4.2\mathrm{K}$ were reported in [60], where use of conjugate gradients [185] was made, as explained elsewhere [55,56]. Discretizations of the $t(r)$ functions into 7001 points covered the distance interval $\left[0\le r/Å\le 70\right]$. Convergence as assessed via ${\displaystyle \frac{\mathsf{\Theta }\left[t\left(r\right)\right]}{{‖\frac{\partial c_2(r)}{\partial {\rho }_N}‖}^2}}$ was rapid and without any problems, the latter quotient reaching values $~{10}^{-9}$ in a few hundred iterations. The explored ranges of wavenumbers for $S_{ET}^{\left(3\right)}(k,k,{\displaystyle \frac{\pi }{3}})$ have been extended, $k>5{Å}^{-1},$ for checking the properties in k-space as k increases: JF3 appears as the limit with increasing k wavenumbers, and the results are consistent. Figure 4 and Figure 5 and Table 1 collect representative results for this system. As seen, for the equilateral components, BHP3 [60] gives a nice representation of the PIMC values, particularly in capturing the negative behavior for $k\lesssim {1Å}^{-1}.$

Nevertheless, some pilot BHP3 minimizations for the QHS fluid (CM3 and ET3) along the crystallization line [79,137] have turned out to be affected from convergence problems. These latter results are therefore far from conclusive from a physical standpoint and only a brief comment follows. The key observation is that the convergence in this case (tested with both gradients and conjugate gradients) appears to be strongly dependent on the value of the norms ${‖\partial c_2/\partial {\rho }_N‖}^2,$ in such a way that the larger this value is, the slower and more physically doubtful the convergence path turns out to be. These norm values are clearly associated with the depths at the origin of the pair direct correlation function $c_2\left(r=0\right)$. To get a better feeling of this situation some data related to the pair direct correlation functions analyzed in this work are worth quoting here:

(i) For helium-3 at ($4.2\mathrm{K};$ $0.02286713{Å}^{-3})$the five instantaneous values $c_{ET2}\left(r=0\right)$, corresponding to the five state points involved, are in $-9.196\lesssim c_{ET2}\left(r=0\right)\lesssim -5.051$ [60], and no BHP3 problems with the convergence were observed.

(ii) For the QHS fluid at $({\lambda }_B^{*}=0.6;{\rho }_N^{*}=0.589)$: the three instantaneous values involved $(\mathsf{\Delta }{\rho }_N^{*}=\pm 0.01715)$ at the origin of the interparticle distances are in $-36.764\lesssim c_{ET2}\left(r^{*}=0\right)\lesssim -30.916,$ while the three centroid values are in $-58.720\lesssim c_{CM2}\left(r^{*}=0\right)\lesssim -46.598.$

(iii) The reference values given by the squared norms ${‖\partial c_2/\partial {\rho }_N‖}^2$ turn out to be as follows:

-

For helium-3, using Richardson four-point derivative$\left({\displaystyle \frac{\mathsf{\Delta }{\rho }_N}{{Å}^{-3}}}=\pm 0.002,\pm 0.004\right),$ the instantaneous reference is ${‖\partial c_{ET2}/\partial {\rho }_N‖}^2\approx 0.4845\times {10}^7{Å}^9$.

-

For the QHS fluid along the crystallization line [79,137], covering the conditions mentioned above compatible with $0.2\le {\lambda }_B^{*}\le 1.9832$ $\left({\sigma }_{HS}=3.5Å\right),$and using Richardson four-point derivative $(\mathsf{\Delta }{\rho }_N^{*}=\pm 0.0042875,\pm 0.008575)$, one finds values for the instantaneous reference in ${0.3050\times {10}^9{Å}^9\lesssim ‖\partial c_{ET2}/\partial {\rho }_N‖}^2\lesssim 0.2598\times {10}^{10}{Å}^9,$ and for the centroid reference in ${0.2881\times {10}^{10}{Å}^9\lesssim ‖\partial c_{CM2}/\partial {\rho }_N‖}^2\lesssim 0.2094\times {10}^{11}{Å}^9.$

As seen, there is a significant BHP3 difference between the foregoing results for helium-3 and the QHS fluid. Note that the comparison is better made using the actual values in $Å$ for distances, since the modelling via hard spheres for helium-3 is not consistent [55] for the current structural purposes. In relation to this, the distinct features of the interparticle potentials $u\left(r\right)$ involved are to be remarked: for helium-3, continuous for $r>0$ and weakly attractive at the minimum $(\approx -11\mathrm{K})$ [106], whilst for QHS there is an infinite discontinuity at $r={\sigma }_{HS+}$. These potentials are obviously connected to the shapes of the pair direct correlation functions, and hence to their density derivatives. Consequently, although applications of BHP3 are clearly dependent on the latter derivatives by reason of their construction Equation (25b), in intricate regions such as the freezing transition this dependence becomes extremely sensitive to the quality of the numerical derivatives. (Recall the basic density increments $\mathsf{\Delta }{\rho }_N{(Å}^{-3})$ about the target state points: $0.002{Å}^{-3}$ for helium-3, and $0.0004{Å}^{-3}$ for QHS). Therefore, a very careful fixing of the QHS intermediate $c_{ET2}$ and $c_{CM2},$ by diminishing $\mathsf{\Delta }{\rho }_N$ for increasing the accuracy in the density derivative calculations, could fix this BHP3 situation. However, apart from the computational overload, possibly including an augmented precision in the calculations, one can still foresee some difficulties. Firstly, the slow minimization of a large quantity $\mathsf{\Theta }\left[t\left(r\right)\right]$ which, together with the fine details of the density derivatives of the pair direct correlation function involved, might misdirect the convergence path and lead to unphysical solutions. Secondly, although OZ2 for centroids is formally exact, OZ2 is an approximation for the instantaneous case, which means for the latter that extra problems can be encountered in the freezing region when applying simple corrections such as BHw. Furthermore, there is the question of the complete calculation of the triplet structure factors by including its general components. Obviously, detailed work to clarify the ranges of applicability of BHP3 is needed.

(v) On the contrary the current DAS3 calculations of $S^{\left(3\right)}(k,k,{\displaystyle \frac{\pi }{3}})$ are directly worked out in k-space and, once the necessary pair information $c_2\left(k\right)$ is available, they are straightforward. Apart from the numerical derivative problems, no intrinsic procedural problems appear. As shown in Figure 5 and Figure 7 and Table 1, DAS3 also gives good representations of the equilateral PIMC results, albeit apparently not so good as BHP3 in the instantaneous low-k region$.$ These applications use the two-point derivative Equation (118a) for the QHS-fluid CM3 and ET3 evaluations, and the four-point derivative Equation (118b) for the helium-3 ET3 evaluations. For this closure the underlying OZ2 problems with the instantaneous class also exist, and its basic assumptions Equations (26) remain to be fully tested. Again, detailed work to clarify the ranges of applicability of DAS3 is needed.

7. Final Remarks and Future Directions

The investigation of the static triplet structures in quantum condensed matter (fluids and solids) via path integrals is a promising and challenging avenue. Although it may be said that an initial background to this topic is already established, key theoretical and computational aspects in both the real (r-space) and the reciprocal Fourier (k-space) remain to be explored and understood. The theoretical issues should focus on: (a) expanding further the limits set by the usual pair-level approach (solid phases [59]), and (b) the treatment of systems composed of nonzero-spin particles, particularly fermionic systems. The computational issues should mainly address the questions of cost-effectiveness posed by the required accuracy when determining the variety of different quantum triplet structures (dimensionality $\ge$ 4-D), with a view to increasing the simulation $N_S$ sample size and, also, the possible integration of proper two- and three-body potentials (e.g., [107,108]) in the structure computational schemes. In this regard, exascale computations, when widely available, could alleviate the effort to be made. (Just as a high hope, if quantum computers became a reality in practice some day, the quantum triplet problems might provide an excellent test bed). All these developments can provide very useful guidance for materials research and future higher-order structural studies.

Some specific theoretical goals within the path integral framework can be the following. (i) The consideration of general external fields including magnetic interactions and/or situations involving more than one spin state in monatomic fluids; the latter can give rise to a special set of treatments for fermionic fluids when studied with the insightful Wigner path integral schemes [98,99]. (ii) A formulation of the actual triplet structures in molecular fluids (e.g., hydrogen fluids) that, going beyond the overall one-site picture [57,58], gives full meaning to the direct atom-atom constructions [11], or yields more global molecular approaches.

As regards the specific computational goals, one can mention the following pending path integral tasks in homogeneous and isotropic monatomic fluids (4-D problems). In the diffraction regime: (iii) the determination of the total continuous linear response functions TLR3, which remain virtually untouched [54]; and (iv) the completion of the ET3 and CM3 structures by fixing their behaviors over significant ranges of the three significant variables which characterize them (i.e, three independent interparticle distances, or two wavenumbers and an angle). In the quantum exchange regimes, there is nothing yet related to calculations of triplet structures that has been undertaken and work on this should be welcome.

As complements to the exact path integral calculations of triplets, one may try a variety of closures. These closures are not as efficient and clear-cut defined as those at the pair level. Despite this remark, there are some results indicating that, when studying the instantaneous and centroid $g_3$ and $S^{(3)}$ functions within the quantum diffraction regime, closures may be a great source of information. Triplets $g_3$ in r-space have been discussed at length elsewhere [54,55,56,57,59,60], and suffice it to say that the results obtained through the AV3 closure Equation (24) [59,60] point to the possibility of achieving significant improvements based on Abe’s expansion [38]. In this regard, it is expected that the negative effects of increasing densities on the long-range isosceles correlations could be fixed to a significant extent, which may reflect in better k-space results. Nonetheless, closure $S^{(3)}$ results for quantum triplets in k-space remain comparatively an unknown issue. For the time being, the closures BHP3 Equations (25) and DAS3 Equations (26) give results $(S^{(3)}$ equilateral components) that make them worth exploring in depth. They may capture traits of the sensitive low-k (negative) region of the equilateral ET3 and/or CM3 structure factors and, also, the behavior of the important main amplitude region. In connection with the low-k region, BHP3 seems better adapted than DAS3. However, and based on pilot calculations, BHP3 variational calculations seem less appropriate for investigating intricate regions, such as the quantum hard-sphere crystallization line, whereas DAS3 calculations appear here as a more robust option for both ET3 and CM3. For large-k wavenumbers the basic JF3 closure (Equation (15) with $c_{JF}^{\left(3\right)}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)=0)$, which is part of both closures, dominates the triplet $S^{(3)}$ behavior.

The k-space closure results obtained so far are just a beginning (e.g., [58], and Figure 5 and Figure 7 herein), and more work is needed to substantiate their relevance, particularly regarding the non-equilateral components. Searching for the limits of applicability of closures, even for qualitative purposes related to the CM3 and ET3 cases in the diffraction regime, seems a rewarding task. The case TLR3 does not seem amenable to closure studies because of self correlations, although there is nothing against trying closure treatments of the pure triplet-atom part ([54] and Equations (113)).