# Building stellar bulges and halo cores from massive clumps observed in the
DYNAMO-HST sample

Mahmoud Hashim,1 Amr A. El-Zant,1 and Antonino Del Popolo2,3,4
1Centre for Theoretical Physics, The British University in Egypt, Sherouk
City, 11837 Cairo, Egypt
2Dipartimento di Fisica e Astronomia, University Of Catania, Viale Andrea
Doria 6, 95125, Catania, Italy
3 Institute of Astronomy and National Astronomical Observatory, Bulgarian
Academy of Sciences, 72, Tsarigradsko Shosse Blvd., 1784 Sofia, Bulgaria
4 Institute of Astronomy, Russian Academy of Sciences, 119017, Pyatnitskaya
str., 48 , Moscow, Russia mahmoud.hashim<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present N-body simulations of the process of bulge formation in disc
galaxies due to inward migration of massive stellar clumps. The process is
accompanied by dark halo heating, with a quasi-isothermal core replacing the
initial central density cusp, transforming an initially dark matter dominated
central region into a baryon dominated one. The characteristics of the clumps
are chosen to be compatible with low redshift observations of stellar clumps
in DYNAMO-HST galaxies, which may be relatively long lived in terms of being
robust against internal starburst-instigated disruption. We thus test for
disruption due to tidal stripping using different clump internal radial
profiles; Plummer, Hernquist and Jaffe, in ascending order of steeper central
density profile. Our calculations predict that in order for clump migration to
be effective in building galactic bulges and dark halo cores, steeply
increasing central clump profiles, or a less massive or less concentrated
haloes, are preferred. The dependence on such factors may contribute to the
diversity in observed total mass distributions and resulting rotation curves
in galaxies. When the process is most efficient, a ’bulge-halo conspiracy’,
with a singular isothermal total density akin to that observed bright
galaxies, results.

###### keywords:

galaxies: evolution – galaxies: haloes – galaxies: bulges – dark matter

††pubyear: 2023††pagerange: Building stellar bulges and halo cores from
massive clumps observed in the DYNAMO-HST sample–B

##  1 Introduction

As a model for cosmic structure formation, the standard cold dark matter (CDM)
based scenario has developed into a particularly effective paradigm describing
larger scales, but is still significantly lacking on smaller (sub)galactic
ones (e.g. Frenk & White, 2012; Primack, 2012). On these scales there has been
no shortage in controversy and apparent challenges, including longstanding
problems relating to the central densities of dark matter haloes and the
properties of subhaloes in the local universe (see Primack, 2012; Del Popolo &
Le Delliou, 2017; Bullock & Boylan-Kolchin, 2017; Salucci, 2019 for reviews).
This consists of a triplet of possibly related issues.

The first of these concerns the apparent discrepancy between the number of
subhaloes predicted by CDM-only N-body simulations compared to the number of
satellite galaxies observed (the so called ‘missing satellite problem"; Klypin
et al., 1999; Moore et al., 1999). The subhaloes obtained in CDM-only
cosmological simulations also appear too dense relative to those observed
around the Milky Way (Read et al., 2006; Garrison-Kimmel et al., 2013, 2014;
but see however Kravtsov & Wu, 2023), with observed subhaloes lacking
simulated counterparts. This has been labelled the ‘too-big-to-fail" problem
(Boylan-Kolchin et al., 2011). The third issue is concerned with the
discrepancy between the centrally divergent dark halo density profiles
obtained in the simulations ((Navarro et al., 1996b, 1997, 2010)) with the
observations of dwarf spirals, dwarf spheroidals, and low surface brightness
galaxies, showing cored profiles (e.g., Moore, 1994; Flores & Primack, 1994;
Burkert, 1995; de Blok et al., 2003; Swaters et al., 2003; Del Popolo, 2009;
Del Popolo & Kroupa, 2009; Del Popolo, 2012; de Blok et al., 2003; Swaters et
al., 2003; Del Popolo, 2009; Walker & Peñarrubia, 2011; Del Popolo & Hiotelis,
2014). This is habitually referred to as the Cusp/Core problem, and
constitutes a central aspect to be studied in the present work. It may also
lie at the heart of all three aforementioned problems; as these issues may in
fact be connected to the overdense self gravitating structures that the CDM
scenario produces (Ogiya & Burkert, 2015).

As the small scale problems affecting CDM-based structure formation have been
known for decades, several sets of solutions have been proposed, tackling the
problems from different fronts. One set relates to the nature of the dark
matter; popular alternatives to CDM, in this context, include warm dark matter
(Colín et al., 2000; Sommer-Larsen & Dolgov, 2001; Bode et al., 2001), self
interacting dark matter (Spergel & Steinhardt, 2000; Goodman, 2000), and fuzzy
dark matter (Peebles, 2000; Hu et al., 2000; Hui et al., 2017), as well as
less studied scenarios invoking direct interaction between the dark matter and
baryonic components in the Universe (Famaey et al., 2020; Salucci et al.,
2020). Some improvements may also be obtained by modifying the power spectrum
of primordial fluctuations (Kamionkowski & Liddle, 2000; Zentner & Bullock,
2003). Another possibility entails fundamental modification to gravitational
theory ((Buchdahl, 1970; Starobinsky, 1980; Milgrom, 1983a, b; Ferraro,
2012)), which naturally implies essential revision of the picture of the
Universe envisioned through the currently standard $\Lambda$CDM model.

As the problems arise at nonlinear scales, where complex baryonic physics
becomes important, it was natural to consider the effect of the physics of the
galaxy formation on the structure of the dark matter haloes. Solutions to the
cusp/core problem, in that context, propose that the central densities of
haloes may be reduced as the baryons settle in the centres of haloes. For
example, during a starburst supernovae driven gas can be blown out of the
central regions, and the central potential may consequently decrease. CDM
particles will then move to orbits of higher energies that take them farther
away from the centre, giving rise to a decrease in the central density, and a
flattening in the inner profile (as originally proposed by Navarro et al.,
1996a and more recently studied by Freundlich et al., 2020; Li et al., 2023).
Generally, a single blowout is insufficient to explain halo cores, and the gas
may in any case re-enter the central region, but the process of dynamically
’heating’ the central halo can be shown to be irreversible; leading to the
production of a core due to repeated blowout (e.g., Del Popolo, 2009; Pontzen
& Governato, 2012; a wider description of the process can be found in section
2.2 of Del Popolo & Pace, 2016). As supernovae shock waves (and gravitational
instabilities) may lead to fully developed turbulence in the feedback-driven
gas, the process of energy transfer to CDM particles can also be envisaged as
arising from potential fluctuations. connected to the density fluctuations of
compressible turbulence. Repeated inflows and outflows can then be
incorporated as part of a complete spectrum of perturbations associated with
stochastic fluctuations heating the dark halo (El-Zant et al., 2016; Hashim et
al., 2023). Full hydrodynamic simulations have indeed shown that potential
fluctuations from feedback driven gas can indeed lead to core formation in the
central region of CDM haloes, at least in the case of dwarf galaxies (Gelato &
Sommer-Larsen, 1999; Read & Gilmore, 2005; Governato et al., 2010; Di Cintio
et al., 2014; Teyssier et al., 2013), but the efficacy of the process depends
on the feedback recipe, and in particular requires a high star formation
threshold to maintain the sustained fluctuations leading to the dynamical
heating of dark halo centres (e.g., Pontzen & Governato, 2014; Benítez-Llambay
et al., 2019; Dutton et al., 2019, 2020).

Another proposed mechanism, which originally inspired the repeated feedback-
driven bursts models (Mashchenko et al., 2006), and is indeed closely related
in terms of dynamics (Hashim et al., 2023), invokes dynamical friction-
mediated energy transfer to the dark halo particles from massive baryonic
clumps during the galaxy formation process (El-Zant et al., 2001, 2004a). In
the standard CDM-based galaxy formation, gas condenses and cools inside dark
matter haloes (White & Rees, 1978). If the gas distribution remains smooth
during the process, adiabatic contraction ensues (Blumenthal et al., 1986;
Gnedin et al., 2004), further exacerbating the problem of central halo
overdensity. This effect can nevertheless be overcome if the baryonic medium
is clumpy (Del Popolo, 2009), and the mechanism has been repeatedly shown to
be effective, in principle, in reducing the central density of cuspy haloes
into cores (e.g., Romano-Díaz et al., 2008; Cole et al., 2011; Del Popolo,
2012; Saburova & Del Popolo, 2014; Inoue & Saitoh, 2011; Nipoti & Binney,
2015).

The dynamical friction coupling mechanism may especially be called for in
massive galaxies, as supernova feedback is generally insufficient for
producing cores in larger galaxy haloes; indeed for massive galaxies dark
matter contraction due to baryonic settling wins out against supernova
feedback driven core formation, even in simulations that produce cores in
small galaxies (Tollet et al., 2016. Even if AGN feedback can partially
mitigate the effect; Macciò et al., 2020; Forouhar Moreno et al., 2022). On
the other hand, non-singular dark matter halo cores are also inferred in
massive low brightness galaxies (de Blok et al., 2001, 2003) and may possibly
also exist in massive luminous ones (Salucci, 2019). The dynamics of such
galaxies seem especially inconsistent with observations if the sole effects of
baryons on an NFW halo is modelled is conceived in terms of adiabatic
contraction (Li et al., 2022). Although, in general, both feedback and
dynamical friction induced halo heating may act in concert during galaxy
formation (Orkney et al., 2021), dynamical friction mediated halo heating may
be especially relevant in way of explaining these results, as well as those of
apparently cored haloes of higher redshift massive galaxies (Dekel et al.,
2021; Ogiya & Nagai, 2022).

When that mechanism was first proposed there were no observations of baryonic
clumps of the required mass. For example, for a Milky Way like galaxy, the
required clumps were predicted to have masses $\gtrsim 10^{8}M_{\odot}$, while
the largest molecular clouds observed in local galaxies were an order of
magnitude smaller in mass (and also not believed to be sufficiently long
lived). There were observations of irregularly shaped and clustered early
galaxies, but these were primarily interpreted in terms of merging and
interacting systems of similar scale rather than collections of compact clumps
embedded in larger hosts (van den Bergh et al., 1996). However, since then, an
increasing number of galaxies, termed ’clump clusters’ or ’chain galaxies’,
have been discovered to host massive and compact clumpy structures that are
more abundant at higher redshifts, and which may form in situ (Elmegreen et
al., 2004, 2007, 2009; Genzel et al., 2011; Förster Schreiber et al., 2011;
Guo et al., 2015, 2018; Messa et al., 2022; Martin et al., 2023; Sattari et
al., 2023). The existence of the clumps may be related to disc formation and
can potentially have important consequences regarding the the structure of the
host galaxy; when gas is infalling towards a forming disc, radiative cooling
induces self-gravitating instabilities that lead to the formation of massive
clumps, which may in turn migrate — due to dynamical friction, as well as to
torques arising from scale perturbations in the disc — to the centres of their
host galaxies (Noguchi, 1999; Agertz et al., 2009; Ceverino et al., 2010;
Mandelker et al., 2014, 2017; Dekel et al., 2022), possibly removing dark
matter from the vicinity of a nascent bulge in the process (Elmegreen et al.,
2008). It has also been suggested that instabilities in purely stellar discs
may lead to clumpy structures that are sufficiently long lived to likewise
form a central bulge (Saha & Cortesi, 2018). Though there are varying views on
the matter (cf. Section 2.2.3 of Del Popolo & Pace, 2016 for a discussion),
some observational studies are consistent with the expectations of an
evolutionary scenario invoking sufficiently long lived (> 100 Myr) clumps
forming in situ and migrating to the centres of galaxies (Genzel et al., 2008,
2011; Shibuya et al., 2016; Soto et al., 2017; Zanella et al., 2019).

Nevertheless, the lifetime of high redshift clumpsis still uncertain; they
may, in theory, indeed disrupt over timespans much shorter than the of inward
migration timescale (Murray et al., 2010; Oklopčić et al., 2017). The
fragmentation mass scale found in simulations may furthermore be subject to
numerical resolution effects (Tamburello et al., 2015). Resolution may also
affect the observationally inferred effective radii and stellr masses of the
clumps (e.g., Dessauges-Zavadsky et al., 2017; Fisher et al., 2017a; Cava et
al., 2018), revising them downwards. Selection effects related to the limited
rest frame wavelength range available to the HST may further affect the
inferred masses. Thus, although early studies using the JWST seem consistent
with clumps surviving for up to a Gyr and centrally migrating (Claeyssens et
al., 2023; Kalita et al., 2023), it is still not obvious that observed higher
redshift clumps have the right properties (mass, size and robustness) for
dynamical friction coupling with the host system to be sufficient for bulge
building or halo core formation.

On the other hand, local anologues of early clumpy galaxies, albeit rare, do
exist and have been studied (e.g. Fisher et al., 2017a, b; Messa et al.,
2019). The present work is an attempt to make use of such observations to pin
down and isolate the potential importance and consequences of dynamical
friction mediated coupling in clumpy galaxies, To this end, we undertake
N-body simulations of stellar clumps embedded in a disc-halo system, with
clump masses and sizes consistent with those presented the relatively well
resolved DYNAMO-HST sample (Ambachew et al., 2022). If sufficiently long
lived, such observed clumps would be in the process of migrating towards the
centres of their host galaxies, a process that may lead to the simultaneous
formation of a central bulge (from the remnants of the clumps), in addition to
a halo core replacing the cusp. Thus, even galaxies that are initially dark
matter dominated at all radii may end up with much less dark matter in their
central regions; as the dark matter is replaced by a nascent baryonic bulge.
Such a scenario appears particularly appealing as it may contribute to the
’conspiring’ of dark halo and baryonic components to maintain flat rotation
curves, a phenomenon sometimes termed as the disk halo or bulge halo
conspiracy. And, as our results here also suggest, to the wide diversity of
observed galaxy rotation curves.

Whether this scenario actually materializes in galaxies will depend on the
timescale associated with the clump migration relative to the lifetime of the
clumps. The latter timescale in turn entails two processes: disruption of the
clumps as a result of further star formation, and tidal stripping due to the
motion of clumps in the disk-halo potential of the host galaxy. The internal
composition and density profile of the DYNAMO-HST clumps, for which only
stellar masses are available, is still uncertain. However, the work of Lenkić
et al. (2021) suggests that the clumps have a complex substructure of stellar
populations with different ages, with the centres of the clumps being
populated with young stars while older ones occupy the outer parts.
Accordingly, and given their age distribution as a function of distance from
their galactic centres, the DYNAMO-HST clumps were interpreted as relatively
long-lived structures; in the sense of being robust against sudden, internal
starburst-induced disruption. We thus focus here on the second effect, namely
that of disruption due to tidal stripping.

We thus conduct numerical simulations of disc galaxies containing clumps of
similar masses and sizes, and radial distribution in the host system, to those
found in the DYNAMO systems. As the internal density profiles of the DYNAMO
clumps is uncertain, we try three different theoretical ones (Plummer,
Hernquist and Jaffe profiles), to test under what conditions the clumps may
survive long enough to centrally migrate to form a stellar bulge and induce a
dark matter core. We also examine the variety of rotation curve shapes that
arise from the process, and how it depends on host halo mass and
concentration. These results are described in Section 4. In the next section
we first outline the properties of the simulated clumps, and how they are
related to the observations, while in Section 3 we describe the simulation
setup, including the effective initial conditions for our disc-halo-clump
system. We summarise and discuss our results in Section 5.

##  2 The DYNAMO-HST stellar clumps

To obtain observationally constrained models of stellar clumps embedded in
galactic discs, we use the recent DYNAMO-HST observations of the stellar
masses of baryonic clumps, as presented in Ambachew et al. (2022). These
involve six galaxies from the DYNAMO sample (Green et al., 2014), which are
selected from Sloan Digital Sky Survey (SDSS) DR4 (Adelman-McCarthy et al.,
2006). The galaxies are local ones, with redshift range $z=0.075-0.13$,
exhibiting a turbulent gas rich bulge-less disk with high star formation rate.
The inferred stellar masses of the galaxies range from $1.7\;{\rm
to}\;6.4\times 10^{10}M_{\odot}$, the upper limit being close to that of the
Milky Way (Licquia & Newman, 2015).

###  2.1 Clump sizes and masses

In Ambachew et al. (2022), the size of the identified clumps are measured by
fitting an elliptical Gaussian function to the 2D brightness distributions of
the clumps around the peak. The effective radius of each clump is then defined
as the mean standard deviation of the major and minor axis of the best fit 2-D
Gaussian functions. This is the region inside which a young stellar population
of the clumps is exceedingly bright compared to the background discs they are
embedded in. The stellar masses of the clumps are then inferred by fitting
stellar population synthesis models to the observed spectral energy
distributions, distinguishing between ’raw’ and ’disc subtracted’ clump
masses.

As the clump sizes are of order of the disc thickness, and as we will assume
the clumps to be self-gravitating and dynamically distinct objects, we
consider only the raw masses, which we refer to here as $M_{\rm obs}$. These
range from $4.71\times 10^{6}$ to $3.74\times 10^{9}M_{\odot}$, with average
mass $2.53\times 10^{8}M_{\odot}$. We note that they are about an order of
magnitude smaller than stellar clump masses quoted in many higher redshift
studies, an effect that may arise from the lower resolution available in these
studies (Fisher et al., 2017a) The observed clump mass-size relation may be
fitted with a power law, $\rm M_{\rm obs}=r^{n}_{\rm obs}$, with $n=2.4$, as
illustrated in Fig. 1.

As mentioned, the masses of the clumps are observationally evaluated by
fitting their light curves with Gaussians and imposing a cutoff at one
standard deviation. We will want to construct physical self gravitating mass
models for the observed clumps. As a 2-D standard deviation (with same
deviation in each dimension) corresponds to a fraction of 0.39 times the light
comprised under the full 2-D Gaussian curve, we will generally assume that
this cutoff — and the associated observationally quoted effective radius for
each clump — corresponds to the half mass and half mass radius of the clumps.
For which we will adopt various mass models to describe (Section 3.2).

###  2.2 Clumps mass function

Figure 1: The DYNAMO-HST clumps mass-size relation (from the data of Ambachew
et al., 2022). The solid line is the power law fit to the mass-size scaling.
The horizontal black-dashed line represents the observational mass limit of
$10^{7}M_{\odot}$. The masses are inferred by fitting the radial light
distribution by 2-D Gaussians, with cuts on orthogonal axes corresponding to
one standard deviation (which corresponds to a fraction of about 0.39 of the
region under such a curve, assuming equal standard deviations on both axes).
To construct physical self gravitating mass models, we associate the effective
radius, thus inferred, with the half mass radius of clumps with total masses
$M_{c}=2M_{\rm obs}=2M_{1/2}$, and with internal mass distributions modelled
through three different physical density profiles (Section 3.2).

According to the data given in Ambachew et al. (2022), the average number of
clumps per galaxy in the DYNAMO sample is about 13, which is the number we
adopt in each of the simulations described below. The number of clumps as a
function of mass may be fit by the following mass function (Ambachew et al.,
2022):

$\phi(x)=\ln(10)\times\left(\frac{M_{\rm obs}}{M_{0}}\right)^{\alpha+1.0},$
(1)

where $x=\log(M_{\rm obs})$, $\alpha=-1.4$ and $M_{0}=10^{9.5}M_{\odot}$. As a
realisation of this mass faction, we take four equal logarithmic mass bins in
which the clumps are distributed. In accordance to the discussion above, the
total mass of each clump $M_{c}$ is assumed to correspond to twice the masses
$M_{\rm obs}$, which are taken to correspond to the half mass of the clumps.
In practice, in order to arrive to a quasi-equilibrium of the simulated halo-
disc plus clump system, while approximately matching the masses $M_{c}$ and
radial location $R$ of the clumps in the observed sample, we start the our
simulated clumps at large $R$ (relative to the observed ones), with masses
$2M_{c}$, and let the system relax, with clumps migrating inwards through
dynamical friction and losing mass due to stripping, for about a rotation
period. This procedure, which sets our effective initial conditions, is
described in Section 3.3.

##  3 Simulation Setup

We advance the dynamics of the disc, halo and clump system, with properties
and parameters described below, with the publicly available GADGET-4 (Springel
et al., 2021) code for N-body simulation of isolated systems. The code employs
fully adaptive time-stepping and uses a hierarchical tree algorithm, in
combination with a particle-mesh scheme, to adaptively compute the
gravitational forces.

The simulations involve two types of particles; dark matter and stellar,
composing the halo and the disc with its embedded clumps, respectively. The
gravitational softening of both species is set equal to $\epsilon=0.035\;\rm
kpc$. The maximum simulation time is set to $\rm 3\;Gyr$. In the runs
described below, the total number of particles composing the halo is
$N=256^{3}$, with mass resolution of $m_{p}=5.9\times 10^{4}M_{\odot}$.
Convergence and initial equilibrium tests are presented in Appendix A

In the next subsections we describe the structure of the components of our
disc-halo-clump system. Then we detail our starting initial conditions, and
the relaxation towards a quasi-equilibrium, with clump positions and mass
distribution consistent with the DYNAMO observations, which will represent the
effective initial conditions for our runs.

###  3.1 The halo and the disc

Figure 2: Available Observed rotation curves of DYNAMO-HST galaxies used in
this work (Sweet et al., 2019), versus simulated ones, taken at the effective
start of the simulations (at $T=0$, after a 500 Myr relaxation period as
described in Section 3.3). The theoretical lines show models with halo
concentrations $c=20.4$ and $c=10.2$.

The N-Body simulations presented here are of clumps embedded in an isolated
galaxy composed of:

  1. 1.

A dark matter halo following a Navarro,Frenk and White (NFW) density profile
(Navarro et al., 1997),

$\rho=\frac{\rho_{s}}{\frac{r}{r_{s}}\left(1+\frac{r}{r_{s}}\right)^{2}},$ (2)

where $\rho_{s}$ and $r_{s}$ are density and radial scales.

In our fiducial model the total host (disc + halo) mass is $10^{12}M_{\odot}$,
so that the halo mass is $(10^{12}-f_{d})M_{\odot}$, where $f_{d}$ is the disc
mass fraction (taken as 0.02 in our fiducial model). We assume a
flat-$\rm\Lambda CDM$ cosmology (with Hubble constant $H_{0}=70.0$ and matter
density ($\Omega_{m}=0.3$), so that the halo viral radius is $r_{\rm
vir}\equiv cr_{s}=204.0~{}{\rm kpc}$, where $c=20.4$ is the halo
concentration. This is larger than the average found in simulations at this
mass scale (which is closer to $c=10$ at low redshift; e.g., Prada et al.,
2012), but it ensures that the model galaxy is initially entirely dark matter
dominated at all radii, including the central region. One of our goals will be
to probe the transformation from such a situation into one whereby the centre
becomes baryon dominated, which occurs concurrently with the formation of the
halo core. We also run models with halo concentration parameters $c=10.2$ for
comparison. In all cases, to avoid a sharp cutoff, our halos are simulated up
to radius larger than the assumed virial radius (namely up to $r=230~{}{\rm
kpc}$). In the absence of the clumps, halo profiles have been checked to,
remain unchanged (after reaching equilibrium state) during the simulation
time, as shown in Fig. 18.

  2. 2.

A stellar disc, with density profile given by:

$\rho_{d}(r)=\rho_{d}\exp\left(\frac{-l}{l_{d}}\right)\operatorname{sech}^{2}\left(\frac{z}{z_{d}}\right),$
(3)

where $l=\sqrt{x^{2}+y^{2}}$, $l_{d}$ is the radial disc density scale and
$z_{d}$ is the vertical scale height. To ensure stability against significant
non-axisymmetric instabilities, such as strong bars, we choose a large value
for the Toomre parameter, namely $Q=2.0$.

The values of the profile parameters of each component of the host system in
our fiducial system are listed in Table 1. For the radial extent of interest
to us (less than 10 kpc), our starting mass distributions, for $c=20.4$ and
$10.2$, approximately correspond to the range spanned by the three more
massive galaxies of the four galaxies in the DYNAMO samples used used here for
which rotation curves are available, as may be seen from from Fig. 2. We also
ran simulations with less massive halos and more massive discs (such as the
one corresponding to the rotation curves in Fig. 16).

| Halo | Disc
---|---|---
$M~{}[M_{\odot}]$ | $(1-f_{d})\times 10^{12}$ | $f_{*}\times 10^{12}$
$r_{s}~{}~{}{\rm[kpc]}$ | 10.0 | –
$c$ | 20.4 | –
$l_{d}~{}~{}{\rm[kpc]}$ | – | 3.5
$z_{d}~{}{\rm[kpc]}$ | – | 0.5
Table 1: The fiducial host (disc + halo) mass distribution parameters. The
total mass of the disc and halo combined is $\rm M_{\rm
host}=10^{12}M_{\odot}$, and the disc mass fraction is $f_{d}=0.02$. The
exponential disc has vertical and radial scale length $l_{z}$ and $l_{d}$, and
the NFW halo has concentration $c$ and radial scale length $r_{s}$. In
addition to those fiducial parameters we ran simulations with less
concentrated halo, with $c=10.2$, and also ones where the halo mass is halved
and disc mass doubled.

###  3.2 The internal profiles of simulated clumps

We add to the halo-disc system stellar clumps with properties constrained by
the DYNAMO-HST observations of Ambachew et al. (2022) (as described in
Section2). The resolution of the observations did not permit detailed
inference of the clump density profiles. Indeed, as mentioned, the observed
radial light profiles of the clumps were simply fit by a Gaussian and cut at
one-sigma, which we take to constrain the half mass.

A primary goal of the current study is to determine the robustness of the
stellar clumps against striping, which may determine the possibility of the
inward migration of the clumps, the building of bulges from their remnants,
and the formation of halo cores from dynamical friction mediated energy
transfer to dark matter particles. As these processes are expected to strongly
depend on the internal structure of the clumps, we try three different
profiles, differing in the steepness of their radial central density
distribution.

The three spherical density distributions we consider consist of: i) a Plummer
profile, characterised by a central nearly constant density core; ii) a
Hernquist profile, with central density rising as $1/r$; iii) a still more
cuspy Jaffe profile, with density going as $1/r^{2}$ in the central region.
The associated mass distributions are

$M(<r)=4\pi\rho_{c}r_{c}^{3}\begin{cases}\frac{\left(\frac{r}{r_{c}}\right)^{2}}{2\left(1+\frac{r}{r_{c}}\right)^{2}},&\text{Herquist}\\\
\frac{\frac{r}{r_{c}}}{1+\frac{r}{r_{c}}},&\text{Jaffe}\\\
\frac{1}{3}\frac{r^{3}}{2^{-\frac{3}{2}}\left(r^{2}+r^{2}_{c}\right)^{\frac{3}{2}}},&\text{Plummer}\end{cases}$
(4)

where $\rho_{c},r_{c}$ are density and radial scales characteristics of the
clumps. These are fixed through the observationally inferred half mass radius
and the tidal truncation radius described below.

###  3.3 Tidal truncation of the clumps

Figure 3: Internal clump profile sketch. In order to determine the scale
parameters of the profile, we set constraints on the clump core mass and
effective radius from DYNAMO-HST observations, and define the clump boundary
through the tidal (Jacobi) radius at the clump’s orbital position in the
host’s disc plane. Figure 4: The different mass profiles of a clump with mass
$M_{c}=10^{8.8}M_{\odot}$. The vertical black-dashed line is at the half mass
radius, which we take to be fixed by observations. This, along with the Jacobi
tidal truncation radius completely fixes the internal clump density
distribution, given the profile.

None of the profiles just described have finite truncation radii, as the
associated densities are finite at any finite radius from the centre. Clumps
embedded in our disc-halo system, however, will be tidally truncated. We use
the Jacobi radius to approximate the expected tidal truncation scale. For a
clump of mass $M_{c}$ and orbital radius $R$, this is defined as (Binney &
Tremaine, 2008):

$\rm r_{J}=\left(\frac{M_{c}}{3M_{\rm host}(<R)}\right)^{1/3}R,$ (5)

where $M_{\rm host}(<R)$ is the enclosed mass of the host (that is disk plus
halo mass) at radius $R$ in the disc plane. The clump profile parameters are
then determined from the condition that the masses of the clumps be completely
contained within $r_{J}$ and that half the mass be contained within the
observed radius $r_{\rm obs}=r_{1/2}$.

Fig. 3 shows a sketch illustrating the two radii that thus fix the internal
clump mass distribution, namely the half mass radius $r_{1/2}$ and the
truncation radius $r_{J}$. As a quantitative example of the resulting mass
profiles, Fig. 4 shows those corresponding to a clump of half size
$r_{1/2}=0.27\;{\rm kpc}$ and total mass $M_{c}=10^{8.8}M_{\odot}$, located at
an orbital radius $R=9.64~{}{\rm kpc}$ in the disc plane.

Figure 5: Setting the effective initial conditions. The initial (at $t=0$)
clump masses are distributed in logarithmic mass bins in accordance with
equation (1), as described in Section 2. The masses are then doubled and the
clumps are placed at large radial positions $R$ in the disc (relative to those
observed), corresponding to the red crosses on the right of the plots. The
left panel shows the positions of the clumps after the (disc-halo-clumps)
system is left to relax by advancing the combined dynamics for a relaxation
period of 500 Myr (corresponding to $T=0$), with different symbols denoting
runs with the three different internal clump profiles used. The grey round
dots correspond to the actual masses and positions of the clumps (in all
combined) DYNAMO galaxies. On the right hand panel, the lines represent the
estimated remnant mass, when the tidal stripping is estimated using eq. 5. The
vertical black dashed line is at $\rm R=1{\rm kpc}$, where the remnants of
larger clumps are expected to merge into a bulge after completing their
central migration. The scaling radius $R_{85\%}$ refers to the observed disc
radius comprising $85\%$ of the light, which we assume to correspond to the
same percentage in mass of the simulated discs.

###  3.4 Quasi-equilibrium start and relaxation to effective initial
conditions

####  3.4.1 Disc halo system

The quasi-equilibrium starting conditions for the disc-halo system are created
using the Disk Initial Conditions Environment (DICE) code (Perret et al.,
2014; Perret, 2016); the particle distribution is thus created using an N-try
MCMC algorithm, which does not require prior knowledge of the distribution
function. The dynamical equilibrium of each component in the system is
computed by integrating the Jeans equations for each disc and halo particle.

####  3.4.2 Clumps

Figure 6: Left: Simulated galaxy with clumps of Hernquist internal mass
profile after the initial relaxation period comprising $500{\rm Myr}$ of
evolution ($T=0$). Right: DYNAMO-HST galaxy D13-5 using F336W filter. Figure
7: Projected density maps of the stellar disc and clumps, for simulations of
the fiducial model (Table 1), with different clump profiles, as noted. The
clumps move to the centre of the disc-halo host due to dynamical friction
coupling. The resulting bulge size differs depending on the clumps’ internal
profile. The time $t=0.5~{}{\rm Gyr}$, corresponding to $T=0$, denotes the
effective starting conditions of our simulations, when the models relax to a
quasi-steady state with clump positions and masses in approximate
correspondence with the observed ones (cf. Fig 5).

Each clump centre of mass is assumed to initially move at the local circular
speed at its starting position in the disc-halo system. The combined system is
only in initial approximate quasi-equilibrium and must therefore be advanced
for a time of the order of a rotation period in order to settle. This is also
necessary for obtaining stable clump masses, as the initial internal dynamical
equilibrium of the embedded clumps, as well as their initial (tidally
truncated) radius can only be approximate, since the Jacobi relation Eq. (5)
for the the tidal radius is approximate. Indeed, as shown in Appendix B, it is
also a better approximation when the clumps are at larger $R$ than if used
directly to truncate the clumps at smaller radii within the disc.

We adapt the following procedure so as to arrive at a mass and positional
distribution of clumps compatible with that observed in the DYNAMO galaxies
after an initial ’relaxation’ period. We initially set the positions of the
clumps on the plane of the disc at relatively large radii, namely with
$8.0<R<16.0{\rm kpc}$, with a normal distribution for their radii and random
azimuthal angles, and imposing the condition that the Jacobi radius $r_{J}$ of
each clump is bigger than its half mass radius $r_{\rm obs}$ (which naturally
orders them in radii $R$ according to their masses).

As the clumps migrate inwards from their initially large radius (due to
dynamical friction), they also lose mass due to stripping. To match the
observed mass-position relation (Ambachew et al., 2022) of the clumps, as
depicted in Fig. 5, we double the initial mass and evolve the system for 500
Myr. This ensures that the observed distribution is approximately matched, and
also that the system reaches a quasi-steady state, which may be assumed to
correspond to the effective start of the simulation, We denote this time with
$T=0$.

We illustrate the procedure in Fig. 5. 111If anything, the applied procedure
somewhat underestimates the preponderance of massive clumps in the central
region, and in turn the associated effects due to dynamical friction coupling
that we describe below. Fig. 6 exhibits a pictorial comparison of the
effective initial conditions in the simulation with one of the observed
DYNAMO-HST galaxies. Finally, it is important to note that, as we will see
below, the physical changes in the system we are looking for, namely halo core
and baryonic bulge formation, are absent during the relaxation period. This
may be expected, as the relaxation process process involves the migration of
clumps through the outer region of the system. The effects we are interested
in, on the other hand, occur with the subsequent migration to the inner parts,
which is described below.

##  4 Results

Dynamical friction coupling between the massive stellar clumps and the host
system leads to clump migration toward the centre. This is illustrated in Fig.
7, where we show projections of the evolution of the disc-clump system. As the
clumps migrate towards the centre they lose energy to the background, composed
of the disc and the halo, which is thus dynamically ’heated’. As the sub-
system of halo particles is more centrally concentrated, it has larger binding
energy and therefore absorbs most of the energy lost by the clumps. This
results in a central halo core replacing the initial density cusp. The
remnants of the inwardly migrating clumps accumulate to form a central bulge.
The end result of these processes is the once dominant dark matter becoming
insignificant in terms of contribution to the central rotation curve, with the
main contribution there coming from the baryons by the end of the evolution.
We describe these processes in more detail below.

###  4.1 Energy transfer from clumps to cusp and the transformation to quasi-
isothermal core

Figure 8: Change in total energy of subsystems of halo particles enclosed
within radius $r$, as given by equation 6, as they gain energy from the
inwardly migrating stellar clumps. The results are illustrated for the
fiducial disc-halo model of Table 1, and for the three internal clump profiles
considered: Plummer spheres with constant density cores; Hernquist spheres,
with $1/r$ inner cusps; and Jaffe spheres, with $1/r^{2}$ inner cusps.

The energy gained by the inner halo as the clumps migrate inwards is
illustrated in Fig. 8, which shows the energy of systems made of subsets of
halo particles enclosed within spherical radii $r$. This is defined as

$E_{r}=\frac{m_{p}}{2}\left(\sum_{i}v_{i}^{2}+\Phi_{i}\right),$ (6)

where $m_{p}$ refers to the mass of a particle (recall that all particles have
the same mass in our simulations), $v_{i}$ is the speed of particle $i$ and
$\Phi_{i}$ is the Newtonian potential (due to all other particles in the
simulation), at its location. The summation is evaluated over all halo
particles within radius $r$.

One immediately sees that the energy input rate is larger for the more
concentrated internal clump mass distribution; namely for the Jaffe profile,
with its steeply increasing central density cusp. This is due to the fact that
as the clumps spiral to the centre they lose mass due to stripping, a process
that is least efficient for steeper internal clump profiles. As the energy
loss rate is proportional to the square of the mass, this is an important
effect. The competition between energy transfer and stripping in fact
determines the efficiency of the energy transfer from the stellar clumps to
the background particles, and associated the erasure of the halo cusp as we
will see below.

Figure 9: Evolution of the (three dimensional) velocity dispersion profile of
the host halo, for the same runs as in Fig. 8.

The halo central density cusp is particularly sensitive to the energy input
because it suffers from ’temperature inversion’, with inner velocity
dispersion decreasing towards the centre, which renders its structure
susceptible to even small energy input; as even modest gain in kinetic energy
is significant given its initially small value. Thus, as the halo particles
gain energy from the in-spiraling clumps, an isothermal core starts to replace
the temperature inversion connected to the initial cusp (Fig. 9). Again, the
effect is most prominent in the case with the most centrally concentrated
clumps (with internal Jaffe profiles).

Figure 10: Time evolution of Velocity dispersion when the halo concentration
is $c=10.2$, instead of $20.4$. The internal clump profile is of the Jaffe
form.

When a less concentrated halo profile is used (initial NFW $c=10.2$ instead of
20.4), halo particles are less bound, and the central halo structure is thus
more susceptible to energy input. This results in a larger effect in terms of
total energy change, as well as increase in velocity dispersion that leads to
the tendency towards an isothermal core (Fig. 10).

###  4.2 Density profiles

####  4.2.1 Cusp-core transformation

Figure 11: Halo density profile for the fiducial runs (corresponding to those
in Fig. 8).

The energy transfer through dynamical friction, and the accompanying advent of
a quasi-isothermal inner velocity dispersion profile, leads to the replacement
the initial halo density cusp by a core. In Fig. 11, we plot the halo density
profile in our fiducial model for the different internal clump profiles. After
a few hundred Myr of evolution from our initial relaxed quasi-steady state, a
nearly constant density core, of order of a kpc in scale, replaces the central
cusp. Given the considerations regarding the efficiency of the stripping as a
function of clump internal profile, the precise timescale and the core size
depends on that profile. The process being again more effective for more
steeply increasing clump central density.

Figure 12: Comparing the evolution of host halo density profiles when the
initial NFW concentration is $c=10.2$ (blue solid), to the $c=20.4$ case
(black solid), for internal clump densities set to the Hernquist and Jaffe
profiles.

As may be expected from the above discussion relating to the formation of the
isothermal core, the associated effect of density core formation is also
larger when the halo is less concentrated (and so its particles less bound).
This is confirmed in the plots shown in Fig. 12. (It is worth noting, however,
that because the initial halo central density is also smaller in this case,
the dynamical friction coupling is actually smaller, and so the effect is
initially slower than in the $c=20.4$. case.)

####  4.2.2 Total density profiles

Figure 13: Total (Halo+Disc+Clumps) density profile for halo concentrations
$c=20.4$ (black solid) and $c=10.2$ (blue solid), all other parameters
remaining fixed. The dashed lines on late time plots correspond to singular
isothermal ($1/r^{2}$) profiles.

The cusp-core transformation is associated with a decrease in halo density,
and hence in dark matter mass enclosed within a given radius in the region
where the cusp reigned. Nevertheless, despite the decrease in central halo
mass, the total mass in the central region actually increases, as the central
bulge, born of clump remnants, materialises. This result, which we illustrate
in Fig. 13, is in contrast to the case when clumps are distributed initially
in the same way as the smooth background (NFW) halo; in this case the total
density profile remains invariant to a good approximation (El-Zant et al.,
2004b; El-Zant, 2008, 2019). It is also opposite to the case when the cusp-
core transformation is due to fluctuating feedback driven gas; in that case
the total density decreases and the core forms due to the shallowing of the
central potential (Hashim et al., 2023; Appendix B).

When the effect of the dynamical friction is strongest — namely in the case of
the Jaffe internal clump profile — the total inner density distribution in
fact converges towards a logarithmic slope $\sim-2$, after sufficient time.
This is in line with observations of the central parts of early type galaxies,
where an inner ’universal’ profile is found, and leads to interesting
consequences regarding the rotation curves of spiral galaxies as well, as we
describe below. (These issues are also discussed further in the concluding
section).

###  4.3 Rotation curves: diversity and bulge halo conspiracies

Figure 14: Rotation curve of dark halo, stellar disc and clumps, as well as
the total, for all simulations with fiducial disc-halo parameters (Table 1),
shown at time $T=0.7~{}{\rm Gyr}$ (solid lines), and initially at $T=0.0$
(dashed lines).

Figure 15: Initial ($T=0$) and final ($T=0.7{\rm Gyr}$) rotation curves in
simulations with Hernquist (left) and Jaffe-profile clumps with halo
concentrations $c=10.2$. Figure 16: Initial ($T=0$) and final ($T=0.7{\rm
Gyr}$) rotation curves in simulation with Jaffe-profile clumps and halo
concentration $c=20.4$, but with halo mass reduced by half and stellar mass
(including both disc and clumps) doubled. Figure 17: Same as Fig. 7, with
Jaffe profile clumps, and when the halo mass is halved halo and stellar mass
(disc plus clumps) is doubled.

In Fig. 14, we plot the rotation curve (defined simply in terms of the mass
enclosed within a given radius) of the halo, disc and clumps, as well as the
total, for simulations of our fiducial disc-halo galaxy model with embedded
clumps endowed with the different internal profiles considered. Initially, for
the parameters chosen, the contribution of the halo is largest even towards
the centre of our simulated galaxy, which does not initially contain a bulge.
As a result of the clump migration, however, a central bulge forms; a process
that is accompanied by decrease in the central dark mass, as probed in
previous subsections. The net result is that, starting from dark matter
domination in the central region, one arrives at baryon domination in that
region (as is likely the case in bright spiral galaxies such as the Milky
Way). The evolved rotation curve now shows a central feature reflective of the
nascent bulge. Furthermore, as may be expected from the previous subsections,
the effect of decrease of central dark matter mass and increase in inner
baryonic mass is enhanced as the internal clump density profiles become
steeper, as one moves from the cored Plummer sphere to a singular isothermal
centre of a Jaffe model.

The effect of halo mass ejection from the central region, and replacement with
a baryonic bulge, is also more significant when a less concentrated initial
halo is used, as can be seen from Fig. 15. Here, the central baryonic and dark
matter contribution is comparable at the effective start of the simulation at
$T=0$, instead of the halo entirely dominating at all radii. But the baryon
domination is complete at the end of the runs. This also comes with the
emergence of ’bulge-halo conspiracy’, whereby the two components’
contributions approximately match at some transition radius to give a flat
rotation curves, despite one component being dominant in the central kpc or
so, and the other contributing most outside that radius.

We also ran simulations with parameters that further clarify the possible
origin of such a bulge halo conspiracy, including the advent of a central bump
of the type often observed in galaxies. For this purpose, we start from a
stellar component (disk plus clumps) that is doubled in mass (in both the disc
and clump masses) and a halo mass that is halved (while fixing $c=20.4$). The
resulting rotation curves are shown in Fig. 16, for Jaffe profile clumps. The
final configuration — with central bump linking the transition from baryonic
to dark matter domination, resulting in a nearly flat rotation curve despite
the transition — is akin to that of many observed galaxies, including that of
the Milky Way (e.g. Sofue, 2017). The corresponding projected density plots,
shown in Fig. 17 also show a thickening of the disc during the evolution,
resulting from the interaction with the clumps, which finally spiral in and
merge to form a prominent central bulge.

##  5 Conclusion

We performed simulations of model disc galaxies with massive stellar clumps
embedded in the discs. The effective initial conditions were designed so that
the initial masses and sizes of the clumps are compatible with those observed
in the HST Dynamo sample of local galaxies.

The clumps are coupled through dynamical friction to the host disc-halo
system, and so spiral to the centre (within a few hundred Myr to $\la{\rm
Gyr}$), forming a central bulge. In the process, the central halo is
dynamically heated; with the ’temperature inversion’, characteristic of the
initial density cusp, giving way to an isothermal core. The combination of
centre-bound clump migration and outward flow of dark matter leads to the
baryonic component dominating in the inner region, even when starting with a
model galaxy that is dark matter dominated at all radii.

The efficiency of the aforementioned processes depend on the rate at which
clumps are stripped as they move the disc-halo potential. which in turn
depends on the form of their internal density profiles. We tried three
different forms. In order of less to more steep inner density radial
variation, these were the Plummer profile (with inner core), the Hernquist
profile (with $1/r$ cusp), and the Jaffe profile (with $1/r^{2}$ cusp). We
find that the steeper the inner clump profile — and so, the more mass in the
inner region, given the same total mass and tidal radius — the more effective
the process of clump migration and concurrent effect of halo cusp-core
transformation.

The replacement of the inner dark matter mass by the baryonic component of
clump remnants thus effectively depends on the competition between two
effects, namely dynamical friction and stripping. Given that the efficiency of
the stripping also depends on the host halo mass and concentration, being more
efficacious when these are increased, more prominent bulges and larger halo
cores are thus naturally produced if the halo concentration is decreased, or
when the halo mass is decreased relative to the disc-clump system (in the
latter case the effect being particularly prominent due to larger clump
masses).

In this context, the variation of the inner clump profiles, as well as the
halo concentration and mass of the halo relative to the stellar component,
lead to the advent of a diversity in final rotation curves; from minor inner
hump, reflective of the transformation from dark matter to baryon dominance in
the very inner region, to the building of a more prominent bulge (Fig. 14,
15). When the effect is strongest, this bulge ’conspires’ with the nascent
cored halo component to produce produce a flat rotation curve, complete with
the characteristic inner bump heralding the transition from baryon to dark
matter dominance (16).

Such conspiracies, as reflected in the rotation curves of disc galaxies, were
once thought to characterise most such galaxies. But already since the 1990’s
the consensus shifted towards the realisation that there was a diversity of
rotation curves, and that the conspiracy was only present to high surface
brightness galaxies (Ashman, 1992). More recent research has been concerned
with the bulge halo conspiracy in early type galaxies, particularly with the
seemingly universal isothermal total (baryonic plus dark) density profile
characterising of such systems (e.g. Koopmans et al., 2006; Gavazzi et al.,
2007; Dutton & Treu, 2014; Chae et al., 2014; Shankar et al., 2017).

This isothermal profile, which also seems to characterise the total mass
distribution of bright disc galaxies over a large radial span (Zhengra et al.,
2019), is reproduced as the final state in our simulations; particularly those
with parameters chosen such that the dynamical friction coupling is most
efficient and stripping least effective (Fig. 13). On the other hand, the
variety of end states we find, depending on the relative importance of these
two effects, may provide further insight, which may be added to existing
explanations and interpretations of the diversity problem of spiral galaxy
rotation curves (e.g., Zentner et al., 2022).

Thus, to sum up. The migration of stellar clumps leads to the formation of
halo cores that replace the initial inner dark matter density cusps with a
baryonic central bulge. Depending on the inner profile of the clumps, the halo
concentration, and its mass relative to the stellar component, a variety of
final mass distributions may arise, including those associated with a singular
isothermal total mass distribution characteristic of massive disc and early
type systems. The variety of possible end states may also form a contributing
factor to the issue of diversity of rotation curves in spiral galaxies and
total mass distribution of early type galaxies. Indeed, in the present
scenario, the steepening of the inner of the inner profile may be expected to
correlate with surface brightness (e.g., with a massive stellar component and
less concentrated halo) as observed, a contention that may be tested in detail
with future simulations.

In concluding, it is important to note that our simulations do not include
feedback, from either statburst or AGN, which may reduce the inner baryonic
dominance that is the final result of our simulations, and perhaps even lead
to total reduction in central mass distribution, contributing further to the
diversity of observed rotation curves. Another related issue of course
pertains to the survival of the observed stellar clumps, at least on
timescales of the order of disc rotation period at solar radius in a Milky Way
like galaxy. In the present study we have assumed that they may be disrupted
through stripping but not through internal starburst in a major gaseous
component. This is suggested by the observations in Ambachew et al. (2022),
which suggest that the clump sample on which this study is based is relatively
robust against such effects. Further progress in examining these issues, as
well as more information regarding the internal clump profile, which is our
main unknown here, may be provided by future observations.

## Acknowledgements

We would like to thank Sarah Sweet for providing rotation curve data of
DYNAMO-HST galaxies used in this paper.

## Data Availability

Data is available upon request.

## References

  * Adelman-McCarthy et al. (2006) Adelman-McCarthy J. K., et al., 2006, ApJS, 162, 38
  * Agertz et al. (2009) Agertz O., Teyssier R., Moore B., 2009, MNRAS, 397, L64
  * Ambachew et al. (2022) Ambachew L., et al., 2022, MNRAS, 512, 3079
  * Ashman (1992) Ashman K. M., 1992, PASP, 104, 1109
  * Benítez-Llambay et al. (2019) Benítez-Llambay A., Frenk C. S., Ludlow A. D., Navarro J. F., 2019, MNRAS, 488, 2387
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press
  * Blumenthal et al. (1986) Blumenthal G. R., Faber S. M., Flores R., Primack J. R., 1986, ApJ, 301, 27
  * Bode et al. (2001) Bode P., Ostriker J. P., Turok N., 2001, ApJ, 556, 93
  * Boylan-Kolchin et al. (2011) Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2011, MNRAS, 415, L40
  * Buchdahl (1970) Buchdahl H. A., 1970, MNRAS, 150, 1
  * Bullock & Boylan-Kolchin (2017) Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343
  * Burkert (1995) Burkert A., 1995, ApJ, 447, L25
  * Cava et al. (2018) Cava A., Schaerer D., Richard J., Pérez-González P. G., Dessauges-Zavadsky M., Mayer L., Tamburello V., 2018, Nature Astronomy, 2, 76
  * Ceverino et al. (2010) Ceverino D., Dekel A., Bournaud F., 2010, MNRAS, 404, 2151
  * Chae et al. (2014) Chae K.-H., Bernardi M., Kravtsov A. V., 2014, MNRAS, 437, 3670
  * Claeyssens et al. (2023) Claeyssens A., Adamo A., Richard J., Mahler G., Messa M., Dessauges-Zavadsky M., 2023, MNRAS, 520, 2180
  * Cole et al. (2011) Cole D. R., Dehnen W., Wilkinson M. I., 2011, MNRAS, 416, 1118
  * Colín et al. (2000) Colín P., Avila-Reese V., Valenzuela O., 2000, ApJ, 542, 622
  * Dekel et al. (2021) Dekel A., et al., 2021, MNRAS, 508, 999
  * Dekel et al. (2022) Dekel A., Mandelker N., Bournaud F., Ceverino D., Guo Y., Primack J., 2022, MNRAS, 511, 316
  * Del Popolo (2009) Del Popolo A., 2009, ApJ, 698, 2093
  * Del Popolo (2012) Del Popolo A., 2012, MNRAS, 419, 971
  * Del Popolo & Hiotelis (2014) Del Popolo A., Hiotelis N., 2014, J. Cosmology Astropart. Phys., 1, 47
  * Del Popolo & Kroupa (2009) Del Popolo A., Kroupa P., 2009, A&A, 502, 733
  * Del Popolo & Le Delliou (2017) Del Popolo A., Le Delliou M., 2017, Galaxies, 5, 17
  * Del Popolo & Pace (2016) Del Popolo A., Pace F., 2016, Ap&SS, 361, 162
  * Dessauges-Zavadsky et al. (2017) Dessauges-Zavadsky M., Schaerer D., Cava A., Mayer L., Tamburello V., 2017, ApJ, 836, L22
  * Di Cintio et al. (2014) Di Cintio A., Brook C. B., Dutton A. A., Macciò A. V., Stinson G. S., Knebe A., 2014, MNRAS, 441, 2986
  * Dutton & Treu (2014) Dutton A. A., Treu T., 2014, MNRAS, 438, 3594
  * Dutton et al. (2019) Dutton A. A., Macciò A. V., Buck T., Dixon K. L., Blank M., Obreja A., 2019, MNRAS, 486, 655
  * Dutton et al. (2020) Dutton A. A., Buck T., Macciò A. V., Dixon K. L., Blank M., Obreja A., 2020, MNRAS, 499, 2648
  * El-Zant (2008) El-Zant A. A., 2008, ApJ, 681, 1058
  * El-Zant (2019) El-Zant A. A., 2019, in Journal of Physics Conference Series. p. 012007, doi:10.1088/1742-6596/1253/1/012007
  * El-Zant et al. (2001) El-Zant A., Shlosman I., Hoffman Y., 2001, ApJ, 560, 636
  * El-Zant et al. (2004a) El-Zant A. A., Hoffman Y., Primack J., Combes F., Shlosman I., 2004a, ApJ, 607, L75
  * El-Zant et al. (2004b) El-Zant A. A., Hoffman Y., Primack J., Combes F., Shlosman I., 2004b, ApJ, 607, L75
  * El-Zant et al. (2016) El-Zant A. A., Freundlich J., Combes F., 2016, MNRAS, 461, 1745
  * Elmegreen et al. (2004) Elmegreen D. M., Elmegreen B. G., Hirst A. C., 2004, ApJ, 604, L21
  * Elmegreen et al. (2007) Elmegreen D. M., Elmegreen B. G., Ravindranath S., Coe D. A., 2007, ApJ, 658, 763
  * Elmegreen et al. (2008) Elmegreen B. G., Bournaud F., Elmegreen D. M., 2008, ApJ, 688, 67
  * Elmegreen et al. (2009) Elmegreen D. M., Elmegreen B. G., Marcus M. T., Shahinyan K., Yau A., Petersen M., 2009, ApJ, 701, 306
  * Famaey et al. (2020) Famaey B., Khoury J., Penco R., Sharma A., 2020, J. Cosmology Astropart. Phys., 2020, 025
  * Ferraro (2012) Ferraro R., 2012, doi:10.1063/1.4756821, 1471, 103
  * Fisher et al. (2017a) Fisher D. B., et al., 2017a, MNRAS, 464, 491
  * Fisher et al. (2017b) Fisher D. B., et al., 2017b, ApJ, 839, L5
  * Flores & Primack (1994) Flores R. A., Primack J. R., 1994, ApJ, 427, L1
  * Forouhar Moreno et al. (2022) Forouhar Moreno V. J., Benítez-Llambay A., Cole S., Frenk C., 2022, MNRAS, 511, 3910
  * Förster Schreiber et al. (2011) Förster Schreiber N. M., et al., 2011, ApJ, 739, 45
  * Frenk & White (2012) Frenk C. S., White S. D. M., 2012, Annalen der Physik, 524, 507
  * Freundlich et al. (2020) Freundlich J., Dekel A., Jiang F., Ishai G., Cornuault N., Lapiner S., Dutton A. A., Macciò A. V., 2020, MNRAS, 491, 4523
  * Garrison-Kimmel et al. (2013) Garrison-Kimmel S., Rocha M., Boylan-Kolchin M., Bullock J. S., Lally J., 2013, MNRAS, 433, 3539
  * Garrison-Kimmel et al. (2014) Garrison-Kimmel S., Boylan-Kolchin M., Bullock J. S., Kirby E. N., 2014, MNRAS, 444, 222
  * Gavazzi et al. (2007) Gavazzi R., Treu T., Rhodes J. D., Koopmans L. V. E., Bolton A. S., Burles S., Massey R. J., Moustakas L. A., 2007, The Astrophysical Journal, 667, 176
  * Gelato & Sommer-Larsen (1999) Gelato S., Sommer-Larsen J., 1999, MNRAS, 303, 321
  * Genzel et al. (2008) Genzel R., et al., 2008, ApJ, 687, 59
  * Genzel et al. (2011) Genzel R., et al., 2011, ApJ, 733, 101
  * Gnedin et al. (2004) Gnedin O. Y., Kravtsov A. V., Klypin A. A., Nagai D., 2004, ApJ, 616, 16
  * Goodman (2000) Goodman J., 2000, New Astron., 5, 103
  * Governato et al. (2010) Governato F., et al., 2010, Nature, 463, 203
  * Green et al. (2014) Green A. W., et al., 2014, MNRAS, 437, 1070
  * Guo et al. (2015) Guo Y., et al., 2015, ApJ, 800, 39
  * Guo et al. (2018) Guo Y., et al., 2018, ApJ, 853, 108
  * Hashim et al. (2023) Hashim M., El-Zant A. A., Freundlich J., Read J. I., Combes F., 2023, MNRAS, 521, 772
  * Hu et al. (2000) Hu W., Barkana R., Gruzinov A., 2000, Phys. Rev. Lett., 85, 1158
  * Hui et al. (2017) Hui L., Ostriker J. P., Tremaine S., Witten E., 2017, Phys. Rev. D, 95, 043541
  * Inoue & Saitoh (2011) Inoue S., Saitoh T. R., 2011, MNRAS, 418, 2527
  * Kalita et al. (2023) Kalita B. S., Silverman J. D., Daddi E., Bottrell C., Ho L. C., Ding X., Yang L., 2023, arXiv e-prints, p. arXiv:2309.05737
  * Kamionkowski & Liddle (2000) Kamionkowski M., Liddle A. R., 2000, Phys. Rev. Lett., 84, 4525
  * Klypin et al. (1999) Klypin A., Kravtsov A. V., Valenzuela O., Prada F., 1999, ApJ, 522, 82
  * Koopmans et al. (2006) Koopmans L. V. E., Treu T., Bolton A. S., Burles S., Moustakas L. A., 2006, The Astrophysical Journal, 649, 599
  * Kravtsov & Wu (2023) Kravtsov A., Wu Z., 2023, arXiv e-prints, p. arXiv:2306.08674
  * Lenkić et al. (2021) Lenkić L., Bolatto A. D., Fisher D. B., Glazebrook K., Obreschkow D., Abraham R., Ambachew L., 2021, MNRAS, 506, 3916
  * Li et al. (2022) Li P., McGaugh S. S., Lelli F., Schombert J. M., Pawlowski M. S., 2022, A&A, 665, A143
  * Li et al. (2023) Li Z., Dekel A., Mandelker N., Freundlich J., François T. L., 2023, MNRAS, 518, 5356
  * Licquia & Newman (2015) Licquia T. C., Newman J. A., 2015, ApJ, 806, 96
  * Macciò et al. (2020) Macciò A. V., Crespi S., Blank M., Kang X., 2020, MNRAS, 495, L46
  * Mandelker et al. (2014) Mandelker N., Dekel A., Ceverino D., Tweed D., Moody C. E., Primack J., 2014, MNRAS, 443, 3675
  * Mandelker et al. (2017) Mandelker N., Dekel A., Ceverino D., DeGraf C., Guo Y., Primack J., 2017, MNRAS, 464, 635
  * Martin et al. (2023) Martin A., et al., 2023, ApJ, 955, 106
  * Mashchenko et al. (2006) Mashchenko S., Couchman H. M. P., Wadsley J., 2006, Nature, 442, 539
  * Messa et al. (2019) Messa M., Adamo A., Östlin G., Melinder J., Hayes M., Bridge J. S., Cannon J., 2019, MNRAS, 487, 4238
  * Messa et al. (2022) Messa M., Dessauges-Zavadsky M., Richard J., Adamo A., Nagy D., Combes F., Mayer L., Ebeling H., 2022, MNRAS, 516, 2420
  * Milgrom (1983a) Milgrom M., 1983a, ApJ, 270, 365
  * Milgrom (1983b) Milgrom M., 1983b, ApJ, 270, 371
  * Moore (1994) Moore B., 1994, Nature, 370, 629
  * Moore et al. (1999) Moore B., Quinn T., Governato F., Stadel J., Lake G., 1999, MNRAS, 310, 1147
  * Murray et al. (2010) Murray N., Quataert E., Thompson T. A., 2010, ApJ, 709, 191
  * Navarro et al. (1996a) Navarro J. F., Eke V. R., Frenk C. S., 1996a, MNRAS, 283, L72
  * Navarro et al. (1996b) Navarro J. F., Frenk C. S., White S. D. M., 1996b, ApJ, 462, 563
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Navarro et al. (2010) Navarro J. F., et al., 2010, MNRAS, 402, 21
  * Nipoti & Binney (2015) Nipoti C., Binney J., 2015, MNRAS, 446, 1820
  * Noguchi (1999) Noguchi M., 1999, ApJ, 514, 77
  * Ogiya & Burkert (2015) Ogiya G., Burkert A., 2015, MNRAS, 446, 2363
  * Ogiya & Nagai (2022) Ogiya G., Nagai D., 2022, MNRAS, 514, 555
  * Oklopčić et al. (2017) Oklopčić A., Hopkins P. F., Feldmann R., Kereš D., Faucher-Giguère C.-A., Murray N., 2017, MNRAS, 465, 952
  * Orkney et al. (2021) Orkney M. D. A., et al., 2021, MNRAS, 504, 3509
  * Peebles (2000) Peebles P. J. E., 2000, ApJ, 534, L127
  * Perret (2016) Perret V., 2016, DICE: Disk Initial Conditions Environment, Astrophysics Source Code Library, record ascl:1607.002 (ascl:1607.002)
  * Perret et al. (2014) Perret V., Renaud F., Epinat B., Amram P., Bournaud F., Contini T., Teyssier R., Lambert J. C., 2014, A&A, 562, A1
  * Pontzen & Governato (2012) Pontzen A., Governato F., 2012, MNRAS, 421, 3464
  * Pontzen & Governato (2014) Pontzen A., Governato F., 2014, Nature, 506, 171
  * Power et al. (2003) Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14
  * Prada et al. (2012) Prada F., Klypin A. A., Cuesta A. J., Betancort-Rijo J. E., Primack J., 2012, MNRAS, 423, 3018
  * Primack (2012) Primack J. R., 2012, Annalen der Physik, 524, 535
  * Read & Gilmore (2005) Read J. I., Gilmore G., 2005, MNRAS, 356, 107
  * Read et al. (2006) Read J. I., Wilkinson M. I., Evans N. W., Gilmore G., Kleyna J. T., 2006, MNRAS, 367, 387
  * Romano-Díaz et al. (2008) Romano-Díaz E., Shlosman I., Hoffman Y., Heller C., 2008, ApJ, 685, L105
  * Saburova & Del Popolo (2014) Saburova A., Del Popolo A., 2014, MNRAS, 445, 3512
  * Saha & Cortesi (2018) Saha K., Cortesi A., 2018, ApJ, 862, L12
  * Salucci (2019) Salucci P., 2019, A&ARv, 27, 2
  * Salucci et al. (2020) Salucci P., Turini N., di Paolo C., 2020, Universe, 6, 118
  * Sattari et al. (2023) Sattari Z., et al., 2023, ApJ, 951, 147
  * Shankar et al. (2017) Shankar F., et al., 2017, ApJ, 840, 34
  * Shibuya et al. (2016) Shibuya T., Ouchi M., Kubo M., Harikane Y., 2016, ApJ, 821, 72
  * Sofue (2017) Sofue Y., 2017, PASJ, 69, R1
  * Sommer-Larsen & Dolgov (2001) Sommer-Larsen J., Dolgov A., 2001, ApJ, 551, 608
  * Soto et al. (2017) Soto E., et al., 2017, ApJ, 837, 6
  * Spergel & Steinhardt (2000) Spergel D. N., Steinhardt P. J., 2000, Phys. Rev. Lett., 84, 3760
  * Springel et al. (2021) Springel V., Pakmor R., Zier O., Reinecke M., 2021, MNRAS, 506, 2871
  * Starobinsky (1980) Starobinsky A. A., 1980, Physics Letters B, 91, 99
  * Swaters et al. (2003) Swaters R. A., Madore B. F., van den Bosch F. C., Balcells M., 2003, ApJ, 583, 732
  * Sweet et al. (2019) Sweet S. M., et al., 2019, MNRAS, 485, 5700
  * Tamburello et al. (2015) Tamburello V., Mayer L., Shen S., Wadsley J., 2015, MNRAS, 453, 2490
  * Teyssier et al. (2013) Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068
  * Tollet et al. (2016) Tollet E., et al., 2016, MNRAS, 456, 3542
  * Walker & Peñarrubia (2011) Walker M. G., Peñarrubia J., 2011, ApJ, 742, 20
  * White & Rees (1978) White S. D. M., Rees M. J., 1978, MNRAS, 183, 341
  * Zanella et al. (2019) Zanella A., et al., 2019, MNRAS, 489, 2792
  * Zentner & Bullock (2003) Zentner A. R., Bullock J. S., 2003, ApJ, 598, 49
  * Zentner et al. (2022) Zentner A., Dandavate S., Slone O., Lisanti M., 2022, J. Cosmology Astropart. Phys., 2022, 031
  * Zhengra et al. (2019) Zhengra C., Posti L., Koopmans L. V. E., Napolitano N. R., 2019, MNRAS, 489, 5483
  * de Blok et al. (2001) de Blok W. J. G., McGaugh S. S., Rubin V. C., 2001, AJ, 122, 2396
  * de Blok et al. (2003) de Blok W. J. G., Bosma A., McGaugh S., 2003, MNRAS, 340, 657
  * van den Bergh et al. (1996) van den Bergh S., Abraham R. G., Ellis R. S., Tanvir N. R., Santiago B. X., Glazebrook K. G., 1996, AJ, 112, 359

##  Appendix A Host Equilibrium Check and convergence test

Figure 18: The dark matter halo density profile over 3 Gyr of evolution time,
in a system with disc and halo (without massive clumps). The two vertical
dashed lines denote the convergence radius at the maximum simulation time of 3
Gyr, $r_{\rm cov}=0.36\;\rm kpc$ (cf. equation 7), and NFW the scale radius
$r_{s}$. Note that the results shown in this work concern a period of
effective evolution (after interval of relaxation to quasi-equilibrium) of 0.7
Gyr, and thus the relevant convergence radius is smaller than shown here.

We seek to verify that our disc-halo configurations do not evolve
significantly in the absence of the stellar clumps; that they quickly reach a
quasi-steady state on a timescale much smaller than that associated with the
effect of dynamical friction mediation between the clumps and the host system,
and that they remain in that state. In Fig. 18 we show results of a simulation
that is run with disc and halo only, without the added stellar clumps. The
convergence radius $r_{\rm cov}$ is defined such that the maximum simulation
time ($t=3\;\rm Gyr$) is equal to the two-body relation time $t_{\rm relax}$
which is given by (Power et al., 2003),

$t_{\rm relax}(r)=\frac{\pi N(<r)}{4\ln N(<r)}\sqrt{\frac{r^{3}}{GM(<r)}},$
(7)

where $N$ and $M$ are the host particles number and mass within radius $r$. As
the results shown in the rest of this work concern a period of less than 1
Gyr, the the relevant convergence radius is smaller still. Also, as can be
seen, except for initial relaxation lasting about 30 Myr (much less than a
typical clump orbital period), the dynamical structure of the host system
remains stable.

Figure 19: Dark matter halo density profile in simulation with the same host
and clump properties but different mass resolutions. The results are shown at
$t=1.2\;{\rm Gyr}$, with the dashed-black line showing the density
distribution at $t=0.3\;{\rm Gyr}$. The internal clump profile is set to the
Jaffe model and the disc and halo parameters correspond to the fiducial runs
of Table 1.

We have also verified that our results are insensitive to the particle numbers
used in the simulation. Fig. 19 shows an example of such a control run.

##  Appendix B Truncated clump mass as function of orbital radius

Figure 20: Mass versus orbital position of a one clump simulations with mass
$2\times 10^{8.5}M_{\odot}$, overplotted with estimations from the Jacobi
formula (eq. 5). The horizontal black-dashed line is at $2\times
10^{8.5}M_{\odot}$.

In order to examine mass stripping of clumps due to tidal forces in a more
controlled manner, and compare with estimates using the Jacobi radius, we ran
one-clump simulations. An example is shown here, where the clump is set
initially at orbital radius $R=11.9{\rm kpc}$, and has mass $M_{c}=4\times
10^{8.5}M_{\odot}$. In Fig. 20, the evolution of the clump mass as a function
of orbital radius is plotted together with theoretical estimates using the
Jacobi formula (equation 5). The extent of mass preservation depends on the
internal clump mass profile. But it is only sensitive to the particular
profile in the inner region. This enables the advent of effective initial
conditions, after about half a Gyr relaxation period, that are similar for all
internal clump mass distribution used in this work (as described in Section
3.3). Furthermore, the Jacobi estimates of the mass stripping agree relatively
well with the simulation results at larger orbital radii and fail at smaller
ones. This also reinforces the plausibility of the procedure used in Section
3.3, which involves truncating the clumps at the Jacobi radius at large radial
distances $R$ in the disc and letting the system of clumps relax, to the
approximate the observed distribution, from there.