**Mathematical Derivation** of the fact that all planets will orbit on the same plane There is some discussion about this Here. If an isolated system of particles under gravitational force is allowed to decrease its energy by means such as inelastic collision, then eventually all particles will end up on the same plane. Reasons for that are already explained in many other posts, but NONE of them includes any mathematical formulation and quantitative analysis.
But HERE, I am looking for a mathematical formulation of such a system. Here is what I have tried. DO NOT complain that I have made ridiculous assumptions: I know what's actually happening with the solar system, but to model it with math, I have to simplify it, because I am unable to come up with better models.
Let there be $N$ particles with positions $\mathbf x_k$ and mass $m_k$. When particles $i,j$ collide, if they stick together, they lost some kinetic energy. The amount of energy lost is
$$
\frac{1}{2} (m_i\dot{\mathbf x}_i^2+m_j\dot{\mathbf x}_j^2)-\frac{1}{2}(m_i+m_j)\left(\frac{m_i\dot{\mathbf x}_i+m_j\dot{\mathbf x}_j}{m_i+m_j}\right)^2\\
=\frac{m_im_j(\dot{\mathbf x}_i^2+\dot{\mathbf x}_j^2-2\dot{\mathbf x}_j\cdot \dot{\mathbf x}_i)}{2(m_i+m_j)}=\frac{m_im_j|\dot{\mathbf x}_i-\dot{\mathbf x}_j|^2}{2(m_i+m_j)}.
$$
So I get an idea of how much energy is lost. However, I can't see why this makes the planets closer to a plane.
If I use this collision model, can I prove that eventually, after many collisions, all particles will be very "close" to a plane? I know that the plane should be perpendicular to the angular momentum vector, but how can I express the notion of "close" with formulae.
Is there a better model to help me predict how much time elapses before all particles are sufficiently close to the plane? i.e. $\mathbf x_k\cdot \mathbf J<\epsilon$ for some small $\epsilon$ we choose. $\mathbf J$ is angular momentum.
I am an ESL learner, so please ask me to clarify things if there are any errors.
 A: For the first part of your question I believe that global conservation of angular momentum can mostly explain it:
We are assuming that we have a very large number of particles, already the three body problem is not possible to solve analytically. What we can do is look at individual collisions:
From a quick search in Wikipedia I found the equations for an inelastic collision, which are:
$$\mathbf{v}_i=\frac{m_i-C_R m_j}{m_i+m_j}\mathbf{u}_i+\frac{m_j+C_R m_j}{m_i+m_j}\mathbf{u}_j \\
\mathbf{v}_j=\frac{m_j-C_R m_i}{m_i+m_j}\mathbf{u}_j+\frac{m_i+C_R m_i}{m_i+m_j}\mathbf{u}_i$$
where $m_i,m_j$ are the masses of  of particles $i$ and $j$, respectively. $\mathbf{u}_i, \mathbf{u}_j$ are the velocities before the collision. $\mathbf{v}_i, \mathbf{v}_j$ are the velocities after the collision, and $C_R$ is the coefficient of restitution (a quantity that tells us how elastic or inelastic the collision is $1\to$elastic, $0\to$perfectly inelastic. It could depend on the material and geometry of the bodies that collide, their energies and the "angle" of collision, but the exact value here is not important).
As you can check, these equations satisfy conservation of total momentum and total angular momentum (the sum of the momenta and angular momenta of the bodies before and after the collision is the same)
In the system we are considering, there should be a total angular momentum $\mathbf J \neq \mathbf0$ (if we consider the bodies as having been placed at random positions with random momenta, it is impossible (probability zero) that the angular momentum is zero) and as you can see from the equations above, the angular momentum of each of the particles after the collision is an average of their angular momenta before the collision. After a very large number of collisions (a large number of averages) all of the individual angular momenta should have almost the same direction as $\mathbf J$ (and adding them still gives $\mathbf J$ because the total momentum was conserved for each individual collision).

For the second part of your question, i.e., a model that predicts the time required for this to happen. It seems to me like a complicated thing, it would obviously have to be a numerical simulation, but it will also be heavily dependent on the initial distribution of bodies. 
A good model may need a variable $C_R$ (with its complicated dependencies) and will require a very large number $N$ of bodies with $6N$ first order coupled ODE's (or $3N$ second order), with all the usual issues in numerical simulations (stability, consistency, convergence, stiffness, etc.).
So... a model no... but a numerical simulation... maybe. You would have to search for articles if this has been done before and if not you could be the one to do it.
