# Supermassive black hole formation in the initial collapse
of axion dark matter

Pierre Sikivie and Yuxin Zhao  Department of Physics, University of Florida,
Gainesville, FL 32611, USA

(July 14, 2024)

###### Abstract

Axion dark matter thermalizes by gravitational self-interactions and forms a
Bose-Einstein condensate. We show that the rethermalization of the axion fluid
during the initial collapse of large scale overdensities near cosmic dawn
transports angular momentum outward sufficientlly fast that black holes form
with masses ranging from approximately $10^{5}$ to a few times
$10^{10}~{}M_{\odot}$.

###### pacs:

95.33.+d

It has long been established Korm  that most galaxies have supermassive black
holes at their centers with masses ranging from approximately
$10^{5}M_{\odot}$ to a few times $10^{10}M_{\odot}$. The event horizons of the
black holes at the center of the Milky Way and the ceter of the large galaxy
M87 have been imaged using very long baseline interferometry EHT . Active
galactic nuclei (AGN) are understood to be powered by accretion onto
supermasive black holes AGN . Recently, low frequency gravitational waves
consistent with emission from supermassive black hole mergers have been
detected by pulsar timing arrays PTA .

How the supermassive black holes form has been an enduring puzzle. The main
impediement to their formation is conservation of angular momentum. If, for
example, a $10^{8}~{}M_{\odot}$ black hole condenses out of a region of
density say $10^{-23}$ gr/cc, which is a value typical of galactic disks, the
material forming the black hole must shrink in all directions by eight orders
of magnitude. Angular momentum conservation makes this difficlt by introducing
a distance of closest approach to the black hole. Angular momentum can be
transported outward if the material falling toward the black hole has
viscosity but in that case the material heats up and acquires pressure
opposing its compression. See ref. bhfrev  for a review of the issues involved
in supermassive black hole formation.

To sidestep the problems that arise when only conventional physics is
involved, it has been proposed that the supermassive black holes form as a
result of the gravothermal collapse of overdenisties of dark matter with very
strong self-interactions SS  or by accretion onto dark stars, i.e stars that
are powered by dark matter annihilation KF . In such scenarios, seed black
holes form that have masses of order $10^{5}-10^{6}~{}M_{\odot}$. The larger
supermassive black holes with mass $10^{9}-10^{10}~{}M_{\odot}$ are then
supposed to be the result of accretion onto, and mergers of, the seed black
holes. These proposals have been challenged by the recent discovery, using the
James Webb Space Telescope, of powerful AGN near cosmic dawn, i.e at redshifts
$z\sim 10$ JWST . There is very little time to grow the black holes powering
the AGN observed at $z\sim 10$. Finally, one may contemplate the possibility
that the supermassive black holes are primordial in nature, i.e. that they
formed long before cosmic dawn. This proposal runs afoul of constraints from
cosmic microwave background observations, although a recent paper HoSt
indicates how it may still be viable.

The purpose of our paper is to show that supermassive black holes form
naturally near cosmic dawn if the dark matter is axions or axion-like
particles. No additional assumptions are required. The crucial step is to
recognize that cold dark matter axions thermalize by their gravitational self-
interactions CABEC ; Erken . Indeed, as explained in more detail below, their
density fluctuations are generically large $(\delta\rho\sim\rho)$ and
correlated over long distances. When the axions thermalize they form a Bose-
Einstein condensate, meaning that most axions go to the lowest energy state
available to them through the thermalizing interactions. When an axion
overdensity collapses near cosmic dawn, the gravitional self-interactions
among the axions produce a long range viscosity that causes outward transport
of angular momentum. The heat produced in the axion case flows into the
thermal distribution accompanying the condensate whereas the condensate,
containing most of the axions, stays in the lowest energy particle state
available. That state is one of rigid rotation where most of the angular
momentum resides far from the central overdensity which may then perhaps
collapse into a black hole.

First let us indicate how axions differ from the other cold dark matter
candidates, such as weakly interacting massive particles (WIMPs) and sterile
neutrinos, which we will refer to collectively as “ordinaray CDM”. Before
density perturbations go non-linear and multi-streaming begins, in particular
during the recombination era, the state of ordinary CDM is customarily given
by a mass density $\rho(\vec{x},t)$ and a velocity field $\vec{v}(\vec{x},t)$.
In linear order of perturbation thory, axions behave in the same way as
ordinary CDM on length scales longer than their Jeans length Zel . When the
axion fluid has mass density $\rho(\vec{x},t)$ and velocity field
$\vec{v}(\vec{x},t)$, all axions are in the particle state of wavefunction

$\Psi(\vec{x},t)=\sqrt{\rho(\vec{x},t)\over Nm}e^{i\beta(\vec{x},t)}$ (1)

where $m$ is the axion mass, $N$ is the number of axions and
$\beta(\vec{x},t)$ is such that $\vec{v}(\vec{x},t)={\hbar\over
m}\vec{\nabla}\beta(\vec{x},t)$. For both ordinary CDM and axions, the above
description is exact only when the dark matter has zero velocity dispersion.
In reality both axions and ordinary CDM have a finite velocity dispersion. In
the case of ordinary CDM this distinction is irrelevant but it is important in
the case of axions Erken . WIMPs of mass 100 GeV have primordial velocity
dispersion $\delta v_{W}\sim 10^{-12}~{}c~{}a(t_{0})/a(t)$ and hence
correlation length $\ell_{W}\sim 2~{}\mu$m $a(t)/a(t_{0})$, where $a(t)$ is
the scale factor, $t$ is cosmic time and $t_{0}$ the present age. Sterile
neutrinos with mass a few keV have primordial velocity dispersion $\delta
v_{\nu}\sim 10^{-8}~{}c~{}a(t_{0})/a(t)$ and hence $\ell_{\nu}\sim$ cm
$a(t)/a(t_{0})$. The correlation lenghts of orinary CDM are far too short to
be relevant in large scale structure formation. Axions are different.

The correlation length and velocity dispersion of cold dark matter axions are
set by the horizon when the axion mass turns on during the QCD phase
transition:

$\ell(t)\sim ct_{1}{a(t_{1})\over a(t)}~{}~{}~{}~{},~{}~{}~{}~{}~{}\delta
v(t)\sim{\hbar\over mct_{1}}{a(t_{1})\over a(t)}$ (2)

where $t_{1}\simeq 3.6\cdot 10^{-7}~{}{\rm s}({\mu{\rm eV}\over m})^{1\over
3}$ is the time when axion field oscillations begin axdm . The thermal
relaxation rate of the axion fluid is of order CABEC ; Erken

$\Gamma(t)\sim 4\pi G\rho(t)m\ell(t)^{2}/\hbar~{}~{}\ .$ (3)

This result can be understood by noting that the axion fluid has density
flucuations $\delta\rho\sim\rho$ correlated over distances of order $\ell$;
see Appendix A. The density fluctuations source gravitational fields $\delta
g\sim 4\pi G\rho\ell$. Since the fluid momentum dispersion is $\delta
p\sim\hbar/\ell$, the fluctaating gravitational forces $m\delta g$ modify the
momentum distribution entirely in a time $\tau\sim\delta p/m\delta g$ The
thermal relaxation rate $\Gamma$ is the inverse of $\tau$.

$\Gamma(t)$ exceeds the Hubble rate $H(t)$ at some time well before equality
between matter and radiation. At that time all conditions for Bose-Einstein
condensation are satisfied and almost all axions go to the lowest energy
particle state available to them through the thermalizing interactions. The
axion fluid correlation length consequently grows to be of order the horizon
at the time. In linear order of perturbation theory, the axion condensate
evolves as ordinary CDM. However, in higher orders and in particular when it
acquires angular momentum by tidal torquing, the condensate can and will lower
its energy by acquiring vorticity whereas ordinary CDM remains vorticity free
Banik .

Consider a large overdensity of dark matter as it is about to collapse near
cosmic dawn. We ignore at first all density perturbations on scales smaller
than that of the overdenisty itself, effectivetly smotthing it. Its smoothed
density field has the general form

$\rho(\vec{r},t)=\rho(\vec{0},t)[1-\left({x_{1}\over
R_{1}(t)}\right)^{2}-\left({x_{2}\over R_{2}(t)}\right)^{2}\left({x_{3}\over
R_{3}(t)}\right)^{2}]$ (4)

where $(x_{1},x_{2},x_{3})$ are appropriately chosen Cartesian corrdinates
centered on the peak, and the $R_{i}~{}(i=1,2,3)$ give the peak’s extent in
the three spatial directions. Ref. Bard  derived the properties of such peaks
in Gaussian random fields, e.g. the number of peaks per unit volume of a given
size, and probability distributions for the $R_{i}$.

Fig. 1 shows the evolution of such an overdensity in a 2-dimensional cut
$(x,\dot{x})$ of its 6-dimensional phase-space. $x$ is the spatial coordinate
along an arbitrary direction through the overdensity and $\dot{x}$ the
corresponding velocity. Curve (b) shows the distribution of particles at time
$t_{\rm in}$ defined as the time when the central part of the overdensity is
at turnaround, and curve (d) the distribution of particles at the time $t_{\rm
coll}$ of collapse of the central overdensity. At $t_{\rm coll}$ the density
is very large near the center but no black hole forms in case of ordinary CDM
because the particles acquired angular momentum through tidal torquing.

The amount of angular momentum acquired by galaxies through tidal torquing is
commonly given by a dimensionless number $\lambda$ called the galactic spin
parameter Peeb . Spin parameters in the range $0.01\lesssim\lambda\lesssim
0.18$ are predicted Efst  and found to be consistent with the amount of
angular momentum baryons are observed to carry in galaxies Hern . We take this
range to be a guide to the amount of angular monetum that the overdesnity
under consideration acquired through tidal torquing. If the dark matter is
ordinary CDM, the angular momentum introduces an average distance of closest
approach to the center of the overdenisty of order $\lambda^{2}R$, far too
large for a black hole to form.

If axions are the dark matter, the axions at a distance $r$ from the center of
the overdensity thermalize at the rate given in Eq. (3) with $\ell\sim r$
since the gravitational fields due to the axion fluid outside the region of
radius $r$ do not help the thermalization within that region Erken . Since
$\rho(\vec{r},t)$ exceeds the average cosmological energy density
$\bar{\rho}(t)=1/6\pi Gt^{2}$,

$\Gamma(t)\gtrsim 2\cdot 10^{16}H(t)\left({m\over\mu{\rm
eV}}\right)\left({r\over 10^{22}~{}{\rm cm}}\right)\left({220~{}{\rm Myr}\over
t}\right)$ (5)

where $H(t)={2\over 3t}$ near cosmic dawn. The thermalization rate is
therefore very large compared to the dynamical evolution rate at time $t_{\rm
in}$. One may readlily verify that it remains much larger than the dynamical
evolution rate during the collapse. However, thermalizaton does not suffice to
justify the angular momentum transport necessary for black hole formation.
Since the axions can only change their momenta by an amount $\delta
p\sim{\hbar\over\ell}$ in a time $\tau\sim 1/\Gamma$, they can only change
their specific angular momentum $L$, i.e. their angular momentum per unit
mass, by an amount $\delta L\sim\hbar/m$ in that time. Hence

$\dot{L}_{\rm max}(r,t)\sim 4\pi G\rho(r,t)r^{2}$ (6)

is the maximum rate at which axions at radius $r$ can gain or lose specific
angular momentum.

The initial oversensity is not spherically symmetric, and hence its collapse
is not isotropic. Instead, its asphericity grows during the collpase anis  as
the overdensity tends to become a pancake or spindle. In our treatment below,
we ignore the fact that the infall is highy anistropic because anisotropy only
produces velocities in the angular (i.e. non radial) directions that are at
most of order $\sqrt{GM/r}$ where $M$ is the mass enclosed by the shell of
radius $r$. Since this is always much less than $c$ during the infall, the
angular velocities do not greatly affect whether a shell falls within the
black hole horizon of the mass $M$ it contains. Although we ignore the effect
of anisotropy on the infall motion we do not ignore its effect on the tidal
torquing experienced by the infalling overdensity. Tidal torquing is important
because it produces angular momentum and hence a minimum radius

$r_{\rm min}={L^{2}\over 2GM}+{\cal O}(L^{4})$ (7)

for each shell.

In case of isotropic radial infall the radius of the shell containing mass $M$
is given at time $t$ by the parametric equations

$\displaystyle r(M,t)$ $\displaystyle=$ $\displaystyle r_{0}(M)\sin^{2}\sigma$
$\displaystyle t$ $\displaystyle=$ $\displaystyle\sqrt{r_{0}(M)^{3}\over
2GM}[\sigma-{1\over 2}\sin(2\sigma)]$ (8)

where $r_{0}(M)$ is the shell’s turnaround radius. The shell’s turnaround time
is

$t_{0}(M)={\pi\over 2}\sqrt{r_{0}(M)^{3}\over 2GM}~{}~{}\ .$ (9)

The axion mass density at shell $M$ is

$\rho(M,t)={1\over 4\pi r(M,t)^{2}}{dM\over dr}(M,t)~{}~{}\ .$ (10)

By using Eqs. (8) to describe the collapse of an axion dark matter overdensity
four approximations are made in addition to ignoring the fact that the infall
is highly anisotropic. First, we ignore all forms or matter and energy
(baryons, photons,neutrinos, dark energy) other than the axionic dark matter.
Second, we use Newtonian gravity and Newton’s laws of motion to describe the
infal. Relativistic corrections are unimportant until a shell approaches the
black hole horizon and are unlikely to change the final black hole mass by
more than a factor two, or so, which is within the uncertainty that we
tolerate throughout. Third, although we keep track of the specific angular
momentum $L(M,t)$ of each shell, we ignore its effect on the radial motion of
the shell during its infall. This is justified because a shell can only fall
into a black hole if its angular momentum is extremely small. Fourth, we
ignore the corrections to the radial motion of the shells due to the wave
nature of the axion fluid. This is justified provided the axion Compton
wavelength is much less than the black hole size. For our smallest black
holes, which have mass of order $10^{6}~{}M_{\odot}$, the axion mass must be
more than of order $10^{-16}~{}{\rm eV}/c^{2}$. The condition is amply
satisfied by QCD axions.

Eqs. (8) describe the evolution of an overdensity entirely in terms of the
function $r_{0}(M)$ that gives the turnaround radius of each shell. We define
our initial time $t_{\rm in}$ to be the turnaround time of the innermost
shells:

$t_{\rm in}=\lim_{M\rightarrow 0}t_{0}(M)\ .$

At that time, the density near the center has the form

$\rho(r,t_{\rm in})=\rho(0,t_{\rm in})[1-\left({r\over R}\right)^{2}+...]$
(11)

where the dots represent additional terms in an expansion in powers of $r/R$.
We require $r_{0}(M)$ to be such that Eq. (11) is reproduced for $r<<R$, and
such that $\rho(r,t_{\rm in})$ approaches the contemporary average
cosmological energy density $\bar{\rho}(t_{\rm in})$ for $r>>R$. The
requirements are met by the choice:

$r_{0}(M)=R[\left({M\over M_{f}}\right)^{1\over 3}+{1\over 5}{M\over
M_{f}}+{1\over 2}({M\over M_{f}})^{2}]$ (12)

with

$M_{f}={4\pi\over 3}\rho(0,t_{\rm in})R^{3}~{}~{}\ .$ (13)

The resulting initial density profile $\rho(r,t_{\rm in})$ is shown in Fig. 2.

The central overdensity collapses at time

$t_{\rm coll}=2t_{\rm in}=\pi\sqrt{3\over 8\pi G\rho(0,t_{\rm in)}}~{}~{}\ ,$
(14)

with the outer shells collapsing later. We define $M_{\star}(t)$ such that all
axions within shell $M_{\star}(t)$ (but outside any black hole that may have
formed) thermalize sufficiently fast that they rotate rigidly, in the three
dimensional sense. Let $\omega(t)$ be their angular rotation frequency at that
time. For $M<M_{\star}(t)$ the specific angular momentum of the axions at the
equator of shell $M$ is therefore

$L(M,t)=\omega(t)r(M,t)^{2}~{}~{}\ .$ (15)

Its rate of change following the motion is

$\dot{L}(M,t)={d\omega\over dt}(t)r(M,t)^{2}+2\omega(t)r(M,t)v_{r}(M,t)$ (16)

where $v_{r}(M,t)=\dot{r}(M,t)$ is the radial velocity of the shell. The
maximum rate at which the shell can shed its specific angular monentum
$L(M,t)$ while collapsing is given by Eq. (6) with $r=r(M,t)$.

If there were no relaxation, $\dot{L}(M,t)=0$ and the angular frequency of
each shell would increase as $r(M,t)^{-2}$. Instead, relaxation allows the
shells within $M_{\star}(t)$ to collapse without harrdly incerasing their
angular rotation frequency $\omega(t)$ because their angular momentum is
transported outward. The condition for a shell to collapse without hardly
increasing its angular rotation frequency is that the RHS of Eq. (6) is larger
than the second term on the RHS of Eq. (16), or equivalently that

$\omega(t)\lesssim 2\pi G\rho(M,t)r(M,t){1\over v_{r}(M,t)}\equiv\omega_{\rm
max}(M,t)~{}~{}\ .$ (17)

$M_{\star}(t)$ is therefore the largest shell such that

$\omega(t)<\omega_{\rm max}(M,t)$ (18)

for all $M<M_{\star}(t)$. Because $\omega_{\rm max}(M,t)$ is, at all times, a
decreasing function of $M$ near $M=0$, $M_{\star}(t)$ is the smallest solution
of

$\omega_{\rm max}(M_{\star}(t),t)=\omega(t)~{}~{}\ .$ (19)

All axions within the volume enclosed by $M_{\star}(t)$ exchange angular
momentum suffiiently fast that they can, and therefore do, rotate with the
common angular frequency $\omega(t)$. On the other hand all axions outside
shell $M_{\star}(t)$ conserve their angular momentum. Of course, the
transition at $M_{\star}(t)$ is not sudden as we take it to be but smoothing
it out is not expected to change the final outcomes significantly.

$\omega(t)$ increases with time for two distinct reasons:

${d\omega\over dt}={d\omega\over dt}\Big{|}_{L}+{d\omega\over
dt}\Big{|}_{T}~{}~{}\ .$ (20)

The first term is due to the conservation of angular momentum within the
volume of axions that rotate rigidly, and the second term is due to tidal
torquing. The first term is

${d\omega\over dt}\Big{|}_{L}=-\omega(t){\dot{I}(t)\over I(t)}$ (21)

with $I(t)$ the moment of inertia of all axions between shells $M_{\rm bh}(t)$
and $M_{\star}(t)$, where $M_{\rm bh}(t)$ is the black hole mass at time $t$

$I(t)={2\over 3}\int_{M_{\rm bh}(t)}^{M_{\star}(t)}dM~{}r(M,t)^{2}$ (22)

and

$\dot{I}(t)={4\over 3}\int_{M_{\rm
bh}(t)}^{M_{\star}(t)}dM~{}r(M,t)v_{r}(M,t)$ (23)

its time derivative following the motion. The second term on the RHS of Eq.
(20) is estimated in Appendix B. Let us write the initial value of $\omega(t)$
as

$\omega(t_{\rm in})\equiv\omega_{\rm in}={j\over t_{\rm in}}~{}~{}\ .$ (24)

In case of rigid rotation, the relationship between $j$ and the spin parameter
is

$\lambda={4\over 5\pi}\sqrt{6\over 5}j=0.279~{}j+{\cal O}(j^{2})~{}~{}\ .$
(25)

We therefore expect $j$ to be in the approximate range 0.03 to 0.8. We find in
Appendix B that

${d\omega\over dt}\Big{|}_{T}\simeq 2.2~{}j~{}{1\over t_{\rm
in}^{2}}\left({t_{\rm in}\over t}\right)^{4\over 3}$ (26)

during the interval $t_{\rm in}<t<t_{\rm coll}$.

Eqs. (17),(19), (20), (21) and (26) were solved numerically. Initially,
$M_{\star}$ is of order $M_{f}$ or larger. As $\omega(t)$ increases, it
becomes more and more difficult to maintain rigid rotation, $M_{\star}(t)$
starts to decrease and then accelerates towards $M=0$. Provided $\omega$ is
sufficently small at time $t_{\rm coll}$, a black hole forms and grows. Soon
thereafter, at a time $t_{f}$, $r_{(}M_{\star}(t),t)$ reaches zero and
relaxation stops. After $t_{f}$, $L(M,t)$ of each shell is conserved. The
black hole mass is the largest $M$ for which $r_{\rm min}(M,t_{f})<2GM/c^{2}$
or equivalently $L(M,t_{f})<2GM/c$.

For given $j$ the black hole mass $M_{\rm bh}$ is found to be very nearly
proportional to $M_{f}$. It has only a small dependence on $t_{\rm in}$ for
given $M_{f}$ and $j$. Fig. 3 shows $M_{\rm bh}$ as a function of $j$ for a)
$M_{f}=3\cdot 10^{10}M_{\odot}$, b) $M_{f}=10^{9}M_{\odot}$, and c)
$M_{f}=3\cdot 10^{7}M_{\odot}$. For the sake of definiteness, collapse was
assumed to occur at redshift $z=10$, so that $t_{\rm coll}\simeq$ 0.5 Gyr, and
hence $t_{\rm in}\simeq$ 0.25 Gyr, implying $\rho(0,t_{\rm in})\simeq 7.0\cdot
10^{-26}$ gr/cc. The three cases were chosen to correspond to overdensities
that will evolve later into a) very large galaxies, b) galaxies similar to the
Milky Way c) small galaxies. Since the ratio $M_{\rm bh}/M_{f}$, a function of
$j$, has little dependence on the collpase time, it is straightforward to
generate results for all plausible cases.

The black hole masses that result from the inital collpase of axion
overdensities near cosmic dawn are in surprisingly good agreement with
observations. First, the range of black hole masses formed, from approximately
$10^{5}M_{\odot}$ to a few times $10^{10}M_{\odot}$ is the mass range of
observed supermassive black holes. A cutoff is predicted at low masses.
Specifially the theory predicts that a black hole mass less than
$10^{6}~{}M_{\odot}$ is unlikely on galaxy scales. Also, masses larger than
$10^{11}~{}M_{\odot}$ are unlikey. They occur only in the largest
overdensities and only if $\omega_{\rm in}$ is unexpectedly small. Second, the
predicted supermassive black holes form near cosmic dawn. The puzzle of why
supermassive black holes appear at high redshifts is removed. Although they
may merge and accrete later, mergers and accretion are not necessary to
explain their size. Third, there is a strong correlation between black hole
mass and galaxy size. On the other hand, for a given galaxy size the black
hole mass ranges widely, by a factor hundred or so. Both the correlation with
galaxy size and the intrinsic variability are in qualitative agreement with
observation. Fourth, the black holes form for approximately the range of
galactic spin parameters expected from tidal torquing.Since there are galaxies
with black holes close to the cutoff in $j$, e.g. the Milky Way whose central
black hole has mass $4.3\cdot 10^{6}M_{\odot}$, the theory suggests that there
are galaxies without supermassive black hole because their $j$ happens to be
larger than the cutoff. This too appears consistent with observation. For
example, the nearby small galaxy M33 appears not to have a supermassive black
hole M33 ..

Our description of supermassive black hole formation does not make any ad-hoc
assumption. No new particle is postulated, other than the standard QCD axion
or an axion-like particle with similar properties. The invoked processes of
thermalization and angular momentum trnasport in the cold axion dark matter
fluid are the same as were described previously in refs. CABEC ; Erken ; Banik
. The theory is predicitve and can be tested further. It should be possible to
predict the distribution of supermassive black hole masses from the
distribution of overdensities of a given size at comsic dawn and the
distribution of galactic spin parameters. The amplitude and spectrum of
gravitational waves produced can be calculated, to be compared with
observation. By including baryons, it may be possible to determine whether the
theory is consistent with the observed relation Msigma  between the stellar
velocity disperion in a galactic bulge and the mass of the supermassive black
hole mass at its center.

###### Acknowledgements.

We are grateful to Laura Blecha for critical feedback and to Jeff Andrews,
Abhishek Chattaraj, Jeff Dror, Antonios Kyriazis, Wei Xue and Fengwei Yang for
useful discussions. This work was supported in part by the U.S. Department of
Energy under grant DE-SC0022148 at the University of Florida.

Appendix A: Axion fluid fluctuations

It has been shown in a variety of contexts that degenerate Bosonic systems
have generically large flucutations in intensity or density Pathria . In this
appendix we show this for the cold dark matter axion fluid.

In the non-relativisitc limit, the scalar quantum field $\phi(\vec{x},t)$
describing axions in a volume $V$ may be written

$\phi(\vec{x},t)={1\over\sqrt{2m}}[\psi(\vec{x},t)e^{-imt}+h.c.]$ (27)

and expanded

$\psi(\vec{x},t)=\sum_{\vec{\alpha}}u^{\vec{\alpha}}(\vec{x},t)a_{\vec{\alpha}}(t)$
(28)

where the $u^{\vec{\alpha}}(\vec{x},t)$ are any set of orthonormal and
complete (ONC)wavefunctions in that volume:

$\int_{V}d^{3}x~{}u^{\vec{\alpha}}(\vec{x},t)^{*}u^{\vec{\beta}}(\vec{x},t)=\delta_{\vec{\alpha}}^{\vec{\beta}}~{}~{}~{},~{}~{}~{}~{}\sum_{\vec{\alpha}}u^{\vec{\alpha}}(\vec{x},t)u^{\vec{\alpha}}(\vec{y},t)^{*}=\delta(\vec{x}-\vec{y})~{}~{}\
.$ (29)

The $a_{\vec{\alpha}}(t)$ and their Hermitean conjugates
$a_{\vec{\alpha}}(t)^{\dagger}$ satisfy canonical equal time commutation
relations. The most general axion system state is given by a linear
combination

$|c_{\\{\cal N\\}}\rangle=\sum_{\\{\cal N\\}}c_{\\{\cal N\\}}|\\{\cal
N\\}\rangle$ (30)

of all possible particle state occupation number eigenstates

$|\\{{\cal N}\\}\rangle=\prod_{\vec{\alpha}}{1\over\sqrt{{\cal
N}_{\vec{\alpha}}!}}(a_{\vec{\alpha}}^{\dagger})^{{\cal
N}_{\vec{\alpha}}}|0\rangle~{}~{}\ .$ (31)

Here $|0\rangle$ is the empty state and $\\{{\cal N}\\}=\\{{\cal
N}_{\vec{\alpha}}:\forall\vec{\alpha}\\}$ is the set of integers giving the
occupation number of each particle state. $N=\sum_{\vec{\alpha}}{\cal
N}_{\vec{\alpha}}$ is the total number of axions in volume $V$.

For the purpose of describing the cold dark matter axion fluid with aveage
number density ${1\over m}\rho(\vec{x},t)$ and average velocity field
$\vec{v}(\vec{x},t)$ we choose as ONC set of wavefunctions the wavefucntion of
Eq. (1) and spatial modulations thereof with waveector $\vec{k}$:

$u^{\vec{k}}(\vec{x},t)=e^{i\vec{k}\cdot\vec{\chi}(\vec{x},t)}\Psi(\vec{x},t)\sim
e^{i\vec{k}\cdot\vec{x}}\Psi(\vec{x},t)~{}~{}\ .$ (32)

If the axion fluid is homogeneous and at rest, this ONC set of wavefunctions
would simply be

$u^{\vec{k}}(\vec{x})={1\over\sqrt{V}}e^{i\vec{k}\cdot\vec{x}}~{}~{}\ .$ (33)

When the axion fluid is inhomgeneous and/or moving, there are many ways to
construct suitable $u^{\vec{k}}(\vec{x},t)$. A particular method is described
in ref. Chak . Another is to start with all the wavefuntions
$e^{i\vec{k}\cdot\vec{x}}\Psi(\vec{x},t)$ and orthonormalize them in
succession. In any such basis, the state of the axion fluid is one where most
particle states have low occupation numbers but those with wavevector
magnitude $k\lesssim m\delta v/\hbar$ are hugely occupied. The occupation
numbers of those states that are occupied are of order $10^{61}$ or larger
CABEC ; Erken .

The number denisty of axions is the operator

$n(\vec{x},t)=\psi(\vec{x},t)^{\dagger}\psi(\vec{x},t)~{}~{}\ .$ (34)

In occupation number eignestates, it has quantum mechanical average

$\langle\\{{\cal N}\\}|n(\vec{x},t)|\\{{\cal N}\\}\rangle=\sum_{\vec{k}}{\cal
N}_{\vec{k}}|u^{\vec{k}}(\vec{x},t)|^{2}={1\over m}\rho(\vec{x},t)~{}~{}\ .$
(35)

Let us define the operator

$\delta n(\vec{x},t)=n(\vec{x},t)-{1\over m}\rho(\vec{x},t)~{}~{}\ .$ (36)

One readily finds

$\langle\\{{\cal N}\\}|\delta n(\vec{x},t)\delta n(\vec{y},t)|\\{{\cal
N}\\}\rangle=|D(\vec{x},\vec{y};t)|^{2}~{}[1+{\cal O}\left({1\over{\cal
N}}\right)]$ (37)

where

$D(\vec{x},\vec{y};t)=\sum_{\vec{k}}{\cal
N}_{\vec{k}}~{}u^{\vec{k}}(\vec{x},t)^{*}u^{\vec{k}}(\vec{y},t)~{}~{}\ .$ (38)

We have therefre

$\langle\\{{\cal N}\\}|(\delta n(\vec{x},t))^{2}\\{{\cal
N}\\}\rangle=\left({1\over m}\rho(\vec{x},t)\right)^{2}[1+{\cal
O}\left({1\over{\cal N}}\right)]~{}~{}.$ (39)

Thus, in eignestates of the occupation numbers, the root mean square deviation
from the average density at every space-time point equals the average density
there. Moreover these deviations are correlated over distances of order
$\ell={1\over\delta k}={\hbar\over m\delta v}$ since $D(\vec{x},\vec{y};t)$
cannot vary much over distances shorter than $1/\delta k$.

In the general system states of Eq. (30) it is not possible to make such
strong statements as Eq (39) because they include system states in which all
axions are in a single state whose wavefunction is a linear combination of
several $u^{\vec{k}}(\vec{x},t)$. In these very special system statss, the
quantum mechanical uncertainty in measuring $n(\vec{x},t)$ vanishes. Consider
nonetheless the quamtum mechanical average in a general state of any operator
$\Omega$

$\langle c_{\\{{\cal N}\\}}|\Omega|c_{\\{{\cal N}\\}}\rangle=\sum_{\\{{\cal
N}\\}}\sum_{\\{{\cal N}^{\prime}\\}}c_{\\{\cal N}\\}^{*}c_{\\{{\cal
N}^{\prime}\\}}\langle\\{{\cal N}\\}|\Omega|\\{{\cal
N}^{\prime}\\}\rangle~{}~{}\ .$ (40)

Even in case the initial values of the complex coefficients $c_{\\{{\cal
N}\\}}$ are very special, they acquire random phases after some time,
especially as a result of thermal relaxation of the fluid. We have

$\cdot\langle c_{\\{\cal N}\\}|\Omega|c_{\\{{\cal
N}\\}}\rangle\cdot=\sum_{\\{{\cal N}\\}}|c_{\\{\cal N}\\}|^{2}\langle\\{{\cal
N}\\}|\Omega|\\{{\cal N}\\}\rangle$ (41)

where $~{}~{}\cdot\langle...\rangle\cdot~{}~{}$ means quatnum mechanical
average followed by an average over the phases of the coefficient $c_{\\{{\cal
N}\\}}$. Using Eq. (39) we have

$\cdot\langle c_{\\{{\cal N}\\}}|(\delta n(\vec{x},t))^{2}|c_{\\{{\cal
N}\\}}\rangle\cdot=({1\over m}\rho(\vec{x},t))^{2}~{}~{}\ .$ (42)

Appendix B: Tidal torquing of the ovdedensity

Consider the tidal torque on all the particles between shels $M_{\rm bh}(t)$
and $M_{\star}(t)$

$\vec{\tau}(t)=\int d\Omega\int_{r(M_{\rm
bh}(t),t)}^{r(M_{\star}(t),t)}~{}r^{2}dr~{}\rho(\vec{r},t)~{}\vec{r}\times(-\vec{\nabla}\Phi(\vec{r},t)+\vec{\nabla}\Phi(\vec{0},t))$
(43)

where $\Phi(\vec{r},t)$ is the gravitational potential due to density
perturbations outside the overdensity of interest. The region surrounding the
overdensity near cosmic dawn is an expanding Einstein-de Sitter space-time
where $\vec{r}=\vec{x}\left({t\over t_{\rm in}}\right)^{2\over 3}$, $\vec{x}$
are comoving coordinates, and $\Phi(\vec{r},t)=\Phi(\vec{x})$ SW . We expand
the gravitational potential in Taylor series about the center of the
overdensity

$\Phi(\vec{x})=\Phi(\vec{0})-\vec{g}\cdot\vec{x}+{1\over
2}\vec{x}^{T}\mathbb{T}\vec{x}+...$ (44)

where $\vec{g}$ is a constant vector and $\mathbb{T}$ a constant symmetric
matrix. Eq. (43) becomes

$\vec{\tau}(t)\simeq-\left({t\over t_{\rm in}}\right)^{4\over 3}\int
d\Omega\int_{r(M_{\rm
bh}(t),t)}^{r(M_{\star}(t),t)}~{}r^{4}dr~{}\rho(\vec{r},t)~{}\hat{n}(\theta,\phi)\times\mathbb{T}~{}\hat{n}(\theta,\phi)~{}~{}\
,$ (45)

where $\hat{n}(\theta,\phi)={1\over r}\vec{r}$ is the unit vector in the
direction of spherical coordinates $(\theta,\phi)$. We expect the asphericity
of the smoothed overdensity, as in Eq. (4), to doninate the integral in Eq.
(45). Let us rewrite Eq. (4) as

$\rho(\vec{r},t_{\rm in})=\rho(\vec{0},t_{\rm in})[1-\left({r\over
R}\right)^{2}A(\theta,\phi)]$ (46)

where $A(\theta,\phi)$ is a linear combination of the spherical harmonics up
to second order.

For $0<t<t_{\rm in}$, we estimate the integral in Eq. (45) using linear
perturbation theory in a homogeneous expanding universe of denisty
$\rho(\vec{0},t)$. In such a universe

$r(M,t)=b(t)r(M,t_{\rm in})$ (47)

where the scale factor $b(t)$ is given by

$\displaystyle b(t)$ $\displaystyle=$ $\displaystyle\sin^{2}\sigma$
$\displaystyle t$ $\displaystyle=$ $\displaystyle\sqrt{3\over 8\pi
G\rho(\vec{0},t_{\rm in})}[\sigma-{1\over 2}\sin(2\sigma)]~{}~{}\ .$ (48)

We have therefore

${\delta\rho(r(M,t)\hat{n},t)\over\rho(\vec{0},t)}=\delta_{+}(t){\delta\rho(r(M,t_{\rm
in})\hat{n},t_{\rm in})\over\rho(\vec{0},t_{\rm in})}$ (49)

where $\delta_{+}(t)$ is the appropriate growth fuction SW  normallized so
that $\delta_{+}(t_{\rm in})=1$:

$\delta_{+}(\sigma)={1\over
2}\left({-3\sigma\cos\sigma\over\sin^{3}\sigma}+{3\over\sin^{2}\sigma}-1\right)~{}~{}\
.$ (50)

Since

$\delta\rho(r(M,t_{\rm in})\hat{n},t_{\rm in})=-\rho(\vec{0},t_{\rm
in})\left({r(M,t_{\rm in})\over R}\right)^{2}A(\theta,\phi)$ (51)

we have

$\delta\rho(r(M,t)\hat{n},t)=-\delta_{+}(t)\rho(\vec{0},t)\left({r(M,t)\over
R}\right)^{2}A(\theta,\phi)~{}~{}\ .$ (52)

We substitute this for $\rho(\vec{r},t)$ in Eq. (45) and define

${3\over 8\pi}\int d\Omega
A(\theta,\phi)\hat{n}(\theta,\phi)\times\mathbb{T}~{}\hat{n}(\theta,\phi)\equiv\Sigma\hat{z}~{}~{}\
.$ (53)

Upon integrating the radial coordinate $r$ from 0 to $R(t)=b(t)R$, Eq.‘(45)
becomes

$\vec{\tau}(t)\simeq{8\pi\over 3}\Sigma\hat{z}\left({t_{\rm in}\over
t}\right)^{4\over 3}{1\over 7}R^{5}b(t)^{7}\rho(\vec{0},t)\delta_{+}(t)~{}~{}\
.$ (54)

Since

$\vec{\tau}(t)=I(t)\hat{z}{d\omega\over dt}\Bigg{|}_{T}$ (55)

with

$I(t)={2\over 3}\int_{0}^{R(t)}4\pi r^{4}dr~{}\rho(\vec{0},t)=4\pi{2\over
15}R^{5}b(t)^{5}\rho(\vec{0},t)$ (56)

we have

${d\omega\over dt}\Bigg{|}_{T}\simeq{5\over 7}\Sigma\left({t_{\rm in}\over
t}\right)^{4\over 3}b(t)^{2}\delta_{+}(t)$ (57)

during the time interval $0<t<t_{\rm in}$. Hence

$\omega(t_{\rm in})\simeq{5\over 7}\Sigma\int_{0}^{t_{\rm in}}dt\left({t_{\rm
in}\over t}\right)^{4\over 3}b(t)^{2}\delta_{+}(t)\simeq
0.448~{}\Sigma~{}t_{\rm in}~{}~{}\ .$ (58)

During the time interval $t_{\rm in}<t<t_{\rm coll}$, tidal torquing continues
but perturbation theory breaks down. The mass distribution becomes dominated
by its anisotropic component. We estimat the integral in Eq. (45) by setting

$\rho(r(M,t)\hat{n}(\theta,\phi),t)=\rho(M,t)A(\theta,\phi)$ (59)

which yields

$\vec{\tau}(t)\simeq\left({t_{\rm in}\over t}\right)^{4\over
3}\Sigma~{}I(t)~{}\hat{z}$ (60)

and hence

${d\omega\over dt}\Bigg{|}_{T}\simeq\left({t_{\rm in}\over t}\right)^{4\over
3}\Sigma~{}~{}\ .$ (61)

Combining this with Eq. (58) yields

${d\omega\over dt}\Bigg{|}_{T}\simeq 2.2~{}{\omega(t_{\rm in})\over t_{\rm
in}}\left({t_{\rm in}\over t}\right)^{4\over 3}$ (62)

during the time interval $t_{\rm in}<t<t_{\rm coll}$.

## References

  * (1) J. Kormendy and L.C. Ho, Ann. Rev. of Astron. and Astroph., 51 (2013) 511.
  * (2) A. Kazunori et al. (the EHT Collaboration), Ap. J. Lett. 875 (2019) L5, and Ap. J. Lett. 930 (2022) L12.
  * (3) J. Frank, A. King and D.J. Raine, Accretion Power in Astrophysics, Third Edition, Cambridge University Press, 2002.
  * (4) G. Agazie et al. (the NANOGrav. Collaboration), Ap. J. Lett. 951 (2023) L8; D.J. Reardon et al., Ap. J. Lett. 951 (2023) L6; H. Xu et al., Res. Astron. Astroph. 23 (2023) 075024; J. Antoniadis et al. (the EPTA Collaboration), arXiv: 2306.16227.
  * (5) K. Inayoshi, E. Visbal and Z. Haiman, Ann. Rev. of Astron. and Astroph., 58 (2020) 27.
  * (6) S. Balberg and S.L. Shapiro, Phys. Rev. Lett. 88 (2002) 101301; J. Pollack, D.N. Spergel and P. Steinhardt, Ap. J. 804 (2015) 2, 131; W.-X. Feng, H.-B. Yu and Y.-M. Zhong, Ap. J. Lett. 914 (2021) 2, L26.
  * (7) T. Rindler-Daller, K. Freese, M.H. Montgomery, D. Winget and B. Paxton, Ap. J. 799 (2015) 210.
  * (8) R. Larson et al., Ap. J. Lett. 953 (2023) L29; A. Bogdan et al., Nature Atron. 8 (2024) 126; R. Maiolino et al., arXiv:2308.01230; L.J. Furtak et al., arXiv: 2308.05735; R. Miaolino et al., Nature 627 (2024) 59;
  * (9) D. Hooper, A. Ireland, G. Krujaic and A Stebbins, JCAP 04 (2024) 021.
  * (10) P. Sikivie and Q. Yang, Phys. Rev. Lett. 103 (2009) 111301.
  * (11) O. Erken, P. Sikivie, H. Tam and Q. Yang, Phys. Rev. D 85 (2012) 063520.
  * (12) M.Y. Khlopov, B.A. Malomed and Y.B. Zel’dovich, MNRAS 215 (1985) 575.
  * (13) J. Preskill, F. Wilczek and M. Wise, Phys. Lett. B120 (1983) 127; L. Abbott and P. Sikivie, Phys. Lett. B120 (1983) 133; M. Dine and W. Fischler, Phys. Lett. B120 (1983) 137.
  * (14) N. Banik and P. Sikivie, Phys. Rev. D 88 (2013) 123517.
  * (15) J.M. Bardeen, J.R. Bond, N. Kaiser and A.S. Szalay, Ap. J. 304 (1986) 15.
  * (16) P.J.E. Peebles, Ap. J. 155 (1969) 393.
  * (17) G. Efstatathiou and B.J.T. Jones, MNRAS 186 (1979) 133; J. Barnes and G. Efstathiou, Ap. J. 319 (1987) 575.
  * (18) X. Hernanadez, C. Park, B. Cervantes-Sodi and Y.-Y. Choi, MNRAS 375 (2007) 163.
  * (19) C.C. Lin, L. Meistel and F.H. Shu, Ap. J. 142 (1965) 1431; Y.B. Zel’dovich, Astron. and Astroph. 5 (1970) 84; J. Binney, Ap. J. 215 (1977) 492.
  * (20) K. Gebhardt et al., Astron. J. 122 (2001) 2469.
  * (21) F. Ferrarese and D. Merritt, Ap. J. 538 (2000) L9; K. Gebhardt et al., Ap. J. (2000) L13.
  * (22) R.K. Pathria and P.D. Beale, Statistical Mechanics, 3rd edition, Elsevier 2011, and references therein.
  * (23) S. Chakrabarty et al., Phys. Rev. D 97 (2018) 043531.
  * (24) S. Weinberg, Gravitation and Cosmology, J. Wiley and Sons, 1972.

Figure 1: Phase space distribution of cold collisionless particles during the
collpase of a large smooth overdensity near cosmic dawn, at four different
times: a) just after the Big Bang, b) when the central part of the overdeisty
is at turnaround, c) just before, and d) at the time $t_{\rm coll}$ of
collapse of the central overdensity. $x$ is the spatial coordinate along an
arbitrary direction through the overdensity. An actual overdensity has small
scale structure which has been smoothed out in the figure. Figure 2: Density
profile at time $t_{\rm in}$, when the central part of the overdensity is at
turnaround, in units where $M_{f}=1$ and $R=1$. In these units the central
density $\rho(\vec{0},t_{\rm in})=3/4\pi$ and the contemporary average
cosmological enegy density $\bar{\rho}(t_{\rm in})=4/3\pi^{3}$. The average
cosmological enegy density is indicated by the horizontal dashed line. The
vertical solid line indicates the radius that contains mass $M_{f}$. Figure 3:
Black hole mass as a function of $j\equiv\omega_{\rm in}t_{\rm in}$ for three
values of $M_{f}$. For given $M_{f}$ there is only a very slight dependence of
the black hole mass on $t_{\rm in}$. The values shown were computed for
$z_{\rm coll}$ = 10 and hence $t_{\rm in}$ = 240 Myr.


# Uber Stable: Formulating the Rideshare System as a Stable Matching Problem

Rhea Acharya1<EMAIL_ADDRESS>Jessica Chen1
<EMAIL_ADDRESS>Helen Xiao1<EMAIL_ADDRESS>
###### Abstract

Peer-to-peer ride-sharing platforms like Uber, Lyft, and DiDi have
revolutionized the transportation industry and labor market. At its essence,
these systems tackle the bipartite matching problem between two populations:
riders and drivers. This research paper comprises two main components: an
initial literature review of existing ride-sharing platforms and efforts to
enhance driver satisfaction, and the development of a novel algorithm
implemented through simulation testing to allow us to make our own
observations. The core algorithm utilized is the Gale-Shapley deferred
acceptance algorithm, applied to a static matching problem over multiple time
periods. In this simulation, we construct a preference-aware task assignment
model, considering both overall revenue maximization and individual preference
satisfaction. Specifically, the algorithm design incorporates factors such as
passenger willingness-to-pay, driver preferences, and location attractiveness,
with an overarching goal of achieving equitable income distribution for
drivers while maintaining overall system efficiency.

Through simulation, the paper compares the performance of the proposed
algorithm with random matching and closest neighbor algorithms, looking at
metrics such as total revenue, revenue per ride, and standard deviation to
identify trends and impacts of shifting priorities. Additionally, the DA
algorithm is compared to the Boston algorithm, and the paper explores the
effect of prioritizing proximity to passengers versus distance from city
center. Ultimately, the research underscores the importance of continued
exploration in areas such as dynamic pricing models and additional modeling
for unconventional driving times to further enhance the findings on the
effectiveness and fairness of ride-sharing platforms.

Our simulation code can be found here.

##  I Introduction

The subway breaks down on our way to work, it’s 2:30 a.m. and we need to get
home after a night of clubbing, we’re in a new city and need to get from the
airport to a museum across town — what’s our first instinct? ”Call an Uber.”
In this digital age, we rely on ride-sharing apps as an integral part of daily
life. Passengers turn to driving apps, trusting that they will be able to get
from their current destination to the next in a fast and safe way, without
breaking their wallets. Ride-sharing apps like Uber, Lyft, and Didi, eager to
gain more revenue, prioritize meeting these needs to provide a positive
experience to passengers, hoping they will use the apps more and spend more
money.

This means that apps have historically maximized rider satisfaction, ensuring
that riders are likely to continue to use their services, as without them, the
apps would not be able to survive. However, research has shown that in
matching problems like this one, when the needs of one party is prioritized,
those of the other can be sometimes hurt. Here, by prioritizing riders, the
matching mechanisms might be creating pessimal outcomes for the drivers. If
the negative effects are too large, this could have significant consequences
on driver satisfaction and equity.

Potential concerns related to this have already been seen by the public eye.
In recent years, these ridesharing apps have come under fire for the driver
experience that they provided. For example, a CNBC article found that only 4%
of Uber drivers are still driving a year later. Additionally, Uber and Lyft
have bought been sued for withholding money from drivers, and drivers are
frequently commenting on the unsustainability of the role in outlets like The
New York Times. Why is this the case? Is it inherent to the nature of the
role, or is it at all impacted by the Uber matching algorithm itself?

Motivated by these questions, we want to dive deeper into the issue of driver
equity, exploring through the driver’s perspective how the design of matching
algorithms can affect the efficiency and equity. Our deliverable will have two
main components: 1) a literature review of existing ridesharing platforms and
efforts at increasing driver satisfaction; 2) a novel algorithm implementation
with simulation testing. We will spend significant time defining an
application of the Gale-Shapley deferred acceptance algorithm to this
rideshare system matching problem, and we will compare the results of this to
other algorithms, like random matching and closest neighbours.

##  II Previous Work

To begin, we want to better understand how the Uber system works and the
general ride-sharing matching algorithm. Unfortunately, ridesharing services
such as Uber and Lyft are not required to share data like taxi services, so
exact algorithms can not be found through literature, yet we will review work
discussing the current research of these algorithms and what they aim to
achieve (in most detail). In order to ground our knowledge of current work in
the field to inform the design of our simulation, we conducted a preliminary
literature review. Additionally, we will focus on the problems that these
algorithms pose to the driver-side of the system: in particular, we will do
brief literature review on the trade-offs between maximizing objectives in
ride-share matching algorithms and driver revenue distribution. Next, we will
look at how Uber and other ride-sharing companies generate matching algorithms
to best generate aggregate revenue for the company, and how this may affect
consumers of the app. Ultimately, we will use these insights to inform the
design of our simulation and will aim to work towards a solution to various
problems addressed, through our preference ordering and nuanced matching
algorithm.

###  II-A Two-Sided Fairness

Most ride-sharing platforms prioritize maximizing passenger satisfaction, as
this is directly correlated with maximizing revenue. There is a trade-off
between income distribution and maximizing the overall objective, which in
this case is to maximize revenue of app and minimize the total distance
covered by drivers. While existing platforms primarily prioritize customer
satisfaction, Sühr et al [1] emphasize the importance of considering fairness
for drivers to ensure their well-being and long-term stability. The authors
propose a framework for fair income distribution by amortizing equality over
time, allowing drivers to receive benefits proportional to their active time
on the platform. Specifically, the paper suggests algorithms that minimize
income inequality while fairly distributing increases in customer waiting
times. As we construct our own algorithm and analyze our results, we try to
strike a balance between these two conflicting ideals, as we applied its
concepts to a static matching problem rather than dynamic. [1]

###  II-B Preference-Aware Task Assignment

Previous work in this field includes OSM-KIID [2], a model whose long-form
name is Online Stable Matching under Known Identical Independent
Distributions. In many matching problems, especially those constructed by
commercial platforms like Uber, the aim is simply to maximize the profit of
the system as a whole. This means that frequently, the individual preferences
of parties on both sides are neglected, which in the short term might seem
inconsequential but in the long term could have more drastic impacts, as
participants may not longer feel that it is worthwhile to themselves to engage
in the matching process. For example, in the rideshare system, if drivers feel
as though their time and effort are not being compensated appropriately or
they always get the least desired rides (like those that start or end in
remote areas), then they may quit their job which is not good for the system
as a whole. The authors of OSM-KIID seek to create a new model that takes into
account both the overall revenue minimization objective as well as the
preference satisfaction one, by maximizing revenue and minimizing blocking
pairs. We incorporate this focus on preferences within our deferred acceptance
algorithm, by taking into account factors like prior income distribution and
desirability of starting and ending locations within the preference ordering
generation. Additionally, we account for individual diversity of preferences
by having unique cost coefficients among other nuances. [2]

##  III Algorithm Design

###  III-A Overview

At the most fundamental level, we want to model the rideshare system as a two-
sided matching problem. Then, we will explore how we can design our matching
algorithm to this bipartite problem with driver equity in mind. Specifically,
we want to achieve equitable income distribution among drivers while
maintaining comparable total welfare, as measured by total revenue and
aggregate wait time. Although, in reality, the algorithms used by Uber, Lyft,
and other ride-sharing apps are dynamic, given the shortened timeframe of this
project, we will limit our algorithm to a static algorithm run over multiple
time periods.

###  III-B Approach

  1. 1.

Reformulate this problem into a deferred-acceptance two-sided matching
problem.

  2. 2.

Driver-proposing: Create cost functions for drivers and (fixed) payment
functions. Drivers will propose a ride if utility is positive for a given
request.

  3. 3.

Both drivers and riders have preference orders, and of note is the fact that
we will incorporate driver income into preference ordering to look into how
driver equity is prioritized

###  III-C Some Relevant Assumptions

  * •

We do not have equilibrium pricing. Instead, price is determined by riders’
fixed willingness to pay given location with their WTP is fixed and known in a
given round.

  * •

We constrain goal to driver-side only, so we assume riders simply want to get
picked up.

  * •

We have our ride-sharing platform as a benevolently matching central
allocation mechanism, with no revenue generation abilities.

  * •

We model this as a static problem, where the algorithm only aims to solve
problem for all drivers and riders currently on the map. Drivers are constant
through all the rounds, in order to allow us to most helpfully keep track of
driver income, but we can generate new passengers at each round.

###  III-D Algorithm Flow

The Gale-Shapley algorithm can be applied to this situation as follows:

  1. 1.

Each rider has a willingness-to-pay(WTP) and both a pickup and drop-off
location.

  2. 2.

Next, drivers see if utility, deemed to be price minus cost, is positive, then
proposes rides to riders based on their preference orderings

  3. 3.

Rider decides to accept/reject depending on wait time, and will do so in an
order based on their own preference orderings. The preference orderings for
both drivers and riders will take into consideration certain
heuristics/tradeoffs as detailed below.

###  III-E Algorithm Discussion

The Gale-Shapley algorithm, also known as deferred acceptance, results in a
stable matching. This means that there is no blocking pair consisting of a
passenger and driver who would prefer to be matched with each other than their
current pairing, based on their preference ordering. Additionally, the Gale-
Shapley algorithm is strategy-proof and optimal for the proposing side, which
in this case consists of the drivers. This means that it is in the driver’s
best interest to report their honest attributes and preferences. Typically,
the Gale-Shapley algorithm is not strategy-proof for the non-proposing side,
and so there can be useful deviations for those entities. However, some
concerns here are less relevant in this specific case, given the name of
ridesharing requests. Since riders are currently in one specific location and
want to get to another, then it would not make sense to misreport either of
these values as then

The area of greater concern could be in the underreporting of willingness-to-
pay (WTP) by passengers, as if they are able to be still matched with this
lower WTP then they would prefer to pay this lower amount if they are still
completing their ride, and thus having the same utility. We could potentially
fix this by having the payment amount be a critical value rather than a fixed
WTP, such that below this point, the ride would not be completed and the
passenger would not derive any utility.

We can make our algorithm more equitable if drivers care less about the drop-
off location’s distance from the center. If drivers are willing to prioritize
passengers that are farther away, then this would help overall driver income
equity but it would hurt the company’s overall revenue.

###  III-F Utilities Functions

Each passenger has a fixed willingness-to-pay (WTP), calculated as a linear
function of the Manhattan ride distance and the Euclidean distance from city
center, wrapped in a Gaussian distribution. In other words,

$\displaystyle WTP=\alpha((\mathtt{rideLength}+\mathtt{cityCenterDist}))/3$

where $\alpha$ is a coefficient randomly sampled from a Gaussian distribution
with mean = 1 and standard deviation = 0.1, $\mathtt{rideLength}$ is the
Manhattan distance from the passenger’s requested pickup and dropoff
locations, and $\mathtt{cityCenterDist}$ is the mean of the Euclidean
distances of the pickup and dropoff locations from the city center. The
stochasticity of $\alpha$ denotes that different riders may have different
tradeoffs and demand factors that are outside of our current scope, naively
capturing different people’s frugality. We calculated the Passenger
willingness ot pay like this to demonstrate that passengers are willing to pay
more for longer rides, as well as rides that make the driver go more out of
the way (rides that begin and end farther from the city center).

Similarly, each driver’s cost is a function of their total trip length and
distance from the city center, multiplied by a personal cost coefficient.

$\displaystyle\mathtt{cost}=\gamma(\mathtt{tripLength}+\mathtt{cityCenterDist})$

where $\gamma\in[0,1]$ is a cost coefficient randomly sampled from a Gaussian
distribution, $\mathtt{rideLength}$ is the sum of the Manhattan distances from
the driver’s current location to the passenger pickup location and the
passenger’s pickup to dropoff locations, and $\mathtt{cityCenterDist}$ is the
Euclidean distance of the dropoff locations from the city center. The cost
coefficient of $\gamma$ denotes that different drivers may have cost factors
that are outside of our current scope.

While our riders simply post their WTP, our drivers must take into
consideration the difference between a rider’s WTP and the cost of the trip.
In other words, we have simplified market equilibrium dynamics to be
constrained to just driver-side profit and utility. Our drivers have higher
costs for longer drives, as well as drives that result in them being father
away from the city center that would limit their ability to pick up new rides.

###  III-G Preference Orders

Both drivers and passengers take into consideration multiple weighted factors
when determining their preference orders. In our project, we experiment with
these hyperparameters for how much weight should be placed on these different
factors to achieve the best results. Here, lower preference scores are
prioritized.

Drivers take into account two factors: profit and distance from the city
center. Normalizing either factor, driver $d$ calculates the following
preference score for a passenger $p$:

$\displaystyle pref_{d}(p)=$
$\displaystyle-\mathtt{p.WTP}+w_{p}\mathtt{d.cost(p)}$
$\displaystyle\mathtt{d.cost(p)}=$
$\displaystyle\alpha(\mathtt{d.pickupDist(p)})$
$\displaystyle+0.1*(\mathtt{d.cityCenterDist(p)})^{2}$

where $w_{p}$ is a randomly generated gaussian random variable with mean = 1
and SD = 0.2 to demonstrate that different drivers value cost of the drive
different from each other. $\mathtt{d.pickupDist(p)}$ is the manhattan
distance between the driver and passenger’s pickup location and
$\mathtt{d.cityCenterDist(p)}$ is the euclidean distance between the
passenger’s dropoff location and the city center, normalized. Drivers
prioritize rides that earn them more money, which is why we add the WTP of the
passenger (the whole function is flipped by the negative since sorting in
python sorts smallest to largest), while they also want rides with the
smallest cost to them possible.

Passengers also take into account two factors when ordering drivers: wait time
and the driver’s current aggregate income. This is where our system attempts
to implement income equity among the drivers. By adding this factor into the
preference orders, the system allows drivers with lower aggregate income in
the current period to get higher priority in matching, thus hypothetically
promoting greater ride equity. We will test whether this is a plausible
solution or not in the simulation. Normalizing either factor, passenger $p$
calculates the following preference score for a passenger $d$:

$\displaystyle pref_{p}(d)=$ $\displaystyle w_{t}(\mathtt{waitLength(p,d)})$
$\displaystyle+w_{i}(\mathtt{aggIncome(d)})$

where $w_{t}$ is the weight applied to the pickup wait time (which, due to our
unit time per distance, is just the Manhattan distance between the driver’s
current location and the passenger’s requested pickup location) and $w_{i}$ is
the weight applied to the driver’s aggregate income thus far. When the system
does not take into consider income distribution, we set $w_{t}=1$ and
$w_{i}=0$. We hypothesize that the greater the $w_{i}$, the more equitable our
overall income distribution, at the tradeoff of our aggregate system revenue.

##  IV Simulation Design

###  IV-A Code

The Github repository containing our original simulation code can be found
here.

###  IV-B City Framework

We will simulate our city using a 100 x 100 Euclidean Grid, with the idea
being that it takes 1 time unit to travel a single Euclidean distance. An
important characteristic of many cities is the city center, which will assume
to be the direct center of the grid. This is the heart of the city where most
people and places of interest exist.

Thus, we generated locations to be Euclidean points where the x,y coordinates
span from -50 to 50. To account for the population densities, they are
simulated randomly from a normal distribution (with mean = 0 and standard
deviation = 20) to simulate that there are more people traveling to and from
the city center.

This is not only relevant for generating what the initial function of each
passenger and driver is, but it is only relevant because then drivers will
prefer rides that have drop-off locations closer to the city center rather
than those that end up in more remote locations. This is because then they
will be more likely to be closer to their own location of interest or to a new
ride.

###  IV-C Driver Generation

Each driver will be instantiated at the start of the first round. They will
have the following attributes, each initialized to 0, None, or the empty list:
$\mathtt{totalincome}$, $\mathtt{currentlocation}$, $\mathtt{totalrides}$,
$\mathtt{preferenceordering}$, $\mathtt{matchedpassengername}$, and
$\mathtt{locations}$, which is a list of locations that they have visited.
With each driver, we will also generate a unique cost coefficient and a
starting location, determined through the function of Euclidean distances from
the center touched upon above. Our set of drivers will stay constant
throughout the rounds, so that we can best analyze the effects on cumulative
driver income equity. In each subsequent round, we will update the attributes
intuitively based on the matching and completed ride. The cost coefficient
will stay the same for a particular driver between rounds and the new starting
location will just be the ending location of the previous ride if one was
completed.

###  IV-D Passenger Generation

Unlike drivers, we will generate a new set of passengers in each round. Each
will have the following attributes, initialized to 0, None, or the empty list
as best fits: $\mathtt{preferenceordering}$, $\mathtt{WTP}$,
$\mathtt{currentlocation}$, and $\mathtt{dropofflocation}$. Once the matching
and rides has been completed, we assume the passengers are no longer involved
in the next round and generate an entirely new batch, so no updating is
necessary directly to the passengers.

###  IV-E Design

For the sake of having enough participants but keeping our computational time
reasonable, we chose to have 15 drivers and 15 passengers in any given round,
and we iterated over 50 total rounds to see the long-term distribution of
income. In addition to specifying the number of rounds and the number of
agents per round, the simulation can also vary the Gaussian distribution of
passenger instantiation away from the city center, standard deviations of
driver cost functions and passenger WTPs, as well as the weight each agent may
assign its different priorities.

###  IV-F Random Matching

As a baseline, we explore what occurs when passengers and drivers are matched
randomly, which we implement by randomly shuffling both lists and then
creating matches in order.

###  IV-G Closest Neighbours

Next, we will apply an algorithm that first matches the pairing that are
closest together. More concretely, we implement this by constructing a
Euclidean distance matrix between passengers and drivers and then solving the
matching problem using a Hungarian algorithm. In this case, this is a central
allocation mechanism, in that it is neither rider-proposing nor driver-
proposing and instead Uber is acting as the third party to create the matches.

###  IV-H Gale-Shapley

The heart of our exploration focuses on the results that occur when we use the
Gale-Shapley deferred acceptance algorithm.

###  IV-I Altering Weights

Within this Gale-Shapley algorithm, we encapsulate various factors, such as
Manhattan distance, cost, distance from city center, and driver income
distribution, into one equation under preference orderings for passengers’
ranking of drivers. In this process, we assign different weights to these
factors, corresponding to how much we care about each factor. Part of our
results will be in modifying these weights and seeing how this leads to
different overall distributions and what tradeoffs surface.

##  V Results

For the most thorough understanding of the performance of the deferred
acceptance algorithm, we will compare the results against that of Random
Matching and Closest Neighbors algorithms. First, we will output a chart of
results of the different algorithms against each other. First, let us look at
the total cumulative incomes of the four different algorithms

DA Driver prop | Random | Closest | Boston Algo
---|---|---|---
14644.953 | 10084.181 | 12839.953 | 8196.113
15934.399 | 9555.306 | 12956.324 | 11372.515
17567.636 | 9247.306 | 13789.982 | 7963.459
18106.353 | 10164.92 | 13058.475 | 12930.511
17168.541 | 9536.810 | 14708.658 | 12343.745
17131.983 | 9469.655 | 12671.108 | 11733.939
13544.27 | 9469.239 | 11482.561 | 9644.839
16419.518 | 8971.500 | 12421.540 | 12778.096
15864.458 | 10368.144 | 13716.894 | 12983.909

As we can see from the chart, the driver proposing DA algorithm always
produced the greatest aggregate income, followed by Closest (closets drivers
and passengers are matched without any side proposing), then the Boston
algorithm where drivers simply propose without passengers able to switch,
followed by random. This makes sense since drives are rejected by passengers
when wait time is too long, which results a lot more often in random and
Boston than Closest or DA.

###  V-A Graph comparing Different Algorithms

These graphs plot the title quantity in terms of number of agents (total
drivers plus passengers) in a given time period.

Figure 1: Mean income per driver vs. number of agents. In the randomly matched
algorithm, we see that mean income stays around the same at $600 per driver,
which is as expected, since the ratio of drivers and passengers stay 1-to-1.
In contrast, we see that the Boston matching algorithm demonstrates greater
variability with increased number of agents, since given the instability,
there is increased capacity for the algorithm to output suboptimal matchings
between drivers and passengers, where passengers may not give input on better
wait times, which may imply better allocation. Then, we see that the closeness
and DA matching algorithms are better equipped to capture greater revenue as
agents increase. These results demonstrate that proximity is a driver of
revenue, and the DA algorithm builds on this revenue driver by further
considering of profit and city-center-clustering. Figure 2: Revenue per ride
vs. number of agents. Here, we see that, despite significantly better revenue
performance for closeness matching and DA, we see that the the closeness
matching and random matching algorithms maintain relatively similar per-ride
revenue. Figure 3: Standard deviation of income per driver vs. number of
agents. As expected, standard deviation of driver income of the non-random
algorithms increase with the increased number of agents, while the random
algorithms maintains around the same variability. Figure 4: Total revenue of
system vs. number of agents. This graph further supports our findings in
previous figures, which demonstrate that the DA algorithm best capitalizes on
more agents’ preferences to maximize revenue for drivers. We also continue to
see that the Boston matching algorithm demonstrates higher variability due to
instability.

###  V-B Graph comparing DA with Fairness

Figure 5: Total revenue based on fairness vs. number of agents. Interestingly,
we find that our hypothetically fair algorithm (DA with driver income taken
into consideration for preference ordering) achieves higher aggregate revenue
compared to the vanilla DA algorithm. This result would suggest that equitable
income distribution not only has positive implications on the labor market and
sustainable business practices, but also increases welfare and revenue
generation. Still, we may have to do further work to confirm the drivers for
this revenue-driving equity. Figure 6: Revenue per ride based on fairness vs.
number of agents. Again, we see our hypothetically fair algorithm performing
better in increasing revenue per ride. Figure 7: SD based on fairness vs.
number of agents. The hypothetically fair algorithm sees slightly greater
standard deviation in income per driver.

###  V-C Graphs with new fairness weights

Figure 8: Revenue total vs. number of agents. Drivers better optimize revenue
by prioritizing an optimal dropoff location (one closer to the city center,
where they are likely to book the next passenger with less travel cost) over
naive proximity to passenger. This result suggests that the clustering of
drivers around passenger hot spots is significant in generating greater
revenue. Therefore, in order to incentivize drivers to drop off passengers in
more remote locations, the system may have to take further measures at
distributing income. Figure 9: Revenue per ride vs. number of agents This
graph supports the one above, showing a strong difference between the revenue
per ride for different weights. In particular, total revenu is increased when
drivers care more about the dropoff location’s distance from the center, but
this might be at the tradeoff of less equity across drivers. Figure 10:
Standard deviation vs. number of agents As we see, the standard deviation of
the driver incomes where drivers put more weight onto dropoff location of
passengers is higher than the sandard deviation of the driver incomes that put
more weight in passenger proximity. This demonstrates that putting more focus
on reaching driver with farther locations to travel can bring more income
equity, lowing the variance of driver incomes and flattening driver income
distribution.

##  VI Reflections and Conclusion

In conclusion, we found that the Gale-Shapley deferred acceptance (DA)
algorithm was better equipped to optimize driver revenue than other matching
algorithms, such as closeness matching or the Boston algorithm. We also found
that stability is important for better income equity, as the Boston algorithm
showed high variability for most metrics of results.

Surprisingly, we found that our hypothetically fair algorithm (which is the DA
algorithm modified to prioritize drivers with lower aggregate income) was able
to generate higher revenue than the vanilla DA algorithm. This may imply that
pursuing income equitability is not only a labor-sustainable practice, but
also something that may achieve revenue upsides.

The ridesharing system is complex and hard to comprehensively evaluate in any
one research endeavor, much less one for an undergraduate class. As such,
there are several further points of study. Of interest to us is implementing a
dynamic or competitive pricing model, where prices increase if demand is
greater than the supply, such as through origin-based surge pricing.
Additionally, we would like to model behavior at abnormal times of driving,
such as in the later hours of the night or in rush hours like post-work or
around holidays like Thanksgiving. This could give insights into how changes
in demand and driver density could affect driver equity and treatment.

##  VII Acknowledgements

## References

  * [1] Tom Sühr, Asia J. Biega, Meike Zehlike, Krishna P. Gummadi, and Abhijnan Chakraborty. 2019. Two-Sided Fairness for Repeated Matchings in Two-Sided Markets: A Case Study of a Ride-Hailing Platform. In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining (KDD ’19). Association for Computing Machinery, New York, NY, USA, 3082–3092.
  * [2] Zhao, B., Xu, P., Shi, Y., Tong, Y., Zhou, Z., & Zeng, Y. (2019). Preference-Aware Task Assignment in On-Demand Taxi Dispatching: An Online Stable Matching Approach. Proceedings of the AAAI Conference on Artificial Intelligence, 33(01), 2245-2252.