The atomic orbitals (wavefunctions) $\psi(r,\theta,\phi)$ are comprised of a radial component $R_{n,l}(r)$, as well as an angular component $Y_{l,m}(\theta,\phi)$. These are obtained by separately solving the radial and angular parts of the Schrödinger equation, the details of which can be found in any text. Often these are referred to as the radial and angular wavefunctions, and the angular wavefunctions are also called spherical harmonics.
As a concrete example, the $2\mathrm p_z$ orbital ($\{n,l,m\} = \{2,1,0\}$) has the radial component $R_{2,1}$ and the angular component $Y_{1,0}$:
$$R_{2,1} = \sqrt{\frac{1}{24}}\left(\frac{Z}{a}\right)^{3/2}\rho\mathrm e^{-\rho/2} \qquad Y_{1,0} = \sqrt{\frac{3}{4\pi}}\cos\theta \tag{1}$$
where $a$ is a collection of constants ($a \equiv 4\pi\varepsilon_0\hbar^2/\mu e^2$), and $\rho$ is related to $r$ by another collection of constants ($\rho \equiv 2Zr/na$). The product of these two expressions is the wavefunction corresponding to the $2\mathrm p_z$ orbital.
The atomic orbitals themselves are normalised in the sense that
$$\int_{\text{all space}} |\psi(r,\theta,\phi)|^2 \mathrm dV = \int_0^{2\pi} \!\!\!\!\int_0^\pi\!\!\!\int_0^\infty |R(r)|^2|Y(\theta,\phi)|^2 r^2\sin\theta\,\mathrm dr\,\mathrm d\theta\,\mathrm d\phi = 1 \tag{2}$$
Conventionally, both the radial and angular wavefunctions are separately normalised. This means that
$$\int_0^\infty |R(r)|^2 r^2\mathrm dr = 1; \qquad \int_0^{2\pi}\!\!\!\!\int_0^\pi |Y(\theta,\phi)|^2 \sin\theta\,\mathrm d\theta\,\mathrm d\phi = 1 \tag{3}$$
(by multiplying these two together, you recover the original normalisation of the entire atomic orbital.)
Normally, since $\psi = \psi(r,\theta,\phi)$, the probability density $|\psi|^2$ depends on all three coordinates. The probability density (technically, $|\psi|^2\,\mathrm dV$) gives you the probability of finding an electron in an infinitesimal cube frustum-like shape,† located at the point $(r,\theta,\phi)$. By summing these $|\psi|^2\,\mathrm dV$ over every single point in space, i.e. by integrating over all three variables, we get the total probability of finding an electron in all space, which is $1$.
However, often we are not interested in the angular dependency, and so what we can do is to sum up all these bits which have one value of $r$, but different values of $\theta$ and $\phi$. That would give us the probability of finding the electron in a tiny spherical shell of radius $r$. Mathematically, we need to integrate over $\theta$ and $\phi$. By virtue of the angular wavefunction being normalised, the integral over $\theta$ and $\phi$ evaluates to $1$, and hence we have
$$|\psi|^2\,\mathrm dV = (|R|^2r^2\,\mathrm dr)\cdot (|Y|^2\sin\theta\,\mathrm d\theta\,\mathrm d\phi) \overset{\text{integrate over }\theta,\phi}{\longrightarrow} |R|^2r^2\,\mathrm dr \tag{4}$$
So, this quantity $|R|^2r^2\,\mathrm dr$ gives us the probability of finding the electron in a spherical shell of radius $r$. The corresponding probability density function is the function which controls how this probability varies with $r$, and we give it a special name, the radial distribution function (RDF); hence Atkins defines the RDF to be equal to $|R|^2r^2$.
This expression is valid regardless of the form of the angular wavefunction, as no matter what it is, when integrating over $\theta$ and $\phi$, it comes out to be $1$. The widely-cited alternative definition $4\pi r^2 |\psi|^2$ is only really relevant for spherically symmetric orbitals (in other words, s-orbitals). In this case, we have $\{l, m\} = \{0,0\}$ and
$$Y_{0,0} = \sqrt{\frac{1}{4\pi}} \tag{5}$$
(you may check that this is normalised according to equation $(3)$.) Hence
$$|\psi|^2 = |R|^2|Y^2| = \frac{1}{4\pi}|R^2| \tag{6}$$
and it therefore follows that $|R|^2r^2 = 4\pi r^2|\psi|^2$, so both definitions are equivalent in the case of spherically symmetric orbitals (s-orbitals). Notice that in the other linked question it is an s-orbital under consideration, and so the second definition $4\pi r^2|\psi|^2$ is used (or $4\pi r^2 \rho$, where the probability density $\rho$ is simply $|\psi|^2$). However, for anything that is not spherically symmetric, these are not the same and the more general expression $|R|^2r^2$ should be used.
Based on the above arguments this is my take on the matter:
- $|R^2|r^2$ is the best definition of the RDF;
- $4\pi r^2|\psi|^2$ is equivalent to that for spherically symmetric orbitals and agrees with a-cyclohexane-molecule's answer here, but does not work for other orbitals (p, d, ...);
- I still don't fully understand why Housecroft has an additional factor of $4\pi$ – something which I've noticed many years ago. There is probably some motivation behind it, but I don't know what it is. In Feodoran's answer $Y_{0,0}$ is taken to be equal to $1$, which leads to an extra factor of $4\pi$. However, I think it is more conventional to have $Y_{0,0}$ normalised as $\sqrt{1/(4\pi)}$.
† Thanks to Ruslan for pointing this out! If you are curious as to what this looks like, there is a good diagram on LibreTexts, which I have re-uploaded to SE here. Note, however, that the convention in this diagram is $(\rho,\phi,\theta)$ instead of $(r,\theta,\phi)$.
