# Flattening of Dark Matter Cusps during Mergers: Model of M31

Pierre Boldrini, Roya Mohayaee, Joseph Silk

<figure><img src="image_1.png" /></figure>

## To cite this version:

Pierre Boldrini, Roya Mohayaee, Joseph Silk. Flattening of Dark Matter Cusps during Mergers: Model of M31. The Astrophysical Journal, 2021, 919 (2), pp.86. 10.3847/1538-4357/ac12d3. hal-02507751

HAL Id: hal-02507751

https://hal.science/hal-02507751v1

Submitted on 3 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Flattenning of Dark Matter Cusps during Mergers: Model of M31

Pierre Boldrini¹, Roya Mohayaee¹, and Joe Silk¹,²,³

¹ Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, F-75014 Paris, France

² Department of Physics and Astronomy, The Johns Hopkins University, Baltimore MD 21218, USA

³ Beecroft Institute for Particle Astrophysics and Cosmology, Department of Physics, University of Oxford, Oxford OX1 3RH, UK

Received 2020 November 13; revised 2021 June 18; accepted 2021 July 7; published 2021 September 28

## Abstract

We run high-resolution fully GPU N-body simulations to demonstrate that the dark matter (DM) distribution in M31 is well fitted by a core-like profile. Rich observational data especially on the giant stellar stream provides stringent constraints on the initial conditions of our simulations. We demonstrate that the accretion of a satellite on a highly eccentric orbit heats up the central parts of M31, causes an outward migration of DM particles, flattens the central cusp over more than a decade in scale, and generates a new model-independent DM profile that is well fitted by a core. Our results imply that core-like central profiles could be a common feature of massive galaxies that have been initially cuspy but have accreted satellites on nearly radial orbits.

*Unified Astronomy Thesaurus concepts: Dark matter density (354); Galaxy dynamics (591); Galaxy mergers (608); Andromeda Galaxy (39)*

## 1. Introduction

Due to its proximity, the Andromeda galaxy (M31) provides a wealth of high-precision observational data for understanding the history of M31 and of the Local Group, and more generally, galaxy formation models in a cold-dark-matter-dominated universe. M31 exhibits challenging features on different scales, ranging from the double nucleus at its center on a scale of a few parsecs (Tremaine 1995; Kormendy & Bender 1999) to the giant stellar stream (GSS) in its outskirts, which extends to tens of kiloparsecs (Ibata et al. 2001).

It is widely believed that the phase structure of M31, namely, its GSS and its shell-like features, result from the accretion of a satellite galaxy (Ibata et al. 2004; Fardal et al. 2006, 2007; Font et al. 2006). The mass, radial velocity, and distance to the stream are observed with good accuracy. The accretion scenario has been simulated with high resolution, and the high-quality data on the stream and the shell-like features have been used to strongly limit the initial parameter space. In previous work, we ran high-resolution simulations of a live M31 and an infalling satellite and showed that in the cosmological scenario, the infalling satellite traces a highly eccentric orbit after reaching its turnaround radius, and then falls toward the center of M31 (Sadoun et al. 2014). The satellite was disrupted by M31, forming the GSS, and its subsequent passages through the center of M31 led to the formation of the shell-like features that we observe today. Our work showed that the infalling satellite was dark matter (DM) rich and in the same plane as most of the present-day satellites of M31. Also, minor merger scenarios, as discussed above, have been successful in explaining the phase structures of M31. Major merger scenarios have also been proposed in the past (Hammer et al. 2010), reviewed by Hammer et al. (2018), and received recent attention, seeming to fit well the kinematic data (e.g., D'Souza & Bell 2018).

The density profiles of galaxy halos provide important constraints on the nature of DM. In a cold-DM-dominated universe, the halos have profiles that diverge at the centers, i.e., cuspy profiles (Fukushige & Makino 1997; Navarro et al. 1997; Moore et al. 1998; Navarro et al. 2010). In contrast, measurements of galaxy rotation curves and dynamical models

of dSph galaxies are often claimed to require shallower slopes that are consistent with a central density core (e.g., Flores & Primack 1994; Moore 1994; Battaglia et al. 2008; Walker & Peñarrubia 2011; Agnello & Evans 2012; Amorisco & Evans 2012; Adams et al. 2014, and references therein). However, it has been pointed out that disequilibrium, inclination, and noncircular motions can mimic the presence of cores in galaxies (Hayashi & Navarro 2006; Read et al. 2016; Oman et al. 2019). The core profiles in dwarf galaxies (Oman et al. 2015; Santos-Santos et al. 2020) generally present a challenge to the present model of cosmology (ΛCDM). This discrepancy has become known as the core-cusp problem. In order to reconcile observations and simulations, many mechanisms involving baryons have been proposed, transforming cusps into cores via changes in the gravitational potential caused by stellar feedback redistributing gas clouds and generating bulk motions and galactic winds (e.g., Navarro et al. 1996; El-Zant et al. 2001; Goerdt et al. 2010; Ogiya & Mori 2011; Pontzen & Governato 2012; Teyssier et al. 2013; Benítez-Llambay et al. 2019). The existence of shallower than Navarro-Frenk-White (NFW) density profiles in low-mass galaxies has sometimes been interpreted as the manifestation of a new DM-specific feature, such as DM self-interactions. The multiple scattering events between DM particles can result in the formation of constant-density cores by removing particles from the centers of halos (e.g., Vogelsberger et al. 2012; Rocha et al. 2013; Zavala et al. 2013; Robertson et al. 2018). Alternatively, in a recent study, we have shown that DM candidates in the form of primordial black holes can induce a cusp-to-core transition in low-mass dwarf galaxies via dynamical friction by DM particles (Boldrini et al. 2020). Yet another mechanism could be due to sinking massive objects, such as subhalos, gas clumps, or globular clusters, that transfer energy and angular momentum to the DM field via dynamical friction, creating a DM core from an initially cuspy density distribution. The DM density profile of the galaxy flattens due to heating via dynamical friction (El-Zant et al. 2001; Ma & Boylan-Kolchin 2004; Tonini et al. 2006; Romano-Díaz et al. 2008; Jardel & Sellwood 2009; Goerdt et al. 2010; Inoue & Saitoh 2011).
In this work, we describe another consequence of accretion of a satellite on a highly elongated orbit in a $\Lambda$CDM universe and show that the density profile of M31 could have been strongly influenced by this mechanism. Such an orbit is motivated, and indeed required, by dynamical modeling of features in M31. We demonstrate that the passage of the satellite, which is at the origin of the GSS, the shells, and the warp of the disk of M31, should have also caused a cusp-to-core transition at the center of the DM distribution in M31. A more general consequence of our work is that accreting satellites on highly eccentric orbits can induce a cusp-to-core transition in CDM halos. These cores are a common feature of many DM halos that have been initially cuspy but have accreted subhalos on highly eccentric orbits. In most cosmological simulations, we expect to see this effect only for DM halos with masses higher than $10^{12} M_{\odot}$ by assuming that they have sufficient resolution to determine the shapes of DM density profiles (Ma & Boylan-Kolchin 2004; Gao & White 2007; Laporte & White 2015). The DM density profile of the Aquarius simulation revealed that density profiles become shallower inwards down to the innermost resolved radius (Navarro et al. 2010). This slight deviation from the NFW model could be evidence for the proposed cusp-to-core mechanism (El-Zant et al. 2001; Ma & Boylan-Kolchin 2004). Indeed, most halos have suffered multiple subhalo mergers, especially at early epochs ($z \sim 2-3$) (Neistein & Dekel 2008). However, the cusp can regenerate itself. As such, cuspy profiles are more common at recent epochs, possibly explaining the presence of transient cores (see Dekel et al. 2003; Laporte & Penarrubia 2015; Boldrini et al. 2020).

In this work, we perform high-resolution N-body simulations with GPUs, which allow parsec resolution, to study the effect of the accretion of a satellite on the central density profile of the DM halo of M31. The initial conditions of our simulations are determined by observations of the mass, density profile, radial velocities, and distances of the GSS, which provide high-precision tests of our model. We consider infalling galaxy scenarios with different halo-to-stellar mass ratios ranging from 0 (Fardal et al. 2007) to 20 (Sadoun et al. 2014) and to 100 in three new models. By analyzing the density profile of the halo of M31, we see that the initial cuspy profile becomes shallower as the DM particles are heated, during the passages of the infalling galaxy, and some DM particles migrate outwards from the central region of M31 in all our five models. As the flattening of the cusp occurs for all of these models, we propose that the central DM core-like profile of M31 is model independent.

The paper is organized as follows. In Section 2 we present a brief summary of observational data. Section 3 provides a description of the N-body modeling of M31 and its satellite, along with details of our numerical simulations. In Section 4, we present the results from simulations and discuss the implications for the cusp-core problem. Section 5 presents our conclusions.

## 2. Observation: Giant Stellar Stream and Shell-like Features of M31

Phase structures, such as streams and shells, are ubiquitous in the universe where galaxies are interacting through gravity and frequently undergo mergers (Malin & Carter 1980, 1983; Hernquist & Quinn 1988, 1989). Our nearest neighbor M31 provides an excellent laboratory for the study of such tidal

features. M31 exhibits a GSS and shells and rings (Ibata et al. 2001; Ferguson et al. 2002; Bellazzini et al. 2003; Zucker et al. 2004; Ibata et al. 2005; Block et al. 2006; Brown et al. 2006; Ibata et al. 2007; Richardson et al. 2008; McConnachie et al. 2009). The main focus of this paper is the GSS, which provides stringent constraints on our models. The GSS is a faint stellar tail with a mass of $2.4 \times 10^8 M_{\odot}$ (corresponding to a luminosity of $3.4 \times 10^7 L_{\odot}$ and a mass to light ratio of 7) and extends to a projected radius of about 68 kpc on the sky (Ibata et al. 2001; Fardal et al. 2006) and is further accompanied by two stellar shells (Ferguson et al. 2002, 2005; Fardal et al. 2007; Richardson et al. 2008; Tanaka et al. 2010; Fardal et al. 2012). The similarity in the stars of the GSS and the shells has been the main reason why many models have taken these features to be of similar origin (Ibata et al. 2004; Font et al. 2006; Fardal et al. 2007).

An empirical minor merger scenario in which a satellite galaxy with no DM falls from a very close distance to M31 very recently (less than 1 Gyr ago) has been extensively used in the literature to model the GSS and the shells (Fardal et al. 2006, 2007). This model however suffered from simplifications. First, the dynamical friction effect was ignored as M31 was not treated as a live galaxy; second, the satellite had no DM in spite of most dwarfs being DM rich; and third, the satellite started its infall from a very close-in radius. The latter is puzzling as the satellite would have been disrupted on such a highly radial orbit (see Sadoun et al. 2014 for full details).

In a recent work, we proposed a new model for the origin of the GSS, the shells, and also the warped structure of the M31 disk itself (Sadoun et al. 2014). In our model, a DM-rich satellite is accreted and falls from its first turnaround radius on an eccentric orbit onto M31. The best agreement with the observational data is obtained when the satellite falls in about 2 Gyr ago and remains on the same plane that presently contains many of the dwarfs of M31 (Conn et al. 2013; Ibata et al. 2013). Unlike the previous model, in which the satellite made no impact on the disk as it was DM poor, the disk of M31 is perturbed by the infall of the massive satellite in our model and becomes warped.

In this paper, we use this cosmologically motivated scenario to set up our simulations. To show that our results are universal and hold for a wide range of initial conditions, we also run simulations for the model proposed by Fardal et al. (2007) in which the satellite is DM poor and starts its infall from a much closer distance to the center of M31. In addition, we also ran simulations for a higher halo-to-stellar mass ratio $M_{DM}/M_{*}$ of 100, which seems to be favored by some models of galaxy formation and evolution (e.g., Guo et al. 2011; Sawala et al. 2014; Henriques et al. 2015). For the latter, we run three simulations for satellites with different DM density profiles. Our simulations are summarized in Table 1.

A major improvement to our previous work (Sadoun et al. 2014) is achieved here because we have used a fully GPU-scaled code and gain in mass resolution by a factor of 100. This allows us to study the impact of the infalling satellite not just on the outer parts of M31 but also on the DM distribution at its center. The rich observational data on the GSS and shells of M31 provide rather demanding tests of our model. In the following section, we discuss the details of our numerical simulations.
<table><caption>**Table 1**
Simulation Parameters</caption><thead><tr><th colspan="7">The Host Galaxy</th></tr><tr><th></th><th>Component</th><th>Profile</th><th>a<br>(kpc)</th><th>r<sub>200</sub><br>(kpc)</th><th>Mass<br>(10<sup>10</sup> M<sub>☉</sub>)</th><th>(x<sub>0</sub>,y<sub>0</sub>,z<sub>0</sub>)<br>(kpc)</th><th>(v<sub>x0</sub>, v<sub>y0</sub>, v<sub>z0</sub>)<br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="6">M31</td><td>DM halo</td><td>NFW</td><td>7.63</td><td>195</td><td>88</td><td>0</td><td>0</td></tr><tr><td>Bulge</td><td>Hernquist</td><td>0.61</td><td>...</td><td>3.24</td><td>0</td><td>0</td></tr><tr><td rowspan="2">Disk</td><td>Exponential</td><td>R<sub>d</sub> = 5.4</td><td>...</td><td>3.66</td><td>0</td><td>0</td></tr><tr><td>disk</td><td>z<sub>d</sub> = 0.6</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td colspan="7">The infalling<br>satellite</td></tr><tr><td colspan="7">Scenario</td></tr><tr><td>Fardal et al. (2007)<br>(M<sub>DM</sub>/M<sub>*</sub>)<sub>sat</sub> = 0</td><td>Stars</td><td>Plummer</td><td>1.03</td><td>...</td><td>0.22</td><td>(-34.75,19.37,-13.99)</td><td>(67.34,-26.12,13.5)</td></tr><tr><td rowspan="2">Sadoun et al. (2014)<br>(M<sub>DM</sub>/M<sub>*</sub>)<sub>sat</sub> = 20</td><td>DM halo</td><td>Hernquist</td><td>12.5</td><td>20</td><td>4.18</td><td>(-84.41,152.47,-97.08)</td><td>0</td></tr><tr><td>Stars</td><td>Plummer</td><td>1.03</td><td>...</td><td>0.22</td><td>(-84.41,152.47,-97.08)</td><td>0</td></tr><tr><td rowspan="2">Satellite A<br>(M<sub>DM</sub>/M<sub>*</sub>)<sub>sat</sub> = 100</td><td>DM halo</td><td>Hernquist</td><td>25</td><td>41</td><td>22</td><td>(-84.41,152.47,-97.08)</td><td>0</td></tr><tr><td>Stars</td><td>Plummer</td><td>1.03</td><td>...</td><td>0.22</td><td>(-84.41,152.47,-97.08)</td><td>0</td></tr><tr><td rowspan="2">Satellite B<br>(M<sub>DM</sub>/M<sub>*</sub>)<sub>sat</sub> = 100</td><td>DM halo</td><td>Hernquist</td><td>25</td><td>123</td><td>22</td><td>(-84.41,152.47,-97.08)</td><td>0</td></tr><tr><td>Stars</td><td>Plummer</td><td>1.03</td><td>...</td><td>0.22</td><td>(-84.41,152.47,-97.08)</td><td>0</td></tr><tr><td rowspan="2">Satellite C<br>(M<sub>DM</sub>/M<sub>*</sub>)<sub>sat</sub> = 100</td><td>DM halo</td><td>Hernquist</td><td>25</td><td>41</td><td>22</td><td>(-168.82,304.94,-194.16)</td><td>0</td></tr><tr><td>Stars</td><td>Plummer</td><td>1.03</td><td>...</td><td>0.22</td><td>(-168.82,304.94,-194.16)</td><td>0</td></tr></tbody></table>

Note. From left to right, the columns give for each component: the initial density profile; the scale length; the virial radius; the mass; the initial positions in a reference frame centered on M31 with the x-axis pointing east, the y-axis pointing north, and the z-axis corresponding to the line-of-sight direction; and the velocities in this reference frame. We set the particle resolution of all the live objects to $4.4 \times 10^4 M_{\odot}$ and the gravitational softening length to 10 pc for all components.

### 3. Simulation: High-resolution Fully GPU code

The initial conditions for the M31 satellite are taken from Sadoun et al. (2014) and Fardal et al. (2007), and furthermore, we simulate three further models for a DM-rich satellite with $M_{DM}/M_{*} = 100$ (see details in Table. 1). To generate our live objects, we use the initial condition code MAGI (Miki & Umemura 2018). Adopting a distribution-function-based method ensures that the final realization of the galaxy is in dynamical equilibrium (Miki & Umemura 2018). We perform our simulations with the high-performance collisionless N-body code GOTHIC (Miki & Umemura 2017). This gravitational octree code runs entirely on GPU and is accelerated by the use of hierarchical time steps in which a group of particles has the same time step (Miki & Umemura 2017). We evolve the M31 galaxy-satellite system over a few Gyr depending on the scenario. We set the particle resolution of all the live objects to $4.4 \times 10^4 M_{\odot}$ and the gravitational softening length to 10 pc. The softening value was estimated using the following criterion: $\epsilon \sim a/N^{1/3}$, where N and a are the number of particles and the scale length of the component, respectively. As all components in our simulations share the same softening, we adopt the smaller value, which is given by the stellar component of the satellite. As in GADGET-2 (Springel 2005), the acceleration Multipole Acceptance Criterion is employed in our GPU code (Miki & Umemura 2017). The time-step parameter $\eta$ is set to 0.5.

To compare with previous simulations, we comment that Sadoun et al. (2014) used NBODYGEN for initial conditions and

GADGET-2, here we test our models using our fully GPU code GOTHIC and MAGI as the initial condition generator, which allows us to achieve a 100 times higher mass resolution. This higher resolution enables us to study the impact of the satellite on the central DM density profile of M31.

### 4. Tests of the Models with M31 observations

We consider five scenarios for the formation of the GSS of M31. The original empirical scenario is that in which the satellite has no DM (Fardal et al. 2007) and the cosmologically motivated scenario in which the satellite is DM rich (Sadoun et al. 2014) (see Table 1). Many models of formation and evolution of galaxies favor a larger halo-to-stellar mass ratio $M_{DM}/M_{*}$ of about 100 whereas the model proposed by Sadoun et al. (2014) favors a lower $M_{DM}/M_{*}$ of 20. Interestingly, this implies that the satellite lies below the halo-to-stellar mass relation as galaxies in this mass range should have a mean halo mass of $\simeq 2 \times 10^{11} M_{\odot}$ (Behroozi et al. 2010; Moster et al. 2013; Henriques et al. 2015).

To conform with these observations of dwarf galaxies, we further ran three simulations with $M_{DM}/M_{*} = 100$ and with different density profiles and initial infall radii, but keeping $10^9 M_{\odot}$ for the stellar component as the mass of the GSS is well constrained by the observations. In these simulations, the satellite has a cuspy profile as summarized in Table 1.

The orbital evolution of satellites in these five scenarios up to their passages through the center of M31 and their subsequent disruptions is shown in Figure 1. The orbital radius is color-
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Orbital evolution of the satellite: orbital radius as a function of time in all scenarios (see Table 1). The orbital radius is color-coded according to the bound stellar mass in the satellite at each time. Satellites with ($M_{\text{DM}}/M_{*\text{sat}}$) = 0 and 20 and satellite A have reached their pericenter at 0.18, 1.6, and 1.47 Gyr and were completely tidally stripped after 0.37, 1.93, and 1.87 Gyr, respectively. However, models B and C exhibit the presence of a remnant of the satellite after 1.86 and 4.54 Gyr, respectively.</figcaption></figure>

coded according to the bound stellar mass in the satellite at each time. We follow the iterative method of Baumgardt & Makino (2003) to determine the number of bound stars over time. Because all our models are stringently constrained by the observational data on the GSS and shells of M31, we stop the simulations when the best agreements with these kinematic data are achieved. Hence, some of these mergers are constrained to be more ancient than others as clearly depicted both in Figures 1 and 2. Satellite galaxies with $M_{\text{DM}}/M_{*} = 0$ and 20 have reached their pericenter at 0.18 and 1.6 Gyr and were completely tidally stripped within 0.37 and 1.93 Gyr, respectively. Satellites with ($M_{\text{DM}}/M_{*\text{sat}}$) = 0 and 20 and satellite A have reached their pericenters at 0.18, 1.6, and 1.47 Gyr and were completely tidally stripped after 0.37, 1.93, and 1.87 Gyr, respectively. However, models B and C exhibit the presence of a remnant of the satellite after 1.86 and 4.54 Gyr, respectively.

In our best-fit model, Sadoun et al. (2014), the simulation goes back to 2.1 Gyr ago and that is why we cannot exclude the cosmological turnaround of the infalling satellite as part of its real orbit. Our $r_{\text{peri}}/r_{\text{apo}}$ matches quite well the result of Lux et al. (2010) if we evaluate this ratio for the follow-up apocenters as is done in Lux et al. (2010) and is demonstrated for our model in Figure 1. We find that $r_{\text{peri}}/r_{\text{apo}}$ is about 0.1, which is exactly within the domain of values given by Lux et al. (2010).

For models with $M_{\text{DM}}/M_{*} = 100$, cases A and C (see Table 1) do not match well the kinematic data. We find that the DM in the satellite prevents the stellar component from being tidally stripped after the first pericenter passage and consequently the resulting GSS is far less massive than required by the observations. Subsequent passages will be necessary to completely strip the satellite but once again fail to reproduce the observations (see Figure 2). For model B in which the satellite has the same initial orbital conditions as in Sadoun et al. (2014) but has a much larger concentration parameter, the mass of the GSS is slightly overestimated but the other data are well matched (see Figure 2). For the Sadoun model, we find $M_{\text{GSS}} = 2.15 \times 10^8 M_{\odot}$ in excellent agreement with the value

of $2.4 \times 10^8 M_{\odot}$ derived from observations (Ibata et al. 2001; Fardal et al. 2006).

In addition to fitting the mass of the GSS, we select our best-fit model by making a detailed comparison with other observations of the stellar stream for all five scenarios (see Table 1). Figure 2 shows the stellar density maps in standard sky coordinates corresponding to stars of the satellite. We represent the observed stream fields as solid rectangles with proper scaling. We note that the simulated stream is in good agreement with the observations regarding the morphology and spatial extent of the GSS only for the Sadoun model and to a lesser extent for model B (see Figure 2). Furthermore, we test the five models against kinematic data. Figure 2 shows the radial velocities of satellite particles as a function of distance along the stream. We obtain good agreement with observations of the radial velocity measurement in the five fields for the Sadoun model and model B (Ibata et al. 2004; Fardal et al. 2006).

The small halo-to-stellar mass ratio adopted in our best-fit model (Sadoun et al. 2014) can be puzzling in the general context of formation and evolution of galaxies, which in general prefers a larger DM mass for the satellite of stellar mass in the range used in this model ($M_{*} = 2 \times 10^9 M_{\odot}$). However, there are many incidences where a small halo-to-stellar mass ratio as such is favored. It is widely believed that dwarfs formed in tidal tails have a low DM content (Guo et al. 2020). It has also been proposed that dwarfs that form in high-velocity dwarf mergers are DM deficient (Shin et al. 2020). There are also cases of dwarf galaxies that grow passively (Cowley et al. 2019). It has also been argued that many of the dwarfs of the Milky Way are DM poor (Hammer et al. 2020). These are just a few arguments in favor of the low halo-to-stellar mass ratio used in the Sadoun model. Here, we emphasize that a large halo-to-stellar mass ratio was dismissed by Sadoun et al. (2014) because it would thicken the disk of M31 by more than is allowed by observations. Here we have run three further simulations with $M_{\text{DM}}/M_{*} = 100$ and have shown that one could not satisfy all observational constraints set by the GSS as well as in the Sadoun model. Indeed, model B, which uses identical initial orbital conditions to Sadoun et al. (2014) and which the merger dates back to about 2 Gyr as in Sadoun et al. (2014), fits well the observations but slightly overestimates the mass of the GSS.

# 5. Results: A Model-independent Cusp-to-core Transition in M31

## 5.1. Heating and Migration of Dark Matter in M31

In this subsection, we study the heating and migration of DM in M31 for our best-fit model (Sadoun et al. 2014). Unlike the previous work of Sadoun et al. (2014), which lacked resolution and computational power, our GPU simulations enable us to study the impact of the accretion of a DM-rich satellite on the DM distribution in the central regions of M31. In order to determine properly the M31 center, we use the method of Power et al. (2003) by applying the shrinking sphere method to M31 (DM halo and the satellite), the M31 halo, and the M31 bulge and have found similar centers at each time. As discussed in Power et al. (2003), this iterative technique localized efficiently the densest region in halos, which are distorted by mergers. We have further checked our result against this method by finding particles with the lowest potential as the
<figure><img src="image_5.png" /></figure>

**Figure 2.** Satellite models constrained by the GSS. Top panels: simulated stellar density maps in standard sky coordinates corresponding to stars of the satellite at different times. We represent the observed stellar fields with black boxes (McConnachie et al. 2003). We find $M_{GSS} = 2.15 \times 10^8 M_{\odot}$ for the Sadoun model in good agreement with the value of $2.4 \times 10^8 M_{\odot}$ derived from observations with a mass-to-light ratio of 7 (Ibata et al. 2001; Fardal et al. 2006). However, for models with large $M_{DM}/M_{*} = 100$, the mass of the GSS is overestimated in model B, and models A and C do not match any of the observations. Bottom panels: simulated radial velocity of satellite particles as a function of the distance along the stream at different times when the best observations are obtained. We represent the radial velocity measurements in five fields by red points with error bars (Ibata et al. 2004; Fardal et al. 2006). A good agreement with observations for the radial velocity measurement is shown for the Sadoun model and to a lesser extent for model B.

<figure><img src="image_6.png" /></figure>

**Figure 3.** Energy transfer via dynamical friction. In the left panel, we show the evolution of the energy distribution for DM particles of M31 within a fixed radius of 500 pc in the Sadoun scenario (see Table 1). The times for the snapshots are chosen because the pericenter passage occurs at about 1.6 Gyr and the simulation is stopped at 2.1 Gyr when the best agreements with observations are achieved. The figure confirms that not only the DM particles that remain within 500 pc have gained energy but also some of them have migrated out of this region. In the right panel, we follow the DM particles that were initially inside the 500 pc radius and show the evolution of the distribution of the energy at the same time intervals as the left panel. Hence, in this figure, unlike the left panel, the number of DM particles is fixed. The histogram clearly shows that most of the particles that were initially inside the 500 pc radius, whether migrated or not, have heated up, which could plausibly be a mechanism for core formation.

center of the halo in order to have a robust determination of the center. In the end, all our simulations are centered on the densest region of the M31 halo.

As the satellite has a radial orbit, its crossings near the center perturb the M31 halo by heating its DM particles via dynamical friction. The pericenter passage of the Sadoun satellite occurs at 1.6 Gyr. In our merger scenarios, the two DM halos interact gravitationally via dynamical friction. At the particle scale, this means that energy transfers occur between DM particles. This can be seen as a heating of DM particles of M31 from the DM particles of the satellite. Indeed, the energy and angular momentum of the satellite galaxy are consumed by dynamical friction. The result is that these energy transfers are going to expand the orbits of DM particles of M31. However,

at kiloparsec scales, we expect to observe the emergence of collective effects and potential fluctuations, which erode the central density cusp of the M31 halo: the transfer of the gravitational energy from the satellite to M31 perturbs the DM distribution of the host halo.

Figure 3 depicts the total energy distribution for DM particles of M31, which are within the 500 pc (left panel), and in the right panel, we follow the particles that were initially within 500 pc at different times in the Sadoun scenario. As the number of particles is not conserved in the left panel of Figure 3, it confirms that some DM particles have gained energy and have migrated outwards from the central region of M31. The right panel demonstrates that the particles that were initially within 500 pc have gained energy.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Migration and heating of DM in the central region of M31 during the passage of a satellite. Density-weighted projected mass (top panel) and total energy maps (bottom panel) of a slice of thickness 200 pc for M31 DM particles at different times in the Sadoun scenario (see Table 1). The top panel shows the migration of DM particles in the central region of M31, especially after the pericenter passage of the satellite (1.6 Gyr). The bottom panel shows the heating of the DM particles that remains within the 500 pc radius central region.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Change in the potential profile and flattening of the density profile of M31. The profile of the DM potential of M31 normalized by the norm of its initial central value $|\Phi_0|$ at different times in the Sadoun scenario (see Table 1) is shown in the left panel. The spherically averaged DM density profile of M31 in the Sadoun scenario in 26 pc thick radial shells at different times is shown in the right panel. Initially, the M31 DM halo assumes an NFW profile (black dashed curve). The convergence radius of 73 pc for the M31 DM halo is shown by the vertical dashed green line below which the simulations cannot be considered to be fully converged according to the criterion of Power et al. (2003). Beyond this region, the flattening of the cusp over almost one decade is evident. In the absence of a satellite, the DM density profile remains cuspy. This scenario ensures the stability of the M31 profile against numerical effects (red dotted curve). The times for the snapshots are chosen because the pericenter passage occurs at about 1.6 Gyr and the simulation is stopped at 2.1 Gyr when the best agreements with observations are achieved. Both panels demonstrate the gradual flattening of the initial cuspy profile.</figcaption></figure>

To further confirm the mechanism of migration and heating behind the core formation, we also show density-weighted projected mass maps in Figure 4. The top panel of this figure clearly demonstrates the migration from the central region as the mass within the 500 pc radius decreases after the pericenter passage of the satellite at 1.6 Gyr. Not only the migrated DM particles but also those that remain within the central 500 pc have gained energy through dynamical friction, which slows down the satellite.

The heating and particle migration in the central region of M31 are expected to lead to the flattening of the DM density profile. In Figure 5, we demonstrate the change in the potential (left panel) and the flattening of the cuspy initial density profile of M31 after the passage of the satellite (at 1.6 Gyr). Initially, the M31 DM halo assumes an NFW profile. We consider DM particles from both M31 and satellite halos to determine the spherically averaged DM density profile of M31 in 26 pc thick radial shells over time. Below a radius of 73 pc, the simulations are affected by two-body relaxation (Power et al. 2003), and
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Impact of softening on the DM profile of M31. The DM density profile of M31 halo for $\epsilon = 10$ and 5 pc in the Sadoun scenario (see Table 1). As the softening length $\epsilon$ does not affect the DM density profile, our simulation results are thus free from such numerical artifacts. The convergence radius of 73 pc for the M31 DM halo is shown by the vertical dashed green line below which the simulations cannot be considered to be fully converged according to the criterion of Power et al. (2003)</figcaption></figure>

beyond this region, Figure 5 clearly demonstrates the flattening of the cusp over more than one decade in radius.

The flattening of the initially cuspy profile can also be seen in the change of the potential energy of M31. The profile of the potential of M31 plotted in the left panel of Figure 5 is determined by computing the gravitational potential at different radii using all DM particles. The potential at a specified position is obtained by summing all interactions with the surrounding point masses with PNBODY (Revaz 2013).

The migration of DM particles from the central regions of M31 is expected to lead to the shallowing of the potential well, which is clearly demonstrated in Figure 5. The two panels of Figure 5 together with Figures 3 and 4 show clearly that energy transfer from the satellite to M31 and subsequent outward migrations of DM particles are responsible for the flattening of the density profile.

## 5.2. Convergence Tests: Two-body Relaxation, Softening, and Mass Resolution

In this subsection, we study the numerical artifacts that could render our result nonphysical. We have conservatively applied the Power et al. (2003) criterion to our DM halos to estimate the radius within which the two-body relaxation time is shorter than the simulation time. We find a convergence radius of 73 pc for the M31 DM halo. The green vertical dashed line in Figures 5, 6, and 8 marks the radius below which the simulations are not fully converged according to the criterion of Power et al. (2003).

Next, we study the effect of the softening parameter, which models the interaction between two Plummer point masses. In order to test and ensure that our simulation results are robust and reliable, we have performed softening tests. Figure 6 shows the DM density profile for a set of simulations that only differ in the value of the softening length $\epsilon$. This figure demonstrates that doubling the softening radius leaves the density profile intact.

Furthermore, we study the effect of mass resolution and run our simulations for three different cases. We see in Figure 7

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Impact of mass resolution on the DM profile of M31. DM density profile of the M31 halo for three different mass resolutions is shown for the Sadoun scenario (see Table 1). The convergence radii below which the simulations cannot be considered to be fully converged (Power et al. 2003) change with the mass resolution and are shown by the vertical dashed lines. We see that convergence is achieved for a mass resolution smaller than $8.8 \times 10^4 M_{\odot}$.</figcaption></figure>

that for a very large mass resolution of $4.4 \times 10^5 M_{\odot}$ the convergence radius of Power et al. (2003) shown by the vertical orange line moves to the right and hence the change in the density is only reliable over a very small range of radii. However, we see a convergence between the mass resolution of $8.8 \times 10^4 M_{\odot}$ and $4.4 \times 10^4 M_{\odot}$.

We emphasize that the flattening of the density profile is observed here over more than a decade between 73 pc and 1 kpc. A few previous studies (e.g., Boylan-Kolchin & Ma 2004) using simple models of head-on collisions have concluded that the cusps remain cusps in mergers. However, such studies could not resolve scales below 1 kpc (e.g., Figure 4 in Boylan-Kolchin & Ma 2004), which is precisely the scales at which the cusp profile flattens. This study was also not carried out in the original work of Sadoun et al. (2014) again because of lack of numerical resolution.

In the next subsection, we provide fits to the newly flattened profiles and study the same migration and heating mechanism in other models of the formation of GSS in M31.

## 5.3. Model-independent Core-like Dark Matter Profile in M31

In order to affirm the model-independent nature of the cusp-to-core transition in M31, we have studied all of our five models (Table 1).

We consider DM particles from both M31 and satellite halos to determine the spherically averaged DM density profile for all our five models. As shown in Figure 8, we find that our profile is well fitted by the following four-parameter formula (Hernquist 1990; Saha 1992; Dehnen 1993; Zhao 1996):

$$ \rho(r) = \frac{\rho_c}{(1 + (r/r_c)^{1/\beta})^{\gamma}}, \quad (1) $$

where $\rho_c$ is the core constant density and $r_c$ is the core radius we find to be about 1.1 kpc for the M31 halo in the Sadoun scenario (see Table 1). We stress that our best-fit core radii are larger than the numerical convergence radius (73 pc) and persist over almost one decade. However, as the value of the core size depends on the fitted DM profile, we defer from
<figure><img src="image_11.png" /></figure>

**Figure 8.** Model-independent core-like profile in M31. Spherically averaged DM density profile for all our five scenarios in 26 pc thick radial shells (see Table 1). In all scenarios, the satellite has a stellar mass $M_* = 2.2 \times 10^9 M_\odot$. Initially, the M31 DM halo assumes an NFW profile (black dashed line). We consider DM particles from both the M31 and satellite halos to determine the DM density profile of M31. The fitting function described by Equation 1 reproduces the simulated density structures and captures the rapid transition from the cusp to the core. We set Poissonian errors for fitting weights. We stress that our best-fit core radii are larger than the numerical convergence radius marked by the vertical dashed green line. We observe a DM core of about 1.1 kpc for the M31 halo in our best-fit model (Sadoun et al. 2014).

imposing the precise value of the radius as a constraint on the M31 core, and here we only demonstrate that a merger with a wide range of halo-to-stellar mass ratio can indeed flatten the central density profile.

A central question in our work is whether or not the initial density profile has become so much shallower that it can be called a core. There are various definitions for a core. The simplest evidence for a core is the visual change in the density profile. Figure 8 clearly shows that the density profile above the convergence radius of 73 pc has been flattened and can no longer be fitted by the same NFW profile. More in-depth evidence is given by the value of $r_c$ in Equation 1, which is larger than the convergence radius and defined by the radius at which the constant density is approximately divided by 2. The third piece of evidence is provided by direct measurement of the slopes of the density profiles as shown in Figure 9. Different ranges of values of these slopes have been used to define a central core (de Blok et al. 2001; Oh et al. 2008, 2011, 2015). Figure 9 clearly shows a flattening of the slope. For the Sadoun model and larger $M_{DM}/M_*$, we clearly see a slope larger than -0.5 spanning over a wide range. Whether this flattening continues to smaller radii below our convergence radius is a question that will be pursued in forthcoming work where we will exploit more powerful GPU resources.

For the Fardal model (Fardal et al. 2007), the density profile at the end of the simulation, which is at 0.85 Gyr when the best match with observations is obtained, in spite of the satellite being far less massive in this scenario, is well fit by a core-like profile (see Figure 8). However, the core size is relatively smaller as expected because the heating and migration mechanisms are less effective for the low-mass satellite. As the satellite in the Sadoun model is 20 times more massive than that in the Fardal model, the core size is larger in this scenario.

<figure><img src="image_12.png" /></figure>

**Figure 9.** Core-like profiles of M31. Inner slope of the DM density profile of M31 is plotted against the radii. In all scenarios, the satellite has a stellar mass $M_* = 2.2 \times 10^9 M_\odot$. It clearly shows a flattening of the slope. For the Sadoun model and larger $M_{DM}/M_*$, we clearly see a slope larger than -0.5 spanning over a wide range, which is often used as a definition of a core profile (de Blok et al. 2001; Oh et al. 2008, 2011, 2015).

In the other scenarios with a larger halo-to-stellar mass ratio, a cusp-to-core transition also occurs (see Figure 8).

It has been proposed that a major merger scenario could explain the kinematic data and stellar abundance data of M31 (Hammer et al. 2010, 2018; D'Souza & Bell 2018). Here, we demonstrate that a merger of satellites with $M_{DM}/M_* = 100$, which can be considered as major mergers, can also lead to the flattening of the initially cuspy profiles and yield new profiles, which are well fitted by core-like functions (see Figure 8). Consequently, a major merger scenario with a satellite on a highly eccentric orbit should also trigger a cusp-to-core formation for the M31 halo. Besides, the resulting DM core is expected to be larger as larger perturber masses lead to a larger constant-density central region (Goerdt et al. 2010), which is shown by model C (green curve in Figure 8). Model B, which conforms best with observations for $M_{DM}/M_* = 100$, also shows a clear flattening of the DM density profile of M31 (see Figure 8). In spite of the fact that models A and C cannot reproduce the observational data on the GSS, they disturb the density profile of the host galaxy and flatten the central cusp, which hints that such a transition can be a generic outcome of most close-encounter merger scenarios.

In addition to the satellite mass and orbit, another dynamical parameter, which impacts the flattening of the cusp, is the concentration of the DM and star distributions within the satellite. Even if the satellite in models A, B, and C is more massive than Fardal's satellite, the satellite DM distribution in these three models is more diffuse compared to the one of Fardal, which is only formed by stars. Indeed, it is harder to slow down a satellite with an extended DM distribution as the mass loss of the satellite can easily compensate for the dynamical friction from the DM background. Consequently, the energy transferred to the DM central region of M31 is reduced, and this explains why the DM density profiles only differ by a factor of 2 at the convergence radius between these models (see Figure 8). Even if dynamical friction increases in proportion to the satellite mass, the DM halo at a large radius is not sufficiently dense compared to the central region to significantly slow down the satellite. The energy transferred via dynamical friction by the satellite to the DM outer regions
of M31 is found to be negligible as depicted by the DM density profiles in Figure 8.

## 5.4. Comparison with Cosmological Simulations

The fact that all of our models with $M_{DM}/M_{*}$ ranging from 0 to 100, i.e., spanning minor to major mergers, redistribute DM in the central part of M31 and are well fit by a core-like profile is strong evidence for the universality of our result. We expect the flattening of the initial cuspy profile to occur universally at least in nearly radial close-encounter merger scenarios as discussed here. However, there is widespread consensus that the hierarchical assembly of CDM halos yields a cuspy DM profile in DM-only simulations (Navarro et al. 1996, 1997; Moore et al. 1999; Fukushima & Makino 2001; Navarro et al. 2004; Diemand et al. 2004). These cosmological simulations find cuspy halos as we do when we run simulations with low mass resolution (see Figure 7). In such cases, there is enough resolution to resolve the energy transfers via dynamical friction between particles in the central region of galaxies. In the absence of particles in the central regions, the energy transfer induced by the dynamical friction is not properly modeled. The mass resolution is directly related to the convergence radius, below which we cannot resolve the physics in the simulation. As demonstrated by Navarro et al. (2010) and our Figure 7, it is necessary to have mass resolution lower than $10^5 M_{\odot}$ to investigate properly the physics below the kiloparsec scale. For instance, in the Illustris TNG simulations, the central regions of galaxies are not resolved as the best resolution in these simulations is $4.5 \times 10^5 M_{\odot}$, achieved in the latest TNG 50 simulation (Pillepich et al. 2019).

Indeed, the DM density profile of resimulated halos in zoom simulations reveals that the density profiles become shallower inwards down to the innermost resolved point (Navarro et al. 2010). Here we are demonstrating that the subhalo accretion could be at the origin of such flattening of the NFW profile. Even if subhalo orbits become more radial and plunge deeper into their host halo at higher host halo mass (Wetzel 2011), the cuspy profile persists if the mass of the infalling satellite is at the lower end of the mass range. In our M31 scenario, however, the flattening of the cusp is more easily detected because the satellite is massive and falls in a highly radial orbit.

Moreover, we highlight also that the stellar component of the satellite plays a major role in core formation. Indeed, as this component is more concentrated compared to the DM of the satellite (see Table 1), it will further slow down the satellite during its infall and thus disturb the central region of M31 even more prominently. Figure 10 compares the DM density profile of M31 after the infall of a satellite with and without a stellar component in the Sadoun scenario (see Table 1). The stellar component of the satellite has a nonnegligible contribution to the flattening of the DM cusp, and this clearly demonstrates that the stellar component cannot be ignored in the cosmological simulations which aim at determining the density profiles of DM halos.

# 6. Conclusion

The GSS and the shell-like features of M31 are likely outcomes of accretion of satellites on highly eccentric orbits as shown by numerical simulations in the past and reaffirmed here. In this work, we have mainly studied the model proposed

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Impact of stars from the satellite. DM density profile of M31 after the infall of satellite with (without) stars in 26 pc thick radial shells in the Sadoun scenario (see Table 1). Initially, the M31 DM halo assumes an NFW profile (black dashed line). We stress that the stellar component of the satellite has a nonnegligible contribution to the flattening of the DM cusp.</figcaption></figure>

by Sadoun in 2014 in which the Satellite is DM rich and falls from its turnaround radius onto M31.

However, we have increased the resolution of the previous simulations by using a fully GPU code. The substantial gain of 100 in mass resolution has enabled us to here study the impact of such an accretion event on the spatial distribution of DM in the central regions of M31. Our simulations show that as the satellite falls onto M31, it is slowed down by dynamical friction and its energy is transferred to the host halo. The DM particles in the central regions of M31 are heated and migrate outwards. Here we have shown that in this process, the initial cusp becomes shallower over almost a decade in radius and is well fit by a core-like profile.

To explore the model-independent aspects of our results, we have studied different models with halo-to-stellar mass ratios ranging from 0 to 100 corresponding to minor and major merger scenarios. An ancient major merger scenario has been proposed for M31 (Hammer et al. 2010), which however did not reproduce the observed features of the GSS. A new version of this model in which the merger occurs 2 Gyr ago (as in the Sadoun model) has been proposed recently (Hammer et al. 2018), which has been taken further to also account for the origin of M32 (D’Souza & Bell 2018) but has also been proposed in the past (Block et al. 2006). Major mergers are quite rare, especially recently, and it has been shown that they could thicken the disk of galaxies by more than is allowed by observations. It is customary to study the stars from the merger but the more massive DM component could churn up the host galaxy and destroy any thin disk (e.g., Toth & Ostriker 1992).

We infer that merger events in which galaxies fall on highly eccentric orbits into their host halos can provide a general mechanism for the flattening of the density profile in a $\Lambda$CDM universe where halos are expected to be cuspy. It has been reported that satellites in host halos with larger mass ratios have slightly more eccentric orbits with lower angular momentum, and moreover, satellites around more massive halos seem to be on more radial orbits at fixed mass ratio (Tormen 1997; Wetzel 2011; Jiang et al. 2015). Hence, we expect a noticeable fraction of galaxies in a $\Lambda$CDM universe to harbor core-like profiles that have been formed during merger and accretion events. In this work, we have ignored the presence of a central
black hole and have assumed an initially NFW profile for M31. The initially cuspy profile must have been even cuspier, due to the redistribution of DM around black holes (Bahcall & Wolf 1976; Gondolo & Silk 1999). We however expect that the mechanism discussed here will not be influenced by a central black hole. A study of the fate of the black hole and the change of the central profile during the merger event studied here will be presented in a forthcoming work. Present observational data from stellar kinematic observations and the disk's H I rotation curve have been inconclusive about the DM distribution in M31 (e.g., Tamm et al. 2012; Blaña Díaz et al. 2018). Here, we have shown that GSS-constrained modeling of M31 favors a shallow cusp or a core-like inner DM distribution in M31. Future observations will shed light on this interesting question for our nearest neighbor.

Similarly, it was demonstrated that the recent passage of the Large Magellanic Cloud (LMC) has perturbed the MW disk (Petersen & Peñarrubia 2020) and also distorted the structure of the DM halo (Belokurov et al. 2019; Garavito-Camargo et al. 2019; Erkal et al. 2020). In contrast with our M31 merger scenario, the heating from the LMC seems to particularly affect the outer regions of the MW halo as the LMC has a much larger pericenter (~50 kpc) than our M31 satellite. With more powerful GPU resources, we could investigate the impact of the LMC infall on the inner regions of the DM halo of the MW.

We thank the reviewer for constructive feedback, which helped to improve the quality of the manuscript. We thank Henry McCracken for constructive comments and discussions. We also thank Miki Yohei for providing us with the nonpublic N-body code, GOTHIC.

## ORCID iDs

Joe Silk https://orcid.org/0000-0002-1566-8148

## References

Adams, J. J., et al. 2014, ApJ, 789, 63

Agnello, A., & Evans, N. W. 2012, ApJL, 754, L39

Amorisco, N. C., & Evans, N. W. 2012, MNRAS, 419, 184

Bahcall, J. N., & Wolf, R. A. 1976, ApJ, 209, 214

Battaglia, G., Helmi, A., Tolstoy, E., et al. 2008, ApJL, 681, L13

Baumgardt, H., & Makino, J. 2003, MNRAS, 340, 227

Behroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, ApJ, 717, 379

Bellazzini, M., Cacciari, C., Federici, L., Fusi Pecci, F., & Rich, M. 2003, A&A, 405, 867

Belokurov, V., Deason, A. J., Erkal, D., et al. 2019, MNRAS, 488, L47

Benítez-Llambay, A., Frenk, C. S., Ludlow, A. D., & Navarro, J. F. 2019, MNRAS, 488, 2387

Blaña Díaz, M., Gerhard, O., Wegg, C., et al. 2018, MNRAS, 481, 3210

Block, D. L., et al. 2006, Natur, 443, 832

Boldrini, P., Miki, Y., Wagner, A. Y., et al. 2020, MNRAS, 492, 5218

Boldrini, P., Mohayaee, R., & Silk, J. 2020, MNRAS, 492, 3169

Boylan-Kolchin, M., & Ma, C.-P. 2004, MNRAS, 349, 1117

Brown, T. M., et al. 2006, ApJ, 652, 323

Conn, A. R., et al. 2013, ApJ, 766, 120

Cowley, W. I., Caputi, K. I., Deshmukh, S., et al. 2019, ApJ, 874, 114

D'Souza, R., & Bell, E. F. 2018, NatAs, 2, 737

de Blok, W. J. G., McGaugh, S. S., Bosma, A., & Rubin, V. C. 2001, ApJL, 552, L23

Dehnen, W. 1993, MNRAS, 265, 250

Dekel, A., Devor, J., & Hetzroni, G. 2003, MNRAS, 341, 326

Diemand, J., Moore, B., & Stadel, J. 2004, MNRAS, 353, 624

El-Zant, A., Shlosman, I., & Hoffman, Y. 2001, ApJ, 560, 636

Erkal, D., Belokurov, V. A., & Parkin, D. L. 2020, MNRAS, 498, 5574

Fardal, M. A., Babul, A., Geehan, J. J., & Guhathakurta, P. 2006, MNRAS, 366, 1012

Fardal, M. A., Guhathakurta, P., Babul, A., & McConnachie, A. W. 2007, MNRAS, 380, 15

Fardal, M. A., et al. 2012, MNRAS, 423, 3134

Ferguson, A. M. N., Irwin, M. J., Ibata, R. A., Lewis, G. F., & Tanvir, N. R. 2002, AJ, 124, 1452

Ferguson, A. M. N., et al. 2005, ApJL, 622, L109

Flores, R. A., & Primack, J. R. 1994, ApJL, 427, L1

Font, A. S., Johnston, K. V., Guhathakurta, P., Majewski, S. R., & Rich, R. M. 2006, AJ, 131, 1436

Fukushige, T., & Makino, J. 1997, ApJL, 477, L9

Fukushige, T., & Makino, J. 2001, ApJ, 557, 533

Gao, L., & White, S. D. M. 2007, MNRAS, 377, L5

Garavito-Camargo, N., Besla, G., Laporte, C. F. P., et al. 2019, ApJ, 884, 51

Goerdt, T., Moore, B., Read, J. I., & Stadel, J. 2010, ApJ, 725, 1707

Gondolo, P., & Silk, J. 1999, PhRvL, 83, 1719

Guo, Q., Hu, H., Zheng, Z., et al. 2020, NatAs, 4, 246

Guo, Q., White, S., Boylan-Kolchin, M., et al. 2011, MNRAS, 413, 101

Hammer, F., Yang, Y., Arenou, F., et al. 2020, ApJ, 892, 3

Hammer, F., Yang, Y. B., Wang, J. L., et al. 2010, ApJ, 725, 542

Hammer, F., Yang, Y. B., Wang, J. L., et al. 2018, MNRAS, 475, 2754

Hayashi, E., & Navarro, J. F. 2006, MNRAS, 373, 1117

Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2015, MNRAS, 451, 2663

Hernquist, L. 1990, ApJ, 356, 359

Hernquist, L., & Quinn, P. J. 1988, ApJ, 331, 682

Hernquist, L., & Quinn, P. J. 1989, ApJ, 342, 1

Ibata, R., Chapman, S., Ferguson, A. M. N., et al. 2004, MNRAS, 351, 117

Ibata, R., Chapman, S., Ferguson, A. M. N., et al. 2005, ApJ, 634, 287

Ibata, R., Irwin, M., Lewis, G., Ferguson, A. M. N., & Tanvir, N. 2001, Natur, 412, 49

Ibata, R., Martin, N. F., Irwin, M., et al. 2007, ApJ, 671, 1591

Ibata, R. A., et al. 2013, Natur, 493, 62

Inoue, S., & Saitoh, T. R. 2011, MNRAS, 418, 2527

Jardel, J. R., & Sellwood, J. A. 2009, ApJ, 691, 1300

Jiang, L., Cole, S., Sawala, T., & Frenk, C. S. 2015, MNRAS, 448, 1674

Kormendy, J., & Bender, R. 1999, ApJ, 522, 772

Laporte, C. F. P., & Penarrubia, J. 2015, MNRAS, 449, L90

Laporte, C. F. P., & White, S. D. M. 2015, MNRAS, 451, 1177

Lux, H., Read, J. I., & Lake, G. 2010, MNRAS, 406, 2312

Ma, C.-P., & Boylan-Kolchin, M. 2004, PhRvL, 93, 021301

Malin, D. F., & Carter, D. 1980, Natur, 285, 643

Malin, D. F., & Carter, D. 1983, ApJ, 274, 534

McConnachie, A. W., Irwin, M. J., Ibata, R. A., et al. 2003, MNRAS, 343, 1335

McConnachie, A. W., et al. 2009, Natur, 461, 66

Miki, Y., & Umemura, M. 2018, MNRAS, 475, 2269

Miki, Y., & Umemura, M. 2017, NewA, 52, 65

Moore, B. 1994, Natur, 370, 629

Moore, B., Ghigna, S., Governato, F., et al. 1999, ApJL, 524, L19

Moore, B., Governato, F., Quinn, T., Stadel, J., & Lake, G. 1998, ApJL, 499, L5

Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Navarro, J. F., Hayashi, E., Power, C., et al. 2004, MNRAS, 349, 1039

Navarro, J. F., et al. 2010, MNRAS, 402, 21

Neistein, E., & Dekel, A. 2008, MNRAS, 388, 1792

Ogiya, G., & Mori, M. 2011, ApJL, 736, L2

Oh, S.-H., de Blok, W. J. G., Brinks, E., Walter, F., & Kennicutt, R. C. 2011, AJ, 141, 193

Oh, S.-H., de Blok, W. J. G., Walter, F., Brinks, E., & Kennicutt, R. C. 2008, AJ, 136, 2761

Oh, S.-H., Hunter, D. A., Brinks, E., et al. 2015, AJ, 149, 180

Oman, K. A., Marasco, A., Navarro, J. F., et al. 2019, MNRAS, 482, 821

Oman, K. A., et al. 2015, MNRAS, 452, 3650

Petersen, M. S., & Peñarrubia, J. 2020, NatAs, 5, 251

Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196

Pontzen, A., & Governato, F. 2012, MNRAS, 421, 3464

Power, C., et al. 2003, MNRAS, 338, 14

Read, J. I., Iorio, G., Agertz, O., & Fraternali, F. 2016, MNRAS, 462, 3628

Revaz, Y. 2013, pNbody: A python parallelized N-body reduction toolbox, Astrophysics Source Code Library, ascl:1302.004

Richardson, J. C., et al. 2008, AJ, 135, 1998

Robertson, A., et al. 2018, MNRAS, 476, L20

Rocha, M., Peter, A. H. G., Bullock, J. S., et al. 2013, MNRAS, 430, 81
Romano-Díaz, E., Shlosman, I., Hoffman, Y., & Heller, C. 2008, ApJL, 685, L105

Sadoun, R., Mohayaee, R., & Colin, J. 2014, MNRAS, 442, 160

Saha, P. 1992, MNRAS, 254, 132

Santos-Santos, I. M. E., et al. 2020, MNRAS, 495, 58

Sawala, T., Frenk, C. S., Fattahi, A., et al. 2014, arXiv:1412.2748

Shin, E.-j., Jung, M., Kwon, G., et al. 2020, ApJ, 899, 25

Springel, V. 2005, MNRAS, 364, 1105

Tamm, A., Tempel, E., Tenjes, P., Tihhonova, O., & Tuvikene, T. 2012, A&A, 546, A4

Tanaka, M., Chiba, M., Komiyama, Y., et al. 2010, ApJ, 708, 1168

Teyssier, R., Pontzen, A., Dubois, Y., & Read, J. I. 2013, MNRAS, 429, 3068

Tonini, C., Lapi, A., & Salucci, P. 2006, ApJ, 649, 591

Tormen, G. 1997, MNRAS, 290, 411

Toth, G., & Ostriker, J. P. 1992, ApJ, 389, 5

Tremaine, S. 1995, AJ, 110, 628

Vogelsberger, M., Zavala, J., & Loeb, A. 2012, MNRAS, 423, 3740

Walker, M. G., & Peñarrubia, J. 2011, ApJ, 742, 20

Wetzel, A. R. 2011, MNRAS, 412, 49

Zavala, J., Vogelsberger, M., & Walker, M. G. 2013, MNRAS, 431, L20

Zhao, H. 1996, MNRAS, 278, 488

Zucker, D. B., et al. 2004, ApJL, 612, L117