be prepared once for the nominal analysis, even if this is extended as
described in Appendix H.

Unlike the preparation of the CB and WB libraries, the actual analysis is
certain to be repeated with different modelling assumptions (Section 6),
and/or with the same ones to check for dependence on the starting point. The
amount of information required in the random access memory should not pose a
problem (Table 25). It should take a few hours to conduct a gradient descent
taking 1000 iterations with $\approx 10$ evaluations each time (one trial
evaluation and one for each model parameter to obtain the gradient by finite
differencing, see Table 19).

This would essentially avoid the burn-in period with the MCMC stage. Even so,
the large number of free parameters means this stage might require $10^{5}$
iterations. In this case, the whole analysis should take about a day. This
would allow variations to be explored, for instance changing the prior on
$q_{\text{ext}}$ and $q_{\text{int}}$, possibly in distinct ways (Section 6).
Some tests are required to ensure the minimum found is global, though these
could be based on gradient descent starting at different points rather than on
MCMC. This may well be possible because the cost function is envisaged to be
fully deterministic, so it should be differentiable. It is likely that the
MCMC analysis will also be more resilient with a deterministic cost function.
In addition, relatively secure modelling assumptions could be tested purely by
gradient descent, especially if the minimum moves very little when these
assumptions are changed in some way.

Our overall assessment is that if local WBs are sufficiently well observed by
Gaia and if the main source of contamination to the WBT is undetected CBs with
properties that are independent of the WB (except for the upper limit to
$a_{\text{int}}$), then the WBT is possible using about 60 cores
simultaneously. Therefore, the WBT is expensive but probably feasible using
currently available technology.

## Acknowledgements

IB is supported by Science and Technology Facilities Council grant
ST/V000861/1 led by Hongsheng Zhao, the primary objective of which is to
implement the WBT. IB is grateful for support from a “Pathways to Research”
fellowship from the University of Bonn. The authors acknowledge helpful
discussions with Benoit Famaey and Elena Asencio regarding this plan. They
also thank Rodrigo Ibata for demonstrating operation of the 2CB preparation
algorithm on the computing cluster at Strasbourg after some numerical
modifications he prepared.

## References

  * Asencio et al. (2021) Asencio E., Banik I., Kroupa P., 2021, MNRAS, 500, 5249
  * Banik (2019) Banik I., 2019, MNRAS, 487, 5291
  * Banik & Kroupa (2019) Banik I., Kroupa P., 2019, MNRAS, 487, 2665
  * Banik & Zhao (2018a) Banik I., Zhao H., 2018a, MNRAS, 473, 419
  * Banik & Zhao (2018b) Banik I., Zhao H., 2018b, MNRAS, 480, 2660
  * Beech (2009) Beech M., 2009, MNRAS, 399, L21
  * Begeman et al. (1991) Begeman K. G., Broeils A. H., Sanders R. H., 1991, MNRAS, 249, 523
  * Bekenstein & Milgrom (1984) Bekenstein J., Milgrom M., 1984, ApJ, 286, 7
  * Belokurov et al. (2020) Belokurov V., et al., 2020, MNRAS, 496, 1922
  * Clarke (2020) Clarke C. J., 2020, MNRAS, 491, L72
  * El-Badry (2019) El-Badry K., 2019, MNRAS, 482, 5018
  * El-Badry et al. (2019) El-Badry K., Rix H.-W., Tian H., Duchêne G., Moe M., 2019, MNRAS, 489, 5822
  * Faber & Gallagher (1979) Faber S. M., Gallagher J. S., 1979, ARA&A, 17, 135
  * Fletcher & Powell (1963) Fletcher R., Powell M. J. D., 1963, The Computer Journal, 6, 163
  * Francis & Anderson (2014) Francis C., Anderson E., 2014, Celestial Mechanics and Dynamical Astronomy, 118, 399
  * Gaia Collaboration (2018) Gaia Collaboration 2018, A&A, 616, A1
  * Gaia Collaboration (2021a) Gaia Collaboration 2021a, A&A, 649, A1
  * Gaia Collaboration (2021b) Gaia Collaboration 2021b, A&A, 649, A9
  * Gentile et al. (2011) Gentile G., Famaey B., de Blok W. J. G., 2011, A&A, 527, A76
  * Hagen et al. (2019) Hagen J. H. J., Helmi A., de Zeeuw P. T., Posti L., 2019, A&A, 629, A70
  * Hernandez et al. (2012) Hernandez X., Jiménez M. A., Allen C., 2012, European Physical Journal C, 72, 1884
  * Lelli et al. (2017) Lelli F., McGaugh S. S., Schombert J. M., Pawlowski M. S., 2017, ApJ, 836, 152
  * McMillan (2017) McMillan P. J., 2017, MNRAS, 465, 76
  * Milgrom (1983) Milgrom M., 1983, ApJ, 270, 365
  * Milgrom (1986) Milgrom M., 1986, ApJ, 302, 617
  * Milgrom (2010) Milgrom M., 2010, MNRAS, 403, 886
  * Öpik (1924) Öpik E., 1924, Publications of the Tartu Astrofizica Observatory, 25, 1
  * Ostriker & Peebles (1973) Ostriker J. P., Peebles P. J. E., 1973, ApJ, 186, 467
  * Pecaut & Mamajek (2013) Pecaut M. J., Mamajek E. E., 2013, ApJS, 208, 9
  * Pittordis (2019) Pittordis C., 2019, PhD thesis, Queen Mary University of London
  * Pittordis & Sutherland (2018) Pittordis C., Sutherland W., 2018, MNRAS, 480, 1778
  * Pittordis & Sutherland (2019) Pittordis C., Sutherland W., 2019, MNRAS, 488, 4740
  * Riello et al. (2021) Riello M., et al., 2021, A&A, 649, A3

##  Appendix A Required input data

The following observables will be required for each star in a WB within 250 pc
of the Sun that satisfies the quality criteria discussed in Section 4 (see
also Pittordis & Sutherland, 2019):

  1. 1.

position, parallax, and proper motion, and their uncertainties (5 variables
$\times\,2$),

  2. 2.

correlation coefficients between the above quantities (10 variables),

  3. 3.

Gaia-band apparent magnitude and uncertainty (1 variable $\times\,2$),

  4. 4.

quality flags, especially the goodness of fit of the Gaia single-source
astrometric solution.

These are required for WBs with $r_{\text{sky}}=1.5-40$ kAU and
$\widetilde{v}_{\text{sky}}\leq 8$ to capture WBs that might be scattered into
the photographic range in at least one MC trial (Section 4). Further
restrictions can be applied before the main analysis based on the MC trials
discussed in Section 4.1.

##  Appendix B Restricting the MCMC analysis to physical parameter values

Some of the model parameters (Table 19) have physical restrictions on the
allowed range, e.g. the slope of $P\left(e\right)$ must lie in the range
$\left(-2,2\right)$ to avoid negative probabilities (Equation 8). To implement
such restrictions, we will repeat the part of the MCMC analysis where the
proposal is set using Gaussian random numbers. Normally, only one proposal is
drawn per MCMC step. However, we will draw as many proposals as are necessary
to obtain one where every parameter is within its physically allowed range.
The cost of obtaining a new proposal without evaluating its binomial
likelihood is assumed to be negligible in comparison to the other costs, even
if several hundred unphysical proposals must first be discarded at every MCMC
step.

##  Appendix C Derivation of stellar mass from luminosity

Gaia astrometry determines the trigonometric parallax and thus the
heliocentric distance of a star. Uncertainties in this will be propagated
using the $5\times 5$ covariance matrix, as discussed in Section 4.1. For a
known parallax, the apparent Gaia-band magnitude $G$ can be converted into an
absolute magnitude $m_{{}_{G}}$, where we use the same notation as Pittordis &
Sutherland (2019) to distinguish apparent from absolute magnitudes in
different colour bands.

The stellar mass $M$ is then found from $m_{{}_{G}}$ using a Newton-Raphson
algorithm based on the tabulated mass-luminosity relationship in table 5 of
Pecaut & Mamajek (2013), which is available at https://www.pas.rochester.edu/
emamajek/EEM_dwarf_UBVIJHK_colors_Teff.txt with updated values. Since this is
given in the $V$-band, the first step is to guess the mass and thereby obtain
$m_{{}_{V}}$ and the $V-I$ colour. The latter is used to obtain the $G-V$
colour using the first Johnson-Cousins relation in table C2 of Riello et al.
(2021). This allows the determination of $m_{{}_{G}}$, whose difference with
the true value is the error of the Newton-Raphson algorithm.

In this way, we will obtain a relation between $\ln M/M_{\odot}$ and
$m_{{}_{G}}$. We will then fit a cubic to it using similar techniques to those
verified in Asencio et al. (2021). This cubic fit can readily be inverted
using the Newton-Raphson algorithm to obtain the mass from $m_{{}_{G}}$ during
each MC trial (Section 4.1). The use of a fitting function also avoids
problems of multiple masses having the same luminosity due to small-scale
structure in the mass-luminosity relation, which could otherwise make an
inversion very difficult.

##  Appendix D Finding WBs in the Gaia catalogue

Starting with a list of $\approx 10^{7}$ point-like sources with reliable
parallaxes $>4$ mas, we will search for pairs with a small angle on the sky.
This will require 2D binning in the Galactic latitude $b$ and longitude $l$.
Since the distribution of stars is close to isotropic within 250 pc, each star
only needs to be compared with $\approx 10^{3}$ stars because the maximum
allowed angular separation corresponds to an angular area $\approx
10^{4}\times$ smaller than the whole sky. Assuming each pairwise comparison
can be done in under 10,000 ns, the entire pair-finding operation should have
a manageable cost.

If the latitude $b$ is such that the star is close to either Galactic pole,
the coordinate system becomes nearly singular. Thus, a revised coordinate
system will be used if $\left|b\right|>45^{\circ}$. The basic idea is to use
the $x$-axis as the polar direction. To keep the system right-handed and
minimize the adjustments, this pole switch will be implemented using

$\displaystyle x\leftrightarrow z\,,~{}y\leftrightarrow-y\,.$ (37)

In the new coordinate system, the latitude is

$\displaystyle b~{}=~{}\arcsin z\,.$ (38)

To find the longitude, we will first find

$\displaystyle r~{}\equiv~{}\sqrt{x^{2}+y^{2}}\,.$ (39)

If $r$ is very small, the longitude $l$ is irrelevant, and can be set to 0 for
safety. However, this case should not arise due to the above-mentioned safety
catch. As a result, $l$ will be found using a previously verified algorithm
that implements

$\displaystyle l~{}=~{}\arcsin\left(\frac{y}{r}\right)$ (40)
$\displaystyle\text{if}~{}x<0$ $\displaystyle\quad l~{}=~{}\mathrm{\pi}-l$
(41) $\displaystyle\text{elseif}~{}y<0$ $\displaystyle\quad
l~{}=~{}l+2\mathrm{\pi}$ (42) end

To speed up the calculations, the coordinates for every star will be stored in
both systems. If a star is close to either pole in one coordinate system, then
it and the few degree region around it will be close to the equator in the
other coordinate system. Thus, only the coordinates in the latter system are
necessary. We will exploit this by storing the coordinates of a star in a
system only if the star has $\left|b\right|<80^{\circ}$ in that system. Since
the search region is much smaller than $35^{\circ}$, the binary finding
algorithm can be run entirely within one or other coordinate system.

If it becomes necessary to transform not only sky positions but also proper
motions, then the directions of the polar and azimuthal vectors are required
in both systems. These can be defined independently of the coordinate system.
For some axis $\widehat{\bm{n}}$, the polar vector is the vector on the
surface of the sphere which most rapidly goes towards the pole. It is
therefore

$\displaystyle\widehat{\bm{n}}_{{}_{b}}~{}\propto~{}\widehat{\bm{n}}-\left(\widehat{\bm{n}}\cdot\widehat{\bm{r}}\right)\widehat{\bm{r}}\,.$
(43)

In a right-handed system, the azimuthal vector is

$\displaystyle\widehat{\bm{n}}_{{}_{l}}~{}\propto~{}\widehat{\bm{n}}\times\widehat{\bm{r}}\,.$
(44)

In both cases, the proportionality constant is $1/\sin b$, which can be found
most efficiently using
$1/\sqrt{1-\left(\widehat{\bm{n}}\cdot\widehat{\bm{r}}\right)^{2}}$.

Once the above-mentioned basic 2D tree structure is used to identify relevant
pairs for comparison, the WB selection criteria will have to be imposed. There
should be exactly one other source within the relevant 3D volume, so the
parallax information must be considered as well. Since parallaxes are
relatively less precise than $r_{\text{sky}}$, it will be sufficient for the
two parallaxes to be consistent within the $\approx 2\sigma$ combined
uncertainty. Importantly, there should be no third star close by on the sky
with a parallax consistent within $\approx 3\sigma$ regardless of its proper
motion, since even an unbound star passing close to a WB would still perturb
it.

The precise selection criteria must be discussed with observers. We expect to
consider only WBs with $\widetilde{v}_{\text{sky}}<8$. The uncertainty on
$\widetilde{v}_{\text{sky}}$ is also important $-$ we expect to only consider
systems where this is $<0.5$, possibly with a tighter tolerance at low
$\widetilde{v}$ (Section 6.3).

##  Appendix E Assigning a systemic RV

The systemic RV of a WB affects $\widetilde{v}_{\text{sky}}$ by inducing an
apparent extra motion along the sky-projected separation of the WB (El-Badry,
2019; Banik, 2019). For small angular separation WBs, the induced apparent
velocity has magnitude $v_{r}\,r_{\text{sky}}/d$, where $v_{r}$ is the
heliocentric RV, and $d$ is the heliocentric distance to the geometric centre
of the WB. For a typical WB separated by $r_{\text{sky}}=10$ kAU at $d=100$
pc, this contribution is 25 m/s for $v_{r}=50$ km/s, corresponding to
$\Delta\widetilde{\bm{v}}_{\text{sky}}=0.06$ along the line between the WBs if
both stars are Sun-like. For systems with other properties,
$\Delta\widetilde{\bm{v}}_{\text{sky}}$ could be larger. Thus, the systemic RV
cannot be neglected. However, the difference in RV between the stars in a WB
can safely be neglected if it is bound with a gravity law not too different
from Newtonian. The possibility of LOS contamination is handled in the
statistical analysis (Section 5), and is expected to be rare (see Section 1.2
and Pittordis & Sutherland, 2019).

Despite its importance, the systemic RV of each WB will generally be unknown
prior to dedicated follow-up of the WBs. To estimate $v_{r}$, the basic
principle will be to determine the relative likelihood that the WB is part of
the Galactic thin or thick disk. A random number will be used to select the
component, and then to assign $v_{r}$ based on its mean and dispersion at the
position of the WB. In Galactocentric cylindrical polar coordinates, the
velocity ellipsoid will be assumed constant over the survey volume as this is
only a small fraction of the Sun’s distance from the Galactic centre. As an
example, the velocity ellipsoid tilts by $\approx 6^{\circ}$ per kpc in a
cylindrically aligned system, so the tilt over 250 pc should be very small
(Hagen et al., 2019).

For each Galactic component, the mean $v_{r}$ must be found based on the
assumption of the disk components rotating around the Galactic centre at a
speed differing from that of the Sun by the appropriate combination of the
Oort constants. The non-circular motion of the Sun also needs to be considered
(Francis & Anderson, 2014).

The dispersion in $v_{r}$ for each component is

$\displaystyle{\sigma_{{}_{\text{LOS}}}}^{2}=\sum_{i=r,\phi,z}{\sigma_{i}}^{2}\left(\widehat{\bm{\sigma}}_{i}\cdot\widehat{\bm{r}}\right)^{2}\,,$
(45)

where $\sigma_{i}$ is the velocity dispersion along Galactocentric cylindrical
polar direction $\widehat{\bm{\sigma}}_{i}$, and $\widehat{\bm{r}}$ is the sky
position of the WB’s geometric centre.

###  E.1 Selecting the Galactic component

The main difficulty is in finding the appropriate membership probability
$P_{j}$ to Galactic component $j$. We will start with a position-based prior

$\displaystyle\rho_{j}~{}\propto~{}\rho_{j,\text{mid}}\exp\left(-\frac{\left|z\right|}{h_{j}}\right)\,,$
(46)

where $\rho_{j,\text{mid}}$ is the mid-plane density of component $j$ with
exponential scale height $h_{j}$ in the vertical direction $z$. Other
$z$-dependencies can also be considered, for instance
$\operatorname{sech}^{2}$.

To find how this prior is altered by the observed $\bm{v}_{t}$, it is
necessary to find the sky-projected velocity dispersion ellipsoid. The
velocity dispersion along any direction can be found by generalising Equation
45. Defining two orthogonal directions within the sky plane ($\equiv$ the
plane orthogonal to $\widehat{\bm{r}}$), the angle $\psi$ with respect to the
first then defines any direction within the sky plane at the geometric centre
of the WB. By varying $\psi$ using 1D gradient descent, we will find the
direction $\widehat{\bm{\sigma}}_{\text{min}}$ along which the velocity
dispersion is minimal, with value $\sigma_{\text{min}}$. The direction
$\psi\pm\mathrm{\pi}/2$ is denoted $\widehat{\bm{\sigma}}_{\text{max}}$, with
associated velocity dispersion $\sigma_{\text{max}}$ also found using Equation
45. The relative probability of Galactic component $j$ based on the observed
position and $\bm{v}_{t}$ is then

$\displaystyle
P_{j}\propto\rho_{j}\exp\left(-\frac{1}{2}\overbrace{\sum_{i}\left(\frac{\Delta\bm{v}_{t}\cdot\widehat{\bm{\sigma}}_{i}}{\sigma_{i}}\right)^{2}}^{{\chi_{j}}^{2}}\right)\,,$
(47)

where $i$ takes two values covering the directions
$\widehat{\bm{\sigma}}_{\text{min}}$ and $\widehat{\bm{\sigma}}_{\text{max}}$,
while $\Delta\bm{v}_{t}$ is the difference between the observed $\bm{v}_{t}$
and the central expectation for Galactic component $j$. The Gaussian factor of
$2\mathrm{\pi}\sigma_{\text{min}}\sigma_{\text{max}}$ has been omitted from
the denominator as only relative probabilities are considered.

Once results are available for both Galactic disk components, the $P_{j}$ will
be normalised to a sum of 1. The WB will be rejected if it is inconsistent
with disk kinematics in the sense that the minimum of the two ${\chi_{j}}^{2}$
exceeds some threshold enclosing e.g. 99% of the total probability of a 2D
Gaussian. A logical flag will be used to record if this is the case.

For each WB, we will need to store whether it can reliably be considered part
of either Galactic disk, the two $P_{j}$, the mean and dispersion in $v_{r}$
if it belongs to Galactic disk component $j$, and the uncertainties in
$\widetilde{v}_{\text{sky}}$ and $\widetilde{v}_{\text{line}}$ (Section 6.3)
for quality control purposes only. These uncertainties will be simple
$1/\sqrt{N-1}$ standard deviations found by MC error propagation (Section
4.1). Note that since $\bm{v}_{t}$ varies between MC trials, the membership
probabilities to different disk components must be recalculated for each
trial.

##  Appendix F Inferring LOS distances

The method described in Banik (2019) will be used to infer the distribution of
the LOS separation of each WB, which in general is only known statistically
based on the sky-projected separation, which differs slightly between MC
trials. The actual separation of each WB along the LOS will be allocated using
an MC approach sampling its distribution. To avoid any artificial breaks in
power-law dependencies, the logarithmic slope of the $a$ distribution will be
fixed at the high-end slope of the WB population as published in Banik & Zhao
(2018b) $-$ this should be very accurate over the vast majority of the
relevant range. A broken power law can also be assumed, but for consistency
this should be the same as the $P\left(a\right)$ assumed for the WB
population, which is one of the model parameters to be fitted (Table 19). The
extra complexity means this is not envisaged in the nominal analysis.

##  Appendix G Pixel area corrections

The velocity transformation described in Section 5.8.2 involves binning in the
space of $\left(\widetilde{v}_{\perp},\widetilde{v}_{\parallel}\right)$ to
arrive at a total $\widetilde{v}$. This makes use of Pythagoras’ Theorem based
on the mid-point of each pixel. To improve the accuracy, we will first
determine the area of each bin in total $\widetilde{v}$ in the space of
$\left(\widetilde{v}_{\perp},\widetilde{v}_{\parallel}\right)$, which involves
finding the area of an annulus. This necessarily involves a factor of
$\mathrm{\pi}$. We will then find the total area of the pixels in
$\left(\widetilde{v}_{\perp},\widetilde{v}_{\parallel}\right)$ assigned to the
corresponding annulus, which involves simply counting pixels and multiplying
by the area of each pixel. Since $\mathrm{\pi}$ is irrational, there is a
small mismatch between these areas. For instance, the innermost annulus should
have an area of $\mathrm{\pi}$, but instead only 3 pixels with unit area are
assigned to this annulus. To approximately correct for this, we will determine
the ratio between the area of the annulus and the area of the pixels assigned
to it. These pixel area corrections will then be used to improve the accuracy
of the algorithm. The corrections become very small at high $\widetilde{v}$,
and are in any case never more than a few percent. These corrections will not
be applied to the DI array described in Section 5.8.2 since it does not rely
on annular binning.

##  Appendix H CBs with exactly equal mass

The WB population has an excess of systems where the stars have an exactly
equal mass (El-Badry et al., 2019). Regardless of why this is, it should be
considered in the WBT. Presumably, a similar excess is also present in the
undetected CB population.

Since $\Delta\widetilde{\bm{v}}=\bm{0}$ for a WB star contaminated by an equal
mass CB, we will use a special protocol based on the DI array to ensure that
absolutely no velocity perturbation is applied for the induced recoil effect
arising from equal mass CBs. The only variable on the CB side is therefore
$\Delta\widetilde{M}$, which cannot be neglected especially in the equal mass
CB case. Therefore, we will set up an array called $\bm{P}_{\text{0dv}}$
holding the likelihood of different $\Delta\widetilde{M}$ values arising from
the case that $q_{\text{int}}$ is exactly $1/2$. Thanks to the 1CB case, many
of the $\Delta\widetilde{M}$ bins receive a contribution as $q_{\text{ext}}$
can be $\ll 1/2$.

###  H.1 Distribution of the WB mass ratio

To have an affordable computational cost, the resolution in $q_{\text{ext}}$
will be 0.05, so 19 steps will be required to cover the range $0.025-0.975$.
Each pixel will be represented by a calculation at the central value of
$q_{\text{ext}}$. Therefore, a flat prior on $q_{\text{ext}}$ would normally
be imposed by setting all bins to have

$\displaystyle P\left(q_{\text{ext}}\right)~{}=~{}\frac{1}{19}\,.$ (48)

To allow an extra $\delta$-function in the distribution of $q_{\text{ext}}$ at
exactly $1/2$ such that the likelihood of an exactly equal mass WB is
$P_{\text{eqm}}$, Equation 48 will be adjusted to

$\displaystyle
P\left(q_{\text{ext}}\right)~{}=~{}\frac{1-P_{\text{eqm}}}{19}\,+\,P_{\text{eqm}}\delta\left(q_{\text{ext}}-\frac{1}{2}\right)\,.$
(49)

The $\delta$-function will be implemented by a boost to the probability in the
central bin, which for 19 pixels is bin 10.

If some slope is introduced on the $q_{\text{ext}}$ distribution, then the
idea will be for $P\left(q_{\text{ext}}\right)$ to rise linearly starting from
some non-zero $q_{\text{ext}}$ until $q_{\text{ext}}=1/2$, before decreasing
symmetrically down to zero such that there is exact symmetry between
$q_{\text{ext}}\to 1-q_{\text{ext}}$. The two triangular regions defined by
$P\left(q_{\text{ext}}\right)$ should have a total area of $1-P_{\text{eqm}}$,
since we also need to add an extra $\delta$-function at $1/2$. Preliminary
results suggest that it will be best to assume $P_{\text{eqm}}=0.02$ and no
WBs in our sample with $q_{\text{ext}}\leq 0.125$, as a constant slope over
the range $0.125-0.5$ provides a good fit.

In general, the basic principle is to combine a smooth and symmetric
$q_{\text{ext}}$ distribution with a $\delta$-function at $1/2$. The smooth
distribution should be integrated over the $q_{\text{ext}}$ range of each
pixel to find the probability in it, with the pixel containing
$q_{\text{ext}}=1/2$ receiving a further contribution of $P_{\text{eqm}}$.
Note that since we need the $q_{\text{ext}}$ distribution of Gaia WBs but the
CB contamination depends on the $q_{\text{int}}$ distribution of all CBs, we
expect to use different distributions for $q_{\text{int}}$ and
$q_{\text{ext}}$. Relative to the underlying WB population, there should be a
deficit of Gaia WBs where $q_{\text{ext}}$ differs greatly from $1/2$ because
a very faint WB companion might remain undetected.

###  H.2 The 1CB case

The $\bm{P}_{\text{0dv}}$ array is required mainly to handle the 1CB case, for
which $\Delta\widetilde{M}$ depends on how much of the total detected system
mass is in the contaminated star. This can be worked out by first obtaining
$x$ in Equation 16 for $q_{\text{int}}=1/2$, and then scaling it according to
Equation 18. Most bins in $\bm{P}_{\text{0dv}}$ will receive a contribution
once the full range of $q_{\text{ext}}$ is considered.

To improve the accuracy of this step, we will use a backward binning protocol
in which we find the range of $q_{\text{ext}}$ that yields a contribution to
each element of $\bm{P}_{\text{0dv}}$. By considering also the minimum and
maximum allowed $q_{\text{ext}}$ (Section H.1), we will find the overlap
range, and thus the integral of $P\left(q_{\text{ext}}\right)$ across this. A
special allowance will be required in case this overlap range includes
$q_{\text{ext}}=1/2$, since then we will need to add $P_{\text{eqm}}$. The
clause for this special case must arise for only one of the
$\bm{P}_{\text{0dv}}$ bins, providing an important cross check. Another check
is that before the MCMC stage, the $\bm{P}_{\text{0dv}}$ array should sum to 1
as factors related to $f_{\text{CB}}$ and $P_{\text{eqm}}$ will be put in as
part of each MCMC trial (Appendix H.4).

Since part of the 1CB probability will be handled through the
$\bm{P}_{\text{0dv}}$ array, the normalisation of the
$\textbf{{P}}_{\text{1CB}}$ array should be reduced before the MCMC stage. We
will arrange this by implementing

$\displaystyle\textbf{{P}}_{\text{1CB}}~{}\to~{}\left(1-P_{\text{eqm}}\right)\textbf{{P}}_{\text{1CB}}\,.$
(50)

It will become clear that a similar adjustment is not required for the 2CB
case because a separate library will be prepared rated for the desired
$P_{\text{eqm}}$ (Appendix H.3)

Note that the pixel with $q_{\text{int}}=1/2$ should receive some probability
even if a uniform prior is assumed on $q_{\text{int}}$. As a result, the value
of $P_{\text{eqm}}$ used here should exceed that in Appendix H.1, where the
true value should be used. Equation 50 should use a value which is higher by
$1/999$ because the values of $q_{\text{int}}$ assumed in Table 8 cover the
range $\left(0.5-499.5\right)\text{\textperthousand}$, so an extra pixel with
half as much weight should be added at $q_{\text{int}}=1/2$ if the goal is a
uniform prior over the range $0.0005-0.5$. Since this argument also applies to
the 2CB case and the scalings are implemented in the main algorithm, at this
point we will set

$\displaystyle P_{\text{eqm}}~{}\to~{}P_{\text{eqm}}+\frac{1}{999}\,.$ (51)

The $\bm{P}_{\text{0dv}}$ array will then be scaled to the desired
$f_{\text{CB}}$ and $P_{\text{eqm}}$ as discussed in Section H.4.

###  H.3 The 2CB case

For the 2CB case, the total $\Delta\widetilde{\bm{v}}$ can only be assumed to
vanish if $q_{\text{int}}=1/2$ for both CBs. If so, only the largest
$\Delta\widetilde{M}$ bin would receive a contribution. While the likelihood
may seem small, a significant amount of hidden mass with no induced recoil on
either WB star could mimic the effect of a gravity law stronger than
Newtonian. This is one major motivation to handle the equal mass CB case with
some care. In the 2CB case where both CBs have an exactly equal mass, the
contribution to $\bm{P}_{\text{0dv}}$ is
${P_{\text{eqm}}}^{2}{f_{\text{CB}}}^{2}$ into the bin with the highest
$\Delta\widetilde{M}$. Since $f_{\text{CB}}$ is a model parameter (Table 19),
this contribution must be applied during the MCMC stage. To optimise this,
${P_{\text{eqm}}}^{2}$ will be determined beforehand, but after Equation 51
has been applied.

It is possible that only one of the two CBs is an exactly equal mass CB. To
handle this case, we will set up an extension to the standard
$\textbf{{P}}_{\text{2CB}}$ library, with the extension denoted
$\textbf{{P}}_{\text{2CB,e}}$ in what follows. This will hold the CB
contamination pattern for each $q_{\text{ext}}$ on the assumption that star 1
is contaminated by an exactly equal mass CB, and that it comprises a fraction
$q_{\text{ext}}$ of the detected WB mass. Star 2 is assumed to be contaminated
by a CB with a uniform $q_{\text{int}}$ distribution, so the CBT will be used
(Section 3.3). Since there is an asymmetry between $q_{\text{ext}}\to
1-q_{\text{ext}}$, the library cannot be restricted to $q_{\text{ext}}\geq
1/2$, and will thus have a similar size to the 1CB library (Table 11). The
computational cost should also be similar, and thus vastly smaller than for
the 2CB library (Section 3.5).

The $\textbf{{P}}_{\text{2CB,e}}$ library will be used to create a revised
$\textbf{{P}}_{\text{2CB}}$ library where none or one of the CBs has
$q_{\text{int}}=1/2$. This is because the case where both CBs have
$q_{\text{int}}=1/2$ leads to no photocentre-barycentre offset, so we will
handle it through the $\bm{P}_{\text{0dv}}$ array. Not considering this case
means that the normalisation of the $\textbf{{P}}_{\text{2CB}}$ array will
drop by a fraction close to ${P_{\text{eqm}}}^{2}$.

The revised $\textbf{{P}}_{\text{2CB}}$ library for each $q_{\text{ext}}$ will
be found using

$\displaystyle\textbf{{P}}_{\text{2CB}}\left(q_{\text{ext}}\right)~{}\to~{}\left(1-P_{\text{eqm}}\right)^{2}\textbf{{P}}_{\text{2CB}}\left(q_{\text{ext}}\right)$
(52)
$\displaystyle+P_{\text{eqm}}\left(1-P_{\text{eqm}}\right)\left[\textbf{{P}}_{\text{2CB,e}}\left(q_{\text{ext}}\right)+\textbf{{P}}_{\text{2CB,e}}\left(1-q_{\text{ext}}\right)\right]\,.$

The basic principle is that $q_{\text{ext}}$ tells us the WB mass fraction in
star 1, not whether star 1 or star 2 is contaminated by an equal mass CB. Both
options are possible, albeit not simultaneously $-$ as explained above, this
possibility is handled later. We will safeguard the original
$\textbf{{P}}_{\text{2CB}}$ library by setting up another version rated for
the desired value of $P_{\text{eqm}}$. Note that this is required even for a
uniform prior on $q_{\text{int}}$ (Equation 51).

###  H.4 Adjustments to the main algorithm

In the main part of the MCMC analysis, the $\bm{P}_{\text{0dv}}$ array must be
scaled to the desired $f_{\text{CB}}$ and $P_{\text{eqm}}$. A separate copy of
the array will be required during the MCMC stage so that the original is still
available after each MCMC trial. Since the original $\bm{P}_{\text{0dv}}$
array is entirely about the 1CB case (Section H.2), it must be scaled by
$2P_{\text{eqm}}f_{\text{CB}}\left(1-f_{\text{CB}}\right)$. To handle the 2CB
case, ${P_{\text{eqm}}}^{2}{f_{\text{CB}}}^{2}$ must then be added to the bin
with the highest $\Delta\widetilde{M}$ (Appendix H.3).

To include probability contributions from cases with
$\Delta\widetilde{\bm{v}}=\bm{0}$, the DI array will receive contributions not
only from the ‘raw’ WB library to handle the 0CB case, but also a version for
every $\Delta\widetilde{M}$ where a stretch is applied in the velocity for the
hidden mass effect (Equation 19). Each time, the contribution will be scaled
by the corresponding element in the $\bm{P}_{\text{0dv}}$ array. As with the
nominal analysis, the DI array will then be directly injected into the
simulated photo on the last $\Delta\widetilde{M}$ index, by which time the DI
array must have received all the contributions it is supposed to receive. Note
that the range of velocity indices must be that determined for the highest
$\Delta\widetilde{M}$, since the range also scales according to Equation 19.


# Comparing Implementations of Self-Interacting Dark Matter in the Gizmo and
Arepo Codes

Helen Meskhidze,1 Francisco J. Mercado,2 Omid Sameie,3 Victor H. Robles,4
James S. Bullock,2 Manoj Kaplinghat,2 James O. Weatherall1
1Department of Logic and Philosophy of Science, 779 Social Science Tower,
University of California Irvine, CA 92697, USA
2Center for Cosmology, Department of Physics and Astronomy, 4129 Reines Hall,
University of California Irvine, CA 92697, USA
3Department of Astronomy, The University of Texas Austin, 2515 Speedway, Stop
C1400, Austin, TX 78712 USA
4Physics Department, Yale Center for Astronomy and Astrophysics, New Haven, CT
06520, USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Self-interacting dark matter (SIDM) models have received great attention over
the past decade as solutions to the small-scale puzzles of astrophysics.
Though there are different implementations of dark matter (DM) self-
interactions in N-body codes of structure formation, there has not been a
systematic study to compare the predictions of these different
implementations. We investigate the implementation of dark matter self-
interactions in two simulation codes: Gizmo and Arepo. We begin with identical
initial conditions for an isolated $10^{10}$ M⊙ dark matter halo and
investigate the evolution of the density and velocity dispersion profiles in
Gizmo and Arepo for SIDM cross-section over mass of 1, 5, and 50 $\rm
cm^{2}g^{-1}$. Our tests are restricted to the core expansion phase where the
core density decreases and core radius increases with time. We find better
than 30% agreement between the codes for the density profile in this phase of
evolution, with the agreement improving at higher resolution. We find that
varying code-specific SIDM parameters changes the central halo density by less
than 10% outside of the convergence radius. We argue that SIDM core formation
is robust across the two different schemes and conclude that these codes can
reliably differentiate between cross-sections of 1, 5, and 50 $\rm
cm^{2}g^{-1}$ but finer distinctions would require further investigation.

###### keywords:

dark matter– software: simulations – cosmology: theory

††pubyear: 2022††pagerange: Comparing Implementations of Self-Interacting Dark
Matter in the Gizmo and Arepo Codes–Comparing Implementations of Self-
Interacting Dark Matter in the Gizmo and Arepo Codes

##  1 Introduction

Self-interacting dark matter (SIDM) is a generic prediction of dark sector
models for physics beyond the Standard Model (Spergel & Steinhardt, 2000; Ahn
& Shapiro, 2005; Ackerman et al., 2009; Arkani-Hamed et al., 2009; Feng et
al., 2009; Loeb & Weiner, 2011; Tulin et al., 2013), and it is a possible
explanation for small-scale structure formation puzzles (for a comprehensive
review of small-scale challenges, see Bullock & Boylan-Kolchin, 2017). Several
groups have confirmed that SIDM with a cross-section over mass of order
$1~{}\rm cm^{2}g^{-1}$ or larger can alleviate small-scale issues (Davé et
al., 2001; Colín et al., 2002; Vogelsberger et al., 2012; Rocha et al., 2013;
Kamada et al., 2017). This is due to the fact that scattering can effectively
transfer kinetic energy within galactic halos and change the dark matter
distribution (see Tulin & Yu, 2018, for a review of SIDM phenomenology). In
particular, the changes introduced by the heat transfer have been shown to
provide an economical way to explain the diverse range of rotation curves of
galaxies (Oman et al., 2015; Ren et al., 2019a) and the “too-big-to-fail"
problem (Vogelsberger et al., 2012; Kaplinghat et al., 2019; Turner et al.,
2021).

Simulations of SIDM with baryons have found that while baryonic feedback can
reduce the central density of a cuspy halo, if an SIDM halo already has a core
present, the feedback does not make a significant difference (Robles et al.,
2017; Fitts et al., 2019; Sameie et al., 2021). This result suggests that SIDM
predictions are fairly robust to feedback implementation and provides further
motivation to use the observed properties of galaxies to test SIDM models.

The growing use of simulation-derived predictions to constrain the
microphysics of SIDM motivates us to compare predictions from different codes.
Since it is impossible to model DM-DM particle scattering directly in a galaxy
formation simulation, the macroscopic effects must be modeled in an
approximate way. As discussed by Tulin & Yu (2018), there are various methods
for implementing DM self-interactions but the differences that may arise from
each have yet to be studied in detail. Here, for the first time, we present a
code comparison of two implementations of simple elastic SIDM, specifically
focusing on the popular Gizmo and Arepo codes. We begin with identical initial
conditions for an isolated $10^{10}$ M⊙ dwarf halo, which is a mass regime of
particular interest for small-scale structure tests. We investigate the
effects of the SIDM implementations within and between the codes by varying
the SIDM cross-section per mass $\sigma$/m = 1, 5, and 50 cm2 g-1.

Our work is structured as follows: section 2 presents our initial conditions
in more detail and outlines the SIDM implementations in each code. Section 3
presents the effects of changing the SIDM parameters, the concentration of the
halo, and the resolution of the simulations. Section 4 outlines our
conclusions.

##  2 The Simulations

###  2.1 Code Descriptions

In the following, we introduce our two simulation codes – Gizmo and Arepo –
and describe the methods each uses for implementing DM self-interactions.

Halo name Mass [M⊙] DM particles Run time ($\alpha,\beta,\gamma$) Scale radius
Force softening Smoothing length Neighbours searched $\sigma$/m Colour $N_{p}$
[Gyr] $r_{s}$ [kpc] $\epsilon$ [pc] $h_{si}$ (Gizmo) $k$ (Arepo) [
$\text{cm}^{2}\,\text{g}^{-1}$] m10HR $10^{10}$ $5\times 10^{6}$ 10
(1,3,1.536) 14 10 – – CDM (High Resolution) 0.25 $\epsilon$ 32 $\pm$ 5 1 0.25
$\epsilon$ 32 $\pm$ 5 5 0.25 $\epsilon$ 32 $\pm$ 5 50 m10 $10^{10}$ $1\times
10^{6}$ 10 (1,3,1.536) 14 10 – – CDM (Baseline; 0.25 $\epsilon$ 32 $\pm$ 5 1
Fiducial Resolution) 0.25 $\epsilon$ 32 $\pm$ 5 5 0.25 $\epsilon$ 32 $\pm$ 5
50 m10SIDM- $10^{10}$ $1\times 10^{6}$ 10 (1,3,1.536) 14 10 0.125 $\epsilon$
16 $\pm$ 5 5 m10SIDM+ 0.5$\epsilon$ 64 $\pm$ 5 5

Table 1: Global parameters of the halos. All simulations described above have
been carried out in Gizmo and Arepo.

####  2.1.1 Gizmo

Gizmo is a massively parallel, multi-physics simulation code that uses a
meshless Lagrangian Godunov-type method (“Meshless finite-mass” or MFM;
Hopkins 2015; Hopkins et al. 2018).111Gizmo, including the SIDM module, is
publicly available here: https://bitbucket.org/phopkins/gizmo-public. Given
that our simulations are DM only, we rely on only the N-body tree-gravity
solver which is derived from Gadget-3 (Springel, 2010).

Gizmo’s implementation of elastic self-interactions uses the methodology
introduced by Rocha et al. (2013), which is based on the rate of scattering of
the DM macro-particles in phase space. The probability of an interaction is
calculated as:

$P_{ij}=(\sigma/m)m_{i}v_{ij}g_{ij}\delta t,$ (1)

where $m_{i}$ is the mass of the macroparticle, $v_{ij}$ is the relative
velocity difference between the two macroparticles, and $g_{ij}$ is the number
density factor that accounts for the overlap of the two macroparticles’
smoothing kernels. A random number is drawn from a uniform distribution in the
interval $(0,1)$ to determine whether an interaction takes place. If an
interaction takes place, the particles are given outgoing velocities
consistent with elastic scattering. The outgoing velocities are calculated in
terms of the center-of-mass velocity of the two particles, their masses, and
their relative speed. The direction of the scatter is randomly chosen (such
that the scatter is isotropic in the center-of-mass frame) and each particle
moves opposite the other (for details and tests of this implementation against
analytic problems see Rocha et al., 2013).

####  2.1.2 Arepo

Arepo is a massively parallel, multi-physics simulation code that employs a
finite-volume method on a moving Voronoi mesh and a tree-particle-mesh method
for gravitational interactions. Details of the underlying method can be found
in (Springel, 2010) while the most recent release of the code is described in
(Weinberger et al., 2020).222Arepo has recently been publicly released (see
Weinberger et al. 2020) and is available here:
https://gitlab.mpcdf.mpg.de/vrs/arepo. Note, however, that the public release
does not include the SIDM implementation.

Arepo estimates the probability of an elastic self-interaction at each time
step by calculating the scattering probability for each particle $i$ with each
of its $k$ nearest neighbours (32 $\pm$ 5 by default and in our baseline
model) $j$ as:

$P_{ij}=(\sigma/m)m_{i}v_{ij}W(r_{ij},h_{i})\delta t,$ (2)

where $W(r_{ij},h_{i})$ is the cubic spline Kernel function and $h_{i}$ is the
smoothing length enclosing the $k$ nearest neighbours of particle $i$ (for
details of this SIDM implementation, see Vogelsberger et al., 2012).

Once the probability of interaction is calculated, a random number is drawn
from a uniform distribution in the interval $(0,1)$ to determine if an
interaction takes place. If it does, a neighbor $j$ must be selected to
scatter with. The set of neighbors is sorted by distance to the original
particle $i$ and the first neighbor $l$ whose pairwise probability to scatter
with the original particle satisfies the inequality $x\leq\sum_{i}^{l}P_{ij}$
is chosen. Once a pair is matched, each particle is given a new velocity that
reflects the original center-of-mass velocity and the two particles’ relative
velocity. The direction is chosen randomly but each particle moves opposite
the other.

Given the different approaches taken by each code in its underlying gravity
solver, the similarities between their SIDM calculations and implementations
are striking. With these similarities in mind, we ask: what differences, if
any, arise between the two codes?

###  2.2 Initial Conditions

Figure 1: The density profile and velocity dispersion profile generated from
SpherIC for the high and fiducial resolution simulations ($5\times 10^{6}$ and
$1\times 10^{6}$ particles respectively).

We use SpherIC, an initial conditions generator for spherically symmetric
systems in equilibrium first presented in Garrison-Kimmel et al. (2013).
SpherIC is based on HALOGEN4MUSE (Zemp et al., 2008)333SpherIC was available
on Bitbucket, but since that site stopped supporting Mercurial repositories,
it is no longer publicly available. However HALOGEN4MUSE is available here:
https://github.com/mzemp/halogen. Our initial conditions are also available
upon request. and generates profiles in the $(\alpha,\beta,\gamma)$-model
family, a generalization of the Navarro-Frenk-White (NFW) model:

$\rho(r)=\frac{\rho_{s}}{(\frac{r}{r_{s}})^{\gamma}[1+(\frac{r}{r_{s}})^{\alpha}]^{(\beta-\gamma)/\alpha}}\,,$
(3)

where $r_{s}$ is the scale radius and $\rho_{s}$ is the scale density (Navarro
et al., 1997; Navarro et al., 2010). The parameters $\gamma$ and $\beta$
characterize the inner and outer power-law slope of the halo, respectively.
The quantity $\alpha$ determines the sharpness of transition between the inner
and outer slope. An NFW profile corresponds to $(\alpha,\beta,\gamma)=(1,3,1)$
and the scale radius in this case is equal to the radius where the logarithmic
slope of the density profile is -2 (i.e., $r_{s}=r_{-2}$). For more general
($\alpha,\beta,\gamma$)-models, we have the following relation (Di Cintio et
al., 2014):

$r_{-2}=\left(\frac{2-\gamma}{\beta-2}\right)^{1/\alpha}r_{s}\,.$ (4)

In this work we utilize a halo model with $(\alpha,\beta,\gamma)=(1,3,1.536)$,
$\rho_{s}=1.998\times 10^{5}\mathrm{\ M_{\odot}kpc^{-3}}$, and $r_{s}=14.14$
kpc as these parameter values lie within the ranges that we would expect for a
typical halo of this mass (see Lazar et al., 2020). We set our halos to have a
total mass of $10^{10}M_{\odot}$. The virial radius (defined as the radius
within which the average density is 100$\rho_{\textrm{crit}}$) is 52.5 kpc and
the mass enclosed is $7.61\times 10^{9}M_{\odot}$. The values we use for
$(\alpha,\beta,\gamma)$ correspond to $r_{-2}=6.56$ kpc - where $r_{-2}$ is
the radius at which the profile has a slope of $-2$.

We use the parameters described above to create initial conditions at two
different resolution levels: a baseline (fiducial) resolution with $10^{6}$
particles and a higher resolution of $5\times 10^{6}$ particles. We show the
radial density and velocity dispersion profiles of these two initial
conditions in Fig. 1.

###  2.3 Runs

We evolve the initial conditions at two resolutions in Gizmo and Arepo. For
each resolution, we run one CDM simulation and three SIDM simulations. The
results for the CDM models are shown in black in the figures throughout this
work. For the code-to-code SIDM implementation comparison, we evolve both
initial conditions using the respective SIDM implementations of Gizmo and
Arepo for 3 different cross-sections: $\sigma$/m $=1,5,50$ cm2 g-1. For the
Gizmo simulations, we adopt the default SIDM smoothing length ($h_{si}$) of
25% of the force softening. Since the force softening we use is 10 pc, our
smoothing length is 2.5 pc. Likewise, for the Arepo simulations, we adopt the
default value for the neighbors searched ($k$) which is 32 $\pm$ 5.

Finally, we test the effects of varying code-specific SIDM parameters at the
fiducial resolution. We set the SIDM cross-section to 5 cm2 g-1 and vary the
smoothing length (in Gizmo) and the number of neighbours searched (in Arepo).
For Gizmo we adopt smoothing lengths of 1.25 pc and 5 pc which we refer to as
SIDM- and SIDM+ respectively. For Arepo we set the number of neighbours
searched to 16 $\pm$ 5 and 64 $\pm$ 5 and refer to these again as SIDM- and
SIDM+. Note that these values are less than and greater than the default
values set for these parameters in the respective codes. In total, we present
a suite of 20 simulations (see Table 1 for a detailed list of all the
simulations).

##  3 Results

Figure 2: A resolution comparison of the two codes. The left panel depicts m10
high resolution and m10 fiducial resolution in Gizmo evolved to 2, 5, and 10
Gyr. The right panel depicts the same for Arepo. The vertical dotted line
represents the radius at which the majority of our simulations are converged
($300$ pc).

###  3.1 Convergence of each code

Before comparing the results from the two codes, we first discuss the
sensitivity of each code internally to resolution. In Fig. 2, we present 3
snapshots (2, 5, and 10 Gyr) of each code (Gizmo in the left panel and Arepo
in the right). The figure shows the results of the high and fiducial
resolution simulations for CDM in black and SIDM $\sigma$/m $=1,5$ and $50$
cm2 g-1 in dark blue, light blue, and green respectively.

All of our high resolution simulations contain at least 200 particles within a
radius of 250 pc throughout the 10 Gyr simulations. Nearly all the fiducial
resolution simulations contain more than 200 particles within 300 pc.444The
only exception is the fiducial resolution Arepo simulation with $\sigma$/m
$=50$ cm2g-1 at late times which contains 200 particles at 400 pc. For
consistency, we mark this radius on all our plots with a black dotted line and
conduct all our comparisons at 300 pc (e.g., the comparison done in Fig. 5).

The error plotted in Fig. 2 and on all plots in this paper is the Poisson
error (calculated as the density at the bin divided by the square root of the
number of particles in the bin) and is represented with the shaded region
around the curve. Some of the profiles are cut off because the simulations do
not have any particles in the innermost region of the halo.

In sum, the simulations plotted in Fig. 2 demonstrate remarkable agreement
between the fiducial and high resolution simulations. Indeed, outside of 250
pc, nearly all the simulations are within one another’s error margins (the
only exception being the $50$ cm2 g-1 run in Gizmo).

###  3.2 Comparison of the codes across resolutions

Having tested the convergence of the codes internally as the resolution is
increased, we next compared the two resolutions across the two codes. For this
comparison, we choose the most extreme SIDM cross-section ($\sigma$/m = 50
cm2g-1). Fig. 3 shows the results of Gizmo and Arepo for this cross-section at
2, 5, and 10 Gyr. In this figure, we plot the fiducial resolution simulations
in green and the high resolution in black. The Gizmo simulations are plotted
with solid lines while the Arepo simulations are plotted with dashed lines. We
find that increasing the resolution of the simulations brings the simulation
results into better agreement. In other words, the density profiles from Gizmo
and Arepo exhibit better agreement at the higher resolution than at the
fiducial resolution.

Figure 3: From top to bottom, these figures show the halo density profiles at
2, 5, and 10 Gyr for a high resolution halo with $5\times 10^{6}$ particles
and $\sigma$/m = 50 cm2g-1.

###  3.3 Sensitivity to SIDM Cross-section

Figure 4: The top panel shows a comparison of the high resolution simulations
of Gizmo (solid) and Arepo (dashed) evolved to 10 Gyr at various SIDM cross-
sections. The bottom panel shows the same for the fiducial resolution
simulation

Next, we consider the sensitivity of each code to different SIDM cross-
sections at high and fiducial resolution. We compare our results for the SIDM
cross-section $\sigma/m$ = 1, 5, and 50 $\text{cm}^{2}\,\text{g}^{-1}$ in Fig.
4. We have plotted the Gizmo simulations with solid lines and the Arepo
simulations with dashed lines and use the same colour conventions as in Fig.
2.

Consistent with Rocha et al. (2013), Elbert et al. (2015), and Fitts et al.
(2019), we found that an increase in the SIDM cross-section results in density
profiles that are shallower and more cored. One can better understand why the
density profiles become more cored by looking at the velocity dispersion
profiles. As Rocha et al. (2013) describe, the core is created by the heat
transport (characterized by the DM velocity dispersion) from large radii to
the inner region. The velocity dispersion curves flattened as the halo was
evolved in our simulations (see the right panel of Figure 4). The flat
velocity dispersion profiles of the SIDM simulations indicate that the SIDM
halos are thermalized within the core, which is a necessary condition for
establishing a cored DM density profile.

Analytic models of gravothermal evolution predict that halo evolution in the
long-mean-free path is self-similar and determined by the product of
scattering time and age (Koda & Shapiro, 2011; Nishikawa et al., 2020). Given
previous work (see, e.g., Ren et al. (2019b) and Robertson et al. (2021)), we
expect the central density for a given halo to be a function of age times
$\sigma/m$ since the scattering rate is proportional to $\sigma/m$. We can see
this behavior clearly in the left panel of Fig. 5 where we have plotted
$\rho_{300}$ (proxy for the core density) vs. age multiplied by $\sigma/m$.
Both Arepo and Gizmo follow the general trend of decreasing core density
(increasing core size) with increasing age times $\sigma/m$. Note that the
initial lack of evolution in $\rho_{300}$ is because the core size is smaller
than $300~{}\rm pc$ and the density at $300~{}\rm pc$ is close to its initial
value.

Fig. 5 provides a generally encouraging picture of agreement between Arepo and
Gizmo. It is worth noting that the Arepo results seem to provide a more
seamless curve when the three different cross-sections are plotted together as
in the left panel of Fig. 5, as expected from analytic models.

To better compare the results of the two codes, we also tracked the difference
between the density profiles at $r=300$ pc through 10 Gyr. We chose $r=300$ pc
since our resolution tests suggest the simulations are well-converged at this
radius (see §3.1 for more). For the fiducial resolution simulations, the
difference in density profiles was approximately 20-40% and the velocity
dispersion profiles were within 5-10% of one another. For the high resolution
simulations, the halo profiles were largely within 30% of one another (see the
right panel of Fig. 5) and the velocity dispersion profiles were within 5%.

We find that the difference between the two codes is much smaller than the
difference between the various SIDM cross-sections we tested. In other words,
we are confident these simulations can be used to distinguish amongst SIDM
cross-sections of 1, 5, and 50 cm2 g-1. However, the code-to-code variation is
large enough that we would argue against using such simulations to
differentiate between effects due to SIDM cross-sections of, e.g., 1 vs. 1.5
cm2 g-1.

Finally, as seen in both Fig. 4 and the right panel of Fig. 5, there is an
inversion of the DM density profiles, computed by the Arepo and Gizmo codes,
for the SIDM cross-section $\sigma$/m $=50$ cm2g-1. For cross-sections of
$\sigma$/m $=1,5$ cm2g-1, the halos simulated with Arepo have higher densities
than those simulated with Gizmo. However, around 4-5 Gyr, the halos with
cross-sections of $50$ cm2g-1 evolved with Gizmo become denser than the those
evolved with Arepo. The inversion is more obvious in the fiducial resolution
simulation (bottom panel of Fig. 4) but is also seen in the higher resolution
simulations (top panel of Fig. 4 and right panel of Fig. 5).

Figure 5: For the left panel, we plot the core density (at 300 pc) of each
cross-section simulated as a function of the time of the simulation multiplied
by the cross-section itself. For the right panel, we show the difference
between the Gizmo and Arepo high resolution density profiles at $r=300$ pc
again tracked through 10 Gyr.

To better understand the results of each code, we also tracked the number of
self-interactions amongst the dark matter particles. As expected, the number
of self-interactions increases with increasing cross-section. The scaling is
nearly linear in both codes, with the interaction cross-section of 50 cm2g-1
exhibiting 8 times as many interactions as the $5$ cm2g-1 simulation and the
$5$ cm2g-1 exhibiting about 4 times as many self-interactions as the $1$
cm2g-1 simulation. We additionally find that the differences in the number of
DM self-interactions per time step are set early on in the simulation and stay
consistent throughout a 10 Gyr run. For example, the Gizmo $1$ cm2g-1 run
consistently exhibits about 10,000 DM self-interactions per Gyr throughout the
entire 10 Gyr simulation and the Arepo $1$ cm2g-1 run consistently exhibits
about 7,000 DM self-interactions per Gyr. Likewise, the Gizmo $5$ cm2g-1 run
consistently has around 40,000 DM self-interactions per Gyr, and the
corresponding Arepo run has around 30,000.

When comparing the self-interactions between the two codes, we find that for
$\sigma$/m $=1,5$ cm2g-1, Gizmo exhibits a greater number of DM self-
interactions ($\sim$25% more at each time step and overall), which results in
the density profiles being more cored (i.e., less dense in the inner regions).
For $\sigma$/m $=50$ cm2g-1, Gizmo begins with slightly more self-interactions
than Arepo ($\sim$5% more). However, consistent with the above discussion,
there is an inversion in this trend between 3 and 5 Gyr where the number of
self-interactions in Arepo overtakes Gizmo. By 6 Gyr, Arepo exhibits $\sim
8\%$ more self-interactions than Gizmo which is maintained across the 10 Gyr.

###  3.4 Sensitivity to code-specific SIDM parameters

We also tested the sensitivity of each code to various code-specific SIDM
parameters. In Figure 6, we plot the difference between the baseline halo
density at 5 Gyr and the result of varying the smoothing factor and the number
of neighbours in Gizmo and Arepo respectively. For Gizmo, a smoothing factor
of 25% of the force softening is the default and has been used in, e.g.,
Elbert et al. (2015). For Arepo, the default number of neighbours searched is
32 $\pm$ 5 but, e.g. Despali et al. (2019) search 64 of the nearest
neighbours. We therefore check the sensitivity of both codes to increasing and
decreasing the default values for these parameters by a factor of 2. As seen
in that figure, these code-specific SIDM choices make up to 10% difference to
the halo density at our innermost converged radius ($\sim$300 pc). However,
there is no general trend with respect to the halo density apparent when the
values are increased or decreased.

Figure 6: Sensitivity of the halo density to code-specific SIDM parameters at
5 Gyr and $r=300$ pc for a fiducial resolution halo and $\sigma$/m = 5 cm2g-1.
We vary the smoothing factor in Gizmo with SIDM- corresponding to a smoothing
factor of 0.125 pc and SIDM+ corresponding to a smoothing factor of 0.5 pc.
Correspondingly, in Arepo, we vary the number of neighbours searched. The
SIDM- simulation corresponds to 16$\pm$5 and the SIDM+ simulation corresponds
to 64$\pm$5\. We then plot the difference between each of these simulations
and the fiducial simulation for each code (e.g., the dark blue solid line
corresponds to (Gizmo SIDM-) - (Gizmo fiducial))/(Gizmo fiducial).

##  4 Discussion and Conclusion

Using an N-body simulation suite of an isolated dwarf dark matter halo, we
compare the SIDM implementations of two simulation codes: Gizmo and Arepo. We
use SpherIC to generate an isolated halo with halo mass of $10^{10}\,{\rm
M}_{\odot}$ as our initial conditions. We adopt constant elastic SIDM cross-
sections $\sigma$/m $=1,5,{\rm\ and\ }50$ $\text{cm}^{2}\,\text{g}^{-1}$
throughout our analysis. Our main conclusions are summarized in the following.

We find that the core density of the halo tracks the number of self-
interactions predicted by the code. We find that Gizmo predicts a greater
number of self-interactions among the DM particles overall and that the
density profiles are generally less dense in the inner region than Arepo.
These general results are consistent with previous results that have found
that DM self-interactions make the halo density profile more cored (see, e.g.,
Rocha et al. 2013; Elbert et al. 2015).

Increasing the resolution of the runs (from the baseline of $1\times 10^{6}$
to $5\times 10^{6}$ DM particles) brings the codes into better agreement. At
this increased resolution, the two codes predict densities that are similar,
to within 10-30%, with Arepo consistently slightly denser than Gizmo. Code-
specific SIDM parameters made less than a 10% difference to the final density.
Generally, it seems that one can reliably predict differences between SIDM
cross-sections of $\sigma$/m = 1, 5, and 50 cm2g-1 based on the agreement
between these two codes. However, the inferred 30% difference in density is
large enough to preclude any finer-grained distinctions (e.g., between results
obtained with with $\sigma$/m = 1.0 vs. 1.5 cm2g-1). Measuring the densities
of halos may pose a larger issue than the uncertainties identified in this
work.

We find that our limited scope, i.e. only comparing the SIDM implementations
in each code, enables us to conduct a fruitful code comparison. In particular,
instead of requiring a large comparison infrastructure and network of
collaborators, we were able to isolate and discuss the effects of SIDM in the
simulations. Naturally, this comparison can be extended to other codes
implementing SIDM as well as to include baryonic physics and/or inelastic
interactions. Having this underlying understanding of the SIDM module will be
invaluable in turning to these larger projects.

## Acknowledgements

The authors all thank the anonymous referee for their careful review of this
manuscript. This article is partially based on work done while HM and FJM were
graduate student researchers under the John Templeton Foundation grant “New
Directions in Philosophy of Cosmology” (grant no. 61048). JSB and JOW also
acknowledge partial support from the same grant. The opinions expressed here
are those of the authors and not necessarily those of the John Templeton
Foundation. FJM and JSB were also supported by National Science Foundation
(NSF) grant AST-1910965. VHR acknowledges support by the Yale Center for
Astronomy and Astrophysics Prize postdoctoral fellowship.

## Data Availability

The data supporting the plots within this article are available on reasonable
request to the corresponding author. A public version of the Gizmo code is
available at http://www.tapir.caltech.edu/ phopkins/Site/GIZMO.html. A public
version of the Arepo code is available at
https://gitlab.mpcdf.mpg.de/vrs/arepo but we obtained a copy with the SIDM
module by directly contacting the code developers.

## References

  * Ackerman et al. (2009) Ackerman L., Buckley M. R., Carroll S. M., Kamionkowski M., 2009, Phys. Rev. D, 79, 023519
  * Ahn & Shapiro (2005) Ahn K., Shapiro P. R., 2005, MNRAS, 363, 1092
  * Arkani-Hamed et al. (2009) Arkani-Hamed N., Finkbeiner D. P., Slatyer T. R., Weiner N., 2009, Phys. Rev. D, 79, 015014
  * Bullock & Boylan-Kolchin (2017) Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343
  * Colín et al. (2002) Colín P., Avila-Reese V., Valenzuela O., Firmani C., 2002, ApJ, 581, 777
  * Davé et al. (2001) Davé R., Spergel D. N., Steinhardt P. J., Wandelt B. D., 2001, ApJ, 547, 574
  * Despali et al. (2019) Despali G., Sparre M., Vegetti S., Vogelsberger M., Zavala J., Marinacci F., 2019, MNRAS, 484, 4563
  * Di Cintio et al. (2014) Di Cintio A., Brook C. B., Dutton A. A., Macciò A. V., Stinson G. S., Knebe A., 2014, MNRAS, 441, 2986
  * Elbert et al. (2015) Elbert O. D., Bullock J. S., Garrison-Kimmel S., Rocha M., Oñorbe J., Peter A. H. G., 2015, MNRAS, 453, 29
  * Feng et al. (2009) Feng J. L., Kaplinghat M., Tu H., Yu H.-B., 2009, J. Cosmology Astropart. Phys., 2009, 004
  * Fitts et al. (2019) Fitts A., et al., 2019, MNRAS, 490, 962
  * Garrison-Kimmel et al. (2013) Garrison-Kimmel S., Rocha M., Boylan-Kolchin M., Bullock J. S., Lally J., 2013, MNRAS, 433, 3539
  * Hopkins (2015) Hopkins P. F., 2015, MNRAS, 450, 53
  * Hopkins et al. (2018) Hopkins P. F., et al., 2018, MNRAS, 480, 800
  * Kamada et al. (2017) Kamada A., Kaplinghat M., Pace A. B., Yu H.-B., 2017, Phys. Rev. Lett., 119, 111102
  * Kaplinghat et al. (2019) Kaplinghat M., Valli M., Yu H.-B., 2019, MNRAS, 490, 231
  * Koda & Shapiro (2011) Koda J., Shapiro P. R., 2011, MNRAS, 415, 1125
  * Lazar et al. (2020) Lazar A., et al., 2020, MNRAS, 497, 2393
  * Loeb & Weiner (2011) Loeb A., Weiner N., 2011, Phys. Rev. Lett., 106, 171302
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Navarro et al. (2010) Navarro J. F., et al., 2010, MNRAS, 402, 21
  * Nishikawa et al. (2020) Nishikawa H., Boddy K. K., Kaplinghat M., 2020, Phys. Rev. D, 101, 063009
  * Oman et al. (2015) Oman K. A., et al., 2015, MNRAS, 452, 3650
  * Ren et al. (2019a) Ren T., Kwa A., Kaplinghat M., Yu H.-B., 2019a, Phys. Rev. X, 9, 031020
  * Ren et al. (2019b) Ren T., Kwa A., Kaplinghat M., Yu H.-B., 2019b, Physical Review X, 9, 031020
  * Robertson et al. (2021) Robertson A., Massey R., Eke V., Schaye J., Theuns T., 2021, MNRAS, 501, 4610
  * Robles et al. (2017) Robles V. H., et al., 2017, MNRAS, 472, 2945
  * Rocha et al. (2013) Rocha M., Peter A. H. G., Bullock J. S., Kaplinghat M., Garrison-Kimmel S., Oñorbe J., Moustakas L. A., 2013, MNRAS, 430, 81
  * Sameie et al. (2021) Sameie O., et al., 2021, MNRAS, 507, 720
  * Spergel & Steinhardt (2000) Spergel D. N., Steinhardt P. J., 2000, Phys. Rev. Lett., 84, 3760
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Tulin & Yu (2018) Tulin S., Yu H.-B., 2018, Phys. Rep., 730, 1
  * Tulin et al. (2013) Tulin S., Yu H.-B., Zurek K. M., 2013, Phys. Rev. D, 87, 115007
  * Turner et al. (2021) Turner H. C., Lovell M. R., Zavala J., Vogelsberger M., 2021, MNRAS, 505, 5327
  * Vogelsberger et al. (2012) Vogelsberger M., Zavala J., Loeb A., 2012, MNRAS, 423, 3740
  * Weinberger et al. (2020) Weinberger R., Springel V., Pakmor R., 2020, ApJS, 248, 32
  * Zemp et al. (2008) Zemp M., Moore B., Stadel J., Carollo C. M., Madau P., 2008, MNRAS, 386, 1543