Intuition for the definition of the Gamma function?

The geometric approach works.

Let’s compute the volume of the $2n$ dimensional ball, $D^{2n}$, in two ways. One way is extremely clever but has been known for centuries and provides interesting insights: it’s based on Liouville’s trick. Specifically, we will compute two integrals in polar coordinates, one of which is the volume of the ball and the other of which reduces to a product of one-dimensional integrals. Both integrands will depend (at most) on the radial coordinate $r$, which lets us separate out the surface area of the boundary of the ball as a common factor. Write this surface area as $S_{2n-1}$.

There’s essentially just one way to do this trick: integrate $\exp(-r^2)$. Its integral over $\mathbb{R}^{2n}$ equals

$$S_{2n-1} \int_0^\infty {\exp\left(- r^2 \right) r^{2n-1} dr}.$$

However, because $r^2 = x_1^2 + x_2^2 + \ldots + x_{2n}^2$, the integrand (in Cartesian coordinates $\left( x_1, x_2, \ldots, x_{2n} \right)$) factors as $\exp\left(-r^2 \right) = \exp\left(-x_1^2 \right) \cdots \exp\left(-x_{2n}^2 \right)$, each of which must be integrated from $-\infty$ to $+\infty$. Whence

$$S_{2n-1} \int_0^\infty {\exp \left(- r^2 \right) r^{2n-1} dr} = \left( \int_{- \infty}^ \infty {\exp \left( -x^2 \right) dx} \right) ^{2n}.$$

I will call the left hand integral $\tfrac{1}{2} \Gamma \left(n \right)$, because that is what it is (as a simple change of variables shows). In the same notation, $\Gamma \left(1/2 \right) = \int_{-\infty}^\infty {\exp\left(-x^2 \right) dx}$. Algebraic re-arrangement of the foregoing yields the volume of $D^{2n}$ as

$$|D^{2n} | = S_{2n - 1} \int_0^1 {r^{2n - 1} dr} = \frac{{S_{2n - 1} }} {{2n}} = \frac {\Gamma \left(1/2 \right)^{2n}} { n \Gamma \left(n \right) }.$$

That was the first method: the result is a familiar one, but has been left expressed in a way that better reveals its origins in polar and Cartesian integration.

The next way to compute the ball's volume is, I believe, new. It is inspired by Cavalieri’s Principle: the idea that you can shift slices of a solid around without changing the volume of the solid. The generalization is to move two-dimensional slices around and to change their shapes while you do so, but in a way that does not change their areas. It follows that the new solid has the same (hypervolume) as the original, although it might have a completely different shape.

We will compute the volume of a region $Q_n$ in $\mathbb{R}^{2n}$. It is conveniently described by identifying $\mathbb{R}^{2n}$ with $\mathbb{C}^{n}$, using coordinates $z_i = \left( x_{2i - 1}, x_{2i} \right)$, in terms of which

$$Q_n = \{ \mathbf{z} \in \mathbb{C}^n :0 \leq \left| {z_1 } \right| \leq \left| {z_2 } \right| \leq \cdots \leq \left| {z_n } \right| \leq 1 \}.$$

If these were real variables, we could make the volume-preserving transformation $w_1 = z_1, w_2 = z_2 – z_1, \ldots , w_i = z_i – z_{i-1}, \ldots, w_n = z_n – z_{n-1}$, with the sole restriction that the sum of the $w_i$ (all of which are nonnegative) not exceed 1. Because they are complex variables, though, we have to consider the area of an annulus bounded by $z_{i-1}$ and $z_i$: it is proportional to $z_i^2 – z_{i – 1}^2$. The circle of the same area has radius $w_i$ for which $w_i^2 = z_i^2 – z_{i – 1}^2$. Therefore, if we define new variables $w_i$ according to this formula, we obtain a new region- - one of substantially different shape- - having the same volume. This region is defined by $\left| {w_1 }^2 \right| + \cdots + \left| {w_n }^2 \right| \le 1$: that is, it’s our old friend $D^{2n}$. Therefore, the volume of $Q_n$ equals the volume of $D^{2n}$ .

Now for the punch line: $Q_n$ is a fundamental domain for the action of $S[n]$, the symmetric group, on the product of $n$ disks $T^{2n} = \left( D^2 \right) ^n$; $S[n]$ acts by permuting the Complex coordinates $z_1, \ldots, z_n$. The volume of $T^{2n}$ equals $|D^2|^n = \pi ^n$. Writing $|S[n]|$ for the number of permutations and equating our two completely different calculations of the volume of the $2n$ ball gives

$$\pi ^ n / |S[n]| = \frac {\Gamma \left(1/2 \right)^{2n}} { n \Gamma \left(n \right) },$$

whence

$$|S[n]| = \frac{{\pi ^n n\Gamma \left( n \right)}}{{\Gamma \left( {1/2} \right)^{2n} }}.$$

This simplifies: the volume formula for $n = 2$ must give the area of the unit circle, equal to $\pi$, whence $\Gamma \left( 1/2 \right)^2 = \pi$. Finally, then,

$$|S[n]| = n\Gamma \left( n \right).$$

I will finish by remarking that Liouville's method is a perfectly natural thing to encounter when working with the multivariate Normal distribution, so it's not really an isolated trick, but is rather a pretty basic result expressing a defining property of Normal (Gaussian) variates. There are, of course, many other ways to compute the volume of $D^{2n}$, but this one gives us the Gamma function directly.

Geometric approach

Note that $\frac{t^n}{n!}$ is the volume of the set $S_t=\{(t_1,t_2,\dots,t_n)\in\mathbb R^{n}\mid t_i\geq 0\text{ and } t_1+t_2+\cdots+t_n\leq t\}$.

So we can perform a change of variables in the integral, replacing $t=t_1+\dots + t_{n+1}$, as:

$$\begin{align}\int_{0}^\infty \frac{t^n}{n!}e^{-t}\,dt &= \int_{0}^{\infty}\left(\int_{(t_i)\in S_t}\,dt_1\dots dt_n\right)e^{-t}\,dt \\&= \int_{t_1,t_2,\dots,t_n,t_{n+1}} e^{-(t_1+\dots +t_{n+1})}\,dt_1\dots dt_{n+1} \end{align}$$

Where all variables go from $0$ to $\infty$. But this can clearly be factored as the $$\left(\int_{0}^\infty e^{-t}\,dt\right)^{n+1}=1$$

So, what's happening is really in $n+1$ dimensions - the linear map $\mathbb R^{n+1}\to\mathbb R^{n+1}$ defined as $$(t_1,\dots,t_n,t_{n+1})\mapsto\left(t_1,\dots,t_n,\sum_{i=1}^{n+1}t_i\right)$$ preserves volumes.

And then using the property $e^{x+y}=e^xe^y$.

---

Moment generating function approach

The exponential random variable $T$ with $P(T<t)=1-e^{-t}$ for $t\geq 0$ has:

$$E(T^n)=\int_{0}^{\infty} t^ne^{-t}\,dt$$

So we get the moment generating function: $$M_T(x) = \sum \frac{E(T^n)}{n!}x^n$$ But we also have:

$$M_T(x) = E(e^{Tx}) = \int_0^{\infty} e^{tx}e^{-t}\,dt = \int_0^{\infty} e^{-t(1-x)}\,dt = \frac{1}{1-x}$$

This means that $E(T^n)=n!$, which is the result you want.

I suppose this could be written without reference to probability as: If $I_n$ is your integral, then, since everything is positive and convergent, then for $0<x<1$ you can show:

$$\sum \frac{I_n}{n!}x^n = \int_0^{\infty} \left(\sum \frac{(xt)^n}{n!}\right)e^{-t}\,dt = \int_0^{\infty} e^{-t(1-x)}\,dt = \frac{1}{1-x}$$

Sorry to "revive," but here's something I noticed while writing 2012 Fall OMO Problem 25. IMO it gives a kind of neat perspective on the gamma function (as a particular case of a "continuous" generating function), so hopefully this is not too off-topic. It may be somewhat related to George Lowther's answer above, but I don't have the background to fully understand/appreciate his post. Also, there might be a bit of seemingly irrelevant setup here.

Anyways, first consider the following "discrete" problem:

Find the number of integer solutions to $a+b+c+d=18$ with $0\le a,b,c,d\le 6$.

This is fairly standard via PIE or generating functions: for a decent discussion, see this AoPS thread.

Now consider a close "continuous" variant:

Let $N$ be a positive integer and $M$ be a positive real number. Find the probability that $y_1,y_2,\ldots,y_{N-1}\le M$ for a random solution in nonnegative reals to $y_0+y_1+\cdots+y_N=1$.

The direct generalization would have $y_0,y_N\le M$ as well, but I'll keep it like this since the OMO problem uses this version (and both versions allow essentially the same solution).

It's easy to generalize the PIE solution to the discrete version, but here's what I got when I tried to generalize the generating functions solution (unfortunately it's not really rigorous, but I feel like it should all be correct).

To extend the discrete generating functions solution, suppose we work with formal power "integrals" indexed by nonnegative reals instead of formal power series (indexed by nonnegative integers). As a basic example, the integral $\int_0^\infty x^t\; dt$ has a $1$ coefficient for each $x^t$ term, and thus corresponds to $y_0,y_N$---the discrete analog would be something like $\sum_{t=0}^{\infty}x^t$. For $y_1,\ldots,y_{N_1}$, which are bounded above by $M$, we instead have $\int_0^M x^t\; dt$. By "convolution" the desired probability is then $$\frac{[x^1](\int_0^M x^t\; dt)^{N-1}(\int_0^\infty x^t\; dt)^2}{[x^1](\int_0^\infty x^t\; dt)^{N+1}}.$$ But $\int_0^M x^t\; dt = \frac{x^M - 1}{\ln{x}}$ and for positive integers $L$, $(-\ln{x})^{-L} = \frac{1}{(L-1)!}\int_0^\infty x^t t^{L-1}\; dt$. Note that this is essentially the gamma function when $x=e^{-1}$, and for $L=1$ we have $\int_0^\infty x^t\; dt = (-\ln{x})^{-1}$. This is not hard to prove by differentiating (w.r.t. $x$) under the integral sign, but for integer $L$ there's a simple combinatorial proof of $$\frac{1}{(L-1)!}\int_0^\infty x^t t^{L-1}\; dt = (\int_0^\infty x^t\; dt)^L$$ which may provide a bit of intuition for the gamma function. Indeed, suppose we choose $L-1$ points at random from the interval $[0,t]$; the resulting sequence is nondecreasing with probability $\frac{1}{(L-1)!}$. These $L-1$ "dividers" split the interval $[0,t]$ into $L$ nonnegative reals adding up to $t$, which corresponds to the convolution on the RHS. It's not hard to check that this yields a bijection, which explains the $[x^t]$ coefficient of the LHS ($\frac{t^{L-1}}{(L-1)!}$). (Compare with the classical "balls and urns" or "stars and bars" argument for the discrete analog. In fact, it shouldn't be hard to interpret the $[x^t]$ coefficient as a limiting case of the discrete version.)

Now back to the problem: by "convolution" again (but used in a substantive way this time) we find $$\begin{align*} [x^1]\left(\int_0^M x^t\; dt\right)^{N-1}\left(\int_0^\infty x^t\; dt\right)^2 &=[x^1](1-x^M)^{N-1}(-\ln{x})^{-(N-1)}(-\ln{x})^{-2} \\ &=[x^1](1-x^M)^{N-1}\frac{1}{N!}\int_0^\infty x^t t^N\; dt \\ &=\frac{1}{N!}\sum_{0\le k\le N-1,\frac{1}{M}}(-1)^k\binom{N-1}{k}(1-kM)^N \end{align*}$$ and $$[x^1](\int_0^\infty x^t\; dt)^{N+1} =[x^1](-\ln{x})^{-(N+1)} =[x^1]\frac{1}{N!}\int_0^\infty x^t t^N\; dt =\frac{1}{N!}.$$ The desired probability thus comes out to $\sum_{0\le k\le N-1,\frac{1}{M}}(-1)^k\binom{N-1}{k}(1-kM)^N$.

You can see how this explicitly applies to the OMO problem in my post here. Everything I wrote above is based off of the "Official 'Continuous' Generating Function Solution" there.

You are really asking for a direct connection between some property of $n!$ and the integral. This can be done from the recursion

$S[n] = nS[n-1], S[0] = 1$

where $S[n]$ is the order of the symmetric group on $n$ elements. The exponential generating function for this series equals $1/(1-t)$. As in your favorite proof, replace that by $\int_{0}^{\infty} e^{(x-1)t} dt$, expand $e^{xt}$ as a series, reverse the order of summation and integration, and you recover a power series in $t$ whose terms are precisely $\int_{0}^{\infty} t^n e^{-t} dt/n!$. The result follows upon term-by-term comparison.

I realize that's not a whole lot more insight, but it does show explicitly a connection between a defining property of $n!$ and the integral for the $\Gamma$ function.

As far as your subsequent musings go: the integrand is the probability density of a Gamma variate, of course! One useful relation is that the Gamma distribution with parameter $n$ is the sum (i.e., convolution) of $n$ independent exponential variates (i.e., Gamma variates with $n=1$). The expectation of an exponential variate is 1 (this is easy), whence the expectation of a Gamma variate must be $n$ (because expectations add), strongly suggesting the mode of its pdf should be near $n$ (as justified for large $n$ by the Central Limit Theorem).

This is another answer in terms of Poisson processes and the Gamma distribution, and it still uses a bit of calculus which you might call a trick, but I think at least it does add another bit of intuition:

Consider the homogeneous Poisson process with rate parameter 1; this means we are counting the number of occurrences of an event that happens with rate 1. Let's calculate the probabilities $p_k(t)$ that we are in state $k$ at time $t$, i.e., that the event occurs $k$ times in the interval $[0,t]$.

Since the event happens with rate 1, probability mass flows from $p_k(t)$ to $p_{k+1}(t)$ with rate 1. This means that ${p'}_0(t) = -p_0(t)$ and ${p'}_{k+1}(t) = p_k(t) - p_{k+1}(t)$. Also, $p_0(0) = 1$ and $p_{k+1}(0) = 0$. Here comes the bit of calculus: these equations have the solution $p_k(t) = \frac{t^k}{k!} e^{-t}$. (A bit fuzzily, we can read this in two parts: $\frac{t^k}{k!}$ is 1 integrated $k$ times, and $e^{-t}$ represents probability mass being lost at rate 1 to states further down the line. See below for yet another fuzzy explanation.)

Now consider the waiting time $T_k$ until the $k$'th occurrence. Clearly, $T_k = t$ means that the transition from state $k-1$ to state $k$ happens at time $t$, so the probability of $T_k \le t$ is the probability that the transition happens before time $t$, and the density is the derivative of this, i.e., the rate at which probability mass flows from state $k-1$ to $k$. This equals the occurrence rate (i.e., 1) times $p_{k-1}(t)$.

So the probability density of the random variable $T_k$ is $1 \cdot p_{k-1}(t) = \frac{t^{k-1}}{(k-1)!} e^{-t}$ (for $t \ge 0$). Since the probability that there is no occurrence ever is obviously zero, $\int_{0}^{\infty} \frac{t^{k-1}}{(k-1)!} e^{-t} dt = 1$.

Incidentally, this is related to a way of thinking about why $\sum_{k=0}^\infty \frac{t^k}{k!} = \lim_{n\to\infty} (1 + \frac{t}{n})^n$. Suppose you start with one unit of money in an account and get 100% interest, continuously compounded. However, the interest from the original account (number 0) is paid not the original account, but to account #1; interest from account #1 is paid to account #2, and so on.

Then the money $m_k(t)$ in account $k$ equals $1$ integrated $k$ times, and the total money $m(t)$ is $\sum_{k=0}^\infty m_k(t) = \sum_{k=0}^\infty \frac{t^k}{k!}$. But on the other hand, $m(t)$ is continuously compounded at 100% interest, so $m(t) = \lim_{n\to\infty} (1 + \frac{t}{n})^n$ by the usual reasoning.

This gives another fuzzy argument why we should have $p_k(t) = \frac{t^k}{k!} e^{-t}$. The change in the $p_k(t)$ over time consists of two parts: on the other hand, each $p_{k+1}(t)$ increases at rate $p_k(t)$; on the other hand, each $p_k(t)$ decreases at rate $p_k(t)$. If we view the $p_k(t)$ as accounts, since we are taking money out of every account at the constant rate 1, the effect is to decrease the total amount at the constant rate 1, i.e. by a factor of $e^{-t}$, which cancels out the increase of $e^t$ due to the accruing interest. It makes some intuitive sense that we can model this effect by simply rescaling the amount of money in each of the accounts by $e^{-t}$.

The Taylor remainder in integral form is a general framework that relates many important identities and ideas. The relation to $\mathrm{B}$ and $\Gamma$ will be clear as follows.

Taylor expand $f\colon\mathbf R\to \mathbf R$ around $0$ with explicit remainder. The expansion reads $$f(x)=\sum_{k=0}^n\frac{x^k}{k!}f^{(k)}(0)+\int_0^x f^{(n+1)}(x-t)\frac{t^n}{n!}\,dt.$$

Consider expanding $x^m$. It is not hard to see that $\dfrac{d^{n+1}}{dx^{n+1}}(x^m)=\dfrac{m!}{(m-n-1)!}x^{m-n-1}$. The crucial property of $f$ that makes it so special to Taylor expand is that $$f^{(k)}(0)=m!\delta_{k, m}=\begin{cases}m!,\, k=m \\ 0, \,\,\,\,\text{otherwise}\end{cases}$$

The Taylor sum is simply zero if $m\gt n$, so $$x^{m} = \frac{m!}{n! (m-n-1)!} \int_0^x t^{m-n-1}(x-t)^n\, dt. \tag{$m\gt n$}$$ Setting $x=1$ and mapping $m-n-1\mapsto p$ immediately gives the beta integral: $$\frac{n!\,p!}{(n+p+1)!}=\int_0^1 t^p(1-t)^n\, dt. \tag{$p, n \gt 0$}$$

---

To find the gamma function, we will need to recognize the Gauss factorial product for $p$ natural, a fairly intuitive one. Manipulating the expansion of $x^{p+n+1}$ yields $$x^{p+n+1}=\frac{(p+n+1)!}{p!\, n!}\int_0^x t^p(x-t)^n\, dt=\frac{(p+n+1)!}{p!\, n!}x^n \int_0^x t^p \left(1-\frac tx\right)^n\, dt,$$ so the finite version of the gamma function is $$\int_0^x t^p \left(1-\frac tx\right)^n\, dt = \left(\frac{x^{p+1}n!}{(p+n+1)!}\right) p!. \tag{$\star$}$$

Taking limits as $x, n \to \infty$, assuring $x/n \to 1$ and recognizing the Gauss product gives the well known gamma function: $$\int_0^\infty t^p e^{-t}\, dt=p!.$$

This is ready to be generalized: we didn't prove it for real $p$, but we can use the standard investigation to extend it from now on.

---

A number of other concepts are also strongly related. To find Cauchy formula for repeated integration, we map $n+1\mapsto 0$. We are left with $$I^{n+1}f(x)=\int_0^x f(t) \frac{(x-t)^n}{n!}\,dt,$$ the usual formula. The Taylor expansion is again zero, but by induction: $If(0)=0$ and $I^kf(0)\implies I^{k+1}f(0)=0.$ We can now describe the Laplace transform and its limit version.

---

From Cauchy, $$\int_0^x f(t) \left(1-\frac tx\right)^n\,dt= \frac{(I^{n+1}f)(x)\,n!}{x^n}.$$ We want to investigate the limiting behaviour as $n\to \infty$ and $x\approx n/s$. Using limits immediately yields the Laplace transform in a limit fashion: $$\mathcal Lf(s)=\int_0^\infty f(t)e^{-st}\, dt=\lim_{n\to\infty}\frac{s^n\, n!(I^{n+1}f)(n/s) }{n^n}.$$ Unfortunately, this still doesn't make the inverse transform obvious, for example. It seems all these processes have an interpretation in terms of approximating a function by Taylor, I don't know how exactly yet.

Here is an elementary observation,

For any integer $ n \geq 1$ and a given real number $\mu > 0$ (this $\mu$ would be the parameter of a Poisson random variable $X$.)

Then by integration by parts, one can obtain the following identity:

$$ \int_0^\mu e^{-x} \cdot x^n \, dx = n! - e^{-\mu} \cdot n! \left( 1 + \mu + \frac{\mu^2}{2!} + \ldots + \frac{\mu^n}{n!} \,\, \right) \quad(\star) $$

The left side is a Riemann integral while the right side can be written as the following $$n! \cdot \sum_{k = n +1}^\infty e^{-\mu} \frac{\mu^k}{k!} \quad \text{which is equal to } \quad n! \cdot \mathbb{P}\big( X \geq n+1 \big) $$

When the parameter $\mu$ goes to infinity, $X$ becomes almost surely infinite and the improper integral $\int_0^\infty e^{-x} \cdot x^n \, dx $ converges, that is, $$\lim_{\mu\to+\infty} \int_0^\mu e^{-x} \cdot x^n \, dx = \lim_{\mu\to+\infty} n! \cdot \mathbb{P}\big( X \geq n+1 \big) = n! \,\, .$$

(the last limit-equality is derived better from ($\star$), since $n$ is fixed.)

I learned the identity $(\star)$ from some book, if I remember correct, may be Hardy's a course of pure math.

The identity

$${ \int _0 ^{\infty} e ^{- t} t ^n \, dt = n! }$$

naturally arises when finding the density of the random variable ${ Z _n = T _1 + \ldots + T _n }$ where ${ T _i \overset{\text{i.i.d}}{\sim} \text{Exp}(\lambda) . }$

Consider radioactive decay of an element, with memoryless waiting time ${ T \sim \text{Exp}(\lambda) . }$ Let ${ n \in \mathbb{Z} _{> 0} . }$ Consider the time taken for ${ n }$ decays to occur, namely the random variable ${ Z _n = T _1 + \ldots + T _n }$ where ${ T _i \overset{\text{i.i.d}}{\sim} \text{Exp}(\lambda) . }$

For ${ n = 2 }$:

The density of ${ Z _2 }$ is the convolution

$${ {\begin{align} &\, f _{Z _2} (x) \\ = &\, (f \ast f)(x) \\ = &\, \int _0 ^x f(t) f(x - t) \, dt \\ = &\, \int _0 ^x \lambda e ^{-\lambda t} \lambda e ^{-\lambda (x - t) } \, dt \\ = &\, \lambda ^2 x e ^{- \lambda x } \end{align}} }$$

that is

$${ f _{Z _2} (x) = \lambda ^2 x e ^{- \lambda x } \quad \text{ for } x > 0. }$$

For ${ n = 3 }$:

The density of ${ Z _3 }$ is the convolution

$${ {\begin{align} &\, f _{Z _3} (x) \\ = &\, (f \ast f _{Z _2})(x) \\ = &\, \int _0 ^x f(t) f _{Z _2} (x - t) \, dt \\ = &\, \int _0 ^x \lambda e ^{-\lambda t} \lambda ^2 (x - t) e ^{-\lambda (x - t) } \, dt \\ = &\, \frac{\lambda ^3 x ^2}{2} e ^{- \lambda x } \end{align}} }$$

that is

$${ f _{Z _3} (x) = \frac{\lambda ^3 x ^2}{2} e ^{- \lambda x } \quad \text{ for } x > 0 . }$$

Now by induction hypothesis

$${ \text{Induction hypothesis :} \quad f _{Z _k} (x) = \lambda \frac{(\lambda x) ^{k - 1}}{(k-1)!} e ^{- \lambda x } \quad \text{ for } x > 0 }$$

we see

$${ {\begin{align} &\, f _{Z _{k+1}} (x) \\ = &\, (f \ast f _{Z _k})(x) \\ = &\, \int _0 ^x f(t) f _{Z _k} (x - t) \, dt \\ = &\, \int _0 ^x \lambda e ^{-\lambda t} \lambda \frac{(\lambda (x - t)) ^{k - 1}}{(k-1)!} e ^{-\lambda (x - t) } \, dt \\ = &\, \lambda \frac{(\lambda x) ^{k}}{k!} e ^{- \lambda x } \end{align}} }$$

that is

$${ f _{Z _{k + 1}} (x) = \lambda \frac{(\lambda x) ^k}{k!} e ^{- \lambda x } \quad \text{ for } x > 0 . }$$

Hence

$${ \boxed{ f _{Z _n} (x) = \lambda \frac{(\lambda x) ^{n - 1}}{(n-1)!} e ^{- \lambda x } \quad \text{ for } x > 0 } }$$

for all ${ n \in \mathbb{Z} _{> 0} . }$

Especially

$${ \int _0 ^{\infty} \lambda \frac{(\lambda x) ^{n - 1}}{(n-1)!} e ^{- \lambda x } \, dx = 1 }$$

that is

$${ \boxed{ \int _0 ^{\infty} x ^{n - 1} e ^{- x} \, dx = (n - 1)! } . }$$

---

Aside: We can also find the density of ${ Z _n }$ in an alternate way. This is from Introduction to mathematical statistics by Hogg, McKean, Craig, and in the context of the answer here.

Note that the time for ${ n }$ decays is

$${ Z _n = \inf \lbrace \tau > 0 : N(\tau) = n \rbrace . }$$

Its cdf is given by

$${ {\begin{aligned} &\, F _{Z _n} (z) \\ = &\, \mathbb{P}(Z _n \leq z) \\ = &\, 1 - \mathbb{P}(Z _n > z) \\ = &\, 1 - \mathbb{P}(\text{There are less than } n \text{ decays in } [0, z]) \\ = &\, 1 - \sum _{k = 0} ^{n - 1} \mathbb{P}(N(z) = k) \\ = &\, 1 - \sum _{ k = 0} ^{n - 1} e ^{- \lambda z } \frac{(\lambda z) ^k}{k!}. \end{aligned}} }$$

From the already computed density

$${ {\begin{aligned} &\, F _{Z _n} (z) \\ = &\, 1 - \int _z ^{\infty} \lambda \frac{(\lambda x) ^{n - 1}}{(n-1)!} e ^{- \lambda x } \, dx \\= &\, 1 - \int _{\lambda z} ^{\infty} \frac{t ^{n - 1}}{(n - 1)!} e ^{- t} \, dt . \end{aligned}} }$$

Equating both the expressions we have

$${ \boxed{ \mathbb{P}(Z _n > z) = \int _{\lambda z} ^{\infty} \frac{t ^{n - 1} e ^{- t} }{(n - 1)!} \, dt = \sum _{ k = 0} ^{n - 1} \frac{(\lambda z) ^k e ^{- \lambda z}}{k!}} . }$$

So roughly speaking, for

$${ \frac{t ^j e ^{- t}}{j!} , }$$

setting ${ j = (n-1) }$ and integrating from ${ t = \lambda z }$ to infinity is same as setting ${ t = \lambda z }$ and adding from ${ j = 0 }$ to ${ (n - 1) . }$