$(p_1^2 + p_2^2) = (E_1 + E_2)^2 - (\vec p_1 + \vec p_2)^2$
This equation should hold in any inertial frame, so let's boost into the COM frame. In the center of mass frame (by definition), $\vec p_1 + \vec p_2 =0$ (since each are the regular 3-momenta of the particles), so
$=(E_1 + E_2)^2 = (2E)^2$ since the energy of each particle is the same in the center of mass frame we just express each energy as $E$.
Note that all the calculations so far have been in the COM frame. So we need the energies of the particles, in this frame. We know their masses, so we need the momenta of the particles, in this frame.
Momentum of one of the particles in the COM frame $P^1_{cm} \equiv p_1 - \dfrac{\vec{p}_1 + \vec{p}_2}{2}$ rather than momentum of center of mass. The other particle will have the opposite 3-momentum, which means that they'll both have the same energy in this frame:
$$E_{cm}^2 = P_{cm}^2 + m^2$$
But then the step I don't understand is the following:
$E^2 = P_{cm}^2 + m^2$
which gets plugged in to $(2E)^2$ to get the final answer.
Quantum field theory does not offer a description of "how" its processes work, just like Newtonian mechanics doesn't offer an explanation of "how" forces impart acceleration or general relativity an explanation of "how" the spacetime metric obeys the Einstein equations.
The predictions of quantum field theory, and quantum electrodynamics (QED) in particular, are well-tested. Given two photons of sufficient energy to yield at least the rest mass of an electron-positron pair, one finds that QED predicts a non-zero amplitude for the process $\gamma\gamma \to e^+ e^-$ to happen. That is all the theory tells us. No "fluctuation", no "virtual particles", nothing. Just a cold, hard, quantitative prediction of how likely such an event is.
All other things - for instance the laughable description in the Wikipedia article you quote - are stories, in this case a human-readable interpretation of the Feynman diagrams used to compute the probability of the event, but should not be taken as the actual statement the quantitative theory makes.
There is no "how", what happens between the input and the output of a quantum field theoretic process is a black box called "time evolution" that has no direct, human-readable interpretation. If we resolve it perturbatively with Feynman diagrams, people like to tell stories of virtual particles, but no one forces us to do that - one may organize the series in another way, may be even forced to do so (e.g. at strong coupling), or one may not use a series at all to compute the probability. The only non-approximative answer to "how" the scattering processes happen in quantum field theory that QFT has to offer is to sit down and derive the LSZ formula for scattering amplitudes from scratch, as it is done in most QFT books. Which, as you may already see from the Wikipedia article, is not what passes as a good story in most circles.
But neither nature nor our models of it are required to yield good stories. Our models are required to yield accurate predictions, and that is what quantum field theory does.
Best Answer
As luksen and Marek point out, when the incoming particle is massless this reaction can't go, unless there's a third particle around to soak up the momentum. If we allow for the existence of this particle, then I don't think there's enough information to solve.
Let $p_{Xi},p_{Xf}$ be the initial and final 4-momenta of that particle. Then conservation of 4-momentum says $$ p_\gamma+p_{Xi}=p_1+p_2+p_{Xf} $$ Assume that the extra particle was initially at rest: $p_{Xi}=(M,0,0,0)$. ($M$ is the mass of $X$, and I'm setting $c=1$.) We don't know anything about any of the components of $p_\gamma$ or $p_{Xf}$, so there are eight unknowns to solve for. The above is 4 equations, so we need 4 additional equations. Unfortunately, we have only two, $$ E_\gamma=|\vec p_\gamma|, $$ $$ E_{Xf}^2=|\vec p_{Xf}|^2+M^2 $$ (Here $E$ means energy and $\vec p$ means 3-momentum in the frame we're working in, of course. These two equations are just $p_\mu p^\mu=m^2$, if you prefer.) So there's not enough information to solve.
The typical assumption is that the extra particle is much more massive than the others. In this case, the energy it absorbs is typically negligible, but the 3-momentum isn't. So if the 3-momentum is what you're interested in, I don't think you can get by with just pretending that particle isn't there.