Stochastic tidal heating by random interactions with extended substructures

TL;DR

This paper models how random gravitational interactions with substructures cause stochastic tidal heating, leading to orbital diffusion and potential ejection of objects, with analytical and Monte Carlo methods applied to planetary systems.

Contribution

It introduces a stochastic calculus framework for tidal heating, deriving analytical solutions and a Monte Carlo method for modeling orbital scattering in complex potentials.

Findings

Stochastic heating causes orbital diffusion and evaporation of planetary bodies.

Analytical expressions for escape rates and orbital distributions are derived.

Dark microhaloes may influence Solar system objects through repeated interactions.

Abstract

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.

Equations

$$\frac{{\mathrm{d}}^2(\text{R}+{\text{r}}_s)}{\mathrm{d}{t}^2}=-\nabla {\Phi }_s(\text{R})-\nabla {\Phi }_g(\text{R}+{\text{r}}_s)+\sum _{i=1}^N{\text{f}}_i(\text{R}+{\text{r}}_s),$$

$$\text{F}\equiv \sum _{i=1}^N{\text{f}}_i,$$

$$\frac{{\mathrm{d}}^2{\text{r}}_s}{\mathrm{d}{t}^2}=-\nabla {\Phi }_g({\text{r}}_s)+\sum _{i=1}^N{\text{f}}_i({\text{r}}_s).$$

$$\frac{{\mathrm{d}}^2\text{R}}{\mathrm{d}{t}^2}=-\nabla {\Phi }_s(\text{R})+T_g\cdot \text{R}+\sum _{i=1}^N{t}_i\cdot \text{R}+\mathcal{O}(R/r_s),$$

$$T_g^{jk}\equiv -\frac{{\partial }^2{\Phi }_g}{\partial x_j\partial x_k},$$

$$T^{jk}\equiv \sum _{i=1}^N{t}_i^{jk}=\sum _{i=1}^N\frac{\partial f_i^k}{\partial x^j}=\frac{\partial }{\partial x^j}\sum _{i=1}^N{f}_i^k=\frac{\partial F^k}{\partial x^j},$$

$$\displaystyle{\textbf{{F}}}_{t}\equiv\sum_{i=1}^{N}t_{i}\cdot{\textbf{{R}}}\approx R\sum_{i=1}^{N}{\mbox{\boldmath$\lambda$}}_{i},$$

$$\displaystyle{\mbox{\boldmath$\Lambda$}}\equiv\sum_{i=1}^{N}{\mbox{\boldmath$\lambda$}}_{i},$$

$$\displaystyle p({\mbox{\boldmath$\Lambda$}})$$

$$\displaystyle\approx\frac{1}{(2\pi)^{3}}\int{\rm d}^{3}k\exp\big{[}-i{\textbf{{k}}}\cdot{\mbox{\boldmath$\Lambda$}}-\phi({\textbf{{k}}})\big{]}~{}~{}~{}~{}{\rm for}~{}~{}~{}~{}N\gg 1,$$

$$\displaystyle\phi({\textbf{{k}}})\equiv n\int_{V^{\prime}}{\rm d}^{3}r\big{(}1-e^{i{\textbf{{k}}}\cdot{\mbox{\boldmath$\lambda$}}({\textbf{{r}}})}\big{)},$$

$$\displaystyle\phi_{\rm in}({\textbf{{k}}})=\int{\rm d}^{3}r\big{(}1-e^{i{\textbf{{k}}}\cdot\mbox{\boldmath$\lambda$}({\textbf{{r}}})}\big{)}n({\textbf{{r}}}),$$

$$\displaystyle p({\textbf{{r}}}){\rm d}^{3}r\sim\exp\big{(}-\frac{4}{3}\pi r^{3}n\big{)}4\pi r^{2}n{\rm d}r,$$

$$\rho (r)=\frac{M}{2\pi c^3}\frac{1}{(r/c)(1+r/c{)}^3},$$

$$\text{f}=-\frac{GM}{(r+c{)}^2}\hat{\text{r}},$$

$$\displaystyle\mbox{\boldmath$\lambda$}=\frac{2GM}{(r+c)^{3}}\hat{\textbf{{u}}},$$

$$\displaystyle p({\mbox{\boldmath$\Lambda$}})\simeq\frac{C}{\pi^{2}q^{3}}\frac{1}{(1+\xi^{2})^{2}}\big{[}1-\big{(}q/\lambda_{0}\big{)}^{1/3}\xi^{1/3}\big{]}^{2}~{}~{}~{}~{}~{}~{}{\rm for}~{}~{}\Lambda<\lambda_{0},$$

$$⟨{\Lambda }^2⟩$$

$$\displaystyle=\frac{4C}{\pi}q^{2}\int_{0}^{\lambda_{0}/q}{\rm d}\xi\frac{\xi^{4}}{(1+\xi^{2})^{2}}\big{[}1-\big{(}q/\lambda_{0}\big{)}^{1/3}\xi^{1/3}\big{]}^{2}$$

$$\simeq \frac{8\pi }{15}\frac{(GM{)}^{2}n}{c^3},$$

$$\displaystyle\Delta{\textbf{{V}}}=\int_{0}^{t}{\rm d}s\,{\textbf{{F}}}_{t}(s)\approx R\int_{0}^{t}{\rm d}s\,{\mbox{\boldmath$\Lambda$}}_{s}.$$

$$⟨\Delta \text{V}⟩$$

$$\displaystyle=R\int_{0}^{t}{\rm d}s\,\int{\rm d}^{3}\Lambda_{s}p(\Lambda_{s}){\mbox{\boldmath$\Lambda$}}_{s}=0.$$

$$\displaystyle|\Delta{\textbf{{V}}}|^{2}=R^{2}\int_{0}^{t}{\rm d}s\int_{0}^{t}{\rm d}s^{\prime}{\mbox{\boldmath$\Lambda$}}_{s}\cdot{\mbox{\boldmath$\Lambda$}}_{s^{\prime}},$$

$$\displaystyle\langle|\Delta{\textbf{{V}}}|^{2}\rangle\simeq R^{2}2t\int_{0}^{\infty}{\rm d}\tau\,\langle{\mbox{\boldmath$\Lambda$}}_{0}\cdot{\mbox{\boldmath$\Lambda$}}_{\tau}\rangle.$$

$$\displaystyle W({\mbox{\boldmath$\Lambda$}}_{0},{\mbox{\boldmath$\Lambda$}}_{\tau})=\int\frac{{\rm d}^{3}k_{0}}{(2\pi)^{3}}\int\frac{{\rm d}^{3}k_{\tau}}{(2\pi)^{3}}\exp\bigg{[}-i({\textbf{{k}}}_{0}\cdot\mbox{\boldmath$\Lambda$}_{0}+{\textbf{{k}}}_{\tau}\cdot\mbox{\boldmath$\Lambda$}_{\tau})-\phi({\textbf{{k}}}_{0},{\textbf{{k}}}_{\tau})\bigg{]},$$

$$\displaystyle\phi({\textbf{{k}}}_{0},{\textbf{{k}}}_{\tau})=n\int{\rm d}^{3}vf({\textbf{{v}}})\int_{V^{\prime}}{\rm d}^{3}r\bigg{[}1-\exp\bigg{(}i{\textbf{{k}}}_{0}\cdot{\mbox{\boldmath$\lambda$}}({\textbf{{r}}})+i{\textbf{{k}}}_{\tau}\cdot{\mbox{\boldmath$\lambda$}}({\textbf{{r}}}+{\textbf{{v}}}\tau)\bigg{)}\bigg{]}.$$

$$\displaystyle\langle{\mbox{\boldmath$\Lambda$}}_{0}\cdot{\mbox{\boldmath$\Lambda$}}_{\tau}\rangle$$

$$\displaystyle=-\bigg{(}\frac{\partial\phi}{\partial{\textbf{{k}}}_{0}}\cdot\frac{\partial\phi}{\partial{\textbf{{k}}}_{\tau}}+\frac{\partial^{2}\phi}{\partial{\textbf{{k}}}_{0}\cdot\partial{\textbf{{k}}}_{\tau}}\bigg{)}_{|{\textbf{{k}}}_{0}|\to 0,|{\textbf{{k}}}_{t}|\to 0}$$

$$\displaystyle=n\int{\rm d}^{3}vf({\textbf{{v}}})\int_{V^{\prime}}{\rm d}^{3}r\,\mbox{\boldmath$\lambda$}({\textbf{{r}}})\cdot\mbox{\boldmath$\lambda$}({\textbf{{r}}}+{\textbf{{v}}}\tau),$$

$$⟨|\Delta \text{V}{|}^2⟩$$

$$=t{R}^{2}8(GM{)}^{2}n\int {\mathrm{d}}^{3}vf(\text{v})\int \frac{{\mathrm{d}}^{3}r}{(r+c{)}^3}{\int }_0^{\infty }\frac{\mathrm{d}\tau }{(r+v\tau +c{)}^3}$$

$$=t{R}^{2}16\pi (GM{)}^{2}n\int {\mathrm{d}}^{3}v\frac{f(\text{v})}{v}{\int }_0^{\infty }\mathrm{d}r\frac{r^2}{(r+c{)}^5}$$

$$=t{R}^2\frac{4\pi }{3}\frac{(GM{)}^2}{c^2}n\int {\mathrm{d}}^{3}v\frac{f(\text{v})}{v},$$

$$T(r)=\sqrt{\frac{2\pi }{3⟨v^2⟩}}\frac{r}{\frac{4\pi }{3}{r}^{3}n+1}.$$

Full text

1 Introduction

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 $\sum\Delta E^{2}$, where $\Delta E$ 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, $D$, 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 $b_{\rm min}$ and $b_{\rm max}$, respectively, which results in an energy variation that is proportional to an ill-defined Coulomb logarithm $\ln(b_{\rm max}/b_{\rm min})$, 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, $\textbf{{F}}=\sum_{i=1}^{N}\textbf{{f}}_{i}$, which depends on the instantaneous relative position of $N\gg 1$ particles and is therefore subject to fluctuations. In this framework, the average acceleration experienced by a test star during a time interval $t>0$ leads to a net variation of the velocity $\langle\Delta\textbf{{v}}\rangle=\langle\textbf{{F}}\,T(F)\rangle$, and a variance $\langle|\Delta\textbf{{v}}|^{2}\rangle=t\,\langle F^{2}\,T(F)\rangle$, where brackets denote averages over the probability function $p(\textbf{{F}})$ to experience a force in the interval $\textbf{{F}},\textbf{{F}}+{\rm d}\textbf{{F}}$, and $T(F)$ is the mean-life of a fluctuation. Chandrasekhar (1941b) derives $p(\textbf{{F}})$ using a method originally devised by Holtsmark (1919) to study the motion of charged particles in a plasma, and suggests that $T(F)$ 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 $T(F)=|\textbf{{F}}|/\sqrt{\langle|{\rm d}\textbf{{F}}/{\rm d}t|^{2}\rangle_{F}}$, where brackets denote an average over the bivariate distribution $W(\textbf{{F}},{\rm d}\textbf{{F}}/{\rm d}t)$, which determines the simultaneous probability to experience a force F and an associated rate of change ${\rm d}\textbf{{F}}/{\rm d}t$. Kandrup (1980) shows that Smoluchowski’s and Chandrasekhar & von Neumann’s derivations of $T(F)$ 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, $\langle\Delta{\textbf{{v}}}\rangle=0$, while the diffusion coefficient is computed as $\langle|\Delta{\textbf{{v}}}|^{2}\rangle=2t\int_{0}^{\infty}{\rm d}\tau\,\langle{\textbf{{F}}}_{0}\cdot{\textbf{{F}}}_{\tau}\rangle$. Here brackets denote averages over the autocorrelation function $W(\textbf{{F}}_{0},\textbf{{F}}_{t})$, which determines the probability to feeling a force $\textbf{{F}}_{0}$ at an initial time $t_{0}=0$, and a force $\textbf{{F}}_{t}$ at a later time $t$. Chandrasekhar (1944) finds that when background objects are assumed to move on straight lines during a short time interval ($\tau$) the autocorrelation function is such that $\langle{\textbf{{F}}}_{0}\cdot{\textbf{{F}}}_{\tau}\rangle\sim\tau^{-1}$, 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 $r\ll D$ (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 $f\sim 1/r^{2}$ at distances $r\gg c$, where the $c$ is the substructure size

This paper extends the work of Chandrasekhar to tracer particles moving in a smooth potential $\Phi_{s}$ subject to random fluctuations of an external tidal field, ${\rm d}\textbf{{F}}/{\rm d}\textbf{{r}}$. The random component is generated by the combined gravitational field of a large number of extended substructures orbiting within a host potential $\Phi_{g}$ 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 ($\langle\Delta{\textbf{{v}}}\rangle$ and $\langle|\Delta{\textbf{{v}}}|^{2}\rangle$, 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 ($E$), 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 $E=0$, such that particles with $E\geq 0$ 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, $L_{c}(E)$.

Section 5 tests our analytical framework by running $N$-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) $N$-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 $t_{0}=0$ in a Keplerian potential. The time-dependent distribution of tracer particles in the 4-dimensional integral-of-motion space, $N(E,\textbf{{L}},t)$, is calculated as a convolution of the Green’s functions with the initial distribution $N(E_{0},\textbf{{L}}_{0},t_{0})$. Tidal evaporation rates, and the rate of production of retrograde orbits are derived analytically from the flux of particles crossing the boundaries $E=0$ and $L_{z}=0$, respectively. The results are compared against Monte-Carlo $N$-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 $\Phi_{s}$ located at a galactocentric radius ${\textbf{{r}}}_{s}$ 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

[TABLE]

where $\Phi_{g}$ is the mean-field gravitational potential of host galaxy, and

[TABLE]

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

[TABLE]

Taylor-expanding (1) at first order and subtracting (3) yields the well-known tidal approximation

[TABLE]

where $T_{g}$ and ${t}_{i}$ are 3$\times$3 tidal tensors evaluated at the centre of the self-gravitating potential $\Phi_{s}$. The smooth component has a form

[TABLE]

while the stochastic tidal tensor

[TABLE]

arises from the gradient in the combined tidal force induced by a set of $N-$substructures distributed across the host galaxy.

Paper I shows that the stochastic term of the tidal force can be approximately written as

[TABLE]

where $\lambda={\rm Trace}({t}_{e})$ is the sum of eigenvalues of the effective tidal tensor $t_{e}=t+\nabla f_{c}$, and {\textbf{{f}}}_{c}={\mbox{\boldmath\Omega}}\times({\mbox{\boldmath\Omega}}\times{\textbf{{R}}}) is the centrifugal force component in a frame that co-rotates with the angular velocity of individual substructures, $\Omega$ (see also Renaud et al. 2011). Note that the vectors {\mbox{\boldmath\lambda}}_{i} 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 $\Omega$ 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

[TABLE]

such that the combined tidal force (7) becomes {\textbf{{F}}}_{t}=\mbox{\boldmath\Lambda}R. A large population of extended substructures homogeneously distributed within a volume $V^{\prime}=4\pi d^{3}/3$ around a test particle generates a stochastic tidal field that is fully specified by the probability density p(\mbox{\boldmath\Lambda}), which defines the probability of experiencing a tidal vector in the interval {\mbox{\boldmath\Lambda}},{\mbox{\boldmath\Lambda}}+{\rm d}{\mbox{\boldmath\Lambda}}. 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

[TABLE]

where $\delta$ is the Dirac’s delta function and

[TABLE]

with $n\equiv N/V^{\prime}$ denoting the number density of substructures. The above derivation can be generalized to inhomogeneous substructure distributions as (Chandrasekhar 1941b; Kandrup 1980; Chavanis 2009)

[TABLE]

where $n({\textbf{{r}}})$ 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, $n({\textbf{{R}}}+{\textbf{{r}}}_{s})\approx n({\textbf{{r}}}_{s})=n$, 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. $|\nabla n/n|^{-1}\gg D$, where the distance $D$ can be measured from the probability of finding the closest substructure within the volume $V^{\prime}$

[TABLE]

which peaks at $D\equiv(2\pi n)^{-1/3}$.

In general, the probability density p({\mbox{\boldmath\Lambda}}) 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

[TABLE]

where $M$ and $c$ are the substructure mass and scale length, respectively. Individual substructures induce a specific tidal force

[TABLE]

which approaches the Keplerian (‘particle’) limit as $c\to 0$. The sum of eigenvalues associated with the force (13) can be straightforwardly calculated using Renaud et al. (2011) formalism

[TABLE]

where $\hat{\textbf{{u}}}$ is a unit vector pointing in a random direction. The maximum value of Equation (14) is $\lambda_{0}\equiv 2GM/c^{3}$ at $r=0$.

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 ($c\ll D$). In this regime, which Paper I characterizes as ’rarefied’ in analogy with the kinetic theory of gases, the combined tidal vector $\Lambda$ fluctuates stochastically with a probability density (9) that can be expressed analytically as

[TABLE]

where $\xi\equiv\Lambda/q$ is a dimension-less quantity, $q\equiv\frac{2\pi^{2}}{3}GMn=(\pi/3)(GM/D^{3})$, and $C$ is a normalization constant that guarantees \int{\rm d}^{3}\Lambda p(\mbox{\boldmath\Lambda})=1. The distribution (15) obeys the isotropic condition p(\mbox{\boldmath\Lambda})=p(\Lambda), which implies that the tidal vectors $\Lambda$ point in random directions. Note that in the weak-force limit ($\Lambda\ll q)$ 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 p(\mbox{\boldmath\Lambda}) becomes flat at small accelerations. The strong-force limit ($\Lambda\gg q)$ 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 p(\mbox{\boldmath\Lambda})=0 for $\Lambda>\lambda_{0}=2GM/c^{3}$.

As we shall see below, the second moment of p({\mbox{\boldmath\Lambda}}) 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

[TABLE]

which diverges in the particle (point-masss) limit $c\to 0$. In practice, the divergence of $\langle\Lambda^{2}\rangle$ 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, $|\Delta{\textbf{{V}}}|<<|{\textbf{{V}}}|$, 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, $\langle\Delta{\textbf{{V}}}\rangle$ and $\langle|\Delta{\textbf{{V}}}|^{2}\rangle$ – 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 $\Phi_{s}$, and (ii) external substructures move on linear trajectories within the host potential $\Phi_{g}$. 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 $\langle\Delta{\textbf{{V}}}\rangle$ and $\langle|\Delta{\textbf{{V}}}|^{2}\rangle$ 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 $\Phi_{s}$. 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 $\Phi_{s}(\textbf{{R}})$ with an orbital frequency, $w\equiv V/R$, where $R$ and $V$ are the moduli of the position and velocity vectors, respectively. During a short time interval, $t\ll w^{-1}$, 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

[TABLE]

where {\mbox{\boldmath\Lambda}}_{s}=\mbox{\boldmath\Lambda}(s) is the tidal vector acting on the tracer particle at the time $s$. 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, p(\mbox{\boldmath\Lambda})=p(\Lambda). Hence, by symmetry, the average velocity increment is

[TABLE]

The computation of the squared velocity increment is considerably more involved. From Equation (17)

[TABLE]

where {\mbox{\boldmath\Lambda}}_{s} and {\mbox{\boldmath\Lambda}}_{s^{\prime}} are tidal vectors acting on a tracer particle at the times $s$ and $s^{\prime}$, 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 $t$, 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 $\tau=s^{\prime}-s$, such that {\mbox{\boldmath\Lambda}}_{s}\cdot{\mbox{\boldmath\Lambda}}_{s+\tau}={\mbox{\boldmath\Lambda}}_{0}\cdot{\mbox{\boldmath\Lambda}}_{\tau}. If the function {\mbox{\boldmath\Lambda}}_{0}\cdot{\mbox{\boldmath\Lambda}}_{\tau} decreases more rapidly than $\tau^{-1}$ for large $\tau$, as it happens in most cases of astrophysical interest, then one can take the limit $t\to\infty$. Hence, the averaged squared velocity increment (19) can be written as a single integral (see Appendix A of Lee 1968 for a detailed derivation)

[TABLE]

To compute \langle{\mbox{\boldmath\Lambda}}_{0}\cdot{\mbox{\boldmath\Lambda}}_{\tau}\rangle we must take into account that the tidal vectors {\mbox{\boldmath\Lambda}}_{0} and {\mbox{\boldmath\Lambda}}_{\tau} 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 W({\mbox{\boldmath\Lambda}}_{0},{\mbox{\boldmath\Lambda}}_{\tau}), which gives the probability that a tracer particle experiences a tidal vector {\mbox{\boldmath\Lambda}}_{0} at $t=0$, and {\mbox{\boldmath\Lambda}}_{\tau} at a later time $t=\tau$. In the straight-line approximation, a substructure with an initial position vector r and a relative velocity v will be located at ${\textbf{{r}}}^{\prime}={\textbf{{r}}}+{\textbf{{v}}}\tau$ at $t=\tau$. Hence, for a population of substructures with a distribution function $n({\textbf{{r}}},{\textbf{{v}}})=nf({\textbf{{v}}})$, where $n\int_{r<d}{\rm d}^{3}r\int^{3}v\,f({\textbf{{v}}})=N$, the autocorrelation function can be written as (see Appendix A)

[TABLE]

where

[TABLE]

The explicit evaluation of W({\mbox{\boldmath\Lambda}}_{0},{\mbox{\boldmath\Lambda}}_{\tau}) is difficult, but the second moment is readily found from its Fourier transform $\tilde{W}({\textbf{{k}}}_{0},{\textbf{{k}}}_{\tau})$ (see Appendix A) as

[TABLE]

with $\partial\phi/\partial{\textbf{{k}}}_{0}=\partial\phi/\partial{\textbf{{k}}}_{\tau}=0$ by symmetry. For most cases of astrophysical interest, $\lambda(r)$ decreases more rapidly than $r^{-3/2}$ at large radii, and one can let $V^{\prime}\to\infty$ 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 $t=0$. During a short time interval, $\tau\ll w^{-1}$, the norm of the orbital plane, \hat{\mbox{\boldmath\Omega}}, remains approximately constant, whereas the substructure moves to a new location ${\textbf{{r}}}^{\prime}={\textbf{{r}}}+{\textbf{{v}}}\tau=(r+v\tau)\hat{\textbf{{r}}}$. Using the framework of §2.1, it straightforward to show that in this configuration the direction of the centrifugal acceleration remains invariant \hat{\textbf{{f}}}_{c}^{\prime}=\hat{\mbox{\boldmath\Omega}}^{\prime}\times[\hat{\mbox{\boldmath\Omega}}^{\prime}\times\hat{\textbf{{r}}}^{\prime}]=\hat{\mbox{\boldmath\Omega}}\times[\hat{\mbox{\boldmath\Omega}}\times\hat{\textbf{{r}}}], and one can set \hat{\mbox{\boldmath\lambda}}\cdot\hat{\mbox{\boldmath\lambda}}^{\prime}\approx 1 in (23). Combination of (23), (20) and (14) then yields

[TABLE]

which also diverges in the particle limit $c\to 0$. The quantity $\langle 1/v\rangle=\int{\rm d}^{3}vf({\textbf{{v}}})/v$ 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 $t=0$ there are $N$ particles distributed within a volume $V^{\prime}=4\pi r^{3}/3$ around the test particle. The number of neighbours inside $V^{\prime}$ will change either because one of the substructures exists this volume, or because another substructure enters from outside $V^{\prime}$. Smoluchowski (1916) shows that the probability $P_{N}(t)$ that at some later time there are still $N$ substructures inside $V^{\prime}$ may be written as $P_{N}(t){\rm d}t=e^{-t/T}{\rm d}t/T$. Here, the time interval $T$ corresponds to the mean life of a state in which the number of substructures within $V^{\prime}$ remains constant. For a uniform distribution of substructures with a number density $n=N/V^{\prime}$ and a mean squared (relative) velocity $\langle v^{2}\rangle$, Smoluchowski (1916) finds (see also Chandrasekhar 1941b)

[TABLE]

Equation (25) has simple asymptotic behaviours. At small distances, $4\pi r^{3}n/3=(2/3)(r/D)^{3}\ll 1$, it recovers the time-scale derived from the straight-line approximation, $T\sim r/\sqrt{\langle v^{2}\rangle}$, which simply corresponds to the time that a substructure moving at a constant speed $\sqrt{\langle v^{2}\rangle}$ would take to cross the radius $r$. At large radii, however, the probability that a substructure leaves/enters $V^{\prime}$ becomes proportional to the number of substructures within this volume. Accordingly, in the limit $4\pi r^{3}n/3\gg 1$ Smoluchowski’s timescale becomes inversely proportional to $N=nV^{\prime}$, such that $T\sim r/\sqrt{\langle v^{2}\rangle}N^{-1}\sim r^{-2}$, which vanishes in the limit $r\to\infty$.

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 $r+c=(2GM/\Lambda)^{1/3}$. The duration of the tidal fluctuation (25) then becomes

[TABLE]

where the time-scale T_{0}\equiv D\big{(}\frac{6}{\pi}\big{)}^{1/3}\sqrt{\frac{2\pi}{3\langle v^{2}\rangle}}\approx 1.8\frac{D}{\sqrt{\langle v^{2}\rangle}} approximately corresponds to the average time that it takes a substructure to travel twice the distance $D$.

In the straight-line limit, $4\pi nr^{3}/3\ll 1$, Equation (26) has a simple form

[TABLE]

which diverges as $T_{sl}\sim\Lambda^{-1/3}$ in the limit $\Lambda\to 0$.

For point-mass particles, $c\to 0\,(\lambda_{0}\to\infty)$, Equation (26) reduces to

[TABLE]

Thus, in the strong-force regime $\Lambda\gg q$ the duration of the fluctuations scales as a power-law $T_{c=0}\sim\Lambda^{-1/3}$, which matches the straight-line asymptotic behaviour of Equation (27). In contrast, for weak forces, $\Lambda\ll q$, Equation (28) scales as $T_{c=0}\sim\Lambda^{2/3}$, thus vanishing in the limit $\Lambda\to 0$, 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 $G=M=1$ and an average separation $D=1$, moving with a relative velocity dispersion $\langle v^{2}\rangle=1$. As expected, comparison between the duration expected in the straight-line approximation (dotted lines) and Smoluchowski’s time-scale $T(\Lambda)$ shows good agreement at large forces ($\Lambda\gg q$), but strongly disagree in the weak-force regime $\Lambda\ll q$, where the straight-line assumption leads to a divergent $T(\Lambda)$. For extended substructures, the duration of tidal fluctuations vanishes in the limits $\Lambda\to 0$ and $\Lambda\to\lambda_{0}$. In the particle-limit $\lambda_{0}\to\infty$, it approaches a power-law behaviour $T\sim\Lambda^{-1/3}$ at $\Lambda\gg q$.

The maximum of $T(\Lambda)$ defines a characteristic time span, $T_{ch}$, associated with the longest, and therefore most likely, tidal fluctuation, $\Lambda_{ch}$. After some algebra, one can show that the solution to $\frac{{\rm d}T}{{\rm d}\Lambda}(\Lambda_{ch})=0$, with $T(\Lambda)$ given by (26), is

[TABLE]

Note that the characteristic amplitude (29) shifts towards smaller values as the ratio $c/D$ increases. Inserting (29) into (26) returns the characteristic duration of tidal fluctuations

[TABLE]

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)

[TABLE]

where $\Psi({\textbf{{V}}},\Delta{\textbf{{V}}},t)$ denotes the probability that a test particle with a velocity V will experience a velocity impulse $\Delta{\textbf{{V}}}$ within a time interval $t$. From Equation (18), the average velocity increment vanishes by symmetry, $\langle\Delta{\textbf{{V}}}\rangle=0$, whereas the variance given by Equation (20) can be re-written as111See Equation (62) of Chandrasekhar (1941b).

[TABLE]

Note that Equation (32) obeys the ergodic property, as the contribution of a tidal fluctuation with a magnitude $\Lambda$ to the averaged velocity impulse is weighted by the duration of the fluctuation, $T(\Lambda)$.

In order to find an analytical solution to Equation (32) it is convenient to approximate $T(\Lambda)\simeq T_{sl}(\Lambda)$. As shown in Fig. 1, this approximation is accurate at large forces, $\Lambda\gg q$, which dominate the behaviour of the variance $\langle\Lambda^{2}\rangle$ (see Paper I). In addition, it is useful to introduce the dimensionless quantity $\chi=q/\lambda_{0}=(\pi/6)(c/D)^{3}$, which approaches asymptotically $\chi\to 0$ as $c\to 0$. Combination of (15) and (27) then yields

[TABLE]

In the limit $\chi\ll 1$ one has $C\to 1$, hence at leading order one can write (33) as

[TABLE]

Thus, the variance of the velocity increments (32)‘ becomes

[TABLE]

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, $\sqrt{\langle v^{2}\rangle}$, corresponds to the inverse of the averaged reciprocal velocity $\langle 1/v\rangle^{-1}$, 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 $\langle v^{2}\rangle^{-1/2}=5\langle 1/v\rangle/\sqrt{6\pi}\simeq 1.15\langle 1/v\rangle$..

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

[TABLE]

where V is the velocity of a tracer particle measured from the centre of the self-gravitating potential $\Phi_{s}$. The average of the reciprocal velocity is

[TABLE]

here ${\rm erf}(x)$ is the error function, which scales as ${\rm erf}(x)\simeq 2x/\sqrt{\pi}$ for $x\ll 1$, and converges towards ${\rm erf}(x)\simeq 1$ at $x\gg 1$. On the other hand, the averaged squared velocity is

[TABLE]

It is straightforward to show that $\langle 1/v\rangle^{-1}$ and $\sqrt{\langle v^{2}\rangle}$ have practically the same asymptotic limits. Indeed, if the tracer particle moves slowly with respect to the substructure population, $V\ll\sigma$, the inverse of the averaged reciprocal velocity becomes $\langle 1/v\rangle^{-1}\simeq\sqrt{\pi}\sigma$, whereas Smoluchowski’s characteristic velocity scales as $\sqrt{\langle v^{2}\rangle}\simeq\sqrt{3}\sigma$. In contrast, if the tracer particle moves with a large velocity, $V\gg\sigma$, both theories asymptotically approach $\sqrt{\langle v^{2}\rangle}=\langle 1/v\rangle^{-1}=V$.

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 $\Phi_{s}$, 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)

[TABLE]

where $A$ corresponds to the ratio between the actual energy change and the value found under the impulse approximation, $w=V/R$ is the orbital frequency of the tracer particle in the potential $\Phi_{s}$, and $T$ is the duration of tidal fluctuations. As expected, the correction becomes increasingly insignificant ($A\approx 1$) in the outskirts of the system, where orbital frequencies are short, $wT\ll 1$. In contrast, the correction is strong in the inner regions of the potential, where orbital frequencies become much higher, $wT\gg 1$, and Equation (39) approaches a power-law asymptotic form $A\sim(wT)^{-3}$.

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

[TABLE]

where $A(\Lambda)=A[T(\Lambda)]$. For $c\ll D$, combination of (15) and (26) yields

[TABLE]

which must be solved numerically.

The strong dependence of the adiabatic correction (39) with the orbital frequency of the tracer particle introduces two relevant regimes for the asymptotic behaviour of (41):

• •

Impulsive regime arises when the characteristic duration of tidal fluctuations is much shorter than the orbital period of the tracer particle in the self-gravitating potential $(T_{ch}\ll w^{-1})$. In this case $A\approx 1$, Equation (41) approaches asymptotically (34), and $\langle\Lambda^{2}TA\rangle_{\rm imp}\simeq\langle\Lambda^{2}T_{sl}\rangle$. Hence, inserting (33) into (40) yields

[TABLE]

• •

Adiabatic regime applies to tidal fluctuations with a characteristic time-span that is much longer than the orbital period of the tracer particle $(T_{ch}\gg w^{-1})$. In this limit, the adiabatic correction scales as $A\sim(wT)^{-3}\approx(wT_{sl})^{-3}$, and (41) becomes

[TABLE]

Inserting (43) into (40) and taking $\chi\ll 1$ yields

[TABLE]

Comparison of (42) and (44) reveals key differences between the impulsive and adiabatic regimes. First, the impulsive velocity variance scales as $\langle|\Delta{\textbf{{V}}}|^{2}\rangle_{\rm imp}\sim\langle v^{2}\rangle^{-1/2}$, whereas in the adiabatic regime $\langle|\Delta{\textbf{{V}}}|^{2}\rangle_{\rm ad}\sim\langle v^{2}\rangle$. 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 $\langle|\Delta{\textbf{{V}}}|^{2}\rangle_{\rm ad}\sim R^{5}/V^{3}$, 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, $c/D$. 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 $T(\Lambda)$ is particularly important for ‘rarefied‘ distributions of substructures ($D\gg c$). Indeed, combination of (34), (30) and (16) shows that adopting a constant value for the duration of tidal fluctuations, $T\approx T_{ch}$, overestimates the magnitude of tidal heating by a factor $\langle\Lambda^{2}\rangle T_{ch}/\langle\Lambda^{2}TA\rangle_{\rm imp}=(6^{1/3}2/9)(D/c)\simeq 0.4D/c$. Importantly, Fig. 2 shows that the impulsive range of frequencies extends towards systematically higher frequencies as the size ratio $c/D$ decreases. The accuracy of the straight-line approximation also improves in this limit, giving a progressively better match between the numerical values of (41) and its asymptotic limits $\langle\Lambda^{2}TA\rangle_{\rm imp}$ and $\langle\Lambda^{2}TA\rangle_{\rm ad}$ as $c/D\to 0$.

4 Evolution driven by tidal heating

In this Section we use the probability theory introduced by Peñarrubia (2015), hereafter P15, to describe the non-equilibrium state of particle ensembles acted on by a stochastic tidal field. For simplicity, we assume that the potential $\Phi_{s}$ is time-independent and has spherical symmetry, such that $\Phi_{s}({\textbf{{R}}},t)=\Phi_{s}(R)$. Extending the analysis to self-gravitating, time-dependent potentials is more involved and will be presented in a separate contribution.

4.1 Diffusion in a static potential

The presence of a stochastic tidal field in the equations of motion (4) induces random fluctuations of the integrals $\{E,{\textbf{{L}}}\}$. In an impulse regime (see §3), the location of a tracer particle is assumed to remain constant. Hence, the variation of orbital energy is equal to the change of kinetic energy

[TABLE]

while the angular momentum varies by an amount

[TABLE]

If one assumes that the population of substructures are isotropically distributed around the tracer particle, and that the velocity impulses are small, $|\Delta{\textbf{{V}}}|\ll V$, then it is straightforward to show that the average variation of these quantities can be written as

[TABLE]

whereas at leading order $\mathcal{O}(|\Delta{\textbf{{V}}}|/V)$ the variance of the integral fluctuations can be approximately written as

[TABLE]

Recall that in our notation brackets denote averages over the spectrum of tidal fluctuations, i.e. $\langle X\rangle=\int{\rm d}^{3}\Lambda\,p(\Lambda)\,X$, see §2.

4.1.1 Free-diffusion in an infinite domain

Following P15, let us define a statistical ensemble of tracer particles as the collection of a large number of individual mass-less particles with energy $E=E_{0}$ and angular momentum ${\textbf{{L}}}={\textbf{{L}}}_{0}$ at the time $t_{0}=0$. The probability that these particles have integrals in the range $E,E+{\rm d}E$ and ${\textbf{{L}}},{\textbf{{L}}}+{\rm d}{\textbf{{L}}}$ at a later time $t>t_{0}$ is given by the function $p(E,{\textbf{{L}}},t|E_{0},{\textbf{{L}}}_{0},t_{0})$, which a solution to a 4-dimensional diffusion equation (see P15 for details)

[TABLE]

with initial conditions $p=\delta(E-E_{0})\delta({\textbf{{L}}}-{\textbf{{L}}}_{0})$ at $t=t_{0}$. The simplest solution corresponds to particles that diffuse freely (i.e. in a domain with no boundaries) in the 4 dimensions of the integral-of-motion space, and with a probability density that is separable in energy and angular momentum, such that

[TABLE]

Under these conditions, P15 shows that the probability functions $p$ are Green’s functions (or propagators) with a Gaussian form

[TABLE]

and

[TABLE]

where $\{\tilde{C}_{E},\tilde{C}_{L},\tilde{D}_{E},\tilde{D}_{L}\}$ are coefficients evaluated at $(E_{0},{\textbf{{L}}}_{0},t-t_{0})$. If the fluctuations of the external tidal field are much weaker than the binding forces all coefficients approach zero, and the probability functions (51) and (52) become sharply peaked about $E=E_{0}$ and ${\textbf{{L}}}={\textbf{{L}}}_{0}$, respectively. As expected, in the limit $\langle|\Delta{\textbf{{V}}}|^{2}\rangle\to 0$, one has that $p(E,t|E_{0},{\textbf{{L}}}_{0},t_{0})\to\delta(E-E_{0})$, and $p({\textbf{{L}}},t|E_{0},{\textbf{{L}}}_{0},t_{0})\to\delta({\textbf{{L}}}-{\textbf{{L}}}_{0})$, which recovers the initial conditions.

The diffusion coefficients are derived from phase-space averages of (47) and (48) over a collection of tracer particles with the same combination of energy and angular momentum at $t_{0}=0$. Using the adiabatically-corrected velocity impulses (40) leads to drift coefficients

[TABLE]

and diffusion coefficients

[TABLE]

with upper bars denoting a phase-space average over tracer particle ensembles (see Appendix C). Note that the coefficients are isotropic in the angular momentum space, i.e. $C_{i}(E,{\textbf{{L}}})=C_{i}(E,L)$ and $D_{i}(E,{\textbf{{L}}})=D_{i}(E,L)$, with a subindex $i=\{E,L\}$. In addition, the cross-coefficients vanish by symmetry, $\overline{\Delta E\Delta{\textbf{{L}}}}=\tilde{D}_{EL}=0$, which justifies the separability of $p(E,{\textbf{{L}}},t|E_{0},{\textbf{{L}}}_{0},t_{0})$ implicitly assumed in Equation (50).

In an impulsive regime, averages over tracer-particle (upper bars) and substructure (brackets) ensembles in (53) and (54) become statistically independent, which simplifies our mathematical treatment greatly. Indeed, for tidal fluctuations with a characteristic duration (30) that is much shorter than the orbital frequency, i.e. $wT_{ch}\ll 1$, one can approximate $\langle\Lambda^{2}TA\rangle\approx\langle\Lambda^{2}T\rangle$, which is independent of the phase-space distribution of tracer particles. Hence, Equations (53) and (54) become

[TABLE]

4.1.2 Diffusion in a confined region

In a self-gravitating potential $\Phi_{s}$, gravitationally-bound particles can only diffuse in a limited volume of the integral-of-motion space. In particular, energies must be negative definite, $E<0$, while the angular momentum is confined within the interval $0\leq L\leq L_{c}(E)$, where $L_{c}(E)$ corresponds to the angular momentum of a circular orbit with energy $E$. In order to confine particles within a given domain we must solve the diffusion equation (49) with appropriate boundary conditions. To this aim, first we set an absorbing boundary at $E=0$, such that $p(E\geq 0)=0$, which implies that particles with $E\geq 0$ escape from a gravitating system to never return333This approximation is only approximately correct. In reality, particles can diffuse in and out of the boundary $E=0$. Although this calls for introducing a permeable boundary, the mathematical treatment significantly more involved (e.g. Carslaw & Jaeger 1986) and must be left for follow-up work.. The angular momentum space is treated as a cubic box with reflecting surfaces at the limits of each dimension444Note that diffusion processes depend on the space geometry. Here we adopt Cartesian symmetry for mathematical convenience, see Appendix B for details., such that $\partial_{L}p(L_{i}=\pm L_{c})=0$, where $i=x,y,z$, which limits the values of the components within $-L_{c}\leq L_{i}\leq L_{c}$. In order to obey these boundary conditions, Appendix B shows that Equations (51) and (52) must be respectively replaced by

[TABLE]

and

[TABLE]

where $L_{c}=L_{c}(E_{0})$, $\lambda_{nml}=\pi^{2}(n^{2}+m^{2}+l^{2})/(4L_{c}^{2})$ is a separation constant, and $\alpha_{000}=1/8$, $\alpha_{n00}=\alpha_{0m0}=\alpha_{00l}=1/4$, $\alpha_{nm0}=\alpha_{n0l}=\alpha_{0ml}=1/2$, and $\alpha_{nml}=1$  for $n,m,l\geq 1$. As demonstrated in Appendix B, the function $p({\textbf{{L}}},t|E_{0},{\textbf{{L}}}_{0},t_{0})$ recovers the Gaussian propagator (52) on short time-scales, $\pi^{2}D_{L}/(4L_{c}^{2})\,t\ll 1$, where $D_{L}=\tilde{D}_{L}/t$ is a ‘static’ coefficient defined in §4.2. As the time progresses, however, all modes with $n,m,l\geq 1$ decay exponentially, and the probability function converges asymptotically towards a constant value $p({\textbf{{L}}},t|E_{0},{\textbf{{L}}}_{0},t_{0})\to 1/(8L_{c}^{3})$ in the limit $t\to\infty$. This is an important result, as it shows that stochastic tidal heating tends to isotropize the initial angular momentum distribution on a time-scale $\pi^{2}D_{L}/(4L_{c}^{2})\,t\gtrsim 1$. This calls for defining the isotropization time-scale as

[TABLE]

where $D_{L}=\tilde{D}_{L}/t$ is given by (54). In §4.3 we discuss the properties of $t_{\rm iso}$ in some detail.

According to the Jeans theorem, the initial equilibrium state of a gravitational system is fully defined by the distribution of particles in the integral-of-motion space $N(E_{0},{\textbf{{L}}}_{0},t_{0})$, which determines the probability to find a particle in the 4-dimensional volume element ${\rm d}E_{0}{\rm d}^{3}L_{0}$ centred at a specific energy $E_{0}=V^{2}/2+\Phi_{s}(R)$ and an angular momentum vector ${\textbf{{L}}}_{0}={\textbf{{R}}}\times{\textbf{{V}}}$ at the time $t_{0}=0$. The (non-equilibrium) state of a gravitating system at a later time $t>t_{0}$ is found by convolving the initial distribution $N(E_{0},{\textbf{{L}}}_{0},t_{0})$ with (56) and (57), which yields

[TABLE]

hence exploiting the fact that $p(E,{\textbf{{L}}},t|E_{0},{\textbf{{L}}}_{0},t_{0})$ is a Green’s function (see P15 for details). In Section 6 we discuss the evolution of planetary orbits driven by stochastic tidal heating as an illustration of the Green’s convolution (59).

4.2 Fokker-Planck equations

Let us write down the Fokker-Planck equations in the integral-of-motion space in order to gain insight onto the physical meaning on the diffusion coefficients. To this end let us assume that the transformation $N(E_{0},L_{0},t_{0})\to N(E,L,t)$ takes place during a short, but non-zero, time interval $t-t_{0}$, and that the variation of energy and angular momentum is small, such that $|\Delta E|/|E|\ll 1$ and $|\Delta{\textbf{{L}}}|/L\ll 1$. Thus, expanding the left-hand side of Equation (59) up to the first order in $t$, and the right-hand side as a Taylor series to the second order in $\Delta E$ and $\Delta{\textbf{{L}}}$ yields (Spitzer 1987; P15)

[TABLE]

The Fokker-Planck equations (60) simplify considerably in an isotropic tidal field, where $\overline{\langle\Delta{\textbf{{L}}}\rangle}=\overline{\langle\Delta E\Delta{\textbf{{L}}}\rangle}=0$ and $2\tilde{D}_{L}=\overline{\langle|\Delta{\textbf{{L}}}|^{2}}\rangle$ by symmetry. Also, from Equations (47) and (48) it is clear that $\langle(\Delta E)^{2}\rangle\sim V^{2}\langle|\Delta\textbf{{V}}|^{2}\rangle$, whereas $\langle\Delta E\rangle^{2}\sim\langle|\Delta\textbf{{V}}|^{2}\rangle^{2}$. Since the Brownian approach requires $\langle|\Delta\textbf{{V}}|^{2}\rangle\ll V^{2}$, one can safely approximate $2\tilde{D}_{E}=\overline{\langle(\Delta E)^{2}\rangle}-\overline{\langle\Delta E\rangle}^{2}\approx\overline{\langle(\Delta E)^{2}\rangle}$, such that

[TABLE]

The linear dependence of the coefficients $\tilde{C}$ and $\tilde{D}$ with the length of the time-interval $t$ in (53) and (54), respectively, calls for the definition of static drift and diffusion coefficients

[TABLE]

such that Equation (61) becomes

[TABLE]

Random fluctuations of the tidal field induce a bulk motion of tracer particles in the integral-of-motion space. At leading order, the main variation corresponds to a flux of particles crossing the energy layer $E$ at a fixed angular momentum L, which can be estimated from Equations (63) and (62) as

[TABLE]

the negative sign of $J$ implies that stochastic tidal heating leads to a steady flow of particles drifting from bound energies $E<0$ towards $E\to 0$. This realization has important consequences for the survival of self-gravitating objects in a fluctuating tidal field, as discussed below. Note also that in an impulsive regime Equation (64) has a simple form

[TABLE]

The fact that $J_{\rm imp}\propto\overline{R^{2}}$ shows that the particle flow is particularly strong in the outskirts of the potential $\Phi_{s}(R)$.

Over a sufficiently long interval of time, particles with a combination of integrals $(E_{0},{\textbf{{L}}}_{0})$ at $t_{0}=0$ will gain sufficient kinetic energy as to escape the potential $\Phi_{s}(R)$, which is known as “tidal evaporation” (Spitzer 1958 and references therein). The unbinding time, $t_{\rm esc}(E_{0},L_{0})$, can be estimated from Equation (62) by expressing $\overline{\langle\Delta E\rangle}=-C_{E}(E_{0},L_{0})t=\frac{1}{2}\overline{R^{2}\langle\Lambda^{2}TA\rangle}t$. Hence, setting $E(t_{\rm esc})=0$ yields an average escape time

[TABLE]

In the impulse approximation ($A\approx 1$), Equation (66) reduces to

[TABLE]

which is independent of the velocity distribution of the tracer particles.

4.3 Example: Keplerian potential

It is worth illustrating the above results with an analytical example. To this end, consider an ensemble of tracer particles orbiting in a Keplerian potential

[TABLE]

As in previous Sections, we assume that the system is surrounded by an homogeneous distribution of Hernquist (1990) spheres with a mass $M$, scale-length $c$ and number density $n$ (see §2.1 for details).

The drift and diffusion coefficients appearing in (53) and (54), respectively, can be calculated analytically under the assumption that (i) random tidal interactions occur in an impulsive regime ($A\approx 1$), and (ii) substructures do not spatially overlap with each other ($c\ll D$). Hence, combination of (55), (163), (164), (166) and (34) yields

[TABLE]

Here, $a$ and $e$ are the semi-major axis and eccentricity of and orbit with energy $E$ and angular momentum $L$

[TABLE]

where $e=0$ corresponds to a circular orbit with angular momentum

[TABLE]

Combination of (66),  (69) and (70) yields the mean impulsive time-scale required to unbind a tracer particle as a function of semi-major axis and orbital eccentricity

[TABLE]

This expression reveals a number of interesting features. Notice first that the unbinding time-scale is mainly determined by the ratio between the semi-major axis of the orbit and the average separation between substructures, $t_{\rm esc}\sim(D/a)^{3}$, and depends only weakly on eccentricity. In particular circular orbits ($e=0$) remain bound a factor $5/2$ longer than radial orbits $(e=1)$. In addition, Equation (72) shows that the escape time-scale decreases linearly with the velocity dispersion of the substructure population, $t_{\rm esc}\sim\langle v^{2}\rangle^{1/2}$, which highlights the fact that impulsive energy injections are dominated by nearby substructures moving with a small relative velocity (see §3.3).

It is interesting to compare the time-scale required to isotropize the angular momentum distribution, Equation (58), against the unbinding time-scale (66). Adopting the impulse approximation, and inserting (71), (164) and (34) into (58) yields

[TABLE]

which indicates that isotropization and unbinding time-scales share the same power-law behaviour, albeit with the former being more sensitive to orbital eccentricity than the latter. In particular, dividing (73) by (72) yields a ratio

[TABLE]

which has values between $12/\pi^{2}\simeq 1.22$ and $80/(21\pi^{2})\simeq 0.386$ for $e=0$ and $e=1$, respectively, with $t_{\rm iso}=t_{\rm esc}$ occurring at $e\simeq 0.256$. This result implies that orbits with low orbital eccentricity tend to escape from the potential $\Phi_{s}$ before the momentum distribution is fully randomized. We will return to this issue below.

5 $N$-body tests

In this Section, we run a number of restricted N-body experiments that follow the dynamical evolution of tracer particles subject to random tidal interactions with a large population of extended substructures. Our goal is to test the main assumptions introduced in previous Sections in order to describe tidal heating as a random walk (a.k.a. Brownian motion) of velocities that leads to a diffusion process a 4-dimensional integral-of-motion space.

5.1 Numerical set-up

To simplify our numerical analysis and speed up $N$-body computations we adopt an analytical Dehnen (1993) model for the host galaxy potential, where

[TABLE]

In what follows, we use $N$-body units with $G=\rho_{0}=r_{0}=1$.

Substructures are assumed to be spherically distributed about the host centre following a number density profile

[TABLE]

which is normalized to $4\pi\int_{0}^{\infty}{\rm d}r\,r^{2}\,n(r)=N$. The scale on which the number density profile varies can be calculated analytically from (76) as

[TABLE]

The local approximation is justified at radii where $d$ is larger than the mean separation between substructures, i.e. $d\gg D=(2\pi n(r)]^{-1/3}$ (see Appendix A of Paper I). For $\alpha>1$, the distance $d$ diverges as $r\to 0$, hence one can safely approximate $n\approx n_{0}$ and $D\approx(2\pi n_{0})^{-1/3}$ at $r\ll r_{0}$.

If one further assumes that substructures move isotropically about the centre of the potential, then the distribution function can be written as (Eddington 1916)

[TABLE]

where $n[r(\Phi_{g})]$ corresponds to the profile (76) expressed as a function of the potential (75). Dynamical equilibrium is guaranteed if and only if the distribution function is definite positive, $f\geq 0$, within the energy range sampled by the orbital distribution. For the models described above this condition requires a steep slope at large radii. In our experiments we set $\gamma=0$, $\alpha=3$ and $\beta=60$. Although not shown here, we have explicitly checked that the particular choice of indices $\gamma$, $\alpha$ and $\beta$ does not significantly change the results shown below. The initial radii and velocities of $N$ substructures are drawn randomly from the distribution function $f(E)$ using a rejection algorithm, while the directions of the position and velocity vectors are isotropically distributed over the surface of a sphere. The velocity dispersion $\langle v^{2}\rangle$ is measured from the relative velocities of the substructure ensemble. For the experiments shown below we find $\langle v^{2}\rangle^{1/2}\simeq 0.45$. The orbits of individual substructures in the potential $\Phi_{g}$ are integrated forward in time using a leap-frog algorithm whose time-step is chosen such that the energy is conserved at least at a $10^{-9}$ level in isolation.

For simplicity, tracer particles orbit an isolated Keplerian potential (68) and follow circular orbits at $t=0$, such that the initial tangential velocity is set to $V=V_{c}(R)=|R\nabla\Phi_{s}(R)|^{1/2}=(Gm/R)^{1/2}$, and the orbital frequency to $w=V/R=(Gm/R^{3})^{1/2}$. The initial energy and angular momentum are therefore $E_{0}=V_{c}^{2}/2+\Phi_{s}(R)$ and $L_{0}=RV_{c}(R)$, while the semi-major axis and eccentricity are $a_{0}=-Gm/(2E_{0})$ and $e_{0}=0$, respectively.

Fluctuations of the external tidal field are added into the equations of motion following two different prescriptions

Direct-force computation of the tidal field. Here the combined tidal tensor induced by a set of $N$-Hernquist (1993) spheres is calculated as a direct summation of the tidal forced induced by individual substructures

[TABLE]

where $\delta_{jk}$ is the Kronecker delta and ${\textbf{{r}}}^{\prime}_{i}=\{x_{i}^{j}\}_{j=1,2,3}={\textbf{{r}}}_{s}-{\textbf{{r}}}_{i}$ is the relative distance between the self-gravitating object and the $i^{th}$ substructure. In the models shown below, the tidal tensor (79) is evaluated at a radius, $r_{s}=0.1r_{0}$, where the number density of substructures (76) is approximately constant, $n\approx n_{0}$. It is worth stressing that this approach comes at a large computational cost. 2. 2.

Monte-Carlo sampling the distribution of velocity increments $\Psi({\textbf{{V}}},\Delta{\textbf{{V}}},\Delta t)$ given by Equation (31). Here we use adiabatic-corrected coefficients $\langle|\Delta{\textbf{{V}}}|^{2}\rangle$ computed from equations (40) and (41), with a time interval $t$ set equal to the time-step of the leap-frog scheme, $\Delta t$. In order to incorporate the effects of a fluctuating tidal field, first we compute the orbit of the test particle in the potential $\Phi_{s}$ from $\textbf{{R}}(t)\to\textbf{{R}}(t+\Delta t)$ by solving the equations of motion

[TABLE]

Subsequently, we add a random velocity ‘kick’ to the velocity vector computed from (80) in isolation, ${\textbf{{V}}}(t+\Delta t)$, such that

[TABLE]

where ${\it Ran}(\Delta{\textbf{{V}}})$ correspond to random velocity increments drawn from the isotropic Gaussian distribution $\Psi({\textbf{{V}}},\Delta{\textbf{{V}}},\Delta t)$.

The potential $\Phi_{s}$ is assumed to be in isolation, i.e. the host galaxy potential is set to $\Phi_{g}=0$, which implies that particles with a specific energy $E\geq 0$ will escape from the system. This choice helps to isolate the effect of stochastic tidal fluctuations by removing the smooth component of the tidal tensor from the equations of motion ($T_{g}=0$). Equations (4) and (80) are solved using a Runge-Kutta scheme (e.g. Press et al. 1992) with a fixed time-step $\Delta t=0.01w^{-1}$. For the models shown below, this yields an energy conservation at a $\sim 10^{-9}$ level in the absence of stochastic tidal fluctuations. One should bear in mind is that at a fixed time-step direct calculations of the force induced by a population of $N-$substructures consumes approximately a factor $N$-times more CPU-power than the Monte-Carlo method. This difference becomes astronomical in the case of Galactic substructures, as discussed in Section 7.

5.2 Tidal force fluctuations

We start by testing the analytical expressions for the distribution of tidal force fluctuations, p(\mbox{\boldmath\Lambda}), derived in Section 2. Recall that our analytical framework assumes that (i) the Euler and Coriolis terms appearing in the non-inertial rest frame of the test particle can be neglected, and that (ii) individual substructures induce tidal forces that point in random directions. In order to quantify the amplitude of tidal fluctuations, we first measure the modulus of the combined tidal vector $\Lambda$ acting on a single tracer particle from a direct summation of the tidal forces (7) induced by individual subtructures

[TABLE]

The average of the variance within an interval of time $t$ is calculated as

[TABLE]

Fig. 3 shows the time evolution of $\langle\Lambda^{2}\rangle_{t}$ computed from $N_{\rm ens}=16$ independent ensembles of $N=5000$ substructures with a mass $M=10^{-7}$ and three different sizes ($c$) given in units of the average separation between substructures ($D$). The ‘saw’-like shape of the curves is caused by impulse variations of $\Lambda^{2}(t)$ over reduced periods of time. As expected, on time-scales $t\gtrsim T_{\rm ch}$, where $T_{\rm ch}$ is the average duration of tidal fluctuations given by (30), ensemble-averaged values of the variance, $\langle\Lambda^{2}\rangle_{\rm ens}(t)=N_{\rm ens}^{-1}\sum_{k=1}^{N_{\rm ens}}\langle\Lambda^{2}\rangle_{t,k}$, converge towards the analytical expression for $\langle\Lambda^{2}\rangle$ derived in the limit $N\to\infty$ (black-dashed lines, Equation 16). However, on time-intervals $t\lesssim T_{\rm ch}$ Equation (16) systematically overestimates the value of $\langle\Lambda^{2}\rangle_{t}$. The reason for this mismatch is simple to understand: as discussed in §2 (see also Paper I), the variance $\langle\Lambda^{2}\rangle$ is completely dominated by the contribution of the nearest object. If the time interval is very short, $t\lesssim T_{\rm ch}$, the positions of substructures do not change appreciably, and the spectrum of tidal fluctuations $p(\Lambda)$ is sharply peaked about the mean, $\langle\Lambda\rangle$. As the time progresses, $t\gg T_{\rm ch}$, stochastic fluctuations begin to sample the large-force tail of the distribution, $\Lambda\sim\lambda_{0}\gg\langle\Lambda\rangle$, which leads to an exponential convergence of $\langle\Lambda^{2}\rangle_{t}$ towards the asymptotic value given by Equation (16).

Comparison between the three panels of Fig. 3 highlights the strong dependence of the magnitude of the tidal fluctuations on the size of substructures. Indeed, according to the analytical expression given by Equation (16), the force variance scales as $\langle\Lambda^{2}\rangle\sim(D/c)^{3}$ for $c/D\ll 1$. The scatter about the ensemble-average value also increases significantly as the size-to-separation ratio $c/D$ decreases, which can be traced back to the fact that $\langle\Lambda^{2}\rangle_{t}$ is largely governed by the few substructures that come closest to the tracer particle within the time interval $t$. Note that for spatially-overlapping populations ($c\gtrsim D$) there is a systematic mismatch between the numerical and analytical values. However, such deviations are to be expected given the analytical expression (16) is only strictly valid for ’rarefied’ populations where $c/D\ll 1$ (for details, see Paper I).

Fig. 3 also reveals that the convergence of $\langle\Lambda^{2}\rangle_{t}$ towards Equation (16) shifts to systematically later times as the size of individual substructures shrinks with respect to their average separation. Given that the time interval $t$ is given in units of the average duration of tidal fluctuations, $T_{\rm ch}$, the x-axis represents an average number of fluctuations within a given time span. Therefore, one may conclude that the average number of fluctuations required to reach the asymptotic behaviour of $\langle\Lambda^{2}\rangle_{\rm ens}(t)$ increases as the ratio $c/D$ decreases. To quantify the convergence time we fit the numerical curves with an exponential function (red lines)

[TABLE]

where $\langle\Lambda^{2}\rangle_{\infty}$ corresponds to the asymptotic value of the variance in the limit $t\to\infty$, $\delta_{s}(t)=1-\exp(-t/\tau_{s})$ is the so-called sampling delay function, and $\tau_{s}$ is defined as the sampling time-scale. On long time-scales, $t\gtrsim\tau_{s}$, the variance of tidal fluctuations converges asymptotically towards $\langle\Lambda^{2}\rangle_{\infty}\approx\langle\Lambda^{2}\rangle$, which suggests that $\tau_{s}$ measures the average time required to sample the large-force tail of $p(\Lambda$). Fig. 4 shows that the sampling time can be approximately written as $\tau_{s}\approx 4.8T_{\rm ch}(D/c)^{2}$, and that this relation is largely independent of the number of substructures, $N$. According to Equation (12), the quadratic slope reflects the fact that the sampling time-scale is inversely proportional to the probability to find the closest substructure at a distance $r\sim c$. Indeed, from Equation (12) we find that $p(r)\sim r^{2}$ for $c/D\ll 1$ ($4\pi c^{3}n/3\ll 1$). For reasons that will become evident in §7, it is useful to define the sampling frequency

[TABLE]

Note that the sampling time becomes arbitrary long in the particle limit, i.e. $\tau_{s}\to\infty$ ($w_{s}\to 0$) as $c\to 0$, which means that the amplitude of tidal fluctuations induced by a background of point-masses never converges, and will therefore reach arbitrarily-large values in the limit $t\to\infty$.

5.3 Diffusion in the integral-of-motion space

In this Section we follow the evolution of binary ’stars’ subject to the tidal fluctuations computed in §5.2. For simplicity, tracer particles are set on circular orbits at $t_{0}=0$, which reduces the dimensionality of the problem to the initial (specific) energy, $E_{0}=-Gm/(2a_{0})$.

Fig. 5 illustrates how random velocity increments acquired by tracer particles over short intervals of time lead to a Brownian motion in the integral-of-motion space. The energy and angular momentum of single tracer particle varies randomly as a response to the fluctuating tidal field generated by ensembles of $N=5000$ Hernquist spheres with a mass $M=10^{-7}$ and a size $c=0.1D$, separated by an average distance $D=(2\pi n_{0})^{-1/3}$. Thin lines correspond to different random realizations of the substructure population. At $t_{0}=0$ tracer particles have a semi-major axis $a_{0}=10D$ and an angular momentum $L_{0}=L_{c}(E_{0})=\sqrt{Gma_{0}}$, which we mark with a red dot for ease of reference. As a response to the fluctuating tidal field, on average particles tend to gain kinetic energy, thus ‘drifting’ towards the energy threshold $E=0$. Simultaneously, a random walk in angular momentum takes the particles away from the initial location on the circular-velocity curve $L_{c}(E)$ (solid black line). Thus, scattered orbits gradually become more radial and less energetically bound.

In order to test the impulse approximation and the adiabatic correction explored in Section 3.3, we integrate the evolution of tracer particles for a time-interval $w\,t_{\rm max}=50000$, where $w=2\pi/P=\sqrt{Gm/a_{0}^{3}}$ is the frequency of a circular orbit. Left, middle and right panels of Fig. 6 correspond to orbits with initial semi-major axes, $a_{0}/D=10,1.0$ and 0.1. These models cover $\sim 3$ orders of magnitude in orbital frequencies, $wT_{\rm ch}\simeq 0.03,1.0$ and 31.6, respectively, where $T_{\rm ch}$ is the typical duration of tidal fluctuations (30). With this choice of parameters, particles with the largest/smallest semi-major axis react impulsively/adiabatically to the fluctuations of the tidal field. In units of the (impulsive) escape time (67), our models cover a much larger range of values, $w\,t_{\rm esc}\simeq 1.3\times 10^{3},4.0\times 10^{7}$ and $1.3\times 10^{12}$ for $a_{0}/D=10.0,1.0$ and 0.1, respectively. Hence, of the three orbits explored here, only particles with $a_{0}=10D$ are evolved for a sufficiently long time to escape from the system (i.e. $t_{\rm esc}<t_{\rm max}$). Recall that on time-intervals $t\ll t_{\rm esc}\simeq t_{\rm iso}$ particles are not “aware” of the finite boundaries of the integral-of-motion space (see Appendix B), hence the evolution of energy and angular momentum is expected to follow the relation exhibited by free-diffusion processes in an unconfined volume, $\sigma^{2}=2Dt$. To test this theoretical expectation, we plot the time-evolution of the energy variance, $\sigma_{E}^{2}=\langle\Delta E^{2}\rangle-\langle\Delta E\rangle^{2}$, and angular momentum variance, $\sigma_{L}^{2}=\langle|\Delta\textbf{{L}}|^{2}\rangle-\langle|\Delta\textbf{{L}}|\rangle^{2}$. Values derived from a direct-force summation (thin coloured lines) are compared against the theoretical predictions derived under the impulse approximation (black-dotted lines) and the adiabatically-corrected values (black-dashed lines).

The first noteworthy result is the convergence of $\sigma_{E}^{2}(t)$ and $\sigma_{L}^{2}(t)$ computed from direct-force $N$-body models towards the free-diffusing relation $\sigma^{2}=2Dt$. As expected from Section 5.2, numerical curves approach Equation (54) on time-scales $t\gtrsim\tau_{s}$, where $\tau_{s}$ is the sampling time-scale plotted in Fig. 4. This can be better seen by looking at the upper axis plots, which measures the time interval in units of the sampling time, $\tau_{s}$. Notice that on very short time-scales, $t\ll\tau_{s}$, the time-evolution of $\sigma_{E}^{2}$ and $\sigma_{L}^{2}$ undergoes periodic oscillations rather than a monotonic increase. These oscillations are particularly prominent in the upper-right panel, which corresponds to orbits with such a high orbital frequency that particles perform $\sim 10^{4}$ revolutions around the potential $\Phi_{s}$ before the location of the closest substructures changes appreciably. Hence, for these objects the external tidal field undergoes a “parallax” oscillation in phase with their orbital motion, which causes cyclic variations of energy and, to a lesser extent angular momentum. Once the spectrum of tidal fluctuations $p(\Lambda)$ is thoroughly sampled and convergence has been reached, we find that the adiabatic-corrected values of $\sigma_{E}^{2}$ and $\sigma_{L}^{2}$ derived from Weinberg (1994a,b,c) formula (black-dashed lines) show an excellent match to all the models explored here. In contrast, the values derived under the impulse approximation largely overestimate the increments of energy and angular momentum. This is particularly evident for ’adiabatic’ models with $w\,T_{\rm ch}\gg 1$ (right panels).

Finally, it is worth highlighting the accelerated increase of energy and angular momentum variance exhibited by the ‘impulsive’ ($wT_{\rm ch}\ll 1$) models at late times. In particular, particles with a large semi-major axis $a_{0}=10D$ show a rapid growth of $\sigma_{E}^{2}(t)$ and $\sigma_{L}^{2}(t)$ at $w\,t\gtrsim 10^{2}$ that cannot be reproduced by Equation (54). Although not shown here, we find that the regime of rapid heating coincides with the unbinding of test particles as the orbital energy approaches the boundary $E=0$. Interestingly, Monte-Carlo models show that sampling velocity ’kicks’ from the probability function (31) (red solid lines) successfully describe the average evolution of $\sigma_{E}^{2}(t)$ and $\sigma_{L}^{2}(t)$ of all our models. Remarkably, these models also reproduce the accelerated heating experienced by test particles as they escape from the system.

6 Application to planetary orbits

This Section illustrates the theoretical framework described in §4 by following the evolution of planetary discs composed of test particles (or ‘comets’) orbiting around a Keplerian potential subject to stochastic fluctuations of an external tidal field. Tidal heating is modelled following two different approaches: (i) Monte-Carlo sampling velocity impulses as described in §5, and (ii) convolving the initial energy–angular momentum distribution with Green’s propagators. For convenience, we choose units where $G=m=D=\langle v^{2}\rangle=1$. In these units, the mass and size of substructures are set to $M=10^{-3}$ and $c=0.1$, respectively.

For simplicity, we adopt idealized conditions at $t_{0}=0$ where particles move on a razor-thin disc with energies homogeneously distributed over an interval $E_{0}\in(E_{1},E_{2})$, and an angular momentum vector that points along the (positive) $z$-axis, $\textbf{{L}}=(0,0,L_{c}[E])$. The initial state of the ensemble is thus fully specified by the probability function

[TABLE]

where $N_{0}=1/(E_{2}-E_{1})$ is a normalization constant which guarantees that $\int_{E_{1}}^{E_{2}}{\rm d}E_{0}\int{\rm d}^{3}L\,N(E_{0},{\textbf{{L}}}_{0},t_{0})=1$. Equation (85) implies that at $t=t_{0}$ comets move on circular orbits on the plane $z=0$ within a radial range $R\in[R_{1},R_{2}]=[Gm/(-2E_{1}),Gm/(-2E_{2})]$.

6.1 Evolution in the integral-of-motion space

The non-equilibrium distribution of integrals at $t>t_{0}$ can be directly derived from the convolution of the initial distribution (85) with the Green’s functions (56) and (57), which yields

[TABLE]

with coefficients $\tilde{C}_{E}=\tilde{C}_{E}(E_{0},L_{c},t-t_{0})$, $\tilde{D}_{E}=\tilde{D}_{E}(E_{0},L_{c},t-t_{0})$ and $\tilde{D}_{L}=\tilde{D}_{L}(E_{0},L_{c},t-t_{0})$ and an integration constant $\lambda_{nml}=\pi^{2}(n^{2}+m^{2}+l^{2})/(4L_{c}^{2})$. Equation (86) simplifies considerably if we assume that the distribution of angular momentum in the $x$ and $y-$directions are tightly peaked around the initial null value, and thus far from the reflecting barriers at $\pm L_{c}$. Under this approximation, propagators in the $L_{x}$ and $L_{y}$ dimensions have a Gaussian form (52) (see Appendix B), and Equation (86) reduces to

[TABLE]

where $\lambda_{\ell}=\pi^{2}\ell^{2}/(4L_{c}^{2})$ and $\cos(\ell\pi)=(-1)^{\ell}$. In order to gain physical intuition into the physical behaviour of $N(E,{\textbf{{L}}},t)$, it is useful to reduce the dimensionality of the probability function by marginalizing over each integral dimension.

6.1.1 Energy space

Let us first calculate the distribution of energy $N(E,t)$ by integrating each side of Equation (87) over ${\rm d}^{3}L=2\pi L_{R}\,{\rm d}L_{R}\,{\rm d}L_{z}$, where $L_{R}=\sqrt{L_{x}^{2}+L_{y}^{2}}$ is the planar component of the angular momentum, which yields

[TABLE]

If the variation of the radial component of the angular momentum is small, $|\Delta L_{R}|\ll L_{c}$, the integral over ${\rm d}L_{R}$ takes the value of unity

[TABLE]

whereas, by construction, the integral over ${\rm d}L_{z}$ is normalized to unity

[TABLE]

Thus, Equation (88) reduces to a single integral

[TABLE]

where $p(E,t|E_{0},L_{c},t_{0})$ is the energy propagator (56) with coefficients computed for circular orbits, $\tilde{C}_{E}=\tilde{C}_{E}(E_{0},L_{c},t-t_{0})$ and $\tilde{D}_{E}=\tilde{D}_{E}(E_{0},L_{c},t-t_{0})$. The role of the negative right-hand term of (90) must be highlighted, as it causes the removal of comets that cross the absorbing barrier at $E=0$ (see Appendix B). This means that the total number of bound comets $\int_{-\infty}^{0}{\rm d}E\,N(E,t)$ decreases steadily with time, as it is shown in more detail below.

6.1.2 Angular momentum space

Similarly, the probability to find comets with an angular momentum component between $L_{R},L_{R}+{\rm d}L_{R}$ at the time $t$ can be calculated as

[TABLE]

whereas the vertical component of angular momentum has a distribution

[TABLE]

here

[TABLE]

is the one-dimensional angular momentum propagator (144) evaluated at $L_{z,0}=L_{c}$ with a coefficient $\tilde{D}_{L}=\tilde{D}_{L}(E_{0},L_{c},t_{0})$.

Unfortunately, the distributions $N(E,t),N(L_{R},t)$ and $N(L_{z},t)$ must be solved numerically owing to the non-trivial dependence of the coefficients $\tilde{C}_{E},\tilde{D}_{E}$ and $\tilde{D}_{L}$ on $E_{0}$. Below we illustrate the time-evolution of angular momentum and energy via comparison against Monte-Carlo $N$-body models.

6.1.3 Monte-Carlo $N$-body experiments

Fig. 7 plots three snap-shots in the evolution of $10^{5}$ tracer particles initially distributed on a disc at $z=0$ within a radial interval $R\in(R_{1},R_{s})=(0.1,4)$. As time progresses, the effects of tidal fluctuations become particularly strong in the outskirts of the disc, where orbits can be scattered off to large distances from the disc plane, $|z|/D\gtrsim 100$. Interestingly, the resulting spatial distribution at late times is reminiscent to that of the Oort cloud in the Solar system (e.g. Oort 1950; Hills 1981). Note also that the number of test particles at small radii $R\ll R_{1}$ also increases with time. As we will see below, these comets have been scattered onto radial, or even retrograde motions.

A sharp drop in the number of particles beyond $R\gtrsim 100$ is also noticeable at late times. Indeed, at such large distances most comets are gravitationally unbound from the central potential, thus moving on parabolic trajectories away from the system. This is clearly visible in the distribution of orbital energies, $N(E,t)$, plotted in the second row of Fig. 7, which shows a long probability tail that extends across $E=0$ (marked with vertical dotted lines) towards positive energies at intermediate and late times. The fraction of particles crossing the energy layer $E=0$ increases very rapidly on short time scales, $t\lesssim 100$, and slowly decline thereafter. A detailed study of the time-evolution of the escape rate will be presented in Section 6.2. Comparison of the energy distribution with the analytical expression (90) (red lines) shows good agreement at low energies ($E\lesssim 0$), which progressively worsens in the proximity of the barrier at $E=0$, a region dubbed the ‘fringe’ by Spitzer & Shapiro (1972). As pointed out by these authors, “Fokker-Plank equations are valid throughout most of the [system], but not in the fringe”. Indeed, both Fokker-Planck and diffusion equations rely on a Taylor-expansion of Einstein’s master equation for small increments of energy $|\Delta E|\ll|E|$ within a finite time interval $t>0$ (e.g. see Sections 2.2 and 3.3 of P15), a condition that does not hold at $E\sim 0$. For the same reason, diffusion equations do not reproduce the ‘accelerated’ unbinding of particles found in Monte-Carlo models (see left panel of Fig. 6), which in turn leads to an overprediction of the number of comets in the fringe. Yet, the largest discrepancy between diffusion and Monte-Carlo models is found in the region of positive energies. Here, Monte-Carlo models show a non-zero probability to find particles at $E\gtrsim 0$ at intermediate and late times. By construction, the tail at positive energies cannot be reproduced by diffusion models with an absorbing barrier at $E=0$, which sets $N(E,t)=0$ for $E\geq 0$ (see Appendix B).

The third row of Fig. 7 shows that the radial angular momentum distributions $N(L_{R},t)$ widens progressively as the time interval increases. Given that at $t_{0}=0$ all particles move on a disc $z=0$, such that $N(L_{R},t_{0}=0)=\delta(L_{R})$, non-zero values of $L_{R}$ must be caused by particles scattered off the initial plane. To see this more clearly, let us define the angle $\cos\theta=\hat{\textbf{{L}}}\cdot\hat{z}$, such that $\theta_{0}=0$ at $t=t_{0}$. Given that the maximum vertical offset from the disc plane goes as $z_{\rm max}=a\sin\theta$ for $e\approx 0$, we find that $z_{\rm max}/a=L_{R}/L$, hence showing that the thickening of the outer disc is related to the widening of the distribution $N(L_{R},t)$. It is also interesting to observe that the excess of particles above the analytical red curve coincides with highest evaporation rates. Indeed, particles that experience large velocity kicks also tend to gain angular momentum, which explains the correlation between the fringe population and the excess of probability in the large-angular momentum tail of the distribution.

The distribution of the rotational angular momentum, $N(L_{z},t)$, is plotted in the bottom panels of Fig. 7, which show a number of remarkable features. Note first that the initial truncations at $L_{z,min}=\sqrt{GmR_{1}}\approx 0.32$ and $L_{z,max}=\sqrt{GmR_{2}}=2$ are quickly erased by tidal perturbations. In particular, the tail extending towards positive angular momentum corresponds to particles scattered into the fringe region, which is poorly described by the analytical expression (92). More interestingly, at intermediate and late times we find a second tail at low, or even negative values of $L_{z}$. These are comets that move on nearly radial or counter-rotating orbits which periodically plummet into the inner regions of the potential from the outskirts of the planetary system, a process discussed in detail by Hills (1981). The fraction of comets on retrograde orbits correlates with the number of particles escaping from the system. We return to this issue in §6.3. Finally, the experiments shown in Fig. 7 confirm that particles become energetically unbound before the isotropization of the angular momentum distribution is complete. This is in agreement with Equation (74), which shows that the isotropization of circular orbits occurs on longer time-scales than the average time that these particles remain bound to the Keplerian potential, i.e. $t_{\rm iso}(e=0)>t_{\rm esc}(e=0)$.

6.2 Evaporation rate

Given enough time, the fraction of energetically-bound comets decays to zero. This is a consequence of the ergodic theorem, which states that one-dimensional Brownian motion will visit every point of space at least once, including the absorbing wall at $E=0$. We can use the energy distribution (90) to derive an analytical expression for the loss of comets as a function of time. As a first step, it is useful to calculate the fraction of bound comets as

[TABLE]

The escape, or “evaporation” rate is defined as the time-derivative of the fraction of comets that have escaped from the potential $\Phi_{s}$ in a time interval $t$, i.e. $f_{\rm esc}(t)\equiv 1-f_{b}(t)$. From (94)

[TABLE]

The right-hand integral in Equation (95) represents a flux of particles crossing the energy layer $E=0$ at the time $t$ which started from $\delta(E-E_{0})\delta(L_{x,0})\delta(L_{y,0})\delta(L_{z,0}-L_{c})$ at $t_{0}=0$. One can thus use the energy term of the diffusion equation (49) and the propagator (56) to measure the flux through the barrier $E=0$ as

[TABLE]

where we have used the relation between static and time-dependent coefficients $\tilde{C}_{E}=C_{E}t$ and $\tilde{D}_{E}=D_{E}t$ (see §4.2).

Combination of (95) and (96) shows that the escape rate corresponds to the average flux of probability crossing $E=0$ at the time $t$ which originates from the energy interval $(E_{1},E_{2})$ at $t_{0}=0$

[TABLE]

One can show that the rate at which comets leave a Keplerian potential approaches a power-law function ${\rm d}f_{\rm esc}/{\rm d}t\sim t^{-2/3}$ on time-scales $t_{\rm esc}(E_{2})\ll t\ll t_{\rm esc}(E_{1})$, where $t_{\rm esc}(E_{0})=E_{0}/C_{E}(E_{0},L_{c},t_{0})$ corresponds to the average unbinding time of particles moving on circular orbits with an energy $E_{0}$ at $t=t_{0}$. To demonstrate this property, it is useful to introduce a dimension-less variable $\beta=t_{\rm esc}/t$. As shown in §4.3, the impulsive approximation is accurate on time-scales $t\sim t_{\rm esc}$, hence one can use the analytical coefficients (69) to solve Equation (97), which yields

[TABLE]

and

[TABLE]

After some algebra, Equation (97) can be written as

[TABLE]

with

[TABLE]

On a time-scale $\beta_{2}\ll\beta\ll\beta_{1}$ Equation (99) becomes approximately constant, $K(E_{1},E_{2},t)\simeq K$. To calculate the value of $K$ one must replace the lower and upper limits of the integral by zero and infinity, respectively, such that

[TABLE]

where $K_{n}(z)$ corresponds to the modified Bessel function555This function satisfies the equation $z^{2}K_{n}^{\prime\prime}+zK_{n}^{\prime}-(z^{2}+n^{2})K_{n}=0$ (e.g. Press et al. 1992). Hence, inserting (100) into (98) yields a escape rate that approaches asymptotically a power-law

[TABLE]

Fig. 8 plots the escape rate measured from Monte-Carlo runs (symbols) and the analytical values derived from diffusion equations (97) (solid lines). Lines and symbols are colour-coded according to the energy of the outer-most circular orbit at $t_{0}=0$, $E_{2}=-Gm/(2R_{2})$. This Figure shows several noteworthy results. Note first that the rate at which particles leave the potential roughly peaks at the escape time of the outer-most radius, $t_{\rm esc}(E_{2})$, which is marked with vertical dotted lines for reference. Within short time intervals, $t\ll t_{\rm esc}(E_{2})$ the flux of particles that become energetically unbound is negligible in all models, except for very extended discs ($R_{2}\gtrsim 10D$) at $t_{0}$. On longer time-scales, $t\gtrsim t_{\rm esc}(E_{2})$, the escape rate converges asymptotically to the power-law curve (101) (black-dashed line), which scales as $\mathcal{R}_{\rm esc}\sim t^{-2/3}$. Comparison between Monte-Carlo models and analytical curves shows that diffusion equations systematically underestimate the fraction of particles crossing $E=0$. This is particularly visible at peak time, $t\sim t_{\rm esc}(E_{2})$. The origin of this mismatch can be traced back to the accelerated unbinding of particles as they escape from the system (see left panels of Fig. 6), which cannot be properly described by diffusion or Fokker-Planck equations, as discussed in §6.1.3. The agreement between Monte-Carlo and diffusion models improves as the fraction of particles in the fringe decreases and the escape rate approaches $\mathcal{R}_{\rm esc}\to 0$.

6.3 Production of retrograde orbits

As demonstrated in Appendix B, stochastic tidal fluctuations tend to isotropize the initial distribution of angular momentum. This implies a non-zero probability that at any given time a tracer particle will scatter onto a retrograde motion. Here we use Equation (92) to calculate the fraction of comets whose rotational direction is reversed as a function of time. Recall that by convection $L_{z,0}>0$ at $t=t_{0}$, hence the fraction of comets on retrograde orbits can be written as

[TABLE]

By definition, the rate at which retrograde orbits are produced corresponds to the time-derivative of (102)

[TABLE]

where $j_{E}(0,t|E_{0},L_{c},t_{0})$ is the flux of particles crossing the energy boundary $E=0$, Equation (96), and

[TABLE]

is the flux of probability crossing $L_{z}=0$ at the time $t$ which originated from $\delta(E-E_{0})\delta(L_{R})\delta(L_{z}-L_{c})$ at the time $t_{0}=0$. With these definitions in place, it is straightforward to derive a physical interpretation of the two terms appearing in Equation (103). In particular, the left-hand term has a negative value and accounts for the loss of retrograde orbits as they cross the absorbing barrier at $E=0$, while the right-hand one is positive and describes a flux of particles that diffuse into the retrograde interval $L_{z}<0$ at a fixed bound fraction.

One can follow similar steps as in Section 6.2 to find analytical solutions for ${\rm d}f_{\rm ret}/{\rm d}t$ on time-scales $t_{\rm esc}(E_{2})\ll t\ll t_{\rm esc}(E_{1})$. Let us first re-write Equation (103) as

[TABLE]

where $\mathcal{R}_{\rm esc}$ is the escape rate defined by Equation (97), and

[TABLE]

After some algebra, one can show that the series in Equation (105) converges to a power-law function

[TABLE]

Inserting (107) into (105) leads to a scale-free rate

[TABLE]

where $\mathcal{R}_{\rm esc}(t)$ is the escape rate is given by (101). Thus, we find that on time-scales $t_{\rm esc}(E_{2})\ll t\ll t_{\rm esc}(E_{1})$ the fraction of bound comets scattered onto retrograde orbits per unit of time is approximately $\sim 3\%$ of those escaping from the central potential.

Upper panel of Fig. 9 shows the rate of production of retrograde orbits ($\mathcal{R}_{\rm ret}$) as a function of time for the models plotted in Fig. 8. There is an overall good agreement between Monte-Carlo simulations and the analytical expression (103) (thin-colour lines), which improves as the time interval increases. As in Fig. 8, we find that the maximum rate $\mathcal{R}_{\rm ret}(t)$ approximately occurs at the average escape time of the initial outer-most disc radius, $t_{\rm esc}(E_{2})$ (marked with vertical dotted lines), which is similar to the isotropization time-scale (vertical dashed lines), as expected from (74). On longer time intervals, $t\gg t_{\rm esc}(E_{2})$, the scale rate approaches the power-law behaviour $\mathcal{R}_{\rm ret}(t)\sim t^{-2/3}$ given by Equation (108). The ratio $\mathcal{R}_{\rm ret}/\mathcal{R}_{\rm esc}$ is plotted in the lower panel of Fig. 9. We find that the rate of comets attaining a retrograde motion reaches a maximum $\mathcal{R}_{\rm ret}/\mathcal{R}_{\rm esc}\sim 0.1$ shortly before the escape rate peaks. Thereafter, the rate ratio slowly converges towards the value predicted by Equation (108), i.e. $\mathcal{R}_{\rm ret}/\mathcal{R}_{\rm esc}\sim 0.03$ (marked with a horizontal dashed line).

7 Discussion: Tidal heating by dark microhaloes

Previous Sections assume that substructures have a fixed mass $(M)$ and scale-radius $(c)$. Here, we extend our statistical formalism to ensembles of objects that cover a wide range of masses and sizes. For illustration, we discuss the effect of dark matter subhaloes on weakly-bound objects in the Solar System, as this provides an interesting, and rather special case of study.

Following Paper I (see §4.1 for details), let us consider statistical ensembles of dark matter subhaloes whose mass function and number density profile match those found in the Aquarius simulations of Milky Way-sized haloes

[TABLE]

with $B_{0}=2.02\times 10^{-13}M_{\odot}^{-1}\,{\rm kpc}^{-3}$ and $M_{0}=2.52\times 10^{7}M_{\odot}$. The function $g(r)=\exp\{-(2/\gamma)[(r/r_{-2})^{\gamma}-1]\}$ is an Einasto profile with parameters $\gamma=0.678$, $r_{-2}=199\,{\rm kpc}$ and $\alpha=-1.9$ (Springel et al. 2008; see also Han et al. 2016, Erkal et al. 2016). The subhalo mass function is believed to be truncated from below on scales comparable to the free-streaming length of the DM-particle candidates. For DM made of WIMPs with a mass $\sim 1-1000$ GeV the cut-off of the mass spectrum lies on sub-solar mass scales, $M_{1}/M_{\odot}\sim 10^{-12}$–$10^{-3}$ (e.g. Schmid et al. 1999; Hofmann et al. 2001; Green et al. 2005; Loeb & Zaldarriaga 2005; Diemand et al. 2005). As a result, the total number of galactic substructures predicted by Equation (109) is extremely large

[TABLE]

and diverges in the perfect-fluid limit $M_{1}\to 0$.

The velocity distribution of DM substructures in the Solar neighbourhood is assumed to be Maxwellian, Equation (36), with a dispersion $\sigma$ that can be computed from the isotropic Jeans equations as

[TABLE]

here $\Phi_{g}(r)$ is the potential of the Aquarius halo, which can be described by a Navarro, Frenk & White (1997) model with a virial mass and radius $M_{200}=1.84\times 10^{12}M_{\odot}$, $r_{200}=246\,{\rm kpc}$, respectively, and a concentration $c_{200}=16.11$ (Springel et al. 2008). Inserting the number density profile (109) into (111) yields $\sigma(r_{\odot})=296\,{\rm km\,s^{-1}}$. The mean relative velocity between the Sun and the dark matter microhaloes in Equation (36) corresponds to the circular velocity of the Milky Way at the solar radius $V_{\odot}=V(r_{\odot})=230\,{\rm km\,s^{-1}}$ at $r_{\odot}\approx 8\,{\rm kpc}$ (Eilers et al. 2018). Hence, we find that the averaged velocity (38) is

[TABLE]

As the cutoff of the mass spectrum is shifted to lower values, the average separation between subhaloes decreases. From Equation (109) we find

[TABLE]

whereas the characteristic duration of tidal fluctuations shortens in proportion to $D(r_{\odot})$. Combination of (30), (112) and (113) yields

[TABLE]

which is comparable to the orbital period of comets with a semi-major axis

[TABLE]

Given that the mean-life of a tidal fluctuation (114) is comparable to the time-span between subsequent fluctuations, comets with a semi-major axis $a\gg a_{\rm ch}$ will experience multiple ‘encounters’ with dark microhaloes during a single orbital revolution. To illustrate this point, Fig. 10 plots the average number of fluctuations per orbital period, $N_{\rm enc}=P/T_{\rm ch}(r_{\odot})$, as a function of the minimum subhalo mass, $M_{1}$. The upper-limit of the subhalo mass function is set at $M_{2}=10^{11}M_{\odot}$, although the results shown in Fig. 10 are largely independent of this choice insofar as $M_{1}\ll M_{2}$. As expected, orbits with a semi-major axis comparable or smaller to that of Pluto ($a\sim 40\,{\rm AU}$, blue solid line) must perform several orbital revolutions before experiencing a single fluctuation of the local tidal field, and therefore responds adiabatically to the subhalo background (see §3.3). In contrast, comets in the Oort cloud ($a\sim 10^{5}\,{\rm AU}$, red solid line) feel a rapidly fluctuating tidal field, and will react impulsively to gravitational interactions with dark matter substructures.

The magnitude of tidal fluctuations depends very strongly on the relation between the mass and size of dark matter substructures. Following Paper I, it is useful to quantify this dependency by introducing a power-law size function

[TABLE]

For the subhaloes found in the Via Lactea II simulation at redshift $z=0$ (Diemand et al. 2008) the best-fitting parameters are $c_{0}\simeq 0.53\,{\rm kpc}$ and $\beta\simeq 0.5$ (Erkal et al. 2016). The average size of subhaloes in a mass bin $M\in(M_{1},M_{2})$ derived from (116) is

[TABLE]

which is grossly consistent with the size of microhaloes expected from an inflation-produced primeval fluctuation spectrum (Berezinsky et al. 2003), as well as with simulations that resolve the size-mass relation of the first structures formed in the Universe (Ishiyama 2014; see also Sánchez-Conde & Prada 2014). Fig. 10 shows that the gap between the mean size of dark matter substructures (black long-dashed line) and their average separation (black short-dashed line) widens as the value of $M_{1}$ decreases. One can use (117) and (113) to compute the size-to-separation ratio, which goes as $\langle c\rangle/D(r_{\odot})\sim 4\times 10^{-4}[M_{1}/(10^{-6}M_{\odot})]^{0.2}$, thus approaching the point-mass behaviour $\langle c\rangle/D(r_{\odot})\to 0$ in the perfect-fluid limit $M_{1}\to 0$.

For dark matter microhaloes with power-law mass and size functions Equations (109) and (116) can be combined as (see Paper I)

[TABLE]

Note that (118) assumes an exact correspondence between size and mass (i.e. no scatter).

Over a time interval $t$, the average velocity increments experienced by comets located at a fixed distance from the Sun $R$ can be calculated as $\langle|\Delta{\textbf{{V}}}|^{2}\rangle=tR^{2}\langle\Lambda^{2}TA\rangle\delta_{s}(t)$, where $\langle\Lambda^{2}TA\rangle$ is given by Equation (41), and $\delta_{s}(t)=1-\exp(-t/\tau_{s})$ is the sampling-delay function defined in §5.2. Since we are mainly interested in the response of comets in the outer-most regions of the Solar system ($a\gg a_{\rm ch}$), one can safely assume that tidal perturbations occur in an impulsive regime, wherein the adiabatic correction (39) is $A\approx 1$. To average both quantities, $\langle\Lambda^{2}T\rangle$ and $\tau_{s}$, over substructures with a known distribution of masses and sizes one can simply write $n\to\int\int{\rm d}M{\rm d}c\frac{{\rm d}^{2}n}{{\rm d}M{\rm d}c}$. Equation (34) then becomes

[TABLE]

This Equation exhibits two well-defined behaviours

• •

$3+\alpha-2\beta\leq 0$. This regime arises when the size function is sufficiently steep, $\beta\geq(3+\alpha)/2$, leading to an average velocity increment that diverges in the limit $M_{1}\to 0$.

• •

$3+\alpha-2\beta>0$. This case corresponds to relatively shallow size functions, $\beta<(3+\alpha)/2$, leading to an asymptotically convergent value of $\langle\Lambda^{2}T\rangle$ as $M_{1}\to 0$.

The second parameter of relevance is the sampling time-scale, $\tau_{s}$, which measures the average time interval required to sample the large-force tail of probability function $p(\Lambda)$. The ratio $\tau_{s}/T_{\rm ch}$ informs on the average number of fluctuations that must occur before the tidal force variance $\langle\Lambda^{2}\rangle$ converges to the analytical value (16). As shown in Fig. 3, the magnitude of tidal fluctuations is exponentially suppressed on time-scales $t\ll\tau_{s}$ (i.e. $\delta_{s}\ll 1$). Here, it is convenient to average the time reciprocal of the sampling time-scale, i.e. the sampling frequency $w_{s}=1/\tau_{s}$. From Equation (84) we find

[TABLE]

Again, this expression reveals two possible regimes

• •

$1+\alpha+2\beta\leq 0$. In this case the sampling frequency diverges in the perfect fluid limit $M_{1}\to 0$, which implies a progressively faster convergence of $\langle\Lambda^{2}T\rangle$ towards the analytical value (119) as the value of $M_{1}$ decreases. As expected from Figs. 3 and 4, this regime arises in substructure populations with a relatively shallow mass function $\beta\leq-(1+\alpha)/2$.

• •

$1+\alpha+2\beta>0$. This regime is associated with a relatively steep size function, $\beta>-(1+\alpha)/2$, which leads to a convergent sampling frequency in the limit $M_{1}\to 0$.

It is straightforward to show that the solution to the equations $3+\alpha-2\beta=0$ and $1+\alpha+2\beta=0$ is $(\alpha,\beta)=(-2,1/2)$, which corresponds to an ensemble of substructures with a constant mass fraction per logarithmic bin, i.e. $M(\Delta n/\Delta\ln M)\sim{\rm const.}$, and a squared mass-size relation $M\sim c^{2}$. For this particular combination, both quantities $\langle\Lambda^{2}T\rangle$ and $\langle w_{s}\rangle$ show a mild logarithmic divergence $\sim\ln(M_{2}/M_{1})$ as $M_{1}\to 0$. Strikingly, Fig. 11 shows that substructures in the Via Lactea II simulation sit closely to this regime. By a narrow margin, however, the best-fitting indices, $(\alpha,\beta)\simeq(-1.9,0.5)$, meet the condition for slowly-converging values of velocity increments suggested by (119), i.e. $\beta<(3+\alpha)/2\simeq 0.55$, as well as the condition for a slowly-converging sampling frequency (and thus a finite time-scale $\tau_{s}$) given by (120), $\beta>-(1+\alpha)/2\simeq 0.45$.

The above conclusion must be taken with caution, as several issues must be addressed before CDM can make a firm prediction on the convergence of $\langle|\Delta{\textbf{{V}}}|^{2}\rangle$ and $\langle w_{s}\rangle$ at low subhalo masses. For example, our statistical ensembles rely on extrapolations of cosmological $N$-body simulations that do not include baryons and have particle-mass resolutions that lie many orders of magnitude above the mass-scale of dark matter microhaloes. Indeed, the largest CDM simulations of Milky Way-like haloes to date, i.e. Aquarius (Springel et al. 2008), Via Lactea II (Diemand et al. 2007) and GHALO (Stadel et al. 2009), have particle-masses $\sim 10^{4}$–$10^{5}M_{\odot}$, i.e. at least $\sim 10$ orders of magnitude above the WIMP mass-scale. On the other hand, simulations of structure formation that do resolve the free-streaming length of WIMPs (e.g. Ishiyama et al. 2014) typically stop at redshift $z\sim 30$ owing to their large computational cost. To the resolution limitations one must add the formidable numerical challenge of following the dynamical evolution of subhaloes that span $\sim 18$ orders of magnitude in mass to redshift $z=0$ (van den Bosch 2017, van den Bosch & Ogiya 2018).

Semi-analytical algorithms offer a flexible alternative approach to reach small low subhalo masses at a relatively low computational cost (e.g. Peñarrubia et al. 2010; Stref & Lavalle 2017; Hiroshima et al. 2018). However, these methods are limited by our incomplete theoretical understanding of the processes of tidal stripping and mass loss (e.g. Daniel et al. 2017 and references therein).

Although often neglected, baryons may play an important role in shaping the mass & size functions of low-mass subhaloes in the Solar neighbourhood. E.g., fly-by encounters with individual stars (Green & Goodwin 2007) and the smooth tidal field of the Milky Way disc (D’Onghia et al. 2010; Peñarrubia et al. 2010; Errani et al. 2017) will enhance mass loss and modify the size of these objects. Owing to their low-binding energies, micro-haloes experience a wide variety of mass loss histories depending on their individual orbits and accretion times. Combination of the above effects will necessarily magnify the scatter of the mass-size relation of low-mass subhaloes in the Solar neighbourhood. Note that it is not straightforward to predict the effect of scatter on the behaviour of $\langle\Lambda^{2}T\rangle$ and $\langle w_{s}\rangle$ because neither of these quantities scale linearly with subhalo size, $c$. Finally, the suppression of structure formation at low masses leads to a relatively shallow subhalo mass function (Schneider et al. 2012; Angulo et al. 2013; Lovell et al. 2014), which may also affect the magnitude of tidal heating. Addressing these issues goes beyond the main goals of this work and will be studied in a separate contribution.

8 Summary

This paper uses stochastic calculus techniques to describe the dynamical response of tracer particles to stochastic variations of the combined tidal field generated by a large population of extended objects. Over a sufficiently-long interval of time the cumulative effect of tidal fluctuations leads to a random walk of orbital velocities known as “tidal heating”. Under the assumption that velocity impulses are small, $|\Delta{\textbf{{V}}}|<<|{\textbf{{V}}}|$, this paper uses diffusion equations to describe the response of tracer particles to random variations of a gravitational field (e.g. Chandrasekhar 1943; Kandrup 1980), which reduces the problem of tidal heating to the computation of diffusion coefficients, $\langle\Delta{\textbf{{V}}}\rangle$ and $\langle|\Delta{\textbf{{V}}}|^{2}\rangle$, where brackets denote averages over the spectrum of tidal fluctuations, p({\mbox{\boldmath\Lambda}}) (see §3.2). This approach solves a number of shortcomings in Chandrasekhar’s (1941a) classical theory, in which tidal heating is treated as the cumulative effect of isolated encounters between a test particle with an infinite background of individual point-mass particles:

• •

Our theory does not require arbitrary truncations of the force spectrum either at large or weak forces.

• •

External forces are generated by extended objects (not particles), which allows us to combine tidal heating from baryonic and dark matter substructures covering a wide range of masses and sizes. On the other hand, the theory breaks down when applied to objects with divergent forces (e.g. black holes).

• •

Modelling the combined tidal force generated by background substructures avoids a detailed analysis of three-body encounters.

• •

Implementation of Weinberg’s adiabatic corrections leads to diffusion coefficients that describe perturbations in impulsive as well as adiabatic regimes.

• •

Our analysis accounts for the period of time required to fully sample the spectrum of tidal fluctuations generated by a finite number of substructures.

In Section 4, we treat stochastic tidal heating of self-gravitating systems as a diffusion process in a confined region of the 4-dimensional integral-of-motion space, where the energy is definite negative $E<0$, and the angular momentum components lie in the range $-L_{c}(E)\leq L_{i}\leq L_{c}(E)$, with $L_{c}(E)$ being the angular momentum of a circular orbit with fixed energy. To this aim, Appendix B derives Green’s functions with boundary conditions that place an absorbing barrier at $E=0$, and reflecting surfaces at $L_{i}=\pm L_{c}(E)$. The coefficients of the Green’s functions, $\langle\Delta E\rangle$, $\langle\Delta E^{2}\rangle$, $\langle\Delta{\textbf{{L}}}\rangle$ and $\langle|\Delta{\textbf{{L}}}|^{2}\rangle$, are calculated from averages over ensembles of tracer particles and substructures. §4.3 discusses the analytical case of a Keplerian potential in detail. Our results show that tidal heating leads to (i) a steady flow of probability drifting from bound energies towards $E\to 0$, which causes “tidal evaporation” as tracer particles gain sufficient energy as to escape from the system, and (ii) “isotropization”, as the initial angular momentum distribution is randomized.

Diffusion equations are tested in Section 5 with the aid of $N$-body experiments in which the the local tidal field is computed as a direct summation of the forces generated by a large ($N\gg 1$) population of Hernquist (1990) spheres in dynamical equilibrium. We find that (i) sampling the large-force tail of the distribution p(\mbox{\boldmath\Lambda}) requires of the order of $N_{\rm enc}\approx 4.8(D/c)^{2}$ fluctuations, where $D/c\gg 1$ is ratio between the mean separation and the size of individual substructures. This implies that our theory never attains statistical convergence when applied to a background of point-masses ($c\to 0$). (ii) Tracer particles in the “fringe” ($E\sim 0$) experience an rapid increase of energy and angular momentum that cannot be explained by our diffusion theory. In contrast, the process of tidal evaporation can be successfully described Monte-Carlo $N$-body simulations that sample random velocity kicks $\Delta{\textbf{{V}}}$ drawn from a probability function $\Psi({\textbf{{V}}},\Delta{\textbf{{V}}},\Delta t)$, Equation (31).

Section 6 follows the evolution of idealized planetary (mass-less) discs with an homogeneous energy distribution at $t_{0}=0$ subject to stochastic fluctuations of the local tidal field. The tidal evaporation rate ($\mathcal{R}_{\rm esc}$), and the rate of comets moving on retrograde orbits ($\mathcal{R}_{\rm ret}$) are computed analytically from the flux of comets going through the energy layer $E=0$ and the angular momentum boundary $L_{z}=0$, respectively. We find that on long time-scales both rates approach a scale-free relation $\mathcal{R}\sim t^{-2/3}$, and that in this regime the number of comets acquiring a retrograde motion corresponds to a small fraction of those escaping from the system, $\mathcal{R}_{\rm ret}/\mathcal{R}_{\rm esc}\simeq 0.03$. Overall, the energy and angular momentum distributions derived from Monte-Carlo simulations and Green’s convolutions are in excellent agreement except for particles in the fringe ($E\sim 0$), where the condition of small increments demanded by the diffusion theory, $|\Delta E/E|\ll 1$, cannot be guaranteed.

Section 7 inspects the effect of dark matter substructures on the dynamics of weakly-bound objects in the Solar system. To this end, we construct substructure ensembles that mimic the subhalo population found in cosmological $N$-body simulations of Milky Way-sized haloes. Extrapolation of the subhalo mass & size functions found in the Aquarius (Springel et al. 2008) and Via Lactea II (Diemand et al. 2007) models suggests that comets in the Oort cloud may be sensitive to the presence of dark subhaloes with sub-solar masses. In particular, for a standard WIMP free-streaming mass Equation (115) indicates that objects with a semi-major axis $a\gtrsim 150\,{\rm AU}$ may experience a large number of tidal interactions with planet-sized microhaloes on a dynamical time-scale. Interestingly, recent studies of the dynamics of Trans-Neptunian Objects with semi-major axes $a\gtrsim 250\,{\rm AU}$ find evidence for a non-equilibrium configuration, which may be caused by perturbations from a super-Earth ($M\gtrsim 10M_{\oplus}$) object known as “Planet Nine” (Batygin & Brown 2016; Becker et al. 2018). The possibility to put constraints on the dark matter clumpiness using objects in the Solar system deserves detailed examination in a separate contribution.

The diffusion and Monte-Carlo techniques presented here provide statistical tools to study the combined effect of baryonic (e.g. free-floating planets, stars, giant molecular clouds, etc) and dark matter substructures on the dynamical evolution of a large number of weakly-bound objects (e.g. comets, protoplanetary discs, wide-binary stars, stellar clusters and gas clouds) at a minor computational cost. Combining observational constraints from a large variety of visible objects represents one of the best hopes to probe the unknown behaviour of dark matter on subsolar mass scales, a question with important implications for particle physics and cosmology.

Appendix A The autocorrelation function

Here we calculate the autocorrelation function W({\mbox{\boldmath\Lambda}}_{0},{\mbox{\boldmath\Lambda}}_{t}). Our main assumptions are that (i) substructures are randomly distributed over the volume $V^{\prime}$, and (ii) they move on straight lines.

By definition, $W$ is the convolution of two probability densities

[TABLE]

where p(\mbox{\boldmath\Lambda}) is the probability of experiencing a tidal vector in the interval {\mbox{\boldmath\Lambda}},{\mbox{\boldmath\Lambda}}+{\rm d}{\mbox{\boldmath\Lambda}}. The computation of $W$ is more straighforward in Fourier space. Application of the convolution theorem to Equation (52) yields

[TABLE]

where the Fourier transform $\tilde{p}({\textbf{{k}}}_{0})$ can be written as (Holtsmark 1919; see Paper I)

[TABLE]

The derivation of $\tilde{p}({\textbf{{k}}}_{t})$ can be done in a similar way under the assumption that substructures move on straight lines with a relative velocity v. Hence, if ${\textbf{{r}}}_{i}$ is the initial position vector, the probability to find the particle at a later time $t$ at the location ${\textbf{{r}}}^{\prime}_{i}$ is $\delta({\textbf{{r}}}^{\prime}_{i}-{\textbf{{r}}}_{i}-{\textbf{{v}}}t)$. Thus,

[TABLE]

Note that the last equalities in (123) and (124) implicitly assume that the vectors {\mbox{\boldmath\lambda}}_{i} are spatially uncorrelated at all times.

The Fourier transform of the autocorrelation function (121) becomes

[TABLE]

The last integral can be re-written as

[TABLE]

which elevated to the $N$-th power becomes

[TABLE]

where $n\equiv N/V^{\prime}$ is the number density of substructures. Equation (125) thus becomes

[TABLE]

with

[TABLE]

For substructures with a velocity distribution $f({\textbf{{v}}})$, where $\int{\rm d}v^{3}\,f({\textbf{{v}}})=1$, Equation (127) can be generalized as (Chandrasekhar 1943)

[TABLE]

Finally, taking the inverse Fourier transform of Equation (125) we obtain

[TABLE]

thus recovering autocorrelation function derived in Chandrasekhar (1943) under the assumption that the effects of individual interactions can be described as a Markov chain process (see his §3).

Appendix B Diffusion in confined regions

Here we derive Green’s propagators in a confined region of the integral-of-motion space. As in previous Sections, we shall assume that diffusion in energy and angular momentum are statistically independent processes, such that solutions to the diffusion equation (49) can be assumed to have a separable form (50).

B.1 Energy space

Let us assume that particles that attain a positive gravitational energy will leave the system to never come back. This calls for the diffusion equation in energy space

[TABLE]

with an absorbing boundary at $E=0$, such that $p(E=0,t)=0$. To solve the differential equation we use the method of images originally introduced by Lord Kelvin (see Feller, 1971), where one places an image source (or sink) at $-E_{0}$ as mirror image of the original source at $+E_{0}$ with a strength or intensity selected to match the boundary condition at $E=0$. The initial conditions in this case can be written as

[TABLE]

where the dimensionless parameter $\eta$ determines the strength of the mirror source. The solution of (130) with boundaries (131) is a linear combination of two mirrored Gaussians

[TABLE]

with coefficients evaluated at $\tilde{C}=\tilde{C}(E_{0},L_{0},t-t_{0})$ and $\tilde{D}=\tilde{D}(E_{0},L_{0},t-t_{0})$. It straightfoward to show that the condition $p(0,t)=0$ implies

[TABLE]

Left panel of Fig. 12 plots the Green’s propagator (132) at three different snaphosts. Far from the boundary $E=0$, the right-hand side of Equation (132) has a negligible contribution, and the energy propagator can be approximated by a Gaussian distribution (black dots). However, as the time progresses, the Green’s function broadens and shifts systematically towards lower energies. By construction, $p(E,t|E_{0},{\textbf{{L}}}_{0},t)\to 0$ as $E\to 0$, which strongly deviates from the behaviour of the free-diffusing Gaussian solution in an infinite volume.

B.2 Angular momentum space

The angular momentum of gravitationally-bound particles cannot exceed the maximum set by circular orbits. Hence, $L\leq L_{c}(E)$, where $L_{c}(E)$ is the angular moment of a circular orbit with energy $E$. E.g., for a Keplerian potential (68) it follows from (70) that the maximum angular momentum is $L_{c}(E)=\sqrt{Gma}$, where $a=Gm/(-2E)$ is the semi-major axis of orbits with energy $E$.

Our goal is to solve the diffusion equation

[TABLE]

with reflecting boundaries $\nabla_{L}p=0$ at $L=0,L_{c}$. Given that diffusion in angular momentum space is assumed to be isotropic, one may attempt to solve (134) in spherical coordinates. However, the probability function corresponds to an infinite sum of Bessel functions and spherical harmonics (Carslaw & Jaeger 1986), which complicates an intuitive interpretation of the result. Instead, here we will solve (134) in Cartesian coordinates, where the differential operator $\nabla_{L}^{2}=\partial^{2}/\partial L_{x}^{2}+\partial^{2}/\partial L_{y}^{2}+\partial^{2}/\partial L_{z}^{2}$. The problem at hand, therefore, is that of isotropic diffusion in a cubic box with reflecting boundaries. To gain insight, first we solve the one-dimensional case, and subsequently provide the general 3D solution.

B.2.1 One dimension

Let us re-write the one-dimensional diffusion equation in angular momentum (134) using a simplified notation

[TABLE]

where $x=L_{x}$ is defined within the interval $x\in[-L_{c},L_{c}]$, and $D=D_{L}(E_{0},L_{0},t_{0})$. To limit the region in which particles can diffuse, we set reflecting boundaries at $x=\pm L_{c}$, such that $\partial p/\partial x|_{\pm L_{c}}=0$. Equation (135) admits separable solutions

[TABLE]

Inserting (136) into (135) and re-arranging yields

[TABLE]

where $\lambda$ is a separation constant. Since the coefficient $D$ does not contain an explicit dependence on time (see §4.2), the function $T(t)$ can be directly integrated from (137), which yields

[TABLE]

Let us now define the variable $x^{\prime}=x+L_{c}$ and expand $X(x^{\prime})$ as a Fourier series with a fundamental period $4L_{c}$

[TABLE]

with a wave number $k_{n}=n\pi/(2L_{c})$. It is clear from (139) that reflecting boundaries ${\rm d}X/{\rm d}x^{\prime}|_{0,2L_{c}}=0$ require $b_{n}=0$ for integers $n\geq 1$. Note also that although the Fourier series (139) extends to $[-2L_{c},0]$, here we are only interested in the solution inside the interval $[0,2L_{c}]$. Inserting (139) into (137) yields a separation constant

[TABLE]

To find the coefficients of the Fourier series we must use the initial condition $p=\delta(x-x_{0})=\delta(x^{\prime}-x_{0}-L_{c})$ at $t=0$. To this aim, we multiply both sides of (139) by $\cos(k_{m}x^{\prime})$ and integrate from 0 to $2L_{c}$

[TABLE]

Using the orthogonality relations

[TABLE]

we find

[TABLE]

Combination of (138) and (139) then yields

[TABLE]

Compare Equation (144) with the Fourier expansion of a Gaussian function

[TABLE]

within the interval $x\in[-L_{c},L_{c}]$. Following the same steps as above, we introduce a variable $x^{\prime}=x+L_{c}$, and choose a fundamental period $4L_{c}$, neglecting the even terms of the series $x^{\prime}\in[-2L_{c},0]$

[TABLE]

where $y=x^{\prime}-x_{0}-L_{c}$. If the mean of the Gaussian function is located far from the boundaries, $|L_{c}-x_{0}|/\sigma\gg 1$, one can safely shift the limits of the integral (145) to infinity. Hence, Equation (145) becomes

[TABLE]

which leads to a Fourier series

[TABLE]

It is clear that Equation (147) recovers (144) if

[TABLE]

Applying the separation constant (140) yields the celebrated behaviour of the free-diffusing Gaussian propagators (52)

[TABLE]

thus demonstrating that far from the boundaries of the confined region particles diffuse as if they were moving in an infinite medium, $p(x,t)\approx p_{G}(x,t)$.

Right panel of Fig. 12 shows the solution (144) at three different snapshots. As expected from the above results, the confined Green’s function (144) has a Gaussian form (black dots) on short time scales $\pi^{2}/(2L_{c})^{2}Dt\ll 1$. However, as the time progresses all modes with $n>0$ decay exponentially, and the propagator converges asymptotically towards a constant value $p(x,\infty)=1/(2L_{c})$. Hence, on long time-scales $\pi^{2}/(2L_{c})^{2}Dt\gg 1$ the angular momentum distribution progressively becomes isotropic. This behaviour is in stark contrast with the free-diffusing solution in an infinite domain, which vanishes $p_{G}\to 0$ in the limit $t\to\infty$.

B.2.2 Three dimensions

The above derivation can be generalized to a three-dimensional space by replacing the diffusion equation (135) with

[TABLE]

and (136) with

[TABLE]

Inserting (150) into (148) and re-arranging yields

[TABLE]

with a separation constant

[TABLE]

Following similar steps as in the 1D case, we find that the separable solution (150) with reflecting boundaries at $x=y=z=\pm L_{c}$ can be expressed as an infinite cosine series

[TABLE]

with $x^{\prime}=x+L_{c}$, $y^{\prime}=y+L_{c}$ and $z^{\prime}=z+L_{c}$. Application of the initial conditions $p(x,y,z,t=0)=\delta(x-x_{0})\delta(y-y_{0})\delta(z-z_{0})$ and the orthogonality relations (142) yields

[TABLE]

with $\alpha_{000}=1/8$, $\alpha_{n00}=\alpha_{0m0}=\alpha_{00l}=1/4$, $\alpha_{nm0}=\alpha_{n0l}=\alpha_{0ml}=1/2$; and $\alpha_{nml}=1$  for $n,m,l\geq 1$.

Appendix C Phase-space averages

Here we briefly discuss how to compute phase-space averages of particle ensembles with a fixed energy and angular momentum (for further details see Appendix A of P15). As a first step, let us define the tangential and radial velocity components as $v_{t}$ and $v_{r}$, respectively, such that the angular momentum becomes $L=rv_{t}$, the velocity volume ${\rm d}^{3}v=2\pi v_{t}{\rm d}v_{t}{\rm d}v_{r}=2\pi{\rm d}v_{r}L{\rm d}L/r^{2}$, and the energy $E=v^{2}/2+\Phi(r)=v_{r}^{2}/2+\Phi(r)+L^{2}/(2r^{2})$, hence at fixed radius ${\rm d}E=v{\rm d}v$. The volume accessible to tracer particles with a particular combination of $E$ and $L$ is a spherical shell with inner (peri-centre) and outer (apo-centre) radii, $R_{p}$ and $R_{a}$, respectively. At these two radii the radial velocity component vanishes, $v_{r}=0$, and the kinetic energy is entirely in the tangential direction $v=v_{t}$. Note that at each radius $v_{r}$ can be both positive and negative. In what follows we take only positive values of $v_{r}$ and double the number of particles with this velocity.

In order to express the probability function $N(E,L,t)$ as a function of phase-space coordinates, one needs to map points in the integral-of-motion space onto phase space through the equation

[TABLE]

where $f=N/\omega$ is the probability of finding a particle with integrals $(E,L)$ in the phase-space volume ${\rm d}^{6}\Omega$ centred at the coordinates $({\textbf{{r}}},{\textbf{{v}}})$ at the time $t$, and $[\omega]\equiv\partial({\textbf{{r}}},{\textbf{{v}}})/\partial(E,L)$ is the matrix of density of states. The Jacobian of this matrix, $\omega(E,L,t)$, defines the maximum phase-space volume that particles with a given combination of energy and angular momentum can access. The simplest derivation of $\omega$ is found by integrating the volume ${\rm d}^{6}\Omega={\rm d}^{3}{\textbf{{r}}}{\rm d}^{3}{\textbf{{v}}}=4\pi r^{2}{\rm d}r\times 2\pi v_{t}{\rm d}v_{t}{\rm d}v_{r}$ over the accessible spherical shell, which yields (Spitzer 1987)

[TABLE]

with a density of states

[TABLE]

Here $P(E,L)$ denotes to the orbital period

[TABLE]

with $R_{p}$ and $R_{a}$ being the peri- and apo-centres of the orbit, respectively, which correspond to the radii where the radial velocity vanishes, $v_{r}=\{2[E-\Phi_{s}(R)-L^{2}/(2R^{2})]\}^{1/2}=0$.

Tracer particles with a given combination of energy and angular momentum follow a delta function in the integral-of-motion space, $N(E,L)=\delta(E-H)\delta(L-rv_{t})$. The distribution function $f=N/\omega$ can be used to compute the average of a generic function $X(r,v)$ in phase-space. Combination of (155), (156) and (157) yields

[TABLE]

which weights a region of the orbit ${\rm d}r$ by the fraction of the orbital period that a particle spends in that region, $(2{\rm d}r/v_{r})/P$. Hence, Equation (158) establishes an ergodic correspondence between ensemble and time averages. It is important to bear in mind that the ergodic property only holds for particle distributions that are fully mixed in phase space (see P15).

The averaged quantities appearing in Equation (55) can be calculated analytically using (158) as

[TABLE]

for $m=2,4$, and

[TABLE]

which in general must be integrated numerically. An exception is the Keplerian potential (68), which we inspect in some detail below.

C.1 Keplerian potential

As a first step, it is useful to express the orbital energy and angular momentum in terms of a semi-major axis ($a$) and eccentricity ($e$) using Equation (70). The peri- and apocentres of the orbit correspond to the radii where $v_{r}=\{2[E-\Phi_{s}(R)-L^{2}/(2R^{2})]\}^{1/2}=\{2[-Gm/(2a)+Gm/R-Gm\,a(1-e^{2})/(2R^{2})]\}^{1/2}=0$, which admits two solutions

[TABLE]

The orbital period can be derived analytically from (157) by changing the integration variable to $\nu=r/a$ and using the roots of the radial velocity (161), which yields

[TABLE]

hence, the orbital period is independent of the orbital angular momentum.

Similarly, the averaged second and fourth power of the radius (159) can be written as

[TABLE]

and

[TABLE]

The averaged squared velocity is

[TABLE]

while the integral (160) has an analytical solution

[TABLE]

Note that $\overline{R^{2}V^{2}}\neq\overline{R^{2}}\times\overline{V^{2}}$ for $e\neq 0$ (non-circular orbits).