Cross-section in relativistic limit: Fermi's golden rule still valid? In order to calculate the cross-section of an interaction process the following formula is often used for first approximations:
$$
\sigma = \frac {2\pi} {\hbar\,v_i} \left| M_{fi}\right|^2\varrho\left(E_f\right)\,V
$$
$$
M_{fi} = \langle\psi_f|H_{int}|\psi_i\rangle
$$
Very often plane waves are assumed for the final state and therefore the density of states is given by
$$
\varrho\left(E_f\right) = \frac{\mathrm d n\left(E_f\right)}{\mathrm d E_f} = \frac{4\pi {p_f}^2}{\left(2\pi\hbar\right)^3}\frac V {v_f}
$$
I understand the derivation of this equation in the context of the non relativistic Schrödinger equation. But why can I continue to use this formula in the relativistic limit: $v_i, v_f \to c\,,\quad p_f\approx E_f/c$. 
Very often books simply use this equation with matrix element derived from some relativistic theory, e.g. coupling factors and propagators from the Dirac equation or Electroweak interaction. How is this justified?
Specific concerns:


*

*Is Fermi's golden rule still valid in the relativistic limit?

*Doesn't the density of final states has to be adapted in the relativistic limit?
 A: Fermi's golden rule still applies in the relativistic limit, and can be rewritten in a Lorentz invariant fashion.  Starting with the transition probability
$$ W_{i\rightarrow f} = \frac{2\pi}{\hbar} |m_{if}|^2 \rho(E) \,,$$ to have $W$ Lorentz invariant we'd like both the matrix element $|m_{if}|^2$ and the density of final states $\rho(E)$ to be invariant.
This can be done by shifting a few terms around. A little bit of handwaving to motivate it: The wave function $\psi$ (which is in the matrix element) has to be normalized by $\int |\psi|^2 dV = 1$, which gives us a density (of probability to encounter a particle) of $1/V$.  Now, a boosted observer experiences length contraction of $1/\gamma$, which changes the density to $\gamma/V$.  To obtain the correct probability again, we should re-normalize the wave function to $\psi' = \sqrt{\gamma}\,\psi $ by pulling the Lorentz factor out.
So we intoduce a new matrix element $$|{\cal M}_{if}|^2 = |m_{if}|^2 \prod_{i=1}^n (2 \gamma_i m_i c^2) =|m_{if}|^2  \prod_{i=1}^n (2E_i)^2 $$ (this is for an $n$-body process).  Now the transition probability (here in differential form) becomes:
$$
  dW = \frac{2\pi}{\hbar} \frac{|{\cal M}_{if}|^2}{ (2E_1)^2  (2E_2)^2  \cdots}
  \cdot  \frac{1}{(2\pi\hbar)^{3n}} 
  \, d^3p_1 \, d^3p_2 \, \cdots 
  \delta({p_1}^\mu + {p_2}^\mu + \ldots - {p}^\mu )
$$
The delta function is there to ensure conservation of momentum and energy.  Now we can regroup the terms:
$$
  \Rightarrow \quad
  dW = \frac{2\pi}{\hbar} \frac{|{\cal M}_{if}|^2}{ 2E_1 2E_2 \cdot \ldots}
  \cdot d_\mathrm{LIPS}
$$
The density of states/"phase space" $d\rho$ is replaced by a relativistic version, sometimes called the Lorentz invariant phase space $d_\mathrm{LIPS}$, which is given by
$$
  d_\mathrm{LIPS} =  \frac{1}{(2\pi\hbar)^{3n}} 
  \prod_{i=1}^n \frac{d^3p_i}{ 2E_i }
  \delta\left(\prod_{i=1}^n {p_i}^\mu - {p}^\mu \right) \,.
$$
The nice thing about the relativistic formula for $dW$ is that, in the case you are scattering particles off one another, it immediately shows us three important contributions: not only the matrix element and phase space, but also the flux factor $1/s$ (where $s = ({p_1}^\mu + {p_2}^\mu)^2$ is the Mandelstam variable, and in case the masses are negligible, $ s \approx 2 E $).  This flux factor is responsible for the general $1/Q^2$ falling slope when you plot cross section over momentum transfer $Q = \sqrt{s}$, which comes entirely from relativistic kinematics.
Hope this answers your questions. Here is a presentation (PDF) that sums it up, with an explicit proof that it is Lorentz invariant.
