Let me present you a derivation which does not use the circular assumption presented in the "dart" proof and uses only the property of the Central limit theorem - i.e. that for any suitable set of i.i.d. variables $\{ X_1,X_2,\ldots,X_n \}$ with $\mathrm{E}\left[X\right]=\mu$ and $\mathrm{Var}\left[X\right]=\sigma^2$, we have

$$Z_n := \frac{\bar{X}-\mu}{\sigma/\sqrt{n}} \overset{n\rightarrow\infty}{\longrightarrow} Z \sim N(0,1)$$

where $N(0,1)$ is the so called **normal distribution** (of the variable $Z$). Central limit theorem provides us the **universality**, i.e. no matter which set we choose, we will allways get $Z\sim N(0,1)$ - the same distribution. From this universality only we will find the formula for the distribution function $f(z)$ of the standard normally distributed $Z$. Any other normal variable $X \sim N(\mu,\sigma)$ is get from $Z$ by scalling and its density function must be $f\left((x-\mu)/\sigma\right)/\sigma$.

**1) Stirling approximation**

Without proof (can be found elsewhere) we state that for large $n$ we have

$$n! \approx \sqrt{2\pi n}\left(\frac{n}{e}\right)^n$$

**2) Poisson distribution**

The key to derive the normal distribution density function is to choose **some** particular set of i.i.d. for which we can find the density function of its sum (average $\bar{X}$ is that sum divided by $n$), i.e. to choose $X_1$. There are multiple choices for such choice, in many derivation of normal distribution function it is common to choose $X_1 \sim \mathrm{Ber}(p)$ Bernoulli, so the sum $S_n = X_1 + X_2 + \ldots + X_n \sim \mathrm{Bin}(n,p)$ is Binomial. In our approach, we set $X_1 \sim \mathrm{Po}(1)$. Random variable $X$ with a Poisson distribution $\mathrm{Po}(\lambda)$ has the known discrete distribution function:

$$\mathbb{P}\left[X = k\right]=\frac{\lambda^k}{k!}e^{-\lambda}, \qquad k=0,1,2,\ldots$$

Moreover $\mathrm{Var}\left[X\right] = \lambda$. For $X_1 \sim \mathrm{Po}(1)$ we have $\mu = \lambda = 1$ and $\sigma = \sqrt{\lambda} = 1$.

Poisson distribution is the distribution of the number $X$ of uniformly distributed events in a given finite time interval with a given $\lambda = \mathrm{E}\left[X\right]$ average number of events.

We have chosed the Poisson distrubition because of one special property - **aditivity**. One can think of the sum of two identical (same $\lambda$) independent poisonian $X_1 + X_2$ as the number of the same events it an interval with twice the duration of the previous one. Since

$$\mathrm{E}\left[X_1+X_2\right]=\mathrm{E}\left[X_1\right]+\mathrm{E}\left[X_2\right]=1 + 1 = 2$$

we immediately get (using the "events" description of Poisson distribution), that also

$$S_n = X_1 + X_2 + \ldots + X_n \sim \mathrm{Po}(n)$$

This distribution is discrete with again the step equal to $1$. However, on the left hand side of CLT we get

$$Z_n = \frac{\bar{X}-\mu}{\sigma/\sqrt{n}} = \frac{S_n - n}{\sqrt{n}}$$

Clearly, $S_n$ is reduced by $n$ and then divided by $\sqrt{n}$. This division causes the $\mathrm{Po}(n)$ in the limit of $n\rightarrow\infty$ be continuous, since now the step is $1/\sqrt{n}$ and this covers eventually all the real numbers.

Consider now the limit of large $n$. Let $z$ be a real number. We can write $z=(k-n)/\sqrt{n}$ (for large $k$ and $n$ natural, in a sense of being close to original real $z$). Then, using the definition of (continuous) distribution function

$$f(z)\approx\frac{1}{1/\sqrt{n}}\left(\mathbb{P}\left[Z_n < z + \frac{0.5}{\sqrt{n}}\right]-\mathbb{P}\left[Z_n < z - \frac{0.5}{\sqrt{n}}\right]\right)=\sqrt{n}\,\mathbb{P}\left[Z_n = z\right]$$

since there is only one point $z$ in the interval $(z-0.5/\sqrt{n},z+0.5\sqrt{n})$. But since

$$\mathbb{P}\left[Z_n = z\right] = \mathbb{P}\left[\frac{S_n - n}{\sqrt{n}} = \frac{k - n}{\sqrt{n}}\right]=\mathbb{P}\left[S_n = k\right]=\frac{n^k}{k!}e^{-n}$$

In the limit $n\rightarrow\infty$ with $k = n + z\sqrt{n}$ we then have explicitely

$$f(z) = \lim_{n\rightarrow \infty} \sqrt{n}\frac{n^{n + z\sqrt{n}}}{(n + z\sqrt{n})!}e^{-n}$$

Using Stirling's approximation

$$f(z) = \lim_{n\rightarrow \infty} \sqrt{n}\frac{e^{n + z\sqrt{n}}}{\sqrt{2\pi (n + z\sqrt{n})}}\left(\frac{n + z\sqrt{n}}{n}\right)^{-n - z\sqrt{n}}e^{-n} = \frac{1}{\sqrt{2\pi}}\lim_{n\rightarrow \infty} e^{z\sqrt{n}}\left(1+\frac{z}{\sqrt{n}}\right)^{-n - z\sqrt{n}}$$

**3) Exponential function limit**

This can be further simplified using the well known

$$e^x = \lim_{n\rightarrow \infty}\left(1+\frac{x}{n}\right)^n$$

i.e. for example

$$\lim_{n\rightarrow \infty} \left(1+\frac{z}{\sqrt{n}}\right)^{- z\sqrt{n}} = e^{-z^2}$$

therefore

$$f(z) = \frac{e^{-z^2}}{\sqrt{2\pi}}\lim_{n\rightarrow \infty} e^{z\sqrt{n}}\left(1+\frac{z}{\sqrt{n}}\right)^{-n}$$

The last limit is performed by Taylor expansion of logarithm since

$$\lim_{n\rightarrow \infty} e^{z\sqrt{n}}\left(1+\frac{z}{\sqrt{n}}\right)^{-n} = \exp\left(\lim_{n\rightarrow \infty} z\sqrt{n}-n\ln\left(1+\frac{z}{\sqrt{n}}\right)\right) = \exp\left(\lim_{n\rightarrow \infty} z\sqrt{n}-z\sqrt{n}+\frac{z^2}{2}\right)$$

imediately recovering

$$f(z) = \frac{1}{\sqrt{2\pi}} e^{-\frac{z^2}{2}}$$