The imprint of dark matter on the Galactic acceleration field

Abstract

Measurements of the accelerations of stars enabled by time-series extreme-precision spectroscopic observations, from pulsar timing, and from eclipsing binary stars in the Solar Neighborhood offer insights into the mass distribution of the Milky Way that do not rely on traditional equilibrium modeling. Given the measured accelerations, we can determine a total mass density, and from this, by accounting for the mass in stars, gas, and dust, we can infer the amount of dark matter. Leveraging the FIRE-2 simulations of Milky Way-mass galaxies, we compare vertical acceleration profiles between cold dark matter (CDM) and self-interacting dark matter (SIDM) with constant cross-section of 1 cm² g^-1 across three halos with diverse assembly histories. Notably, significant asymmetries in vertical acceleration profiles near the midplane at fixed radii are observed in both CDM and SIDM, particularly in halos recently affected by mergers with satellites of Sagittarius/SMC-like masses or greater. These asymmetries offer a unique window into exploring the merger history of a galaxy. We show that SIDM halos consistently exhibit higher local stellar and dark matter densities and steeper vertical acceleration gradients, up to 30% steeper near the Solar Neighborhood. SIDM halos also manifest a more oblate halo shape in the Solar Neighborhood. Furthermore, enhanced precision in acceleration measurements and larger datasets promise to provide better constraints on the local dark matter density, complementing our understanding from kinematic analysis of their distribution within galaxies.

1 Introduction

Dark matter (DM) constitutes approximately 85% of the matter in the universe (Planck et al., 2020), but its nature remains elusive. The standard cold dark matter (CDM) model, which assumes that DM particles are collisionless and non-interacting, has been successful in explaining large-scale structures in the Universe. However, there are several inconsistencies at small scales, such as the too-big-to-fail tension, core-cusp tension, planes-of-satellites, and others (e.g. Flores & Primack, 1994; Moore, 1994; Moore et al., 1999; Bullock & Boylan-Kolchin, 2017; Tulin & Yu, 2018; Sales et al., 2022).

One unresolved issue is the “diversity of rotation curves”, where the observed rotation curves for dwarf satellites and Milky Way-mass galaxies show a stunning diversity, with inferred inner DM profiles ranging from cored to NFW-like and even more concentrated than NFW profiles (e.g. Oman et al., 2015; Zavala et al., 2019). On the other hand, CDM simulations predict a steeply rising central density profile, and reproducing such diversity in inner DM profiles remains a challenge even with additions of baryonic components with active feedback (Walker & Penarrubia, 2011; Santos-Santos et al., 2020; Ebisu et al., 2022).

This discrepancy has spurred exploration into alternative models, such as self-interacting dark matter (SIDM) (Carlson et al., 1992; Spergel & Steinhardt, 2000), which allows DM particles to interact with each other through elastic collisions. Such interactions can lead to a transfer of kinetic energy from the hot outer region to the dense inner region, resulting in a more cored density profile (Spergel & Steinhardt, 2000; Yoshida et al., 2000; Colín et al., 2002; Tulin & Yu, 2018). Moreover, a large self-interaction cross-section ($\sigma/m\leq 10$ cm² g^-1) can induce a core-collapse phase that drives the DM mass towards the center and forms steeper inner DM profiles (Kahlhoefer et al., 2019; Zavala et al., 2019; Sameie et al., 2020).

Recent studies have shown that SIDM can reproduce some observed properties of dwarf galaxies and Milky Way (MW)-mass galaxies better than CDM (Vogelsberger et al., 2012; Rocha et al., 2013; Peter et al., 2013; Vogelsberger et al., 2019), such as their shapes (Sameie et al., 2018, 2021; Vargya et al., 2022) and inner densities (Kaplinghat et al., 2016). One effective test of DM involves measuring the DM mass distribution within a galaxy.

The most straightforward and least assumptive method for probing the mass distribution (including stars, gas and DM) in the Galaxy is through acceleration measurements of stars in the MW. Recent advances in high-precision spectrographs (e.g., Schwab et al., 2016; Fischer et al., 2016; Wright & Robertson, 2017; Pepe et al., 2021), improving precision in pulsar timing data (e.g., Keith et al., 2013; Pennucci, 2019; Goncharov et al., 2021), and high-precision spectroscopic observations of eclipsing binary stars by space-based missions (Hełminiak et al., 2019) now enable direct measurement of Galactic accelerations in the Solar Neighborhood with multiple independent techniques (e.g Silverwood & Easther, 2019; Chakrabarti et al., 2020, 2021; Phillips et al., 2021; Chakrabarti et al., 2022, 2022).

Acceleration measurements can be used to determine the shape of the MW potential in the Solar Neighborhood and constrain the Galactic mid-plane density (Oort limit, Chakrabarti et al., 2020, 2021; Donlon et al., 2024), from which we can determine the local DM density after accounting for the baryon budget. Notably, in comparison to alternative methodologies (reviewed comprehensively by de Salas & Widmark, 2021), acceleration measurements require far fewer assumptions. Using Poisson’s equation, given an acceleration field (i.e., the gradient of the potential), we can straightforwardly determine the local density.

These properties of the MW depend on the DM model (Sameie et al., 2018; Vargya et al., 2022; Sameie et al., 2021), indicating that measurement of the local dark matter density and the shape of the potential can provide discriminating power for constraining the nature of DM. Prior work on determinations of the Oort limit has focused on kinematic analysis (e.g. McKee et al., 2015; Schutz et al., 2018; Guo et al., 2020), which models a snapshot in time of the positions and speeds of stars under simplifying assumptions of equilibrium and symmetry both across the midplane and axisymmetrically, which can lead to inaccurate inferences of Galactic parameters for a time-dependent potential (Haines et al., 2019).

Chakrabarti et al. (2020) demonstrated that perturbations from past mergers such as the Gaia-Sausage-Enceladus merger (Helmi et al., 2018; Belokurov et al., 2018), the Sagittarius dwarf galaxy (Ibata et al., 1994; Johnston et al., 1995; Newberg et al., 2002), the ongoing merger with the LMC (Besla et al., 2007; Kallivayalil et al., 2013) and other possible disturbances in the MW such as formation of a warped disk (Ostriker & Binney, 1989; Evans et al., 1998) or large planar disturbances in the outer gas disk of the Galaxy (Chakrabarti & Blitz, 2009; Chakrabarti et al., 2019) can induce asymmetries in the Galactic acceleration profile. These non-equilibrium effects are naturally taken into account while studying DM through direct measurement of Galactic accelerations, because extreme-precision time-series observations of stars provide the Galactic acceleration today without assuming equilibrium or symmetry as in kinematic analyses. Moreover, kinematic analyses only yield the average acceleration for a bulk population of stars, which offers less constraining power.

Cosmological simulations with a realistic baryonic disk and assembly history can serve as a natural test laboratory for interrogating tools to measure local DM density (e.g. Necib et al., 2019; Ou et al., 2024) and for interpreting Galactic acceleration observations (Loebman et al., 2012, 2014). Moreover, by examining these simulations across different DM models, we can effectively constrain the mass distribution and nature of DM.

In this paper, we use cosmological baryonic simulations with varying merger histories described in Sec. 2 to analyze their Galactic acceleration fields (Sec. 3) under different DM models. In Sec. 2.2, we compare the measured shape of the MW potential from pulsar timing with our simulations. Here, we compare to fundamental Galactic parameters determined from the time-rate of change of the orbital period of binary pulsars (Chakrabarti et al., 2021; Donlon et al., 2024) as using the time-rate of change of the spin-period (Phillips et al., 2021) leads to large uncertainties due to the unknown effect of magnetic braking on the spin periods. Our conclusions are presented in Sec. 4.

2 Simulations of MW-mass galaxies

We use three MW–mass galaxies using cosmological-baryonic simulations from the Latte suite (Wetzel et al., 2023) of FIRE-2 simulations with initial conditions derived from AGORA (Kim et al., 2014) run with GIZMO (Hopkins, 2015). The halos are labelled m12f, m12i, and m12m are some of the most MW-like isolate disks in the suite and span a range of assembly histories. All the halos use identical baryonic FIRE-2 physics (Hopkins et al., 2018) and two distinct DM models: CDM, which employs cold-collisionless DM and SIDM, which adopts self-interacting DM with a cross-section of $\sigma/m=1$ cm² g^-1. The SIDM implementation in GIZMO (introduced in Rocha et al. 2013 and Peter et al. 2013) employs a Monte Carlo approach to determine the scattering probability for the nearest neighbors of each DM particle using a spline kernel with adaptive smoothing length (Monaghan & Lattanzio, 1985). It then assigns velocities isotropically to the scattered particles to conserve energy and momentum. The cross-section of $\sigma/m=1$ cm² g^-1 for SIDM is often used for MW-mass galaxies to achieve a balance between addressing small-scale problems (Tulin & Yu, 2018) and maintaining agreement with larger-scale galactic and cluster observations (Randall et al., 2008).

These halos are all isolated systems with no massive companion halos within 4 Mpc, and have a virial mass of approximately $\sim 1-1.5\times 10^{12}$ $\textrm{M}_{\odot}$. Each simulation uses an initial particle mass of $m_{\mathrm{b}}=7100$ M_$\odot$ for stars and gas, and $m_{\mathrm{DM}}=35000$ M_$\odot$ for DM with a minimum physical spatial resolution $\epsilon_{\textrm{gas}}=1$ pc, $\epsilon_{\star}=4$ pc, and $\epsilon_{\textrm{DM}}=20$ pc. Also, the potential and acceleration values are tracked for every particle in the simulation.

The properties of the halos, including halo shapes and density profiles, along with a description of the simulation methods are detailed in Sameie et al. (2021) and Vargya et al. (2022). These halos are part of a suite of simulations dedicated to exploring alternative DM models which are different from the fiducial FIRE-2 simulations presented in Wetzel et al. (2023). In these simulations, any thermal energy to momentum for stellar winds energy conversion resulting from stellar mass-loss processes is neglected for the sub-resolution regions which leads to a lower (almost a factor of 2) stellar mass than our fiducial FIRE-2 simulations.

Notably, FIRE-2 CDM halos exhibit MW–like stellar-to-halo mass ratios (Hopkins et al., 2018), stellar morphologies and kinematics (Ma et al., 2017; McCluskey et al., 2024), and other properties such as gas fractions, scale radii, scale heights, and satellite populations (e.g. El-Badry et al., 2018; Escala et al., 2018; Sanderson et al., 2018; Samuel et al., 2021), making them ideal for comparing with observational Galactic accelerations in the Solar Neighborhood.

Additionally, the three halos have distinct mergers and assembly histories. The simulation m12i features a thick young disk with an intermediate formation epoch (Sanderson et al., 2020) with a merger with first pericentric passage around 6 Gyr before the present day. m12f had a major merger with a pair of satellites similar in total mass to the progenitor of Sagittarius stream (Jiang & Binney, 2000; Niederste-Ostholt et al., 2010; De Boer et al., 2015; Gibbons et al., 2017) with the first pericentric passage about 3 Gyr before the present day (Arora et al., 2022; Garavito-Camargo et al., 2023). In contrast, m12m is the earliest forming, has a massive disk, and has no massive mergers in the last 7 Gyr of its evolution (Debattista et al., 2019). Table LABEL:tab:gal_prop lists the total stellar mass within 10% of the virial radius, where the cumulative density is 200 times the critical density of the Universe ($\textrm{M}_{\star}$) and the 3D stellar half-mass radius ($\textrm{r}_{1/2}$) for all of the galaxies. While the total virial mass of these systems is the same in the SIDM and CDM runs, the SIDM simulations have systematically higher stellar mass throughout the disk due to increased star formation rates at late times (Sameie et al., 2021) and larger galaxy sizes, quantified by $\textrm{r}_{1/2}$, compared to the CDM simulations.

2.1 The Solar Neighborhood

Figure 1: Relative variation in surface density ($\delta_{\textrm{dens}}$) in the XY plane around the Solar Neighborhood within $\textrm{R}_{\odot}=8.1\pm 0.6$ kpc and $|Z|<5$ kpc for m12f halos; i. SIDM (top), and ii. CDM (bottom) computed for particles species (all–left, baryons–middle, DM–right). Most variation with overdensities (red) and underdensities (blue) are in the baryonic component, while DM densities are relatively uniform and smooth.

Table 1: Properties of simulated galaxies and their azimuthal variation at the Solar Circle, $\textrm{R}_{\odot}=8.1\pm 0.6$ kpc, comparing CDM and SIDM models with the MW at present day. $\textrm{M}_{\star}$: Stellar mass within 10% of the radius where the cumulative density is 200 times the critical density of the Universe. $\textrm{r}_{1/2}$: 3D stellar half-mass radius.

halo property

$\textrm{surface density}^{\textrm{c}}$ [$\textrm{M}_{\odot}/\textrm{pc}^{-2}$]

$\textrm{volume density}^{\textrm{d}}$ [$10^{-3}\textrm{M}_{\odot}/\textrm{pc}^{-3}$]

scale height [pc]

Galaxy

Physics

$\textrm{M}_{\star}^{\textrm{b}}[\textrm{M}_{\odot}]$

$\textrm{r}_{1/2}^{\textrm{b}}[\textrm{kpc}]$

total

baryonic

total

stars

gas

DM

stars thin

stars thick

$\textrm{cold}^{*}$ gas

$\textrm{MW}^{\textrm{a}}$

?

${5^{+0.4}_{-0.5}\times 10^{10}}$

${4.2^{+1}_{-1}}$

70${}^{+5}_{-5}$

43.8${}^{+3.4}_{-3.4}$

99.5${}^{+9.5}_{-9.5}$

42.4${}^{+3}_{-3}$

45.9${}^{+8}_{-8}$

10.5${}^{+2.5}_{-2.5}$

300${}^{+60}_{-60}$

900${}^{+180}_{-180}$

150

$m12i$

CDM

$3.4\times 10^{10}$

4.3

44.7${}^{+24}_{-5}$

25.3${}^{+24}_{-4}$

23.1${}^{+15}_{-11}$

7.8${}^{+6.3}_{-4.5}$

6.6${}^{+8.2}_{-5.8}$

8.7${}^{+0.2}_{-0.2}$

566.4

$2000^{\dagger}$

321.5

SIDM

$5.2\times 10^{10}$

3.8

55.9${}^{+24}_{-17}$

36.8${}^{+23}_{-17}$

30.3${}^{+18}_{-15}$

11.2${}^{+5.2}_{-4.5}$

10.4${}^{+13}_{-10}$

9.1${}^{+0.1}_{-0.1}$

488

1503.4

359.1

$m12f$

CDM

$5.8\times 10^{10}$

3.6

61.1${}^{+20}_{-13}$

39.2${}^{+19}_{-11}$

33.5${}^{+19}_{-11}$

12.9${}^{+4.3}_{-3.1}$

10.4${}^{+15.2}_{-7.6}$

10.1${}^{+0.3}_{-0.3}$

569.9

$2000^{\dagger}$

363.4

SIDM

$6.7\times 10^{10}$

4.7

67.9${}^{+20}_{-16}$

43.4${}^{+20}_{-16}$

35.7${}^{+22}_{-10}$

16.7${}^{+5.5}_{-3.2}$

7.6${}^{+16}_{-6.4}$

11.4${}^{+0.3}_{-0.1}$

437.6

1078.6

322.8

$m12m$

CDM

$5.1\times 10^{10}$

8.2

63.3${}^{+19}_{-15}$

39.6${}^{+18}_{-15}$

35${}^{+20}_{-10}$

16${}^{+5.7}_{-2.8}$

8.3${}^{+13.2}_{-7.3}$

10.6${}^{+0.2}_{-0.2}$

153.3

682

297.6

SIDM

$6.9\times 10^{10}$

8.3

87.9${}^{30}_{-26}$

57${}^{+29}_{-25}$

49.5${}^{+28}_{-13}$

25${}^{+8.6}_{-5.4}$

10.7${}^{+19}_{-7.4}$

13.9${}^{+0.3}_{-0.3}$

218.4

742.9

317.2

Given the cosmological nature of these simulations in an arbitrary comoving box, we establish the galactocentric coordinates for each galaxy through a two-step process. First, we use the “shrinking spheres” method (Power et al., 2003) to determine the center position of each galaxy. We rotate the system to align the total angular momentum of the young stars (age $\leq 1$ Gyr) along the $Z$ direction, which orients the Galactic disk in the XY plane for each simulation. We establish the Solar Circle with a fixed cylindrical radius of $\textrm{R}_{\odot}=8.1\pm 0.6$ kpc (Gravity Collaboration et al., 2018). The rotation curves of these simulations are roughly flat at 8.1 kpc, with values close to $220\pm 20$ km s^-1, which roughly match that of the MW. Our choice is strictly motivated from the measured Solar position, in Sec. 3.2.1, we discuss our results for the DM density as a function of changing cylindrical radius from the center of each halo. While one could scale the cylindrical radius of each halo based on disk size or other parameters, this location would vary for each CDM and SIDM halo independently. Instead, we focus on what happens at a fixed radius away from the center in each halo.

We select a cylindrical region with a fixed Galactocentric cylindrical radius of $\textrm{R}_{\odot}=8.1\pm 0.6$ kpc and a height of $|Z|\leq 5$ kpc, centered around the Solar Circle. This cylindrical region is then divided into small 50 kpc bins in the X and Y direction, covering the $|Z|\leq 5$ kpc. Within each of these bins, we independently compute the surface density ($\Sigma_{X,Y}$) for baryons (stars and gas), DM, and for all (baryons and DM) the species combined together. Following this, we introduce a metric, $\delta_{\textrm{dens}}$, to quantify the relative surface density variation:

$$\delta_{\textrm{dens}}=\frac{\Sigma_{X,Y}-\overline{\Sigma_{R}}}{\overline{\Sigma_{R}}}$$ (1)

where $\overline{\Sigma_{R}}$ is the median surface density at $\mathrm{R={(X^{2}+Y^{2})}^{1/2}}$.

Fig. 1 illustrates relative surface density variation ($\delta_{\textrm{dens}}$) for m12f SIDM (top panel) and CDM (bottom panel) halos at the present day, categorized by species: all (left), baryons (middle), and DM (right). Major density variations are predominantly within the baryonic component, where overdense regions (red) can exhibit densities approximately twice as high ($\delta_{\textrm{dens}}\sim 1$). Conversely, the DM component showcases a uniform distribution, characterized by $|\delta_{\textrm{dens}}|<0.05$. This uniformity in DM serves to counterbalance the fluctuations observed in baryonic matter. However, when considering all species combined, notable relative density variations persist within the Solar Neighborhood. Similar patterns are observed in m12i and m12m simulations. Notably, no discernible systematic trend in densities is observed between CDM and SIDM models. Furthermore, majority of the high density regions exhibit elevated gas density and are actively forming new stars. We note that $\delta_{\textrm{dens}}$ is roughly independent on the choice of $Z$ range.

Next, we categorize the bins into three density regions – low, median, and high – based on their surface density in quartile ranges $\leq 25^{\textrm{th}}$ percentile, $25^{\textrm{th}}-75^{\textrm{th}}$ percentile, and $\geq 75^{\textrm{th}}$ percentile respectively. The three density regimes effectively provide upper and lower limits on the Galactic acceleration profiles.

Table LABEL:tab:gal_prop summarizes the galaxy-wide properties along with average variations in the low, median, and high density regions around the Solar Circle for all of the simulations. While, the median baryonic volume density in the Solar Circle across all halos is generally much lower ($5-8\sigma$) compared to the MW, the total surface density across the halos is relatively close to that of the MW (within $2-3\sigma$), except for m12i CDM that is about $6\sigma$ away. Also, McCluskey et al. (2024) motivated that the MW forms an unusually dynamically cold disk which could explain higher average density, contrasting with other observed galaxies of similar mass ($\sim 10^{12}$ M_$\odot$) through analysis of the Latte disks. Additionally, our estimates for the MW are solely based on the local volume around the Solar Neighborhood. In terms of matching local properties, the m12m SIDM halo is the most similar to the MW, exhibiting comparable thin and thick disk, and cold gas scale heights.

We observe that both the total volume and surface density between regions within a same halo can vary by a factor 2. For example, in the low, median, and high density regions, the volume density for m12m SIDM is approximately (36, 50, 79) $\cdot 10^{-3}\textrm{M}_{\odot}/\textrm{pc}^{-3}$. Only $\sim 6\%$ of our selected bins exhibit total volume densities within $1\sigma$ of the MW value of $99\pm 9.5\ \cdot 10^{-3}\textrm{M}_{\odot}/\textrm{pc}^{-3}$ and 10% bins have total volume density values higher than the MW. The majority of the variations arise from the baryonic component (stars and gas) and the DM component is evenly distributed azimuthally (lower errors in Table LABEL:tab:gal_prop, see Fig. 1.).

Comparing DM models, the SIDM halos exhibit both higher stellar and DM volume and surface densities in the Solar Neighborhood than their CDM counterparts. This is primarily due to SIDM particles’ capacity to exchange energy and momentum in addition to gravity, making them more responsive and sensitive to the baryonic potential (Sameie et al., 2018; Elbert et al., 2018; Despali et al., 2019; Robles et al., 2019; Santos-Santos et al., 2020). The higher DM density facilitates the entrapment of additional gas, forming new stars and establishing a feedback loop, where DM particles can respond faster to the rising stellar density (Despali et al., 2019; Sameie et al., 2021). While at about 8.1 kpc from the Galactic center, we anticipate only about one DM scattering event per Hubble time, these rates are higher in the inner regions of the galaxy (Vargya et al., 2022), which leads to more pronounced differences in DM density between CDM and SIDM in the inner regions (see Fig. 5 in Sec. 3.2.1). Consequently, SIDM halos exhibit higher density in DM and stellar distribution within the baryon-dominated potential, resulting in increased oblateness in the inner regions and around the Solar Circle.

In Sec. 3.2.1 (see Fig. 5), we show the azimuthally average DM density in different 2D cylindrical radii (between 4 to 12 kpc), showing that the trend of denser SIDM holds across a variety of radii. The differences between the CDM and SIDM are more pronounced in the inner regions, and considering that the MW has a thinner disk compared to the simulated galaxies, the expected effect of SIDM versus CDM as a function of vertical distance from the midplane ($Z$) in the MW should be even more significant.

2.2 Comparison with local potential models

Table 2: Best fit potential model and oblateness/shape parameters.

Galaxy	Physics	$\mathrm{M_{MN}}$ [$10^{10}\textrm{M}_{\odot}$]	a [kpc]	b [kpc]	$\mathrm{\log_{10}(-\gamma_{MN}/Gyr^{-2})}$	$\mathrm{\log_{10}(\alpha_{1}/Gyr^{-2})}$	$\mathrm{\log_{10}(-\gamma_{LO}/Gyr^{-2})}$	$\frac{c}{a}\Big{\rfloor}_{\textrm{R}_{\odot}}$
$\textrm{MW}^{\textrm{a}}$	?	10	6.5	0.26	3.94	$3.54\pm 0.16$	$3.3\pm 0.9$	–
m12i	CDM	$1.95\pm 0.04$	$3.91\pm 0.31$	$0.80\pm 0.12$	$2.83\pm 0.02$	$3.00\pm 0.00$	$2.87\pm 0.01$	0.46
	SIDM	$3.11\pm 0.07$	$3.17\pm 0.34$	$0.85\pm 0.23$	$3.01\pm 0.03$	$3.14\pm 0.01$	$2.99\pm 0.01$	0.43
m12f	CDM	$3.54\pm 0.13$	$3.30\pm 0.26$	$0.82\pm 0.33$	$3.08\pm 0.05$	$3.17\pm 0.02$	$3.15\pm 0.00$	0.40
	SIDM	$4.20\pm 0.05$	$3.50\pm 0.18$	$0.79\pm 0.05$	$3.17\pm 0.02$	$3.20\pm 0.00$	$3.16\pm 0.01$	0.38
m12m	CDM	$3.49\pm 0.24$	$6.31\pm 0.44$	$0.78\pm 0.04$	$2.97\pm 0.04$	$3.17\pm 0.01$	$2.92\pm 0.00$	0.30
	SIDM	$5.30\pm 0.09$	$6.81\pm 0.20$	$0.78\pm 0.04$	$3.11\pm 0.05$	$3.28\pm 0.01$	$2.95\pm 0.01$	0.27

Ideally, one would integrate the acceleration field (from e.g., pulsar timing or other acceleration measurements) to infer the potential directly. However, in practice, the scarcity of data points (e.g., 20 pulsars within 2 kpc of the Sun (Donlon et al., 2024)) currently limits our ability to perform such integration. Consequently, in the MW, observes use static parameterize models to fit the available data accurately. Here, we present our simulations with similar parameters through potential model fitting on stars in order for a direct comparison. We consider two axisymmetric potentials previously used in Chakrabarti et al. (2021). Specifically, a low-order (LO) expansion $\alpha\gamma$ potential (eq. 2, Chakrabarti et al., 2021) near the position of the sun (i.e adopted Solar Circle in our case), and a single Miyamoto-Nagai (MN) potential (eq. 3) for the Galactic disk.

The LO potential profile is defined as

$$\Phi(R,Z)=V_{\mathrm{LSR}}^{2}\log\left(R/R_{\odot}\right)+\log\left(R/R_{\odot}\right)\gamma_{\mathrm{LO}}Z^{2}+\frac{1}{2}\alpha_{1}Z^{2}$$ (2)

where $V_{\mathrm{LSR}}$ is the local standard of rest velocity, $\gamma_{\mathrm{LO}}$ determines the shape of the potential affecting the vertical accelerations, and $\alpha_{1}$ represents the squared angular frequency of oscillations in the vertical density profile. A higher value of $\log{-\gamma_{\mathrm{LO}}}$ value indicates a flatter potential profile in the Solar Neighborhood. F

For stars located within $\textrm{R}_{\odot}=8.1\pm 0.6$ kpc and at heights $|Z|\leq 1.1$ kpc, localized around the midplane, we fit the LO potential profile paramters $V_{\mathrm{LSR}}$, $\gamma_{\mathrm{LO}}$, and $\alpha_{1}$ using linear regressions on the potential values stored as a property for the particles in the simulation. To estimate the uncertainties on these parameters, we employ a Markov Chain Monte Carlo (MCMC) bootstrap approach, using 1000 samples. This involves repeatedly resampling the data with replacement and fitting the model to each resampled dataset.

To model the gravitational effect of the simulated galactic disk for stars within $3\textrm{kpc}\leq\textrm{R}\leq 15\textrm{kpc}$ and at heights $|Z|\leq 1.1$ kpc, we use a single MN profile for the limited vertical range around the midplane, given by

$$\Phi_{\mathrm{MN}}(R,Z)=-\frac{GM_{\mathrm{MN}}}{\sqrt{R^{2}+(a+\sqrt{Z^{2}+b^{2}})^{2}}}.$$ (3)

We fit the 2D density of the simulation with the density functional form of the MN profile using the scale mass ($M_{\mathrm{MN}}$), scale radius ($a$), and scale height ($b$) for the Galactic disk. The fitting is performed using a least squares optimization method, specifically the Nelder-Mead method (Lagarias et al., 1998), to find the best fit parameters and estimate the uncertainties on these parameters using MCMC bootstrap using 1000 samples. Arora et al. (2022) showed that such analytic potential models can reconstruct the rotation curve of the FIRE simulated galaxies to about $\pm 2\%$ accuracy in the disk region.

The oblateness parameter for the MN disk is computed by expanding eq. 3 to first order in $R$ and second order in $Z$ around ($R_{\odot},0$):

$$\gamma_{\mathrm{MN}}=-\frac{GM_{\mathrm{MN}}}{b}\frac{a+b}{\left(R_{\odot}^{2}+(a+b)^{2}\right)^{5/2}}\frac{3R_{\odot}^{2}}{2}.$$ (4)

Table 2 shows the best fit parameters and oblateness ($\gamma$) obtained for the simulations using the LO and MN potentials along with the best fit parameters for the MW from Donlon et al. (2024). It is noteworthy that the $\gamma$ values derived from the LO potential ($\gamma_{\mathrm{LO}}$) and the MN disk ($\gamma_{\mathrm{MN}}$) are similar in all simulations, except for the m12m SIDM case. This suggests that the shape of the local potential at the Solar Circle is primary influenced by the stellar disk.

The values obtained in the simulations are within $0.5\sigma$ of the MW best fit parameters from (Donlon et al., 2024). Our simulated SIDM halos match the observed MW values better than the CDM counterparts, particularly with regard to the shape and oblateness of DM halos. Specifically, the SIDM halos at the Solar Circle demonstrate lower minor to major axis ratios (by 0.02-0.03); in other words, they are more oblate than the CDM halos. This happens because SIDM magnifies the variations in SFR that enhances the concentration of both baryons and DM (see Table. LABEL:tab:gal_prop).

3 Vertical acceleration profiles

Figure 2: Vertical acceleration (top row) and change in the vertical velocity over a 10-year baseline (bottom row) as a function of vertical height $Z$ from the Solar Circle for CDM (blue) and SIDM (green) models in the three simulations (columns). The corresponding shaded regions show the upper and lower limits on the profiles, as computed from high and low density bins respectively. It is evident that both the vertical acceleration profile and the change in vertical velocity have consistently larger values in the SIDM simulations compared to their CDM counterparts.

We calculate the vertical accelerations above and below the Galactic midplane using a binned method ($\Delta Z=0.05$ kpc) to present smooth trends as a function of vertical distance ($Z$) from the midplane.

We slice the low, median, and high density regions as defined in Sec. 2.1 into 50 kpc vertical distance $Z$ bins. Next, we average the gravitational accelerations computed for stars and DM particles in each region within each $Z$ bin. With this approach, we obtain a smooth and accurate representation of the median acceleration profile; the high density bin provide an upper limit, and the low density bin provide a lower limit for each profile.

Fig. 2 plots the vertical acceleration $\mathrm{a_{Z}}$ (top row) and change in the vertical velocity over a 10-year baseline ($\mathrm{\Delta v_{Z}\equiv a_{Z}\Delta t}$) (bottom row) at present day for the m12i (left column), m12f (middle column), and m12m (right column) in CDM (blue) and SIDM (green) model. The shaded regions represent the upper and lower limits on the profiles derived from the high and low density regions respectively. SIDM halos have consistently larger magnitude of their accelerations than CDM, thanks to their higher surface densities, especially at larger distances from the mid-plane. We note that the degree of difference in $\mathrm{\Delta v_{Z}}$ between SIDM and CDM profiles varies with the simulation, with the smallest difference observed for m12f (one with a massive merger about 3 Gyr ago before the present day) and the largest difference for m12m (no mergers in the last 10 Gyr). This observation suggests that mergers may have a role in removing or mitigating differences between SIDM and CDM profiles. We discuss this more in Sec. 3.2.

3.1 Median vertical acceleration gradient

Figure 3: Vertical acceleration gradient ($\mathrm{\frac{da_{Z}}{dZ}(Z)}$) (left) and the asymmetric difference between $+Z$ and $-Z$ (right) with respect to vertical height $Z$ for CDM (solid) and SIDM (dotted) halos. SIDM halos have consistently steeper profiles compared to the CDM halos. Notable asymmetries are observed in m12f and m12i for $Z\geq 1$ kpc (dashed line in second panel). The MW’s measured $\mathrm{\frac{da_{Z}}{dZ}}$ at midplane is approximately -3200 $\textrm{Gyr}^{-2}$. Gray shading represents the current measurement uncertainty for pulsar studies in $\mathrm{\frac{da_{Z}}{dZ}}$ (Chakrabarti et al., 2021; Donlon et al., 2024), and the green region indicates the expected precision increase, scaling pulsar timing precision by a power of 5/3 over a 20 year baseline. Both the shaded regions are scaled to match the minimum value in m12m SIDM. Fig. 5 in Sec. 3.2.1 plots the $\mathrm{da_{Z}/dZ(Z=0)}$ at midplane ($Z=0$) as a function of 2D cylindrical radius. SIDM halos show consistently larger value of $\mathrm{\frac{da_{Z}}{dZ}(Z=0)}$ across a range of radii. The differences between CDM and SIDM are more pronounced closer to the galactic center.

Fig. 3 depicts the vertical acceleration gradient ($\mathrm{\frac{da_{Z}}{dZ}}$) (left) and the asymmetric difference around the midplane (right) as a function of vertical height $Z$ for simulations (color-coded) in CDM (solid line) and SIDM (dotted line) models. The gradients are computed using total-variation regularization differentiation algorithm which avoids the noise amplification in finite-difference methods for noisy data (Chartrand, 2011). The asymmetric difference is given by

$$\Bigg{|}\Delta\mathrm{\frac{da_{Z}}{dZ}}\Bigg{|}\equiv\Bigg{|}\mathrm{\frac{da_{Z}}{dZ}\Big{\rfloor}_{+Z}-\frac{da_{Z}}{dZ}\Big{\rfloor}_{-Z}}\Bigg{|}$$ (5)

The gray shaded region indicates the current measurement uncertainty in $\mathrm{\frac{da_{Z}}{dZ}}$ from Donlon et al. (2024), while the green shaded region represents the anticipated measurement uncertainty with a 20-year baseline (Lorimer & Kramer, 2005). Note, the majority of the azimuthal variation stems from the observable baryonic component and the azimuthal distribution of DM is effectively constant in these simulations. In fact, variations arising from the baryonic component can be mitigated with precise measurements, as such we only show the median profiles for $\mathrm{\frac{da_{Z}}{dZ}}$ in Fig. 3,  4, and  5.

The $\mathrm{\frac{da_{Z}}{dZ}}$ at the midplane ($Z=0$ kpc) in our simulations closely aligns with the MW’s best-fit value of $-3200\pm 2560$ Gyr^-2 (Donlon et al., 2024), albeit within $1\sigma$ error of the MW measurement. Discrepancies primarily arise from lower median total matter density in the Solar Neighborhood in the simulations relative to the measured MW values that is locally measured. Chakrabarti et al. (2021) reported a value of $-4900^{+1600}_{-2700}$ Gyr^-2 using data from 14 binary pulsars distributed over $\sim$ 1 kpc from the Sun, which is comparable to the more updated value given in Donlon et al. (2024). This value corresponds to the square of the frequency of low-amplitude vertical oscillations (denoted $\alpha_{1}$ in Chakrabarti et al. (2021)), and is used to determine the mid-plane density, or the Oort limit. This gives a value of the Oort limit of $0.08^{0.05}_{-0.02}~{}M_{\odot}/\rm pc^{3}$, which is close to but lower than kinematic estimates. Due to the limitation of data available at that time, the analysis in Chakrabarti et al. (2021) could not constrain global properties, like the mass. Donlon et al. (2024) find a value for the Oort limit of $0.062\pm 0.017~{}M_{\odot}/\rm pc^{3}$, which is lower still relative to kinematic estimates. Using simplified density profiles, Donlon et al. (2024) estimates the galaxy mass within the Solar Circle to be nearly twice as high as the currently accepted value from Bland-Hawthorn & Gerhard (2016). MW measurements may however be impacted by high uncertainties currently in binary pulsar data. These measurements will continue to improve as more pulsar timing data become available.

Notably, the depth of the $\mathrm{\frac{da_{Z}}{dZ}(Z)}$ as a function of vertical distance from the midplane varies significantly across simulations, with SIDM simulations consistently displaying deeper valleys by 10-30% within 1 kpc of the midplane. The most substantial gradient and the most pronounced difference between SIDM and CDM profiles are observed in m12m because it is the most massive disk (interestingly, it also has the lowest Toomre Q value). In contrast, m12i CDM displays the shallowest valley close to the midplane, aligning with the volume densities of their respective systems. Interestingly, m12f, which experienced a recent merger, shows the least pronounced difference between CDM and SIDM profiles, while m12m, which had the earliest merger, exhibits the most significant difference between the DM profiles close to the midplane. In the outskirts ($|Z|\geq 2$ kpc), the difference between profiles is smaller. This observation suggests that the presence of recent major mergers such as Sag dSph and the LMC in a galaxy might mitigate the distinctive signatures of DM models on the vertical acceleration gradient profile, particularly close to the mid-plane. This could be due to multiple factors, such as an epoch of increased star formation triggered by the influx of gas from the merging satellites (Di Matteo et al., 2007; Hopkins et al., 2010; Pearson et al., 2019), and/or heating of the galactic center due to the energy exchange with the merging satellite (e.g., Barnes & Hernquist, 1992).

Given the uncertainties regarding the exact location of the Solar Circle within our simulation, Sec. 3.2.1, Fig.5 shows $\mathrm{\frac{da_{Z}}{dZ}}$ as a function of 2D cylindrical radius spanning a range of 4–12 kpc. SIDM halos have consistently higher values of $\mathrm{\frac{da_{Z}}{dZ}(Z=0)}$ across all radii. The difference between SIDM and CDM model is most pronounced near the galactic center. Notably, while m12m exhibits the most significant difference, m12f displays relatively minor disparities, aligning with the trends observed in Fig.3. These highlight the consistency of our results across a broad range of radii.

We also note a significant asymmetry around the mid-plane in $\mathrm{\frac{da_{Z}}{dZ}}$ (right panels in Fig. 3), particularly pronounced for m12f, which underwent recent massive mergers, and also for m12i. In contrast, m12m, which had no recent mergers in the CDM case, exhibits a less pronounced asymmetry. For SIDM simulations, a similar trend is observed between m12f and m12i, with noticeable asymmetries in m12m as well.

These asymmetries are influenced by various factors, including the orbits of merging satellites, which can excite long-standing nodes in the Galactic disk (Weinberg, 1998), leading to higher variation in azimuthal direction. As argued by Chakrabarti et al. (2019, 2020), recent mergers leave a discernible imprint on the vertical acceleration profiles persisting long after the complete tidal stripping of satellites. These signatures become more pronounced at greater distances from the mid-plane, highlighting the enduring impact of recent merger events. We observe large vertical asymmetries especially in the gas disk in the outer part of the MW (Levine et al., 2006) that may arise from an interaction with a massive dark matter sub-halo (Chakrabarti & Blitz, 2009). The global asymmetries in the baryonic component, may allow us, together with the observed asymmetry in the total acceleration field, to more comprehensively model past interactions with dwarf galaxies.

3.2 Temporal and spatial variation in the vertical acceleration gradient.

Figure 4: Left: Azimuthally averaged vertical acceleration gradient ($\mathrm{\frac{da_{Z}}{dZ}}$) as a function of vertical height $Z$ for the m12f CDM halo at different epochs: at present day (magenta color) for CDM (solid) and SIDM (dotted), at the time of the satellite’s first pericentric passage (3 Gyr before the present day, indigo color), and 4.2 Gyr before the present day (orange color). Before the merger event, the profile is symmetric across the mid-plane ($Z=0$), but it becomes highly asymmetric during the merger and retains some asymmetry until the present day. Right: The vertical acceleration gradient ($\mathrm{\frac{da_{Z}}{dZ}}$) computed in different quadrants (marked in insets in the bottom left corner) for the m12f CDM halo at the time of the first pericentric passage. Q. I and II, which contain two satellites close to peri, exhibit the most asymmetry across the midplane. In contrast, Q. III and IV, which have no satellites, display relatively symmetric profiles. The inset in bottom left plots the disk plane for the m12f CDM at first pericenter and marks the relevant quadrants.

In previous sections, we observed a consistent trend: azimuthally averaged vertical acceleration gradient profiles in SIDM simulations were systematically steeper, typically by 10-30%, compared to their CDM counterparts. However, it remains whether temporal and azimuthal variations expected in $\mathrm{\frac{da_{Z}}{dZ}(Z)}$ are similar in level and scale as the DM model dependent differences because if this is true, it could imply that the variations we see in the simulations might not be due to actual differences in DM properties, but are actually transient effects due to mergers. In this section, we examine the spatial and temporal variations in the $\mathrm{\frac{da_{Z}}{dZ}(Z)}$ using the m12f CDM simulation as an example. This simulation experienced a major merger with two satellites with their first pericentric passage approximately 3 Gyr ago. With the total mass ratio of 15 and 18, when each satellite was at pericentric distance of about 25 kpc.

The left panel of Fig. 4 illustrates the azimuthally averaged $\mathrm{\frac{da_{Z}}{dZ}(Z)}$ profiles at various epochs: present day (magenta) both for CDM (solid) and SIDM (dotted), during the first pericentric passage (3 Gyr before the present day, indigo), and 4.2 Gyr before the present day (orange). These profiles display significant variations (10-30%) in depth, influenced by the average total matter density. Approximately 1 Gyr before the pericentric passage, the profile is steeper due to enhanced star formation in the region around the Solar Circle and symmetric around the mid-plane. At the first pericentric passage, the profile is shallower and highly asymmetric. By the present day, the increasing baryon density has substantially steepened the acceleration profile, with the asymmetry persisting to some extent. These temporal variations in profile depth are comparable to the highest order variations seen in m12m CDM and SIDM profiles in Fig. 3 and higher than the variations observed in the the m12f CDM and SIDM profiles.

The right panel of Fig. 4 depicts $\mathrm{\frac{da_{Z}}{dZ}(Z)}$ computed in four azimuthal quadrants (marked in insets on the bottom left) during the first pericentric passage in m12f CDM. All quadrants exhibit consistent profile depths within a 10% range, except for Q. III, which shows increased depth due to high-density gas prompting new star formation. The variation at large $|Z|$ stem from limited particle data away from the midplane leading to noisy gradients. Quadrants containing the two satellites (Q. I and II) display more significant asymmetry across the midplane, while the quadrants without any satellites (Q. III and IV) have relatively symmetric profiles, highlighting the impact of satellite orbits on gradient asymmetry.

3.2.1 Radial variation in the local DM density

Figure 5: Left: Azimuthally averaged DM density ($\rho_{\textrm{DM}}$) within a vertical range of $|Z|\leq 0.2$ kpc as a function of 2D cylindrical radius ($R_{\textrm{cylindrical}}$) for: m12f (magenta), m12i (cyan), and m12m (red) for CDM (solid circles), and SIDM (plus) at present day. In general, SIDM halos have a higher DM density closer to the galactic center. Notably, m12m SIDM maintains consistently higher density compared to the CDM, whereas m12i halo, with the lowest stellar mass (see Table LABEL:tab:gal_prop), shows a transition where SIDM becomes less dense (puffier) compared to CDM close to 8 kpc. Right: Vertical acceleration gradient ($\mathrm{\frac{da_{Z}}{dZ}(Z=0)}$) at the midplane ($Z=0$ kpc) as a function of cylindrical radius ($R_{\textrm{cylindrical}}$) for all the halos as in the left panel. The trends in $\mathrm{\frac{da_{Z}}{dZ}(Z=0)}$ are similar to those observed in DM density, with SIDM consistently exhibiting steeper gradients compared to CDM. The differences between SIDM and CDM are more pronounced closer to the center. Additionally, the current MW DM density (Bland-Hawthorn & Gerhard, 2016) and $\mathrm{\frac{da_{Z}}{dZ}(Z=0)}$ (Donlon et al., 2024) are shown with gray crosses with error bars.

Fig.5 plots the azimuthally averaged DM density ($\rho_{\textrm{DM}}$) calculated within a vertical range of $|Z|\leq 0.2$ kpc (left) and vertical acceleration gradient ($\mathrm{\frac{da_{Z}}{dZ}(Z=0)}$) at midplane ($Z=0$ kpc) as a function of 2D cylindrical radius ($R_{\textrm{cylindrical}}$) for the three simulations (color-coded) in both CDM (solid circles) and SIDM (plus) models. In the inner regions of the galaxies, SIDM models, influenced by the baryonic component, exhibit systematically higher DM density and steeper gradients at the mid-plane compared to their CDM counterparts for all radii. These differences decrease further away from the center. These pronounced differences are consistent with the higher scattering rates of DM particles close to the center (Vargya et al., 2022).

In m12m, SIDM maintains more pronounced differences between densities and $\mathrm{\frac{da_{Z}}{dZ}(Z=0)}$ even at larger $R_{\textrm{cylindrical}}$, with no sign of convergence between CDM and SIDM up to 12 kpc. Conversely, in m12i, beyond 8 kpc, the SIDM profile becomes less dense than the CDM profile. This “puffier” profile is due to the absence of a strong baryonic potential in the outskirts, allowing DM particles to self-interact and form less dense, rounder profiles. Additionally, there are minor differences in the $\mathrm{\frac{da_{Z}}{dZ}(Z=0)}$ at larger $R_{\textrm{cylindrical}}$.Furthermore, m12f exhibits the least variation between CDM and SIDM at all radii in both density and $\mathrm{\frac{da_{Z}}{dZ}(Z=0)}$. However, $\mathrm{\frac{da_{Z}}{dZ}(Z=0)}$ remains systematically steeper in the SIDM model.

The galatic disks in these simulated halos are thicker compared to the MW’s disk (Wetzel et al., 2023; McCluskey et al., 2024). This implies that the effects of SIDM versus CDM could be even more pronounced in a thinner disk like that of the MW.

4 Discussion and conclusions

Direct measurement of Galactic accelerations presents a promising avenue for testing DM, bypassing the need to solve the “inverse problem” (Binney & Tremaine, 2011) inherent in traditional methods of inferring the mass distribution within galaxies from observed motions. By comparing observed acceleration profiles with predictions from simulations with different DM models, we can assess the compatibility of these models with observational data. However, testing SIDM in particular with this strategy offers some challenges due to the ability of SIDM to respond efficiently to changes in the baryon distribution. Galaxy formation can alter even the shape and density profile of a CDM halo; in SIDM, which can more easily alter its energy and angular momentum through interactions, the expectation is that the halo’s shape and profile are more tightly correlated with the disk properties (Sameie et al., 2021; Vargya et al., 2022). However, this response from the SIDM component in the galaxy subsequently affects the baryonic disk, since the deepened DM potential increases the gas density and boosts star formation. Thus the long-term co-evolution of the SIDM and baryons leads to a slightly higher stellar mass, a slightly more dense DM halo in the disk plane, and slightly higher star formation rates than in CDM. These differences are evident in our Table LABEL:tab:gal_prop and in Fig. 5, and are the means by which SIDM solves the so-called “diversity problem.”

At the Solar Circle, the distribution and properties of baryonic matter play a role equal to or larger than DM’s in creating the observed Galactic acceleration profiles (see Fig. 4). The steeper vertical acceleration gradients in SIDM relative to the same halo in CDM (Fig. 3), which are calculated at the Solar Circle, thus reflect the response of both species to the change in the DM model. However, the variation in the acceleration profile among CDM galaxies, which is simply due to their different star formation and assembly histories, spans nearly as large a range. Thus, observations of the density distribution in the MW alone are unlikely to distinguish CDM from SIDM. However, a large sample of measurements of the disk-plane density of galaxies could potentially show a signal, since in SIDM one would expect the mean density in such a sample to be statistically higher than expected from CDM. The trend of density with radius in such a sample may also help distinguish the two theories, since the differences in density between CDM and SIDM become more pronounced at smaller radius (where one expects more frequent self-interactions). At the Solar Circle (8.1 kpc), with a cross-section of 1 cm² g^-1, we expect about 1 scattering event per Hubble time ($\sim 10$ Gyr) per DM particle in our simulations, with much higher scattering rates in the inner regions of the galaxy (Vargya et al., 2022).

Measuring Galactic accelerations may help to illuminate the merger history of the MW, since mergers induce asymmetries in the acceleration profile (Chakrabarti et al., 2021) (see Fig. 3, 4 in this paper). This highlights the need for flexible disequilibrium models to fully utilize acceleration data (Donlon et al., 2024), given that merger-induced variations often exceed the differences induced by changing the DM model. Additionally, we only explored the Galactic acceleration profiles in simulations run with standard Monte Carlo implementation for SIDM (Rocha et al., 2013; Peter et al., 2013) using the FIRE-2 prescription for baryons (Hopkins et al., 2018). In the future, one should explore different simulation suites performed with varied baryon prescriptions (e.g. Teyssier, 2002; Menon et al., 2015; Wadsley et al., 2017; Weinberger et al., 2020) and alternative SIDM implementations (e.g. Vogelsberger et al., 2012; Fry et al., 2015; Meskhidze et al., 2022).

Our main conclusions are summarized below:

1. The Solar Neighborhood in the MW is comparatively denser than the median of azimuthally averaged Solar Circle regions in these three FIRE-2 simulations with CDM and SIDM model. (see Sec. 2.1 and Table LABEL:tab:gal_prop). Regions with local density similar to the Solar Neighborhood are relatively limited and make up about 6% of the volume at Solar Circle. The higher density in the MW can be attributed to the relatively thin MW disk even compared to other $\sim 10^{12}$ M_$\odot$ galaxies (McCluskey et al., 2024) and/or our placement in a dense region of the Galaxy. It should be noted that we did not explore all of the FIRE-2 MW-mass simulations, so a systematic comparison to better quantify this in relation to the MW remains to be conducted.

2. The shape of the potential and density in the Solar Neighborhood are predominantly influenced by the Galactic disk (see Sec. 2.2). However, SIDM particles greater responsiveness to the disk potential leads to higher flattening in the DM shape at Solar Circle. This leads to measurable distinctions in the local dark matter density within the Solar Neighborhood. SIDM halos consistently demonstrate denser and more oblate potentials, resulting in vertical acceleration gradient profiles that are steeper by 10-30% compared to CDM. (see Sec. 3.1 and Fig. 3). However, the galaxy-to-galaxy variation in density is broader than the difference between SIDM and CDM at 1 cm²/g.

3. Recent mergers with total mass equal to or greater than the Sagittarius dwarf or the SMC will induce measurable asymmetries across the midplane in the vertical acceleration gradient profiles (Chakrabarti et al., 2019, 2020), which can persist over extended periods to the present day (see Sec. 3.2 and Fig. 4). These asymmetries—influenced by factors such as satellite mass and orbit—can be larger than the systematic difference between SIDM and CDM, which poses challenges for probing the nature of DM in a single galaxy. Nonetheless, these asymmetric profiles are a promising probe of the merging history of a galaxy.