# Evolution of Cosmic Voids in the Schrödinger-Poisson formalism

Aoibhinn Gallagher∗<EMAIL_ADDRESS>Peter Coles  Department of
Theoretical Physics, Science Building, Maynooth University, Maynooth, Co.
Kildare, Ireland.

###### Abstract

We investigate the evolution of cosmic voids in the Schrödinger-Poisson
formalism, finding wave-mechanical solutions for the dynamics in a standard
cosmological background with appropriate boundary conditions. We compare the
results in this model to those obtained using the Zel’dovich approximation. We
discuss the advantages of studying voids in general and the advantages of
Schrödinger-Poisson description over other approaches. In particular
emphasizing the utility of the free-particle approximation. We also discuss a
dimensionless number, similar to the Reynolds number, for this system which
allows our void solutions to be scaled to systems of different physical
dimensions.

##  1 Introduction

Over the past few decades, galaxy redshift surveys, exemplified by those
resulting from the Sloan Digital Sky Survey described by Eisenstein et al.
(2011), have revealed that the present-day Universe displays a rich
hierarchical pattern of structure - the ”Cosmic Web” - that encompasses a vast
range of length scales, from individual galaxies to clusters and filaments
surrounding enormous voids. On the other hand, observations of the cosmic
microwave background, most recently from Planck Collaboration (2020), suggest
that the early Universe was almost homogeneous, with only slight temperature
ﬂuctuations seen in the cosmic microwave background radiation. Models of
structure formation link these observations through the eﬀect of gravity,
relying on the fact that small initially over-dense regions accrete additional
matter as the universe expands (a mechanism known as gravitational
instability). The growth of density perturbations via gravitational
instability is well understood in the linear regime, but the non-linear regime
is much more complicated and generally not amenable to analytic solution.
Numerical N-body simulations such as those discussed in Springel et al. (2005)
have led the way towards a fuller understanding of the structure formation
process, but although such calculations have been priceless in establishing
quantitative predictions of the large-scale structure expected to arise in a
particular cosmology, it remains important to develop as full an understanding
as possible using analytical models. This is not only because simulation is
not the same as comprehension but also because $N$–body methods are
computational expensive which makes it difficult to explore a large parameter
space using them alone.

To complement the ”brute force” approach of purely computational techniques, a
variety of analytic techniques has been developed for describing various
aspects of the evolution of density perturbations to various levels of
approximation; for a good review of many of these approaches,see Bernardeau et
al. (2002). None of these capture all the details of the structure formation
process, but each can provide useful information in particular circumstances.

In this paper we study a particular formulation of the gravitational
instability problem - the Schrödinger-Poisson (SP) approach (explained below)
- which is based on wave mechanics rather than fluid mechanics. This approach
was first proposed in the cosmological context by Widrow & Kaiser (1993), and
has subsequently generated a great deal of interest; see e.g. Coles (2002),
Szapudi & Kaiser (2003), Coles & Spencer (2003), Short & Coles (2006a), Short
& Coles (2006b) and Johnston et al. (2010). It turns out that this approach,
while initially developed as an approximate method for handling the evolution
of cold dark matter, is also applicable to a scenarios in which the behaviour
of dark matter is inherently wave-like rather than particle-like, such as if
the dark matter is a very light bosonic particle, perhaps in the form of a
condensate, or some other form of “fuzzy” dark matter; see for example Brook &
Coles (2009), Schwabe et al. (2020), Dentler et al. (2021), Hui (2021) and
references therein.

The behaviour of cosmic voids provides an interesting test case for the SP
system that complements the case of spherical collapse studied by Johnston et
al. (2010). While gravitational collapse tends to amplify asymmetry, Icke
(1984) proved a “Bubble Theorem” according to which an initially asymmetric
void becomes more spherical as the Universe evolves; see also numerical
computations by Bertschinger (1985). The study of a spherically symmetric void
is therefore arguably more realistic than the collapse of a spherical over-
density. Moreover, because matter tends to become compressed in sheets or
filaments as the Cosmic Web develops, voids actually dominate space in terms
of volume fraction. As well as these physical aspects of void evolution there
is an important statistical property shown by White (1979) that voids contain
information about the entire hierarchy of $n$-point correlation functions at
all orders. For these reasons and others there has been a great deal of recent
interest in the behaviour of cosmic voids, see e.g. Wojtak et al. (2016),
Pisani et al. (2020) and Minoguchi et al. (2021).

Pulling together these chains of thought, in this paper we apply the
Schrödinger-Poisson system to the case of isolated spherical voids. Although
of course in a realistic context void expansion would be hindered by
neighbouring structures, it is important to obtain as full an understanding as
possible of the behaviour of the simplest case within the SP system. The
layout of the paper is as follows: in Section 2, we outline the background
cosmology needed to build the results of this paper; in Section 3, we outline
the basics of the SP system; in Sections 4 and 5, we look directly at the void
and apply what was covered in previous sections to solve for the dynamics of
the void. In Section 6 we discuss scaling properties of the solutions in terms
of an effective Reynolds number. In Section 7 we present our conclusions and
suggest future applications of our results and the approach that leads to
them.

##  2 Background Cosmology

The current standard model of cosmology assumes the Cosmological Principle,
according to which the universe is assumed to be homogeneous and isotropic, at
least on large scales. Space-times consistent with this are described by the
Robertson-Walker metric:

$ds^{2}=c^{2}dt^{2}-a^{2}(t)\left(\frac{dr^{2}}{1-\kappa
r^{2}}+r^{2}d\theta^{2}+r^{2}\sin^{2}(\theta)d\phi^{2}\right),$ (1)

where $\kappa$ is the spatial curvature: $\kappa=0$ represents flat space;
$\kappa=+1$ represents constant positively curved space (closed universe); and
$\kappa=-1$ represents constant negatively curved space (open universe). The
time coordinate $t$ is cosmological proper time and $a(t)$ is the cosmic scale
factor. The time evolution of the cosmic scale factor, $a(t)$, is determined
via Einstein’s gravitational field equations through the Friedman equation,

$3\left(\frac{\dot{a}}{a}\right)^{2}=8\pi G\rho-\frac{3\kappa
c^{2}}{a^{2}}+\Lambda$ (2)

the deceleration equation,

$\frac{\ddot{a}}{a}=-\frac{4\pi
G}{3}\left(\rho+3\frac{p}{c^{2}}\right)+\frac{\Lambda}{3}$ (3)

and the density-pressure relation,

$\dot{\rho}=-3\frac{\dot{a}}{a}\left(\rho+\frac{p}{c^{2}}\right)$ (4)

where the dots denote derivatives with respect to cosmological proper time t,
thus describing the global expansion or contraction of the universe. These
models can be further parameterised by the Hubble parameter $H=\dot{a}/a$ and
the density parameter $\Omega=8\pi G\rho/3H^{2}$. As usual the present epoch
is defined by $t=t_{0}$, when $H=H_{0}$ and $\Omega=\Omega_{0}$.

Throughout this paper we assume Newtonian gravity, since the scale of the
perturbations are much smaller than the effective cosmological horizon
$d_{H}=c/H$. We also do all calculations in co-moving coordinates,

$\textbf{x}\equiv\textbf{r}/a(t),$ (5)

where x is the co-moving spatial coordinate (x is fixed in the frame of the
Hubble expansion), r is the spatial coordinate in our reference frame, and
$a(t)$ is the cosmic scale factor.

###  2.1 The Zel’dovich Approximation

A simple and elegant approximation to the non-linear stage of gravitational
evolution is the Zel’dovich approximation, introduced by Zel’Dovich (1970).
The initial density distribution considered homogeneous and collisionless. The
initial (Lagrangian) coordinates are q, then the Eulerian coordinates at time
t are given by

$\textbf{r}(\textbf{q},t)=a(t)[\textbf{q}+b(t)\textbf{s}(\textbf{q})],$ (6)

where $a(t)$ is the cosmic scale factor, $b(t)$ describes the evolution of a
perturbation in the linear regime, and is defined by

$\ddot{b}+2\Big{(}\frac{\dot{a}}{a}\Big{)}\dot{b}-4\pi G\rho b=0.$ (7)

In a flat matter-dominated universe $b\propto t^{2/3}$.
$\textbf{s}(\textbf{q})$ here is the displacement term, described by the
initial velocity potential $\phi_{0}$ by,

$\nabla\phi_{0}(\textbf{q})=\textbf{s}(\textbf{q}).$ (8)

Mass conservation requires
$\rho(\textbf{r},t)d\textbf{r}=\rho_{0}d\textbf{q}$, so the density field as a
function of Lagrangian coordinates is

$\rho(\textbf{q},t)=\rho_{0}\left|\frac{\partial\textbf{q}}{\partial\textbf{r}}\right|=\overline{\rho}\left|\delta_{ij}-b(t)\frac{\partial
s_{i}}{\partial q_{j}}\right|^{-1}$ (9)

where $\overline{\rho}=(a_{0}/a)^{3}\rho_{0}$. Here
$\frac{\partial\textbf{r}}{\partial\textbf{q}}$ is the deformation tensor,
which accounts for the gravitational evolution of the density field. Due to
the nature of $\textbf{s}(\textbf{q})$, the deformation tensor is a real
symmetric matrix with eigenvectors that define a set of three principal
(orthogonal) axes. Rewriting Equation (9) in terms of the eigenvalues
$-\alpha(\textbf{q})$, $-\beta(\textbf{q})$ and $-\gamma(\textbf{q})$, we get

$\rho(\textbf{q},t)=\frac{\overline{\rho}}{[1-b(t)\alpha(\textbf{q})][1-b(t)\beta(\textbf{q})][1-b(t)\gamma(\textbf{q})]}$
(10)

If we use comoving coordinates (Equation (5)) in 1 dimension, we get the much
simpler expression for $\rho$

$\rho(q,t)=\frac{\rho_{0}}{1-b(t)\frac{\partial s(q)}{\partial q}}$ (11)

The Zel’dovich approximation is an excellent tool for studying structure
formation and is highly accurate up to the point of shell crossing Shandarin &
Zeldovich (1989), at which point it breaks down. To combat this we need a
model that goes beyond shell crossing, this is one reason why we look to the
SP model. Other methods for studying multi-stream flow can be seen in Buchert
& Domínguez (1998) and Gough & Uhlemann (2022).

##  3 Schrödinger Poisson Regime

We start with the Newtonian equations for a self-gravitating perfect fluid.The
Euler Equation

$\frac{\partial\textbf{v}}{\partial
t}+(\textbf{v}\cdot\nabla)\textbf{v}+\frac{1}{\rho}\nabla p+\nabla\phi=0,$
(12)

the Continuity Equation

$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\textbf{v})=0,$ (13)

and the Poisson Equation

$\nabla^{2}V=4\pi G\rho.$ (14)

Here v, $p$ and $\rho$ are the fluid velocity, pressure and density,
respectively, and $V$ here is the gravitational potential. We simplify our
treatment by assuming a cold (pressureless), irrotational fluid moving under
the influence of a general potential $V$. This means that Equation (12) now
becomes the much simpler Bernoulli equation

$\frac{\partial\phi}{\partial t}+\frac{1}{2}(\nabla\phi)^{2}=-V,$ (15)

where $\textbf{v}=\nabla\phi,$ where $\phi$ is the velocity potential.

We now make the Madelung Transformation $\psi=\alpha e^{i\phi/\nu}$, where
$\rho=\psi\psi^{*}=\alpha^{2}$. This gives

$i\nu\frac{\partial\psi}{\partial
t}=-\frac{\nu^{2}}{2}\nabla^{2}\psi+V\psi+P\psi$ (16)

where

$P=\frac{\nu^{2}}{2}\frac{\nabla^{2}\alpha}{\alpha},$ (17)

from Equations (15) and (13). Equation (16) is Schrödinger Equation, with
potential $V$, $\nu$ acting as $\hbar$ with the addition of non-linear term
$P$. Although $\psi$ is governed by the same equation as the evolution of a
single-particle wave function, that is not how it should be interpreted. In
particular, $|\psi|^{2}$ represents a physical density not a probability
density, and its evolution is completely unitary - there is nothing like the
wave-function collapse that occurs in standard quantum mechanics in this
system.

It’s important to note a crucial advantage of this description, namely that
because $\rho=|\psi|^{2}$ the condition that $\rho\geq 0$ is automatically
enforced if one applies, e.g., perturbation theory to $\psi$. This is not the
case for approaches based on standard Eulerian perturbation theory applied to
$\delta=(\rho-\rho_{0})/\rho_{0}$ which predict $\rho<0$ when $\delta<-1$ and
are therefore very unsuitable for describing voids. Note also that because the
wavefunction describes a delocalized particle there are no singularities
analogous to the caustics that form in the Zel’dovich approximation.

In the cosmological context we take $V$ to be the gravitational potential
determined via the Poisson Equation and we get a system of coupled Partial
Differential Equations (PDEs). We set $\nu$ to be an adjustable parameter so
we have some control over the frequency of any oscillatory solutions that
could arise from this system. $\nu$ has dimensions such that $\phi/\nu$ is
dimensionless, $\phi$ is a velocity potential so has dimensions $L^{2}T^{-1}$,
these are also the dimensions of kinematic viscosity, so $\nu$ can be thought
of as a viscosity parameter for the model; we return to this later.

The original intentions of this model were to find a fluid interpretation of a
quantum system, Widrow & Kaiser (1993) first proposed this be used to simulate
cold dark matter (CDM), and it has been used predominantly in a cosmological
context since. The quantum nature of this model is also an advantage as it
gives us a new perspective on the behaviour of large scale structure. This is
especially useful when one seeks to look beyond shell crossing, as this model
can handle multiple streams.

##  4 Cosmological Solutions

Analytic solutions to the SP system can be obtained straightforwardly for the
simplest case of a homogeneous and isotropic fluid. For this $P\equiv 0$,
since $|\psi|^{2}=\alpha^{2}=\rho$ is defined to be homogeneous. We can start
our study of a void by using solutions for an open and flat universe for the
interior and exterior of the void in order to make a snapshot of a simple void
configuration as follows.

###  4.1 Spatially flat universe

For the following we assume a flat universe (with $\kappa=0$). We begin here
by taking the well known cosmological result for the evolution of the density
field under a globally defined cosmic time, $t$, given by

$\rho(t)\equiv|\psi|^{2}=\frac{\Lambda c^{2}}{8\pi
G}\text{cosech}^{2}\Big{(}\frac{3}{2}\sqrt{\frac{\Lambda c^{2}}{3}}t\Big{)}.$
(18)

Under spherical symmetry the Laplacian operator is simply

$\nabla^{2}\equiv\frac{\partial^{2}}{\partial
r^{2}}+\frac{2}{r}\frac{\partial}{\partial r}.$

This density evolution obtained previously can be substituted into Equation
(14) to obtain an equation for the gravitational potential $V$. This allows
the calculation of $\phi$ by substitution into Equation (16), resulting in the
wave-function $\psi$, as shown in Johnston et al. (2010):

$V=\frac{\Lambda c^{2}}{12}(\text{coth}^{2}(\lambda t)-3)r^{2}$ (19)
$\psi=\sqrt{\frac{\Lambda c^{2}}{8\pi G}}\mathrm{cosech}\,(\lambda
t)\exp\Big{(}\frac{i}{v}\sqrt{\frac{\Lambda c^{2}}{12}}\coth\left[\lambda
t\right]r^{2}\Big{)}$ (20)

where $\lambda\equiv\frac{3}{2}\sqrt{\Lambda c^{2}/3}$. We therefore also find

$\phi=\sqrt{\frac{\Lambda c^{2}}{12}}\coth(\lambda t)r^{2}.$ (21)

Also following Johnston et al. (2010), a series expansion of $\rho=\alpha^{2}$
and $\phi/r^{2}$ about $\Lambda=0$ provides the well known results of
$\rho\propto t^{2}$ and $V=r^{2}/9t^{2}$ for this particular case.

###  4.2 Open universe

Working in a similar way to Section 4.1, we find solutions for a universe with
negative curvature, $\kappa=-1$, an open universe, to be

$\rho=\frac{\gamma}{A^{3}(\cosh(\eta)-1)^{3}},$ (22)
$V=\left\\{\frac{G\pi\gamma}{12A^{3}\left[\sinh\left(\frac{\eta}{2}\right)\right]^{3}}\right\\}r^{2},$
(23)

$\psi=\frac{C}{\sqrt{2A}\left[\sinh\left(\frac{\eta}{2}\right)\right]^{3}}\exp\left[\frac{i}{v}C\coth\left(\frac{\eta}{2}\right)\mathrm{cosech}\,\left(\frac{\eta}{2}\right)r^{2}\right],$
(24)

where $A\equiv 4\pi G\rho_{0}/3$, $C=\gamma^{1/2}/4A$, $\mu=3M/4\pi$ and
$\gamma$ is a constant. For a simpler derivation we have converted to
conformal time by using $ad\eta=dt$ where $a$ is the scale factor of the open
universe.

Having established this simple solution we can now proceed to study the
dynamics of a void, taking into account the fact that the density surrounding
the void will not remain homogeneous as the interior expands.

##  5 Modelling an Underdensity

In this Section we look at simple cases of a void growing under its own
gravity. We will use the solutions found in Section 4 to form an analytic
model of a void, in analogous fashion to the spherical collapse model studied
by Johnston et al. (2010). We will also discuss the results of the numerical
solutions to the SP system of equations for a simple one dimensional void
model, both ‘free-particle’ and including the gravitational potential in the
evolution. These solutions will be compared to the Zel’dovich approximation
for the same initial conditions, which is exact up to the moment of shell
crossing, Shandarin & Zeldovich (1989). We see that shell crossing happens
rather quickly in the case of a void though the formation of multi-stream
regions forming at the edge of the void. This means that the Zel’dovich
approximation applies for a shorter portion of the evolution in the void case
than in the collapse case.

###  5.1 Analytic solutions to the void problem

Using the results from Section 4 and following the same method as described by
Johnston et al. (2010), we can construct a simple model of a void. We
construct the model of the void with an inner region and an outer fluid
separated by a shell. The inner fluid is modeled by the open universe
solutions and the outer fluid is modeled by the flat universe solutions. Both
fluids are homogeneous but the fluid in the inner region is less dense than
the outer fluid so expands more quickly, pushing into the surrounding medium
and evacuating the void. This creates a small region on the cusp where the
density is higher than the outer region. An illustration of how this would
look is shown in Figure 1. This figure shows the inside of the void in white,
this is the area that gets emptied out. The yellow area is the shell of the
void, this is where the denser region is created at the edge of the void. The
pink area is the the rest of the universe which extends to infinity.

Figure 1: An illustration of the void, shell and flat space in white, yellow
and pink respectively. It shows $R_{o}$ as the radius to the inner shell
radius and $R_{f}$ as the radius of the outer shell.

Using the cosmological solutions shown in Section 4, and following the same
procedure as Johnston et al. (2010), we compute an expression for the
expansion of the radii of both sections of the void, in conformal time $\eta$.
We have $R_{o}$ for the open-universe section of the void and $R_{f}$ for flat
universe section of the void.

$R_{o}=\mu^{1/3}\frac{2A}{\gamma^{1/3}}(\sinh^{2}(\eta/2)),$ (25)

and

$R_{f}=B\left(\sinh(\eta)-\eta\right)^{2/3},$ (26)

where $B$ is a free parameter corresponding to the initial conformal time
$\eta_{0}$ when the configuration is set up, at which point $R_{o}=R_{f}$. The
evolution is illustrated by Figure 2 which shows the evolution of each of
these radii with respect to conformal time with an arbitrary scaling on the
vertical axis just to show that the inner region of the void expands much more
quickly than the outer region, and an overlapping region develops after some
time as the inner void expands into a region that was previously in the space
exterior to the void.

Figure 2: Evolution of the radii of both the inner section and outer section
of the void, with respect to conformal time $\eta$.

###  5.2 The ‘free particle’ void

In this Section we will present numerical solutions of the void evolution in
one dimension using the SP system. We start by studying the approximate case
of a void where the gravitational potential is calculated at each time-step,
but not included in the Schrödinger equation for the evolution of the void. In
other words, the gravitational potential is updated as the fluid moves but the
effect of this evolution on the fluid motion is ignored. This has the effect
that the fluid moves as a collection of ‘free particles’, a shortcut that has
been shown to be remarkably accurate for some applications; see Short & Coles
(2006a). The Schrödinger - Poisson system in this approach simplifies to

$\displaystyle\frac{\partial\psi}{\partial b}$
$\displaystyle=-\frac{\nu^{2}}{2}\frac{\partial^{2}\psi}{\partial x^{2}}$ (27)
$\displaystyle\frac{\partial^{2}V}{\partial x^{2}}$ $\displaystyle=4\pi
G|\psi|^{2}.$ (28)

These equations thus decouple, making the computation much simpler. We begin
with a compensated void, shown in Figure 3. To get this shape we use a
smoothed $\tanh$ function.

Figure 3: Initial fluid density field for compensated void in one dimension.

We need to set up an initial velocity field in which the fluid particles are
moving away from the centre of the void, quickest at the centre and more
slowly towards the edges of the void. This behaviour can be seen in Figure 4a.
To get $\phi$, as shown in Figure 4b, we solve $\nabla v=\phi$. This solution
for $\phi$ is only determined up to a constant, which we fix by setting the
potential at infinity to be zero (the edge of the box).

We performed these calculations using periodic boundary conditions. These
boundary conditions work quite well with these initial conditions as we have
constructed them in such away that both the density field and the velocity
potential go to zero at the boundary. The boundary is also far enough away
from any structure so do not generate any significant cross-boundary
interaction as the boundaries are effectively infinitely far away. The
Schrödinger-Poisson solution is seen in Figure 5; this solution is analogous
to the Zel’dovich approximation if we work in comoving coordinates x, defined
by Equation (5), and time coordinate $b$, as defined by Equation (7), so these
two solutions are comparable.

(a) velocity field

(b) velocity potential

Figure 4: Initial velocity field and its potential.

The Zel’dovich approximation is calculated using Equation (11), for $\rho_{0}$
in Figure 3 and $s(q)=\nabla\phi_{0}(q)$, where $\phi_{0}$ is that in Figure
4b. Shell-crossing occurs soon after $b=1$, where multi-streaming begins and
the fluid density field is no longer well defined under Zel’dovich, as it
cannot account for multi-streaming. Comparing Figures 5 and 6, we see the
evolution of the void at three different time-steps in the Schrödinger-Poisson
and two different time-steps in the Zel’dovich approximation. At early time-
step $b=0.2$, we see the void begins to expand slightly and a small amount of
matter accumulates at the boundary of the void, as seen in Section 5.1. Just
before shell-crossing we see that the two models deviate, with Zel’dovich void
boundaries shooting off to infinity as it struggles to account for multi-
streaming, and the SP void is further expanded with shorter peaks. In the
final time-step $b=1.2$, we see that the interference effect caused by multi-
streaming is in full effect as the void continues to expand. This interference
effect occurs after shell crossing, and is also seen in gravitational collapse
described by SP dynamics in, e.g., Coles & Spencer (2003) Johnston et al.
(2010) and Uhlemann et al. (2019).

(a)

(b)

(c)

Figure 5: Evolution of the void computed numerically shown at three time-
steps. The first is an early time-step, the second is just before shell-
crossing, the third is post shell-crossing.

(a)

(b)

Figure 6: The time-steps shown here are the same as the first two Figure 5:
the first is an early time-step, the second is just before shell crossing.
There is no post shell-crossing plot in this case as a singularity in the
density forms at shell-crossing.

###  5.3 Void Expansion Speed

It can be seen in Figure 5 that the edges of the void move outwards in
comoving coordinates as time goes on. Since the peaks in $\rho$ are sharply
defined at least until the strong interference develops after shell-crossing,
their positions can be used to define an expansion speed that can be extracted
easily from the numerical calculations. Doing this, we found the group
velocity of this peak to be generally constant in time (as defined by $b$),
with some jumps once the interference pattern starts to develop. This can be
seen in Figure 7, which shows the position of the highest peak at the cusp of
the void changes with respect to time. The resulting plot is roughly linear in
the coordinates used, which indicates a free expansion of the void region
discussed by e.g. Bertschinger (1985) and Coles & Barrow (1990)). At very late
times one expects a void embedded in a baryonic fluid to enter an adiabatic
phase described by the Sedov solution for blast waves (Sedov (1946)), but the
physical behaviour of the quantum fluid described here would probably be quite
different and in any case we only look at the initial phase of the expansion
here.

Figure 7: Positions of the peak of the mass at the edge of the void by SP
regime as a function of time; the jumps are caused by the formation of
subsidiary peaks through wave interference.

###  5.4 Including the potential term

In this section we look at including a static potential and a time-varying
potential term in the Schrödinger-Poisson system of equations. For the static
potential, the potential is calculated by Equation (28) from the initial
conditions, and then used as a time-independent constant in

$\frac{\partial\psi}{\partial
b}=-\frac{\nu^{2}}{2}\frac{\partial^{2}\psi}{\partial x^{2}}+V\psi.$ (29)

For the time-varying potential, both Equations (29) and (28) are solved
simultaneously as a coupled system of PDEs for a time-dependent $\psi$ and
$V$. Since the evolution of the gravitational potential is quite small, we did
not expect either of these extensions to have a drastic effect on the results,
but we include them here for completeness. The effect of the potential on the
system is only seen at late time-steps and only has an effect on the
interference pattern seen after shell crossing. These effects can be seen by
comparing Figures 5, 6, 8 and 9.

(a)

(b)

Figure 8: Solution with a static potential, just before shell crossing and
shortly after shell crossing.

(a)

(b)

Figure 9: Solution with a time-dependent potential, just before and just after
shell-crossing.

The initial conditions show in Figures 3 and 4 were used for these
calculations also, so they can be compared. Comparing Figures 5, 8 and 9, we
don’t see much of a difference between the images, although it is worth noting
that for the same time-step in Figures 8 and 9 we see the interference pattern
coming into effect slightly sooner than in Figure 5. This is expected, as we
expect the inclusion of the evolution of the potential to speed up the
gravitational effects of the system. It is also worth noting that the
interference pattern is almost identical between the static and time-varying
potential. This shows us that the potential does not change much over this
time-scale which in turn validates the use of the free-particle approximation.
We do not see a notable difference in the behaviour of the outward peak
velocity, when the potential is included.

##  6 Viscosity in a Schrödinger fluid

It is interesting to remark that the parameter $\nu$ involved in the SP
formalism has the same dimensions as a viscosity. This was discussed for
example in Short & Coles (2006b) in connection with the adhesion model, an
extension of the Zel’dovich approximation; see Gurbatov et al. (1989).

In order to discuss kinematic viscosity in the SP system in more detail one
should compare it not with the Euler Equation (12) but with a Navier-Stokes
Equation

$\frac{\partial\textbf{v}}{\partial
t}+(\textbf{v}\cdot\nabla)\textbf{v}+\frac{1}{\rho}\nabla
p-\eta\frac{\nabla^{2}\textbf{v}}{\rho}+\nabla\phi=0,$ (30)

in which case the kinematic viscosity $\eta$ arises naturally as a consequence
of the quantum potential term discussed above in Section 3, though its
behaviour is however somewhat different from that of a classical viscous
fluid, specifically because the viscosity varies spatially and can even be
negative.

We will defer a fuller discussion of the origin and role of viscosity in the
SP system to future work but, continuing with the simple fluid analogy of the
Schrödinger-Poisson model, one can suggest defining something analogous to a
Reynolds number for this model as

$R=\frac{ul}{\nu},$ (31)

In normal fluid mechanics (e.g. Landau & Lifshitz (1959)), the Reynolds number
is a dimensionless quantity formed from $u$ (a velocity), $l$ (a reasonable
length-scale for the system) and $\nu$ is the kinematic viscosity of the
fluid. We see in Section 3 that $\nu$ is essentially a kinematic viscosity.
Although $\nu$ is well defined, it is not clear how one would define the other
two quantities. There are many velocity quantities present in this model and
it is not exactly clear which one we should choose. There is a momentum
associated with the wave-function, but this is not uniform in space and it is
not clear how one could convert this to a single number. There is also the
velocity of the outward peak of the void, $v_{peak}$, a sort of expansion
velocity. Although this is not well-defined, it is a single number with
dimensions of velocity that captures the behaviour of the void. There is also
the issue of the length-scale to choose for this Reynolds number. It is
obvious to choose the size of the void, however, since the void is expanding,
this length is not constant. For this we have decided to choose the initial
size of the void for our length. This is still somewhat arbitrary as the
initial time and therefore size of the void is arbitrarily chosen.

To test the reliability of this pseudo-Reynolds number, we looked at the
relationships between the quantities used to define it. Since $u$ is the
calculated value $v_{peak}$, we started by fixing $l$, and comparing various
$\nu$ values with the $v_{peak}$ produced by this. As seen in Figure 10a, it
is clear that $\nu$ and $v_{peak}$ are proportional to one another. Then we
fixed $\nu$ and calculated $v_{peak}$ for various $l$ values and as seen in
Figure 10b, it is clear that $l$ and $v_{peak}$ have an inversely proportional
relationship. Due to the nature of how $v_{peak}$ is calculated it is not
possible to fix $u$ and look at the relationship between $\nu$ and $l$, but it
suffices for a quick illustration.

(a) Peak-velocity calculated for fixed $l$ and varying $\nu$.

(b) Peak-velocity calculated for fixed $\nu$ and varying $l$.

Figure 10: Peak-velocity calculated for (a) fixed $l$ and varying $\nu$ and
(b) fixed $\nu$ and varying $l$.

The next obvious step is to compute a Reynolds number for the system described
above. The size of the void can be found by taking the space between the edges
of the void, it is not clear from Figure 3, but the size of the void was set
to be $l=2.0$. The frequency parameter, or viscosity, was set to be
$\nu=0.05$. As seen above the peak velocity was calculated and found to be
$v_{\rm peak}=1.55$. This would give us a Reynolds number $R=62.0$. This is
equivalent to a void $30~{}$Mpc wide, with expansion velocity of $27~{}$ km
s-1 and a viscosity of $4\times 10^{10}~{}$ km${}^{2}s^{-1}$. The actual value
for the effective viscosity is not well known, so we used the value stated and
discussed in Gibson (2006) for illustration.

##  7 Conclusions and Discussion

In this paper we have explored the dynamical evolution of cosmic voids through
the lens of the Schrödinger-Poisson approach. Following Widrow & Kaiser
(1993), we argue that whether or not dark matter does behave quantum-
mechanically, this is an advantageous approximate approach that avoids the
formation of caustics in the Zel’dovich approximation outlined in Section 2.1.
In an era dominated by massive $N$–body simulations there is still space for
analytical approaches such as these because they are much less expensive from
a computational perspective and allow us to explore parameter space much more
quickly.

Combined with other works using the SP model to study gravitational collapse
such as Johnston et al. (2010) and Uhlemann et al. (2019), demonstrates that
this has much greater potential for the evolution of voids than the the
Zel’dovich approximation; see Figure 6. The efficacy of this approach is
further enhanced by the fact that this model has scaling properties so that
computations for voids of different sizes and initial expansion rates can be
computed cheaply.

Obviously isolated voids such as those we have studied here are idealized so
their main function is as a testing ground to enable us to develop insight
into the behaviour of the SP system. In future work we plan to explore the
evolution of more realistic 3D structures. A particularly interesting
potential application of the approximate dynamics deployed here is the problem
of cosmic reconstruction, for example from redshift-space to real space, or
evolving structure back in time from the present epoch to some earlier state,
discussed for example by Brenier et al. (2003).

The above comments apply to our use of the SP dynamics as approximate model
for the behaviour of collisionless cold dark matter. If it turns out that dark
matter is not cold but instead is for example ultra-light axionic matter then
$\nu$ takes on physical relevance as it is inversely proportional to the mass
of the light particle. It is important to stress however that the SP approach
has a much wider range of potential use than this specific scenario.

## Acknowledgments

PC is grateful to William Hunter for discussions during a preliminary
investigation of this topic.

## References

  * Bernardeau et al. (2002) Bernardeau F., Colombi S., Gaztañaga E., Scoccimarro R., 2002, Phys. Rep., 367, 1
  * Bertschinger (1985) Bertschinger E., 1985, ApJS, 58, 1
  * Brenier et al. (2003) Brenier Y., Frisch U., Hénon M., Loeper G., Matarrese S., Mohayaee R., Sobolevskiĭ A., 2003, Monthly Notices of the Royal Astronomical Society, 346, 501
  * Brook & Coles (2009) Brook M. N., Coles P., 2009, Gravitational Stability of Vortices in Bose-Einstein Condensate Dark Matter (arXiv:0902.0605)
  * Buchert & Domínguez (1998) Buchert T., Domínguez A., 1998, A&A, 335, 395
  * Coles (2002) Coles P., 2002, The Wave Mechanics of Large-scale Structure (arXiv:astro-ph/0209576)
  * Coles & Barrow (1990) Coles P., Barrow J. D., 1990, MNRAS, 244, 557
  * Coles & Spencer (2003) Coles P., Spencer K., 2003, Monthly Notices of the Royal Astronomical Society, 342, 176–184
  * Dentler et al. (2021) Dentler M., Marsh D. J. E., Hložek R., Laguë A., Rogers K. K., Grin D., 2021, Fuzzy Dark Matter and the Dark Energy Survey Year 1 Data, doi:10.48550/arXiv.2111.01199, https://arxiv.org/abs/2111.01199
  * Eisenstein et al. (2011) Eisenstein D. J., et al., 2011, AJ, 142, 72
  * Gibson (2006) Gibson C. H., 2006, Journal of Applied Fluid Mechanics
  * Gough & Uhlemann (2022) Gough A., Uhlemann C., 2022, arXiv e-prints, p. arXiv:2206.11918
  * Gurbatov et al. (1989) Gurbatov S. N., Saichev A. I., Shandarin S. F., 1989, MNRAS, 236, 385
  * Hui (2021) Hui L., 2021, ARA&A, 59, 247
  * Icke (1984) Icke V., 1984, MNRAS, 206, 1P
  * Johnston et al. (2010) Johnston R., Lasenby A. N., Hobson M. P., 2010, Monthly Notices of the Royal Astronomical Society, 402, 2491
  * Landau & Lifshitz (1959) Landau L. D., Lifshitz E. M., cop. 1959, Fluid mechanics. A Course of Theorical Physics, Pergamon Press, Oxford
  * Minoguchi et al. (2021) Minoguchi M., Nishizawa A. J., Takeuchi T. T., Sugiyama N., 2021, MNRAS, 503, 2804
  * Pisani et al. (2020) Pisani A., et al., 2020, Astro2020: Decadal Survey on Astronomy and Astrophysics, science white papers, no. 40; Bulletin of the American Astronomical Society, Vol. 51, Issue 3, id. 40 (2019)
  * Planck Collaboration (2020) Planck Collaboration 2020, A&A, 641, A1
  * Schwabe et al. (2020) Schwabe B., Gosenca M., Behrens C., Niemeyer J. C., Easther R., 2020, Physical Review D, 102
  * Sedov (1946) Sedov L. I., 1946, Journal of Applied Mathematics and Mechanics, 10, 241
  * Shandarin & Zeldovich (1989) Shandarin S. F., Zeldovich Y. B., 1989, Rev. Mod. Phys., 61:2
  * Short & Coles (2006a) Short C. J., Coles P., 2006a, J. Cosmology Astropart. Phys, 2006, 012
  * Short & Coles (2006b) Short C. J., Coles P., 2006b, J. Cosmology Astropart. Phys, 2006, 016
  * Springel et al. (2005) Springel V., et al., 2005, Nature, 435, 629
  * Szapudi & Kaiser (2003) Szapudi I., Kaiser N., 2003, ApJ, 583, L1
  * Uhlemann et al. (2019) Uhlemann C., Rampf C., Gosenca M., Hahn O., 2019, Physical Review D, 99
  * White (1979) White S. D. M., 1979, MNRAS, 186, 145
  * Widrow & Kaiser (1993) Widrow L. M., Kaiser N., 1993, ApJ, 416, L71
  * Wojtak et al. (2016) Wojtak R., Powell D., Abel T., 2016, MNRAS, 458, 4431
  * Zel’Dovich (1970) Zel’Dovich Y. B., 1970, A&A, 500, 13


# Complete corrected formula for generating functions of the hypergeometric
distribution

Ken Yamamoto  Faculty of Science, University of the Ryukyus, Nishihara,
Okinawa 903-0213, Japan

###### Abstract

The hypergeometric distribution is a popular distribution, whose properties
have been extensively investigated. Generating functions of this distribution,
such as the probability-generating function, the moment-generating function,
and the characteristic function, are known to involve the Gauss hypergeometric
function. We elucidate that the existing formula of generating functions is
incomplete and provide the corrected formula. In view of the utility and
applicability of the hypergeometric distribution, the exact and correct
foundation is crucially important.

Keywords: hypergeometric distribution; Gauss hypergeometric function;
probability-generating function; moment-generating function

MSC: 60E10; 33C05; 62E15

##  1 Introduction

The hypergeometric distribution is a useful and popular discrete probability
distribution in probability theory and statistics. This distribution describes
the number of successes in a sample drawn from a finite population without
replacement. The hypergeometric distribution has been applied to biology [4],
complex networks [3], and blockchain technology [5]. Furthermore, the
hypergeometric distribution provides the theoretical basis for Fisher’s exact
test [11] and the Mantel–Haenszel test [13] used in medical statistics, and
mark–recapture methods [12] in ecology. A simple example of the hypergeometric
distribution is picking colored balls from an urn [6, 1]. Let an urn contain
$K$ white balls and $N-K$ black balls. We let $X$ denote the number of white
balls in $n$ balls picked randomly. The probability that $k$ balls are white
is given by

$\mathrm{Pr}\left(X=k\right)=\left.\binom{K}{k}\binom{N-K}{n-k}\middle/\binom{N}{n}\right..$
(1)

This is the probability mass function of the hypergeometric distribution. The
support of this distribution, i.e., the set of $k$ with
$\mathrm{Pr}\left(X=k\right)>0$, consists of integers within
$\max\\{0,n+K-N\\}\leq k\leq\min\\{n,K\\}$. If we consider the factorial of
negative integers as infinity, Eq. (1) is correct for all integers $k$.

The hypergeometric distribution is related to several distributions. When
picking balls _with replacement_ , the number of white balls follows the
binomial distribution [6], which is simpler than the hypergeometric
distribution. When we pick balls one by one without replacement until the
number of black balls reaches a given value, the number of white balls follows
the negative hypergeometric distribution [8]. The multivariate hypergeometric
distribution [2] appears when the balls have more than two colors. In the
Pólya urn model [7], at each step, one ball is drawn, its color is observed,
and the ball is returned in the urn with additional balls of the same color.
In this model, the number of white balls follows a distribution referred to as
the beta-binomial distribution [6].

Many publications, e.g., [6, 1, 2, 10], contain the moment-generating
function, the characteristic function, and the probability-generating function
of the hypergeometric distribution. The most fundamental function among them
is the probability-generating function

$G_{X}(z)\coloneqq
E\left[z^{X}\right]=\sum_{k=0}^{\infty}\mathrm{Pr}\left(X=k\right)z^{k}.$

Using $G_{X}$, the moment-generating function is expressed as
$M_{X}(t)=G_{X}(e^{t})$, the characteristic function as
$\varphi_{X}(t)=G_{X}(e^{it})$ ($i=\sqrt{-1}$), and the cumulant-generating
function as $K_{X}(t)=\ln G_{X}(e^{t})$. The existing formula for the
probability-generating function of the hypergeometric distribution is

$G_{X}(z)=\frac{(N-n)!(N-K)!}{N!(N-K-n)!}{}_{2}F_{1}\left(-n,-K;N-K-n+1;z\right).$
(2)

Here, ${}_{2}F_{1}$ is the Gauss hypergeometric function defined by

${}_{2}F_{1}\left(a,b;c;z\right)\coloneqq\sum_{k=0}^{\infty}\frac{(a)_{k}(b)_{k}}{(c)_{k}}\frac{z^{k}}{k!},$

where $(a)_{k}$ is the Pochhammer symbol:

$(a)_{0}=1,\quad(a)_{k}=a(a+1)\cdots(a+k-1).$

The Gauss hypergeometric function ${}_{2}F_{1}\left(a,b;c;z\right)$ diverges
when $c$ is zero or negative integer because $(c)_{k}=0$ for $k\geq|c|+1$.
Therefore, Eq. (2) is not properly defined when $n\geq N-K+1$. The factorial
$(N-K-n)!$ also diverges for such $n$, and Eq. (2) takes an indeterminate form
“$\infty/\infty$.” The condition $n\geq N-K+1$ indicates that the number $n$
of balls drawn is greater than the number of black balls. In this study, we
aim to establish the complete formula for $G_{X}(z)$ by providing the
determinate expression for $n\geq N-K+1$.

######  Theorem 1.1 (Main result).

The complete probability-generating function of $X$ becomes

$\displaystyle G_{X}(z)=$
$\displaystyle\frac{(N-n)!(N-K)!}{N!(N-K-n)!}{}_{2}F_{1}\left(-n,-K;N-K-n+1;z\right)$
$n\leq N-K$, (3a) $\displaystyle G_{X}(z)=$
$\displaystyle\frac{n!K!}{N!(n+K-N)!}z^{n+K-N}{}_{2}F_{1}\left(n-N,K-N;n+K-N+1;z\right)$
$n\geq N-K+1$. (3b)

Equation (3b) is the correct determinate formula for $n\geq N-K+1$. Although
this result might be trivial for expert mathematicians, the explicit
expression will be useful for students, engineers, practitioners, and
programmers who are not necessarily familiar with mathematics.

######  Remark 1.1.

As the first two arguments of the hypergeometric function, $-n$ and $-K$ in
Eq. (3a), and $n-N$ and $K-N$ in Eq. (3b), are negative integers, the infinite
sum of ${}_{2}F_{1}$ becomes finite. Therefore, ${}_{2}F_{1}$ in Eq. (3b)
becomes a polynomial of $z$ for any $N$, $K$, and $n$, whose degree is
$\min\\{n,K\\}$ for Eq. (3a) and $\min\\{N-n,N-K\\}$ for Eq. (3b). In both
Eqs. (3a) and (3b), $G_{X}(z)$ has degree $\min\\{n,K\\}$; the degree of Eq.
(3b) is calculated as $n+K-N+\min\\{N-n,N-K\\}=\min\\{n,K\\}$.

######  Remark 1.2.

We set the range of Eq. (3b) as $n\geq N-K+1$ in order to supplement the range
$n\leq N-K$ of Eq. (3a). However, Eq. (3b) is defined even for $n=N-K$, and
both Eqs. (3a) and (3b) for $n=N-K$ have the same form

$G_{X}(z)=\frac{n!K!}{N!}{}_{2}F_{1}\left(-n,-K;1;z\right).$

See Remark 3.1 and Fig. 1 for this overlap in detail.

The author believes that it is worth presenting the moment-generating function
and characteristic function of the hypergeometric random variable $X$
explicitly, because they appear more frequently than the probability-
generating function in many publications.

######  Corollary 1.1 (Moment-generating function and characteristic
function).

The moment-generating function $M_{X}(t)=G_{X}(e^{t})$ and characteristic
function $\varphi_{X}(t)=G_{X}(e^{it})$ of $X$ are

$M_{X}(t)=\begin{dcases}\frac{(N-n)!(N-K)!}{N!(N-K-n)!}{}_{2}F_{1}\left(-n,-K;N-K-n+1;e^{t}\right)&n\leq
N-K,\\\
\frac{n!K!}{N!(n+K-N)!}z^{n+K-N}{}_{2}F_{1}\left(n-N,K-N;n+K-N+1;e^{t}\right)&n\geq
N-K+1,\end{dcases}$

and

$\varphi_{X}(t)=\begin{dcases}\frac{(N-n)!(N-K)!}{N!(N-K-n)!}{}_{2}F_{1}\left(-n,-K;N-K-n+1;e^{it}\right)&n\leq
N-K,\\\
\frac{n!K!}{N!(n+K-N)!}z^{n+K-N}{}_{2}F_{1}\left(n-N,K-N;n+K-N+1;e^{it}\right)&n\geq
N-K+1,\end{dcases}$

respectively.

##  2 Proofs of Theorem 1.1

We do not prove Eq. (3a), which is a well-known existing result and present
only the proof of Eq. (3b) for the $n\geq N-K+1$ case. We provide two proofs.

###### First proof.

When $n\geq N-K+1$, the lower bound of the support of $X$ becomes
$\max\\{0,n+K-N\\}=n+K-N$. The probability-generating function can be computed
as follows.

$\displaystyle G_{X}(z)$
$\displaystyle=\frac{n!(N-n)!}{N!}\sum_{k=n+K-N}^{\infty}\frac{K!}{k!(K-k)!}\frac{(N-K)!}{(n-k)!(N-K-n+k)!}z^{k}$
$\displaystyle=\frac{n!K!}{N!(n+K-N)!}z^{n+K-N}\sum_{j=0}^{\infty}\frac{(n+K-N)!}{(n+K-N+j)!}\frac{(N-n)!}{(N-n-j)!}\frac{(N-K)!}{(N-K-j)!}\frac{z^{j}}{j!}$
$\displaystyle=\frac{n!K!}{N!(n+K-N)!}z^{n+K-N}\sum_{j=0}^{\infty}\frac{(n-N)_{j}(K-N)_{j}}{(n+K-N+1)_{j}}\frac{z^{j}}{j!}$
$\displaystyle=\frac{n!K!}{N!(n+K-N)!}z^{n+K-N}{}_{2}F_{1}\left(n-N,K-N;n+K-N+1;z\right),$

where $j=k-n-K+N$ is introduced in the second equality, and relations

$\displaystyle\frac{(M+j)!}{M!}$
$\displaystyle=(M+1)(M+2)\cdots(M+j)=(M+1)_{j},$
$\displaystyle\frac{M!}{(M-j)!}$
$\displaystyle=M(M-1)\cdots(M-j+1)=(-1)^{j}(-M)(-M+1)\cdots(-M+j-1)=(-1)^{j}(-M)_{j}$

are used in the third equality.

∎

###### Second proof.

We introduce the random variable $X^{\prime}$ as the number of _black_ balls
left in the urn. The number of white balls left becomes $N-n-X^{\prime}$, and
the relation $K=X+(N-n-X^{\prime})$ holds for the total number of white balls.
Therefore, we obtain

$X=X^{\prime}-N+K+n.$

Picking $n$ balls from the urn is equivalent to choosing $N-n\eqqcolon
n^{\prime}$ balls to be left in the urn. That is, $X^{\prime}$ follows a
hypergeometric distribution representing the number of black balls in randomly
selected $n^{\prime}$ balls, out of total $K^{\prime}=N-K$ black balls. When
$n\geq N-K+1$, $n^{\prime}$ satisfies $n^{\prime}\leq
N-K^{\prime}-1<N-K^{\prime}$. Hence, the probability-generating function of
$X^{\prime}$ is expressed in the form of Eq. (3a) with $(n,K)$ replaced by
$(n^{\prime},K^{\prime})$:

$\displaystyle G_{X^{\prime}}(z)$
$\displaystyle=\frac{(N-n^{\prime})!(N-K^{\prime})!}{N!(N-K^{\prime}-n^{\prime})!}{}_{2}F_{1}\left(-n^{\prime},-K^{\prime};N-K^{\prime}-n^{\prime}+1;z\right)$
$\displaystyle=\frac{n!K!}{N!(n+K-N)!}{}_{2}F_{1}\left(n-N,K-N;n+K-N+1;z\right).$

In conclusion,

$\displaystyle G_{X}(z)$
$\displaystyle=E\left[z^{X}\right]=E\left[z^{X^{\prime}+n+K-N}\right]=z^{n+K-N}G_{X^{\prime}}(z)$
$\displaystyle=\frac{n!K!}{N!(n+K-N)!}z^{n+K-N}{}_{2}F_{1}\left(n-N,K-N;n+K-N+1;z\right).$

∎

######  Remark 2.1.

The Gauss hypergeometric function ${}_{2}F_{1}\left(a,b;c;z\right)$ cannot be
defined for $c=0,-1,-2,\ldots$, but the limit

$\lim_{c\to-m}\frac{{}_{2}F_{1}\left(a,b;c;z\right)}{\Gamma(c)}=\frac{(a)_{m+1}(b)_{m+1}}{(m+1)!}z^{m+1}{}_{2}F_{1}\left(a+m+1,b+m+1;m+2;z\right)$
(4)

holds for nonnegative integer $m$ (see [9]), where $\Gamma$ is the gamma
function. In this sense, ${}_{2}F_{1}\left(a,b;c;z\right)/\Gamma(c)$ exists
for any $c$. Using this property, we can attain Eq. (3b) from Eq. (3a) as a
formality by taking a suitable limit. Equation (3a) can be rewritten as

$G_{X}(z)=\frac{(N-n)!(N-K)!}{N!}\lim_{\varepsilon\to
0}\frac{{}_{2}F_{1}\left(-n,-K;N-K-n+1+\varepsilon;z\right)}{\Gamma(N-K-n+1+\varepsilon)}.$

The right-hand side remains finite even for $n\geq N-K+1$, which becomes

$\frac{(N-n)!(N-K)!}{N!}\frac{(-n)_{n+K-N}(-K)_{n+K-N}}{(n+K-N)!}z^{n+K-N}{}_{2}F_{1}\left(K-N,n-N;n+K-N+1;z\right),$

owing to Eq. (4). By using

$(-n)_{n+K-N}=(-n)(-n+1)\cdots(K-N-1)=(-1)^{n+K-N}\frac{n!}{(N-K)!}$

and

$(-K)_{n+K-N}=(-K)(-K+1)\cdots(n-N-1)=(-1)^{n+K-N}\frac{K!}{(N-n)!},$

Eq. (3b) is obtained.

##  3 Alternative expression

As stated in the second proof in the previous section, picking $n$ balls is
equivalent to selecting $N-n$ balls to be left in the urn. It is expected that
a kind of symmetry for $n$ and $N-n$ exists in $G_{X}(z)$. However, Eq. (3b)
does not directly display this symmetry. In fact, Eq. (3b) switches its form
at $n=N-K$, not $n=N/2$. We seek an alternative expression of $G_{X}(z)$ in
this section, focusing on this symmetry.

######  Corollary 3.1.

The probability-generating function $G_{X}(z)$ in Eq. (3b) is rewritten in the
following form.

  1. $\mathrm{(a)}$

When $K\leq N/2$,

$\displaystyle G_{X}(z)=$
$\displaystyle\frac{(N-n)!(N-K)!}{N!(N-K-n)!}{}_{2}F_{1}\left(-n,-K;N-K-n+1;z\right)$
$n\leq\dfrac{N}{2}$, (5a) $\displaystyle G_{X}(z)=$
$\displaystyle\frac{n!(N-K)!}{N!(n-K)!}z^{K}{}_{2}F_{1}\left(n-N,-K;n-K+1;z^{-1}\right)$
$n\geq\dfrac{N}{2}$. (5b)

  2. $\mathrm{(b)}$

When $K\geq N/2$,

$\displaystyle G_{X}(z)=$
$\displaystyle\frac{(N-n)!K!}{N!(K-n)!}z^{n}{}_{2}F_{1}\left(-n,K-N;K-n+1;z^{-1}\right)$
$n\leq\dfrac{N}{2}$, (6a) $\displaystyle G_{X}(z)=$
$\displaystyle\frac{n!K!}{N!(n+K-N)!}z^{n+K-N}{}_{2}F_{1}\left(n-N,K-N;n+K-N+1;z\right)$
$n\geq\dfrac{N}{2}$. (6b)

As a formality, Eq. (5b) is formed from Eq. (5a) by replacing $n$ with $N-n$
and $z$ with $z^{-1}$, and multiplying $z^{K}$. The same replacement rule is
applicable to Eq. (6b). The factor $z^{K}$ represents that the symmetry is
incomplete.

###### Proof.

We use the formula [9]

${}_{2}F_{1}\left(-m,b;c;z\right)=\frac{(b)_{m}}{(c)_{m}}(-z)^{m}{}_{2}F_{1}\left(-m,1-c-m;1-b-m;z^{-1}\right),$
(7)

where $m$ is a nonnegative integer.

$\mathrm{(a)}$ Equation (5a) is the same as Eq. (3a) except for the range of
$n$. When $K\leq N/2$, the range $n\leq N/2$ is included in $n\leq N-K$, and
Eq. (5a) is proven.

To obtain Eq. (5b), we apply the formula (7) to Eq. (3a) with $m=K$, $b=-n$,
and $c=N-K-n+1$:

$\displaystyle\frac{(N-n)!(N-K)!}{N!(N-K-n)!}{}_{2}F_{1}\left(-n,-K;N-K-n+1;z\right)$
$\displaystyle=\frac{(N-n)!)N-K)!}{N!(N-K-n)!}\frac{(-n)_{K}}{(N-K-n+1)_{K}}(-z)^{K}{}_{2}F_{1}\left(n-N,-K;n-K+1;z^{-1}\right).$
(8)

By substituting relations

$\displaystyle(-n)_{K}$
$\displaystyle=(-n)(-n+1)\cdots(-n+K-1)=(-1)^{K}\frac{n!}{(n-K)!},$
$\displaystyle(N-K-n+1)_{K}$
$\displaystyle=(N-K-n+1)(N-K-n+2)\cdots(N-n)=\frac{(N-n)!}{(N-K-n)!},$

we obtain

$G_{X}(z)=\frac{n!(N-K)!}{N!(n-K)!}z^{K}{}_{2}F_{1}\left(n-N,-K;n-K+1;z^{-1}\right).$

The range of $n$ where hypergeometric functions in Eq. (8) are properly
defined is $K\leq n\leq N-K$. Therefore, Eq. (5b) for $N/2\leq n\leq N-K$ is
proven.

Equation (5b) for $n\geq N-K$ can be derived from Eq. (3b). By using Eq. (7)
with $m=N-n$, $b=K-N$, and $c=n+K-N+1$, Eq. (3b) becomes

$\displaystyle\frac{n!K!}{N!(n+K-N)!}\frac{(K-N)_{N-n}}{(n+K-N+1)_{N-n}}(-1)^{N-n}z^{K}{}_{2}F_{1}\left(n-N,-K;n-K+1;z^{-1}\right)$
$\displaystyle=\frac{n!(N-K)!}{N!(n-K)!}z^{K}{}_{2}F_{1}\left(n-N,-K;n-K+1;z^{-1}\right),$

where

$(K-N)_{N-n}=(-1)^{N-n}\frac{(N-K)!}{(n-K)!},\quad(n+K-N+1)_{N-n}=\frac{K!}{(n+K-N)!}.$

$\mathrm{(b)}$ Equation (6b) can be directly obtained from Eq. (3b), as Eq.
(5a) was derived in $\mathrm{(a)}$ from Eq. (3a). As above, Equation (6a) is
derived from Eq. (3b) by using the formula (7), where Eqs. (3a) and (3b)
covers $n\leq N-K$ and $N-K\leq n\leq N/2$, respectively. ∎

(a) (b)

Figure 1:  Covering regions of Eqs. (3b), (5b), and (6b) for $N=4$ on the
$n$–$K$ plane.

######  Remark 3.1.

In Corollary 3.1, “$n\leq N/2$” and “$n\geq N/2$” are used to highlight the
symmetry of $n$ and $N-n$. These ranges can be extended as long as the third
argument of the hypergeometric function is a positive integer. For example,
Eqs. (5a) and (5b) are valid for $n\leq N-K$ and $n\geq K$, respectively; both
equations include $K\leq n\leq N-K$, and they become identical. Similarly, the
ranges of $n$ in Eqs. (6a) and (6b) can be extended to $n\leq K$ and $n\geq
N-K$, respectively. Hence, Eqs. (5a) and (6b) are essentially identical to
Eqs. (3a) and (3b), respectively. Figure 1 shows regions where the equations
in Theorem 1.1 and Corollary 3.1 are defined. In Fig. 1(a), the overlap of two
regions consists of the points on the descending diagonal $n=N-K$. This
overlap is mentioned in Remark 1.2. In Fig. 1(b), the overlap of Eqs. (5b) and
(6a) occurs along the ascending diagonal $n=K$. Furthermore, when $N$ is even,
the $n=K=N/2$ case is included in all four regions, and $G_{X}(z)$ is the same
regardless of which equation is used. Specifically, when $N=2m$ and $N=K=m$,
the probability-generating function $G_{X}(z)$ can be written using the
Legendre polynomial $P_{m}$:

$G_{X}(z)=\frac{(m!)^{2}}{(2m)!}{}_{2}F_{1}\left(-m,-m;1;z\right)=\frac{(m!)^{2}}{(2m)!}(z-1)^{m}{}_{2}F_{1}\left(-m,m+1;1;(1-z)^{-1}\right)=\frac{(m!)^{2}}{(2m)!}(z-1)^{m}P_{m}\left(\frac{z+1}{z-1}\right),$

where the formula [9]

${}_{2}F_{1}\left(-m,b;c;z\right)=\frac{(b)_{m}}{(c)_{m}}(1-z)^{m}{}_{2}F_{1}\left(-m,c-b;1-b-m;(1-z)^{-1}\right)$

is used in the third equality, and the hypergeometric representation [9] of
the Legendre polynomial

$P_{m}(x)={}_{2}F_{1}\left(-m,m+1;1;\frac{1-x}{2}\right)$

is used in the last equality.

##  4 Moments

A typical application of the probability- and moment-generating functions is
the calculation of moments. In many books, e.g., [6, 1, 8, 2, 10], the mean
and variance of the hypergeometric distribution become

$E\left[X\right]=\frac{nK}{N},\quad
V\left[X\right]=\frac{n(N-n)K(N-K)}{N^{2}(N-1)},$

respectively. More generally, factorial moments,
$E\left[X(X-1)\cdots(X-r+1)\right]$ for $r=1,2,\ldots$, can be exactly
calculated for the hypergeometric distribution:

$E\left[X(X-1)\cdots(X-r+1)\right]=\frac{n!K!(N-r)!}{N!(K-r)!(N-r)!}.$ (9)

Although direct computation of Eq. (9) based on the probability mass function
(1) is possible, the factorial moments can be systematically calculated using
the probability-generating function $G_{X}$ as

$E\left[X(X-1)\cdots(X-r+1)\right]=\left.\frac{d^{r}G_{X}(z)}{dz^{r}}\right|_{z=1}.$

Using the probability-generating function in Eq. (3a), the factorial moment in
Eq. (9) can be easily obtained using formulas (see [9]):

$\frac{d^{r}}{dz^{r}}{}_{2}F_{1}\left(a,b;c;z\right)=\frac{(a)_{r}(b)_{r}}{(c)_{r}}{}_{2}F_{1}\left(a+r,b+r;c+r;z\right)$

and

${}_{2}F_{1}\left(a,b;c;1\right)=\frac{\Gamma(c)\Gamma(c-a-b)}{\Gamma(c-a)\Gamma(c-b)}.$

What computation is performed for $n\geq N-K+1$ in which Eq. (3b) is applied?
The effective formula [9] in this case is

$\frac{d^{r}}{dz^{r}}(z^{c-1}{}_{2}F_{1}\left(a,b;c;z\right))=(c-r)_{r}z^{c-r-1}{}_{2}F_{1}\left(a,b;c-r;z\right).$
(10)

For ${}_{2}F_{1}\left(n-N,K-N;n+K-N+1;z\right)$ ($a=n-N$, $b=K-N$, and
$c=n+K-N+1$) in Eq. (3b), $c-r=n+K-N+1-r$ becomes $0$ or a negative integer
when $r\geq n+K-N+1$. Therefore, ${}_{2}F_{1}\left(a,b;c-r;z\right)$ cannot be
defined for such $r$, and a special treatment is necessary.

######  Proposition 4.1 (Supplement formula to Eq. (10)).

For positive integers $m$ and $r$ which satisfy $m\leq r$,

$\frac{d^{r}}{dz^{r}}(z^{m-1}{}_{2}F_{1}\left(a,b;m;z\right))=\frac{(m-1)!(a)_{r-m+1}(b)_{r-m+1}}{(r-m+1)!}{}_{2}F_{1}\left(a+r-m+1,b+r-m+1;r-m+2;z\right).$

###### Proof.

Using the gamma function instead of the Pochhammer symbol, Eq. (10) can be
written as

$\frac{d^{r}}{dz^{r}}(z^{c-1}{}_{2}F_{1}\left(a,b;c;z\right))=\frac{\Gamma(c)}{\Gamma(c-r)}z^{c-r-1}{}_{2}F_{1}\left(a,b;c-r;z\right).$

In taking the $c\to m$ limit, Eq. (4) is used because $m-r$ is a nonpositive
integer:

$\lim_{c\to
m}\frac{{}_{2}F_{1}\left(a,b;c-r;z\right)}{\Gamma(c-r)}=\frac{(a)_{r-m+1}(b)_{r-m+1}}{(r-m+1)!}z^{r-m+1}{}_{2}F_{1}\left(a+r-m+1,b+r-m+1;r-m+2;z\right).$

∎

In the computation of the $r$th factorial moment of $X$ for $n\geq N-K+1$, Eq.
(10) can be directly applied for $r\leq n+K-N$, and Proposition 4.1 is needed
for $r\geq n+K-N+1$. In both cases, the final formula of the factorial moment
becomes identical with the existing formula (9). Therefore, the moments are
not affected when starting from the new formula (3b) for $n\geq N-K+1$. The
author presumes that Eq. (3b) has not focused thus far because the moments
become the same as the formulas obtained from the existing expression (3a).
Moreover, individual values of generating functions have probably been less
frequently required, which has kept researchers away from careful
investigation of the generating functions of the hypergeometric distribution.

## References

  * [1] Balakrishnan, N. and Nevzorov, V. B. (2004) A Primer on Statistical Distributions.  New York: Wiley.
  * [2] Bishop, Y. M. M., Fineberg, S. E., and Holland, P. W. (1975) Discrete Multivariate Analysis: Theory and Practice.  Cambridge, MA: MIT Press.
  * [3] Casiraghi, G., Nanumyan, V., Scholtes, I., and Schweitzer, F. (2017) From relational data to graphs: inferring significant links using Generalized hypergeometric ensembles, volume 10540 of Lecture Notes in Computer Science, (pp. 111–120).  Springer: Cham.
  * [4] Fury, W., Batliwalla, F., Gregersen, P. K., and Li, W. (2006) Overlapping probabilities of top ranking gene lists, hypergeometric distribution, and stringency of gene selection criterion.  In 2006 International Conference of the IEEE Engineering in Medicine and Biology Society (pp. 5531–5534).
  * [5] Hafid, A., Hafid, A. S., and Samih, M. (2020) A novel methodology-based joint hypergeometric distribution to analyze the security of sharded blockchains.  In IEEE Access, volume 8 (pp. 179389–179899).
  * [6] Johnson, N. L., Kotz, S., and Kemp, A. W. (1992) Univariate Discrete Distributions.  New York: Wiley, 2nd edition.
  * [7] Mahmoud, H. M. (2009) Pólya Urn Models.  Boca Raton: CRC Press.
  * [8] Mukhopadhyay, P. (2012) An Introduction to the Theory of Probability.  Singapore: World Scientific.
  * [9] Olver, F. W., Lozier, D. W., Boisvert, R. F., and Clark, C. W. (2010) NIST Handbook of Mathematical Functions.  Cambridge: Cambridge University Press.
  * [10] Patel, J. K., Kapadia, C. H., and Owen, D. B. (1976) Handbook of Statistical Distributions.  New York: Dekker.
  * [11] Rivals, I., Personnaz, L., Taing, L., and Potier, M.-C. (2007) Enrichment or depletion of a GO category within a class of genes: which test?. Bioinformatics 23 401–407.
  * [12] Seber, G. A. F. (2013) Capture–Recapture methods—I, (pp. 35–42).  Wiley: New York.
  * [13] Wittes, J. and Wallenstein, S. (1987) The power of the Mantel–Haenszel test. Journal of the American Statistical Association 82 1104–1109.