Density of $|t_1 - t_2|$ where $t_1$ and $t_2$ are iid with $P(t) = \alpha e^{-t\alpha}$ I am trying to answer the following question from my quantum mechanics textbook and my probability theory is admittedly rusty (this is not schoolwork as should be clear from my post history on Phys and Math SE):

The probability density for the decay of a radioactive nucleus is $P(t) = \alpha e^{-t\alpha}$ where $t > 0$ is the (unpredictable) lifetime of the nucleus, and $\alpha^{-1}$ is the mean lifetime for such a decay process. Calculate the probability density for $|t_1 — t_2|$, where $t_1$ and $t_2$ are the lifetimes of two such identical independent nuclei.

My attempt at a solution is below:

We begin by computing the probability that the random variable $z = |{t_1 - t_2}|$ is less than some value (the cumulative distribution function). We have
$$\textrm{P}(z\leq x) = \textrm{P}(|{t_1 - t_2}|\leq x) =  \textrm{P}(-x \leq {t_1 - t_2}\leq x) = \textrm{P}(t_2-x \leq {t_1}\leq t_2+ x)$$
Consider the $t_1 - t_2$ plane and the two lines $t_2={t_1}+x$ and $t_2= t_1 - x$. For $z\leq x$ we require that the event $(t_1,t_2)$ be in the region between these two lines. Then since we are given that the two nuclei are (1) independent (so we can form the joint probability distribution function for the $t_1$ and $t_2$ random variables as the product of the two individual probability distributions) and (2) identical (so the two individual probability distributions are identical), we find that the joint distribution is $f(t_1,t_2) = \alpha^2 e^{-\alpha(t_1 + t_2)}$. Then to find the probability noted above we compute the integral of this distribution over the relevant region (call it $R$)
\begin{gather*}
    \textrm{P}(z\leq x) = \textrm{P}(t_2-x \leq {t_1}\leq t_2+ x) = \int_R dt_1dt_2 \ f(t_1,t_2) 
    \\
    =\int_0^x \int_0^{t_1+x} \alpha^2 e^{-\alpha(t_1 + t_2)} dt_1dt_2+ \int_x^\infty \int_{t_1-x}^{t_1+x} \alpha^2 e^{-\alpha(t_1 + t_2)}dt_1dt_2 
    \\
    = \frac{e^{-3x\alpha}-3e^{-x\alpha}+2}{2} -\frac{e^{-3x\alpha}-e^{-x\alpha}}{2} = 1- e^{-x\alpha}
\end{gather*}
Now to obtain the pdf from the cdf we take the derivative at the given value in the range of the random variable $z$:
$$f_z(x) = \frac{d\textrm{P}(z\leq x)}{dx} = \alpha e^{-x\alpha}$$
which is the same as the distributions for the individual decays.

But I'm sceptical that I've arrived at the conclusion that $z$ has the same distribution as the individual decay time RVs which doesn't seem to match my intuition on physical grounds. Am I doing something wrong here?
 A: It's surprising but correct.  The exponential distribution is memoryless, meaning that the distribution of time until a decay is the same whenever you start.  It's easy to show it's the same if you pick any fixed time to start, and you've shown it's also the same if you pick a random time to start. In particular, if you pick the time nucleus 1 decayed as the start and nucleus 2 hasn't decayed yet it's the same, and if you pick the time nucleus 2 decayed as the start and nucleus 1 hasn't decayed yet it's the same.
A: I would handle with a different approach which, imo, is more insightful:
Consider $T_i\overset{\textrm{iid}}{\sim}\mathrm{Exp}(\alpha), ~i\in\{1, 2\}; Z:=|T_1-T_2|.$ Now
\begin{align}\mathbb P(Z\leq z) &= \mathbb P(|T_1-T_2|\leq z) \\&= \mathbb P(T_1\leq T_2+z)-\mathbb P(T_1\leq T_2-z)\\&=\int_0^\infty \mathbb P(T_1\leq t_2+z) \alpha\exp(-\alpha t_2)~\mathrm dt_2-\int_z^\infty\mathbb P(T_1\leq t_2-z)\alpha\exp(-\alpha t_2)~\mathrm dt_2\\&= \int_0^\infty\alpha(1-\exp(-\alpha(t_2+z)))\exp(-\alpha t_2)~\mathrm dt_2-\int_z^\infty \alpha (1-\exp(-\alpha(t_2-z)))\exp(-\alpha t_2)~\mathrm dt_2\\&= 1-\exp(-\alpha z) \int_0^\infty\alpha\exp(-2\alpha t_2)~\mathrm dt_2-\alpha\left[\int_z^\infty\exp(-\alpha t_2)~\mathrm dt_2-\exp(\alpha z) \int_z^\infty \exp(-2\alpha t_2)~\mathrm dt_2\right]\\&= 1-\frac{\exp(-\alpha z) }2-\alpha\left[\frac{\exp(-\alpha z) }{\alpha}+\frac{\exp(\alpha z-2\alpha z) }{2\alpha}\right]\\&= 1-\exp(-\alpha z),\tag 1\end{align}
which is the cdf of $\mathrm{Exp}(\alpha). $
A: For exponential distribution family, operating its survival function $S(t) = e^{-\alpha t}$ is usually slightly more convenient than treating the distribution function $F(t) = 1 - e^{-\alpha t}$.  Except it, my answer below is almost the same as User1865345's.
For any $t > 0$, by the independence of $T_1$ and $T_2$:
\begin{align}
 & P[|T_1 - T_2| > t] = P[T_1 > T_2 + t] + P[T_2 > T_1 + t] \\
=& \int_0^{\infty}P[T_1 > s + t]f_{T_2}(s)ds + \int_0^\infty P[T_2 > s + t]f_{T_1}(s)ds \\
=& 2\int_0^{\infty}e^{-\alpha(s + t)}\alpha e^{-\alpha s}ds  \\
=& e^{-\alpha t} \int_0^\infty 2\alpha e^{-2\alpha s}ds \\
=& e^{-\alpha t}, 
\end{align}
which is the survival function of an $\exp(\alpha)$ r.v., hence the result.
