Fluctuations of the gravitational field generated by a random population of extended substructures
Abstract
A large population of extended substructures generates a stochastic gravitational field that is fully specified by the function , which defines the probability that a tracer particle experiences a force within the interval . This paper presents a statistical technique for deriving the spectrum of random fluctuations directly from the number density of substructures with known mass and size functions. Application to the subhalo population found in cold dark matter simulations of Milky Waysized haloes shows that, while the combined force distribution is governed by the most massive satellites, the fluctuations of the tidal field are completely dominated by the smallest and most abundant subhaloes. In light of this result we discuss observational experiments that may be sufficiently sensitive to Galactic tidal fluctuations to probe the “dark” lowend of the subhalo mass function and constrain the particle mass of warm and ultralight axion dark matter models.
keywords:
Galaxy: kinematics and dynamics; galaxies: evolution.References
1 Introduction
Observations of the cosmic microwave background provide robust evidence for the existence of large amounts of nonbaryonic matter (e.g Planck collaboration 2014) that behaves like a perfect fluid on large scales (e.g. Peacock 1999). In theory, deviations from the perfectfluid model are expected to arise on scales comparable to the freestreaming length of the DM particle candidates. Below this scale, fluctuations of the power spectrum are heavily suppressed, de facto imposing a truncation at the low end of the halo mass function (e.g. Benson 2017 and references therein). However, for most DMparticle candidates the truncation falls well below the threshold of galaxy formation (White & Rees 1978; Bullock et al. 2000), which greatly complicates observational tests. For example, in cold dark matter (CDM) models made of WIMPs with masses GeV (e.g. Bertone et al. 2005) the cutoff is placed on planetsize mass scales, (e.g. Schmid et al. 1999; Hofmann et al. 2001; Green et al. 2005; Loeb & Zaldarriaga 2005; Diemand et al. 2005), which leads to one of the most striking predictions from the cold dark matter (CDM) paradigm, namely the existence of a very large number of selfgravitating “microhaloes” with subsolar masses devoid of baryons (i.e. ‘dark’).
The lack of visible matter in these objects together with their tiny masses makes the detection of DM microhaloes extremely challenging. Current observational efforts to test the presence of dark subhaloes in the Milky Way range from searching for gammaray annihilation signals (e.g. Ackermann et al. 2014; Bringmann et al. 2014) to detecting gaps in narrow tidal streams induced by close encounters with individual subhaloes (Ibata et al. 2002; Johnston et al. 2002; Yoon et al. 2011; Carlberg 2013; Erkal & Belokurov 2015; Ngan et al. 2016, Erkal et al. 2016; Bovy et al. 2017), with no unambiguous results to date. Stronglylensed galaxies provide complementary constraints on the clumpiness of dark matter haloes in the innermost region of galaxies (Koopmans 2005; Vegetti & Koopmans 2009; Li et al. 2013; Vegetti et al. 2014).
In general, current observational tests tend to probe the highmass end of the subhalo mass function. For example, the sizes of stream gaps produced by encounters with DM substructure in Milky Waylike galaxies are too small to be detectable for subhalo masses – (Carlberg 2012; Erkal& Belokurov 2015; Erkal et al. 2016; Bovy et al. 2017), with the added complexity that tidal heating by giant molecular clouds (GMCs) becomes nonnegligible on similar mass scales, (Amorisco et al. 2016). Perturbations of Einstein rings around lensed galaxies are also expected to be dominated by relatively massive subhaloes with (Li et al. 2016).
To the observational challenge of finding tiny invisible objects one must add the extreme dynamical range of CDM subhaloes. Collisionless body cosmological simulations that explore the properties of CDM substructures within Milky Waysize systems, e.g. the Aquarius Project (Springel et al. 2008), Via Lactea (Diemand et al. 2007) and GHALO (Stadel et al. 2009), have particlemass resolution –, which lies orders of magnitude above the mass scale of microhaloes.
To date it remains unclear to what extent the dearth of unresolved minihaloes may affect the predictions on the gammaray annihilation signals (e.g. Koushiappas 2009), and/or the dynamics of selfgravitating systems orbiting in the host galaxy. Similarly, the survivability of microhaloes in the tidal field of the parent galaxy is still poorly understood. On the one hand, flyby encounters with individual stars (Green & Goodwin 2007) and the tidal field of the Milky Way disc (D’Onghia et al. 2010; Errani et al. 2017) may wipe out a large fraction of the subhalo population with small orbital pericentres. On the other hand, the steep density profile of microhaloes (Anderhalden & Diemand 2013; Angulo et al. 2016) greatly increases their resilience to tidal mass stripping (Goerdt et al. 2007; Peñarrubia et al. 2010).
For those interested in the history of science the search for dark matter microhaloes bears unmistakable resemblance to the discussions on the ’reality of molecules’ during the 19th century. The existence of astronomical numbers of invisible molecules had been theoretically postulated by Avogrado in 1811 in order to explain the observed thermodynamical properties of gases. A heated epistemologic debate on whether molecules were real physical entities or just a useful mathematical construct lasted for nearly a century until Einstein (1905) introduced a revolutionary method for inferring the number and size of these objects from observations of the Brownian motion of particles suspended on the surface of a liquid. Given the very large number of degrees of freedom involved in the problem, Einstein entirely abandoned the classical Newtonian approach of solving the phasespace trajectories of individual molecules from deterministic equations of motion, focusing instead on a statistical description of the response of macroscopic objects to repeated interactions with a very large number of small particles.
This contribution follows a similar line of reasoning. Rather than following the dynamical evolution and the orbits of individual DM microhaloes in order to simulate their gravitational interactions with tracer particles, our goal is to construct a probability function that fully specifies the combined gravitational field generated by subhalo ensembles. As an illustration, in Section 2 we run numerical experiments that illustrate the stochastic nature of the combined gravitational force generated by a large subhalo population. Following up the work of Holtsmark (1919) we develop a simple statistical technique for deriving the spectrum of force fluctuations directly from the number density and the mass & size functions of extended substructures. In Section 3 we extend the analysis to the combined tidal field generated by these objects. Interestingly, we find that the spectrum of tidal forces induced by DM subhaloes is particularly sensitive to the lowmass truncation of the halo mass function, suggesting that the existence of ‘dark’ satellites could in principle be tested with observational experiments that measure fluctuations in the Galactic tidal field. Section 4 compares the force and tidal distributions expected in Milky Waylike galaxies for DM haloes made of cold/warm particles and ultralight axions. Section 5 summarizes our finding and discusses further applications of our method.
2 Stochastic fluctuations of the field
Let us assume that the gravitational acceleration experienced by a tracer particle at a distance from the galaxy centre can be expressed as
(1) 
where is the meanfield gravitational potential of the galaxy, and
(2) 
is the specific force induced by a set of selfgravitating substructures orbiting in the host galaxy. In systems where we expect to fluctuate stochastically along the phasespace trajectory of the particle, .
Solving (1) for an arbitrarilylarge population of substructures appears an impossibly difficult task given that the trajectory of the tracer particle is coupled with the equations of motion of each individual substructure. To attack this problem we shall follow a statistical method originally introduced by Holtsmark (1919) to describe the motion of charged particles in a plasma.
The crux of our analysis lies in the function , which defines the probability density that the test particle experiences a force in the interval . The simplest derivation assumes that the substructures are homogeneously distributed within a volume around the test particle, such that
(3) 
where is the Dirac’s delta function. Fourier transforming yields
Note that the last equality in Equation (2) implicitly assumes that the particles are randomly distributed within the volume or, equivalently, that the forces are spatially uncorrelated. The last integral can be rewritten as
which elevated to the th power becomes
where is the number density of substructures. It is useful to define the function
(5) 
such that the inverse Fourier transform of (2) becomes
(6) 
The above derivation can be easily generalized to inhomogeneous distributions of substructures by taking into account that the number density varies with radius, Equation (5) then becomes (Chandrasekhar 1941; Kandrup 1980; Chavanis 2009)
(7) 
where is centred at the location of the test particle, and is the number density profile. On scales the number density can be assumed to be roughly constant, , such that , which is typically known as the local approximation (see Appendix A).
The function contains all physical information on the number, distribution and masses of the substructure population. Unfortunately, it can be rarely expressed in an analytical form. Below we inspect a few notable exceptions with broad applications in Astronomy.
2.1 Pointmass particles
As a first step it is useful to review the analysis of fluctuations in the gravitational field induced by a random distribution of equalmass particles (see Kandrup 1980 for a formal description of this method).
The gravitational force exerted by a pointmass particle can be written as
(8) 
where is a unit vector. For simplicity, particles are assumed to be homogeneously distributed across the galaxy. Inserting Equation (8) into (5) and writing yields
where . Note that the resulting function is isotropically oriented in Fourier space, i.e. , which implies that the force distribution must also be isotropic, i.e. . Indeed, combination of Equations (6) and (2.1) leads to the wellknown Holtsmark (1919) distribution
The asymptotic behaviour of this function is studied in Chandrasekhar (1943) and reproduced here for completeness. In the weakforce limit we can approximate , such that
(11) 
We thus find that the effect of having an increasing number of particles at large distances combines with the declining Keplerian force such that in the weakforce limit becomes flat. Deriving the strongforce limit requires a number of nontrivial steps. Chandrasekhar (1943) finds
(12) 
It is important to stress that the largeforce behaviour is entirely dominated by the contribution of the nearest particle. This can be straightforwardly shown by writing the probability of finding the closest particle within the volume as (see Appendix A of Chavanis 2009)
(13) 
which peaks at , and making the transformation , where is the force induced by the nearest particle. For nearby objects, , Equation (13) reduces to , whereas . Hence,
(14) 
as in Equation (12).
It is useful to illustrate the above results by means of two simple numerical experiments. In experiment A we randomly generate independent samples of equalmass particles distributed homogeneously within a volume . For each particle ensemble we measure the total force . In experiment B we assign isotropic velocities to one random set of pointmasses generated in experiment A, such that the resulting orbital distribution is in equilibrium within a a harmonic spherical potential . The orbits of individual particles are integrated in the smooth potential for 10 dynamical times, where . Clearly, experiment A is designed to mimic the mathematical framework devised in Section 2, whereas experiment B provides a simplistic representation of the dynamics of pointmass particles in a galaxy.
Fig. 1 shows the combined force induced by equalmass particles with homogeneously distributed over a sphere of radius . Comparison between the upper and lowerleft panels shows clear differences between the force measured from random statistical samples and that generated by an equilibrium ensemble at different timesteps. The reason behind the mismatch can be traced to the different speed of fluctuations, , in experiments A and B. In dynamicalequilibrium systems such as experiment B this quantity is modulated by timescale of the orbital motions (see Chandrasekhar & von Neumann 1942 for a detailed discussion), whereas in experiment A the speed of fluctuations is random by construction. Yet, the right panel shows that the probability of measuring a force in the interval are statistically indistinguishable. Equation (2.1) (longdashed line) matches the numerical distributions, in agreement with the early body results of Ahmad & Cohen (1973) and del Popolo (1996). This is an important result, as it shows that the probability density fully specifies the stochastic perturbations induced by an ensemble of equalmass particles moving in an equilibrium orbital configuration in the host potential.
The specific force induced by a random distribution of particles upon a test star fluctuates around a mean value , which for ease of reference is marked with horizontal dottedlines in the left panels of Fig. 1. The mean force can be calculated analytically from Equation (2.1) by changing the integration variable to and integrating over the volume , which yields
where the third equality follows from the change of variable and the definition of the Gamma function, , whereas the fourth equality results from carrying the integration in complex space ]. For the experiments shown in Fig. 1 with , and we find . Note that the typical force has a magnitude , where is the average distance between particles.
A key feature of the Holtsmark distribution is the divergence of the moments for . This is due to the dominant contribution of the nearest particles to the largeforce, powerlaw tail of the distribution. Using Equation (12) it is straightforward to show that the variance of diverges in the strongforce limit as
(16) 
In practice this means that the maximum force experienced by the test particle can grow up to arbitrarilylarge values as the number realizations and timesteps increase in experiment A and B, respectively.
2.2 Extended substructures
Selfgravitating substructures with an extended matter distribution generate forces that (i) in general do not diverge at arbitrarilyclose distances, and (ii) asymptotically approach the Keplerian limit on scales larger than the size of these systems. The simplest modification of Equation (8) that accounts for this behaviour is
(17) 
The density profile associated to (17) can be readily found by solving Poisson’s equation, which yields
This profile corresponds to a Hernquist (1990) model with a scale length .
Let us follow the same steps as in Section 2.1 to derive the spectrum of perturbations induced by a random distribution of Hernquist (1990) spheres. Inserting Equation (17) into (5) and writing yields
where is a dimensionless quantity. The last equality follows from the change of variable , with the function being defined as
(19)  
The special functions FS, FC correspond to Fresnel sine and cosine integrals, respectively, whereas Si is the sine integral (see e.g. Press et al. 1992 for details). Note that in the limit one finds , , and , which reduces , thus recovering Equation (2.1).
Inserting (2.2) and (19) into (6), and following similar steps as in Equation (2.1) yields
(20) 
Fig. 2 shows numerical solutions to the integral (20) for different values of the scalelength . The first noteworthy result is the presence of a truncation of the force distribution at , such that for , which was to be expected given that (17) is not centrally divergent for .
Interestingly, the truncation does not generally correspond to the maximum force exerted by a single substructure, . To understand this result let us calculate the contribution of the single nearest substructure following the same steps as in Section 2.1. Recall that for nearby objects , so that Equation (13) reduces to , whereas from Equation (17) we have . It is straightforward to show that the transformation , where , yields
(21) 
which recovers Equation (14) in the limit . The function is plotted in Fig. 2 with thindotted lines. Note first that for Equation (21) greatly overestimates the probability to experience weak forces. The mismatch is caused by the break down of the closestparticle approximation at large distances, where .
The behaviour of in the strongforce limit is less clearcut. Comparison of Equations (20) and (21) at shows that the distribution of forces generated by the nearest substructure either over/underestimates the collective distribution (20) depending on whether the scalelength is smaller/larger than the typical separation between substructures, . Hence, there are two relevant regimes that define the largeforce behaviour of .

Rarefied regime arises when the separation between substructures is much larger than their individual sizes (). In this case large forces are dominated by the nearest substructure, such that .

Saturated regime arises when substructures overlap with each other (). Here large forces are dominated by the collective contribution of the ensemble rather than by individual objects, hence
The results shown in Fig. 2 indicate that the transition between these behaviours occurs on a scale . In what follows we limit our analysis to rarefied distributions, as an ensemble of selfgravitating, overlapping substructures is not dynamically stable.
Fig. 2 also shows that weak interactions become more likely as the substructure size increases, indicating that the mean specific force and the substructure scalelength are inversely correlated. To inspect this issue in more detail let us calculate by inserting (20) into (2.1). In contrast to the average force exerted by a random distribution of pointmasses, here cannot be expressed analytically owing to the nontrivial functional form of in (19). In the upper panel of Fig. 3 we plot numericallycomputed values of as a function of the substructure scalelength , which is given in units of the typical separation between substructures (). The mean force induced by substructures decreases monotonically as grows. For we find that that is a factor smaller than an ensemble of pointmasses. As expected, for the mean force converges asymptotically to the value given by Equation (2.1).
Crucially, the truncation of the force spectrum (20) at removes the divergence of the force variance discussed in Section 2.1. This implies that the maximum force experienced by a test particle orbiting in a clumpy medium does not grow up to arbitrarilylarge values insofar as the force generated by individual clumps does not diverge at arbitrarilyclose distances. To illustrate this result let us compute in the rarefied regime, where nearby particles dominate the largeforce tail of Equation (20), such that . From Equation (21)
which shows that the truncation of the largeforce spectrum implies a finite variance. As expected, the divergence found for pointmasses arises in the limit (). The close agreement between the variance derived from Equations (20) and (21) visible in the lower panel of Fig. 3 suggests that, as in the case of pointmass particles, nearby substructures dominate the largeforce, powerlaw tail of the force distribution.
2.3 Substructures with a mass function
The above results can be straightforwardly extended to ensembles of substructures with a mass distribution by writing the number density as
Powerlaw mass functions are of particular interest
(23) 
where is our mass unit, a normalization parameter and the powerlaw index. In general, the mass function is defined within a range , with . To normalize the distribution we impose the condition , where is the total number of substructures within a volume . Hence,
(24) 
where . Below we examine the impact of the mass function (23) on the substructure models outlined in Sections 2.1 and 2.2 separately.
2.3.1 Pointmass particles
Let us first calculate the function by taking into account that particles in the mass range contribute to the number density by an amount . Hence, from Equation (2.1) the total contribution of particles with masses between is
Thus, varying the parameters of the mass function simply modifies the numerical value of .
2.3.2 Extended substructures
Selfgravitating substructures in dynamical equilibrium have sizes that correlate with their masses. It is useful to quantify this dependency by introducing a size function , which determines the variation of the scaleradius as a function of . The powerlaw mass function seen in Section 2.3.1 can be readily modified in order to account for the masssize relation by writing
and
(27) 
where denotes Dirac’s delta function. Equation (2.2) then becomes
where , and is defined as
(29) 
with the functions given by Equation (19).
The relation between mass & scaleradius of extended substructures, , plays a key role in defining the largeforce behaviour of . In particular, it determines whether the strongest forces arise from the most massive objects or, in contrast, from the smallest, most abundant substructures. To study this issue in some detail let us adopt a powerlaw relation
(30) 
Given that in the rarefied regime the variance of the force distribution is dominated by the contribution of the nearest substructure, we can use Equations (2.2), (23) and (30) to derive an analytic expression for the force variance, which yields
Note that for ensembles with a broad mass range, , Equation (2.3.2) exhibits two welldefined regimes

Macrostructuredominated force distribution. This regime corresponds to a size function with index , such that , yielding a variance dominated by the largest objects.

Microstructuredominated force distribution. In this case the leastmassive substructures govern the largeforce behaviour of . This regime arises when the size function of substructures is sufficiently steep, , which leads to a force distribution whose variance scales as , thus diverging in the limit . This is an important result, as it shows that the probability of experiencing large force fluctuations strongly depends on the shape of the size function on mass scales .
2.3.3 Example: cosmological index
It is useful to illustrate the above results with an idealized representation of the cosmological body models studied in Section 4.1, which follow a powerlaw mass function (23) with an index (e.g. Springel et al. 2008), a fixed normalization parameter , and a lower mass limit that contains information on the quantum attributes of the dark matter particle candidates.
Notice first that a mass function with index within the interval contains a divergent number of substructures in the limit . Indeed, integration of (23) yields , which grows to arbitrarilylarge values as the minimum mass of substructures decreases.
The average force exerted by an ensemble of pointmass particles with is completely dominated by the most massive objects. This can be straightforwardly shown from Equation (2.3.1)
where we have substituted . Note that the mean force scales as independently of the lower limit of the mass function, . This is noteworthy given that, in principle, the number of light particles can become arbitrarily large.
As discussed in §2.3.2, the variance of the force distribution depends on the shape of the size function . Substituting in Equation (2.3.2) yields
(33) 
Hence, for the numerator in Equation (33) can be approximated as , which shows that the variance is dominated by the most massive objects, whereas for the variance diverges in the limit , indicating that the smallest substructures induce the largest fluctuations.
To gain further insight onto this result we show in Fig. 4 the force distribution arising from a population of extended substructures with a powerlaw mass function . To calculate we insert Equations (29) and (30) in (20) and integrate over . For simplicity we adopt . Also, for ease of comparison with previous Sections the volume is set to with , which leaves as the only free parameter. Each panel corresponds to a different value of the index in the size distribution (30). Adopting in the left panels reveals a scant sensitivity of to the lowermass limit, . Indeed, this is a macrostructuredominated distribution, where the mean and the variance of are largely determined by the most massive substructures in the ensemble. In contrast, for powerlaw indices and 1.5 (middle and right panels respectively) the largeforce distributions are microstructuredominated, exhibiting a strong dependence on the lower mass limit, . Notice in particular that for steep size functions () the largeforce tail of the force distribution approaches the pointmass behaviour in the limit , which according to Equation (16) leads to a divergent variance.
Interestingly, CDM subhaloes found in collisionless body simulations of structure formation exhibit powerlaw size functions with . For example, using the Via Lactea II models (Diemand et al. 2008) Erkal et al. (2016) find that provide the best fit to the masssize relation at redshift . Note, however, that cosmological simulations of Milky Waysize haloes typically resolve substructures with masses , with the lowmass limit only partially imposed by the finite mass resolution of the body runs (see van den Bosch 2017). Thus, (boldly) extrapolating these results orders in magnitude in mass, down to the mass scale associated with the streaming length of CDM particles (, e.g. SánchezConde & Prada 2014; Moliné et al. 2017 and references therein) suggests a negligible contribution of the smallest dark matter minihaloes to the force distribution . We will return to this issue in §4.
3 Stochastic tidal forces
The equations of motion that govern the phasetrajectory of a tracer particle within a selfgravitating system that moves through a clumpy medium can be written as
(34) 
where and are 33 tidal tensors evaluated at the centre of the selfgravitating potential . The smooth component has a form
(35) 
while the stochastic tidal tensor
(36) 
arises from the gradient in the combined tidal force induced by a set of substructures distributed across the external galaxy.
The tidal force induced by a single substructure can be written as
(37) 
To simplify the analysis it is useful to diagonalize the tensor by rotating the coordinates to a new frame , which corotates with the angular velocity of the substructure, . In the noninertial frame the effective tensor has a diagonal form, where is the centrifugal force component (see Renaud et al. 2011 for details). Unfortunately, the resulting field is not spherically isotropic, which greatly complicates the mathematical analysis of random forces outlined in Section 2. Yet, at leading order one can isotropize the tidal field by writing Equation (37) as
(38) 
where is the sum of eigenvalues of the effective tidal tensor , and is a unit vector irrelevant for our analysis (see §3.1). Equation (38) also neglects the Euler and Coriolis terms appearing in the noninertial rest frame, implicitly assuming that during the duration of the encounter the angular frequency remains constant and in a parallel direction to the galactocentric velocity , respectively. This approximation is less inaccurate in an impulsive regime, where the timescale of the tidal fluctuations is much shorter than the orbital period of the tracer particle in the potential .
With these approximations in place it is relatively straightforward to extend the analysis of Section 2 to systems experiencing stochastic tidal perturbations. In analogy with Equation (2) let us first define the vector
(39) 
which arises from the combined tidal field of the external substructures. From (38) and (39) the combined tidal force can be written as
(40) 
A tracer particle embedded in a sea of substructures isotropically distributed in space experiences a force that fluctuates in random directions according to a probability function , which is fully specified by the tidal distribution . Below we derive this function for isotropic ensembles of pointmass particles and Hernquist (1990) spheres.
3.1 Pointmass particles
Consider a vector
(41) 
where corresponds to the trace of the tidal tensor induced by an external Keplerian force (8) (see Equation (18) of Renaud et al. 2011).
The distribution can be expressed in an analytical form. Inserting (3.1) into (6) and integrating over yields
where is a dimensionless quantity.
The distribution (3.1) exhibits two welldefined asymptotic behaviours depending on the magnitude of . In the weaktides limit () Equation (3.1) becomes flat
(44) 
whereas in the strongtides limit () it follows a powerlaw tail
(45) 
As in Section 2.1, the largeforce behaviour is entirely dominated by the nearest particle. Indeed, from Equation (13) the probability to find the closest particle at a distance is , while (41) leads to . Hence, from (14) one has that
(46) 
which matches the strongtide limit (45) exactly.
It is straightforward to show that both the mean and the variance of the distribution diverge in the limit , that is when the distance to the closest particle becomes arbitrarily small. Thus, the divergence of the tidal moments is caused by the singular behaviour of (41) at .
3.2 Extended substructures
The trace vector associated with a Hernquist (1990) sphere can be straightforwardly derived from Equation (17) of Renaud et al. (2011), which yields
(47) 
As expected, Equation (47) approaches asymptotically the pointmass behaviour (41) at distances much larger than the size of the substructure, .
We can now follow the same procedure as in Section 2.2 to derive the distribution of tidal forces generated by random ensembles of Hernquist (1990) spheres. Inserting Equation (47) into (5) and writing yields
where is a dimensionless quantity, and the function is defined as