Emergence of a stellar cusp by a dark matter cusp in a low-mass compact ultra-faint dwarf galaxy
Shigeki Inoue


firstpage–lastpage 2014

Relating the finite-volume spectrum and the two-and-three-particle S matrix
 for relativistic systems of identical scalar particles
    Stephen R. Sharpe
    December 30, 2023
===================================================================================================================================







  Recent observations have been discovering new ultra-faint dwarf galaxies as small as ∼20  pc in half-light radius and ∼3  km s^-1 in line-of-sight velocity dispersion. In these galaxies, dynamical friction on a star against dark matter can be significant and alter their stellar density distribution. The effect can strongly depend on a central density profile of dark matter, i.e. cusp or core. In this study, I perform computations using a classical and a modern analytic formulae and N-body simulations to study how dynamical friction changes a stellar density profile and how different it is between a cuspy and a cored dark matter haloes. This study shows that, if a dark matter halo has a cusp, dynamical friction can cause shrivelling instability which results in emergence of a stellar cusp in the central region 2  pc. On the other hand, if it has a constant-density core, dynamical friction is significantly weaker and does not generate a stellar cusp even if the galaxy has the same line-of-sight velocity dispersion. In such a compact and low-mass galaxy, since the shrivelling instability by dynamical friction is inevitable if it has a dark matter cusp, absence of a stellar cusp implies that the galaxy has a dark-matter core. I expect that this could be used to diagnose a dark matter density profile in these compact ultra-faint dwarf galaxies.




instabilities – methods: numerical – methods: analytical – galaxies: dwarf – galaxies: kinematics and dynamics.





§ INTRODUCTION


Dark matter (DM) density profiles in dwarf galaxies have long been debated. Theoretical studies such as cosmological N-body simulations have demonstrated that DM density increases toward the galactic centre independent of a halo mass <cit.>. On the other hand, observations proposed that dwarf galaxies seem to have nearly constant densities of DM at their central regions <cit.>.[Low surface brightness galaxies have also been observed to have DM cores <cit.> although I do not discuss these galaxies.]

As a possible solution, if DM haloes consist of warm or self-interacting particles, all dwarf galaxies are expected to have central DM cores. It has also been proposed, alternatively, that (recursive) baryonic feedback can turn a cusp into a core by flattening the inner slopes of the primordial DM density profiles in dwarf galaxies as massive as M_ DM∼ 10^10–10^11  M_⊙ <cit.>. If the latter scenario is the case, since dynamical masses of some ultra-faint dwarf galaxies (UFDs) in the local group have been observed to be significantly smaller than the mass threshold above which the baryonic effect is influential to their central DM densities, they could be expected to preserve the primordial DM density profiles which may be cuspy. Accordingly, it is interesting to try to determine DM density profiles of such low-mass UFDs. It is, however, still impossible to know whether their DM haloes have cusps or cores because of only a handful of stars observable by spectroscopy to measure their line-of-sight velocities (LOSVs) and model their DM haloes. Hence, it is worthwhile looking for an alternative method to deduce which type of DM the low-mass galaxies have, cusp or core. For example, <cit.> have proposed a method using a fraction of wide binaries which can be disrupted by tidal force depending on their DM potential in UFDs.

The smallest UFDs are as tiny as R_ h30  pc in half-light radius and L ∼10^2-3  L_⊙ in luminosity <cit.> although current observations still cannot reject the possibility that some of them are extended globular clusters. Recently, <cit.> has analytically discussed that dynamical friction (DF) on a star against dark matter may be marginally effective on the timescale of ∼10  Gyr in Draco II — observed physical properties of which are R_ h=19^+8_-6  pc, brightness M_ v=-2.9±0.8, LOSV dispersion σ_ h=2.9±2.1  km s^-1 measured within R_ h <cit.> — by adopting the Chandrasekhar DF formula <cit.> to his singular isothermal DM halo model. His result implies that DF against DM could significantly change stellar distribution in UFDs more compact and/or less massive than Draco II, which will be discovered by future observations.

The effect of DF strongly depends on the DM density profile. It has been known that DF drag force becomes significantly weaker in cored density distribution than in cuspy one, once a massive particle enters the core <cit.>. Studies using N-body simulations have demonstrated that drag force by DF does cease practically in a constant-density core, probably by non-linear effects <cit.>. Therefore, if an extremely low-mass UFD has a constant-density core of DM, DF could be too weak to affect the stellar distribution. On the other hand, if such a UFD has a DM cusp, DF against DM could be strong enough to make alterations to its stellar distribution, such as emergence of a stellar cusp or formation of a nucleus cluster as a remnant of stars fallen into the galactic centre. Current observations of low-mass compact UFDs are limited to a close distance of d30  kpc from the sun because of their faintness. At this distance, each star in a UFD can be resolved since the typical size of observational smearing is smaller than the mean separation of stars even at the galactic centres. Therefore, the expected stellar cusp and  the nucleus cluster would be observed as a dense group of stars at the galactic centre if it exists.

This study addresses the effect of DF by DM on stellar distribution in an extremely low-mass and compact UFD and focus on how different it is between cuspy and cored DM density profiles. In Section <ref>, I perform analytical estimation of stellar shrivelling due to DF based on a classical and a modern DF formulae. In Section <ref>, I perform N-body simulations resolving every single star and demonstrate the same as the analytical estimation presented in Section <ref>. Finally, I present discussion and summary of this work in Section <ref>.



§ ANALYSIS USING DYNAMICAL FRICTION FORMULAE


In this study, I describe a density profile of a DM halo by a Dehnen model <cit.>,

    ρ_ DM(r) = ρ_ DM,0r_ s^4/r^γ(r+r_ s)^4-γ,

where ρ_ DM,0 and r_ s are scale density and radius, and γ is an inner density slope. I assume a cuspy DM halo to be represented by setting γ=1, which corresponds to the Hernquist model profile <cit.>, and a cored DM halo is represented by γ=0. Local velocity dispersion of DM is generally computed by solving Jeans equation, 

    σ_ DM^2(r) = 1/ρ_ DM∫^∞_rρ_ DMGM_ DM(r')/r'^2 dr',

where G is the gravitational constant, and M_ DM(r) is mass of DM enclosed within r. Here, I assume that gravity of a baryon component is negligible and that velocity distribution is isotropic. The analytic solutions of σ_ DM for γ=0 and 1 can be found in <cit.> and <cit.>.

This study discusses how the DF against DM affects stellar distribution and how different it is between cuspy and cored DM haloes. I use a Plummer's model for stellar distribution of a compact UFD,

    ρ_⋆(r) = 3M_⋆/4π r_⋆^3(1+r^2/r_⋆^2)^-5/2,

where M_⋆ and r_⋆ are the total mass and scale radius of stars. This model has a core of stars in r≪ r_⋆ where the density is nearly constant. Two-dimensional half-light radius R_ h and integrated mass M_ h inside R_ h are obtained by assuming a constant mass-to-luminosity ratio and integrating equation (<ref>).  In this study, I assume r_⋆=20  pc (R_ h=r_⋆ in a Plummer's model). Luminosity-weighted LOSV dispersion inside R_ h is given as 

    σ^2_ h = 4π/M_ h∫^∞_0dz∫^R_ h_0ρ_⋆σ^2_⋆(R',z)R' dR',

where stellar velocity dispersion σ_⋆ is computed from equation (<ref>) in which ρ_⋆ is substituted for ρ_ DM. Since stellar gravity is now assumed to be negligible, setting σ_ h gives ρ_ DM,0 when the other parameters in equation (<ref>) are fixed:[When σ_ h=1.5  km s^-1, the values of ρ_ DM,0 in the cored DM models are 9.8, 7.0, 5.8 and 5.1×10^-1  M_⊙ pc^-3 for r_ s=125, 250, 500  pc and 1  kpc. Those in the cuspy models are 1.1, 0.45, 0.20 and 0.096×10^-1  M_⊙ pc^-3, respectively.] σ_ h^2∝ρ_ DM,0. In what follows, I discuss the two cases of σ_ h to 1.5 and 3.0  km s^-1.[The cusp and the core models of equation (<ref>) have the finite total masses, M_ cusp,tot=2πρ_ DM,0r_ s^3 and M_ core,tot=(4/3)πρ_ DM,0r_ s^3, respectively. When σ_ h=1.5  km s^-1, for r_ s=125  pc, the total masses of cuspy and cored haloes are M_ cusp,tot=1.3×10^6  M_⊙ and M_ core,tot=8.0×10^6  M_⊙. For r_ s=1  kpc, M_ cusp,tot=6.0×10^7  M_⊙ and M_ core,tot=2.2×10^9  M_⊙.] Fig. <ref> illustrates radial profiles of σ_ DM, circular velocities v_ circ≡√(GM_ DM(r)/r) and σ_⋆ normalised by σ_ h in my cuspy and cored halo models with r_ s=125  pc and 1  kpc.




 §.§ Analytic formulae of dynamical friction



  §.§.§ The Chandrasekhar formula



I consider DF against DM on a star. The Chandrasekhar DF formula under Maxwellian velocity distribution[Although <cit.> has also proposed more general forms of analytic DF not relying on Maxwell distribution, I refer to equation (<ref>) as Chandrasekhar formula in this paper.] is given as 

    F_ DF = -4πlnΛ G^2ρ_ DMm_⋆^2/v_⋆^2[ erf(X)-2X/√(π)exp(-X^2)],

where m_⋆ and v_⋆ are a mass and a velocity of a star, X≡ v_⋆/(√(2)σ_ DM), and Λ is a parameter, whose proper value is still under debate <cit.>[Basically, Λ is defined to be a ratio between the minimum and maximum impact parameters of two-body gravitational interaction, i.e., Λ≡ b_ max/b_ min, where b_ min∼ Gm_⋆/v_⋆^2, b_ max∼ r_ s in the classical formula <cit.>.]. The direction of F_ DF is presumed to be opposite to the velocity vector v_⋆. By assuming a circular orbit, i.e. v_⋆=v_ circ, equation (<ref>) can be solved with σ_ DM from equation (<ref>). In the case considered here, equation (<ref>) is independent of ρ_ DM,0 (i.e. σ_ h) since ρ_ DM, v_⋆^2 and σ_ DM^2 are proportional to ρ_ DM,0 although the DF timescale, ∼ m_⋆σ_ h/F_ DF, depends on ρ_ DM,0. In addition, because the parentheses in equation (<ref>) is independent of m_⋆, F_ DF/m_⋆^2 is independent of m_⋆.

The left panel of Fig. <ref> shows the Chandrasekhar DF force of equation (<ref>) with lnΛ=15 normalised by m_⋆^2 in the cuspy and the cored haloes with various r_ s. The Figure indicates that the strength of DF is remarkably different between the cuspy and the cored DM haloes; DF increases monotonically towards the centre in a DM cusp (the blue lines), whereas it is approximately constant or gently decreases in a core (the red lines). The behaviour of DF is almost independent of r_ s in the cuspy haloes, whereas DF in a core becomes weaker when r_ s is larger. Although the difference of DF inside r_ s between the cusp and the core becomes smaller with decreasing r_ s, it is still quite large in r10  pc even in the case of r_ s=125  pc. This means that, if a stellar component of a UFD is deeply embedded in a DM halo (i.e. R_ h≪ r_ s), one can expect that DF strongly depends on a density profile of DM and could change the stellar distribution if the compact UFD has a cusp of DM. Moreover, in a cuspy halo, a star undergoing DF migrates into an inner radius, at which DF is even stronger (see Section <ref>).



  §.§.§ Petts et al. formula


The Chandrasekhar formula of equation (<ref>) is, however, based on several assumptions such as the Maxwellian velocity distribution and the invariable parameter of Λ. Therefore, improved formulae have been invented by previous studies. Recently, <cit.> proposed more sophisticated DF modelling based of the general Chandrasekhar formula <cit.>, which uses a distribution function instead of the Maxwellian distribution and takes high-velocity encounters into account. They demonstrated that their improved DF model can reproduce orbits of infalling particles in cuspy and cored density fields better than the classical formula. They formulate DF force[<cit.> proposed two models of DF: `P16' and `P16f'. I use their P16f model in this study since they concluded that P16f is more accurate than P16.] as 

    F_ DF = -2π^2 G^2ρ_ DMm_⋆^2/v_⋆^2∫^v_ esc_0J(v_ DM)f(v_ DM)v_ DM dv_ DM,



    J= ∫^v_⋆+v_ DM_|v_⋆-v_ DM|(1+v_⋆^2-v_ DM^2/V^2)log(1+b_ max^2V^4/G^2m_⋆^2) dV,

where f(v_ DM) represents a distribution function of DM, which is defined so that 4π∫ f(v_ DM)v_ DM^2 dv_ DM=1, and escape velocity v_ esc=√(-2Φ), and V corresponds to relative velocity of encounter. In equation (<ref>), the maximum impact parameter

    b_ max = min( ρ_ DM(r)/dρ_ DM/dr, r).

The value of ρ_ DM/(dρ_ DM/dr) can be taken as the distance within which the density field can be considered to be homogeneous <cit.>. However, since it can diverge in a constant-density core, b_ max is limited to be ≤ r <cit.>. Equation (<ref>) can be solved analytically (see Appendix <ref>). 

In addition, <cit.> also introduced a `tidal-stalling' radius, r_ TS, at which the tidal radius of a massive particle is equal to its orbital radius. The tidal radius is 

    r_ t = Gm_⋆/Ω^2-d^2Φ/dr^2,

where Ω^2=GM_ DM/r^3. They argued that DF ceases within the tidal-stalling radius because of non-linear effects, and showed that the radius of r_ TS matches well results of N-body simulations (a radius of DF cessation in a core, see Section <ref>). In the  DF model, F_ DF=0 in r<r_ TS although equation (<ref>) still returns a non-zero value.  

The right panel of Fig. <ref> shows the  DF force of equation (<ref>) normalised by m_⋆^2. Unlike the Chandrasekhar formula, now F_ DF/m_⋆^2 weakly depends on σ_ h and m_⋆. Here, I assume σ_ h=1.5  km s^-1 and m_⋆=0.5  M_⊙, however the results hardly change between σ_ h=1.5 and 3.0  km s^-1. Radii of r_ TS become about 1.4 times smaller when σ_ h=3.0  km s^-1. In the right panel of Fig. <ref>, although the differences between the cuspy and the cored haloes are still quite large, it is remarkable that the  formula predicts DF significantly stronger than the classical formula in the central regions of the cored haloes <cit.>. On the other hand, the DF in the cuspy haloes is similar to that given by the Chandrasekhar formula.



 §.§ Orbital integration with the formulae


Using the models and the DF formulae described above, I perform orbital integration of stars under the potential given by the DM distribution of equation (<ref>). With the initial spatial distribution of equation (<ref>), the velocity distribution of stars is given by Eddington's formula <cit.> with isotropy. I do not take into account mutual interactions between the stars, therefore the result is independent of the number of stars. For the sake of statistics, I use a random sample of ten million stars in each run. While integrating their orbits with respect to time, the stars are decelerated by DF represented by the analytic formulae every timestep. I assume m_⋆=0.5  M_⊙ as a typical mass of a star as old as ∼10  Gyr <cit.>. The analytic DF is considered to work until a star reaches the radius r_ limit at which M_ DM(r_ limit)=m_⋆. In my models, r_ limit≃0.1 and 0.5  pc in the cuspy and the cored DM haloes. When a star enters r_ limit with a velocity slower than v_ circ|_r=r_ limit, the star is stopped there and considered to be fallen into the galactic centre by DF. When the Chandrasekhar formula is applied, I set lnΛ=15. When the  formula is applied, DF ceases within r_ TS (i.e. F_ DF=0). I use a second-order leap-frog integrator for the orbital computations with a constant and shared timestep of Δ t=0.01× r_ limit/v_ circ|_r=r_ limit. I confirmed the convergence of my results with respect to Δ t. Orbits of stars are time-integrated until t=10  Gyr, and I obtain stellar surface densities as functions of radius in the runs.


Fig. <ref> shows the results of the orbital integration using the Chandrasekhar (top) and  (bottom) formulae in the cuspy and the cored haloes with r_ s=125  pc and σ_ h=1.5  km s^-1. Although the stellar surface density is nearly constant in R5  pc in the initial state (black line), the density profiles are significantly steepened after t∼4  Gyr in the cuspy DM halo (blue lines). In addition, sharp stellar cusps like nucleus clusters emerge at the galactic centres, which have five and four per cents of the total number of stellar particles within R<0.5  pc in the top and the bottom panels. The cusps mainly consist of stars fallen into the centres by DF. The steepened stellar distribution profiles are nearly exponential outside the stellar cusps. Because the Chandrasekhar and the  formulae are not significantly different (Fig. <ref>), the results of Fig. <ref> are similar in the cuspy halo. These results corroborate the expectation that, as <cit.> proposed,  a stellar distribution in a low-mass compact UFD can be affected by DF against DM if it has a cusp. 

If a DM halo has a core, however, the DF approximated by the analytic formulae is significantly less efficient to steepen the stellar profile  (the red lines), in spite of the same σ_ h meaning similar DM masses within R_ h. When the Chandrasekhar formula is applied (the top panel), the stellar density hardly changes even at t=10  Gyr in the cored DM halo. Although the  formula (the bottom panel) steepens the stellar density profile more than the classical formula, the density slope is clearly shallower than that in the cuspy DM, and a stellar cusp does not form. No stars are fallen into the centre by either DF modellings. The absence of the stellar cusp is due to the week DF in the DM core and the DF cessation assumed in r<r_ TS in the  model.


Fig. <ref> shows the same results but for different settings for  r_ s and σ_ h. The effect of DF becomes weaker for larger r_ s and σ_ h (i.e. higher ρ_ DM,0), but the steepening of a surface density profile and formation of a stellar cusp by a DM cusp can be seen even when r_ s=1  kpc and σ_ h=3.0  km s^-1. In the case of cored DM haloes, on the other hand, the stellar density profiles are almost intact even though σ_ h and r_ s are the same as in the cuspy halo models. In the case of the cored halo with σ_ h=3.0  km s^-1 and r_ s=125  pc, the  formula predicts weak steepening, but a stellar cusp does not emerge. Thus, significance of the DF effect strongly depends on whether the DM halo has a cusp or a core even if the r_ s and σ_ h are the same. The most noticeable difference is the emergence of a stellar cusp in a DM cusp.

The total mass of the nucleus remnants consisting of stars fallen into the centre can depend not only on DM density but also stellar distribution. If R_ h is larger, stars have more extended distribution, therefore DF timescale becomes longer on average. Thus, a larger R_ h leads a smaller fraction of stars to fall into the centre by DF. As a result, a less prominent stellar cusp would form in such an extended galaxy.



§ N-BODY SIMULATIONS


As I showed in Section <ref>, the analytic formulae are useful to estimate the magnitude of DF. The formulae, however, still ignore non-linear effects. For example, they assume DF as a corrective effect of two-body interactions and do not take into account orbital periodicity of particles or reaction of field particles. To address further the effect of DF using more realistic models, I perform N-body simulations in which the models are fully self-consistent, and DF drag force naturally arises as mutual interactions between particles.



 §.§ Settings

The initial conditions of my N-body simulations are the same as the DM and stellar models (equation <ref> and <ref>) with the parameters used in Fig. <ref>: r_ s=125  pc, σ_ h=1.5  km s^-1 (ignoring stellar potential) and b=20  pc. The total stellar mass is set to M_⋆=2500  M_⊙, and a mass of a single stellar particle is m_⋆=0.5  M_⊙, i.e. the number of stellar particles is N_⋆=5000. A stellar particle has a softening length of ϵ_⋆=0.1  pc. Velocity distribution is given by Eddington's formula taking into account the total potential of the DM and the stars. Although the actual LOSV dispersion of stars inside R_ h is slightly higher than 1.5  km s^-1 because of self-gravity of the stars, the increase is only a few per cent. Although every single star is resolved with a point-mass particle, the interactions between stars in the simulations are still collisionless (see Appendix <ref>).

DF can arise if m_⋆≫ m_ DM, where m_ DM is a mass of a DM particle in simulations. Since m_⋆=0.5  M_⊙ in my simulations, m_ DM should be 0.05  M_⊙. Achieving such a high resolution requires approximately 2.6 and 16.0×10^7 particles for the cuspy and the cored DM haloes. To lighten the heavy burden of the N-body computations, I employ an orbit-dependent refinement method for a multi-mass spherical model proposed by <cit.>. This method divides a DM halo into i shells and the central sphere (the zeroth shell). Basically, each shell is resolved into DM particles with each mass resolution m_ DM,i and softening length ϵ_i (see Table <ref>). After assigning a DM particle in the i-th shell its initial position and velocity and computing its pericentre distance in the fixed potential, if the pericentre intrudes into the inner j-th shell, the particle is split into m_ DM,i/m_ DM,j particles with the mass m_ DM,j and the softening length ϵ_j.[Therefore, the mass ratio m_ DM,i/m_ DM,j has to be a natural number.] The split particles are distributed on random positions while keeping the initial radius of the parent particle, and directions of their tangential velocities are randomly reassigned while keeping the initial radial velocity and the kinematic energy of the parent particle. This refinement method can, by a substantial factor, reduce the computational run time by decreasing the number of DM particles in outer regions that are not important to this study, while preventing the outer particles with larger masses from entering the innermost region resolved with the smallest particle mass. After the refinement, 1.53 and 6.88×10^7 particles are required to represent the cuspy and the cored DM haloes.


I use a simulation code ASURA <cit.>,[ASURA is an N-body/smoothed particle hydrodynamics (SPH) code although this study only uses the N-body part.] in which a symmetric form of a Plummer softening kernel <cit.>, a parallel tree method with an computational accelerator GRAPE <cit.> and the second-order leap-frog integrator with individual timesteps are used. The number of stellar particles, N_⋆=5000, in my simulations may be too small to obtain a statistically certain density profile. To reinforce this point, I perform ten runs with the same initial condition but different random-number seeds.



 §.§ Results




  §.§.§ Evolution of the stellar density profiles


I obtain three surface density profiles observed from perpendicular angles for each of the ten runs. Then, I compute a stacking of the thirty profiles of stellar surface density at the same time t for each case of the cuspy and the cored halo. The centre of the stellar distribution is defined to be the median position among all stellar particles in each snapshot.



Fig. <ref> shows the stackings of stellar surface density profiles in the cuspy (blue) and the cored (red) DM haloes at t=4, 7 and 10  Gyr. The shaded regions indicate the ranges of upper and lower 1σ-deviations of the stackings. In the DM cusp, the stellar density profile clearly demonstrates the emergence of a stellar cusp after t=7  Gyr; the density slope becomes remarkably steeper in R2  pc than that in R2  pc. On the other hand, such a stellar cusp does not emerge in the cored DM halo although the outer density slope of the stars in R2  pc is similar to that in the cuspy DM halo, which is nearly exponential with radius. From this result, it can  be seen that a DM cusp can generate a stellar cusp by DF even if the stellar density profile is initially flat. In the N-body simulations, the stellar cusps have masses of 33.8^+4.5_-5.3  M_⊙ within R<2  pc, which corresponds to 1.4 per cent of the total stellar mass. Additionally, the difference between the two cases is significant in spite of the same σ_ h, which means that the two halo models are considered to be similar in observations. The emergence and the absence of stellar cusps in the cuspy and cored DM are approximately consistent with the results of my orbital integration models using the analytic DF formulae (Fig. <ref> and <ref>).



Fig. <ref> shows the evolution of R_ h and σ_ h during the N-body simulations. Stellar half-mass radii R_ h decrease slightly; by ≃0.5 and 1  pc in the cored and cuspy DM. LOSV dispersions σ_ h with in R_ h are almost constant even after the emergence of the stellar cusps. This means that DF is not effective for most of stars around the half-mass radius although the central regions in r≪ R_ h are significantly affected.



  §.§.§ Evolution of the DM density profiles

As I showed above, DM exerts DF on stars and can cause a low-mass compact UFD to have a stellar cusp if it has a cuspy DM halo. On the other hand, the DM particles can be kinematically heated by the stars spiraling into the centre, as the reaction of DF. Previous studies have shown that a DM cusp can be disrupted or made shallower by objects spiraling into the centre <cit.>. <cit.> also demonstrated that a DM cusp disrupted by infalling objects can be revived if a central remnant of the infalling objects is sufficiently massive. Hence, it is also interesting to look into evolution of the dark matter density profiles in my N-body simulations, i.e. whether the cuspy halo is still cuspy or cored after the creation of the stellar cusp.
 
Fig. <ref> shows DM density profiles in my N-body simulations of the cuspy halo model, in which I make a stacking of the ten runs. Here, the halo centre is defined to be a position of the particle that has the highest DM density in each snapshot. I use a method like SPH to compute the local DM densities for the centering; a cubic spline kernel is applied to 128 neighbouring DM particles. The Figure indicates that the DM cusp in the initial state is significantly weakened in r1  pc at t=4  Gyr (orange). Eventually, the initial DM cusp is turned into a core extending to r≃2  pc at t=10  Gyr (green). This result means that the central DM is kinematically heated by infalling stars, and the DM density is decreased in r2  pc. The size of this region where DM is affected is consistent with the size of the stellar cusp in the N-body simulations (Fig. <ref>). Since the softening length of DM particles is ϵ_0=0.1  pc, the peaks of DM densities at r≃0.1  pc may be transient fluctuation.

From the consistency of the sizes between the stellar cusps and the DM cores created, the size of a stellar cusp may be regulated by a DM density profile flattened by infalling stars. If this is the case, a larger number of stars falling into the centre can create a larger DM core and a broader stellar cusp. In the central region of a cuspy DM halo, the significance of DF basically depends on σ_ h[F_ DF is almost independent from r_ s in the  formula for cuspy haloes (Fig. <ref>).]. In addition, the number of stars in the central region where stars can reach the centre by DF within ∼10  Gyr depends on the initial stellar distribution, i.e. M_⋆ and R_ h.






§ DISCUSSION AND SUMMARY




 §.§ Summary and interpretation of the results


As I showed in Section <ref> and <ref>, DF on stars against DM can largely alter the stellar density distribution if the galaxy has cuspy DM distribution and a compact stellar component having R_ h≃20  pc and σ_ h3  km s^-1. The most important result obtained from my N-body simulations is that the DF by the DM cusp can arouse emergence of a stellar cusp in the galactic centre 2  pc. On the other hand, if a DM halo has a core, DF is not efficient to generate such a stellar cusp, in spite of the same σ_ h, although stellar density can be affected and increase slightly in a wide range 10  pc.

The results mentioned above can be explained by the differences of density and velocity dispersion between a cuspy and a cored DM haloes. According to the analytic formulae, DF becomes stronger when a background density is higher and a velocity dispersion is lower. In a cuspy halo, DM density increases toward the centre, and velocity dispersion decreases (see Fig. <ref>), therefore DF becomes stronger towards the centre. In this case, orbital shrinkage by DF brings a star to an inner region where DF is even stronger: `DF shrivelling instability' <cit.>. On the other hand, in a cored halo, density and velocity dispersion are nearly constant in the central region, therefore DF drag force is approximately independent of radius; it actually decreases gently towards the centre (Fig. <ref>). This means that a cored halo is relatively stable against the DF shrivelling of stars.

Interpretation of the above results should be considered carefully. It should be noted that presence of a stellar cusp in a low-mass compact UFD is not necessarily evidence to prove a DM cusp. It is because we do not know the initial condition of the stellar density; a galaxy can create a stellar cusp at its birth even if its DM halo has a core. It could be said, however, that it is inevitable to have a stellar cusp if a low-mass compact UFD has a DM cusp. In other words, if a low-mass compact UFD is observed to have no stellar cusp and be old enough, it suggests that its DM halo would have a large core with a size of r_ s100  pc. I discuss UFDs in current observations on this point in Section <ref>.



 §.§ The analytic formulae vs N-body simulations

It is interesting to compare the results of the analytic formulae with those of N-body simulations although it is not the main purpose of this study. Early studies using numerical simulations have discussed that the Chandrasekhar formula assuming Maxwellian velocity distribution and a invariant Λ can give a quite accurate estimate of DF force in various cases <cit.>. It was also reported, however, that the formula can be inaccurate in some specific cases; analyses and N-body simulations have shown that DF can be enhanced around a constant-density core, and then suppressed in the core <cit.>. These phenomena are inconsistent with predictions by the simplified Chandrasekhar formula. Various physical mechanisms of the deviation from the analytic formula have been proposed: orbital resonance between a massive and field particles <cit.>, coherent velocity field among particles <cit.>, a non-Maxwellian velocity distribution <cit.>, decrease of low-velocity particles <cit.> and inhomogeneity of background density and a variable Λ <cit.>. 

In the case of a DM cusp, the two analytic formulae predict similar DF force in Fig. <ref>, and the results of my orbital integration models are qualitatively consistent with the N-body simulations. However, the stellar cusp in my N-body simulations have the size of R≃2  pc, and it could be attributed to the weakened DM cusp shown in Fig. <ref>, where DF is weakened. In addition, the DM density centre is not necessarily be fixed onto the stellar centre in the simulations, and the slippage of the centres can broaden the stellar cusp. Therefore, the broadness of the stellar cusp in the simulations does not necessarily mean inaccuracy of the DF modellings. However, the stellar density slopes outside the stellar cusp is steeper in the orbital integration models. Moreover,  in my N-body simulations, the total mass of the stellar cusps within R=2  pc are about four times smaller than those predicted by the analytic DF modellings. On these points, both analytic formulae would be overestimating DF in the cuspy haloes.

In the case of a DM core, on the other hand, DF cannot generate a stellar cusp in either the orbital integration or the N-body models. However, there is a difference worthy of special mention: the simplified Chandrasekhar formula hardly changes the stellar density slopes, whereas the N-body simulations show significant increase of the stellar densities in a wide range of R10  pc (Fig. <ref>). This result shows that the Chandrasekhar formula assuming Maxwellian distribution and a invariable Λ underestimates DF in the cored halo despite that the classical formula overestimates in the cuspy halo. It is noteworthy that using the Petts et al. formula with their tidal-stalling model can dramatically improve the reproducibility of DF effect in the cored haloes, and the result of the stellar density profile is almost consistent with the N-body simulations (see the bottom panels of  Fig. <ref> and <ref>). Thus, the DF modelling proposed by <cit.> seems to be more accurate than the simplified Chandrasekhar formula although it may not be perfect yet in a DM cusp. Although it is beyond the scope of this study to investigate the physical reasons of the differences between the analytic formulae and my N-body simulations, I consider that the N-body simulations would be physically more credible than the analytic models.



 §.§ Comparison with observations


Here, I discuss the validity of my models of low-mass compact UFDs and the results in comparison with current observations. First, it is still very difficult or impossible to determine masses and sizes of DM haloes of UFDs with accuracy in current observations. I have to note, therefore, that the parameters in my DM halo models might be arbitrary. Recent observational studies have argued that galaxies have the universal DM surface density, μ_ DM≡ρ_ DM,0r_ s, over quite a wide range of luminosity when they are assumed to have cored DM haloes <cit.>.[This universality can be explained by assuming the Faber-Jackson law for DM haloes <cit.>.] <cit.> and <cit.> derived μ_ DM=140^+80_-30 and 70±4  M_⊙ pc^-2 from their galaxy samples including some satellite galaxies of the Milky Way. Although it has to be noted that they assumed different models for their cored haloes, my cored DM models of equation (<ref>) have μ_ DM=122 and 176  M_⊙ pc^-2 for r_ s=125 and 250  pc, respectively, when σ_ h=1.5  km s^-1. Accordingly, my cored halo model with r_ s=125–250  pc and σ_ h=1.5  km s^-1 would be consistent with the observed universality of μ_ DM if it is extrapolated to extremely low-mass galaxies. 

For my stellar model, I assume the uniform stellar mass of m_⋆=0.5  M_⊙. However, of course, stars generally have different masses according to their initial mass function and stellar evolution. Although stellar scattering by massive stars would not be efficient since encounters between stars are expected to be rare in a low-mass compact UFD (see Appendix <ref>), mass segregation does occur on the same timescale as DF because of mass-dependence of DF. Because more massive stars sink faster into the centre of a DM cusp, a stellar cusp would mainly consist of massive stars.[The massive objects are
stellar remnants. Even the typical mass of White Dwarfs is larger than 0.5  M_⊙.] The massive objects in the stellar cusp could be a heating source of less massive stars around it and prevent the less massive stars from falling into the centre. Thus, I have to note that my N-body simulations lack this effect.

The best UFD that is the most similar to my N-body model (R_ h≃20  pc, σ_ h=1.5  km s^-1 and log(M_⋆/M_⊙)=3.4) would be Draco II, which has R_ h∼19^+8_-6  pc, σ_ h=2.9±2.1  km s^-1, the total  luminosity log(L_ V/L_⊙)=3.1±0.3 and an age  ∼12  Gyr <cit.>. Although the most probable value of σ_ h observed is nearly twice higher than that in my model, it is within the error range. If I adopt a stellar mass-to-luminosity ratio M/L_ V=2–3 M_⊙/L_⊙ for a metal-poor system of 12  Gyr from a simple stellar population model of <cit.>, the stellar mass of Draco II is approximately log(M_⋆/M_⊙)=3.1–3.9. At the Heliocentric distance of Draco II, 20±3  kpc <cit.>, the size of a stellar cusp expected from my N-body simulations, ≃2  pc, corresponds to ≃0.3  arcmin. Unfortunately, the cusp region is smaller than the size of the innermost bin of a stellar surface density profile shown in fig. 3 of <cit.>, therefore it might be still challenging for the current observations to detect a stellar cusp in Draco II even if it is present. In addition, the observed number of stars belonging to Draco II may be too small to excavate a stellar cusp in the current observations. Although DM density profiles of UFDs may differ from one another even if they have the same sizes and LOSV dispersions, it could improve statistics for proving absence of stellar cusps to make a stacking like Fig. <ref> among UFDs having similar and sufficiently small R_ h and σ_ h. Because of the faintness of UFDs as small as Draco II, observations are limited to the close distance from the solar system: d30  kpc. It could be expected, however, that future observations will explore vaster regions to discover such faint UFDs. 



§ ACKNOWLEDGMENTS

The author thanks the referee for his/her useful comments that helped improve
the article greatly, and Takayuki R. Saitoh for kindly providing the simulation code ASURA. This study was supported by World Premier International Research Center Initiative (WPI), MEXT, Japan and CREST, JST. The numerical computations presented in this paper were carried out on Cray XC30 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.




[Antonini & MerrittAntonini &
  Merritt2012]am:12
Antonini F.,  Merritt D.,  2012, , 745, 83

[Arca-Sedda &
  Capuzzo-DolcettaArca-Sedda & Capuzzo-Dolcetta2014]ac:14
Arca-Sedda M.,  Capuzzo-Dolcetta R.,  2014, , 785, 51

[Arca-Sedda &
  Capuzzo-DolcettaArca-Sedda & Capuzzo-Dolcetta2017]ac:17
Arca-Sedda M.,  Capuzzo-Dolcetta R.,  2017, , 464, 3060

[Belokurov et al.,Belokurov
  et al.2009]bel:09
Belokurov V.,  et al., 2009, , 397, 1748

[Binney & TremaineBinney &
  Tremaine2008]bt:08
Binney J.,  Tremaine S.,  2008, Galactic Dynamics Second Edition.
Princeton Univ. Press, Princeton

[Bontekoe & van AlbadaBontekoe & van
  Albada1987]bva:87
Bontekoe T. R.,  van Albada T. S.,  1987, , 224, 349

[ChandrasekharChandrasekhar1943]c:43
Chandrasekhar S.,  1943, ApJ, 97, 255

[Cole, Dehnen & WilkinsonCole
  et al.2011]cdw:11
Cole D. R.,  Dehnen W.,    Wilkinson M. I.,  2011, , 416, 1118

[de Blokde Blok2010]d:10
de Blok W. J. G.,  2010, Adv. Astron., 2010, 789293

[DehnenDehnen1993]d:93
Dehnen W.,  1993, , 265, 250

[Di Cintio, Brook, Dutton, Macciò,
  Obreja & DekelDi Cintio et al.2017]dbd:16
Di Cintio A.,  Brook C. B.,  Dutton A. A.,  Macciò A. V.,  Obreja
  A.,    Dekel A.,  2017, , 466, L1

[Donato et al.,Donato
  et al.2009]d:09
Donato F.,  et al., 2009, , 397, 1169

[Dosopoulou & AntoniniDosopoulou &
  Antonini2016]da:16
Dosopoulou F.,  Antonini F.,  2016, preprint (astro-ph/1611.06573)

[Dubinski & CarlbergDubinski &
  Carlberg1991]dc:91
Dubinski J.,  Carlberg R. G.,  1991, , 378, 496

[El-Badry, Wetzel, Geha, Hopkins,
  Kereš, Chan & Faucher-GiguèreEl-Badry
  et al.2016]ewg:16
El-Badry K.,  Wetzel A.,  Geha M.,  Hopkins P. F.,  Kereš D.,
  Chan T. K.,    Faucher-Giguère C.-A.,  2016, , 820, 131

[Gilmore, Wilkinson, Wyse, Kleyna, Koch, Evans
  & GrebelGilmore et al.2007]gww:07
Gilmore G.,  Wilkinson M. I.,  Wyse R. F. G.,  Kleyna J. T.,  Koch A.,  Evans
  N. W.,    Grebel E. K.,  2007, ApJ, 663, 948

[Goerdt, Moore, Read &
  StadelGoerdt et al.2010]gmr:10
Goerdt T.,  Moore B.,  Read J. I.,    Stadel J.,  2010, , 725, 1707

[Goerdt, Moore, Read, Stadel & ZempGoerdt
  et al.2006]gmr:06
Goerdt T.,  Moore B.,  Read J. I.,  Stadel J.,    Zemp M.,  2006, MNRAS, 368,
  1073

[Governato et al.,Governato
  et al.2010]g:10
Governato F.,  et al., 2010, Nat, 463, 203

[Gradshteyn, Ryzhik, Jeffrey &
  ZwillingerGradshteyn et al.2007]grj:07
Gradshteyn I. S.,  Ryzhik I. M.,  Jeffrey A.,    Zwillinger D.,  2007,
  Table of Integrals, Series, and Products

[Hayashi & ChibaHayashi &
  Chiba2012]hc:12
Hayashi K.,  Chiba M.,  2012, , 755, 145

[Hayashi & ChibaHayashi &
  Chiba2015]hc:15
Hayashi K.,  Chiba M.,  2015, , 803, L11

[HernandezHernandez2016]h:16
Hernandez X.,  2016, , 462, 2734

[Hernandez & GilmoreHernandez &
  Gilmore1998]hg:98
Hernandez X.,  Gilmore G.,  1998, , 297, 517

[HernquistHernquist1990]h:90
Hernquist L.,  1990, , 356, 359

[Homma et al.,Homma
  et al.2016]hc:16
Homma D.,  et al., 2016, , 832, 21

[InoueInoue2009]i:09
Inoue S.,  2009, MNRAS, 397, 709

[InoueInoue2011]i:11
Inoue S.,  2011, , 416, 1181

[Inoue & SaitohInoue &
  Saitoh2011]is:11
Inoue S.,  Saitoh T. R.,  2011, , 418, 2527

[Ishiyama et al.,Ishiyama
  et al.2013]imp:11
Ishiyama T.,  et al., 2013, , 767, 146

[Just, Khan, Berczik, Ernst &
  SpurzemJust et al.2011]jkb:11
Just A.,  Khan F. M.,  Berczik P.,  Ernst A.,    Spurzem R.,  2011,
  , 411, 653

[Just & PeñarrubiaJust &
  Peñarrubia2005]jp:05
Just A.,  Peñarrubia J.,  2005, , 431, 861

[Klypin, Kravtsov, Bullock &
  PrimackKlypin et al.2001]kkb:01
Klypin A.,  Kravtsov A. V.,  Bullock J. S.,    Primack J. R.,  2001,
  ApJ, 554, 903

[Kormendy & FreemanKormendy &
  Freeman2016]kf:16
Kormendy J.,  Freeman K. C.,  2016, , 817, 84

[KroupaKroupa2002]k:02
Kroupa P.,  2002, Science, 295, 82

[Laevens et al.,Laevens
  et al.2015]l:15
Laevens B. P. M.,  et al., 2015, , 813, 44

[Lin & TremaineLin &
  Tremaine1983]lh:83
Lin D. N. C.,  Tremaine S.,  1983, , 264, 364

[MakinoMakino2004]m:04
Makino J.,  2004, PASJ, 56, 521

[MarastonMaraston2005]m:05
Maraston C.,  2005, , 362, 799

[Martin et al.,Martin
  et al.2016]m:16
Martin N. F.,  et al., 2016, , 458, L59

[Navarro, Frenk & WhiteNavarro
  et al.1997]nfw:97
Navarro J. F.,  Frenk C. S.,    White S. D. M.,  1997, ApJ, 490, 493

[Ogiya & MoriOgiya &
  Mori2014]om:14
Ogiya G.,  Mori M.,  2014, , 793, 46

[Oh, de Blok, Brinks, Walter &
  Kennicutt Jr.Oh et al.2011]obb:10
Oh S.-H.,  de Blok W. J. G.,  Brinks E.,  Walter F.,    Kennicutt Jr.
  R. C.,  2011, AJ, 141, 193

[Peñarrubia, Ludlow, Chanamé &
  WalkerPeñarrubia et al.2016]plc:16
Peñarrubia J.,  Ludlow A. D.,  Chanamé J.,    Walker M. G.,
  2016, , 461, L72

[Petts, Gualandris & ReadPetts
  et al.2015]pgr:15
Petts J. A.,  Gualandris A.,    Read J. I.,  2015, , 454, 3778

[Petts, Read & GualandrisPetts
  et al.2016]prg:16
Petts J. A.,  Read J. I.,    Gualandris A.,  2016, , 463, 858

[Pontzen & GovernatoPontzen &
  Governato2012]pg:12
Pontzen A.,  Governato F.,  2012, , 421, 3464

[Read, Goerdt, Moore, pontzen & StadalRead
  et al.2006]rgm:06
Read J. I.,  Goerdt T.,  Moore B.,  pontzen A. P.,    Stadal J.,  2006, MNRAS,
  373, 1451

[Saitoh, Daisaka, Kokubo, Makino,
  Okamoto, Tomisaka, Wada & YoshidaSaitoh et al.2008]sdk:08
Saitoh T. R.,  Daisaka H.,  Kokubo E.,  Makino J.,  Okamoto T.,
  Tomisaka K.,  Wada K.,    Yoshida N.,  2008, PASJ, 60, 667

[Saitoh, Daisaka, Kokubo, Makino,
  Okamoto, Tomisaka, Wada & YoshidaSaitoh et al.2009]sdk:09
Saitoh T. R.,  Daisaka H.,  Kokubo E.,  Makino J.,  Okamoto T.,
  Tomisaka K.,  Wada K.,    Yoshida N.,  2009, PASJ, 61, 481

[Saitoh & MakinoSaitoh &
  Makino2009]sm:09
Saitoh T. R.,  Makino J.,  2009, ApJL, 697, L99

[Saitoh & MakinoSaitoh &
  Makino2010]sm:10
Saitoh T. R.,  Makino J.,  2010, PASJ, 62, 301

[Saitoh & MakinoSaitoh &
  Makino2012]sm:12
Saitoh T. R.,  Makino J.,  2012, , 17, 76

[Saitoh & MakinoSaitoh &
  Makino2013]sm:13
Saitoh T. R.,  Makino J.,  2013, , 768, 44

[Silva, Lima, de Souza, Del Popolo, Le
  Delliou & LeeSilva et al.2016]sls:16
Silva J. M.,  Lima J. A. S.,  de Souza R. E.,  Del Popolo A.,  Le
  Delliou M.,    Lee X.-G.,  2016, , 5, 021

[Simon et al.,Simon  et al.2016]s:16
Simon J. D.,  et al., 2016, preprint (astro-ph/1610.05301)

[Spano, Marcelin, Amram, Carignan,
  Epinat & HernandezSpano et al.2008]sma:08
Spano M.,  Marcelin M.,  Amram P.,  Carignan C.,  Epinat B.,
  Hernandez O.,  2008, MNRAS, 383, 297

[Springel et al.,Springel
  et al.2008]swv:08
Springel V.,  et al., 2008, MNRAS, 391, 1685

[Tanikawa, Yoshikawa, Nitadori &
  OkamotoTanikawa et al.2013]tyn:13
Tanikawa A.,  Yoshikawa K.,  Nitadori K.,    Okamoto T.,  2013, ,
  19, 74

[Willman et al.,Willman
  et al.2005]w:05
Willman B.,  et al., 2005, , 129, 2692

[Zelnikov & KuskovZelnikov &
  Kuskov2016]zk:16
Zelnikov M. I.,  Kuskov D. S.,  2016, , 455, 3597

[Zemp, Moore, Stadel, Carollo &
  MadauZemp et al.2008]zms:08
Zemp M.,  Moore B.,  Stadel J.,  Carollo C. M.,    Madau P.,  2008,
  , 386, 1543









§ ANALYTIC SOLUTION OF INTERACTION INTENSITY


Equation (<ref>) integrates intensity of interactions with DM particles over possible relative velocities and impact parameters. <cit.> have mentioned that their formula (equation <ref>) containing J(v_ DM) requires a double integral which is quite expensive in numerical computations. For the sake of practical use of their formula, here I find the analytical solution of equation (<ref>).

By letting A≡ v_⋆^2-v_ DM^2 and B=b_ max/G^2m_⋆^2, the indefinite integral of equation (<ref>) can be obtained as,


    j(V) = ∫(1+A/V^2)log(1+BV^4) dV


    = (V-A/V)log(BV^4+1) -4V + I_ 1 +  Im(I_ 2)/√(2)B^1/4 +C,

where C is an integration constant, and 

    I_ 1(V) = (A√(B)-1)log(√(B)V^2+1-√(2)B^1/4V/√(B)V^2+1+√(2)B^1/4V)


    I_ 2(V) = (A√(B)+1)log(√(B)V^2-1-√(2)iB^1/4V/√(B)V^2-1+√(2)iB^1/4V).

Furthermore, Im(I_ 2) is transformed as follows,[Im[log(x+iy)]=arctan(y/x) <cit.>.]

    Im(I_ 2) = -2(A√(B)+1)arctan(√(2)B^1/4V/√(B)V^2-1).

Eventually, the definite integral, equation (<ref>), is

    J(v_ DM) = j(v_⋆+v_ DM) - j(|v_⋆-v_ DM|)

Thus, since equation (<ref>) can be solved analytically, the equation (<ref>) actually does not require a double integration but a single integration in the numerical computations.



§ COLLISIONLESSNESS OF THE SIMULATIONS


In my N-body simulations, every single star is resolved with a point-mass particle, and I should ascertain whether the gravitational interactions between the stellar particles are collisional or collisionless. Perturbations on a star by the others can be approximated as,

    Δ v_⊥^2 ≃ 8N_⋆(Gm_⋆/r_⋆ v_⋆)^2lnΛ_⋆,

where 

    Λ_⋆≡b_⋆, max/b_⋆, min∼v_⋆^2r_⋆/Gm_⋆,

where b_⋆, min∼ Gm_⋆/v_⋆^2,  b_⋆, max∼ r_⋆. Here, using a stellar mass fraction f_⋆, v_⋆^2∼ GN_⋆ m_⋆/(r_⋆ f_⋆), and the number of crossings that are required for dynamical relaxation is  

    n_ relax≡v_⋆^2/Δ v_⊥^2=N_⋆/f_⋆^2/8ln(N_⋆/f_⋆).

f_⋆≃0.03 within r_⋆ in both of the cuspy and the cored haloes. In the initial settings of my simulations, n_ relax=5.8×10^5. The crossing time is t_ cross∼ r_⋆/v_⋆≃ r_⋆/σ_ h=13.0  Myr. Eventually, I estimate the relaxation timescales to be approximately t_ relax≡ n_ relax× t_ cross∼10^3  Gyr which are significantly longer then the age of the Universe: ∼10  Gyr. Therefore, the stellar interactions in my N-body simulations are regarded to be collisionless.

