If I pick any random star in our galaxy, what is the mean & std dev of distances from that star to its, say, 5 nearest neighbors? As this distribution likely depends on the distance from the galactic center, is there a useful function of said radial position that can be applied?

  • 1
    $\begingroup$ The nearest neighbor distribution function isn't (anything like) a normal distribution, so can you say exactly what you want with the standard deviation and mean. $\endgroup$
    – ProfRob
    Jun 6, 2022 at 14:15
  • $\begingroup$ public.nrao.edu/ask/… says the average distance to the nearest star is 5 ly, but I don't know if anyone has specifically studied 5 nearest stars. As @ProfRob notes, this distribution is unlikely to be normal. You might try to find a general distribution from the galactic center and use that to derive what you want or use a catalog like GAIA to compute it yourself. $\endgroup$ Jun 6, 2022 at 14:51
  • $\begingroup$ Our galaxy has "arms" that spiral outwards, so the density is higher in the arms than in between them. So there's not really anything useful to derive from averaging them all together. $\endgroup$ Jun 6, 2022 at 15:51
  • $\begingroup$ The nearest neighbor function is a reasonably well known statistical thing and then you just scale that with the density. At any point, the average distance to the nearest neighbor is $n^{-1/3}$ where $n$ is the number of star per cubic distance unit. $\endgroup$
    – ProfRob
    Jun 6, 2022 at 17:37

1 Answer 1


We can treat this as a Poisson point process: the probability that there are $N$ points inside a volume $V$ with point density is $\rho$ $$\Pr[N|V,\rho]=\frac{(\rho V)^N}{N!}e^{-\rho V}.$$ This assumes independence of where the points are, which to a first approximation is reasonable for star systems.

The probability that there is no star within distance $d$ of the origin if stars have density $\rho$ is $\Pr[N=0|d,\rho]=e^{-4\pi \rho d^3/3}$. This means that the cumulative distribution function of the closest star $F_1(d)=1-\Pr[N=0|d]=1-e^{-4\pi \rho d^3/3}$. If we take the derivative we get the probability density $$f_1(d)=4\pi \rho d^2 e^{-4\pi \rho d^3/3}.$$ The average becomes $\int_0^\infty x f_1(x) dx = [\Gamma(1/3)/(36\pi)^{1/3}]\rho^{-1/3}\approx 0.5540 \rho^{-1/3}.$

The median is simple: $1/2 = F_1(d)$ gives $[3\ln(2)/4\pi\rho]^{1/3}\approx 0.5490 \rho^{-1/3}$.

The variance can be found from $\int_0^\infty x^2 f_1(x) dx$ but the result is a messy expression $V=\left([2\Gamma(2/3)/(48\pi^2)^{1/3}] - [\Gamma(1/3)/(36\pi)^{1/3}]^2\right )\rho^{-2/3} $.

Distribution of <span class=$f_n(d)$ found using Monte Carlo estimation." /> Distribution of $f_n(d), n=1$ to 5 found using Monte Carlo estimation. $\rho=1$, $10^5$ samples.

Calculating the corresponding distributions for n stars can be done using order statistics applied to $f_1(d)$. The expressions become rather complicated, though. I generally favor just doing a Monte Carlo simulation.

Local variation The star density $\rho(x)$ varies in different places, but the variation is small over distances comparable to the nearest 5 stars.

A common model of the Milky Way is a combination of exponential disks $\rho \propto e^{-z/h}e^{-R/r}$ where $h$ is the scale height of the disk (220 pc for the thin disk and 450 pc for the thick disk population) and $r$ a scale length $\sim 2.6$ kpc along the disk. This makes the above means increase exponentially as one moves outward, but the scale factors are far larger than the typical means unless we care about very large $N$. This is also true inside globular and open clusters.

Multiplicity One can make a more complex model taking stellar multiplicity into account. This makes $f_n(d)$ a mixture distribution of $f_i(d), i\leq n$ with weights set by $p_i$, the probability of i-ariness, and the number of integer partitions of n. For $f_2(d)$ this is $$\hat{f}_2(d)= (1-p_1) f_1(d) + p_1 f_2(d)$$ (either the closest star is is multiple and $f_1(d)$ applies, or it is single and we use $f_2(d)$). For $f_3(d)$ things become far more complex: $$\hat{f}_3(d)=(1-p_1-p_2)f_1(d) + (1-p_1)^2f_2(d) + (p_1^2+2p_1p_2)f_3(d)$$ (closest system is trinary, second closest system is binary or higher order while first is at most binary, closest two systems are single or one of them binary). It gets opaque very fast, and numerical simulation is likely more useful.

Monte Carlo simulation of distances to nearby stars, with multiplicities from (Eggleton & Tokovinin 2008)

Monte Carlo simulation of distances to nearby stars, with multiplicities from (Eggleton & Tokovinin 2008)

  • 2
    $\begingroup$ I think this approach requires the assumption that the star locations are locally stochastically independent, which I don't think is justified since stars often form near one another in collapsing gas clouds. Our own Sun is in the minority of stars in the Milky Way in that it isn't in a gravitationally bound orbit with other star(s). How does this answer account for these (binary, trinary, etc...) systems and clumps like globular clusters in the galactic halo? If not, is this approach a way to compute an "upper-bound"? $\endgroup$
    – Connor Garcia
    Jun 6, 2022 at 22:45
  • 1
    $\begingroup$ @ConnorGarcia - Yes, the Poisson process model assumes independence. If we have multiple star systems this doesn't change $f_1(d)$, but the higher order distributions become mixture distributions with weights set by the multiplicity distribution. The clumpiness scale is generally spatially varying more slowly than star distances, so unless you look for $f_n(d)$ for $n\gg 5$ this will not matter, you can use the local mean. $\endgroup$ Jun 7, 2022 at 5:31

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .