Stochastic tidal heating by random interactions with extended substructures
Gravitating systems surrounded by a dynamic sea of substructures experience fluctuations of the local tidal field which inject kinetic energy into the internal motions. This paper uses stochastic calculus techniques to describe ‘tidal heating’ as a random walk of orbital velocities that leads to diffusion in a 4-dimensional energy–angular momentum space. In spherical, static potentials we derive analytical solutions for the Green’s propagators directly from the number density and velocity distribution of substructures with known mass & size functions without arbitrary cuts in forces or impact parameters. Furthermore, a Monte-Carlo method is presented, which samples velocity ’kicks’ from a probability function and can be used to model orbital scattering in fully generic potentials. For illustration, we follow the evolution of planetary orbits in a clumpy environment. We show that stochastic heating of (mass-less) discs in a Keplerian potential leads to the formation, and subsequent “evaporation” of Oort-like clouds, and derive analytical expressions for the escape rate and the fraction of comets on retrograde orbits as a function of time. Extrapolation of the subhalo mass function of Milky Way-like haloes down to the WIMP free-streaming length suggests that objects in the outer Solar system experience repeated interactions with dark microhaloes on dynamical time-scales.
keywords:Cosmology: dark matter; kinematics and dynamics; Oort cloud; methods: statistical.
Understanding the evolution of gravitating systems embedded in a clumpy medium is both theoretically and numerically challenging. First, because gravitational forces cannot be shielded, which thwarts any clear-cut definition of a physical ‘boundary’ between the system and the surrounding background. Second, because dynamical equilibrium can never be attained in regions where the dynamical time is comparable to the time-scale on which the external field fluctuates. The analytical hurdles involved in the description of non-equilibrium systems subject to long-range forces (see Padmanabhan 1990; Lynden-Bell 1999 for reviews) often calls for the heuristic assumption of “isolation”, whereby the contribution of distant objects to the local acceleration is ignored. However, this approximation breaks down on long time-scales, as the cumulative effect of repeated interactions with background objects dominates over the secular evolution of the system.
The first attempt to construct a statistical theory for the response of a single, free-moving, mass-less particle to a fluctuating external force is due to Chandrasekhar. He proposed three different approaches to tackle this problem. In his first classical paper, Chandrasekhar (1941a) divides the force acting on a tracer particle into two components: a force that changes very slowly and can be expressed as the gradient of a smooth potential, plus a random contribution of “chance stellar encounters” of short duration . The computation of the random part relies on two assumptions (i) an infinite homogeneous medium, and (ii) independent encounters, such that the average effect of the clumpy background will be the sum of the effects of separate two-body encounters. Alas, Chandrasekhar finds that these simplifications lead to a total energy variation , where is the energy exchanged during a single encounter, that diverges when the integral over the minimum two-body separation (the so-called “impact parameter”) is taken to distances comparable to the intra-particle separation, , as well as when it tends to infinity. Both issues are typically resolved by truncating the range of impact parameters at small and large distances and , respectively, which results in an energy variation that is proportional to an ill-defined Coulomb logarithm , whose computation is a long-standing matter of debate (e.g. Just & Peñarrubia 2005 and references therein).
Due to these shortcomings, Chandrasekhar (1941b) proposes to abandon the two-body approximation altogether, and argues in favour of stochastic methods that describe the combined force exerted by the background onto a test star, , which depends on the instantaneous relative position of particles and is therefore subject to fluctuations. In this framework, the average acceleration experienced by a test star during a time interval leads to a net variation of the velocity , and a variance , where brackets denote averages over the probability function to experience a force in the interval , and is the mean-life of a fluctuation. Chandrasekhar (1941b) derives using a method originally devised by Holtsmark (1919) to study the motion of charged particles in a plasma, and suggests that may correspond to Smoluchowski’s (1916) time-scale associated with a stochastic system, yet also warning that “Smoluchowski’s ideas cannot be applied without further deep generalizations of them”. To this aim, Chandrasekhar & von Neumann (1942; 1943) define the mean life of fluctuations as , where brackets denote an average over the bivariate distribution , which determines the simultaneous probability to experience a force F and an associated rate of change . Kandrup (1980) shows that Smoluchowski’s and Chandrasekhar & von Neumann’s derivations of yield consistent results modulo numerical constants of order unity. We return to this issue in §3.2.
A third approach was proposed by Chandrasekhar (1944), in which the response of tracer particles to fluctuating forces is treated as a Brownian motion in velocity space, where the first and second moments of the velocity increments correspond to the drift and diffusion coefficients, respectively. In an isotropic medium the drift coefficient vanishes by symmetry, , while the diffusion coefficient is computed as . Here brackets denote averages over the autocorrelation function , which determines the probability to feeling a force at an initial time , and a force at a later time . Chandrasekhar (1944) finds that when background objects are assumed to move on straight lines during a short time interval () the autocorrelation function is such that , which leads to a diffusion coefficient that diverges logarithmically. Lee (1968) shows that the divergence arises from the assumption of an infinite medium, and argues in favour of truncating the range of impact parameters at large and small distances, which re-introduces the Coulomb logarithm found by Chandrasekhar (1941a) in the theory.
An appealing aspect of the stochastic analysis of Chandrasekhar (1941b) is the cancelation of the net force contribution of particles separated by large distances, which eliminates the need of an arbitrary cut-off at weak forces. Yet, this approach still requires a truncation at short separations, or strong forces, due to the divergent force generated by point-mass particles at (Chandrasekhar 1941b; Cohen et al. 1950; Kandrup 1980). Recently, Peñarrubia (2018) shows that this last hindrance can be removed from the theory by considering a background of extended substructures with individual forces f that (i) are not centrally divergent, and (ii) approach a Keplerian limit at distances , where the is the substructure size
This paper extends the work of Chandrasekhar to tracer particles moving in a smooth potential subject to random fluctuations of an external tidal field, . The random component is generated by the combined gravitational field of a large number of extended substructures orbiting within a host potential in dynamical equilibrium. Section 2 summarizes the results of Peñarrubia (2018), who computes the spectrum of tidal fluctuations generated by an homogeneous distribution of extended objects using Holtsmark (1919) statistical technique. Section 3 applies autocorrelation (Chandrasekhar 1944) and stochastic (Chandrasekhar 1941b) methods to derive the drift and diffusion coefficients ( and , respectively) from the number density and velocity distribution of substructures with fixed mass and size without arbitrary cuts in forces or impact parameters. Both derivations are shown to provide consistent results modulus a numerical factor of order unity that arises from the different assumptions on which the two methods rest. However, while the derivation of the autocorrelation function assumes that background objects move on straight-line trajectories, the stochastic approach requires no a priori information on the motion of substructures in the host potential. The method proposed by Chandrasekhar (1941a), which treats individual encounters separately, is not explored here, as it requires a non-trivial treatment of the 3-body problem (e.g. Heggie & Rasio 1996). In Chandrasekhar’s theories, the computation of diffusion coefficients relies on the impulsive approximation, which assumes that the location of a tracer particle does not vary appreciably during an encounter. In Section 3.3 we apply adiabatic corrections computed by Weinberg (1994a,b,c), which allow for the motion of tracer particles during the duration of a tidal fluctuation.
Section 4 uses the probability theory presented by Peñarrubia (2015) to describe the non-equilibrium state of gravitating systems subject to random tidal interactions as a diffusion process in the integral-of-motion space. Analytical expressions for Green’s functions are given for spherical, static potentials, where the integrals correspond to the energy (), and the three components of the angular momentum (L). The derivation of the Green’s functions is described in detail in Appendix B, and accounts for the fact that gravitationally-bound particles can only diffuse in a confined region of the integral-of-motion space. To this end, we set an absorbing boundary at , such that particles with escape from a gravitating system, and treat the angular momentum volume as a cubic box with reflecting surfaces placed at the angular momentum of circular orbits with a fixed energy, .
Section 5 tests our analytical framework by running -body experiments where (i) the tidal tensor is computed directly from the relative positions of a population of extended substructures in dynamical equilibrium within a host potential, and (ii) -body particles experience velocity “kicks” which are sampled from a probability function using Monte-Carlo techniques.
As an application of the theoretical methods, Section 6 follows the evolution of idealized mass-less discs at in a Keplerian potential. The time-dependent distribution of tracer particles in the 4-dimensional integral-of-motion space, , is calculated as a convolution of the Green’s functions with the initial distribution . Tidal evaporation rates, and the rate of production of retrograde orbits are derived analytically from the flux of particles crossing the boundaries and , respectively. The results are compared against Monte-Carlo -body models.
Following up on the results of Peñarrubia (2018), Section 7 discusses the effects of dark matter microhaloes on weakly-bound objects in the outskirts of the Solar system. Our analysis relies on bold extrapolations of the subhalo properties found in Milky Way-like haloes down to free-streaming mass-scales, and must be therefore taken with caution. At face value, the results suggest that objects in the Oort cloud may experience a very large number of interactions with microhaloes during a single orbital period, which opens up the interesting possibility to use Trans-Neptunian Objects (TNO’s) to probe the (local) subhalo mass function on sub-solar mass scales.
2 Stochastic fluctuations of the tidal field
This Section briefly summarizes the analytical framework of Peñarrubia (2018), hereafter Paper I, who introduce a statistical technique for deriving the spectrum of random fluctuations of the tidal field generated by a large population of extended substructures.
2.1 Probability theory
As a starting point, consider a tracer particle at a distance R from the centre of a small, system with a self-gravitating potential located at a galactocentric radius from a larger host galaxy. For simplicity, we shall work in the collision-less, mean-field limit, where the granularity in the system may be ignored. Hence, in the reference frame of the host galaxy the gravitational acceleration experienced by a tracer particle can be written as
where is the mean-field gravitational potential of host galaxy, and
is the specific force induced by a set of extended substructures in a state of dynamical equilibrium within the host galaxy.
Similarly, the equations of motion that define the trajectory of the smaller system about the parent galaxy can be written as
where and are 33 tidal tensors evaluated at the centre of the self-gravitating potential . The smooth component has a form
while the stochastic tidal tensor
arises from the gradient in the combined tidal force induced by a set of substructures distributed across the host galaxy.
Paper I shows that the stochastic term of the tidal force can be approximately written as
where is the sum of eigenvalues of the effective tidal tensor , and is the centrifugal force component in a frame that co-rotates with the angular velocity of individual substructures, (see also Renaud et al. 2011). Note that the vectors have random directions if substructures are isotropically distributed around the test particle. As discussed in Paper I, Equation (7) neglects the Euler and Coriolis terms appearing in the non-inertial rest frame, which is a reasonable approximation in an impulsive regime, where one can assume that the angular frequency remains constant during the encounter duration. We will return to this issue in Section 3.
In analogy with Equation (2), it is useful to define the tidal vector
such that the combined tidal force (7) becomes . A large population of extended substructures homogeneously distributed within a volume around a test particle generates a stochastic tidal field that is fully specified by the probability density , which defines the probability of experiencing a tidal vector in the interval . Following up the method of Holtsmark (1919), which was originally devised to study the motion of charged particles in a plasma, Paper I derives the spectrum of tidal fluctuations as
where is the Dirac’s delta function and
with denoting the number density of substructures. The above derivation can be generalized to inhomogeneous substructure distributions as (Chandrasekhar 1941b; Kandrup 1980; Chavanis 2009)
where is the number density profile, and r is centred at the location of the test particle. In the local approximation the number density can be assumed to be roughly constant, , thus recovering (10). Paper I shows that this approximation holds in regions where the number density profile varies on scales that are much larger than the averaged separation between subhaloes, i.e. , where the distance can be measured from the probability of finding the closest substructure within the volume
which peaks at .
In general, the probability density can be rarely expressed in an analytical form. An interesting exception with broad applications in cosmology corresponds to a large population of Hernquist (1990) spheres with a density profile
where and are the substructure mass and scale length, respectively. Individual substructures induce a specific tidal force
which approaches the Keplerian (‘particle’) limit as . The sum of eigenvalues associated with the force (13) can be straightforwardly calculated using Renaud et al. (2011) formalism
where is a unit vector pointing in a random direction. The maximum value of Equation (14) is at .
Given that an ensemble of self-gravitating, overlapping objects is not dynamically stable, in what follows we limit our analysis to populations of Hernquist (1990) spheres with sizes much smaller than their typical separation (). In this regime, which Paper I characterizes as ’rarefied’ in analogy with the kinetic theory of gases, the combined tidal vector fluctuates stochastically with a probability density (9) that can be expressed analytically as
where is a dimension-less quantity, , and is a normalization constant that guarantees . The distribution (15) obeys the isotropic condition , which implies that the tidal vectors point in random directions. Note that in the weak-force limit ( the probability function approaches asymptotically a constant value, which implies that the effect of having an increasing number of particles at large distances exactly balances with the declining force, such that becomes flat at small accelerations. The strong-force limit ( is dominated by the contribution of the nearest object (see Paper I), thus the probability density is truncated at the maximum force derivative exerted by an individual substructure (14), such that for .
As we shall see below, the second moment of is particularly relevant for understanding the dynamical response of a tracer particle subject to a stochastic tidal field. After some algebra, Paper I finds that
which diverges in the particle (point-masss) limit . In practice, the divergence of means that as the time progresses the maximum tidal force experienced by a tracer particle can grow up to arbitrarily-large values.
3 Tidal heating
A test particle surrounded by a large population of moving substructures experiences random fluctuations of the local tidal force due to the rapid change of the (relative) position of nearby objects. Over a sufficiently long interval of time the cumulative effect of multiple encounters contributes to the randomization of the peculiar velocity. If the velocity impulses are small, , the fluctuations can be thought to occur independently, and the average effect of a series of fluctuations will be the sum of the expectation values of the effect of a single fluctuation. With these simplifications in place, the formalism of Brownian motion can be used to describe the response of self-gravitating objects to stochastic variations of a gravitational field (e.g. Chandrasekhar 1943), which reduces the problem of tidal heating to the computation of diffusion coefficients.
The averaged velocity increments acquired by tracer particles over short intervals of time, and – here brackets denote averages over multiple interactions– play a key role in the Brownian motion theory. Below we explore two independent methods to calculate these quantities. The first method (§3.1) follows the usual derivation of velocity increments, where one assumes that (i) the speed of fluctuations is much faster than the orbital velocity of the tracer particle in the potential , and (ii) external substructures move on linear trajectories within the host potential . The second approach is presented in §3.2 and follows up on the arguments of Chandrasekhar (1941b) and Kandrup (1980) to derive the speed of fluctuations directly from the phase-space distribution of substructures, without making assumptions on the orbital motion of these objects. Section 3.3 extends this formalism to encounters in a non-impulsive regime by applying Weinberg (1994a,b,c) adiabatic corrections to the coefficients obtained in §3.2.
3.1 Impulse/Straight-line approximation
The classical derivation of and relies on two key assumptions: (i) substructures move on straight lines, and (ii) the typical duration of tidal fluctuations is much shorter than the orbital period of the tracer particle about the potential . Under these conditions it is relatively straightforward to extend the analysis of Lee (1968) to velocity increments arising from stochastic fluctuations of an external tidal field.
Consider first a test particle orbiting in a potential with an orbital frequency, , where and are the moduli of the position and velocity vectors, respectively. During a short time interval, , the particle position does not change appreciably. In contrast, from Equation (7) the net contribution of the stochastic tidal forces leads to a velocity variation
where is the tidal vector acting on the tracer particle at the time . Equation (17) is typically known as the impulse approximation. For a distribution of substructures uniformly distributed around the tracer particle the distribution of tidal fluctuations is isotropic, . Hence, by symmetry, the average velocity increment is
The computation of the squared velocity increment is considerably more involved. From Equation (17)
where and are tidal vectors acting on a tracer particle at the times and , respectively. Recall that our working assumptions are that (i) test particles suffer a large number of individual encounters with neighbour substructures during a time interval , and (ii) velocity increments induced by subsequent encounters are statistically uncorrelated. Under those conditions, the product of tidal vectors only depends on the length of the time interval , such that . If the function decreases more rapidly than for large , as it happens in most cases of astrophysical interest, then one can take the limit . Hence, the averaged squared velocity increment (19) can be written as a single integral (see Appendix A of Lee 1968 for a detailed derivation)
To compute we must take into account that the tidal vectors and are not statistically independent. Indeed, these quantities are related through the trajectories of individual substructures in the host potential. Following Chandrasekhar (1944), let us define the autocorrelation function , which gives the probability that a tracer particle experiences a tidal vector at , and at a later time . In the straight-line approximation, a substructure with an initial position vector r and a relative velocity v will be located at at . Hence, for a population of substructures with a distribution function , where , the autocorrelation function can be written as (see Appendix A)
The explicit evaluation of is difficult, but the second moment is readily found from its Fourier transform (see Appendix A) as
with by symmetry. For most cases of astrophysical interest, decreases more rapidly than at large radii, and one can let for simplicity.
For a random population of Hernquist (1990) spheres, the averaged squared velocity increment (20) can be expressed analytically under a proper choice of coordinates. Let us adopt a coordinate frame in which the tangential velocity of the substructure is parallel to the position vector r at . During a short time interval, , the norm of the orbital plane, , remains approximately constant, whereas the substructure moves to a new location . Using the framework of §2.1, it straightforward to show that in this configuration the direction of the centrifugal acceleration remains invariant , and one can set in (23). Combination of (23), (20) and (14) then yields
which also diverges in the particle limit . The quantity is the reciprocal of a characteristic velocity whose meaning is discussed in §3.2.
3.2 Stochastic approach
The autocorrelation method described in §3.1 rests upon the assumption substructures move on straight-line trajectories. This condition is not required in a stochastic approach, where fluctuations of the combined tidal field arise from the presence at some point of time of substructures in the vicinity of the test particle.
Suppose that at there are particles distributed within a volume around the test particle. The number of neighbours inside will change either because one of the substructures exists this volume, or because another substructure enters from outside . Smoluchowski (1916) shows that the probability that at some later time there are still substructures inside may be written as . Here, the time interval corresponds to the mean life of a state in which the number of substructures within remains constant. For a uniform distribution of substructures with a number density and a mean squared (relative) velocity , Smoluchowski (1916) finds (see also Chandrasekhar 1941b)
Equation (25) has simple asymptotic behaviours. At small distances, , it recovers the time-scale derived from the straight-line approximation, , which simply corresponds to the time that a substructure moving at a constant speed would take to cross the radius . At large radii, however, the probability that a substructure leaves/enters becomes proportional to the number of substructures within this volume. Accordingly, in the limit Smoluchowski’s timescale becomes inversely proportional to , such that , which vanishes in the limit .
Following Chandrasekhar (1941b) and Paper I, let us assume that the tidal field acting on a tracer particle is entirely dominated by the nearest substructure. Hence, from Equation (14) one can identify . The duration of the tidal fluctuation (25) then becomes
where the time-scale approximately corresponds to the average time that it takes a substructure to travel twice the distance .
In the straight-line limit, , Equation (26) has a simple form
which diverges as in the limit .
For point-mass particles, , Equation (26) reduces to
Thus, in the strong-force regime the duration of the fluctuations scales as a power-law , which matches the straight-line asymptotic behaviour of Equation (27). In contrast, for weak forces, , Equation (28) scales as , thus vanishing in the limit , which strongly deviates from the divergent duration one would expect for substructures moving on straight-line trajectories. This result was first noticed by Chandrasekhar & von Neumann (1942), and is discussed in detail by Kandrup (1980), who warns that adopting the straight-line approximation leads to a serious overestimation of the contribution of distant encounters to the velocity increments derived in §3.1.
To illustrate these results, Fig. 1 shows the duration of fluctuations induced by substructures with a mass and an average separation , moving with a relative velocity dispersion . As expected, comparison between the duration expected in the straight-line approximation (dotted lines) and Smoluchowski’s time-scale shows good agreement at large forces (), but strongly disagree in the weak-force regime , where the straight-line assumption leads to a divergent . For extended substructures, the duration of tidal fluctuations vanishes in the limits and . In the particle-limit , it approaches a power-law behaviour at .
The maximum of defines a characteristic time span, , associated with the longest, and therefore most likely, tidal fluctuation, . After some algebra, one can show that the solution to , with given by (26), is
which does not depend on substructure mass or size.
Following Chandrasekhar (1941b), let us model impulsive velocity increments induced by a fluctuating tidal field as random-walk process in velocity space. For a substructure population distributed homogeneously around the tracer particle, the distribution of velocity impulses is isotropic and has a Gaussian form (Chandrasekhar 1943; Kandrup 1980)
where denotes the probability that a test particle with a velocity V will experience a velocity impulse within a time interval . From Equation (18), the average velocity increment vanishes by symmetry, , whereas the variance given by Equation (20) can be re-written as111See Equation (62) of Chandrasekhar (1941b).
Note that Equation (32) obeys the ergodic property, as the contribution of a tidal fluctuation with a magnitude to the averaged velocity impulse is weighted by the duration of the fluctuation, .
In order to find an analytical solution to Equation (32) it is convenient to approximate . As shown in Fig. 1, this approximation is accurate at large forces, , which dominate the behaviour of the variance (see Paper I). In addition, it is useful to introduce the dimensionless quantity , which approaches asymptotically as . Combination of (15) and (27) then yields
In the limit one has , hence at leading order one can write (33) as
Thus, the variance of the velocity increments (32)‘ becomes
Expression (35) is remarkably similar to the result derived from the autocorrelation function, Equation (24), which suggests that the characteristic speed in Smoluchowski’s equation, , corresponds to the inverse of the averaged reciprocal velocity , modulo a numerical factor of order unity that arises from the different assumptions on which the two methods rest222To be precise, equating (24) to (35) yields ..
To gain further physical insight onto the physical meaning of these quantities, let us for example consider a population of extended substructures with an isotropic Maxwellian velocity distribution displaced by a velocity V
where V is the velocity of a tracer particle measured from the centre of the self-gravitating potential . The average of the reciprocal velocity is
here is the error function, which scales as for , and converges towards at . On the other hand, the averaged squared velocity is
It is straightforward to show that and have practically the same asymptotic limits. Indeed, if the tracer particle moves slowly with respect to the substructure population, , the inverse of the averaged reciprocal velocity becomes , whereas Smoluchowski’s characteristic velocity scales as . In contrast, if the tracer particle moves with a large velocity, , both theories asymptotically approach .
3.3 Adiabatic corrections
The results obtained Sections 3.1 and 3.2 assume that tidal interactions occur in an impulsive regime, wherein the location of a tracer particle does not vary appreciably during the duration of a flyby encounter/tidal fluctuation. This approximation is generally valid in the outskirts of self-gravitating systems, where the orbital period can be much longer than the time-scale on which the tidal field varies. However, it is bound to fail in the inner regions of the self-gravitating potential , where particles complete their orbits on short time-scales and thus react adiabatically to external perturbations. In this Section we introduce adiabatic corrections that allow for the motion of stars during tidal field fluctuations. In particular, we use the corrections found by Weinberg (1994a,b,c) (see also Gnedin & Ostriker 1999)
where corresponds to the ratio between the actual energy change and the value found under the impulse approximation, is the orbital frequency of the tracer particle in the potential , and is the duration of tidal fluctuations. As expected, the correction becomes increasingly insignificant () in the outskirts of the system, where orbital frequencies are short, . In contrast, the correction is strong in the inner regions of the potential, where orbital frequencies become much higher, , and Equation (39) approaches a power-law asymptotic form .
In order to account for the adiabatic response of tracer particles to relatively long tidal fluctuations, we ‘correct’ the impulsive variation of the kinetic energy (32) as
which must be solved numerically.
Comparison of (42) and (44) reveals key differences between the impulsive and adiabatic regimes. First, the impulsive velocity variance scales as , whereas in the adiabatic regime . The opposite behaviour means that while the outskirts of self-gravitating objects are predominantly heated by nearby substructures with small relative velocities, only fast encounters can perturb the internal regions of the potential. Furthermore, we find that , which implies that particles close to the centre of the potential moving with high peculiar velocities will be barely affected by fluctuations of the external tidal field.
Fig. 2 shows numerical solutions to Equation (41) for different substructure size-to-average separation ratios, . Blue-dotted and orange-dashed lines show the impulsive and adiabatic limits derived from Equations (34) and (43), respectively. Note first that accounting for the force dependence of is particularly important for ‘rarefied‘ distributions of substructures (). Indeed, combination of (34), (30) and (16) shows that adopting a constant value for the duration of tidal fluctuations, , overestimates the magnitude of tidal heating by a factor