¹¹institutetext: Centre for Theoretical Physics, The British University in Egypt, Sherouk City 11837, Cairo, Egypt

Gravitational Brownian motion as inhomogeneous diffusion:
Black hole populations in globular clusters

Zacharias Roupas Gravitational Brownian motion as inhomogeneous diffusion: Black hole populations in globular clustersGravitational Brownian motion as inhomogeneous diffusion: Black hole populations in globular clusters

Recent theoretical and numerical developments supported by observational evidence strongly suggest that many globular clusters host a black hole (BH) population in their centers. This stands in contrast to the prior long-standing belief that a BH subcluster would evaporate after undergoing core collapse and decoupling from the cluster. In this work, we propose that the inhomogeneous Brownian motion generated by fluctuations of the stellar gravitational field may act as a mechanism adding a stabilizing pressure to a BH population. We argue that the diffusion equation for Brownian motion in an inhomogeneous medium with spatially varying diffusion coefficient and temperature, which was first discovered by Van Kampen, also applies to self-gravitating systems. Applying the stationary phase space probability distribution to a single BH immersed in a Plummer globular cluster, we infer that it may wander as far as $\sim 0.05,\,0.1,\,0.5{\rm pc}$ for a mass of $m_{\rm b}\sim 10^{3},\,10^{2},\,10{\rm M}_{\odot}$ , respectively. Furthermore, we find that the fluctuations of a fixed stellar mean gravitational field are sufficient to stabilize a BH population above the Spitzer instability threshold. Nevertheless, we identify an instability whose onset depends on the Spitzer parameter, $S=(M_{\rm b}/M_{\star})(m_{\rm b}/m_{\star})^{3/2},$ and parameter $B=\rho_{\rm b}(0)(4\pi r_{c}^{3}/M_{b})(m_{\star}/m_{\rm b})^{3/2}$ , where $\rho_{\rm b}(0)$ is the Brownian population central density. For a Plummer sphere, the instability occurs at $(B,S)=(140,0.25)$ . For $B>140,$ we get very cuspy BH subcluster profiles that are unstable with regard to the support of fluctuations alone. For $S>0.25,$ there is no evidence of any stationary states for the BH population based on the inhomogeneous diffusion equation.

Key Words.:

stellar clusters, gravitational kinetic theory, black holes

1 Introduction

The study of Brownian motion was introduced in astrophysics by Chandrasekhar (1943b), who used it to establish, among other things, the concept of dynamical friction (Chandrasekhar 1943a). Brownian motion has also been used to model or estimate the motion of a massive black hole in the center of a stellar cluster or galaxy (Chatterjee et al. 2002b, a, 2003; Merritt 2005; Merritt et al. 2007; Merritt 2013; Bortolas et al. 2016; Lingam 2018; Di Cintio et al. 2020). While significant developments have been achieved regarding the statistical mechanics and a general kinetic theory of orbit-averaged motion in action-angle variables (Tremaine & Weinberg 1984; Binney & Lacey 1988; Chavanis 2012, 2013; Vasiliev 2017; Roupas et al. 2017; Fouvry & Bar-Or 2018; Tremaine 2020) or other canonical variables (Roupas 2020), a diffusion equation in configuration space for inhomogeneous Brownian motion – including varying velocity dispersion and diffusion coefficient– has not been applied as of yet. The paradigm of inhomogeneous Brownian motion is directly relevant to a self-gravitating system consisting of a population of heavier, fewer bodies immersed in a bigger, highly populated cluster of lighter bodies, whose profile is typically inhomogeneous regarding the density, velocity dispersion, and diffusion coefficient.

Such a population of a significant number of black holes (BHs) is believed to exist in the center of many globular clusters. Its existence is supported by numerical and theoretical developments (Merritt et al. 2004; Mackey et al. 2008; Morscher et al. 2013; Breen & Heggie 2013; Morscher et al. 2015; Wang et al. 2016; Arca-Sedda 2016; Rodriguez et al. 2016; Chatterjee et al. 2017; Kremer et al. 2018; Askar et al. 2018; Arca Sedda et al. 2018; Weatherford et al. 2018, 2019; Kremer et al. 2019) as well as observational evidence (Maccarone et al. 2007; Barnard et al. 2008; Strader et al. 2012; Irwin et al. 2010; Roberts et al. 2012; Chomiuk et al. 2013; Miller-Jones et al. 2015; Taylor et al. 2015; Minniti et al. 2015; Bahramian et al. 2017; Giesers et al. 2018; Shishkovsky et al. 2018; Abbate et al. 2019). These recent advances stand in contrast to a long-held belief that globular clusters cannot retain their BHs (Spitzer 1969; Kulkarni et al. 1993; Sigurdsson & Hernquist 1993). Spitzer instability (Spitzer 1969) seems unavoidable for massive stars which are subject to mass segregation, decoupling from the cluster and undergoing gravothermal collapse. The resulting black hole subcluster, according to a traditional view, would itself undergo core collapse and become so dense that would evaporate due to two-body and three-body encounters, apart from one or two BHs that may be retained in the centre (Kulkarni et al. 1993; Sigurdsson & Hernquist 1993). However, examining this picture more carefully, we may realize that since Spitzer instability depends on total mass and the BH-population has significantly less total mass than the progenitor subcluster of massive stars, the former might not undergo gravothermal collapse and might, instead, become stabilized. It has also been suggested that the BH-population does not stay decoupled from the cluster for as long a time as initially thought (Breen & Heggie 2013; Morscher et al. 2013).

Here, we propose an additional theoretical component to this picture. We investigate whether a BH-population that is well within the Spitzer instability regime may be supported by Brownian motion induced by the fluctuations of the cluster’s gravitational field. To this end, we describe the Brownian motion of the BHs due to fluctuations of the field by an inhomogenous diffusion equation that takes into account the density, velocity dispersion, and dynamical friction coefficient spatial variations. This model for the inhomogeneous diffusion of Brownian particles was discovered by van Kampen (1988) in a general setting and we argue that it also applies to self-gravitating systems.

In the next section, we review the Van Kampen inhomogeneous diffusion equation for Brownian particles and its stationary solution. In Section 2.1, we study a single massive BH and in Section 2.2, we consider a BH population immersed in a fixed Plummer profile. In the final section, we discuss our conclusions. In Appendix A, we re-derive the Van Kampen diffusion equation.

2 Inhomogeneous diffusion of black holes in stellar clusters

A standard approach in gravitational kinetic theory is to apply the orbit-averaged Fokker-Planck equation (e.g. Merritt 2013; Binney & Tremaine 2008). This approach comes with a drawback that Merritt (2013)¹¹1In page 253, section 5.5. refers to as a “kludge.” He remarks that while the diffusion coefficients are averaged over the volume filled by an orbit, no account is taken of the fact that the perturbations acting on a test star come from stars in regions far from the test star, which may have very different properties than those near to the test star. For this reason he concludes there is no sense in which the orbit-averaged Fokker-Planck equation might correctly describe an inhomogeneous system. We wish to describe the motion of heavier – that is, Brownian – bodies inside an inhomogeneous (in phase-space) medium of lighter bodies and, in addition, we wish to focus on the distribution of position. Therefore, we shall not use the standard orbit-averaged approach with $E$ , $L$ degrees of freedom, but the older approach from Chandrasekhar (1943a, b) in the phase-space $x$ , $v$ (see also section 5.4 of Merritt 2013). We nevertheless generalize Chandrasekhar (1943b) approach here in order to account for inhomogeneity in space and temperature.

This may be achieved by considering the Kramers equation of the Fokker-Planck type (Kramers 1940) in a local sense (Equation 34) and derive a diffusion equation in a way that is similar to the derivation of the Smolukowski equation from the Kramers equation, as done by Chandrasekhar (1943b). This problem was already resolved by van Kampen (1988) and we provide a review in Appendix A.

Therefore, we suggest that the Van Kampen diffusion equation for Brownian motion in inhomogeneous media (van Kampen 1988) applies to self-gravitating systems, such as stellar clusters, which is the focus of this work, as well as dark matter haloes. These are typically spatially inhomogeneous and non-isothermal and they acquire a velocity dispersion that is also varying. If we denote $\mu(x)\equiv 1/m\eta(x),$ with $\eta(x)$ as the diffusion coefficient and $m$ the mass of a Brownian particle, $\Phi(x)$ and $T(x)$ the potential and local temperature of the system in which the Brownian particle is embedded and $p(x)$ the probability density of the Brownian particle in one dimenion $x$ , the Van Kampen diffusion equation reads:

\frac{\partial p(x,t)}{\partial t}=\frac{\partial}{\partial x}\left\{\mu(x)\left(m\Phi(x)^{\prime}p(x,t)+\frac{\partial}{\partial x}(T(x)p(x,t))\right)\right\}.

(1)

As we noted above, we reproduce this equation in Appendix A. To be more specific, it is derived by expanding the Fokker-Planck equation (34 of Kramers with respect to $\eta^{-1}$ , defined in the Langevin equation (36). This quantity is given in units of time and corresponds to the timescale of the diffusion, as was suggested by Chandrasekhar (1943a). Therefore, the Van-Kampen equation is strictly valid for times $t\gg\eta^{-1}$ .

The stationary distribution

\frac{\partial p}{\partial t}=0

(2)

of Equation 1 is

p(x)=\frac{p(0)}{\beta(0)}\beta(x)e^{-\int_{0}^{x}d\tilde{x}\,\beta(\tilde{x})m\Phi(\tilde{x})^{\prime}},

(3)

where it is assumed that $p(0)^{\prime}=T(0)^{\prime}=\Phi(0)^{\prime}=0$ . We denoted $\beta=1/T$ . In three dimensions, the diffusion Equation 1 becomes

\frac{\partial p(\bf{r},t)}{\partial t}=\nabla\left\{\mu({\bf{r}})\left(m(\nabla\Phi({\bf{r}}))p({\bf{r}},t)+\nabla(T({\bf{r}})p({\bf{r}},t))\right)\right\}.

(4)

In the spherically symmetric case, the stationary radial probability density can therefore be written as

p(r)=\frac{\beta(r)e^{-\int_{0}^{r}d\tilde{r}\,\beta(\tilde{r})m\Phi(\tilde{r})^{\prime}}}{\int_{0}^{\infty}dr\,4\pi r^{2}\beta(r)e^{-\int_{0}^{r}d\tilde{r}\,\beta(\tilde{r})m\Phi(\tilde{r})^{\prime}}}.

(5)

The result is that at equilibrium, the probability density does not depend on $\mu(r)$ , but only on $\beta(r)$ and $\Phi(r)$ . We note that the stationary phase space distribution function may be written at first order in $\eta(x)^{-1}$ by use of Equations (52), (54), and (58):

	$\displaystyle f(x,v)$	$\displaystyle=p(x)\sqrt{\frac{m}{2\pi T(x)}}e^{-\frac{mv^{2}}{2T(x)}}$
		$\displaystyle\left\{1-\eta(x)^{-1}\left\{v\left(\frac{1}{2}\frac{T(x)^{\prime}}{T(x)}p(x)+p(x)^{\prime}\right)+v^{3}\frac{mT(x)^{\prime}}{T(x)^{2}}p(x)\right\}\right\}.$		(6)

The temperature, $T(x),$ refers to the host cluster. There a small correction term arises to the Maxwellian distribution. This requires further investigation and may be the subject of interesting future developments, but it is not discussed in this work.

2.1 Single black hole

Refer to caption — Figure 1: Stationary probability distribution $p_{\rm b}$ of inhomogeneous diffusion (5) with respect to distance, $r,$ for a single Brownian body, $m_{\rm b}$ , inside a Plummer external gravitational potential. The Brownian body may be considered to be a massive BH and the gravitational potential may be generated by a globular cluster of individual average stellar mass, $m_{\star}$ . We denote $r_{c}$ the softening radius of the potential.

Here, we consider a single BH of mass, $m_{\rm b}$ , as a Brownian particle immersed inside a stellar cluster of total mass, $M_{\star},$ and average individual stellar mass, $m_{\star}$ .

We model the distribution of the stellar cluster with a Plummer sphere,

	$\displaystyle\rho_{\star}(x)$	$\displaystyle=\frac{M_{\star}}{\frac{4}{3}\pi r_{c}^{3}}\frac{1}{(1+x^{2})^{5/2}},\;\mathcal{M}_{\star}(x)=M_{\star}\frac{x^{3}}{(1+x^{2})^{3/2}},$		(7)
	$\displaystyle\sigma_{\star}(x)^{2}$	$\displaystyle=\frac{GM_{\star}}{6r_{c}}\frac{1}{(1+x^{2})^{1/2}},$		(8)

where $\rho_{\star}(x)$ , $\sigma_{\star}(x)$ denote the mass density and velocity dispersion of the host cluster, respectively, while $\mathcal{M}_{\star}(x)$ denotes the total stellar mass contained within $x$ . The half-mass radius is related to the softening radius $r_{c}$ by

r_{\rm hm}=(2^{2/3}-1)^{-1/2}r_{c}\simeq 1.3r_{c}.

(9)

We identify the temperature as

T_{\star}(r)=m_{\star}\sigma_{\star}(r)^{2}.

(10)

The probability density of the BH position is given in a straightforward way from Equation (5) which, in the case of Plummer external potential, may be written as

p_{\rm b}=p_{\rm b}(0)(1+x^{2})^{-\frac{1}{2}-3\frac{m_{\rm b}}{m_{\star}}},

(11)

where

p_{\rm b}(0)=\frac{1}{4\pi r_{c}^{3}}\frac{4\Gamma\left(3\frac{m_{\rm b}}{m_{\star}}+\frac{1}{2}\right)}{\sqrt{\pi}\Gamma\left(3\frac{m_{\rm b}}{m_{\star}}-1\right)},

(12)

and $\Gamma$ is the gamma-function. In Figure 1, we plot the probability density with respect to $x$ for several different values of the individual mass ratio $m_{\rm b}/m_{\star}$ . For a dense stellar cluster, such as a globular cluster or a nuclear star cluster, it is $r_{c}\sim 1-5{\rm pc}$ and $m_{\star}\sim 0.5{\rm M}_{\odot}$ . It is evident that an intermediate-mass BH with $m_{\rm b}\sim 10^{3}{\rm M}_{\odot}$ inside a globular cluster wanders as far as $\sim 0.05{\rm pc}$ from the center, and in a nuclear star cluster, as far as $\sim 0.2{\rm pc}$ , although a Plummer sphere might not be a fairly good approximation for the density of the latter. A steeper external profile would induce smaller diffusion and, therefore, stricter boundaries for the BH fluctuation.

In order to verify, in a straightforward manner, the validity of the inhomogeneous diffusion Equation (1), we performed direct numerical $N$ -body simulations using the AMUSE 13.2 framework (Portegies Zwart et al. 2009, 2013; Pelupessy et al. 2013; Portegies Zwart & McMillan 2018). We considered a stellar cluster with $N_{\star}=10^{3}$ , $m_{\star}=0.5{\rm M}_{\odot}$ , $r_{c}=0.5{\rm pc}$ and a BH with $m_{\rm b}=50{\rm M}_{\odot}$ placed in the center of the cluster at distance $r_{\rm b,ini}=1{\rm AU}$ and zero velocity. We use this initial condition for the BH to demonstrate that even if it is left sitting still in the center, it will eventually be kicked out due to gravitational fluctuations. In any case, we find that the distribution function seems to not be affected by the BH initial conditions, at least when it is initially bound inside the core $\sim 0.6r_{c}$ . We performed 1000 runs, initially sampling the stellar cluster from an isotropic Plummer phase-space distribution function $f_{\rm Pl}(x,v)\sim(-E(x,v))^{7/2}$ . In our previous calculation (Equation 11), we assumed the Plummer potential to be fixed. Therefore, our sampling of the BH position should occur at a time is shorter than the relaxation time of the cluster $t_{\rm rh}=0.138(N_{\star}/\ln\Lambda)(G\rho_{\star})^{-1/2}$ (e.g., Heggie & Hut 2003; Portegies Zwart & McMillan 2018) so that the deviation from Plummer distribution is small. The relaxation time for a Plummer sphere is (Section 8, Table 1 of Heggie & Hut 2003):

t_{\rm rh}=\frac{0.206}{\ln\Lambda}\frac{N_{\star}r_{c}^{3/2}}{\sqrt{GM_{\star}}}=10{\rm Myr},

(13)

where we used $\Lambda=0.11N_{\star}$ (Portegies Zwart & McMillan 2018). The BH diffusion timescale $t_{\rm b,rel}$ is given from the two-body relaxation timescale (in a closely related context see El-Zant et al. 2020, 2016):

t_{\rm b,rel}=\frac{1}{8\pi\ln\Lambda}\frac{\sigma_{\star}^{3}}{G^{2}\rho_{\star}m_{\rm b}}.

(14)

Dividing the two timescales, we get at the center $r=0:$

\frac{t_{\rm b,rel}}{t_{\rm rh}}=0.05\frac{m_{\star}}{m_{\rm b}}.

(15)

This verifies that the diffusion Equation (1) is valid at times that are comparable to the relaxation timescale of a cluster for $m_{\star}<m_{\rm b}$ . Now, for the parameters in our simulation, the last equation gives $t_{\rm b,rel}=5\cdot 10^{-4}t_{\rm rh}=5\cdot 10^{-3}{\rm Myr}$ . We acquire the BH position at each simulation run within the time window, $t=0.3\pm 0.05{\rm Myr}$ , which is sufficiently bigger than $t_{\rm b,rel}$ , so that the BH has enough time to diffuse, and sufficiently smaller than $t_{\rm rh}$ , so that the Plummer potential does not get drastically altered.

We calculate the probability density by assuming the ergodic hypothesis holds (at each run, we record several positions in time) and performing an ensemble average. The fit of the theory of inhomogeneous diffusion to the simulations is very good, as depicted in Figure 2. In Figure 2(a), we demonstrate in a single simulation run that the BH rapidly (at $\sim t_{\rm b,rel}=0.005{\rm Myr}$ ) diffuses in the medium, while in Figure 2(b), we present the BH probability density of 1000 simulation runs. The RMS position $r_{\rm RMS}\equiv<r_{\rm b}^{2}>^{1/2}$ of the simulations $r_{\rm RMS}^{\rm sim}=0.037{\rm pc}$ is well-matched by the theoretical prediction $r_{\rm RMS}^{\rm theory}=0.036{\rm pc}$ in the chosen time window.

2.2 Black hole subcluster

Here, we consider a two-component model, namely a population of $N_{\rm b}$ Brownian particles with average mass, $m_{\rm b}$ embedded in a cluster of $N_{\star}$ bodies with average mass, $m_{\star}$ . Our description wishes to describe a BH population immersed in a globular cluster, but applies generally to any self-gravitating system which hosts a subsystem of significantly less bodies. The total mass of the host is $M_{\star}=N_{\star}m_{\star}$ and the total mass of the Brownian bodies is $M_{\rm b}=N_{\rm b}m_{\rm b}$ .

The potential $\Phi$ includes the potential of the host cluster, but also can account for the self-gravity of the Brownian population if not negligible. Given the distribution of the host cluster and neglecting the feedback of the Brownian population onto the distribution of the host, the diffusion equation (1) together with the Poisson equation form a system of equations that determines the equilibrium distribution of the Brownian particles. This system may be formulated as follows. In the spherically symmetric case, the gravitational field at any point $r$ may be decomposed according to Poisson equation as

\frac{d\Phi}{dr}=G\frac{\mathcal{M}_{\star}(r)+\mathcal{M}_{\rm b}(r)}{r^{2}},

(16)

where $\mathcal{M}_{\star}$ , $\mathcal{M}_{\rm b}$ denote the mass contained in $r$ of the host cluster and the subcluster respectively. The equilibrium corresponds to the stationary solution of the inhomogeneous diffusion equation (1). We therefore get the system for the Brownian population:

	$\displaystyle\frac{dp_{\rm b}(r)}{dr}$	$\displaystyle=-p_{\rm b}(r)\left(Gm_{\rm b}\frac{\mathcal{M}_{\star}(r)+\mathcal{M}_{\rm b}(r)}{r^{2}T_{\star}(r)}+\frac{T_{\star}(r)^{\prime}}{T_{\star}(r)}\right),$		(17)
	$\displaystyle\frac{d\mathcal{M}_{\rm b}(r)}{dr}$	$\displaystyle=4\pi r^{2}M_{\rm b}p_{\rm b}(r).$		(18)

This is a system of the unknown distributions $\{p_{b},\mathcal{M}_{\rm b}\}$ given the host cluster distributions $\mathcal{M}_{\star}(r)$ and $T_{\star}(r)$ and subject to the constraint $\mathcal{M}_{b}(\infty)=M_{\rm b}$ . The constraint suggests that equilibria may exist only for certain parameter values. This formulation, where the host cluster distribution is fixed, does not take into account the feedback of the Brownian population onto the host. This may be significant in the proximate regions close to the population if it is sufficiently compact, but we do not be consider this point in this work.

Supposing that the system is characterised by a length scale, $r_{c}$ , we may introduce the dimensionless variables:

	$\displaystyle x$	$\displaystyle=\frac{r}{r_{c}},\;\tilde{M}_{\rm b}(x)=\frac{\mathcal{M}_{\rm b}(x)}{M_{\rm b}},\;\tilde{M}_{\star}(x)=\frac{\mathcal{M}_{\star}(x)}{M_{\star}},$		(19)
	$\displaystyle\tilde{p}_{\rm b}(x)$	$\displaystyle=4\pi r_{c}^{3}p_{\rm b},\;y(x)=-\ln\frac{\tilde{p}_{\rm b}(x)}{\tilde{p}_{\rm b}(0)}.$		(20)

The system (17)-(18) becomes

	$\displaystyle\frac{dy(x)}{dx}$	$\displaystyle=\frac{Gm_{\rm b}M_{\star}}{r_{c}T_{\star}(x)}\frac{1}{x^{2}}\left(\tilde{M}_{\star}(x)+\frac{M_{b}}{M_{\star}}\tilde{M}_{\rm b}(x)\right)+\frac{T_{\star}(x)^{\prime}}{T_{\star}(x)},$		(21)
	$\displaystyle\frac{d\tilde{M}_{\rm b}(x)}{dx}$	$\displaystyle=\tilde{p}_{\rm b}(0)x^{2}e^{-y(x)},$		(22)

with initial conditions

y(0)=0,\quad\tilde{M}_{\rm b}(0)=0,\quad y(0)^{\prime}=0,\quad\tilde{M}_{\rm b}(0)^{\prime}=0

(23)

and subject to the boundary constraint

\tilde{M}_{\rm b}(\infty)=1\Leftrightarrow\tilde{p}_{\rm b}(0)=\left(\int_{0}^{\infty}dx\,x^{2}e^{-y(x)}\right)^{-1}.

(24)

For given host distributions $\tilde{M}_{\star}(x)$ , $T(x)$ the system (21)-(24) may be solved numerically.

We consider again the case where the stellar distribution is that of a Plummer sphere, as in Equation 7. The temperature is identified as $T_{\star}(r)=m_{\star}\sigma_{\star}(r)^{2}$ . The system (21)-(22) becomes

	$\displaystyle\frac{dy(x)}{dx}$	$\displaystyle=6\frac{m_{\rm b}}{m_{\star}}\left(\frac{x}{1+x^{2}}+\frac{M_{\rm b}}{M_{\star}}\frac{\sqrt{1+x^{2}}}{x^{2}}\tilde{M}_{b}(x)\right)-\frac{x}{1+x^{2}},$		(25)
	$\displaystyle\frac{d\tilde{M}_{\rm b}(x)}{dx}$	$\displaystyle=\tilde{p}_{\rm b}(0)x^{2}e^{-y(x)},$		(26)

subject again to the conditions (23)-(24).

The reach (or otherwise) of an equilibrium and the specific form of the equilibrium distribution of the Brownian particles, which we consider hereof to be a BH population, depend on both the number of bodies and their individual mass. In Figure3(a), we calculate the series of equilibria of BH populations immersed in a Plummer stellar profile for various ratios $m_{\rm b}/m_{\star}$ . Different points of each curve correspond to the equilibrium state of a different BH population with total mass, $M_{\rm b}$ , and corresponding central density $\tilde{p}(0)$ . We identify an instability that sets in at the maximum of the series of equilibria curve. No equilibria exist above the maximum, whereas, according to the Poincaré theorem of linear series of equilibria, the branch beyond the turning point (dotted lines in Figure 3(a)) are unstable.

We further discover numerically that when the BH subcluster mass is scaled to become the Spitzer parameter,

S\equiv\frac{M_{\rm b}m_{\rm b}^{3/2}}{M_{\star}m_{\star}^{3/2}}

(27)

and the central probability density is scaled as

B\equiv\tilde{p}_{\rm b}(0)\left(\frac{m_{\star}}{m_{\rm b}}\right)^{3/2},

(28)

then all curves with different ratios, $m_{\rm b}/m_{\star},$ merge to form a single one, as in Figure 3(b) (we estimate the exact value of the exponent to be $1.51$ and not $3/2$ , but we consider this small deviation a numerical precision effect without any physical significance). Thus, the curve 3(b) is global and applies to all Brownian populations immersed in a Plummer profile. The instability sets in at

B_{\rm I}=140,\;S_{\rm I}=0.25.

(29)

Any stationary equilibrium with

B>B_{\rm I}

(30)

is unstable within the current framework. In addition no stationary equilinbrium states exist above

S>S_{\rm I}.

(31)

This value $S_{\rm I}$ is significantly larger than the Spitzer value (Spitzer 1969), $S_{\rm SP}=0.16,$ calculated for isothermal equipartition.

We stress at this point that while we assume here that the BH population to be immersed inside a fixed density profile of the cluster, it should, in practice, affect the cluster profile at least in the vicinity of BH subcluster’s denser regions. This effect, not taken into account in the current work, may influence the cluster’s ability to support the BH population. Nevertheless, it seems possible that the BH population will locally heat up and inflate the population of lighter bodies in an effect that is similar to osmosis. This could possibly enhance, rather than reduce, a cluster’s ability to support the BH population via gravitational fluctuations and could lead to formation of a core (as e.g., proposed by Merritt et al. 2004). Such BH population feedback onto the cluster density profile requires further investigation and is not studied here.

The density profile of a BH population at the onset of the instability appears in Figure 4, where the Plummer profile, $\rho_{\star}$ , is also plotted . It is evident that the BH subcluster may extend up to $0.2{\rm pc}$ inside the host cluster and it may be almost twice as dense in the centre.

In Figure 6, we plot the maximum number of Brownian bodies, that is, BHs, in our context and with respect to their individual average mass for inhomogeneous diffusion and Spitzer instability. We assume $m_{\star}=0.5{\rm M}_{\odot}$ and consider two cases of $N_{\star}=10^{6},\,10^{8}$ . In the typical case of a globular cluster with $N_{\star}=10^{6}$ and $m_{\rm b}=10{\rm M}_{\odot},$ for an inhomogeneous diffusion we get $N_{\rm b}\sim 1400,$ while the Spitzer threshold is significantly lower at $N_{\rm b}\sim 890$ .

In Figure 7, we compare the probability distribution of a single intermediate-mass BH with BH populations of lighter stellar mass BHs but with equal total mass. It is evident from Figure 7(a) that in order to be able to discriminate a BH population from an intermediate-mass BH with a mass of $500m_{\star}$ , it may be required to probe inside the inner $0.1{\rm pc}$ . In Figure 7(b), we show that for a BH population with $N_{\rm b}=1090$ , $m_{\rm b}=10{\rm M}_{\odot}$ the first unstable profile is very similar to that of an intermediate-mass BH of equal total mass, while in order to discover observationally the second unstable profile, if stabilized by other processes, one has to probe the inner $\sim 0.005{\rm pc}$ .

3 Conclusions

In this work, we study inhomogeneous diffusion, that is, diffusion in a medium with varying mass density and temperature and, hence, also varying damping and diffusion coefficients. This is the typical case for self-gravitating systems that are spatially inhomogeneous and trapped in states with varying velocity dispersions. We argue that the inhomogeneous diffusion equation (4) applies to self-gravitating systems that involve a sub-population of fewer and heavier bodies and describes gravitational Brownian motion. The corresponding stationary states are given in Equations 5 and 2. We calculated the spatial probability distribution function of a single Brownian particle immersed in a Plummer profile, as depicted in Figure 1. We estimate that a single intermediate-mass BH may wander as far as $\sim 0.05{\rm pc}$ in a typical Globular cluster, while a single typical stellar BH may wander even ten times farther. We validate the reality of inhomogeneous diffusion of BHs immersed in stellar clusters, expressed by Equation 1, with the use of $N$ -body simulations, as in Figure 2. The theoretical prediction is well-matched by the simulation results.

Applying our framework to a Brownian population of massive bodies (focusing on BHs) inside a stellar cluster, which follows a Plummer density profile, we identify an instability that sets in for Brownian populations with Spitzer parameter, $S_{\rm I}=0.25$ , and a new global parameter, $B_{\rm I}=140$ , which depends on the central density of the Brownian population. This is depicted in Figure 3. For $B>B_{\rm I}$ , any stationary equilibrium state is unstable despite the fluctuations of the gravitational field. For $S>S_{\rm I},$ no stationary states exist. This is a manifestation of the Spitzer instability, reinterpreted as the inability of the cluster to support the sub-population of heavier bodies by gravitational fluctuations. The dependence of the onset of the instability on the individual mass ratios in such a framework arises naturally. Furthermore, since the ordinary Spitzer instability occurs at a lower value, $S_{\rm SP}=0.16$ , the inhomogeneous diffusion allows more massive BH populations to reach stationary states than one would expect from isothermal equipartition.

An important limitation of our model for BH populations in Section 2.2 regarding its physical applicability is that it does not take into account the feedback of the BH population to the gravitational potential of the host cluster. Such a limitation for analyses on the Brownian motion of a single massive BH was also noted in (Merritt et al. 2007). The situation is very similar since our BH population has the same mass of our assumed single massive BH and it turns out to have about the same spatial probability distribution. Still, it was our intention to neglect this feedback in order to quantify the effect solely of the gravitational fluctuations to the BH population and inspect their stabilizing efficiency. If anything, it seems plausible that the BH population will locally heat up and inflate the population of lighter bodies, as was also suggested by Merritt et al. (2004). This will enhance, and not reduce, the cluster’s ability to support the BH population via gravitational fluctuations and lead to the formation of a core. We remark that in his seminal work on Brownian motion, Einstein (1905, 1956) interpreted it precisely as a response to osmotic pressure. Osmosis involves the diffusion of the solvent particles (in our case the stars) in the region of the solute (in our case, the BH population), with result the inflation of the region containing the solute. Such BH population feedback onto the cluster density profile, including the possible reality of such a phenomenon as ”gravitational osmosis,” requires further investigation. Another future improvement could include following the time evolution of the diffusion equation (4) for the BH population.

In conclusion, the fact that a BH subcluster can be supported by random fluctuations of the gravitational field beyond the limit of Spitzer instability threshold supports the idea that globular clusters can retain a significant BH population.

References

Abbate et al. (2019) Abbate, F., Possenti, A., Colpi, M., & Spera, M. 2019, ApJL, 884, L9
Arca-Sedda (2016) Arca-Sedda, M. 2016, MNRAS, 455, 35
Arca Sedda et al. (2018) Arca Sedda, M., Askar, A., & Giersz, M. 2018, MNRAS, 479, 4652
Askar et al. (2018) Askar, A., Arca Sedda, M., & Giersz, M. 2018, MNRAS, 478, 1844
Bahramian et al. (2017) Bahramian, A., Heinke, C. O., Tudor, V., et al. 2017, MNRAS, 467, 2199
Barnard et al. (2008) Barnard, R., Stiele, H., Hatzidimitriou, D., et al. 2008, ApJ, 689, 1215
Binney & Lacey (1988) Binney, J. & Lacey, C. 1988, MNRAS, 230, 597
Binney & Tremaine (2008) Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press)
Bortolas et al. (2016) Bortolas, E., Gualandris, A., Dotti, M., Spera, M., & Mapelli, M. 2016, MNRAS, 461, 1023
Breen & Heggie (2013) Breen, P. G. & Heggie, D. C. 2013, MNRAS, 432, 2779
Chandrasekhar (1943a) Chandrasekhar, S. 1943a, ApJ, 97, 255
Chandrasekhar (1943b) Chandrasekhar, S. 1943b, Reviews of Modern Physics, 15, 1
Chatterjee et al. (2002a) Chatterjee, P., Hernquist, L., & Loeb, A. 2002a, Phys. Rev. Lett., 88, 121103
Chatterjee et al. (2002b) Chatterjee, P., Hernquist, L., & Loeb, A. 2002b, ApJ, 572, 371
Chatterjee et al. (2003) Chatterjee, P., Hernquist, L., & Loeb, A. 2003, ApJ, 592, 32
Chatterjee et al. (2017) Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2017, ApJ, 834, 68
Chavanis (2012) Chavanis, P.-H. 2012, Physica A: Statistical Mechanics and its Applications, 391, 3680–3701
Chavanis (2013) Chavanis, P. H. 2013, Astron. Astrophys., 556, A93
Chomiuk et al. (2013) Chomiuk, L., Strader, J., Maccarone, T. J., et al. 2013, ApJ, 777, 69
Di Cintio et al. (2020) Di Cintio, P., Ciotti, L., & Nipoti, C. 2020, in IAU Symposium, Vol. 351, IAU Symposium, ed. A. Bragaglia, M. Davies, A. Sills, & E. Vesperini, 93–96
Einstein (1905) Einstein, A. 1905, Annalen der Physik, 322, 549
Einstein (1956) Einstein, A. 1956, Investigations on the theory of Brownian movement (Dover Publications, Inc.)
El-Zant et al. (2016) El-Zant, A. A., Freundlich, J., & Combes, F. 2016, MNRAS, 461, 1745
El-Zant et al. (2020) El-Zant, A. A., Freundlich, J., Combes, F., & Halle, A. 2020, MNRAS, 492, 877
Fouvry & Bar-Or (2018) Fouvry, J.-B. & Bar-Or, B. 2018, MNRAS, 481, 4566
Giesers et al. (2018) Giesers, B., Dreizler, S., Husser, T.-O., et al. 2018, MNRAS, 475, L15
Heggie & Hut (2003) Heggie, D. & Hut, P. 2003, The Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics (Cambridge University Press)
Irwin et al. (2010) Irwin, J. A., Brink, T. G., Bregman, J. N., & Roberts, T. P. 2010, ApJL, 712, L1
Kramers (1940) Kramers, H. A. 1940, Physica, 7, 284
Kremer et al. (2019) Kremer, K., Chatterjee, S., Ye, C. S., Rodriguez, C. L., & Rasio, F. A. 2019, ApJ, 871, 38
Kremer et al. (2018) Kremer, K., Ye, C. S., Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2018, ApJL, 855, L15
Kulkarni et al. (1993) Kulkarni, S. R., Hut, P., & McMillan, S. 1993, Nature, 364, 421
Lau & Binney (2019) Lau, J. Y. & Binney, J. 2019, MNRAS, 490, 478
Lingam (2018) Lingam, M. 2018, MNRAS, 473, 1719
Maccarone et al. (2007) Maccarone, T. J., Kundu, A., Zepf, S. E., & Rhode, K. L. 2007, Nature, 445, 183
Mackey et al. (2008) Mackey, A. D., Wilkinson, M. I., Davies, M. B., & Gilmore, G. F. 2008, MNRAS, 386, 65
Merritt (2005) Merritt, D. 2005, ApJ, 628, 673
Merritt (2013) Merritt, D. 2013, Dynamics and evolution of galactic nuclei (Princeton, N.J: Princeton University Press)
Merritt et al. (2007) Merritt, D., Berczik, P., & Laun, F. 2007, AJ, 133, 553
Merritt et al. (2004) Merritt, D., Piatek, S., Portegies Zwart, S., & Hemsendorf, M. 2004, ApJL, 608, L25
Miller-Jones et al. (2015) Miller-Jones, J. C. A., Strader, J., Heinke, C. O., et al. 2015, MNRAS, 453, 3918
Minniti et al. (2015) Minniti, D., Contreras Ramos, R., Alonso-García, J., et al. 2015, ApJL, 810, L20
Morscher et al. (2015) Morscher, M., Pattabiraman, B., Rodriguez, C., Rasio, F. A., & Umbreit, S. 2015, ApJ, 800, 9
Morscher et al. (2013) Morscher, M., Umbreit, S., Farr, W. M., & Rasio, F. A. 2013, ApJL, 763, L15
Pelupessy et al. (2013) Pelupessy, F. I., van Elteren, A., de Vries, N., et al. 2013, Astron. Astrophys., 557, A84
Portegies Zwart & McMillan (2018) Portegies Zwart, S. & McMillan, S. 2018, Astrophysical Recipes, 2514-3433 (IOP Publishing)
Portegies Zwart et al. (2009) Portegies Zwart, S., McMillan, S., Harfst, S., et al. 2009, NewA, 14, 369
Portegies Zwart et al. (2013) Portegies Zwart, S., McMillan, S. L. W., van Elteren, E., Pelupessy, I., & de Vries, N. 2013, Computer Physics Communications, 184, 456
Rauch & Tremaine (1996) Rauch, K. P. & Tremaine, S. 1996, NewA, 1, 149
Roberts et al. (2012) Roberts, T. P., Fabbiano, G., Luo, B., et al. 2012, ApJ, 760, 135
Rodriguez et al. (2016) Rodriguez, C. L., Morscher, M., Wang, L., et al. 2016, MNRAS, 463, 2109
Roupas (2020) Roupas, Z. 2020, Journal of Physics A Mathematical General, 53, 045002
Roupas et al. (2017) Roupas, Z., Kocsis, B., & Tremaine, S. 2017, ApJ, 842, 90
Shishkovsky et al. (2018) Shishkovsky, L., Strader, J., Chomiuk, L., et al. 2018, ApJ, 855, 55
Sigurdsson & Hernquist (1993) Sigurdsson, S. & Hernquist, L. 1993, Nature, 364, 423
Spitzer (1969) Spitzer, Lyman, J. 1969, ApJL, 158, L139
Strader et al. (2012) Strader, J., Chomiuk, L., Maccarone, T. J., Miller-Jones, J. C. A., & Seth, A. C. 2012, Nature, 490, 71
Taylor et al. (2015) Taylor, M. A., Puzia, T. H., Gomez, M., & Woodley, K. A. 2015, ApJ, 805, 65
Tremaine (2020) Tremaine, S. 2020, MNRAS, 493, 2632
Tremaine & Weinberg (1984) Tremaine, S. & Weinberg, M. D. 1984, MNRAS, 209, 729
van Kampen (1988) van Kampen, N. 1988, Journal of Physics and Chemistry of Solids, 49, 673
Van Kampen (1985) Van Kampen, N. G. 1985, Phys. Rep. , 124, 69
Vasiliev (2017) Vasiliev, E. 2017, The Astrophysical Journal, 848, 10
Wang et al. (2016) Wang, L., Spurzem, R., Aarseth, S., et al. 2016, MNRAS, 458, 1450
Weatherford et al. (2019) Weatherford, N. C., Chatterjee, S., Kremer, K., & Rasio, F. A. 2019, arXiv e-prints, arXiv:1911.09125
Weatherford et al. (2018) Weatherford, N. C., Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2018, ApJ, 864, 13

Appendix A Diffusion equation for gravitational Brownian motion

We reproduce here the diffusion equation of Brownian motion inside an inhomogeneous medium with varying temperature, as was first proposed by van Kampen (1988). We further argue that this diffusion equation describes gravitational Brownian motion. It applies to self-gravitating systems, which are not only inhomogeneous, but also are typically non-isothermal, in the sense that the velocity dispersion varies with position.

As suggested by Chandrasekhar (1943b), the gravitational field, $\bm{g}(\bm{r},t),$ at a point, $\bm{r},$ and time, $t,$ of an $N$ -body self-gravitating system may be decomposed to the sum of a mean field $\bm{g}_{\rm m}$ and a fluctuating field $\bm{g}_{\rm f:}$

\bm{g}(\bm{r},t)=\bm{g}_{\rm m}(\bm{r},t)+\bm{g}_{\rm f}(\bm{r},t).

(32)

The mean field represents the effect of the system as a whole to each point of space at any instant of time through the smoothed out continuous mass density function $\rho(\bm{r},t):$

\bm{g}_{\rm m}=-\nabla\Phi,\quad\Phi(\bm{r},t)=-\int d\tilde{\bm{r}}\,G\frac{\rho(\bm{r},t)}{|\bm{r}-\tilde{\bm{r}}|}.

(33)

The fluctuating field accounts for the deviations from this smoothed out field that arise due to the granularity of the system. Chandrasekhar (1943a) showed further that each body is subject to a dynamical friction force with a damping coefficient, $\eta,$ which is reciprocal to the relaxation timescale of the system. Central to his derivation is a Fokker-Planck type of equation discovered by Kramers, namely,

\left(\frac{\partial}{\partial t}+v\frac{\partial}{\partial x}-\Phi^{\prime}\frac{\partial}{\partial v}\right)f(x,v,t)=\eta\frac{\partial}{\partial v}\left(v+\frac{T}{m}\frac{\partial}{\partial v}\right)f(x,v,t),

(34)

which we consider here in one-dimension for simplicity without loss of generality. We denote $f$ as the phase-space probability distribution function, $m$ the mass of the Brownian particle, $T$ the temperature in units $k_{\rm B}=1,$ and $\eta$ the damping coefficient. Equation 34 was derived by Kramers (1940) according to very general considerations. Chandrasekhar (1943b) rederived alternatively Kramers equation and used it in self-gravitating systems to establish the concept of dynamical friction (Chandrasekhar 1943a). Here, following van Kampen (1988), we consider Kramers equation, but with varying damping (dynamical friction in our case) coefficient and temperature,

\eta=\eta(x),\quad T=T(x).

(35)

From Equation 34, we derive the diffusion equation, which is to be a generalization of the Smolukowski equation and we argue that it describes gravitational Brownian motion.

First, for the purposes of completeness, we reproduce the derivation of (34) following Chandrasekhar (1943b). We assume that the motion of a Brownian particle is described by a Langevin equation,

\frac{dv}{dt}=-\eta v+g_{\rm m}+g_{\rm f},

(36)

where $\eta$ is the dynamical friction coefficient (Chandrasekhar 1943a), $g_{\rm m}=-\Phi^{\prime}$ is the mean field and $g_{\rm f}$ the fluctuating field. The latter term has statistically defined properties rendering the Langevin equation a stochastic differential equation that cannot be solved as an ordinary differential equation might be. Instead, we would aim to derive a distribution function $f(x,v,t+\Delta t)$ governing the probability of occurrence of $x$ , $v$ at time $t+\Delta t$ given the distribution, $f(x,v,t),$ at time, $t$ .

We let $\Delta t$ be short with respect to the relaxation time $\Delta t\ll\eta^{-1}$ but long compared to the period of $g_{\rm f}$ fluctuations. We define the velocity variation due to field fluctuation within $\Delta t$ ,

{\rm B}(\Delta t)=\int_{t}^{t+\Delta t}g_{\rm f}(\xi)d\xi.

(37)

Physically, it describes the net acceleration exerted on the Brownian particle arising from field fluctuations in time, $\Delta t$ . We impose the physical requirement that for $t\rightarrow\infty$ , i.e. $t\gg\eta^{-1}$ that is at zeroth order of $\eta^{-1}$ , the distribution function, $f,$ is locally Maxwellian. Chandrasekhar (1943b) has proven, as in his Chapter II.2, that this requirement implies that ${\rm B}$ satisfies the Brownian probability distribution:

\psi({\rm B})=\left(\frac{m}{4\pi\eta T\Delta t}\right)^{1/2}e^{-\frac{m}{4\eta T\Delta t}{\rm B}^{2}}.

(38)

The spatial and velocity displacements, $\Delta x$ , $\Delta v,$ at $\Delta t$ are derived by the Langevin equation (36), so that we get

\Delta x=v\Delta t,\quad{\rm B}=\Delta v+(\eta v+\Phi^{\prime})\Delta t.

(39)

The latter equation follows by integration of the Langevin equation (36) from $t$ to $t+\Delta t$ , introducing ${\rm B}$ by use of (37), and finally solving with respect to ${\rm B}$ .

Given that Brownian motion is governed by two-body relaxation and it proceeds in a short timescale, it is reasonably expected to be modelled as a Markoff process²²2We caution the reader that the Markoff assumption may not be valid for any gravitational system. In the case of regular orbits, e.g. in clusters whose self-gravity is dominated by a central potential, the relaxation of orbital angular momentum may proceed in a different timescale than the energy and encounters between stars may be coherent, as in resonant relaxation (Rauch & Tremaine 1996). There is recent evidence of resonant relaxation occuring in globular clusters, questioning also the local-scattering relaxation theory (Lau & Binney 2019). as is pointed out also by Chandrasekhar (1943a). In this case the probability distribution at any time, $t+\Delta t,$ can be derived from the distribution at previous time, $t$ . Using Equation 39, the transition probability, $\psi(v-\Delta v;\Delta v),$ that $v$ changes by $\Delta v$ may be expressed with respect only to $\Delta v$ . Then we have:

f(x+v\Delta t,v,t+\Delta t)=\int d(\Delta v)\,f(x,v-\Delta v,t)\psi(x,v-\Delta v;\Delta v),

(40)

where $\psi$ is given from (38) for ${\rm B}(\Delta v)$ given in (39). Taylor expanding $f(x+v\Delta t,v,t+\Delta t)$ , $f(x,v-\Delta v,t)$ and $\psi(x,v-\Delta v;\Delta v)$ we get the Fokker-Planck equation:

\left(\frac{\partial f}{\partial t}+v\frac{\partial f}{\partial x}\right)\Delta t+\mathcal{O}((\Delta t)^{2})=-\frac{\partial(f\left\langle\Delta v\right\rangle)}{\partial v}+\frac{1}{2}\frac{\partial(f\left\langle\Delta v^{2}\right\rangle)}{\partial v^{2}}+\mathcal{O}((\Delta t)^{2}),

(41)

where the mean values $\left\langle\bullet\right\rangle$ are calculated with the distribution $\psi(x,v;\Delta v)$ . We get $\left\langle\Delta v\right\rangle=-(\eta v+\Phi^{\prime})\Delta t$ and $\left\langle(\Delta v)^{2}\right\rangle=(2\eta T/m)\Delta t+\mathcal{O}(\Delta t^{2})$ . Making a substitution in (41), we get finally Kramers equation (34), generalized for an inhomogeneous medium:

\left(\frac{\partial}{\partial t}+v\frac{\partial}{\partial x}-\Phi(x)^{\prime}\frac{\partial}{\partial v}\right)f(x,v,t)=\eta(x)\frac{\partial}{\partial v}\left(v+\frac{T(x)}{m}\frac{\partial}{\partial v}\right)f(x,v,t).

(42)

Now, we derive the corresponding diffusion equation from Kramers equation following van Kampen (1988), who uses a well-established method that was also used to derive the Smolukowski equation from the Kramers equation (Chandrasekhar 1943b), but modified in a straightforward manner to account for the dependence of $\eta$ and $T$ on $x$ . That is a method of elimination of fast variables (Van Kampen 1985) in which we expand the solution to (42) in powers of

\varepsilon(x)\equiv\eta(x)^{-1},

(43)

which is the local relaxation time (Chandrasekhar 1943b). We assume

f=f^{(0)}+f^{(1)}+f^{(2)}+\mathcal{O}(\varepsilon^{3})

(44)

and get up to the second order the equations

		$\displaystyle\frac{\partial}{\partial v}\left(vf^{(0)}+\frac{T}{m}\frac{\partial f^{(0)}}{\partial v}\right)=0$		(45)
		$\displaystyle\varepsilon\left(\frac{\partial f^{(0)}}{\partial t}+v\frac{\partial f^{(0)}}{\partial x}-\Phi^{\prime}\frac{\partial f^{(0)}}{\partial v}\right)=\frac{\partial}{\partial v}\left(vf^{(1)}+\frac{T}{m}\frac{\partial f^{(1)}}{\partial v}\right)$		(46)
		$\displaystyle\varepsilon\left(\frac{\partial f^{(1)}}{\partial t}+v\frac{\partial f^{(1)}}{\partial x}-\Phi^{\prime}\frac{\partial f^{(1)}}{\partial v}\right)=\frac{\partial}{\partial v}\left(vf^{(2)}+\frac{T}{m}\frac{\partial f^{(2)}}{\partial v}\right).$		(47)

Requiring that $f^{(0)}$ vanishes for $|v|\rightarrow\infty$ , the first equation gives

f^{(0)}=s(x,t)e^{-\frac{mv^{2}}{2T(x)}},

(48)

for some function, $s(x,t),$ to be determined. The integral of the right-hand side of Eq. (46) is zero:

\int_{-\infty}^{+\infty}dv\left(vf^{(1)}+\frac{T}{m}\frac{\partial f^{(1)}}{\partial v}\right)=0

(49)

and therefore the integral on the left-hand side should also be zero:

\int_{-\infty}^{+\infty}dv\,\frac{\partial f^{(0)}}{\partial t}+\int_{-\infty}^{+\infty}dv\,v\frac{\partial f^{(0)}}{\partial x}-\int_{-\infty}^{+\infty}dv\,\Phi^{\prime}\frac{\partial f^{(0)}}{\partial v}=0

(50)

Substituting Equation 48 we have

	$\displaystyle\int_{-\infty}^{+\infty}dv\,v\frac{\partial f^{(0)}}{\partial x}=s^{\prime}\int_{-\infty}^{+\infty}dv\,ve^{-\frac{mv^{2}}{2T}}+s\frac{mT^{\prime}}{2T^{2}}\int_{-\infty}^{+\infty}dv\,v^{3}e^{-\frac{mv^{2}}{2T}}=0$
	$\displaystyle\int_{-\infty}^{+\infty}dv\,\Phi^{\prime}\frac{\partial f^{(0)}}{\partial v}=s\Phi^{\prime}\int_{-\infty}^{+\infty}dv\,\frac{\partial}{dv}\left(e^{-\frac{mv^{2}}{2T}}\right)=0$

and therefore Equation 50 gives the integrability condition:

\frac{\partial f^{(0)}}{\partial t}=0,

(51)

which results finally in

f^{(0)}=f^{(0)}(x,v)=s(x)e^{-\frac{mv^{2}}{2T(x)}}.

(52)

We substitute this in Equation 46 and get

	$\displaystyle\varepsilon e^{-\frac{mv^{2}}{2T}}$	$\displaystyle\left(vs^{\prime}+\frac{v}{T}m\Phi^{\prime}s+\frac{v^{3}}{2T^{2}}mT^{\prime}s\right)=$
		$\displaystyle=\frac{\partial}{\partial v}\left(vf^{(1)}+\frac{T}{m}\frac{\partial f^{(1)}}{\partial v}\right)=\frac{T}{m}\frac{\partial}{\partial v}\left(e^{-\frac{mv^{2}}{2T}}\frac{\partial}{\partial v}e^{\frac{mv^{2}}{2T}}f^{(1)}\right).$		(53)

The ansatz

f^{(1)}(x,v,t)=\left(vh(x,t)+v^{3}q(x,t)\right)e^{-\frac{mv^{2}}{2T}}

(54)

gives

	$\displaystyle h=h(x)=-\frac{(T(x)s(x))^{\prime}+m\Phi(x)^{\prime}s(x)}{T(x)}\varepsilon(x)$		(55)
	$\displaystyle q=q(x)=-\frac{mT(x)^{\prime}}{6T(x)^{2}}\varepsilon(x)s(x).$		(56)

Both $h$ and $q$ do not depend on time. The general solution may be obtained by adding a solution of the homogeneous problem

f^{(1)}(x,v,t)=\varepsilon\left(-v\frac{T^{\prime}}{T}s-vs^{\prime}-v\frac{m\Phi^{\prime}}{T}s-v^{3}\frac{mT^{\prime}}{6T^{2}}s+w(t)\right)e^{-\frac{mv^{2}}{2T}},

(57)

where $\varepsilon$ , $T$ , $\Phi$ , $s$ depend only on $x$ . We have

p(x,t)\equiv\int_{-\infty}^{+\infty}dv\,f(x,v,t)=\sqrt{2\pi\frac{T(x)}{m}}\left(s(x)+\varepsilon(x)w(t)\right)+\mathcal{O}(\varepsilon^{2}),

(58)

where $p$ is the probability density in space. For $N$ Brownian particles, $n=Np$ can be identified as their average number density.

Equation (47) gives the integrability condition

\int_{-\infty}^{+\infty}dv\,\frac{\partial f^{(1)}}{\partial t}=-\frac{\partial}{\partial x}\int_{-\infty}^{+\infty}dv\,vf^{(1)}.

(59)

Since $\partial f^{(0)}/\partial t=0,$ we have

$\displaystyle\frac{\partial p(x,t)}{\partial t}$	$\displaystyle=-\frac{\partial}{\partial x}\int_{-\infty}^{+\infty}dv\,vf^{(1)}+\mathcal{O}(\varepsilon^{2})$
	$\displaystyle=\frac{\partial}{\partial x}\left\{\varepsilon\left(\frac{\partial}{\partial x}(\frac{T}{m}\sqrt{2\pi\frac{T}{m}}(s+\varepsilon w))+\Phi^{\prime}\sqrt{2\pi\frac{T}{m}}(s+\varepsilon w)\right)\right\}$
	$\displaystyle+\mathcal{O}(\varepsilon^{2}),$	(60)

which gives, by substitution of (58), to first order in $\varepsilon,$ the diffusion equation

\frac{\partial p(x,t)}{\partial t}=\frac{\partial}{\partial x}\left\{\mu(x)\left(m\Phi(x)^{\prime}p(x,t)+\frac{\partial}{\partial x}(T(x)p(x,t))\right)\right\},

(61)

where we denote $\mu(x)=\varepsilon(x)/m=1/m\eta(x)$ . This is the Van Kampen inhomogeneous diffusion equation that we suggest applies to gravitational Brownian motion.