A square contains many random points. From each point, a disc grows until it hits another disc. What proportion of the square is covered by the discs?

Question

A square lamina contains $n$ independent uniformly random points. At a given time, each point becomes the centre of a disc whose radius grows from $0$, at say $1$ cm per second, and stops growing when the disc hits another disc.

Consider the discs after all growth has stopped. Here is an example with $n=20$.

Let $P=$ proportion of the square that is covered (i.e. occupied) by the discs. In the example above, $P\approx0.383$.

What does $P$ approach as $n\to\infty$ ?

The shape of the lamina (for example, square) does not matter, since we are taking $n\to\infty$. The rate of growth (for example, $1$ cm per second) does not matter, as long as all the discs start growing at the same time and grow at the same rate.

My thoughts

I tried to find the probability that a new random point in the square lies in one of the existing discs. I also tried to find the average area of a disc. But I haven't succeeded. These questions seem to be complicated by the fact that the size of a point's disc is determined not only by the point's distance to its neighbors, but also its neighbors' distances to their neighbors, and so on.

Context

This question was inspired by another question, "A disc contains $n$ random points. Each point is connected to its nearest neighbor. What does the average cluster size approach as $n\to\infty$?". Both of these questions are about inherent properties of the $2D$ Poisson process.

Answer 1

I did some literature search, and found that:

1. OP's question exactly corresponds to what is called lilypond model introduced by Häggström and Meester (1996). Daley et al. (1999) independently studied the same model and reported a Monte-Carlo estimate of the limit proportion with $$ \text{mean} = 0.3487 \qquad\text{and}\qquad \text{SD} = 0.00045 $$ using $10^6$ samples.

2. It seems that OP's question has not been answered in the literature yet, which is not surprising considering that statistical properties of a geometric functional of a continuum percolation model is notoriously hard to answer, if not impossible.

3. Interestingly, 1D version of OP's question has a definite answer: $$P_{\text{1D,limit}} = 1 - e^{-1}$$ See Section 6 of Daley et al. (1999), for instance.

---

References.

• O. Häggström, and R. Meester. “Nearest Neighbor and Hard Sphere Models in Continuum Percolation.” Random Structures & Algorithms 9, no. 3 (1996): 295–315.

• D. J. Daley, H. Stoyan, and D. Stoyan. “The Volume Fraction of a Poisson Germ Model with Maximally Non-Overlapping Spherical Grains.” Advances in Applied Probability 31, no. 3 (1999): 610–24.

Answer 2

Attached a MATHEMATICA script which generates a set of random circles as needed. The calculation process could be more efficient but it would also be less transparent.

data = {{1.0032, 0.1304},{0.8713, 0.0446}, {0.8458,0.2497}, 
{0.6789, 0.2751}, {0.6137, 0.2179}, 
{0.5437, 0.1209}, {1.0398, 0.5311}, 
{1.0525, 0.6202}, {0.8347, 0.5502}, 
{0.8061, 0.6297}, {0.7202, 0.6170}, 
{0.6487, 0.7219}, {0.6964, 0.7839}, 
{0.8442, 0.9461}, {0.0763, 0.2990}, 
{0.1160, 0.5915}, {0.2941, 0.5979}, 
{0.2305, 0.7871},{0.1717, 0.8078}, 
{0.0556,0.9334}};
boundary = {{0.0270, 0.0250}, {1.0764, 0.0250}, 
{1.0780, 1.0744}, {0.0238,1.0744}, {0.0270,0.0250}};
(**********************************)
SeedRandom[1];
data = RandomReal[{-1, 1}, {200, 2}];
boundary = {{-1, -1}, {1, -1}, {1, 1}, {-1, 1}, {-1, -1}};
(**********************************)
{X, Y} = Transpose[boundary];
xmin = Min[X];
xmax = Max[X];
ymin = Min[Y];
ymax = Max[Y];
n = Length[data];
gr1 = ListPlot[data, PlotStyle -> Black];
gr2 = ListLinePlot[boundary, PlotStyle -> Thick];
circs = Table[{data[[k]], 0, 0}, {k, 1, n}];
For[i = 1, i <= n, i++, 
 For[j = i + 1, j <= n, j++, 
  dist[i, j] = Norm[circs[[i, 1]] - circs[[j, 1]]]]]
minDist[circs_] := Module[{n = Length[circs], d, dmin = Infinity, pi, pj, i, j, i0, j0},
  For[i = 1, i <= n, i++,
   For[j = i + 1, j <= n, j++, 
     If[circs[[i, 3]] == 1 && circs[[j, 3]] == 1, Continue[]];
     d = dist[i, j] - circs[[i, 2]] - circs[[j, 2]];
     If[d <= dmin,
      dmin = d;
      pi = circs[[i]];
      pj = circs[[j]];
      i0 = i;
      j0 = j];
     ];
   ]; Return[{dmin, i0, j0, pi, pj}]
  ]
dmin = 0;
grs = {};
While[dmin < Infinity,
  {dmin, i0, j0, pi, pj} = minDist[circs];
  If[dmin == Infinity, Continue[]];
  If[pi[[3]] == 0 && pj[[3]] == 0,
   circs[[i0, 2]] += dmin/2;
   circs[[j0, 2]] += dmin/2;
   circs[[i0, 3]] = 1;
   circs[[j0, 3]] = 1;
   For[k = 1, k <= n, k++,
    If[circs[[k, 3]] != 1, circs[[k, 2]] += dmin/2]
    ]];
  If[circs[[i0, 3]] == 1 && circs[[j0, 3]] == 0,
   circs[[j0, 2]] += dmin;
   circs[[j0, 3]] = 1;
   For[k = 1, k <= n, k++,
    If[circs[[k, 3]] != 1, circs[[k, 2]] += dmin]]
   ];
  If[circs[[j0, 3]] == 1 && circs[[i0, 3]] == 0,
   circs[[i0, 2]] += dmin;
   circs[[i0, 3]] = 1;
   For[k = 1, k <= n, k++,
    If[circs[[k, 3]] != 1, circs[[k, 2]] += dmin]]
   ];
  grcircs = Table[Graphics[{Blue, Disk[circs[[k, 1]], circs[[k, 2]]]}, PlotRange -> {{xmin, xmax}, {ymin, ymax}}, PlotRangeClipping -> True], {k, 1, n}];
  AppendTo[grs, grcircs];
  ];
Show[grs, gr1, gr2, PlotRange -> {{xmin, xmax}, {ymin, ymax}}, PlotRangeClipping -> True]

The proposed setup.

A random setup

My guess is that the maximum ratio is achieved with a honeycomb structure, going to the limit. From Lagrange $\frac{\pi}{2\sqrt{3}}$.