Properties of fermionic systems with the Path-integral ground state method

We investigate strongly correlated many-body systems composed of bosons and fermions with a fully quantum treatment using the path-integral ground state method, PIGS. To account for the Fermi-Dirac statistics, we implement the fixed-node approximation into PIGS, which we then call FN-PIGS. In great detail, we discuss the pair density matrices we use to construct the full density operator in coordinate representation, a vital ingredient of the method. We consider the harmonic oscillator as a proof-of-concept and, as a platform representing quantum many-body systems, we explore helium atoms. Pure $^4$He systems demonstrate most of the features of the method. Complementarily, for pure $^3$He, the fixed-node approximation resolves the ubiquitous sign problem stemming from anti-symmetric wave functions. Finally, we investigate $^3$He-$^4$He mixtures, demonstrating the method's robustness. One of the main features of FN-PIGS is its ability to estimate any property at temperature $T=0$ without any additional bias apart from the FN approximation; biases from long simulations are also excluded. In particular, we calculate the correlation function of pairs of equal and opposite spins and precise values of the $^3$He kinetic energy in the mixture.


Introduction
Strongly interacting fermions are a building block of matter in scales ranging from nuclear physics to neutron stars [1][2][3][4]. Over the past several decades, efforts to unravel theories to describe these systems have raised several issues that attest to how challenging it can be to approach them comprehensively [5].
Mixtures of bosons and fermions are an active field of research with scientific interest renewed by the recent observation of quantum gas degeneracies [6][7][8]. A prominent example is helium systems, which have, for a long time now, been the bedrock of strongly interacting many-body physics. Despite that, they exhibit properties and features that are not entirely understood yet, such as the kinetic energy of isotopic 3 He-4 He mixtures [9,10]. One of the goals of this work is to investigate the kinetic energy of mixtures of 3 He in the normal phase and 4 He in the superfluid one, with both species receiving a fully quantum many-body treatment concerning the statistics they obey. For this purpose, we have revisited a quantum Monte Carlo (QMC) method at zero temperature [11], demonstrating its accuracy. It allows the quantum many-body treatment to be subjected exclusively to an approximation of the nodal structure of the wave function for the projection is long enough in imaginary time. In this work, we review in depth the PIGS method, including its extension to treat fermions using the fixed-node approximation. With a different focus, an excellent review of the method was put forward by Yan and Blume [37].
PIGS estimates any property without biases, regardless of whether or not the property commutes with the system's Hamiltonian. In our approach for fermionic systems, all properties are estimated within the fixed-node approximation, where only a single given nodal region is considered, avoiding changes of sign in the sampled probability. The fixed-node approximation is straightforward to implement within PIGS. FN-PIGS is a state-of-the-art method for treating fermions, paving the way for further studies of local, non-local, and dynamical properties free from extrapolations.
In the following sections, we present details of the PIGS method and results for a proof-ofconcept system. We then show results for pure 4 He systems before reaching the primary purpose of our work, which is the investigation of fermionic systems. We investigate systems composed of 3 He atoms and the much more complex case of the 3 He- 4 He mixture.
PIGS and PIMC methods are similar. Many essential aspects of the method can be found in the magnificent review by David Ceperley about path integrals [19]. Here, we give technical details and particularities of the numerical treatment involved in PIGS. We also provide a detailed description of how estimators are constructed, particularly for the radial distribution function.
Finally, we mention that we have limited the study of 3 He- 4 He mixtures to some of their basic properties, leaving aside other important and interesting characteristics of these quantum liquids. For instance, we have yet to attempt to apply the worm algorithm to evaluate the one-body density matrix and the condensate fraction of the 4 He isotope. With our work as a stepping stone, we expect these properties to be investigated in the future.

The PIGS method
The formal solution of Schrödinger's equation is the time evolution operatorK(t) = e −i tĤ , where t is the time andĤ is the Hamiltonian of the quantum system. This operator is intrinsically connected to the density operator in statistical mechanics,ρ(τ) = e −τĤ , where τ = 1/k B T is the inverse temperature, k B is the Boltzmann constant, and T is the temperature. One can transform ρ intoK by applying the so-called Wick rotation τ → i t. Therefore, τ is also commonly referred to as the imaginary time. This transformation leads to several approaches that explore a quantumto-classical mapping to calculate statistical properties and expectation values in quantum systems without dealing with integration measures that are not positive-definite or even real, such as the renowned PIMC method [19].
In the PIGS method the density matrixρ(τ) = exp(−τĤ) projects a trial state |Ψ T 〉 onto a quantum state |Ψ τ 〉 where τ is a real number. We assume that all quantum states are not normalized. If |Ψ T 〉 is nonorthogonal to the ground state |Ψ 0 〉, the projection is guaranteed to exponentially approach |Ψ 0 〉 as τ increases. This can be easily verified expanding |Ψ T 〉 into eigenstates of the Hamiltonian, where E i < E i+1 are eigenvalues ofĤ and c i are expansion coefficients. For a large enough τ, the term with the smallest decay rate, namely the ground state, entirely dominates the expansion.
In practice, to carry on computer simulations, these projections are performed in coordinate space representation, where the density operator has the matrix elements where R ≡ {r i | i = 1, . . . N } is the set of coordinates r i of all N particles in the system, which can be of any spatial dimension. The PIGS method takes advantage of the convolution property of the density matrix, also called imaginary time composition, which reads and which becomes an explicit convolution operation in coordinate space: The last equation is an exact expression of the density matrix for any imaginary time τ. Apart from normalization constants, the ground state expectation value of any quantity O associated with the operatorÔ, can be estimated in the following way. Consider initially O(τ) given by We can insert completeness relations1 = dR i |R i 〉〈R i | to the left side of |Ψ T 〉 and to the right side of 〈Ψ T | such that with Ψ * T (R −M ) = 〈Ψ T |R −M 〉 and Ψ T (R M +1 ) = 〈R M +1 |Ψ T 〉. Next, we use imaginary time composition to slice the density operator on the right side ofÔ into M operators e −δτH , with δτ = τ/M , such that Insert another completeness relation to the left of the right-most short-imaginary time operator, where we can see the formation of a "chain" of density matrix elements ρ(R i , R i+1 ; δτ) with endpoints R 1 and R M +1 , that are then integrated over the entire space. The same procedure can be repeated to slice e −τH on the left side ofÔ into M e −δτH operators and then insert intermediate This equation has a visual representation that is shown in Fig. 1. Each of the 2M + 2 configurations R i , represented as a bead, is linked to R i+1 through the matrix element ρ(R i , R i+1 ; δτ), which is represented by a spring. The exception is the central configurations R 0 and R 1 , connected by the operatorÔ, which is represented by a horizontal bar. The endpoint configurations R −M and R M +1 carry the trial wave function Ψ T .

Non-local operators
The probabilistic nature of ρ and Ψ T offers the possibility of using Monte Carlo techniques to sample the configurations R i , with i = 0, 1. To sample the central configurations R 0 and R 1 , we are left with a choice to use any suitable probability distribution p(R 0 , R 1 ), since Eq. (13) can also be written as where such that R 0 and R 1 are amenable to importance sampling. The most convenient choice is to take p(R 0 , R 1 ) = ρ(R 0 , R 1 ; δτ), such that R 0 and R 1 are also connected by a spring, and the same algorithm can be used to sample the configurations of all the beads R −M , ..., R M +1 in the chain. Despite its convenience, this may not be the best choice depending on the nature ofÔ, and a more appropriate p can be selected case-by-case. Finally, O(τ) can be estimated by sampling the probability distribution

Local operators
A particular relevant case of the PIGS method happens when the operatorÔ is local in coordinate representation, such that: forcing the configurations R 1 and R 0 to be the same. The integration in one of these coordinates can therefore be analytically performed, effectively collapsing the two central beads into a single one. This results in a chain of 2M + 1 beads, all connected by density operators ρ(R i , R i+1 ; δτ), as shown in Fig. 2.
In this particular case, one does not need to worry about choosing p, and estimates for O(τ) can be obtained by sampling the probability distribution Ground state converged configurations must be used when evaluating the local quantities O(R 0 , R 1 ).
In particular, it is imperative that the configurations (R 0 , R 1 ), points at whichÔ is effectively evaluated, satisfy this ground state convergence condition. Another interesting point here, which contrasts this method to others that share the same imaginary time composition approach enclosed in Eq. (7) through Eq. (13), is that, despite having a total of (2M + 2) beads (or intermediary configuration), we are slicing the imaginary time component into M chunks. Configurations from Figure 2: Visual representation of the PIGS method for the particular case whenÔ is local in configuration space. The total number of beads is 2M + 1, and O is calculated in the central bead R 0 . Density matrix elements connect all internal beads. Therefore, the same algorithm can be used to sample all configurations R i .
Ψ T (endpoints) can be sampled via the Metropolis-Hastings algorithm. Those of the remaining beads can be efficiently sampled by a multi-level Metropolis algorithm (see Sec. 5).
In summary, the PIGS method estimates ground state properties without introducing any biases. This calculation is done through the application of a local operatorÔ in converged configurations of the ground state, which are obtained by projecting a trial quantum state represented by the wave function Ψ T at both extremities of a linear polymer using the density matrix of the system. Estimates can be refined by adding more beads and/or applying the estimator O to all beads with converged configurations when [Ô,ρ] = 0.
To fully account for quantum fluctuations, the probability distribution P must include Bose or Fermi statistics, whereP is the operator that permutes P particles of the system. In the PIGS method, the appropriate quantum statistics of the system can be implemented through the trial state Ψ T . In this case, then and all terms in the sum of Eq. (19) are the same under the integration after relabelling particle indexes. Thus, it is possible to consider the density operator elements ρ(R, R ; δτ) for distinguishable particles [11,19]. For fermionic systems, Ψ T must be anti-symmetric under particle exchange. For this reason, P may change sign, which is the source of the sign problem in the PIGS method. The fixed-node approximation, which will be discussed in the next section, can be used to circumvent this problem.

Fixed-node approximation
The fixed-node approximation is a common way to address the sign problem in QMC methods. It was envisioned by Anderson [38], and its first application to the investigation of the electron gas ground state attained results that became an essential ingredient for the density functional theory [39,40]. In a somewhat more intricate fashion, similar ideas were employed to restrict path integrals from sign changes [28]. In the PIGS method, the fixed-node approximation can be implemented straightforwardly.
In fermionic systems, the nodes of the ground state wave function Ψ 0 carve up the configuration space into nodal regions. Ceperley demonstrated that two distinct nodal regions of the ground state differ only by the permutations of particle indexes, i.e. by the sign of Ψ 0 [28]. Consequently, properties associated with a given nodal region Ω are identical to those of other nodal regions. Thus, one can solve Schrödinger's equation inside one of the Ω regions and still obtain expectation values pertinent to the entire system.
In practice, the nodal hypersurfaces of a given ansatz Ψ T for the ground state are converted into infinite potential walls, such that all configurations R i , with i = {−M , . . . , 0, 1, . . . , M + 1}, are forced to remain inside the initial, arbitrary nodal region. In general, the nodal regions of Ψ 0 are unknown and are likely different than those of Ψ T . Nevertheless, it has been shown that, by using a single Ω, one obtains the lowest energy estimate of the system, consistent with the nodal structure of the trial state [23,28].
In a typical simulation, Monte Carlo moves are rejected if a configuration of the trial wave function or that of any bead crosses the nodal hyper-surface. A test to determine if the surface has been crossed can be established by verifying if the sign of Ψ T (R i ) has changed after the trial displacement. Cross-recross errors happen when a configuration R i leaves a nodal region with a particular sign and ends in another region with the same sign as the initial one. To avoid this error, we use sufficiently small values of the imaginary time step δτ, which results in acceptance ratios greater than 90%.
To account for the infinite potential walls at the nodal hyper-surface, the density matrix ρ(R, R ; δτ) must smoothly vanish at these boundaries. This is accomplished by implementing the image action approximation [41]. For instance, in a short-time δτ approximation, one can multiply ρ(R, R ; δτ) by e S I , where λ = ħ h 2 /2m, and d(R) is the distance between configuration R and the closest nodal region. The use of this approximation is strongly recommended. Without that, there would be a non-vanishing probability of finding beads in the proximity of the nodes. An exact analytical expression for d(R) is difficult to obtain, but one can rely on the Newton-Raphson method, which gives the approximation

The density matrix
A central ingredient of the PIGS method is the representation of the density operator of the system e −τĤ in coordinate space or, in other words, the elements ρ(R, R ; τ) of the density matrix [19,42].
A path integral representation of this object can be obtained by formally bringing the number of imaginary time compositions in Eq. (5) to infinity [43], which yields where the integration measure is and is called the total action of path X (t), with X i = X (i × δτ) and δτ = τ/M . Notice that all paths start at R and end at R after an imaginary time τ, which is effectively a boundary condition on this functional integration. These paths are indeed continuous, driven random walks in the M → ∞, or conversely δτ → 0 limit. In general, the challenge of computing these matrix elements stems from the fact that the Hamiltonian is almost always a sum of two or more non-commuting terms, such that it is impossible to calculate the matrix elements from each contribution individually. For the typical case whereĤ is a sum of a kinetic term plus a potential term,Ĥ =K +V , this can be seen by noticing that with higher-order terms following the general recipe from the Baker-Campbell-Hausdorff formula [44][45][46]. However, this relation indicates that these contributions can be neglected in the limit τ → 0. In fact, the Trotter-Suzuki decomposition shows that, precisely in that limit, the contributions coming from the commutators vanish [47]. This is a rigorous mathematical result and makes complete sense physically. It represents the observed nature that, in the limit of very short τ, which translates into either short time intervals or high temperatures, the system behaves classically [48]. Since the Trotter-Suzuki limit is intrinsically enclosed in each element of the total action in Eq. (25), we can rigorously separate K and V under the path integral exponentiation. The kinetic term can be analytically solved, having a well-known Gaussian distribution form, such that the elements become where d is the number of spatial dimensions. The Gaussian envelope from the kinetic contribution implicitly sets a relevant length scale over which the random walks have finite, non-vanishing probabilities, centered around the point where the interaction potential is being calculated and which is controlled in size by δτ. These envelopes can be formally combined into the integration measure of Eq. (23), such that the density matrix elements are represented by a Wiener process through where D W [X (t)] is the Wiener measure [49], which attributes statistical weights to each freeparticle path starting at R and ending at R after an imaginary time τ. The Wiener process expression of the density matrix describes a Brownian random walk (BRW) driven by the interaction term [50]. The functional integration can be written as an average over these random walks, in what is known as the Feynman-Kac formula [51]: where ρ 0 is the density matrix element for the non-interacting system, .
Equation (30) defines a framework amenable to several approximations guided by physical intuition that will be discussed next. However, more common mathematical approaches can also be used, such as manipulating high-order commutators in the Baker-Campbell-Haussdorff expression [52][53][54], or general classical mechanics methods. One example is the use of van Vleck determinants to find trajectories with high stability against initial conditions [55,56].

Primitive approximation
The most common, widely used expression for the density matrix considers that τ is small enough that the Brownian motion is effectively captured almost entirely by its initial and final states at R and R . This leads to the following expression, which is formally correct up to O(τ 2 ), The symmetrized form of the potential action U(R, R , τ) in the primitive approximation is simply the interatomic potential. The primitive approximation is easy to compute during simulations, not requiring any operations related to the interaction potential and, therefore, being available to virtually every many-body system, apart from singular cases [57]. It performs exceptionally well for smooth, slowly varying potentials. The pitfall here is that representing the entire range of the Brownian random walks exclusively by the endpoints is a pretty drastic restriction that can only be reasonable for very short imaginary time steps τ. Consequently, it requires slicing the density operator into a sometimes prohibitively large number of chunks. This renders sampling the distribution of Eq. (16) a computationally demanding task in terms of time and resources. In practical terms, this means having to equilibrate, in a Markov process sense, many interacting polymers, each composed of a large number of beads. More accurate approximations are highly welcome for simulations to be efficient and, despite not being very common, can also be easily achieved.

Semiclassical approximation
In many cases, it is possible to obtain a substantially better result by considering a WKB approach, which consists of finding the path amongst the Brownian random walks that has the largest statistical weight, denoted X WKB (t), and evaluating the interaction functional in Eq. (30) over that path. Aside from being physically intuitive, it can be rigorously shown that the path with the largest weight is exactly the classical path connecting R to R , justifying the name of the approximation. Moreover, since the BRWs are paths of the non-interacting system, X WKB (t) is simply a straight line where The integration can be translated into a kinematics problem with a definite, analytical solution for some potentials. When that is not possible, direct numerical integration, with tabulation of the resulting terms, is a good alternative. This table can then be interpolated during simulations. In several situations, this is computationally faster than dealing with analytical forms, where intricate potentials must be evaluated repeatedly at every iteration.
By considering the interaction along the entire path connecting R to R , the semiclassical approximation can capture the physics of the system better than the primitive approach and should always be preferred whenever possible. However, since it is not derived directly from Baker-Campbell-Hausdorff terms, it is uncertain to what order of τ this approximation is correct, although we certainly know that it is at least O(τ 2 ). Despite that, given a certain τ, one can calculate the relative weight of X WKB to other possible random walks and thus have insights on how large τ should be to obtain the desired accuracy. These additional BRWs can be sampled with a simple Monte Carlo algorithm. The downside of this approximation is that, as we shall see, some estimators require derivatives of the density matrix that may not be easily computed from this expression.

Pair product approximation
In a variety of quantum systems, the many-body physics is encapsulated into a two-body potential v, such that the interaction term of the Hamiltonian is where (i, j) denotes a pair of particles. With that, we can write the many-body density matrix of Eq. (30) as The pair product approximation, ρ PP , consists of neglecting product correlations of orders higher than two in the Wiener process, such that and where SciPost Physics Core With that, the critical ingredient for the many-body density matrix is now the interaction term of one pair of particles, a much simpler object. This approximation has several advantages. First, it is exact for a pair of particles by construction. Second, for homogeneous systems, the neglected correlation effects are largely attenuated since they tend to be canceled out by the presence of other particles in other directions that will have opposite correlations. Third, this approximation is well-defined for all sorts of interactions. Moreover, finally, for large imaginary times, this expression will approach the solution of a twobody Schrödinger equation (weighted by the corresponding energy eigenvalue), resulting in a many-body density matrix equivalent to a Jastrow-type wave function. These wave functions are well-known to accurately capture most ground state short-range correlations in systems that exhibit collective phenomena [30,33,34,58,59]. For systems composed of helium atoms, as stated by Leggett [60], the Jastrow function ansatz is the archetypal form of a variational ground state wave function. Corrections due to long-range correlations can be made based on quantizing the classical sound field and considering the zero-point motion of longitudinal phonons [61,62]. Although this correction improves the accuracy of the static structure factor at small wave vectors, contributions to the energy are small because the interatomic potential is not long-range. Further improvements based on parametrical expansions are also possible [63]. In other systems, it may be very well the case where higher-order and long-range correlations are important, which must be addressed case-by-case.

Central potentials
When the interaction is isotropic, depending only upon the magnitude r of the relative coordinate of the pair of particles (i, j), r = |r i − r j |, it is possible to employ a center of mass transformation such that the pair density matrix elements become where r CM = (r i + r j )/2 and r CM = (r i + r j )/2 are the center of mass coordinates, and r = r i − r j and r = r i − r j the relative coordinates at the endpoints. The center of mass density matrix ρ CM is effectively a free-particle one-body term on these coordinates, with the resulting mass m CM = m i + m j . The relative density matrix ρ rel is the solution to the problem of a single particle of reduced mass m rel = m i m j /(m i + m j ) under the action of an external potential v(r).
In this scenario, one can immediately make use of one of the short imaginary time approximations to obtain an expression for ρ rel , then reconstruct the pair density matrix ρ pair , and finally obtain the many-body object with the pair product approximation. Alternatively, since this problem is much simpler than the many-body one, one can attempt to solve the one-body Schrödinger equation and construct an exact density matrix from where φ i are the eigenfunctions and E i the eigenvalues of the corresponding single-particle Hamiltonian.

Numerical convolutions
Accurate density matrices for imaginary times τ larger than the ones for which short-time approximations hold well are essential to reduce the total required number of beads in a simulation. They allow the PIGS method to project the system's ground state faster or with more particles with the same computational resources. Here, once again, the convolution property is quite convenient since we can use the approximations for short imaginary times and then self-compose them to obtain the expression for a larger imaginary time. This is not practical for the full many-body matrix because it requires numerical integration of a highly multi-dimensional object. However, the central potential case discussed in the last section offers a great advantage: the relative coordinates term can be expanded in partial waves. In three dimensions, ρ rel becomes where θ is the angle between r and r and P are Legendre polynomials. It can be shown that the partial waves ρ are solutions of the one-body problem of a particle under an external potential v(r) and a centrifugal barrier λ rel ( + 1)/r 2 , with λ rel = ħ h 2 /2m rel . Each partial wave ρ (r, r , τ) in Eq. (41) is the solution to the following Bloch equation: with boundary conditions ρ (r, r , 0) = δ(r, r ) and ρ (0, r , t) = 0. For the free particle case, v(r) = 0, the solution is the free-particle contribution Finally, the partial waves ρ for the three dimensional case can be written as where i is the modified spherical Bessel function, and CBRW stands for centrifugal Brownian random walks. In the short imaginary time limit, the CBRW converges to the usual BRWs, which can be shown by employing an asymptotic expansion of i l .
Since the partial waves are completely independent, each one of them satisfies its own convolution property, which is now a one-dimensional integration. One can then consider imaginary times short enough such that an approximation (primitive or semiclassical) is reasonable, and then compute the partial waves with Eq. (44) and perform numerical integrations, each one effectively doubling the initial imaginary time and yielding the exact density matrix elements apart from numerical errors [64].
The starting temperature can be relatively high, such as T ∼ 10 3 K for Helium systems. The integration is repeated until the desired temperature is reached. In the calculation of ρ rel using Eq. (41), the sum of partial wave components at the desired temperature is truncated when contributions to the sum are smaller than a chosen threshold ε, typically ε ∼ 10 −8 . As a bonus, the larger the imaginary time interval, the smaller the number of effectively contributing waves. For fixed r and r , contributions from partial waves with large angular momenta become increasingly irrelevant with increasing δτ, which is a direct consequence of the spherical Bessel function weights.
Attempts to sample configurations that enclose larger relative distances become less rare at low temperatures since the free particle Gaussian distribution mainly dictates sampling. This distribution broadens as δτ increases. Such subtlety is at the heart of implementing FN-PIGS and other path-integral methods. Even if one can find an incredibly accurate density matrix at low temperatures, one must still employ a substantial number of beads in the simulation. Using few beads results in large displacements in the bisection algorithm, which tend to be rejected by the repulsive interaction part of the density matrix (particles tend to fall too close to others). The ideal average displacement, and therefore the associated value of δτ and the number of beads, is primarily controlled by the density of the system.

Construction of the many-body action
In practice, evaluating the action for a pair of particles, viz.
is a costly task to be performed during simulation time, even when analytical expressions are available. A possible procedure is to store the u i j values in a table and form the many-body action via the pair product approximation, which yields and the full density matrix elements are ρ(R, R ; δτ) = ρ 0 (R, R ; τ)e −U(R,R ,δτ) . Even then, simply tabulating u i j as a function of distances r, r and the angle θ = arccos(r · r /r r ) for the required time intervals used in sampling P can generate very large four-dimensional arrays.
A much more efficient representation can be achieved by decomposing the pair action into diagonal and off-diagonal contributions [19], and writing u i j as where u ep (r, r , δτ) = [u i j (r, r, δτ) + u i j (r , r , δτ)]/2 is the so-called end-point term, and u od is the remaining contribution. The end-point term can be efficiently stored in tables since it is essentially a one-dimensional object.
The off-diagonal terms can be expanded by noticing that the density matrix elements are centered around the diagonal r = r , a consequence of the Gaussian envelope that we discussed earlier in the path integral formalism. A suitable change of variables that encloses the proximity to the diagonal is given by (r, r , θ ) → (q, s, z), with q = (r + r )/2, s = |r−r | and z = r − r . The variable q is the only one that is not restricted by the Gaussian term [19]. With that, u od can be written as the following polynomial representation: This expression can be truncated, and the coefficients c i j are obtained by fitting a surface to u od for a set of values of q that are relevant to the problem. The coefficients are then stored in one-dimensional tables that can be efficiently accessed during simulations, allowing for quick construction of the many-body action. For Helium systems, i = 2 is an adequate truncation point. Figure 3 shows the resulting coefficients for a pair of 3 He atoms interacting via the Aziz potential [65] at temperature T = 50 K.

Sampling algorithm
In the PIGS method, sampling the relevant probability distributions can be done with stochastic techniques. However, for polymers with more than 16 beads [19], detailed balance is challenging to achieve when one relays on the canonical Metropolis algorithm and its variants based on configuration-by-configuration sampling [66,67]. The multi-level Metropolis algorithm speeds up the sampling process and is, therefore, preferred [19]. This algorithm divides the sampling into stages (levels). In the initial stages, a crude but fast approximation for the probability density is used to decide if the proposed moves are accepted. In the final stage, the most accurate expression of the density matrix is employed, and its associated probability density is sampled for all beads.
The idea is to use more simple forms of the density matrix during the early stages to filter trial moves with a reasonable probability of being accepted in the final stage. It is undesirable to spend simulation time evaluating sophisticated expressions of the density matrix before the last stage. For example, one can use the primitive approximation to discard situations where hardcore particles overlap. The more accurate pair density matrix expressions can be used at the final stage.
The multi-level Metropolis algorithm implementation we employ is known as the bisection algorithm. In a bisection of level L, a segment R = {R i | i, . . . , i + N } with 2 L + 1 beads of the open necklace is randomly selected. If (i + N ) > M , the segment is discarded, and another one is selected. At the first stage, = 1, a trial configuration R m is proposed for the central bead R m of the segment delimited by (R i , R i+N ), N = 2 L . If the proposed move in this stage or any subsequent one is rejected, a different segment R is selected, and the iteration is restarted. When the trial configuration R m is accepted, the algorithm proceeds to the next stage. At stage = 2, the chosen segment is split in two with extremities (R i , R m ) and (R m , R i+N ). Each one includes the updated configuration R m of the previous stage. A new configuration for the central bead of each one of these two new segments is then proposed. Although the proposed moves are independently generated, the criterion for their acceptance, as we will see in Eq. (54), considers all of them simultaneously. If the proposed moves of the centers in this stage are accepted, the algorithm proceeds to the next stage. The same recipe is repeated recursively in the following stages, considering bisections of segments constructed in the prior stages. 2 −1 segments are to be considered at each stage. In Fig. 4, we show a graphical representation of the bisection algorithm with L = 2.
A segment with extremities (R e1 , R e2 ) has its central bead R m separated by an imaginary time δτ = 2 L− δτ from the endpoints. The probability density associated with this segment is with the action U given by a single sum over all the segments S = R e1 − R m − R e2 considered at the present stage, The product of the free-particle density matrices in Eq. (50) is proportional to and therefore can be exactly sampled with the trial displacement where η is a vector of random numbers generated from a normal distribution of unitary variance and zero mean. This trial configuration is accepted or rejected by the Metropolis algorithm according to the acceptance probability where U trial is obtained by taking R m = R trial m in U , with U 0 = U trial 0 = 0 by definition [67]. When calculating U trial at any stage, we consider the updated configurations of all previous stages. The original configurations of the segment are used to compute U . For < L, U is calculated from the primitive approximation.
The algorithm proceeds to the next stage if the trial move is accepted. Otherwise, the entire segment is left unchanged, a new one is selected, and the iteration is restarted. In the last stage, = L, each segment is composed of three beads. Except for the endpoints R i and R i+N and those at the center of all segments with three beads, all others have configurations from previously accepted moves. The imaginary time interval to be employed in this stage is the one used for calculating the pair density matrix, and U is calculated using the more refined pair product expression. If accepted, bead configurations are updated as {R i = R i | i + 1, . . . , i + N − 1}; otherwise, R remains unchanged. Since configurations of beads R i and R i+N are never changed, the procedure of choosing a segment R needs to be repeated several times for the same polymer.
In the final stage, the denominator in Eq. (54) cancels the approximated probability densities used in the previous stages. The final acceptance probability A = A is how the most accurate expression of the many-body density matrix for accepting moves of all bead configurations is obtained. For this reason, the movements at all levels must be either all accepted or all rejected. In this way, the detailed balance required for an unbiased estimate of properties is satisfied.
Generally, acceptance is controlled by the imaginary time interval δτ and the total number of stages L. For bosonic systems, the usual acceptance is ∼ 20%. However, a finite probability of having a cross-recross error exists for fermionic systems, where the fixed-node approximation is employed. This error occurs when a trial movement crosses two nodal surfaces, reaching pockets with the same wave function sign. To avoid this error, we monitor the number of attempted moves that result in a sign change of the trial wave function and are discarded due to the fixed node approximation. This number, for acceptance ratios above 90%, represents less than 0.4% of the total attempted moves in a simulation, and it is even smaller, 0.02%, for acceptance ratios over 99%. In this scenario, having a cross-recross error due to a large displacement is improbable. This unlikelihood is reflected in the fact that ground state energies obtained with acceptance ratios above 90% statistically agree among themselves. However, higher acceptance rations imply much longer simulation times. Thus, 90% is a recommended value.
A Monte Carlo sweep is finally completed with attempts to move the configurations that carry wave functions Ψ T at the extremities of the open necklace. Trial configurations are generated according to the Metropolis algorithm [66,67]. The canonical algorithm is associated with a step size that regulates the average displacement of the particles and the acceptance ratios. This step size is also adjusted to give an overall acceptance ratio of 90% when the fixed node approximation is used.
Over the course of decades, extensive variational studies provide a collection of suitable trial functions related to various physical systems. Nonetheless, the art of constructing these functions is still an active field, including machine learning techniques [68][69][70][71][72][73][74] and quantum computing [75]. It is essential to mention that accurate trial functions expedite simulations in the sense that shorter projections, resulting in fewer beads, will be necessary to reach convergence.
Most of our simulations were performed with an L = 3 bisection algorithm. In this case, the number of contiguous beads selected to be treated independently is 7, the two extremities of the segment are kept fixed. The total number of beads in each polymer usually considered for the density matrix is typically a few times larger. Therefore, in our simulations, a Monte Carlo sweep consists of several applications of the bisection algorithm until, on average, all beads are given a chance to have updated configurations.

Estimators
The PIGS method paves the way to evaluate the expected values of several properties of physical systems at T = 0. Unbiased estimates for local properties, such as the total, kinetic and potential energies, and the radial distribution function, can be obtained by constructing estimators from the corresponding operators. In particular, operators in coordinate representation have a direct estimator expression that follows from Eqs. (14) and (15). Once the sampling algorithm has achieved an equilibrium state, the configurations of the polymers are used to calculate a Monte Carlo ground state estimate O(τ) associated with operatorÔ following the expression where 〈. . . 〉 P denotes the average over uncorrelated configurations sampled from the probability distribution of Eq. (16) and O(R 0 , R 1 ) is given by Eq. (15).

Densities and radial distribution function
To obtain an operator expression for the radial distribution function, we initially consider the one-body number density operator, which is defined in terms of field operators aŝ Any state vector in the coordinate representation of a system of N particles can be written as such that the action of field operators in these states is given bŷ where |R + r〉 correspond to the configuration R with the addition of one particle at r, and where |R−r〉 correspond to the configuration R with the annihilation of one particle at r. Naturally, this only makes sense if there is one particle at r in the configuration R, which is where the delta functions stem from. We can see then thatd (1) (r) is a local operator, sincê therefore its estimator is given by where r j is the position of particle j in R 0 , and the average is over the probability distribution P loc from Eq. (18). For a uniform liquid, this final expression is simplified to d (1) (r; τ) = N /V , where V is the volume of the system. We can use the same approach to calculate the two-body number density operator, which is written in terms of field operators aŝ with the resulting estimator given by Finally, the radial distribution function operator, defined aŝ has, for a uniform liquid with spherical symmetry, the following estimator: where r = |r 1 − r 2 |. For systems with two spin components, g(r) can be further separated in different manners. A simple one is in up-and-down spins where the plus and minus signs correspond to g ↑↑ and g ↑↓ , respectively.

Total energy
Total energies can be calculated in several ways by identifyingÔ with the Hamiltonian of the system. If we follow the same recipe as the one used for the radial distribution function, we can show thatĤ is local in coordinate representation. The estimator expression, also known as the direct estimator, is given by where the derivative represents a Laplacian over all particle coordinates of R 0 . This estimator is not convenient since it involves second derivatives of density matrix elements.

Thermodynamic estimator
It can be more convenient to consider a different approach to estimate the total energy. Sincê ρ(δτ) = exp(−δτĤ), we can writeĤρ Therefore, the expected value of the total energy can be estimated via the so-called thermodynamic estimator This estimator has the advantage of requiring only a first-order, non-positional derivative of the density matrix.

Mixed estimator
Given the fact that [Ĥ,ρ] = 0,Ĥ does not need to be estimated in the central bead configuration R 0 . It is possible to estimate the system's total energy by applying the energy estimator to each one of the internal beads of the open polymer and, as a consequence, attenuate the statistical error associated with the estimation process. It is possible to moveĤ to one of the endpoints of the open necklace, which carries the trial wave function Ψ T , such that This leads to the ubiquitous mixed estimator of the Diffusion Monte Carlo (DMC) method [76].
The mixed estimator allows for the total energy of the system for a given configuration R M to be estimated without bias through the so-called local energy where Despite not being essential as in the DMC method, the total energy mixed estimator is beneficial in the PIGS method. The total energy can be estimated with half the number of beads necessary to estimate quantities that do not commute withĤ because filtering the ground state is only required from one of the extremities of the open necklace, and the Hermiticity ofĤ can be employed. Thus, this estimator offers a practical way to determine how many beads are needed in a given simulation to obtain converged configurations at the central bead, where most quantities are calculated.
It is important to mention some peculiarities of the mixed estimator. Suppose that the trial wave function Ψ T is not orthogonal to any of the true eigenstates of the system. In that case, the estimator will lead to vanishing variances of total energy estimates when Ψ T approaches the true eigenstate. Nevertheless, if Ψ T is not precisely the true eigenstate, a quantity O that does not commute with the Hamiltonian requires an extrapolation, O = 2O mix − O var , where O var is the variational estimate obtained with Ψ T [25]. When a simulation displays a zero-variance situation, the true wave function is known. In this case, the obvious approach is to perform a simple Monte Carlo simulation to sample the probability distribution directly if integrations cannot be performed analytically.

Potential energy
The potential energy can be immediately estimated using a direct estimator. A simple calculation gives us the expression which is straightforward to calculate during simulations.

Kinetic energy
If the potential energy operatorV is independent of the single particle mass m, one can estimate the kinetic energy by considering the expression to construct the thermodynamic estimator ofK, which is then given by A local estimator of the kinetic energy can also be constructed by directly using the quantum mechanical operator associated with this quantity. Nevertheless, it involves the Laplacian of the density matrix ρ(R 0 , R 1 ; δτ), which can be numerically unstable and have a high computational cost. In practice, to estimate the kinetic energy, the most accessible approach is to consider the straightforward estimation of the potential energy subtracted from the total energy estimate.

Results
As a proof-of-concept, the harmonic oscillator helps clarify several aspects involved in implementing the PIGS method for states described by symmetrical and anti-symmetrical wave functions. We then discuss 4 He atoms as a paradigm of strongly correlated many-body systems. It offers the possibility of showing characteristics of the PIGS method in its simplest form. Conversely, results for bulk 3 He show how our extension for Fermi particles can be applied in this prototype manybody system, establishing the basis for the treatment of any fermionic system. Finally, we show results for several 3 He concentrations in the 3 He-4 He mixture.

The harmonic oscillator as a proof-of-concept
As an initial application of the PIGS method, we employ a particle in a one-dimensional harmonic potential well in units of m = ω = ħ h = 1. The ground state of the system under these conditions has energy E 0 = 1/2 and is associated with an even wave function ψ 0 (x), i.e., ψ 0 (−x) = ψ 0 (x), ∀x ∈ ℜ. To verify the properties of the PIGS method, we choose a symmetrical trial state , where we assume b = 1/2. Although the functional form of Ψ T is the same as of the true ground state, b = 1/2 guarantees only a partial superposition of Ψ T (x) with ψ 0 (x) as long as Ψ T (x) is not orthogonal to the ground state. In our calculations, the parameter b was set to 1.1, which yields a variational energy of 0.665 ± 0.003. This value can be obtained in a standard variational Monte Carlo calculation or equivalently by imposing τ = 0 in the PIGS method. Of course, the exact value could also be calculated analytically. We tested two situations for the initial coordinates of the polymers: all beads starting at the origin and each bead starting at a different, random position. No difference in overall performance was observed between these two initial conditions. We attribute this to the simplicity of the system. The imaginary time evolution to the ground state was performed using δτ = 0.1 and the primitive approximation of Eq. (32) for the density matrix. These simulations were performed with the bisection algorithm using an L = 3 level. In a conventional PIGS run, we assigned 5% of the computational time to equilibration or thermalization to ensure detailed balance is achieved, with the remaining time used to estimate quantities of interest. Block averages of all quantities are aggregated to avoid correlations in the calculations of variances, and estimated errors are computed after the system has reached equilibrium.
The energy of the system was calculated using both the thermodynamic and the mixed estimator, Eq. (69) and Eq. (71), respectively. Both estimators give unbiased estimates, subjected only to statistical uncertainties. However, for small δτ values and considering a single bead, the energy standard errors associated with the thermodynamic estimator are larger than those corresponding to the mixed estimator. The larger error bars are a consequence of derivatives of the density matrix with respect to δτ, which includes a narrow Gaussian for small values of δτ. A ground state energy estimate of 0.50 ± 0.01 was obtained with the thermodynamic estimator, in excellent agreement with the analytical result E 0 = 1/2. For the mixed estimator, in The first excited state of the harmonic oscillator is associated with an odd wave function ψ 1 (x), i.e., ψ 1 (−x) = −ψ 1 (x), ∀x ∈ ℜ, which divides the x-axis into two regions of opposite signs, similar to the nodal regions of fermionic systems. For this reason, this state was chosen as a proof of concept to study the behavior of the PIGS method when it is mandatory to handle a state described by an anti-symmetric wave function.
The choice of the trial state Ψ T (x) = x exp(−b x 2 ), orthogonal to the ground state ψ 0 (x), allows the imaginary time evolution to converge to the first excited state. The variational energy obtained by setting the parameter b = 1.1 in this trial function is E(τ = 0) = 1.992 ± 0.006.

SciPost Physics Core
Submission E 1 = 3/2 is the exact value.
In this case, the fixed-node approximation must be considered if we expect a reliable description of the properties of the system, especially its energy. This approximation was performed by restricting the sampling algorithm to positive values of x, and employing the primitive density matrix approximation of Eq. (32) along with the image action correction of Eq. (21). Results are displayed in Fig. 6. We performed simulations without the image action correction to showcase its importance. The simulation details were similar to those of the ground state energy estimates. However, without the image action correction, a much shorter imaginary time step δτ = 0.01 was needed to obtain accurate results.
The imaginary time evolution towards the first excited state using the mixed estimator can be seen in Fig. 6. Results are notably better when the image action correction, represented by the circles, is used. Convergence is obtained for τ = 2. Even though the results obtained with and without the image action correction are indistinguishable within a 99% confidence interval, the introduction of the image action correction systematically improves the results and reduces the computational cost of each simulation. For τ = 5.0, the estimated energy with the image action correction is 1.51 ± 0.02, which is in excellent agreement with the first-excited state energy value. For calculations without the image action correction, the estimated energy is 1.56 ± 0.02. The oscillatory behavior of the results is related to the vanishing of probabilities in the vicinity of the nodal surfaces not being well captured without the correction.  Figure 6: Evolution of estimates of the first-excited state energy for the harmonic oscillator as a function of the imaginary time τ. The circles stand for results when the image action is employed. The squares represent these quantities without using the image action correction in the simulations. Error bars, when not visible, are smaller than the size of the symbol. The mixed estimator was used in both cases. The line shows the analytical value of the first excited-state energy E 1 = 3/2.
Since the density matrix of the harmonic oscillator is known exactly [77], features of the PIGS method related to properties of the thermodynamic and mixed estimators for the total energy can be further investigated. The exact action, allows large τ values to be used without requiring convolutions. The exact action is straightforward for the ground state, where we do not apply the fixed node approximation. A single bead is enough if the mixed estimator is applied to one of the wave functions at the extremities. The thermodynamic estimator was used only at the central bead of a necklace formed by three beads. This choice allowed converged configurations to the ground state at both sides of the central bead. It offered a hint of the difficulties encountered in estimating quantities that do not commute with the Hamiltonian of the system. We projected the ground state from the same even trial wave function used in the calculations with the primitive approximation. Estimates of the total energy using only configurations from the central bead as a function of τ are presented in Fig. 7. As already mentioned, for small τ values, the error bars associated with the thermodynamic estimator of the energy applied only to the central bead are considerably larger than the errors associated with the mixed estimator. Nevertheless, in this case, using the exact density matrix of the harmonic oscillator and with a fixed number of beads, the error bars decrease as τ increases since the derivatives of the density matrix with respect to τ include a broadening Gaussian. After convergence (τ = 5.0), the average energies obtained were E th = 0.4999±0.0002 for the thermodynamic estimator and E mix = 0.499 ± 0.003 for the mixed estimator. When the fixed-node approximation is used to study the first excited state, which requires an infinite positive potential for x ≤ 0, the action of Eq. (75) is no longer the exact action of the problem. As a result, a small value of δτ and convolutions of the density matrix is necessary for convergence. We projected the first excited state from the same odd trial wave function used in the calculations with the primitive approximation and performed simulations with and without the image action correction. In the former case, after convergence, the average energy obtained with the mixed estimator was 1.56 ± 0.02. When applied to the single central bead, a larger error bar is associated with the thermodynamic estimator, resulting in 1.38 ± 0.17. These results were improved when the image action correction was considered, with the thermodynamic estimator giving the energy 1.496±0.007 and the mixed estimator giving 1.502±0.006. In this simple model, the location of the node of the trial function is precisely known. Therefore, the true distance to the nodal surface, without approximation, was used to calculate the image action. For this reason, the convergence is to the true ground state energy.

Liquid helium
A complete description of mixtures of helium atoms is given by the Hamiltonian where m α is the isotope mass and N α is the number of atoms of each species. In most computer simulations, the two-body inter-atomic potential is proposed by Aziz and collaborators [78]. In our simulations, we use a Hartree-Fock damped form that mimics the entire configuration interaction in the intermediate-range, the HFD-B3-FCI1 potential [65], which gives excellent results in the short-, medium-, and long-range regions. The i < j sum in the interacting potential is over the total number of atoms in the system. Simulations of pure isotopes are performed by omitting the sum in the α index.
In studying isotope mixtures of helium atoms, it is interesting to have additional results for systems made from pure 3 He and 4 He atoms. Moreover, the PIGS method applied to a many-body bosonic system, where the difficulties associated with fermionic systems are not present, helps clarify practical aspects of the method.

Pure 4 He
The 4 He trial function adopted in the extremities of the necklace was chosen to be of the Jastrow form, The factor f (r i j ) = exp[−u(r i j )/2] explicitly correlates pairs of particles through a pseudo-potential of the McMillan form u(r i j ) = (b/r i j ) 5 , where b is a parameter [79]. Certainly, more sophisticated wave functions will favor a reduction in the number of beads considered for converged results and, consequently, make for a faster simulation. Our simulations were performed with 108 4 He atoms at density ρ = 0.02186 Å −3 . A convolution of the density matrix with 20 beads would produce converged configurations of the system ground state when the many-body action developed in Section (4.4) is considered. In other words, results obtained from a necklace of 41 beads guarantee that the central bead will have converged ground state configurations to estimate any property. Similar results can be obtained if the primitive approximation is used instead. However, more beads will be necessary to reach convergence. Table 1: Ground state total energy per particle of each pure system (first column) and densities in units of Å −3 (second column). Results were obtained using the thermodynamic E th (third column) and the mixed E mix estimators (fourth column). In the last columns, the number of particles N considered at each simulation is reported along with the experimental value. All the energies are in units of K. The equilibration stage was performed using 5% of the computational time of the run, with the remaining time employed in production stages. Again, block averages of all the quantities of interest are formed to avoid statistical correlations. The total ground state energy of the system obtained through the thermodynamic estimator of Eq.  Table 1. These results agree with other implementations of the PIGS method [30,33]. In the literature, different two-body interatomic potentials that may or may not consider three-body effects in an effective manner [78] are commonly used. Consequently, the predicted total energy per particle values have slight variations of the order of 2.6% [81][82][83][84][85].
In the upper panel of Fig. 8, we display the evolution of the ground state energy estimates with mixed and thermodynamic estimators as a function of the imaginary projection time. Although both estimates are statistically indistinguishable, those made with the thermodynamic estimator systematically show a lower value. This behavior is related to the fact that, in contrast to the mixed estimator, the thermodynamic estimator suffers from the finite character of the imaginary time steps since it depends on derivatives of the density matrix. Nonetheless, a strict agreement between different estimators is only expected for δτ → 0. Results obtained with both estimators can be combined, reducing the statistical uncertainty of a final estimate.
Estimates of the potential energy can easily be made at the central bead of the necklace following Eq. (72). As a result, we find V = −21.61 ± 0.01 K. The kinetic energy can then be obtained by subtracting the potential energy from the total energy. For this purpose, we can use either the thermodynamic or the mixed estimator. We found K 4 = 14.31 ± 0.01 K, which is in excellent agreement with experimental data [80,88,89] and with theory [30,90]. We can also calculate the kinetic energy through the thermodynamic estimator of Eq. (74), K 4 = 14.2 ± 0.1 K, a result in agreement with the previous estimate within the statistical errors.

Pure 3 He
We now focus on a strongly correlated Fermi system formed from 3 He atoms, which also has a long history in the literature [60,91]. For the liquid 3 He system, the trial function at the extremities of where with and η(r) given by Here λ B , s B , ω B , and λ B are variational parameters. The ↑ (↓) symbol corresponds to the spin up (down) configuration.
To study finite size effects in our simulations, we have considered pure unpolarized systems of 54, 66, and 114 3 He atoms (corresponding to closed Fermi shells) with the fixed-node approximation and the image action correction. Results for the total energy per particle of the system are presented in Table 1. At density ρ = 0.01635 Å −3 , the results obtained with the thermodynamic estimator applied to the central bead are in excellent agreement with the mixed ones. For the N = 114 case, we obtained E th = −2.34 ± 0.02 K and E mix = −2.35 ± 0.01 K. These results agree, within statistical uncertainties, with the ones obtained with N = 66 particles, suggesting N = 66 a suitable number for pure 3 He simulations. However, these results do not agree with the experimental value of −2.47 ± 0.01 K [87]. As predicted in Refs. [95,96], we found finite size effects in the total energy of approximately 0.1 K between the N = 54 and N = 114 simulations. Nevertheless, the kinetic energy for a simulation with N = 114 atoms, displayed in Table 2, K 3 = 12.34 ± 0.02 K, is in excellent agreement with the experimental result from Ref. [94]. The estimated potential energy is V = −14.68 ± 0.01 K. Finally, the lower panel of Fig. 8 shows the evolution of the binding energy per 3 He atom as a function of the projection time for the N = 54 simulation. As expected, fluctuations in the computed values of the total energy of 54 3 He atoms are bigger than those observed for the system of 108 4 He atoms (upper panel).
Continued efforts spanning the last four decades to measure the kinetic energy of pure 3 He face several experimental challenges [93,94,[97][98][99][100][101][102]. In Fig. 9, our estimate is compared with some of the most recent experimental and theoretical [92,[103][104][105] values from the literature. The first attempt to investigate a fermionic system using the PIGS method took an approach that resembles a "released-node" simulation, resulting in a lower kinetic energy [11]. Conversely, considering the Fermi statistics through all beads, as in our approach, yields more accurate results.   [94].
When studying the 3 He liquid phase, the total pair correlation function and its decomposition into parallel and antiparallel spin components, Eq. (66), is a quantity of interest. In Fig. 10, we displayed our results for the pure system obtained using the converged configurations of the central bead of the necklace. The curve that considers atoms of antiparallel spins has a more pronounced peak since such atoms do not experience the Pauli exclusion principle. r (Å) Figure 10: Total pair correlation functions (solid red line) and its components, the antiparallel (dashed blue line) and the parallel spins functions (dashed-dotted magenta line) versus the radial distance. The yellow shaded area shows the region where the pair correlation of antiparallel spin is larger than that of parallel spins. Between r = 8 Å and r = 10 Å there is a region where the antiparallel pair function is smaller than the one for parallel spins. This region is not easily seen in the figure scale, even though it is depicted in green.

The 3 He-4 He mixture
A paradigm of strongly correlated quantum many-body systems composed of both bosons and fermions are mixtures of 3 He and 4 He atoms. Despite properties of bulk 4 He and 3 He being well understood, in the 3 He-4 He mixture, there is still disagreement between experimental data and theoretical results. Such inconsistencies are demonstrated by a recent study of the kinetic energies of each one of the components in the liquid phase at finite temperature [10]. The observed discrepancy for relative concentrations above x = 0.20 was partly attributed to the lack of Fermi statistics to describe the 3 He atoms; they were treated as distinguishable particles. In this context, it is crucial to perform a study where the appropriate quantum statistics treat both components to observe if improvements in the description of the system can be achieved. With this goal, in simulations using the FN-PIGS method, we describe the liquid helium mixture by a wave function of the Jastrow-Slater form, Eq. (78), which explicitly includes Fermi statistics and back-flow correlations together with the Bose statistics for the 4 He atoms. However, in the present case, the pair function of the Jastrow form correlates all particles in the mixture 3 He-4 He, with each species and inter-species pair making use of its own variational parameter b.
Furthermore, the Slater determinant will now also depends on the 4 He coordinates through the back-flow correlation, increasing the richness of the nodal surface of the system. The parameters we considered for the back-flow correlation were the same as in the description of bulk 3 He.
We performed simulations with 108 atoms in the mixture and considered different concentrations of 3 He atoms, x = N3 He /(N3 He + N4 He ). The equilibrium density of the mixture for each concentration x was calculated using the relation ρ = [x m 3 + (1 − x)m 4 ]/V , where V is the volume of the simulation cell.
Results for the kinetic energy of both species can be seen in Fig. 11, and the values are displayed in Table 2. The kinetic energies of 3 He show an improvement in the direction of experimental results, ranging from about 1.5 to 0.4 K from the lower to the higher concentrations, compared with results from the literature [10]. Nevertheless, agreement with experimental data is still lacking because, while these seem to fluctuate around 11 K, theory shows values that decrease with increasing concentrations. However, it is puzzling that, for pure 3 He, there is an excellent agreement between theoretical and experimental values. Estimated values of the kinetic energy for the 4 He component not only improve literature results but are also in excellent agreement with experimental data. The result for pure 4 He corroborates that estimates of kinetic energies made with PIGS are very accurate for bosonic particles. x Figure 11: Kinetic energies as a function of 3 He concentration x. The solid squares and triangles represent our 4 He and 3 He results, respectively. The empty symbols refer to the most recent experimental data. For 3 He: empty triangles are from [93]; the upside down triangle is from [94]; the pentagons are from [102]. At x = 1.0 (pure 3 He), our result and an experimental value from [102] are indistinguishable in the figure. For 4 He: diamonds are from [93]; the empty square is from [88]. The diamond at x = 0.35 and the empty square at x = 0.0 (pure 4 He) are indistinguishable from our results; lines are guides to the eye.

Conclusions
We propose an extension of the PIGS method by incorporating the fixed-node approximation, calling it FN-PIGS. As in the bosonic case, FN-PIGS allows estimations without the need for variational results for the extrapolation procedure of any property, regardless of whether it commutes with the Hamiltonian of a fermionic system. Moreover, further developments, such as those that could be achieved through quantum computing, offer a possibility to remove the bias coming from the fixed-node approximation in hybrid quantum-classical algorithms [75], which would make any FN-PIGS estimate numerically exact. It is also possible that other attempts to reduce or even eliminate the sign problem can be easily incorporated into the method [9,27,37,106]. An essential feature of FN-PIGS is that it does not rely on importance sampling transformations [107]. Moreover, the method does not require the introduction of a periodically updated energy shift parameter E T to stabilize a population of configurations or walkers. This situation introduces bias and is known as the population-control error [108,109]. Additionally, one can directly distribute parallel processes and periodically collect averages to obtain precise results with increasingly smaller statistical errors. This approach avoids difficulties associated with strategies where thousands of walkers are submitted to branching algorithms, eventually resulting in the necessity of load balance among processors for efficiency.
We have presented results for a proof-of-concept scenario with the harmonic oscillator to clarify several aspects of implementing the FN-PIGS method. Key features of the method were conveyed by examining properties of mixtures of normal 3 He in superfluid 4 He, where both the Fermi-Dirac and the Bose-Einstein statistics intervene, and both species receive a quantum treatment. Our results for different 3 He concentrations in 4 He show an improvement in the direction of the experimental values. However, a significant difference between our results and the experimental ones remains for the 3 He kinetic energies. Indeed, a better knowledge of how the nodal structure of the 3 He isotope evolves in a mixture with 4 He is necessary to obtain more accurate results. As put forward in the literature, this disagreement could be related to an underestimation of the 3 He kinetic energy contribution associated with the tail of the measured momentum distribution.
Nonetheless, the kinetic, potential, and total energies obtained with FN-PIGS for pure liquid 4 He and 3 He systems at their equilibrium densities agree with the reported experimental and theoretical data. The investigation of the characteristics of fermionic nodes in wave functions is an active research area [110], and hopefully, advances in this direction will contribute to the analysis of 3 He-4 He mixtures.
Our results for pure 3 He and 3 He-4 He mixtures raise the interesting question of how the dilution of 3 He in superfluid 4 He is responsible for shaping the 3 He nodal structure. After all, results for pure 3 He and pure 4 He are more accurate than those obtained for the mixture in the sense that theory and experiment are in excellent agreement. Indeed, the 3 He-4 He mixture is fascinating because, despite being one of the simplest combinations of bosons and fermions, it still offers several challenges for experiments and theory. Results from quantum Monte Carlo methods point to the direction of a 4 He kinetic energy that decreases with an increasing concentration of 3 He. Experimental data of 3 He for this quantity still leave some margin for interpreting its behavior. This scenario makes explicit the need for further studies to understand these systems comprehensively.
FN-PIGS is a robust method that gives accurate results for several physical properties in strongly correlated quantum many-body systems. The method does not rely on variational results to estimate properties like the contact parameter in ultracold quantum cases [111][112][113], making it a valuable tool in the toolbox of quantum Monte Carlo methods. Potential applications to other systems are not hard to find. For instance, FN-PIGS can be applied in the investigation of neutron matter and ultracold Fermi gases. In nuclear matter at densities lower than 0.003 fm −3 , where the complexities of the asymmetric nuclear Hamiltonian can be neglected, a central potential can be used. For quantum Fermi gases, the equation of state, or the spectral weight when an impurity is introduced in a polarized medium, and other properties of interest can be estimated without any extrapolation. For reviews on these topics, see Refs. [114,115].
As a final remark, in naming the method fixed-node path-integral ground state (FN-PIGS), we follow the standard nomenclature found in the literature when the fixed-node approximation is incorporated into a given existing method. However, we believe that Density Matrix Projection (DMP) better represents how the method operates and reflects its main capabilities. In particular, results can converge to an excited state if, in the extremities of the necklace, the wave functions are orthogonal to the true ground state of the Hamiltonian.