Stochastic substructure forces

# Fluctuations of the gravitational field generated by a random population of extended substructures

Jorge Peñarrubia
Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
jorpega@roe.ac.uk
###### 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 Way-sized 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” low-end of the subhalo mass function and constrain the particle mass of warm and ultra-light axion dark matter models.

###### keywords:
Galaxy: kinematics and dynamics; galaxies: evolution.

## 1 Introduction

Observations of the cosmic microwave background provide robust evidence for the existence of large amounts of non-baryonic matter (e.g Planck collaboration 2014) that behaves like a perfect fluid on large scales (e.g. Peacock 1999). In theory, deviations from the perfect-fluid model are expected to arise on scales comparable to the free-streaming 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 DM-particle 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 cut-off is placed on planet-size 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 self-gravitating “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 gamma-ray 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. Strongly-lensed galaxies provide complementary constraints on the clumpiness of dark matter haloes in the inner-most 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 high-mass end of the subhalo mass function. For example, the sizes of stream gaps produced by encounters with DM substructure in Milky Way-like 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 non-negligible 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. Collision-less -body cosmological simulations that explore the properties of CDM substructures within Milky Way-size systems, e.g. the Aquarius Project (Springel et al. 2008), Via Lactea (Diemand et al. 2007) and GHALO (Stadel et al. 2009), have particle-mass 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 gamma-ray annihilation signals (e.g. Koushiappas 2009), and/or the dynamics of self-gravitating 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, fly-by 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 phase-space 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 low-mass 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 Way-like galaxies for DM haloes made of cold/warm particles and ultra-light 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

 d2Rdt2=−∇Φg(R)+N∑i=1fi(R), (1)

where is the mean-field gravitational potential of the galaxy, and

 F≡N∑i=1fi, (2)

is the specific force induced by a set of self-gravitating substructures orbiting in the host galaxy. In systems where we expect to fluctuate stochastically along the phase-space trajectory of the particle, .

Solving (1) for an arbitrarily-large 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

 p(F)=1V∫d3r1×...×1V∫d3rNδ(F−∑ifi), (3)

where is the Dirac’s delta function. Fourier transforming yields

 ~p(k) = ∫d3Feik⋅Fp(F) = 1V∫d3r1×...×1V∫d3rN∫d3Feik⋅Fδ(F−∑jfj) = 1V∫d3r1×...×1V∫d3rNeik∑jfj = [1V∫d3reik⋅f]N.

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 re-written as

 1V∫d3reik⋅f=1V∫Vd3r[1−(1−eik⋅f)]=1−1V∫Vd3r(1−eik⋅f),

which elevated to the -th power becomes

where is the number density of substructures. It is useful to define the function

 ϕ(k)≡n∫Vd3r(1−eik⋅f), (5)

such that the inverse Fourier transform of (2) becomes

 p(F)=1(2π)3∫d3kexp[−ik⋅F−ϕ(k)]. (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)

 ϕin(k)=∫d3r(1−eik⋅f)n(r), (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 Point-mass particles

As a first step it is useful to review the analysis of fluctuations in the gravitational field induced by a random distribution of equal-mass particles (see Kandrup 1980 for a formal description of this method).

The gravitational force exerted by a point-mass particle can be written as

 f=−GMr2^r, (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

 ϕ(k) = 2πn∫∞0drr2∫+1−1d(cosθ)(1−exp[−ikGMcosθ/r2]) = 4πn∫∞0drr2(1−sin[kGM/r2]kGM/r2) = 415(2πGM)3/2nk3/2≡ak3/2,

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 well-known Holtsmark (1919) distribution

 p(F) = 1(2π)2∫∞0dkk2exp(−ak3/2)∫+1−1dxexp[−ikFx] = 12π2∫∞0dkk2exp(−ak3/2)sin(kF)kF.

The asymptotic behaviour of this function is studied in Chandrasekhar (1943) and reproduced here for completeness. In the weak-force limit we can approximate , such that

 limF→0p(F)=12π2∫∞0dkk2exp(−ak3/2)=13π2a2. (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 weak-force limit becomes flat. Deriving the strong-force limit requires a number of non-trivial steps. Chandrasekhar (1943) finds

 limF→∞p(F)≈158a(2π)3/2F−9/2=12(GM)3/2nF−9/2. (12)

It is important to stress that the large-force 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,

 p0(F)=nr2drF2dF=12(GM)3/2nF−9/2, (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 equal-mass 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 point-masses 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 point-mass particles in a galaxy.

Fig. 1 shows the combined force induced by equal-mass particles with homogeneously distributed over a sphere of radius . Comparison between the upper- and lower-left panels shows clear differences between the force measured from random statistical samples and that generated by an equilibrium ensemble at different time-steps. The reason behind the mismatch can be traced to the different speed of fluctuations, , in experiments A and B. In dynamical-equilibrium systems such as experiment B this quantity is modulated by time-scale 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) (long-dashed 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 equal-mass 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 dotted-lines 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

 ⟨F⟩ = ∫∞0d3Fp(F)F = 2π∫∞0dzsin(z)z∫∞0dFexp(−az3/2F−3/2) = 2πΓ(1/3)a2/3∫∞0dzsin(z)z2 = 2πΓ(1/3)Γ(3)a2/3 ≃ 8.879GMn2/3,

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 large-force, power-law tail of the distribution. Using Equation (12) it is straightforward to show that the variance of diverges in the strong-force limit as

 ⟨F2⟩=∫∞0d3Fp0(F)F2∝∫∞0dFF1/2→∞. (16)

In practice this means that the maximum force experienced by the test particle can grow up to arbitrarily-large values as the number realizations and time-steps increase in experiment A and B, respectively.

### 2.2 Extended substructures

Self-gravitating substructures with an extended matter distribution generate forces that (i) in general do not diverge at arbitrarily-close 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

 f=−GM(r+c)2^r. (17)

The density profile associated to (17) can be readily found by solving Poisson’s equation, which yields

 ρ(r)=M2πc31(r/c)(1+r/c)3.

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

 ϕ(k) = = 4πn∫∞0drr2(1−sin[kGM/(r+c)2]kGM/(r+c)2) = ≡ A(k)k3/2,

where is a dimension-less quantity. The last equality follows from the change of variable , with the function being defined as

 A(k) = a(I1+I2+I3)    with (19) I1 = 12√2π{ψ3[−5+2cos(1ψ2)+4√2πFC(√2π1ψ) +ψ[−4+3b4sin(1ψ2)]} I2 = +ψ[−3+2cos(1ψ2)+ψ2sin(1ψ2)]} I3 =

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

 p(F)=12π2F∫∞0dkksin(kF)exp[−A(k)k3/2]. (20)

Fig. 2 shows numerical solutions to the integral (20) for different values of the scale-length . 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

 p0(F)=12(GM)3/2nF−9/2(1−√F/f0)2    for    F

which recovers Equation (14) in the limit . The function is plotted in Fig. 2 with thin-dotted 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 closest-particle approximation at large distances, where .

The behaviour of in the strong-force limit is less clear-cut. Comparison of Equations (20) and (21) at shows that the distribution of forces generated by the nearest substructure either over/under-estimates the collective distribution (20) depending on whether the scale-length is smaller/larger than the typical separation between substructures, . Hence, there are two relevant regimes that define the large-force 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 self-gravitating, 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 scale-length 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 point-masses, here cannot be expressed analytically owing to the non-trivial functional form of in (19). In the upper panel of Fig. 3 we plot numerically-computed values of as a function of the substructure scale-length , 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 point-masses. 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 arbitrarily-large values insofar as the force generated by individual clumps does not diverge at arbitrarily-close distances. To illustrate this result let us compute in the rarefied regime, where nearby particles dominate the large-force tail of Equation (20), such that . From Equation (21)

 ⟨F2⟩ = ∫∞0d3Fp0(F)F2 = 2π(GM)3/2n∫f00dFF1/2(1−√F/f0)2 = 4π3(GM)3/2n√f0 = 4π3(GM)2nc,

which shows that the truncation of the large-force spectrum implies a finite variance. As expected, the divergence found for point-masses 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 point-mass particles, nearby substructures dominate the large-force, power-law 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

 n→∫dndMdM.

Power-law mass functions are of particular interest

 dndM=B0(MM0)α, (23)

where is our mass unit, a normalization parameter and the power-law 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,

 B0=n0Mα0×⎧⎪⎨⎪⎩1+αM1+α2−M1+α1,α≠−11ln(M2/M1),α=−1, (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 Point-mass 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

 ϕ(k) = 415(2πG)3/2B0k3/2∫M2M1dMM3/2(MM0)α = k3/2415(2πG)3/2B0Mα0×⎧⎪⎨⎪⎩Mα+5/22−M5/2+α1α+5/2,α≠−5/2ln(M2/M1),α=−5/2, ≡ k3/2a′.

Thus, varying the parameters of the mass function simply modifies the numerical value of .

Inserting and (24) into (2.1) yields

 ⟨F⟩ = 2πΓ(1/3)Γ(3)a′2/3 ≃ 8.879Gn2/30(α+1α+5/2Mα+5/22−Mα+5/21Mα+12−Mα+11)2/3,

where the indices and have been excluded for clarity. Hence, for we find that the average force scales as

• for (shallow mass profile).

• for .

• for (steep mass profile).

Similarly, replacing in Equation (12) and inserting the resulting distribution into (16) yields a divergent independently of the value of . As discussed in Section 2.1, the divergence of the second moment of can be traced back to the singular behaviour of the Keplerian force (8) at .

#### 2.3.2 Extended substructures

Self-gravitating 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 scale-radius as a function of . The power-law mass function seen in Section 2.3.1 can be readily modified in order to account for the mass-size relation by writing

 n→∫∫d2ndMdcdMdc

and

 d2ndMdc=B0(MM0)αδ[c−c(M)]. (27)

where denotes Dirac’s delta function. Equation (2.2) then becomes

 ϕ(k) = ×∫∞ψM(k)dz(z−ψM)2[1−z2sin(1z2)] ≡ A′(k)k3/2,

where , and is defined as

 A′(k)=415(2πG)3/2B0Mα0∫M2M1dMMα+3/23∑j=1Ij(k,M). (29)

with the functions given by Equation (19).

The relation between mass & scale-radius of extended substructures, , plays a key role in defining the large-force 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 power-law relation

 c(M)=c0(MM0)β. (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

 ⟨F2⟩ = 4π3∫M2M1dMdndMG2M2c(M) = 4π3Mβ−α0G2B0c0×⎧⎨⎩M3+α−β2−M3+α−β13+α−β, 3+α−β≠0ln(M2/M1), 3+α−β=0,

Note that for ensembles with a broad mass range, , Equation (2.3.2) exhibits two well-defined regimes

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

• Microstructure-dominated force distribution. In this case the least-massive substructures govern the large-force 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 α=−2

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 power-law 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 arbitrarily-large values as the minimum mass of substructures decreases.

The average force exerted by an ensemble of point-mass particles with is completely dominated by the most massive objects. This can be straightforwardly shown from Equation (2.3.1)

 ⟨F⟩ ≃ 8.879Gn2/30(2M1/22−M1/21M−11−M−12)2/3 ≃ 14.095GB2/30M4/30M1/32,

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

 ⟨F2⟩=4π3M2+β0G2B0c0M1−β2−M1−β11−β. (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 power-law 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 lower-mass limit, . Indeed, this is a macrostructure-dominated distribution, where the mean and the variance of are largely determined by the most massive substructures in the ensemble. In contrast, for power-law indices and 1.5 (middle and right panels respectively) the large-force distributions are microstructure-dominated, exhibiting a strong dependence on the lower mass limit, . Notice in particular that for steep size functions () the large-force tail of the force distribution approaches the point-mass behaviour in the limit , which according to Equation (16) leads to a divergent variance.

Interestingly, CDM subhaloes found in collision-less -body simulations of structure formation exhibit power-law 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 mass-size relation at redshift . Note, however, that cosmological simulations of Milky Way-size haloes typically resolve substructures with masses , with the low-mass 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ánchez-Conde & Prada 2014; Moliné et al. 2017 and references therein) suggests a negligible contribution of the smallest dark matter mini-haloes to the force distribution . We will return to this issue in §4.

## 3 Stochastic tidal forces

The equations of motion that govern the phase-trajectory of a tracer particle within a self-gravitating system that moves through a clumpy medium can be written as

 d2R′dt2=−∇Φs(R′)+Tg⋅R′+N∑i=1ti⋅R′, (34)

where and are 33 tidal tensors evaluated at the centre of the self-gravitating potential . The smooth component has a form

 Tjkg≡−∂2Φg∂xj∂xk, (35)

while the stochastic tidal tensor

 Tjk≡N∑i=1ti=N∑i=1∂fki∂xj=∂∂xjN∑i=1fki=∂Fk∂xj, (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

 ft≡t⋅R′. (37)

To simplify the analysis it is useful to diagonalize the tensor by rotating the coordinates to a new frame , which co-rotates with the angular velocity of the substructure, . In the non-inertial 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

 ft≈Rλ^u≡R\boldmathλ, (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 non-inertial 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 time-scale 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

 \boldmathΛ≡N∑i=1\boldmathλi. (39)

which arises from the combined tidal field of the external substructures. From (38) and (39) the combined tidal force can be written as

 Ft=N∑i=1\boldmathλiR=% \boldmathΛR. (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 point-mass particles and Hernquist (1990) spheres.

### 3.1 Point-mass particles

Consider a vector

 \boldmathλ=2GMr3^u, (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).

Inserting Equation (41) into (5) and integrating over , with , yields

 ϕ(k) = = = 2π23GMnk≡qk,

with .

The distribution can be expressed in an analytical form. Inserting (3.1) into (6) and integrating over yields

 p(\boldmathΛ) = 1(2π)2∫∞0dkk2exp(−qk)∫+1−1dxexp[−ikΛx] = 12π2∫∞0dkk2exp(−qk)sin(kΛ)kΛ = 1π2q31(1+ξ2)2,

where is a dimension-less quantity.

The distribution (3.1) exhibits two well-defined asymptotic behaviours depending on the magnitude of . In the weak-tides limit () Equation (3.1) becomes flat

 limΛ→0p(\boldmathΛ)=278π8(GMn)−3, (44)

whereas in the strong-tides limit () it follows a power-law tail

 limΛ→∞p(\boldmathΛ)=23GMnΛ4. (45)

As in Section 2.1, the large-force 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

 p0(\boldmathΛ)=nr2drΛ2dΛ=23GMnΛ4, (46)

which matches the strong-tide 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

 \boldmathλ=2GM(r+c)3^u. (47)

As expected, Equation (47) approaches asymptotically the point-mass 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

 ϕ(k) = = 2πn∫∞0drr2(2−sin[k2GM/(r+c)3]kGM/(r+c)3) = ≡ Q(k)k,

where is a dimension-less quantity, and the function is defined as

 Q(k)