<p>We can treat this as a <a href="https://en.wikipedia.org/wiki/Poisson_point_process" rel="nofollow noreferrer">Poisson point process</a>: the probability that there are <span class="math-container">$N$</span> points inside a volume <span class="math-container">$V$</span> with point density is <span class="math-container">$\rho$</span> <span class="math-container">$$\Pr[N|V,\rho]=\frac{(\rho V)^N}{N!}e^{-\rho V}.$$</span> This assumes independence of where the points are, which to a first approximation is reasonable for star systems.</p>
<p>The probability that there is <em>no</em> star within distance <span class="math-container">$d$</span> of the origin if stars have density <span class="math-container">$\rho$</span> is <span class="math-container">$\Pr[N=0|d,\rho]=e^{-4\pi \rho d^3/3}$</span>. This means that the cumulative distribution function of the closest star <span class="math-container">$F_1(d)=1-\Pr[N=0|d]=1-e^{-4\pi \rho d^3/3}$</span>. If we take the derivative we get the probability density <span class="math-container">$$f_1(d)=4\pi \rho d^2 e^{-4\pi \rho d^3/3}.$$</span> The average becomes <span class="math-container">$\int_0^\infty x f_1(x) dx = [\Gamma(1/3)/(36\pi)^{1/3}]\rho^{-1/3}\approx 0.5540 \rho^{-1/3}.$</span></p>
<p>The median is simple: <span class="math-container">$1/2 = F_1(d)$</span> gives <span class="math-container">$[3\ln(2)/4\pi\rho]^{1/3}\approx  0.5490 \rho^{-1/3}$</span>.</p>
<p>The variance can be found from <span class="math-container">$\int_0^\infty x^2 f_1(x) dx$</span> but the result is a messy expression <span class="math-container">$V=\left([2\Gamma(2/3)/(48\pi^2)^{1/3}] - [\Gamma(1/3)/(36\pi)^{1/3}]^2\right )\rho^{-2/3} $</span>.</p>
<p><a href="https://i.sstatic.net/Holp9.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Holp9.png" alt="Distribution of <span class="math-container">$f_n(d)$</span> found using Monte Carlo estimation." /></a>
<em>Distribution of <span class="math-container">$f_n(d), n=1$</span> to 5 found using Monte Carlo estimation. <span class="math-container">$\rho=1$</span>, <span class="math-container">$10^5$</span> samples.</em></p>
<p>Calculating the corresponding distributions for n stars can be done using <a href="https://en.wikipedia.org/wiki/Order_statistic" rel="nofollow noreferrer">order statistics</a> applied to <span class="math-container">$f_1(d)$</span>. The expressions become rather complicated, though. I generally favor just doing a Monte Carlo simulation.</p>
<p><strong>Local variation</strong> The star density <span class="math-container">$\rho(x)$</span> varies in different places, but the variation is small over distances comparable to the nearest 5 stars.</p>
<p>A common model of the Milky Way is a combination of exponential disks <span class="math-container">$\rho \propto e^{-z/h}e^{-R/r}$</span> where <span class="math-container">$h$</span> is the scale height of the disk (220 pc for the thin disk and 450 pc for the thick disk population) and <span class="math-container">$r$</span> a scale length <span class="math-container">$\sim 2.6$</span> 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 <span class="math-container">$N$</span>. This is also true inside globular and open clusters.</p>
<p><strong>Multiplicity</strong>
One can make a more complex model taking stellar multiplicity into account. This makes <span class="math-container">$f_n(d)$</span> a mixture distribution of <span class="math-container">$f_i(d), i\leq n$</span> with weights set by <span class="math-container">$p_i$</span>, the probability of i-ariness, and the number of integer partitions of n. For <span class="math-container">$f_2(d)$</span> this is <span class="math-container">$$\hat{f}_2(d)= (1-p_1) f_1(d) + p_1 f_2(d)$$</span> (either the closest star is is multiple and <span class="math-container">$f_1(d)$</span> applies, or it is single and we use <span class="math-container">$f_2(d)$</span>). For <span class="math-container">$f_3(d)$</span> things become far more complex: <span class="math-container">$$\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)$$</span> (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.</p>
<p><a href="https://i.sstatic.net/CDA33.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/CDA33.png" alt="Monte Carlo simulation of distances to nearby stars, with multiplicities from (Eggleton &amp; Tokovinin 2008)" /></a></p>
<p><em>Monte Carlo simulation of distances to nearby stars, with multiplicities from (<a href="https://academic.oup.com/mnras/article/389/2/869/974100" rel="nofollow noreferrer">Eggleton &amp; Tokovinin 2008</a>)</em></p>

CC BY-SA 4.0
2022-06-06T22:45:20.387Z
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"?
Connor Garcia
CC BY-SA 4.0
2022-06-07T05:31:08.790Z
@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.
Anders Sandberg
CC BY-SA 4.0
2022-06-06T20:36:59.300Z
2022-06-07T08:42:32.177Z
2022-06-07T08:42:32.177Z
Anders Sandberg
<p>If I pick any random star in our galaxy,  what is the mean &amp; 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?</p>

CC BY-SA 4.0
2022-06-06T14:15:58.590Z
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.
ProfRob
CC BY-SA 4.0
2022-06-06T14:51:45.220Z
https://public.nrao.edu/ask/what-is-the-average-distance-between-stars-in-our-galaxy/ 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.
Barry Carter
CC BY-SA 4.0
2022-06-06T15:51:10.017Z
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.
Greg Miller
CC BY-SA 4.0
2022-06-06T17:37:24.497Z
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.
ProfRob
<galaxy><distances>
Distribution function for distance to nearest stars from any given star in our galaxy
Carl Witthoft
# Distribution function for distance to nearest stars from any given star in our galaxy

(Asked by: username_0 on 2022-06-06)

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?

**Tags:** galaxy, distances
> **Comments:**
> **username_1**: 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.  
> **username_2**: https://public.nrao.edu/ask/what-is-the-average-distance-between-stars-in-our-galaxy/ 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.  
> **username_3**: 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.  
> **username_1**: 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.  

---

## 1 answers

### Answer by username_4 on 2022-06-06. Score: 4

We can treat this as a [Poisson point process](https://en.wikipedia.org/wiki/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=](https://i.sstatic.net/Holp9.png)$f\_n(d)$ found using Monte Carlo estimation." />](https://i.sstatic.net/Holp9.png)
*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](https://en.wikipedia.org/wiki/Order_statistic) 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)](https://i.sstatic.net/CDA33.png)](https://i.sstatic.net/CDA33.png)

*Monte Carlo simulation of distances to nearby stars, with multiplicities from ([Eggleton & Tokovinin 2008](https://academic.oup.com/mnras/article/389/2/869/974100))*

> **Comments:**
> **username_5**: 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"?  
> **username_4**: @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.  

---

