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

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.


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(Navarro et al. , 1997, , ★ mahmoud.hashim@bue.edu.eg† amr.elzant@bue.edu.eg2010)) 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 Λ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 feedbackdriven 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. 2019Dutton et al. , 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(de Blok et al. , 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 ≳ 10 8  ⊙ , 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 dark matter hosts (van den Bergh et al. 1996; even if the latter situation may also, in principle, materialise in the context of hierarchical structure formation).However, since then, an increasing number of galaxies have been discovered to host massive and compact clumpy structures of the right mass.These clumps, which are more abundant at higher redshifts, may collapse in situ within a forming galaxy (Elmegreen et al. 2004(Elmegreen et al. , 2007(Elmegreen et al. , 2009;;Genzel et al. 2011;Förster Schreiber et al. 2011;Guo et al. 2015Guo et al. , 2018;;Messa et al. 2022;Martin et al. 2023;Sattari et al. 2023).Their existence 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. 2014Mandelker et al. , 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(Genzel et al. , 2011;;Shibuya et al. 2016;Soto et al. 2017;Zanella et al. 2019).
Nevertheless, the lifetime of of high redshift clumps, particularly primarily gaseous ones, is still uncertain; they may 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 stellar masses of the clumps (e.g., Dessauges-Zavadsky et al. 2017;Fisher et al. 2017a;Cava et al. 2018;Faure et al. 2021).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 analogues 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 materialises 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 stellar 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.

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 DY-NAMO 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  = 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 to 6.4 × 10 10  ⊙ , the upper limit being close to that of the Milky Way (Licquia & Newman 2015).

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  obs .These range from 4.71 × 10 6 to 3.74 × 10 9  ⊙ , with average mass 2.53 × 10 8  ⊙ .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, M obs = r n obs , with  = 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 clumpcorresponds to the half mass and half mass radius of the clumps.For which we will adopt various mass models to describe (Section 3.2).

Clumps mass function
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): where  = log( obs ),  = −1.4 and  0 = 10 9.5  ⊙ .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   is assumed to correspond to twice the masses  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   and radial location  of the clumps in the observed sample, we start the our simulated clumps at large  (relative to the observed ones), with masses 2  , 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.4.2.

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 particlemesh 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  = 0.035 kpc.The maximum simulation time is set to 3 Gyr.In the runs described below, the total number of particles composing the halo is  = 256 3 , with mass resolution of   = 5.9 × 10 4  ⊙ .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.

The halo and the disc
The N-Body simulations presented here are of clumps embedded in an isolated galaxy composed of: (i) A dark matter halo following a Navarro,Frenk and White (NFW) density profile (Navarro et al. 1997), where   and   are density and radial scales.
In our fiducial model the total host (disc + halo) mass is 10 12  ⊙ , so that the halo mass is (10 12 −   )  ⊙ , where   is the disc mass fraction (taken as 0.02 in our fiducial model).We assume a flat-ΛCDM cosmology (with Hubble constant  0 = 70.0 and matter density (Ω  = 0.3), so that the halo viral radius is  vir ≡   = 204.0kpc, where  = 20.4 is the halo concentration.This is larger than the average found in simulations at this mass scale (which is closer to  = 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  = 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  = 230 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. A1.
(ii) A stellar disc, with density profile given by: Table 1.The fiducial host (disc + halo) mass distribution parameters.The total mass of the disc and halo combined is M host = 10 12 M ⊙ , and the disc mass fraction is   = 0.02.The exponential disc has vertical and radial scale length   and   , and the NFW halo has concentration  and radial scale length   .In addition to those fiducial parameters we ran simulations with less concentrated halo, with  = 10.2, and also ones where the halo mass is halved and disc mass doubled.
where  = √︁  2 +  2 ,   is the radial disc density scale and   is the vertical scale height.To ensure stability against significant nonaxisymmetric instabilities, such as strong bars, we choose a large value for the Toomre parameter, namely  = 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  = 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).

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/; iii) a still more cuspy Jaffe profile, with density going as 1/ 2 in the central region.The associated mass distributions are where   ,   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.The (isotropic) internal velocities are set by the spherical Jeans equation.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.

Tidal truncation of the clumps
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   and orbital radius , this is defined as (Binney & Tremaine 2008): where  host (< ) is the enclosed mass of the host (that is disk plus halo mass) at radius  in the disc plane.The clump profile parameters are then determined from the condition that the masses of the clumps be completely contained within   and that half the mass be contained within the observed radius  obs =  1/2 .Fig. 3 shows a sketch illustrating the two radii that thus fix the in- 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 R = 1kpc, where the remnants of larger clumps are expected to merge into a bulge after completing their central migration.The scaling radius  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.
ternal clump mass distribution, namely the half mass radius  1/2 and the truncation radius   .As a quantitative example of the resulting mass profiles, Fig. 4 shows those corresponding to a clump of half size  1/2 = 0.27 kpc and total mass   = 10 8.8  ⊙ , located at an orbital radius  = 9.64 kpc in the disc plane.

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.

Clumps
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  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 <  < 16.0kpc, with a normal distribution for their radii and random azimuthal angles, and imposing the condition that the Jacobi radius   of each clump is bigger than its half mass radius  obs (which naturally orders them in radii  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 quasisteady state, which may be assumed to correspond to the effective start of the simulation, We denote this time with  = 0.
We illustrate the procedure in Fig. 5.1 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.

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.

Energy transfer from clumps to cusp and the transformation to quasi-isothermal core
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 .This is derance of massive clumps in the central region, and in turn the associated effects due to dynamical friction coupling that we describe below.
Figure 8. Change in total energy of subsystems of halo particles enclosed within radius , as given by equation 6, as they gain energy from the inwardly migrating stellar clumps.The results are illustrated for the fiducial dischalo model of Table 1, and for the three internal clump profiles considered: Plummer spheres with constant density cores; Hernquist spheres, with 1/ inner cusps; and Jaffe spheres, with 1/ 2 inner cusps.
defined as where   refers to the mass of a particle (recall that all particles have the same mass in our simulations),   is the speed of particle  and Φ  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 .
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.
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).
When a less concentrated halo profile is used (initial NFW  = 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).

Cusp-core transformation
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.
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  = 20.4.case.)

Total density 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 ∼ −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).

Rotation curves: diversity and bulge halo conspiracies
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 pre-vious 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 re-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  = 0.7 Gyr (solid lines), and initially at  = 0.0 (dashed lines).placement 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  = 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  = 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.

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 < ∼ 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/ cusp), and the Jaffe profile (with 1/ 2 cusp).We find that the steeper the inner clump profileand 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 the 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 purely N-body simulations of the observed stellar clumps naturally do not include feed-back, from either starburst 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.They are also supported by simulations showing that the relative preponderance of the most massive stellar clumps, which are primarily behind the processes described here, may actually increase as a result of feedback (Faure et al. 2021).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 in tandem with detailed simulations.Myr (much less than a typical clump orbital period), the dynamical structure of the host system remains stable.
We have also verified that our results are insensitive to the particle numbers used in the simulation.Fig. A2 shows an example of such a control run.

Figure 1 .
Figure 1.The DYNAMO-HST clumps mass-size relation (from the data ofAmbachew et al. 2022).The solid line is the power law fit to the masssize scaling.The horizontal black-dashed line represents the observational mass limit of 10 7  ⊙ .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   = 2 obs = 2 1/2 , and with internal mass distributions modelled through three different physical density profiles (Section 3.2).

Figure 2 .
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  = 0, after a 500 Myr relaxation period as described in Section 3.4.2).The theoretical lines show models with halo concentrations  = 20.4 and  = 10.2.

Figure 3 .
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 .
Figure 4.The different mass profiles of a clump with mass   = 10 8.8  ⊙ .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.

Figure 5 .
Figure 5. Setting the effective initial conditions.The initial (at  = 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  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  = 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 R = 1kpc, where the remnants of larger clumps are expected to merge into a bulge after completing their central migration.The scaling radius  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.

Figure 7 .
Figure 7. Projected density maps of the stellar disc and clumps, for simulations of the fiducial model (Table1), 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  = 0.5 Gyr, corresponding to  = 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).

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

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

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

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

Figure 16 .
Figure 16.Initial ( = 0) and final ( = 0.7Gyr) rotation curves in simulation with Jaffe-profile clumps and halo concentration  = 20.4,but with halo mass reduced by half and stellar mass (including both disc and clumps) doubled.

Figure 17 .
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.

Figure A2 .
Figure A2.Dark matter halo density profile in simulation with the same host and clump properties but different mass resolutions.The results are shown at  = 1.2 Gyr, with the dashed-black line showing the density distribution at  = 0.3 Gyr.The internal clump profile is set to the Jaffe model and the disc and halo parameters correspond to the fiducial runs of Table1.