k0 = Sqrt[s] - m;
k1 = 0;
k2 = 0;
k3 = (Sqrt[s] - m);

p0 = m;
p1 = 0;
p2 = 0;
p3 = 0;

s = Simplify[2*(p*k) + m^2]

で、s=m^2になってしまいます。

すると、k0 = wとか、新しい変数を使う必要があるのでしょうか?

すると、実験室系では、重心系のように

Simplify[y3 //. {p1 -> 0, p2 -> 0,
k0 -> p3, k1 -> 0, k2 -> 0, k3 -> -p3,
q0 -> p0, q1 -> p3*Sqrt[1 - z^2], q2 -> 0, q3 -> p3*z,
j0 -> p3, j1 -> -p3*Sqrt[1 - z^2], j2 -> 0, j3 -> -p3*z,
p0 -> (s + m^2)/(2 Sqrt[s]), p3 -> (s - m^2)/(2 Sqrt[s]),
z -> 1 + t/(2 p3^2), t -> 2 m^2 - s - u}]

とか、纏められるのでしょうか?