You have two incoming slow momenta (\(P_{1}\) and \(P_{2}\)) and two outgoing slow momenta (\(P_{3}\) and \(P_{4}\)). In order to be consistent with the non-relativistic result, I will assume elasticity:
$$ \Vert P_{1} \Vert^{2} = \Vert P_{3} \Vert^{2}, \qquad \Vert P_{2} \Vert^{2} = \Vert P_{4} \Vert^{2}. $$
Energy-Momentum Vectors
The four external slow momenta can be written as
$$ P_{1} = \begin{pmatrix} m_{1} \cosh{(\xi_{1})} & m_{1} \sinh{(\xi_{1})} \end{pmatrix}; $$
$$ P_{2} = \begin{pmatrix} m_{2} \cosh{(\xi_{2})} & m_{2} \sinh{(\xi_{2})} \end{pmatrix}; $$
$$ P_{3} = \begin{pmatrix} m_{1} \cosh{(\xi_{3})} & m_{1} \sinh{(\xi_{3})} \end{pmatrix}; $$
$$ P_{4} = \begin{pmatrix} m_{2} \cosh{(\xi_{4})} & m_{2} \sinh{(\xi_{4})} \end{pmatrix}. $$
Conservation
The conservation of energy leads to
$$ m_{1} \cosh{(\xi_{1})} + m_{2} \cosh{(\xi_{2})} = m_{1} \cosh{(\xi_{3})} + m_{2} \cosh{(\xi_{4})}. $$
Similarly, the conservation of momentum leads to
$$ m_{1} \sinh{(\xi_{1})} + m_{2} \sinh{(\xi_{2})} = m_{1} \sinh{(\xi_{3})} + m_{2} \sinh{(\xi_{4})}. $$
These two expressions can be combined to yield two equivalent expressions:
$$ m_{1} \exp{(\xi_{1})} + m_{2} \exp{(\xi_{2})} = m_{1} \exp{(\xi_{3})} + m_{2} \exp{(\xi_{4})}; $$
$$ \frac{m_{1}}{\exp{(\xi_{1})}} + \frac{m_{2}}{\exp{(\xi_{2})}} = \frac{m_{1}}{\exp{(\xi_{3})}} + \frac{m_{2}}{\exp{(\xi_{4})}}. $$
2-Mandelstam Invariants
You can introduce three 2-Mandelstam invariants:
$$ s = -(P_{1} + P_{2})^{2} = -(P_{3} + P_{4})^{2}; $$
$$ t = -(P_{1} - P_{3})^{2} = -(P_{2} - P_{4})^{2}; $$
$$ u = -(P_{1} - P_{4})^{2} = -(P_{2} - P_{3})^{2}. $$
Due to the conservation laws, you have
$$ s + t + u = 2 m_{1}^{2} + 2 m_{2}^{2}. $$
Regge-Mandelstam Invariants
You can introduce six Regge-Mandelstam invariants:
$$ r_{12} \equiv -\frac{P_{1} \cdot P_{2}}{m_{1} m_{2}} = \frac{s - m_{1}^{2} - m_{2}^{2}}{2 m_{1} m_{2}} = \cosh{(\xi_{1} - \xi_{2})}; $$
$$ r_{34} \equiv -\frac{P_{3} \cdot P_{4}}{m_{1} m_{2}} = \frac{s - m_{1}^{2} - m_{2}^{2}}{2 m_{1} m_{2}} = \cosh{(\xi_{3} - \xi_{4})}; $$
$$ r_{13} \equiv -\frac{P_{1} \cdot P_{3}}{m_{1}^{2}} = \frac{2m_{1}^{2} - t}{2 m_{1}^{2}} = \cosh{(\xi_{1} - \xi_{3})}; $$
$$ r_{24} \equiv -\frac{P_{2} \cdot P_{4}}{m_{2}^{2}} = \frac{2m_{2}^{2} - t}{2 m_{2}^{2}} = \cosh{(\xi_{2} - \xi_{4})}; $$
$$ r_{14} \equiv -\frac{P_{1} \cdot P_{4}}{m_{1} m_{2}} = \frac{m_{1}^{2} + m_{2}^{2} - u}{2 m_{1} m_{2}} = \cosh{(\xi_{1} - \xi_{4})}; $$
$$ r_{23} \equiv -\frac{P_{2} \cdot P_{3}}{m_{1} m_{2}} = \frac{m_{1}^{2} + m_{2}^{2} - u}{2 m_{1} m_{2}} = \cosh{(\xi_{2} - \xi_{3})}. $$
Since the hyperbolic cosine satisfies
$$ \cosh{(x)} \geq 1; $$
you have the conditions:
$$ s \geq (m_{1} + m_{2})^{2}, \qquad t \leq 0, \qquad u \leq (m_{1} - m_{2})^{2}. $$
These conditions define the range of values allowed for the 2-Mandelstam invariants. The first one says that \(s\) only takes positive values above the threshold value \((m_{1} + m_{2})^{2}\). Thus, the vector \(P_{s} = P_{1} + P_{2} = P_{3} + P_{4}\) is slow and you can write
$$ P_{s} = \begin{pmatrix} \sqrt{s} \cosh{(\xi_{s})} & \sqrt{s} \sinh{(\xi_{s})} \end{pmatrix}. $$
The second condition says that \(t\) only takes negative values. This means that the vector \(P_{t} = P_{1} - P_{3} = P_{4} - P_{2}\) is fast and you can write:
$$ P_{t} = \begin{pmatrix} \sqrt{-t} \sinh{(\rho_{t})} & \sqrt{-t} \cosh{(\rho_{t})} \end{pmatrix}. $$
Note that
$$ P_{s} \cdot P_{t} = \sqrt{s} \sqrt{-t} \sinh{(\xi_{s} - \rho_{t})}. $$
When \(\xi_{s} = 0\) and \(\rho_{t} = 0\), you have the center-of-momentum frame.
Since \(r_{12} = r_{34}\), you have
$$ \cosh{(\xi_{1} - \xi_{2})} = \cosh{(\xi_{3} - \xi_{4})}. $$
This has two solutions:
$$ \xi_{1} - \xi_{2} = \xi_{3} - \xi_{4}; $$
or
$$ \xi_{1} - \xi_{2} = \xi_{4} - \xi_{3}. $$
The first solution implies that
$$ \xi_{3} = \xi_{1} - \xi_{2} + \xi_{4}. $$
Combining this with the conservation law gives
$$ \xi_{3} = \xi_{1}, \qquad \xi_{4} = \xi_{2}. $$
This is the forward solution. The second solution is compatible with \(r_{14} = r_{23}\). Indeed, since you can write the slow rapidity in terms of the slow speed as
$$ \xi = \operatorname{atanh}{(v)}, $$
it follows that
$$ \operatorname{atanh}{(v_{1})} - \operatorname{atanh}{(v_{2})} = \operatorname{atanh}{(v_{4})} - \operatorname{atanh}{(v_{3})}. $$
This is the relativistic analog of the relation
$$ u_{1} - u_{2} = v_{2} - v_{1}; $$
found for the non-relativistic problem.
Although \(r_{13} \neq r_{24}\), you have
$$ t = 2 m_{1}^{2} (1 - r_{13}) = 2 m_{2}^{2} (1 - r_{24}). $$
This means that
$$ m_{1}^{2} (1 - \cosh{(\xi_{1} - \xi_{3})}) = m_{2}^{2} (1 - \cosh{(\xi_{2} - \xi_{4})}). $$
But recall that
$$ \sinh^{2}{\left( \frac{x}{2} \right)} = \frac{\cosh{(x)} - 1}{2}. $$
Thus,
$$ m_{1}^{2} \sinh^{2}{\left( \frac{\xi_{1} - \xi_{3}}{2} \right)} = m_{2}^{2} \sinh^{2}{\left( \frac{\xi_{2} - \xi_{4}}{2} \right)}. $$
Solving for the Outgoing Data
You can first write
$$ \xi_{4} = \xi_{1} - \xi_{2} + \xi_{3}. $$
Then it follows that
$$ \exp{(\xi_{4})} = \frac{\exp{(\xi_{1})} \exp{(\xi_{3})}}{\exp{(\xi_{2})}}. $$
From the conservation law you find
$$ m_{1} \exp{(\xi_{1})} + m_{2} \exp{(\xi_{2})} = m_{1} \exp{(\xi_{3})} + m_{2} \left( \frac{\exp{(\xi_{1})} \exp{(\xi_{3})}}{\exp{(\xi_{2})}} \right). $$
This gives
$$ \exp{(\xi_{3})} = \left[\frac{ m_{1} \exp{(\xi_{1})} + m_{2} \exp{(\xi_{2})} }{m_{1} \exp{(\xi_{2})} + m_{2} \exp{(\xi_{1})} }\right] \exp{(\xi_{2})}. $$
Plugin this back into the first expression gives
$$ \exp{(\xi_{4})} = \left[\frac{ m_{1} \exp{(\xi_{1})} + m_{2} \exp{(\xi_{2})} }{m_{1} \exp{(\xi_{2})} + m_{2} \exp{(\xi_{1})} }\right] \exp{(\xi_{1})}. $$
The speeds can be introduced back via the relation
$$ v = \tanh{(\xi)} \quad \Longrightarrow \quad \exp{(2 \xi)} = \frac{1 + v}{1 - v}. $$
You find
$$ \frac{\sqrt{1 + v_{3}}}{\sqrt{1 - v_{3}}} = \frac{\sqrt{1 + v_{2}}}{\sqrt{1 - v_{2}}} \frac{m_{1} \sqrt{(1 + v_{1})(1 - v_{2})} + m_{2} \sqrt{(1 - v_{1})(1 + v_{2})}}{m_{1} \sqrt{(1 - v_{1})(1 + v_{2})} + m_{2} \sqrt{(1 + v_{1})(1 - v_{2})}}; $$
$$ \frac{\sqrt{1 + v_{4}}}{\sqrt{1 - v_{4}}} = \frac{\sqrt{1 + v_{1}}}{\sqrt{1 - v_{1}}} \frac{m_{1} \sqrt{(1 + v_{1})(1 - v_{2})} + m_{2} \sqrt{(1 - v_{1})(1 + v_{2})}}{m_{1} \sqrt{(1 - v_{1})(1 + v_{2})} + m_{2} \sqrt{(1 + v_{1})(1 - v_{2})}}. $$
This leads to the identity
$$ \frac{(1 + v_{3})(1 - v_{4})}{(1 - v_{3})(1 + v_{4})} = \frac{(1 + v_{2})(1 - v_{1})}{(1 - v_{2})(1 + v_{1})}. $$
When \(m_{2} = m_{1}\), this solution reduces to
$$ \xi_{3} = \xi_{2}, \qquad \xi_{4} = \xi_{1}. $$
This is the exchange solution.
The article in Wikipedia on elastic collisions was helpful and it was fun to verify its results.
Composite Rapidities
The slow rapidity \(\xi_{s}\) associated to \(P_{s}\) can be writen as
$$ \xi_{s} = \log{\left( \frac{m_{1} \exp{(\xi_{1})} + m_{2} \exp{(\xi_{2})}}{\sqrt{s}} \right)} = \log{\left( \frac{m_{1} \exp{(\xi_{3})} + m_{2} \exp{(\xi_{4})}}{\sqrt{s}} \right)}. $$
Similarly, the fast rapidity \(\rho_{t}\) associated to \(P_{t}\) can be written as
$$ \rho_{t} = \log{\left( \frac{m_{1} \exp{(\xi_{1})} - m_{1} \exp{(\xi_{3})}}{\sqrt{-t}} \right)} = \log{\left( \frac{m_{2} \exp{(\xi_{4})} - m_{2} \exp{(\xi_{2})}}{\sqrt{-t}} \right)}. $$
Regimes
The condition on \(u\) does not allow for a specific kind of rapidity for the vector \(P_{u} = P_{1} - P_{4} = P_{3} - P_{2}\). Depending on the values of \(s\) and \(t\), you can have all three cases.
Slow \(u\)-channel
In the region
$$ 0 < u \leq (m_{1} - m_{2})^{2} $$
the vector \(P_{u}\) is slow and given by
$$ P_{u} = \begin{pmatrix} \sqrt{u} \cosh{(\xi_{u})} & \sqrt{u} \sinh{(\xi_{u})} \end{pmatrix}. $$
Then
$$ P_{s} \cdot P_{u} = - \sqrt{s} \sqrt{u} \cosh{(\xi_{s} - \xi_{u})}; $$
and
$$ P_{t} \cdot P_{u} = - \sqrt{-t} \sqrt{u} \sinh{(\rho_{t} - \xi_{u})}. $$
The composite slow rapidity \(\xi_{u}\) is given by
$$ \xi_{u} = \log{\left( \frac{m_{1} \exp{(\xi_{1})} - m_{2} \exp{(\xi_{4})}}{\sqrt{u}} \right)} = \log{\left( \frac{m_{1} \exp{(\xi_{3})} - m_{2} \exp{(\xi_{2})}}{\sqrt{u}} \right)}. $$
Null \(u\)-channel
If \(u = 0\), then the vector \(P_{u}\) is null and can be written as
$$ P_{u} = \begin{pmatrix} k_{u} & k_{u} \end{pmatrix}. $$
Then
$$ P_{s} \cdot P_{u} = \sqrt{s} k_{u} \exp{(-\xi_{s})}; $$
and
$$ P_{t} \cdot P_{u} = \sqrt{-t} k_{u} \exp{(\rho_{t})}. $$
You also have
$$ k_{u} = m_{1} \cosh{(\xi_{1})} - m_{2} \cosh{(\xi_{4})} = m_{1} \sinh{(\xi_{1})} - m_{2} \sinh{(\xi_{4})}. $$
From here it follows that
$$ \xi_{4} = \xi_{1} - \log{\left( \frac{m_{1}}{m_{2}} \right)}. $$
Using the conservation law, you can also write
$$ k_{u} = m_{1} \cosh{(\xi_{3})} - m_{2} \cosh{(\xi_{2})} = m_{1} \sinh{(\xi_{3})} - m_{2} \sinh{(\xi_{2})}. $$
From here it follows that
$$ \xi_{3} = \xi_{2} + \log{\left( \frac{m_{1}}{m_{2}} \right)}. $$
Note that
$$ \xi_{1} + \xi_{2} = \xi_{3} + \xi_{4}. $$
Fast \(u\)-channel
If \(u < 0\), then the vector \(P_{u}\) is fast and given by
$$ P_{u} = \begin{pmatrix} \sqrt{-u} \sinh{(\rho_{u})} & \sqrt{-u} \cosh{(\rho_{u})} \end{pmatrix}. $$
Then
$$ P_{s} \cdot P_{u} = \sqrt{s} \sqrt{-u} \sinh{(\xi_{s} - \rho_{u})}; $$
and
$$ P_{t} \cdot P_{u} = \sqrt{-t} \sqrt{-u} \cosh{(\rho_{t} - \rho_{u})}. $$
The composite fast rapidity \(\rho_{u}\) can be written as
$$ \rho_{u} = \log{\left( \frac{m_{1} \exp{(\xi_{1})} - m_{2} \exp{(\xi_{4})}}{\sqrt{-u}} \right)} = \log{\left( \frac{m_{1} \exp{(\xi_{3})} - m_{2} \exp{(\xi_{2})}}{\sqrt{-u}} \right)}. $$
