Surely someone has found the solutions to the hard sphere collisions (in $n$ dimensions) of two bodies of mass $m_1$ and $m_2$, respectively--that is the resultant velocities (or momenta) of the two bodies after a hard sphere collision given $\vec{p}_1$ and $\vec{p}_2$. I have a solution (both for elastic and inelastic collisions), but I wonder whether or not it's correct.
FYI, here are my solutions for elastic, hard-sphere collisions:
$$ \vec{p}_{1f} = \vec{p}_1 + 2\frac{\left(m_1\vec{p}_2 - m_2\vec{p}_1\right)\circ \hat{r}}{m_1 + m_2}\hat{r}, \\ \vec{p}_{2f} = \vec{p}_2 - 2\frac{\left(m_1\vec{p}_2 - m_2\vec{p}_1\right)\circ \hat{r}}{m_1 + m_2}\hat{r}. $$
wher $\hat{r}$ is in the direction of the "hard sphere" (i.e. the unit vector connecting the two centers of the two spheres). You can see that my formulation trivially conserved momentum--the question is whether or not it conserves energy as well.
If it helps, here is my full explanation (this is quite long and the link will expire in 30 days from 09/23/2014). For these purposes, you should skip to the section on n-dimensional cases (page 6). I wrote this as a reference to myself for doing physics simulations since I could not find any reference on hard sphere collisions (other than references setting up the problem which give no solutions).
Here is an abbreviated work on finding the final momenta from the elastic case:
Just as in the $1$D case, we add and subtract the same momentum from each particle--except now we are adding and subtracting a vector: $\vec{p} = p\hat{r}$, which means that $|\vec{p}| = p$ (such that $\vec{p}_{1f} = \vec{p}_1 + p\hat{r}$ and $\vec{p}_{2f} = \vec{p}_2 - p\hat{r}$). This gives the final energy as:
\begin{align} E_f = \frac{|\vec{p}_1 + p\hat{r}|^2}{2m_1} + \frac{|\vec{p_2} - p\hat{r}|^2} {2m_2} \end{align}
We now write the final energy in terms of the magnitude of $p$ and the dot product between $\hat{r}$ and $\vec{p}_1$ and likewise for $\vec{p}_2$. Note that if you choose the direction for $\vec{p}$ to act a priori (like in the hard-sphere case), the dot products $\vec{p}_1\cdot \hat{r}$ and $\vec{p}_2\cdot\hat{r}$ will be known:
\begin{align*} E_f = \frac{p_1^2 + p^2 + 2p\hat{r}\circ\vec{p}_1}{2m_1} + \frac{p_2^2 + p^2 - 2p\hat{r}\circ\vec{p}_2}{2m_2} \end{align*}
As usual the elastic case greatly simplifies:
\begin{align*} \require{cancel} \cancel{\frac{p_1^2}{2m_1}} + \cancel{\frac{p_2^2}{2m_1}}= \frac{\cancel{p_1^2} + p^2 + 2p\hat{r}\circ\vec{p}_1}{2m_1} + \frac{\cancel{p_2^2} + p^2 - 2p\hat{r}\circ\vec{p}_2}{2m_2} \\ \frac{m_1 + m_2}{2m_1m_2}p^2 - 2\frac{m_1 \vec{p}_2\circ\hat{r} - m_2\vec{p}_1\circ\hat{r}}{2m_1m_2}p = 0 \end{align*}
This leads to a trivial solution of $p = 0$ (no collision occurs) and:
\begin{align} \frac{m_1 + m_2}{2m_1m_2}p - 2\frac{m_1 \vec{p}_2\circ\hat{r} - m_2\vec{p}_1\circ\hat{r}}{2m_1m_2} = 0 \\ p = 2\frac{m_1 \vec{p}_2\circ\hat{r} - m_2\vec{p}_1\circ\hat{r}}{m_1 + m_2} \end{align}
Hence my solutions of:
$$ \vec{p}_{1f} = \vec{p}_1 + 2\frac{\left(m_1\vec{p}_2 - m_2\vec{p}_1\right)\circ \hat{r}}{m_1 + m_2}\hat{r}, \\ \vec{p}_{2f} = \vec{p}_2 - 2\frac{\left(m_1\vec{p}_2 - m_2\vec{p}_1\right)\circ \hat{r}}{m_1 + m_2}\hat{r}. $$
Finding the missing equations
Coordinate transforms just complicate the issue. The heart of the matter is that in n dimensions you have n degrees of freedom for the velocities of the COM of each sphere, and you only have n momentum conservation equations plus one energy conservation equation. That means you need an additional n-1 equations to solve the problem.
The additional equations come from the fact that if you assume a frictionless impact*, then the impulse/force imparted to each object must be normal to the surface. This gives a constraint on the direction of the impulse. Along with the magnitude, this direction would give us n degrees of freedom. So without it, it gives us n-1: exactly the number we need.
Solving for final velocities
Since we have spheres, the direction normal to the impact surface will be the vector that points from one COM to the other. Let's call this normal vector $\hat{j}$.
Now the velocities of balls A and B can be given as $V_{a1}$ and $V_{b1}$ before the collision and $V_{a2}$ and $V_{b2}$ after the collision.
Let's call the impulse that is imparted on ball A $J$. Then we have:
$$J \bullet\hat{j} =|J|$$ $$ V_{a1} + \frac{J}{M_a}=V_{a2}$$ $$ V_{b1} - \frac{J}{M_b}=V_{b2}$$
Then we can write the conservation of energy (multiplied by 2) as: $$M_a V_{a1} \bullet V_{a1} + M_b V_{b1} \bullet V_{b1} = M_a V_{a2} \bullet V_{a2} + M_b V_{b2} \bullet V_{b2}$$
Plugging in for our velocities after to collision gives:
$$M_a V_{a1} \bullet V_{a1} + M_b V_{b1} \bullet V_{b1} = M_a \left( V_{a1} + \frac{J}{M_a}\right) \bullet \left( V_{a1} + \frac{J}{M_a}\right) + M_b \left( V_{b1} - \frac{J}{M_b}\right) \bullet \left( V_{b1} - \frac{J}{M_b}\right)$$
$$M_a V_{a1} \bullet V_{a1} + M_b V_{b1} \bullet V_{b1} = M_a V_{a1}\bullet V_{a1} + 2 J\bullet V_{a1} + \frac{J\bullet J}{M_a} + M_b V_{b1}\bullet V_{b1} - 2 J\bullet V_{b1} + \frac{J\bullet J}{M_b}$$
$$0 = 2 J\bullet V_{a1} + \frac{J\bullet J}{M_a} - 2 J\bullet V_{b1} + \frac{J\bullet J}{M_b}$$
dividing by the magnitude of $J$ yields:
$$0 = 2 \hat{j} \bullet (V_{a1}-V_{b1}) + |J|\left(\frac1{M_a} + \frac1{M_b}\right)$$
$$|J| = 2 \frac{\hat{j} \bullet (V_{a1}-V_{b1})}{\frac1{M_a} + \frac1{M_b}}$$
Final Velocities
$$V_{a2}=V_{a1}+2 \frac{\hat{j} \bullet (V_{a1}-V_{b1})}{M_a\left(\frac1{M_a} + \frac1{M_b}\right)}$$ $$V_{b2}=V_{b1}+2 \frac{\hat{j} \bullet (V_{b1}-V_{a1})}{M_b\left(\frac1{M_a} + \frac1{M_b}\right)}$$
Finding Direction of Impulse
Let's call the vector pointing from the center of ball B to the center of ball A $R$. We know that before the impact $R=R_0+R't$ where $R_0$ is the difference in the balls' initial positions $P_a-P_b$ and $R'$ is the difference in the balls' initial velocities $V_{a1}-V_{b1}$.
At the moment of impact: $$\hat{j}=\frac{R}{|R|}$$ and $$|R|=r_a+r_b$$ Where $r_a$ and $r_b$ are the radii of the balls.
We can transform this second equation to find t:
$$R\bullet R=(r_a+r_b)^2$$ $$R'\bullet R't^2+2R_0\bullet R't+R_0\bullet R_0-(r_a+r_b)^2=0$$ $$t=-\frac{R_0\bullet R'+\sqrt{(R_0\bullet R')^2-R'\bullet R'(R_0\bullet R_0-(r_a+r_b)^2)}}{R'\bullet R'}$$
If this time is negative or imaginary, then the collision will not occur. Otherwise $R$ at the moment of impact can be found by plugging in t:
$$R=R_0-R'\frac{R_0\bullet R'+\sqrt{(R_0\bullet R')^2-R'\bullet R'(R_0\bullet R_0-(r_a+r_b)^2)}}{R'\bullet R'}$$
Wrap Up
So given initial positions $P_a$ and $P_b$, initial velocities $V_{a1}$ and $V_{b1}$, and ball radii $r_a$ and $r_b$ the final velocities $V_{a2}$ and $V_{b2}$ can be found be following the following procedure:
Calculate the relative positions and velocities: $$R_0=P_a-P_b$$ $$R'=V_{a1}-V_{b1}$$
Calculate impact time: $$t=-\frac{R_0\bullet R'+\sqrt{(R_0\bullet R')^2-R'\bullet R'(R_0\bullet R_0-(r_a+r_b)^2)}}{R'\bullet R'}$$ If t is negative or imaginary, then no collision occurs, and the initial velocities are the final velocities. Otherwise:
Calculate the impulse direction: $$\hat{j}=\frac{R_0+R't}{|R_0+R't|}$$ Calculate the impulse: $$J=2 \hat{j}\frac{\hat{j} \bullet R'}{\frac1{M_a} + \frac1{M_b}}$$ Calculate the final velocities: $$V_{a2} = V_{a1} + \frac{J}{M_a}$$ $$V_{b2} = V_{b1} - \frac{J}{M_b}$$
*The collision must be frictionless if we're not considering moments of inertia as any friction would apply torque to the ball which would transfer energy and rotational inertia to the ball. The problem becomes slightly more complicated in that case as one has to assume how the friction will cause the balls to interact, which is not trivial as many of the simple assumptions conflict with energy conservation.
