An analytical description of substructure-induced gravitational perturbations in stellar systems

Perturbations to stellar systems can reﬂect the gravitational inﬂuence of dark matter substructures. On scales much smaller than the size of a stellar system, we point out analytical connections between the stellar and dark matter distributions. In particular, the density and velocity power spectra of the stars are proportional to the density power spectrum of the perturbing dark matter, scaled by k − 4 . This relationship allows easy evaluation of the suitability of a stellar system for detecting dark substructure. As examples, we show that the Galactic stellar halo is expected to be sensitive to cold dark matter substructure at wavenumbers k (cid:2) 0.3 kpc − 1 , and the Galactic disc might be sensitive to substructure at wavenumbers k ∼ 4 kpc − 1 . The perturbations considered in this work are short-lived, being rapidly erased by the stellar velocity dispersion, so it may be possible to attribute a detection to dark matter substructure without ambiguity.


INTRODUCTION
Dark matter has no known nongravitational interaction with Standard Model particles.Consequently, an important experimental strategy is to infer the properties of the dark matter from its gravitational influence alone.If the dark matter is cold, it should cluster on scales much smaller than galaxies, forming collapsed and tightly bound haloes that contain no visible matter at all.Such dark haloes could be detected gravitationally, but unambiguous detection remains elusive.
We use analytic methods to derive how the density and velocity fields of a gravitationally perturbed stellar system, and the statistics ★ E-mail: mdelos@carnegiescience.edu thereof, are related to the perturbing matter distribution.We consider arbitrary systems but focus on scales much smaller than the size of the system.On these length scales, the velocity dispersion of the stars erases perturbations on time scales much shorter than an orbital period.Although their short lifetimes limit the amplitudes of these perturbations, their transience can also be an advantage.Perturbations have been detected in stellar streams, but due to their long lifetimes, it is not clear that these perturbations can be attributed to dark matter substructure (e.g.Capuzzo Dolcetta et al. 2005;Küpper et al. 2008Küpper et al. , 2010Küpper et al. , 2012;;Amorisco et al. 2016;Pearson et al. 2017;Ibata et al. 2020;Qian et al. 2022;Weatherford et al. 2023).Sources of short-lived inhomogeneity are easier to disambiguate.
We find simple analytic relationships between the gravitationally induced perturbations of a stellar system and the properties of the perturbing structures.The power spectra of the stellar density and velocity fields become algebraically related to the power spectrum of the perturber density field.By comparing these induced power spectra to the intrinsic Poisson noise associated with the number of stars, one can easily test the sensitivity of a stellar system as a probe of dark substructure.
This article is organized as follows.In section 2, we explore how perturbations to a stellar distribution are related to the density field of the perturbing matter.In section 3, we specialize to two-point statistics, deriving relationships between the stellar density and velocity power spectra and the nonlinear matter power spectrum.In section 4, we discuss the regimes in which these results are valid, and we explore some applications.We conclude in section 5.

PERTURBATIONS TO THE STELLAR DISTRIBUTION
We begin by exploring how a stellar distribution is perturbed by the gravity of a separate distribution of matter.Let  (x, v, ) be 2 M. S. Delos the distribution function of the stars, which depends on position x, velocity v, and time .The collisionless Boltzmann equation reads

𝜕 𝑓 𝜕𝑡
where a(x, ) is the gravitational acceleration at the position x and time .Now separate  (x, v, ) =  0 (v) +  1 (x, v, ) into the unperturbed distribution function  0 , which we approximate to be spatially uniform and constant in time, and the perturbation  1 .At linear order in the perturbations, we obtain Note that we approximate the acceleration a(x, ) to be entirely perturbative.
To solve this, we Fourier transform from position x to wavenumber k, defining the solution to which is for an arbitrary integration constant  0 , which can be interpreted as the time that perturbations began.We will take  0 → −∞, i.e., we assume that perturbations began in the arbitrarily distant past.We now approximate that the unperturbed stellar velocity distribution is Maxwellian with dispersion  * per dimension, so that Here,  * ≡ ∫ d 3 v  0 (v) is the unperturbed stellar number density.We also specialize to the gravitational acceleration sourced by a mass density field (k, ).Then the perturbation to the distribution function, given by equation ( 4), becomes

Density perturbations
The fractional stellar density perturbation is The interpretation of this equation is simple if we imagine how the stars would respond to an external density field that is only present at one moment  = 0, so that (k,  It can be convenient to eliminate the integral over the matter density field  by temporal Fourier transformation from time  to frequency .In wavenumber-frequency space, the density perturbation in terms of the similarly transformed matter density field (k, ).Substituting Δ ≡  −  ′ separates the time integrals into a product of . With some manipulation, we obtain where the frequency dependence is contained in the function For numerical evaluation, it is useful to note that where  () ≡ e −  2 ∫  0 d e  2 is the Dawson integral.We plot () in figure 1.Note that () ranges from 1 for || ≪ 1 to −1/ 2 for || ≫ 1.Thus, in the low-frequency limit, Recall that ( * ) −1 is the time scale over which the velocity dispersion erases stellar perturbations on the scale  −1 .The interpretation of equation ( 13) is that if the driving density field  changes only on much longer time scales || −1 ≫ ( * ) −1 , then the stellar response  * aligns with the density field in space and time but is suppressed by the factor ( * ) −2 .Conversely, in the high-frequency limit, This means that if the driving density field  changes on much shorter time scales || −1 ≪ ( * ) −1 , then the stellar response  * is opposite from  in phase and is suppressed by the factor || −2 .

Velocity perturbations
Through similar calculations, or by using the continuity equation, the mean velocity of the stars is in wavenumber space and in wavenumber-frequency space.
If the stellar distribution  * (k, ) were known, then the mass distribution (k, ) that is gravitationally driving it could be inferred using equation ( 10), and there would be no need to consider v.However, knowledge of  * (k, ) means knowing the stellar distribution at all times.In practice, we can measure stellar positions and velocities today only.If we define the present day to be  = 0, then at the present time.That is, by measuring the present-day density and velocity fields of the stellar distribution, it is possible in principle to observationally infer the zeroth and first frequency moments of the quantity [/( * )] (k, ).

STELLAR AND DARK MATTER POWER SPECTRA
Further exploration of the possibility of probing the full matter density field (k, ) lies beyond the scope of this article.Instead, we will focus on the simpler possibility of probing its power spectrum.Suppose that the gravitating matter distribution  is statistically homogeneous (in space) and constant (in time).We can define a "spatiotemporal" matter density power spectrum P  (k, ) such that ) Here  * is the complex conjugate of , the angle brackets denote the average, and   and  3  are the one-and three-dimensional Dirac delta functions, respectively.Predictions for P  (k, ) can be estimated directly from numerical simulations or semianalytic substructure models, which is a subject of ongoing work.Here we approximate that the stellar perturbations arise entirely from dark matter, i.e. that  is the dark matter density field, and we discuss how P  (k, ) connects to conventional descriptions of dark matter structure.
As we will see, it is useful to define a function   (k, ) such that Here   () is the spatial matter density power spectrum defined such that where (k) is the spatial Fourier transform of the density field at some fixed time, which we can take to be  = 0.   () is constant in time by assumption.Also, we assume statistical isotropy at fixed times, so   () does not depend on the direction of k.
Note that   () = ρ2 (), where () is the power spectrum of the density contrast with respect to the cosmological average ρ, a commonly discussed quantity in cosmology.However, we emphasize that   () is not the matter power spectrum extrapolated using linearorder cosmological perturbation theory (as is standard in cosmology) but rather the nonlinear matter power spectrum, which describes the nonlinearly evolved matter density field.

Dark matter substructure
Nonlinear dark matter structure is often discussed in terms of subhalo models, and   () can be evaluated from such models in the following way.Suppose that haloes of mass  have differential number density d/d per mass interval and spatial volume and that a halo of mass  is spherical with density profile ( |) as a function of radius .Then if halo positions are uncorrelated, where ( |) is the threedimensional Fourier transform of ( |), i.e., If haloes have a nontrivial spatial distribution, then equation ( 22) also includes a two-halo term (e.g.Scherrer & Bertschinger 1991).However, subhalo positions are expected to be mostly uncorrelated due to phase mixing and the subdominance of their mutual gravitation.
Returning to the spatiotemporal matter power spectrum P  (k, ), the function   = P  /  (equation 20) can be estimated as follows.Suppose the dark matter field were moving coherently at a velocity u.Then as a function of time, (k, ) = (k)e −ik•u .Fourier transforming in time leads to and so For the simple scenario of coherent motion, evidently Notice how the spatiotemporal power spectrum P  =     can be anisotropic even if the constant-time spatial power spectrum   is isotropic.
A more realistic case is where the dark matter substructure has a Maxwellian velocity distribution with dispersion  per dimension.If density fields moving in different directions are uncorrelated, it suffices to average equation ( 26) over this velocity distribution, yielding (see also the more careful discussion in Delos & Schmidt 2022).This description is appropriate for studying perturbations in a stellar distribution that is at rest with respect to the dark matter distribution, such as the stars in the spheroidal component of a galaxy.A further sophistication is to suppose that the stellar distribution moves at velocity u 0 with respect to the centre of the Maxwellian dark matter velocity distribution.Now if u are the velocities in the rest frame of the stellar distribution, then u ′ ≡ u+u 0 are Maxwellian distributed.In this case, This description is appropriate when studying perturbations to stars in a galactic disk or in a satellite galaxy, for example, since these systems exhibit net motion with respect to the dark matter substructure.

Stellar density power
Consider now the spatial power spectrum  * (k) of the stellar density contrast, defined similarly to equation ( 21) but potentially anisotropic.From equations ( 17), ( 19), and (20),  * (k) is related to the dark matter density power spectrum   () as with () again given by equation ( 12).In general, the integral must be evaluated numerically.Note that |()| is numerically well behaved (see figure 1), as is   in realistic cases (equations 27 and 28).However, there are some cases where the relationship between  * and   can be written analytically.For this purpose, it is useful to note that equation (11) implies

Maxwellian dark matter substructure
Suppose that the dark matter substructure has a Maxwellian velocity distribution in the rest frame of the stellar distribution with dispersion  per dimension, so   is given by equation ( 27).Substituting equation (30) into equation ( 29) and carrying out the  integral first yields where we define the integral We note that if the dark matter and the stars have the same velocity dispersion, so  * = , then this quantity evaluates to  * (1) = (9 + 2 √ 3)/27 ≃ 0.74.If on the other hand the stars are much colder than the dark matter, so  * ≪ , then it evaluates to  * (0) = /2 3/2 ≃ 1.11.We plot  * () as the black curve in figure 2.

Shifted Maxwellian dark matter substructure
Now consider dark matter substructure with the same Maxwellian velocity distribution, but suppose that the stars move at net velocity u 0 10 −2 10 −1 10 0 velocity dispersion ratio σ * /σ Integral factor between the stellar density power spectrum  * and the matter density power spectrum   .This is the bracketed factor in equation ( 33), and it depends on the ratio  * / between the stellar and matter velocity dispersions (horizontal axis) and on k • u 0 / (different colours), the ratio between the projection of the mean stellar velocity onto the wavevector and the matter velocity dispersion.The k • u 0 = 0 case (black curve) corresponds to  * () given by equation ( 32).The dotted lines show the approximation of equation ( 34), which is valid in the limit of a cold stellar distribution.
with respect to the dark matter, so that   is given by equation ( 28).
Then equation ( 29) becomes where k ≡ k/ is the unit vector along k.We are not able to evaluate this expression analytically, but we show the numerically evaluated integral in brackets as the coloured solid curves in figure 2. Note that by inspection, it can only depend on the magnitude of k • u 0 , since |()| 2 is an even function of .Nevertheless, we can analytically treat the following limiting scenarios.If the wavevector k is perpendicular to the mean stellar velocity u 0 , then equation (33) simply reduces to equation (31).Conversely, suppose that | k • u 0 | ≫  * , corresponding to a cold stellar distribution with a large mean velocity that is not close to perpendicular with the wavevector k.Then ( + k • u 0 / * ) 2 ≃ ( k • u 0 / * ) 2 , since |()| 2 only has significant support for || ∼ O (1).In this limit where we have used that . As we noted above,  * (0) = /2 3/2 ≃ 1.11.Evidently, the stellar perturbations are exponentially suppressed by a nonzero mean velocity u 0 .This approximation is shown in figure 2 with the dotted curves.

Stellar velocity power
Let us define the fixed-time stellar velocity power spectrum      (k) such that where  and  are the indices of the mean velocity vectors.From equations (18-20), we can obtain (36) Integral factor between the stellar velocity power spectrum     and the matter density power spectrum   .This is the bracketed factor in equation ( 39).As with figure 32, we show how it depends on the stellar and matter velocity dispersions  * and , respectively, and the component k • u 0 of the mean stellar velocity along the wavevector.The k • u 0 = 0 case (black curve) corresponds to   () given by equation ( 38).The dotted lines show the approximation of equation ( 40), which is valid in the limit of a cold stellar distribution.

Maxwellian dark matter substructure
Similarly to section 3.2.1,we can consider a Maxwellian dark matter velocity distribution in the stellar rest frame and obtain in this case where we define the function For the cases of hot ( * = ) and cold ( * ≪ ) stellar distributions,   (1) = 2/3 5/2 ≃ 0.40 and   (0) = 3/2 5/2 ≃ 1.67, respectively.We plot   () as the black curve in figure 3.

Shifted Maxwellian dark matter substructure
Similarly to section 3.2.2,we also consider the case where the dark matter substructure has a Maxwellian velocity distribution but the stars have a net velocity u 0 with respect to the dark matter.In this case where k ≡ k/ is again the unit vector along k.We plot the bracketed integral as the coloured solid curves in figure 3, and note that it again depends only on the magnitude of k • u 0 .In the limit of a cold stellar distribution,  * ≪ | k • u 0 |, and similarly to equation ( 34) we obtain This approximation corresponds to the dotted curves in figure 3.

Density-velocity cross-correlations
Finally, we consider the cross power between the stellar density contrast and the mean velocity of stars.We may define a cross power Integral factor between the stellar density-velocity cross power spectrum  *   and the matter density power spectrum   .This is the negative of the bracketed factor in equation ( 43).It is only large when the stellar velocity dispersion  * , the matter velocity dispersion , and the component k • u 0 of the mean stellar velocity along the wavevector are all comparable.Note also that it is negative when k and by equations (17-20) we obtain In particular, if the dark matter has a Maxwellian velocity distribution and the stars have no net motion, so   is given by equation ( 27), then there is no cross power between v and .
If the stars have net velocity u 0 with respect to the dark matter, then We are not able to integrate this expression analytically, nor is the  * ≪ k • u 0 limit helpful, since it simply leads to  *   (k) = 0, but we show the numerically evaluated integral in figure 4. Note that it is negative for the k • u 0 > 0 that we show, but it is odd in k • u 0 , so it would be positive for k • u 0 < 0. However, figure 4 indicates that this integral is only large in magnitude when  * ∼  ∼ k • u 0 , which is a physically implausible scenario.If the stars are as hot as the dark matter, then we do not also expect a comparable magnitude of coherent (e.g.rotational) motion.Thus, we generally expect  *   (k) ≃ 0.

DISCUSSION
We now discuss the regimes in which the results of this analysis are applicable, and we present an application.

Applicability
The principal approximations made in these calculations are that the perturbing matter fields are the only source of gravitational acceleration and that the unperturbed systems are spatially uniform.Stars reside inside galaxies, of course, and we have neglected both the host galaxy's gravity and its overall density structure.These approximations are valid under the following conditions.
(i) The time scale ( * ) −1 for the stellar velocity dispersion to erase perturbations is much shorter than the stars' orbital time scale.If perturbations can persist over an orbital period, then one cannot neglect orbital dynamics.
(ii) The length scale  −1 under consideration is much smaller than the overall size of the system.Otherwise, the global structure of the system must be accounted for.
The first condition means that the stars must comprise a structure that is supported by its velocity dispersion on the scales that we are studying.The second condition does not impose an additional requirement, since the velocity dispersion sets the system's size.
Thus, the calculations in this work are applicable to spheroidal galaxies, spheroidal components of galaxies, and globular star clusters, on scales much smaller than the size of the system.These calculations also apply to galactic disks or stellar streams, but only if  −1 is much smaller than the system's smallest transverse size (although for larger scales, one can modify the description to include an approximate account of orbital dynamics; see Delos & Schmidt 2022).

Sensitivity estimates
A straightforward application of the power spectra of stellar perturbations discussed in section 3 is to estimate whether gravitational perturbations would be detectable in a given stellar system.The idea is to compare the predicted power spectrum to the intrinsic Poisson noise associated with the finite number of stars.We focus on the stellar density power spectrum  * discussed in section 3.2.If the stars have mean number density n * , then the power spectrum of the Poisson noise is  * = n−1 * .These estimates will be optimistic, because clustering in the stars can arise also due to the history of a stellar system.That is, stars could be clustered because they formed together or accreted together.Although we do not consider this effect further, we remark that historical clustering can be distinguished in principle from clustering associated with gravitational perturbations.When clustering arises from history, density and velocity perturbations must be strongly correlated in space; because otherwise they would rapidly dissipate.On the other hand, we showed that the transient density and velocity perturbations induced by substructure are spatially uncorrelated at a fixed time.

The stellar halo at 20 kpc
As a first example, we note that the Milky Way's stellar halo at a radius of about 20 kpc has a number density of order n * ∼ 10 −5 pc −3 , based on a power-law density profile with index -3.5 and local normalization of about 10 −4 M ⊙ pc −3 (Helmi 2008) and assuming the average star weighs 0.5 M ⊙ .The power spectrum of the Poisson noise is hence  * () = n−1 * ∼ 10 −4 kpc 3 .The radial velocity dispersion is about 120 km s −1 (Helmi 2008), and if we approximate that it is isotropic and the same for stars and substructures, then equation (31) implies that perturber power spectra larger than   () ∼ 10 13 (/kpc −1 ) 4 M 2 ⊙ kpc −3 produce stellar density perturbations that exceed the level of the Poisson noise.The dashed line in figure 5 shows this sensitivity limit.
For comparison, we also show the nonlinear power spectrum that  et al. (2008) of a Milky Way-like dark matter halo.The subhalo mass function is taken to be d/d ∝  −2 normalized so that the number density of 10 6 -10 7 M ⊙ subhaloes is 6 × 10 −4 kpc −3 .Subhaloes have Hernquist (1990) density profiles with scale radii  = 1.04 kpc (/10 8 M ⊙ ) 1/2 .This description is highly approximate, and a more accurate characterization of the local nonlinear matter power spectrum is a subject of ongoing work.Apparently, for  ≲ 0.3 kpc −1 , this substructure model would induce perturbations in the distribution of halo stars at ∼ 20 kpc that exceed the level of the Poisson noise.We show the power spectra with different maximum subhalo masses imposed in order to test which subhaloes are driving this result.The perturbations from sub-10 8 M ⊙ subhaloes would be detectable (orange curve), but those from sub-10 7 M ⊙ subhaloes (green curve) likely are not.

Local disk stars
The dotted lines in figure 5 show the estimated sensitivity of nearby stars in the Galactic disk to perturbations by substructure, plotted only up to  −1 = 0.3 kpc, which is approximately the vertical scale height of the dominant component of the disk (Bland-Hawthorn & Gerhard 2016).We assume the stars comprise a mass density of about 0.04 M ⊙ pc −3 and have a velocity dispersion of 30 km s −1 per dimension (Bland-Hawthorn & Gerhard 2016).We again assume that stars have average mass 0.5 M ⊙ and that the substructure velocity dispersion is 120 km s −1 per dimension.We adopt a rotation velocity of  0 = 250 km s −1 .The sensitivity depends on the angle between the wavevector k and the rotation velocity, which is why we show two dotted curves.Sensitivity is maximized (  is lowest) when k is perpendicular to the rotation velocity and minimized (  is highest) when it is parallel.
Figure 5 suggests that dark substructure around  ∼ 4 kpc −1 could be barely detectable through induced density perturbations of the local stellar field.However, it should be noted that substructure that crosses the disk may be disrupted to a significantly greater degree than is reflected in the dark-matter-only simulation on which the predicted   () in figure 5 is based (e.g.D 'Onghia et al. 2010).Another concern is that the calculations in this work neglect the selfgravity of the stellar perturbations themselves, which can be relevant for galactic disks (e.g.Widrow 2023), although this effect would only amplify the stellar perturbations (and so improve the sensitivity).

CONCLUSION
We derived analytic relationships between the distribution of a gravitationally perturbed stellar system and that of the matter that sources the perturbations.These results are expressed in Fourier space and are valid in the limit that the length scale  −1 under consideration is much smaller than the size of the stellar system.We considered the full stellar density and velocity fields in section 2 and their two-point statistics in section 3.In particular, the stellar density and velocity power spectra are related in a simple way to the nonlinear matter power spectrum.In both cases, the power spectrum of the stars is proportional to that of the perturbers, albeit weighted by  −4 .The cross power between stellar density and velocity is essentially zero.
Although the description is general, the motivation for this work was to explore the degree to which stellar systems can be used as a probe of dark matter substructures.As an example, we found that the stellar halo at ∼ 20 kpc is expected to be sensitive to cold dark matter substructure at  ≲ 0.3 kpc −1 .Perturbations are mostly driven by subhaloes larger than 10 8 M ⊙ , however, which may already contain visible galaxies.We also found that local disk stars may be sensitive to substructure at  ∼ 4 kpc −1 , with perturbations driven by subhaloes as light as ∼ 10 7 M ⊙ .However, this conclusion depends on the extent to which subhaloes survive encounters with the Galactic disk.
These tests of the sensitivity of the stellar disk and halo relied on highly approximate descriptions of dark matter substructure.More accurate predictions of the nonlinear matter power spectrum inside a Milky Way-like halo are a subject of ongoing work.We have demonstrated in particular that it is important to know the "spatiotemporal" two-point statistics of the dark matter field, rather than only the spatial two-point statistics at a fixed time.
′ ) ∝   ( ′ ), where   is the Dirac delta function.The resulting density perturbation is spatially in phase with  but evolves in time as  * ∝  e − (  *  ) 2 /2 .Here  ≡ /(  * ) and  * is the stellar velocity dispersion per dimension.The upper panel shows the magnitude of  (  ), while the lower panel shows its complex argument.For low frequencies |  | ≪   * ,  (  ) ≃ 1, so the stellar and matter fields are in phase.For high frequencies |  | ≫   * ,  (  ) ≃ −1/ 2 , so the stellar and matter fields are maximally out of phase.
Diemand et al. (2008))f the Milky Way's stellar halo and disk to probe dark matter substructure.We plot the nonlinear matter power spectrum   ( ); the dashed and dotted curves show the sensitivity limits for the Milky Way's stellar halo at ∼ 20 kpc and local disk stars, respectively.Specifically, for each wavenumber , the sensitivity limit corresponds to the smallest value of   ( ) that would induce fluctuations in the stellar density field that exceed the level of the intrinsic Poisson noise arising from the finite number of stars.The sensitivity for disk stars depends on the angle (shaded band) and is strongest (lowest   ) for wavevectors k perpendicular to the disk rotation velocity.For comparison, the solid curves show predicted cold dark matter (CDM) density power spectra, with different colours corresponding to different upper limits on the subhalo mass .We use the same subhalo model asBovy et al. (2017)andDelos & Schmidt (2022), which is based on the CDM simulation ofDiemand et al. (2008).
Delos & Schmidt (2022)ons 22 and 23)to the subhalo model used byBovy et al. (2017)andDelos & Schmidt (2022)to study perturbations to stellar streams at a comparable Galactocentric radius.This subhalo model is based on the numerical simulation of Diemand