published: Clas. Quantum. Grav. 39, 225006 (2022),
doi.org/10.1088/1361-6382/ac987f

#  The stability of galaxies in an expanding universe obtained by Newtonian
dynamics

Søren Toxvaerd  Department of Science and Environment, Roskilde University,
Postbox 260, DK-4000 Roskilde, Denmark

###### Abstract

The dynamics of galaxies in an expanding universe is often determined for
gravitational and dark matter in an Einstein-de Sitter universe, or
alternatively by modifying the gravitational long-range attractions in the
Newtonian dynamics (MOND). Here the time evolution of galaxies is determined
by simulations of systems with pure gravitational forces by classical
Molecular Dynamic simulations. A time reversible algorithm for formation and
aging of gravitational systems by self-assembly of baryonic objects, recently
derived (Eur. Phys. J. Plus 2022, 137:99), is extended to include the Hubble
expansion of the space. The algorithm is stable for billions of time steps
without any adjustments. The algorithm is used to simulate simple models of
the Milky Way with the Hubble expansion of the universe, and the galaxies are
simulated for times which corresponds to more than 25 Gyr. The rotating
galaxies lose bound objects from time to time, but they are still stable at
the end of the simulations. The simulations indicate that the explanation for
the dynamics of galaxies may be that the universe is very young in
cosmological times. Although the models of the Milky Way are rather stable at
13-14 Gyr, which corresponds to the cosmological time of the universe, the
Hubble expansion will sooner or later release the objects in the galaxies. But
the simulations indicate that this will first happen in a far away future.

##  I Introduction

The constituents of the universe are often found in clusters of objects, in
planetary systems and in galaxies. The galaxies with stars were established
shortly after the ”Big Bang” Kolb2018 . The two groupings of objects have a
common feature, the objects rotate around their center of gravity. But a
rotating galaxy seems, however, to be unstable unless it contains a
significantly larger mass (dark matter) than given by the known baryonic
matter in the galaxy Martino2020 .

The universe expands with a Hubble velocity Hubble1929 , and the expansion of
the universe affects the dynamics of objects with gravitational forces.
Simulations of galaxies with dark matter in an expanding universe have been
performed for many decades. In Klypin1983 ; Centrella1983  the authors solved
the dynamics of baryonic objects by the Particle-Particle/Particle-Mess (PPPM)
method Hockney1974 , where each mass unit is treated as moving in the
collective field of all others, and the Poison equation for the PPPM grid is
solved numerically. The dynamics are obtained for Zeldovich’s adiabatic
expansion of the Einstein-de Sitter universe Zeldovich1970 . Later,
simulations with large scaled computer packages with PPPM Diemand2008 ;
Alimi2012  are with many billions of mass units. The evolution of galaxies has
also been obtained from hydrodynamical large scaled cosmological simulations
Schaye2010 ; Dubois2014 ; Vogelsberger2014 , and with co-evolving dark matter
gas and stellar objects Dubois2016 ; Ludlow2021 . Cosmological simulations of
galaxy formation is reviewed in Vogelsberger2019 .

An alternative theory for the dynamics of the galaxies, but without dark
matter is the modified Newtonian dynamics (MOND), where Milgrom assumed that
the very weak gravitational force experienced by a star in the outer regions
of a galaxy varies inversely linearly with the distance, as opposed to the
inverse square of the distance in Newton’s law of gravity Milgrom1983 .
Merging galaxies with MOND and with galaxies embedded in dark matter shows
differences in the merging time-scales, but only very small differences in the
final stellar distribution Nipoti2007 . The MOND simulations techniques for
modified gravity is given in Llinares2018 . The MOND theory is extensively
tested and compared with the Lambda cold dark matter model $\rm{\Lambda}$CDM
and in favour of the MOND theory Kroupa2015 ; Banik2022 .

Newtonian cosmology is derived in Ellis2014 ; Ellis2015 . Galactic dynamics
with pure Newtonian gravity forces has been compared with $N$-body MOND
simulations, and in favour of gravitational dynamics Wu2007 ; Derakhshani2014
, and Angus et al. find similarly some shortcomings (mass-to-light ratios) of
MOND at simulations of the Carina Dwarf Spheroidal Galaxy orbiting the Milky
Way Angus2014 . The Newtonian dynamics simulations of galaxy systems described
below complement the existing analyzes of the stability of galaxies. The
$N$-body simulations are for models of celestial objects in the Milky Way, and
the objects interact with pure gravitational forces and with and without the
Hubble expansion of the space.

The dynamics of models for galaxies are based on an algorithm for the
classical dynamics of systems of celestial objects Toxvaerd2022 . Here (see
Appendix A) the algorithm is extended to include the Hubble expansion of the
space. The discrete algorithm is time reversible and stable for billions of
time steps without any adjustments. The algorithm is used to simulate simple
models of the Milky Way with and without the Hubble expansion of the universe.
The different models of galaxies are simulated for times which correspond to
more than twice the age of the Universe and both the galaxies without- as well
as galaxies with the Hubble expansion of the space are stable. The simulations
indicate that an alternative explanation to the dynamics of galaxies with dark
matter is, that the universe is very young in cosmological time units given by
one rotation of the galaxy. A galaxy as old as our universe lose bound
objects, but although the Hubble expansion sooner or later will release the
stars in the Milky Way, the simulations show that it will first happen in the
distant future.

##  II The expanding universe

A galaxy, $l$, far away in the universe moves away from the Earth at a speed
proportional to its ”proper distance”, $Hr_{kl}(t)$, to the Earth, $k$,
measured at the ”cosmological time” $t$, and this behaviour is explained by an
expanding universe. The Hubble constant $H$ is the expansion coefficient in
Hubble’s law Hubble1929 ; Lan2022

$\textbf{v}^{H}(\textbf{r}_{kl}(t))=H\textbf{r}_{kl}(t)$ (1)

for the velocity of the increase of the distance $\textbf{r}_{kl}(t)$ from the
Earth to a galaxy. The Hubble expansion can be obtained by an intensive
expansion of the space independently of the baryonic matter in the universe

$\textbf{v}^{H}(\textbf{r}(t))=H\textbf{r}(t),$ (2)

by which the distance between pairs of positions
$\textbf{r}_{k}(t),\textbf{r}_{l}(t),$ or
$\textbf{r}_{k}(t),\textbf{r}_{k}(t+\delta t)$ increases with the Hubble
velocity. Eq. (2) fulfills the Cosmological principle and the Copernican
principle for expansion of the universe Book . The Cosmological principle,
which was first formulated by Newton Newton1687  demands, that no place in the
universe is preferred, and the Copernican principle demands, that no direction
in the universe is preferred. Eq. 2 will accelerate the expansion of the
universe, and observations of galaxies indicates, that the distances between
the galaxies accelerate (dark energy)Riess1998 ; Perlmutter1999 .

The Hubble constant $H$ is quoted in $\textrm{km s}^{-1}\textrm{Mpc}^{-1}$,
for the velocity in $\textrm{km s}^{-1}$ of a galaxy 1 megaparsec ($3.09\times
10^{19}$ km) away. Its value is Soltis2021

$H=72.1\pm 2.0\ \textrm{km}\ \textrm{s}^{-1}\textrm{Mpc}^{-1}.$ (3)

For the Solar system with distances in au (1 au$\approx$ Earth’s distance to
the Sun =149597870.700 km (definition)) the Hubble constant is $33.93\times
10^{-11}$ $\textrm{km s}^{-1}\textrm{au}^{-1}$, and the Earth’s mean velocity
in the Solar system is $\approx 30$ $\textrm{km s}^{-1}$, so the Hubble
velocity is about $\approx 10^{-11}$ of the Earth’s velocity. This small
Hubble velocity has no effect on the orbits of the planets and the stability
of a planetary system (see later).

The Milky Way is a barred spiral galaxy. The extension of the barred disk is
2.5-3 kpc Rix2013 , and the extension of the halos is $\approx$ 100 000 pc,
but some recent investigations indicate that this might be substantial larger,
$\approx$ 300 kps or even more Deason2020 ; Li2021 . The Milky Way rotates
with a rotational velocity $220\pm 10$ $\textrm{km s}^{-1}$ for stars at a
distance of $8.0\pm 0.310$ kpc from its center Camarillo2018 , but the
rotations at different distances from the center of the galaxy deviate from
calculations of the rotation caused by the gravitational forces between its
constituents, and this deviation is explained by the existence of unknown dark
matter Weinberg1984 ; Chiba2020 ; Martino2020 . The Milky Way contains very
old stars from the period just after the ”Big Bang” Bond2013 , so the galaxy
has presumably exist for more than thirteen billion years. The rotational time
(galactic year) of the Milky Way is estimated to be $\approx$ 220-240 million
years, and this implies that an old star in the Milky Way with this rotational
velocity must have performed of the order $\approx$ 60 cycles.

The Hubble velocity between pairs of stars in the Milky Way at distances of
e.g. $5\times 10^{-2}$ Mpc is 3.5 $\textrm{km s}^{-1}$, and it is not
negligible compared with the rotational velocity of the galaxy. The expansion
of the universe is with a Hubble velocity, which is of the order $\approx 1\%$
or more of the rotational velocity of the Milky Way, and the Hubble expansion
will affect the rotation and the stability of the galaxy.

Figure 1:  Positions of the objects. A: Site view of the system; B: Top view.
The green dots are the positions at the start and the red dots are the
positions at $t=2.5\times 10^{4}(10^{7}$ time steps).

##  III Simulation of galaxies in an expanding universe

A galaxy and the Milky Way contains hundred of billion of stars, and a
substantial amount of baryonic gas Gupta2012 ; Bergma2018  and it is not
possible directly to obtain the Newtonian dynamics of a galaxy with this
number of objects. Instead, we have simulate models of small ”galaxies” of
hundred of objects in orbits around their center of gravity, and in an
expanding space with various strength of the expansion. A recent article
describes how to create a collection of objects with rotation about their
center of gravity by Molecular Dynamics (MD) simulations, obtained by an
extension of Newton’s discrete algorithm Toxvaerd2022 . The algorithm is
further developed to also include the Hubble expansion of space, and the
algorithm with the extensions are given in the Appendix.

The collection of baryonic objects with Newtonian dynamics will, for many
starting positions and velocities of the objects spontaneously form a system
with many of the objects in rotations about a heavy mass in its center. If
these MD systems shall simulate the dynamics of a galaxy in the expanding
universe, then one must relate distances, times, and Hubble expansion in the
MD systems with the corresponding distances, times and Hubble expansion in a
galaxy.

###  III.1 Determination of distances, times and $H$ for the Hubble expansion
in the MD systems

A system $without$ $a$ $Hubble$ $expansion$ of the space is obtained from
ensembles of $N=1000$ objects with a disk-shaped distribution of their start
positions and with $\approx$ zero velocities in the z-direction. The
distribution of the objects at the start and after the creation of a system
with $\approx 300$ of the objects in circulations around a heavy object at the
center of mass is shown in Figure 1.

Figure 2:  $log$ mean distances to the mass center as a function of the $log$
mean velocities of the objects. The bound objects (red spheres) are located on
the lower branch of the distribution and the upper branch shows (blue spheres)
the mean locations of the free objects. Magenta line:
$log(y)=log(Gm_{c}/x^{2})$.

The system with positions shown in Figure 1 is simulated for a long period of
time. The system is evolving and with mergers of objects, and after
$t=2.5\times 10^{6}$ it contains $N=557$ objects, and with 270 of them in
circulations around the heavy mass in the center. The mean distances,
$\bar{r}_{i,c}$, in the succeeding time interval $\Delta t\in[2.5\times
10^{6},3.0\times 10^{6}]$, as a function of the mean velocity,
$\bar{v}_{i,c}$, relative to the center $c$ is shown in Figure 2. A bound
object will have a relative small mean distance to the center of mass in
contrast to a free object, which has a large and monotonic increasing distance
to the center. The lower branch of mean distances in Figure 2 are for the
collection of bound objects, and the the upper branch are the mean distances
of the free objects moving away from the center with constant velocities. An
object, $i$, in an elliptic orbit around a heavy gravitational center with the
mass $m_{c}$ will have a mean velocity $\bar{v}_{i,c}=\sqrt{Gm_{c}/a_{i}}$,
where $a_{i}$ is the semimajor axis in its ellipse. The magenta line is the
line $log(y)=log(Gm_{c}/x^{2})$ with the mass $m_{c}=340$ of the center of the
system (for units in the MD systems see the Appendix). The systematic
deviation of the mean distances of the bound objects from the straight line
for increasing distances to the center is partly caused by, that the
eccentricities of the orbits of the bound objects increase with increasing
distance to the center of mass, but the distribution of objects for greater
distances $2000<\bar{r}_{i,c}$ shows, that the objects are not lined up in a
row of objects in elliptical orbits as in a planetary system.

The mean distances, rotation times and velocities for the bound objects, shown
in Figure 2 are used to relate the MD system with the Milky Way.

The Milky Way has an extension of $\approx 100000$ pc, but the extension of
some of the halos may be more. The lower branch in Figure 2 with the bound
objects is within a MD distance $100000$ MD length units or less from the
center. So a relation between the distances in the Milky Way in pc and the
distances in the MD system is 1 ps= 1 MD length unit.

The galactic year of the Milky Way is about 220 million years Camarillo2018 ,
and the Milky Way contains a very old star Bond2013 , which has been created
shortly after the ”Big Bang”. A star with a rotation time of 220 million
years, must have performed only $\approx$ 60 cycles in the Milky Way after it
is created. The distance $R_{0}$ of our planet system to the center of the
Milky way is $R_{0}=8.34\pm 0.16$ kpc Reid2014 . These data are used to relate
the cosmological time with the time in the MD systems. The rotation time of
objects in the MD system at $\approx$ 8000 MD length units is
$t_{orbit}\approx 1.3\times 10^{5}$, so these objects have performed $\approx
60$ rotations at the MD time $8.0\times 10^{6}$. The Figure 2 is for 0.3125 of
the cosmological time of the universe, or $\approx$ 4.2 Gyr, and the age of
the universe corresponds to $\approx$ $8.0\times 10^{6}$ in the MD systems.

Figure 3: The $log$ mean distances $log(\bar{r}_{i,c})$ of the objects as a
function of their $log$ mean velocities $log(\bar{v}_{i,c})$ with respect to
the center $c$ for different values of the Hubble constant $H$. The mean
positions are obtained for the time interval $\Delta t\in[2.5\times
10^{6},3.0\times 10^{6}]$.

The mean velocity $\bar{v}$ of the bound objects in the MD system in Figure 2
is $\bar{v}\approx 0.15$. The value of the Hubble expansion constant in the MD
system is obtained from the ratio $H/\bar{v}$ between the Hubble constant per
unit length and the rotational velocity of the galaxy as

$(H/\bar{v})_{\textrm{MD}}=(H/\bar{v})_{\textrm{Milky Way}},$ (4)

which determines the value of $H$ in the MD system,

$\displaystyle H=0.15\times 72\times 10^{-6}\textrm{km}\ \textrm{s}^{-1}\
\textrm{pc}^{-1}/220\ \textrm{km}\ \textrm{s}^{-1}$ $\displaystyle\approx
5.\times 10^{-8}\textrm{ per MD length unit}.$ (5)

Let the estimates of the extent of the Milky Way, the extent of the MD model
of a galaxy, the rotation time of the Milky Way and the mean velocity in the
MD system all be with an uncertainty with a factor of two, then the total
uncertainty of the estimate of $H$ in the MD systems is $\approx 2^{4}\approx
10$. The value of the Hubble constant in the MD units for the model of the
Milky Way in the expanding universe is

$H=5.\times 10^{-8\pm 1}.$ (6)

Figure 4: The release of the last bound object $l$ from the system with
$H=5.\times 10^{-7}$. Blue: The orbit for $\bar{r}_{l,c}(t)<15000$ and green:
for $\bar{r}_{l,c}(t)>15000$. The center of gravity with red. The blue sphere
is the position at $t=6.5\times 10^{6}$, and the green sphere is the position
at $t=7.75\times 10^{6}$.

###  III.2 The formation and stability of galaxies obtained by classical
dynamics

The MD system without an expansion of the space, with positions shown in
Figure 1 at $t=0$ (green) and $t=2.5\times 10^{4}$ (red), and Figure 2 at
$t=2.5\times 10^{6}$ is formed spontaneously. This system is used to
investigate the formation and stability of MD systems to various values of $H$
in the Hubble expansion. Systems with $H=0,$ $5.\times 10^{-8},$ $5.\times
10^{-7}$ and $5.\times 10^{-6}$, respectively, are simulated and their
behavior is shown below.

In Figure 3 the distributions of objects at $t=2.5\times 10^{6}$ is shown for
various strength of the Hubble constant. The four systems are started with the
disk formed distribution of 1000 objects with equal masses, shown with green
dots in Figure 1. All four systems spontaneously formed gravitational systems.
The systems contain many bound objects with different masses and in
circulation around their centre of gravity, but the system with $H=5.\times
10^{-6}$ was unstable and with no bound objects already after $t=1.0\times
10^{6}$. The radius of the disk in the Milky way is $\approx 15$ kpc Rix2013
corresponding to a a distance 15000 MD length units in the MD systems. At
$t=2.5\times 10^{6}$, the number of objects within this distance to the center
of the galaxy are 144, 152, 19 and 0, respectively for the four different
values of the expansion. The system with $H=5.\times 10^{-7}$ and only 19
objects with $\bar{r}_{i,c}<15000$ is, however, also unstable and without any
objects with distances less than 15000, but first after $t=6.5\times 10^{6}$
(Figure 4).

Figure 5:  Number of objects $i$ with mean distances $\bar{r}_{i,c}(t)<15000$
to the center of the galaxy $c$ as a function of time. The systems are started
from the ”gas” distribution shown with green dots in Figure 1. The colors
(same as in Figure 3) are for different strength of of the Hubble expansion.
The inset to the left shows the time evolution of the $L_{z}$-components of
the angular momenta of the galaxies with different strength of the Hubble
expansion. In the right inset the number of bound objects with
$\bar{r}_{i,c}(t)<15000$ for $H=5.\times 10^{-8}$ are compared with the total
number of bound objects $\bar{r}_{i,c}(t)<10^{5}$ in the galaxy.

Figure 5 shows the time evolution of the number of objects within a mean
distance $\bar{r}_{i,c}<15000$ from the center $c$. Almost all of the objects
are within this distance at the start of the evolution, but as the self-
assembly and aging of the Newtonian systems evolves, the number declines
slowly for $H\leq 5.\times 10^{-8}$, whereas it rapidly goes to zero for
faster expansions. The arrow to the left marks the time $2.5\times 10^{6}$
with the distribution of objects shown in Figure 3, and the arrow in the
middle is at the time $8.0\times 10^{6}$, which should correspond to the age
of our universe. The systems with $H\leq 5.\times 10^{-8}$ are rather stable
at this time, but with a weak declining number of objects within the distance
$\bar{r}_{i,c}<15000$. The systems are started from the same start
distribution, shown with green spheres in Figure 1, and the self-assembly with
merging of objects results in slightly different masses of the center of the
galaxies and small differences between the numbers for $H=0$ and $H=5.\times
10^{-8}$. The left inset in Figure 5 shows the $L_{z}$ components of the
angular momenta L for the galaxies with the different strength of the
expansion. The angular momentum of a system without a Hubble expansion is
conserved for Newton’s discrete dynamics (blue curve in the inset), but the
Hubble expansion increases the angular momenta, and this will sooner or later
destroy the galaxies (Figure 4). The number of bound object with mean distance
$\bar{r}_{i,c}(t)<15000$ to the center of the galaxy with $H=5.\times 10^{-8}$
is compared with the corresponding total number of all bound objects ( objects
with $\bar{r}_{i,c}(t)<10^{5}$ ) in the galaxy in the right inset of Figure 5.
The two sets of numbers decreases with time, but very slowly and the
destruction of the system with first happen in a far away future for the
galaxy with $H=5.\times 10^{-8}$.

A series of control simulations have been performed in order to ensure, that
the results shown in Figure 5 are representative for the formation and aging
of gravitational systems with Hubble expansion. In one set of test simulations
the three different expansions was first included in the dynamics at
$t=2.5\times 10^{6}$ in the system without expansion. But with the same result
that the system with $H=5.\times 10^{-8}$ was rather stable, but the two
systems with stronger Hubble expansions were unstable. In another set of
simulations the Hubble expansion was added to a Newtonian system with another
distribution of the objects at the start of the formation of the galaxy, but
with the same result: All the systems are simulated for times much longer,
that corresponds to ours cosmological time $\approx$13.4 Gyr, and the systems
are still rather stable for an expansion with $H\leq 5.\times 10^{-8}$, which
should correspond to the Hubble expansion of the universe.

Figure 6:  $log(\bar{r}_{i,c})$ as a function of $log(\bar{v}_{i,c})$ and at
different times. With red is for $t=2.5\times 10^{6}$, and with green is for
$t=8.\times 10^{6}$ ($\approx$ 13.4 Gyr). $\textrm{N}_{\textrm{disk}}$ are for
$\bar{r}_{i,c}<15000$ and $\textrm{N}_{\textrm{halos}}$ are for
$15000<\bar{r}_{i,c}<100000$. $<\bar{r}_{i,c}>$ is the mean positions for the
$\textrm{N}_{\textrm{disk}}+\textrm{N}_{\textrm{halos}}$ objects. The blue
line, $r=m_{c}/v^{2}$, is the line for mean positions of objects in ellipses,
and with the same eccentricities.

Figure 6 shows the $log$ mean distances $log(\bar{r}_{i,c})$ of the objects as
a function of their $log$ mean velocities $log(\bar{v}_{i,c})$ with respect to
the center $c$ of the system with $H=5.\times 10^{-8}$, and at different
times. With red is for $t=2.5\times 10^{6}$, also shown in Figure 3, and with
green is for $t=8.\times 10^{6}$ ($\approx$ 13.4 Gyr). The numbers
$\textrm{N}_{\textrm{disk}}$ are for objects with $\bar{r}_{i,c}<15000$ and
the numbers $\textrm{N}_{\textrm{halos}}$ are for objects with
$15000<\bar{r}_{i,c}<100000$. $<\bar{r}_{i,c}>$ is the mean distances for the
$\textrm{N}_{\textrm{disk}}+\textrm{N}_{\textrm{halos}}$ objects. The blue
line, $log(r)=log(Gm_{c}/v^{2})$ is the line for mean distances of objects in
ellipses with the same eccentricities. The distributions of the objects
deviates from this line for a ”Kepler-like” order with the objects in
unperturbed elliptical orbits around the center of mass. In particular, the
objects with large distances to the center of gravity deviate from the
straight line by having in general a higher velocity than given by the blue
line, and the aging of the system does not remove this tendency. The system at
$t=8.\times 10^{6}$ should correspond to the cosmological time 13.4 Gyr of the
universe, and the systems are still aging with losses of bound objects and
with a monotonic increasing rotational velocities, as shown in the inset of
Figure 5 for $L_{z}(t)$.

##  IV Discussion

The dynamics of galaxies in an expanding universe are often determined for
gravitational and dark matter in an Einstein-de Sitter universe, or
alternatively by the modified Newtonian dynamics (MOND)Milgrom1983 . But here
the dynamics of galaxies is obtained for systems with pure gravitational
matter by classical Molecular Dynamic simulations. The dynamics of the systems
is not explained by dark matter, nor by a modification of the gravitational
attractions at long distances, but by that the universe is a young universe
measured in cosmological times. The Milky Way is believed to be more than 13
Gyr, but it corresponds to only $\approx$ 60 rotations or galactic years, and
it is nothing compared to the age of our solar system with 4.56 Gyr, i.e.
number of rotations of the Earth. The stability of the solar system is still
to debate Batygin2008 ; Mogavero2021 , and in the light of these facts the
universe must be characterized as a very young universe. The dynamics and
aging of the MD systems with gravitational matter shows, however, that the
objects in the small MD systems are not in a stable steady state after only
sixty rotations, so how can hundred of billion of stars and a substantial
amount of baryonic gas then be it?

According to the present investigation, the explanation for the dynamics of
the universe could be, that the galaxies are formed spontaneously at the
relatively high concentration of baryonic matter shortly after the Big Bang.
The universe expands and the galaxies are aging, but the galaxies are not in a
steady state at a time which corresponds to 13.4 Gyr (Figure 5, Figure 6). The
systems with- and without a Hubble expansion still occasionally lose bound
objects at 13.4 Gyr, and the mean distance of the bound objects increase also
over time. The bound objects far away from the center of gravity are not
distributed as in a stable Newtonian system. They have in general a higher
velocity, than if they were ordered in simple Kepler-like orbits (Figure 6).
The hypothesis can perhaps be tested by statistical analyses of the density
and velocity distribution (distribution of light intensity) in galaxies as
functions of the age of the galaxies.

The Hubble expansion increases the angular momentum of a galaxy and this fact
implies, that the galaxy sooner or later will be destroyed, and that all the
objects will be free in an open universe. The end of the system with
$H=5.\times 10^{-7}$ is shown in Figure 4, with the release of the last bound
object in the system. But this system has a strength of the Hubble expansion,
which probably is ten times stronger than in the universe. The simulations
indicate, however, that for the systems with a Hubble expansion similar to the
Hubble expansion of the universe, this destruction might first happen in a far
away future. These systems are simulated for what corresponds to more than 25
Gyr, and they still contain many objects with $\bar{r}<15000$ (Figure 5).

The algorithm is well suited for analyzes of formation and aging of
gravitational systems in an expanding universe. The central difference
algorithm, Eqn. (A11)+ (A3), is time reversible, and the underlying Newtonian
discrete algorithm, which is used in almost all MD simulations in physics and
chemistry is exact in the sense, that it has the same invariances as Newton’s
analytic dynamics Toxvaerd2022 ; Toxvaerd2014 . The algorithm is stable, the
present simulations are with more than $6\times 10^{9}$ discrete time steps
for each systems, and each simulation has taken 4-5 months, but without any
adjustments.

##  Appendix A A time reversible algorithm for discrete classical dynamics in
an expanding space

###  A.1 Newton’s discrete central difference algorithm

In Newton’s classical discrete dynamics Newton1687 ; Toxvaerd2020  a new
position $\textbf{r}_{k}(t+\delta t)$ at time $t+\delta t$ of an object $k$
with the mass $m_{k}$ is determined by the force $\textbf{f}_{k}(t)$ acting on
the object at the discrete positions $\textbf{r}_{k}(t)$ at time $t$, and the
position $\textbf{r}_{k}(t-\delta t)$ at $t-\delta t$ as

$m_{k}\frac{\textbf{r}_{k}(t+\delta t)-\textbf{r}_{k}(t)}{\delta
t}=m_{k}\frac{\textbf{r}_{k}(t)-\textbf{r}_{k}(t-\delta t)}{\delta t}+\delta
t\textbf{f}_{k}(t),$ (7)

where the momenta $\textbf{p}_{k}(t+\delta t/2)=m_{k}(\textbf{r}_{k}(t+\delta
t)-\textbf{r}_{k}(t))/\delta t$ and $\textbf{p}_{k}(t-\delta
t/2)=m_{k}(\textbf{r}_{k}(t)-\textbf{r}_{k}(t-\delta t))/\delta t$ are
constant in the time intervals in between the discrete positions. Newton
postulated Eq. (A1) and obtained his second law, and the analytic dynamics in
the limit $lim_{\delta t\rightarrow 0}$.

The algorithm, Eq. (A1), is usual presented as the ”Leap frog” algorithm for
the velocities

$\textbf{v}_{k}(t+\delta t/2)=\textbf{v}_{k}(t-\delta t/2)+\delta
t/m_{k}\textbf{f}_{k}(t).$ (8)

The positions are determined from the discrete values of the
momenta/velocities as

$\textbf{r}_{k}(t+\delta t)=\textbf{r}_{k}(t)+\delta t\textbf{v}_{k}(t+\delta
t/2).$ (9)

Newton’s discrete algorithm, Eq. (A1), has been re-derived many times and with
different names, and almost all Molecular Dynamics (MD) simulations in physics
and chemistry as well as many simulations in astrophysics are performed with
his discrete algorithm. The rearrangement of Eq. (A1) gives e.g. the ”Verlet”
algorithm Verlet1967 .

The classical discrete dynamics between $N$ spherically symmetrical objects
with masses $m^{N}=m_{1},m_{2},..,m_{N}$ and positions
$\textbf{r}^{N}(t)=\textbf{r}_{1}(t),\textbf{r}_{2}(t),..,\textbf{r}_{N}(t)$
is obtained by Eq. (A2) and Eq. (A3). The momentum and angular momentum for a
conservative system of the $N$ objects is conserved due to Newton’s third law,
and the energy is also conserved Lee1983 ; Toxvaerd1994 ; Toxvaerd2014 , so
Newton’s time reversible and symplectic discrete dynamics has the same
invariances as his analytic dynamics. Newton’s discrete algorithm make it
possible to simulate systems of interacting objects for billions of time steps
without $any$ adjustments for the infinite range of the gravitational forces
or corrections for symplecticity, time reversibility, conservation of
momentum, angular momentum and energy.

###  A.2 Newton’s discrete central difference algorithm with merging of
gravitational objects

The discrete central difference algorithm for merging of objects is derived in
Toxvaerd2022 . Let all the spherically symmetrical objects have the same
(reduced) number density $\rho=(\pi/6)^{-1}$ by which the diameter
$\sigma_{i}$ of the spherical object $i$ is

$\sigma_{i}=m_{i}^{1/3}$ (10)

and the collision diameter

$\sigma_{ij}=\frac{\sigma_{i}+\sigma_{j}}{2}.$ (11)

If the distance $r_{ij}(t)$ at time $t$ between two objects is less than
$\sigma_{ij}$ the two objects merge to one spherical symmetrical object with
mass

$m_{\alpha}=m_{i}+m_{j},$ (12)

and diameter

$\sigma_{\alpha}=(m_{\alpha})^{1/3},$ (13)

and with the new object $\alpha$ at the position

$\textbf{r}_{\alpha}(t)=\frac{m_{i}}{m_{\alpha}}\textbf{r}_{i}(t)+\frac{m_{j}}{m_{\alpha}}\textbf{r}_{j}(t),$
(14)

at the center of mass of the the two objects before the fusion. (The object
$\alpha$ at the center of mass of the two merged objects $i$ and $j$ might
occasionally be near another object $k$ by which more objects merge, but after
the same laws.)

The momenta of the objects in the discrete dynamics just before the fusion are
$\textbf{p}^{N}(t-\delta t/2)$ and the total momentum of the system is
conserved at the fusion if

$\textbf{v}_{\alpha}(t-\delta
t/2)=\frac{m_{i}}{m_{\alpha}}\textbf{v}_{i}(t-\delta
t/2)+\frac{m_{j}}{m_{\alpha}}\textbf{v}_{j}(t-\delta t/2),$ (15)

which determines the velocity $\textbf{v}_{\alpha}(t-\delta t/2)$ of the
merged object.

The algorithm for the system of baryonic objects consists of the equations
(A2)+(A3) for time steps without merging of objects, and the fusion of objects
is given by the equations (A5),(A6), (A7), (A8) and (A9). The discrete
algorithm with merging of objects has the same invariances as Newton’s exact
analytic and his discrete dynamics Toxvaerd2022 .

###  A.3 The algorithm for discrete classical mechanics in the expanding
universe

Newton’s discrete dynamics changes the position of an object $k$. If the space
expands monotonously over time with the Hubble velocity $\textbf{v}^{H}$, then
the expansion also changes the distance between two positions. The new
position $\textbf{r}_{k}(t+\delta t)$ is the sum of the change due to the
gravitational force on $k$ and the contribution from the Hubble expansion. The
mean location of an object changes from $\textbf{r}_{k}(t-\delta
t/2)=(\textbf{r}_{k}(t-\delta t)+\textbf{r}_{k}(t))/2$ at $t\in[t-\delta t,t]$
to $\textbf{r}_{k}(t+\delta t/2)=(\textbf{r}_{k}(t)+\textbf{r}_{k}(t+\delta
t))/2$ at $t\in[t,t+\delta t]$ . The Hubble expansion changes the distance
between the two positions by the Hubble velocity

$\displaystyle\textbf{v}^{H}=H\delta\textbf{r}_{k}=H(\textbf{r}_{k}(t+\delta
t/2)-\textbf{r}_{k}(t-\delta t/2))$ $\displaystyle=\delta
tH\frac{\textbf{r}_{k}(t+\delta t)-\textbf{r}_{k}(t)}{2\delta t}+\delta
tH\frac{\textbf{r}_{k}(t)-\textbf{r}_{k}(t-\delta t)}{2\delta t}$
$\displaystyle=\frac{\delta tH}{2}\textbf{v}_{k}(t+\delta t/2)+\frac{\delta
tH}{2}\textbf{v}_{k}(t-\delta t/2).$ (16)

By including the Hubble velocity, Eq. (A10), in the Newtons algorithm, Eq.
(A2), and after a re-arrangement, one obtains the algorithm for the classical
mechanics with a Hubble expansion of the space included in the Newtonian
dynamics

$\textbf{v}_{k}(t+\delta t/2)=\frac{(1+\delta tH/2)\textbf{v}_{k}(t-\delta
t/2)+\delta t/m_{k}\textbf{f}_{k}(t)}{1-\delta tH/2}.$ (17)

The discrete classical dynamics with Hubble expansion, Eqn. (A11)+(A3), is
still time reversible, but it increases the velocities, the momenta and the
angular momenta.

The small galaxies in the articles are obtained for thousand objects, which at
the start of the simulation, at $t=0$ are separated with a mean distance
$<r_{ij}(0)>\approx 1000$ and with a Maxwell-Boltzmann distributed velocities
with mean velocity $<\mid\textbf{v}_{i}(0)\mid>\approx 1$. The gravitational
strengths in the article are in units of the gravitational constant $G=1$ and
the mass unit $m=m_{i}(0)=1$ and diameters of the objects $\sigma_{i}(0)=1$.
For the set-up of the systems see also Toxvaerd2022 . The systems are followed
more than $6.4\times 10^{9}$ MD time steps, i.e. $t=1.6\times 10^{7}$ time-
units, which corresponds to 25 Gyr of our universe.

Acknowledgements This work was supported by the VILLUM Foundation Matter
project, grant No. 16515. Isaac Newton is gratefully acknowledged.

Data Availability Statement Data will be available on request.

References

## References

  * (1) Kolb, E. W., & Turner, M. S. 2018, The early Universe(CRS Press, 2018)
  * (2) de Martino, I., Chakrabarty, S. S., Gallo, A., Ostorero, L.,& Diaferio, A 2020, Universe 6, 107
  * (3) Hubble, E. 1929, Proc. Nath. Acad. Sci. 15, 168
  * (4) Klypin, A. A.,& Shandarin, S. F. 1983, MNRAS 204, 891
  * (5) Centrella, J., & Melott, A. L. 1983, Nature 305, 196
  * (6) Hockney, R. W., Goel, S. P.,& Eastwood, J.W. 1974, J. Comput. Phys. 14 148
  * (7) Zeldovich, Ya. B. 1970, Astrom. & Astrophys. 5, 84
  * (8) Diemand, J., Kuhlen, M., Madau, P., Zemp, M., Moore, B., Potter, D., & Stadel J. 2008, Nature 454, 735
  * (9) Alimi, J.-M. et al., 2008, SC2012: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE Computer Soc. Press, 2012, CA p. 1
  * (10) Schaye, J., et al. 2010, MNRAS 402 1536
  * (11) Dubois, Y., Volonteri, M., & Silk, J. 2014, MNRAS 440 1590
  * (12) Vogelsberger, M., et al. 2014, MNRAS 444 1518
  * (13) Dubois, Y., et al. 2016, MNRAS 463 3948
  * (14) Ludlow, A. D., Fall, S. M., Schaye, J., & Obreschkow, D. 2021, MNRAS 508 5114
  * (15) Vogelsberger, M., et al. 2014, NATURE REVIEWS PHYSICS 2 42-66
  * (16) Milgrom, M., 1983, Ap. J. 270, 371
  * (17) Nipoti, C., Londrillo, P., & Ciotti, L. 2007, MNRAS 381 L104
  * (18) Llinares, C. 2018, INTERNATIONAL JOURNAL OF MODERN PHYSICS D 27 (15)
  * (19) Kroupa, P. 2015 Can. J. Phys. 93  1
  * (20) Banik, I, & Zhao, H. 2022 Symmetry 14  1331
  * (21) Ellis, G. F. R., & Gibbons, G. W. 2014, Clas. Quantum Grav. 31, 025003; $ibid$ 32055001
  * (22) Ellis, G. F. R., & Gibbons, G. W. 2015, Clas. Quantum Grav. 32055001
  * (23) Wu, X., Famaey, B., Gentile, G., Combes, O. T., Angus, G. W., & Robin, A. C. 2007 Apj 665 L101
  * (24) Derakhshani, K. & Haghi, H. Apj 2014 78348; $ibid.$ 2014 785:166
  * (25) Angus, G. W., Gentile, G., Diaferio, A., Famaey, B., & van der Heyden, K. J. 2014, MNRAS 440 746
  * (26) Toxvaerd, S. 2022, Eur. J. Phys. 137:99
  * (27) Zhang, T-J., Ma, C., & Lan T. 2010, Adv. Astron. 2010, 184284.
  * (28) Carrol, B. W., & Ostlie, D. A. 1996, An Introduction to Modern Astrophysics(Addiso-Wesley Publishing Company, Reading 1996)
  * (29) Newton, I 1687, PHILOSOPHIÆ NATURALIS PRINCIPIA MATHEMATICA. LONDINI, Anno MDCLXXXVII. Second Ed.1713; Third Ed. 1726.
  * (30) Riess, A. G., Filippenko, A. V., & Challis, P. et al. 1998, Astronomical J. 116, 1009
  * (31) Perlmutter, S, Aldering, A., Goldhaber, A et al. 1999, Astrophys. J. 517, 565
  * (32) Soltis, J., Casertano, S., & Riess, A. G. 2021, Astrophys. J. Lett. 908: L5
  * (33) Rix, H-W., & Bovy 2013, J. Astron. Astrophys. Rev. 21:61
  * (34) Daeson, A. J., Fattahi, A., Frenk, C. S., Grand, R. J. J., Oman, K. A., Garrison-Kimmel, S., Simpson, C. M., & Navarro, J. F. 2020, MNRAS 496, 3929
  * (35) Li, Z-Z,& J. Han, J. 2021, Astrophys. J. Lett. 915 L18
  * (36) Camarillo, T., Dredger, P., & Ratra, B 2018, Astrophys. Space Sci. 363:268
  * (37) Tremaine, S,& Weinberg, M. D. 1984, MNRAS 209 729
  * (38) Chiba, R., Friske, J. K., & Schönrich, R. 2014, MNRAS 783 130
  * (39) Bond, H. E., Nelan E. P., VandenBerg, D. A., Schaefer, G. H., & Harmer, D. 2013, Astrophys. J. Lett. 765:L12
  * (40) Gupta, A., Mathur, S., Krongold, Y., Nicastro, F.,& Galeazzi 2012, Astrophys. J. Lett. 756:L8
  * (41) Bergma, J. N., Anderson, M. E., Miller, M. J., Hodges-Kluck, E., Dai, X., Li, J-T., Li, Y., & Qu, Z. 2018, Astrophys. J. 862:3
  * (42) Reid, M. J., Menten, K. M., Brunthaler, A. et al. 2014, Astrophys. J. 783:130
  * (43) Batygin, K., & Laughlin, G. 2008, Astrophys. J. 683 1207
  * (44) Mogavero, F. 2021, & Laskar J. 2021, A&A 655 A1
  * (45) Toxvaerd, S. 2014, J. Chem. Phys. 140, 044102
  * (46) Toxvaerd, S. 2020, Eur. J. Phys. 135:267
  * (47) Verlet, L 1967, Phys. Rev. 159, 98
  * (48) Lee, T. D. 1983, Phys. Lett. B 125B, 217
  * (49) Toxvaerd S., Phys. Rev. E 50, 2271 (1994).


# Realizing a Fake Projective Plane as a Degree 25 Surface in
${\mathbb{P}}^{5}$

Lev Borisov  Hill Center Department of Mathematics, Rutgers University, NJ
08854<EMAIL_ADDRESS>and  Zachary Lihn  Department of Mathematics,
Columbia University, New York, NY 10027, USA<EMAIL_ADDRESS>
###### Abstract.

Fake projective planes are smooth complex surfaces of general type with Betti
numbers equal to that of the usual projective plane. Recent explicit
constructions of fake projective planes embed them via their bicanonical
embedding in $\mathbb{P}^{9}$. In this paper, we study Keum’s fake projective
plane $(a=7,p=2,\\{7\\},D_{3}2_{7})$ and use the equations of [1] to construct
an embedding of fake projective plane in $\mathbb{P}^{5}$. We also simplify
the 84 cubic equations defining the fake projective plane in $\mathbb{P}^{9}$.

##  1\. Introduction

The Enriques–Kodaira classification splits compact complex surfaces $S$ into
10 classes based largely on their Kodaira dimension $k(S)$. While surfaces
with Kodaira dimension $<2$ are better understood, those of general type with
maximum Kodaira dimension $k(S)=2$ still need a detailed classification.

To each minimal model of a surface $S$ one associates a triple of numerical
invariants $(p_{g},q,K_{S}^{2})$, where $p_{g}=h^{0}(S,K_{S})$ is the
geometric genus, $q=h^{1}(S,{\mathcal{O}}_{S})$ is the irregularity, and
$K_{S}^{2}$ is the self-intersection number of the canonical class $K_{S}$.
These determine all the other classical invariants such as the topological
Euler characteristic $e_{top}(S)=12\chi({\mathcal{O}}_{S})-K^{2}_{S}$ and the
plurigenera $P_{m}(S)=h^{0}(S,mK_{S})$ [6]. It turns out that producing
surfaces with low $p_{g}$ and $q$ is quite difficult and a complete
classification appears far away [2]. In the case of $p_{g}=q=0$, one has the
Bogomolov-Miyaoka-Yau inequality $K^{2}_{S}\leq 9$. The focus of this paper is
the extreme case of surfaces with $p_{g}=q=0$ and $K^{2}_{S}=9$. These are the
fake projective planes (often called FPPs for short) which by definition are
complex projective surfaces of general type with Hodge diamond

$\begin{array}[]{ccccc}&&1&&\\\ &0&&0&\\\ 0&&1&&0\\\ &0&&0&\\\
&&1&&\end{array}$

which is the same as that of ${\mathbb{C}\mathbb{P}}^{2}$. The existence of a
fake projective plane was first proved by Mumford [13] by expressing the
surface as a quotient of a 2-adic analog of the complex two-dimensional ball

$\mathcal{B}^{2}=\\{(z_{1},z_{2})\in{\mathbb{C}}^{2}:|z_{1}|^{2}+|z_{2}|^{2}\leq
1\\}$

by a finitely generated group.

The general theory ensures that each fake projective plane is algebraic. By
Noether’s formula we know that $c_{1}^{2}=9$ and so all FPPs have
$c_{1}^{2}=3c_{2}=9$, where $c_{1},c_{2}$ are the Chern numbers. This implies
that each FPP is a quotient of $\mathcal{B}^{2}$ by an infinite discrete group
[15]. These ball quotients are determined by their fundamental group up to
holomorphic or anti-holomorphic isomorphism [12] and come in complex conjugate
pairs [8]. Each of the groups are arithmetic [9] and come in a finite list of
classes [14].

Based on the work of Prasad and Yeung [14], a complete classification was
obtained by Cartwright and Steger [3]. All fake projective planes are
quotients of $\mathcal{B}^{2}$ by explicit co-compact torsion-free arithmetic
subgroups of $\text{PU}(2,1)$. The classification was accomplished with
significant use of computer calculations. There are 50 conjugate pairs of fake
projective planes split among 28 classes. Each FPP is a ball quotient
$\mathcal{B}^{2}/\Gamma$ where $\Gamma$ is the fundamental group, and where
the automorphism group is $N(\Gamma)/\Gamma$ with $N(\Gamma)$ the normalizer
of $\Gamma$ in $\text{PU}(2,1)$. The torsion of the Picard group of
${\mathbb{P}}^{2}_{fake}$ is equal to the abelianization of $\Gamma$. Various
cover relations between related surfaces are also known [3].

###  1.1. The Geometry of Keum’s Fake Projective Plane

In this paper, we will focus on the fake projective plane
$(a=7,p=2,\\{7\\},D_{3}2_{7})$ in Cartwright-Steger classification. First
constructed in [7], it is named Keum’s fake projective plane and we will
denote it by ${\mathbb{P}}^{2}_{Keum}$. Its automorphism group has maximum
order among all FPPs, being equal to the semi-direct product of a normal
cyclic subgroup $C_{7}$ of order 7 and a non-normal cyclic subgroup $C_{3}$ of
order 3. By the Cartwright-Steger classification, there are three other fake
projective planes in its class including Mumford’s first fake projective
plane.

For the rest of this paper, we will let $K$ denote the canonical class of
Keum’s fake projective plane. The minimal resolution $Y$ of the quotient
${\mathbb{P}}^{2}_{Keum}/C_{7}$ by the subgroup $C_{7}$ of its automorphism
group has interesting geometry which we describe briefly.

Recall that a singular point of type $\frac{1}{m}(1,a)$ is a cyclic quotient
singularity given local analytically by the action $(x,y)\mapsto(\zeta
x,\zeta^{a}y)$ on ${\mathbb{C}}^{2}$ for $\zeta$ a primitive $m$th root of
unity. $Y$ has three singular points of type $\frac{1}{7}(1,3)$ permuted by
the residual $C_{3}$ automorphism group of ${\mathbb{P}}^{2}_{Keum}$. It is
also a Dolgachev surface fibered over ${\mathbb{P}}^{1}$, with generic fibers
of genus 1, two multiple fibers, three nodal fibers, and one fiber of type
$I_{9}$. The two multiple fibers are $2F_{3}$ and $3F_{2}$, which have
multiplicity 2 and 3 respectively. The reductions $F_{3}$ and $F_{2}$ are
linearly equivalent to $3K_{Y}$ and $2K_{Y}$. We refer to [7, 1] for more
details.

###  1.2. Explicit Construction of ${\mathbb{P}}^{2}_{Keum}$

In [1], Keum’s fake projective plane was explicitly constructed via its
bicanonical embedding as the vanishing set of 84 cubic equations in
${\mathbb{P}}^{9}$. One first constructs a birational model $Y_{0}$ of $Y$ as
a system of quadrics in 8 variables defined over ${\mathbb{Q}}(\sqrt{-7})$.
Included is a construction of the double and triple fibers and the $C_{3}$
action on $Y_{0}$. A degree 7 extension of the field of rational functions of
$Y_{0}$ gives the sevenfold cover of $Y_{0}$, which is exactly
${\mathbb{P}}^{2}_{Keum}$. Ten sections of ${\mathcal{O}}(2K)$ are constructed
from this description and the embedding in ${\mathbb{P}}^{9}$ is finally given
by 84 cubic equations in the 10 variables $P_{0},\ldots,P_{9}$.

A perennial question is how to simplify the equations of a fake projective
plane, which can have polynomials with coefficients hundreds to thousands of
decimal digits long. In this paper, we give a simplified version of the
equations of Keum’s fake projective plane in [1]. We use the equations to find
an embedding of ${\mathbb{P}}^{2}_{Keum}$ as a degree 25 surface in
${\mathbb{P}}^{5}$. The embedding is given by sections of ${\mathcal{O}}(5H)$,
where $H$ is a divisor such that $3H$ is linearly equivalent to $K$. Finally,
we exhibit the surface as a system of 56 sextics in ${\mathbb{P}}^{5}$ with
coefficients in the field ${\mathbb{Q}}(\sqrt{-7})$.

The paper is organized as follows. In Section 2 we outline the steps to
simplify the 84 cubics defining ${\mathbb{P}}^{2}_{Keum}$ in
${\mathbb{P}}^{9}$. We follow the strategy described in [1] by explicitly
calculating the nonreduced linear cuts on ${\mathbb{P}}^{2}_{Keum}$
corresponding to $2$-torsion in the Picard group. Using these equations, in
Section 3 we describe the steps to embed ${\mathbb{P}}^{2}_{Keum}$ in
${\mathbb{P}}^{5}$. Specifically, we compute global sections of
${\mathcal{O}}(5H)$ as global sections of the divisor $18H-9H-4H$ and explain
the key idea that allowed us to find
$H^{0}\left({\mathbb{P}}^{2}_{fake},{\mathcal{O}}(4H)\right)$. Section 4
concludes with future directions.

######  Remark 1.1.

A defining feature of recent constructions of fake projective planes is their
heavy use of computer algebra software. To that end, this project depended
heavily on the use of the Mathematica software system [10] and the computer
algebra systems Magma [11] and Macaulay2 [4].

######  Remark 1.2.

The 84 cubics in ${\mathbb{P}}^{9}$ and the 56 sextics in ${\mathbb{P}}^{5}$
are still too large to be included in the printed paper.

##  2\. Simplification of Keum’s Fake Projective Plane

We will begin by simplifying the explicit equations of Keum’s fake projective
plane ${\mathbb{P}}^{2}_{Keum}$ found in [1]. This is done by looking for
nonreduced cuves on ${\mathbb{P}}^{2}_{Keum}$ which correspond to $2$-torsion
in the Picard group. We proceed by making a coordinate change that makes the
curve nicer in our new basis.

### Step 1: Finite Field Search for Nonreduced Curves

By the Cartwright-Steger classification, the torsion in the Picard group of
${\mathbb{P}}^{2}_{Keum}$ is $C_{2}^{3}$. In addition, the automorphism group
is $C_{7}\rtimes C_{3}$, the semidirect product of $C_{7}$ and $C_{3}$.

We claim that 2-torsion classes give nonreduced curves in $|2K|$. Let $L$ be a
2-torsion class in the Picard group. By [5], we have
$h^{0}({\mathbb{P}}^{2}_{Keum},K+L)=1$. Hence, up to scaling, there is a
unique section $s_{L}\in H^{0}({\mathbb{P}}^{2}_{Keum},K+L$). The square of
$s_{L}$ is in $H^{0}({\mathbb{P}}^{2}_{Keum},2K)$ and gives rise to a
nonreduced curve.

We will further assume that the nonreduced curve is $C_{3}$ invariant. This
reduces the search to nonreduced curves of the form

$a_{0}P_{0}+a_{1}(P_{1}+P_{2}+P_{3})+a_{2}(P_{4}+P_{5}+P_{6})+a_{3}(P_{7}+P_{8}+P_{9})$

up to scaling (so we subsequently set $a_{0}=1)$. To look for such curves we
look at a finite field reduction of ${\mathbb{P}}^{2}_{Keum}$ over
${\mathbb{F}}_{p}$ for suitable $p$. More precisely, such suitable $p$
contains a square root of $-7$ and has the same Hilbert polynomial for
${\mathbb{P}}^{2}_{Keum}$ modulo $p$. We picked $p=43$ with $\sqrt{-7}\equiv
6\mod{43}$ which was an arbitrary small prime with the aforementioned
conditions. Using Magma, we ran an exhaustive search for all
$a_{1},a_{2},a_{3}$ in ${\mathbb{F}}_{43}$ and checked if the corresponding
curve is nonreduced. We obtained the curve

$P_{0}+24(P_{1}+P_{2}+P_{3})+0(P_{4}+P_{5}+P_{6})+28(P_{7}+P_{8}+P_{9}).$

### Step 2: Lift to Characteristic 0

We lift this curve to ${\mathbb{Q}}(\sqrt{-7})$ as follows. Using Magma, we
calculate some points in ${\mathbb{F}}_{43}$ lying on
${\mathbb{P}}^{2}_{Keum}$ and the nonreduced curve. We then apply a variant of
Hensel lifting to lift the curve to ${\mathbb{Z}}/43^{k}{\mathbb{Z}}$ for
higher $k$ at each step, obtaining a $p$-adic approximation.

The lifting process was done by finding, at each point, two linearly-
independent tangent vectors in ${\mathbb{P}}^{9}({\mathbb{F}}_{43})$ that are
orthogonal to all polynomials defining ${\mathbb{P}}^{2}_{Keum}$ and the
linear cut. We modified the points, tangent vectors, and the linear cut at
each stage to lift them to higher powers of $43$ such that the orthogonality
conditions held; this reduced to solving a system of linear equations modulo
$43$. After a sufficiently high power of $43$ we identify the corresponding
algebraic numbers by applying a lattice reduction algorithm. We obtain the
curve

$\displaystyle
P_{0}+\frac{(-1+\sqrt{-7})}{2}(P_{1}+P_{2}+P_{3})+\frac{(272-848\sqrt{-7})}{7}(P_{4}+P_{5}+P_{6})$
$\displaystyle+\frac{(832-192\sqrt{-7})}{7}(P_{7}+P_{8}+P_{9})$

which we verify is nonreduced numerically.

Thus we have found one nontrivial $C_{3}$-invariant torsion line bundle. It is
not $C_{7}$-invariant because the corresponding nonreduced linear cut is not
$C_{7}$ invariant. Its orbit therefore has $7$ elements, which combined with
knowledge of the torsion of the Picard group as $C_{2}^{3}$ shows that the
action of the automorphism group on the torsion in Picard group is transitive.

### Step 3: Setting Up the Coordinate Change

Finally we set up the coordinate change to find a nicer basis for
$H^{0}({\mathbb{P}}^{2}_{Keum},2K)$ in order to simplify the equations
defining the fake projective plane. We use a coordinate change from $P_{i}$ to
$Q_{i}$ that respects the automorphisms on the surface such that the the
nonreduced cut becomes

$Q_{0}+Q_{1}+Q_{2}+Q_{3}+Q_{7}+Q_{8}+Q_{9}.$

These conditions leave one free parameter in the coordinate change. We fix the
free parameter by choosing it in such a way to set the ”simplest” coefficient
in the equations to 1. This allows us to find a version of the 84 equations
with significantly smaller coefficients.

We simplify the equations further by reducing the number of monomials in the
equations. We take random linear combinations of the seven equations in each
$C_{7}$ weight and select those that span the space and have the fewest
monomials.

##  3\. Embedding of a fake projective plane into ${\mathbb{P}}^{5}$

In this section, we will describe the process that led us to find the
equations of an embedding of ${\mathbb{P}}^{2}_{Keum}$ in ${\mathbb{P}}^{5}$.
Let $H$ be a divisor such that $3H=K$, where $K=K_{{\mathbb{P}}^{2}_{Keum}}$
is the canonical divisor of ${\mathbb{P}}^{2}_{Keum}$. Calculations of
$h^{0}({\mathbb{P}}^{2}_{Keum},nH)$ show that $|5H|$ has the expected
dimension such that the corresponding map to projective space is
${\mathbb{P}}^{5}$. Thus we aim to construct $|5H|$ explicitly, which will
give the desired map ${\mathbb{P}}^{2}_{Keum}\to{\mathbb{P}}^{5}$.

$n$ | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12
---|---|---|---|---|---|---|---|---|---|---
$h^{0}\left({\mathbb{P}}^{2}_{Keum},nH\right)$ | 0 | 3 | 6 | 10 | 15 | 21 | 28 | 36 | 45 | 55

Table 1. Dimensions of $H^{0}({\mathbb{P}}^{2}_{Keum},nH)$ for different
values of $n$, where $3H=K$

.

Recall that $Y$ denotes the quotient ${\mathbb{P}}^{2}_{Keum}/C_{7}$ of Keum’s
fake projective plane by its $C_{7}$ automorphism subgroup. It has residual
automorphism group $C_{3}$ and has a double fiber $F=3K_{Y}$.

We will construct $|5H|$ as the space $|18H-9H-4H|$. We first find
$|9H|=|18H-9H|$ by expressing 112 cubic equations in the $Q_{i}$ (which lie in
$18H=6K$) that vanish on $9H$. Crucial to this construction is the preimage of
the double fiber $F$ of $Y$ which we use to find points on $9H$. We then
compute $|4H|$. This required the use of several important ideas which are
detailed in Step 2 below. Finally, after constructing $4H$ we may find $5H$ as
linear combinations of the equations of $9H$ vanishing on $4H$. We conclude by
using the explicit equations in ${\mathbb{P}}^{5}$ to reconstruct the $C_{3}$
action on ${\mathbb{P}}^{2}_{Keum}$ in its embedding into ${\mathbb{P}}^{5}$.

### Step 1: Constructing $|9H|$

The preimage of the double fiber on $Y$ has divisor class $3K=9H$ [1]. Hence
to construct $|9H|$ we are led to find polynomials on $Y$ vanishing on the
double fiber. Recall that [1] constructs the surface $Y$ as a system of
quadrics in the variables $u_{0},u_{1},w_{1},\ldots,w_{6}$ with the double
fiber given by $\\{u_{1}=0\\}$. We compute a number of random points on the
double fiber of $Y$ and use the equations to construct points on
${\mathbb{P}}^{2}_{Keum}$ lying on the preimage of the double fiber. We then
look for polynomials vanishing on these points to compute
$H^{0}({\mathbb{P}}^{2}_{Keum},9H)$. The search for cubic polynomials gave 112
cubics with 16 in each $C_{7}$ weight.

### Step 2: Constructing $|4H|$

We may attempt to construct $|4H|$ as follows. The action of the $C_{7}$
automorphism subgroup on $H^{0}({\mathbb{P}}^{2}_{Keum},4H)$ gives a
$C_{7}$-representation which splits $H^{0}({\mathbb{P}}^{2}_{Keum},4H)$ into
three one-dimensional $C_{7}$-eigenspaces. The Holomorphic Lefschetz Fixed-
Point formula shows that the eigenvalues are $\xi^{3},\xi^{5},\xi^{6}$, where
$\xi$ is a seventh root of unity. Thus
$H^{0}({\mathbb{P}}^{2}_{Keum},4H)\cong{\mathbb{C}}r_{3}\,\oplus\,{\mathbb{C}}r_{5}\,\oplus\,{\mathbb{C}}r_{6}$
where $r_{3},r_{5},r_{6}$ are sections of $4H$ with $C_{7}$-weights 3, 5, and
6 respectively. In addition, the $C_{3}$-action on the surface implies
$r_{5}=\sigma(r_{3}),r_{6}=\sigma^{2}(r_{3})$ for $\sigma$ an order 3
automorphism on ${\mathbb{P}}^{2}_{Keum}$. The product $d=r_{3}r_{5}r_{6}$ is
therefore a $C_{3}$-invariant section with $C_{7}$-weight 0 in
$H^{0}({\mathbb{P}}^{2}_{Keum},12H)$ (it is then invariant under the whole
automorphism group).

Set $s_{i}=r_{i}^{3}\in H^{0}({\mathbb{P}}^{2}_{Keum},12H)$ for
$i\in\\{3,5,6\\}$. The equation

$s_{3}s_{5}s_{6}=d^{3}$

in $H^{0}({\mathbb{P}}^{2}_{Keum},36H)$ allows us to narrow down parameters in
the search for $r_{3},r_{5},r_{6}$. Since $s_{3},s_{5},s_{6},$ and $d$ lie in
$H^{0}({\mathbb{P}}^{2}_{Keum},12H)$, they are quadratic in the variables
$Q_{0},\ldots,Q_{9}$ for the fake projective plane. It is sufficient to
construct $s_{3}$ since $s_{5}$ and $s_{6}$ may be constructed from $s_{3}$
with the $C_{3}$ action. Additionally, since $s_{3}$ has $C_{7}$ weight
$3\times 3\equiv 2\mod{7}$, we narrow the search down to $C_{7}$-weight 2
quadratics.

We may further reduce the number of parameters with additional data. The curve
$\\{r_{3}=0\\}$ passes through the two $C_{7}$ fixed points

$\displaystyle p_{1}$ $\displaystyle=(0\colon 0\colon 0\colon 0\colon 0\colon
0\colon 0\colon 1\colon 0\colon 0),$ $\displaystyle p_{2}$
$\displaystyle=(0\colon 0\colon 0\colon 0\colon 0\colon 0\colon 0\colon
0\colon 1\colon 0).$

It follows that at these points the curve $\\{s_{3}=0\\}$ vanishes with
multiplicity 3, which place additional conditions on the coefficients of
$s_{3}$.

Now we begin describing the details of the calculation. We first calculate the
order 3 neighborhoods of the points $p_{1}$ and $p_{2}$. This was done by
computing the tangent space and solving for the conditions of the
neighborhoods vanishing on the FPP. We began by solving for the order 2
neighborhood and then for the third order. To speed up calculations, it was
sufficient to take some equations for ${\mathbb{P}}^{2}_{Keum}$ locally
cutting out the point. After computing these neighborhoods, we posit the
general form for $s_{3}$ as weight 2 quadratics in the variables and then
solve for the conditions of being identically 0 at the higher order
neighborhoods. We are able to solve for two of these variables, narrowing down
the general form for $s_{3}$ to 6 variables.

We now want to solve for the sextic equation $s_{3}s_{5}s_{6}-d^{3}=0$. The
requirement that $d$ be invariant under the full automorphism group forces it
to be of the form

$e_{1}Q_{0}^{2}+e_{2}(Q_{1}Q_{6}+Q_{2}Q_{4}+Q_{3}Q_{5})+e_{3}(Q_{1}Q_{9}+Q_{2}Q_{7}+Q_{3}Q_{8})$

for undetermined coefficients $e_{1},e_{2},e_{3}$. We also obtain the general
forms for $s_{5}$ and $s_{6}$ by applying the $C_{3}$ automorphism to $s_{3}$.
To solve for the coefficients, we compute some points of
${\mathbb{P}}^{2}_{Keum}$ with high accuracy and substitute them into
$s_{3}s_{5}s_{6}-d^{3}=0$ to obtain a system of 24 cubics in 6 variables. We
solve this system of equations by applying the trick of [2]. The Hilbert
polynomial of the system of equations modulo 37 with $\sqrt{-7}\equiv
17\mod{37}$ is 3, which suggests that there are 3 solutions for this system.
By applying successive linear conditions on the system and checking the
Hilbert polynomial at each step, we are able to take linear cuts that drop the
Hilbert polynomial eventually to 1. At some point there are 3 different
choices for the linear cuts corresponding to our 3 solutions. We were able to
lift these 3 solutions modulo $37^{200}$ and then use the lattice reduction
algorithm to obtain the corresponding solutions over
${\mathbb{Q}}(\sqrt{-7})$. The three solutions differed by a cube root of
unity. We selected the solution defined over the desired field of definition
to proceed.

The solution for these coefficients allow us to fully determine
$s_{3},s_{5},s_{6}$, and $d$. The equations for $s_{3}$ and $d$ are given
below, with $s_{5}$ and $s_{6}$ found by applying the $C_{3}$ automorphism.
Points on $\\{r_{3}=0\\}$ may then be calculated by solving for the
simultaneous conditions $\\{s_{3}=0,d=0\\}$. These points were used later in
the construction.

$\displaystyle s_{3}$
$\displaystyle=\frac{\left(-212275+26525i\sqrt{7}\right)Q_{0}Q_{5}}{2470336}+\frac{\left(22575+51275i\sqrt{7}\right)Q_{0}Q_{8}}{1235168}$
$\displaystyle+\frac{\left(139475+17575i\sqrt{7}\right)Q_{1}Q_{2}}{9881344}+\frac{\left(196875-91425i\sqrt{7}\right)Q_{3}Q_{4}}{2470336}$
$\displaystyle+\frac{\left(-303625-270725i\sqrt{7}\right)Q_{3}Q_{7}}{4940672}+\frac{\left(139475+17575i\sqrt{7}\right)Q_{6}^{2}}{1235168}$
$\displaystyle+\frac{\left(795725-287175i\sqrt{7}\right)Q_{6}Q_{9}}{4940672}+\frac{\left(-57575-549675i\sqrt{7}\right)Q_{9}^{2}}{9881344}$
$\displaystyle d$
$\displaystyle=\frac{25}{9881344}\Bigl{(}3407\sqrt{-7}Q_{0}^{2}+17045Q_{0}^{2}-2812\sqrt{-7}Q_{1}Q_{6}-22316Q_{1}Q_{6}$
$\displaystyle+329\sqrt{-7}Q_{1}Q_{9}-21987Q_{1}Q_{9}-2812\sqrt{-7}Q_{2}Q_{4}-22316Q_{2}Q_{4}+329\sqrt{-7}Q_{2}Q_{7}$
$\displaystyle-21987Q_{2}Q_{7}-2812\sqrt{-7}Q_{3}Q_{5}-22316Q_{3}Q_{5}+329\sqrt{-7}Q_{3}Q_{8}-21987Q_{3}Q_{8}\Bigr{)}$

### Step 3: The map ${\mathbb{P}}^{2}_{Keum}\to{\mathbb{P}}^{5}$

With the computations of $9H$ and $4H$ we may now find $5H$. We look at
suitable linear combinations of the 112 polynomials vanishing on $18H-9H=9H$
additionally vanishing on $4H$ to obtain $18H-9H-4H=5H$.

We first compute some random points on $4H$ by solving for $\\{d=0,s_{3}=0\\}$
on the FPP. $5H$ is then found by looking for linear combinations of the
cubics defining $9H$ for each weight that vanish on these points. To verify
that they are in $|5H|$ we also check that they do not vanish on the whole
fake projective plane.

The six resulting degree 3 polynomials give us the map
${\mathbb{P}}^{2}_{Keum}\to{\mathbb{P}}^{5}$. We calculate points in the image
of this map in ${\mathbb{P}}^{5}$ and find 56 degree 6 polynomials in new
variables $Z_{1},\ldots,Z_{6}$ that vanish at these points. These give the
desired embedding of the fake projective plane.

######  Remark 3.1.

The $C_{7}$-weights on the variables $Z_{1},Z_{2},\ldots,Z_{6}$ are
$1,2,\ldots,6$. There is no weight 0 variable. The construction required that
we shift the $C_{7}$ weights by $3$. This may be explained by viewing our
construction of $H^{0}({\mathbb{P}}^{2}_{Keum},5H)$ as given by an embedding

$H^{0}({\mathbb{P}}^{2}_{Keum},5H)\hookrightarrow
H^{0}({\mathbb{P}}^{2}_{Keum},18H)$

with the map given by tensoring with $s_{3}\otimes f$ for $s_{3}\in
H^{0}({\mathbb{P}}^{2}_{Keum},4H)$ and $f\in
H^{0}({\mathbb{P}}^{2}_{Keum},9H)$. While $f$ has weight 0, $s_{3}$ has weight
$3$ and therefore shifts the weights of $H^{0}({\mathbb{P}}^{2}_{Keum},5H)$ by
3.

We take care to reconstruct the automorphism group. While the $C_{7}$-action
is preserved under our construction, the non-$C_{3}$-invariance of $s_{3}$
introduces a scaling factor in the $C_{3}$ action. We fix the coefficients of
this scaling factor and recompute the equations with the scaling to find a
better basis for the action. As before, we take random linear combinations of
the equations that span the space and take the simplest ones to further
simplify the equations.

Finally, we use Magma to verify that the Hilbert polynomial is as expected.
The verification process for the FPP is carried out as in [1] working modulo
$p=1327$ with $\sqrt{-7}=103\mod{1327}$. Thus we have constructed Keum’s fake
projective plane as a degree 25 surface in ${\mathbb{P}}^{5}$.

##  4\. Future Directions

One hopes to find a coordinate change to additionally simplify the 56
equations in ${\mathbb{P}}^{5}$.

A related construction of interest is that of Mumford’s original fake
projective plane [13]. This surface has not been explicitly constructed yet.
It lies in the same class as ${\mathbb{P}}^{2}_{Keum}$ and two other fake
projective planes. We are currently attempting to find this surface by
computing a seven-to-one cover of ${\mathbb{P}}^{2}_{Keum}$, after which
several cover relations may yield the surface and the two fake projective
planes in the same class.

##  5\. Acknowledgements

The authors thank the DIMACS REU program at Rutgers University for supporting
this research project. This work was carried out while the second author was
supported by NSF grant CCF-1852215.

## References

  * [1] L. Borisov, On equations of fake projective planes with automorphism group of order $21$, arXiv:2109.02070, submitted.
  * [2] L. Borisov and E. Fatighenti, New explicit constructions of surfaces of general type, arXiv preprint arXiv:2004.02637 (2020).
  * [3] D. Cartwright, T. Steger, Enumeration of the 50 fake projective planes, C. R. Acad. Sci. Paris, Ser. I 348 (2010), 11-13.
  * [4] D. Grayson and M. Stillman, Macaulay2, a software system for research in algebraic geometry, http://www.math.uiuc.edu/Macaulay2.
  * [5] S. Galkin, I. Karzhemanov, E. Shinder, Fake projective planes, automorphic forms, exceptional collections, arXiv:1602.06107.
  * [6] R. Hartshorne, Algebraic Geometry, Graduate Texts in Mathematics (1977).
  * [7] J. Keum, A fake projective plane with an order 7 automorphism, Topology 45 (2006) 919-927.
  * [8] V. S. Kharlamov, V. M. Kulikov, On real structures on rigid surfaces, Izv. Russ. Akad. Nauk. Ser. Mat. 66, no. 1, (2002) 133-152; Izv. Math. 66, no. 1, (2002) 133-150.
  * [9] B. Klingler, Sur la rigidit de certains groupes fondamentaux, l’arithmeticite des reseaux hyperboliques complexes, et les ”faux plans projectifs”, Invent. Math. 153 (2003), 105-143.
  * [10] Mathematica 11.1, Wolfram Research, Inc., 2016.
  * [11] W. Bosma, J. Cannon, and C. Playoust, The Magma algebra system. I. The user language, J. Symbolic Comput., 24 (1997), 235-265; v2.23-2.
  * [12] G. D. Mostow, Strong rigidity of locally symmetric spaces, Annals Math. Studies 78, Princeton Univ. Press, Princeton, N.J.; Univ. Tokyo Press, Tokyo 1973.
  * [13] D. Mumford, An algebraic surface with $K$ ample, $K2=9,p_{g}=q=0$, Amer. J. Math. 101 (1979) 233-244.
  * [14] G. Prasad, S.-K. Yeung, Fake projective planes, Invent. Math. 168 (2007), 321-370; Addendum, 182 (2010), 213-227; Addendum. Invent. Math. 168(2007), 321-370.
  * [15] S.-T. Yau, Calabi’s conjecture and some new results in algebraic geometry, Proc. Nat. Acad. Sci. U.S.A. 74 (1977), no. 5, 1798-1799.