# Dark Matter Deficient Galaxies And Their Member Star Clusters Form
Simultaneously During High-velocity Galaxy Collisions In 1.25 pc Resolution
Simulations

Joohyun Lee Center for Theoretical Physics, Department of Physics and
Astronomy, Seoul National University, Seoul 08826, Korea  Eun-jin Shin Center
for Theoretical Physics, Department of Physics and Astronomy, Seoul National
University, Seoul 08826, Korea  Ji-hoon Kim Center for Theoretical Physics,
Department of Physics and Astronomy, Seoul National University, Seoul 08826,
Korea  Seoul National University Astronomy Research Center, Seoul 08826, Korea

###### Abstract

How diffuse dwarf galaxies that are deficient in dark matter—such as
NGC1052-DF2 and NGC1052-DF4—formed is a mystery. Along with their luminous
member globular clusters (GCs), the so-called dark matter deficient galaxies
(DMDGs) have challenged observers and theorists alike. Here we report a suite
of galaxy collision simulations using the adaptive mesh refinement code Enzo
with 1.25 pc resolution, which demonstrates that high-velocity galaxy
collisions induce the formation of DMDGs and their star clusters (SCs)
simultaneously. With a numerical resolution that is significantly better than
our previous study, we resolve the dynamical structure of the produced DMDGs
and the detailed formation history of their SCs, which are possible
progenitors of the DMDG’s member GCs. In particular, we show that a galaxy
collision with a high relative velocity of $\sim 300\;{\rm km\;s}^{-1}$,
invoking severe shock compression, spawns multiple massive SCs
($M_{\star}\gtrsim 10^{6}\,M_{\odot}$) in $<150\;{\rm Myr}$ after the
collision. At the end of the $\sim 800\;{\rm Myr}$ evolution in our fiducial
run, the resulting DMDG of $M_{\star}\simeq 3.5\times 10^{8}\,M_{\odot}$ hosts
10 luminous ($M_{V}\lesssim-8.5$ mag), gravitationally bound SCs with a line-
of-sight velocity dispersion $11.2\;{\rm km\;s}^{-1}$. Our study suggests that
DMDGs and their luminous member SCs could form simultaneously in high-velocity
galaxy collisions while being in line with the key observed properties of
NGC1052-DF2 and NGC1052-DF4.

galaxies: formation — galaxies: evolution — galaxies: star formation —
galaxies: star clusters: general — cosmology: theory — methods: numerical

††software:  yt (Turk et al., 2011), Enzo (Bryan et al., 2014; Brummel-Smith
et al., 2019), the Grackle chemistry and cooling library (Smith et al., 2017),
numpy (Oliphant, 2006), scipy (Virtanen et al., 2020), matplotlib (Hunter,
2007)

##  1 Introduction

At the center of the intense debate on the formation mechanism that can
address abnormal lack of dark matter and unusually bright populations of
globular clusters (GCs) are two recently discovered ultra-diffuse dwarf
galaxies, NGC1052-DF2 and NGC1052-DF4 (hereafter DF2 and DF4, respectively).
The spectroscopic study of their member GCs has revealed surprisingly small
velocity dispersions compared to what is normally expected by their observed
stellar masses. This implies that DF2 and DF4 have very small dark-matter-to-
stellar-mass ratios (dark matter deficient galaxy (DMDG); van Dokkum et al.,
2018b, 2019).111 This anomalous characteristic is based on the distance
measurement, the accuracy of which is still under debate (see e.g., Trujillo
et al. 2019 and Monelli & Trujillo 2019 for details, but also van Dokkum et
al. 2018a, Danieli et al. 2020 and Shen et al. 2021b). The exceptionally large
population of their luminous member GCs is also a mystery. The GC luminosity
function (GCLF) of DF2 peaks at $M_{V}\sim-9$, approximately 1.5 mag shifted
to the brighter side from the Milky Way’s GCLF (van Dokkum et al., 2018c; Ma
et al., 2020a; Shen et al., 2021a), indicating a large number of massive GCs.
Extensive spectroscopic studies of DF2 have found that their stellar
population is $\sim$9 Gyr old on average and that the metallicity of the main
stellar body differs from that of its member GCs (Fensch et al., 2019; Ruiz-
Lara et al., 2019). As for the kinematics, recent studies hint that DF2 is
likely a prolate rotator (Emsellem et al., 2019; Lewis et al., 2020).

Table 1: A Suite of Idealized Galaxy Gollision simulations listed with their
initial configurations and the stellar mass of the most massive DMDG formed

Run Name | Stellar Feedback Energy | Initial Gas Metallicity | Disk Inclination Angle | Disk Spin Relative Angle | $M_{\star,\,\rm DMDG}$ | $t_{\rm end}$
---|---|---|---|---|---|---
| (1051 erg/150$\,M_{\odot}$ of Stars) | ($Z_{\odot}$) | ($\theta_{1}$,$\phi_{1}$), ($\theta_{2}$,$\phi_{2}$) |  | ($10^{8}\,M_{\odot}$) | (Myr)
(1) | (2) | (3) | (4) | (5) | (6) | (7)
Fiducial††footnotemark:  | 1 | 1 | (0∘, 0∘), (0∘, 0∘) | 0∘ | 3.47 (3.72) | 800 (250)
FB3 | 3.3 | 1 | (0∘, 0∘), (0∘, 0∘) | 0∘ | 2.76 | 250
FB7 | 6.7 | 1 | (0∘, 0∘), (0∘, 0∘) | 0∘ | 1.46 | 250
FB10 | 10 | 1 | (0∘, 0∘), (0∘, 0∘) | 0∘ | - | 50
Low-metal‡‡footnotemark:  | 1 | $10^{-2}$ | (0∘, 0∘), (0∘, 0∘) | 0∘ | 0.90 | 250
Retrograde | 1 | 1 | (0∘, 0∘), (180∘, 0∘) | 180∘ | 1.76 | 250
Tilt-1 | 1 | 1 | ($-22\hbox to0.0pt{.\hss}^{\circ}5$, 0∘), ($22\hbox to0.0pt{.\hss}^{\circ}5$, 0∘) | 45∘ | 2.33 | 250
Tilt-1-FB10 | 10 | 1 | ($-22\hbox to0.0pt{.\hss}^{\circ}5$, 0∘), ($22\hbox to0.0pt{.\hss}^{\circ}5$, 0∘) | 45∘ | - | 50
Tilt-2 | 1 | 1 | (0∘, $-22\hbox to0.0pt{.\hss}^{\circ}5$), (0∘, $22\hbox to0.0pt{.\hss}^{\circ}5$) | 45∘ | - | 50
Tilt-2-FB10 | 10 | 1 | (0∘, $-22\hbox to0.0pt{.\hss}^{\circ}5$), (0∘, $22\hbox to0.0pt{.\hss}^{\circ}5$) | 45∘ | - | 50

  * •

† This is the fiducial run presented in all of the figures in this article.
See Sections 2 and 3 for more information.

  * •

‡ This is the run presented in the right panels of Figure 5 in Section 3.2.

Note. — The simulation configurations and the stellar mass of the most massive
DMDG formed: (1) run name, (2) stellar feedback energy, (3) initial
metallicity of the gas disks, (4) angle by which each of the two colliding
galactic disks is rotated about the +$x$-axis ($\theta_{1}$, $\theta_{2}$) and
about the +$y$-axis ($\phi_{1}$, $\phi_{2}$), starting from its original
position in the $x$-$y$ plane (see the top-left panel of Figure 1), (5)
relative angle between the angular momentum vectors of the two galactic disks,
(6) stellar mass within 12 kpc from the center of the most massive DMDG (“-”
$=$ no DMDG formed), (7) time since the pericentric approach of the two
progenitor disks when we end the simulation.

These observations could provide a testbed for various formation scenarios
proposed for DF2 and DF4. In one of such scenario suggested by Silk (2019), a
“mini-Bullet cluster-like” event on a galaxy scale decouples dark matter from
dissipative baryons.222 Other scenarios include the tidal dwarf galaxies
(Fensch et al., 2019) and the dwarf galaxies with tidally stripped dark matter
halos (Ogiya, 2018; Jackson et al., 2021). See Shin et al. (2020) for more
information. The high-velocity galaxy collision induces a strong shear and
shock compression in the interstellar medium (ISM), instigating the formation
of DMDGs and their bright member GCs (Silk, 2019; Trujillo-Gomez et al.,
2020). Numerical studies have demonstrated the formation of young massive
clusters (YMCs) and/or GC candidates triggered by galaxy mergers in idealized
simulations (Li et al., 2004; Bournaud et al., 2008; Kruijssen et al., 2011,
2012; Hopkins et al., 2013; Renaud et al., 2015; Maji et al., 2017; Moreno et
al., 2019; Lahén et al., 2019, 2020) and in cosmological simulations (Li et
al., 2017; Kim et al., 2018; Ma et al., 2020b). To validate the scenario by
Silk (2019), we previously simulated idealized galaxy collisions using two
different gravito-hydrodynamics codes and showed that
$\sim$10${}^{8}\,M_{\odot}$ DMDGs are produced when two gas-rich, dwarf
galaxies collide with a relative velocity of $\sim 300\;{\rm km\;s}^{-1}$
(Shin et al., 2020). However, the chosen resolution of 80 pc limited our
ability to examine the inner structure of the resulting DMDGs, or to resolve
the formation of their member SCs.

In this Letter, we examine the galaxy collision–induced DMDG formation
scenario using adaptive mesh refinement (AMR) simulations with significantly
improved 1.25 pc resolution. We not only demonstrate the DMDGs’ formation but
also explore the detailed formation histories and the dynamical properties of
the DMDG and their member SCs which could eventually evolve into bright GCs.

##  2 Simulations

The initial condition used for this study is similar to what was adopted in
Shin et al. (2020) but with greatly improved particle resolution: a dark
matter halo of mass $3.88\times 10^{9}\,M_{\odot}$ with a particle mass
$1.55\times 10^{3}\,M_{\odot}$ each, a stellar disk of $4.21\times
10^{8}\,M_{\odot}$ with a particle mass $1.05\times 10^{3}\,M_{\odot}$ each,
and a gas disk of $1.68\times 10^{9}\,M_{\odot}$. The initial gas fraction is
thus $f_{\rm gas}=M_{\rm gas}/M_{\rm galaxy}=0.28$, but because $\sim 30\%$ of
the gas is consumed before the collision, $f_{\rm gas}$ becomes $\sim 0.2$
when two galaxies collide at $t=0$. Our Fiducial run is a co-planar
prograde–prograde collision of the two identical galaxies with a relative
velocity of $\sim 300\;{\rm km\;s}^{-1}$ at a 40 kpc distance, and with a
pericentric distance of 1 kpc (see Table 1).333 Readers should note that this
initial condition models a rare type of hyperbolic galaxy collision events,
but a type that does exist in a simulated universe such as IllustrisTNG (for
details, see Shin et al., 2020). By varying the stellar feedback energy,
initial gas metallicity, disk inclinations and spins, we form a suite of
numerical experiments. The simulation configurations and the stellar masses of
the resulting DMDGs are listed in Table 1.

Figure 1:  $t=-100$, 50, 150, 800 Myr snapshots of the fiducial, idealized
collision of two identical gas-rich dwarf galaxies of $M_{200}=5.95\times
10^{9}\,M_{\odot}$ each, and the resulting DMDGs (the Fiducial run in Table
1). Surface densities of gas (top row), stars (middle row) and dark matter
(bottom row). $t=0$ is set to the moment when the two progenitor galaxies are
at pericentric approach. The co-planar galaxy collision (i.e., both progenitor
disks lying in the $x$-$y$ plane) with a relative velocity of 300 ${\rm
km\;s}^{-1}$ produces multiple DMDGs with stellar mass
$\sim$10${}^{8}\,M_{\odot}$ each. The dotted circles of 10 kpc radius are
centered on the most massive DMDG produced. See Section 3.1 for more
information.

The simulations are performed with the publicly available code Enzo (Bryan et
al., 2014; Brummel-Smith et al., 2019) with the Zeus hydrodynamics solver
(Stone & Norman, 1992a, b). Radiative gas cooling and heating rates along with
the metagalactic ultraviolet (UV) background are computed by the Grackle
library that interpolates the lookup table from Cloudy (Smith et al., 2017).
To describe feedback-regulated star formation, we employ the star-forming
molecular cloud model (SFMC; for details, see Kim et al., 2013, 2019), a
change from a simpler prescription in Shin et al. (2020). A particle
representing a SFMC forms when (1) the density of a gas cell exceeds $4\times
10^{4}$ cm-3, (2) the gas flow is converging, (3) the cooling time of the cell
is shorter than its dynamical time ($t_{\rm dyn}$), and (4) the created SFMC
particle is heavier than $M_{\rm SFMC}=10^{3}\,M_{\odot}$. Depicting
inefficient star formation in molecular clouds (Krumholz & Tan, 2007), the
SFMC particle returns 80% of its original mass in 12 $t_{\rm dyn}$,444 The
local star formation efficiency (SFE) per freefall time ($\epsilon_{\rm
ff}\sim 1\%$) is based on Kim et al. (2011, 2019). Other studies have also
found that varying the SFE parameter by not more than an order of magnitude
does not significantly affect the SC properties (Li et al., 2020; Ma et al.,
2020b). together with thermal energy peaking at 1 $t_{\rm dyn}$ (four
different energy choices for Fiducial to FB10 runs between 1051-52 erg per
150$\,M_{\odot}$ of stars formed) and metals (2% of the returned gas; see also
Cen & Ostriker, 1992; Kim et al., 2011).555 While it is possible that the
momentum imparted by supernovae could be underestimated in our thermal
feedback model, the subgrid physics implementation adopted here have
consistently shown that the the star formation is efficiently self-regulated
and the DMDGs are formed, across multiple resolution (resolution study at 80,
5, 1.25, and 0.625 pc) and feedback strengths (Fiducial, FB3, and FB7 run in
Table 1).

Figure 2:  Left panel: the radial profile of enclosed masses in a logarithmic
scale from the center of the most massive DMDG at $t=800$ Myr in the Fiducial
run. The DMDG consists of mostly stars, lacking in dark matter and gas. The
vertical magenta line denotes the half-mass radius $R_{\rm h}$ of the DMDG.
Middle panel: the dynamical mass within $r$ estimated from the individual SC’s
circular velocity, $M_{<}(r)=r\,v_{\rm circ}^{2}(r)/{\rm G}$, now in a linear
scale plot (circles). Among 42 gravitationally bound SCs the DMDG hosts, the
five most massive SCs are shown with colored markers. For comparison, the true
curve of the enclosed total mass is also included (dotted line with squares).
Right panel: the line-of-sight stellar velocities in a (5 kpc)2 box centered
on the DMDG at $t\,=\,800$ Myr, projected to a plane whose vertical coordinate
($z^{\,\prime}$) is aligned with the total angular momentum of the DMDG
(LDMDG). Colored circles denote the location of SCs and their line-of-sight
velocities with radii of $2\,(M_{\rm SC}/\,M_{\odot})^{1/3}$ kpc. Two
isosurface density contours are also shown: ${\rm
log[\Sigma_{\star}/(\,M_{\odot}\,pc^{-2})]}=0$ and 1. See Section 3.1 for more
information.

A pair of colliding galaxies is initialized in a $(1.311{\,\rm Mpc})^{3}$ box
of $64^{3}$ root grid. We adaptively refine the galaxies up to level $l_{\rm
max}=14$ achieving 1.25 pc resolution.666 We create a relaxed galactic disk by
gradually increasing the resolution for 250 Myr, thereby suppressing
artificial starbursts: maximum spatial resolution $=80$ pc for the first 100
Myr ($l_{\rm max}=8$), 20 pc for the next 100 Myr ($l_{\rm max}=10$), and 5 pc
for the last 50 Myr ($l_{\rm max}=12$). The relaxation step ends at
$t=-150\;{\rm Myr}$, the beginning of our 1.25 pc resolution simulation
discussed in this article ($t=0$ is when the two galaxies are at pericentric
approach; see, e.g., Figure 4). During the relaxation, SFMC formation is
either not allowed (i.e., dense gas clumps are supported only by the Jeans
pressure floor) or allowed with a coarser-resolution recipe. After the two
progenitor galaxies pass each other ($t=100$ Myr), we refine the cells only in
a pre-defined region that only includes the produced DMDGs. This RefineRegion
is initially set to a (16 kpc)3 box containing the most massive DMDG
throughout the run, and is surrounded by static nested boxes of cascading
resolution from level 0 to 8. At level $l$, a cell splits into eight child
cells whenever the cell contains more gas than $M^{l}_{\rm
ref,\,gas}=2^{-0.554(l-14)}\times M^{14}_{\rm ref,\,gas}$ or more particles
than $M^{l}_{\rm ref,\,part}=2^{-0.387(l-14)}\times M^{14}_{\rm ref,\,part}$
where $M^{14}_{\rm ref,\,gas}=2000\;{\rm M_{\odot}}\simeq 2\;M_{\rm
Jeans}^{180\,{\rm K}}$ and $M^{14}_{\rm ref,\,part}=4000\;{\rm
M_{\odot}}=4\;M_{\rm SFMC}$.

##  3 Results

###  3.1 Collision-induced DMDGs: Kinematics and Morphology

Figure 1 shows the collision sequence of two gas-rich dwarf galaxies
approaching with a relative velocity of $\sim 300\;{\rm km\;s}^{-1}$ (the
Fiducial run in the first row of Table 1; similar to what was extensively
tested in Shin et al. 2020). We define $t=0$ as the time when the two galaxies
are at pericentric approach. As shown in the top row, we find that extremely
dense gas clumps form during $t=0-50\;{\rm Myr}$ as the colliding columns of
gas undergo severe shock compression near the first contact point (the center
of each panel) and along the tidally stripped gas streams. The highly
compressed gas clumps of $\Sigma_{\rm gas}\sim 10^{3-4}\,{\rm M}_{\odot}\,{\rm
pc}^{-2}$ (see also Figure 3)—a typical value for extremely dense giant
molecular clouds—invokes the formation of DMDGs and their member SCs (Kim et
al., 2018; Lahén et al., 2020; Ma et al., 2020b). Strong bursts of star
formation continue inside the gas clumps during $t=50-150\;{\rm Myr}$ (middle
row). The newly formed stars remain near the first contact point and begin to
form one or more gravitationally bound structures without the help of pre-
existing dark matter potentials. In contrast, the two dark matter halos simply
pass through each other during this high-velocity collision. As a result, as
seen in the bottom row of Figure 1, hardly any dark matter particle remains in
the resulting DMDGs. Among the three DMDGs identified by the Hop halo finder
(Efstathiou et al., 1985) in the Fiducial run, the most massive one has
$M_{\star}\simeq 3.5\times 10^{8}\,M_{\odot}$ and almost no dark matter (see
also Figure 2). In Figure 1, the most massive DMDG found is marked with a
dotted circle in each panel. The DMDG drifts away slightly from the first
contact point by $t=800\;{\rm Myr}$ (the fourth panels in each row).

In Figure 2, we present the dynamical properties of the most massive DMDG
produced in the Fiducial run at $t=800\;{\rm Myr}$. The left panel displays
the radial profiles of the enclosed masses in different components: stars,
gas, and dark matter. The DMDG primarily consists of stars and is nearly
devoid of not only dark matter but also gas, which is in line with the
reported gas-poor nature of DF2 ($M_{\rm gas,\,DF2}\lesssim 3.15\times
10^{6}\,M_{\odot}$; Chowdhury, 2019). The half-mass radius $R_{\rm h}$ of the
DMDG’s stellar body is 0.026 kpc, which is $>10$ times smaller than that of a
typical ultra-diffuse galaxy ($\sim 1.5\;{\rm kpc}$) or that of DF2 ($\sim
2.2\;{\rm kpc}$; van Dokkum et al., 2015, 2018b). However, it is important to
note that our simulated DMDGs evolved only for $<800\;{\rm Myr}$ in an
idealized, isolated environment. Some morphological properties of a compact
stellar object, such as $M_{\rm gas}$ and $R_{\rm h}$ could be sensitive to
the energetic supernovae feedback (Trujillo-Gomez et al., 2021), or the tidal
interaction with massive neighbors.777For example, the DMDGs formed in an
isolated simulation may harbor more gas than DF2 does. However, Shin et al.
(2020) demonstrated that if the simulated DMDG orbits around another massive
galaxy, the tidal interaction could deprive the DMDG of most of its gas in a
few Gyr. $R_{\rm h}$ also changes frequently at this early stage of evolution
as a number of new SCs in the outskirts leave the DMDG, or because of the
dynamical relaxation or heating.

Figure 3:  Similar to Figure 1, but zoomed in on the most massive SC (“SC1”)
created around the most massive DMDG in the Fiducial run. $t=$ 70, 110, 160,
and 210 Myr snapshots of the gas surface density and stellar distribution in a
(1 kpc)2 box (top and middle rows), and in a (6 kpc)2 box (bottom row). The
red dots in the middle and bottom panels denote the stars in SC1 (the same
color codes and SCs as in Figure 2). The dotted circle of 1 kpc radius marks
the location of the DMDG in each of the bottom panels. See Section 3.2 for
more information.

We now study the kinematics of SCs formed simultaneously with the DMDG that
are identified by the Hop algorithm reconfigured to locate SCs. We fit each
SC’s stellar density profile to a power law,
$\rho(r)\,=\rho_{0}\,(1+r^{2}/a^{2})^{-{(\gamma+1)/2}}$ (Elson et al., 1987),
and obtain its scale radius $a$. $R=10\,a$ is adopted as the extent of each SC
when we determine its physical quantities such as mass, half-mass radius,
orbital velocity, average metallicity, and creation time. In the Fiducial run,
a total of 50 SCs are found within 12 kpc from the center of the DMDG at
$t=800\;{\rm Myr}$; among them, 42 are gravitationally bound to the DMDG, and
eight are unbound. The middle panel of Figure 2 presents the dynamical mass
within radius $r$ estimated from each SC’s circular velocity,
$M_{<}(r)=r\,v_{\rm circ}^{2}(r)/{\rm G}$. Shown for comparison is the true
curve of the DMDG’s enclosed mass at $r$. It is notable that the dynamical
masses estimated by the gravitationally bound SCs lie close to the true
enclosed mass. This implies that the gravitationally bound “member” SCs of the
DMDG are virialized, a key assumption when the dynamical masses of stellar
systems are observationally evaluated.

In the right panel of Figure 2, we inspect the dynamical structure of the DMDG
in the Fiducial run. It displays the line-of-sight velocities of the stars and
SCs in the DMDG, viewed from an angle that is perpendicular to the angular
momentum vector of the DMDG (LDMDG) at $t=800\;{\rm Myr}$. We first find that
the line-of-sight velocity dispersion of the 42 gravitationally bound SCs is
$\sigma=\,12.6\;{\rm km\;s}^{-1}$. The 10 most luminous bound SCs
($M_{V}<-8.5$; see Figure 6 and Section 3.2) have $\sigma=\,11.2\;{\rm
km\;s}^{-1}$. The small spread in velocity—smaller than what is normally
expected by its stellar mass—is an indication that the resulting galaxy has a
small dark-matter-to-stellar-mass ratio. Our finding is in line with the
observed value of $8.4\;{\rm km\;s}^{-1}$ for the luminous SCs in DF2 (van
Dokkum et al., 2018b) and $5.8\;{\rm km\;s}^{-1}$ in DF4 (van Dokkum et al.,
2019), considering that the observation is conducted in an oblique plane.
Second, from the “butterfly-shaped” line-of-sight velocities and the
isosurface density contours, we find that the DMDG is an oblate rotator with a
major-to-minor axis ratio $\sim 2$. The rotation velocity of the stars is
$\sim 15\;{\rm\,km\;s}^{-1}$ at $4\;R_{\rm h}<r<10\;R_{\rm h}$ with a
noticeable velocity gradient along the major axis of the system at $r\lesssim
R_{\rm h}$ ($v_{\rm circ,\,max}\sim 35\;\rm\,km\;s^{-1}$ at $r=25$ pc). For
comparison, the stellar body of DF2 is reported to be a prolate rotator with a
rotational velocity of $10.8-12.4\;{\rm km\;s}^{-1}$ (Emsellem et al. 2019 and
Lewis et al. 2020, respectively).

However, by comparing the specific angular momentum of the DMDG (lDMDG) and
that of the system of 42 bound “member” SCs around the DMDG’s center (lSC), we
discover that the SCs are not in a coherent rotation around the DMDG, nor
preferentially in the disk plane of the DMDG. (i.e., $|$l${}_{\rm
SC}|/|$l${}_{\rm DMDG}|=0.042$). The rotation axis of the system of SCs about
the DMDG—if it exists at all—is not aligned with that of the whole DMDG (i.e.,
lSC and lDMDG makes a 69∘ angle). The discrepancy in the dynamical states of
the DMDG and the system of its member SCs could be attributed to their
different formation processes (see Section 3.2).

Lastly, using the suite of simulations that we performed (Table 1), we remark
on how the formation of the collision-induced DMDGs depends on the simulation
parameters. (1) First, DMDGs form if a sufficiently large amount of gas
participates in the supersonic collision and forms self-gravitating gas
clumps. For example, in the Tilt-1 run, although the collision is no longer
co-planar, a portion of the disk gas along the $x$-axis (see Figure 1) still
clashes into a large column of gas in the other disk, producing two DMDGs. In
contrast, in the Tilt-2 run, the disk gas finds only a small column of gas in
the other disk during its first near head-on collision, so no DMDG forms (for
a quantitative discussion, see Section 3.3 of Shin et al., 2020). (2) The
strength of stellar feedback affects the internal and kinetic energy of the
ISM and dictates the DMDG formation. In the FB10 run where the feedback is 10
times more energetic than in the Fiducial run, clumps are easily dispersed via
stellar feedback, with no DMDG formed. (3) The metallicity of the disk gas
influences the DMDG formation. In the Low-metal run where the initial gas
metallicity is set to $10^{-2}\,Z_{\odot}$ (100 times lower than in the
Fiducial run), the relatively less metal-enriched gas slows down the gas
cooling and suppresses star formation. As a result, the mass of the resulting
DMDG is only one-third of what is found in the Fiducial run, hosting only nine
SCs (as opposed to 50 in the Fiducial run; see also Figure 5). (4) The
velocity shear in the gas at the first contact surface affects the dynamical
properties of the resulting DMDG. In the Fiducial run of a prograde–prograde
collision, the ISM experiences a strong shear as the colliding columns of gas
move in the opposite directions. Therefore, a DMDG is born in a clump that is
in coherent rotation, resulting in a rotation-supported structure. In
contrast, in the Retrograde run with a prograde–retrograde collision, the
velocity vectors in the collapsing gas clumps are highly disordered, resulting
in a DMDG with less rotation (i.e., $|$l${}_{\tt Retrograde}|/|$l${}_{\tt
Fiducial}|\sim 0.1$).

###  3.2 SCs in the Collision-induced DMDGs: Formation History and Key
Observables

We turn our attention to the formation history and the properties of SCs in
the collision-induced DMDGs. Figure 3 again shows the collision sequence of
the Fiducial run, but zooming in on the site where the most massive SC of the
most massive DMDG was formed (“SC1”; red dots). After a large fraction of the
disk gas loses its initial momentum during the first near head-on collision
and stalls at the first contact surface, pockets of dense, shock-compressed
gas emerge ($t=50\;{\rm Myr}$ panels in Figure 1). A few of these very dense
clumps ($\Sigma_{\rm gas}\gtrsim 10^{3-4}\,{\rm M}_{\odot}\,{\rm pc}^{-2}$)
rapidly cools and collapses to spawn clusters of stars before or while the
DMDG’s main stellar body forms ($M_{\star,\,\rm SC}\gtrsim 10^{6}\,M_{\odot}$
each; $t=70\;{\rm Myr}$ panels in Figure 3; Kim et al., 2018; Madau et al.,
2020). As the clumps continue to turn gas into stars in the next $\sim
50\;{\rm Myr}$, they are fed by the accreting cold gas while the frequent
supernova explosions from young SFMCs eject the hot gas out of the new SCs
($t=110\;{\rm Myr}$ panels in Figure 3; see also the bottom panel in Figure
4). The SCs consume or lose most of their gas by $t=160\;{\rm Myr}$, and
continue to orbit around the DMDG until the end of our simulation at
$t=800\;{\rm Myr}$.

Figure 4:  Top panel: the star formation history in the most massive DMDG and
its five most massive member SCs in the Fiducial run (the same SCs and colored
markers as in Figure 2). Bottom panel: zoomed in on the time that the five SCs
are created. The five SCs and most of the stars in the DMDG’s stellar body
form during $t=0-150$ Myr, immediately after the progenitor galaxies’
pericentric approach. See Section 3.2 for more information.  Figure 5:  Top
row: the stellar masses versus the average metallicities of the SCs formed in
the DMDG. Bottom row: the average creation times versus the average
metallicities of the SCs. The left panels are for the Fiducial run (42 bound
$+$ eight unbound SCs) and the right panels are for the Low-metal run (two
bound $+$ seven unbound SCs). Selected massive SCs are shown with the colored
markers (the same SCs and markers as in Figures 2 and 4 for the Fiducial run).
The average metallicities of SCs are correlated with their masses. We also
denote the average creation time and metallicity of the DMDG as a whole with a
black circle in the bottom panels. In the Low-metal run, the stars in the DMDG
is more metal-enriched than the two bound SCs. See Section 3.2 for more
information.

Figure 4 presents the star formation history in the most massive DMDG of the
Fiducial run and its five most massive member SCs. Most of the stars in the
DMDG and the five SCs are born in $<150\;{\rm Myr}$ after the progenitor
galaxies’ pericentric approach ($t=0$–150 Myr), which is consistent with our
earlier studies (Shin et al., 2020). After the initial burst of star
formation, the DMDG maintains to form stars at a modest pace ($\sim
10^{-1}-10^{-2}\,{\rm M}_{\odot}\,{\rm yr}^{-1}$), while none of the massive
SCs form additional stars at a meaningful rate after $t=150$ Myr as they are
nearly devoid of gas. The widths of star formation peaks for these SCs are
$\sim$50 Myr, which are significantly narrower than that of the DMDG.

Figure 6:  Top panel: the absolute magnitudes ($M_{V}$) versus the half-light
radii of the SCs formed in the most massive DMDG in the Fiducial run. $M_{V}$
is estimated from the SC’s mass, assuming the mass-to-light ratio $M/L_{V}=2$.
The five most massive SCs are shown with the colored markers (the same SCs and
markers as in Figures 2, 4 and 5). We overplot the luminous GCs located in DF2
($M_{V}<-8.5$; triangles from van Dokkum et al. 2018c; squares from Ma et al.
2020a); here the half-light radii and $M_{V}$ in the F606W band are used. See
Section 3.2 for more information.

Then in Figure 5, we make scatter plots of all the SCs formed with the DMDG on
the plane of the stellar mass versus the average metallicity, and on the plane
of the average creation time versus the average metallicity. The left panels
are for the Fiducial run and the right panels are for the Low-metal run (note
that the Low-metal run with the initial metallicity of $10^{-2}\,Z_{\odot}$
produces nine SCs compared to 50 in the Fiducial run; see Table 1 and Section
3.1 for more information). In the bottom panels, we also denote the average
creation time and metallicity of the stellar body of the DMDG with a black
circle. Two observations could be made. First, we find that the average
metallicity of a SC correlates positively with its mass (top row) rather than
with its creation time (bottom row). We argue that because a massive SC
creates more stars that consequently injects more metals into the ISM, and
retains a higher fraction of supernova ejecta in its deeper gravitational
potential well, the next generations of stars in the SC inherit higher
metallicities from the metal-enriched ISM (Bailin & Harris, 2009; Fensch et
al., 2014). Therefore, the more massive a SC is, the higher its average
stellar metallicity is. Second, readers may note that in the Low-metal run,
the average metallicity of the DMDG is higher than that of the two bound SCs
(colored markers in the bottom-right panel). A similar discrepancy is observed
in DF2: ${\rm[M/H]}\sim-1.07$ for DF2’s stellar body and $\sim-1.63$ for its
SCs (Fensch et al., 2019). The higher metallicity in the DMDG’s main stellar
body could be attributed to the continuing star formation after the initial
starburst. The later generations of stars formed in the more metal-enriched
ISM raise the average stellar metallicity of the DMDG, while the SCs’
metallicities rarely change once the initial starburst ceases.888 Note that
this trend is found only in the Low-metal run, but not in the Fiducial run
where the initial gas metallicity is already set to $1\,Z_{\odot}$.

Finally, we analyze the individual characteristics of the SCs, such as their
luminosities. Figure 6 shows the scatter plot on the plane of the absolute
magnitude ($M_{V}$) versus the half-light radius ($=$ half-mass radius $R_{\rm
h}$) for all the SCs in the most massive DMDG of the Fiducial run. The
absolute magnitude of a SC is estimated from its mass, assuming a typical
constant mass-to-light ratio $M/L_{V}=2$ (e.g., van Dokkum et al., 2018c). The
SC masses plotted in Figure 6 range from $1.2\times 10^{5}$ to $9.0\times
10^{6}\,M_{\odot}$. Among them, 32 gravitationally bound SCs are with
$-8.5<M_{V}<-6.5$, while 10 bound SCs are with $M_{V}<-8.5$. For comparison,
we include the luminous SCs observed in DF2 ($M_{V}<-8.5$; van Dokkum et al.,
2018c; Ma et al., 2020a). The absolute magnitudes of the luminous SCs are
similar in our Fiducial run and DF2. Note that the peak of our cluster
luminosity function (the bottom panel of Figure 6) lies between the peaks of
the Milky Way’s GCLF and DF2’s (van Dokkum et al., 2018c).999 84 (Milky Way)
$/$ 4 (DF2) GCs are with $-8.5<M_{V}<-6.5$, and 20 (Milky Way) $/$ 11 (DF2)
GCs are with $M_{V}<-8.5$. We argue that the extremely dense, shock-compressed
gas clumps induced by the high-velocity galaxy collision could give rise to
multiple luminous SCs—a scenario considered by Trujillo-Gomez et al. (2020,
2021) and demonstrated with simulations by Kim et al. (2018) and Madau et al.
(2020). Meanwhile, the half-light radii $R_{\rm h}$ of the simulated SCs are
1.5–3.0 times smaller than the estimated radii of the GCs in DF2. As was for
$R_{\rm h}$ of our simulated DMDG, it is important to note that our simulated
SCs evolved only for $<800\;{\rm Myr}$. If evolved for a sufficiently long
time, the radii of SCs could change due to, e.g., the dynamical heating by
their host galaxy.101010 For example, Gnedin et al. (1999) showed that during
orbital passages near the host galaxy, GCs could be dynamically heated by the
host. We also note that the detailed dynamics of supernova blast waves or
remnants is not fully resolved with 1.25 pc resolution (or SFMC particle of
$\gtrsim 200\,{\rm M}_{\odot}$), making it difficult to accurately estimate
the disruption time scale and the inner structure of each SC. This warrants an
even higher resolution simulation.

##  4 Conclusion

The diffuse dwarf galaxies DF2 and DF4, with their alleged absence of dark
matter and multiple luminous member GCs, have been a topic of intense debate
as the astrophysical community tries to comprehend their formation mechanism.
Using a suite of high-velocity galaxy collision simulations with 1.25 pc
resolution, we have demonstrated that such a collision induces the formation
of the DMDGs with $M_{\star}\sim 10^{8}\,M_{\odot}$. With numerical resolution
significantly improved from our previous study (Shin et al., 2020), we resolve
the dynamical structure of the produced DMDGs and find that a co-planar
prograde–prograde galaxy collision can produce a compact, rotation-supported,
oblate DMDG. We have also resolved and studied the SCs forming simultaneously
with the DMDGs, which later could become progenitors of GCs. A galaxy
collision with a high relative velocity of $\sim 300\;{\rm km\;s}^{-1}$
creates pockets of severely shock-compressed gas that rapidly turn into
multiple luminous SCs ($M_{\star}\gtrsim 10^{6}\,M_{\odot}$,
$M_{V}\lesssim-8.5$ mag) that are in line with the observed GCs in DF2. The
majority of stars in the DMDGs and their SCs form less than $150\;{\rm Myr}$
after the progenitor galaxies’ pericentric approach. Many properties of the
observed DMDGs’ GCs are reproduced. In our Fiducial run, 10 luminous,
gravitationally bound, virialized SCs have a line-of-sight velocity dispersion
$11.2\;{\rm km\;s}^{-1}$. In the Low-metal run with the initial gas
metallicity of $10^{-2}\,Z_{\odot}$, the DMDG shows a higher average stellar
metallicity than the SCs.

Our study suggests that a high-velocity galaxy collision is a promising
candidate for the simultaneous formation of DMDGs and their luminous member
GCs. In the forthcoming studies, we aim to examine the galaxy collision-
induced DMDG formation scenario with more sophisticated subgrid physics that
fully takes into account the evolution of interstellar medium due to energetic
stellar feedback with even higher resolution (see footnote 5), in a
cosmological context for $\sim 10$ Gyr, comparable to DF2’s age ($\sim 9$ Gyr;
van Dokkum et al., 2018c; Fensch et al., 2019) designed to follow its later
evolution.

The authors would like to thank Jérémy Fensch, Yongseok Jo, Taysun Kimm, Myung
Gyoon Lee, Boon Kiat Oh, Joseph Silk, Sebastian Trujillo-Gomez, and Sukyoung
Yi for insightful discussions. Ji-hoon Kim acknowledges support by Samsung
Science and Technology Foundation under Project Number SSTF-BA1802-04, and by
the POSCO Science Fellowship of POSCO TJ Park Foundation. His work was also
supported by the National Institute of Supercomputing and Network/Korea
Institute of Science and Technology Information with supercomputing resources
including technical support, grants KSC-2019-CRE-0163 and KSC-2020-CRE-0219.
The publicly available Enzo and yt codes used in this work are the products of
collaborative efforts by many independent scientists from numerous
institutions around the world. Their commitment to open science has helped
make this work possible.

## References

  * Bailin & Harris (2009) Bailin, J., & Harris, W. E. 2009, ApJ, 695, 1082, doi: 10.1088/0004-637X/695/2/1082
  * Bournaud et al. (2008) Bournaud, F., Duc, P. A., & Emsellem, E. 2008, MNRAS, 389, L8, doi: 10.1111/j.1745-3933.2008.00511.x
  * Brummel-Smith et al. (2019) Brummel-Smith, C., Bryan, G., Butsky, I., et al. 2019, JOSS Software, 4, 1636, doi: 10.21105/joss.01636
  * Bryan et al. (2014) Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014, ApJS, 211, 19, doi: 10.1088/0067-0049/211/2/19
  * Cen & Ostriker (1992) Cen, R., & Ostriker, J. P. 1992, ApJ, 399, L113, doi: 10.1086/186620
  * Chowdhury (2019) Chowdhury, A. 2019, MNRAS, 482, L99, doi: 10.1093/mnrasl/sly192
  * Danieli et al. (2020) Danieli, S., van Dokkum, P., Abraham, R., et al. 2020, ApJ, 895, L4, doi: 10.3847/2041-8213/ab8dc4
  * Efstathiou et al. (1985) Efstathiou, G., Davis, M., White, S. D. M., & Frenk, C. S. 1985, ApJS, 57, 241, doi: 10.1086/191003
  * Elson et al. (1987) Elson, R. A. W., Fall, S. M., & Freeman, K. C. 1987, ApJ, 323, 54, doi: 10.1086/165807
  * Emsellem et al. (2019) Emsellem, E., van der Burg, R. F. J., Fensch, J., et al. 2019, A&A, 625, A76, doi: 10.1051/0004-6361/201834909
  * Fensch et al. (2014) Fensch, J., Mieske, S., Müller-Seidlitz, J., & Hilker, M. 2014, A&A, 567, A105, doi: 10.1051/0004-6361/201423971
  * Fensch et al. (2019) Fensch, J., van der Burg, R. F. J., Jeřábková, T., et al. 2019, A&A, 625, A77, doi: 10.1051/0004-6361/201834911
  * Gnedin et al. (1999) Gnedin, O. Y., Lee, H. M., & Ostriker, J. P. 1999, ApJ, 522, 935, doi: 10.1086/307659
  * Hopkins et al. (2013) Hopkins, P. F., Cox, T. J., Hernquist, L., et al. 2013, MNRAS, 430, 1901, doi: 10.1093/mnras/stt017
  * Hunter (2007) Hunter, J. D. 2007, CSE, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jackson et al. (2021) Jackson, R. A., Kaviraj, S., Martin, G., et al. 2021, MNRAS, 502, 1785, doi: 10.1093/mnras/stab093
  * Kim et al. (2013) Kim, J.-h., Krumholz, M. R., Wise, J. H., et al. 2013, ApJ, 775, 109, doi: 10.1088/0004-637X/775/2/109
  * Kim et al. (2019) Kim, J.-h., Wise, J. H., Abel, T., et al. 2019, ApJ, 887, 120, doi: 10.3847/1538-4357/ab510b
  * Kim et al. (2011) Kim, J.-h., Wise, J. H., Alvarez, M. A., & Abel, T. 2011, ApJ, 738, 54, doi: 10.1088/0004-637X/738/1/54
  * Kim et al. (2018) Kim, J.-h., Ma, X., Grudić, M. Y., et al. 2018, MNRAS, 474, 4232, doi: 10.1093/mnras/stx2994
  * Kruijssen et al. (2012) Kruijssen, J. M. D., Pelupessy, F. I., Lamers, H. J. G. L. M., et al. 2012, MNRAS, 421, 1927, doi: 10.1111/j.1365-2966.2012.20322.x
  * Kruijssen et al. (2011) Kruijssen, J. M. D., Pelupessy, F. I., Lamers, H. J. G. L. M., Portegies Zwart, S. F., & Icke, V. 2011, MNRAS, 414, 1339, doi: 10.1111/j.1365-2966.2011.18467.x
  * Krumholz & Tan (2007) Krumholz, M. R., & Tan, J. C. 2007, ApJ, 654, 304, doi: 10.1086/509101
  * Lahén et al. (2019) Lahén, N., Naab, T., Johansson, P. H., et al. 2019, ApJ, 879, L18, doi: 10.3847/2041-8213/ab2a13
  * Lahén et al. (2020) —. 2020, ApJ, 891, 2, doi: 10.3847/1538-4357/ab7190
  * Lewis et al. (2020) Lewis, G. F., Brewer, B. J., & Wan, Z. 2020, MNRAS, 491, L1, doi: 10.1093/mnrasl/slz157
  * Li et al. (2017) Li, H., Gnedin, O. Y., Gnedin, N. Y., et al. 2017, ApJ, 834, 69, doi: 10.3847/1538-4357/834/1/69
  * Li et al. (2004) Li, Y., Mac Low, M.-M., & Klessen, R. S. 2004, ApJ, 614, L29, doi: 10.1086/425320
  * Li et al. (2020) Li, H., Vogelsberger, M., Marinacci, F., Sales, L. V., & Torrey, P. 2020, MNRAS, 499, 5862, doi: 10.1093/mnras/staa3122
  * Ma et al. (2020a) Ma, J., Wang, S., Wang, S., et al. 2020a, MNRAS, 496, 3741, doi: 10.1093/mnras/staa1775
  * Ma et al. (2020b) Ma, X., Grudić, M. Y., Quataert, E., et al. 2020b, MNRAS, 493, 4315, doi: 10.1093/mnras/staa527
  * Madau et al. (2020) Madau, P., Lupi, A., Diemand, J., Burkert, A., & Lin, D. N. C. 2020, ApJ, 890, 18, doi: 10.3847/1538-4357/ab66c6
  * Maji et al. (2017) Maji, M., Zhu, Q., Li, Y., et al. 2017, ApJ, 844, 108, doi: 10.3847/1538-4357/aa7aa1
  * Monelli & Trujillo (2019) Monelli, M., & Trujillo, I. 2019, ApJ, 880, L11, doi: 10.3847/2041-8213/ab2fd2
  * Moreno et al. (2019) Moreno, J., Torrey, P., Ellison, S. L., et al. 2019, MNRAS, 485, 1320, doi: 10.1093/mnras/stz417
  * Ogiya (2018) Ogiya, G. 2018, MNRAS, 480, L106, doi: 10.1093/mnrasl/sly138
  * Oliphant (2006) Oliphant, T. E. 2006, A Guide to NumPy, Vol. 1 (Trelgol Publishing USA)
  * Renaud et al. (2015) Renaud, F., Bournaud, F., & Duc, P.-A. 2015, MNRAS, 446, 2038, doi: 10.1093/mnras/stu2208
  * Ruiz-Lara et al. (2019) Ruiz-Lara, T., Trujillo, I., Beasley, M. A., et al. 2019, MNRAS, 486, 5670, doi: 10.1093/mnras/stz1237
  * Shen et al. (2021b) Shen, Z., van Dokkum, P., & Danieli, S. 2021b, ApJ, 909, 179, doi: 10.3847/1538-4357/abdd29
  * Shen et al. (2021a) Shen, Z., Danieli, S., van Dokkum, P., et al. 2021a, ApJ, 914, L12, doi: 10.3847/2041-8213/ac0335
  * Shin et al. (2020) Shin, E.-j., Jung, M., Kwon, G., et al. 2020, ApJ, 899, 25, doi: 10.3847/1538-4357/aba434
  * Silk (2019) Silk, J. 2019, MNRAS, 488, L24, doi: 10.1093/mnrasl/slz090
  * Smith et al. (2017) Smith, B. D., Bryan, G. L., Glover, S. C. O., et al. 2017, MNRAS, 466, 2217, doi: 10.1093/mnras/stw3291
  * Stone & Norman (1992a) Stone, J. M., & Norman, M. L. 1992a, ApJS, 80, 753, doi: 10.1086/191680
  * Stone & Norman (1992b) —. 1992b, ApJS, 80, 791, doi: 10.1086/191681
  * Trujillo et al. (2019) Trujillo, I., Beasley, M. A., Borlaff, A., et al. 2019, MNRAS, 486, 1192, doi: 10.1093/mnras/stz771
  * Trujillo-Gomez et al. (2020) Trujillo-Gomez, S., Kruijssen, J. M. D., Keller, B. W., & Reina-Campos, M. 2020, MNRAS, 506, 4841, doi: 10.1093/mnras/stab1895
  * Trujillo-Gomez et al. (2021) Trujillo-Gomez, S., Kruijssen, J. M. D., & Reina-Campos, M. 2021, arXiv e-prints, arXiv:2103.08610.  https://arxiv.org/abs/2103.08610
  * Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9, doi: 10.1088/0067-0049/192/1/9
  * van Dokkum et al. (2019) van Dokkum, P., Danieli, S., Abraham, R., Conroy, C., & Romanowsky, A. J. 2019, ApJ, 874, L5, doi: 10.3847/2041-8213/ab0d92
  * van Dokkum et al. (2018a) van Dokkum, P., Danieli, S., Cohen, Y., Romanowsky, A. J., & Conroy, C. 2018a, ApJ, 864, L18, doi: 10.3847/2041-8213/aada4d
  * van Dokkum et al. (2018b) van Dokkum, P., Danieli, S., Cohen, Y., et al. 2018b, Nature, 555, 629, doi: 10.1038/nature25767
  * van Dokkum et al. (2018c) van Dokkum, P., Cohen, Y., Danieli, S., et al. 2018c, ApJ, 856, L30, doi: 10.3847/2041-8213/aab60b
  * van Dokkum et al. (2015) van Dokkum, P. G., Abraham, R., Merritt, A., et al. 2015, ApJ, 798, L45, doi: 10.1088/2041-8205/798/2/L45
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261, doi: 10.1038/s41592-019-0686-2


# Combined Sieve Algorithm for Prime Gaps

Seth Troisi
<EMAIL_ADDRESS>
###### Abstract

A new Combined Sieve algorithm is presented with cost proportional to the
number of enumerated factors over a series of intervals. This algorithm
achieves a significant speedup, over a traditional sieve, when handling many
($[10^{4},10^{7}]$) intervals concurrently. The speedup comes from a space-
time tradeoff and a novel solution to a modular equation. In real world tests,
this new algorithm regularly runs 10,000x faster. This faster sieve paired
with higher sieving limits eliminates more composites and accelerates the
search for large prime gaps by 30-70%. During the development and testing of
this new algorithm, two top-10 record merit prime gaps were discovered.

_K_ eywords primes, prime differences, prime gaps, parallel sieve

##  1 Introduction

A prime gap is the difference between successive primes. A record is kept of
the smallest known prime for each gap. For example the smallest prime followed
by a gap of 4 is 7; a gap of 6 first follows 23; and a gap of 96 first follows
360,653. Until his death in 2019, Thomas R. Nicely maintained a list of record
prime gaps, first in print and later on his faculty website [1]. Since 2019, a
group from MersenneForum.org, including the author, has maintained the records
on GitHub[2].

The search for prime gaps precedes computers but the digital era accelerated
the search. Computer searches for prime gaps start circa 1930 when D.H. Lehmer
constructed a physical computer from bicycle chains, The Lehmer Sieve, and
used it to construct a table of the primes less than 37 million [3]. In 1961
Baker, Gruenberger, And Armerding added to the search in [4] [4] Both listed
first occurrences of prime gaps. Lander and Parkin using a CDC 6600 extended
the search and published a paper exclusively on prime differences in 1967 [5].
At least a dozen papers from as many authors have followed since.

This paper outlines a new algorithm that speeds up searches for large gaps by
combining the sieving of many nonadjacent but regularly spaced intervals at a
substantial speedup. With a faster sieve more small factors can be found
reducing the number of candidates tested per interval. Ready-to-run code,
benchmarks, and validation can be found on
GitHub111https://github.com/sethtroisi/prime-gap.

###  1.1 New Prime Gaps with Merit 38.048 & 37.051

Large gaps become more common as the density of primes decrease. The merit of
a prime gap, $g_{k}$, following the prime, $p_{k}$, is defined as
$\frac{g_{k}}{\log(p_{k})}$. It is the ratio of the gap to the “average” gap
near that point. Appendix 3 shows the merit of each record prime for gaps
under $100,000$. Gaps with merit $>35$ are exceedingly rare, only 30 are
known, so finding two using this new algorithm is both an accomplishment and
involved quite a bit of luck. The two new prime gaps are:

A prime gap of 35,308 follows $\bm{100054841\frac{953\\#}{210}-9670}$ with
merit = 38.048 this is the 4th largest known merit222https://primegap-list-
project.github.io/lists/top20-overall-merits/.

A prime gap of 36,702 follows $\bm{164065661\frac{1019\\#}{210}-18566}$ with
merit = 37.051 this is the 9th largest know merit.

##  2 Searching for Prime Gaps

An exhaustive search for the first occurring prime gap has been performed up
to $2^{64}$[6]. Thomas Nicely explains the search process

> No general method more sophisticated than an exhaustive search is known for
> the determination of first occurrences and maximal prime gaps. As in the
> present study, this is most efficiently done by sieving successive blocks of
> positive integers for primes, recording the successive differences, and thus
> determining directly the first occurrences and maximal gaps. [1]

The largest known first occurring prime gap is the gap of $1550$ following
$18361375334787046697$ with merit
$\frac{1550}{\log(1.836e19)}=\textbf{34.944}$ found by Bertil Nyman in 2014
and independently confirmed in 2018[2]. A history of the upper bound for
exhaustive searches is found on Thomas Nicely’s page [7] and an updated
version on GitHub [6].

##### Larger Gaps

The natural approach to search for larger gaps is similar. First, an interval
around a test number, $N$, is sieved using small primes to mark numbers with
small factors as composite. Then a primality test is run on the remaining
numbers to find the ”surrounding” primes, the next and previous prime. The
prime gap containing $N$ is the difference between the two surrounding primes.
An existing respected program PGsurround.pl[8] use this approach with the
additional optimization of searching only $N$ of a certain forms known to
produce larger gaps.

The speed a program using this approach can find the primes surrounding $N$ is
determined by the speed of the primality test and how many composite numbers
are removed by the sieve. Generally less than 2% of execution time is spent on
the sieve stage; the Combined Sieve doesn’t change this dynamic, even with the
deepest sieves, sieving up to $10^{11}$, and smallest $N$, $N\sim 3,000$ bits,
sieving takes less than $10\%$ of amortized runtime (see Appendix 2).

##### Primorials

It’s possible to construct arbitrarily large gaps by looking near primorials.

Let $K=P\\#=\displaystyle\prod_{p\in\mathbb{P}}^{P}p\ $

$\displaystyle\begin{array}[]{rll}2\mid&(K+2)\implies&K+2\notin\mathbb{P}\\\
3\mid&(K+3)\implies&K+3\notin\mathbb{P}\\\
2\mid&(K+4)\implies&K+4\notin\mathbb{P}\\\ &...&\\\ &(K+X)&\\\ &...&\\\
P\mid&(K+P)\implies&K+P\notin\mathbb{P}\end{array}$

Thus there are no primes in the interval $[K+2,K+P]$ and the gap from
$prev\\_prime(K+2)$ to $next\\_prime(K+2)$ is $\geq P$. If both $K\pm 1$ are
composite (by far the average case) then the gap is $\geq 2P+3$. Note that as
the gap grows larger so does $\log(K)$ so merit doesn’t grow arbitrarily
large.

The average gap can be increased by choose a small high composite divisor $d$
(e.g. $11\\#$) and looking near $K=\frac{P\\#}{d}$. To illustrate why,
consider $K=\frac{P\\#}{2}$ all odd $X$, even those $>P$, are multiples of two
($K$ is odd so $K+X$ is even) and multiples of small primes eliminate all even
$X\leq 2P$ that aren’t powers of two.

Because of the way primorials align their prime factors a greater share of
nearby numbers are composite, meaning larger gaps. For this reason nearly all
of the records ($\approx 97\%$ of the 95,000+ record gaps with gap $>2000$)
are of the form $m\frac{P\\#}{d}-a$, $m$ being a small multiple, $a$ being the
offset to the first prime below $m\frac{P\\#}{d}$.

###  2.1 Sieving

For each $N$, an interval (typically a few thousand numbers corresponding to
average or optimistic prime gap) is sieved around $N$ and numbers without any
factors, heading away from $N$, are tested for primality till the surrounding
primes are found.

The traditional sieve calculates, for all small primes $p$, $N\equiv r\bmod p$
and uses that to mark $N-r\pm i\times p$ as composite. E.g. $N\equiv 5\mod
73\implies 73\mid N-5,\ \ 73\mid N-5-73\text{,\ \ and\ \ }73\mid N-5+73$.

###  2.2 Probabilistic Primality Testing

The best known deterministic primality tests are much slower than their
probabilistic counterparts with respectively runtime of
$\tilde{\mathcal{O}}(\log(n)^{6})$ and $\tilde{\mathcal{O}}(\log(n)^{2})$. As
the probability of a composite being marked as prime can be reduced to an
arbitrary level (BPSW has no know counterexamples) and no primes is ever
marked as composite, it’s advantageous to use a probabilistic primality tests.
Notable gaps, with high merit, often have their endpoints later certified as
prime with a deterministic test.

GMP version 6.2.0 [9] mpz_probab_prime_p uses a BPSW primality test. GMP is
well optimized for numbers of all sizes and is a natural choice for primality
testing333Some of the inspiration for this paper came from analysis done while
improving GMP’s next_prime code. Some testing was done with PrimeForm/GW which
is faster for large numbers (1000+ digits) but more difficult to work with.

##  3 Combined Sieve Algorithm

For reasons described above, the majority of prime gap searches are centered
around numbers of the form $N=m\frac{P\\#}{d}$. Often millions of values of
$m$ are tested for the same $K=\frac{P\\#}{d}$. The Combined Sieve algorithm
sieves $M$ sequential intervals simultaneously at a dramatically speedup to
sieving $M$ intervals serially. The algorithm works on any arithmetic sequence
but for the purpose of prime gap it’s beneficial to focus on numbers of this
deficient primorial form.

Figure 1: Speedup for various values of $M_{c}$ with
$K=\frac{1511\\#}{7\\#},X=18,000,$ and $P_{limit}=10^{11}$

Note that $m$’s that share a factor with $d$ are skipped partially to avoid
duplicate testing (e.g. $m=100,d=35$ would be a repeat of $m=20,d=7$) and
partially because small $d$ tend to generate lower merit gaps. Figure 1
accounts for this by counting the number of $m$ coprime to $d$ less than $M$,
$M_{c}$ = $|\\{m:1\leq m\leq M\land(m,d)=1\\}|$

###  3.1 Space–Time Tradeoff

By caching $K\bmod p$, $N\bmod p$ is reduced to a machine word sized modulo.

$\begin{split}N\bmod p=m\frac{P\\#}{d}\bmod p&=mK\bmod p\\\
&\equiv(\underbrace{(m\bmod p)}_{\text{$<$64bits}}\underbrace{(K\bmod
p)}_{\text{500-20,000 bits}})\bmod p\end{split}$ (1)

After this trade-off only the first $m$ is slow for the other $M_{c}-1$ values
of $m$ each $N\bmod p$ can be computed with a memory lookup, a 64 bit multiply
and a 64 bit modulo irregardless of the size of $K$. This means the max
speedup is $M_{c}$ since $K\bmod p$ is still calculated once for each $p$.
This upper bound is easily seen in Figure 1 up to $M_{c}\sim 10^{5}$.

There are two easy ways to cache $K\bmod p$. $K\bmod p$ can be caching for all
primes than two nested for-loops iterate over $m,p$ this requires storing the
modulo for all primes in memory. The other method inverts the order and
iterating over $p,m$ this requires storing the list of coprime $m$ and a very
large boolean array in memory of the composite status of all offsets for all
$M_{c}$ intervals.

###  3.2 Large Prime Optimization

This section details a complicated trick that allows the algorithm to avoid
calculating $mK\bmod p$ for each $m$ and instead directly compute the set of
$m$ where $p$ divides a number in the intervals around $N=mK$.

Assume $X$ numbers before and after $N$ are sieved for an interval of size
$S=2X+1$, on average $\frac{S}{2}$ numbers will be divisible by $2$,
$\frac{S}{3}$ by $3$, $\frac{S}{5}$ by $5$. In general, each prime $p$ will
divide $\lceil\frac{S}{p}\rceil$ or $\lfloor\frac{S}{p}\rfloor$ numbers in the
interval around $N$. When $p\gg S$ most primes will not divide a number in the
interval. Mertens’ second theorem states

$\sum_{p}^{n}\frac{1}{p}=\log(\log(n))+M_{m}+o(1)\\\ $ (2)

Where $M_{m}\approx\ 0.2615$ is the Meissel-Mertens constant.

With $X=100000,N=\frac{503\\#}{210}$, all of the $17,984$ primes $2\leq p\leq
S$ divide one or more number in the interval. $\frac{24744}{60514}=41\%$ of
the primes $S<p\leq 10^{6}$ divide a number in the interval.
$\frac{81355}{50769036}=0.16\%$ of the primes $10^{6}<p\leq 10^{9}$ divide a
number in the interval. This means that $99.84\%$ of large primes for any
given $m$ don’t divide a number in the interval but still require processing.
As the sieve limit increases, to $10^{10}$, $10^{11}$, $10^{12}$, …, the
fraction of primes dividing a number in the interval decreases further.

####  3.2.1 Testing many intervals

With a fixed prime $p$ and fixed $K$ and $1\leq m\leq M$, sieving can be
restated as trying to finding $m$’s where

$\begin{split}-X\leq\ &mK\leq X\bmod p\\\ &\equiv\\\ 0\leq\ &mK+X\leq 2X\bmod
p\\\ \end{split}$ (3)

For each $p$ larger than some threshold (naively $p>10S$ i.e. $p$ divides
fewer than 1 in 10 intervals), $K\bmod p$ is calculated, and the smallest $m$
satisfying Equation 3 is found. $mK\bmod p$ is computed and one of
$\lfloor\frac{mK}{p}\rfloor p$ or $\lceil\frac{mK}{p}\rceil p$ is marked as a
composite. The algorithm continues finding larger $m$ and marking off
composites till $m>M$.

Because we’re only interested in $m$ which are coprime to $d$, work is wasted
when $m$ shares a factor with $d$ and the search advances to the next $m$
without marking a composite. The smaller and more numerous the factors of $d$,
the more unnecessary work. Very large $d$ are uncommon so a reasonable worst
case is $d=17\\#$ where $81.9\%$ of $m$ share a factor with $d$.

Each solution to Equation 3 (ignoring $m$ which share a factor with $d$ and
the last solution with $m>M$) is another factor of a number in one of the
interval. This is a very powerful transformation, by solving Equation 3 the
Combined Sieve algorithm can enumerate $m$’s where $p$ divides a number in the
interval $[mK-X,mK+X]$. Each prime divides an average $\frac{MS}{p}$ numbers
because $p\gg S$ this is a large reduction from testing all $M_{c}$ intervals
individually. With $S=60,001$ and $P\sim 10^{10}$ only 1 in $160,000$ $m$ will
need to calculate $mK\bmod p$.

Equation 2 gives an estimate of the number of primes $p\leq P_{limit}$ that
divide a number in each interval as $\approx S\ (\log(\log(P_{limit}))+o(1))$.
This means that all prime factors $10S\leq p\leq n$ for all $M$ intervals can
be enumerated by solving Equation 3 $\approx M\times
S(\log(\log(P_{limit}))-\log(\log(\text{10S})))+\text{PrimePi}(P_{limit})-\text{PrimePi}(10S)$
times444This is one of the many theoretical results verified during execution.
The count of all $m$ tested over all $p$ generally matches to $\approx
0.05\%$.

###  3.3 Solving Equation 3

$\text{Let }G(p,K,L,R)=m\text{ be the smallest solution to }L\leq mK\leq
R\bmod p$ (4)

$\displaystyle G(p,K,L,R)=$ $\displaystyle 0$ $K=0\ \|\ L=0$ (5a)
$\displaystyle G(p,K,L,R)=$ $\displaystyle G(p,K\bmod p,L,R)$ $K\geq p$ (5b)
$\displaystyle G(p,K,L,R)=$ $\displaystyle\lceil L/K\rceil$ $\lceil
L/K\rceil\leq\lfloor R/K\rfloor$ (5c) $\displaystyle G(p,K,L,R)=$
$\displaystyle G(p,p-K,p-R,p-L)$ $2K\geq p$ (5d) $\displaystyle G(p,K,L,R)=$
$\displaystyle\lceil(L+\frac{p\times
G(K,-p\bmod{K},L\bmod{K},R\bmod{K}))}{K}\rceil$ (5e)

This solution assumes $p$ and $K$ are coprime which is guaranteed by the
definition of $K$ as a primorial and $p$ a large prime (at minimum $p>S>L,R$).
$p$ and $K$ coprime also guarantees the existence of a solution $(L\times
K^{-1})\bmod p$.

A handwavy proof that 5e returns the minimal answer can be constructed by
inspecting 4 with and without the $\bmod\ p$ condition. If Equation 4 has a
solution without considering $\bmod\ p$, the minimal solution is the first
multiple of $K\geq L\implies m=\lceil L/K\rceil$ , otherwise $L\leq Km-ip\leq
R$ for some $i\geq 1$.

$\displaystyle\begin{array}[]{rlll}L\leq&\ mK-ip&\leq R&=\\\ L-mK\leq&-i\times
p&\leq R-mK&\\\ \\\
\lx@intercol\hfil\text{consider}\bmod{K}\hfil\lx@intercol\\\ \\\
L\leq&i\times-p&\leq R\bmod{K}\\\ \lx@intercol\hfil
i=G(K,-p\bmod{K},L\bmod{K},R\bmod{K})\hfil\lx@intercol\\\ &\implies\\\
\lx@intercol\hfil
m=\lceil\frac{L+pG(K,-p\bmod{K},L\bmod{K},R\bmod{K})}{K}\rceil\hfil\lx@intercol\end{array}$

5d which reduces $K>\frac{p}{2}$ to $K<\frac{p}{2}$ is the result of some
rearrangement.

$\displaystyle\begin{split}\text{Let }G(p,K,L,R)=m&\text{ be the smallest
solution to }\\\ L\leq mK\leq R\bmod{p}&\equiv\\\ -L\geq m(-K)\geq-R&\equiv\\\
-R\leq m(-K)\leq-L&\equiv\\\ \text{Which is equivalent to}&\\\
G(p,-K,-R,-L)&\equiv G(p,p-K,p-R,p-L)\end{split}$

Equation 5e either directly returns $m$ (5a, 5c) or reduces to a new sub-
problem with $K<\frac{p}{2}$ (5d) or $p<K$ (5e). It follows that $G$ is
$\mathcal{O}(\log(p))$.

Given $p<2^{64}$ and $K\bmod p$ is cached, this runs quickly. In practice it
takes $2\log_{2}(\frac{p}{R-L})\leq 2\log_{2}(\frac{10^{13}}{10000})\leq 60$
steps. On an i7-2600K with $p$ a 50 bit prime and $R-L=S=100,000$ an optimized
implementation takes 100-200 ns555See https://github.com/sethtroisi/prime-
gap/blob/main/BENCHMARKS.md.

##  4 Parameter Selection

###  4.1 Sieve Limit

Mertens’ third theorem estimates the fraction of numbers that don’t have a
factor less than $P_{limit}$

$\prod_{p}^{P_{limit}}\frac{p-1}{p}\approx\frac{1}{\log(P_{limit})e^{\gamma}}$
(6)

Where $\gamma\approx 0.577$ is the Euler–Mascheroni constant.

$P_{limit}$ | primes | theoretical | experimental
---|---|---|---
$10^{4}$ | 1229 | 0.060960 | 0.06088
$10^{5}$ | 9592 | 0.048768 | 0.04875
$10^{6}$ | 78498 | 0.040640 | 0.040638
$10^{8}$ | 5761455 | 0.030480 | 0.030480
$4\times 10^{9}$ | 189961812 | 0.025394 | 0.025394
$10^{10}$ | 455052511 | 0.024384 | 0.024384
$10^{11}$ | $4.1\times 10^{9}$ | 0.022167 |
$10^{12}$ | $3.8\times 10^{10}$ | 0.020320 |
$10^{13}$ | $3.5\times 10^{11}$ | 0.018757 |

Joined with the fact that a number near $N$ is prime
$\approx\frac{1}{\log(n)}$ of the time, the expected number of primality tests
that will need to be performed to find a prime after removing numbers with
factors less than $P_{limit}$ is

$E(N,P_{limit})\approx\frac{\log(N)}{\log(P_{limit})e^{\gamma}}$ (7)

A sieve limit can be chosen such that the quasi-linear cost of calculating
$K\bmod p$ and solving Equation 3 multiplied by the expected number of
primality tests reduced is near the cost of running the avoided primality
tests.

In practice combined_sieve computes the derivative of Equation 7 at regular
intervals giving an estimate of primality tests avoided per second of
additional sieving. When that rate falls below the benchmarked rate of
primality tests per second, combined_sieve advises the user to stop sieving.
This allows for better control of sieve limit based on CPU and memory
characteristics.

A counter-intuitively observation is the optimal sieve limits can decrease as
$N$ increases. The optimal limit strongly depends on $M$ and for large $N$ it
doesn’t make sense to sieve a lot of intervals when only a few can be tested
per day.

Anecdotally setting $P_{limit}$ much higher than $5SM$ is sub-optimal as the
ratio of time spent performing $K\bmod p$ starts to dwarf the time spent
solving Equation 3. Some values used during exploratory searches are listed in
Appendix 1.

###  4.2 Sieve Interval

For any $N$, if both of the surrounding primes is within $X$ of $N$ then time
was wasted finding factors of numbers further away. If either of the
surrounding primes is further away than $X$ additional traditional sieving
will need to be done. Choosing $X$ in $[6\log(N),10\log(N)]$ is reasonable. In
practice combined_sieve performs a small test sieve of a few $m$ and issues a
warning if the estimated probability of the surrounding primes exceeding $X$
is greater than 0.5%.

##  5 Runtime

Let

$\displaystyle S=$ $\displaystyle\ 2X+1=\text{the sieve interval}$
$\displaystyle M=$ number of combined sieves to run simultaneously
$\displaystyle M_{c}=$ number of $m$ coprime to $d$ $\displaystyle C_{mult}=$
$\displaystyle\ \text{small prime multiplier constant}\approx
10\frac{M}{M_{c}}\sim 50$ $\displaystyle P_{limit}=$ sieve limit
$\displaystyle|P|=$ $\displaystyle\
\text{PrimePi}(P_{limit})\approx\frac{P_{limit}}{\log(P_{limit})}$
$\displaystyle|P_{small}|=$ $\displaystyle\
\text{PrimePi}(C_{mult}S)\approx\frac{C_{mult}S}{\log(C_{mult}S)}$
$\displaystyle|P_{large}|=$ $\displaystyle\ |P|-|P_{small}|$

Note: This section ignores the small one-time cost of finding primes up to
$P_{limit}$.

Sieving each coprime $m$ without the Space-Time Tradeoff and Large Prime
Optimization involves computing $(m\times K)\bmod p$ then marking off an
average of $\frac{S}{p}$ composites. $N\bmod p$ can be computed in $\log(N)$
time,

$|P|\log(N)+S(\log(\log(P_{limit}))+M_{m})$ (8)

The Combined Sieve algorithm computes $(K\bmod p)$ once and reused this for
each $m$. Small primes, which are likely to divide a number in the sieve
interval (e.g. $\frac{p}{S}<C_{mult}$) are handled similarly to the
traditional sieve but use cached $K\bmod p$ to calculating $N\bmod p$ with 64
bit integers. For larger primes the $(m\times K)\bmod p$ calculation is
avoided for the vast majority of $m$ at the amortized cost of $\frac{MS}{p}+1$
calls to Equation 3. The cost for $M$ intervals is

$\underbrace{|P|\ \log(K)}_{K\bmod p\text{ for all
primes}}+\underbrace{M_{c}S(\log\log(P_{limit})+M_{m}))}_{\text{marking off
all multiples}}+\\\ \underbrace{M_{c}|P_{small}|}_{m\times K\bmod p\text{ for
small
primes}}+\underbrace{(MS(\log\log(P_{limit})-\log\log(C_{mult}S))+|P_{large}|)}_{\text{calls
to Equation
\ref{eq:modeq}}}\times\underbrace{\log_{2}(\frac{P_{limit}}{S})}_{\text{cost
to solve Equation \ref{eq:modeq}}}$ (9)

With

$\displaystyle K=$ $\displaystyle\frac{8887\\#}{13\\#}$
$\displaystyle\log_{2}(K)=$ $\displaystyle 12,651\ \
(\log(P\\#)\approx\log(e^{P})\approx P)$ $\displaystyle M=$ $\displaystyle
10,000$ $\displaystyle M_{c}=$ $\displaystyle 1,917$ $\displaystyle
P_{limit}=$ $\displaystyle 10^{10}$ $\displaystyle|P|=$ $\displaystyle
455,052,511$ Optimizing the searching above and below for up to 10 merit
$\displaystyle S=$ $\displaystyle 2X+1=2\times 10\times 9000+1=180,001$
$\displaystyle C_{eq\ref{eq:modeq}}=$ $\displaystyle 30\ \text{experimental
constant for speed of Equation \ref{eq:modeq}}$ $\displaystyle C_{mod}=$
$\displaystyle\frac{1}{16}\ \text{experimental constant for speed of modulo}$

The experimental speedup of the Combined Sieve over the traditional sieve is

$\begin{split}&|P|\log_{2}(K)C_{mod}+M_{c}S(\log(\log(P_{limit}))+M_{m})+\\\
&M_{c}|P_{small}|+(MS(\log(\log(P_{limit}))-\log\log(C_{mult}S))+|P_{large}|)C_{eq\ref{eq:modeq}}\log_{2}(\frac{P_{limit}}{S})\\\
\approx\ &455,052,511\frac{12,651}{16}+1,917\times 180,001(3.14+0.26)+\\\
&1,917\times 602,492+(10,000\times 180,001\times(3.14-2.86)+452,618,854)\times
30\log_{2}(\frac{10^{10}}{400,001})\\\ \approx\ &3.60\times 10^{11}+1.17\times
10^{9}+1.15\times 10^{9}+5.24\times 10^{11}\\\ \approx\ &8.86\times 10^{11}\\\
\textbf{vs}&\\\
&M_{c}(|P|\log_{2}(K)C_{mod}+S(\log(\log(P_{limit}))+M_{m}))\\\ \approx\
&1,917\times(455,052,511\frac{12,651}{16}+180,001(3.14+0.26))\\\ \approx\
&6.90\times 10^{14}\\\\[3.00003pt] \frac{6.90\times 10^{14}}{8.86\times
10^{11}}=\ &\ 778\text{x speedup}\end{split}$ (10)

This result matches experimental observations of speedups. In general the
speedup mostly depends on the magnitude of $M$ (The speedup is bounded by
$M_{c}$, $1,917x$ in this example), moderately on $P_{limit}$, and a small
amount on the magnitude of $N$.

##  6 Results

Traditionally sieving with primes twice as large, means finding half as many
factors in the same amount of time. Using the Combined Sieve algorithm $K\bmod
p$ takes the same amount of time irrespective of $p$ but larger primes require
solving Equation 3 fewer times. As long as a reasonable percentage of time is
spent solving Equation 3 vs calculating $K\bmod p$, the sieve speeds up as it
processes larger primes. Anecdotally it takes less than twice as long to sieve
to a 10x higher sieve limit (e.g. if it takes 2 hours to handle primes up to
50 million, all primes up to 500 million can be handle in less than 4 hours).

Use of the Combined Sieve allows for substantial higher sieve limits which
reduces the number of primality tests that need to be perform. Appendix 1 has
estimates of new optimal sieve limits. Calculating the resulting speedups
compared to the sieve limits used by PGSurround.pl the overall speedup from
using the Combined Sieve is 30-70%.

##### Real world searching

A number of complex relationships govern overall speed. To maximize the number
of records found a wide search strategy is utilized, various combinations of
$\frac{P\\#}{d}$ are sieved with $M=2,000$. For each $m_{n}$ an auxiliary
program computes $Prob_{n}=P(\text{record gap}|m_{n}\times K,\text{sieved
interval})$, for the top 5% $\\{Prob_{n},m_{n},K,\text{sieved interval}\\}$ is
stored in a database. Several background threads read from the database
testing the highest likelihood intervals. If
$\overline{\frac{Prob}{\text{time}}}$ for a $\frac{P\\#}{d}$ is above average,
combined_sieve is run with the largest $M$ that fits in memory and a higher
sieve limit. A number of additional metrics and stats are collected to verify
assumptions and integrity during this process; see Appendix 4 for plots of
theoretical and experimental distributions.

###  6.1 Future work

combined_sieve contains a number of other significant speedups not discussed
in this paper. Section 6 briefly describes an optimization involving directly
computing the probability of a record gap for each interval using the
remaining numbers in an interval. This can be utilized to only test the most
likely 5-30% of $m$ or to avoid testing the high sides if previous prime is
close to $N$.

##### Acknowledgments

The author would like to thank Craig Citro for advice on publication and Iulia
Turc for regular check-ins on progress.

## References

  * [1] Thomas R. Nicely  “New maximal prime gaps and first occurrences” An updated copy is available at https://faculty.lynchburg.edu/~nicely/gaps/gaps.html In _Math. Comput._ 68, 1999, pp. 1311–1315
  * [2] Graham Higgins  “Prime gap list” https://primegap-list-project.github.io/lists/top20-overall-merits/, https://github.com/primegap-list-project/prime-gap-list, GitHub repository, 2020
  * [3] J. Selfridge and D.. Lehmer  “Tables Concerning the Distribution of Primes up to 37 Millions”, 1959
  * [4] F. Gruenberger and G. Armerding  “Statistics on the First Six Million Prime Numbers”, Paper  Rand Corporation, 1961
  * [5] L.. Lander and T.. Parkin  “On First Appearance of Prime Differences”  In _Mathematics of Computation_ 21, 1967, pp. 483–488
  * [6] Graham Higgins  “Exhaustively analyzed gaps” https://primegap-list-project.github.io/fully-analyzed/, GitHub repository, 2020
  * [7] Thomas R. Nicely  “Some Results of Research in Computational Number Theory”  URL: https://faculty.lynchburg.edu/~nicely/
  * [8] Dana Jacobsen  “PGsurround.pl” A historical version can be found at http://sti15.com/misc/tgaps5-30.txt
A more recent version was provided by Robert W. Smith and is available on
GitHub https://gist.github.com/sethtroisi/8f0ade362690c6bd416f025bc38df2be
A second version of unknown date is available on GitHub
https://gist.github.com/sethtroisi/ac50b89b8041b2c78f3dbcc1f9cc0116

  * [9] Torbjörn Granlund and the GMP development team  “GNU MP: The GNU Multiple Precision Arithmetic Library” http://gmplib.org/, 2020

## Appendix

Table 1: Experimental Sieve Limits $P$ | $\log(P\\#)$ | sieve limit | expected primality tests |  PGSurround.pl sieve limit | Theoretical Speedup
---|---|---|---|---|---
1511 | 1469 | $250\times 10^{9}$ | 31.5 | $5.0\times 10^{6}$ | 70.2%
2221 | 2054 | $800\times 10^{9}$ | 40.9 | $1.4\times 10^{7}$ | 71.8%
4441 | 4353 | $2000\times 10^{9}$ | 86.3 | $1.3\times 10^{8}$ | 51.4%
5333 | 5228 | $2000\times 10^{9}$ | 103.7 | $2.3\times 10^{8}$ | 47.0%
8887 | 8779 | $1000\times 10^{9}$ | 178.4 | $1.1\times 10^{9}$ | 32.5%
Table 2: Combined Sieve Amortized Runtime (time per coprime m) $M$ | $d$ | $X$ | $M$ | sieve limit | Combined Sieve(s) | Testing(s)
---|---|---|---|---|---|---
1511 | $7\\#\times 1009$ | 16,000 | 10,000,000 | $100\times 10^{9}$ | 0.0112 | 0.134
1511 | - | - | - | $250\times 10^{9}$ | 0.0129 | 0.127
2221 | $5\\#\times 2027$ | 23,000 | 10,000,000 | $800\times 10^{9}$ | 0.0252 | 0.568
4441 | $11\\#\times 4001$ | 47,000 | 10,000,000 | $2000\times 10^{9}$ | 0.0605 | 6.04
5333 | $11\\#\times 5003$ | 65,000 | 1,000,000 | $2000\times 10^{9}$ | 0.158 | 11.23
8887 | $17\\#\times 8009$ | 100,000 | 100,000 | $1000\times 10^{9}$ | 0.709 | 72.16
Figure 3: Merit of record gap for gaps up to 100,000 Figure 4: Stats and
distribution for $K=\frac{3001\\#}{2190}\ \ M=1+[0,100000)$

Top row: Probability of gap $X$, Distribution of expected gap, CDF of expected
gap
Middle row: Sum of probability of merit 12+ gap after $M$ tests, histogram of
probabilities, CDF of probability
Bottom row: Sum of probability of record gap after $M$ tests, histogram of
probabilities, CDF of probability


# Reducing Spatial Discretization Error with Linear Discontinuous Source
Tilting in
Iterative Quasi-Monte Carlo for Neutron Transport

Samuel Pasmann1 1University of Notre Dame, Department of Aerospace and
Mechanical Engineering, Fitzpatrick Hall, Notre Dame, IN 46556,
<EMAIL_ADDRESS>2Oregon State University, Department of Nuclear Science and
Engineering, 1791 SW Campus Way, Corvallis, OR 3North Carolina State
University, Department of Mathematics, 3234 SAS Hall, Box 8205, Raleigh NC
27695  Ilham Variansyah2 1University of Notre Dame, Department of Aerospace
and Mechanical Engineering, Fitzpatrick Hall, Notre Dame, IN 46556,
<EMAIL_ADDRESS>2Oregon State University, Department of Nuclear Science and
Engineering, 1791 SW Campus Way, Corvallis, OR 3North Carolina State
University, Department of Mathematics, 3234 SAS Hall, Box 8205, Raleigh NC
27695  C.T. Kelley3 1University of Notre Dame, Department of Aerospace and
Mechanical Engineering, Fitzpatrick Hall, Notre Dame, IN 46556,
<EMAIL_ADDRESS>2Oregon State University, Department of Nuclear Science and
Engineering, 1791 SW Campus Way, Corvallis, OR 3North Carolina State
University, Department of Mathematics, 3234 SAS Hall, Box 8205, Raleigh NC
27695  and Ryan G. McClarren1 1University of Notre Dame, Department of
Aerospace and Mechanical Engineering, Fitzpatrick Hall, Notre Dame, IN 46556,
<EMAIL_ADDRESS>2Oregon State University, Department of Nuclear Science and
Engineering, 1791 SW Campus Way, Corvallis, OR 3North Carolina State
University, Department of Mathematics, 3234 SAS Hall, Box 8205, Raleigh NC
27695

##  1 Introduction

###  1.1 The Iterative Quasi-Monte Carlo Method

Recently, iterative Quasi-Monte Carlo (iQMC) was introduced as a new method of
neutron transport which combines deterministic iterative methods and quasi-
Monte Carlo simulation for more efficient solutions to the neutron transport
equation [1]. Primary advantages of iQMC include a vectorized multigroup
scheme, $O(N^{-1})$ solution convergence, and the use of advanced iterative
Krylov solvers, including GMRES and BiCGSTAB, which converge with far fewer
iterations than the standard source iteration.

iQMC treats the scattering and fission terms as internal fixed sources, and
thus the Monte Carlo transport sweep is reduced to a particle ray trace which
provides a well-suited application for Quasi-Monte Carlo. Quasi-Monte Carlo is
the use of low-discrepancy sequences in place of pseudo-random number
generators and provides a more efficient sampling of the phase space for a
given number of particles. The improved sampling technique provides a
theoretical $O(N^{-1})$ convergence compared to the $O(N^{-1/2})$ of pseudo-
random numbers.

In each iQMC sweep, particles are emitted with an initial position and
direction assigned from the LDS and a statistical weight calculated from a
piecewise-constant source over a given mesh. $N$ particles are created and
traced out of the volume, tallying the scalar flux with a path-length tally
estimator. After a complete sweep, the new scalar flux approximation is sent
to the iterative solver to update the source strength. iQMC was shown to
achieve $O(N^{-1})$ across multiple test problems, however was also observed
to plateau in convergence in some problems [1, 2].

Previous iQMC results utilized a uniform Cartesian grid with a piecewise-
constant source. Similar to “teleportation error” in Implicit Monte Carlo
(IMC) methods [3], the spatial discretization and piecewise-constant source
can lead to a significant spatial error that limits convergence of the overall
method. Taking concepts from IMC, we have developed a history-based
discontinuous piecewise-linear source tilting scheme to reduce spatial error
in iQMC. The source tilting method is described below and afterward we present
results from a fixed-source 2D reactor-like problem adapted from the Takeda-1
Benchmark problem [4].

##  2 METHODOLOGY

###  2.1 Piecewise-Constant Scheme

In a 2-dimensional example, previous iQMC studies utilized a piecewise-
constant flux approximation

$\phi_{\mathrm{Constant}}=a_{i,j},$ (1)

where $i$ and $j$ denote the 2-dimensional spatial indices and $a_{i,j}$ is
the cell-averaged flux

$a_{i,j}=\frac{1}{\Delta x_{i}\,\Delta
y_{j}}\int_{y_{j-1}}^{y_{j}}\int_{x_{i-1}}^{x_{i}}\phi(x,y)\,dx\,dy.$ (2)

Given a path-length tally estimator with path length $S$, total cross section
$\Sigma_{t}$, and continuous particle weight capture with initial weight
$w_{0}$, the resultant tally is

$\int_{0}^{S}w_{0}e^{-\Sigma_{t}S^{\prime}}\,dS^{\prime}=w_{0}\left(\frac{1-e^{-\Sigma_{t}S}}{\Sigma_{t}}\right).$
(3)

###  2.2 Piecewise-Linear Scheme

In our discontinuous linear scheme, the scalar flux is now represented as

$\phi_{\mathrm{Linear}}=a_{i,j}+b_{i,j}(x-x_{\mathrm{mid},i})\\\
+c_{i,j}(y-y_{\mathrm{mid},j})+d_{i,j}(x-x_{\mathrm{mid},i})(y-y_{\mathrm{mid},j}),$
(4)

where $x_{\mathrm{mid},i}$ and $y_{\mathrm{mid},j}$ represent the cell
midpoint. The $a_{i,j}$ term is the same from the piecewise-constant scheme in
Eq. 2, while the $b_{i,j}\,$, $c_{i,j}\,$, $d_{i,j}\,$ terms respectively
represent the linear flux-tilt in the $x$, $y$, and $xy$ directions, where the
linear terms $b_{i,j}\,$ and $c_{i,j}\,$ are

$b_{i,j}=\frac{12}{\Delta x_{i}^{3}\Delta
y_{j}}\int_{y_{j-1}}^{y_{j}}\int_{x_{i-1}}^{x_{i}}\left(x-x_{\mathrm{mid},i}\right)\phi(x,y)\,dx\,dy,$
(5)

and

$c_{i,j}=\frac{12}{\Delta x_{i}\Delta
y_{j}^{3}}\int_{y_{j-1}}^{y_{j}}\int_{x_{i-1}}^{x_{i}}\left(y-y_{\mathrm{mid},j}\right)\phi(x,y)\,dx\,dy,$
(6)

and the bilinear term is

$d_{i,j}=\frac{12}{(\Delta x_{i}\Delta
y_{j})^{3}}\int_{y_{j-1}}^{y_{j}}\int_{x_{i-1}}^{x_{i}}\left(x-x_{\mathrm{mid},i}\right)\\\
\left(y-y_{\mathrm{mid},j}\right)\phi(x,y)\,dx\,dy,$ (7)

The new terms can be similarly tallied as

$\int_{0}^{S}w_{0}e^{-\Sigma_{t}S^{\prime}}\left[\left(x_{0}+\mu_{x}S^{\prime}\right)-x_{\mathrm{mid}}\right]\,dS^{\prime},$
(8)

for the $b_{i,j}\,$ and $c_{i,j}\,$ terms (swapping x-variables for
y-variables respectively) and

$\int_{0}^{S}w_{0}e^{-\Sigma_{t}S^{\prime}}\left[\left(x_{0}+\mu_{x}S^{\prime}\right)-x_{\mathrm{mid}}\right]\\\
\left[\left(y_{0}+\mu_{y}S^{\prime}\right)-y_{\mathrm{mid}}\right]\,dS^{\prime},$
(9)

for $d_{i,j}\,$, where $x_{0}$ and $y_{0}$ denote particle initial position
associated with the initial weight $w_{0}$.

##  3 TEST PROBLEM

###  3.1 2D Fixed-Source Reactor Problem

To evaluate the effect of using the proposed linear source, we’ve designed a
2-dimensional, 2-group, fixed-source reactor problem inspired by the Takeda-1
Benchmark problem [4] and using the same cross sections. The problem features
a core region surrounded by reflective boundary conditions and reflector
material, making the system near critical ($k=0.96883\pm 0.00049$). A fixed
source was placed outside the fuel in the reflector region. Figure 1 depicts
the problem setup. A reference scalar flux result is obtained from a high-
fidelity Monte Carlo simulation generated using the Monte Carlo code MC/DC [5]
with $10^{10}$ particle histories.

The iQMC simulations were run with a $25\mathrm{X}25$ uniform mesh similar to
other benchmark results of the Takeda-1 problem [4]. The Halton Sequence was
used to generate particle positions and angles, while GMRES was used to
iterate to $\Delta\phi/\phi\leq 1\times 10^{-9}$. Figure 2 shows the total
(summed across groups) scalar flux from a resulting iQMC simulation with
source tilting. Figures 3 and 4 show the source strength from a piecewise-
constant simulation, while Figures 5 and 6 depict the source strength from a
piecewise-linear simulation. Finally, Figure 7 displays the $L_{\infty}$-norm
of the scalar flux relative error for iQMC simulations with and without source
tilting as a function of the number of particles per iteration.

Figure 1: Core configuration of the test problem inspired by the Takeda-1
Benchmark [4]. Figure 2: Total scalar flux results from iQMC simulation with
source tilting and $N=2\times 10^{6}$ particles per iteration. Figure 3: Final
piecewise-constant thermal group source results from $N=2\times 10^{6}$
particles per iteration. Figure 4: Final piecewise-constant fast group source
results from $N=2\times 10^{6}$ particles per iteration. Figure 5: Final
piecewise-linear thermal group source results from $N=2\times 10^{6}$
particles per iteration. Figure 6: Final piecewise-linear fast group source
results from $N=2\times 10^{6}$ particles per iteration.

Figure 7: Scalar flux convergence as a function of particle histories per
iteration ($N$). With a piecewise-constant source the error plateaus almost
immediately. However, with a linear approximation iQMC is able to converge at
the expected rate of $O(N^{-1})$.

##  4 CONCLUSIONS

We have developed a history-based linear discontinuous source tilting scheme
for the iterative Quasi-Monte Carlo (iQMC) method. The source tilting
technique was shown to reduce spatial error in a 2-dimensional, 2-group,
reactor problem. Figure 7 shows that the relative error plateaus from
piecewise-constant simulations, while the piecewise-linear results are able to
converge at the expected $O(N^{-1})$. It is unlikely however, that piecewise-
linear source tilting completely eliminates the spatial discretization error.
More tests are needed to further explore the limitations of the method and to
evaluate performance on more difficult 3D problems.

##  5 Acknowledgments

This work was funded by the Center for Exascale Monte- Carlo Neutron Transport
(CEMeNT) a PSAAP-III project funded by the Department of Energy, grant number:
DE- NA003967 and the National Science Foundation, grant number DMS-1906446.

## References

  * [1] S. PASMANN, I. VARIANSYAH, C. T. KELLEY, and R. MCCLARREN, “A Quasi-Monte Carlo Method with Krylov Linear Solvers for Multigroup Neutron Transport Simulations,” _Nuclear Science and Engineering_ (2022).
  * [2] S. PASMANN, I. VARIANSYAH, C. KELLEY, and R. G. MCCLARREN, “iQMC: Iterative Quasi-Monte Carlo for k-Eigenvalue Neutron Transport Simulations,” in “International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering,” Niagara Falls, Ontario, Canada (8 2023).
  * [3] R. T. WOLLAEGER, A. B. WOLLABER, T. J. URBATSCH, and J. D. DENSMORE, “Implicit Monte Carlo with a Linear Discontinuous Finite Element Material Solution and Piecewise Non-Constant Opacity,” _Journal of Computational and Theoretical Transport_ , 45, 123–157 (2 2016).
  * [4] T. TAKEDA and H. IKEDA, “3-D Neutron Transport Benchmarks,” _Journal of Nuclear Science and Technology_ , 28, 656–669 (1991).
  * [5] I. VARIANSYAH, J. P. MORGAN, J. NORTHROP, K. E. NIEMEYER, and R. G. MCCLARREN, “Development of MC/DC: a performant, scalable, and portable Python-based Monte Carlo neutron transport code,” in “International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering,” Niagara Falls, Ontario, Canada (8 2023).