# A Unifying View of Fermionic Neural Network Quantum States: From Neural
Network Backflow to Hidden Fermion Determinant States

Zejun Liu  The Anthony J. Leggett Institute for Condensed Matter Theory and
IQUIST and NCSA Center for Artificial Intelligence Innovation and Department
of Physics, University of Illinois at Urbana-Champaign, IL 61801, USA  Bryan
K. Clark<EMAIL_ADDRESS>The Anthony J. Leggett Institute for Condensed
Matter Theory and IQUIST and NCSA Center for Artificial Intelligence
Innovation and Department of Physics, University of Illinois at Urbana-
Champaign, IL 61801, USA

###### Abstract

Among the variational wave functions for Fermionic Hamiltonians, neural
network backflow (NNBF) and hidden fermion determinant states (HFDS) are two
prominent classes to provide accurate approximations to the ground state. Here
we develop a unifying view of fermionic neural quantum states casting them all
in the framework of NNBF. NNBF wave-functions have configuration-dependent
single-particle orbitals (SPO) which are parameterized by a neural network. We
show that HFDS with $r$ hidden fermions can be written as a NNBF with an
$r\times r$ determinant Jastrow and a restricted low-rank $r$ additive
correction to the SPO. Furthermore, we show that in NNBF wave-functions, such
determinant Jastrow’s can generically be removed at the cost of further
complicating the additive SPO correction increasing its rank by $r$. We
numerically and analytically compare additive SPO corrections generated by the
product of two matrices with inner dimension $r$. We find that larger $r$
wave-functions span a larger space and give evidence that simpler and more
direct updates to the SPO’s tend to be more expressive and better
energetically. These suggest the standard NNBF approach is preferred amongst
other related choices. Finally, we uncover that the row-selection used to
select single-particle orbitals allows significant sign and amplitude
modulation between nearby configurations and is partially responsible for the
quality of NNBF and HFDS wave-functions.

##  I Introduction

Simulating quantum many body systems is difficult often requiring various
approximations to make progress. Since the early history of quantum mechanics,
variational approaches have been one of these core approximations. In the
variational approach, one starts with a parameterized class of wave-functions
and searches amongst this class for the lowest energy state. Hartree Fock, a
variational search over the class of non-interacting wave-functions was
amongst the earliest variational approaches for simulating fermions, and works
by optimizing the parameterized single particle orbitals (SPO) which make up a
Slater Determinant [1]. Even today, many of the state-of-the-art approaches
for fermions and frustrated magnets build on top of this non-interacting class
dressing and complementing it in various ways. For example, a standard class
of Fermionic variational wave-functions is the Slater-Jastrow form,
$\Psi(c)=J(c)\Psi_{\text{SD}}(c)$ , where $\Psi_{\text{SD}}(c)$ is a non-
interacting Slater Determinant and $J(c)$ is a multiplicative Jastrow factor
which adds significant correlation on top of the Slater Determinant [2, 3, 4].
The Jastrow factor $J(c)$ has historically been a strictly positive function
taking the form of an exponential of one and two-body operators
$J(c)=\exp[-u_{1}(c)+u_{2}(c)]$. More recently, with advances in optimization
as well as the progress in machine learning (ML) architectures which
approximate a wide class of functions [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16, 17, 18], Jastrows have become more sophisticated, parameterized, and often
negative spanning wave-function classes such as correlated product states,
RBM, etc [19, 20, 21, 13, 16, 17, 18].

In addition, to further improving the Jastrow factor, recent work has applied
ML architectures to dress the Slater Determinant part of the Slater-Jastrow
wave-function. Two such paradigms include the Neural Network Backflow (NNBF)
[22, 11, 12, 23, 13] and the hidden Fermion determinant state (HFDS) [7, 8].
NNBF replaces the static set of SPO with a configuration-dependent set which
are generated by a neural network. HFDS instead works in the paradigm of
projected hidden fermions using neural networks to replace the standard Slater
Determinant with a larger determinant which includes SPO’s from an additional
projected $r$ hidden fermions.

Figure 1: (a-d) Visualization of wave-functions. Red area are configuration-
dependent matrices output from the NN; other colored area depend on those
outputs; black area are optimized but configuration-independent and white area
is zero. Same color in a given line are identical (up to a trivial added
identity). Brackets with $c$ in the lower-right indicate row-selection of that
matrix with configuration $c$. All wave-functions (including HFDS) can be
written in the style of NNBF as an additive correction to the SPO. Additive
correction to lower wave-functions are less constrained and this additional
freedom is often reflected in better energetics. Note that (c) is formally a
subset of (d) even when $r>L_{s}$. (b) is a subset of (c) when the NN for (c)
has two additional layers and otherwise we see that both (b) and (c) and (a)
and (b) can be optimized towards each other with low infidelity. The
determinant Jastrow in the l.h.s of (b) can be replaced by an exponential
Jastrow with energetics of similar quality. (e) Wave function architecture for
various NNBF: Fermion configuration $c$ as a set of binary numbers is input to
a feed-forward neural network with two hidden layers. The output from the NN
gives the $r\times N_{s}$-size $C_{{\bm{\theta}}}(c)$ and $r\times r$-size
$F_{{\bm{\theta}}}(c)$ backflow matrices, which are then combined with the
static matrices $A$ and $B$ to form SPO for various NNBF wave functions. The
amplitudes from the wave functions are used to sample new fermion
configurations with Monte Carlo.

These seemingly distinct wave-functions both achieve similar significant
improvements to the variational power for fermionic systems. One cannot help
but inquire about the existence of some underlying mechanism that links these
states. In this work, we investigate explicitly the connection between NNBF
and HFDS. Surprisingly, we will find that that HFDS can be cast exactly as a
Jastrow-NNBF with a specific determinant Jastrow and a restricted “low-rank”
SPO correction. This motivates us then to construct a series of Jastrow-NNBF
wave-functions which differ in their Jastrow’s and the way the neural network
output is used to modify their single particle orbitals.

From both analytical and numerical arguments, we identify the relation between
these different wave-functions (see Fig. 1). We consider the effect of the
multiplicative determinant Jastrow on these wave-functions. We find that
energetically, it is similar in quality to exponential Jastrow. Moreover,
instead of explicitly multiplying a NNBF-type wave-function by a determinant
Jastrow, we can absorb the determinant Jastrow into a slightly more
complicated additive correction to the single particle orbitals. We also
compare various ways to additively modify the single particle orbitals.
Generically, we find that simpler, more general and more direct updates to the
configuration-dependent single particle orbitals tend to be more expressive
giving evidence that the type of NNBF induced by HFDS is less expressive than
the standard NNBF when the number hidden fermions is less than the physical
fermion number. Finally, we also give a qualitative understanding for why NNBF
(also HFDS) differs from other wave-functions by rewriting these wave-
functions as a fixed Slater determinant times an effective Jastrow and
quantitatively showing that the effective Jastrow with configuration-dependent
row-selection endows the wave-function ansatz with more flexibility compared
with other Jastrow factors such as determinant Jastrow.

##  II Neural Network Backflow

Both the NNBF and the HFDS build on top of Slater determinants. A Slater
determinant is represented by a set of $N_{s}$ single particle orbitals (where
$N_{s}$ is the total number of fermions) each of size $L_{s}=(2s+1)L$ where
$L$ is the number of lattice sites and $s$ is the spin of fermions.

In this work, we consider fermions on a discrete lattice, and represent the
SPO’s as an $L_{s}\times N_{s}$ matrix $A$. Let $A[c]$ be the $N_{s}\times
N_{s}$ matrix generated by taking the rows of $A$ which correspond to the
locations of fermions in a configuration $c$. The amplitude of a Slater
determinant wave-function for $c$ is then $\Psi_{\text{SD}}(c)=\det A[c]$.

In NNBF, we instead make the SPO configuration-dependent,

$\Psi(c)=J(c)\det(\widetilde{A}_{\bm{\theta}}(c)[c])$ (1)

where the subscript ${\bm{\theta}}$ are parameters to the neural network (NN)
which takes as input a configuration $c$ and outputs a matrix (or a series of
matrices) which is used to construct $\widetilde{A}_{\bm{\theta}}(c)$. Take
$J(c)=1$ when it is not specified.

The NN needs to be permutation invariant with respect to the input; in
standard NNBF, this network is a feed forward neural network (FFNN) which
takes the input $c$ as an $L_{s}$ bit binary number with 1’s at the sites the
fermions occupy. To make comparisons fair, throughout this paper we fix all
the layers (number of layers and layer widths $w$) except for a final linear
layer whose output changes to be the full set of relevant matrices (In Sec.
IV.2, we will discuss the possibility of reconstructing the SPO for one wave-
function with another by adding a few more layers, otherwise, we assume the
NN’s have identical architecture).

There are various ways the new SPO can be generated from the output of the NN.
For example, through an additive correction

$\widetilde{A}_{\theta}(c)=A+M_{\bm{\theta}}(c)$ (2)

where $M_{\bm{\theta}}(c)$ is an $L_{s}\times N_{s}$ matrix. Alternatively, a
left rotation, $\widetilde{A}_{\bm{\theta}}(c)=L_{\bm{\theta}}(c)A$ (shown in
Ref. [11] to be nearly identical to an additive correction) or a right
rotation $\widetilde{A}_{\bm{\theta}}(c)=A\widetilde{F}_{\bm{\theta}}(c)$ can
be used to generate the new SPO. Left and right rotations are qualitatively
distinct with left rotations mixing sites within each individual single
particle orbital and right rotations mixing SPO between each other at fixed
site. We can even compose changes to the SPO doing first an additive
correction and then a right rotation,

$\widetilde{A}_{\bm{\theta}}(c)=(A+M_{\bm{\theta}}(c))\widetilde{F}_{\bm{\theta}}(c)$
(3)

Interestingly, right rotations can be factored out of the Slater determinant
and treated as a Jastrow factor of the form

$J_{\det}(c)=\det\widetilde{F}_{\bm{\theta}}(c).$ (4)

The converse of this is also true: any determinant Jastrow factor $\det
F_{\bm{\theta}}(c)$ for any $F_{\bm{\theta}}(c)$ whose size $r\leq N_{s}$ can
be absorbed into Eqn. (3) by letting $\widetilde{F}$ be $F$ padded with the
identity up to size $N_{s}$. This will lead us to interchangably treat a
determinant Jastrow and a final right rotation of the SPO as equivalent
(paying attention when necessary that $r\leq N_{s}$) throughout. For example,
the wave-function of eqn. (3) will be JNNBF,

$\Psi_{\textrm{TNNBF}}\underset{\mathrm{r\leq
N_{s}}}{=}\Psi_{\textrm{JNNBF}}\equiv J_{\det}\Psi_{\textrm{NNBF}}$ (5)

##  III Hidden Fermion Determinant States

In HFDS, starting with an $L_{s}\times N_{s}$ set of SPO’s $A$, the HFDS
supplements the physical system with the introduction of $r$ hidden fermions
and consequently an additional $r$ single particle orbitals as well as some
number of new sites. The $r$ additional single particle orbitals on the
original ‘real’ sites are represented by a static $L_{s}\times r$ matrix $B$
and $B[c]$ represents the rows selected from $B$ corresponding to the real
Fermion configuration.

The evaluation of both the ‘real’ and ’hidden’ single particle orbitals on the
extra sites is done implicitly. The NN outputs the value of the original
$N_{s}$ and extra $r$ SPO’s with respect to the (never explicitly specified)
locations of the $r$ hidden fermions on the hidden sites. This gives the
$r\times N_{s}$ matrix $C_{\bm{\theta}}(c)$ for the original SPO’s and the
$r\times r$ matrix $F_{\bm{\theta}}(c)$ for the new SPO’s.

This will then give an amplitude, for $r$ hidden fermions of

$\begin{split}\Psi_{r\text{-HF}}(c)&=\det\begin{bmatrix}A[c]&B[c]\\\
C_{\bm{\theta}}(c)&F_{\bm{\theta}}(c)\end{bmatrix}\\\ \end{split}$ (6)

In the case of HFDS (Eqn. (6)), the NN parameters ${\bm{\theta}}$ and $A$ and
$B$ are all optimized when applied to the ground state approximation for
Fermionic Hamiltonian.

##  IV Comparison between Wave Functions

We start by rewriting HFDS with $r$ hidden fermions in a somewhat different
form. Without loss of generality, we assume that $F_{\bm{\theta}}(c)$ is
invertible for all $c$ (see Appendix A). Then we can rewrite Eqn (6) as

$\Psi_{r\text{-HF}}(c)=\det(F_{\bm{\theta}}(c))\det\left((A-BF_{{\bm{\theta}}}(c)^{-1}C_{\bm{\theta}}(c))[c]\right)$
(7)

In the parlance of NNBF, the ${r\text{-HF}}$ wave-function is neural network
backflow with an additive correction of the form

$\widetilde{A}_{\bm{\theta}}(c)=A+B(F_{\bm{\theta}}(c))^{-1}C_{\bm{\theta}}(c)$
(8)

multiplied by a determinant Jastrow (alternatively an additive correction
composed with a right-rotation with $F$).

This formula can be derived either directly by working with determinants of
block matrices; or alternatively by using the matrix determinant lemma and
computing the change in the determinant between the (padded with identity)
matrix $A[c]$ and the matrix $(A+BC_{\bm{\theta}}(c))[c]$.

${r\text{-HF}}$ then is a form of NNBF but does differ from the most standard
Jastrow-NNBF in two ways. First, there is a somewhat non-standard Jastrow
factor consisting of the determinant of a $r\times r$ matrix. We will want to
address two questions about this Jastrow: (a) Is it preferable to use a
determinant Jastrow compared to a more basic (possibly also NN-inspired)
Jastrow and (b) Does the addition of a determinant Jastrow on top of NNBF span
a larger class of wave-functions than no Jastrow at all. Second, we will be
interested in understanding the difference in the additive corrections between
${r\text{-HF}}$ and standard NNBF. Allowing the NNBF to use as its Jastrow an
$r\times r$ determinant and fixing the NN’s to the same depth and width, we
want to ascertain whether using the NN output for matrices $C$ and $F$ which
are multiplied as in Eqn. (8) to generate the additive SPO is better then just
outputting $M$ directly.

###  IV.1 Comparisons of Additive Corrections

We will address these in reverse order starting first with comparing the
additive correction between ${r\text{-HF}}$ and JNNBF.

To make progress on this, it will also be useful to define a closely related
wave-function

$\widetilde{A}_{\bm{\theta}}(c)=A+BC_{\bm{\theta}}(c)$ (9)

as well as its determinant-Jastrow multiplied form ${r\text{-JBF}}$
(equivalently r-TBF). ${r\text{-JBF}}$ differs from ${r\text{-HF}}$ in that it
directly uses the outputted matrix $C_{\theta}(c)$ whereas ${r\text{-HF}}$
uses $(F_{\theta}(c)^{-1})C_{\theta}(c)$.

If we consider Eqn. (8) or Eqn. (9), we see that the additive correction to
the SPO is the product of an $L_{s}\times r$ matrix ($B$) and an $r\times
N_{s}$ matrix, ($F^{-1}_{{\bm{\theta}}}C_{\bm{\theta}}$ or $C_{\bm{\theta}}$),
resulting in a low-rank additive correction of, at most, rank $r$. As
generically the rank of the additive correction to JNNBF will be full rank
(i.e. $N_{s}$), almost none of the JNNBF SPO’s will be representable by either
Eqn. (8) or Eqn. (9) when $r<N_{s}$. On the other hand, we can show every
additive SPO generated by ${r\text{-JBF}}$ can also be generated by JNNBF.
This is done by converting the ${r\text{-JBF}}$ neural network (even for
$r>L_{s}$) into a JNNBF neural network of the same size by first increasing
the ${r\text{-JBF}}$ neural network by one linear layer (encoding $B$ in the
weights) so that it outputs $M_{\bm{\theta}}(c)$ and then compressing the two
final linear layers into a single linear layer; This results in a JNNBF wave-
function with the same initial layers. Notice, it is also straightforward to
see that $L_{s}$-JBF and JNNBF span the same space (let $B$ be the identity)
and that $(r+1)$-JBF contains $r$-JBF (but the converse is generically not
true when $r<N_{s}$). We also see evidence of this numerically as the energy
of ${r\text{-HF}}$ (or ${r\text{-JBF}}$) decrease monotonically with $r$ out
to, at least, $r=L_{s}$ with the most significant decreases at $r<N_{s}$. See
Appendix. B for more explicit arguments for these descriptions.

Analytically proving whether ${r\text{-HF}}$ is strictly included in JNNBF is
made difficult by the $F^{-1}$ term. That said, JNNBF is almost certainly
larger than ${r\text{-HF}}$ (for $r\leq L_{s}$) as we expect on general
grounds that ${r\text{-JBF}}$ and ${r\text{-HF}}$ span a very similar (if not
identical) space. In fact, in the $w\rightarrow\infty$ limit they do span the
same set of configuration-dependent SPO’s under the assumption of universal
approximation theorem for arbitrary-width NN [24]. At finite $w$, we have
numerically considered low-energy (with respect to the $4\times 4$ Hubbard
model) ${r\text{-JBF}}$ and ${r\text{-HF}}$ states and found we can optimize
them towards each other (in both directions) so that their infidelity is lower
then $10^{-2}$. This is consistent also with finding their optimized energies
being very close with ${r\text{-JBF}}$ being slightly lower other than at
$r=2$. Unsurprisingly, we find that JNNBF finds lower-energy variational
states than both ${r\text{-JBF}}$ and ${r\text{-HF}}$ at all $r$.

This leaves us with the following situation:

$\displaystyle{r\text{-HF}}$ $\displaystyle\approx{r\text{-JBF}}$ (10)
$\displaystyle{r\text{-JBF}}$
$\displaystyle\underset{\text{spo}}{\subset}(r+1)\text{-JBF}\underset{\text{spo}}{\subset}L_{s}\text{-JBF}\underset{\text{spo}}{=}\text{JNNBF}$
(11) $\displaystyle r\text{-HF}$
$\displaystyle\underset{\text{spo}}{\subset}\text{JNNBF}:w\rightarrow\infty$
(12)

where we use the notation $\underset{\text{spo}}{\subset}$ to indicate one
class of wave-functions contains the same set of SPO’s as the other class (a
strictly stronger inclusion than the wave-functions being the same). Notice
that Eqn. (11) is also true if we remove the Jastrow from each term.

###  IV.2 Benefit of Jastrow on Neural Network Backflow

In the previous subsection, we argued that when dressed with a determinant
Jastrow, the additive corrections being used by JNNBF were generically a
superset of the additive correction of other wave-functions such as
${r\text{-HF}}$ and ${r\text{-JBF}}$ whose additive corrections had restricted
rank. In this subsection, we will consider whether the determinant Jastrow is
important at all - i.e. is NNBF identical, in some limits, to JNNBF.

To make progress on this, we will actually start by focusing on comparing
${r\text{-JBF}}$ (or ${r\text{-HF}}$) with ${r\text{-BF}}$. To make this
comparison, we will take the ${r\text{-JBF}}$ (${r\text{-HF}}$) wave-functions
and rewrite them as a BF wave-function with a different additive correction.
-i.e. we will have the Jastrow term absorbed into the additive correction.

When $r\leq N_{s}$, we can rewrite (see Appendix C for details) these wave-
function as

$\Psi_{r\text{-HF}/\text{JBF}}=\det\left((A-B^{\prime}C_{\bm{\theta}}^{\prime}(c))[c]\right)$
(13)

where

$B^{\prime}=\begin{bmatrix}-A_{1}&B\end{bmatrix},$ (14)

and

$C_{\bm{\theta}}^{\prime}(c)=\begin{cases}\begin{bmatrix}F-I&\bm{0}\\\
F^{-1}C_{1}F&F^{-1}C_{2}\end{bmatrix},&\text{for }r\text{-HF}\\\
\begin{bmatrix}F-I&\bm{0}\\\ C_{1}F&C_{2}\end{bmatrix},&\text{for
}r\text{-JBF}\end{cases}$ (15)

are of size $L_{s}\times 2r$ and $2r\times N_{s}$, respectively. Notice that,
this means each ${r\text{-JBF}}$ (${r\text{-HF}}$) wave-function can be
written as a slightly more complicated additive correction of rank $2r$. As
both $2r$-BF and ${r\text{-JBF}}$ can be written as a neural network backflow
with rank-$2r$ additive corrections, this suggests that they might span a
similar space of wave-functions. Generically, we might expect that $2r$-BF
actually spans a larger space because the $2r\times L_{s}$ SPO’s are
completely general for $2r$-BF while for ${r\text{-JBF}}$ (${r\text{-HF}}$)
the upper right $r\times(N_{s}-r)$ corner of the additive correction is forced
to be zero. In the $w\rightarrow\infty$ limit it is strictly true that $2r$-BF
contains the (determinant absorbed) SPO’s of both ${r\text{-JBF}}$ and
${r\text{-HF}}$ (but not the converse).

Unsurprisingly, we find that energetics of $2r$-BF and ${r\text{-JBF}}$ are
similar on both $4\times 8$ and $4\times 4$ Hubbard models (see Fig. 3) with
$2r$-BF being slightly lower in energy.

Moreover, through numerical fidelity matching of the low energy states (see
Fig. 4), we find that $2r$-BF can represent both ${r\text{-JBF}}$ and
${r\text{-HF}}$ with infidelities less then $10^{-2}$ for all $r$. In the case
of ${r\text{-JBF}}$, essentially this tells us that $2r$-BF can output the
matrix $C_{1}F$ from ${r\text{-JBF}}$ where previously it was only outputting
$F$ and $C$ separately. The converse (${r\text{-JBF}}$ and ${r\text{-HF}}$
matching $2r$-BF) also achieves infidelities less then $10^{-2}$ except at
$r=1$ where the fidelity matching is worse in that direction.

We even attempt to go beyond just fidelity matching and try to directly
optimize the SPO distance between $2r$-BF and ${r\text{-JBF}}$ by optimizing
the $2r$-BF wave-function (the opposite is not possible because of the zeros),
see Fig. 6. In some limit, optimizing the SPO distance and the fidelity are
closely related (i.e. their respective infideltiies become zero at the same
time), but they do probe the closeness of wave-functions in different ways.
For example, it is possible to have low SPO-distance but not super-high
fidelity because even if the SPO distance is small, its contribution to the
fidelity can be amplified by large values generated from the inverse of the
row-selected SPO (see Appendix. D for a detailed discussion of this). We find
that the $2r$-BF can match the SPO’s quite successfully at all $r$’s for
${r\text{-JBF}}$ and at $r\geq 4$ for ${r\text{-HF}}$.

While the above numerical experiments are performed by fixing the NN for all
the wave functions, we can also ask whether there are ‘minor’ architectural
changes to the FFNN for $2r$-BF which would allow for an exact reconstruction
of other wave-functions. We demonstrate an example of this reconstruction in
Appendix. E where we write a ${r\text{-JBF}}$ wave-function as a
${r\text{-BF}}$ wave-function with two additional layers (with different
activation functions) and widths expanded only by an extra factor of $r$.

It is worth mentioning that it is also possible to reconstruct ${r\text{-HF}}$
from $2r$-BF by adding a few more layers to the FFNN, where we can use
adjugate formula to obtain the inverse of $F$, i.e.,
$F^{-1}=\text{adj}(F)/\det F$. However, this involves the evaluation of
determinants within the NN, and it will require widths on the order of $r!$
number of neurons if implemented directly via the definition of determinant
making this direct approach impractical (we know that the evaluation of
determinant can be done in $O(n^{3})$ steps, but this will require many more
layers within the NN). Nevertheless, the fidelity and the energetics from Fig.
2 indicate that even without any change in architecture, ${r\text{-JBF}}$ or
NNBF are already as expressive as ${r\text{-HF}}$ and effective in
approximating the ground state without the necessity of $F^{-1}$ to be
included in the additive correction.

We now directly consider the relation between NNBF and
${r\text{-JBF}}/{r\text{-HF}}$. On the one hand, when $2r<L_{s}$, almost none
of the NNBF SPO’s will be representable by either of these two wave-functions.
On the other hand, we already know that all the $2r$-BF wave-functions are
representable by NNBF and given the evidence that $2r$-BF is a superset of
${r\text{-JBF}}$, this suggests that all ($r\leq N_{s})$ ${r\text{-JBF}}$ have
(determinant absorbed) SPO’s that are representable by NNBF; this is
rigorously true for $w\rightarrow\infty$. This argument also suggests that
JNNBF and NNBF should span the same space for all $r\leq N_{s}$. Upon
optimization (see Fig. 2), we find that their energetics are extremely close
with very minimal $r$ dependence in JNNBF out to $r=N_{s}$ (surprisingly
actually even for $r>N_{s}$ there is still very little $r$ dependence). We
again attempt to match the fidelity of $r$-JNNBF with NNBF; This optimization
is done through $L_{s}$-BF which is equivalent to NNBF, and we find an
infidelity of essentially $10^{-5}$ for all $r$, suggesting that NNBF really
does contain JNNBF (see Fig. 5).

The conclusion of this section is that the addition of a Jastrow seems to not
increase the space of wave-functions spanned by NNBF , and increases the space
of wave-functions spanned by ${r\text{-BF}}$ up to $2r$-BF when $r\leq N_{s}$.
When the number of hidden fermions is very large - more then the number of
actual fermions, it is less clear whether the Jastrow should have additional
benefit although in practice we see almost no decrease in energy for JNNBF at
$r>N_{s}$ (See Fig. 2).

###  IV.3 Exponential vs Determinant Jastrows

What we have seen in the previous section, is that the Jastrow matters very
little for JNNBF, but does increase the ‘rank’ or accessible matrix size of BF
from $r$ to $2r$. Also, when $r>N_{s}$, formally the Jastrow can be important.
In all the previous cases, we have been explicitly considering a determinant
Jastrow. If one is going to use a Jastrow to enhance the wave-function, it’s
worth considering whether the determinant Jastrow is better than a more
standard alternative. In particular, a reasonable comparison is to use an
exponential Jastrow where the term beging exponentiated also comes from NN; We
call this class of wave-functions expJBF. Like a standard Jastrow, we have
that there is always a positive correction to the determinant backflow. We can
additionally test whether this is important by one more Jastrow, expJsBF which
multiplies by an additional NN-generated scalar that is allowed to be positive
and negative.

In this section, we primarily compare the energetic of these wave-functions,
to numerically check which have lower energy on our protoyptical $4\times 4$
Hubbard model (See Fig. 2(b)). We find that the energy of expJBF and expJsBF
are almost identical although the signed version gives slightly better results
at large $r$. In comparison against JBF, we find that at small rank, we get
lower energy from the exponential Jastrows than the determinant Jastrow; This
relationship seems to switch at large enough $r$.

A plausible mechanism for this difference is that the role of the Jastrow is
to be able to introduce exponentially large separation between different
configurations. An exponential Jastrow naturally has this ability; And the
determinant Jastrow essentially emulates this exponential separation through
the product of the respective matrices eigenvalues. At low $r$, it is
restricted in the number of such terms in the product and has trouble matching
the exponential Jastrow. This all said, the difference between the various
Jastrow’s are tiny especially when $r$ is large, and it is not clear that one
should take those tiny differences at this level seriously. The high-level
result is that a Jastrow marginally improves a restricted form the NNBF wave-
function, but different reasonable Jastrow’s lead to essentially the same
marginal improvement; there seems to be nothing fundamentally special or
powerful about the determinant Jastrow.

###  IV.4 Determinant Jastrows vs Effective Jastrows

Throughout this paper, we’ve been making a somewhat artificial distinction
between the Jastrow and Slater Determinant pieces of the wave-functions. There
is a sense in which one could make a different decomposition of the wave-
functions we’ve been considering, by identifying the components of wave
function amplitude as ”effective Jastrow” and ”Determinant” part. In
particular, we can write all the wave-functions considered so far as

$\Psi(c)=\det(A[c])\det(\bar{F}(c))$ (16)

where $\bar{F}$ absorbs the rest of the wave-function except the static Slater
determinant part $A[c]$. Notice that (1) $\bar{F}$ will depend on the NN in
some way and (2) unlike all the other Jastrow’s we’ve considered, the matrix
$\bar{F}$ for which we are taking the determinant of may involve row-
selection. In particular, we can write NNBF as

$\Psi_{\text{NNBF}}(c)=\det(A[c])\det(I-CA[c]^{-1}B[c])$ (17)

and r-HF as

$\Psi_{{r\text{-HF}}}(c)=\det(A[c])\det(F-CA[c]^{-1}B[c])$ (18)

Naively, it’s not obvious why (or if) these wave-functions should be superior
to a Slater Determinant dressed with a determinant Jastrow - i.e.
$\bar{F}=F_{\bm{\theta}}(c)$ which we call JSD. Numerically, we have examined
this wave-function and find JSD is much worse in energy than any of the NNBF-
type wave-functions. It is naively surprising that the NN in JSD is unable to
learn $F_{\theta}(c)=\det(I-CA[c]^{-1}B[c])$ while it should in some
$w\rightarrow\infty$ limit.

To understand the mechanism for this distinction, let us first look at the
determinant Jastrow $\det F(c)$. As NN, as well as determinant, is a smooth
function where small change of the input value will only cause a small change
to the output, the determinant Jastrow will not change drastically as the
input $c$ moves between the nearest configurations. However, the amplitude of
the exact ground state is not necessarily a smooth function of the
configurations, thus is likely to go beyond the modulation ability of a
determinant Jastrow at moderate $w$.

On the other hand, as of the effective Jastrow from NNBF
$\det(I-CA[c]^{-1}B[c])$, the presence of $A[c]^{-1}B[c]$, which is
essentially not a smooth function of input configuration due to row-selection,
will impose a drastic change on top of the smooth change from NN, making it
possible to capture the pattern from the true ground state, while still
maintaining the generalization ability of the NN even at small $w$.

We can explicitly test and quantify whether the effective Jastrow with row-
selection generates significantly stronger modulation than directly using the
NN output as determinant Jastrow. To accomplish this, we want to look at the
difference in amplitudes for configurations that are ”nearest neighbors”.
First, we look at a prototypical example for $\bar{F}$ of both ${r\text{-HF}}$
and $r$-JSD at $r=N_{s}$ where we measure the ratio of nearest-neighbor
amplitudes of the effective Jastrow,
$f_{\text{ratio}}(c_{j})=\det\bar{F}(c_{j+1})/\det\bar{F}(c_{j})$ for random
configurations sampled from $|\Psi(c)|^{2}$. In Fig. 7(a), we can readily see
an obvious distinction between the $f_{\text{ratio}}$, from ${r\text{-HF}}$
and ${r\text{-JSD}}$. More specifically, the ratios from ${r\text{-JSD}}$ show
that the determinant Jastrow, $\det\bar{F}=\det F_{\bm{\theta}}$, modulates
the wave function amplitude in a weak way for both the sign structure and the
magnitude. In contrast, the effective Jastrow from ${r\text{-HF}}$ modulates
both the magnitude and the sign structure in a strong way.

We can further quantify this modulation by computing the real and imaginary
components of the log-ratio $f^{\prime}(c_{j})=\ln f_{\text{ratio}}(c_{j})$
which act as an effective derivative (see Sec. V.4 for details), and the norms
of them measure the modulation on the magnitude and sign structure
respectively. Using the log importantly lets us see rapid relative changes in
the amplitude even when $|f_{\text{ratio}}|\ll 1$. In Fig. 7(b,c), we again
observe an sharp distinction between the effective Jastrow and determinant
Jastrow from all types of the wave functions considered in this work. Note
that the larger the norm is, the stronger the modulation is, so that the
modulation from effective Jastrow is nearly $\exp(0.75)\approx 2.1$ times
stronger, based on the differences of the average norms in Fig. 7(b,c), than
the determinant Jastrow on both magnitude and sign structure. Interestingly,
we can also observe from Fig. 7(b,c) that at lower rank, the effective Jastrow
modulations from ${r\text{-HF}}$, ${r\text{-JBF}}$ and ${r\text{-BF}}$ are
weaker than their counterparts at higher rank, which is consistent with their
relatively higher energies in Fig. 2.

Alternatively, we compute the log-ratio of the matrix elements from SPO, and
show the mean and variance over these $L_{s}\times N_{s}$ elements in Fig. 8.
From them, we can extract some information about the modulation on the SPO
from various types of wave function ansatz. Although the construction of SPO
does not involve the row-selection as the wave function amplitude, we still
observe weaker modulations from JSD. The mean for JSD is always lower than
those from NNBF-type wave functions, in both magnitude (real part of log-
ratio) and sign (imaginary part of log-ratio) modulation, indicating that on
each SPO element, JSD changes the values in a relatively weaker way. Moreover,
the smaller variance from JSD shows that JSD impose a more uniform modulation
over the entire SPO, suggesting that JSD is less capable of identifying which
part of SPO is worth stronger modulation than the rest.

Considering the geometry for the SPO, it is easy to see that, in the cases of
$r\leq N_{s}$, the SPO of JSD is always within the linear space spanned by the
column vectors from the static matrix $A$, while NNBF is dynamically changing
the space for the SPO given different configuration inputs. In Appendix. G, we
demonstrate that it is impossible for JSD to reproduce the SPO from NNBF.

Overall we see that the effective Jastrow with row-selection imposes stronger
modulation on the wave-function amplitudes, achieving better performances on
ground state approximation. The investigation on SPO also suggests that NNBF
produces a broader class of SPO, that goes beyond the representability of
Jastrow-Determinant type wave functions, indicating that NNBF is likely to be
a better ansatz candidate for Fermionic systems.

##  V Numerical Methods

To further support our statements about the relation among the NN wave
functions, we provide some numerical evidences with the methods given in this
section, including ground state approximation for Hubbard model, mutual
learning of wave functions and log-ratio for wave function modulation.

###  V.1 Model and Wave Function Architecture

The variational wave functions are optimized to approximate the ground state
of Hubbard model

$\hat{H}=-t\sum_{\langle
i,j\rangle,\sigma}(\hat{c}_{i\sigma}^{\dagger}\hat{c}_{j\sigma}+h.c.)+U\sum_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}$
(19)

where $\hat{c}_{i\sigma}$ ($\hat{c}_{i\sigma}^{\dagger}$) is annihilation
(creation) operator for spin-$\sigma$ ($\sigma=\uparrow,\downarrow$) electron
on site $i$, and
$\hat{n}_{i\sigma}=\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}$ is the number
operator. We will consider this model on an $L=4\times 4$ square lattice with
periodic boundary conditions on both spatial directions, interaction strength
$U/t=8$, and $\frac{1}{8}$ hole doping, i.e., $N_{\uparrow}=N_{\downarrow}=7$
[25].

The neural network architecture is designed to have two hidden layers, each
with 2,048 neurons and ReLU activation function, parametrized by
${\bm{\theta}}$, the set of weights and bias inside the NN. The output is an
oversized array, which is reshaped into an $r\times N_{s}$ matrix
$C_{\bm{\theta}}(c)$ and an $r\times r$ matrix $F_{\bm{\theta}}(c)$, which,
together with the static matrices (with constant entries independent of $c$)
$A$ and $B$, are used to construct the SPO and corresponding determinant
matrices for various types of wave functions (see Fig. 1(e)).

Figure 2: Variational energy results from ground state approximation on
Hubbard model ($L=4\times 4$, PBC, $U/t=8$, $N_{\uparrow}=N_{\downarrow}=7$).
The two vertical dotted line show the locations of $r=N_{s}$(blue) and
$r=L_{s}$(pink), and the horizontal dashed line in (a) is the lowest energy
from ${r\text{-BF}}$ at $r=L_{s}$ from Fig.5.  Figure 3: Replot of the energy
data from Fig. 2 for comparison between ${r\text{-JBF}}$ (blue) and $2r$-BF
(green). In addition to the results for previous model ($L=4\times 4$, PBC,
$U/t=8$, $N_{\uparrow}=N_{\downarrow}=7$)(circle markers), the results from a
larger lattice ($L=4\times 8$, PBC, $U/t=8$,
$N_{\uparrow}=N_{\downarrow}=14$)(square markers), and results for small-size
NN (width=10)(empty circle markers) are displayed as well.

###  V.2 Ground State Approximation

To optimize the wave functions towards the ground state, in this work, we use
supervised wave-function optimization with Adam [26, 27] (SWO-Adam, see
Appendix H). Within SWO, at each so-called imaginary-time step $t$, we set a
target state
$|\Phi_{t}\rangle=e^{-\tau\hat{H}}|\Psi_{t-1}\rangle\approx(1-\tau\hat{H})|\Psi_{t-1}\rangle$
for the present variational wave function $|\Psi_{t}\rangle$. And a loss
function to minimize in this procedure is defined as the logarithm of the
wave-function fidelity:

$\mathcal{L}=-\ln\frac{\langle\Psi_{t}|\Phi_{t}\rangle\langle\Phi_{t}|\Psi_{t}\rangle}{\langle\Psi_{t}|\Psi_{t}\rangle\langle\Phi_{t}|\Phi_{t}\rangle}$
(20)

The gradient of Eqn.(20) is estimated with Monte Carlo sampling from
probability distribution $p(c)=|\Psi_{t-1}(c)|^{2}$,

$\begin{split}\frac{\partial\mathcal{L}}{\partial{\bm{\theta}}}&\approx
2\sum_{\\{c|p(c)\\}}\Re\left[\alpha(c)\operatorname{Tr}\left(\widetilde{A}_{\bm{\theta}}(c)^{-1}\frac{\partial\widetilde{A}_{\bm{\theta}}(c)}{\partial{\bm{\theta}}}\right)\right]\end{split}$
(21)

where the coefficient $\alpha(c)$ is given by

$\alpha(c)=\frac{|\psi(c)|^{2}}{\sum_{\\{c\\}}|\psi(c)|^{2}}-\frac{\psi(c)^{*}\phi(c)}{\sum_{\\{c\\}}\psi(c)^{*}\phi(c)}$
(22)
$\psi(c)=\frac{\Psi_{t}(c)}{\Psi_{t-1}(c)},\quad\phi(c)=\frac{\Phi_{t}(c)}{\Psi_{t-1}(c)}$
(23)

By iteratively feeding the estimated gradient Eqn. (21) into Adam optimizer to
minimize the loss function Eqn. (20) at each step $t$, the variational state
$|\Psi_{t}\rangle$ is undergoing imaginary time evolution stochastically, with
energy being decreased until convergence.

The results, in terms of relative energy error to the exact ground energy for
the chosen model $E_{0}=-11.868$ [25], are shown in Fig. 2.

###  V.3 Mutual Learning of Wave Functions

The mutual learning of wave functions is the optimization of one wave function
$|\Psi_{1}\rangle$, e.g., ${r\text{-JBF}}$, towards the other one
$|\Psi_{2}\rangle$, e.g., ${r\text{-HF}}$. Two metrics are defined for this
purpose: one is the wave function fidelity

$\mathcal{F}_{\text{wf}}=\frac{\langle\Psi_{1}|\Psi_{2}\rangle\langle\Psi_{2}|\Psi_{1}\rangle}{\langle\Psi_{1}|\Psi_{1}\rangle\langle\Psi_{2}|\Psi_{2}\rangle}$
(24)

and the other is the square distance between the SPO’s.

$\mathcal{L}_{\text{spo}}=\frac{1}{S}\sum_{\\{c\\}}\left|\widetilde{A}_{1}(c)-\widetilde{A}_{2}(c)\right|^{2}$
(25)

where $\widetilde{A}_{1/2}(c)$ are the $L_{s}\times N_{s}$-size SPO matrices
generated from $|\Psi_{1/2}\rangle$ for configuration $c$, and $S$ is number
of samples in $\\{c\\}$. These two metrics are used as the loss functions for
the optimization, and the target states are those from the ground state
approximation for Hubbard model.

In Fig. 4, we show the fidelity after mutual learning among ${r\text{-HF}}$,
${r\text{-JBF}}$ and $2r$-BF ($r\leq N_{s}$) by maximizing the wave function
fidelity Eqn. (24). In Fig. 5, we show the same results from learning
${r\text{-JNNBF}}$ ($r\leq N_{s}$) with $L_{s}$-BF, and the comparison between
energy after the optimization.

In the Fig. 6, we show the results from mutual learning by minimizing the SPO
distance Eqn. (25). Notice that optimizing SPO distance is more strict than
optimizing the wave function fidelity, as there is a significant amount of
degree of freedom for the SPO to give the same amplitude during the operations
of the row selection and determinant evaluation. Hence, in Fig. 6(b,c) we
further show the SPO fidelity and the wave function fidelity that are
consistent with the SPO distance in Fig. 6(a). Here the SPO fidelity, defined
as

$\mathcal{F}_{\text{spo}}=\sum_{\\{c\\}}\frac{\det(\mathcal{S}_{1}(c)^{\dagger}\mathcal{S}_{2}(c))\det(\mathcal{S}_{2}(c)^{\dagger}\mathcal{S}_{1}(c))}{\det(\mathcal{S}_{1}(c)^{\dagger}\mathcal{S}_{1}(c))\det(\mathcal{S}_{2}(c)^{\dagger}\mathcal{S}_{2}(c))}$
(26)

measures the average overlap between the SPO’s for each configuration $c$.

Figure 4: Wave function fidelity from mutual learning of different wave
functions by maximizing the wave function fidelity: Optimization of $r$-HF
towards $r$-JBF (red circle), and $2r$-BF (red square); $r$-JBF towards $r$-HF
(blue circle), and $2r$-BF (blue diamond); $2r$-BF towards $r$-HF (green
square), and $r$-JBF (green diamond). The target states are from the low
energy space of Hubbard model ($L=4\times 4$, PBC, $U/t=8$,
$N_{\uparrow}=N_{\downarrow}=7$).

.

Figure 5: Results of mutual learning of $L_{s}$-BF towards $r$-JNNBF by
maximizing the wave function fidelity. (a) Wave function fidelity (b)
variational energy after the optimization. The horizontal lines in (b) is the
lowest energy of $L_{s}$-BF after the optimization.

.

Figure 6: (a) Square distance of SPO’s per configuration and matrix element,
(b) SPO fidelity ($1-\mathcal{F}_{\text{spo}}$ in the plot) and (c) wave
function fidelity ($1-\mathcal{F}_{\text{wf}}$ in the plot) from mutual
learning of different wave functions by minimizing the SPO distance:
Optimization of $r$-HF towards $r$-JBF (red circle); $r$-JBF towards $r$-HF
(blue circle); $2r$-BF towards $r$-HF (green square), and $r$-JBF (green
diamond). The target states are from the ground state approximation for
Hubbard model ($L=4\times 4$, PBC, $U/t=8$, $N_{\uparrow}=N_{\downarrow}=7$).

###  V.4 Amplitude Modulation

To look into the modulation from the neural networks on the wave function , we
compute the log-ratio $f^{\prime}(c)$ for some quantity $f(c)$ as follows:

Given a Markov chain for Monte Carlo sampling from wave function
$|\Psi\rangle$, we collect all the configuration samples from this chain but
remove those successively duplicated ones, such that the remaining ones have
the neighboring samples differing by the position of exactly one Fermion,
i.e., $|c_{j}-c_{j+1}|=2$ for all $c_{j}$ from the collected sample set
$\\{c_{j}\\}_{j=1}^{J}$. Then, between neighboring samples, we compute the
ratio of the quantity of interest $f(c)$ from them, i.e.,
$f_{\text{ratio}}(j)=\frac{f(c_{j+1})}{f(c_{j})}$.

Here the following quantities will be considered: (1) $f(c)=\det
F_{{\bm{\theta}}}(c)$, where $\det F_{{\bm{\theta}}}(c)$ is the determinant
Jastrow from the definition of wave function ansatz, e.g., the one from Eqn.
(7); (2) $f(c)=\det\bar{F}(c)$, the effective Jastrow factor obtained by
factoring out the static Slater determinant part of the wave function, e.g.,
$\det\bar{F}(c)=\det\left(F_{\bm{\theta}}(c)-C_{\bm{\theta}}(c)A[c]^{-1}B[c]\right)$
for ${r\text{-HF}}$; (3) $f(c)=\mathcal{S}_{kj}(c)$, the matrix element from
the SPO. In Fig. 8(a), we show a snapshot of the ratios for the effective
Jastrow from ${r\text{-HF}}$ and $r$-JSD.

Next we take the logarithm of $f_{\text{ratio}}(j)$, i.e., $f^{\prime}(j)=\ln
f_{\text{ratio}}(j)$, because of the following reasons: (1) It is independent
of the norm of wave function, allowing for comparison among different un-
normalized wave functions; (2) The modulation on the magnitude and phase
structure can be studied separately, by treating the real and imaginary part
of $f^{\prime}(j)$ separately.

To have a global estimation of the modulation, we calculate the average square
norm of the real and imaginary part from $f^{\prime}(j)$ over the sample set
$\\{c_{j}\\}_{j=1}^{J}$,

$\mathcal{F}_{\text{magn}}=\frac{1}{J}\sum_{j=1}^{J}|\Re
f^{\prime}(j)|^{2},\quad\mathcal{F}_{\text{phase}}=\frac{1}{J}\sum_{j=1}^{J}|\Im
f^{\prime}(j)|^{2}$ (27)

In Fig. 7(b,c), we show these two estimations for the Jastrow factors from all
types of wave functions considered in this work. Moreover, in Fig. 8, we
estimate the modulation on each element of SPO matrix, and then show the mean
and the variance of the modulation over the SPO elements for the selected
types of wave functions.

Figure 7: (a) Samples of $f_{\text{ratio}}$ for the effective Jastrow factor
from ${r\text{-HF}}$ and $r$-JSD ($r=N_{s}$) along a Markov chain. (b) The
average square norm of the real part of $f^{\prime}$. (c) The average square
norm of the imaginary part of $f^{\prime}$. In (b) and (c), the results for
the defined Jastrow factors ($F$) and the effective Jastrow factors
($F^{\prime}$) are plotted for all the wave functions considered in this
paper. The wave functions are from the ground state approximation for Hubbard
model ($L=4\times 4$, PBC, $U/t=8$, $N_{\uparrow}=N_{\downarrow}=7$). Figure
8: The average norm of effective derivative on the SPO elements for
${r\text{-HF}}$, ${r\text{-JBF}}$, $2r$-BF and $r$-JSD. (a) The mean of
$\mathcal{F}_{\text{magn}}$ over all the SPO elements; (b) The mean of
$\mathcal{F}_{\text{phase}}$; (c) The variance of $\mathcal{F}_{\text{magn}}$
over all the SPO elements; (d) The variance of $\mathcal{F}_{\text{phase}}$;
The wave functions are from the ground state approximation for Hubbard model
($L=4\times 4$, PBC, $U/t=8$, $N_{\uparrow}=N_{\downarrow}=7$).

##  VI Summary and Overview

In this work, we have considered a series of different Jastrow-NNBF wave-
functions and their relation with HFDS. We show that HFDS can be written in
the Jastrow-NNBF form with an $r\times r$ determinant Jastrow, $\det
F_{\bm{\theta}}(c)$, and an additive correction to the SPO of the form
$BF_{\bm{\theta}}(c)^{-1}C_{\bm{\theta}}(c)$, where $F^{-1}C$ is an $r\times
N_{s}$ matrix (see Fig. 1(a) bottom and Eqn. (7)). This is not of full-rank
for $r\leq N_{s}$ and is to be contrasted with the standard NNBF where the
full-rank additive correction (an $L_{s}\times N_{s}$ matrix) is generated
directly.

The determinant Jastrow piece turns out to not be critical for NNBF wave-
functions: it can be replaced by an exponential Jastrow which achieves similar
(and often better) energetics (see Fig. 2(b)) or alternatively when $r\leq
N_{s}$ it can be absorbed into the additive correction to the SPO resulting in
a new more complicated additive correction $BC’$ where $C’$ is now “taller” by
size $r$ (see r.h.s of Fig. 1(a,b)) and Eqns. (13)-(15)).

We are then left comparing various different additive corrections of the form
$BC’$ where $C’$ is generated from various combinations of neural network
output. Two general principles appear here. First and unsurprisingly, larger
$r$ is better until $r=L_{s}$. Larger $r$ allows access to more flexibility in
the matrices that can be generated including, out to $r=N_{s}$, access to
higher rank matrices. Beyond $r\geq L_{s}$, we can reduce additive corrections
to $r=L_{s}$. Secondly, the different ways we’ve considered to use the NN
output to produce a $C’$ of size $2r$ (i.e. $r$-JBF, $r$-HF, $2r$-BF, etc) all
are roughly similar in energetics and expressibility and, where they differ,
favor the states which are simpler, more direct, and more generic updates
(i.e. $2r$-BF). In particular, we see that $r$-HF and $r$-JBF are empirically
very close energetically (with a slight edge to $r$-JBF, see Fig. 2(a)) and
that we can optimize them towards each other with relatively low infidelity in
each direction(see Fig. 4). $2r$-BF is consistently lower in energy than both
$r$-JBF and $r$-HF(see Fig. 3) and optimizing toward and from $2r$-BF can be
done with low infidelity (see Fig. 4). In the limit of large width, $2r$-BF
can represent any SPO in $r$-JBF and $r$-HF; and with two extra layers can
represent any SPO in $r$-JBF at finite width. Amongst the various NNBF wave-
functions, these two principals suggest that using standard NNBF is better
energetically and likely spans the space of all the other wave-functions we’ve
considered when $r\leq N_{s}$. For $r>N_{s}$, it likely makes sense to use
Jastrow-NNBF with either an exponential or determinant Jastrow (although the
improvement in energy even then is somewhat marginal).

NNBF and HFDS have both been successful in capturing highly accurate ground
states. We compare a wave-function JSD which is similar on the surface, a
determinant Jastrow times a Slater Determinant, and find that it is
significantly worse. To help understand what aspects are endowing NNBF with
its particular advantage, we rewrite all our wave-function as a (non-
configuration independent) Slater Determinant times an “effective Jastrow” and
then compare properties of these effective Jastrows. In cases, such as NNBF
and HFDS where the effective Jastrow involves row-selection, we see that both
the sign and amplitude are rapidly oscillating as a function of configuration,
something which does not happen at moderate width for JSD. This suggests that
some of the advantage of NNBF is coming from the ability to rapidly change the
wave-function.

Unlike the tensor network ansatz, which are well understood in terms of
entanglement [28, 29, 30, 31, 32, 33, 34], the development of NQS largely
relies on heuristic experiments, and lacks a comprehensive understanding for
the underlying principles [35, 36, 37, 38, 39, 40]. We hope our work paves the
way for unifying various types of NQS and identifying the crucial components
for an efficient wave-function ansatz, which, in turn, may lead to the
discovery of more accurate and scalable NQS.

###### Acknowledgements.

We acknowledge the helpful conversations with Andrew Liu, Di Luo and Kieran
Loehr. This work made use of the Illinois Campus Cluster, a computing resource
that is operated by the Illinois Campus Cluster Program (ICCP) in conjunction
with the National Center for Supercomputing Applications (NCSA) and which is
supported by funds from the University of Illinois at Urbana-Champaign.

##  Appendix A Validity of Eqn. (7)

In Sec. IV, we write down the transformed wave function expression for HFDS,
Eqn. (7), by assuming that the backflow matrix $F_{\bm{\theta}}(c)$ is
invertible. Here we are proving that even if originally $F_{\bm{\theta}}(c)$
is not invertible for all the configuration input $c$’s, we can always
construct an exact transformation (without changing the amplitude output or
the neural network architecture) on HFDS, such that

  1. 1.

$F_{\bm{\theta}}(c)$ is invertible for all $c$’s when the backflow matrix
$\begin{bmatrix}C_{\bm{\theta}}(c)&F_{\bm{\theta}}(c)\end{bmatrix}$ for
${r\text{-HF}}$ is at least rank $r$.

  2. 2.

Eqn. (7) is still valid for all $c$’s when there exists backflow matrix
$\begin{bmatrix}C_{\bm{\theta}}(c)&F_{\bm{\theta}}(c)\end{bmatrix}$ for
${r\text{-HF}}$ with rank less than $r$.

Given the original HFDS wave function form

$\begin{split}\Psi_{r\text{-HF}}(c)&=\det\begin{bmatrix}A[c]&B[c]\\\
C_{\bm{\theta}}(c)&F_{\bm{\theta}}(c)\end{bmatrix}\end{split}$ (28)

We can divide the configuration basis sets into two disjoint sets: invertible
set $V$ whose elements $c$ give invertible $F_{\bm{\theta}}(c)$, and
uninvertible set $\bar{V}$ with uninvertible $F_{\bm{\theta}}(c)$, or
equivalently, $\det F_{\bm{\theta}}(c)=0$. The underlying reason for $\det
F_{\bm{\theta}}(c)=0$ is that the $r\times r$ matrix $F_{\bm{\theta}}(c)$ is
made up of less than $r$, say, $r_{c}$, linearly independent vectors.
Performing Gaussian elimination on the columns of $F_{\bm{\theta}}(c)$, we
will end up with

$F_{\bm{\theta}}(c)\rightarrow\begin{bmatrix}\widetilde{f}_{1,{\bm{\theta}}}(c)&\widetilde{f}_{2,{\bm{\theta}}}(c)&\cdots&\widetilde{f}_{r_{c},{\bm{\theta}}}(c)&0&\cdots&0\end{bmatrix}$
(29)

where $\\{\widetilde{f}_{j,{\bm{\theta}}}(c)\\}_{j=1}^{r_{c}}$ are
$r$-dimensional column vectors. Because of the presence of zero vectors in
Eqn. (29), we get the unwanted $\det F_{\bm{\theta}}(c)=0$ and uninvertible
$F_{\bm{\theta}}(c)$. Nevertheless, suppose that a $c\in\bar{V}$ has backflow
matrix $\begin{bmatrix}C_{\bm{\theta}}(c)&F_{\bm{\theta}}(c)\end{bmatrix}$ is
of rank at least $r$, but unfortunately has rank-deficient
$F_{\bm{\theta}}(c)$ as given by Eqn. (29), we can mix it with new linearly
independent vectors from $C_{\bm{\theta}}(c)$ to remove those zero vectors,
such that $F_{\bm{\theta}}(c)$ is full-rank and invertible. To do so, we
introduce an $(N_{s}+r)\times(N_{s}+r)$ mixing matrix $K$, which is upper-
triangular with 1’s on the diagonal:

$K=\begin{bmatrix}1&k_{12}&\cdots&k_{1,N_{s}+r}\\\ 0&1&\cdots&k_{2,N_{s}+r}\\\
0&0&\ddots&\vdots\\\ 0&0&\cdots&1\end{bmatrix}$ (30)

We right multiply the ${r\text{-HF}}$ matrix with this mixing matrix $K$, the
result of which will be transparent if we rewrite the ${r\text{-HF}}$ matrix
in terms of column vectors:

$\begin{split}&\begin{bmatrix}A[c]&B[c]\\\
C_{\bm{\theta}}(c)&F_{\bm{\theta}}(c)\end{bmatrix}K\\\
&=\begin{bmatrix}h_{1}&\cdots&h_{N_{s}}&f_{1}&\cdots&f_{r}\end{bmatrix}K\\\
&=\begin{bmatrix}\overline{h}_{1}&\cdots&\overline{h}_{N_{s}}&\overline{f}_{1}&\cdots&\overline{f}_{r}\end{bmatrix}\end{split}$
(31)

where

$\overline{h}_{m}=\sum_{j=1}^{m-1}k_{j,m}h_{j}+h_{m}$ (32)
$\overline{f}_{m}=\sum_{j=1}^{N_{s}}k_{j,N_{s}+m}h_{m}+\sum_{j=1}^{m-1}k_{N_{s}+j,N_{s}+m}f_{j}+f_{m}$
(33)

Effectively, the $K$ matrix adds to each column vector a linear combination of
previous ones. By choosing appropriate matrix entries for $K$, the right
multiplication of $K$ would increase the rank of $F_{\bm{\theta}}(c)$ up to
$r$ (full-rank) via adding new column vectors from $C_{\bm{\theta}}(c)$. As
long as those $c\in\bar{V}$ satisfy the condition that the backflow matrix
$\begin{bmatrix}C_{\bm{\theta}}(c)&F_{\bm{\theta}}(c)\end{bmatrix}$ is of rank
at least $r$, there always exists some column vector from $C_{\bm{\theta}}(c)$
to mix in to make $F_{\bm{\theta}}(c)$ full-rank.

Moreover, as $K$ is an upper-triangular square matrix with constant entries
(not dependent on $c$), and $\det K=1$, this mixing will not change the wave
function amplitude of HFDS. The practical realization on HFDS for this mixing
is simply to add one more linear layer with weights constructed from $K$ to
the output layer of the NN that generates the backflow matrix. This new linear
layer can just be absorbed into the last layer of the original NN, as a
result, this mixing can be achieved by simply adjusting the weights and bias
on the original output layer even without changing the NN architecture.

Since the effect of the mixing matrix $K$ is to mix new vectors into
$F_{\bm{\theta}}(c)$, the choices for the non-zero entries of $K$ are almost
generic, as long as it circumvents the situation where the mixing accidentally
decreases the rank for those originally full-ranked $F_{\bm{\theta}}(c)$,
$c\in V$. But this only happens for a zero-measure set in the parameter space
defining $K$, which can be excluded from the determination of $K$.

At last, we consider the situation where the backflow matrix has rank less
than $r$, which indicates that there are not sufficient column vectors from
$C_{\bm{\theta}}(c)$ to make $F_{\bm{\theta}}(c)$ full-rank, such that
$F_{\bm{\theta}}(c)$ is always rank-deficient, i.e., $\det
F_{\bm{\theta}}(c)=0$ no matter what mixing matrix $K$ we have chosen. In this
case, not only the columns of the $r\times(N_{s}+r)$ backflow matrix are
linearly dependent, but the rows are also linearly dependent. As a
consequence, the ${r\text{-HF}}$ matrix for HFDS is rank-deficient, and will
give zero amplitude from its determinant. Nevertheless, Eqn. (7) is still
valid, since it also gives zero amplitude as the original HFDS amplitude, due
to $\det F_{\bm{\theta}}(c)=0$. Although it brings with an ill-defined
$F^{-1}$ in the expression, in practice, this none-valued output can be
replaced with 0 during programming without affecting the performance.

In fact, this mixing procedure is also able to bring $F_{\bm{\theta}}(c)$ away
from rank-deficient regimes of ${\bm{\theta}}$, where $F_{\bm{\theta}}(c)$ has
extremely small but nonzero singular values and $F^{-1}$ has extremely large
entries, via, for example, optimizing the entries of matrix $K$ with some
proper loss function. By doing so, the learning of ${r\text{-HF}}$ with
${r\text{-JBF}}$ or ${r\text{-BF}}$ will be practically easier.

##  Appendix B Explicit Construction for Subset SPO

In Sec. IV.1, we give statements on the relation between SPO’s from different
classes of wave function ansatz. Here we provide explicit construction of SPO
to support some of these statements.

  1. 1.

We start with ${r\text{-HF}}\underset{\text{spo}}{\subset}(r+1)\text{-HF}$ for
any $r$. Recall that the SPO for ${r\text{-HF}}$ is given by $A-BF^{-1}C$,
which is an $L_{s}\times N_{s}$ matrix; it is straightforward to see that it
can be reproduced in the class of $(r+1)$-HF, by letting

$A^{\prime}=A,\quad B^{\prime}=\begin{bmatrix}B&\vline&\begin{matrix}0\\\ 0\\\
\vdots\\\ 0\end{matrix}\end{bmatrix}$ (34) $C^{\prime}=\begin{bmatrix}C\\\
\hline\cr\begin{matrix}0&0&\cdots&0\end{matrix}\end{bmatrix},\quad
F^{\prime}=\begin{bmatrix}F&\vline&\begin{matrix}0\\\ 0\\\ \vdots\\\
0\end{matrix}\\\
\hline\cr\begin{matrix}0&0&\cdots&0\end{matrix}&\vline&1\end{bmatrix}$ (35)

such that $A^{\prime}-B^{\prime}F^{\prime-1}C^{\prime}=A-BF^{-1}C$ for any
$c$. Note that the constant elements on the backflow matrices $C$ and $F$ can
be constructed by assigning the weights on the final linear layer to those
positions as 0 and the bias as the corresponding constant values (0 or 1
here).

Similarly, we can also show that
${r\text{-JBF}}\underset{\text{spo}}{\subset}(r+1)\text{-JBF}$ and
${r\text{-BF}}\underset{\text{spo}}{\subset}(r+1)\text{-BF}$, using the same
construction as ${r\text{-HF}}$ (Eqn. (34) and (35)), such that
$A^{\prime}-B^{\prime}C^{\prime}=A-BC$ for any $c$.

  2. 2.

Next, we show that for any $r$, including the cases of $r\geq L_{s}$, we have
${r\text{-JBF}}\underset{\text{spo}}{\subset}L_{s}\text{-JBF}$ and
${r\text{-BF}}\underset{\text{spo}}{\subset}L_{s}\text{-BF}$, both of which
amount to showing $A^{\prime}-B^{\prime}C^{\prime}=A-BC$, where $B^{\prime}$,
$C^{\prime}$ are now matrices of size $L_{s}\times L_{s}$ and $L_{s}\times
N_{s}$, respectively. Copy the neural network generating $C$ for
${r\text{-BF}}$ or ${r\text{-JBF}}$ and add one more linear layer on top of
it. The weights and bias of this additional layer are given by

$W_{(mn)(kl)}=\delta_{nl}B_{mk},\quad b_{(mn)}=0$ (36)

where we use the matrix index $(mn),(kl)$ to index the weights and bias. The
output of this new linear layer will directly give us the $L_{s}\times N_{s}$
matrix $C^{\prime}=BC$, because

$\begin{split}C^{\prime}_{mn}=&\sum_{k,l}W_{(mn)(kl)}C_{(kl)}+b_{(mn)}\\\
=&\sum_{k}B_{mk}C_{kn}\equiv(BC)_{mn}\end{split}$ (37)

Note that this additional linear layer $(W^{\prime},b^{\prime})$ can be
absorbed into previous linear layer $(W,b)$ for $C$ simply by $W\leftarrow
W^{\prime}W$, $b\leftarrow W^{\prime}b+b^{\prime}$. As we have gotten
$C^{\prime}=BC$, this construction is completed by the assignments of
$A^{\prime}=A$ and $B^{\prime}=I_{L_{s}}$.

  3. 3.

Here we demonstrate that, given an NNBF that outputs matrix $M$ (see Eqn. 2),
we can generate an $L_{s}$-BF, that outputs matrix $C$, such that $BC=M$.

We start by rewriting a general $L_{s}$-BF where we absorb $B$ into the NN,
which amounts to adding one more linear layer with weights as

$W^{\prime}=\begin{bmatrix}B&\bm{0}&\cdots&\bm{0}\\\ \bm{0}&B&\cdots&\bm{0}\\\
\vdots&\vdots&\ddots&\vdots\\\ \bm{0}&\bm{0}&\cdots&B\\\ \end{bmatrix}$ (38)

where this is an $N_{s}\times N_{s}$ block matrix; We can combine it with the
weight matrix from previous linear layer $W$ which is given as

$W=\begin{bmatrix}W_{1}\\\ W_{2}\\\ \vdots\\\ W_{N_{s}}\end{bmatrix}$ (39)

where $W_{h}$, $h=1,2,...,N_{s}$ ($H$ is the width of last hidden layer) are
arbitrary $r\times H$ matrices. Now the resulting weight
$W^{\prime\prime}=W^{\prime}W$ is still a generic $(L_{s}N_{s})\times H$
matrix with no constraints, because each block of $W^{\prime\prime}$ is given
by

$W^{\prime\prime}_{h}=BW_{h}$ (40)

and is of full-rank $L_{s}$ as $B$ and $W_{h}$ are of size $L_{s}\times r$ and
$r\times H$, respectively, with $r\geq L_{s}$. We can assume $B$ is of full
rank $L_{s}$ (if not, it can be brought to full rank by mixing with $A$, see
Appendix. A), such that there are complete sets of linearly independent
vectors within it. Finally, the arbitrary choice of $W_{h}$ enables
$W^{\prime\prime}_{h}$ to be any $L_{s}\times H$ matrix.

The reverse direction of the relation between NNBF and ${r\text{-BF}}$ ($r\geq
L_{s}$) is also true by simply choosing a common basis set as $B$ for each
$W_{h}$.

##  Appendix C Absorption of Determinant Jastrow into SPO

In this section, we show the derivation for the absorption of determinant
Jastrow factor into SPO for ${r\text{-HF}}$ and ${r\text{-JBF}}$.

For $r\leq N_{s}$, we embed the $r\times r$ matrix $F_{{\bm{\theta}}}(c)$ in
an $N_{s}\times N_{s}$ block-diagonal matrix with identity on the extra
diagonal block

$[F_{\bm{\theta}}]\rightarrow\begin{bmatrix}F_{\bm{\theta}}&\vline&\bm{0}\\\
\hline\cr\bm{0}&\vline&I_{N_{s}-r}\end{bmatrix}\equiv\widetilde{F}_{\bm{\theta}}$
(41)

This leaves the wave function amplitude unchanged, since $\det
F=\det\widetilde{F}$, but allows us to multiply the two square matrices for
determinant Jastrow and from SPO before the determinants are computed.

$\begin{split}\Psi_{r\text{-HF}}(c)&=\det
F_{\bm{\theta}}(c)\det((A-BF^{-1}_{\bm{\theta}}(c)C_{\bm{\theta}}(c))[c])\\\
&=\det\widetilde{F}_{\bm{\theta}}(c)\det((A-BF^{-1}_{\bm{\theta}}(c)C_{\bm{\theta}}(c))[c])\\\
&=\det\left[(A-BF^{-1}C)\widetilde{F}[c]\right]\\\
&=\det\left[(A+A(\widetilde{F}-I)-BF^{-1}C\widetilde{F})[c]\right]\end{split}$
(42)

From Eqn. (42), we see a new SPO matrix of size $L_{s}\times N_{s}$. Next we
transform it into BF.

After writing $C$ in terms of block matrices as well, i.e.,

$C=\begin{bmatrix}C_{1}&\vline&C_{2}\end{bmatrix}$ (43)

where $C_{1}$ and $C_{2}$ are of size $r\times r$ and $r\times(N_{s}-r)$,
respectively, we have

$\begin{split}BF^{-1}C\widetilde{F}&=BF^{-1}\begin{bmatrix}C_{1}&C_{2}\end{bmatrix}\begin{bmatrix}F&\bm{0}\\\
\bm{0}&I\end{bmatrix}\\\
&=B\begin{bmatrix}F^{-1}C_{1}F&F^{-1}C_{2}\end{bmatrix}\end{split}$ (44)

Meanwhile, with $A$ also in block matrix form, i.e.,

$A=\begin{bmatrix}A_{11}&\vline&A_{12}\\\ \hline\cr
A_{21}&\vline&A_{22}\end{bmatrix}=\begin{bmatrix}A_{1}&\vline&A_{2}\end{bmatrix}$
(45)

where $A_{1}$ and $A_{2}$ are of size $L_{s}\times r$ and
$L_{s}\times(N_{s}-r)$, respectively, $A_{11}$ and $A_{21}$ are of size of
$r\times r$ and $(L_{s}-r)\times r$, respectively, etc, we have

$\begin{split}A(\widetilde{F}-I)&=\begin{bmatrix}A_{11}&A_{12}\\\
A_{21}&A_{22}\end{bmatrix}\begin{bmatrix}F-I&\bm{0}\\\
\bm{0}&\bm{0}\end{bmatrix}\\\ &=\begin{bmatrix}A_{11}(F-I)&\bm{0}\\\
A_{21}(F-I)&\bm{0}\end{bmatrix}\\\
&=A_{1}\begin{bmatrix}(F-I)&\bm{0}\end{bmatrix}\end{split}$ (46)

Combining Eqn. (44) and Eqn. (46) together, we eventually obtain

$\begin{split}\Psi_{\text{HF}}(c)&=\det\begin{bmatrix}A[c]&B[c]\\\
C_{\bm{\theta}}(c)&F_{\bm{\theta}}(c)\end{bmatrix}\\\
&=\det\left((A-B^{\prime}C_{\bm{\theta}}^{\prime}(c))[c]\right)\end{split}$
(47)

where

$B^{\prime}=\begin{bmatrix}-A_{1}&B\end{bmatrix},\quad
C_{\bm{\theta}}^{\prime}(c)=\begin{bmatrix}F-I&\bm{0}\\\
F^{-1}C_{1}F&F^{-1}C_{2}\end{bmatrix}$ (48)

are of size $L_{s}\times 2r$ and $2r\times N_{s}$, respectively.

Likewise, for ${r\text{-JBF}}$, we have

$\begin{split}\Psi_{\text{JBF}}(c)&=\det\left(F_{\bm{\theta}}(c)\right)\det\left((A-BC_{\bm{\theta}}(c))[c]\right)\\\
&=\det\left((A-B^{\prime}C_{\bm{\theta}}^{\prime}(c))[c]\right)\end{split}$
(49)

where

$B^{\prime}=\begin{bmatrix}-A_{1}&B\end{bmatrix},\quad
C_{\bm{\theta}}^{\prime}(c)=\begin{bmatrix}F-I&\bm{0}\\\
C_{1}F&C_{2}\end{bmatrix}$ (50)

are of size $L_{s}\times 2r$ and $2r\times N_{s}$, respectively.

##  Appendix D Relation between Wave Function Fidelity and SPO Distance

In Fig. 6, we observe that at $r=1,2$, though the SPO distance Eqn. (25) is
small (Fig. 6(a)), the wave function fidelity Eqn. (24) is relatively worse
(Fig. 6(b)), compared with larger $r$ cases. In this section we will look into
this observation by showing the relation between the wave function fidelity
and the SPO distance.

Within VMC, the wave function fidelity can be calculated as

$\mathcal{F}_{\text{wf}}=\frac{\left(\frac{1}{S}\sum_{\\{c\\}}\frac{\Phi(c)}{\Psi(c)}\right)^{2}}{\frac{1}{S}\sum_{\\{c\\}}\left(\frac{\Phi(c)}{\Psi(c)}\right)^{2}}$
(51)

where $S$ is the size of sample set $\\{c\\}$ sampled from $|\Psi(c)|^{2}$. As
$\Phi(c)/\Psi(c)$ is close to 1 here, we rewrite it as

$\frac{\Phi(c)}{\Psi(c)}\equiv 1+x(c)$ (52)

For most of $c$, we should have $|x(c)|\ll 1$. We relate $x(c)$ to the the SPO
distance component with Taylor expansion around $\Psi(c)$ up to 1st order,

$x(c)=\frac{\det\bar{A}(c)}{\det\bar{A}_{0}(c)}-1\approx\sum_{jk}\bar{A}_{0}(c)^{-1}_{kj}\left(\bar{A}(c)-\bar{A}_{0}(c)\right)_{jk}$
(53)

where $\bar{A}(c)$ and $\bar{A}_{0}(c)$ are the square matrices in the
determinant wave function form of $|\Phi\rangle$ and $|\Psi\rangle$,
respectively. Note that the truncation in Taylor expansion is one of the
sources of error for our final approximation formula, which can nonetheless be
systemically improved by including higher-order terms.

In the following, we will perform Taylor expansion up to 2nd order on wave
function fidelity as well. The first-order term will vanish as
$\mathcal{F}_{\text{wf}}=1$ is a stationary point. After defining

$\overline{x}=\frac{1}{S}\sum_{\\{c\\}}x(c),\quad\overline{x^{2}}=\frac{1}{S}\sum_{\\{c\\}}x(c)^{2}$
(54)

and $|\overline{x}|\ll 1$, $|\overline{x^{2}}|\ll 1$, we can either use chain
rule to obtain the Hessian as

$\frac{\delta}{\delta\bar{A}(c^{\prime})_{j^{\prime}k^{\prime}}}\frac{\delta\ln\mathcal{F}_{\text{wf}}}{\delta\bar{A}(c)_{jk}}=\frac{1}{S^{2}}\bar{A}_{0}(c)^{-1}_{kj}\bar{A}_{0}(c^{\prime})^{-1}_{k^{\prime}j^{\prime}}(1-S\delta_{c,c^{\prime}})$
(55)

or directly expand on $\ln\mathcal{F}_{\text{wf}}$ in terms of $\overline{x}$
and $\overline{x^{2}}$, we obtain the approximation formula for
$\ln\mathcal{F}_{\text{wf}}$ , the deviation from perfect fidelity
$\mathcal{F}_{\text{wf}}=1$, as the variance of $x(c)$,

$\ln\mathcal{F}_{\text{wf}}\approx\overline{x}^{2}-\overline{x^{2}}$ (56)

In Fig. 9, we examine the accuracy of Eqn. (56) by comparing with the exact
results from VMC.

Figure 9: Comparison between the approximate results Eqn. (56) for the wave
function fidelity and the exact results from VMC. The target wave functions
are ${r\text{-JBF}}$ from ground state approximation, and $2r$-BF and
${r\text{-HF}}$ are optimized towards the target wave functions by minimizing
the SPO distances.

We also found that $\overline{x}^{2}\ll\overline{x^{2}}$. Then by just
comparing $\overline{x^{2}}$ with the SPO distance Eqn. (25), which we rewrite
as

$\mathcal{L}_{\text{spo}}=\frac{1}{S}\sum_{\\{c\\}}\|x(c)\|^{2}=\frac{1}{S}\sum_{\\{c\\}}\sum_{jk}|\bar{A}(c)-\bar{A}_{0}(c)|_{jk}^{2}$
(57)

we can see that wave function fidelity is effectively a weighted average of
SPO distance on each element, with the weights given by $\bar{A}_{0}(c)^{-1}$
from the target state. Further evidence is shown from the comparison between
Fig. 10(a,c) and Fig. 10(b,d), where the difference of wave function fidelity
is affected not only by the SPO distance, but also the target wave function.

Figure 10: Probability distributions of $\|x(c)\|^{2}$ and $|x(c)|^{2}$ at
$r=1$ and $14(N_{s})$. The corresponding SPO distance and wave function
fidelity are given in the legends.

##  Appendix E Reconstruction of $r$-JBF with $2r$-BF

Despite the existence of matrix multiplication within ${r\text{-JBF}}$, it is
possible to reconstruct ${r\text{-JBF}}$ rigorously from $2r$-BF with feed-
forward neural networks still. In this section we show it explicitly.

First, with the feed-forward neural network that generates backflow matrices
$C_{\bm{\theta}}(c)$ and $F_{\bm{\theta}}(c)$ for ${r\text{-JBF}}$, we add an
activation function $a(x)=\ln(x)$ to this layer (note that originally for
${r\text{-JBF}}$, this layer does not have activation function on it. Also
note that for the case where $x=0$, we can nonetheless replace it with
$x+\epsilon$, $\epsilon\ll 1$, which will not affect the subsequent steps).
Now, we have the output from this modified layer as $\ln C_{jk}$
$(j=1,2,...,r$; $k=1,2,...,N_{s})$, and $\ln F_{lm}$ $(l,m=1,2,...,r)$. In the
following, we will work on the first $r$ columns of $\ln C$, i.e., $\ln
C_{1}$, in order to reproduce $C_{1}F$ with feed-forward neural network. As
for the remaining $\ln C_{2}$, we construct the weights and bias as an
identity map, and the non-linear activation function introduced later will
bring $\ln C_{2}$ back to the original $C_{2}$ from r-JBF.

Next, we add one more layer to this neural network, whose activation function
is $a(x)=\exp(x)$, bias are 0’s, and weights are

$W^{[C]}_{(xyz)(jk)}=\delta_{xj}\delta_{zk},\quad
W^{[F]}_{(xyz)(lm)}=\delta_{zl}\delta_{ym}$ (58)

where $x,y,z=1,2,...,r$, such that the output from this layer is

$\begin{split}D_{xyz}=&\exp\left(\sum_{jk}W^{[C]}_{(xyz)(jk)}\ln
C_{1,jk}+\sum_{lm}W^{[F]}_{(xyz)(lm)}\ln F_{lm}\right)\\\
=&C_{1,xz}F_{zy}\end{split}$ (59)

At the end, we add one more linear layer with non-zero weights as

$W_{(jm)(xyz)}=\delta_{jx}\delta_{my}$ (60)

such that the output is

$\widetilde{C}_{1}=\sum_{xyz}W_{(jm)(xyz)}D_{xyz}=\sum_{z}C_{1,jz}F_{zm}\equiv(C_{1}F)_{jm}$
(61)

Overall, by introducing two more layers of widths $r^{2}(N_{s}+r)$ and
$r(N_{s}+r)$ to the original neural network for r-JBF, we are able to
reproduce the backflow matrices with $2r$-BF. Therefore, we have analytically
proved that r-JBF is a subset of r-BFat finite $w$ if we allow ${r\text{-BF}}$
to have additional layers.

##  Appendix F Computational Cost

In this section, we give the computational cost for obtaining the amplitude
from various types of NNBF wave functions.

Let us start with NN part. The input size is $L_{s}$, and the output size is
$rN_{s}+r^{2}$, and assume the width of NN is $w$, then the computational cost
due to the matrix-vector multiplication is given by
$O(w^{2}+wL_{s}+wr^{2}+wrN_{s})$.

Next, the cost of the evaluation of determinant is

$\begin{split}r\text{-HF}:\quad&O((N_{s}+r)^{3})=O(N_{s}^{3}+r^{3})\\\
r\text{-BF}:\quad&O(rN_{s}^{2}+N_{s}^{3})\\\
r\text{-JBF}:\quad&O(r^{3}+rN_{s}^{2}+N_{s}^{3})\end{split}$ (62)

In Appendix. E, we show that to reproduce r-JBF from $2r$-BF, we can add two
more layers of width $r^{3}+r^{2}N_{s}$ and $rN_{s}+r^{2}$, this will give an
extra cost as $O(r^{3}(N_{s}+r)^{2})=O(r^{3}N_{s}^{2}+r^{5})$.

##  Appendix G Failure of SPO Matching with JSD towards NNBF

In this section, we illustrate the impossibility of matching the SPO of NNBF
with JSD for rank $r\leq N_{s}$, thus, combining with the numerical evidences
Fig. 2(b) and Fig. 8, we conclude that NNBF represents a broader space of SPO.

For simplicity, we consider the case of $N_{s}$-JSD and $N_{s}$-BF. Since the
NN generating the back flow matrices has linear layer as the last layer, there
are constant matrices $b$ in the backflow matrices, then we rewrite the SPO as

$\begin{split}\mathcal{S}_{\text{JSD}}(c)&=A_{\text{JSD}}F_{\text{JSD},{\bm{\theta}}}(c)\\\
&=A_{\text{JSD}}b_{\text{JSD}}+A_{\text{JSD}}\widetilde{F}_{\text{JSD},{\bm{\theta}}}(c)\\\
\mathcal{S}_{\text{BF}}(c)&=A_{\text{BF}}-B_{\text{BF}}C_{\text{BF},{\bm{\theta}}}(c)\\\
&=A_{\text{BF}}-B_{\text{BF}}b_{\text{BF}}-B_{\text{BF}}\widetilde{C}_{\text{BF},{\bm{\theta}}}(c)\end{split}$
(63)

In order to match the SPO,
$\mathcal{S}_{\text{JSD}}(c)=\mathcal{S}_{\text{BF}}(c)$, we have to solve the
following equations for static matrices

$\begin{split}A_{\text{JSD}}&=-B_{\text{BF}}\\\
A_{\text{JSD}}b_{\text{JSD}}&=A_{\text{BF}}-B_{\text{BF}}b_{\text{BF}}\end{split}$
(64)

However, the second equation above is an over-determinant equation, so there
is no exact solution for $b_{\text{JSD}}$, unless
$A_{\text{BF}}-B_{\text{BF}}b_{\text{BF}}$ falls into the linear space spanned
by the column vectors of $A_{\text{JSD}}$. But this special case is not
generically true for the BF wave function we are considering.

In conclusion, we have showed that it is impossible for JSD to reproduce the
same SPO from BF.

##  Appendix H Supervised Wave-function Optimization

Figure 11: (a) Energy and (b) fidelity (associated with the loss function Eqn.
(20)) from ${r\text{-HF}}$ at $r=64$ during SWO optimization ($\tau=0.018$)
with optimizer SGD (red), Momentum (blue) and Adam (green) for Hubbard model
($L=4\times 8$, PBC, $U/t=8$, $N_{\uparrow}=N_{\downarrow}=14$).

In this section, we provide detailed information on the SWO method, including
the derivation for the gradient formula Eqn. (21) we use for NNBF (and its
variants) and comparison between various optimizers.

To begin with, we note that SWO is a first-order optimization method for
variational wave function to undergo imaginary time evolution in a stochastic
way. As already mentioned in Sec. V, at each imaginary-time step $t$, we set
the target state as the state evolved by imaginary time $\tau$, i.e,
$|\Phi_{t}\rangle=e^{-\tau\hat{H}}|\Psi_{t-1}\rangle\approx(1-\tau\hat{H})|\Psi_{t-1}\rangle$,
where $\tau$ needs to be small enough for the latter approximation to be
valid. In practice, we notice that the energy can be minimized successfully
with $\tau<1/E_{t=0}$, where $E_{t=0}$ is the variational energy of randomly
initialized state.

To derive the gradient formula Eqn. (21), we first take derivative of the loss
function Eqn. (20) with respect to the parameter ${\bm{\theta}}$,

$\begin{split}\frac{\partial\mathcal{L}}{\partial{\bm{\theta}}}&=-\frac{\partial_{{\bm{\theta}}}\langle\Psi_{t}|\Phi_{t}\rangle}{\langle\Psi_{t}|\Phi_{t}\rangle}-\frac{\partial_{{\bm{\theta}}}\langle\Phi_{t}|\Psi_{t}\rangle}{\langle\Phi_{t}|\Psi_{t}\rangle}+\frac{\partial_{{\bm{\theta}}}\langle\Psi_{t}|\Psi_{t}\rangle}{\langle\Psi_{t}|\Psi_{t}\rangle}\\\
&=-\frac{\sum_{c}(\partial_{\bm{\theta}}\langle\Psi_{t}|c\rangle)\langle
c|\Phi_{t}\rangle}{\sum_{c}\langle\Psi_{t}|c\rangle\langle
c|\Phi_{t}\rangle}-\frac{\sum_{c}\langle\Phi_{t}|c\rangle(\partial_{\bm{\theta}}\langle
c|\Psi_{t}\rangle)}{\sum_{c}\langle\Phi_{t}|c\rangle\langle
c|\Psi_{t}\rangle}\\\
&+\frac{\sum_{c}(\partial_{\bm{\theta}}\langle\Psi_{t}|c\rangle)\langle
c|\Psi_{t}\rangle+\sum_{c}\langle\Psi_{t}|c\rangle(\partial_{\bm{\theta}}\langle
c|\Psi_{t}\rangle)}{\sum_{c}\langle\Psi_{t}|c\rangle\langle
c|\Psi_{t}\rangle}\end{split}$ (65)

This expression is exact since we are summing over all the configurations $c$.
Next, we use the sample set $\\{c|p(c)\\}$ from probability distribution
$p(c)=|\Psi_{t-1}(c)|^{2}$ to replace the exact summation. Then, for example,
the first term of Eqn. (65) ends up as

$\begin{split}&\frac{\sum_{c}(\partial_{\bm{\theta}}\langle\Psi_{t}|c\rangle)\langle
c|\Phi_{t}\rangle}{\sum_{c}\langle\Psi_{t}|c\rangle\langle
c|\Phi_{t}\rangle}=\frac{\sum_{c}\frac{(\partial_{\bm{\theta}}\langle\Psi_{t}|c\rangle)\langle
c|\Phi_{t}\rangle}{p(c)}p(c)}{\sum_{c}\frac{\langle\Psi_{t}|c\rangle\langle
c|\Phi_{t}\rangle}{p(c)}p(c)}\\\
\approx&\frac{\sum_{\\{c\\}}\frac{(\partial_{\bm{\theta}}\langle\Psi_{t}|c\rangle)\langle
c|\Phi_{t}\rangle}{p(c)}}{\sum_{\\{c\\}}\frac{\langle\Psi_{t}|c\rangle\langle
c|\Phi_{t}\rangle}{p(c)}}=\frac{\sum_{\\{c\\}}\frac{(\partial_{\bm{\theta}}\langle\Psi_{t}|c\rangle)\langle
c|\Phi_{t}\rangle}{\langle\Psi_{t-1}|c\rangle\langle
c|\Psi_{t-1}\rangle}}{\sum_{\\{c\\}}\frac{\langle\Psi_{t}|c\rangle\langle
c|\Phi_{t}\rangle}{\langle\Psi_{t-1}|c\rangle\langle
c|\Psi_{t-1}\rangle}}\end{split}$ (66)

As we are optimizing over determinant wave function form, i.e.,
$\Psi_{\bm{\theta}}(c)=\det\widetilde{A}_{\bm{\theta}}(c)$, we can further
simplify the derivative term in Eqn. (65) as

$\begin{split}\frac{\partial\Psi(c)}{\partial{\bm{\theta}}}&=\frac{\partial\det\widetilde{A}_{\bm{\theta}}(c)}{\partial{\bm{\theta}}}\\\
&=\det\widetilde{A}_{\bm{\theta}}(c)\operatorname{Tr}\left(\widetilde{A}_{\bm{\theta}}(c)^{-1}\frac{\partial\widetilde{A}_{\bm{\theta}}(c)}{\partial{\bm{\theta}}}\right)\\\
&=\Psi(c)\operatorname{Tr}\left(\widetilde{A}_{\bm{\theta}}(c)^{-1}\frac{\partial\widetilde{A}_{\bm{\theta}}(c)}{\partial{\bm{\theta}}}\right)\end{split}$
(67)

After defining the amplitude ratio $\phi(c)$ and $\psi(c)$ as Eqn. (23) to
simplify the notations, we factor out the derivative part, and put together
all the other factors as a coefficient $\alpha(c)$ given by Eqn.(22), then we
end up with the gradient estimation formula Eqn. (21). That only the real part
of gradient is taken is because the parameter set ${\bm{\theta}}$ are real
numbers and each term in Eqn. (65) is accompanied with its complex conjugate.

Finally, we note that, as SWO is a first-order optimization method, it can be
further improved with ML optimizers, such as Momentum-based gradient descent,
Adam, etc. In Fig. 11, we compare the performance from gradient descent (GD),
momentum and Adam within SWO on the ground state approximation with
${r\text{-HF}}$. It turns out that SWO-Adam is more efficient than other two
optimizers, as it optimizes the fidelity to a higher value at each step, so as
to minimize the energy faster than the other two optimizers. Henceforth, we
use SWO-Adam for optimization in our numerical tests.

##  Appendix I Notations

In Table 1, we give a summary for the notations used in this paper.

Table 1:  Table of Notations $c$ |  $L_{s}$-size binary vector
---|---
${\bm{\theta}}$ | Parameters (weights and bias) of neural networks
$A$ |  $L_{s}\times N_{s}$ constant matrix
$A[c]$ |  $N_{s}\times N_{s}$ matrix from row selection on $A$ based on $c$
$B$ |  $L_{s}\times r$ constant matrix
$B[c]$ |  $N_{s}\times r$ matrix from row selection on $B$ based on $c$
$C_{{\bm{\theta}}}(c)$ |  $r\times N_{s}$ backflow matrix from neural networks,
$F_{{\bm{\theta}}}(c)$ |  $r\times r$ backflow matrix from neural networks,
${r\text{-HF}}$ | $\Psi(c)=\det(F_{\bm{\theta}}(c))\det\left((A-BF_{{\bm{\theta}}}(c)^{-1}C_{\bm{\theta}}(c))[c]\right)$
${r\text{-JBF}}$ | $\Psi(c)=\det(F_{\bm{\theta}}(c))\det\left((A-BC_{\bm{\theta}}(c))[c]\right)$
${r\text{-BF}}$ | $\Psi(c)=\det\left((A-BC_{\bm{\theta}}(c))[c]\right)$
${r\text{-JNNBF}}$ | $\Psi(c)=\det(F_{\bm{\theta}}(c))\Psi_{L_{s}\text{-BF}}(c)$
$r$-JSD | $\Psi(c)=\det(F_{\bm{\theta}}(c))\det A[c]$
$r$-expJBF | $\Psi(c)=\exp\left({f_{{\bm{\theta}}}(c)}\right)\Psi_{r\text{-BF}}(c)$
$r$-expJsBF | $\Psi(c)=\exp\left({f_{1,{\bm{\theta}}}(c)}\right)\cos\left({f_{2,{\bm{\theta}}}(c)}\right)\Psi_{r\text{-BF}}(c)$

## References

  * Slater [1929] J. C. Slater, The theory of complex spectra, Phys. Rev. 34, 1293 (1929).
  * Jastrow [1955] R. Jastrow, Many-body problem with strong forces, Phys. Rev. 98, 1479 (1955).
  * Mitáš and Martin [1994] L. c. v. Mitáš and R. M. Martin, Quantum monte carlo of nitrogen: Atom, dimer, atomic, and molecular solids, Phys. Rev. Lett. 72, 2438 (1994).
  * Foulkes  _et al._ [2001] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Quantum monte carlo simulations of solids, Rev. Mod. Phys. 73, 33 (2001).
  * Carleo and Troyer [2017] G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science 355, 602 (2017).
  * Hermann  _et al._ [2023] J. Hermann, J. Spencer, K. Choo, A. Mezzacapo, W. M. Foulkes, D. Pfau, G. Carleo, and F. Noé, Ab initio quantum chemistry with neural-network wavefunctions (2023).
  * Moreno  _et al._ [2022] J. R. Moreno, G. Carleo, A. Georges, and J. Stokes, Fermionic wave functions from neural-network constrained hidden states, Proceedings of the National Academy of Sciences 119, e2122059119 (2022).
  * Lovato  _et al._ [2022] A. Lovato, C. Adams, G. Carleo, and N. Rocco, Hidden-nucleons neural-network quantum states for the nuclear many-body problem, Phys. Rev. Res. 4, 043178 (2022).
  * Choo  _et al._ [2020] K. Choo, A. Mezzacapo, and G. Carleo, Fermionic neural-network states for ab-initio electronic structure, Nature Communications 11, 10.1038/s41467-020-15724-9 (2020).
  * Yoshioka  _et al._ [2021] N. Yoshioka, W. Mizukami, and F. Nori, Solving quasiparticle band spectra of real solids using neural-network quantum states, Communications Physics 4, 10.1038/s42005-021-00609-0 (2021).
  * Luo and Clark [2019] D. Luo and B. K. Clark, Backflow transformations via neural networks for quantum many-body wave functions, Phys. Rev. Lett. 122, 226401 (2019).
  * Pfau  _et al._ [2020] D. Pfau, J. S. Spencer, A. G. D. G. Matthews, and W. M. C. Foulkes, Ab initio solution of the many-electron schrödinger equation with deep neural networks, Phys. Rev. Res. 2, 033429 (2020).
  * Hermann  _et al._ [2020] J. Hermann, Z. Schätzle, and F. Noé, Deep-neural-network solution of the electronic schrödinger equation, Nature Chemistry 12, 10.1038/s41557-020-0544-y (2020).
  * Spencer  _et al._ [2020] J. S. Spencer, D. Pfau, A. Botev, and W. M. C. Foulkes, Better, faster fermionic neural networks, arXiv preprint arXiv:2011.07125 (2020).
  * Inui  _et al._ [2021] K. Inui, Y. Kato, and Y. Motome, Determinant-free fermionic wave function using feed-forward neural networks, Phys. Rev. Res. 3, 043126 (2021).
  * Nomura  _et al._ [2017] Y. Nomura, A. S. Darmawan, Y. Yamaji, and M. Imada, Restricted boltzmann machine learning for solving strongly correlated quantum systems, Phys. Rev. B 96, 205152 (2017).
  * Stokes  _et al._ [2020] J. Stokes, J. R. Moreno, E. A. Pnevmatikakis, and G. Carleo, Phases of two-dimensional spinless lattice fermions with first-quantized deep neural-network quantum states, Phys. Rev. B 102, 205122 (2020).
  * Ferrari  _et al._ [2019] F. Ferrari, F. Becca, and J. Carrasquilla, Neural gutzwiller-projected variational wave functions, Phys. Rev. B 100, 125131 (2019).
  * Huse and Elser [1988] D. A. Huse and V. Elser, Simple variational wave functions for two-dimensional heisenberg spin-½ antiferromagnets, Phys. Rev. Lett. 60, 2531 (1988).
  * Changlani  _et al._ [2009] H. J. Changlani, J. M. Kinder, C. J. Umrigar, and G. K.-L. Chan, Approximating strongly correlated wave functions with correlator product states, Phys. Rev. B 80, 245116 (2009).
  * Mezzacapo  _et al._ [2009] F. Mezzacapo, N. Schuch, M. Boninsegni, and J. I. Cirac, Ground-state properties of quantum many-body systems: Entangled-plaquette states and variational monte carlo, New Journal of Physics 11, 10.1088/1367-2630/11/8/083026 (2009).
  * Ruggeri  _et al._ [2018] M. Ruggeri, S. Moroni, and M. Holzmann, Nonlinear network description for many-body quantum systems in continuous space, Phys. Rev. Lett. 120, 205302 (2018).
  * Chen  _et al._ [2022] Z. Chen, D. Luo, K. Hu, and B. K. Clark, Simulating 2+1d lattice quantum electrodynamics at finite density with neural flow wavefunctions (2022), arXiv:2212.06835 [hep-lat] .
  * Pinkus [1999] A. Pinkus, Approximation theory of the mlp model in neural networks, Acta Numerica 8, 143–195 (1999).
  * Dagotto  _et al._ [1992] E. Dagotto, A. Moreo, F. Ortolani, D. Poilblanc, and J. Riera, Static and dynamical properties of doped hubbard clusters, Phys. Rev. B 45, 10741 (1992).
  * Kochkov and Clark [2018] D. Kochkov and B. K. Clark, Variational optimization in the ai era: Computational graph states and supervised wave-function optimization (2018), arXiv:1811.12423 [cond-mat.str-el] .
  * Kingma and Ba [2017] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization (2017), arXiv:1412.6980 [cs.LG] .
  * White [1992] S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).
  * Verstraete  _et al._ [2006a] F. Verstraete, M. M. Wolf, D. Perez-Garcia, and J. I. Cirac, Criticality, the area law, and the computational power of projected entangled pair states, Phys. Rev. Lett. 96, 220601 (2006a).
  * Vidal [2007] G. Vidal, Entanglement renormalization, Phys. Rev. Lett. 99, 220405 (2007).
  * Orús [2019] R. Orús, Tensor networks for complex quantum systems (2019).
  * Verstraete  _et al._ [2006b] F. Verstraete, M. M. Wolf, D. Perez-Garcia, and J. I. Cirac, Criticality, the area law, and the computational power of projected entangled pair states, Phys. Rev. Lett. 96, 220601 (2006b).
  * Orús [2014] R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Annals of Physics 349, 117 (2014).
  * Vidal [2008] G. Vidal, Class of quantum many-body states that can be efficiently simulated, Phys. Rev. Lett. 101, 110501 (2008).
  * Deng  _et al._ [2017] D.-L. Deng, X. Li, and S. Das Sarma, Quantum entanglement in neural network states, Phys. Rev. X 7, 021021 (2017).
  * Passetti  _et al._ [2023] G. Passetti, D. Hofmann, P. Neitemeier, L. Grunwald, M. A. Sentef, and D. M. Kennes, Can neural quantum states learn volume-law ground states?, Phys. Rev. Lett. 131, 036502 (2023).
  * Denis  _et al._ [2023] Z. Denis, A. Sinibaldi, and G. Carleo, Comment on ”can neural quantum states learn volume-law ground states?” (2023), arXiv:2309.11534 [quant-ph] .
  * Trigueros  _et al._ [2023] F. B. Trigueros, T. Mendes-Santos, and M. Heyl, Mean-field theories are simple for neural quantum states (2023), arXiv:2308.10934 [quant-ph] .
  * Sharir  _et al._ [2022] O. Sharir, A. Shashua, and G. Carleo, Neural tensor contractions and the expressive power of deep neural quantum states, Phys. Rev. B 106, 205136 (2022).
  * Clark [2018] S. R. Clark, Unifying neural-network quantum states and correlator product states via tensor networks, Journal of Physics A: Mathematical and Theoretical 51, 10.1088/1751-8121/aaaaf2 (2018).


Oscillations in Modified Combinants of Hadronic Multiplicity Distributions

P. Agarwal1$\star$, H. W. Ang1, Z. Ong1, A. H. Chan1 and C. H. Oh1

1 Department of Physics, National University of Singapore

*<EMAIL_ADDRESS>
50th International Symposium on Multiparticle Dynamics
(ISMD2021)
12-16 July 2021 10.21468/SciPostPhysProc.?

## Abstract

Oscillations in modified combinants ($C_{j}$s) have been of interest to
multiparticle production mechanisms since the 1990s [1]. Recently, there has
been a discussion on how these oscillations can be reproduced by compounding a
binomial distribution with a negative binomial distribution [2, 3]. In this
work, we explore a stochastic branching model based on a simple interaction
term $\lambda\overline{\psi}\phi\psi$ for partons and propose a hadronization
scheme to arrive at the final multiplicity distribution. We study the effects
that compounding our model with a binomial distribution has on $C_{j}$s and
explore its physical implications. We find that there is a significant
difference in the oscillations in $C_{j}$s between high energy $pp$ and
$p\bar{p}$ scattering that our model can reproduce.

##  1 Introduction

In high energy particle production, one of the characteristic measurements is
that of the multiplicity distribution. It assigns a probability $P(n)$ to
producing $n$ hadrons for a given interaction. Over the years, many
interesting features of these distributions have been studied [4]. Modified
combinants can be defined as the coefficients of recursion for the
multiplicity distribution. That is

$(n+1)P(n+1)=\left<n\right>\sum_{j=0}^{n}\ C_{j}P(n-j).$ (1)

These quantities are known to exhibit oscillations w.r.t. their rank $j$. The
plots (1) shown showcase these oscillations across three variables:
pseudorapidity window, collision type and centre of mass energy. By looking at
the plots themselves, there are no obvious candidate explanations for what
causes these oscillations or what physical information might be gleaned from
their period or amplitude.

Figure 1: Oscillations in Modified Combinants of hadronic multiplicity
distributions with rank.

Herein lies the motivation of this work: we propose a physical model of
multiple production that has clear physical meaning to each of its variables.
By reproducing these oscillations, we are able to shed some light on what
exactly causes this phenomena. We shall focus here mostly on the difference
across collision types.

##  2 Multiple Cluster Formation

It was pointed out recently by Wilk et al. [2, 5] that clan/cluster models
result in oscillations in $C_{j}$ if the cluster production is explained using
a binomial distribution. In such a model, the generating function of the
"daughter" distribution that produces particles within a cluster is compounded
with the generating function of the "parent" distribution describing cluster
production itself. The final multiplicity can be calculated using chain rule
on the generating function, as summarised in the Faà di Bruno’s formula:

$P(n)=\frac{1}{n!}\sum_{k=1}^{n}P_{\text{BD}}(r)\cdot B_{n,r}(p_{1},2!\
p_{2},...,(n-r+1)!\ p_{n-r+1})$ (2)

where $B_{n,r}$ are Bell’s polynomials and $p_{n}$ is the daughter
distribution.

##  3 Stochastic Branching Models

Hadron production is a QCD phenomena which has inherent non-perturbative
effects. One way to circumvent this problem is to employ branching models
derived from perturbation theories as the resulting stochastic equations are
independent of the coupling [6]. We use one such branching model with the
interaction term $\lambda\overline{\psi}\phi\psi$ for our daughter
distribution [7]. There are three fundamental processes possible dictated by
the same vertex: $\psi$ Bremsstrahlung, $\overline{\psi}$ Bremsstrahlung and
$\overline{\psi}\psi$ pair production.

Figure 2: $\psi$ Bremsstrahlung, $\overline{\psi}$ Bremsstrahlung and
$\overline{\psi}\psi$ pair production.

The evolution equation of this system can be easily derived. We assign a
probability to each of the three kinds of branching processes: $a$ for $\psi$
Bremsstrahlung, $b$ for $\overline{\psi}$ Bremsstrahlung and $c$ for
$\overline{\psi}\psi$ pair production. A state with $n$ "effective" quarks,
$m$ "effective" antiquarks and $o$ "effective" scalar glueballs has
probability given by $P_{mno}$ and evolves as

$\displaystyle\frac{\text{d}P_{mno}(t)}{\text{d}t}=$ $\displaystyle\
na(P_{mn,o-1}-P_{mno})+mb(P_{mn,o-1}-P_{mno})$
$\displaystyle+c((o+1)\,P_{m-1,n-1,o+1}-o\,P_{mno})\text{.}$ (3)

##  4 Hadronisation & Modes of Interaction

###  4.1 Hadronisation

Having generated our partons via branching, we propose a hadronisation scheme
that translates this partonic distribution to a final state hadronic
distribution. Suppose we have parton multiplicity distributions $p(n)$ for
quarks and $\bar{p}(n)$ for antiquarks.

Assumptions:

  * •

We ignore "exotic" particles such as tetra/penta/hexa-quark structures and
focus on quark-antiquark mesons and tri quark/antiquark baryons.

  * •

For a given partonic state, all possible physical states have a priori the
same probability of occurring. For example, if we have 3 quarks and 3
antiquarks, the final state is equally likely to contain either 3 mesons or 1
baryon-antibaryon pair.

Under these assumptions, the hadron multiplicity distribution $P(n)$ is given
by:

$P(n)=\sum_{i=0}^{n}\sum_{j=0}^{n-i}\ \frac{1}{m_{3(n-i-j)+i,3j+i}}\
p(3(n-i-j)+i)\ \bar{p}(3j+i).$ (4)

Here, $m_{n,\bar{n}}$ are the number of final possible hadronic states for
given $n$ quarks and $\bar{n}$ antiquarks. It is symmetric matrix calculated
for $n>\bar{n}$ using

$m_{n,\bar{n}}=1+\left\lfloor\frac{\bar{n}}{3}\right\rfloor.$

###  4.2 Modes of Interaction

We assume a 3-(anti)quark structure for (anti)proton (3). When a $p\bar{p}$
interaction happens, we get a fast moving di(anti)quark and a held-back
(anti)quark for each (anti)proton. For hadronisation to result in a final
state colour singlet, the held-back quark-antiquark chain hadronise together
while the fast moving diquark-diantiquark chain hadronise together [8].

Figure 3: Hadron chains form to keep the final state a colour singlet. (Figure
borrowed from [8])

The story is different for $pp$ interactions. The held-back quark pair no
longer form a colour singlet and therefore each hadronises with a fast moving
diquark instead. This is a fundamental difference in the nature of interaction
between $pp$ and $p\bar{p}$ that we will be probing.

##  5 Behaviour of $\bf{C_{j}}$ in our Model

There are two main ways in our model to distinguish between $pp$ and
$p\bar{p}$. First is the difference in branching probabilities of quarks ($a$)
and antiquarks (b) while the second is the difference between the
hadronisation mechanisms of the two interactions. The second of these can be
studied simply by constraining the model such that $a=b$. Sample results are
shown in Fig. 4.

Figure 4: Oscillations in $C_{j}$ from our model. The plot on the left is for
$a\neq b$ while the one on the right is for $a=b$

The effect of $a\neq b$ is brought out by imposing $b\sim 2a$. Notice that the
period of oscillation is highly dependent on the branching probabilities of
our partons. Another interesting feature of our model is the seemingly linear
dependence of the amplitude of oscillations on the maximum number of clusters
allowed. The period of oscillation also depends on the probability of cluster
formation. These details are currently being investigated.

##  6 Conclusion

The phenomena of oscillations in $C_{j}$s is an interesting mystery of
multiplicity phenomenology. By studying a physical model of these
oscillations, there is significant information to be drawn out from the period
and amplitude of these oscillations. We find that there is a significant
difference between $pp$ and $p\bar{p}$ which our model is able to capture.
Upon preliminary fitting, our model is able to describe the multiplicity
distribution across $pp$ and $p\bar{p}$ interactions with $\chi^{2}/d.o.f.$
better than 1.5 in all cases tested (not shown here). The period of
oscillation depends on the branching probability of partons as well as the
probability of cluster formation.

## Acknowledgements

#### Funding information

PA, HWA and ZO are supported by the National University of Singapore Research
Scholarship.

## References

  * [1] I. M. Dremin,  _QCD and hybrid NBD on oscillating moments of multiplicity distributions in lepton and hadron initiated reactions_ ,  Phys. Atom. Nucl. 68, 119 (2005),  10.1134/1.1858565,  hep-ph/0401013.
  * [2] G. Wilk and Z. Wlodarczyk,  _Counting distributions from the perspective of combinants_ ,  Int. J. Mod. Phys A 36(13), 2150072 (2021).
  * [3] H. W. Ang, A. H. Chan, M. Ghaffar, M. Rybczyński, G. Wilk and Z. Włodarczyk,  _A look at multiparticle production via modified combinants_ ,  Eur. Phys. J. A 56(4), 117 (2020),  10.1140/epja/s10050-020-00140-w,  1908.11062.
  * [4] J. F. Grosse-Oetringhaus and K. Reygers,  _Charged-Particle Multiplicity in Proton-Proton Collisions_ ,  J. Phys. G 37, 083001 (2010),  10.1088/0954-3899/37/8/083001,  0912.0023.
  * [5] M. Rybczyński, G. Wilk and Z. Włodarczyk,  _A look at multiplicity distributions via modified combinants_ ,  Ukr. J. Phys. 64(8), 738 (2019),  10.15407/ujpe64.8.738,  1906.11531.
  * [6] R. C. Hwa,  _BRANCHING PROCESSES IN MULTIPARTICLE PRODUCTION_ , pp. 556–613,  10.1142/9789814503259_0014,  https://www.worldscientific.com/doi/pdf/10.1142/9789814503259_0014.
  * [7] P. Agarwal, H.W. Ang, Z. Ong, A.H. Chan and C.H. Oh,  _On multiparticle production in very high energy scattering_ ,  EPJ Web Conf. 240, 07001 (2020),  10.1051/epjconf/202024007001.
  * [8] A. Capella, U. Sukhatme, C.-I. Tan and J. Tran Thanh Van,  _Dual parton model_ ,  Physics Reports 236(4), 225 (1994),  https://doi.org/10.1016/0370-1573(94)90064-7.