Doing Olber's paradox with specific systems is simple enough. Take stars on a grid of length $\ell$, then the total luminosity is just
\begin{equation} L = \sum_{i,j,k} \frac{L_0}{4\pi\ell^2(i^2+j^2+k^2)} \end{equation}
which is enough to show divergent. I have attempted to do it for a more realistic case, with a uniform distribution of stars. For any region of space $\Omega$, there is a probability of having $n$ stars such that
\begin{equation} E[P_\Omega(N_\odot)] = \rho_\odot V(\Omega) \end{equation}
so that we have a constant density of stars.
Now let's consider a split of our universe into shells. The exact shape of the shell doesn't matter too much, I have attempted both a cubical one and a spherical one. In both cases, for a shell at a distance of order $d$ and thickness $\ell$, we have that the volume is roughly $\approx d^2 \ell$, and the luminosity is $\approx d^{-2}$.
What I need to prove then is the bounds of the total luminosity.
\begin{equation} P(L_{\mathrm{tot}} < K) = P(\sum_{d = 1}^\infty L_d < K) \end{equation}
Now via some joint distribution theorems, I obtained that
\begin{equation} P(L_{\mathrm{tot}} < K) \leq \prod_{d = 1}^\infty P(L_d < K) \end{equation}
This may not be the best upper bound, but it is certainly more pleasant to work with. The total luminosity of the shell is defined by the distance $d$ and a random variable $N_d$. There are many little issues with respect to the bounds of the distance, but roughly it is equal to
\begin{equation} P(L_d < K) = P(\frac{L_0}{4\pi \ell^2 d^2} N_d < K) \end{equation}
By the law of large number, we have that
\begin{equation} \lim_{d \to \infty} \frac{N_d}{4\pi \ell d^2} \approx \rho_\odot \end{equation}
or more specifically,
\begin{equation} \lim_{d \to \infty} P(\frac{L_0}{4\pi \ell^2 d^2} N_d - \frac{\rho_\odot}{4\pi \ell d^2} < K - \frac{\rho_\odot}{4\pi \ell d^2}) \end{equation}
The issue here is that this means that, if $K > E[N]$, then the limit will in fact be $1$ and not $0$. If we want our luminosity to be (almost surely) divergent, our upper bound should be zero.
Does this stem from the choice of the joint distribution upper bound selected previously, or is something else at play here? And what would be a better way to get the correct result, ideally without having to perform an infinite chain of convolutions?
Edit : For more details on what I tried :
Consider $\mathbb{R}^3$ split into cubes of sides $\ell$, with the cube $(0,0,0)$ at $(-\ell/2, \ell/2)$. The distance of any star in cube $(i,j,k)$ is in the range
\begin{equation} r \in [\frac{\ell}{\sqrt{2}} \sqrt{i^2 + j^2 + k^2}, \frac{\ell}{\sqrt{2}} \sqrt{(|i| + 1)^2 + (|j| + 1)^2 + (|k| + 1)^2}] \end{equation}
and the luminosity is within
\begin{eqnarray} L_{\min(i,j,k)} &=& \left[ \frac{L}{2\pi \ell^2 r^2_{\max}}, \frac{L}{2\pi \ell^2 r^2_\min} \right] \end{eqnarray}
If we consider a cube of those cubes of side $2d + 1$ (so that we have all the cubes at coordinates $(\pm d, j, k)$ and so forth), the boundary of that cube has a total number of
\begin{equation} (2d + 1)^3 - (2(d-1) + 1)^3 = 24 d^2 + 2 \end{equation}
Any star in that boundary is at a distance
\begin{equation} r \in [\frac{\ell}{\sqrt{2}} d, \sqrt{3}\frac{\ell}{\sqrt{2}} (d+1)] \end{equation}
And the luminosity is within
\begin{eqnarray} L_{d} &=& \left[ \frac{L}{6\pi \ell^2 (d+1)^2}, \frac{L}{2\pi \ell^2 d^2} \right] \end{eqnarray}
Now we compute the probability of the total luminosity :
\begin{eqnarray} P(L_{\mathrm{tot}} < K) &=& P(\sum_{d = 1}^\infty L_{d} < K) \\ \end{eqnarray}
We're summing over the luminosity of every boundaries of cubes of distance $d$ here. The probability of that sum being inferior to $K$ is some integral of the joint probability distribution over $[0,K]^d$, to be precise on the subset where $\sum L < K$. But this will itself be inferior to the integral on the whole domain, and the joint probability is independent, so that
\begin{eqnarray} P(L_{\mathrm{tot}} < K) &\leq& \prod_{d = 1}^\infty P(L_{d} < K) \\ \end{eqnarray}
We can always find upper bounds by considering a random variable stricly superior to the previous one. By considering that every star has the maximal luminosity in each shell, the luminosity of each shell is computed by
\begin{eqnarray} P(L_{\mathrm{tot}} < K) &\leq& \prod_{d = 1}^\infty P(\frac{L}{2\pi \ell^2 d^2} \sum_a^{24 d^2 + 2} N_a < K) \\ \end{eqnarray}
From this, we can attempt to apply the law of large numbers, but we can already see that at $d \to \infty$, this is gonna be $\approx N < K$, which is true for every large enough values of $K$, while we would like the probability to converge to zero so that the product does as well. I suspect the approximation of the joint distribution is at fault here.
Edit : Another attempt :
Take the joint probability distribution
\begin{eqnarray} P(L_{\mathrm{tot}} < K) &=& P(\sum_{d = 1}^\infty \sum_{a}^{24d^2 + 2} L_{a,d} < K)\\ &=& \int_{\sum k_{a,d} < K} \prod_{d = 1}^\infty \prod_{a = 1}^{24d^2 + 2} f(k_{a,d}) dk_{a,d} \end{eqnarray}
Using independence and Holder's inequality,
\begin{eqnarray} P(L_{\mathrm{tot}} < K) &\leq& \lim_{D \to \infty} \prod_{d = 1}^D \prod_{a = 1}^{24d^2 + 2} \left( \int_0^K f^{2^{D(24d^2 + 2)}}(k) dk \right)^{2^{-D(24d^2 + 2)}}\\ &\leq& \lim_{D \to \infty} \prod_{d = 1}^D \prod_{a = 1}^{24d^2 + 2} \left( P(L \leq K) \right)^{2^{-D(24d^2 + 2)}}\\ \end{eqnarray}
Short of assuming that $P(L \leq K) < 1$ for a single cube (which isn't necessarily a great assumption), still can't prove that it converges to $0$.
