- Split View
-
Views
-
Cite
Cite
David R. Cole, Walter Dehnen, Justin I. Read, Mark I. Wilkinson, The mass distribution of the Fornax dSph: constraints from its globular cluster distribution, Monthly Notices of the Royal Astronomical Society, Volume 426, Issue 1, 11 October 2012, Pages 601–613, https://doi.org/10.1111/j.1365-2966.2012.21885.x
- Share Icon Share
Abstract
Uniquely among the dwarf spheroidal (dSph) satellite galaxies of the Milky Way, Fornax hosts globular clusters. It remains a puzzle as to why dynamical friction has not yet dragged any of Fornax's five globular clusters to the centre, and also why there is no evidence that any similar star cluster has been in the past (for Fornax or any other tidally undisrupted dSph). We set up a suite of 2800 N-body simulations that sample the full range of globular cluster orbits and mass models consistent with all existing observational constraints for Fornax. In agreement with previous work, we find that if Fornax has a large dark matter core, then its globular clusters remain close to their currently observed locations for long times. Furthermore, we find previously unreported behaviour for clusters that start inside the core region. These are pushed out of the core and gain orbital energy, a process we call ‘dynamical buoyancy’. Thus, a cored mass distribution in Fornax will naturally lead to a shell-like globular cluster distribution near the core radius, independent of the initial conditions. By contrast, cold dark matter-type cusped mass distributions lead to the rapid infall of at least one cluster within Δt = 1–2 Gyr, except when picking unlikely initial conditions for the cluster orbits (∼2 per cent probability), and almost all clusters within Δt = 10 Gyr. Alternatively, if Fornax has only a weakly cusped mass distribution, then dynamical friction is much reduced. While over Δt = 10 Gyr this still leads to the infall of one to four clusters from their present orbits, the infall of any cluster within Δt = 1–2 Gyr is much less likely (with probability 0–70 per cent, depending on Δt and the strength of the cusp). Such a solution to the timing problem requires (in addition to a shallow dark matter cusp) that in the past the globular clusters were somewhat further from Fornax than today; they most likely did not form within Fornax, but were accreted.
1 Introduction
The Fornax dwarf spheroidal (dSph) galaxy is the most massive undisrupted dSph satellite of the Milky Way (Walker et al. 2009). Like all dSphs it is dark matter dominated even in its central regions. It is unique among the undisrupted dSphs in having globular clusters (GCs); it has five, with three of them at a projected distance outside of the half-light radius (see Table 1). There is also evidence of two shell-like structures, which may be the remnants of a merger occurring more than 2 Gyr ago (Coleman et al. 2004, 2005).
One apparent paradox about these clusters is that, because they appear to orbit in a massive background of dark matter, they should be affected by dynamical friction which will cause their orbits to decay. Fornax's GCs are metal poor and very old, comparable to the oldest GCs in the Milky Way with ages of the order of a Hubble time (Buonanno et al. 1998, 1999; Mackey & Gilmore 2003b; Greco et al. 2007). During their lifetime it would be expected that they fall to the centre of Fornax and form a nuclear star cluster (Tremaine, Ostriker & Spitzer 1975; Tremaine 1976). However, no bright stellar nucleus is observed in Fornax, or in fact any other undisrupted dSph. This is known as the timing problem for Fornax's clusters because it seems highly improbable that Fornax's GCs would be observed just briefly before they fall into the core.
Several solutions to the timing problem have been proposed. Oh, Lin & Richer (2000) suggested two ideas: first, that a population of black holes transferred energy to the clusters through close encounters and secondly, that a strong tidal interaction between the Milky Way and Fornax could inject energy into their orbits. There is no observational evidence for a population of black holes in the centre of Fornax, while the currently observed proper motion indicates that the orbit of Fornax around the Milky Way never takes it closer than at present (Dinescu et al. 2004; Lux, Read & Lake 2010). Thus, both ideas appear to be disfavoured. Angus & Diaferio (2009) proposed that all but the most massive cluster could avoid sinking into the centre of Fornax if their current distance is much larger than projected but still within the tidal radius of Fornax of ∼1.9 kpc. However, this (i) requires special arrangement of the current projected positions and (ii) stills leaves a timing problem for the most massive cluster and is therefore not a complete solution. (Moreover, their analysis was based on Chandrasekhar's simple dynamical friction formula, which is not suitable for accurate estimates.)
Using numerical simulations and analytic arguments Goerdt et al. (2006) proposed that the current distribution of the Fornax clusters can be explained by the diminution of dynamical friction on the edge of a cored matter distribution which causes the clusters to remain outside the dark matter core radius. Dynamical reasons for this ‘core-stalling’ effect have been explored in Read et al. (2006b), Inoue (2009) and Cole, Dehnen & Wilkinson (2011). Support for this result was provided by Sánchez-Salcedo, Reyes-Iturbide & Hernandez (2006) who showed that a cored matter distribution in dwarf galaxies can significantly delay the infall times of the GCs (even if Chandrasekhar's simple dynamical friction formula is used).
Measuring and/or constraining the dark matter distribution in Fornax is interesting as a test of our current cosmological model. Collisionless cosmological simulations (which ignore the effects of baryons) predict a universal density distribution for dark matter haloes, with a central density cusp ρ ∝ r−γ where γ ∼ 1 (Dubinski & Carlberg 1991; Navarro, Frenk & White 1996b). If the dark matter distribution in Fornax is found to deviate strongly from this prediction, this could imply that baryons have an important dynamical role in shaping the central dark matter distribution in dwarf galaxies (e.g. Navarro, Eke & Frenk 1996a; El-Zant, Shlosman & Hoffman 2001, Read & Gilmore 2005; Goerdt et al. 2010; Cole et al. 2011; Pontzen & Governato 2012), or that we must turn to more exotic cosmological models (e.g. Tremaine & Gunn 1979; Hogan & Dalcanton 2000; Kochanek & White 2000; Strigari et al. 2006; Villaescusa-Navarro & Dalal 2011; Macciò et al. 2012).
The first evidence that dSphs may have a constant density core came from Kleyna et al. (2003) who found indirect evidence for a core in the Ursa Minor dSph. The Milky Way dSphs have been observed intensively in recent years, primarily because these systems are the most dark matter dominated known. They contain mostly intermediate or old stellar populations which are likely to be well mixed in the dark matter potential because star formation ceased many dynamical times ago. This in turn implies that they are ideal laboratories for studying the mass structure of their dark matter haloes. The intense observational effort means that there is a wealth of kinematical data available to form the basis for theoretical models of these systems. One line of approach has been based on the Jeans equations where a parametric light profile for the stars is assumed and a velocity dispersion profile is derived based on an underlying parametrized dark matter profile (Peñarrubia, McConnachie & Navarro 2008; Strigari et al. 2008; Walker et al. 2009). This approach has its drawbacks (see e.g. Amorisco & Evans 2011b). Most importantly, the degeneracy between the mass profile and the velocity anisotropy (which is poorly constrained), means that both cusped and cored density distributions are consistent with even the latest data. However, modelling the dSphs as two chemically distinct populations with different scalelengths, it appears that this degeneracy can be broken (Battaglia et al. 2008; Amorisco & Evans 2011a,b; Walker & Peñarrubia 2011); the results favour a cored density distribution in the two dSph galaxies best analysed to date: Fornax and Sculptor (but see Breddels et al. 2012).
In this paper, we follow the work of Goerdt et al. (2006) by examining what the current location of the GCs can tell us about Fornax's mass distribution. Our work improves on this previous analysis in several key respects: (i) we use several mass models for the underlying potential in Fornax that sample the full range consistent with the latest data, (ii) we use the latest data for Fornax's GCs as constraints on their phase space distribution and (iii) we run thousands of N-body models to sample the uncertainties in the cluster distribution and Fornax mass model.
This large search of the available parameter space allows us to address whether or not there are multiple solutions to Fornax's timing problem. To focus the discussion, we phrase the timing problem as a contradiction with either of the following two hypotheses.
- (i)
The Fornax GC system is in a near steady-state; consequently, none of the GCs should fall into the core of Fornax within a Hubble time.
- (ii)
Our present cosmic epoch of observing Fornax is not special; consequently, the system does not evolve significantly on a time-scale short compared to a Hubble time. In particular, within 1–2 Gyr none of the clusters should fall into the core of Fornax with high probability.
The first hypothesis can be justified by the assumption that the present state of Fornax's GC system is also representative of its past. This assumption, however, is not necessary, as one can easily think of alternative scenarios. The second hypothesis, on the other hand, is harder to avoid and is similar to the cosmological principle, though here expressed in terms of the epoch of observation rather than its vantage point and orientation. We will refer to contradictions with these hypotheses as the long-term and immediate timing problem, respectively.
This paper is organized as follows. Section 2 reviews the properties of the Fornax system relevant for our study, Section 3 details our modelling approach and Sections 4 and 5 present the simulations results and assess the probability that none of the five clusters will sink into the core of Fornax within either 1–2 Gyr or a Hubble time. Finally, in Section 6 we discuss the implications of our results and draw our conclusions.
2 The Fornax system
We summarize in Table 1 the most relevant data for our study and their origin.
2.1 The dSph
As already discussed above, Fornax is the most massive of the Milky Way's dSph [except possibly for disrupted objects such as Sagittarius (Sgr) dwarf] and unique amongst (undisrupted) dSph to host a GC system. The very fact that Fornax can hold on to a GC system implies that Galactic tides cannot be strong enough to pull these clusters off. This in turn requires that its Galactic orbit never carries Fornax too close to the Milky Way, where the tidal forces become exceedingly strong.
Lux et al. (2010) estimate the peri-galactic radius of Fornax, based on its observed position, distance, radial velocity and proper motion to be 100–130 kpc, which indeed is only slightly smaller than its current distance of about 140 kpc (see Table 1).
Based on this result, we estimate the tidal radius for Fornax using the method of Read et al. (2006a), where the dSph is modelled as a spherical satellite orbiting the Milky Way represented by a Hernquist (1990) model. We then solve equation 7 of Read et al. (2006a), which accounts for the orbit of the satellite about the host and the orbits of the stars within the satellite. We find a tidal radius of 1.8–2.8 kpc, based on the range of masses for Fornax given in Table 2 and using the extremal values for the orbital data taken from Lux et al. (2010) and a total (extended) mass for the Milky Way of 1 - 2 × 1012 M⊙.
2.2 The globular clusters
Our principal sources for GC data are those published by Mackey & Gilmore 2003a,b and Greco et al. (2007) who have carried out thorough surveys of the Fornax globular clusters. For our purposes the main data required are the clusters’ masses, sizes, three-dimensional positions and velocities.
The best estimates for these quantities are given in Table 1. The values for the core radius rc of each cluster are based on the surface brightness profiles calculated in Mackey & Gilmore (2003a). These are Elson, Fall & Freeman (1987, EEF hereafter) models and the King (1962) model core radius rc is related to the EFF scale parameter a by rc = a(22/γ −1 1)1/2, where γ is the power-law slope of the surface brightness at large radii.
Fig. 1 plots the distribution of the five Fornax GCs in mass M and projected radius Rp from the dSph. There are two interesting observations to be made from this figure. First, the radial distribution of clusters is consistent with that of the stars within Fornax: there is about half of the total cluster light within the stellar half-light radius of 668 pc. There are two possible interpretations of this. Either it is a coincidence, or the formation histories of the clusters and the Fornax galaxy are closely related, in particular, the clusters formed within the same entities as the stars.
Secondly, there is a weak correlation between M and Rp. In particular, the lightest cluster is furthest away from Fornax and the heaviest is the second closest. The remaining three are about equally massive and cover a spread of Rp. This correlation is in the sense expected from mass segregation such as driven by dynamical friction.
3 Modelling Approach
The basis for our approach is to take the most up-to-date observations of Fornax's GCs and combine these with plausible mass models consistent with the latest kinematic data for Fornax's stars. We then create, for each Fornax mass model, initial conditions for the Fornax GC system, which are consistent with the relevant observations, and evolve them for 10 Gyr into the future.
3.1 Mass models for Fornax
The stellar distribution functions are calculated numerically following the approach of Gerhard (1991) and Saha (1992) allowing various velocity anisotropy profiles. As in the earlier work (Kleyna et al. 2002; Wilkinson et al. 2002), the models are compared to the data on a star-by-star basis.
The fourth model, large core (LC), was motivated by the recent work of Walker & Peñarrubia (2011). These authors applied a non-parametric statistical modelling technique to two chemically distinct stellar populations within Fornax to define the enclosed mass at the half-light radii of the two populations. The resulting model possesses an LC with near-constant density and γ ≈ 0.1 for r ≲ 500 pc.
The radial profiles of density, enclosed mass and logarithmic density slope of the four mass models are shown in Fig. 2 for comparison. These models cover a wide range of inner density slopes, including shallow profiles, such as those suggested by Gilmore et al. (2007) based on observations of dSph galaxies, but also a SC, such as those predicted by cosmological simulations (Dubinski & Carlberg 1991; Navarro et al. 1996b). Note, however, that our models represent the overall mass distribution of Fornax including both the stars and the dark matter.
Fig. 3 shows the stellar velocity dispersion for our three MCMC-based Fornax models (curves) plotted together with the observed velocity dispersion as measured by Walker et al. (2007). The model velocity dispersions were calculated assuming an ergodic distribution function with a Plummer (1911) density profile with core radius 668 pc for the stellar component, but the relevant halo model for the underlying mass distribution. In view of the fact that these simple ergodic models have not been fitted to the velocity-dispersion data (apart from the assumed mass models), they provide a surprisingly good description of these data. (This and the similarity between the model predictions are exactly the reason why inferring the mass profile from data like these is hardly possible.)
Modelling the stellar velocity dispersion of a single stellar population as described above is not appropriate for the LC model. This model has been normalized to roughly agree with the estimates for the enclosed mass derived from two different tracer populations by Walker & Peñarrubia (2011) and plotted in Fig. 2.
3.2 Modelling the globular cluster system
This is done in such a way that for clusters GC2–5 the line-of-sight velocity matches the observed value (which may be considered a prior for our sampling).
As the spatial distribution of the clusters is consistent with that of the stars in Fornax, it seems reasonable to base our kinematical parameters for the clusters on the observed stellar kinematics. The measured stellar velocity dispersion is approximately flat over the range of radii observed (Walker et al. 2007; Łokas 2009), and we use σ = 10.5 km s−1. For the velocity anisotropy we assume β = −0.33 as suggested for the stars (Łokas 2009), i.e. a mild tangential bias, which gives σr ≈ 9.5 km s−1 and σt ≈ 11 km s−1.
With this modelling approach, in particular the wide range of initial radii sampled, we generate many different orbits for the GCs, covering the complete range of all possible orbits for them. By allowing such a wide distribution of cluster orbits, we can explore the effects of possible narrow choices for the cluster distribution function, such as preferably circular orbits, afterwards by restricting our analysis correspondingly.
3.3 N-body simulations
For each Fornax mass model, we ran 700 N-body simulations, each with different initial conditions for the five clusters.
To generate the N-body initial conditions for the Fornax mass models, we sample positions from the density (2) and velocities from self-consistent ergodic distribution functions ƒ(E), which only depend on the specific orbital energy E, thus giving everywhere isotropic velocity distributions. The forces between particles representing Fornax are softened with softening length ε = 10 pc.
The mass ratio of the lightest background particle in all our models to the lightest cluster is 1:1680 (in model WC) and that of the heaviest particle to the lightest cluster 1:62 (in model LC). The clusters orbit mainly within the high-resolution region of our models defined by the volume where the resolution is better than that would have been achieved with the same number of particles of constant mass (within approximately 2 kpc).
The simulations are performed using the publicly available N-body code gyrfalcon, which uses Dehnen's (2000, 2002) algorithm for force approximation. The total energy conservation was typically a few parts in 104.
4 Raw Simulation Results
As detailed at the end of Section 3.2, our simulations cover a wide range of initial cluster orbits, some of which may not be very realistic. In this section, we ignore any implications of our sampling of the initial cluster orbits and simply consider the individual simulations on their own merit.
For each of the four mass models of Table 2, Fig. 4 plots the distributions of each cluster orbit in initial and final rapo at t = 2 and 10 Gyr in the left-hand and right-hand subpanels, respectively.
4.1 Orbital decay after 2 Gyr
After 2 Gyr (left-hand subpanels in Fig. 4) rapo is significantly less than initially for most orbits with initial rapo < 2–3 kpc, except for model LC, which we discuss separately in Section 4.3.
However, both GC1 (red) and GC5 (cyan) rarely sink substantially, and GC1 only ever falls into the centre of Fornax for the most cusped model SC. This is not surprising since GC1 is not only by far the lightest of the five clusters (see Table 1), but also the one furthest away (in projection) from the centre of Fornax and hence less likely to pass through high-density regions. As such it requires most dynamical friction to fall into Fornax, but will suffer least, since drag force ∝mass2ρ. Though GC5 is five times more massive, it has the second greatest projected distance from the centre of Fornax and hence also suffers significantly less friction than the other clusters for most orbits sampled.
The remaining three clusters GC2–4 are all dragged inwards when initially rapo < 2–3 kpc, presenting the Fornax timing problem. For all of these clusters, orbits with initial rapo ≲ 1 kpc show similar effects of dynamical friction, presumably because orbits with apo-centres as small as that spend sufficient time in high-density regions to suffer substantially from dynamical friction.
Since Rp ≤ rapo, an initial rapo ∼ 1 kpc is the absolute possible minimum for GC2, which is observed at Rp = 1.05 kpc. Consequently, GC2 only rarely falls in as much as GC3 and GC4.
Apart from the initial rapo, the infall of these clusters depends most strongly on the inner density profile of the background mass distribution. Model SC shows the greatest effect of dynamical friction on the clusters. For GC3 and GC4, a significant proportion of the simulations find their instantaneous apo-centres inside 30 pc. GC2 does not show such a marked effect but a number of simulations already have an apo-centre inside 100 pc.
By contrast, model WC shows significantly reduced dynamical friction. GC3 again is most affected. However, even in the most extreme cases, the apo-centres have decayed to no less than ∼100 pc from the centre. In this model, the logarithmic density slope γ (equation 3) was initially only 0.1 at 100 pc, which implies a near-flat density profile within this radius. As expected, model IC is intermediate between models WC and SC.
4.2 Orbital decay after 10 Gyr
After 10 Gyr, the trends shown at 2 Gyr are amplified. Four of the GCs fall into Fornax for most simulations. Only GC1 still shows not much evidence of migrating to the centre of Fornax. All clusters have a small fraction of simulations where they remain at large rapo. This requires the initial rapo to be large as well. In general, rapo decreases at least slightly, even when initially large (the few significant increases of rapo are caused by cluster encounters).
For the most cusped model SC and to some degree for model IC too, GC3, GC4 and even GC2 can obtain values for rapo down to close to their softening length of 5 pc, i.e. they are essentially at the very centre of the galaxy. For model WC, this is not the case, i.e. none of the clusters can reach the very centre of the galaxy in this case. However, the vast majority of simulations do not obtain such very small rapo.
We observe that the distributions of rapo after 10 Gyr are quite similar between models WC, IC and SC (in particular the latter two), when GC3 has settled at rapo < 100 pc for most of our simulations, in fact all simulations with initial rapo < 2.8 kpc, and the distributions for rapo of GC4 are remarkably similar. For the most massive GC3, this similarity reflects the fact that the cluster has essentially sunken to the very centre of the galaxy.
However, this cannot explain the similarities for GC4, which we think is caused by a similarity of the background mass profiles after 10 Gyr. The action of cluster dynamical friction causes dynamical heating of the background particles dragging the clusters inwards. This in turn erases the initial central density cusp (Read et al. 2006b; Cole et al. 2011). Goerdt et al. (2010) discuss the formation of a density core in this way and provide an empirical formula for the expected core size created as the clusters fall to the centre. For model WC this predicts a core radius of 262 pc which is larger than our stalling radii but agrees within a factor of a few. However, this formula fails to predict the behaviour of the clusters in the IC and SC models as it gives stalling radii of 176 and 100 pc, respectively.
In Fig. 5, we plot the density profiles of the final models for cases with and without cluster infall. For the SC and IC models, the central density profiles are significantly reduced and in fact have become much more similar to each other. However, only for model SC is this reduction slightly stronger when clusters have reached the core of Fornax. In all cases, the density still keeps on increasing inwards, and hence does not really correspond to a constant-density core in the classical sense1.
4.3 The large core model
The LC model shows very unusual behaviour. After 2 Gyr, all GCs have apo-centres closely clustered together with a strong peak at ∼1 kpc. This becomes even more pronounced at 10 Gyr. Detailed examination of the cluster orbits in individual simulations shows two interesting behaviours. First, orbits with initial rapo ≳ 900 pc decay move in quite rapidly (mostly in less than 2 Gyr) to rapo ∼ 900 pc, where they stall. This confirms the work of Goerdt et al. (2006) and Read et al. (2006b) which showed that massive satellites orbiting outside of an harmonic density core stall at the edge of the core. This behaviour is believed to be due to the reduction of dynamical friction due to the resonant effects of particles in the harmonic core.
Secondly, any cluster which has an initial orbit within the harmonic core moves out to the edge of the core. Fig. 6 shows the evolution of the orbit of GC4 in one of our simulations. As can be seen the cluster orbit absorbs energy as it expands. The overall conservation of energy for the simulation is not affected by this change and energy is conserved to approximately three parts in 104 during the simulation. This behaviour is unexpected and we do not believe that it has been reported previously, though there is some evidence for orbital radii expanding again after falling in at the edge of harmonic cores in what has been called the ‘kickback effect’ (Inoue 2009; Goerdt et al. 2010). We will discuss this further in section 6.
4.4 Summary
The raw empirical results from our simulations can be summarized as follows.
- (i)
Cluster orbits with large initial rapo are not significantly affected by dynamical friction, because these orbits spend no or too little time in high-density regions, where frictional drag is exerted. However, most simulations with initial rapo less than the current tidal radius of Fornax suffer from dynamical friction.
- (ii)
Cluster GC3 is most likely affected by dynamical friction, followed by GC4 and GC2, while clusters GC1 and GC5 are least likely affected after 2 Gyr. This ordering is expected from the masses and initial projected radii of the clusters.
- (iii)
For all except model LC, cluster GC3 always reaches the core of Fornax within 10 Gyr (unless its initial orbit was unrealistically large with rapo > 2.8 kpc beyond the present-day tidal radius of Fornax), constituting the long-term timing problem.
- (iv)
The effect of dynamical friction at 2 Gyr is increasing with the central mass density from model WC to SC, as expected from Chandrasekhar's dynamical friction formula.
- (v)
The effect of dynamical friction after 10 Gyr is more similar for the three halo models with WC to SC than after 2 Gyr. This similarity can be understood, at least qualitatively, by the stalling of dynamical frictions as consequence of core formation.
- (vi)
Model LC shows no effect of dynamical friction, but rather the opposite: ‘dynamical buoyancy’, when the clusters are pushed out of the core. Like dynamical friction, this effect appears strongest for the most massive cluster.
In particular, for a cold dark matter (CDM)-type SC model, almost all simulations with initial rapo < 2.8 kpc (Fornax's tidal radius) suffer significantly from dynamical friction within 10 Gyr, or even within only 2 Gyr. This reflects the Fornax timing problem and is in contradiction with the claims of Angus & Diaferio (2010) (as discussed in Section 1), possibly because our sampling discourages very nearly circular orbits, but the usage of Chandrasekhar's simple formula by these authors certainly plays a role too.
5 The probability of cluster sinking
Our results presented in the previous section and Fig. 4 demonstrate that the more massive clusters GC3 and GC4 will fall into the core of Fornax within 2 Gyr or less, unless either the mass distribution of Fornax has a core or the cluster has an initial orbit with large apo-centric radius rapo. This latter scenario requires that the observed cluster position be either near peri-centre (special orbital phase) or has intrinsic radius much larger than projected (special projection geometry). Either is atypical, implying that this scenario is inherently unlikely. In this section, we try to quantify just how unlikely.
Since we know neither the current orbital phase nor projection geometry for any of the GCs, our most sensible approach is to marginalize over both, assuming uniform distributions. However, as we will see, our sampling of initial cluster position and velocity did not generate a uniform sampling in these quantities. As a consequence, any quantitative statements based on the raw results about the probability of cluster infall will be biased. In order to eliminate this bias, we must emulate a uniform sampling of initial orbital phases and projection geometries.
In Fig. 7, we plot for each simulated cluster and each halo mass model the instantaneous apo-centric radius rapo after 2 and 10 Gyr against p(R ≤ Rp|orbit). For all mass models, there is a clear correlation between p(R ≤ Rp|orbit) and rapo at later times in the sense that larger rapo(t > 0) are achieved only when the initial projected radius was relatively small for the initial orbit. This makes perfect sense: for small p(R ≤ Rp|orbit) the orbit spends most of its time at r > Rp with little dynamical friction.
As is evident from Fig. 7, the distributions of p(R ≤ Rp|orbit) from our simulations are not uniform: there are many more simulations with small p(R ≤ Rp|orbit), in particular for clusters with small Rp, such as GC3 and GC4.3 This non-uniformity is simply a consequence of our sampling procedure, which favours orbits with rapo ≫ Rp (in particular for clusters with small Rp, such as GC3 and GC4) resulting in non-uniform sampling of orbital phase and projection geometry and hence introducing a bias.
Since rapo(t > 0) is strongly correlated with p(R ≤ Rp|orbit), we can read off Fig. 7 for the unbiased probabilities of infall into Fornax. For example, in model SC rapo(t = 2 Gyr) > 200 pc for GC4 (magenta) in almost all simulations with p(R ≤ Rp|orbit) < 0.1 but hardly for any simulation at larger p(R ≤ Rp|orbit), implying that cluster GC4 will have rapo > 200 pc after 2Gyr with ∼10 per cent probability in halo model SC.
A more quantitative procedure for removing the bias of our non-uniform sampling of orbital phase and projection is to give each simulated cluster orbit a weight such that the weighted distributions of simulated orbits have, or at least are consistent with, a uniform sampling. Obviously, the weight to choose is simply the reciprocal of the actual frequency ƒ(p(R ≤ Rp|orbit)) of fractions p(R ≤ Rp|orbit) for all simulations of the same cluster in a particular halo model. The probability for, say rapo(t = 2 Gyr) > 200 pc for GC4, is then obtained as the weighted fraction of simulations that obtain rapo(t = 2 Gyr) > 200 pc for this cluster.
Combining this for all clusters, we estimate the probability that none of the five clusters when initially on orbit with rapo < 2.8 kpc (our upper limit for the tidal radius of Fornax) will fall into Fornax, in the sense that rapo < 100 pc or <200 pc, within 1, 2 or 10 Gyr. The results for each halo model are given in Table 3 (except for model LC when no cluster ever falls into Fornax). Columns 2–4 give the probabilities based on all of our simulations. If we restrict the analysis to high or low eccentricities (columns 5–7 and 8–10, respectively) or to orbits with initial rapo < 1.8 kpc (not given in Table 3), the results are hardly altered.
In other words, these probabilities depend mainly on the halo model and only weakly on the distribution function of the cluster orbits. This is a direct consequence of the tight relation, seen in Fig. 7, between p(R ≤ Rp|orbit) and rapo at later times. This result implies that we do not need to investigate further the implications of different distribution functions for the initial cluster orbits.
6 Conclusions
The Fornax galaxy is unique among the Milky Way dSphs in having five surviving GCs. These clusters are metal poor and very old – comparable with the oldest GCs in the Milky Way. There is a well-known timing problem for these clusters: over a Hubble time or even a much shorter time-scale dynamical friction should drag them to the centre of Fornax, where they would form a nuclear star cluster. Yet no such stellar nucleus is observed in Fornax or, in fact, in any other dSph galaxy. The one possible exception to this is the Sgr dSph which has the GC M54 at its core (Bellazzini et al. 2008). This dSph is currently only ≈ 16 kpc from the galactic centre and strongly affected by the tidal field of the Milky Way. It seems likely that Sgr has lost much of its mass to the Milky Way and was once much more massive, of the order of the mass of the Large Magellanic Cloud (Łokas et al. 2010). If so the effects of dynamical friction would be correspondingly greater in a more massive galaxy and the likelihood of a GC falling to its centre is much greater. Recent work has also indicated that there is at least a possibility that M54 may not be at the core of Sgr but at ≈ 2 kpc in the foreground (Siegel et al. 2011).
In this study, we have extended previous work on what the current location of Fornax's GCs can tell us about Fornax's mass distribution. We explored four different mass models for the underlying potential in Fornax, we used the latest data for Fornax's GCs as constraints on their phase space distribution and we ran thousands of N-body models to sample the uncertainties in the cluster position and velocity distribution for each of our five mass models. This large search of the available parameter space allowed us to hunt for viable solutions to Fornax's timing problem.
6.1 Caveats
Our models all assume a spherical mass distribution for Fornax, yet the stars are observed to be ellipsoidal in projection, implying that the true intrinsic distribution is aspherical, most likely triaxial. The parameter space of such configurations is considerably larger (and the freedom in cluster orbital projections is smaller), potentially allowing more solutions to the timing problem. However, previous studies suggest that, while dynamical friction is on average stronger on box orbits than on loop orbits(Capuzzo-Dolcetta & Vicari 2005), triaxiality has no significant overall effect on the strength of dynamical friction (Cora, Vergne & Muzzio 2001; Sachania 2009). This can be understood as cancellation of two opposing effects. First, because angular momentum is not conserved, there is no barrier for eccentric orbits to reach arbitrary small radii and hence high densities and strong drag. Secondly, also because angular momentum is not conserved, successive peri-centric passages have different radii and densities such that an orbiting cluster only rarely suffers very high drags.
Another minor caveat of our modelling is our ignorance of the tidal field of the Milky Way potential. In principle, it would be trivial to have our N-body models orbit in a static model for the Galactic potential. However, such a procedure would add even more only weakly constrained parameters without adding significant benefits. With our existing models, we take the Galactic tidal effects into account by discounting any simulated cluster orbits with apo-centric radius beyond (best estimates of) Fornax's tidal radius.
Finally, by modelling each GC as a massive particle, we have ignored their inner dynamics and tidal interaction with Fornax. Mass-loss rates for these clusters, however, are likely to be too small to significantly reduce their mass (Goerdt et al. 2010). Tidal disruption is arguably only relevant for GC1, as discussed in Section 6.3 below.
6.2 Solutions to the timing problem
Our simulations demonstrate that for our normal mass models (WC, IC or SC) for Fornax the infall of clusters GC3 and GC4 within a Hubble time is unavoidable and the infall of all clusters except GC1 most likely. This constitutes the long-term timing problem in the sense of hypothesis (i) from Section 1, which is only really a problem if one assumes that the present GC distribution of Fornax is representative of the distant past.
Only our LC model avoids this problem completely. Rather than dynamical friction, this model shows ‘dynamical buoyancy’, when clusters are pushed out of the core. Following Tremaine & Weinberg (1984), this is likely caused by this model's inverted distribution function, when over a significant range of orbital energies E. Such models4 are generally thought to be unstable (Binney & Tremaine 2008, section 5.5). Hence, it seems unlikely that such a model emerges from any reasonable formation mechanism. However, our N-body models appear to be stable, and some further research is needed to clarify the cause for ‘dynamical buoyancy’, its relation to the core stalling effect (reported by Read et al. 2006b) and its prospects in natural systems. We have considered the possibility that this effect is due to numerical noise. It is possible that the LC model may not be thoroughly tested by our numerical convergence work (Appendix A) due to its unusual nature. We have subsequently run two simulations with identical initial conditions using the LC model and one GC which shows the ‘dynamical buoyancy’ behaviour. One has the fiducial resolution of our main work, 2 million particles, and the other has five times better resolution, 10 million particles (see Fig. 8). We find the evolution of the GC's orbit to be almost identical in the two cases and the energy of the GC is the same to within 1 per cent as the orbit moves out. This strongly suggests that we see a real effect, which we will investigate further in future work.
For cusped mass models (IC and SC), such as predicted by CDM cosmogony via dissipationless simulations, clusters GC3 or GC4 will sink into the centre of Fornax within 1–2 Gyr with ∼90 per cent probability (in the sense that solutions where this does not occur cover only ∼10 per cent of the possible orbital phases and projections). In fact, we estimate the probability that no cluster obtains rapo ≲ 100 pc within 2 Gyr to be no more than 2 per cent in Table 3. This is largely independent of the assumed GC orbital distribution function and constitutes the more severe immediate timing problem in the sense of hypothesis (ii) from Section 1. It implies that if Fornax indeed has a cusped density profile, then our cosmic epoch of observation is necessarily very special.
For a shallow cusp model (WC), only the most massive cluster GC3 has a 90 per cent probability to reach, within 1–2 Gyr, rapo ≲ 200 pc significantly less than its current projected radius of Rp = 430 pc. For this model, no cluster can reach rapo ≲ 100 pc within 2 Gyr and the chances that no cluster will reach rapo ≲ 200 pc within 1 or 2 Gyr are, respectively, 92 and 32 per cent (see Table 3). These numbers are not unlikely and hence avoid the long-term timing problem.
6.2.1 A steady-state solution
Our finding thus suggest two possible solutions to the timing problem. First, Fornax has an LC (perhaps between models LC and WC) and dynamical friction is slow or has stalled a long time ago. In this case, Fornax may have been on its current orbit for a Hubble time with its GC system hardly evolving. In this case, the consistency of the cluster distribution with the stellar distribution (as discussed in Section 2.2) cannot be a coincidence, but hints at a common formation scenario.
6.2.2 An evolving solution
In the second solution Fornax has a small core or shallow cusp (as model WC) and dynamical friction is still ongoing, albeit slowly enough that the absence of a central nucleus in Fornax (or in fact any other undisrupted dSph) is perfectly plausible. In this case, the clusters must have been further away from Fornax in the past than today, and the current (weak) consistency of their distribution with that of the stars is just a coincidence. Also, a Hubble time ago the clusters most likely were more than the current tidal radius of 1.8–2.8 kpc away from Fornax. This in turn required that Fornax did not orbit the Milky Way for a Hubble time on its present orbit. However, even a simple adiabatic evolution of Fornax's orbit may be sufficient to solve this problem. For example, for a slowly growing Galaxy with always flat rotation curve the peri-centric tidal radius of Fornax evolves like (even when neglecting mass-loss of Fornax due to Galactic tides).
This second solution appears more natural and also fits with the weak indication of mass segregation, as would be induced by dynamical friction, in the current mass–radius relation (see Fig. 1). However, this model implies that the GCs have not formed within Fornax, but are most likely accreted. One may, of course, consider these two solutions as the extreme ends of a continuity of solutions with various degrees of cusp strengths and hence dynamical friction effects.
6.3 The case of GC1
An interesting aspect relates to the cluster GC1. Peñarrubia, Walker & Gilmore (2009) demonstrated that, uniquely of all Fornax's GCs, GC1 would be tidally disrupted if it fell to the centre of Fornax. While we find that GC1 does not sink to the centre of Fornax in almost all of our models, this still leaves us with a puzzle. Why should the one cluster vulnerable to tides be on an orbit where it would hardly ever suffer disruption? For our steady-state solution to the timing problem above, this puzzle can be resolved by the postulation that Fornax once had a richer GC system and we only see the survivors. Such survivors are either massive enough or on remote orbits to avoid tidal disruption. However, a high (initial) frequency of GCs (of ∼104-5 M⊙) appears rather implausible for a small galaxy like Fornax.
For the evolving solution to the timing problem, on the other hand, the fact that GC1 would be disrupted poses no problem at all. In this picture, low-mass clusters, such as GC1, would not be dragged down much, and there is no need to postulate a large early population of clusters.
Acknowledgments
Research in Theoretical Astrophysics at Leicester is supported by an STFC rolling grant. JIR would like to acknowledge support from SNF grant PP00P2_128540/1.
This research used the ALICE High Performance Computing Facility at the University of Leicester. Some resources on ALICE form part of the DiRAC Facility jointly funded by STFC and the Large Facilities Capital Fund of BIS. We also thank Matt Walker for providing the stellar velocity dispersion data to compare to our mass models in Fig. 3.
References
For model WC, the density actually increases slightly at r < 100 pc. This is presumably caused by a slight instability of this model, which has for some range of specific energies E and hence may be unstable (see Binney & Tremaine 2008, section 5.5).
In addition to p(R ≤ Rp|orbit) one may also utilize the fraction p(υz ≤ υlos|Rp, orbit) of the line-of-sight velocity to be less than observed. This fraction should also be uniformly distributed independently of p(R ≤ Rp|orbit), and thus allow an additional independent constraint. However, its computation is more involved and it seems unlikely that much would be gained from it, mainly because of the relatively large uncertainty of the observed υlos.
For model LC, no simulation has small p(R ≤ Rp|orbit) for clusters with Rp > 1 kpc. This is a consequence of the larger mass of this model which implies that our sampling did not obtain any nearly unbound orbits with large initial rapo.
For self-gravitating systems, an inverted distribution function typically occurs if the transition between a near-constant density core and steep power-law decay is fast, i.e. when η in equation (2) is larger than ∼2.
This follows from the probability density for the polar angle θ ∈ [0, π] and with R = r sin θ.
Appendix
Appendix A: Numerical convergence
In order to ensure that our simulations do not suffer from numerical noise we ran simulations with two different mass models, one cusped and one cored; each with two different orbits, one circular and one eccentric; and each of these at three different resolutions: with N = 4 × 105, 106 and 4 × 106. We then compared the evolution in each case of one cluster over 10 Gyr.
The evolution of the orbital radius of a single cluster moving on an eccentric orbit in model WC (see Table 2) is shown in Fig. A1. It can be seen that orbital evolution is very similar for all resolutions. In particular, the decay of the orbit follows the same time-scale, with the time and radius of the first-stalling of the cluster being the same. It has been shown that the two-body noise in a simulation can cause the cluster orbit to precess and cause artificial decay of the orbit once in the core (Read et al. 2006b). We selected the above combination of orbit and density profile precisely because Read et al. (2006b) showed that convergence is most difficult for an eccentric orbit in a cored halo. This is the case where numerical friction caused by orbit precession has the largest effect on orbital decay. The simulations shown above give a strong indication that such effects are not significant at even lower resolutions than the one used for the main body of this work. Our simulations are well converged.
Appendix B: The probability for R ≤ Rp
The fraction of orbital phases and projections for which R ≤ Rp on a given orbit is identical to the probability for R ≤ Rp, when R is computed for that orbit at random orbital phase and projection.