Why is $\gamma \to 2\gamma$ kinematically forbidden?
It's not strictly forbidden for a massless particle to split into two, but the "volume of available phase space" is zero. Remember that the rates for processes are computed by finding a so-called matrix element, and then integrating it over phase space. The allowed phase space here is of measure zero, because you have the additional constraint that the two final particles be exactly collinear.
More concretely, suppose we applied some regulator like a discrete lattice. Depending on how the regulator is set up, the rate for the process you cited may be nonzero, or exactly zero. But for a reasonable regulator, the rate should limit to zero as it is removed, by the argument above.
I just want to compile here a curated list of all the things I've seen said about this question here on this site, both right and wrong. The only answer I find truly satisfactory is "Correct Answer 3", as given by knzhou above.
Question: Can a massless particle spontaneously split in two?
Correct Answer 1:(if the particle is a photon) No, because this would violate conservation of angular momentum. The photon must have helicity $\pm 1$ (for reasons I don't fully understand). There's no way to have 3 numbers, all of them $\pm 1$, add up to zero. Of course, this argument relies on the photon having spin 1.
Correct Answer 2:(if the particle is a photon) No, by Furry's Theorem in QED. This is a dynamical argument rather than a kinematic one.
Correct Answer 3:(for any massless particle) Perhaps in some sense yes, but only with probability zero. Consider 1 incoming particle with momentum $P_1$ and 2 outgoing particles with momenta $P_2,P_3$ and masses $m_1,m_2,m_3$, in $D+1$ dimensions. The momenta $P_a^\mu$ have $3(D+1)$ degrees of freedom. Conservation of momentum says that $P_1^\mu = P_2^\mu + P_3^\mu$, which is $D+1$ constraints. The mass constraint says that $\eta_{\mu\nu}P^\mu_a P^\nu_a = (m_a)^2$, 3 more constraints. These constraints are "generically" independent, as when $m_2 + m_3 < m_1$, so they result in an "expected dimension" of the phase space of allowed configurations of $3(D+1) - (D+1) - 3 = 2D-1$ (i.e. this is the actual dimension when $m_2 + m_3 < m_1$). But when $m_1 = m_2 = m_3 = 0$, the dimension of phase space drops to $D+1$ because of positivity considerations -- the $P_a$ must be null and collinear [1]. So the dimension has dropped by $2D-1 - (D+1) = D-2$. Thus when $D \geq 3$, the dimension is smaller than "expected". Now, any reasonable integral over this phase space will be integrating a quantity of dimension $2D-1$ (the "expected dimension"), but in the massless case with $D \geq 3$, the integration will only be over a space of dimension $D + 1 < 2D-1$, so it will come out zero. As the probability of such an event happening will always be calculated by such an integral, the probability will always be zero.
Marginally Correct Answer 4:(if the particle is a photon) No; this is not observed experimentally. This is maybe an answer, though not an explanation. Even as an answer, it's marginal, because to be a serious answer, one would need to do some theoretical analysis to figure out where one should look to find such a strange effect and then actually search for it. Maybe it turns out to be as simple as looking at the lamp next to me and observing that its light stays yellow rather than becoming infrared as it cross the room -- but maybe not -- it depends on the theoretical analysis. This begs the question -- what goes into such an analysis?
Possibly Correct Answer 5:(for any massless particle) No; for any "decay" process there must be a nonzero drop in total invariant mass from the incoming particles to the outgoing ones. I'm not quite sure whether this answer is correct -- certainly some physical processes occurring in nature lead to a gain of invariant mass -- witness the existence of elements heavier than iron! So it seems this argument hinges on what exactly a "decay" is, and whether the hypothetical $\gamma \to 2\gamma$ would count as a "decay". At any rate, this argument doesn't seem "kinematical" to me, because kinematics is always time-reversible, and this argument rests on a fundamentally non-time-reversible premise.
Incorrect Answer 1:(for any massless particle) The momenta given in the counterexample above don't represent massless particles because their speed is not 1. In fact, the speed of $(\frac 1 2, -\frac 1 2, 0, 0)$ is $\frac{\frac 1 2}{\frac 1 2} = 1$. Remember that unlike massive particles, there is no canonical normalization for the 4-momentum vector of a massless particle.
Incorrect Answer 2:(for any massless particle) No; the momenta given in the counterexample above don't "really" represent 3 distinct particles. In fact they do -- there's a different between two photons of energy $\frac 1 2$ and one photon of energy $1$. If that doesn't convince you, I hope you sort yourself out eventually, but you can also trivially modify the counterexample to have $P_2 = (-\frac 1 3, \frac 1 3, 0,0)$ and $P_3 = (-\frac 2 3, \frac 2 3, 0, 0)$.
Incorrect Answer 3:(for any massless particle) No, because we could pass to a reference frame where $P_2$ and $P_3$ are going in opposite directions; this would put $P_1$ at rest, which is impossible because it is massless. This argument overlooks the possibility of $P_1,P_2,P_3$ being collinear, which as we've seen is perfectly consistent with momentum conservation.
Incorrect Answer 4:(for any massless particle) No, because the Lorentz group is noncompact. This argument was given as an answer to the linked question of which this one is a duplicate. This argument is not explicit enough to say what's wrong it, but the above counterexample shows it must be wrong.
[1] Let us carefully work out the analysis of Correct Answer 3 above. Choose coordinates where
$P_1 = (1,-1,0,\dots,0)$.
Write $P_2 = (P^0,P^1,P^2,\dots,P^D)$.
Then by momentum conservation ($P_1 = P_2 + P_3$) we have
$P_3 = (1-P^0,-1-P^1,-P^2,\dots,-P^D)$.
The mass constraints say that
$(1) \qquad 0 = \eta_{\mu\nu}P_2^\mu P_2^\nu = (P^0)^2 - (P^1)^2 - (P^2)^2 - \dots -(P^D)^2$ and
$(2) \qquad 0 = \eta_{\mu\nu}P_3^\mu P_3^\nu = (1-P^0)^2 - (-1-P^1)^2 - (P^2)^2 - \dots -(P^D)^2$.
Equating the two expressions, we obtain
$0 = -2P^0 - 2P^1$, i.e. $P^1 = -P^0$.
Therefore, from the first equation we have
$0 = (P^2)^2 + \dots + (P^D)^2$, which forces $P^2 = \cdots = P^D = 0$ because the $P^\mu$ are real.
So we find that $P_2,P_3$ are each scalar multiples of $P_1$. So in coordinate-independent language, the solutions to our equation come when the $P_a$ are collinear, null, and satisfy $P_1 = P_2 + P_3$. So the space of solutions is given by choosing
the null vector $P_1$ (that's $D$ dimensions), plus
a scalar $c$ such that $P_2 = cP_1$ (one more dimension)
(we then take $P_3 = (1-c)P_1$). So the space of solutions has dimension $D+1$. The "dimension deficit" is
$(\text{"expected dimension"}) - (\text{actual dimension}) = (2D-1) - (D+1) = D-2$.
For those keeping track of home, what has happened is that in this case we effectively get the $D-1$ additional constraints $P^\mu = 0$ for $\mu \geq 2$, but one dimension of degeneracy has also crept into the system of constraints, which is why we end up with a "dimension deficit" of $D-2$. In fact, when $D = 1$, the positivity considerations are trivial and it's only the degeneracy that matters -- in this case we end up with a larger phase space than expected!