# Spikes and accretion of unbound, collisionless matter around black holes

Stuart L. Shapiro  Departments of Physics and Astronomy, University of
Illinois at Urbana-Champaign, Urbana, IL 61801, USA

###### Abstract

We consider the steady-state density and velocity dispersion profiles of
collisionless matter around a Schwarzschild black hole (BH) and its associated
rate of accretion onto the BH. We treat matter, which could be stars or dark
matter particles, whose orbits are unbound to the BH, but still governed by
its gravitational field. We consider two opposite spatial geometries for the
matter distributions: an infinite, 3D cluster and a 2D razor-thin disk, both
with zero net angular momentum. We demonstrate that the results depend
critically on the adopted geometry, even in the absence of angular momentum.
We adopt a simple monoenergetic, isotropic, phase-space distribution function
for the matter as a convenient example to illustrate this dependence. The
effect of breaking strict isotropy by incorporating an unreplenished loss cone
due to BH capture of low-angular momentum matter is also considered.
Calculations are all analytic and performed in full general relativity, though
key results are also evaluated in the Newtonian limit. We consider one
application to show that the rate of BH accretion from an ambient cluster can
be significantly less than that from a thin disk to which it may collapse,
although both rates are considerably smaller than Bondi accretion for a
(collisional) fluid with a similar asymptotic particle density and velocity
dispersion (i.e., sound speed).

##  I Introduction

A black hole (BH), and especially a supermassive black hole (SMBH), typically
will steepen the density profile of stars and/or dark matter (DM) within the
hole’s sphere of influence, i.e., within a radius $r_{h}\approx
GM/v_{\infty}^{2}$. Here, $M$ is the mass of the hole and $v_{\infty}$ is the
velocity dispersion in the cluster or galaxy core outside $r_{h}$. The density
profile of this spike depends on the properties of the matter, as well as the
formation and evolution history of the BH. Most of the treatments to date
focus on matter bound to the BH and hence moving on orbits that are confined
inside $r_{h}$. In some analyses the matter is treated as completely
collisionless, as was assumed for stars in typical galaxies in Refs Peebles
(1972); Young (1977), and for cold DM in galaxies and halos in Ref Gondolo and
Silk (1999), where in both of these cases the orbiting species were assumed
bound to adiabatically growing BHs. In other cases the matter is assumed to be
weakly collisional and slowly evolving on long relaxation timescales. Such
relaxation may be governed by distant, cumulative, two-body, gravitational
encounters (Coulomb scattering) for stars bound to a central, massive BH in a
globular clusters or to a SMBH in a dense galaxy core Bahcall and Wolf (1976);
Frank and Rees (1976); Lightman and Shapiro (1977); Shapiro and Marchant
(1978); Cohn and Kulsrud (1978). For collisionless DM bound to a SMBH in a
galaxy core, relaxation may occur due to gravitational encounters with ambient
stars Gnedin and Primack (2004); Merritt (2004); Shapiro and Heggie (2022),
or, in the case of self-interacting DM (i.e. SIDM) to particle self-
interactions Shapiro and Paschalidis (2014); Shapiro (2018); Luca and Khoury
(2023), and/or annihilations Vasiliev (2007); Shapiro and Shelton (2016). In
all of these instances the orbiting matter forms power-law density spikes
inside $r_{h}$.

Here we consider collisionless matter, either dark matter or stars (both of
which we shall refer to as “particles”) that move on orbits unbound to the BH.
These unbound orbits can take particles infinitely far away from the BH, but
particles that plunge inside $r_{h}$ on such orbits also may, in some
instances, generate a density spike around the BH. Those particles traveling
inward with sufficiently low angular momentum about the BH are eventually
captured. By collisionless here we mean that the particles interact solely
with the gravitational field of the central BH and that any scattering due to,
or collisions with, their neighbors, themselves or other intruders remain
generally unimportant over the age of the system. We take the BH to be a
nonspinning Schwarzschild BH with a mass $M$ fixed in time and we determine
the density profile and accretion rate of the unbound ambient matter.
Collsionless particles move on geodesic orbits about the BH and for this
analysis we assume their self-gravity is unimportant. We take all of the
particles to have the same mass and to describe them we adopt a simple,
monoenergetic, phase-space distribution function of the form $f=f(E)$, where
$E$ is the conserved “energy at infinity” per unit mass of a particle. Such a
distribution function has been adopted in several previous investigations of
unbound particles about a BH Zeldovich and Novikov (1971); Shapiro and
Teukolsky (1983, 1985); Gamboa et al. (2021) and, though idealized, it is
sufficient to illustrate our main conclusions. Any distribution function that
depends solely on conserved integrals of the motion, such as $E$,
automatically satisfies the time-independent, collisionless Boltzmann (Vlasov)
equation. Hence the profiles and accretion rates we determine from our adopted
distribution function yield steady-state solutions.

Given our adopted distribution function we compare two opposite spatial
geometries for the unbound matter – infinite 3-D clusters and 2-D razor-thin
disks. The adopted disks, though extreme, serve to highlight the differences
between extended, spherical-like vs. thin disk-like density distributions for
unbound, collisionless particles orbiting a central BH and their associated
rates of accretion. We note that collisionless dark matter in the early
universe, and even the first generation stars, may in fact form thin sheets or
“pancakes” Zel’dovich (1970); Kulkarni et al. (2022), so our extreme disks may
mimic some of their features, should massive BHs reside within them. In both
our cluster and disk cases we analyze the case where the orbital velocities
are isotropic at every point. We then consider the anisotropic velocity
distribution associated with an unreplenished loss cone that can arise from BH
capture of particles with sufficiently low angular momentum. In all cases the
net angular momentum of the orbiting matter is assumed to be zero, so we are
only examining the effect of spatial geometry on the profiles and accretion
rates.

Infinite, unbound collisionless clusters around Schwarzschild BHs have been
previously investigated with a different approach in Refs Rioseco and Sarbach
(2017, 2017), who adopted a Maxwell-Jünter phase-space distribution function
as an application. A treatment involving thin (equatorial plane) disks around
Kerr black holes was recently provided by Ref Cieślik et al. (2022), again
with a Maxwell-Jünter distribution function and a different approach from the
one adopted here. We compare the accretion rates found by these studies with
the ones calculated here.

The significance of determining the density profiles and accretion rates of
stars and DM near BHs is that it provides clues as to the nature, formation
history and evolution of the systems in which they are found. The presence of
a sufficiently steep density spike not only conveys the presence of a massive,
central BH but, in the case of DM, may lead to an observable excess of gamma
rays or other form of radiation, if particle annihilation occurs at a
sufficient rate in the innermost regions where the density is highest [see
Daylan et al. (2016); Cholis et al. (2022); Leane et al. (2022); Gondolo and
Silk (1999); Fields et al. (2014); Wanders et al. (2015); Chiang et al. (2020)
and references therein]. These applications motivate our adopting a simple
illustration to point out the importance of unbound collisionless particles
and their global spatial geometry in determining the existence of a stellar,
or DM, spike around a BH and their consumption rate by the BH.

In Section II we consider an infinite, unbound, monoenergetic, 3D cluster of
collisionless particles that orbit a Schwarzschild BH. We derive the particle
density and velocity profiles and the associated particle accretion rate onto
the BH. We consider both an isotropic velocity distribution and one with an
unreplenished loss cone. The derivations and quoted quantities are analytic
and fully general relativistic, but the final expressions are also quoted for
the slow-velocity, weak-field (i.e. Newtonian) limit. In Section III we repeat
the analysis in Section I, but now for particles confined to a 2D razor-thin
disk. In Section IV we compare the accretion rates found for the clusters and
disks. In the Appendix we provide two alternative derivations for the cluster
accretion rate, all arriving at the same answer. We adopt geometric units
throughout, setting $G=1=c$, unless otherwise noted.

##  II Unbound, Monoenergetic Cluster

###  II.1 Isotropic Distribution

The gravitational field of our system is governed by the Schwarzschild BH and
we adopt the following familiar form for the spacetime metric:

$ds^{2}=-(1-2M/r)dt^{2}+\frac{dr^{2}}{1-2M/r}+r^{2}d\theta^{2}+r^{2}{\rm
sin}^{2}\theta d\phi^{2}.$ (1)

The monoenergetic phase-space distribution function we shall employ everywhere
in space is given by

$f(E)=K\delta(E-E_{\infty}),\ \ \ K,E_{\infty}\ {\rm constant},$ (2)

where $E=-{\bf p\cdot e_{t}}$ is the “energy at infinity” of a particle per
unit mass, ${\bf p}$ is the particle energy-momentum 4-vector, divided by its
rest mass, and ${\bf e_{t}}=\partial/\partial t$ is the time coordinate basis
vector. All particles are assumed to have the same rest mass $m$. The constant
$E_{\infty}$ satisfies $E_{\infty}>1$, appropriate for an unbound orbit, while
the constant $K$, given $E_{\infty}$, is determined by the asymptotic density
at infinity [see Eq. 9 below]. This distribution function, specifically chosen
to describe a homogeneous density of unbound, collisionless particles moving
randomly at the same speed far from the BH, in fact applies everywhere as a
consequence of Liouville’s theorem. Our adopted function, though simple, may
also be thought of as a Green’s function for other energy distribution
functions (e.g., Maxwellians or power laws) and our results below in which
$E_{\infty}$ appears can always be integrated over $E_{\infty}$ for these
alternative energy distributions.

####  II.1.1 Density

We now repeat the derivation of the number density of particles at all radii
$r$ we provided previously (see Appendix A in Ref Shapiro and Teukolsky
(1985)), so that we may refer to and modify some of the equations in
subsequent sections where we change the boundary conditions or contrast the
results with other cases. The particle energy per unit mass measured by a
local, static, orthonormal observer with 4-velocity ${\bf u}={\bf
e_{\hat{t}}}$ at radius $r$ is given by

$E_{local}\equiv p^{\hat{t}}=-{\bf p\cdot
e_{\hat{t}}}=\frac{E}{(1-2M/r)^{1/2}},$ (3)

where a caret on a variable denotes an orthonormal component. The number
density of particles at $r$ measured by this observer is then

$\displaystyle n(r)$ $\displaystyle=$ $\displaystyle\int f(E)d^{3}\hat{p}$ (4)
$\displaystyle=$ $\displaystyle 4\pi\int f(E)\hat{p}^{2}d\hat{p},$ (5)

where

$\hat{p}=[(p^{\hat{t}})^{2}-1]^{1/2}$ (6)

is the particle 3-momentum per unit mass. Using Eqns. (3) and (6) we obtain

$\hat{p}^{2}d\hat{p}=\frac{\hat{p}EdE}{1-2M/r}=\frac{(E^{2}-1+2M/r)^{1/2}}{(1-2M/r)^{3/2}}EdE.$
(7)

Substituting Eqs. (2) and (7) into (5) gives

$n(r)=4\pi K\frac{(E_{\infty}^{2}-1+2M/r)^{1/2}E_{\infty}}{(1-2M/r)^{3/2}}.$
(8)

Evaluating Eq. (8) at $r=\infty$ determines $K$:

$K=\frac{n_{\infty}}{4\pi(E_{\infty}^{2}-1)^{1/2}E_{\infty}}.$ (9)

So in general the density profile is

$\frac{n(r)}{n_{\infty}}=\frac{(E_{\infty}^{2}-1+2M/r)^{1/2}}{(E_{\infty}^{2}-1)^{1/2}(1-2M/r)^{3/2}}.$
(10)

An important special case applies to nonrelativistic particles at infinity
moving with velocity $\hat{v}_{\infty}\ll 1$, for which

$E_{\infty}\approx 1+\frac{1}{2}\hat{v}_{\infty}^{2},\ \ \
E_{\infty}^{2}-1\approx 2(E_{\infty}-1)\approx\hat{v}_{\infty}^{2}.$ (11)

In this limit Eq. (9) reduces to

$K\approx\frac{n_{\infty}}{4\pi\hat{v}_{\infty}},\ \ \ \ (\hat{v}_{\infty}\ll
1),$ (12)

and Eq. (8) becomes

$\frac{n(r)}{n_{\infty}}\approx\frac{[1+2M/(r\hat{v}_{\infty}^{2})]^{1/2}}{(1-2M/r)^{3/2}},\
\ \ \ (\hat{v}_{\infty}\ll 1).$ (13)

We note that Eq. (13) in the weak-field domain, wherein $M/r\ll 1$, is
consistent with earlier Newtonian derivations in Refs. Zeldovich and Novikov
(1971); Shapiro and Teukolsky (1983).

####  II.1.2 Velocity Dispersion

The velocity dispersion may be calculated formally from

$\langle\hat{v}^{2}(r)\rangle=\frac{\int\hat{v}^{2}f(E)d^{3}\hat{p}}{\int
f(E)d^{3}\hat{p}},$ (14)

where the locally measured velocity is $\hat{v}=\hat{p}/E_{local}$ or, using
Eqs. (3) and (6),

$\hat{v}=\frac{(E^{2}-1+2M/r)^{1/2}}{E}.$ (15)

Substituting Eqs. (2) and (15) into (14) and integrating yields

$\langle\hat{v}^{2}(r)\rangle=\hat{v}^{2}(E_{\infty})=\frac{(E_{\infty}^{2}-1+2M/r)}{E_{\infty}^{2}}.$
(16)

The first equality in Eq. (16) also may be arrived at trivially by noting that
all particles move on geodesics with the same specific energy $E_{\infty}$,
and the second equality also can be obtained by invoking
$E_{local}=\gamma=1/(1-\hat{v}^{2})^{1/2}$, replacing $E_{local}$ by $E$ using
Eq. (3), and inverting for $\hat{v}$.

####  II.1.3 Accretion Rate

Here we generalize our earlier Newtonian derivation of the DM accretion rate
onto the BH [see Ref. Shapiro and Teukolsky (1983), Eqs. (14.2.13)-(14.2.20)]
to one that is fully general relativistic. The phase-space momentum element
$d^{3}\hat{p}$ may be expressed as

$d^{3}\hat{p}=2\pi p^{\hat{\bot}}dp^{\hat{\bot}}dp^{\hat{r}}=\frac{4\pi
JdJdE}{(1-2M/r)^{1/2}|v^{\hat{r}}|r^{2}},$ (17)

where $\bot$ denotes directions perpendicular to the radial direction (whereby
$J=rp^{\hat{\bot}}$ is the particle angular momentum per unit mass) and where
Eqs. (12.4.9), (12.4.13) and (12.4.16) of Ref. Shapiro and Teukolsky (1983)
were employed to relate $dp^{\hat{r}}$ to $dE$:

$dp^{\hat{r}}=\frac{dE}{(1-2M/r)^{1/2}|v^{\hat{r}}|}.$ (18)

An additional factor of 2 arises in Eq. (17) since for a given $E$,
$p^{\hat{r}}$ can be either positive or negative. Let $N^{-}(r,E,J)$ be the
number of particles per interval $dr,dE$, and $dJ$ with inward-directed radial
velocity:

$\displaystyle N^{-}(r,E,J)drdEdJ$ $\displaystyle=$
$\displaystyle\frac{1}{2}f(E,J)d^{3}rd^{3}\hat{p}$ $\displaystyle=$
$\displaystyle\frac{8\pi^{2}Jf(E,J)}{|v^{\hat{r}|}(1-2M/r)^{1/2}}drdEdJ.$

The resulting total capture rate for particles onto the central BH as measured
by the local observer is then

$\displaystyle\frac{dN_{tot}}{d\tau}$ $\displaystyle=$
$\displaystyle\int^{\infty}_{(1-\frac{2M}{r})^{\frac{1}{2}}}dE\int_{0}^{J_{min}(E)}dJ|v^{\hat{r}}|N^{-}(r,E,J)\
\ \ \ \ \ \ \ $ (20) $\displaystyle=$ $\displaystyle
8\pi^{2}\int^{\infty}_{(1-\frac{2M}{r})^{\frac{1}{2}}}dE\int_{0}^{J_{min}(E)}dJ\frac{Jf}{(1-\frac{2M}{r})^{1/2}},$

where particles moving with specific angular momentum less than a critical
value $J_{min}$ will be captured by BH as they approach pericenter. The region
$J<J_{min}$ thus defines a capture loss cone. For nonrelativistic (NR)
particles at large distances with $E-1\ll 1$, which is the case for typical
stars in clusters and galaxies, and cold dark matter, we note that

$J_{min}(E)=4M\ \ \ \ ({\rm NR\ particles},\ \hat{v}_{\infty}\ll 1).$ (21)

Substituting Eqs. (2) and (9) into Eq. (20) and evaluating the integral yields

$\frac{dN_{tot}}{d\tau}=\frac{\pi
n_{\infty}J_{min}^{2}(E_{\infty})}{E_{\infty}(E_{\infty}^{2}-1)^{1/2}}\frac{1}{(1-2M/r)^{1/2}}.$
(22)

The particle capture rate as measured by a static observer at infinity is then
given by

$\frac{dN_{tot}}{dt}=\frac{dN_{tot}}{d\tau}\frac{d\tau}{dt}=\frac{dN_{tot}}{d\tau}(1-\frac{2M}{r})^{1/2},$
(23)

or

$\frac{dN_{tot}}{dt}=\frac{\pi
n_{\infty}J_{min}^{2}(E_{\infty})}{E_{\infty}(E_{\infty}^{2}-1)^{1/2}}.$ (24)

It is this observer who measures the steady depletion of particles from the
ambient cluster and their acquisition by the BH.

The corresponding rest-mass accretion rate measured by this observer,
$dM_{0}/dt=mdN_{tot}/dt$, is then given by

$\frac{dM_{0}}{dt}=\frac{\pi\rho_{\infty}J_{min}^{2}(E_{\infty})}{E_{\infty}(E_{\infty}^{2}-1)^{1/2}},$
(25)

where $\rho_{\infty}\equiv mn_{\infty}$ is the asymptotic rest-mass density.
Equation (25) agrees with the result obtained by a different route in Ref
Gamboa et al. (2021), [see their Eq. (76), setting
$\gamma_{\infty}=E_{\infty}$ and $L_{c}(E_{\infty})=mJ_{min}(E_{\infty})$].
The rate of accretion of total mass-energy onto the black hole is then given
by $dM/dt=mE_{\infty}dN/dt$ or

$\frac{dM}{dt}=\frac{\pi\rho_{\infty}J_{min}^{2}(E_{\infty})}{(E_{\infty}^{2}-1)^{1/2}}.$
(26)

Evaluating Eqs. (25) and (26) for NR particles at infinity, using Eqs. (11)
and (21), yields

$\frac{dM}{dt}\approx\frac{dM_{0}}{dt}\approx\frac{16\pi
M^{2}\rho_{\infty}}{\hat{v}_{\infty}}\ \ \ \ ({\rm NR\ particles},\
\hat{v}_{\infty}\ll 1),$ (27)

which agrees with the previous Newtonian result we derived in Ref. Shapiro and
Teukolsky (1983) [see Eq. (14.2.20)]. Since the steady-state accretion rate
can be evaluated at large $r\gg M$, as we did in Ref. Shapiro and Teukolsky
(1983), it is no surprise that the Newtonian derivation for the rate provided
there for nonrelativistic particles at infinity yields the exact same result
derived in GR for the limiting case given by Eq. (27). Note that to convert
back from geometric to physical units, one simply multiplies the right-hand
sides of Eqs. (25)–(27) by $G^{2}/c^{2}$.

It is interesting to note that the particle density as measured by a locally
static observer blows up when the observer is stationed arbitrarily close to
the BH horizon [see Eq. (10) as $r\rightarrow 2M$]. (However, recall that
static observers cannot exist at or inside the horizon.) Using Eq. (12.4.17)
of Ref. Shapiro and Teukolsky (1983), for the particle radial velocity, i.e.,

$|v^{\hat{r}}|=\left[1-\frac{1}{E^{2}}(1-\frac{2M}{r})(1+\frac{J^{2}}{r^{2}})\right]^{1/2},$
(28)

it is seen that the particles all have velocities which are measured by this
observer to approach the speed of light and move in the radial direction as
they approach the horizon. Based on these observations one might naively take
the radial matter flux to be $\sim n(r)\hat{v}$, and then estimate the
accretion rate by multiplying this product by the invariant area $4\pi r^{2}$.
Doing so yields a rate that blows up near the horizon. However, this estimate
is too naive, as the actual steady-state depletion rate of unbound particles
from the ambient gas must be obtained by more careful considerations, such as
those we incorporated above, leading to the net accretion rates onto the BH as
measured by a distant, static observer, Eqs. (24)-(26). The latter remain
perfectly finite. This contrast is sufficiently striking that we are motivated
to provide two alternative derivations in the Appendix for the net accretion
rate, leading to the same equations obtained above.

The integration of an isotropic distribution function $f(E)$ over all phase
space allows for so-called “white hole” orbits, i.e., outgoing trajectories at
the horizon. Eliminating these orbits leads to a density profile of the form
$n(r)\sim 1/(1-2M/r)^{1/2}$ as $r\rightarrow 2M$. Thus, we still have $n(r)$
blowing up near the horizon, but only with a mild (“redshift”) factor instead
of the $1/(1-2M/r)^{3/2}$ behavior exhibited by Eq. (10). Using this
“modified” density near the horizon to measure the inward flux then yields the
correct locally measured accretion rate, Eq. (22). Furthermore, choosing to
measure the density by an observer comoving with the net flow should lead to a
finite density at the horizon Rioseco and Sarbach (2017, 2017). However, as
these modifications do not affect our computations of the accretion rates,
which are our primary targets, we do not implement them here.

####  II.1.4 Massless and Extremely Relativistic Particles

To evaluate Eq. (26) for massless particles (e.g. photons), we may first set
$E_{\infty}=(1-\hat{v}_{\infty}^{2})^{-1/2}$ and then note that
$J_{min}=b(E_{\infty}^{2}-1)^{1/2}=\hat{v}_{\infty}(1-\hat{v}_{\infty}^{2})^{-1/2}$,
where $b$ is the critical impact parameter for particle capture. Now
$b=3\sqrt{3}$ is the critical impact parameter for massless particles [see,
e.g., Ref. Shapiro and Teukolsky (1983), Eqs. (12.4.36) and (12.5.11)]. Next
let $n_{\infty}mE_{\infty}\equiv\epsilon_{\infty}^{r}$ be the energy density
of particles far from the BH. Then taking the limit as
$\hat{v}_{\infty}\rightarrow 1$ with these substitutions gives

$\frac{dM}{dt}=27\pi M^{2}\epsilon_{\infty}^{r}\ \ \ (m=0),$ (29)

which agrees with the result quoted in Exercise 14.4 in Ref. Shapiro and
Teukolsky (1983) that was obtained by an alternative approach.

We also note that Eq. (29) also applies to extremely relativistic (ER)
particles with nonzero rest mass $m$. For such particles the accretion rates
may also be written as

$\frac{dM}{dt}=E_{\infty}\frac{dM_{0}}{dt}\approx 27\pi
M^{2}E_{\infty}\rho_{\infty}\ \ \ ({\rm ER\ particles},\hat{v}_{\infty}\gg
1).$ (30)

###  II.2 Loss-Cone Effect

Given that particles having nearly radial, inward velocities result in their
being captured by the BH, it is interesting to examine a scenario in which,
after a sufficient time has passed and steady state is achieved, there is a
depletion of particles in a low-angular momentum capture loss cone about the
BH. For perfectly collisionless gas in the complete absence of perturbations,
the loss cone cannot be replenished, since particle self-interactions and
gravitational scattering by stars or other perturbers are assumed absent. (For
scenarios in which self-interactions of dark matter particles around BHs may
be important, see, e.g., Refs. (Shapiro and Paschalidis (2014); Shapiro
(2018); Luca and Khoury (2023)) and where their gravitational scattering off
stars may be significant, see, e.g., Refs. (Merritt (2004); Gnedin and Primack
(2004); Shapiro and Heggie (2022)), and references therein). In reality, the
slightest gravitational perturbations, at large distances, due, for example,
to small density anisotropies or to intruders, or the weakest self-
interactions, will likely be sufficient to replenish the narrow loss cone. Our
analysis above of the matter profiles and accretion rate will thus apply in
these situations. But immediately below we will consider a perturbation-free,
perfectly collisionless cluster that has depleted its loss cone and cannot
refill it.

To describe this idealized situation a minimal modification to our adopted
distribution function will suffice:

$\displaystyle f(E)$ $\displaystyle=$ $\displaystyle K\delta(E-E_{\infty}),\ \
J_{min}(E)\leq J\leq J_{max}(E),$ (31) $\displaystyle=$ $\displaystyle 0,\ \
0\leq J\leq J_{min}(E),$

where again $K$ and $E_{\infty}>1$ are constants.

####  II.2.1 Density

To determine the density profile we again evaluate Eq. (4), employing Eq.
(17), which yields

$n(r)=\int^{\infty}_{(1-\frac{2M}{r})^{\frac{1}{2}}}f(E)dE\int_{J_{min}(E)}^{J_{max}(E)}\frac{4\pi
JdJ}{r^{2}|v^{\hat{r}}|(1-\frac{2M}{r})^{1/2}}.$ (32)

where $|v^{\hat{r}}|$ is given by Eq. (28) and $J_{max}(E)$ is the maximum
specific angular momentum that a particle moving on a geodesic with specific
energy $E$ can have,

$J_{max}(E)=r\left[\frac{(E^{2}-1+\frac{2M}{r})}{(1-\frac{2M}{r})}\right]^{1/2}$
(33)

[set $dr/d\tau=0$ and $\tilde{l}=J_{max}(E)$ in Eq. (12.4.13) in Ref. Shapiro
and Teukolsky (1983)]. Substituting Eq. (31) into (32) and integrating yields

$\frac{n(r)}{n_{\infty}}=\frac{[E_{\infty}^{2}-(1-2M/r)(1+\left(J_{min}(E)/r\right)^{2})]^{1/2}}{(E_{\infty}^{2}-1)^{1/2}(1-2M/r)^{3/2}},$
(34)

where the normalization constant $K$ is again given by Eq. (9). Comparing Eqs
(10) and (34) shows that the spike density is somewhat lower everywhere, but
mostly around $r\lesssim J_{min}$, if a loss cone is established and
unreplenished.

Equation (34) reduces to Eq. (10) in the absence of a loss cone, i.e. when the
captured particle distribution is assumed to be continually replenished so
that $J_{min}(E)=0$ in Eq. (31). For nonrelativistic particles with
$\hat{v}_{\infty}\ll 1$ and an empty loss cone with $J_{min}=4M$ the density
becomes

$\displaystyle\frac{n(r)}{n_{\infty}}$ $\displaystyle\approx$
$\displaystyle\frac{[1+2M/(r\hat{v}_{\infty}^{2})(1-\frac{8M}{r}+\left(\frac{4M}{r}\right)^{2})]^{1/2}}{(1-2M/r)^{3/2}},$
(35) $\displaystyle\ \ \ \ ({\rm NR\ particles},\ \hat{v}_{\infty}\ll 1).$

For sufficiently small $\hat{v}_{\infty}\ll 1$ the density exhibits a minimum
at $r\approx 4M$ outside the event horizon.

####  II.2.2 Velocity Dispersion

As all particles have the same energy $E_{\infty}$, they have the same
velocity and velocity dispersion profile as in the absence of a loss cone,
i.e., Eq. (16).

####  II.2.3 Accretion Rate

Given a loss cone that is unreplenished with particles once they have been
captured, the accretion rate becomes zero in steady state.

##  III Unbound, Monoenergetic, Razor-Thin Disk

###  III.1 Isotropic Distribution

We now derive the surface density profile around the BH for unbound,
collisionless particles assumed to reside in a razor-thin disk. Since the
distribution function is zero outside of the disk, we will simply work in the
plane of the disk and drop the vertical dependence of any quantity. For
comparison purposes we adopt the same distribution function adopted for the
large 3D cluster, but now confined to the 2D disk plane, i.e., we again use
Eq. 2 to describe unbound particles moving isotropically in the plane of the
disk.

####  III.1.1 Density

The surface number density $\Sigma{\text{*}}(r)$ measured by a locally static
observer is then given by

$\Sigma{\text{*}}(r)=\int f(E)d^{2}\hat{p}=2\pi\int f(E)\frac{EdE}{1-2M/r},$
(36)

where we used Eq. (7) to evaluate $d^{2}\hat{p}=2\pi\hat{p}d\hat{p}$.
Substituting Eq. (2) in (36) and integrating yields

$\frac{\Sigma{\text{*}}(r)}{\ \ \Sigma{\text{*}}_{\infty}}=\frac{\Sigma(r)}{\
\ \Sigma_{\infty}}=\frac{1}{1-2M/r},$ (37)

where $\Sigma(r)=m\Sigma{\text{*}}(r)$ is the surface rest-mass density, $m$
is the particle rest mass and $K=\Sigma_{\infty}/(2\pi mE_{\infty})$. Far from
the BH the surface density is thus seen to be flat. Once again, eliminating
the “white hole” orbits softens the blowup near the horizon, whereby the
modified surface density scales as $\Sigma(r)\sim 1/(1-2M/r)^{1/2}$ as
$r\rightarrow 2M$.

####  III.1.2 Velocity Dispersion

As all particles again have the same energy $E_{\infty}$, they have the same
velocity and velocity dispersion profile given by Eq. (16).

####  III.1.3 Accretion Rate

To derive the accretion rate we adapt Eq. (II.1.3) to a plane, yielding

$\displaystyle N^{-}(r,E,J)drdEdJ$ $\displaystyle=$
$\displaystyle\frac{1}{2}f(E,J)d^{2}rd^{2}\hat{p}$ (38) $\displaystyle=$
$\displaystyle 4\pi\frac{dJdEdr}{(1-2M/r)^{1/2}}\frac{f(E,J)}{|v^{\hat{r}|}},\
\ \ \ $

where we used $d^{2}r=2\pi rdr$, $J=p^{\hat{\bot}}r$, Eq. (18) and

$d^{2}\hat{p}=dp^{\hat{\bot}}dp^{\hat{r}}=4\frac{dJ}{r}\frac{dE}{(1-2M/r)^{1/2}|v^{\hat{r}}|}.$
(39)

The factor of $4$ is inserted since for a given $E$ and $J$, both
$p^{\hat{r}}$ and $p^{\hat{\bot}}$ can be either positive or negative.

The resulting capture rate measured by the local observer is then given by

$\displaystyle\frac{dN_{tot}}{d\tau}$ $\displaystyle=$
$\displaystyle\int^{\infty}_{(1-\frac{2M}{r})^{\frac{1}{2}}}dE\int_{0}^{J_{min}(E)}dJ|v^{\hat{r}}|N^{-}(r,E,J)\
\ \ \ \ \ \ \ $ (40) $\displaystyle=$ $\displaystyle
4\pi\int^{\infty}_{(1-\frac{2M}{r})^{\frac{1}{2}}}dE\int_{0}^{J_{min}(E)}dJ\frac{f}{(1-\frac{2M}{r})^{1/2}},$

yielding

$\frac{dN_{tot}}{d\tau}=\frac{2\Sigma{\text{*}}_{\infty}J_{min}(E_{\infty})}{E_{\infty}}\frac{1}{(1-2M/r)^{1/2}}.$
(41)

The depletion rate of the disk as measured by a distant static observer is
then obtained using Eq. (23), which gives

$\frac{dN_{tot}}{dt}=\frac{2\Sigma{\text{*}}_{\infty}J_{min}(E_{\infty})}{E_{\infty}}.$
(42)

The corresponding rates of accretion of rest mass and total mass-energy onto
the BH are then

$\frac{dM_{0}}{dt}=\frac{2\Sigma_{\infty}J_{min}(E_{\infty})}{E_{\infty}}$
(43)

and

$\frac{dM}{dt}=2\Sigma_{\infty}J_{min}(E_{\infty}),$ (44)

respectively. The rate for NR particles at infinity, using Eq. (21), is then

$dM/dt\approx 8M\Sigma_{\infty}\ \ ({\rm NR\ particles},\ \hat{v}_{\infty}\ll
1).$ (45)

We note that Eq. (45) agrees with Eq. (93) in Ref. Cieślik et al. (2022) for
Maxwell-Jüttner particle temperatures approaching zero and Schwarzschild BHs
[set $\alpha\equiv a/M=0$ and note $\rho_{s,\infty}\equiv\Sigma_{\infty}$ in
Eq. (93)].

It is interesting to observe that in this slow-velocity limit, the above rate
of mass-energy accretion for a 2D razor-thin disk depends only on the surface
density of distant particles and not on their velocity dispersion, in contrast
to accretion from a large 3D cluster, which depends on both the asymptotic
density and the velocity dispersion [see Eq. (27)].

We note that to convert back from geometric to physical units, one simply
multiplies the right-hand sides of Eqs. (43)-(45) by $G^{2}/c$.

####  III.1.4 Massless and Extremely Relativistic Particles

Defining the asymptotic particle surface energy density to be
$\mathcal{E}_{\infty}^{r}\equiv\Sigma_{\infty}E_{\infty}$ and noting that
$J_{min}\approx 3\sqrt{3}ME_{\infty}$ when $E_{\infty}\gg 1$, we can evaluate
Eq. (44) for massless particles, yielding

$\frac{dM}{dt}=6\sqrt{3}M\mathcal{E}^{r}_{\infty},\ \ \ (m=0).$ (46)

We again note that Eq. (46) also applies to ER particles with nonzero $m$. For
such particles the accretion rates may also be written as

$\frac{dM}{dt}=E_{\infty}\frac{dM_{0}}{dt}\approx
6\sqrt{3}M\Sigma_{\infty}E_{\infty}\ \ \ ({\rm ER\
particles},\hat{v}_{\infty}\gg 1).$ (47)

We point out that Eq. (47 agrees with Eqs. (124) and (125) in Ref. Cieślik et
al. (2022) for Maxwell-Jüttner particle temperatures approaching infinity and
Schwarzschild BHs [again set $\alpha\equiv a/M=0$ and note
$\rho_{s,\infty}\equiv\Sigma_{\infty}$ in Eq. (124) and
$\epsilon_{s,\infty}=\Sigma_{\infty}E_{\infty}$ in Eq. (125)].

###  III.2 Loss-Cone Effect

####  III.2.1 Density

Here we treat the scenario whereby collisionless particles in the loss cone
are never replenished once captured, whereby the distribution function may
again be represented by Eq. (31). Using Eq. (39), the surface density is then
given by

$\displaystyle\Sigma{\text{*}}(r)$ $\displaystyle=$ $\displaystyle\int
f(E)d^{2}\hat{p}$ $\displaystyle=$
$\displaystyle\int^{\infty}_{(1-\frac{2M}{r})^{\frac{1}{2}}}f(E)dE\int_{J_{min}(E)}^{J_{max}(E)}\frac{4dJ}{r|v^{\hat{r}}|(1-\frac{2M}{r})^{1/2}}$

Substituting Eqs. (31) and (33) and integrating yields

$\frac{\Sigma(r)}{\Sigma_{\infty}}=\frac{1-\frac{2}{\pi}\arctan\left(\frac{J_{min}/r}{\left[\frac{E_{\infty}^{2}-1+2M/r}{1-2M/r}-(J_{min}/r)^{2})\right]^{1/2}}\right)}{1-2M/r}.$
(49)

For NR particles in Newtonian gravitation, Eq. (49) reduces to

$\frac{\Sigma(r)}{\Sigma_{\infty}}\approx
1-\frac{2}{\pi}\arcsin\left(\frac{4M}{r\hat{v}_{\infty}}\right),\\\ \ \
(\hat{v}_{\infty}\ll 1,\ M/r\ll 1).$ (50)

####  III.2.2 Velocity Dispersion

As all particles again have the same energy $E_{\infty}$, they have the same
velocity and velocity dispersion profile given by Eq. (16).

####  III.2.3 Accretion Rate

Given a loss cone that is unreplenished with particles once they have been
captured, the accretion rate becomes zero in steady state.

##  IV Comparison of Accretion Rates

Comparing the accretion rates of 3D clusters and 2D thin-disks with filled
loss cones is not entirely straightforward, given their different geometries
and defining parameters. The results will depend on the different physical
systems that are compared. To give one example, let us consider a large,
homogeneous spherical cluster of particles with density $\rho_{\infty}$,
isotropic velocity dispersion $v_{\infty}$, radius $R_{c}\gg r_{h}$, and total
mass $M_{c}$. Imagine that it undergoes collapse parallel to the $z$-axis to a
thin pancake in the $x-y$ plane, preserving its surface density along
cylinders. This scenario might mimic one way that thin sheets of collisionless
particles form in the early universe. We will compare the unbound,
collisionless particle accretion rates for the large spherical cluster and the
pancake for NR particles with $v_{\infty}\ll 1$.

The surface density in the pancake is

$\Sigma_{\infty}(r_{\perp})=2\rho_{\infty}R_{c}\left(1-\frac{r^{2}_{\perp}}{R^{2}_{c}}\right)^{1/2},$
(51)

where $r_{\perp}$ is the radius in the $x$-$y$ plane measured from the center
of the pancake. Taking the surface density at $r_{\perp}$ in the central core,
requiring $r_{h}\ll r_{\perp}\ll R_{c}$, whereby the particles remain largely
unperturbed by the central BH, gives

$\Sigma_{\infty}\approx 2\rho_{\infty}R_{c}.$ (52)

The ratio of BH accretion rates in the spherical cluster vs. the thin disk is
then approximated by

$\frac{\dot{M}_{c}}{\dot{M}_{d}}\approx\pi\frac{M}{R_{c}}\frac{1}{v_{\infty}},$
(53)

where $\dot{M}_{c}$ is given by Eq. (27) and $\dot{M}_{d}$ is given by Eq.
(45), substituting Eq. (52) for $\Sigma_{\infty}$. We can evaluate
$v_{\infty}$ if we assume that it is some fraction of the virial value,
whereby

$v_{\infty}^{2}\approx\frac{M_{c}}{R_{c}}.$ (54)

Equation (53) then yields

$\frac{\dot{M}_{c}}{\dot{M}_{d}}\approx\pi\left(\frac{M}{M_{c}}\right)^{1/2}\left(\frac{M}{R_{c}}\right)^{1/2}\ll
1.$ (55)

The strong inequality above for this one extreme example shows the dominance
of disk vs. cluster accretion, thereby demonstrating that the rate of unbound,
collisionless particles can depend significantly on the geometry of the
ambient particle cloud, other parameters being equal.

It is also interesting to compare the collisionless particle rates to the
standard spherical Bondi accretion rate for particles comprising a true fluid.
The Bondi rate for an adiabatic gas with $1\leq\Gamma\leq 5/3$ Bondi (1952)
(see also Shapiro and Teukolsky (1983) for a relativistic treatment) is given
by

$\frac{dM}{dt}=\dot{M}_{b}=4\pi\lambda\frac{M^{2}\rho_{\infty}}{a^{3}_{\infty}},$
(56)

where $a_{\infty}$ is the asymptotic sound speed of the fluid and $\lambda$ is
a constant of order unity that depends on $\Gamma$. We will equate
$a_{\infty}$ to $v_{\infty}$ for comparison below, obtaining

$\frac{\dot{M}_{b}}{\dot{M}_{c}}\approx\frac{\lambda}{4}\frac{1}{v_{\infty}^{2}}\gg
1,$ (57)

and

$\frac{\dot{M}_{b}}{\dot{M}_{d}}\approx\frac{\pi\lambda}{4}\frac{M}{R_{c}}\frac{1}{v_{\infty}^{3}}\approx\frac{\pi\lambda}{4}\left(\frac{M}{M_{c}}\right)\left(\frac{R_{c}}{M_{c}}\right)^{1/2},$
(58)

where we used Eq. (54) to obtain the second equality in Eq. (58). Equation
(57) shows that Bondi accretion dominates collisionless particle accretion for
a large 3D cluster. However Eq. (58) suggests that the ratio for an extended
2D thin disk depends on the particular system, since the first factor in
parentheses on the right-hand side may be much smaller than unity, but the
second factor is much bigger than unity. Once again, spatial geometry counts.

###  IV.1 Applications

There may exist a near “universal” value of the surface density for DM that
spans, within a factor of two, over at least nine (and possibly more) galaxy
magnitudes and across several different Hubble types Donato et al. (2009);
Kormendy and Freeman (2004):

$\Sigma_{\infty}^{u}\approx 140~{}{\rm M_{\odot}~{}pc^{-2}}.$ (59)

The cosmological implications of this observation are not yet resolved, but
Eq. (45) suggests that a thin disk with the surface density
$\Sigma_{\infty}^{u}$ and a central BH would have an accretion rate given by

$\frac{dM^{u}_{d}}{dt}\approx 3.3\times 10^{-2}{\rm
M_{\odot}~{}yr^{-1}}\left(\frac{M}{10^{6}{\rm
M_{\odot}}}\right)\left(\frac{\Sigma_{\infty}^{u}}{140~{}{\rm
M_{\odot}~{}yr^{-1}}}\right).$ (60)

The surface density $\Sigma_{\infty}^{u}$ is also within a factor of two of
estimates of DM in the Galactic neighborhood, where $\rho_{D}\sim 0.008~{}{\rm
M_{\odot}~{}pc^{-3}}$ and $D\sim 8.5~{}{\rm kpc}$ Bovy and Tremaine (2012),
yielding $\Sigma_{\infty}^{Gal}\sim\rho_{D}D\approx 68~{}{\rm
M_{\odot}~{}pc^{-2}}$. Since $\rho_{D}$ may scale as $r^{-1}$ should it obey
an NFW profile Navarro et al. (1997), the Galactic DM surface density
estimated as $\sim\rho_{D}r$ would be constant all the way down to the BH
sphere of influence at $r_{h}$. If it were to reside in a thin disk near the
Galactic Center its accretion rate would be comparable to the universal value
given by Eq. (60) for a BH of mass $4.3\times 10^{6}$ Genzel et al. (2010);
Ghez et al. (2008). If instead it were to occupy a large spherical cluster
outside $r_{h}$ and move with a velocity dispersion of $\sim 100~{}{\rm
km~{}s^{-1}}$ then Eq. (27) suggests it will accrete at a much smaller rate of
$dM^{c}/dt\approx 2.4\times 10^{-7}~{}{\rm M_{\odot}~{}yr^{-1}}$.

By comparison the Bondi value for the baryon accretion rate onto Sgr A* at the
Galactic Center, which is determined from the gas density and temperature
inferred from the diffuse X-ray emission observed by Chandra at $\sim 2$
arcsec ($\sim 0.1$ pc) from the black hole, is $dM^{b}/dt\sim 2\times
10^{-5}~{}{\rm M_{\odot}~{}yr^{-1}}$. In fact, the baryon accretion rate is
believed to be $\sim 10^{-8}~{}{\rm M_{\odot}~{}yr^{-1}}$, or roughly three
orders of magnitude below the Bondi value as determined from polarization
measurements Marrone et al. (2007) and models of the near-horizon accretion
flow and emitted luminosity Shcherbakov and Baganoff (2010); Ressler et al.
(2017). This difference may be due to the angular momentum of the stellar
winds that may be supplying the gas, or possibly to more exotic effects such
as the heating of the gas by DM annihilation in the spike about the BH
Bennewitz et al. (2019).

##  V Summary and Conclusions

We have examined the steady-state density and velocity profiles, and the
associated accretion rates, of collisionless particles (e.g. stars or DM)
moving around a central Schwarzschild black hole in unbound orbits. We
considered two distinct spatial geometries for the particle: an infinite 3D
cluster and a 2D razor-thin disk, both without net angular momentum. We
adopted the same simple monoenergetic, phase-space distribution function for
the particles for both cases, arguing that, though idealized, this assignment
was sufficient to illustrate the features that might distinguish nonrotating
spherical-like and disklike collisionless systems orbiting a black hole. We
treated both a totally isotropic velocity profile at each point and one in
which an empty loss cone is present due to the capture of low-angular momentum
particles that are captured by the BH and not replenished. In all cases the
net angular momentum of the systems was assumed to be zero so that the only
differences were due to the different spatial geometries and velocity
anisotropies adopted.

We found that even in the weak-field region, where $r\gg M$, a mild spike
arises in the locally measured particle density $n(r)$ for the 3D cluster but
that the surface density $\Sigma^{*}$ remains constant with $r$ for the 2D
razor-thin disk. We also found that, at least for one simple application, the
rate of accretion of the disk was much larger than that of the cluster.
However, both rates were much lower than the Bondi accretion rate for a fluid
with a comparable particle density and velocity dispersion (i.e., sound speed)
far from the BH.

While these differences may not be so striking when more realistic phase-space
distribution functions and geometries are considered, the results do suggest
that the spatial distribution of particles around a black hole is a feature
that affects the resulting steady-state particle profiles and accretion rates
significantly. So this is just one other factor that must be accounted for, in
addition to knowing what the nature of the particles are (e.g., collisionless
or collisional fluid matter, or mildly collisional by virtue of self-
interactions and/or annihilations) and the global properties of their
distributions (bound or unbound, with or without net angular momentum, subject
or not to gravitational intruders, etc) in assessing their profiles and
capture rates about a black hole.

Acknowledgments: It is a pleasure to thank Yuk Tung Liu and Olivier Sarbach
for several useful discussions. This paper was supported in part by NSF Grants
No. PHY-2006066 and No. PHY-2308242 and NASA Grant 80NSSC17K0070 to the
University of Illinois at Urbana-Champaign.

##  Appendix A Cluster Accretion Rates: Alternative Derivations

###  A.1 Alternative Derivation 1

Here we obtain the accretion rate for an unbound, monoenergetic, nonrotating
cluster as measured by a locally static observer by calculating the total
inward particle flux across a sphere of radius $r$ $\times$ the area of the
sphere $\times$ the fraction of these particles that move within the loss cone
and are thus captured:

$\frac{dN_{tot}}{d\tau}=\left(\frac{1}{4}n\hat{v}\right)\times\left(4\pi
r^{2}\right)\times\cal{P},$ (61)

where $\hat{v}$ is the magnitude of the isotropic 3-velocity at $r$ given by
Eq. (16) and $\cal{P}$ is the fraction of particles captured. For $\cal{P}$ we
have

${\cal{P}}=\frac{\int_{0}^{J_{min}(E_{\infty})}dJJ}{\int_{0}^{J_{max}(E_{\infty})}dJJ}=\frac{J_{min}^{2}(E_{\infty})}{J_{max}^{2}(E_{\infty})},$
(62)

while the density $n$ is given by Eq. (10) and $J_{max}(E)$ is given by Eq.
(33). Assembling the factors in Eq. (61) then yields Eq. (22) for the locally
measured accretion rate at $r$, $dN_{tot}/d\tau$, from which, using Eq. (23),
the rates measured by a distant observer, Eqs. (24) - (26) for $dN_{tot}/dt$,
$dM_{0}/dt$ and $dM/dt$, respectively, follow immediately.

###  A.2 Alternative Derivation 2

Here we provide yet another derivation of the accretion rate found above. The
maximum impact parameter $b_{max}$ for a particle of energy $E_{\infty}$
falling inward from infinity to be captured by the BH is given by

$b_{max}^{2}=\frac{J_{min}^{2}(E_{\infty})}{E_{\infty}^{2}-1},$ (63)

[see, e.g. Ref. Shapiro and Teukolsky (1983), Eq. (12.4.35)], with the typo
corrected for the missing square on $E_{\infty}$), whereby the capture cross
section is

$\sigma_{cap}=\pi b_{max}^{2}=\frac{\pi
J_{min}^{2}(E_{\infty})}{E_{\infty}^{2}-1}.$ (64)

So far from the BH the accretion rate is obtained as the intensity of
particles for an isotropic distribution $\times$ the area of a large sphere
about the BH $\times$ the solid angle within which a particle is captured:

$\frac{dN_{tot}}{dt}=\left(\frac{n_{\infty}\hat{v}_{\infty}}{4\pi}\right)(4\pi
r^{2})(\Delta\Omega_{cap}),$ (65)

where using Eq. (16) at $r\rightarrow\infty$, we have

$\hat{v}_{\infty}=\frac{(E_{\infty}^{2}-1)^{1/2}}{E_{\infty}}$ (66)

and where

$\Delta\Omega_{cap}=\frac{\sigma_{cap}}{r^{2}}.$ (67)

Assembling the factors in Eq. (65) again yields Eq. (24) for $dN_{tot}/dt$,
from which Eq. (22) for $dN_{tot}/d\tau$ and Eqs. (25) - (26) for $dM_{0}/dt$
and $dM/dt$, respectively, again follow immediately.

## References

  * Peebles (1972) P. J. E. Peebles, General Relativity and Gravitation 3, 63 (1972).
  * Young (1977) P. J. Young, Astrophys. J. 217, 287 (1977).
  * Gondolo and Silk (1999) P. Gondolo and J. Silk, Physical Review Letters 83, 1719 (1999).
  * Bahcall and Wolf (1976) J. N. Bahcall and R. A. Wolf, Astrophys. J. 209, 214 (1976).
  * Frank and Rees (1976) J. Frank and M. J. Rees, Mon. Not. R. Astro. Soc. 176, 633 (1976).
  * Lightman and Shapiro (1977) A. P. Lightman and S. L. Shapiro, Astrophys. J. 211, 244 (1977).
  * Shapiro and Marchant (1978) S. L. Shapiro and A. B. Marchant, Astrophys. J. 225, 603 (1978).
  * Cohn and Kulsrud (1978) H. Cohn and R. M. Kulsrud, Astrophys. J. 226, 1087 (1978).
  * Gnedin and Primack (2004) O. Y. Gnedin and J. R. Primack, Physical Review Letters 93, 061302 (2004).
  * Merritt (2004) D. Merritt, Physical Review Letters 92, 201304 (2004).
  * Shapiro and Heggie (2022) S. L. Shapiro and D. C. Heggie, Physical Review D 106, 043018 (2022).
  * Shapiro and Paschalidis (2014) S. L. Shapiro and V. Paschalidis, Phys. Rev. D 89, 023506 (2014).
  * Shapiro (2018) S. L. Shapiro, Phys. Rev. D 98, 023021 (2018), eprint 1809.02618.
  * Luca and Khoury (2023) V. D. Luca and J. Khoury, Journal of Cosmology and Astroparticle Physics 2023, 048 (2023).
  * Vasiliev (2007) E. Vasiliev, Phys. Rev. D 76, 103532 (2007).
  * Shapiro and Shelton (2016) S. L. Shapiro and J. Shelton, Phys. Rev. D 93, 123510 (2016).
  * Zeldovich and Novikov (1971) Y. B. Zeldovich and I. D. Novikov, _Relativistic Astrophysics. Vol.1: Stars and Relativity_ (University of Chicago Press, Chicago, 1971).
  * Shapiro and Teukolsky (1983) S. L. Shapiro and S. A. Teukolsky, _Black Holes, White Dwarfs, and Neutron Stars: The Physics of Compact Objects_ (Wiley, New York, 1983).
  * Shapiro and Teukolsky (1985) S. L. Shapiro and S. A. Teukolsky, Astrophys. J. 298, 34 (1985).
  * Gamboa et al. (2021) A. Gamboa, C. Gabarrete, P. Dominguez-Fernandez, D. Nunez, and O. Sarbach, Physical Review D 104, 083001 (2021).
  * Zel’dovich (1970) Y. B. Zel’dovich, Astron.Astrophys. 5, 84 (1970).
  * Kulkarni et al. (2022) M. Kulkarni, E. Visbal, G. L. Bryan, and X. Li, Astrophys. J. Lett. 941, L18 (2022).
  * Rioseco and Sarbach (2017) P. Rioseco and O. Sarbach, Journal of Physics: Conference Series 831, 012009 (2017).
  * Rioseco and Sarbach (2017) P. Rioseco and O. Sarbach, Classical and Quantum Gravity 34, 095007 (2017), eprint 1611.02389.
  * Cieślik et al. (2022) A. Cieślik, P. Mach, and A. Odrzywołek, Phys. Rev. D 106, 104056 (2022), eprint 2208.04218.
  * Daylan et al. (2016) T. Daylan, D. P. Finkbeiner, D. Hooper, T. Linden, S. K. N. Portillo, N. L. Rodd, and T. R. Slatyer, Physics of the Dark Universe 12, 1 (2016).
  * Cholis et al. (2022) I. Cholis, Y.-M. Zhong, S. D. McDermott, and J. P. Surdutovich, Phys. Rev. D 105, 103023 (2022), eprint 2112.09706.
  * Leane et al. (2022) R. K. Leane, S. Shin, L. Yang, G. Adhikari, H. Alhazmi, T. Aramaki, D. Baxter, F. Calore, R. Caputo, I. Cholis, et al., arXiv:2203.06859 (2022).
  * Fields et al. (2014) B. D. Fields, S. L. Shapiro, and J. Shelton, Physical Review Letters 113, 151302 (2014).
  * Wanders et al. (2015) M. Wanders, G. Bertone, M. Volonteri, and C. Weniger, JCAP 4, 004 (2015).
  * Chiang et al. (2020) B. T. Chiang, S. L. Shapiro, and J. Shelton, Phys. Rev. D 102, 023030 (2020), URL https://link.aps.org/doi/10.1103/PhysRevD.102.023030.
  * Bondi (1952) H. Bondi, Mon. Not. R. Astro. Soc. 112, 195 (1952).
  * Donato et al. (2009) F. Donato, G. Gentile, P. Salucci, C. Frigerio Martins, M. I. Wilkinson, G. Gilmore, E. K. Grebel, A. Koch, and R. Wyse, Monthly Notices of the Royal Astronomical Society 397, 1169 (2009).
  * Kormendy and Freeman (2004) J. Kormendy and K. C. Freeman, in _Dark Matter in Galaxies_ , edited by S. Ryder, D. Pisano, M. Walker, and K. Freeman (ASP, San Francisco, 2004), vol. 220, p. 377.
  * Bovy and Tremaine (2012) J. Bovy and S. Tremaine, Astrophys. J. 756, 89 (2012).
  * Navarro et al. (1997) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 490, 493 (1997).
  * Genzel et al. (2010) R. Genzel, F. Eisenhauer, and S. Gillessen, Reviews of Modern Physics 82, 3121 (2010).
  * Ghez et al. (2008) A. M. Ghez, S. Salim, N. N. Weinberg, J. R. Lu, T. Do, J. K. Dunn, K. Matthews, M. R. Morris, S. Yelda, E. E. Becklin, et al., Astrophys. J. 689, 1044 (2008).
  * Marrone et al. (2007) D. P. Marrone, J. M. Moran, J.-H. Zhao, and R. Rao, Astrophys. J. 654, L57 (2007), eprint astro-ph/0611791.
  * Shcherbakov and Baganoff (2010) R. V. Shcherbakov and F. K. Baganoff, Astrophys. J. 716, 504 (2010), eprint 1004.0702.
  * Ressler et al. (2017) S. M. Ressler, A. Tchekhovskoy, E. Quataert, and C. F. Gammie, Mon. Not. R. Astro. Soc. 467, 3604 (2017), eprint 1611.09365.
  * Bennewitz et al. (2019) E. R. Bennewitz, C. Gaidau, T. W. Baumgarte, and S. L. Shapiro, Mon. Not. R. Astro. Soc. 490, 3414 (2019), eprint 1907.00015.


# Stochastic Code Generation

Swapnil Sharma  Nikita Anand  Kranthi Kiran GV

###### Abstract

Large language models pre-trained for code generation can generate high-
quality short code, but often struggle with generating coherent long code and
understanding higher-level or system-level specifications. This issue is also
observed in language modeling for long text generation, and one proposed
solution is the use of a latent stochastic process. This approach involves
generating a document plan and then producing text that is consistent with it.

In this study, we investigate whether this technique can be applied to code
generation to improve coherence. We base our proposed encoder and decoder on
the pre-trained GPT-2 based CodeParrot model and utilize the APPS dataset for
training. We evaluate our results using the HumanEval benchmark and observe
that the modified Time Control model performs similarly to CodeParrot on this
evaluation.

Machine Learning, ICML

##  1 Introduction

Large pre-trained language models have gained immense popularity recently.
Large language models have achieved impressive results in various natural
language processing tasks, including language translation, summarization, and
text generation. These models are trained on a large corpus of text data and
are able to capture complex patterns and relationships in the data. One
impactful application of these models is code generation.

Code generation is the process of automatically generating computer programs
based on natural language descriptions or higher-level specifications.
Multiple large scale models like PyMT5 (Clement et al., 2020), Codebert (Feng
et al., 2020a) and Codex (Chen et al., 2021) have been trained for the task
and are used to aid software developers. They are used in integrated
development environments for code completion, to convert natural language
documentation (docstrings) into code and vice versa.

These models produce high quality short pieces of code but as pointed out by
Chen et al., 2021, they tend to struggle with long and higher-level or system-
level specifications. An experiment was run on Codex that indicates
performance degradation as docstring length increases.

As stated by Wang et al., 2022, this problem also occurs in text generation
where models tend to meander or lose coherency when producing long text. They
propose Time Control, a language model that learns a mapping of how text
changes in a document with the changes in the document plan generated via a
stochastic process. Text is then generated to be consistent with this document
plan. This helps the model learn structure and leads to locally and globally
more coherent text.

In this work, we adapt the stochastic process for the code generation task. We
do this by basing the proposed encoder-decoder architecture on the pre-trained
CodeParrot model (Tunstall et al., 2022) which was trained using 180 gigabytes
of Python code. We use the Automated Programming Progress Standard dataset
Hendrycks et al., 2021 to generate the latent plan and train the encoder and
fine-tune the decoder.

We compare our model’s (Modified Time Control) performance to the regular pre-
trained CodeParrot model on the HumanEval (Chen et al., 2021) Evaluation
framework and see that they perform similarly. We also compare the results
manually to get a better understanding of when it works well and when it
doesn’t.

##  2 Related Work

Autoregressive models, which are commonly used for generating text, may
struggle to produce long and coherent text because they lack the ability to
model text structure and dynamics (Wang et al., 2022; Lin et al., 2021). This
can result in globally incoherent text as the models are unable to effectively
plan and anticipate long context in the text. When forced to generate longer
texts, this incoherence is often worsened as the models struggle to extend
beyond their expected text end point. To address this issue, previous research
has employed planning-based approaches in an attempt to generate globally
coherent text (Puduppully et al., 2019; Kiddon et al., 2016).

Program induction is a machine learning technique in which a model generates
program outputs from a latent program representation (Zaremba & Sutskever,
2014; Pierrot et al., 2019). Program synthesis involves generating a program
from a natural language specification (Feng et al., 2020b). Both approaches
have been improved upon through the incorporation of inductive biases and the
use of recurrence and abstract syntax trees.

The Codex paper by Chen et al., 2021 introduces a large language model called
Codex that is finetuned on code from GitHub and is able to generate Python
code. The model was evaluated on the HumanEval dataset (Chen et al., 2021),
which measures functional correctness for synthesizing programs from
docstrings, and was found to solve 28.8% of the problems, outperforming the
GPT-3 model which solved 0%. The limitations of Codex include difficulty with
docstrings describing long chains of operations and with binding operations to
variables.

##  3 Method

In this section we describe the parts of Language Modeling via Stochastic
Processes (Wang et al., 2022) that we use to run our experiments. In the
following Experiments section, we describe the changes made to adapt this
method to code generation in detail.

###  3.1 Architecture

We adapt the same structure as Time Control but make modifications to run for
code generation. The proposed solution is a modified encoder-decoder model
built on top of CodeParrot.

#### Encoder

The encoder in this study is a nonlinear mapping that transforms raw input
data into a low-dimensional latent space. The goal of the encoder is to map
high-dimensional sequential data into a latent space that follows a stochastic
process of interest, similar to the Brownian bridge process described by Wang
et al., 2022.

As demonstrated by Wang et al., 2022, the density of a Brownian bridge process
between an arbitrary start point, $z_{0}$, at time $t=0$ and end point,
$z_{T}$, at time $t=T$, is given by:

$p(z_{t}|z_{0},z_{T})=\mathcal{N}\left(\left(1-\frac{t}{T}\right)z_{0}+\frac{t}{T}z_{T},\frac{t(T-t)}{T}\right)$

This equation acts as a noisy interpolation between the start and end points
of the trajectory, with the center being the most noisy. To ensure that every
triplet of observations, ($x_{1},x_{2},x_{3}$), follows the Brownian bridge
transition density, we use a contrastive objective.

#### Decoder

The decoding component of the system is responsible for generating code from
latent plans. To do this, we first map all lines in the training dataset to
our learned latent space using the pretrained encoder, $f_{\theta}$. This
produces a Brownian bridge trajectory of sentence-level latent codes,
$(z_{0},\cdots,z_{t},\cdots,z_{T})$, for each entry in the dataset. Rather
than training a decoder from scratch, we fine-tune the existing model,
CodeParrot, to generate code conditioned on past context and the latent
structure.

To fine-tune CodeParrot, we use a standard auto-regressive language model that
has been modified as follows. Given a document with $W$ tokens and $T$
sentences used to train the decoder, $x_{1}\cdots x_{W}$, we obtain
embeddings, $z_{1}\cdots z_{3}$, for each sentence using the encoder,
$f_{\theta}$. At time $t$, the decoder must predict $x_{t}$ using all past
tokens, $x_{<t}$, as well as the sentence embedding, $z_{s_{t}}$, where the
index $s_{t}\in[T]$ maps each token to its corresponding sentence. This
process can be viewed as a reconstruction objective, as the identity of
$x_{t}$ is encoded in $z_{s_{t}}$.

⬇

\-----Question\-----

Polycarp analyzes the prices of the new berPhone. At his disposal are the
prices for n last days: a_1, a_2, ... , a_n, where a_i is the price of
berPhone on the day i. Polycarp considers the price on the day i to be bad if
later (that is, a day with a greater number) berPhone was sold at a lower
price. For example, if n=6 and a=[3, 9, 4, 6, 7, 5], then the number of days
with a bad price is 3 - these are days 2 (a_2=9), 4 (a_4=6) and 5 (a_5=7).
Print the number of days with a bad price. You have to answer t independent
data sets.

\-----Input\-----

The first line contains an integer t (1 <= t <= 10000$) - the number of sets
of input data in the test. Input data sets must be processed independently,
one after another. Each input data set consists of two lines. The first line
contains an integer n (1 <= n <= 150000) - the number of days. The second line
contains n integers a_1, a_2, ... , a_n (1 <= a_i <= 10^6), where a_i is the
price on the i-th day. It is guaranteed that the sum of n over all data sets
in the test does not exceed 150000.

\-----Output\-----

Print t integers, the j-th of which should be equal to the number of days with
a bad price in the j-th input data set.

\-----Example\-----

Input 5 6 3 9 4 6 7 5 1 1000000 2 2 1 10 31 41 59 26 53 58 97 93 23 84 7 3 2 1
2 3 4 5

Output 3 0 1 8 2

Listing 1: Sample Introductory Level Question from APPS
Level | Total Count | Test Set
---|---|---
Introductory | 3639 | 1000
Interview | 5000 | 3000
Competition | 1361 | 1000
Table 1: Problems of varying difficulty levels in the APPS dataset

###  3.2 Dataset

The Automated Programming Progress Standard (APPS) dataset introduced by
Hendrycks et al., 2021 consists of coding problems collected from multiple
open access websites like Codeforces, Kattis etc. It contains 10,000 problems
of varying levels; introductory, interview level and competition level. The
introductory problems are ones whose solutions do not require any complicated
algorithms, and can be solved by people with a year or two of programming
experience. The interview level questions are ones that are used to test
technical proficiency and require non-trivial solutions. They can typically
involve the use of data structures like trees, linked lists etc. The
competition level questions are ones that require advanced collegiate level
knowledge.

The latent structure was labelled using an open-source tool called tree-
sitter111https://github.com/tree-sitter/tree-sitter. This tool is a parser
that builds a syntax tree for a source file. It works with any programming
language, is robust enough that it can work with minor errors and can be
embedded in any application because it is dependency free.

###  3.3 Choice of Model

Some of the most popular language models like Codex which powers the industry
product Github Copilot, are not open-source and thus unusable in this work.
Our choice of pre-trained model for code generation is largely based on what
was open sourced and popular on the Hugging Face library (Wolf et al., 2019).

For our experiments we use CodeParrot (Tunstall et al., 2022) which is a GPT-2
(Radford et al., 2019) based model trained to generate Python code. It was
downloaded 22,416 times from Hugging Face in the last month. CodeParrot is an
auto-regressive Left-to-Right Language Model that predict the probability of a
token given the previous tokens. It has 1.5 billion parameters and was trained
on 180 GB of Python code. The Time Control model (Wang et al., 2022) we
reference is also based on GPT-2 which lead us to believe this model would be
a good fit for our experiments.

##  4 Experiments

The purpose of this study is to evaluate the capabilities of Time Control in
capturing code structure and generating more robust code. Specifically, we aim
to determine whether the generated code is able to compile and solve the
problem stated in the prompt. To assess Time Control, we use a modified
version with latent dimensions (d = 32).

The encoder architecture for this evaluation consists of a frozen, pretrained
CodeParrot and a trainable multi-layer perceptron (MLP). We extract the last
layer hidden state of CodeParrot that corresponds to the end-of-sentence
token, and train the four-layer MLP on top of this hidden state. The MLP has
intermediate rectified linear unit activations and is trained using stochastic
gradient descent with a learning rate of 1e-4 and momentum of 0.9.

For the decoder, we also use CodeParrot as the base model. Similar to previous
work (Wang et al., 2022), we modify CodeParrot to pass the encoder output
after appending it to positional embeddings in the transformer layer. You can
see one of the good samples generated in Figure 1 and some of the bad ones 2.

Table 2: Model accuracy with baseline CodeParrot.
Model | Pass@1 | Pass@10 | Pass@100
---|---|---|---
TimeControl | 3.7% | 6.58% | 12.73%
CodeParrot 110M | 3.80% | 6.57% | 12.78%
CodeParrot 1.5B | 3.58% | 8.03% | 14.96%

#### Datasets

As the dataset for this study, we use APPS, which provides sufficient length
code to fully utilize our available resources. We modify the dataset so that
each entry contains five sections (question, solution, class statement, def
statement, import statement) marked with section identifier tokens.
Specifically, each entry is represented as,

        [QUESTION]
            <text>

        [SOLUTION]
            [CLASS_STATEMENT]

            [DEF_STATEMENT]

            [IMPORT_STATEMENT]

No additional tokens are added in this dataset.

To label the code in the dataset with the appropriate section identifiers, we
utilize tree-sitter to mark each line of code with tree-sitter specific
labels, and then retain only the relevant labels while discarding the rest.

⬇

[QUESTION]:

In this problem we assume the Earth to be a completely round ball and its
surface a perfect sphere. The length of the equator and any meridian is
considered to be exactly 40 000 kilometers. Thus, travelling from North Pole
to South Pole or vice versa takes exactly 20 000 kilometers.

Limak, a polar bear, lives on the North Pole. Close to the New Year, he helps
somebody with delivering packages all around the world. Instead of coordinates
of places to visit, Limak got a description how he should move, assuming that
he starts from the North Pole. The description consists of n parts. In the
i-th part of his journey, Limak should move t_{i} kilometers in the direction
represented by a string dir_{i} that is one of: "North", "South", "West",
"East".

Limak isn’t sure whether the description is valid. You must help him to check
the following conditions: If at any moment of time (before any of the
instructions or while performing one of them) Limak is on the North Pole, he
can move only to the South . If at any moment of time (before any of the
instructions or while performing one of them) Limak is on the South Pole, he
can move only to the North .

[SOLUTION]:

def main():

latitude = 0

n, = read()

for i in range(n):

l, d = read(str)

l = int(l)

if latitude == 0:

if d != "South":

return "NO"

if latitude == 20000:

if d != "North":

return "NO"

if d == "South":

latitude += l

if not (0 <= latitude <= 20000):

return "NO"

if latitude != 0:

return "NO"

return "YES"

def read(callback=int):

return list(map(callback, input().strip().split()))

def write(value, end="\n"):

if value is None: return

try:

if not isinstance(value, str):

value = " ".join(map(str, value))

except:

pass

print(value, end=end)

write(main())

Figure 1: Problem Statement and Solution generated by Modified Time Control

⬇

[QUESTION]

Given an array of is, you can find a array where each element has equal
letters.

[SOLUTION]

a = list(str(input()))

for i in range(1,len(a)):

sh = a[i].swapcase()

tell = sh[’what’]

...

x,n = list(a.swapcase()) a.sort() pivot = """" pivot_pos =
pivot_pos__starting_point__starting_point__starting_point__starting_point__len()
if pivot_pos: pivot_pos += 1 pivot_pos =
pivot_pos__starting_point__starting_point__end() gcdrom = pivot ocdr = pivot
while True: a,r = pivot_pos__starting_point__starting_point__starting_point...

Figure 2: Sample from Modified Time Control with incomprehensible question and
output

###  4.1 Evaluation

To evaluate the model, we use the HumanEval framework. HumanEval contains 164
handwritten programming problems; this is done to ensure that the models could
not have seen the problems before as most are trained on open access code from
Github which contain solutions to problems from the web like the APPS dataset.

We don’t use metrics like BLEU Score that are used for natural language
because it neglects syntactic and semantic features of code (Ren et al., 2020)
and doesn’t work well for evaluatinh code generation.

Our results, shown in Table 2, indicate that modified Time Control performs
similarly to CodeParrot on the HumanEval dataset. As described by Chen et al.,
2021, ”Pass@k can be interpreted as the result of evaluating the best out of k
samples, where the best sample is picked by an oracle with prior knowledge of
the unit tests”. These findings suggest that Modified Time Control is a
promising approach, as we were able to maintain code quality with our limited
computational resources. We hope that with longer per-sentence lengths and
additional computing resources, we will be able to further improve the
performance of our model.

##  5 Conclusion

The purpose of this study was to investigate the potential of Time Control to
improve code generation through the use of a latent stochastic process that
generates output with learned implicit structure. Our results suggest that
modified Time Control is capable of learning to map programming code to smooth
Brownian bridge trajectories, and in some cases, we observed improved
performance in the code generation process. These findings are encouraging for
the use of Time Control in code generation.

To further explore the potential of Time Control, we suggest using larger
datasets containing multiple files in a project as dataset entries. We believe
that this approach may allow for the learning of code-level architecture. In
addition, the use of enhanced datasets and additional computing resources may
lead to even more improved results with modified Time Control.

##  6 Future Work

Due to computational constraints, we were unable to utilize approximately 70%
of our training set and had to limit the decoder to less than 750 tokens. It
is possible that the use of the full dataset may have produced different
results.

In an ideal scenario, we would have preferred to conduct these experiments
using a more structured dataset, such as programs with the model-view-
controller software architecture pattern. Such a dataset would provide higher-
level specifications and have a well-defined structure that may be
particularly well-suited to benefit from a latent process like the one used in
this study.

## References

  * Chen et al. (2021) Chen, M., Tworek, J., Jun, H., Yuan, Q., de Oliveira Pinto, H. P., Kaplan, J., Edwards, H., Burda, Y., Joseph, N., Brockman, G., Ray, A., Puri, R., Krueger, G., Petrov, M., Khlaaf, H., Sastry, G., Mishkin, P., Chan, B., Gray, S., Ryder, N., Pavlov, M., Power, A., Kaiser, L., Bavarian, M., Winter, C., Tillet, P., Such, F. P., Cummings, D., Plappert, M., Chantzis, F., Barnes, E., Herbert-Voss, A., Guss, W. H., Nichol, A., Paino, A., Tezak, N., Tang, J., Babuschkin, I., Balaji, S., Jain, S., Saunders, W., Hesse, C., Carr, A. N., Leike, J., Achiam, J., Misra, V., Morikawa, E., Radford, A., Knight, M., Brundage, M., Murati, M., Mayer, K., Welinder, P., McGrew, B., Amodei, D., McCandlish, S., Sutskever, I., and Zaremba, W.  Evaluating large language models trained on code.  2021\.
  * Clement et al. (2020) Clement, C., Drain, D., Timcheck, J., Svyatkovskiy, A., and Sundaresan, N.  PyMT5: multi-mode translation of natural language and python code with transformers.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pp. 9052–9065, Online, November 2020\. Association for Computational Linguistics.  doi: 10.18653/v1/2020.emnlp-main.728.  URL https://aclanthology.org/2020.emnlp-main.728.
  * Feng et al. (2020a) Feng, Z., Guo, D., Tang, D., Duan, N., Feng, X., Gong, M., Shou, L., Qin, B., Liu, T., Jiang, D., and Zhou, M.  Codebert: A pre-trained model for programming and natural languages.  _CoRR_ , abs/2002.08155, 2020a.  URL https://arxiv.org/abs/2002.08155.
  * Feng et al. (2020b) Feng, Z., Guo, D., Tang, D., Duan, N., Feng, X., Gong, M., Shou, L., Qin, B., Liu, T., Jiang, D., et al.  Codebert: A pre-trained model for programming and natural languages.  In _Findings of the Association for Computational Linguistics: EMNLP 2020_ , pp. 1536–1547, 2020b.
  * Hendrycks et al. (2021) Hendrycks, D., Basart, S., Kadavath, S., Mazeika, M., Arora, A., Guo, E., Burns, C., Puranik, S., He, H., Song, D., and Steinhardt, J.  Measuring coding challenge competence with apps.  _NeurIPS_ , 2021.
  * Kiddon et al. (2016) Kiddon, C., Zettlemoyer, L., and Choi, Y.  Globally coherent text generation with neural checklist models.  In _Proceedings of the 2016 conference on empirical methods in natural language processing_ , pp. 329–339, 2016.
  * Lin et al. (2021) Lin, C.-C., Jaech, A., Li, X., Gormley, M. R., and Eisner, J.  Limitations of autoregressive models and their alternatives.  In _Proceedings of the 2021 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pp. 5147–5173, 2021.
  * Pierrot et al. (2019) Pierrot, T., Ligner, G., Reed, S. E., Sigaud, O., Perrin, N., Laterre, A., Kas, D., Beguir, K., and de Freitas, N.  Learning compositional neural programs with recursive tree search and planning.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Puduppully et al. (2019) Puduppully, R., Dong, L., and Lapata, M.  Data-to-text generation with content selection and planning.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 33, pp. 6908–6915, 2019.
  * Radford et al. (2019) Radford, A., Wu, J., Child, R., Luan, D., Amodei, D., and Sutskever, I.  Language models are unsupervised multitask learners.  2019\.
  * Ren et al. (2020) Ren, S., Guo, D., Lu, S., Zhou, L., Liu, S., Tang, D., Sundaresan, N., Zhou, M., Blanco, A., and Ma, S.  Codebleu: a method for automatic evaluation of code synthesis.  _CoRR_ , abs/2009.10297, 2020.  URL https://arxiv.org/abs/2009.10297.
  * Tunstall et al. (2022) Tunstall, L., von Werra, L., and Wolf, T.  _Natural Language Processing with Transformers: Building Language Applications with Hugging Face_.  O’Reilly Media, 2022.  ISBN 9781098103248.  URL https://books.google.com/books?id=_0uezgEACAAJ.
  * Wang et al. (2022) Wang, R. E., Durmus, E., Goodman, N., and Hashimoto, T.  Language modeling via stochastic processes.  In _International Conference on Learning Representations_ , 2022.  URL https://openreview.net/forum?id=pMQwKL1yctf.
  * Wolf et al. (2019) Wolf, T., Debut, L., Sanh, V., Chaumond, J., Delangue, C., Moi, A., Cistac, P., Rault, T., Louf, R., Funtowicz, M., and Brew, J.  Huggingface’s transformers: State-of-the-art natural language processing.  _CoRR_ , abs/1910.03771, 2019.  URL http://arxiv.org/abs/1910.03771.
  * Zaremba & Sutskever (2014) Zaremba, W. and Sutskever, I.  Learning to execute.  _arXiv preprint arXiv:1410.4615_ , 2014.