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)?
"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.
$\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.
