Probability density in Fermi golden rule Consider Fermi golden rule
$$\Gamma _{{i\rightarrow f}}={\frac  {2\pi }{\hbar }}\left|\langle f|H'|i\rangle \right|^{{2}}\rho $$
I don't understand why $\left|\langle f|H'|i\rangle \right|^{{2}}$ is defined as a probability density. 
Since it is an integral over the volume of normalization, in my view, it should be a probability itself and not a density. 
What kind of probability density is it? In other words, w.r.t. which variable should I integrate it to get the total probability (i.e. 1)?
 A: "Probability density" is nowhere to be found in the WP article you are link-quoting. I'll try to detail @rob 's recommendation that should have sufficed to answer your question, together with the article linked. Set $\hbar=1$ here for simplicity, easily reinstatable by elementary dimensional analysis. I sense the confusion might reside in the bra-ket notation used. 
Recall states such as $|i\rangle$, etc are dimensionless... they are not dimensionful like $|x\rangle$ or $\psi(x)=\langle x|\psi\rangle$. So thinking about spacial probabilities and volume integrals is a canard: you never considered space integrals at the level of this treatment (but you may need them to produce the final matrix element). You just do time-dependent perturbation theory w.r.t. time-varying dimensionless coefficients $a(t)$.
In these units, then, the transition amplitude $\langle f|H'|i  \rangle $ has dimensions of frequency, so, inverse time. Squared, it produces the decay rate to a single state k, with energy very, very close to that of i:
$$\Gamma_{i\rightarrow k }= \frac{\mathrm{d}}{\mathrm{d}t} \left|a_k(t)\right|^2 \sim    {2|\langle k| H'|i\rangle |^2} ~  t. $$ 
Integrated over a small time interval, it will give you the probability of your single state decay $i \rightarrow k$ during that interval. You know the drill of accounting for a collection of i s...
Putting lots of final states into the decay channels, their uneven contributions to the relevant integral  will give you
$$  
\Gamma_{i \rightarrow f}=  {2 \pi}   \left | \langle f|H'|i  \rangle \right |^{2} \rho ,
$$
where ρ, the density of states, has dimensions of time, since integration over its argument, frequency or energy, must net one (dimensionless). In total,  the decay rate $\Gamma_{i \rightarrow f}$ has 
dimensions of inverse time, or energy (the width of plots versus energy). 
"Golden" refers to the rather surprising outcome that it is constant.
A: $\Gamma_{i \to f}$ refers in some sense to the probability per unit time of an initial state $i$, being found in the particular final state, $f$. The probability of $i$ transitioning to $\textit{some}$ state will be (hopefully) $1$, so $|\langle f|H'|i\rangle|^2$ is a "function" of the final state, such that summing (or integrating) over the final states will give you unity. 
For example, you can have an incoming plane wave solution scattering off of a localized potential which is $0$ outside some radius (like a hard sphere). If a ping pong ball bounces off a bowling ball, it has a probability density of being directed in any direction $(\theta, \phi)$ (these are the final states - the momentum is usually uniquely determined by the direction for simple cases so the state itself can be thought of as just the direction). Integrating over a region of $\theta$ and $\phi$ values (final states) gives the probability (not the density) of the ping pong going in that direction. 
