Assessing the Fornax globular cluster timing problem in different models of dark matter
KKCL-2021-07DESY 21-023
Assessing the Fornax globular cluster timing problem in different models of darkmatter
Nitsan Bar, Diego Blas, Kfir Blum, and Hyungjin Kim
3, 1 Department of Particle Physics and Astrophysics,Weizmann Institute of Science, Rehovot 7610001, Israel Theoretical Particle Physics and Cosmology Group, Department of Physics,King’s College London, Strand, London WC2R 2LS, UK DESY, Notkestrasse 85, 22607 Hamburg, Germany (Dated: February 24, 2021)We investigate what the orbits of globular clusters (GCs) in the Fornax dwarf spheroidal (dSph)galaxy can teach us about dark matter (DM). This problem was recently studied for ultralightdark matter (ULDM). We consider two additional models: (i) fermionic degenerate dark matter(DDM), where Pauli blocking should be taken into account in the dynamical friction computation;and (ii) self-interacting dark matter (SIDM). We give a simple and direct Fokker-Planck derivationof dynamical friction, new in the case of DDM and reproducing previous results in the literaturefor ULDM and cold DM. ULDM, DDM and SIDM were considered in the past as leading to coresin dSphs, a feature that acts to suppress dynamical friction and prolong GC orbits. For DDM wederive a version of the cosmological free streaming limit that is independent of the DM productionmechanism, finding that DDM cannot produce an appreciable core in Fornax without violating Ly- α limits. If the Ly- α limit is discounted for some reason, then stellar kinematics data does allowa DDM core which could prolong GC orbits. For SIDM we find that significant prolongation ofGC orbits could be obtained for values of the self-interaction cross section considered in previousworks. In addition to reassessing the inspiral time using updated observational data, we give a newperspective on the so-called GC timing problem, demonstrating that for a cuspy cold DM profiledynamical friction predicts a z = 0 radial distribution for the innermost GCs that is independentof initial conditions. The observed orbits of Fornax GCs are consistent with this expectation witha mild apparent fine-tuning at the level of ∼ CONTENTS
I. Introduction 1II. Dynamical friction: microphysics 3A. Classical gas 3B. Degenerate Fermi gas 3C. Bose gas 4III. Dynamical friction in a CDM halo: core vs.cusp 4IV. Degenerate dark matter (DDM) 6A. Stellar LOSVD constraints on DDM inFornax 7B. Structure formation constraints on DDM 8C. Orbital decay time in DDM 9V. Self-interacting dark matter (SIDM) 10VI. Discussion 11A. Comments regarding GC data used in thiswork (GC-by-GC discussion) 11B. Statistical discussion 13VII. Summary 14Acknowledgments 15 A. Dynamical friction in exotic media: derivationfrom the Boltzman equation 151. A classical gas medium 162. Degenerate fermionic dark matter 163. Ultralight dark matter 17B. Maximum entropy DDM halos 17C. Jeans modeling 19D. CDM velocity dispersion in a cored profile 19E. Orbits under dynamical friction 19F. The radial and projected CDF of GCs 201. CDF of GCs in a cuspy halo 202. CDF of GCs in a cored halo 213. Accounting for distribution of GC masses 214. Projected radius distribution 21G. Exploration of initial conditions 22References 24
I. INTRODUCTION
The Milky Way dwarf spheroidal (dSph) satellitegalaxies are broadly believed to be dominated by dark a r X i v : . [ a s t r o - ph . GA ] F e b matter (DM) [1, 2] , and this fact combined with theirsmall sizes and near-by locations makes them interest-ing test beds of the small-scale behaviour of DM [4–10].In fact, some of the basic predictions of the most com-monly considered paradigm of DM — collisionless colddark matter (CDM) — may be in tension with observa-tions (see, e.g., [11]). Conclusive kinematic data for adecisive test of CDM in dSphs is difficult to obtain, butupcoming observatories may supply it [12].One intriguing puzzle about the dSph galaxies con-cerns the globular clusters (GCs) of the Fornax dSph[13]: some of Fornax’s six known GCs [2, 14] have orbitaldecay times due to dynamical friction (DF) which seemto fall significantly short of their age [13]. If estimatedna¨ıvely based on the Chandrasekhar formula [15], assum-ing the usual CDM cusp density profile (see, e.g. [16]),one obtains an instantaneous DF time of less than 1 Gyrfor the most troublesome GC4. On the other hand, thestellar content of the GCs is old, >
10 Gyr [17, 18], asis much of the stellar content of Fornax itself [19, 20].It may seem unlikely then, that we observe some of theGCs just a short time before they fall to the center of thegalaxy. We show a visualization in Fig. 1.
NFW GC4 - - - - - - [ Gyr ] r [ k p c ] FIG. 1. Orbital radius vs. time, calculated for the FornaxGC4 assuming a slightly eccentric orbit. t = 0 represents to-day. The orbit calculation assumes the CDM Navarro-Frenk-White (NFW) [21] profile of [16]. Part of the scope of this work is to give an analyticalperspective on DF, allowing us to sharpen the GC tim-ing puzzle. When the dust settles (in Sec. VI B, usingtools developed throughout the paper) we obtain reason-ably robust predictions for the late-time distribution ofGCs: for a cuspy halo, the cumulative number of GCscontained within radius r has the form F ∆ t ( r ) ∝ τ ( r )∆ t , (1)where τ ( r ) is the instantaneous DF time (defined pre-cisely in Sec. III), ∆ t is the age of the system, and the See, however, a contrary claim in [3]. pre-factor is proportional to the initial number of GCscontained inside r ∼ τ (cid:28) ∆ t is of order τ / ∆ t .We should say in advance that although the timingpuzzle is very interesting, our analysis suggests that thepossible tension it entails is not very severe. For a cuspyCDM halo, when one takes into account projection ef-fects and the fact that Fornax hosts not just one, buta collection of GCs, then the timing puzzle may be as-cribed to a mild (but quite persistent) chance fluctuationwith a probability of 25% or so. The lack of a ∼ M (cid:12) nuclear star cluster in Fornax, the remnant of old tidallydisrupted GCs, may exacerbate the tension, but is noteasy to assess [13] (and we do not attempt to analysethis issue here).Many explanations were suggested for the GC timingpuzzle [22–35], of which an exciting class of ideas entails amodification to the nature of dark matter, going beyondCDM [29, 32–35]. In particular, Refs. [29, 34, 35] stud-ied ultralight dark matter (ULDM) and showed that inthe particle mass window m (cid:46) − eV, ULDM wouldsuppress DF enough to eliminate the timing puzzle. How-ever, most of this mass range for ULDM has been scru-tinised in the last few years, resulting with disfavouringevidence [36–39]. Motivated by the fact that the com-bination of GC age and orbit measurements probes thedetails of the DM halo and microphysics, we extend theDF analysis to additional DM models. The first model isdegenerate dark matter (DDM), in which the phase spacedistribution of DM in dSph cores is affected by Pauliblocking [40, 41]. The second model is self-interactingDM (SIDM), in which self interactions between DM par-ticles produce a cored isothermal distribution.In Sec. II we focus on the microphysics and calculateDF for CDM, DDM, and ULDM. Our results for DDMare new; for ULDM, we make contact with a differentderivation in the literature; while for SIDM the micro-physics of the DF calculation is argued to be similar tothat in CDM.All three DM models can, in principle, naturally pro-duce cored isothermal halos. As we show in Sec. III,a cored isothermal distribution of DM suppresses DFin part due to a phase-space effect (associated with the“core stalling” [42] identified in past numerical work), asthe velocity of the inspiraling GC can become paramet-rically lower than the DM velocity dispersion.For ULDM, DF and the Fornax GC timing puzzle werestudied in recent works [29, 34, 35] and we do not reviewthem again. As noted above, constraints from galaxydynamics and from cosmological Ly- α analyses suggest asimilar behaviour to CDM.For DDM (Sec. IV), we formulate a robust version ofthe Ly- α bound that is insensitive to DM model buildingand cosmological history, finding that it disfavours an ap-preciable core. If one chooses to discount the Ly- α bound(see, e.g. [29] for a qualitative discussion of concerns re-garding systematic uncertainties), then stellar kinemat-ics does allow a considerable DDM core which could leadto significant suppression of DF and prolong the settlingtime of the innermost GCs.For SIDM (Sec. V), stellar kinematics allows a con-siderable core. If the SIDM cross section is as large asthat considered in [43], then the DF settling time for theinnermost GCs can be significantly longer than in thecuspy halo CDM model.Our approach is mostly analytical. Of course, this haslimitations and one may be justified in expecting thatmore progress would require numerical simulations. Inmost current simulations, however, both GC formationand DF correspond to sub-grid phenomena. Accordingto such recent simulation work targeting the GC tim-ing puzzle [44], reasonable initial conditions for the GCs(derived from the simulations) can lead to the observedconfiguration in a standard cuspy CDM halo. The tim-ing puzzle may thus be even less significant than the mild25% that we find with analytical tools.We summarize in Sec. VII. Many details of the calcu-lations are deferred to the Appendices. II. DYNAMICAL FRICTION: MICROPHYSICS
Dynamical friction can be described in terms of theFokker-Planck theory for the motion of a probe particle(a GC in our case) traveling through a gas of spectatorparticles (DM particles in our case). In App. A we de-rive the Fokker-Planck equation as the small-momentum-exchange limit of the Boltzmann equation, governing themotion of a probe object in different background media,accounting for the gravitational interaction between theprobe and the medium particles. Our calculation is di-rect, in the sense that it simply amounts to computingthe collision integral while taking care to account for thequantum statistics of spectator gas particles. Here webypass the details of the calculation, while utilising themain results.The Fokker-Planck equation is characterised by a setof momentum space diffusion coefficients, calculated inApp. A 1, A 2, and A 3 for the case of a medium com-posed of a classical gas, degenerate Fermi gas, and Bosegas, respectively. Of particular importance for our anal-ysis is the diffusion coefficient D || , corresponding to thediffusion in momentum parallel to the probe object’s in-stantaneous velocity. The DF deceleration acting on aprobe with mass m (cid:63) moving with instantaneous velocity V w.r.t. the medium is computed as [45] d V dt = D || m (cid:63) ˆ V = − πG m (cid:63) ρV C V . (2)In the second line, to compare the DF arising in different types of medium we define the dimensionless coefficient C as follows [29], C = − V D || πG m (cid:63) ρ , (3)where ρ is the mass density of the medium.Different microphysics properties of the medium (inour case, the DM galactic halo) predict different resultsfor C . In the next subsection we discuss three scenarios. A. Classical gas
This is the appropriate limit for a halo composed ofa gas of classical particles. We will adopt this limit todescribe DF in the ordinary CDM model, as well as forthe SIDM model . For a homogeneous classical gas withan isotropic distribution function f v ( v ), DF is describedby the Chandrasekhar formula [15] (see also App. A 1), C class = 4 π ln Λ (cid:90) V dv m v m f v ( v m ) , (4)where ln Λ is the Coulomb logarithm. If the gas distri-bution function is a Maxwellian with velocity dispersion σ , f v ( v ) = (2 πσ ) − / exp( − v / (2 σ )), we have C Max = ln Λ (cid:18) erf( X ) − X √ π e − X (cid:19) (5) → ln Λ (cid:40) V (cid:29) σ √ √ π V σ V (cid:28) σ , where X ≡ V / ( √ σ ) and where in the second line weshow the asymptotic scaling of C at large and small X . B. Degenerate Fermi gas
This is the relevant limit for DF at the core of a halosupported by the degeneracy pressure of light fermionicDM (DDM model [40, 41]). In the high-degeneracy limitwe have f v ( v ) = 3 / (4 πv F ) θ ( v F − v ), where θ ( x ) is theHeaviside function, the Fermi velocity v F is related tothe medium density via ρ = gm v F π , (6) m is the mass of the particles and g the number of de-grees of freedom (e.g. g = 2 for Weyl fermions). The This is a good approximation for the SIDM cross-sections of in-terest, which are small enough such that SIDM particles travelacross distances larger than the size of the system without col-liding with each other. See Sec. V. calculation in App. A 2 gives the following limiting be-haviour, C DDM → ln Λ (cid:40) V (cid:29) v FV v F V (cid:28) v F . (7)Thus, in both limits V (cid:29) v F and v (cid:28) v F , we findthat DF in a degenerate medium is equivalent to DF in aclassical medium with the replacement σ → (cid:0) π (cid:1) v F ≈ . v F . Note that the 3-dimensional velocity dispersionassociated with the classical isotropic Maxwellian distri-bution is (cid:104) v x + v y + v z (cid:105) = (cid:104) v (cid:105) = 3 σ , while the disper-sion for the degenerate distribution is (cid:104) v (cid:105) = (3 / v F .Therefore, the pressure in the different types of mediummatches when v F ≈ . σ . Similarly, Eqs. (7) and (5) tellus that DF in these media match when v F ≈ . σ .As an aside, it is interesting to note that to leading or-der in m/m (cid:63) , the diffusion coefficient of a classical gas hasthe same functional form with respect to the distributionfunction as the diffusion coefficient of a degenerate gas(c.f. Eq. (A12) and Eq. (A21)). This is somewhat sur-prising, because the Fokker-Planck calculation took intoaccount Pauli exclusion in the medium whereas Eq. (4)does not. Moreover, according to Eq. (4), only particleswith velocities smaller than the probe object’s contributeto the DF. For the case of degenerate matter, one couldhave expected that the opposite should happen: onlyparticles close to the Fermi surface contribute to DF. Werefer the reader again to App. A 2 for the detailed com-putation that leads us to Eq. (7).Finally, note that above we evaluated DF in the zero-temperature limit and not in the finite-temperature limit.In Sec. IV we consider a finite-temperature density pro-file, so we should keep this caveat in mind. We havenot explored DF of degenerate matter within the moresophisticated treatment of Refs. [46, 47]. C. Bose gas
This is the relevant limit for the case where halo parti-cles follow the Bose-Einstein statistics, as in the ULDMmodel. The diffusion coefficients can be obtained eitherby solving a Langevin equation with stochastic fluctua-tions of the gravitational potential [34] or, as we do inApp. A 3, by using a kinetic equation . Both approachesprovide identical results.Up to a slight modification of Coulomb logarithm, DFfor the bosonic gas includes a contribution to the C termthat is identical to that of the classical gas in Eq. (5). Inaddition to this, ULDM large-scale density fluctuations(manifested by Bose-enhancement terms in the kinetic While this paper was being prepared for publication, Ref. [48]appeared which also presents a kinetic theory derivation of theULDM diffusion coefficients. theory computation) cause additional velocity drift thatcan be characterised by an extra term to C → C + ∆ C ,with ∆ C = ln Λ (cid:18) m eff m (cid:63) (cid:19) (cid:18) erf( X eff ) − X eff √ π e − X (cid:19) , (8)where m eff = π / ρ/ ( mσ ) is the ULDM massenclosed in an effective de Broglie volume and X eff ≡ v/ √ σ eff with σ eff = σ/ √
2. Nu-merically, m eff ≈ . × (cid:0) − eV /m (cid:1) [ ρ/ (3 × M (cid:12) / kpc )][(10 km / s) /σ ] M (cid:12) . With these numbersand keeping in mind a typical GC mass m ∗ ∼ M (cid:12) ,the ∆ C effect becomes quantitatively important in For-nax for m (cid:46) × − eV.The kinetic theory result summarised above assumedthat the scale size of the system – e.g., the radius r of aGC orbit – is much larger than the effective de Brogliewavelength of the ULDM particles, r dB ≈ πmσ ≈ (cid:18)
10 km / s σ (cid:19) (cid:18) − eV m (cid:19) pc , (9)and thus much larger than ULDM quasi-particle excita-tions or than the soliton core that is ubiquitously foundin ULDM simulations (see [29] for a review). For r < r dB ,the treatment above breaks down and must be modifiedby taking into account large-scale coherence effects of theULDM. This can be done via solving the Schroedingerequation, as shown in [29, 35], which indeed found thatDF becomes suppressed at r (cid:46) r dB . We refer the readerto [29, 35] for more details on DF and the Fornax GCpuzzle in the context of ULDM. Here we only note thatfor m (cid:38) − eV, where r (cid:29) r dB and m eff (cid:28) m ∗ for theFornax GCs, DF in the ULDM medium becomes quanti-tatively similar to DF in a classical medium. III. DYNAMICAL FRICTION IN A CDM HALO:CORE VS. CUSP
It is natural to define an instantaneous DF time, τ , via τ = V πG m (cid:63) ρC , (10)such that (including here only the DF effect)˙ V = − τ V . (11)A crude estimate of the time scale it would take a GCto settle down to the dynamical center of a halo canbe obtained by computing τ , using the current instanta-neous position and velocity of the GC. Assuming a CDM Formally, the ∆ C term is there also for standard CDM but isnegligible unless the individual DM particles are extremely mas-sive. NFW distribution, and plugging an estimate of the darkmatter density and velocity dispersion corresponding tothe present observed position of each GC into Eqs. (5)and (10), the result we find is summarised in the columnmarked τ CDM (highlighted in blue) in Tab. I. For GC3and GC4 the DF time estimated in this way is 2.6 and0.9 Gyr, respectively: much shorter than the age of thesystem. However, estimating an orbital decay time from theinstantaneous value of τ can be misleading. In a real-istic galaxy, the DM phase space distribution and withit the instantaneous value of τ could change along theorbit of the GC. To obtain a better estimate of the ac-tual settling time one could track the orbit of the GCsemi-analytically, using the phase space-dependent valueof τ along the orbit [28, 42, 55]. Some details of thiscalculation are given in App. E.The semi-analytic integration reproduces results fromN-body simulations [16]. To demonstrate this, we useEqs. (5) and (10) while reading the CDM density andvelocity dispersion from the N-body simulations of [16]to integrate the orbit of a GC. In Fig. 2 we compare ourresults to two different scenarios from [16].The first scenario, denoted NFW, contains a cuspyNFW-like halo (the density profile of this model is shownin Fig. 11). The orbit of a GC in this halo is shown bythe blue dashed line for the simulation of [16] and by ablue solid line for the semi-analytic tracking. The secondscenario, denoted ISO, contains an isothermal core halo(also shown in Fig. 11). The GC orbits are shown bythe red lines. Again, the semi-analytic method (solid)compares reasonably well with the simulation (dotted).The results have a mild dependence on the choice ofthe Coulomb log, and we make slightly different choicesfor the different scenarios. For NFW we follow [25] insetting ln Λ NFW = ln b max σ Gm (cid:63) . (12)However, instead of the b max = 0 .
25 kpc used in [25], weadopt b max = 0 . ISO = ln 2 V rGm (cid:63) . (13)We have checked that changing the definition of theCoulomb log according to different prescriptions in theliterature changes the predicted infall time of GCs at thelevel of a few tens of percent. This would not be crucialfor our main results.We can gain some insight on the difference between theDF settling time in the cusp vs. the core profiles. In thecentral part of a cuspy NFW halo, the density scales as ρ ∝ /r and the circular velocity scales as V circ ∝ r / . Our estimates are larger than those previously obtained in [29];we discuss the differences in Sec. VI.
NFWISO Semi - AnalyticSimulation [ Gyr ] r [ k p c ] FIG. 2. Radius of an infalling GC with mass m (cid:63) = 3 × M (cid:12) , based on the simulations of [16]. In dotted blue(thick dotted red) we plot the simulation result (Fig. 3 in[16]) for the NFW (ISO) halo. In solid lines, we plot oursemi-analytic integration. Horizontal dot-dashed lines showthe radii r f where M halo ( r f ) = m (cid:63) , in which the semi-analytictreatment should break down. Let us simplify matters by assuming (as was often donein previous works) that the GC moves on an approxi-mately circular orbit, V = V circ . With this, consider-ing the NFW halo of [16] and using Eq. (5), we find C Max ≈ . r/ kpc) . ln Λ. The DF time τ defined inEq. (10) then scales as τ ∝ r . This is a rough estimate:if we use the simulation data of [16] for σ and ρ we find asimilar but slightly different scaling, τ ∝ r . , plotted insolid blue in Fig. 4. The important point is the approx-imately power law decline of τ towards small r . This isthe cause of the fast orbital decay of the GC in the cuspyhalo model.The situation is different in a cored halo. In a core, thedensity ρ ≈ ρ = const., the circular velocity V circ ∝ r ,while a Jeans analysis shows that for an isotropic ve-locity distribution the velocity dispersion is constant σ ≈ √ Gρ r c [42], where r c is the core radius (seeApp. D). This implies V circ / ( √ σ ) ∼ r/r c , as corrobo-rated in Fig. 3 by comparing to the simulation data from[16]. At r < r c the low-velocity approximation in Eq. (5)gives C ∝ ( r/r c ) ln Λ. The ( r/r c ) factor can be thoughtof as a phase-space suppression of DF: it arises from thefactor (cid:82) V dv m v m f v ( v m ) in Eq. (4), because the veloc-ity dispersion inside an isotropic core is greater than thecircular velocity (which we assumed to match the instan-taneous GC velocity).Altogether, referring to Eq. (10), an isotropic core pre- Ref. [29] assumed that V = V circ and also took V circ equal tothe velocity dispersion σ , which would lead to a constant C Max .While this is roughly correct, for the NFW halo of [16] we findmild radial dependence of V circ /σ ∝ r . , as shown in Fig. 3. TABLE I. Some details of Fornax GCs. For the galactic center of Fornax we use an updated measurement [20], based on surfacebrightness modelling. This estimate is ≈
160 pc off relative to the center defined by previous works [2, 16, 18, 29, 31, 44],leading to different projected radii of GCs. We set the distance to Fornax as 147 ± r ⊥ bypropagating the distance error, added in quadrature with a 13 pc [20] uncertainty on the center. For relative radial velocities∆ v r , we use the galactic radial velocity RV Fornax = 55 . ± .
63 km/s [49] and set ∆ v r = RV GC − RV Fornax , adding errors inquadrature. For GC6, the values correspond to a small sample of stars, likely contaminated by background [14]. r c/h refers toKing radius for GC1-GC5 and half-light radius for GC6. The CDM instantaneous DF time (Eq. (10)) estimates are based onthe NFW profile of [16]. The instantaneous DF time of DDM and SIDM are based on Secs. IV and V. m (cid:63) [10 M (cid:12) ] r ⊥ [kpc] ∆ v r [km / s] r c/h [pc] Refs. τ CDM [Gyr] τ (135)DDM [Gyr] τ SIDM [Gyr]GC1 0 . ± .
10 1 . ± .
05 3 . ± .
18 10 . ± . . . ± .
28 0 . ± .
03 3 . ± . . ± . . .
12 8 . . ± .
84 0 . ± .
02 4 . ± .
66 1 . ± . .
63 1 .
48 2 . . ± .
15 0 . ± . − . ± .
64 1 . ± . .
91 10 . . . ± .
24 1 . ± .
05 3 . ± .
77 1 . ± . . . ∼ .
29 0 . ± . − . ± .
36 12 . ± . .
45 16 . dicts an approximately constant τ (see also [22]), τ ≈ √ π √ σ πG m (cid:63) ρ ln Λ (14) ≈ π √ r c √ ρ √ Gm (cid:63) ln Λ ≈ .
95 4ln Λ (cid:18) r c (cid:19) × M (cid:12) m (cid:63) (cid:32) ρ × M (cid:12) kpc (cid:33) Gyr . In the second line we used Eq. (D4) and in the thirdline we used values relevant for Fornax GCs. Again, wecan compare this estimate to numerical simulations. Thedashed red line in Fig. 4 shows τ as calculated by usingthe velocity dispersion and density read off the cored ISOmodel of [16]. τ is approximately constant in the core,exceeding the value of τ found for the NFW halo. NFW ISO1.2 r [ kpc ] V c i r c / σ r FIG. 3. The ratio of circular velocity V circ to the radial veloc-ity dispersion σ r , reproduced from [16] for NFW (thick red)and isothermal (ISO, blue) halos. The result that DF is suppressed in a cored halo, incomparison to a cusp, is consistent with the finding of
NFWISO [ kpc ] τ ( r ) [ G y r ] FIG. 4. Orbital decay time calculated using ρ and σ fromFig. 1 of [16] and assuming a test object on a circular orbit.NFW denotes a cuspy profile, and ISO denotes an isothermalcored profile. We used here GC mass m (cid:63) = 3 × M (cid:12) . Refs. [2, 25, 56], further confirmed in [16, 28, 42, 55]. In the next two sections we consider this perspective inexploring DM models that predict a cored halo.
IV. DEGENERATE DARK MATTER (DDM)
Ref. [40] (see also [41]) made the interesting obser-vation that light fermionic DM would produce a core indSphs, if the DM particle mass m is light enough to placethe halo in the degenerate regime. We call this model de-generate DM (DDM). The DDM core scale radius r c can The core stalling observed in N-body simulations was initiallyascribed in [56] to a failure of the Chandrasekhar formula. How-ever, as we explained here (see also [42]), semi-analytic trackingusing the Chandrasekhar formula along the orbit reproduces thisresult. be estimated via r c = AG ρ ( g m ) (15) ≈ (cid:18) ρ M (cid:12) / kpc (cid:19) − (cid:18) g m × (120 eV) (cid:19) − pc , where ρ is the core central density and where A = (cid:0) π/ (cid:1) / ≈ .
78. In App. B we give a derivation ofEq. (15), modelling the dSph halo by a maximum entropyconfiguration (at fixed total mass and energy, similarlyto [57]) . The maximum entropy halo is isothermal, scal-ing as ρ ∝ /r at large r . Between the degenerate coreand the 1 /r regime there are intermediate features thatdepend on the temperature.Inside r (cid:46) r c DDM particles are described by a degen-erate distribution function with Fermi velocity v F relatedto their mass density via Eq. (6). DF for this system ischaracterised by Eq. (7), so inside the DDM core, where ρ and v F are constant, Eq. (10) yields a constant DFtime, τ DDM ≈ π G gm m (cid:63) ln Λ (16)= 4 . M (cid:12) m (cid:63) × (150 eV) gm Gyr . Note that if one inserts the DDM halo core radiusEq. (15) into Eq. (14), one obtains the same paramet-ric dependence as in Eq. (16). This is a result of thesimilarity between Eqs. (5) and (7). The interesting fea-ture of the DDM model is that it produces the core dueto Pauli blocking.Naively, Eq. (16) suggests that DF in a DDM corecould be arbitrarily suppressed by decreasing m . Thishappens because decreasing m at fixed ρ is tied to in-creasing v F . The DF effect on Fornax GCs is thus aninteresting test bed of DDM, and in Sec. IV C we explorethis point in more detail. Before entering that discussion,however, we first consider observational limits on m .First, the Fermi velocity cannot be arbitrarily high ina gravitationally bound halo [58]. In Sec. IV A we makethis analysis more precise by fitting stellar line-of-sightvelocity distribution (LOSVD) data to the DDM halomodel; we find that while the Fornax LOSVD data in-deed constrains m (cid:38)
100 eV or so, this constraint byitself would still allow a significant modification of DFcompared to the CDM prediction.A second and much tighter constraint comes from cos-mological structure formation as observed through Ly- α forest statistics. We show in Sec. IV B that this con-straint directly affects the same combination, gm , that Our approach in App. B is similar to that of [40], but differsin that we include also non-zero temperature solutions. Suchsolutions were noted but not implemented in [40]. We find thatthese solutions could expand the range of applicability of theDDM model in dSphs. appears in Eqs. (16) and (15). Imposing the Ly- α con-straint excludes DDM from making an appreciable corein Fornax on scales r (cid:38)
100 pc, meaning that DDMcould not significantly affect the orbits of GCs. Whileearlier work on DDM argued that a non-thermal produc-tion mechanism for DDM could avoid the cosmologicalconstraint, we formulate a rather robust version of thebound which appears difficult to evade.
A. Stellar LOSVD constraints on DDM in Fornax
In this section we summarise the results of a Jeansanalysis for the DDM model in Fornax . The DDM pro-file is described in App. B and the details of the Jeansanalysis are given in App. C.In Fig. 5 we plot LOSVD data of Fornax [10] along-side fits of the density profile presented in App. B. Ourfitting procedure is based on a simple χ minimization,where χ ≡ (cid:80) N data i =1 ( σ LOS ,i − σ LOS ( r i )) /σ i and σ i is thereported uncertainty for radial bin i . At a given particlemass m , our fit has three free parameters: the degen-eracy parameter µ /T , the central core density ρ , andthe stellar velocity anisotropy parameter β , taken to beconstant in r . DDM, β =- =
135 eV, χ / d.o.f = [ pc ] σ L O S [ k m / s ] FIG. 5. LOSVD compared with data [10] for DDM with g =2. The best-fit parameters of the profile are µ /T = 3 and ρ = 3 . × M (cid:12) / kpc . In Fig. 6 we plot the circular velocity compared tothe Fermi velocity, which we define using v F ( r ) ≡ (cid:112) µ ( r ) /m , for the m = 135 eV fit. As can be expected,the circular velocity scales as V ∝ r , whereas v F remainsconstant at small radii. This leads to a suppression ofDF, as explained in Sec. III. For previous analyses, see [40, 41, 59–62]. V circ ( r )= GM ( r ) (cid:0) v F ( r ) [ pc ] V e l o c i t i e s [ k m / s ] FIG. 6. The circular velocity V circ and the Fermi velocity v F for the DDM halo with m = 135 eV, ρ = 3 . × M (cid:12) / kpc , µ /T = 3 and g = 2. The vertical lines show the estimatedorbital radii ( r = r ⊥ × / √
3) of the three GCs closest to thedynamical center of Fornax, c.f. Table I.
B. Structure formation constraints on DDM
The free streaming of light DM suppresses the matterpower spectrum [63], notably constrained via Ly- α foreststatistics [64–66], with details depending on the cosmo-logical DM production mechanism. Refs. [40, 41] consid-ered non-thermal mechanisms for cosmological produc-tion of DDM, aiming to bypass the Ly- α bounds. Withsuch mechanisms in mind, values of m in the ballparkof 100 eV were considered in these works. We now re-visit the cosmological bound and formulate a conserva-tive limit that is insensitive to the cosmological produc-tion mechanism of DM. Our results suggest that mech-anisms of the kind proposed in [40, 41] should not beable to produce m < . α bound.The instantaneous DM free streaming wavelength k FS depends on the DM velocity dispersion , k FS ( z ) ≡ (cid:114) H ( z ) c s ( z ) (cid:39) (cid:114) H ( z ) σ ( z ) . (17)Here z is the redshift, H ( z ) is the Hubble rate, we haveset σ ( z ) = (cid:112) (cid:104) v (cid:105) and (cid:10) v (cid:11) = (cid:82) d p p m + p f ( p ) (cid:82) d pf ( p ) . (18)Bounds on warm DM (WDM) [66] effectively constrain k FS ( z ). Specifically, they apply to z (cid:46) , where densityperturbations on comoving scales of the order of λ ≈ (1 + z ) / H ( z ) ≈ . This formula is correct up to an order unity factor relating thespeed of sound c s and the velocity dispersion σ [67]. This factoris unimportant for our analysis. We can convert the WDM limit of [66] into a boundon DDM by the following prescription. At a given en-ergy density, the coldest possible distribution functionof DDM is the fully degenerate distribution f ( p ) = θ ( p F − p ), where the Fermi momentum p F is related tothe energy density via ρ = g (2 π ) (cid:90) d p (cid:112) m + p f (19)= g π (cid:20) p F (cid:113) m + p F ( m + 2 p F ) − m sinh − (cid:16) p F m (cid:17)(cid:21) ≈ (cid:40) gmp F π p F (cid:28) m gp F π p F (cid:29) m . The p F parameter redshifts as p F ∝ (1+ z ). The velocitydispersion for this distribution is (cid:10) v (cid:11) = 1 − (cid:18) mp F (cid:19) + 3 (cid:18) mp F (cid:19) arctan (cid:16) p F m (cid:17) ≈ (cid:40) p F m p F (cid:28) m p F (cid:29) m . (20)Using Eq. (20) we can calculate k FS ( z ), compare thisto the k FS ( z ) of WDM, and cast the bounds of [66] intothe most conservative, maximally-cold DDM model bymatching the k FS ( z ) curves of the two models. To re-call, WDM was defined [66] by the distribution function f WDM = (exp( p/T ) + 1) − , where T ∝ (1 + z ). Com-paring the DM mass density for DDM and WDM in thenonrelativistic regime, we have ρ = (cid:40) gmp F π DDM ζ (3) gmT π WDM . (21)Matching the density in the two models implies:( T /m ) WDM ( p F /m ) DDM = (cid:34) ζ (3) (cid:0) gm (cid:1) DDM ( gm ) WDM (cid:35) . (22)On the other hand, still in the nonrelativistic regime wecan compare the velocity dispersion in the two models, (cid:104) v (cid:105) DDM (cid:104) v (cid:105) WDM ≈ ζ (3)75 ζ (5) (cid:20) ( p F /m ) DDM ( T /m ) WDM (cid:21) = 3 ζ (3)75 ζ (5) (cid:18) ζ (3)2 (cid:19) (cid:34) (cid:0) gm (cid:1) WDM ( gm ) DDM (cid:35) , (23)where in the second line we used Eq. (22). We illus-trate the comparison in Fig. 7. For model parametersof interest to this discussion, the ratio of velocity dis-persions is redshift-independent at z (cid:46) . We cantherefore conclude that the WDM constraint of [66], m W DM > .
96 keV at the 95%CL for g = 2, impliesthe constraint: gm DDM > × (1 . . (24) WDM 2.96 keVWDM 1 keV DDM 150 eVDDM 500 eVDDM 1450 eV - [ z ] V e l o c i t y d i s pe r s i on < v > FIG. 7. A comparison of the velocity dispersions (in naturalunits) of WDM and DDM for various particle masses. Includ-ing CMB data, Ref. [66] puts a bound of m > .
96 keV onWDM (95% CL), for which we plot the velocity dispersion ingreen solid line.
We expect that model-building around the bound ofEq. (24) would be quite difficult. No production mecha-nism should be able to create a colder distribution func-tion for DDM. In particular, the skewed momentum dis-tribution scenarios of [41] and the scalar decay models of[40, 68] should all satisfy this bound. Using Eq. (15), wefind that Eq. (24) constrains the DDM core in Fornax to r c (cid:46)
20 pc, irrelevant for the orbits of GCs.The Ly- α analyses may be affected by systematic un-certainties related, among other things, to the thermalhistory of the inter-cluster medium (IGM) and otherbaryonic effects [66]. Keeping this caveat in mind, itseems sensible to take Fig. (7) with a grain of salt. If weallow k FS of the (coldest possible) DDM model to exceedthe nominal bound of [66] by, say, a factor of ∼ .
5, wecould relax Eq. (7) to gm > × (500 eV) . With such(rather ad-hoc) relaxed bound we could allow a DDMcore r c (cid:46)
80 pc, still irrelevant for GC orbits in Fornax.Going down to m = 150 eV (still at g = 2), which wouldallow a DDM core radius of r c ≈
385 pc encompassingsome GC orbits, would amount to k FS smaller by a factorof 20 than the nominal WDM bound. C. Orbital decay time in DDM
As we have seen, Ly- α analyses exclude DDM fromproducing a core extending to the observed orbital po-sitions of GCs in Fornax. This means that DDM wouldnot change the standard CDM predictions for the DF set-tling time of the GCs. Nevertheless, given that the Ly- α bound is subject to some debate, it is interesting to seewhat DDM could do to DF subject only to the LOSVDconstraints of Sec. IV A.In order to estimate the instantaneous DF time scale τ ,we use Eq. (10), with a modified Eq. (7). In order to in-terpolate between a quasi-degenerate core and a classical gas in the outskirts of the halo, we adopt C DDM = 1 . + v F V ln Λ , (25)such that in the regime V (cid:29) v F , we retrieve CDM-likebehavior, c.f. Sec. III. For ln Λ, we adopt the choice forln Λ ISO , c.f. Eq. (13).We use the GC masses and projected radii collectedin Table I, combined with the density profile derived inthe LOSVD fits. We correct for the projection effect byrelating the assumed true orbital radius to the observedprojected radius of the GC via r true /r ⊥ = 2 / √
3. Wealso assume that the GCs are on circular orbits, setting V true /V circ ( r true ) = 1. (This de-projection procedure is,of course, simplistic: we will shortly report a more com-prehensive treatment.) The results are summarised inTable I. For m = 135 eV we find na¨ıve orbital decay time-scales of 1 .
48 Gyr and 10 . .
63 Gyr and 0 .
99 Gyr. There-fore, while the na¨ıve DF time of GC4 in DDM is muchlonger than in cuspy CDM, for GC3 the na¨ıve time incuspy CDM actually exceeds that of DDM.However, as discussed in Sec. III, the instantaneous τ can be misleading when comparing different halo mor-phologies: a realistic estimate of the GC settling timerequires orbit integration. We turn to this analysis next,finding that the real orbit settling times of both GC4 andGC3 are in fact longer in DDM compared with cuspyCDM.To obtain a more comprehensive estimate of the DFsettling time and the impact of projection effects, weuse the orbit integration explained in App. E with ini-tial conditions that we vary as follows. For each GC,we scan the range r true ∈ [1 , r ⊥ (the logic behind thisrange is explained in App. G). For each r true we scan over V true ∈ [0 . , . V circ ( r true ). For each V true we test pos-itive and negative cos θ . Finally, we test the two cases,∆ v y = (cid:112) v − ∆ v r , ∆ v z = 0 and ∆ v y = 0 , ∆ v z = (cid:112) v − ∆ v r .For each starting point in phase-space, we integratethe equations of motion, stopping the integration when( r apoenter + r pericenter ) / (cid:46) . r initial ≡ f r r initial , or af-ter 10 Gyr (the earlier of the two). We then report theintegration time as τ inspiral .In Fig. 8 we plot the result of this procedure for GC3,comparing the DDM halo for m = 135 eV (top panel)and the cuspy CDM halo from [16] (bottom panel). Forthe representative phase-space point r true /r ⊥ = 2 / √ V true /V circ ( r true ) = 1, highlighted in Fig. 8 by a red dot,we find that the inspiral time in the DDM halo is in factlonger ( ∼ ∼ . f r (set as00 . τ ( r ) ∝ r , therefore the inspiraltime is mostly sensitive to the initial radius. In a coredmodel, τ ( r ) is a weak function of r , potentially even non-monotonuous. Therefore, the definition of the inspiraltime in the cored model becomes sensitive to the radiusat which the orbit integration is stopped. . . . . . . . . . . . . . . . . . . . . . . . . . FIG. 8. Contours of the inspiral time of GC3, defined inApp. G, for DDM ( top ) and cuspy CDM model ( bottom ).The “na¨ıve” estimates written on top are those given in Ta-ble I based on an evaluation of the instantaneous DF time at r true /r ⊥ = 2 / √ V true /V circ ( r true ) = 1. The different linetypes correspond to different discrete choices in our scan ofthe initial conditions in phase-space, explained in more detailin App. G. V. SELF-INTERACTING DARK MATTER(SIDM)
Self-interacting DM (SIDM) is a simple modification ofCDM, that could arise in many models [43, 69–76]. Theself interactions can be expressed in terms of the crosssection per unit mass, σ/m , which could be velocity-dependent [74]. The scattering mean free path is l = mρσ = 48 10 M (cid:12) / kpc ρ / gr σ/m kpc , (26)and the time between scatterings l/v is t scat = 2 .
35 20 km / s v M (cid:12) / kpc ρ / gr σ/m Gyr . (27)When l is larger than the distance across the halo, weexpect that the microphysics of DF in the SIDM modelwill be similar to that of CDM. On the other hand, themorphology of an SIDM halo could be different to that inCDM as long as t scat is smaller than the age of the system.Given a large enough cross section, SIDM produces coredhalos which affect the orbital settling time of GCs asdiscussed in Sec. III.We follow Ref. [43] in modelling the SIDM halo profile.Inside some radius r , we assume an hydrostatic profilewith central density ρ c and pressure P = σ ρ . For aself-gravitating spherical halo, the density profile obeys1 r ∂ r ( r ∂ r ln ρ iso ) = − r ρ iso ρ c , (28)where r ic ≡ σ / √ πGρ c (similar to the King radius [45]).Beyond r we match the density to the NFW profile, ρ NFW = ρ s ( r/r s ) − (1 + r/r s ) − , fixing ρ s and r s by im-posing continuity of ρ iso ( r ) = ρ NFW ( r ) and of the en-closed mass M iso ( r ) = M NFW ( r ). This procedure isconsistent with an initially NFW-like cusp profile thatwas deformed into a cored isothermal profile due to theSIDM scatterings. Altogether, the halo model has threefree parameters, ρ c , r and σ , that we can constrain withLOSVD data.In Fig. 9 we plot a LOSVD fit, following the same χ procedure as in Sec. IV. The model is taken to illustratea large core solution that would include the orbit of GC3.In Fig. 10 we show the inspiral time of GC3, using thesame procedure as in Fig. 8. We use Eq. (5) with σ fromthe LOSVD fit and adopt ln Λ = ln Λ ISO as in Eq. (13).We find that the large core SIDM model significantlyincreases the inspiral time of GC3 compared to the CDMprediction.Ref. [43] pointed out that in the SIDM core region, DMparticles have undergone about a single collision duringthe age of the system, i.e. (cid:104) σv (cid:105) m ρ c t age = t age t scat ≈ . (29)1 SIDM, β χ / d.o.f =
500 1000 1500 200005101520 r [ pc ] σ L O S [ k m / s ] FIG. 9.
Top:
LOSVD data of Fornax dSph modeled bydifferent SIDM profiles. The central density is ρ c = 2 . × M (cid:12) / kpc , the velocity dispersion is σ = 17 km / sec and r /r ic = 6. With this assumption, we can estimate the cross sectionimplied by the LOSVD fit (cid:104) σv (cid:105) m ≈
18 10 Gyr t age . × M (cid:12) / kpc ρ c cm g kms . (30)This result is compatible with the baseline model of [43],which predicted (cid:104) σv (cid:105) /m ∼
25 cm g − km s − . It wouldsignificantly increase the DF settling time of the inner-most Fornax GCs. FIG. 10. DF time for SIDM halo, analogously to Fig. 8. It has been pointed out [72, 77] that this assumption may besimplistic.
VI. DISCUSSION
In Fig. 11 we summarise the key features of differ-ent models of DM discussed in this work, including bothcuspy and core halo models. In the spirit of [2], we alsoadd the hybrid coreNFW model of [78] with density ρ coreNFW = (cid:101) f n ρ NFW + n (cid:101) f n − (1 − (cid:101) f )4 πr r c M NFW , (31)where (cid:101) f = tanh( r/r c ) and M NFW = (cid:82) r d r (cid:48) ρ NFW . Thismodel aims to describe a CDM-dominated halo modifiedby baryonic feedback.The LOSVD data (top left panel of Fig. 11) is de-scribed reasonably well in all models, with the fit of theISO model of [16] being slightly worse.The instantaneous DF time for a GC with m ∗ =3 × M (cid:12) is shown in the bottom right panel of Fig. 11.It illustrates the fact that the main impact of the mi-crophysics of DM (as in DDM and SIDM) on DF comesfrom their prediction of a cored halo morphology, andnot from the exotic microphysics per-se. The formationof a core due to baryonic feedback in CDM [25, 56] couldtherefore have similar consequences.The density profiles (top right panel of Fig. 11) demon-strate the cusp for NFW, large-cores and the intermedi-ate coreNFW . We also plot an estimate of the stellar den-sity. This may become important for large-core models,whose density is only larger by a factor of two or so thanthe stellar density at small r . In these cases, accountingfor the stellar-induced potential could slightly change ourresults numerically, but not qualitatively: as far as theDF microphysics is concerned, background stars wouldcontribute to the DF of a GC just like CDM particles,and since the total mass density is constrained by theLOSVD fit, the separation into DM and stellar compo-nents is not essential for the DF computation.It is interesting to compare the circular velocity profilesof different models (bottom left panel of Fig. 11) to thephase-space parameters of GCs.In Sec. VI A we briefly consider each of the six GCs,noting the implications w.r.t. the timing puzzle. Next,in Sec. VI B we consider the GC system as a whole anduse the tools we have developed throughout this work tore-evaluate the problem. A. Comments regarding GC data used in this work(GC-by-GC discussion)
Let us re-evaluate the role of each GC in the timingpuzzle, in light of the observational data used in this work For the coreNFW model we adopt the best-fit of [10], with ρ ≈ . M (cid:12) , r s ≈ . n ≈ . r c ≈ .
52 kpc. We derivethe velocity-dispersion of the halo by solving the Jeans equation,Eq. (D1), assuming isotropy. Fornax: data ( ) and models, = const. M(cid:0)(cid:1) N(cid:2)(cid:3) (cid:4)(cid:5)(cid:6)
ISO DDMSIDM CoreNFW
L(cid:7)(cid:8)k(cid:9) s [ pc ] V c(cid:10)(cid:11)(cid:12) ( r ) [ (cid:13)(cid:14) s ] Stars5 × × × (cid:0) ×
8M (cid:1) / k(cid:2)(cid:3)3 ] (cid:4) = (cid:5) * M s(cid:6)(cid:7)(cid:8)(cid:9)(cid:10)(cid:11)(cid:12) N(cid:13)(cid:14) (cid:15)(cid:16)(cid:17) I(cid:18)(cid:19) D(cid:20)(cid:21)S(cid:22)(cid:23)(cid:24) C(cid:25)(cid:26)(cid:27)(cid:28)(cid:29)(cid:30)
200 500 1000 2000 [ pc ] G(cid:31) ] FIG. 11. Comparison of models.
Top left:
LOSVD data and fits. The M19 NFW and M19 ISO models refer to the halos of[16], for which we only fit the velocity anisotropy. The DDM and SIDM models are based on Sec. IV and Sec. V. The χ /d.o.fis ≈ . ≈ . β NFW = − . β ISO = 0 . β DDM = − . β SIDM = 0 . β coreNFW = − . Top right:
Density profiles. In addition to DM,we also plot an estimate of the stellar density, assuming a Plummer profile with scale r p = 851 pc [20] and mass 4 × M (cid:12) [17]. Bottom left:
Circular velocity. (Note how cored models require some tuning to explain the large radial velocity of GC4, | ∆ v r | = 8 . ± .
64 km/s, at its small projected radius r ⊥ ≈ .
154 pc).
Bottom right:
Instantaneous DF time, evaluated for m (cid:63) = 3 × M (cid:12) . (see Tab. I) including projected radii that are differentthan those of earlier analyses [2, 18, 29, 31] and LOSvelocities that were often ignored in past analyses, butare in fact known fairly well from observations.1. GC1 and GC5 – these GCs do not seem to pose atiming problem as they are located at fairly largeprojected radii, r ⊥ ≈ . . ÷ . ÷
30 km/s, the measured GCLOS velocities are smaller than V circ by a factor of5 or so. This could hint that GC1 and GC5 areclose to the apocenter of fairly eccentric orbits. Ifthis is indeed the case, then the na¨ıve instantaneousDF time overestimates the true orbital settling timebecause the GCs typically experience stronger DFwhen they venture into smaller radii, as expectedif the orbit is eccentric.2. GC2 – does not seem to present a timing problem.3. GC3 – for a cuspy profile, our “na¨ıve” instanta-neous DF time of 2 . τ ≈ . r ⊥ and the updated LOSVD data [10] comparedto the older data [1]. Of more relevance, however,is the actual physical inspiral time obtained withthe orbit integration of App. G. We find that botha simple cusp profile and the intermediate coreNFW profile predict an inspiral time ≈ . ÷ ≈ C ≈ . coreNFW profile predicts inspiral time of ≈ (cid:38) r ⊥ ≈
150 pc forthis GC), GC4 needs to be on a circular orbit with r true /r ⊥ (cid:38) r true /r ⊥ (cid:38) .
6. These possibilitiesare somewhat tuned, either w.r.t. to the projectionangle or w.r.t. to the orbital phase. In comparison,an NFW profile can comfortably accommodate theradial velocity of GC4. For the coreNFW model,it is marginally possible to have GC4 on a circularorbit without tuning in radius.5. GC6 – the newly rediscovered GC [14] probablyhas a smaller mass ( ≈ . × M (cid:12) [44]) than theother five GCs. It does not seem to reinforce thetiming puzzle.Ref. [14] notes that GC6 has an elongated shapeand may be undergoing tidal disruption. This maycomprise some evidence in favour of a cuspy halo.To summarise, the usual suspects for a GC timing puz-zle, GC3 and GC4, are found here to have somewhatlonger settling times than previously thought [29], butnevertheless much shorter than their age. B. Statistical discussion
Let us finally use the tools we developed to re-evaluatethe timing puzzle. Consider the cuspy NFW and thecored ISO profiles in Fig. 11. Using the results developedin App. F, we can map a distribution of GC initial radiiinto the cumulative distribution function (CDF) of GCprojected radii today.An example of such a calculation is shown in Fig. 12,with the NFW result on the left panel and the ISO onthe right. For concreteness, in making Fig. 12 we usedan initial distribution of GC radii of the form f ( r ) ∝ r exp ( − a r ), with r given in kpc and the parameter a in kpc − . We give this initial distribution ∆ t = 10 Gyrto evolve. For the NFW example we set a − = 0 . a − = 0 . f is used here mainly for illus-tration. Physically, the scaling f ∝ r at small r (cid:28) a − could arise naturally if the initial 3D distribution of GCsis constant in radius, consistent with the current stellardistribution in Fornax. The peak of the distribution is at r = 2 a − , comparable to the current half-light radius.For the NFW halo, the derivation in App. F showsthat before projection effects are taken into account thetime-evolved CDF of GC radii at r (cid:28) r cr ≈ F ∆ t ( r ) ≈ A τ ( r )∆ t , (32)independent on initial conditions, with A an O (1) coef- ficient. This small- r approximation is shown by thedashed red line in the left panel of Fig. 12. The dashedgreen line shows the full unprojected GC CDF, consistentwith Eq. (32) for r <
500 pc.Projection has a significant effect, meaning that typ-ically, a considerable number of GCs observed at pro-jected radius smaller or equal to r ⊥ are in fact located at r > r ⊥ . The solid blue line shows the projected radiusCDF: for r ⊥ (cid:46) . r ⊥ .The GC timing puzzle is reflected in the left panel ofFig. 12, by the presence of two GCs (GC4 and GC6)inside of r ⊥ <
250 pc or so, where the CDF shown by thesolid blue line predicts that no more than one GC shouldbe expected. While Fig. 12 shows the CDF resultingfrom just one example of initial conditions, we could notfind initial conditions that would fit the innermost GC4and GC6, while not at the same time overshooting GC3and GC2 further out; nor could we find initial conditionswhich fit GC3, while not at the same time undershootingGC4 and GC6. The reason for this (mild) inconsistencyis the growth in r of F ∆ t in Eq. (32), which is model-independently (that is, irrespective of initial conditions)predicted to possess a strong slope F ∆ t ( r ) ∝ τ ( r ) / ∆ t ∝ r . . We stress, however, that the inconsistency is indeedquite mild: the Poisson probability to see two or moreGCs where only one GC is expected is about 25%. Thisdoes not seem like severe fine-tuning.As the right panel of Fig. 12 shows, the cored ISO pro-file can easily provide a time-evolved projected CDF inexcellent agreement with the data. While here we per-formed the calculation for the cored ISO profile, our anal-ysis in this paper makes clear that essentially any modelof DM microphysics would lead to similar conclusions aslong as it produces a sizeable core in Fornax.Fig. 12 suggests that even for a cuspy halo (left panel),the GC timing puzzle does not invoke extreme fine-tuningin the sense that it is not difficult to find apparently rea-sonable initial conditions that evolve to the observed con-figuration of GCs. Another point to consider is relatedto GCs that inspiral down and are tidally disrupt, pre-sumably forming a stellar nucleus [13]. In the NFW casein Fig. 12, about 50% of GCs present initially in the haloarrived at the dynamical center within ∆ t . This may beexpected to produce a nucleus for Fornax with a stellarmass in the ballpark of 10 M (cid:12) , which – as far as we areaware – is not observed. However, Ref. [13] suggestedthat because of the small number of GCs involved, 3-body interactions between accreted GCs could precludethe formation of the nucleus.More insight could come from numerical simulations,even though both DF and GC formation involve sub-grid physics in most existing simulations. Recently, For an NFW halo, A ≈ . N cr where N cr specifies (approx-imately) the initial number of GCs located inside the criticalradius r cr ≈ GC4GC6 GC3 GC2 GC5GC1
NFW
500 1000 1500 2000 250001234567 r ⟂ [ pc ] CD F GC4GC6 GC3 GC2 GC5GC1
ISO
500 1000 1500 2000 250001234567 r ⟂ [ pc ] CD F FIG. 12. Example of a calculated CDF of projected radii using τ ( r ) from Fig. 11. Left:
NFW halo.
Right:
Cored ISOhalo. Solid blue line shows the CDF after projection effects are taken into account. Dashed green line shows the result beforeprojection. For the NFW case, the small- r prediction of Eq. (32) is shown by the red dashed line. Observed Fornax GCs arealso shown. The initial radial GC PDFs used to make the plot are explained in the text. Ref. [44] explored the survival of GCs in cuspy Fornax-like halos using hydrodynamical cosmological simulations[84, 85], to which sub-grid formation of GCs was addedand their orbital evolution under DF was tracked in post-processing. Ref. [44] calculated the projected radius CDFof GCs at z = 0, finding that their simulated CDF is con-sistent with the observed positions of the Fornax GCs.According to [44], only around 33% of GCs are tidallydisrupted in the simulations; a somewhat smaller num-ber than the 50% in our example in Fig. 12. Fornax isfound to be special, in that only about 3% of Fornax-like galaxies in the simulation ended up with 5 or moresurviving GCs today. VII. SUMMARY
We revisit the calculation of globular cluster (GC) or-bits under dynamical friction (DF), considering differentmicroscopic models of dark matter (DM) and differenthalo morphologies that they predict. We focus on theFornax dwarf spheroidal (dSph) galaxy, which hosts 6GCs and which has been noted in previous literature topose a GC timing problem, that is, the future orbits ofsome of its GCs are much shorter than their current age.We presented semi-analytical computations of DF andof GC orbits under DF. For a cuspy DM halo, we showedthat the current cumulative distribution function (CDF)of GC radii takes an approximately power law form thatcan be deduced from stellar kinematics and age mea-surements. Including projection effects, we demonstrated Note that [44] considered GC6 as a candidate for a tidally dis-rupted GC, assuming that it’s projected distance from the centerof Fornax is just 30 pc. In comparison, using the Fornax centerof [20], we find a projected distance of 254 pc for GC6. that the GC timing problem does not appear very severe:the existence of the innermost GCs could be accountedfor at the cost of moderate fluctuation with a Poissonprobability of about 25%. A cuspy halo, in conjunctionwith GC orbits, does place interesting constraints on theinitial distribution of GCs in dSphs, as it suggests thatan O (1) of initially present GCs must be tidally disruptby the halo, making GC formation quite prolific. Testingthese constraints would likely require high resolution nu-merical simulations including baryonic effects, where GCformation can be modelled from first principles.The recent numerical simulations of Ref. [44] made aninteresting step in this direction. According to Ref. [44]a reasonable distribution of initial conditions for the GCsmay naturally lead to the observed present configuration.However, both GC formation and DF were treated in [44]at the sub-grid and post-processing level, and it is notclear (to us) if the resolution of the simulation was highenough to resolve the inner region of Fornax containingthe innermost GCs. Moreover, theoretical insight aboutthe role that initial conditions play in shaping the presentdistribution of GCs is important. We thus believe thatour analytical approach remains useful.The fact remains that the combination of GC age andorbit measurements could probe the Fornax DM halo andmicrophysics. At the level of the microphysics, in Sec. IIwe calculated DF for three models of DM: fermionic de-generate DM (DDM), where Pauli blocking affects theDF derivation; bosonic ultralight DM (ULDM), where anastronomical de-Broglie scale comes into play; and self-interacting DM (SIDM), where – in the limit we weremostly interested in – the microphysics of DF shouldmostly follow that of CDM, but the halo morphologyis different.For ULDM, DF and specifically the Fornax GC prob-lem were studied in a number of works. Constraints fromgalaxy dynamics and cosmological Ly- α analyses exclude5a soliton core reaching out to the orbits of Fornax GCs,and lead to a similar behaviour as in CDM.For DDM, we gave a new derivation of DF. We thenformulated a robust (in terms of DM model building)version of the Ly- α bound, showing that it excludes anappreciable core, leading again to CDM-like behaviourat the scale of GC orbits. At the same time, stellar kine-matics in Fornax could still allow a considerable DDMcore. If the Ly- α bound is discounted, for some reason,then DDM could lead to significant suppression of DFand prolongation of the settling time of GC3 and GC4.For SIDM, stellar kinematics allows a considerablecore. If the SIDM cross section is as large as that consid-ered in [43], then the DF settling time for GC3 and GC4can be significantly longer than in the cuspy halo CDMmodel. ACKNOWLEDGMENTS
We thank Justin Read for clarifications on the kine-matic data and Valerie Domcke, Joshua Eby and DanielKaplan for useful discussions. NB is grateful for the sup-port of the Clore scholarship of the Clore Israel Founda-tion. KB is incumbent of the Dewey David Stone andHarry Levine career development chair at the WeizmannInstitute of Science, and was supported by grant 1784/20from the Israel Science Foundation. HK is supported bythe Deutsche Forschungsgemeinschaft under Germany’sExcellence Strategy - EXC 2121 Quantum Universe -390833306.
Appendix A: Dynamical friction in exotic media:derivation from the Boltzman equation
In this section we provide an economical derivation ofgravitational DF acting on a nonrelativistic probe objectmoving in a medium, with different medium microphysicsincluding a classical gas as well as quantum Fermi andBose gases. We neglect interactions apart from minimalgravity. We start with a quick recap of the derivation ofthe Fokker-Planck equation, governing the phase spacedistribution functions of the probe and medium particles.We consider the following elastic scattering process oftwo particle species,1( p ) + 2( k ) → p (cid:48) ) + 2( k (cid:48) ) . The phase space distribution function for the particlespecies 1 evolves according to the Boltzmann equation, df dt = C [ f ] . (A1)The collision integral C [ f ] contains information about the elastic scattering processe, and is written as C [ f ] = (2 π ) E p (cid:90) d Π k d Π p (cid:48) d Π k (cid:48) δ (4) ( p + k − p (cid:48) − k (cid:48) ) |M| × (cid:104) f ( p (cid:48) ) f ( k (cid:48) )(1 ± f ( p ))(1 ± f ( k )) − f ( p ) f ( k )(1 ± f ( p (cid:48) ))(1 ± f ( k (cid:48) )) (cid:105) , (A2)where |M| is a squared matrix element averaged overinitial and final spins, and d Π k = g E k d k (2 π ) is the Lorentzinvariant phase element with the number of internal de-grees of freedom g . The sign in 1 ± f i refers to bosons (+)or fermions ( − ), respectively. It is convenient to write theabove Boltzmann equation in the following form, df dt = (cid:90) d p (cid:48) (2 π ) (cid:104) S ( p (cid:48) , p ) f ( p (cid:48) )(1 ± f ( p )) − S ( p , p (cid:48) ) f ( p )(1 ± f ( p (cid:48) )) (cid:105) , (A3)where the function S encodes the response of themedium, and is defined as S ( p , p (cid:48) ) ≡ (2 π ) E p E p (cid:48) (cid:90) d Π k d Π k (cid:48) δ (4) ( p + k − p (cid:48) − k (cid:48) ) ×|M| f ( k )(1 ± f ( k (cid:48) )) . (A4)The function S ( p , p (cid:48) ) can be interpreted as a differentialrate at which a particle of momentum p is converted intoa particle with momentum p (cid:48) .The Boltzmann equation can be greatly simplified ifthe momentum exchange q = p (cid:48) − p (A5)is smaller than the typical momentum given by the distri-bution function f . In such cases, the Boltzmann equa-tion is reduced to the nonlinear Fokker-Planck equation, df dt = − ∂∂p i [ f (1 ± f ) D i ] (A6)+ 12 ∂∂p i (cid:20) ∂∂p j ( D ij f ) ± f ∂∂p j D ij (cid:21) , where the diffusion coefficients are defined as D i ( p ) = (cid:90) d q (2 π ) q i S ( p , p + q ) , (A7) D ij ( p ) = (cid:90) d q (2 π ) q i q j S ( p , p + q ) . (A8)are the diffusion coefficients.The gravitational scattering of a probe particle of mass M and a particle in the medium with mass m is describedby the spin-averaged matrix element |M| = 12 s + 1 (16 πG ) m M [( q ) − q ] , (A9)entering Eq. (A4). In the nonrelativistic limit, we canneglect q and maintain only q in Eq. (A9).6The problem of calculating the diffusion coefficientsfor different types of medium amounts to evaluatingEqs. (A7) and (A8), where in the response functionEq. (A4) we can select the appropriate sign in 1 ± f corresponding to the medium particle’s spin-statistics (orsetting 1 ± f → D || , the first diffusion coefficient corresponding to mo-tion parallel to the probe object’s instantaneous velocity. D || is simply given by Eq. (A7) when we select q i to alignwith the direction of p .
1. A classical gas medium
We first re-derive the relaxation of massive classicalobjects, such as supermassive black holes or GCs, in abackground medium consisting of other classical objectssuch as stars or CDM particles. In the nonrelativisticlimit, the function S ( p , p (cid:48) ) is simplified as S ( p , p (cid:48) ) (cid:39) g χ (4 πGmM ) q (cid:90) d k (2 π ) d k (cid:48) (2 π ) × (2 π ) δ (4) ( p + k − p (cid:48) − k (cid:48) ) f ( k ) (A10)where g χ is the number of internal degrees of freedom ofdark matter. Here, M and m are the mass of the particlespecies 1 and 2, respectively. In the small momentumexchange limit, the δ -function for the energy conservationcan be expanded as δ ( E p + E k − E p (cid:48) − E k (cid:48) ) (cid:39) q (cid:18) M µ r q · ∂∂ p (cid:19) δ (cid:104) ˆ q · (cid:18) k m − p M (cid:19) (cid:105) , (A11)where µ r = mM/ ( m + M ) is the reduced mass. Us-ing these approximate expressions in nonrelativistic andsmall momentum exchange limit, we obtain the diffusioncoefficients as D i ( p ) = (cid:90) d q (2 π ) q i S ( p , p + q ) (A12)= 4 πG m M (cid:18) Mm (cid:19) ln Λ ∂∂p i h ( p ; f )and D ij ( p ) = (cid:90) d q (2 π ) q i q j S ( p , p + q )= 4 πG m M ln Λ ∂ ∂p i ∂p j g ( p ; f ) (A13)where ln Λ = (cid:82) q max q min dq/q is the Coulomb logarithm, andwe have used identities Eqs. (26)–(27) of [34] to per-form the angular integration at the second step in eachequation. The Rosenbluth potentials h ( p ) and g ( p ) are defined as [86] h ( p ; f ) = g χ (cid:90) d k (2 π ) f ( k ) (cid:12)(cid:12) k m − p M (cid:12)(cid:12) , (A14) g ( p ; f ) = g χ (cid:90) d k (2 π ) (cid:12)(cid:12)(cid:12) k m − p M (cid:12)(cid:12)(cid:12) f ( k ) (A15)This reproduces the well known diffusion coefficientsin a classical system, see Eq. (7.83) in Binney &Tremaine [87]. For the Maxwell-Boltzmann distribution f ( k ) = (2 π ) / n / [ g χ ( mσ ) ] e − v k / σ , it is straightfor-ward to find ∂h∂p i = − v i v n M σ X (cid:20) erf( X ) − X √ π e − X (cid:21) ≡ − v i v n M σ G ( X ) , (A16) ∂g∂p i ∂p j = √ σ M (cid:104) X i X j X (cid:16) G ( X ) −
13 erf( X ) (cid:17) + δ ij X ) − G ( X ) X (cid:105) , (A17)where v = p /M , v = | v | and X = v √ σ . (A18)
2. Degenerate fermionic dark matter
We now consider the diffusion of astrophysical objectssuch as GCs in a halo of fermionic dark matter. In thiscase, the response function S becomes S ( p , p (cid:48) ) (cid:39) g χ (4 πGmM ) q (cid:90) d k (2 π ) d k (cid:48) (2 π ) × (2 π ) δ (4) ( p + k − p (cid:48) − k (cid:48) ) f ( k )(1 − f ( k (cid:48) ))(A19)Expanding f ( k (cid:48) ) around k , we find additional contribu-tion to the function S due to quantum statistics as∆ S ( p , p (cid:48) ) (cid:39) − πg χ (4 πGmM ) q × (A20) (cid:90) d k (2 π ) (cid:18) q · ∂∂ p (cid:19) δ (cid:104) ˆ q · (cid:18) k m − p M (cid:19)(cid:105) f ( k ) , which is the same as Eq. (A10) upon substituting M/µ → f → f . We find D i ( p ) = 4 πG m M ln Λ µ r ∂∂p i (cid:104) h ( p ; f ) − µ r M h ( p ; f ) (cid:105) (A21) D ij ( p ) = 4 πG m M ln Λ ∂ ∂p i ∂p j (cid:104) g ( p ; f ) − g ( p ; f ) (cid:105) (A22)For the degenerate case, one can perform the k and k (cid:48) integrations without expanding f ( k (cid:48) ). This computation7was already done in the context of neutrino transport ina hot and dense medium [88] and dark matter thermal-ization in neutron stars [89]. We find: S ( p , p (cid:48) ) = g χ (4 πGmM ) q m T πq z − e − z (cid:18) ξ − z (cid:19) (A23)where z = − q /T , E − = m + k − , k − = ( m /q )( q + q / m ) , and ξ − = ln (cid:20) e ( E − − µ ) /T e ( E − − µ ) /T e z (cid:21) . (A24)Integrating this response function with respect to q ,one obtains the diffusion coefficients for a degeneratemedium.
3. Ultralight dark matter
It is discussed in [90] that the dynamical relaxationof stars in a ULDM halo proceeds as stars scatter offULDM quasi-particles whose size is of the order of thede Broglie wavelength, λ dB ∼ π/mv . This observationis confirmed by Bar-Or et al [34], where the dynamicalrelaxation time scale as well as diffusion coefficients arecomputed in more rigorous way by using Fokker-Planckequation and stochastic gravitational potential.The Boltzmann equation approach can also reproducethe dynamical relaxation time scale and diffusion coef-ficients. The gravitational scattering between ultralightdark matter and star can be described by the same ma-trix element, Eq. (A9), where M and m are the mass ofstar and ultralight dark matter, respectively. We treatthe star as a pointlike particle, and this can be justi-fied since the maximum momentum exchange q ∼ mv ismuch smaller than 1 /r with a typical radius of star r .The function S ( p , p (cid:48) ) is S ( p , p (cid:48) ) (cid:39) g χ (4 πGmM ) q (cid:90) d k (2 π ) d k (cid:48) (2 π ) × (2 π ) δ (4) ( p + k − p (cid:48) − k (cid:48) ) f ( k )(1 + f ( k (cid:48) )) . (A25)The quantum correction is the same as Eq. (A21) withan opposite sign. Therefore, the diffusion coefficients are D i ( p ) = 4 πG m M ln Λ µ r ∂∂p i (cid:104) h ( p ; f ) + µ r M h ( p ; f ) (cid:105) (A26) D ij ( p ) = 4 πG m M ln Λ ∂ ∂p i ∂p j (cid:104) g ( p ; f ) + g ( p ; f ) (cid:105) (A27)This reproduces the results of [34]. Appendix B: Maximum entropy DDM halos
In the derivation of the quasi-degenerate density pro-file, we adopt the assumption that a galactic structuremay be described as a statistical ensemble close to equi-librium, in the sense of a maximal Boltzmann-Gibbs en-tropy. A similar approach to ours can be found in anumber of earlier works [40, 57, 91–94].The phase space distribution function f ( r , p ) and dif-ferential particle number density dN are related via(2 π ) dNd xd p = gf ( r , p ) . (B1)The entropy of the gas is then given by the functional S = − g (cid:90) d pd r (2 π ) [ f ln f + (1 − f ) ln(1 − f )] . (B2)Supplemented with Lagrange multipliers for the totalenergy and the total number of particles, the variationproblem can be carried out along the lines of [57]. Themaximum entropy result is f ( r , p ) = 11 + exp[ z ( r , p )] , (B3)where z ( r , p ) = β p m + βm Φ( r ) + α , (B4)with β and α the energy and particle number Lagrangemultipliers. The gravitational potential is given byΦ( r ) = − Gmg (cid:90) d p (cid:48) d r (cid:48) (2 π ) f ( r (cid:48) , p (cid:48) ) | r − r (cid:48) | . (B5)By construction, the gravitational potential solves thePoisson equation, ∇ Φ = 4 πGρ , where the density is ρ ( r ) = mg (cid:90) d p (2 π ) f ( r , p ) . (B6)We can make progress by evaluating the density, ρ ( r ) = − mg (cid:18) m πβ (cid:19) / PolyLog (cid:20) , − e ϕ (cid:21) , (B7)where ϕ ≡ − βm Φ( r ) − α ≡ βµ = µ/T , defining also the“chemical potential” µ ( r ) = − m Φ( r ) − (cid:101) α , with (cid:101) α ≡ α/β .Using the Poisson equation and the definition of ϕ , wehave ∇ Φ = − βm ∇ ϕ = 4 πGρ . (B8)In the degenerate limit, µ/T = ϕ (cid:29)
1, The PolyLogfunction asymptotes to − PolyLog (cid:20) , − e ϕ (cid:21) → ϕ / √ π = 4( µ/T ) / √ π . (B9)8It’s therefore useful to rewrite Eq. (B8) as ∇ µ ( r ) (B10)= 4 √ π gGm / µ / PolyLog (cid:2) , − e βµ f (cid:3) √ π ( βµ ) / , where we defined µ (0) ≡ µ and µ ( r ) = µ h ( r ). Let usalso define a dimensionless radius x via r = r x , with r given by r = (cid:115) π √ gGm / µ / . (B11)Then, one finds the equation ∂ x ( x ∂ x h ) = x PolyLog (cid:2) , − e ( µ /T ) h (cid:3) √ π ( µ /T ) / . (B12)Given the solution for h , the density is simply ρ = √ π gm µ − PolyLog (cid:104) , − e µ T h (cid:105) √ π ( µ /T ) / . (B13)In the limit µ /T (cid:29)
1, the right-hand side ofEq. (B12) becomes − x h / , which is a Lane-Emdenequation with a scale r . The central density be-comes ρ (0) ≡ ρ ≈ ( √ / π ) gm / µ / , i.e. µ =(3 π / √ / ( ρ /g ) / /m / , which implies r ≈ (cid:18) π G ρ g m (cid:19) . (B14)The solution for ρ is constant near the origin and falls as ρ ∝ /r at r (cid:29) r , with a “wriggle” feature near r ∼ r .An example with µ /T = 10 is shown in Fig. 13. Our natural unit notation, which includes (cid:126) →
1, may mask thefact that the characteristic radius Eq. (B14) is determined byquantum degeneracy pressure. This is easy to unmask by restor-ing (2 π ) → (2 π (cid:126) ) on the left-hand side of Eq. (B1). Tracking (cid:126) through the computation gives a factor of (cid:126) on the right-handside of Eq. (B14). Exact solution, μ ( ) T + μ ( ) T x - - - (cid:0) ρ [ a r b . un i t s ] FIG. 13. The density profile found by solving Eq. (B12) andinserting into Eq. (B13) (blue line), compared with an analyt-ical ansatz (orange) that demonstrates asymptotic behaviorof the solution. Evidently, for x (cid:46) /x for large x . It is interesting to compare this class of solutions to thesolutions obtained from the prescription of [41]. Rescal-ing the equations of [41] by r from Eq. (B14), we plotthe solution we find (at constant ρ ) in Fig. 14 (namedRSU), along with different solutions corresponding dif-ferent µ /T . Evidently, the profile used by [41] bears astrong resemblance to the µ /T ∼ µ /T → ∞ , the solution is a corewith finite radius. We show it both by solving the den-sity profile for µ /T (cid:29) μ ( )/ T1 310 20
RSULE - - - - - - x ρ [ k g / m ] FIG. 14. The density profile found by solving Eq. (B12) andinserting into Eq. (B13) for different values of µ /T , keeping ρ constant. The dashed blue line (RSU) is based on thetreatment of [41]. The dotted red line (LE) is based on solvingthe Lane-Emden equation, which is the µ /T → ∞ limit ofthe equations, as described in the text. μ ( )/ T1 310 20
RSULE
500 1000 1500 20000510152025 r [ pc ] V e l o c i t i e s [ k m / s ] FIG. 15. The circular velocities (cid:112) GM ( r ) /r for the densityprofiles that appear in Fig. 14. Appendix C: Jeans modeling
Following [45], the equation for the second velocity mo-ments of a static spherical system of particles under theinfluence of a gravitational potential with enclosed mass M ( r ) is 1 ν ddr ( ν ¯ v r ) + 2 β ¯ v r r = − GMr , (C1)where ν ( r ) is the particles’ density, ¯ v r ( r ) is the radialsecond velocity moment, ¯ v θ is the angular second velocitymoment and β ≡ − ¯ v θ / ¯ v r is the velocity anisotropy. Forconstant β , Eq. (C1) is solved by ν ¯ v r ( r ) = Gr β ∞ (cid:90) r r (cid:48) β − ν ( r (cid:48) ) M ( r (cid:48) ) dr (cid:48) . (C2)This can be related to the line-of-sight (LOS) velocity, σ ( r ) = 2 I ( r ) ∞ (cid:90) r (cid:18) − β r r (cid:48) (cid:19) ν ¯ v r ( r (cid:48) ) r (cid:48) √ r (cid:48) − r dr (cid:48) . (C3)In modelling stellar kinematics in Fornax, we assume aPlummer profile with density and surface density givenby, respectively, ν ( r ) = 1(1 + r /r p ) / L πr p (C4) I ( r ) = 1(1 + r /r p ) Lπr p . (C5)We use the radius parameter r p = 710 pc [10], consis-tent with the stellar sample on which the kinematics datais based. (This radius parameter is about 20% smallerthan the r p ≈
851 pc reported in a new morphologicalstudy [20]. The difference is not crucial for our analy-sis. Moreover, we prefer to consider the photometry andspectroscopy of the same data set.)
Appendix D: CDM velocity dispersion in a coredprofile
In this appendix we discuss some features of coredCDM halos, notably DF, partially following [42].The velocity dispersion of dark matter is important tothe discussion, therefore let us write the Jeans equationfor the second radial velocity moment with a constantvelocity anisotropy β , which has the solution [45]¯ v r ( r ) = Gr β ρ ( r ) ∞ (cid:90) r r (cid:48) β ρ ( r (cid:48) ) M ( r (cid:48) ) r (cid:48) dr (cid:48) , (D1)which is the same as Eq. (C2) but in slightly differentnotation.Consider a finite-core toy model, where the density is ρ (0) for r < r c and 0 for r > r c . Then, the solution ofEq. (D1) is [42]¯ v r ( r ) = 2 πGρ (0)3( β + 1) 1 r β (cid:0) r β +2c − r β +2 (cid:1) . (D2)For isotropic velocity dispersion β = 0, this reduces to¯ v r ( r ; β = 0) = 2 πGρ (0)3 (cid:0) r − r (cid:1) . (D3)Thus, for r (cid:28) r c , σ r ≡ (cid:113) ¯ v r ( r (cid:28) r c ; β = 0) ≈ (cid:114) πGρ (0)3 r c (D4) ≈ (cid:32) ρ (0)10 M (cid:12) kpc (cid:33) r c . We can also note the ratio, X ≡ V circ √ σ r ≈ (cid:114) GM ( r )2 r (cid:46)(cid:114) πGρ (0) r c rr c , (D5)indicating that the low-velocity approximation of theChandrasekhar deceleration may apply inside a core, seeEq. (5). This implies a “phase-space suppression” to DF,as discussed in the main text (Sec. III). Appendix E: Orbits under dynamical friction
In this appendix we review the solution of an orbitunder the influence of DF. We write the equations ofmotion (EoM) in circular coordinates,¨ r = (¨ r − r ˙ ϕ )ˆ r + (2 ˙ r ˙ ϕ + r ¨ ϕ ) ˆ ϕ (E1)= − GM ( r ) r ˆ r − (cid:12)(cid:12)(cid:12) d ˙ r dt (cid:12)(cid:12)(cid:12) DF ˙ r | ˙ r | . (E2)We express the deceleration | d ˙ r /dt | DF as | ˙ r | /τ , where τ appears in Eq. (10).0Defining r = R x , t = T ¯ t , T = R /GM ( R ), we find x (cid:48)(cid:48) − xϕ (cid:48) = − x M ( R x ) M ( R ) − x (cid:48) τ /T (E3)2 x (cid:48) ϕ (cid:48) + xϕ (cid:48)(cid:48) = − xϕ (cid:48) τ /T (E4)where (cid:48) is differentiation with respect to ¯ t . Note, τ candepend on r and | ˙ r | . For a circular orbit, for example,the initial conditions can be set as x (0) = 1, x (cid:48) (0) = 0, ϕ (0) = 0 and ϕ (cid:48) (0) = 1 /x (0) = 1, which has a revolutiontime of ∆ (cid:101) t = 2 π .Solving the orbit of a decelerating test object generallyrequires numerical integration. We can understand somefeatures of the solution analytically, however. Defining v ϕ ≡ r ˙ ϕ , the ˆ ϕ part of the EoM has the solution rv ϕ = ( rv ϕ ) exp − t (cid:90) dt (cid:48) τ . (E5)This solution expresses the decay of angular momentumof the test object. Using the circular velocity v = GM ( r ) /r , we can express the ˆ r part of the EoM as v ϕ − v = r (cid:18) ¨ r + ˙ rτ (cid:19) . (E6)We can gain more analytical intuition by consideringnearly-circular orbits, assuming that the inspiral rate ismuch smaller than the circular velocity, r/τ (cid:28) v circ . As-suming that ˙ r ∼ r/τ , ¨ r ∼ r/τ and r/τ (cid:28) v circ , Eq. (E6)implies v ϕ ≈ v circ . We can use this to write − rv ϕ τ = ˙ rv ϕ + r ˙ v ϕ ≈ ˙ rv circ + r ˙ v circ (E7)= 12 v circ ˙ r (cid:18) d ln Md ln r (cid:19) . (E8)Rearranging, we find˙ rr ≈ − (cid:0) d ln Md ln r (cid:1) τ . (E9)Using this, we can estimate the time it takes a test objectto fall from r down to r < r : t ( r ; r ) = r (cid:90) r dr r (cid:18) d ln Md ln r (cid:19) τ ( r, v circ ( r )) . (E10)Given the mass profile of the halo, M ( r ), and a DF modelencapsulated by τ , Eq. (E10) is a simple and quick esti-mate of the inspiral time of a test object.For eccentric orbits, the approximation above isless justified. Defining eccentricity as e ≡ ( r apo − r peri ) / ( r apo + r peri ) with apocenter radius r apo and peri-center radius r peri , we numerically tested Eq. (E10) for e >
0. In these calculations we defined r and r via( r apo + r peri ) /
2, where r apo and r peri are obtained per cy-cle of the orbital phase. With these definition, in numer-ical experiments representative of Fornax GCs we findthat Eq. (E10) holds to better than 30% accuracy for e (cid:46) . Appendix F: The radial and projected CDF of GCs
Consider a population of identical GCs (all with thesame mass), that start off their life at some initial time t = 0 on approximately circular orbits with a radial prob-ability distribution function (PDF) f ( r ) w.r.t. to initialradial coordinate r . The CDF of initial GC positionsis F ( r ) = (cid:82) r dyf ( y ). We are interested in comput-ing the PDF and CDF of GC radial positions today, at t = ∆ t ; call these f ∆ t ( r ) and F ∆ t ( r ).DF causes GC orbits to inspiral inwards, and inte-grating along the orbit we can compute the function r = r ( r ; ∆ t ) and invert it to obtain r = r ( r ; ∆ t ). Neglecting tidal disruption, we have F ∆ t ( r ) = F ( r ( r ; ∆ t )) . (F1)Now we can use explicit results for r ( r ; ∆ t ) to connect F ∆ t ( r ) with F ( r ) in different halo models. To this endwe can use Eq. (E10),∆ t = r (cid:90) r dr (cid:48) r (cid:48) (1 + α ( r (cid:48) )) τ ( r (cid:48) ) , (F2)where α ( r ) ≡ d ln M/d ln r . Let us consider the generalfeatures of F ∆ t for different halo shapes.
1. CDF of GCs in a cuspy halo
We have seen in the main analysis that a cuspy halo(i.e. the inner region of an NFW halo, where α ≈ τ . For an approximately constant α and power law τ = ¯ τ ( r/ ¯ r ) β , it is useful to define the critical radius r cr via τ ( r cr ) = 2 β α ∆ t. (F3)The physical meaning of r cr is that GCs that start theirlife at r ≤ r cr arrive at the origin within t ≤ ∆ t . Usingour power law form for τ , we have r cr = ¯ r (cid:18) β α ∆ t ¯ τ (cid:19) /β . (F4)In terms of r cr , the solution of Eq. (F2) evaluates to r ( r ; ∆ t ) = r cr (cid:32) (cid:18) rr cr (cid:19) β (cid:33) /β . (F5) The monotonous decrease of r with time, that allowed this in-version, is lost for non-circular orbits. We could accommodateelliptical orbits approximately, by letting r represent the averagebetween peri- and apo-center per cycle. r (cid:28) r cr today, we can expandtheir starting point: r ( r ; ∆ t ) = r cr + r cr β (cid:18) rr cr (cid:19) β + ... (F6)= r cr (cid:18) α β τ ( r )∆ t + ... (cid:19) . In other words, for cuspy CDM halos, GCs that are cur-rently seen at r (cid:28) r cr must have originated near r ≈ r cr .This means that for GCs with r (cid:28) r cr today, the ra-dial distribution today is not very sensitive to the (dif-ficult to predict) initial distribution. We can make thispoint manifest by expanding Eq. (F1), using Eqs. (F4)and (F6): F ∆ t ( r ) ≈ F ( r cr ) + (1 + α )2 β f ( r cr ) r cr τ ( r )∆ t + ... , (F7)where the ... refer to higher powers of the small ra-tio τ ( r ) / ∆ t . Above, the r -independent constant F ( r cr )counts GCs that have already settled to the center ofthe halo. These GCs at r ≈ F ( r cr ) on theRHS of Eq. (F7). We thus have F ∆ t ( r ) ≈ A ( τ ( r ) / ∆ t ),where A = (1+ α )2 β f ( r cr ) r cr is an order unity constant ( r -independent) coefficient. We can summarise this section with two importantconclusions. First, for a cuspy halo, all GCs that areborn at r < r cr have arrived at r ≈ f ( r ), that differ from each otheronly at r < r cr ; unless stellar age and metallicity mea-surements can identify the remnants and approximatelycount tidally disrupted GCs, on time scales of Gyrs afterthe disruption. Second, the radial CDF of GCs at smallradii r (cid:28) r cr should follow F ∆ t ( r ) ≈ A ( τ ( r ) / ∆ t ), withorder unity A , irrespective of initial conditions.Fig. 16 illustrates both of these two points, by showingtwo examples of f and the resulting F ∆ t . The DF time τ ( r ) and the critical radius r cr are measurable given amodel of the DM halo, fitted to stellar kinematics, andgiven GC age measurements that define ∆ t . This makesthe predicted shape of F ∆ t measurable, in principle. Inpractice, however, projection effects (explained below)complicate the interpretation. In addition, the collectionof GCs in Fornax seems too sparse to draw robust con-clusions. For the inner region of an NFW profile, we have seen that α ≈ β ≈
2, so A ≈ . N cr , where N cr = f ( r cr ) r cr counts thenumber of GCs that were located in a region of order r cr around r cr . Predicting the actual value of N cr would require understand-ing of the initial cosmological formation of GCs, which is still notunder full theory control.
2. CDF of GCs in a cored halo
Inside a core we expect α ≈ τ . Eq. (F2) is evaluated to∆ t = 1 + α τ ln r r (F8)and the radial CDF today is F ∆ t ( r ) ≈ F ( re t (1+ α ) τ ) ≈ F ( re ∆ t τ ) . (F9)The distribution of GCs inside a core reflects a stretchedversion of the initial conditions. Because of this sensitiv-ity to initial conditions, the degree of possible fine-tuningin the current positions of GCs may be difficult to assess.
3. Accounting for distribution of GC masses
The masses of GCs in Fornax vary over about an orderof magnitude around 10 M (cid:12) , and the instantaneous DFtime satisfies τ ∝ /m ∗ up to logarithmic correctionsthat we neglect here. It is therefore necessary to revisethe prediction of the radial CDF of GCs to account fordifferent GC masses.In the case of a cuspy profile, where τ ∝ r β , the criticalradius scales as r cr ∝ m /β ∗ . For example, using the NFWfit of Fornax (which gives β ≈ . t = 12 Gyr we have: r cr ≈ . (cid:18) m (cid:63) M GC4 (cid:19) . kpc (F10) ≈ . (cid:18) m (cid:63) M GC3 (cid:19) . kpc . (F11)Suppose we have a set of GC masses m ∗ i with initialradial distribution functions f ,i ( r ). Summing over allGC masses we find that as long as r (cid:28) r cr ,i , Eq. (F7)predicts that the total radial CDF today is (again omit-ting GCs that have already settled to the center of thehalo) (cid:88) i F ∆ t,i ( r ) ≈ F ∆ t, ( r ) (cid:88) i f ,i ( r cr ,i ) f , ( r cr , ) (cid:18) m ∗ i m ∗ (cid:19) β − , (F12)where F ∆ t, ( r ) is the radial CDF of GCs of mass m ∗ .We see that Eq. (F12) simply reproduces Eq. (F7) up toa modified over-all multiplicative constant.
4. Projected radius distribution
In reality we only know the projected distance of GCsfrom the center of Fornax, r ⊥ , and not the true radial2 FIG. 16. Initial radial PDF f ( left ) and resulting current radial CDF F ∆ t ( right ) for a cuspy halo with τ ( r ) ∝ r . and r cr = 1 .
12 (units on x-axis are arbitrary). Two different examples for f are shown, leading to nearly identical F ∆ t . The twoversions of f are normalised to yield F ∆ t ( ∞ ) = 6. On both panels, r cr is marked with a vertical black line. On the left, the r (cid:28) r cr approximation F ∆ t ≈ A ( τ ( r ) / ∆ t ) is shown by the red dotted curve. distance r . To obtain the CDF of projected radii, we canstart with the surface density of GCs,Σ ∆ t ( r ⊥ ) = ∞ (cid:90) −∞ dzn ( r ) = 2 ∞ (cid:90) r ⊥ dr rn ( r ) (cid:112) r − r ⊥ , (F13)where the 3D number density n ( r ) is related to the radialPDF via n = f ∆ t ( r ) / πr . Using this relation we haveΣ ∆ t ( r ⊥ ) = 12 π ∞ (cid:90) r ⊥ dr f ∆ t ( r ) r (cid:112) r − r ⊥ . (F14)The CDF in r ⊥ , that we define by F ⊥ ∆ t ( r ⊥ ), is given by: F ⊥ ∆ t ( r ⊥ ) = 2 π r ⊥ (cid:90) dRR Σ ∆ t ( R ) (F15)= F ∆ t ( r ⊥ ) + ∞ (cid:90) r ⊥ drf ∆ t ( r ) (cid:32) − (cid:114) − r ⊥ r (cid:33) . The CDF of projected radii contains the CDF of trueradii, evaluated at r = r ⊥ , plus another term that countsGCs at r > r ⊥ which projection casts into LOS insideof r ⊥ . The added projection term can exceed the un-projected term, meaning that most GCs seen inside r
The goal of this section is to explore the implicationsof uncertainties due to the line-of-sight projection in thetrue positions and velocities of GCs. Different projection
FIG. 17. The effect of projection. Solid blue: un-projected ra-dial CDF computed at the projected radius, F ∆ t ( r ⊥ ). Dashedblack: CDF of projected radii F ⊥ ∆ t ( r ⊥ ). angles lead to different true positions and velocities of aGC, affecting the orbital settling time under DF.Consider the orbits of test bodies in a spherically-symmetric gravitational potential Φ( r ). A given orbitlies on an orbital plane. In the coordinate system inFig. 18, one can parameterize the orbital plane with theunit vector ˆ n = (cos α sin β, sin α sin β, cos β ). On the or-bital plane, the radius r ( ϕ ) and phase ϕ ( t ) completely de-fine the orbit. As in App. E, defining r = r true x , t = T ¯ t , T = r /GM ( r true ), one finds x (cid:48)(cid:48) − xϕ (cid:48) ≈ − x M ( r true x ) M ( r true ) (G1) (cid:0) x ϕ (cid:48) (cid:1) (cid:48) ≈ (cid:48) is differentiation with respect to ¯ t . (Herewe tentatively neglect DF.) The initial conditions are3 x (0) = 1, ϕ (0) = 0, x (cid:48) (0) = ( − ∆ v r cos θ +∆ v y sin θ ) / ( r true /T ), θ (0) = 0 and | ϕ (cid:48) (0) | = (cid:113) ∆ v z + (∆ v r sin θ + ∆ v y cos θ ) / ( x (0) r true /T ). Ev-idently, x (cid:48) (0) + ( x (0) ϕ (cid:48) (0)) = (∆ v r + ∆ v y +∆ v z ) / ( r true /T ) .Given measured r ⊥ , ∆ v r and a model Φ( r ), we ex-plore the remaining orbital parameters which affect theinspiral time. We start with the true radius r true . Given r true and, for simplicity, assuming a circular orbit, theprobability of observing r ⊥ < x r true is P ( r ⊥ /r true 2) = 1 / 3. We therefore explore orbitswith true radius in the range r true = r ⊥ × [1 , v true ∈ [max(∆ v r , . V circ ( r )) , max(∆ v r , . V circ ( r ))]. Definingthe eccentricity e ≡ ( r apo − r peri ) / ( r apo + r peri ), we find amaximal e ∼ . − . v true in the central (cid:46) v true we expect smaller inspiral times. Increas-ing the upper bound, however, results in larger inspiraltimes – but also in more tuning. An test object spendsrelatively little time near the pericenter. Specifically, T peri /T apo ≡ ( r peri /v peri ) / ( r apo /v apo ) ≈ (1 − e ) / (1 + e ) ,yielding about 1 / e = 0 . v z > v y is,however, important: under the transformation ∆ v y →− ∆ v y , cos θ → − cos θ , so | ϕ (cid:48) (0) | remains constant, but x (cid:48) (0) → − x (cid:48) (0). Since the specific energy is (cid:15) ≈ v r ) = ˙ r l r + Φ( r ) , (G3)where l is the specific angular momentum, this trans-formation returns the same orbit. The inspiral time istherefore invariant under this transformation. We shallexplore then ∆ v y > θ positive or negative.To sum up, for each GC (and a given model of thehalo), we scan the range r true ∈ [1 , r proj . For each r true we scan over V true ∈ [0 . , . V circ ( r true ). For eachtrue velocity we scan positive and negative cos θ . Finally,we test the two cases, ∆ v y = (cid:112) v − ∆ v r , ∆ v z = 0and ∆ v y = 0 , ∆ v z = (cid:112) v − ∆ v r . For each point inphase-space, we integrate the full equations of motion asin App. E. For each integration, we stop when ( r apo + r peri ) / (cid:46) . r initial or after 10 Gyr (the first of the two).We then denote the integration time as τ inspiral . FIG. 18. The coordinate system that we adopt to analysea given GC. The galactic dynamical center is in the origin.The observer is located at a very large X . The GC is locatedsomewhere on the line Z = 0 , Y = r ⊥ . The true radius istherefore r true = r ⊥ / sin α . We assume ∆ v r , the componentof velocity in the X direction, can be measured. We assumethat the rest of the components cannot be measured for now.The dotted line is the quasi-stable orbit of the GC. [1] M. G. Walker, M. Mateo, E. W. Olszewski,J. Penarrubia, N. W. Evans, and G. Gilmore, “Auniversal mass profile for dwarf spheroidal galaxies?,”The Astrophysical Journal no. 2, (2009) 1274.[2] D. R. Cole, W. Dehnen, J. I. Read, and M. I.Wilkinson, “The mass distribution of the Fornax dSph:constraints from its globular cluster distribution,” Mon.Not. Roy. Astron. Soc. (2012) 601, arXiv:1205.6327 [astro-ph.CO] .[3] F. Hammer, Y. Yang, F. Arenou, J. Wang, H. Li,P. Bonifacio, and C. Babusiaux, “Orbital evidences fordark-matter-free Milky Way dwarf spheroidal galaxies,” arXiv:2002.09493 [astro-ph.GA] .[4] R. A. Flores and J. R. Primack, “Observational andtheoretical constraints on singular dark matter halos,”Astrophys. J. Lett. (1994) L1–4, arXiv:astro-ph/9402004 .[5] B. Moore, “Evidence against dissipation-less darkmatter from observations of galaxy haloes,” Nature no. 6491, (Aug., 1994) 629–631.[6] W. De Blok, “The core-cusp problem,” Advances inAstronomy (2010) .[7] A. Fattahi, J. F. Navarro, T. Sawala, C. S. Frenk, L. V.Sales, K. Oman, M. Schaller, and J. Wang, “The colddark matter content of Galactic dwarf spheroidals: nocores, no failures, no problem,” arXiv:1607.06479[astro-ph.GA] .[8] J. I. Read, G. Iorio, O. Agertz, and F. Fraternali, “Thestellar mass–halo mass relation of isolated field dwarfs:a critical test of ΛCDM at the edge of galaxyformation,” Mon. Not. Roy. Astron. Soc. no. 2,(2017) 2019–2038, arXiv:1607.03127 [astro-ph.GA] .[9] L. J. Chang and L. Necib, “Dark Matter DensityProfiles in Dwarf Galaxies: Linking Jeans ModelingSystematics and Observation,” arXiv:2009.00613[astro-ph.CO] .[10] J. Read, M. Walker, and P. Steger, “Dark matter heatsup in dwarf galaxies,” Mon. Not. Roy. Astron. Soc. no. 1, (2019) 1401–1420, arXiv:1808.06634[astro-ph.GA] .[11] J. S. Bullock and M. Boylan-Kolchin, “Small-ScaleChallenges to the ΛCDM Paradigm,” Ann. Rev.Astron. Astrophys. (2017) 343–387, arXiv:1707.04256 [astro-ph.CO] .[12] J. Simon, K. Bechtol, S. Birrer, S. Chakrabarti, F.-Y.Cyr-Racine, I. Dell’Antonio, A. Drlica-Wagner,C. Fassnacht, M. Geha, D. Gilman, Y. Hezaveh,D. Kim, T. Li, L. Strigari, and T. Treu, “Testing thenature of dark matter with extremely large telescopes.”Snowmass 2021 - letter of interest, 2020. .[13] S. Tremaine, “The formation of the nuclei of galaxies.ii-the local group,” The Astrophysical Journal (1976) 345–351.[14] M.-Y. Wang, S. Koposov, A. Drlica-Wagner, A. Pieres,T. Li, T. de Boer, K. Bechtol, V. Belokurov, A. Pace,D. Bacon, et al., “Rediscovery of the sixth star clusterin the fornax dwarf spheroidal galaxy,” TheAstrophysical Journal Letters no. 2, (2019) L13. [15] S. Chandrasekhar, “Dynamical Friction. I. GeneralConsiderations: the Coefficient of Dynamical Friction.,”Astrophys. J. (Mar, 1943) 255.[16] N. Meadows, J. F. Navarro, I. Santos-Santos,A. Ben´ıtez-Llambay, and C. Frenk, “Cusp or core?Revisiting the globular cluster timing problem inFornax,” Mon. Not. Roy. Astron. Soc. no. 3, (Jan,2020) 3336–3342, arXiv:1910.11887 [astro-ph.GA] .[17] T. de Boer and M. Fraser, “Four and one more: Theformation history and total mass of globular clusters inthe fornax dsph,” Astronomy & Astrophysics (2016) A35.[18] A. D. Mackey and G. F. Gilmore, “Surface brightnessprofiles and structural parameters for globular clustersin the Fornax and Sagittarius dwarf spheroidalgalaxies,” Mon. Not. Roy. Astron. Soc. (2003) 175, arXiv:astro-ph/0211396 [astro-ph] .[19] A. del Pino, S. L. Hidalgo, A. Aparicio, C. Gallart,R. Carrera, M. Monelli, R. Buonanno, and G. Marconi,“Spatial dependence of the star formation history in thecentral regions of the fornax dwarf spheroidal galaxy,”Monthly Notices of the Royal Astronomical Society no. 2, (2013) 1505–1516.[20] M.-Y. Wang, T. de Boer, A. Pieres, T. Li,A. Drlica-Wagner, S. Koposov, A. Vivas, A. Pace,B. Santiago, A. Walker, et al., “The morphology andstructure of stellar populations in the fornax dwarfspheroidal galaxy from dark energy survey data,” TheAstrophysical Journal no. 2, (2019) 118.[21] J. F. Navarro, C. S. Frenk, and S. D. M. White, “AUniversal density profile from hierarchical clustering,”Astrophys. J. (1997) 493–508, arXiv:astro-ph/9611107 [astro-ph] .[22] X. Hernandez and G. Gilmore, “Dynamical friction indwarf galaxies,” Mon. Not. Roy. Astron. Soc. (1998) 517, arXiv:astro-ph/9802261 .[23] K. S. Oh, D. N. C. Lin, and H. B. Richer, “GlobularClusters in the Fornax Dwarf Spheroidal Galaxy,”Astrophys. J. no. 2, (Mar., 2000) 727–738.[24] J. M. Lotz, R. Telford, H. C. Ferguson, B. W. Miller,M. Stiavelli, and J. Mack, “Dynamical friction in deglobular cluster systems,” Astrophys. J. (2001) 572, arXiv:astro-ph/0102079 .[25] T. Goerdt, B. Moore, J. I. Read, J. Stadel, andM. Zemp, “Does the Fornax dwarf spheroidal have acentral cusp or core?,” Mon. Not. Roy. Astron. Soc. no. 3, (May, 2006) 1073–1077, arXiv:astro-ph/0601404[astro-ph] .[26] G. W. Angus and A. Diaferio, “Resolving the timingproblem of the globular clusters orbiting the fornaxdwarf galaxy,” Monthly Notices of the RoyalAstronomical Society no. 2, (2009) 887–893.[27] R. Cowsik, K. Wagoner, E. Berti, and A. Sircar,“Internal dynamics and dynamical friction effects in thedwarf spheroidal galaxy in Fornax,” Astrophys. J. (2009) 1389–1394, arXiv:0904.0451 [astro-ph.CO] .[28] K. Kaur and S. Sridhar, “Stalling of Globular ClusterOrbits in Dwarf Galaxies,” Astrophys. J. no. 2,(Dec., 2018) 134, arXiv:1810.00369 [astro-ph.GA] .[29] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten,“Ultralight scalars as cosmological dark matter,” Phys. Rev. D95 no. 4, (2017) 043541, arXiv:1610.08297[astro-ph.CO] .[30] G. Y. C. Leung, R. Leaman, G. van de Ven, andG. Battaglia, “A dwarf–dwarf merger and dark mattercore as a solution to the globular cluster problems inthe fornax dsph,” Monthly Notices of the RoyalAstronomical Society no. 1, (Dec, 2019) 320–336. http://dx.doi.org/10.1093/mnras/stz3017 .[31] P. Boldrini, R. Mohayaee, and J. Silk, “Embeddingglobular clusters in dark matter minihaloes solves thecusp–core and timing problems in the fornax dwarfgalaxy,” Monthly Notices of the Royal AstronomicalSociety no. 3, (2020) 3169–3178.[32] L. Berezhiani, B. Elder, and J. Khoury, “DynamicalFriction in Superfluids,” JCAP no. 10, (2019) 074, arXiv:1905.09297 [hep-ph] .[33] S. Hartman, H. Winther, and D. Mota, “Dynamicalfriction in finite temperature superfluids, and theFornax dwarf spheroidal,” arXiv:2011.00116[astro-ph.CO] .[34] B. Bar-Or, J.-B. Fouvry, and S. Tremaine, “Relaxationin a Fuzzy Dark Matter Halo,” Astrophys. J. (2019) 28, arXiv:1809.07673 [astro-ph.GA] .[35] L. Lancaster, C. Giovanetti, P. Mocz, Y. Kahn,M. Lisanti, and D. N. Spergel, “Dynamical Friction in aFuzzy Dark Matter Universe,” JCAP no. 01,(2020) 001, arXiv:1909.06381 [astro-ph.CO] .[36] V. Irˇsiˇc, M. Viel, M. G. Haehnelt, J. S. Bolton, andG. D. Becker, “First constraints on fuzzy dark matterfrom Lyman- α forest data and hydrodynamicalsimulations,” Phys. Rev. Lett. no. 3, (2017)031302, arXiv:1703.04683 [astro-ph.CO] .[37] N. Bar, D. Blas, K. Blum, and S. Sibiryakov, “Galacticrotation curves versus ultralight dark matter:Implications of the soliton-host halo relation,” Phys.Rev. D98 no. 8, (2018) 083027, arXiv:1805.00122[astro-ph.CO] .[38] D. J. Marsh and J. C. Niemeyer, “Strong Constraintson Fuzzy Dark Matter from Ultrafaint Dwarf GalaxyEridanus II,” Phys. Rev. Lett. no. 5, (2019)051103, arXiv:1810.08543 [astro-ph.CO] .[39] M. Safarzadeh and D. N. Spergel, “Ultra-light DarkMatter is Incompatible with the Milky Way’s DwarfSatellites,” arXiv:1906.11848 [astro-ph.CO] .[40] V. Domcke and A. Urbano, “Dwarf spheroidal galaxiesas degenerate gas of free fermions,” JCAP no. 01,(2015) 002, arXiv:1409.3167 [hep-ph] .[41] L. Randall, J. Scholtz, and J. Unwin, “Cores in dwarfgalaxies from fermi repulsion,” Monthly Notices of theRoyal Astronomical Society (Jan, 2017) stx161. http://dx.doi.org/10.1093/mnras/stx161 .[42] J. A. Petts, A. Gualandris, and J. I. Read, “Asemi-analytic dynamical friction model that reproducescore stalling,” Mon. Not. Roy. Astron. Soc. no. 4,(Dec., 2015) 3778–3791, arXiv:1509.07871[astro-ph.GA] .[43] M. Kaplinghat, S. Tulin, and H.-B. Yu, “Dark MatterHalos as Particle Colliders: Unified Solution toSmall-Scale Structure Puzzles from Dwarfs to Clusters,”Phys. Rev. Lett. no. 4, (2016) 041302, arXiv:1508.03339 [astro-ph.CO] .[44] S. Shao, M. Cautun, C. S. Frenk, M. Reina-Campos,A. J. Deason, R. A. Crain, J. D. Kruijssen, andJ. Pfeffer, “The survival of globular clusters in a cuspy Fornax,” arXiv:2012.08058 [astro-ph.GA] .[45] J. Binney and S. Tremaine,Galactic Dynamics: Second Edition. 2008.[46] S. Tremaine and M. D. Weinberg, “Dynamical frictionin spherical systems.,” Mon. Not. Roy. Astron. Soc. (Aug., 1984) 729–757.[47] M. D. Weinberg, “Orbital Decay of Satellite Galaxies inSpherical Systems,” Astrophys. J. (Jan., 1986) 93.[48] B. Bar-Or, J.-B. Fouvry, and S. Tremaine, “Relaxationin a Fuzzy Dark Matter Halo. II. Self-consistent kineticequations,” arXiv:2010.10212 [astro-ph.GA] .[49] B. Hendricks, A. Koch, M. Walker, C. I. Johnson,J. Pe˜narrubia, and G. Gilmore, “Insights from theoutskirts: Chemical and dynamical properties in theouter parts of the fornax dwarf spheroidal galaxy,”Astronomy & Astrophysics (Dec, 2014) A82. http://dx.doi.org/10.1051/0004-6361/201424645 .[50] A. Lauberts, ESO/Uppsala survey of the ESO(B) atlas.1982.[51] B. Letarte, V. Hill, P. Jablonka, E. Tolstoy, P. Fran¸cois,and G. Meylan, “VLT/UVES spectroscopy of individualstars in three globular clusters in the Fornax dwarfspheroidal galaxy,” Astron. Astrophys. no. 2, (July,2006) 547–554, arXiv:astro-ph/0603315 [astro-ph] .[52] J. E. Morrison, S. R¨oser, B. McLean, B. Bucciarelli, andB. Lasker, “The Guide Star Catalog, Version 1.2: AnAstrometric Recalibration and Other Refinements,”Astrophys. J. no. 3, (Mar., 2001) 1752–1763.[53] M. F. Skrutskie, R. M. Cutri, R. Stiening, M. D.Weinberg, S. Schneider, J. M. Carpenter, C. Beichman,R. Capps, T. Chester, J. Elias, J. Huchra, J. Liebert,C. Lonsdale, D. G. Monet, S. Price, P. Seitzer,T. Jarrett, J. D. Kirkpatrick, J. E. Gizis, E. Howard,T. Evans, J. Fowler, L. Fullmer, R. Hurt, R. Light,E. L. Kopan, K. A. Marsh, H. L. McCallon, R. Tam,S. Van Dyk, and S. Wheelock, “The Two Micron AllSky Survey (2MASS),” Astrophys. J. no. 2, (Feb.,2006) 1163–1183.[54] S. S. Larsen, J. P. Brodie, and J. Strader, “Detailedabundance analysis from integrated high-dispersionspectroscopy: globular clusters in the fornax dwarfspheroidal,” Astronomy & Astrophysics (Oct,2012) A53. http://dx.doi.org/10.1051/0004-6361/201219895 .[55] S. M. Inoue, “The test for suppressed dynamical frictionin a constant density core of dwarf galaxies,” Mon. Not.Roy. Astron. Soc. (2009) 709, arXiv:0901.4861[astro-ph.CO] .[56] J. I. Read, T. Goerdt, B. Moore, A. P. Pontzen,J. Stadel, and G. Lake, “Dynamical friction in constantdensity cores: A failure of the Chandrasekhar formula,”Mon. Not. Roy. Astron. Soc. (2006) 1451–1460, arXiv:astro-ph/0606636 [astro-ph] .[57] D. Lynden-Bell and R. Wood, “The gravo-thermalcatastrophe in isothermal spheres and the onset ofred-giant structure for stellar systems,” Mon. Not. Roy.Astron. Soc. (Jan., 1968) 495.[58] S. Tremaine and J. E. Gunn, “Dynamical Role of LightNeutral Leptons in Cosmology,” Phys. Rev. Lett. (1979) 407–410. [,66(1979)].[59] C. Di Paolo, F. Nesti, and F. L. Villante, “Phase spacemass bound for fermionic dark matter from dwarfspheroidal galaxies,” Mon. Not. Roy. Astron. Soc. no. 4, (2018) 5385–5397, arXiv:1704.06644 [astro-ph.GA] .[60] D. Savchenko and A. Rudakovskyi, “New mass boundon fermionic dark matter from a combined analysis ofclassical dSphs,” Mon. Not. Roy. Astron. Soc. no. 4,(2019) 5711–5720, arXiv:1903.01862 [astro-ph.CO] .[61] A. Boyarsky, O. Ruchayskiy, and D. Iakubovskyi, “ALower bound on the mass of Dark Matter particles,”JCAP (2009) 005, arXiv:0808.3902 [hep-ph] .[62] J. Alvey, N. Sabti, V. Tiki, D. Blas, K. Bondarenko,A. Boyarsky, M. Escudero, M. Fairbairn, M. Orkney,and J. I. Read, “New Constraints on the Mass ofFermionic Dark Matter from Dwarf SpheroidalGalaxies,” arXiv:2010.03572 [hep-ph] .[63] J. Bond, G. Efstathiou, and J. Silk, “Massive Neutrinosand the Large Scale Structure of the Universe,” Phys.Rev. Lett. (1980) 1980–1984.[64] M. Viel, J. Lesgourgues, M. G. Haehnelt, S. Matarrese,and A. Riotto, “Constraining warm dark mattercandidates including sterile neutrinos and lightgravitinos with WMAP and the Lyman-alpha forest,”Phys. Rev. D (2005) 063534, arXiv:astro-ph/0501562 .[65] M. Viel, G. D. Becker, J. S. Bolton, and M. G.Haehnelt, “Warm dark matter as a solution to the smallscale crisis: New constraints from high redshift Lyman-?forest data,” Phys. Rev. D88 (2013) 043502, arXiv:1306.2314 [astro-ph.CO] .[66] J. Baur, N. Palanque-Delabrouille, C. Y`eche,C. Magneville, and M. Viel, “Lyman-alpha Forests coolWarm Dark Matter,” JCAP (2016) 012, arXiv:1512.01981 [astro-ph.CO] .[67] M. Shoji and E. Komatsu, “Massive Neutrinos inCosmology: Analytic Solutions and FluidApproximation,” Phys. Rev. D (2010) 123516, arXiv:1003.0942 [astro-ph.CO] . [Erratum:Phys.Rev.D 82, 089901 (2010)].[68] G. Choi, M. Suzuki, and T. Yanagida, “Degeneratefermion dark matter from a broken U (1) B − L gaugesymmetry,” Phys. Rev. D no. 3, (2020) 035022, arXiv:2004.07863 [hep-ph] .[69] D. N. Spergel and P. J. Steinhardt, “Observationalevidence for selfinteracting cold dark matter,” Phys.Rev. Lett. (2000) 3760–3763, arXiv:astro-ph/9909386 .[70] M. Kaplinghat, R. E. Keeley, T. Linden, and H.-B. Yu,“Tying Dark Matter to Baryons with Self-interactions,”Phys. Rev. Lett. (2014) 021302, arXiv:1311.6524[astro-ph.CO] .[71] M. Rocha, A. H. Peter, J. S. Bullock, M. Kaplinghat,S. Garrison-Kimmel, J. Onorbe, and L. A. Moustakas,“Cosmological Simulations with Self-Interacting DarkMatter I: Constant Density Cores and Substructure,”Mon. Not. Roy. Astron. Soc. (2013) 81–104, arXiv:1208.3025 [astro-ph.CO] .[72] A. Sokolenko, K. Bondarenko, T. Brinckmann,J. Zavala, M. Vogelsberger, T. Bringmann, andA. Boyarsky, “Towards an improved model ofself-interacting dark matter haloes,” JCAP (2018)038, arXiv:1806.11539 [astro-ph.CO] .[73] J. Zavala, M. Vogelsberger, and M. G. Walker,“Constraining self-interacting dark matter with themilky way’s dwarf spheroidals,” Monthly Notices of theRoyal Astronomical Society: Letters no. 1, (Feb,2013) L20–L24. http://dx.doi.org/10.1093/mnrasl/sls053 .[74] S. Tulin and H.-B. Yu, “Dark Matter Self-interactionsand Small Scale Structure,” Phys. Rept. (2018)1–57, arXiv:1705.02358 [hep-ph] .[75] A. Fitts et al., “Dwarf Galaxies in CDM, WDM, andSIDM: Disentangling Baryons and Dark MatterPhysics,” Mon. Not. Roy. Astron. Soc. no. 1,(2019) 962–977, arXiv:1811.11791 [astro-ph.GA] .[76] A. Robertson, R. Massey, V. Eke, J. Schaye, andT. Theuns, “The surprising accuracy of isothermalJeans modelling of self-interacting dark matter densityprofiles,” arXiv:2009.07844 [astro-ph.CO] .[77] K. Bondarenko, A. Sokolenko, A. Boyarsky,A. Robertson, D. Harvey, and Y. Revaz, “From dwarfgalaxies to galaxy clusters: Self-Interacting DarkMatter over 7 orders of magnitude in halo mass,” arXiv:2006.06623 [astro-ph.CO] .[78] J. I. Read, O. Agertz, and M. L. M. Collins, “Darkmatter cores all the way down,” Mon. Not. Roy. Astron.Soc. no. 3, (2016) 2573–2590, arXiv:1508.04143[astro-ph.GA] .[79] J. I. Read, M. Wilkinson, N. Evans, G. Gilmore, andJ. T. Kleyna, “The Tidal stripping of satellites,” Mon.Not. Roy. Astron. Soc. (2006) 429–437, arXiv:astro-ph/0506687 .[80] E. Hardy in IAU Symposium. 206, Extragalactic StarClusters, D. Geisler, E. Grebel, and D. Minniti, eds.(San Fransisco: ASP), 62, 2002.[81] J. Strader, J. P. Brodie, D. A. Forbes, M. A. Beasley,and J. P. Huchra, “Spectroscopy of globular clusters inthe fornax dwarf galaxy,” The Astronomical Journal no. 3, (Mar, 2003) 1291–1297. https://doi.org/10.1086%2F367599 .[82] S. Martocchia, E. Dalessandro, M. Salaris, S. Larsen,and M. Rejkuba, “Is fornax 4 the nuclear star cluster ofthe fornax dwarf spheroidal galaxy?,” Monthly Noticesof the Royal Astronomical Society no. 4, (Jun,2020) 4518–4528. http://dx.doi.org/10.1093/mnras/staa1482 .[83] B. Hendricks, C. Boeche, C. I. Johnson, M. J. Frank,A. Koch, M. Mateo, and J. I. Bailey, “Evidence for achemical enrichment coupling of globular clusters andfield stars in the fornax dsph,” Astronomy &Astrophysics (Dec, 2015) A86. http://dx.doi.org/10.1051/0004-6361/201526996 .[84] J. Pfeffer, J. M. D. Kruijssen, R. A. Crain, andN. Bastian, “The E-MOSAICS project: simulating theformation and co-evolution of galaxies and their starcluster populations,” Mon. Not. Roy. Astron. Soc. no. 4, (Apr., 2018) 4309–4346, arXiv:1712.00019[astro-ph.GA] .[85] J. M. D. Kruijssen, J. L. Pfeffer, R. A. Crain, andN. Bastian, “The E-MOSAICS project: tracing galaxyformation and assembly with the age-metallicitydistribution of globular clusters,” Mon. Not. Roy.Astron. Soc. no. 3, (July, 2019) 3134–3179, arXiv:1904.04261 [astro-ph.GA] .[86] M. N. Rosenbluth, W. M. MacDonald, and D. L. Judd,“Fokker-Planck Equation for an Inverse-Square Force,”Phys. Rev. (1957) 1–6.[87] J. Binney and S. Tremaine,Galactic Dynamics: Second Edition. 2008.[88] S. Reddy, M. Prakash, and J. M. Lattimer, “Neutrinointeractions in hot and dense matter,” Phys. Rev. D (1998) 013009, arXiv:astro-ph/9710115 .[89] B. Bertoni, A. E. Nelson, and S. Reddy, “Dark MatterThermalization in Neutron Stars,” Phys. Rev. D (2013) 123505, arXiv:1309.1721 [hep-ph] .[90] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten,“Ultralight scalars as cosmological dark matter,” Phys.Rev. D95 no. 4, (2017) 043541, arXiv:1610.08297[astro-ph.CO] .[91] D. Lynden-Bell, “Statistical mechanics of violentrelaxation in stellar systems,” Mon. Not. Roy. Astron.Soc. (1967) 101. [92] P.-H. Chavanis, “Phase transitions in self-gravitatingsystems: Self-gravitating fermions and hard spheresmodels,” Phys. Rev. E (2002) 056123, arXiv:cond-mat/0109294 .[93] P. Chavanis, “Statistical mechanics of violent relaxationin stellar systems,” arXiv:astro-ph/0212205 .[94] P.-H. Chavanis, “Statistical mechanics andthermodynamic limit of self-gravitating fermionsinddimensions,” Physical Review E no. 6, (Jun,2004) . http://dx.doi.org/10.1103/PhysRevE.69.066126http://dx.doi.org/10.1103/PhysRevE.69.066126