# Shapes of Milky-Way-Mass Galaxies with Self-Interacting Dark Matter

Drona Vargya,1 Robyn Sanderson,1,2 Omid Sameie,3 Michael Boylan-Kolchin,3
Philip F. Hopkins,4 Andrew Wetzel,5 Andrew Graus3
1Department of Physics and Astronomy, University of Pennsylvania, 209 South
33rd Street, Philadelphia, PA 19104, USA
2Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue,
New York, NY 10010, USA
3Department of Astronomy, The University of Texas at Austin, 2515 Speedway,
Stop C1400, Austin, TX 78712, USA
4TAPIR, California Institute of Technology, MC 350-17, Pasadena, CA 91125, USA
5Department of Physics and Astronomy, University of California, Davis, One
Shields Avenue, Davis, CA 95616, USA
E-mail<EMAIL_ADDRESS>(DV)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Self-interacting dark matter (SIDM) models offer one way to reconcile
inconsistencies between observations and predictions from collisionless cold
dark matter (CDM) models on dwarf-galaxy scales. In order to incorporate the
effects of both baryonic and SIDM interactions, we study a suite of
cosmological-baryonic simulations of Milky-Way (MW)-mass galaxies from the
Feedback in Realistic Environments (FIRE-2) project where we vary the SIDM
self-interaction cross-section $\sigma/m$. We compare the shape of the main
dark matter (DM) halo at redshift $z=0$ predicted by SIDM simulations (at
$\sigma/m=0.1$, 1, and 10 cm2 g-1) with CDM simulations using the same initial
conditions. In the presence of baryonic feedback effects, we find that SIDM
models do not produce the large differences in the inner structure of MW-mass
galaxies predicted by SIDM-only models. However, we do find that the radius
where the shape of the total mass distribution begins to differ from that of
the stellar mass distribution _is_ dependent on $\sigma/m$. This transition
could potentially be used to set limits on the SIDM cross-section in the MW.

###### keywords:

dark matter – methods: numerical – galaxies: structure – galaxies: haloes

††pubyear: 2021††pagerange: Shapes of Milky-Way-Mass Galaxies with Self-
Interacting Dark Matter–B

##  1 Introduction

The cold dark matter (CDM) plus dark energy $\left(\Lambda\text{CDM}\right)$
cosmological model has been the most successful model for understanding the
large-scale structure of the Universe. However, on length scales smaller than
$\sim 1$ Mpc and masses smaller than $\sim 10^{11}$ M⊙, there are challenges
to this model from discrepancies between predictions from CDM-only simulations
and observational data. On dwarf galaxy ($M_{\star}\lesssim 10^{9}$ M⊙)
scales, these discrepancies include the core-cusp, diversity, missing
satellites, too-big-to-fail (TBFT), and planes-of-satellites “problems”
(Bullock & Boylan-Kolchin, 2017; Tulin & Yu, 2018, and references therein).
The core-cusp problem (Moore, 1994; Kuzio de Naray et al., 2008; Oh et al.,
2008; Walker & Peñarrubia, 2011; Oh et al., 2015) arises from observational
evidence that the cores of dark matter (DM) dominated galaxies are less dense
and less cuspy (in inner density profile slope) than predicted by CDM-only
simulations. The diversity problem, both in the field and among the Milky
Way’s (MW) satellites (Oman et al., 2015; Kaplinghat et al., 2019), arises
from observational evidence that galaxies’ rotation curves appear to be more
varied than in CDM simulations. The missing satellites problem (Klypin et al.,
1999; Moore et al., 1999) refers to the smaller number of observed satellite
(dwarf) galaxies around the MW and the Local Group than the number of
subhaloes predicted by CDM-only simulations. The too-big-to-fail problem
(Boylan-Kolchin et al., 2011, 2012; Garrison-Kimmel et al., 2014) arises from
a mismatch in the central densities of satellite and field galaxies, which are
smaller than predicted by CDM-only simulations; it is expected that haloes so
massive are resistant to star formation suppression from re-ionization (thus
“too big to fail” in forming stars). Finally, the planes-of-satellites
phenomenon refers to the apparent alignment of the orbital planes of satellite
galaxies in the MW (e.g. Lynden-Bell 1976; Kroupa et al. 2005; Pawlowski et
al. 2012; Fritz et al. 2018; Pawlowski & Kroupa 2020; though see also Riley &
Strigari 2020), M31 (e.g. Ibata et al., 2013; Conn et al., 2013), and
Centaurus A (Müller et al., 2018), which is not commonly seen in CDM-only
simulations (e.g. Libeskind et al., 2009; Ibata et al., 2014); however,
transient coplanar configurations are seen when baryons are included (e.g.
Ahmed et al., 2017; Shao et al., 2019; Samuel et al., 2021).

In order to solve these challenges on small scales without affecting large-
scale structure, solutions within the CDM framework have been proposed to
reconcile observations with predictions through a more complete incorporation
of the baryonic physics (Navarro et al., 1996; Mashchenko et al., 2008;
Peñarrubia et al., 2012; Pontzen & Governato, 2012; Governato et al., 2012;
Brooks & Zolotov, 2014a; Oñorbe et al., 2015). Alternatively, the
aforementioned discrepancies may hint toward a theory of DM beyond the CDM
model.

Self-interacting dark matter (SIDM) models (Spergel & Steinhardt, 2000) assume
DM particles can exchange energy and momentum via a force mediator with a
cross-section close to the regime of the strong interaction (Ahn & Shapiro,
2005; Ackerman et al., 2009; Feng et al., 2009; Arkani-Hamed et al., 2009;
Loeb & Weiner, 2011; Tulin et al., 2013). On galactic scales, the DM
interaction rate becomes comparable to Hubble time-scale only deep inside of
the gravitational potential, leaving the outer radii and extragalactic scales
intact (Rocha et al., 2013; Vogelsberger et al., 2016; Sameie et al., 2019;
Bondarenko et al., 2020). In the inner regions of DM-dominated systems,
energy-exchange through self-interactions will result in an isothermal density
profile if the SIDM local collision rate
$\Gamma_{\mathrm{scatter}}\propto\rho_{\mathrm{DM}}\,\sigma_{\mathrm{x}}/m_{\mathrm{x}}\geq
1$, where $\rho_{\mathrm{DM}}$ is the DM density and
$\sigma_{\mathrm{x}}/m_{\mathrm{x}}$ is the self-interaction cross-section per
unit mass. The value of $\sigma_{\mathrm{x}}/m_{\mathrm{x}}$ is constrained by
observations of galaxy clusters (Yoshida et al., 2000; Randall et al., 2008;
Peter et al., 2013a) and of the internal stellar kinematics of MW satellites
(Koda & Shapiro, 2011; Zavala et al., 2013; Valli & Yu, 2018; Correa, 2020;
Hayashi et al., 2020).

Cosmological DM-only simulations of SIDM models have made baseline predictions
for their velocity profiles, density profiles, halo shapes, and substructures
(Rocha et al., 2013; Peter et al., 2013b; Vogelsberger et al., 2012). These
simulations predict isothermal density profiles and spherical shapes for DM
haloes and their substructures. Introducing baryonic components in SIDM haloes
couples the central DM densities to the baryonic potential (Kaplinghat et al.,
2014; Elbert et al., 2018; Sameie et al., 2018), leading to substantial
differences from DM-only predictions when baryons are dynamically important.
This suggests that a plausible explanation of the observed diversity in the DM
distribution in field galaxies and the MW’s satellites could be that it is a
byproduct of baryonic mass assembly and DM self-interactions (Creasey et al.,
2017; Kamada et al., 2017; Ren et al., 2019; Despali et al., 2019; Sameie et
al., 2020b; Sameie et al., 2020a).

DM self-interactions also lead to more spherical halo shapes in SIDM than CDM
(Peter et al., 2013a). Cosmological mass assembly in CDM creates triaxial DM
haloes: since angular momentum exchange is inefficient, the DM particles
retain substantial memory of their initial in-fall directions, resulting in
haloes with ellipsoidal minor-to-major axis ratios as low as $c/a\sim 0.5$
(Vera-Ciro et al., 2011; Butsky et al., 2016). In pure SIDM, particles can
more efficiently exchange angular momentum through interactions, isotropizing
their orbits until $c/a\sim 1$. However, if baryons dominate the gravitational
potential, DM self-interactions tie the DM halo shapes to the baryonic
distribution. Semi-analytic modeling suggests that the SIDM density profile
should scale with the _total_ gravitational potential (Kaplinghat et al.,
2014). If baryons dominate the central density of galaxies, the shape of the
SIDM distribution should then follow that of the baryons. $N$-body SIDM
simulations of MW-mass systems with embedded baryonic discs support these
predictions (Sameie et al., 2018), as do the SIDM cosmological-baryonic
simulations of slightly more massive disc galaxies (at $z\sim 0.2$) by Despali
et al. (2019).

In this work and a companion paper (Sameie et al., 2021), we examine high-
resolution cosmological-baryonic simulations of MW-mass galaxies from the
“Feedback In Realistic Environments" (FIRE) project. The initial conditions
and CDM simulations are part of the second generation of simulations, the
FIRE-2 suite (Hopkins et al., 2018); we also study the same initial conditions
resimulated with several different SIDM cross-sections. As in the original
FIRE-2 suite, gravity and hydrodynamics are solved with GIZMO and baryonic
feedback is implemented with the FIRE-2 feedback recipes, which are held
constant across all runs (for exact details, see Sameie et al., 2021, and
Section 3 of this paper). Simulations with SIDM use the Monte Carlo approach
to scattering described in Rocha et al. (2013). We also resimulate some haloes
without baryons, in both CDM and SIDM, to isolate feedback effects. Our goal
in this work is to gauge the extent to which halo shapes can serve as a
discriminator between CDM and SIDM and the extent to which this depends on the
self-consistent inclusion of baryonic physics.

The CDM cosmological-baryonic versions of these simulations have previously
been shown to produce a realistic population of satellite galaxies around MW-
mass hosts that does not suffer from the missing satellites or TBFT problems
of small-scale structure formation (Wetzel et al., 2016; Garrison-Kimmel et
al., 2019; Samuel et al., 2020b, 2021). Furthermore, studies across mass
scales have shown that the core-cusp (Chan et al., 2015; Oñorbe et al., 2015;
El-Badry et al., 2016) and diversity (El-Badry et al., 2018) problems are also
mitigated with this feedback implementation. Other groups find similar results
with different physics implementations (e.g. Brooks et al., 2013; Brooks &
Zolotov, 2014b; Tollet et al., 2016; Dutton et al., 2016). Since baryonic
physics can thus at least partially reconcile observations with the standard
$\Lambda$CDM cosmological model, we must also carefully gauge whether SIDM,
combined with baryonic feedback, _over_ -corrects the potential small-scale
problems for $\Lambda$CDM.

This work is organized as follows. In §2, we outline the method used to
determine the shapes of haloes. In §3, we detail the suite of simulations used
in this study. In §4, we compare the results for the densities, velocities,
scattering rates, shapes, and triaxiality between FIRE-2 MW-mass CDM and SIDM
with previous results from simulations and observations. In §5, we discuss the
results. In §6, we give a summary of our results and conclusions.

##  2 Methods of Determining Shapes

To determine halo shapes, we use the iterative algorithm introduced in
Dubinski & Carlberg (1991) (also see Allgood et al., 2006; Vera-Ciro et al.,
2011; Bullock & Boylan-Kolchin, 2017; Sameie et al., 2018; Robertson et al.,
2019). This procedure fits a triaxial ellipsoid to the approximate isodensity
surface of particles starting from a series of spherical radii $\\{r\\}$ from
the galactic-centre by determining the weighted inertia tensor for particles
inside each $r$. We determine the axis ratios of these ellipsoids for each
separate species in the simulations (DM, stars, and gas) and for the total
mass distribution (which includes all particles from each species). The axis
lengths of the ellipsoids are labelled as $a(r)\geq b(r)\geq c(r)$, where
$a(r)$, $b(r)$, and $c(r)$ are lengths of the major, intermediate, and minor
semi-axes, respectively. The axis lengths are defined as functions of $r$ to
allow for changing shapes at different radii. The axis ratios are then defined
as:

$s(r)\equiv\frac{c(r)}{a(r)}\,,\quad p(r)\equiv\frac{c(r)}{b(r)}\,,\quad
q(r)\equiv\frac{b(r)}{a(r)}\,.$ (1)

We begin the iterative algorithm by calculating the weighted (or “reduced”)
inertia tensor, which is a symmetric matrix defined as

$\tilde{I}_{ij}(r)=\sum_{n=1}^{N_{c}}\frac{m_{n}\,x_{i,n}\,x_{j,n}}{d_{n}^{2}(r)}\bigg{/}\sum_{n=1}^{N_{c}}m_{n}\,,\quad
i,j\in\left\\{1,2,3\right\\}\,,$ (2)

where $N_{c}$ is the number of particles within the ellipsoid of each
component (or species), $m_{n}$ is the $n$th particle mass, and $x_{i,n}$ is
the $i$th coordinate of the $n$th particle for each component with respect to
a Cartesian coordinate system. In our final, best-fitting ellipsoid coordinate
system, $x_{1}$ ($x_{2}$, $x_{3}$) corresponds to the distance along the major
(intermediate, minor) axis. The tensor is “reduced” by normalizing the
particle positions $\\{x\\}$ by the ellipsoidal distance $d_{n}(r)$111The
unweighted inertia tensor $I_{ij}(r)$ (without the tilde), does not “reduce”
the matrix with the ellipsoidal normalization distance measure $d_{n}(r)$.
(which is measured in the eigenvector coordinate system from the centre to the
$n$th particle), where

$d_{n}^{2}(r)=x_{1,n}^{2}+\frac{x_{2,n}^{2}}{q^{2}(r)}+\frac{x_{3,n}^{2}}{s^{2}(r)}\,.$
(3)

We then find the three eigenvalues
$\left(\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}\right)$ of the matrix
$\tilde{I}_{ij}$ and set the ellipsoidal orientation to the corresponding
orthogonal eigenvectors
$\left\\{\mathbf{e_{1}},\mathbf{e_{2}},\mathbf{e_{3}}\right\\}$ (i.e. the
principal axes). The square roots of the eigenvalues are used to find the axis
ratios: $s=\left(\lambda_{3}/\lambda_{1}\right)^{1/2}$,
$p=\left(\lambda_{3}/\lambda_{2}\right)^{1/2}$, and
$q=\left(\lambda_{2}/\lambda_{1}\right)^{1/2}$. The axis lengths are then
computed with these axis ratios: $a(r)=r$, $b(r)=r\,q(r)$, and $c(r)=r\,s(r)$.
This ensures that the triaxial ellipsoid is contiguous to the bounding sphere
of radius $r$ at two points.

For every $r$, the ellipsoid is initialized as a sphere, i.e.
$s(r)=p(r)=q(r)=1$, and the inertia matrix, eigenvalues, and eigenvectors are
computed. In the second (and every subsequent) iteration, the inertia matrix
is recomputed using particles that fall inside the reshaped and reoriented
ellipsoid from the previous iteration. This method keeps the largest axis
length $a(r)$ constant, and thus, constrains this semi-major axis of the
ellipsoid to lie on the surface of the bounding sphere. Therefore, particles
are added and removed to the set only within spherical radius $r$. We continue
the iterations until either $\Delta s=\left|s_{k}-s_{k-1}\right|$, $\Delta
p=\left|p_{k}-p_{k-1}\right|$, and $\Delta q=\left|q_{k}-q_{k-1}\right|$ are
all $<10^{-3}$, or until a maximum of $k=1000$ iterations.

##  3 Simulations of Milky-Way-Mass Galaxies

This work compares different resimulations of three zoomed-in, cosmological-
baryonic simulations of MW-mass haloes from the Latte suite of FIRE-2
simulations (see Hopkins et al., 2018). The initial conditions for the haloes
are all drawn from the same low-resolution cosmological box (AGORA; Kim et
al., 2014) and are labelled m12f, m12i, and m12m. The size of the zoomed-in
region varies between $2$–$5$ Mpc, depending on the formation history of each
halo. The haloes are selected to have present-day virial masses between
$1.2$–$1.6\times 10^{12}$ M⊙, similar to that of the MW, and to have no
massive companions within at least $5R_{200\mathrm{m}}\sim 1.5$
Mpc.222$R_{200\mathrm{m}}$ is the radius within which the total mass density,
$\bar{\rho}\equiv 3M(<R_{200\mathrm{m}})/4\pi R_{200\mathrm{m}}^{3}$, is 200
times the average matter density.

The full FIRE-2 suite of 18 cosmological-baryonic zooms that have been run at
this mass scale includes simulations that form thin discs similar to the MW’s,
as well as some that form spheroids (El-Badry et al., 2018; Garrison-Kimmel et
al., 2018). We select these particular three systems for resimulation because
the properties of their CDM versions have been extensively compared to the MW.
Their disc structure (Sanderson et al., 2020), their satellite galaxy systems
(Samuel et al., 2020a), and their stellar haloes (Sanderson et al., 2018) all
have reasonable similarity to the MW. These are not the only three systems
within the suite for which this is true, but they are the ones for which the
similarities and differences have been most thoroughly quantified. As in the
CDM versions, we use an initial mass resolution of 7100 M⊙ for the star and
gas particles, and 35000 M⊙ for DM particles, for all resimulations.

The primary difference between the three haloes chosen for resimulation is
their formation histories. m12m forms earliest and has the largest disc scale
radius of the three (Debattista et al., 2019). m12f forms latest and includes
a late interaction with a Small Magellanic Cloud (SMC)-mass galaxy that
disrupts the disc (Sanderson et al., 2018). m12i has an intermediate formation
epoch, the largest proportion of accretion from low-mass satellites, and a
thicker, younger disc with a significant outer warp (Wetzel et al., 2016;
Sanderson et al., 2020).

We compare resimulations of these three haloes with the following set of
variations at redshift $z=0$, summarized in Table 1:

Table 1: Summary of simulated MW-mass galaxy properties. All simulations have baryonic particle mass $\sim 7100$ M⊙ and DM particle mass 35000 M⊙. An interaction cross-section of $\sigma/m=0$ cm2 g-1 indicates standard collisionless CDM. $M_{\mathrm{vir}},\ r_{\mathrm{vir}}$: Bryan & Norman (1998) spherical virial quantities. $r_{-2}$: spherical radius where log-slope of DM density profile is $-2$. $d_{1}$: scattering radius, determined as shown in Figure 7 (§4.1); DNI indicates that $\Gamma_{\mathrm{scatter}}<t_{z=0}^{-1}$ for all $d$ with sufficient particles to determine $\rho_{\mathrm{DM}}(d)$ and $v_{\mathrm{rel}}(d)\approx 1.3\,v_{\mathrm{rms}}(d)$ (see Figure 7); that is, the scattering rate profile “does not intersect” this characteristic value, so $d_{1}$ is undefined. $M_{\star,90},\ r_{\star,90}$: Mass and spherical radius of 90% of stellar mass within 30 kpc of the central galaxy. ${\dagger}$: $d_{1}$ is determined using ellipsoidally-averaged DM local collision rate profiles from ellipsoidally-averaged density and velocity profiles; all other quantities in this table are determined with spherical averaging. ${\ddagger}$: Values for all simulations are taken at redshift $z=0$ except for m12i SIDM+Baryon $\sigma/m=10$ cm2 g-1, which is evaluated at $z=0.1$. Initial Conditions | $\sigma/m$ | $M_{\mathrm{vir}}$ | $r_{\mathrm{vir}}$ | $r_{-2}$ | $d_{1}\,^{\dagger}$ | $M_{\star,90}$ | $r_{\star,90}$ | Reference
---|---|---|---|---|---|---|---|---
| $\left[\text{cm}^{2}\text{ g}^{-1}\right]$ | $\left[10^{12}\ \text{M}_{\odot}\right]$ | $\left[\text{kpc}\right]$ | $\left[\text{kpc}\right]$ | $\left[\text{kpc}\right]$ | $\left[10^{10}\ \text{M}_{\odot}\right]$ | $\left[\text{kpc}\right]$ |
CDM-only |  |  |  |  |  |  |  |
m12f | 0 | 1.28 | 284.2 | 64.6 | – | – | – | Garrison-Kimmel et al. (2017)
m12i | 0 | 0.90 | 252.8 | 24.5 | – | – | – | Wetzel et al. (2016)
m12m | 0 | 1.14 | 273.9 | 35.5 | – | – | – | Garrison-Kimmel et al. (2019)
SIDM-only |  |  |  |  |  |  |  |
m12f | 1 | 1.28 | 284.0 | 38.9 | DNI | – | – | Sameie et al. (2021)
m12f | 10 | 1.25 | 282.0 | 28.2 | 24. | – | – | Sameie et al. (2021)
CDM+Baryon |  |  |  |  |  |  |  |
m12f | 0 | 1.33 | 287.9 | 19.5 | – | 5.3 | 12.4 | Garrison-Kimmel et al. (2017)
m12i | 0 | 0.96 | 258.4 | 17.0 | – | 3.2 | 16.7 | Wetzel et al. (2016)
m12m | 0 | 1.23 | 280.8 | 18.6 | – | 4.9 | 21.2 | Hopkins et al. (2018)
SIDM+Baryon |  |  |  |  |  |  |  |
m12m | 0.1 | 1.22 | 279.9 | 17.0 | DNI | 6.1 | 23.0 | Sameie et al. (2021)
m12f | 1 | 1.36 | 289.8 | 13.5 |  8.8 | 6.2 | 15.7 | Sameie et al. (2021)
m12i | 1 | 0.98 | 260.0 | 10.7 |  7.4 | 5.0 | 13.9 | Sameie et al. (2021)
m12m | 1 | 1.24 | 281.5 |  9.8 |  9.8 | 6.6 | 20.2 | Sameie et al. (2021)
m12f | 10 | 1.27 | 283.2 | 28.2 | 23. | 5.1 | 15.9 | Sameie et al. (2021)
m12i$\,{}^{\ddagger}$ | 10 | 0.92 | 237.8 | 24.5 | 19. | 4.5 | 12.3 | Sameie et al. (2021)
m12m | 10 | 1.22 | 279.6 |  5.4 | 22. | 8.0 | 20.4 | Sameie et al. (2021)

  1. 1.

CDM-only – Collisionless CDM without baryons, for all three haloes;

  2. 2.

SIDM-only – Collisional self-interacting DM without baryons, at $\sigma/m=$ 1
and 10 cm2 g-1 for m12f;

  3. 3.

CDM+Baryon – Collisionless CDM with baryons and full hydrodynamics, using
FIRE-2 feedback recipes, for all three haloes;

  4. 4.

SIDM+Baryon – Collisional self-interacting DM with baryons and full
hydrodynamics, with identical baryonic physics to the fiducial suite, at
$\sigma/m=0.1$ cm2 g-1 for m12m, and 1 and 10 cm2 g-1 for all three haloes,
but at redshift $z=0.1\ \left(t=12.5\ \textrm{Gyr}\right)$ for m12i at the
latter cross-section.

All DM self-interactions are realized using Monte Carlo elastic (non-
dissipative) scattering, as described in Rocha et al. (2013). We evaluate the
m12i SIDM+Baryon $\sigma/m=$ 10 cm2 g-1 at $z=0.1$, the latest epoch currently
available for this resimulation. Based on the behavior observed in the other
runs, we expect the radial density and velocity profile of this simulated
galaxy to be relatively stable between $z=0.1$ and $z=0$. The shape profiles
of the various species continue to evolve to $z=0$ in the central parts of the
galaxies, but this effect is least pronounced for the DM component.

The baryonic runs listed above and in Table 1 use the standard set of FIRE-2
feedback recipes (Hopkins et al., 2018) with one exception, which is to ignore
the conversion of thermal to kinetic energy in the unresolved Sedov-Taylor
phase of the expansion of shocks generated by mass loss from massive stars. As
discussed in (Sameie et al., 2021), the default prescription in FIRE-2 had the
effect of converting nearly all the thermal energy into momentum, giving the
stellar winds a similar effect to a miniature supernova and resulting (perhaps
counter-intuitively) in higher star formation rate (SFR) and stellar mass in
the simulated galaxies, and subsequently less diversity among density
profiles. However, for this study we use the resimulations of the CDM haloes
that ignore this “$PdV$” work for the sub-res regions, for consistency with
the existing SIDM runs.

Figure 1 shows Spectral Energy Distribution (SED)-modeled and ray-traced
images (integrated over filter band-passes similar to Hubble Space Telescope
visible light) of the simulated galaxies in a wide-angle view focused on the
stellar halo. It is apparent that although there are some small differences,
every simulation contains a reasonable-looking spiral galaxy, system of dwarf
satellite galaxies, and stellar halo. One obvious difference between
resimulations of the same initial conditions is in the particular dwarf
satellites that appear in each simulation, which vary between even identical
runs due to the stochastic impact of supernovae. The most massive satellites
are stable to this effect but can have small phase differences in their orbits
from run to run, as is evident in m12f. The other obvious difference is in the
star formation rate (as evident in m12i). Star formation is a highly non-
linear process, and also varies stochastically from run to run even for
identical initial conditions, again mainly because of the random occurrence
and clustering of supernovae and also (as clearly seen in m12f) from
differences in the orbital phase of mergers. However, these differences do not
significantly change the global properties of the haloes and their central
galaxies, as we next demonstrate.

Figure 1: Face-on halo view of cosmological-baryonic simulated MW-mass
galaxies. SED-modeled, ray-traced images of starlight in three sets of MW-mass
galaxy simulations: m12f (left), m12i (centre), and m12m (right); all
DM+Baryon simulations are shown except m12m SIDM+Baryon $\sigma/m=0.1$ cm2 g-1
(which is quite similar to the m12m CDM+Baryon case). Each panel is 300 kpc
across and the galaxy has been rotated to show the disc face-on. All
simulations except m12i SIDM+Baryon $\sigma/m=10$ cm2 g-1 (discussed in §3)
are depicted at $z=0$.

Figure 2 also shows SED-modeled and ray-traced images, but with close-up views
of the stellar discs in the simulated systems. There is remarkable uniformity
in the structure and size of the discs across all DM simulations. The m12f
SIDM simulation has a slightly higher SFR in its outskirts, probably due to
the timing of a merger with a roughly $\sim$SMC-mass object visible in Figure
1, while m12i and m12m have no noticeable increase in star formation with
$\sigma/m$. Generally, the discs of the CDM galaxies tend to be slightly more
massive and compact than in SIDM.

Figure 2: Face-on disc view of cosmological-baryonic simulated MW-mass
galaxies. SED-modeled, ray-traced images of starlight in three sets of MW-mass
galaxy simulations: m12f (left), m12i (centre), and m12m (right); all
DM+Baryon simulations are shown except m12m SIDM+Baryon $\sigma/m=0.1$ cm2 g-1
(which is quite similar to the m12m CDM+Baryon case). Each panel is 30 kpc
across and the galaxy has been rotated to show the disc face-on. All
simulations except m12i SIDM+Baryon $\sigma/m=10$ cm2 g-1 (discussed in §3)
are depicted at $z=0$.

Figure 3 summarizes the characteristic masses and radii of the resimulations,
quantifying the impressions given by examining Figures 1 and 2. The left-hand
panel shows the DM virial mass $M_{\mathrm{vir}}$ versus DM scale radius
$r_{-2}$; the right-hand panel shows the stellar-to-halo-mass ratio
$M_{\star,90}/M_{\mathrm{vir}}$ versus the radius enclosing 90% of the stellar
mass $r_{\star,90}$. These values are all computed using spherical volumes.
$M_{\mathrm{vir}}$ values are roughly the same over scattered domains of
$r_{-2}$ for each set of simulations, while the
$M_{\star,90}/M_{\mathrm{vir}}$ ratios are generally more scattered over
$r_{\star,90}$. The CDM-only and SIDM-only simulations have larger $r_{-2}$
than their CDM+Baryon and SIDM+Baryon counterparts due to the increased
concentration produced by the central baryonic component. For CDM+Baryon and
SIDM+Baryon, the m12f and m12i simulations have the smallest and largest
$M_{\mathrm{vir}}$, respectively, while the m12m simulations fall in-between.
The interaction cross-section does not otherwise seem to produce any clear
trends in the global DM distribution; m12m’s scale radius decreases as
$\sigma/m$ increases, while m12f and m12i show no clear trend. In all cases
the virial mass varies by less than 10% across all cross-sections. Finally,
steady growth in $M_{\rm vir}$ and $M_{\star,90}$ for m12i at 10 cm2 g-1 from
$0.1$ (plotted) to $z=0$ (all other runs) would bring this halo into
consistency with the other m12i simulations, as would further contraction of
the DM distribution due to the baryonic component (leading to a decrease in
$r_{-2}$).

In terms of the stellar distributions, there is significant variation in the
stellar-to-halo mass ratio across the different resimulations, while
$r_{\star,90}$ appears roughly independent of $\sigma/m$ for most cases
(although $r_{\star,50}$ does have a trend with cross section; see Sameie et
al. 2021). Again m12m shows the opposite trend from m12i and m12f, as well as
the largest variation in $M_{\star,90}/M_{\mathrm{vir}}$. While in m12f and
m12i $M_{\star,90}/M_{\mathrm{vir}}$ shows no trend for larger $\sigma/m$, for
m12m a larger $\sigma/m$ gives rise to a relatively more massive central
galaxy (recall that the DM halo does not change appreciably in mass between
runs).

Figure 3: Global DM and stellar properties of simulated MW-mass galaxies. DM
virial mass $M_{\mathrm{vir}}$ as a function of DM scale radius $r_{-2}$
(left) and stellar-to-halo-mass ratio $M_{\star,90}/M_{\mathrm{vir}}$ as a
function of radius enclosing 90% of stellar mass (right) for m12f (circles),
m12i (pentagons), and m12m (squares). In the left-hand panel (DM properties),
the CDM-only, SIDM-only, CDM+Baryon, and SIDM+Baryon simulations are
represented with hollow cyan, solid cyan, hollow blue, and solid blue markers,
respectively; in the right-hand panel (stellar properties), the CDM+Baryon and
SIDM+Baryon simulations are represented with hollow green and solid green
markers, respectively. In both panels, the SIDM+Baryon simulations with
$\sigma/m=0.1$, 1, and 10 cm2 g-1 are shown with these respective numbers
inside the markers. The simulation of m12i SIDM+Baryon at $\sigma/m=10$ cm2
g-1 is evaluated at $z=0.1$ (instead of $z=0$ like all other simulations),
indicated with a decrease in alpha (lighter shade).

##  4 Results

In this §, we compare the density, velocity, and shape profiles of the
different simulations described in §3. Throughout the rest of this work, we
use the same series of line-styles to denote different DM cross-sections,
different colours to show the different species, and gradient shaded-areas to
indicate the SIDM local collision regions (LCR) for different $\sigma/m$;
these are given in Figure 4 and apply to Figures 5 through 11, as well as
Figures 12, 13, and 14 in the Appendix.

Figure 4: Legend for comparisons of simulated MW-mass galaxies in CDM and
SIDM. Line-styles denote different DM cross-sections listed above in units of
cm2 g-1, colours show the different species, and shaded-areas indicate the
SIDM local collision regions (LCR). These line-styles, colours, and shaded-
areas are used throughout the rest of this work and the Appendix.

###  4.1 Densities, Velocities, and Scattering Rates

Predictions for density and velocity distributions in SIDM are closely related
to those for the halo’s shape, since the same interactions that heat the inner
regions also make the outer halo more spherical, by preferentially scattering
DM particles on plunging, radial orbits. Before discussing the shapes, we will
review the density and velocity profiles for DM (discussed in full in Sameie
et al., 2021) and discuss the profiles for the stellar and gas components as
well. While Sameie et al. (2021) presents spherically-averaged profiles, here
we show profiles computed using ellipsoids fit to isodensity contours for each
species (§2). The difference in density, and therefore DM scattering rate, can
vary by up to 60% from the spherically-averaged value, depending on the
flattening (Appendix A).

To compute the density $\rho$, we use the mass enclosed in shells that follow
the triaxial ellipsoidal surfaces calculated using the method described in §2,
which approximately follow isodensity surfaces. However, we use much broader
spacing in distance for the density calculation than for the ellipsoid fits,
spacing shells by roughly every 10th point in $r$ for which a fit is carried
out, to allow enough space between shells to get sufficient numbers of
particles and to mitigate problems caused by the twisting of the ellipsoid
axes between shells. To estimate the density at ellipsoidal distance
$d_{\kappa}$, we select all particles $N_{\kappa}$ within a shell $\kappa$
centered on the ellipsoid with semi-major axis $a_{\kappa}$ (discarding shells
with $N_{\kappa}<100$). The shell half-thickness $\Delta d_{\kappa}$ is set by
the difference in semi-major axis between the isodensity surface at
$d_{\kappa}$ and its inner neighbor, such that $\Delta d_{\kappa}\equiv
a_{\kappa}-a_{\kappa-1}$. Then the density $\rho(d_{\kappa})$ is computed as

$\rho(d_{\kappa})\equiv\frac{\sum_{n=1}^{N_{\kappa}}m_{n}}{\frac{4}{3}\pi
q_{\kappa}s_{\kappa}\left[\left(d_{\kappa}+\Delta
d_{\kappa}\right)^{3}-\left(d_{\kappa}-\Delta
d_{\kappa}\right)^{3}\right]}\,,$ (4)

where $q_{\kappa}$ and $s_{\kappa}$ are the axis ratios of the ellipsoid used
to calculate $d_{\kappa}$. Likewise the root-mean-square (RMS) velocity
$v_{\mathrm{rms}}(d_{\kappa})$ is

$v_{\mathrm{rms}}(d_{\kappa})\equiv\sqrt{\frac{1}{N_{\kappa}}\sum_{n=1}^{N_{\kappa}}\mathbf{v}_{n}\cdot\mathbf{v}_{n}}\,,$
(5)

where $\mathbf{v}_{n}$ is the velocity vector of particle $n$ inside shell
$\kappa$.

Figure 5: Ellipsoidally-averaged density profiles of the simulated MW-mass
galaxies. Density $\rho$ (top) and $\rho d^{2}$ (bottom; an ideal isothermal
profile is flat in this view) as a function of ellipsoidal distance $d$ for
m12f (left), m12i (centre), and m12m (right). Line-styles, colours, and
shaded-areas follow the legend in Figure 4. Figure 6: Ellipsoidally-averaged
velocity profiles of the simulated MW-mass galaxies. Root-mean-square velocity
$v_{\mathrm{rms}}$ as a function of ellipsoidal distance $d$ for m12f (left),
m12i (centre), and m12m (right). Line-styles, colours, and shaded-areas follow
the legend in Figure 4.

We also calculate the DM local collision rate:

$\displaystyle\Gamma_{\mathrm{scatter}}(d)$
$\displaystyle=\rho_{\mathrm{DM}}(d)\,v_{\mathrm{rel}}(d)\,\sigma/m$
$\displaystyle\approx\frac{4}{\sqrt{3\pi}}\,\rho_{\mathrm{DM}}(d)\,v_{\mathrm{rms}}(d)\,\sigma/m\,,$
(6)

where $\rho_{\mathrm{DM}}(d)$ is the DM density and $v_{\mathrm{rel}}(d)$ is
the relative DM velocity. We approximate $v_{\mathrm{rel}}$ by
$4/\sqrt{3\pi}\,v_{\mathrm{rms}}\approx 1.3\,v_{\mathrm{rms}}$, which holds
exactly for a Maxwellian velocity distribution. We define the local collision
region (LCR) scattering radius $d_{1}$ by requiring
$\Gamma_{\mathrm{scatter}}(d_{1})=t_{z=0}^{-1}$, where $t_{z}$ is the time at
redshift $z$, thus $t_{z=0}=13.8$ Gyr is age of the Universe (Planck
Collaboration et al., 2014, Table A.1: Planck+WP+BAO).333m12i SIDM+Baryon
$\sigma/m=10$ cm2 g-1 is evaluated at $z=0.1$, which gives
$\Gamma_{\mathrm{scatter}}(d_{1})=t_{z=0.1}^{-1}$ where $t_{z=0.1}=12.5$ Gyr.
Thus, inside $d_{1}$, DM particles have experienced at least one self-scatter
within the age of the Universe (based on the density distribution at $z=0$),
giving an approximate volume inside which we expect the SIDM differences to be
the largest. The $d_{1}$ for different simulations are given in Table 1, and
the LCR for different SIDM+Baryon $\sigma/m$ are represented as gray shaded-
areas (as shown by the legend in Figure 4) and used in Figures 5, 6, 8, 9, 11,
13, and 14.

As expected, the SIDM-only density profiles for m12f (Figure 5, left) have
much lower densities in the central region compared to all the other DM
profiles, while the DM profiles in all the CDM-only, CDM+Baryon, and
SIDM+Baryon simulations are remarkably similar. The SIDM-only densities are
too low to produce smooth curves for $r\lesssim 2$ kpc.

In the simulations with baryons, the density of stars dominates the central
region, while the DM and gas approximately follow each other about a magnitude
below the stellar component. The stellar mass density of m12m is lower in the
central region and higher at larger radii compared to that of m12f and m12i.
Notably, when examining the bottom row showing $\rho d^{2}$, there does not
appear to be a transition at $d_{1}$ from an isothermal (flat in this view)
profile at $r<d_{1}$ to a Navarro-Frenk-White (NFW)-like profile at $r>d_{1}$,
as posited in Spergel & Steinhardt (2000). This is another illustration of the
effect of the growing galaxy in the centre in altering the density profile
well beyond the baryon-dominated region, thanks to the significant radial
anisotropy in the DM velocity distribution.

Root-mean-square (RMS) velocity profiles are shown in Figure 6. The CDM-only
and SIDM-only simulations have much lower central RMS velocities than the
corresponding DM curves in the CDM+Baryon and SIDM+Baryon simulations, except
in the extreme outskirts ($r\gtrsim 100$ kpc). The additional component added
by the baryons deepens the central potential and subsequently increases the DM
density, leading to higher RMS velocities. The biggest differences are between
CDM-only and CDM+Baryon in the central region for m12f and m12i. The velocity
dispersion in the CDM+Baryon case is much higher than for CDM-only in both
these simulations due to baryonic deepening of the potential. In the SIDM
case, the effects of baryonic contraction are offset by the ability of the DM
to thermalize; that is, to transfer some of the energy outward that would
otherwise go into raising the central velocity dispersion. Surprisingly, for
m12m, the SIDM+Baryon runs have significantly larger central velocity
dispersion than the CDM+Baryon run does. This is probably related to the fact
that the central galaxy in m12m is substantially more massive as $\sigma/m$
increases without becoming significantly larger in extent (Figure 3), implying
that the average stellar density is larger in the SIDM+Baryon runs than the
CDM+Baryon one. Indeed, the central galaxy in m12i shows the same tendency for
the SIDM systems to be slightly more massive in stellar mass, but has slightly
less variation between DM models than m12m. This variety illustrates how the
central density and velocity dispersion of the DM are shaped by interplay with
the growing galaxy in the center of the halo.

Figure 7 shows the local collision rate $\Gamma_{\mathrm{scatter}}$ as a
function of ellipsoidal distance $d$. The scattering radius $d_{1}$ is smaller
for $\sigma/m=1$ cm2 g-1 than 10 cm2 g-1, while the rate for 0.1 cm2 g-1 is
less than $t_{z=0}^{-1}$ at all radii for the $z=0$ density distribution
(though this is not necessarily true at all $z$). The profile for m12f SIDM-
only $\sigma/m=1$ cm2 g-1 also does not reach $t_{z=0}^{-1}$, while 10 cm2 g-1
does. This indicates that while self-interactions may have occurred earlier,
the resulting heating reduces the scattering rate by $z=0$ to less than
$t_{z=0}^{-1}$ everywhere. On the other hand, all SIDM haloes with baryons
have substantially higher interaction rates at the present day, likely as a
result of the additional depth in the gravitational potential created by the
central galaxy. The growing galaxy can thus amplify the effect of a nonzero
SIDM cross-section in the central portion of the halo by keeping the
scattering rate higher over time.

Figure 7: Ellipsoidally-averaged DM local collision rate profiles of the
simulated MW-mass SIDM galaxies. The local DM self-scattering rate
$\Gamma_{\mathrm{scatter}}$ (Equation 6) is shown as a function of ellipsoidal
distance $d$ for SIDM simulations. Line-styles and colours follow the legend
in Figure 4. The point where $\Gamma_{\mathrm{scatter}}$ intersects the gray
horizontal lines at $t_{z=0}^{-1}$ and $t_{z=0.1}^{-1}$ marks the scattering
radius $d_{1}$, shown as the shaded-areas of Figures 5, 6, 8, 9, 11, 13, and
14.

###  4.2 Shape Profiles

The profiles for the minor-to-major axis ratio, $s=c/a$, for the different
haloes in the simulations are shown in Figure 8. As summarized in Figure 4,
different colours distinguish between simulations with and without baryons and
among species (DM, stars, and gas), while line-styles show different DM
interaction cross-sections $\sigma/m$. Axis ratios are calculated using the
method described in §2, and we remove all data points where the triaxial
ellipsoids do not enclose at least 5000 particles (see Appendix A of Vera-Ciro
et al., 2011, where at least 3000 particles are used).

Figure 8: Shape profiles of simulated MW-mass galaxies. Minor-to-major axis
ratio $s=c/a$ is shown as a function of semi-major axis $r=a$ for different
species (DM, stars, and gas) in three sets of MW-mass galaxy simulations: m12f
(left), m12i (centre), and m12m (right). Line-styles, colours, and shaded-
areas follow the legend in Figure 4. An alternate version of this plot using
the geometrical mean of the axis lengths, $R=\left(abc\right)^{1/3}$, can be
found in Appendix B.

Axis ratio differences between the CDM+Baryon and SIDM+Baryon simulations are
small compared to the differences between the CDM-only and SIDM-only
simulations over the same range of cross-sections (Figure 8). SIDM-only
creates the most spherical DM haloes, obtained for the highest $\sigma/m$.
Adding baryons to these simulations increases the overall roundness, but
increasing $\sigma/m$ does not produce the large changes seen in the non-
baryon cases. The effect of the growth of the central galaxy is to standardize
the shape in roughly the same range of axis ratios, between approximately 0.6
and 0.8. In m12f (left) and m12i (centre), most of the scaling of the shape
with cross-section is also erased. In m12m (right) there are still substantial
differences between haloes with different cross-sections, but the variation in
shape has been “recentered” around the CDM case, while with DM-only the shape
just gets progressively rounder for larger cross-sections.

In the central region of these MW-mass galaxies ($r\lesssim 10$ kpc),
CDM+Baryon produces a more spherical DM distribution than SIDM+Baryon in
nearly all cases. This is consistent with, though less pronounced than, the
effects discussed in Sameie et al. (2018), but contrary to expectations from
analytic predictions considering only DM (Tulin & Yu, 2018). The _stellar_
distribution is also frequently flatter for SIDM than CDM across all three
galaxies. These differences are greatest in the region where SIDM is
collisional (gray shaded-areas), and the degree of flattening in the DM
distribution parallels the flattening in the stars, indicating ongoing
dynamical coupling between the stellar and DM distributions in the inner
galaxy.

It is also apparent from the shape curves, which stop when the ellipsoid no
longer encloses at least 5000 particles, that the SIDM haloes with baryons
remain much denser in their centers at late times than those without baryons
(as discussed in §4.1). The LCR gray shaded-areas are calculated from the
_present-day_ DM densities (Figure 6), indicating that while the SIDM-only
haloes have reached an equilibrium where even their innermost regions have a
very low collision rate relative to the age of the Universe, the gradually
increasing potential depth due to the central galaxy counteracts collisional
heating and maintains a much higher central collision rate at late times. This
also supports the idea that there is ongoing information exchange between the
DM and stars in the inner galaxy, and suggests that this region may not be in
equilibrium between the two species.

The m12m series of simulations contains most of the exceptions to these
generalizations and is thus worth discussing in more detail. This galaxy is
the earliest of the three to form and has the largest scale radius (Figure 3),
meaning that its baryonic component has had the longest time to shape the DM
distribution (and vice-versa) over the largest range of radii. As pointed out
in §3, its global properties have the opposite trend with $\sigma/m$ from the
other two haloes. It is also the only halo whose outer shape ($\gtrsim 20$
kpc) is consistent between the CDM-only and CDM+Baryon cases, and the only one
where there is significant variation of the shape with cross-section across
all radii. Its particular assembly history (early accretion of many small
galaxies) thus appears particularly sensitive to SIDM effects. This could be
because it simply has more time to establish equilibrium between the SIDM and
stellar components, and is driven there more rapidly by a central relatively
large galaxy that forms early.

A full comparison of the central DM shapes in these simulations is
challenging, since many of the profiles are noisy due to low particle number.
This is especially true for the m12f SIDM-only simulations, whose shape
profiles inside $5$–$10$ kpc are ambiguous thanks to their low central
densities. In all the DM+Baryon simulations, the shape of the stellar
component is well-resolved (at all radii) and closely follows the shape of the
_total_ mass distribution (at lower radii), which is accessible through
dynamical modeling. We will examine in more detail whether the differences in
the three-dimensional (3D) shape of the _stellar_ distribution are observable
in future work.

###  4.3 Triaxiality

To better understand the shapes, we also calculate the triaxiality parameter
$T$ (Franx et al., 1991):

$T\equiv\frac{a^{2}-b^{2}}{a^{2}-c^{2}}=\frac{1-q^{2}}{1-s^{2}}\,,$ (7)

where an ellipsoid is oblate if $0<T<1/3$ ($T=0$ is a perfect oblate
distribution, $c\ll a=b$), triaxial if $1/3<T<2/3$ ($T=0.5$ is a maximally
triaxial distribution), and prolate if $2/3<T<1$ ($T=1$ is a perfect prolate
distribution, $b=c\ll a$). Figure 9 shows the triaxiality for the stellar and
DM components of all simulations, as well as the total mass distribution.
There are a wide variety of behaviors on display. Consistent with our other
results, m12m has significantly different behavior than m12f and m12i.

Figure 9: Triaxiality profiles of simulated MW-mass galaxies. Triaxiality $T$
is shown as a function of semi-major axis $r=a$ for different species (DM and
stars) in three sets of MW-mass galaxy simulations: m12f (left), m12i
(centre), and m12m (right). Line-styles, colours, and shaded-areas follow the
legend in Figure 4. An alternate version of this plot using the geometrical
mean of the axis lengths, $R=\left(abc\right)^{1/3}$, can be found in Appendix
B.

The DM-only haloes for all three set of simulations are highly prolate or
triaxial in the centre, especially m12f. For CDM, this triaxiality/prolateness
is well understood. Interestingly, though, at large radii m12i and m12f stay
relatively triaxial while the m12m CDM-only halo transitions to an oblate
shape, which is unusual for a typical CDM-only simulated halo at this mass
scale. The DM+Baryon haloes m12f and m12i tend to be quite oblate at
intermediate radii ($2\lesssim r\lesssim 50$ kpc), with very little difference
between DM models. At large radii ($r\gtrsim 50$ kpc) most haloes (even m12m)
transition to a somewhat more triaxial shape. This is also expected since the
SIDM interaction rates (Figure 7) are quite low at these distances, so the DM
behavior should not differ appreciably from CDM. The degree of triaxiality at
large radii varies substantially, however, with no clear trend with
$\sigma/m$.

In the inner regions of the haloes, m12f shows significant triaxiality in the
CDM+Baryon case and more oblate structure in the SIDM+Baryon cases, while in
m12m there is a huge variation in the degree of triaxiality in the central
part of the halo. Interestingly, the triaxiality of the _total_ mass
distribution does not universally follow the stellar distribution in the inner
galaxy the way that the $s=c/a$ axis ratio does.

###  4.4 Comparison to Previous Work

We first compare the DM shape profiles to previous results from DM-only
simulations, those with an analytic disc model, and CDM+Baryon simulations
(see Sameie et al., 2018, Figure 7 for a summary). In Figure 10, we plot the
shape of the DM component for all the FIRE-2 MW-mass galaxy simulations of
Table 1. This comparison uses log-scale for the $y$-axis axis ratio $s=c/a$
and geometric mean of the axis lengths $R=\left(abc\right)^{1/3}$ for the
$x$-axis, for direct comparison with the summary in Figure 7 of Sameie et al.
(2018). Plots of the axis ratio $s=c/a$ and triaxiality $T$ versus geometric
mean radius $R$ for all individual species (as in Figures 8 and 9) are given
in Appendix B.

Figure 10: Shape profiles of DM in simulated MW-mass galaxies compared to
estimates for the MW. Axis ratio $s=c/a$ is shown as a function of geometrical
mean radius $R=\left(abc\right)^{1/3}$ for the DM component in all
simulations. Line-styles and colours follow the legend in Figure 4 (but with
thicker lines) for the simulations. The dark blue dashed-line follows a
previous SIDM+Baryon semi-analytic model for $\sigma/m=1$ cm2 g-1 (see Figure
7 of Sameie et al., 2018). The SIDM+Baryon simulations are generally more
spherical than CDM-only, but not notably different than CDM+Baryon, and show
less variation with $\sigma/m$ than in SIDM-only. All simulations except CDM-
only are loosely consistent with the various estimates for the MW. For
simulations with baryons, halo-to-halo variation is comparable to the
difference between DM models. The estimated values of the halo density
flattening parameter $q_{\rho}^{\mathrm{DM}}=\left(c/a\right)_{\rho}$ are
summarized in Table 2, and are plotted over their approximate regions of
validity (i.e. the range spanned by the data used for the constraints). The
coloured boxes and lines indicate the various estimates for
$q_{\rho}^{\mathrm{DM}}=\left(c/a\right)_{\rho}$ summarized in Table 2.
Constraints derived from models of individual tidal streams are shown in
orange (Sagittarius stream), red (Palomar 5 stream), green (Grillmair-
Dionatos-1 stream), and yellow (statistically detected streams around NGC3201
and M68). Constraints based on Jeans modeling are shown in purple, and
measurements combining equilibrium assumptions with the disc rotation curve
and/or other disc data are shown in brown.

We see a clear trend towards more spherical haloes at larger $\sigma/m$ for
the DM-only simulations, but find that the SIDM+Baryon simulations are not as
spherical at larger radii as expected from the SIDM+Baryon semi-analytic model
of Sameie et al. (2018). Also, instead of the concave-up shape predicted by
this model for SIDM+Baryon (with the innermost and outermost regions the most
spherical), we see a concave-down trend for all the curves (where the
intermediate radii are most spherical). This appears consistent with the idea
that SIDM can respond more quickly to the influence of the central galaxy than
CDM (which generally has a more spherical inner halo when compared across
resimulations of the same initial conditions). It is also consistent with the
picture that CDM and SIDM should behave similarly in the halo outskirts, where
the shape is driven mainly by the connection with the local filaments (see
e.g. Vera-Ciro et al., 2011) and therefore tends to be less spherical than at
intermediate radii. However, the variation in formation histories across the
different haloes dominates over the variation with $\sigma/m$.

We also compare our results to estimates of the halo density flattening
parameter (minor-to-major axis ratio) $q_{\rho}^{\mathrm{DM}}=(c/a)_{\rho}$
from the literature, summarized in Table 2. This quantity is sometimes
referred to in the literature (see e.g. Hattori et al., 2020) as simply
parameter $q$, but is changed here to distinguish from our intermediate-to-
major axis ratio $q(r)=b(r)/a(r)$. In modeling the kinematics of various MW
tracers of the potential such as tidal streams, globular clusters, or “field”
halo stars unassigned to a given stream, the parameter
$q_{\rho}^{\mathrm{DM}}$ usually represents the flattening of the best-fitting
axisymmetric NFW model for the DM density, and is comparable to our minor-to-
major axis ratio $s(r)=c(r)/a(r)$. In the case of Law & Majewski (2010), where
a triaxial halo is used, we cite the value of $(c/a)_{\rho}$ quoted by the
authors within 20 kpc, which was determined by fitting ellipsoids to the
contours of the Laplacian of the potential.

Table 2: Measurements of the shape of the MW halo from previous work. Uncertainties on $q_{\rho}^{\mathrm{DM}}=\left(c/a\right)_{\rho}$ parameter values are quoted as given in the various works and plotted in Figure 10, which usually correspond to 90% confidence intervals or equivalent. Values with no uncertainties have no easily interpretable range given in the corresponding paper, or are lower limits. $r_{\mathrm{min}}$ and $r_{\mathrm{max}}$ denote the approximate range of galactocentric radii over which the estimates are made, either as given by the authors or as specified for the data-set used. _Notes_ : (a) This paper includes the rotation curve, velocity dispersion, and vertical force profile of the disc as additional constraints. (b) In these two cases, a flattening of $q_{\rho}^{\mathrm{DM}}=0.8$ in the density was assumed to detect the stream statistically before using it to fit a parametrized mass model. (c) As pointed out by Hattori et al. (2020), this paper finds a prolate aspect ratio using an action finder that has known difficulties for orbits in prolate mass distributions. (d) As Law & Majewski (2010) point out and Debattista et al. (2013) confirm, this model does not admit a stable galactic disc; Vera-Ciro & Helmi (2013) shows that the discrepancy can be explained by the influence of the LMC. Reference | $q_{\rho}^{\mathrm{DM}}=\left(c/a\right)_{\rho}$ | $r_{\mathrm{min}}$ | $r_{\mathrm{max}}$ | Data used | Colour
---|---|---|---|---|---
|  | $\left[\text{kpc}\right]$ | $\left[\text{kpc}\right]$ |  |
Hattori et al. (2020) | $0.993^{+0.01}_{-0.005}$ | 1 | 30 | RR Lyrae + constraintsa | brown
Palau & Miralda-Escudé (2020) | $0.87^{+0.02}_{-0.02}$ | 8 | 38 | NGC3201b | yellow
Malhan & Ibata (2019) | $0.82^{+0.25}_{-0.13}$ | 13.5 | 15 | Grillmair-Dionatos-1 (GD-1) | green
Posti & Helmi (2019) | $1.30^{+0.25}_{-0.25}$ | 30 | 150 | Globular clustersc | purple
Wegg et al. (2019) | $1.00^{+0.09}_{-0.09}$ | 1 | 30 | RR Lyrae | purple
Palau & Miralda-Escudé (2019) | $0.87^{+0.06}_{-0.06}$ | 7 | 43 | M68b | yellow
Bovy et al. (2016) | $1.3^{+0.5}_{-0.3}$ | 13.5 | 15 | Grillmair-Dionatos-1 (GD-1) | green
Bovy et al. (2016) | $0.93^{+0.16}_{-0.16}$ | 5 | 19 | Palomar 5 (Pal 5) | red
Küpper et al. (2015) | $0.95^{+0.16}_{-0.12}$ | 5 | 19 | Palomar 5 (Pal 5) | red
Loebman et al. (2014) | $0.4^{+0.1}_{-0.1}$ | 10 | 50 | SDSS halo stars | purple
Koposov et al. (2010) | $>0.68$ | 13.5 | 15 | Grillmair-Dionatos-1 (GD-1) | green
Law & Majewski (2010) | $0.44$ | 20 | 40 | Sagittarius (Sgr)d | orange
Johnston et al. (2005) | $0.88^{+0.04}_{-0.05}$ | 13 | 50 | Sagittarius (Sgr) | orange
Helmi (2004) | $1.30^{+0.20}_{-0.05}$ | 13 | 50 | Sagittarius (Sgr) | orange
Ibata et al. (2001) | $>0.7$ | 16 | 60 | Sagittarius (Sgr) | orange
Olling & Merrifield (2000) | $0.80^{+0.45}_{-0.30}$ | 1 | 25 | HI gas + disc rotation curve | brown

The various estimates for the shape of the MW vary as widely as shapes of the
simulated galaxies, underlining the difficulty of the measurement. One
positive development from this work is that in most DM models the value of
$s=c/a$ for the DM haloes of the simulated systems appears to be fairly
constant ($\Delta s\lesssim 0.2$) over a wide range of radii ($3$–$100$ kpc),
which should in principle simplify efforts to model the dark halo. We caution,
however, that (1) none of our models include a Large Magellanic Cloud
(LMC)-like companion, which is likely to affect this assertion (e.g. Law &
Majewski, 2010; Vera-Ciro & Helmi, 2013; Vasiliev et al., 2020) and (2) that
this statement assumes that any rotation of the principal axes with radius is
precisely incorporated into the model.

##  5 Discussion

The presence of baryons, and their resulting effect on the shapes of MW-mass
galaxies, shows far wider variety than expected from DM-only and semi-analytic
models. Importantly, MW-mass galaxies in SIDM haloes at the preferred values
of $\sigma/m$ based on studies of dwarf galaxies and galaxy clusters still
have density, velocity, and shape profiles that are consistent with
observations, as well as producing a realistic-looking disc galaxy at the
centre. Thus, there is no immediate discrepancy produced by the introduction
of a nonzero self-interaction cross-section (at least in the velocity-
independent, elastic collision model considered here) that can rule out this
type of SIDM.

More interesting is the question of whether the variation in shape due to a
nonzero cross-section could be constrained well enough to differentiate SIDM
from CDM. From Figure 10 the hope of doing this seems fairly dim, since there
is as much variation in shape from different assembly histories as from
different DM models. However, the importance of the question merits a closer
look at comparisons involving the particular radii where we expect the
differences to be largest.

The key region for looking at shape variations produced by SIDM is likely to
be between about $2$–$20$ kpc, still inside $d_{1}$ (so the self-interactions
have a chance to shape the system), but outside the region where the shape is
utterly dominated by the central galaxy’s baryons. We have the advantage that
in this region we can still tightly constrain the shape of the stellar and gas
components from observations, and look instead at where the shape of the total
mass distribution (constrained using dynamical modeling) departs from the
shape of the stars.

We do this by plotting the difference between the total and stellar axis
ratios $\Delta s=s_{\mathrm{tot}}-s_{\star}$ in Figure 11. The region $r<5$
kpc is dominated by the bulge dynamics and hence varies substantially between
galaxies. However, outside this region we see that there are clear transitions
in the slope for all the curves between the flattened, stellar-disc-dominated
regime (out to around 20 kpc) and the region where the total shape is
determined by the more spherical DM halo. For CDM+Baryon simulations, the
transitions roughly occur at $r_{\star,90}$, represented by the green vertical
lines. For the SIDM+Baryon simulations, the transitions roughly occur at
$d_{1}$, represented by the shaded-areas, following the legend in Figure 4.
While the CDM+Baryon transitions are largely dependent on the stellar mass,
the SIDM+Baryon transitions are instead dependent on the DM self-interaction
cross-section. This holds for $\sigma/m=1$ and 10 cm2 g-1, but not for the
m12m SIDM+Baryon $\sigma/m=0.1$ cm2 g-1, which has no detectable $d_{1}$ at
$z=0$. Since the density for $\sigma/m=0.1$ cm2 g-1 does not reach the levels
needed for local collisions to occur within the age of the halo, this cross-
section behaves similarly to CDM and the shape transition occurs at
$r_{\star,90}$. Likewise, constructing such a test would be difficult near
$\sigma/m=10$ cm2 g-1 in m12m since $d_{1}$ and $r_{\star,90}$ are very close
in this case. Otherwise, we see that for SIDM+Baryon, across all three
simulations, _the shape of the total mass distribution departs from that of
the stellar mass distribution at steadily increasing radius as $\sigma/m$
increases_. Thus, for a given galaxy, constraining the radius of the
transition from where stars dominate its shape to where DM is the dominant
influence and comparing this to (1) the galaxy scale length and (2) the
predicted $d_{1}$ as a function of $\sigma/m$ provides a way to constrain the
SIDM cross-section.

Figure 11: Difference in shape profiles of simulated MW-mass galaxies. Minor-
to-major axis ratio difference $\Delta s=s_{\mathrm{tot}}-s_{\star}$ is shown
as a function of semi-major axis $r=a$ in three sets of MW-mass galaxy
simulations: m12f (left), m12i (centre), and m12m (right). Line-styles and
shaded-areas follow the legend in Figure 4. The black curves represent the
axis ratio difference $\Delta s$ and the vertical lines represent
$r_{\star,90}$, the spherical radius of 90% of stellar mass within 30 kpc of
the central galaxy, given in Table 1.

Currently, there are several methods in use for constraining the shape of the
MW’s total mass distribution using dynamical modeling. A non-exhaustive sample
of measurements using these different methods is listed in Table 2. One
strategy is to model the precession of the orbital planes of tidal streams.
Several attempts to do this for the Sagittarius stream alone (Helmi, 2004;
Johnston et al., 2005; Law et al., 2009; Vera-Ciro & Helmi, 2013; Vasiliev et
al., 2020) have come to inconsistent conclusions, mostly due to differing
parametrizations of the potential and the treatment of the influence of the
LMC on the outer portion of the stream (as pointed out in Vera-Ciro & Helmi
2013 and Erkal et al. 2020). However, since we need mainly to constrain the
region $2\lesssim r\lesssim 20$ kpc for these tests, simultaneous fits of
several streams with smaller apocenters may provide a way forward (Bovy et
al., 2016; Bonaca & Hogg, 2018; Reino et al., 2020).

Another strategy is to constrain the total flattening using Jeans or
Schwarzschild modeling of equilibrium stellar populations as in Loebman et al.
(2014); Hattori et al. (2020), and many other works. However, this is limited
to regions where we have sufficient stellar tracers observed to derive the
distribution functions used in the model. Currently such efforts have been
made in the bulge (to $r\sim 5$ kpc) and in the space observatory Gaia
satellite’s six-dimensional (6D) volume ($5\lesssim r\lesssim 11$ kpc).
However, additional data from ground-based spectroscopic surveys, future Gaia
data releases, and new distance estimators (e.g Das & Sanders, 2019; Auge et
al., 2020) promise to expand the volume accessible to this technique
appreciably in the near future. These new data will also provide excellent
constraints on the _stellar_ shape profile, an equally important quantity in
this approach.

A third strategy, as employed by e.g. Blitz & Spergel (1991); Olling &
Merrifield (2000); Narayan et al. (2005) in the MW and Das et al. (2020) for
external galaxies, is to use the flaring of the HI disc to constrain the
flattening, under the assumption that the gas is in dynamical equilibrium.
Results seem to favor a relatively spherical halo with
$q_{\rho}^{\mathrm{DM}}\sim 0.8$. This technique could provide an independent
assessment with different data and systematics than methods using stellar
kinematics.

Finally, a lower limit on $s=c/a$ as a function of radius may be obtained by
searching for evidence of the truncation or scattering of tidal streams by
orbital resonances, which are more common in more highly flattened potentials
(Valluri et al., 2012; Hesp & Helmi, 2018; Pearson et al., 2015; Valluri et
al., 2020). These effects are quite pronounced at $c/a\sim 0.3$–$0.4$, where
most of our shape profiles begin their transition toward the more spherical
halo, and are much less dominant by $c/a\sim 0.6$–$0.7$, where most of the
total mass profiles end up at large $r$. Looking for an abrupt transition in
the prevalence of streams as a function of galactocentric distance, perhaps
even in a stacked sample of external galaxies where 6D information is not
available, could be an additional way to find constraints on the halo
flattening transition and hence on SIDM.

##  6 Conclusion

We perform a suite of cosmological-baryonic zoom simulations of Milky Way
(MW)-mass galaxies for several different models with self-interacting dark
matter (SIDM), one proposed solution to the challenges of the cold dark matter
(CDM) plus dark energy $\left(\Lambda\text{CDM}\right)$ cosmological model at
small-scales. These dark matter (DM) simulations are compared between CDM and
SIDM (with interaction cross-sections $\sigma/m=0.1$, 1, and 10 cm2 g-1) and
with expectations from the literature.

For the SIDM+Baryon simulations, the variation in axis ratio with SIDM cross-
section is not as large as expected in the literature. The assembly history of
the central galaxy is the dominant influence inside the local collision region
(LCR). Variations in the assembly and evolution of the galaxy thus dominate
the resulting shape. At larger radii, greater differences between CDM and SIDM
axis ratios are also expected according to DM-only models, but again the halo-
to-halo variation in the assembly history of DM+Baryon causes larger
differences in halo shape than the variation in $\sigma/m$ does. In general,
the flattening profile with radius tends to be concave-down (most spherical at
intermediate radii) rather than concave-up (most spherical at small and large
radii) as predicted by previous work.

Although we find that halo-to-halo variation is larger than variation due to
different $\sigma/m$, a possible test for SIDM could lie in predictions of the
difference between the shape of the stellar distribution and that of the total
mass distribution, both of which can be constrained by different methods. The
overall shape of the total distribution (inferred from dynamics) is dominated
by the flattened stellar component inside the LCR and by the more spherical DM
component outside the LCR. The radius of this transition occurs at radius of
90% of stellar mass $r_{\star,90}$ for CDM+Baryon simulations, but increases
with increasing $\sigma/m$ and occurs at the LCR scattering radius $d_{1}$ for
SIDM+Baryon simulations, as the response of the DM to the growing galaxy
becomes more important. There are several promising possibilities for
measuring this radius with new survey data in the MW, and perhaps in other
galaxies, in the coming decade.

Our results at very small radius ($\lesssim 2$ kpc) are limited by the need
for better DM particle resolution for the SIDM-only simulations particularly,
since as predicted the self-interactions significantly heat the central halo
and reduce the particle density, limiting our ability to measure the shape in
this region with our Lagrangian approach to simulations. In the case of all
DM+Baryon simulations, this limitation is mitigated by the presence of the
central galaxy, which deepens the gravitational potential and boosts the DM
density in nearly all cases. However, better resolution in the central part
would enable us to better study the transport of energy and angular momentum
between the stellar and DM components, which will be the subject of future
work.

Finally, we note that the galaxies formed in SIDM haloes with baryonic
feedback differ from those in CDM in mostly subtle ways, and are generally
similar (and consistent with observations) in their large-scale properties
such as mass and scale radius. Therefore, SIDM remains a valid possibility for
new DM physics.

## Acknowledgements

RES acknowledges support from NASA grant 19-ATP19-0068 and HST-AR-15809 from
the Space Telescope Science Institute (STScI), which is operated by AURA,
Inc., under NASA contract NAS5-26555. MBK acknowledges support from NSF CAREER
award AST-1752913, NSF grant AST-1910346, NASA grant NNX17AG29G, and HST-
AR-15006, HST-AR-15809, HST-GO-15658, HST-GO-15901, and HST-GO-15902 from
STScI. AW received support from NASA through ATP grants 80NSSC18K1097 and
80NSSC20K0513; HST grants GO-14734, AR-15057, AR-15809, and GO-15902 from
STScI; a Scialog Award from the Heising-Simons Foundation; and a Hellman
Fellowship. ASG is supported by the Harlan J. Smith postdoctoral fellowship.
This research is part of the Frontera computing project at the Texas Advanced
Computing Center (TACC). Frontera is made possible by National Science
Foundation award OAC-1818253. Simulations in this project were run using Early
Science Allocation 1923870, and analysed using computing resources supported
by the Scientific Computing Core at the Flatiron Institute. This work used
additional computational resources of the University of Texas at Austin and
TACC, the NASA Advanced Supercomputing (NAS) Division and the NASA Center for
Climate Simulation (NCCS), and the Extreme Science and Engineering Discovery
Environment (XSEDE), which is supported by National Science Foundation grant
number OCI-1053575.

## Data Availability

The simulations used for this study are currently proprietary to members of
the FIRE collaboration. Please contact the authors if interested.

## References

  * Ackerman et al. (2009) Ackerman L., Buckley M. R., Carroll S. M., Kamionkowski M., 2009, Phys. Rev. D, 79, 023519
  * Ahmed et al. (2017) Ahmed S. H., Brooks A. M., Christensen C. R., 2017, MNRAS, 466, 3119
  * Ahn & Shapiro (2005) Ahn K., Shapiro P. R., 2005, MNRAS, 363, 1092
  * Allgood et al. (2006) Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS, 367, 1781
  * Arkani-Hamed et al. (2009) Arkani-Hamed N., Finkbeiner D. P., Slatyer T. R., Weiner N., 2009, Phys. Rev. D, 79, 015014
  * Auge et al. (2020) Auge C., et al., 2020, AJ, 160, 18
  * Blitz & Spergel (1991) Blitz L., Spergel D. N., 1991, ApJ, 370, 205
  * Bonaca & Hogg (2018) Bonaca A., Hogg D. W., 2018, ApJ, 867, 101
  * Bondarenko et al. (2020) Bondarenko K., Sokolenko A., Boyarsky A., Robertson A., Harvey D., Revaz Y., 2020, arXiv e-prints, p. arXiv:2006.06623
  * Bovy et al. (2016) Bovy J., Bahmanyar A., Fritz T. K., Kallivayalil N., 2016, ApJ, 833, 31
  * Boylan-Kolchin et al. (2011) Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2011, MNRAS, 415, L40
  * Boylan-Kolchin et al. (2012) Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2012, MNRAS, 422, 1203
  * Brooks & Zolotov (2014a) Brooks A. M., Zolotov A., 2014a, ApJ, 786, 87
  * Brooks & Zolotov (2014b) Brooks A. M., Zolotov A., 2014b, ApJ, 786, 87
  * Brooks et al. (2013) Brooks A. M., Kuhlen M., Zolotov A., Hooper D., 2013, ApJ, 765, 22
  * Bryan & Norman (1998) Bryan G. L., Norman M. L., 1998, ApJ, 495, 80
  * Bullock & Boylan-Kolchin (2017) Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343
  * Butsky et al. (2016) Butsky I., et al., 2016, MNRAS, 462, 663
  * Chan et al. (2015) Chan T. K., Kereš D., Oñorbe J., Hopkins P. F., Muratov A. L., Faucher-Giguère C. A., Quataert E., 2015, MNRAS, 454, 2981
  * Conn et al. (2013) Conn A. R., et al., 2013, ApJ, 766, 120
  * Correa (2020) Correa C. A., 2020, arXiv e-prints, p. arXiv:2007.02958
  * Creasey et al. (2017) Creasey P., Sameie O., Sales L. V., Yu H.-B., Vogelsberger M., Zavala J., 2017, MNRAS, 468, 2283
  * Das & Sanders (2019) Das P., Sanders J. L., 2019, MNRAS, 484, 294
  * Das et al. (2020) Das M., McGaugh S. S., Ianjamasimanana R., Schombert J., Dwarakanath K. S., 2020, ApJ, 889, 10
  * Debattista et al. (2013) Debattista V. P., Roškar R., Valluri M., Quinn T., Moore B., Wadsley J., 2013, MNRAS, 434, 2971
  * Debattista et al. (2019) Debattista V. P., Gonzalez O. A., Sand erson R. E., El-Badry K., Garrison-Kimmel S., Wetzel A., Faucher-Giguère C.-A., Hopkins P. F., 2019, MNRAS, 485, 5073
  * Despali et al. (2019) Despali G., Sparre M., Vegetti S., Vogelsberger M., Zavala J., Marinacci F., 2019, MNRAS, 484, 4563
  * Dubinski & Carlberg (1991) Dubinski J., Carlberg R. G., 1991, ApJ, 378, 496
  * Dutton et al. (2016) Dutton A. A., Macciò A. V., Frings J., Wang L., Stinson G. S., Penzo C., Kang X., 2016, MNRAS, 457, L74
  * El-Badry et al. (2016) El-Badry K., Wetzel A., Geha M., Hopkins P. F., Kereš D., Chan T. K., Faucher-Giguère C.-A., 2016, ApJ, 820, 131
  * El-Badry et al. (2018) El-Badry K., et al., 2018, MNRAS, 473, 1930
  * Elbert et al. (2018) Elbert O. D., Bullock J. S., Kaplinghat M., Garrison-Kimmel S., Graus A. S., Rocha M., 2018, ApJ, 853, 109
  * Erkal et al. (2020) Erkal D., Belokurov V. A., Parkin D. L., 2020, MNRAS, 498, 5574
  * Feng et al. (2009) Feng J. L., Kaplinghat M., Tu H., Yu H.-B., 2009, J. Cosmology Astropart. Phys., 2009, 004
  * Franx et al. (1991) Franx M., Illingworth G., de Zeeuw T., 1991, ApJ, 383, 112
  * Fritz et al. (2018) Fritz T. K., Battaglia G., Pawlowski M. S., Kallivayalil N., van der Marel R., Sohn S. T., Brook C., Besla G., 2018, A&A, 619, A103
  * Garrison-Kimmel et al. (2014) Garrison-Kimmel S., Boylan-Kolchin M., Bullock J. S., Kirby E. N., 2014, MNRAS, 444, 222
  * Garrison-Kimmel et al. (2017) Garrison-Kimmel S., et al., 2017, MNRAS, 471, 1709
  * Garrison-Kimmel et al. (2018) Garrison-Kimmel S., et al., 2018, MNRAS, 481, 4133
  * Garrison-Kimmel et al. (2019) Garrison-Kimmel S., et al., 2019, MNRAS, 487, 1380
  * Governato et al. (2012) Governato F., et al., 2012, MNRAS, 422, 1231
  * Hattori et al. (2020) Hattori K., Valluri M., Vasiliev E., 2020, arXiv e-prints, p. arXiv:2012.03908
  * Hayashi et al. (2020) Hayashi K., Ibe M., Kobayashi S., Nakayama Y., Shirai S., 2020, arXiv e-prints, p. arXiv:2008.02529
  * Helmi (2004) Helmi A., 2004, ApJ, 610, L97
  * Hesp & Helmi (2018) Hesp C., Helmi A., 2018, arXiv e-prints, p. arXiv:1804.03670
  * Hopkins et al. (2018) Hopkins P. F., et al., 2018, MNRAS, 480, 800
  * Ibata et al. (2001) Ibata R., Lewis G. F., Irwin M., Totten E., Quinn T., 2001, ApJ, 551, 294
  * Ibata et al. (2013) Ibata R. A., et al., 2013, Nature, 493, 62
  * Ibata et al. (2014) Ibata R. A., Ibata N. G., Lewis G. F., Martin N. F., Conn A., Elahi P., Arias V., Fernando N., 2014, ApJ, 784, L6
  * Jin et al. (2020) Jin Y., Zhu L., Long R. J., Mao S., Wang L., van de Ven G., 2020, MNRAS, 491, 1690
  * Johnston et al. (2005) Johnston K. V., Law D. R., Majewski S. R., 2005, ApJ, 619, 800
  * Kado-Fong et al. (2020) Kado-Fong E., Greene J. E., Huang S., Beaton R., Goulding A. D., Komiyama Y., 2020, ApJ, 900, 163
  * Kamada et al. (2017) Kamada A., Kaplinghat M., Pace A. B., Yu H.-B., 2017, Phys. Rev. Lett., 119, 111102
  * Kaplinghat et al. (2014) Kaplinghat M., Keeley R. E., Linden T., Yu H.-B., 2014, Phys. Rev. Lett., 113, 021302
  * Kaplinghat et al. (2019) Kaplinghat M., Valli M., Yu H.-B., 2019, MNRAS, 490, 231
  * Kim et al. (2014) Kim J.-h., et al., 2014, ApJS, 210, 14
  * Klypin et al. (1999) Klypin A., Kravtsov A. V., Valenzuela O., Prada F., 1999, ApJ, 522, 82
  * Koda & Shapiro (2011) Koda J., Shapiro P. R., 2011, MNRAS, 415, 1125
  * Koposov et al. (2010) Koposov S. E., Rix H.-W., Hogg D. W., 2010, ApJ, 712, 260
  * Kroupa et al. (2005) Kroupa P., Theis C., Boily C. M., 2005, A&A, 431, 517
  * Küpper et al. (2015) Küpper A. H. W., Balbinot E., Bonaca A., Johnston K. V., Hogg D. W., Kroupa P., Santiago B. X., 2015, ApJ, 803, 80
  * Kuzio de Naray et al. (2008) Kuzio de Naray R., McGaugh S. S., de Blok W. J. G., 2008, ApJ, 676, 920
  * Law & Majewski (2010) Law D. R., Majewski S. R., 2010, ApJ, 714, 229
  * Law et al. (2009) Law D. R., Majewski S. R., Johnston K. V., 2009, ApJ, 703, L67
  * Li et al. (2018) Li H., Mao S., Cappellari M., Graham M. T., Emsellem E., Long R. J., 2018, ApJ, 863, L19
  * Libeskind et al. (2009) Libeskind N. I., Frenk C. S., Cole S., Jenkins A., Helly J. C., 2009, MNRAS, 399, 550
  * Loeb & Weiner (2011) Loeb A., Weiner N., 2011, Phys. Rev. Lett., 106, 171302
  * Loebman et al. (2014) Loebman S. R., et al., 2014, ApJ, 794, 151
  * Lynden-Bell (1976) Lynden-Bell D., 1976, MNRAS, 174, 695
  * Malhan & Ibata (2019) Malhan K., Ibata R. A., 2019, MNRAS, 486, 2995
  * Mashchenko et al. (2008) Mashchenko S., Wadsley J., Couchman H. M. P., 2008, Science, 319, 174
  * Moore (1994) Moore B., 1994, Nature, 370, 629
  * Moore et al. (1999) Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P., 1999, ApJ, 524, L19
  * Müller et al. (2018) Müller O., Pawlowski M. S., Jerjen H., Lelli F., 2018, Science, 359, 534
  * Narayan et al. (2005) Narayan C. A., Saha K., Jog C. J., 2005, A&A, 440, 523
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  * Oñorbe et al. (2015) Oñorbe J., Boylan-Kolchin M., Bullock J. S., Hopkins P. F., Kereš D., Faucher-Giguère C.-A., Quataert E., Murray N., 2015, MNRAS, 454, 2092
  * Oh et al. (2008) Oh S.-H., de Blok W. J. G., Walter F., Brinks E., Kennicutt Robert C. J., 2008, AJ, 136, 2761
  * Oh et al. (2015) Oh S.-H., et al., 2015, AJ, 149, 180
  * Olling & Merrifield (2000) Olling R. P., Merrifield M. R., 2000, MNRAS, 311, 361
  * Oman et al. (2015) Oman K. A., et al., 2015, MNRAS, 452, 3650
  * Palau & Miralda-Escudé (2019) Palau C. G., Miralda-Escudé J., 2019, MNRAS, 488, 1535
  * Palau & Miralda-Escudé (2020) Palau C. G., Miralda-Escudé J., 2020, arXiv e-prints, p. arXiv:2010.14381
  * Pawlowski & Kroupa (2020) Pawlowski M. S., Kroupa P., 2020, MNRAS, 491, 3042
  * Pawlowski et al. (2012) Pawlowski M. S., Pflamm-Altenburg J., Kroupa P., 2012, MNRAS, 423, 1109
  * Peñarrubia et al. (2012) Peñarrubia J., Pontzen A., Walker M. G., Koposov S. E., 2012, ApJ, 759, L42
  * Pearson et al. (2015) Pearson S., Küpper A. H. W., Johnston K. V., Price-Whelan A. M., 2015, ApJ, 799, 28
  * Peter et al. (2013a) Peter A. H. G., Rocha M., Bullock J. S., Kaplinghat M., 2013a, MNRAS, 430, 105
  * Peter et al. (2013b) Peter A. H. G., Rocha M., Bullock J. S., Kaplinghat M., 2013b, MNRAS, 430, 105
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A16
  * Pontzen & Governato (2012) Pontzen A., Governato F., 2012, MNRAS, 421, 3464
  * Posti & Helmi (2019) Posti L., Helmi A., 2019, A&A, 621, A56
  * Randall et al. (2008) Randall S. W., Markevitch M., Clowe D., Gonzalez A. H., Bradač M., 2008, ApJ, 679, 1173
  * Reino et al. (2020) Reino S., Rossi E. M., Sanderson R. E., Sellentin E., Helmi A., Koppelman H. H., Sharma S., 2020, arXiv e-prints, p. arXiv:2007.00356
  * Ren et al. (2019) Ren T., Kwa A., Kaplinghat M., Yu H.-B., 2019, Physical Review X, 9, 031020
  * Riley & Strigari (2020) Riley A. H., Strigari L. E., 2020, MNRAS, 494, 983
  * Robertson et al. (2019) Robertson A., Harvey D., Massey R., Eke V., McCarthy I. G., Jauzac M., Li B., Schaye J., 2019, MNRAS, 488, 3646
  * 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. (2018) Sameie O., Creasey P., Yu H.-B., Sales L. V., Vogelsberger M., Zavala J., 2018, MNRAS, 479, 359
  * Sameie et al. (2019) Sameie O., Benson A. J., Sales L. V., Yu H.-b., Moustakas L. A., Creasey P., 2019, ApJ, 874, 101
  * Sameie et al. (2020a) Sameie O., Chakrabarti S., Yu H.-B., Boylan-Kolchin M., Vogelsberger M., Zavala J., Hernquist L., 2020a, arXiv e-prints, p. arXiv:2006.06681
  * Sameie et al. (2020b) Sameie O., Yu H.-B., Sales L. V., Vogelsberger M., Zavala J., 2020b, Phys. Rev. Lett., 124, 141102
  * Sameie et al. (2021) Sameie O., Boylan-Kolchin M., Sanderson R., Vargya D., Hopkins P., Wetzel A., Bullock J., Graus A., 2021, arXiv e-prints, p. arXiv:2102.12480
  * Samuel et al. (2020a) Samuel J., Wetzel A., Chapman S., Tollerud E., Hopkins P. F., Boylan-Kolchin M., Bailin J., Faucher-Giguère C.-A., 2020a, arXiv e-prints, p. arXiv:2010.08571
  * Samuel et al. (2020b) Samuel J., et al., 2020b, MNRAS, 491, 1471
  * Samuel et al. (2021) Samuel J., Wetzel A., Chapman S., Tollerud E., Hopkins P. F., Boylan-Kolchin M., Bailin J., Faucher-Giguére C.-A., 2021, MNRAS,
  * Sanderson et al. (2018) Sanderson R. E., et al., 2018, ApJ, 869, 12
  * Sanderson et al. (2020) Sanderson R. E., et al., 2020, ApJS, 246, 6
  * Shao et al. (2019) Shao S., Cautun M., Frenk C. S., 2019, MNRAS, 488, 1166
  * Spergel & Steinhardt (2000) Spergel D. N., Steinhardt P. J., 2000, Phys. Rev. Lett., 84, 3760
  * Tollet et al. (2016) Tollet E., et al., 2016, MNRAS, 456, 3542
  * 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
  * Valli & Yu (2018) Valli M., Yu H.-B., 2018, Nature Astronomy, 2, 907
  * Valluri et al. (2012) Valluri M., Debattista V. P., Quinn T. R., Roškar R., Wadsley J., 2012, MNRAS, 419, 1951
  * Valluri et al. (2020) Valluri M., Price-Whelan A. M., Snyder S. J., 2020, arXiv e-prints, p. arXiv:2009.09004
  * Vasiliev et al. (2020) Vasiliev E., Belokurov V., Erkal D., 2020, arXiv e-prints, p. arXiv:2009.10726
  * Vera-Ciro & Helmi (2013) Vera-Ciro C., Helmi A., 2013, ApJ, 773, L4
  * Vera-Ciro et al. (2011) Vera-Ciro C. A., Sales L. V., Helmi A., Frenk C. S., Navarro J. F., Springel V., Vogelsberger M., White S. D. M., 2011, MNRAS, 416, 1377
  * Vogelsberger et al. (2012) Vogelsberger M., Zavala J., Loeb A., 2012, MNRAS, 423, 3740
  * Vogelsberger et al. (2016) Vogelsberger M., Zavala J., Cyr-Racine F.-Y., Pfrommer C., Bringmann T., Sigurdson K., 2016, MNRAS, 460, 1399
  * Walker & Peñarrubia (2011) Walker M. G., Peñarrubia J., 2011, ApJ, 742, 20
  * Wegg et al. (2019) Wegg C., Gerhard O., Bieth M., 2019, MNRAS, 485, 3296
  * Wetzel et al. (2016) Wetzel A. R., Hopkins P. F., Kim J.-h., Faucher-Giguère C.-A., Kereš D., Quataert E., 2016, ApJ, 827, L23
  * Yoshida et al. (2000) Yoshida N., Springel V., White S. D. M., Tormen G., 2000, ApJ, 544, L87
  * Zavala et al. (2013) Zavala J., Vogelsberger M., Walker M. G., 2013, MNRAS, 431, L20

##  Appendix A Ellipsoidal Versus Spherical Shells

As discussed in §4, we use the triaxial ellipsoidal shape profiles to
calculate density and velocity profiles in ellipsoidal shells. However, past
work has largely relied on spherical shells to calculate these profiles.
Figure 12 shows the difference in DM local collision rate
$\left[\Gamma_{\mathrm{scatter,E}}(d)-\Gamma_{\mathrm{scatter,S}}(r)\right]/\Gamma_{\mathrm{scatter,E}}(d)$
between ellipsoidal shells and spherical shells, where
$\Gamma_{\mathrm{scatter,E}}(d)$ and $\Gamma_{\mathrm{scatter,S}}(r)$ are the
$\Gamma_{\mathrm{scatter}}$ for ellipsoidal and spherical methods,
respectively. We note that the ellipsoidal shells have larger
$\Gamma_{\mathrm{scatter}}$ due to having larger estimates of the DM density,
particularly at larger $d$. This is the effect of “smearing” across isodensity
contours when using spherical shells to compute densities. The DM velocity
profiles have negligible differences between the ellipsoidal and spherical
methods, suggesting that the local velocity ellipsoid is relatively isotropic.

We calculate the difference in scattering radius when using ellipsoidal shells
($d_{1}$) rather than spherical shells ($r_{1}$) as $\Delta_{1}=d_{1}-r_{1}$,
where $\Delta_{1}=2.$ kpc for SIDM-only $\sigma/m=10$ cm2 g-1. For SIDM+Baryon
$\sigma/m=1$ and 10 cm2 g-1, the differences are $\Delta_{1}=0.6$–$1.6$ kpc
and $\Delta_{1}=1.$–$2.$ kpc, respectively. The $\Delta_{1}\sim 5$–$20$%
difference demonstrates the importance of using ellipsoidal shells fit to the
isodensity contours to estimate density and velocity profiles. The use of
ellipsoidal shells is even more important for the stars, gas, and total mass
distributions, since these have less spherical shapes, and thus larger density
and velocity differences between the two methods compared to the DM component.

Figure 12: Difference in DM local collision rate of simulated MW-mass
galaxies. The difference when averaging over shells of ellipsoidal distances
$d$ versus spherical radii $r=a$,
$\left[\Gamma_{\mathrm{scatter,E}}(d)-\Gamma_{\mathrm{scatter,S}}(r)\right]/\Gamma_{\mathrm{scatter,E}}(d)$
(Equation 6) as a function of $d=r=a$ for SIDM simulations. Shells are matched
so that spherical radius $r$ is equal to the semi-major axis $a$ of the
ellipsoid used to compute $d$. Line-styles and colours follow the legend in
Figure 4. Using ellipsoidal shells leads to larger
$\Gamma_{\mathrm{scatter}}$, mostly because of an increase in the computed DM
density. The difference in the local collision rate increases for the
ellipsoidal shells method away from the centre of the galaxies.

##  Appendix B Geometrical mean of axis lengths

Figure 13 shows the axis ratio $s=c/a$ versus the geometrical mean of the axis
ratios $R=\left(abc\right)^{1/3}$, instead of semi-major axis distance $r=a$
as in Figure 8. This plot is included to facilitate comparison with previous
work on halo shapes. The most important difference to note relative to Figure
8 is that the axis ratio curves have shifted towards smaller radii, since the
ellipsoidal geometrical mean radius $R$ is always less than or equal to semi-
major axis $r=a$ by definition. Due to the relatively spherical shape of the
DM component, its curve is shifted the least, while the profile of the gas
component, which is the most flattened, has shifted the most. Rapid changes in
the flattening of neighboring ellipsoids have the effect of producing non-
functional curves, such as observed for the gas in m12m SIDM+Baryon
$\sigma/m=1$ cm2 g-1 near $R\sim 1$ kpc. This is a result of the fact that the
definition of $R$ does not guarantee that it must always increase with
increasing ellipsoid semi-major axis $r=a$.

Figure 13: Shape profiles of simulated MW-mass galaxies as a function of $R$.
Minor-to-major axis ratio $s=c/a$ as a function of geometrical mean radius $R$
for different species (DM, stars, and gas) in three sets of MW-mass galaxy
simulations: m12f (left), m12i (centre), and m12m (right). Line-styles,
colours, and shaded-areas follow the legend in Figure 4. An alternate version
of this plot using the semi-major axis $r=a$ can be found in the main body of
the work in Figure 8. The literature often plots axis ratios versus this
geometrical mean of the axis lengths $R$, which is therefore given here for
comparison purposes. We note that the main difference is in how the curves
have shifted leftward towards smaller radii.

Similarly, Figure 14 shows the triaxiality $T$ versus the geometrical mean of
the axis lengths $R$, instead of semi-major axis $r=a$. The plot shows a
similar shift leftward for all the triaxiality curves, but the shift is less
pronounced since the triaxiality depends on all three axis lengths $a$, $b$,
and $c$. This shift does in general reduce the apparent prolateness of the
CDM+Baryon and SIDM+Baryon MW-mass galaxies at small radii, especially for
m12f and m12i.

Figure 14: Triaxiality profiles of simulated MW-mass galaxies as a function of
$R$. Triaxiality $T$ as a function of geometrical mean radius $R$ for
different species (DM and stars) in three sets of MW-mass galaxy simulations:
m12f (left), m12i (centre), and m12m (right). Line-styles, colours, and
shaded-areas follow the legend in Figure 4. An alternate version of this plot
using the semi-major axis $r=a$ can be found in the main body of the work in
Figure 9. The literature often plots triaxiality versus this geometrical mean
of the axis lengths $R$, which is therefore given here for comparison
purposes. We note that the main difference is in how the curves have shifted
leftward towards smaller radii.

Overall, we consider that the use of semi-major axis length $r=a$ is
preferable to the use of the geometrical mean of the axis lengths $R$, since
observations of galaxies are 3D light and velocity (redshift) distributions
projected as two-dimensional (2D) light and velocity distributions on to the
celestial sphere. Therefore, measuring a realistic $R$ for any galaxy is a
difficult and degenerate task. On the other hand, deprojecting _only_ the
semi-major and semi-minor axis-lengths of nested ellipsoids is a more feasible
task, since measurements of the line-of-sight velocity distribution (e.g. Li
et al., 2018; Jin et al., 2020) or stacked observations (e.g. Kado-Fong et
al., 2020) can be used to constrain or marginalize over the inclination angle.
This motivates the use of $r=a$, inside which density, velocity, and shape
profiles of the baryons (and thus also DM) are then more easily estimated.