This is an old post, and you have probably figured this out, but I will provide an answer to make the post complete.
Consider four-momentum of two objects, object 1 in frame 1, and object 2 in frame 2, with the relative velocity $v_{\rm rel}$ between them. Recall from special relativity in natural untis (c=1) that
$$
p_{1}^\mu = (E_1, \vec{p}_1) = (\gamma_1 m_1, \gamma_1m_1\vec{v}_{1}), \quad \gamma_1 = \frac{1}{\sqrt{1 - v_{1}^2}}
$$
where $v_{1}$ is the relative velocity between frame 1 and some other inertial reference frame (e.g. when boosting from one frame to another).
Consider the inner product between the four-vectors $p_1$ and $p_2$ in the rest-frame of $p_2$, i.e. $p_2=(E_2, \vec{0})$, when boosting between the two frames:
$$
p_1\cdot p_2 = \frac{m_1m_2}{\sqrt{1-v_{\rm rel}^2}}.
$$
Here, $v_{\rm rel}$ is the relative velocity between frame 1 and frame 2. Solving for $v_{\rm rel}$ gives
$$
v_{\rm rel} = \frac{\sqrt{(p_1\cdot p_2)^2 - m_1^2m_2^2}}{p_1\cdot p_2}
$$
which is manifestly Lorentz-invariant, as the product $p_1\cdot p_2$ and masses $m_1,m_2$ are Lorentz-invariant.
Writing the poduct $p_1\cdot p_2$ in an arbitrary frame, we have
$$
p_1\cdot p_2 = E_1E_2 - \vec{p}_1\cdot \vec{p}_2 = \gamma_1\gamma_2m_1m_2 - \gamma_1\gamma_2m_1m_2\vec{v}_1\cdot \vec{v}_2 = \gamma_1\gamma_2m_1m_2(1 - \vec{v}_1\cdot \vec{v}_2)
$$
where we substituted $E_1=\gamma_1m_1$, $\vec{p}_1=\gamma_1m_1\vec{v}_1$ and similary for $E_2$, $\vec{p}_2$. We would like to substitute this into the expression for $v_{\rm rel}$ above, but first we recall some vector-identities:
$$
\begin{align}
(\vec{a}\times \vec{b})^2 &= (\vec{a}\times\vec{b})_i(\vec{a}\times\vec{b})_i \\
&= \epsilon_{ijk}\epsilon_{ilm}a_jb_ka_lb_m \\
&= (\delta_{jl}\delta_{km} - \delta_{jm}\delta_{kl})a_jb_ka_lb_m \\
&= a_ja_jb_kb_k - a_jb_ja_kb_k \\
&= a^2b^2 - (\vec{a}\cdot\vec{b})^2
\end{align}
$$
so we can write
$$
\begin{align}
(p_1\cdot p_2)^2 &= (m_1m_2\gamma_1\gamma_2)^2(1 - 2(\vec{v}_1\cdot \vec{v}_2) + (\vec{v}_1\cdot\vec{v}_2)^2) \\
&= (m_1m_2\gamma_1\gamma_2)^2(1 - 2(\vec{v}_1\cdot \vec{v}_2) + v_1^2v_2^2 - (\vec{v}_1\times\vec{v}_2)^2) \\
&= (m_1m_2\gamma_1\gamma_2)^2((\vec{v}_1 - \vec{v}_2)^2 + 1 - v_1^2 - v_2^2 + v_1^2v_2^2 - (\vec{v}_1\times\vec{v}_2)^2) \\
&=(m_1m_2\gamma_1\gamma_2)^2((\vec{v}_1 - \vec{v}_2)^2 + (1-v_1^2)(1-v_2^2) - (\vec{v}_1\times\vec{v}_2)^2) \\
&= (m_1m_2\gamma_1\gamma_2)^2\left((\vec{v}_1 - \vec{v}_2)^2 + \frac{1}{\gamma_1^2\gamma_2^2}- (\vec{v}_1\times\vec{v}_2)^2\right)
\end{align}
$$
Then we get
$$
(p_1\cdot p_2)^2 - m_1^2m_2^2 = (m_1m_2\gamma_1\gamma_2)^2\left((\vec{v}_1 - \vec{v}_2)^2 - (\vec{v}_1\times\vec{v}_2)^2\right)
$$
Inserting into $v_{\rm rel}$ gives
$$
v_{\rm rel} = \frac{\sqrt{(p_1\cdot p_2)^2 - m_1^2m_2^2}}{p_1\cdot p_2} = \frac{\sqrt{(\vec{v}_1-\vec{v}_2)^2 - (\vec{v}_1\times\vec{v}_2)^2}}{1-\vec{v}_1\cdot\vec{v}_2}
$$
Now we derived two common expressions for the relative velocity in natural units. For colinear motion, i.e. $\vec{v}_1\| \vec{v}_2$, $v_{\rm rel}$ trivially reduces to
$$
v_{\rm rel} = \frac{|v_1-v_2|}{1-v_1v_2} = \frac{E_1E_2|v_1-v_2|}{p_1\cdot p_2} = \frac{|p_1E_2-p_2E_1|}{p_1\cdot p_2}
$$
For the scattering cross-section, we recognize the numerator of $v_{\rm rel}$ from the flux in the cross-section, but we still have the denominator. It is not the relative velocity that is used for the cross-section, it is the Møller velocity, which is defined as the numerator of the relative velocity.
Source:
Landau, L.D.; Lifshitz, E.M. (2002) [1939]. The Classical Theory of Fields. Course of Theoretical Physics. Vol. 2 (4th ed.).