Neural Wave Functions for Superfluids

Understanding superfluidity remains a major goal of condensed matter physics. Here we tackle this challenge utilizing the recently developed Fermionic neural network (FermiNet) wave function Ansatz [D. Pfau et al., Phys. Rev. Res. 2, 033429 (2020).] for variational Monte Carlo calculations. We study the unitary Fermi gas, a system with strong, short-range, two-body interactions known to possess a superfluid ground state but difficult to describe quantitatively. We demonstrate key limitations of the FermiNet Ansatz in studying the unitary Fermi gas and propose a simple modification based on the idea of an antisymmetric geminal power singlet (AGPs) wave function. The new AGPs FermiNet outperforms the original FermiNet significantly in paired systems, giving results which are more accurate than fixed-node diffusion Monte Carlo and are consistent with experiment. We prove mathematically that the new Ansatz, which only differs from the original Ansatz by the method of antisymmetrization, is a strict generalization of the original FermiNet architecture, despite the use of fewer parameters. Our approach shares several advantages with the original FermiNet: the use of a neural network removes the need for an underlying basis set; and the flexibility of the network yields extremely accurate results within a variational quantum Monte Carlo framework that provides access to unbiased estimates of arbitrary ground-state expectation values. We discuss how the method can be extended to study other superfluids.


I. INTRODUCTION
The unitary Fermi gas (UFG) is a paradigmatic example of a strongly interacting system of two-component fermions that possesses a superfluid ground state and lies in the crossover region between a Bardeen-Cooper-Schrieffer (BCS) superconductor and a Bose-Einstein condensate [2,3].The effective range of the interaction is zero and the s-wave scattering length diverges (the "unitarity limit"), so the UFG has no intrinsic length scale.The only remaining length is the inverse of the Fermi wavevector 1/k F , on which all thermodynamic quantities depend.For example, regardless of the particle density, the ground-state energy per particle of a unitary Fermi gas can be written as where E F G is the energy per particle of a non-interacting Fermi gas of the same density.The dimensionless constant ξ is known as the Bertsch parameter [4].
Because of the universality of the UFG model, it can be used to describe many real physical systems at different scales, such as the neutron matter in the inner crust of a neutron star [5] or the quantum criticality of an s-wave atomic superfluid [6,7].The size of the pairs in the UFG is comparable to the inter-particle spacing, which is also a feature of many high-T c superconductors [8][9][10].As a result, the UFG has been studied extensively [11].Although the UFG is an idealized model, it can be accurately realized in the laboratory using ultracold atomic gases in which the interactions have been tuned by using an external magnetic field to drive the system across a Feshbach resonance [12].
The UFG has been studied for decades, but it remains difficult to calculate its ground-state properties accurately using analytic methods.Mean-field treatments such as BCS theory [13] give good results for systems with weak interactions, but there is no guarantee of success in the strongly interacting regime.As a result, various quantum Monte Carlo (QMC) methods [14,15] have been used to simulate the properties of the UFG to high accuracy at zero and finite temperature.Methods used include variational Monte Carlo (VMC), fixed-node diffusion Monte Carlo (FN-DMC), fixed-node Green function Monte Carlo, auxiliary field Monte Carlo and diagrammatic Monte Carlo [16][17][18][19][20][21][22][23][24][25][26].However, a full quantitative description remains an open and challenging problem.
Recent advances in machine learning algorithms and the growing availability of inexpensive GPU-based computational resources have allowed neural-network-based approaches to permeate many areas of computational physics, including lattice [27][28][29][30] and continuum [1,[31][32][33] QMC simulations.Here we employ a neural network Ansatz within a VMC approach to study the unitary Fermi gas.The Ansatz we use, the Fermionic Neural Network (FermiNet) [1], gives very accurate results for atoms and molecules [1,[34][35][36] and has recently been applied to periodic solids and the homogeneous electron gas (HEG) with comparable success [37].In the case of the HEG, the variational optimization of the FermiNet Ansatz discovered the quantum phase transition between the Fermi liquid and Wigner crystal ground states without external guidance [38].In contrast, previous approaches required different Ansätze to be used for the two different phases.The FermiNet has not previously been applied to fermionic superfluids such as the UFG.
The paper is organized as follows.Section II describes the architecture of the FermiNet.We find that the original FermiNet Ansatz is insufficient to capture the twoparticle correlations of superfluids.Although a FermiNet wave function with one determinant and a sufficiently large neural network is in principle able to represent any fermionic state [1], it is often advantageous to use a network of a fixed size and a small linear combination of Fer-miNet determinants.In the case of the unitary Fermi gas, however, we find that the number of block-diagonal determinants required to describe the ground state accurately scales exponentially with the system size.This is the first example in which the FermiNet has been seen to fail both quantitatively and qualitatively, and suggests that the FermiNet wave function may not be able to represent arbitrary fermionic wave functions in practice.To remedy the problem, we utilize the neural-network part of the FermiNet architecture to build a different type of wave function based on the idea of an antisymmetric geminal power singlet wave function (AGPs) [30,[39][40][41][42], which we discuss in detail in Section III.This leads to substantial improvements, even though the neural-network part of the wave function remains unchanged.The implementation of the AGPs wave function using the FermiNet, as well as its relation to the original block-diagonal multideterminant FermiNet, are discussed in Section IV.Our computational results are presented in Section V, followed by a summary and discussion in Section VI.The Appendix includes detailed explanations and derivations of important formulae, as well as implementation and training details.

II. FERMINET
The Fermionic Neural Network, or FermiNet [1], is a neural network that can be used to approximate the ground-state wave function of any system of interacting fermions.The inputs to the network are the positions r 1 , r 2 , . . ., r N and spin coordinates σ 1 , σ 2 , . . ., σ N of the N particles, and the output is the value of the wave function Ψ(r 1 , σ 1 , r 2 , σ 2 , . . ., r N , σ N ) corresponding to those inputs.The network is trained using the variational Monte Carlo (VMC) method [14]: the weights and biases that define the network are varied at each training iteration to minimize the energy expectation value according to the variational principle.If the network is flexible enough, the approximate wave function obtained after training may be very close to the true ground state.The FermiNet provides a more general and accurate alternative to the conventional Slater-Jastrow (SJ) and Slater-Jastrow-backflow (SJB) Ansätze that have been used in most VMC and FN-DMC calculations to date, and may improve VMC and FN-DMC results for strongly correlated systems.
In conventional SJ Ansätze, the antisymmetry of the N -electron wave function is represented using Slater determinants, which are antisymmetrized products of single-particle orbitals.For simulations of solids, it is common to use one determinant only; for molecules, a linear combination of determinants is usually employed.In both cases, the presence of determinants guarantees that the wave function has the correct exchange antisymmetry.To improve the representation of electronic correlations, especially the correlations that chemists call "dynamic", the determinants are multiplied by a totally symmetric non-negative function of the electron coordinates known as a Jastrow factor.This acts to decrease the value of the wave function as pairs of electrons approach each other, reducing the total Coulomb repulsion energy.
If the Hamiltonian is independent of spin and all the single-particle orbitals are eigenfunctions of total S z , one can assign spins to the electrons and every Slater determinant can be factored into a product of spin-up and spin-down Slater determinants [14,43].The wave function is no longer antisymmetric under the exchange of electrons of opposite spin, but expectation values of spinindependent operators are unaltered.Including a spinassigned Jastrow factor expressed in the form e J , a onedeterminant SJ Ansatz becomes: where {r ↑ } and {r ↓ } are the sets of position coordinates of the N ↑ electrons assigned to be spin up and the N ↓ electrons assigned to be spin down, respectively.One can improve the SJ Ansatz by transforming the electron coordinates as where α and ᾱ are the two possible spin components of an electron, r βα ij = |r β i − r α j |, and η ∥ (r) and η ∦ (r) are parameterized functions of a single distance argument.The coordinate-transformed SJ Ansatz is called a Slater-Jastrow-backflow (SJB) wave function, and the new coordinates are called quasiparticle coordinates.Note that the quasiparticle coordinate x α j is invariant under the exchange of any two position vectors in {r α /j } = {r α 1 , . . ., r α j−1 , r α j+1 , . . ., r α N α } or in {r ᾱ} = {r ᾱ 1 , . . ., r ᾱ N ᾱ } [14].
The backflow transformation replaces every singleparticle orbital ϕ α i (r α j ) by a transformed orbital ϕ α i (x α j ), which depends on the position of every electron in the system.Exchanging the coordinates of any two spinparallel electrons still exchanges two rows of the Slater determinant, so the antisymmetry is preserved.The downside is that moving one electron now changes every element of the Slater matrix, preventing the use of efficient rank-1 update formulae and increasing the cost of re-evaluating the determinant by a factor of N .Despite the extra cost, however, the enrichment of the description of correlations between electrons makes SJB wave functions significantly better than SJ wave functions and they are frequently used in VMC and FN-DMC simulations.
The FermiNet [1] takes the idea of permutation equivariant backflow much further, replacing the orbitals ϕ α i (r α j ) entirely by neural networks.The orbitals represented by these networks differ from SJB orbitals because they are not functions of a single three-dimensional vector x α j but depend in a very general way on r α j and all the elements of the sets {r α /j } and {r ᾱ}.They are best written as ϕ α i (r α j ; {r α /j }; {r ᾱ}).The exchange antisymmetry is maintained because ϕ α i (r α j ; {r α /j }; {r ᾱ}) is totally symmetric on exchange of any pair of coordinates in {r α /j } or {r ᾱ}.Furthermore, because they are represented as neural networks, the FermiNet orbitals need not be expanded in terms of an explicit basis set, widening the class of functions they can represent [44].In order to build functions with the correct exchange symmetry properties, a carefully constructed neural network architecture is used, which is described below.
The FermiNet architecture consists of two parts: the one-electron stream, which takes electron-nucleus separation vectors r α i − R I and distances |r α i − R I | as inputs, and the two-electron stream, which takes electronelectron separations r α i − r β j and distances |r α i − r β j | as inputs, with i, j ∈ {1, 2, . . ., N α } and α, β ∈ {↑, ↓}.The inputs to the one-electron stream are concatenated to form one input vector for each electron, and the inputs to the two-electron stream are concatenated to form one input vector for each pair of electrons: where the superscript 0 means that the vectors are the inputs to the first layer of the network.The distances between particles are passed into the network to help it to model the wave function cusps, i.e., the discontinuities in the derivatives of the wave function when two electrons or an electron and a nucleus coincide.These discontinuities create divergences in the kinetic energy that exactly cancel the divergences in the potential energy as pairs of charged particles approach each other [1].Each electron stream consists of several layers.At each layer l ∈ {0, . . ., L−1}, the outputs h lα i and h lαβ ij from the streams are averaged and concatenated in the following way: The concatenated one-electron vectors are then passed into the next layer, as are the two-electron vectors: where V l and W l are matrices, b l and c l are vectors, and all of them are optimizable.We denote the number of hidden units in each layer in the one-electron stream by n l such that h lα i ∈ R n l , l ∈ {0, 1, 2, . . ., L} [45].The outputs from the final layer L of the one-electron streams are used to build the many-particle FermiNet orbitals: where w kα i is an optimizable vector and g kα i an optimizable scalar.The χ kα i (r α j ) factor is an envelope function to ensure that the wave function satisfies the relevant boundary conditions.For example, in a system which requires the wave function to tend to zero as |r α j − R m | → ∞, exponential envelopes are used: where π kα im and σ kα im are variational parameters.No attempt is made to ensure that the FermiNet orbitals are normalized or orthogonal to each other.
As mentioned earlier, FermiNet orbitals are not functions of one electron position r α i only, but also depend on the positions of all the other electrons in the system in an appropriately permutation invariant way.No Jastrow factor is needed as the electron-electron correlations are included in the network.The full wave function is thus a block-diagonal determinant (BD) of the Fer-miNet orbitals ϕ kα i (r α j ; {r α /j }; {r ᾱ}).Multiple determinants may also be used, in which case the wave function is a weighted linear combination where the superscript D specifies the number of determinants of FermiNet orbitals in the linear combination of determinants that makes up the full wave function, and the "Slater FermiNet" subscript serves to specify this specific wave function Ansatz and is discussed in more detail below.In practice, the weights ω k are absorbed into the orbitals, which are not normalized.The VMC method is then applied to the FermiNet Ansatz and the parameters of the network are optimized using a second-order method known as the Kroneckerfactored approximate curvature algorithm [46].The aim is to minimize the expectation value of the Hamiltonian ⟨H⟩, which acts as our loss function.For a more detailed explanation of the FermiNet architecture, see Pfau et al. [1] and the discussion of the improved JAX implementation [47] in Spencer et al. [36].
Consider the basis {a 1 , a 2 , a 3 } of the Bravais lattice generated by repeating the finite simulation cell periodically.Any position vector may be written as r = s 1 a 1 +s 2 a 2 +s 3 a 3 .To ensure that the FermiNet represents a periodic function, the position coordinates s i are replaced in the FermiNet inputs by pairs of periodic functions, s i → (sin(2πs i ), cos(2πs i )).Thus, if any electron is moved by any simulation-cell Bravais lattice vector, the inputs to the network are unchanged.It follows that the output, the value of the wave function, is also unchanged.A periodic envelope function is used to improve the speed of convergence [38]: where the k m are simulation-cell reciprocal lattice vectors up to the Fermi wavevector of the non-interacting Fermi gas.This specific way of adapting the FermiNet to periodic systems was proposed by Cassella et al. [38], although other similar methods exist [37,48].
The FermiNet has only been used to study systems of electrons interacting via Coulomb forces to date, but can easily be adapted to systems of other spin-1/2 particles simply by changing the Hamiltonian.Here we use the periodic FermiNet Ansatz to approximate the ground state of the UFG Hamiltonian in a cubic box subject to periodic boundary conditions.Since there are no atomic nuclei and the wave function has no electron-nuclear cusps, the inputs to the one-electron streams are simpler than shown in Eq. ( 4), containing only the particle coordinates [49]: h 0α i = r α i , with respect to an origin placed at one corner of the simulation cell.A detailed discussion of translational symmetry of the wave function can be found in section E of the Appendix.
As will be demonstrated below, the Slater FermiNet is sufficient to learn the superfluid ground state for small systems but fails for large systems.Hence, we propose a modification to the method of building orbitals.The motivation for this modification comes from earlier work using antisymmetrized products of two-particles orbitals known as antisymmetrized geminal power (AGP) wave functions [16,30,39,40,42,50,51].We describe the an-tisymmetrized geminal power singlet (AGPs) wave function in the next section.
The authors of Ref. [1] used the term "FermiNet wave function" to refer to all wave functions constructed using a FermiNet neural network.Now that we are going to use almost the same neural network to generate AGPsbased pairing wave functions in addition, more precise terminology is required.Wave functions of the type introduced in Ref. [1], which contain many-particle generalizations of the one-particle orbitals that appear in Slater determinants, will be called one-determinant or multideterminant Slater FermiNets.Wave functions built using determinants of many-particle generalizations of pairing functions will be called one-determinant or multideterminant AGPs FermiNets.Since every AGPs Fer-miNet determinant is built using one pairing function or "geminal", we also refer to one-geminal or multi-geminal AGPs FermiNets.

III. ANTISYMMETRIZED GEMINAL POWER WAVE FUNCTION
The FermiNet and other Ansätze that expand the ground state as a linear combination of Slater determinants give very accurate results for many molecules and solids, but may still fail to capture strong two-particle correlations in superfluids.An alternative starting point, which is better at capturing two-particle correlations, is the antisymmetric geminal power (AGP) wave function [40][41][42]52].This uses an antisymmetrized product of two-particle functions known as pairing orbitals or geminals instead of an antisymmetrized product of singleparticle orbitals.
Although one can build a general AGP wave function with pairings between arbitrary particles, the UFG Hamiltonian only contains interactions between particles of opposite spin.It is therefore sufficient to consider pairing orbitals involving particles of opposite spin only.In this case, the wave function is called an antisymmetrized geminal power singlet (AGPs).The rest of this section summarizes the main features of the AGPs Ansatz and explains how the FermiNet architecture can be modified to produce many-particle generalizations of AGPs pairing orbitals.Detailed discussions of AGPs wave functions, including derivations of the equations, can be found in Refs.[30,[40][41][42]52] and the Appendix.

A. AGP Singlet Wave Functions
It is helpful to start by considering an unpolarized system with an even number (N = 2p) of particles and total spin S z = 0.An AGPs wave function for such a system is constructed using a singlet pairing function of the form where φ(r i , r j ) is a symmetric function of its arguments.We work with spin-assigned wave functions, so we set the spins of particles 1, 2, . . ., p to ↑ and the spins of particles p + 1, p + 2, . . ., 2p to ↓. If, for example, i ≤ p and j > p, so that particle i is spin-up and particle j is spin-down, the spin-assigned pairing function is The spin-assigned singlet pairing function is equal to zero if the spins of particles i and j are the same.The spin-assigned AGPs wave function is a determinant of spatial pairing functions [16,39]: Like all spin-assigned wave functions, it depends on position coordinates only.For convenience, we have changed the particle labeling scheme: i and j now run from 1 to p and arrow superscripts have been added to distinguish up-spin from down-spin particles.Note that the AGPs wave function coincides with the BCS wave function projected onto a fixed particle number subspace (see Ref. 39 and Appendix A 1).It is therefore suitable for describing singlet-paired systems, including s-wave superfluids.

B. AGPs with Unpaired States
We can generalize the spin-assigned AGPs wave function to allow for unpaired particles.Consider a system with N = 2p + u + d particles, where p is the number of pairs, u is the number of unpaired spin-up particles, and d is the number of unpaired spin-down particles.The total number of spin-up particles is p + u and the total number of spin-down particles is p + d.The AGPs wave function can be written as a determinant of pairing functions and single-particle orbitals [39,42,52] as shown in Eq. ( 16) where φ(r ↑ i , r ↓ j ) is an arbitrary singlet pairing function and ϕ σi i (r σi j ) are arbitrary single-particle functions.For the UFG considered in this paper, we only need the case where u = 1 and d = 0 or vice versa.This represents a fully paired 2p-particle system to which one particle has been added.

IV. AGP SINGLET FERMINET
Having discussed the form of the AGPs wave function, we now discuss how it can be implemented using FermiNet.In the original Slater FermiNet architecture, the outputs of the one-electron stream are used to build FermiNet orbitals ϕ kα i (r α j ; {r α /j }; {r ᾱ}).The full manyparticle wave function is a weighted sum of terms, each of which is the product of one up-spin and one down-spin determinant of the FermiNet orbital matrices, as shown in Eq. (11).
To build a many-particle pairing function using the neural-network part of FermiNet, one can make use of its outputs h Lα i from the last layer L of the one-electron stream.Instead of using these outputs to build FermiNet orbitals, as in Eq. ( 9), they can be used to build FermiNet pairing orbitals, also known as FermiNet geminals: where χ k (r) are the envelope functions, w k are vectors, g k a scalar, and ⊙ denotes the element-wise product.Note that the same FermiNet geminal is used for all pairs of particles, so the envelope functions in Eq. ( 17) do not require the particle and spin indices that appear in the envelope functions of the FermiNet orbitals defined in Eq. ( 9).This construction generates a manyparticle pairing function between particles r α i and r ᾱ j , retaining the permutation invariant property possessed by FermiNet orbitals.Depending on the number of Fer-miNet geminals generated, the wave function can be written as one or a weighted sum of multiple determinants of FermiNet geminals, where the superscript D in Ψ D AGPs FermiNet specifies the number of determinants (and thus FermiNet geminals) appearing in the linear combination that makes up the wave function.This is analogous to a weighted sum of conventional single-determinant AGPs wave functions of the type defined in Eq. ( 15), but the replacement of the two-particle pairing orbitals by FermiNet geminals that depend on the positions of all the particles makes it much more general.
Although using the outputs from the one-electron stream is sufficient to build an AGPs, one can also include the outputs from the two-electron stream: where w k 1 and w k 2 are vectors.Note that Eqs. ( 17) and ( 19) are two possible ways of building a many-particle pairing function.There are many others ways and they are all valid as long as the appropriate symmetries are preserved.An alternative method is given by Xie et al. [53].
The benefit of building AGPs-like wave functions using the FermiNet is that the many-particle pairing function φ(r α i , r ᾱ j ; {r α /i }; {r ᾱ /j }) now depends not only on r α i and r ᾱ j but also on the positions of the other particles in the system [54].Correlations between the singlet pair and the other particles can thus be captured.In a similar way, the original Slater FermiNet replaced Hartree-Fock-like single-particle orbitals ϕ α i (r α j ) by many-particle orbitals (FermiNet orbitals) ϕ α i (r α j ; {r α /j }; {r ᾱ}), helping to capture correlations between the particle at r α j and all other particles.

A. Relations between the Slater FermiNet and the AGPs FermiNet
Next, we clarify the relation between the Slater Fer-miNet with block-diagonal determinants and the AGPs FermiNet, showing that the AGPs FermiNet is the more general of the two.A FermiNet geminal with a twoparticle stream term is even more general than a Fer-miNet geminal without, so it is sufficient for this purpose to omit the two-particle stream term.We also neglect the envelope functions and the bias term, g k , which is set to 0 in all results presented here.The use of envelope functions circumvents numerical difficulties in finite systems and can speed up the network optimization, but does not affect the generality of the Ansatz.
Let us first define a many-particle pairing function in the following way: where h

Lα(k) i
= [h Lα i ] k are the outputs from the final layer of the one-electron stream for particle i of spin α.As we explain below, Eq. ( 20) is equivalent to the simpler FermiNet geminal described above: We choose to write the many-particle pairing function in the form of Eq. ( 20) only because this makes it easier to relate to FermiNet orbitals.Since Eqs. ( 20) and ( 21) are equivalent, the choice does not affect the conclusions of the argument.In the rest of this section, for the sake of simplicity, we omit the sets {r ↑ /i } and {r ↓ /j } from the arguments of the many-particle pairing functions and orbitals.
To explain the equivalence of Eqs. ( 20) and ( 21), it is helpful to represent the matrix W kl as its singular-value decomposition (SVD): where U ∈ R n L ×n L and V ∈ R n L ×n L are orthogonal matrices and n L is the size of the vectors h Lα i output by the final layer L of the one-electron stream.This is also known as the number of hidden units in layer L. The many-particle pairing function in Eq. (20) becomes Given the universal approximation theorem [44], and the fact that every layer of the network contains an arbitrary linear transformation, it is reasonable to assume that the functions h Lα i and Oh Lα i , where O is U or V , have the same variational freedom and information content.In other words, we assume that any network capable of representing h Lα i can also represent Oh Lα i , since this is merely a rotation of the vectors h Lα i in the last layer.We thus define hL↑ i = U h L↑ i and hL↓ i = V h L↓ i , such that the many-particle pairing function becomes which is equivalent to Eq. ( 21).
To relate the AGPs FermiNet and the Slater FermiNet, we expand an AGPs determinant constructed using the many-particle pairing function from Eq. ( 20) as a sum of block-diagonal determinants of FermiNet orbitals.It will be sufficient to consider matrices W ∈ R n L ×n L of rank M , with p ≤ M ≤ n L .We can decompose any such matrix using rank factorization, where F and G are matrices in R M ×n L .Equation ( 20) then becomes where the last line defines the functions ϕ ↑ γ and ϕ ↓ γ .In the case when M = p, where p = N/2 is the number of pairs in the system, the determinant of the many-particle pairing function can be written as a block-diagonal determinant of FermiNet orbitals: where [ϕ α ] ij = ϕ α i (r α j ) are matrices in R M ×p with M = p.The product of two p × p determinants can be written as the determinant of a single 2p × 2p matrix, with the p × p spin-up and spin-down blocks on the diagonal.Therefore, a single-geminal AGPs FermiNet wave function constructed using the many-particle pairing function from Eq. ( 20) with a rank-p matrix W kl is equivalent to a 2p × 2p block-diagonal determinant of FermiNet orbitals.The equivalence is already well known [51] for AGPs wave functions constructed using conventional two-particle orbitals.
Now consider the more general case where p ≤ M ≤ n L .The Cauchy-Binet formula states that where the sum is over all M p distinct choices of p rows from the two M × p matrices ϕ ↑ γ (r ↑ i ) and ϕ ↓ γ (r ↓ i ).The products of the determinants of the two p × p matrices associated with each such choice are summed to reproduce the AGPs.This is similar to the linear combination of multiple block-diagonal-determinants of FermiNet orbitals without weights given by Eq. ( 11) and in the original FermiNet paper [1] [55].
Note that the intermediate layers, i.e., the one and two-electron streams, are identical in the Slater FermiNet and the AGPs FermiNet.The only modifications are made at the orbital shaping layer, or, equivalently, the method of antisymmetrization has changed.Thus, the representational power of the intermediate layers of the AGPs FermiNet remains the same as for the Slater Fer-miNet.Thus, it must be the method of antisymmetrization that limits the performance of the Slater FermiNet when applied to the UFG.
We have shown that a single AGPs determinant con-structed using the many-particle pairing function from Eq. ( 20) with a matrix W kl of rank greater than p contains multiple block-diagonal determinants of FermiNet orbitals.If the rank of W kl is equal to p, the AGPs is equivalent to a single block-diagonal FermiNet determinant.Conversely, any single-determinant FermiNet wave function can be written as an AGPs of rank p.Therefore, the AGPs FermiNet provides a more powerful Ansatz with fewer variational parameters than the Slater Fer-miNet, since the former contains the latter.
It is worth mentioning another advantage of using the FermiNet to build geminals.By generating more sets of independent parameters w k i in Eq. ( 19), one can easily construct an arbitrary number N det of FermiNet geminals φ k (r α i , r ᾱ j ; {r α /i }; {r ᾱ /j }) with k ∈ {1, 2, . . ., N det }, all without the use of a basis set.This allows one to use weighted sums of AGPs determinants as trial wave functions, similar to the weighted sum of conventional Fer-miNet determinants seen in Eq. (11).

B. AGPs FermiNet with Unpaired States
To extend the AGPs FermiNet to systems with unpaired states, such as an odd-number of particle system, we use FermiNet geminals and orbitals to replace both the pairing orbitals and the single-particle orbitals in Eq. (16).In this work we consider systems with equal numbers of up-spin and down-spin particles, which are assumed to be fully paired, and systems containing one additional unpaired particle, which may have spin up or spin down.For example, the AGPs FermiNet with an extra spin-up particle is given by where φ k (r ↑ i , r ↓ j ; {r ↑ /i }; {r ↓ /j }) can either be defined as Eq. ( 17) or Eq. ( 19), and ϕ k↑ i (r ↑ j ; {r ↑ /j }; {r ↓ }) is defined in Eq. ( 9).In practice, we generate the required number of FermiNet geminal and orbital parameters in batch at the orbital projection layer.For example, one FermiNet geminal and one FermiNet orbital are generated per determinant for a system with one extra spin-up particle.

V. RESULTS
The power of the AGPs FermiNet Ansatz may be demonstrated by studying the UFG.The Hamiltonian is where is the modified Pöschl-Teller potential, which is widely used in variational and diffusion QMC simulations [16][17][18][19][20][21][22][23].It would be preferable to use a delta function interaction with an infinite s-wave scattering length, but it is difficult to simulate systems with delta-like potentials using QMC methods.Thus, a finite but short-ranged interaction is typically used.The s-wave scattering length of the Pöschl-Teller potential diverges when v 0 = 1.By changing the value of µ at fixed v 0 = 1, it is possible to vary the effective range of the interaction, r e = 2/µ, whilst holding the s-wave scattering length infinite.
We choose to study a system with density parameter r s = 1, where r s , the radius of a sphere that contains one particle on average, provides a convenient measure of the inter-particle distance.Throughout this work, we employ the dimensionless system based on Hartree atomic units: the unit of length is the Bohr radius, a 0 , and the unit of energy is the Hartree.To ensure that the range of the interaction is small compared with the inter-particle separation, we set µ = 12 (r e = 1/6), keeping v 0 = 1 to ensure that the scattering length remains infinite [56].We have also simulated the system with k F = 1 (equivalent to r s = (9π/4) 1/3 ≈ 1.92) and µ = 12 to compare with the fixed-node diffusion Monte Carlo (FN-DMC) result from Forbes et al. [21].
We use both the Slater FermiNet with multiple blockdiagonal determinants and the AGPs FermiNet with multiple geminals to study the unitary Fermi gas from N = 4 to N = 38 particles in a cubic box subject to periodic boundary conditions, as well as AGPs FermiNet on the N = 66 system [57].The same network size, number of determinants, and number of training iterations are used for both Ansätze.The FermiNet orbitals are given by Eq. ( 9) without the bias term.The FermiNet geminal used for systems containing from 4 to 28 particles is the one defined in Eq. ( 17).Unless otherwise stated, all calculations used a linear combination of 32 determinants or 32 geminals.Including contributions from the twoelectron stream improves the optimization rate and can achieve a slightly lower variational energy in larger systems, so Eq. ( 19) was used for systems of N ≥ 29.The inclusion of plane-wave envelopes as defined in Eq. ( 12) also improves the optimization rate.For molecular and electron gas systems, we have found that the bias term in the FermiNet orbital projection (Eq.9) does not affect the accuracy or optimization of the model.We hence set the bias term, g kα i and g k as appropriate, to zero for all calculations presented here.
A comparison of the ground-state energy expectation values given by the two Ansätze is shown in Fig. 1a.The Slater FermiNet, which consists of a linear combination of block-diagonal determinants of FermiNet orbitals, performs well when the number of particles N is smaller than around 10, but the AGPs FermiNet is much superior in larger systems.It is clear that the Slater FermiNet Ansatz has difficulties learning the ground states of large paired systems [58].
In systems containing an odd number of particles, one must be left unpaired.This raises the energy a little and explains the zigzag shape of Fig. 1a.The odd-even staggering is lost for larger systems with the Slater FermiNet Ansatz, indicating the absence of pair formation [16,59].The Slater FermiNet fails to learn the superfluid state.For the AGPs FermiNet, by contrast, the amplitude of the odd-even zigzag remains constant, superposed on the linear increase with N expected of any extensive quantity.
Another comparison between the two Ansätze is shown in Fig. 1b, which depicts the ratio of the interacting and non-interacting energies per particle, known as the Bertsch parameter [4] and defined in Eq. ( 1), as a function of N .All FermiNet energies are variational and the non-interacting energies are exact, so the AGPs Fer-miNet, for which the Bertsch parameter is lower by up to around 30%, is much the better of the two Ansätze.
We next compare our results with the state-of-the-art FN-DMC results of Forbes et al. [21], shown in Fig. 2 for the case k F = 1 and µ = 12.The AGPs FermiNet achieves a lower energy per particle than FN-DMC for all system sizes except for N = 4 and N = 6.The dependence of the Bertsch parameter on system size is also smoother when calculated with the AGPs FermiNet [60].A full training curve of N = 66 with comparison to the FN-DMC energy can be found in Appendix F.
The pairing gap may be found using the approximation formula [15,59] where N is the total number of particles in the box.The results from N = 4 to N = 36 are shown in Fig. 3. Also shown is the thermodynamic (N → ∞) limit of the BCS pairing gap including Gorkov's polarization correction [61]: Here a is the scattering length of the interaction, which is infinite in the UFG.In this limit, ∆ BCS = 1.804EFG and ∆ Gorkov = 0.815E FG , where 2m is the average energy per particle of an unpolarized non-interacting Fermi gas and e is Euler's number [62].The UFG is a strongly coupled system, so the BCS and Gorkov estimates of the gap need not be accurate.
The striking collapse of the pairing gap with increasing system size shows that the Slater FermiNet Ansatz struggles to describe paired states in systems of more than 10 particles.The AGPs FermiNet Ansatz behaves much better, although the oscillations with system size suggest that significant finite-size errors remain even for the largest systems simulated.
Another signature of fermionic superfluidity is the presence of off-diagonal long-ranged order in the twobody density matrix (TBDM), ρ the largest eigenvalue of which diverges as the number of particles N tends to infinity [63].The superfluid condensate fraction c may be obtained by evaluating [64] where Ω is the volume of the simulation cell, N ↑ is the number of spin-up particles, and ρ (2)T R ↑↓ (r) is the rotational and translational average of the TBDM The one-body density matrix, by contrast, tends to zero in the r → ∞ limit [63].A full discussion of the methods used to evaluate the condensate fraction in QMC simulations can be found in the Appendix and the CASINO manual [64].
After fully training both the Slater FermiNet and the AGPs FermiNet for the N = 38 particle system, we used the resulting neural wave functions to compute the quantity Ω 2 N ↑ ρ (2)T R ↑↓ (r).The results are shown in Fig. 4, which provide further evidence that the Slater FermiNet fails to converge to the superfluid ground state; the quantity (r) appears to be approaching zero in the large pair-separation limit, implying that the condensate fraction is also zero.The same quantity for the AGPs Fer-miNet approaches a finite value which we estimated to be roughly c = 0.44(1) using the eight data points with separations r/r s ≥ 2.0.This value is consistent with previous estimations from experiments and the most recent AFMC value from [22] (Table I).
In addition, we also estimated the condensate fractions for the N = 66 UFG at a fixed µ = 12 with two different densities: r s = 1 (k F r e = 0.32) and k F = 1 (k F r e = 0.17), respectively.We compute the quantity Ω 2 N ↑ ρ (2)T R ↑↓ (r) at five sequentially-spaced separations r near r = L/2, where the quantity has approached its asymptotic value.We then take the average of the five data points to get estimated values of the condensate fraction.[65].Our estimate of the condensate fraction for N = 66 is c = 0.42(1) at k F r e = 0.32, and c = 0.52(1) at k F r e = 0.17, which are both consistent with the experiments.The results are summarized in Table I.FIG.1: Comparison between results obtained using the AGPs FermiNet and the Slater FermiNet for different numbers of particles, N , with r s = 1 and µ = 12.All simulations used 32 determinants, 300,000 optimization steps, and the same hyperparameters, which are detailed in the Appendix.
Although VMC methods are generally considered to be less accurate than FN-DMC methods, an important advantage of VMC methods is that almost any expectation value, including any reduced density matrix, may be estimated without bias.The same is not true of FN-DMC simulations, which sample the wave function instead of its square modulus and produce biased "one-sided" estimates of the expectation values of operators that do not commute with the Hamiltonian [14].Thus, there are very few unbiased and accurate first-principles calculations of the condensate fraction.Our approach, having both the advantages of VMC and surpassing the accuracy of DMC, provides solutions to these problems and a more accurate way to estimate general expectation values.
Finally, we study how the number of block-diagonal determinants required to achieve a given accuracy scales with the number of particles in the system.We choose six even-particle systems from N = 4 to N = 14 and compare the energies obtained using linear combinations of multiple block-diagonal FermiNet Slater determinants against energies obtained using a single-determinant (and thus single-geminal) FermiNet AGPs wave function.All other hyperparameters are as given in Table II (Appendix B 1).In Fig. 5, we show that the number of blockdiagonal determinants required to achieve a given percentage accuracy increases approximately exponentially with the number of particles.Plots for each individual system, along with a more detailed discussion, can be found in Appendix C.These results suggest that multideterminant Slater FermiNet wave functions constructed using a neural network of fixed size are incapable of describing the ground state of the UFG accurately unless  the number of block-diagonal determinants rises exponentially with system size.Hence, in practice, the AGPs FermiNet is required for studying paired systems.

VI. DISCUSSION
In this work, we used neural wave functions to study the benchmark superfluid system known as the UFG [70].We showed that the Slater FermiNet Ansatz has difficulties in describing paired systems with strong, shortranged attractive interactions between particles of opposite spin.Hence, we proposed a way to improve the variational Ansatz by using determinants of FermiNet geminals, similar to an AGPs or a BCS wave function.We showed mathematically that the Slater FermiNet is a limiting case of the AGPs FermiNet despite the use of fewer parameters in the latter.It follows that any Fer-miNet wave function can in principle be written as an AGPs FermiNet wave function.
We compared the total energies and energies per particle of the UFG as calculated using the Slater FermiNet and the AGPs FermiNet.The former fails to produce a paired state when the number of particles, N , is greater than around 10, while the AGPs FermiNet works very well.
As the UFG has a superfluid ground-state, we computed the pairing gap and condensate fraction for the N = 38 system and compared estimates made with the Slater FermiNet and the AGPs FermiNet.There is a clear qualitative difference between the pairing gap obtained using the AGPs FermiNet and the Slater FermiNet, with the latter approaching zero as the number of particles N increases.Calculations of the superfluid condensate fraction show a similar behavior: the AGPs FermiNet gives an accurate finite result, while the value obtained using the Slater FermiNet tends to zero in the limit of large system size.Although the AGPs pairing gap shows significant finite size errors, it lies close to the mean-field BCS result with Gorkov-Melik-Barkhudarov corrections [61].Taken together, these results show that the Slater FermiNet is unable to represent large systems with superfluid ground states.The AGPs FermiNet is much more suitable for studying paired systems such as the UFG.
To demonstrate the success of the AGPs FermiNet, we also compared our calculated total energies with state-ofthe-art fixed-node diffusion QMC energies obtained using a Jastrow-BCS Ansatz [21].For all systems with more than a few particles, the AGPs FermiNet achieves lower (i.e., better) variational energies than FN-DMC using the same model interaction and system parameters.
The inability of the Slater FermiNet Ansatz to accurately describe the UFG ground state comes as a surprise because the original FermiNet paper [1] argued that any many-body fermionic wave function can be represented as a single determinant of FermiNet orbitals.However, the mathematical argument relies on the construction of FermiNet orbitals with unphysical discontinuities.Whether or not any wave function can be represented as a single determinant of FermiNet orbitals of the type used in practice, which are differentiable everywhere except at electron-electron and electron-nuclear coalescence points, remains an open question.Another limitation is that the architecture of the FermiNet neural network, which is rather simple, may not be able to represent an arbitrary many-electron FermiNet orbital.Even if a single-determinant Slater FermiNet wave function is general in principle, there is no guarantee that it is equally easy to represent all wave functions.It may be that producing an accurate representation of a paired wave function requires the width and number of layers in the neural network to increase rapidly with system size.Furthermore, if a network of fixed size is used, it may be necessary to increase the number of Slater FermiNet determinants rapidly as the system size increases.The observation that the Slater FermiNet works well when N ⪅ 10 but that the quality of the results degrades rapidly for larger systems, along with the scaling study presented in the final part of the Results section, suggests that this is, in fact, the case.
The AGPs FermiNet introduced in this paper shares many of the strengths of the Slater FermiNet.In particular, there is no need to construct and optimize a new basis set for every new system or particle type.If the AGPs FermiNet proves equally successful in other paired systems, it may now be relatively easy to investigate the importance of pairing in molecules, electron-positron systems, electron-hole liquids, and other s-wave superfluids.Another strength of the AGPs FermiNet is the ease with which it is possible to optimize linear combinations of determinants of FermiNet pairing orbitals, such as the one in Eq. ( 31).This is much more difficult to accomplish with conventional wave functions based on explicit twoelectron pairing orbitals or pairing orbitals represented as outer products of single-particle orbitals or basis functions.Just as the many-particle orbitals in a Slater FermiNet radically generalize single-particle orbitals by incorporating electron-electron terms in a permutationequivariant fashion, so the pairing functions in an AGPs FermiNet generalize BCS-style pairing functions by incorporating the effects of the remaining electrons in a permutation-equivariant fashion.
The AGPs FermiNet introduced here has a straightforward Pfaffian extension and can thus be applied to nons-wave and triplet pairing.Therefore, we expect it to become a powerful tool for understanding strongly correlated non-s-wave superfluid and superconducting systems such as Helium-3 or high-T c and p-wave superconductors.Finally, our approach is not limited to the FermiNet neural network and can be readily adapted to use more recent architectures such as the Psiformer [71], GLOBE and MOON [72], and DeepErwin [73].this project by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC).JK is part of the Munich Quantum Valley, which is supported by the Bavarian state government with funds from the Hightech Agenda Bayern Plus.WTL is supported by an Imperial College President's PhD Scholarship; HS is supported by the Aker Scholarship; and GC is supported by the UK Engineering and Physical Sciences Research Council (EP/T51780X/1).We also acknowledge the support of the Imperial-TUM flagship partnership.The two-body density matrix (TBDM) in first quantized notation can be written as where α and β denote the spin or particle species.The superfluid condensate fraction in a finite and periodic system is defined as where Ω is the volume of the simulation cell, N α is the number of particles with spin α, and ρ (2)T R αβ (r) is the translational and rotational average of the TBDM given in Eq. (39).
The one-body density matrix (OBDM) is expected to tend to zero as r → ∞.However, because of finite-size effects, the OBDM is not necessarily zero within our simulation cell.We therefore use an improved estimator in Eq. (B3) that removes the one-body contribution explicitly [64]:  This quantity can then be estimated using Monte Carlo sampling.
Appendix C: How many block-diagonal determinants does the Slater FermiNet need to achieve the same accuracy as the AGPs FermiNet with one determinant?
In the main text, we have demonstrated that the Slater FermiNet with 32 block-diagonal determinants is able to capture superfluidity in small systems but fails at larger systems.Therefore, we are interested in the scaling of the original block-diagonal determinant FermiNet wave function with respect to the system size, i.e. how many block-diagonal determinants do we need in order for the Slater FermiNet to converge to the ground state at each system size?
To answer this question, first we set the AGPs Fer-miNet with one determinant energies as baselines and plot the percentage difference of the Slater FermiNet from the baseline against the number of block-diagonal determinants used in the Slater FermiNet wave function, repeated at different system sizes from 4 to 14 even particles systems.The results are shown in Fig. 6.
As the results suggested, it becomes more difficult for the Slater FermiNet to get close to the AGPs FermiNet baseline as the number of determinant increases, especially for larger systems.This is due to the limited performance of the optimizer as the number of determinant increases.Hence, due to the constraints on time and resources, it is not feasible to continue to increase the number of determinants until the Slater FermiNet achieve the same accuracy as the AGPs FermiNet.Instead, we de-cided to set two thresholds for the percentage difference between the two results.Here, we have used 5% and 10% as the thresholds, plotted as the two horizontal lines in Fig. 6.By plotting the x-intercept of the curves with the two threshold lines against system size, we can determine the relationship between the two, as shown in Fig. 5.As the y-axis in Fig. 5 is in logarithmic scale, a roughly linear relationship suggests an exponential scaling of the number of block-diagonal determinants as system size increases.
The result indicates that, in theory, the Slater Fer-miNet is capable of converging to the ground state given the number of block-diagonal determinants is sufficiently large.In practice, the number of block-diagonal determinants required to learn the ground state increase exponentially as system size gets bigger and the number get inaccessible rapidly.Although the use of dense determinants in the UFG does provide lower energies comparing to the blockdiagonal determinants, they are still significantly higher than the AGPs FermiNet energies.In addition, the qualitative behaviors, such as the absence of odd-even staggering, are still similar to the block-diagonal FermiNet, which are qualitatively different from the AGPs Fer-miNet.
The total energy of the UFG simulation cell, measured in units of the free Fermi gas energy EFG.The Slater FermiNet Ansatz begins to fail when N ⪆ 10.
The Bertsch parameter ξ (the ratio of the interacting and non-interacting ground-state energies per particle) as a function of the number of particles N .

FIG. 2 :
FIG. 2: Comparison of the system-size dependent values of the Bertsch parameter, ξ, as calculated using the AGPs FermiNet and FN-DMC, with k F = 1 and µ = 12.According to the variational principle, lower values are better.The error bars on the AGPs FermiNet results are smaller than the sizes of the crosses.Inset: difference between the AGPs and FN-DMC values of the Bertsch parameter.The errors in the inset are obtained adding the standard errors of and AGPs FermiNet results in quadrature.The latter are obtained by computing the standard error of the MCMC-averaged Bertsch parameter accumulated over 50,000 inference steps.

FIG. 3 :
FIG. 3: Pairing gaps calculated with the Slater FermiNet and the AGPs FermiNet for different numbers of particles N , with r s = 1 and µ = 12.

FIG. 4 :FIG. 5 :
FIG. 4: Comparison of the TBDM estimators calculated using the AGPs FermiNet and the Slater FermiNet with N = 38, r s = 1 and µ = 12.The error bars show the standard error of the TBDM estimator, accumulated over 2,000 inference steps.Most of the error bars are so small that they are obscured the symbols. ρ

FIG. 8 :
FIG.8: Training curve of the 66 particle UFG with 32 determinants and using the geminals in Eq. 19.Red dashed line is the FN-DMC result from[21] under the same density and interaction width.

TABLE II :
Hyperparameters used in all simulations.

TABLE III :
Network sizes and number of determinants used in all simulations.The corresponding mathematical symbols mentioned in the main text of the paper, where available, are also listed.
2. Estimation of the Two-Body Density Matrix Percentage difference between the Slater FermiNet with various number of block-diagonal determinants and the AGPs FermiNet with one determinant at different system sizes, from 4 to 14 particles.

TABLE IV :
Total energy of the UFG with 36 to 38 particles and their corresponding pairing gap using different wave functions.