A Glimpse into Quantum Triplet Structures in Supercritical 3He

A methodological study of triplet structures in quantum matter is presented. The focus is on helium-3 under supercritical conditions (4 < T/K < 9; 0.022 < ρN/Å−3 < 0.028), for which strong quantum diffraction effects dominate the behavior. Computational results for the triplet instantaneous structures are reported. Path integral Monte Carlo (PIMC) and several closures are utilized to obtain structure information in the real and the Fourier spaces. PIMC involves the fourth-order propagator and the SAPT2 pair interaction potential. The main triplet closures are: AV3, built as the average of the Kirkwood superposition and the Jackson–Feenberg convolution, and the Barrat–Hansen–Pastore variational approach. The results illustrate the main characteristics of the procedures employed by concentrating on the salient equilateral and isosceles features of the computed structures. Finally, the valuable interpretive role of closures in the triplet context is highlighted.


Introduction
Quantum triplet structure studies are key to developing the statistical mechanics of equilibrium many-body systems at low (nonzero) temperatures. Despite the fact that triplet structures are not determined experimentally [1], the computational approach to this topic not only justifies itself by the gaining of knowledge about quantum matter, but it also is crucial for further applications of quantum reasoning (e.g., phase transitions and design of materials, phonon-phonon interactions in superfluids, time-dependent phenomena, etc.) [2][3][4][5].
As stressed elsewhere [6], exact quantum triplet calculations are presently an extremely demanding task. This contrasts sharply with their counterparts in the classical domain where the calculations are far more affordable [7][8][9][10]. Putting aside the high dimensionality that the triplet functions can reach (e.g., 10-D for spatial triplets in a monatomic solid), one should note the quantum variety of physically significant n-particle structures that a system can exhibit in both the real space (r-space) and the reciprocal Fourier space (k-space) (see [11,12] for a basic description). For the reader to grasp the overall situation, suffice it to consider a monatomic homogenous and isotropic fluid with substantial quantum effects that make a classical description meaningless. This fluid shows six basic triplet functions that can be classified into three types, namely, instantaneous, total continuous linear response, and centroids (the classical counterpart only has two basic functions in a single class) [5,6,11,12]. Each of these three types contains one generalized triplet correlation function H 3 (r 12 , r 13 , r 23 ) together with its related Fourier transform S (3) (k 1 , k 2 , cos(k 1 , k 2 )), where r jm = q j − q m is the distance between particles j and m, and k 1 and k 2 are two wavevectors of moduli k 1 and k 2 , respectively. As seen, these triplet functions are already 4-D.
Accordingly, one can obtain an initial impression of the magnitude and expected cost of the related quantum computations, which would increase if higher-order structures were to be dealt with. However, this impression becomes more acute when the inherent features of the actual computations enter the discussion. Very powerful methods to calculate most of the properties of equilibrium quantum condensed matter are based on Feynman's path integrals (PI) [13]. Within PI two main simulation techniques are available: path integral Monte Carlo (PIMC) and path integral molecular dynamics (PIMD) [11,[14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30]. In a way similar to the classical MC and MD techniques, the quantum PIMC and PIMD are highly accurate, their procedural "errors" (e.g., statistical, numerical) being diminished by increasing the simulation run lengths and/or the precision of the computations. By focusing on definiteness for the PIMC simulations of quantum equilibrium structures, they may also be regarded as "exact" in that they provide self-contained solutions unattainable via basic frameworks in statistical mechanics. In relation to this, one should recall the case of the exact (and theoretically revealing) Bogolyubov-Born-Green-Kirkwood-Yvon hierarchy (BBGKY), which needs the knowledge of higher-order structures to define lower-order structures [31]. Thus, although BBGKY is an exact formulation, it leads to calculational schemes (not in use) that cannot be "exact", since they need extra information to break/close the working equations (e.g., hypothesis about the form of triplet correlations -closures-to calculate the pair correlations). This does not occur with PIMC, which is self-contained (e.g., pair and triplet correlation functions can be calculated independently of each other) and whose accuracy can in principle be arbitrarily increased so as to reproduce the targeted theoretical values of the model selected [14,15,18]. Note that the foregoing polysemous use of the term "exact" is independent of another use referring to the computational effort required to deal with an increasing number of particles in a system.
The PI difficulty lies in the extended simulation samples N S × P that PIMC and PIMD use; N S stands for the conventional number of actual particles, and P > 1 is an integer number to be optimized that serves to represent the thermal quantum delocalization of an actual particle (theoretical accuracy is reached in the Trotter's limit P → ∞) [18]. As a rule, P increases with the quantum effects and there are ways to soften its impact on the number of calculations (e.g., pair actions [18,19], fourth-order propagators [20,[25][26][27][28][29][30], parallel computing [20], etc.). All in all, when very strong quantum effects (including bosonic exchange [18,23,28]) are to be studied, the PI undertaking of the triplet task in its entirety, and within reasonable CPU times (and electric power consumption), may remain today out of the reach of most interested researchers. In this regard, note that triplet r-space information, though expensive, is still affordable, whereas its counterpart in k-space is highly demanding because of the necessity to scan appropriate sets of k wavevectors commensurate with the simulation box. It thus seems that the related numerical evaluations of detailed quantum r-and k-structures for triplets (and beyond) could be appealing targets for the PI implementations in the coming exascale computers [32]. For completeness, note that, apart from the present author's works, PI work on significant aspects of the equilibrium structures in quantum matter can also be found in the general references [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30], more specifically [14,16,18,20,[27][28][29][30]. However, none of the latter deals with all the aforementioned types, nor goes beyond the pair level.
For the time being, the main structural features of triplets in fluids with quantum behavior can be extracted by means of PIMC (or PIMD). Recent works by the present author have shown how to tackle this problem when dealing with quantum diffraction effects [2,6,[33][34][35][36]. In these works, one can find comprehensive descriptions of the methods employed, and also positive identifications (or well-grounded indicators) of physically significant triplet patterns in the statistical distributions of the actual particles. The systems studied were the quantum hard-sphere fluid, bare [2,33] and with Yukawa attractions [35], helium-3 at very low densities [34], liquid neon [6], and liquid para-hydrogen [6,36]. By defining the parameter γ = ρ N λ 3 B as a (rough) measure of the magnitude of the quantum effects, where ρ N is the number density, and λ B = h/ √ 2πmk B T is the thermal de Broglie wavelength, the studied conditions covered up to γ 2.7.
As complementary tools, one can also employ the so-called closures, which for triplets intend to infer their characteristics from the available information at the pair level: g 2 (r) correlations, c 2 (r) direct correlation functions, S (2) (k) structure factors, and other auxiliary functions [3,5,37]. Closures are certainly approximations and imply far less expensive calculations than the exact PI techniques. Furthermore, closures may turn out to be highly accurate, as shown at the pair level [11,38], or, alternatively, very useful as interpretive tools for analyzing complex structural problems [2,6,36]. From the point of view of the present author, any theoretical object allowing further reasoning and a deeper understanding always deserves careful study, and this turns out to be the case of closures in the quantum domain.
In the hope of shedding some more light on quantum triplet structures, and as a preliminary part of a larger project, this article addresses some relevant issues. The focus is on the instantaneous triplet structures in the r and the k spaces of supercritical helium-3.
The triplet behavior is explored for conditions 4 < T(K) < 9; 0.022 < ρ N Å −3 < 0.028 (critical point: [39]. Why helium-3 under these conditions? One obvious reason is to give a service by communicating more experience on the, as yet, unexplored triplet topic. Another reason is that quantum diffraction effects (γ 3.2) beyond those mentioned above can be analyzed in (a model of) a real system as important as helium-3. Furthermore, the rationalization of the quantum triplet structures in terms of closures built from the underlying pair structures is worth pursuing [2,6,36]. The supercritical conditions are far from fermionic exchange, known to be present for T ≤ 1K and affected by the computational "sign problem", which precludes practical applications of PIMC (see [17] for a pioneering PIMC approach to this problem). In this article, the equilateral and isosceles correlations are determined in r-space with PIMC and closures, and in k-space with closures. PIMC involves the fourth-order propagator put forward in [25][26][27] (compare with the early application of the primitive propagator that was utilized in [34] for the study of helium-3 at very low densities). The triplet closures employed are Jackson-Feenberg convolution JF3 [3], Kirkwood superposition KS3 [37], the intermediate AV3 = (KS3 + JF3)/2 [2], and the variational Barrat-Hansen-Pastore approach (BHP) [5]. KS3, JF3, and AV3 are utilized for r-space and k-space, whereas BHP is utilized only for k-space. The effects in r-space arising from changes in temperature and in density are discussed, and the significant role of the closures is highlighted.
The outline of this article is as follows. Section 2 contains a summary description of the underlying theory and methods. Section 3 is devoted to the main computational details, and Section 4 gives the results and their discussion. Finally, Section 5 collates the conclusions of this work.

Path Integral Monte Carlo
For a monatomic homogeneous and isotropic fluid at equilibrium, in which quantum exchange can be neglected, the general form of the PI-canonical (N, V, T) partition function reads as [14,18] where m is the particle mass, β = 1 k B T is the inverse temperature, the actual N particles (j = 1, 2, . . . , N) are represented by N necklaces with P beads apiece (t = 1, 2, . . . , P), r N,t = r t 1 , . . . , r t N where r t j denotes the coordinates of bead t belonging to necklace j, dr N,t = dr t 1 . . . dr t N , and W NP is the effective potential for the whole set of N × P beads. Equation (1) resembles the form of a semiclassical partition function [11,13], which gives rise to the so-called (semi-)classical isomorphism [14] that is at the root of the great PI success. The number P is assumed to be optimized in what follows.
The specific expression of W NP depends on the selection of the thermal propagator. In this work this propagator is the fourth-order choice in the final form given by Voth et al. to the Suzuki-Chin developments [25][26][27]. In this case P is an even integer, and W NP is given by where the * in the first t-summation implies the cyclic property t + 1 = P + 1 ≡ 1, α is a real number in the interval [0, 1], v r t jm is the continuous pair potential v(r) acting between equal-t beads in different necklaces, r t jm = r t j − r t m , η t jm = r t jm /r t jm , and the asymmetry in the bead contributions coming from odd-t (t = 1, 3, . . . , P − 1) and even-t (t = 2, 4, . . . , P) is to be noticed. This asymmetry has an obvious impact on the thermodynamic evaluations, but it is of critical importance to the structural computations, for only the odd-numbered beads are physically significant. The reader is referred to [11,27,28] for the related discussions.
By focusing only on the instantaneous structures (ET for equal-t) at the fluid number density ρ N = N/V, one finds at the pair and triplet levels the following direct canonical averages . . . in the r-and k-spaces In connection with the foregoing formulas, some remarks are to be made: (a) q j stands for the absolute position vector of the actual particle j, so that r jm = q j − q m is the distance between the actual particles j and m; (b) in Equations (5) and (6) terms of the δ(k)− type are omitted, as is customary [40,41]; (c) γ is the angle between the wavevectors k 1 and k 2 ; and ET is the Fourier transform of an involved r-space function, H ET3 (r 12 , r 13 , r 23 ) [6,11]. The latter not only contains g ET3 (r 12 , r 13 , r 23 ), but also additional terms which include two-and one-body contributions. It is worth stressing that Equations (3)-(6) are the common quantities in the structural study of fluids at equilibrium. Actually, S ET (k) is directly related to the system response that can be obtained in X-ray or elastic neutron diffraction experiments [18,31]. However, for many purposes, the grand canonical extensions are essential, as will be considered later on (see [6,11] for the other functions and their physical meanings).
There are well-known problems with the simulation of Equations (5) and (6) in the low-k regions, caused by the finite particle sample size N S (or volume V S ) and the required commensurability of the wavevectors with the central box [40,41]. For a cubic box of length L the allowed wavevectors must be of the form k = 2π L k x , k y , k z with the k x , k y , k z components being integer numbers, and a summary description of these problems seems worth including here (see [11] for a review and complete discussion).
First, for example, one may focus just on the pair level for simplicity and also for its immediate repercussions. Note that the central thermodynamic connection S cannot be strictly determined via a single direct simulation (χ T = isothermal compressibility), since k and N S extrapolations should be carried out [41]. (Due to the lack of a δ(k)-term, the case k = 0 in Equation (5) yields an expression far from the true physical meaning of this component [11,42]). Further, the long-range oscillations about zero of h ET2 (r 12 ) cannot be determined with reasonably workable finite sample sizes (unless one had access to computational resources whose use could be hardly justifiable). Although, these oscillations are small in magnitude, they are decisive to obtain a significant value of S ET (k = 0) via the Fourier transform of h ET2 (r 12 ). Nothing of this is new, as this sort of problem was already known in classical statistical mechanics, and was circumvented via the use of the Ornstein-Zernike pair scheme (OZ2) with the introduction of the short-ranged direct correlation function c(r 12 ) [40,41]. One might argue that the finite sample size affects the calculation of every property in a simulation and that there is nothing about the S ET (k = 0) evaluation, via the volume integral of h ET2 (r 12 ), that is not shared by other properties, such as the energy or the pressure, which can also be formulated with the use of integrals involving the pair function g ET2 (r 12 ) (and the interparticle potential) [18,41]. However, it is worth recalling that, when evaluating energies or pressures of systems composed of electrically neutral particles, the standard way to cope with the finite size effects is to correct the simulation results with continuum contributions. The latter can be evaluated by defining the pair function involved as g 2 (r 12 ) = 1 for distances longer than half the box length L/2. These corrections account for the missing long-range interactions, thereby being fully dependent on the features of the interparticle potential [41], but not on the complete structure of the actual system studied. As stated above, and apart from well-known asymptotic behavior questions [42], the latter definition, g 2 (r 12 ) = 1, would be meaningless when Fourier transforming h 2 (r 12 ), and the corresponding argument is flawed.
Second, the sampling of commensurate wavevectors implies the definitions of sets {k uv } v=1,2,3,... compatible with the wavenumbers k u = |k uv | to be analyzed (as many sets as wavenumbers). Therefore, and focusing on the use of (6) for triplets, the increased computational cost involved prompts one to look for reliable alternatives. In this regard, there are methods [11] based on the use of Ornstein-Zernike (OZn) schemes [43] that, while in waiting for a widespread availability of exascale computers, can be explored in the quantum domain for dealing with these serious triplet drawbacks. The OZn schemes are deeply rooted in the grand canonical ensemble developments in both the classical [43] and the quantum [11,14] domains and utilize the aforementioned closures. They are known to produce highly accurate quantum frameworks at the pair level for the three types of pair S (2) (k) [12,38] (for PI centroids the chain of OZn schemes is exact [34]). Now, before entering the discussion on closures, three outstanding aspects of the quantum OZ2 treatments deserve to be mentioned. First, the high accuracy obtainable over the whole range of wavenumbers k ≥ 0, even with the use of moderately sized PIMC samples that provide the basic pair structures g 2 (r). Second, the excellent agreement with experiment that they produce, and third, their cost-effective character, since they only require a very low computational cost as compared to PI simulations of the structure factors. For more information on these important conceptual and practical subjects, the reader is referred to [2,6,11,[33][34][35][36]38].

Closure Approximations
The closures in r-space for the instantaneous triplet correlations employed in this work are the following: Kirkwood superposition (KS3) [37], Jackson-Feenberg convolution (JF3) [3], and their average AV3 = (KS3 + JF3)/2 [2]. By making use of the following conventions: g ET2 r jm → g r jm , and h ET2 r jm = g ET2 r jm − 1 → h r jm = g r jm − 1, these closures can be cast as g KS3 ET3 (r 12 , r 13 , r 23 ) = g(r 12 )g(r 13 )g(r 23 ), g JF3 ET3 (r 12 , r 13 , r 23 ) = g KS3 ET3 (r 12 , r 13 , Equation (9) shows the form of a truncated h-expansion in which the triple h-product, subtracted from JF3, is recovered. Surprisingly, this form is more useful than expected, as shown in a recent study of the quantum hard-sphere fluid along the crystallization line [2].
The closures in k-space for the instantaneous calculations contained in this work are based on the OZ direct correlation functions c ET2 (r 12 ) and c ET3 (r 12 , r 13 , r 23 ), together with their corresponding Fourier transforms c  ET (k 1 , k 2 , γ) [11]. These schemes lead to the following classical-like recipes for the pair and triplet structure factors Equations (10) and (11) are approximations to the instantaneous case within a grand canonical ensemble description. However, the pair level Equation (10) is highly accurate for fluids with substantial quantum behavior [11]. Moreover, the latter assertion is true for the fluid phases of helium in which quantum exchange is negligible. This was demonstrated for 3 He and 4 He at 4.2K in [38], where one can observe the almost indistinguishability between the OZ2 and the exact PIMC results. It is worth remarking that the calculations in r-space and in k-space can be interconnected through closures, which allows one to achieve excellent approximations to the grand canonical ensemble structures without performing simulations in this ensemble. Thus, by starting from the PIMC canonical pair structures g ET2 (r 12 ) and subjecting them to an extended OZ2 treatment, one can obtain a set of improved functions, g ET2 (r 12 ), c ET2 (r 12 ), and S (2) ET (k), which incorporate grand canonical corrections. The latter contribute to correct the deficient asymptotic behavior of the canonical pair radial correlations [42]. Hereafter, this method will be referred to as BDH+BHw (Baxter-Dixon-Hutchison procedure plus Baumketner-Hiwatari corrections) [44][45][46], of which there is plenty of information regarding quantum applications (see [11] and references therein). The calculation of S (2) ET (k = 0) is just a by-product that arises directly from this approach and provides very valuable estimates of the isothermal compressibility. (As a matter of fact, the exact value of the latter quantity for the model system under study can be obtained via application of OZ2 to the centroid functions [12]).
In addition, the complete knowledge of S (2) ET (k) allows one to extend, through the inverse Fourier transform, the pair radial correlations g ET2 (r 12 ) up to arbitrarily long distances. Clearly, the latter operation benefits the closure calculations of triplets [11]. Furthermore, the knowledge of c ET2 (r 12 ) over a range of densities, at constant temperature, is key to tackling the far more intricate problem posed by c ET (k 1 , k 2 , γ) depends heavily on the closure selected [5,10,[33][34][35][36]. Apart from the "trivial" JF3 solution c (3) ET = 0, which makes JF3 a compulsory reference, and the KS3 solution [5], the present work makes use of the elaborate variational procedure proposed by Barrat, Hansen, and Pastore (BHP) in the classical domain [5]. Therefore, the calculations below are based on the following equations For simplicity of notation, the dependence of the direct correlation functions on ρ N is included only in c 2 (r). The auxiliary function t(r) defines the BHP closure by cutting Baxter's c n -hierarchy [47] at the triplet level, Equation (12a). The fixing of t(r) is performed through the functional minimization of [t(r)] ( R max is an upper limit for the integrations in Equation (12c)). Once t(r) is obtained, the double Fourier transform of c 3 leads to S (3) (k 1 , k 2 ) [5]. Also, note the exact relationship arising from Equation (12a) The BHP applications to the ET2 case make the replacements of the pair and triplet (BHP classical) quantities by their quantum instantaneous analogs. By taking advantage of Equation (12d), approximate estimates of the double-zero momentum transfer ET (k 1 = 0, k 2 = 0) can be obtained. (In actual fact, the use of the exact centroid function framework would yield the exact value of the latter component) [6,34]. Although Equation (12b) is already an approximation to the exact classical behavior of c 3 , BHP is known to capture some interesting traits of quantum triplet structure factors [36]. Whether or not BHP may yield quantum triplet results of a quality comparable to that achievable at the pair level with BDH+BHw, only further applications and comparison with PI results will settle the question. Therefore, BHP is utilized in this work for continuing the pilot exploration of such topic when strong quantum diffraction effects are involved.

Computational Details
The helium-3 state points studied are: . The atom mass of helium-3 is set to 3.01603 amu. Conditions SP1 to SP3 are taken from [48], allowing the comparison of the r-space triplet instantaneous results under independent variations of temperature and density. The selection of SP4 is made to carry out the triplet instantaneous calculations in k-space with closures. SP4 and its adjacent states along the 4.2 K isotherm were studied in [38], where their c ET2 (r 12 ), c (2) ET (k), and S (2) ET (k) functions were obtained via BDH+BHw. For purposes of interpretation of the results, the SP1-SP4 closeness is an advantage worth exploiting in this work.
The PIMC simulations follow the general lines already described in other works by the present author [33,38], and only a summary is given here. The interatomic potential employed is SAPT2 [49,50], which produces very reliable results for this system [11,38]. Consequently, SAPT2 can be regarded as adequate for the present structural purposes. The canonical ensemble is used for the basic r-space calculations, involving the fourthorder propagator α = 1 3 [27], and with the sample sizes: N S × P = 1372 × 66 (SP1), 1372 × 80 (SP2), and 1372 × 22 (SP3) (state point SP4 was studied with 1024 × 66 in [38]). The necklace normal-mode algorithm [51] is utilized, and the usual Metropolis sampling procedure is applied by setting the acceptance ratio for the different P-moves to 50%. One kpass is defined as 10 3 N S × P attempted bead moves. As stated earlier, the canonical pair instantaneous structures g ET2 (r) are needed to undertake the calculations of the triplet instantaneous closures, and the run lengths to obtain these g ET2 (r) are in between 500 kpasses and 2000 kpasses. The triplet instantaneous structures computed, g ET3 (r, s, s), are fixed with run lengths in between 750 and 3660 kpasses. The sampling of the pair and triplet structures uses a spacing in the interparticle distances set to ∆ g = 0.1Å [33]. The statistical error bars remain controlled: for example, at the first peaks of g ET2 (r) one finds the error bars (one standard deviation) well below 1%, and at the first peak of the equilateral g ET3 (r, r, r) one finds that the error bars remain ≤ 1% (see the Supplementary Material). The current applications using the PI fourth-order propagator employed [27] should be compared to those of the primitive propagator reported more than a decade ago in [34], where gaseous helium-3 was studied at 5.23 K and very low densities <0.0021Å −3 with sample sizes having: N S = 108, 500, P ≤ 130. The reduction in P and the possibilities for increasing N S (or equivalently, for analyzing wider ranges of the r-correlations) are powerful advantages offered by this efficient propagator when studying increasing densities.
Real space triplet calculations with closures use as data input the improved PIMCg ET2 (r) functions, which are extended up to distances longer than half the box-length L/2 with the use of S (2) ET (k). At a given state point, after calculating the PIMC-g ET2 (r), this task is accomplished in three steps: (1) application of Baxter-Dixon-Hutchinson's treatment (BDH) of the OZ2 equation [44,45]; (2)   ET (k). (For details see [11,38,52,53]). The crucial point is that c ET2 (r), which is short ranged, is fixed over a finite range R Z of distances: c ET2 (r ≥ R Z ) = 0. In this regard, there may appear more than one R Z value (hereafter R Z -zeros), for which the main part of c ET2 is kept essentially invariable, but obviously they yield different tails for the c ET2 decay towards zero with increasing r [53]. These tails only have (generally) a small effect on S (2) ET (k) in the region of very low-k values (see [38] for noticeable exceptions), and an effective averaging method has been proposed to deal with this situation [53,54]. However, this tail effect may or may not become important when the density derivatives involving isothermal sets {c ET2 (r; ρ N ; R Z )} must be computed, and the results below illustrate this point. Also, closure results for KS3 Equation (7) are trivial, but those for JF3 and AV3 depend on the convolution integral shown in Equation (8), which contains the total correlation function h(r) ≡ h ET2 (r). It is worthwhile to mention that, in evaluating this convolution, use of a well-known expansion in Legendre polynomials P l (x) is made [5,55], extending the expansion up to l = 30. The final length for these calculations is set to 70Å, which allows one to deal appropriately with the long-range oscillations of h(r) about zero (for more details see [2,33]).
BHP Fourier space triplet calculations at SP4 minimize with respect to t(r) the functional [t(r)] given in Equation (12c). The initial t 0 (r) is taken as t 0 (r) = h(r), and the integration range of distances R max is set to: (a) 70Å, and (b) 100Å. Two r-distance discretizations are studied: ∆r = 0.01Å and 0.005Å, which in defining t(r) imply 7001 points in 0 ≤ r/Å ≤ 70 (with 0.01), or 20,001 points in 0 ≤ r/Å ≤ 100 (with 0.005). Concomitantly, Fourier t(k)-values are treated in the same way by taking in each case equivalent discretizations (e.g., 20,001 points in 0 ≤ k/Å −1 ≤ 100, with ∆k = 0.005 Å −1 ). The numerical method chosen is a combination of conjugate gradient and pure gradient descents [5,10,56], as explained in detail elsewhere [34,35]. Such combination drives in general the minimization further down when the conjugate gradients "run out of steam" [56]. By doing so, a (double) is obtained. It is worth remarking that a usual criterion for convergence in this context [10] is defined in terms of a ratio between [t(r)] and a reference density derivative quantity, by requiring [t(r)] ≤ Once convergence is reached at a given step τ C , the final calculations of c ET (k 1 , k 2 ) [5] can be performed, thereby giving S ET (k 1 , k 2 ) as indicated in Equation (11). As regards the isothermal derivative of the direct correlation function, , the results obtained in [38] give two possibilities for carrying out the numerical treatment, since five states at 4.2 K were OZ2-studied with BDH+BHw: ρ N Å −3 = 0.02286713 ± ∆ρ N ; ∆ρ N = 0, 0.002, 0.004. There is the simple derivative estimate involving the two states adjacent to SP4 ("finer" Stirling), and also the more accurate estimate obtainable with Richardson's extrapolation that involves the four states around SP4 [57]. To visualize the situation, by denoting x = 0.004Å −3 , Stirling estimate is accurate up to terms of order O x 2 2 , while Richardson's extrapolation is accurate up to terms O x 4 . Consequently, for the sake of comparison, these two algorithms are employed in this work. Now, the significant ranges R Z for nonzero pair direct correlation functions [53,54] arising from the whole OZ2 treatment must be considered. As these OZ2 computations show, the significant R Z -zeros of the different c ET2 (r; ρ N ; R Z ) do not coincide with one another, and to calculate the density derivatives the c ET2 (r)-data regions needed are padded with zero-values. This extension with zero-values is also applied to every c ET2 (r; ρ N ) beyond its selected R Z up to the limit fixed for the variational calculation of t(r) R max = 70Å; 100Å . The latter action is consistent with the initial choice t 0 (r) = h(r), thus allowing for a long-range nonzero behavior of t(r) to develop. The longest {R Z − zeros} applications may be expected to perform better, in response to the wider radial c ET2 -behavior that they contain. However, the effect of this operation deserves closer inspection. Therefore, separate computations based on the Richardson extrapolation are carried out with the two sets {c ET2 (r; ρ N ; R Z )} The current minimizations have square norms of the derivative of c ET2 (r; ρ N ) at state point SP4 that are ≈ 5 × 10 6 . Note that: (a) rapid convergences are achieved, e.g., τ C =500-700 iterations; and (b) final [t(r)] values are ≈ 10 −2 − 10 −3 , which make ≈ 10 −9 − 10 −10 . At the final stages of the different minimizations, the convergences in the auxiliary function t(r) yield typically t τ (r) − t τ+1 (r) 2 10 −9 . Fourier transforms are performed via Fourier sums over the discretizations mentioned above. The behaviors of the auxiliary functions t(r) and t(k) are consistent with significant applications of the Fourier transform: the two functions tend effectively to zero as both r and k increase. In fact, by focusing the attention on the basic quantity t(r) it is worth remarking that the onset of its quick decay shows up in the region defined by the zeros {R Z } employed. To illustrate this point, some representative results, once convergence is reached, are quoted. Thus, for the case R max = 70Å; ∆r = 0.01Å Thus, there appears a t(r)-nonzero tail for distances greater than the longest maximal-R Z that defines the density-derivative nonzero range (r < R Z ), albeit the related features are rather small. The further enlargement to R max = 100Å only brings about very slight changes in the t(r) absolute values, thereby producing final S ET (k 1 , k 2 ) results in close agreement with the former at R max = 70Å (see Section 4). Among these changes in t(r), two may be worth mentioning: (i) within the essential region of nonzero density derivatives the values remain stable up to four/five decimal places; and (ii) for both R max = 70Å and 100Å and, roughly speaking, for distances in intervals within the range 17.6 r/Å 22.3, one finds a relative increase in the |t(r)| values with respect to their decay trends |t(r)| ∼ 10 −4 . As stated above, these changes do not alter the physics of the results obtained, although the second item might be subjected to a closer numerical inspection in future work. The foregoing general features are maintained under the different integration conditions used for [t(r)]. Figure 1 contains the final results at the pair level for state points SP1, SP2, and SP3. One observes in Figure 1a the rightward shift and smoothing of the instantaneous radial structure g ET2 (r) when the temperature is lowered (SP3-8.99 K vs. SP1-4.21 K) and the density is held fixed (0.022872 Å −3 ). Additionally, in Figure 1a the "compression" (or closer packing) of the instantaneous radial structure is apparent when the density is increased at constant temperature (4.21 K: SP1-0.022872Å −3 vs. SP2-0.027299Å −3 ). These are the expected features of g ET2 (r). In stark contrast with the results in the classical domain, the smoothing is consistent with the increasing quantum delocalization of the particles under diminishing temperatures. For visualization purposes, this general smoothing of structures is perhaps easier to grasp via the semiclassical Feynman-Hibbs Gaussian picture [13]. In this semiclassical approximation, each quantum particle is described by a thermal Gaussian packet. The pair radial function arises from a convolution, involving two of these packets, which smears out the sharper features of the pair radial function between the centers of the packets. With diminishing temperatures the width of the Gaussian packet increases, and through the convolution so does the smoothing effect (see [11] for further references and some related calculations). foregoing values are estimates obtained: (1) via three-point quadratic interpolation experimental (pressure, volume) data; (2) by selecting one of the representative str factors arising from the BDH+BHw iterative procedure (five BHw iterations, and th estcase) [38,53,54]. Although other estimates are possible, for the current pur their overall influence on the r-and k-results is slight. In relation to this, the av over -results for the ET2-computed isothermal compressibilities are worth qu / 0.005491 (SP1), 0.002867 (SP2), 0.005090 (SP3). Accordingly, no avera structure functions over the significant -ranges [53,54] are carried out in this work calculating ;

r-Space Triplet Correlation Functions
. Even though OZ2 applications to the instantaneous a proximations [11], not only these computations yield -values close to the experim values, but also keep their correct ordering. Given the high sensitivity of this ther namic quantity, the current structures obtained at the pair level are to be regarde very good representation of the actual helium-3 structures under the conditions st Their use as data input to triplet closure calculations is then fully reliable. Further su for the previous statement is provided by the PIMC calculation of pressures, which For the sake of comparison, it is worthwhile to mention that, after rounding three decimal places, the salient features of at SP4(4.2 K; 0.022867 Å tained in [38] coincide with those shown in Figure 1a for SP1 (see the Supplementa terial). Typical instantaneous estimates of the isothermal compressibility at SP4 / 0.0062 as ET2-computed [38], with the experimental estimate(s) be ET k = 0, or equivalently the corresponding isothermal compressibilities, which compare very well to the experimental values [48]. Thus, one finds: SP1 (χ T /bar −1 = 0.005715 (experimental), 0.005202 (computed)); SP2 (χ T /bar −1 = 0.002383 (experimental), 0.002844 (computed)); and SP3 (χ T /bar −1 = 0.005003 (experimental), 0.005071 (computed)). The foregoing values are estimates obtained: (1) via three-point quadratic interpolation of the experimental (pressure, volume)−data; (2) by selecting one of the representative structure factors arising from the BDH+BHw iterative procedure (five BHw iterations, and the largest-R Z case) [38,53,54]. Although other estimates are possible, for the current purposes, their overall influence on the r-and k-results is slight. In relation to this, the averages over R Z -results for the ET2-computed isothermal compressibilities are worth quoting: χ T /bar −1 = 0.005491 (SP1), 0.002867 (SP2), 0.005090 (SP3). Accordingly, no averages of structure functions over the significant R Z -ranges [53,54] are carried out in this work when calculating . Even though OZ2 applications to the instantaneous g ET2 (r) are approximations [11], not only these computations yield χ T -values close to the experimental values, but also keep their correct ordering. Given the high sensitivity of this thermodynamic quantity, the current structures obtained at the pair level are to be regarded as a very good representation of the actual helium-3 structures under the conditions studied. Their use as data input to triplet closure calculations is then fully reliable. Further support for the previous statement is provided by the PIMC calculation of pressures, which yields: SP1( For the sake of comparison, it is worthwhile to mention that, after rounding-off to three decimal places, the salient features of g ET2 (r) at SP4 (4.2 K; 0.022867Å −3 ) as obtained in [38] coincide with those shown in Figure 1a for SP1 (see the Supplementary Material). Typical instantaneous estimates of the isothermal compressibility at SP4 were χ T /bar −1 ≈ 0.0062 as ET2-computed [38], with the experimental estimate(s) being ≈ 0.00575 bar −1 [48]. (As regards the pressure at SP4, PIMC gave p bar =32.11 [38]). To grasp the SP1-SP4 differences in the ET2 computed χ T , one must remark the combined effect of the substantially larger N S (1372 vs. 1024) plus the more extensive statistics of the pair structures in the present work than what could be achieved in [38]. The current sampling yields highly accurate SP1 structures, hence it is the very low-k region of S (2) ET (e.g., χ T -computation at k = 0) which benefits greatly from this, whereas the rest of the structure factor remains practically invariable. The previous discussion lends additional quantitative support to the SP1-SP4 great proximity expected for their triplet structures. Figures 2 and 3 display the triplet instantaneous correlations computed with PIMC and closures and, for definiteness, Table 1 contains the PIMC salient features. As regards the PIMC results, the equilateral correlations g ET3 (r, r, r) follow closely the pattern set by the pair structures g ET2 (r): (1) although the triplet salient features are more pronounced, the positions of their peaks and valleys at a given state point are close to those at the pair level (Figures 2a and 3a,c); and (2) their relative heights are in correspondence with those at the pair level. In addition, their expected asymptotic behavior for long distances is that they tend to unity [33]. The isosceles correlations g ET3 (r, s, s) are plotted at the r-slice for the equilateral main peak positions r M (Figures 2b and 3b,d). A relevant isosceles trait is that for long s-distances these correlations do not tend to unity in general, but to the limiting pair value g ET2 (r M ) [33]. Owing to the range of distances L 4 scanned in the simulations, the foregoing asymptotic features are only hinted at by the oscillations in the related graphs. Nevertheless, PIMC runs at very low densities and the current closure applications (see below) indicate that the calculations are correct and would show these behaviors if allowances for much longer L and computational time were made. As regards the closure applications, three main points that agree with previous experience are to be commented [6,36]. First, the JF3 (AV3) bad behavior for short-range distances. Second, the remarkably good fitting provided by the closures when considered globally. Third, the closure departure shown in Figure 2c from the exact PIMC behavior. This departure is apparent as one inspects longer r-distances in the isosceles correlations measured with g ET3 (r, s M , s M ), where s M is a distance in the close vicinity of the main peak of g ET3 (r, r, r), i.e., Figure 2c shows the heights in the proximity of the main peaks belonging to increasing-r slices. One also notices that the use of the enlarged pair structures, incorporating grandcanonical corrections for the calculations of the closure values, is consistent with the PIMC canonical triplet structures.
Once again, a combined use of the closures KS3, JF3, and AV3 yields surprising results [2,6,36]. One can observe the excellent overall fit they give of the PIMC results for distances 2.5Å. The short-range distance behavior is correctly described by KS3, although JF3 clearly fails. Closure AV3 reduces the magnitude of such JF3-failure, but it is not sufficient to fix it. In the cases analyzed here, a combination of KS3 for the shortrange distances, say r 0 2.75Å, and AV3 for the distances beyond yields a very good representation of the triplet correlations in r-space (recall the in-built correctness of the asymptotic behavior of closures). The question of amending the closure failures remains intriguing and, because of the current applications to helium-3, may be more focused now. The instantaneous results in r-space obtained so far point decidedly to a prominent role played by the underlying pair correlations in shaping the triplet correlations when strong quantum diffraction effects are present. In relation to the quantum hierarchy of (instantaneous) correlation functions, we owe this theoretical pair-based picture to closures. Table 1. Detail of the salient features of the PIMC triplet calculations in the canonical ensemble. The reported positions and correlation values are in the close vicinities of: fp = first peak, fv = first valley, sp = second peak, sv = second valley. Distance r M is in the vicinity of the position of the equilateral first peak. For notational convenience: r = r 12 , s = r 13 , u = r 23 . Statistical error bars in the g ET3 values are ≤ 1%.      Figure 4a displays the instantaneous direct correlation functions c ET2 (r) associated with the instantaneous g ET2 (r) at the four state points about the target state point selected for the triplet calculations in k-space, SP4 4.2 K; 0.022867Å −3 [38]. It is interesting to note the fast decay to zero of these functions. Figure 4b contains representative results for the closure function t(r), Equation (12b), as arising from the variational Barrat-Hansen-Pastore (BHP) R max = 70Å; ∆r = 0.01Å calculations. One observes that, apart from the deeper bowl within short range, t(r) maintains its oscillations somewhat close to those of the starting choice t 0 (r) = h ET2 (r) = g ET2 (r) − 1. The effect of the derivative computation ∂c 2 (r 12 ;ρ N ) ∂ρ N on the final t(r) seems unimportant on the scale of the graph, as seen in the curves of the applications of Richardson's extrapolation (RM: four-point derivative) and Stirling (SM: two-closest-point derivative). Other details such as the effects of varying the ranges {R Z } considered or using a larger R max with a finer discretization, e.g., R max = 100Å; ∆r = 0.005Å , will be addressed below. In interpreting the k-space results, one must keep in mind that , ,

k-Space Triplet Structure Factors
sentially the Fourier transform of the involved correlation function , , which , , forms a distinctive part [11]. To draw significant conc from this line of thought, one would need an exhaustive -knowledge that g yond the equilateral/isosceles computations [6,11]. Consequently, Figure 5 conta results based on the alternative given by direct correlation functions. The panel (a) the equilateral instantaneous components , , /3 obtained with the four cl KS3, JF3, AV3, and BHP. One observes good agreement among the different appr past the main peak ≳ 2.1 Å ), which is related to the decreasing influence (Equation (11) In addition to the equilateral features, Figure 5b gives the isosceles compon the close vicinity 2.1 Å of the maximum of the equilateral components. D ancies between JF3 and AV3 are small (KS3 details can be inferred from the graph results follow very closely the JF3/AV3 trend within 1.32 ≲ , but they depart ably from it for 0 ≲ 1.32. In interpreting the k-space results, one must keep in mind that S ET (k 1 , k 2 , γ) is essentially the Fourier transform of the involved correlation function H ET3 (r 12 , r 13 , r 23 ), of which g ET3 (r 12 , r 13 , r 23 ) forms a distinctive part [11]. To draw significant conclusions from this line of thought, one would need an exhaustive g ET3 -knowledge that goes beyond the equilateral/isosceles computations [6,11]. Consequently, Figure 5 contains the results based on the alternative given by direct correlation functions. The panel (a) shows the equilateral instantaneous components S ET 2.7, 2.7, π 3 = 0.9814, S ET 3.1, 3.1, π 3 = 0.8744. tion with this, as was to be expected, relatively appreciable effects occur for low-k wavenumbers, e.g., the double-zero momentum transfer for which 70-SM yields 0,0 0.0275. The discrepancies between 70-SM and 70-RM can affect the third decimal places as observed in the equilateral components, but the relative importance diminishes with increasing k-wavenumbers due to the BHP-properties of , , (e.g., it monotonically increases for 0.9 ≲ /Å ≲ 2.1, remaining above 0.874 for larger wavenumbers).  In addition to the equilateral features, Figure 5b gives the isosceles components in the close vicinity k M = 2.1Å −1 of the maximum of the equilateral components. Discrepancies between JF3 and AV3 are small (KS3 details can be inferred from the graph). BHP results follow very closely the JF3/AV3 trend within 1.32 γ ≤ π, but they depart noticeably from it for 0 ≤ γ 1.32.   Figure 5c relates to the little overall importance of employing Stirling (70-SM) or Richardson (70-RM) derivatives in the current case, a fact that could change dramatically under conditions involving higher densities or more structured {c 2 (r; R Z )} functions. In connection with this, as was to be expected, relatively appreciable effects occur for low-k wavenumbers, e.g., the double-zero momentum transfer for which 70-SM yields S ET (0, 0) = −0.0275. The discrepancies be-tween 70-SM and 70-RM can affect the third decimal places as observed in the equilateral components, but the relative importance diminishes with increasing k-wavenumbers due to the BHP-properties of S (3) ET k, k, π 3 (e.g., it monotonically increases for 0.9 k/Å 2.1, remaining above 0.874 for larger wavenumbers).
At this point, one must also comment the effect of increasing the variational range of distances R max for obtaining t(r), together with the decrease in the r-discretization interval.
In this regard, BHP calculations of S  ET (k 1 , k 2 , γ) precludes a discussion of the relative merits of the closures employed in this work. However, some points allow one to make an educated guess. These points are: (a) the exactness of Equation (12d); (b) the diminishing influence of the triplet direct correlation function as wavenumbers k increase; and (c) the previous experience gained when studying liquid para-hydrogen [36]. On these grounds, it can be assumed that BHP gives a more acceptable description than the other three closures. Besides, pilot PIMC results at SP4 (2400 kpasses; N S × P = 128 × 66) indicate that [58]: (d) an equilateral negative-region is expected to be a possible genuine feature of the triplet instantaneous structure factor (compare with [36]); and (e) there is a better behavior of BHP in such region, as deduced from the following partial median S ET k, k, π 3estimates obtained with PIMC: (i) at k = 0.5Å −1 , median ≈ −0.007, and (ii) at k = 1Å −1 , median ≈ −0.045. The latter magnitudes are not final and differ quantitatively from BHP's, but they signal the negative trend. Accordingly, the pronounced negative region found in KS3 (and AV3) for the equilateral components is likely unphysical, which should be consistent with the discrepancies between PIMC and KS3 for the salient r-isosceles features shown (for SP1 ≈ SP4) in Figure 2c. More relevant quantities, such as the actual intensities of the S ET main equilateral peak and its associated isosceles features, are also expected to be better captured by BHP, although only final PIMC results will settle this question. In relation to this, note that the individual S ET -values fluctuate strongly when PIMC-sampling the configurations [36]. Therefore, the computation of precise average values entails very long run lengths in this type of calculations.

Conclusions
The current computational structure study has dealt with supercritical helium-3. Real space and Fourier space properties of the triplet instantaneous correlations have been investigated.
As regards the real space results, the exact PIMC equilateral and isosceles features show the influence of independent variations in temperature and in density. Thus, the structure is smoothed by the decrease in temperature, whereas it is sharpened by the increase in density. When inspecting the equilateral correlations, one observes that the salient features, i.e., positions and heights of the peaks and valleys, follow the patterns set by the pair correlations, albeit the equilateral triplets show far more pronounced first peaks and valleys. These traits are in accordance with what was obtained for the hard-sphere fluid [2] and liquid para-hydrogen [6].
From the comparison with PIMC one finds, once again, that the triplet closures used (KS3, JF3, AV3) reveal themselves as a great help in providing real-space physical pictures of triplet correlations in fluids with quantum behavior. This closure usefulness is thus consistent with that found in other applications [2,6]. One might expect this positive working of closures to occur when studying fluid helium (far from quantum exchange), because of its similarity to a quantum hard-sphere fluid [2,16]. Speculations are fine when they are consistent with known related facts, but there is nothing like direct proof. The current applications to helium-3 fill the conditions to be analyzed in that a continuous interparticle potential (as complementary to the hard-sphere potential) plus very strong quantum diffraction effects (as an extension to para-hydrogen conditions) are dealt with at the same time. The main conclusion is that pair correlations contribute decisively to shape the correlations at the triplet level. A combination of KS3 for small ranges of distances plus the use of the intermediate AV3 = (KS3 + JF3)/2 beyond these ranges yields a very good representation of the exact correlations. On the one hand, this (KS3+AV3)-representation is surprisingly accurate for: (1) the equilateral correlations over long ranges of distances, and (2) the isosceles correlations over short-medium ranges of distances. On the other hand, (KS3+AV3) loses predictive power for isosceles correlations beyond the ranges mentioned. Even though ranges of distances more general than those considered herein remain to be computed, the variety of quantum systems and conditions studied so far indicates that this usefulness of the closures utilized is not a fortunate coincidence. Therefore, it should be regarded as a general fact for fluids with strong quantum behavior. Although it is not fully clear how to continue expanding Equation (9), because of the convergence properties, this goal seems now to deserve a try.
The triplet instantaneous structure factor computations with closures at (T = 4.2 K, ρ N = 0.022867Å −3 ) show the important role played by the triplet direct correlation function. Because of Equation (11), beyond the region of the main peak k 2.1Å −1 , the smallness of the latter function makes KS3, JF3, AV3, and BHP be in close agreement with one another on the equilateral component values, and also on the isosceles components for large angles γ 75 • . At γ = 60 • BHP separates from the rest at the main peak, and far more noticeably within 0 ≤ k/Å −1 < 1.5. The same sort of discrepancies are found when inspecting the isosceles components, e.g., at k = 2.1Å −1 , for angles γ 75 • . By construction JF3 takes positive values all over the possible ranges of wavevectors. However, KS3 (and AV3) and BHP take negative equilateral values for low-k wavenumbers. BHP (absolute) negative values are small, but this does not happen to KS3 which near 1Å −1 , displays a pronounced dip below zero. Given the KS3 defect in reproducing the spatial triplet isosceles correlations as the ranges of distances are enlarged, its overall k-space behavior cannot be regarded as correct. The same may be said of AV3. BHP seems, however, better adapted to the task of giving better estimates of triplet structure factors. A reason for this is Equation (12b), an exact relationship which is a BHP's in-built feature. The negative-valued region for the equilateral components obtained in the BHP calculations finds qualitative support in pilot PIMC results. Accordingly, one is led to surmise that such negative region should be a genuine triplet fact. Furthermore, the calculation of estimates of the double-zero momentum transfer component, which is a response property of a quantum fluid independent of the structural description [6], can be better carried out with the centroid correlations due to the exact OZ3 framework they provide. There is ongoing work, involving PIMC and closures, which is focused mainly on the delicate problem of determining facts as accurately as possible about triplets in k-space for fluids with quantum behavior. The results, covering the instantaneous and the centroid structures, will be the subject of a future article.
Funding: This research received no external funding.