Random N-body problem Suppose there are $N$ unit-mass particles whose initial positions
are uniformly distributed in a unit-radius disk.
Each particle is assigned a randomly oriented initial velocity vector $v_i$ of length $1$
(green below).
Added: Robert Israel's incisive comment suggests that I should also stipulate
that $\sum_i v_i = 0$.
Then the particles act upon one another via inverse-square gravity.

          

     

Dots show initial positions inside unit disk.
Green vectors: initial velocity; red vectors: final velocity.



Q. What is the probability that $k \ge 1$ of the $N$ particles remain within
  a disk of some radius $R \ge 1$ forever?

In the illustration above, $N=8$ and $R=3$.
(But I do not trust my crude simulations.)
I am wondering if the answer is: zero, independent of $k$ and $R$ and the 
gravitational constant?
Perhaps the answer differs for points in $\mathbb{R}^2$ vs. points in $\mathbb{R}^3$ (or $\mathbb{R}^d$, $d > 3$)?
 A: 
What is the probability that $k \ge 1$ of the $N$ particles remain within
  a disk of some radius $R \ge 1$ forever?

If the gravitational constant is small enough, then the probability is definitely less than one, since it's easy for one of the particles to be initially headed away from the others at a high enough speed to escape.
For $N\le 3$, the probability appears to be nonzero given a high enough gravitational constant. For $N=1$, the probability is 1. For $N=2$, we have Keplerian orbits. For $N=3$, there are figure-eight configurations known as 3-body choreographies that are believed to be stable in the KAM sense, and there seem to be solutions of this type for a variety of potentials, including $-1/r$. Stability against small perturbations means that you have bound states that occupy a region of phase space with nonvanishing volume, which implies nonzero probability.
For $N\ge4$, this appears to be an open problem. In general, we expect such systems to be unbound for thermodynamic reasons. This is because for large $N$ the probability of a given state for a particular particle falls off like $e^{-E/T}$, where $E$ is the kinetic energy and $T$ is the temperature. Since there is no bound on $E$, any system like this is generically expected to evaporate its particles off into the surrounding space. This is really just a phase space argument. There is infinite phase space out there at large distances. Such evaporation is in fact observed in galaxies and globular clusters, where we see stars shooting off with anomalously high velocities.
So all such systems should be considered guilty until proven innocent, i.e., we expect them to be unbound by default unless we have some tricky way of constructing a specific example and proving that it's bound. (Numerical simulation doesn't work, because these systems are normally chaotic.) To prove that the probability of being bound is nonzero, we also need to prove that it remains bound under a small perturbation. For $N\ge 4$, I did not come across any mention of any examples that are known to be stable, but I could certainly be missing relevant work.
By the way, physically the "right" way to state the problem in two dimensions is probably to use a $1/r$ force, not a $1/r^2$ force, since we expect Gauss's law to hold for fundamental reasons.
Gerhard Paseman says:

I bet cosmologists thought about this problem in considering origins of the universe.

This would be more about the far future than the distant past. The relevant tool is general relativity, not Newtonian mechanics. The thermodynamic equilibrium state of a gravitationally interacting system according to classical GR is basically a black hole, although it can be more complicated than that depending on the large-scale geometry and topology. You don't get stable orbiting states with a purely gravitational interaction, because gravitational radiation sucks energy away.
A: The probability of a bounded motion is neither 0 nor 1. (I am assuming that the initial total momentum is 0). Indeed, for the case of two bodies, we have an explicit solution, and it shows that there is an inequality
for the initial data which distinguished bounded motion from unbounded one.
For arbitrary number of bodies, there is an open set of initial positions when the motion is bounded (think of something like a solar system), and also an open set
of initial positions which give unbounded motion. 
A: Update: There's nothing obviously preventing (negative) energy from bleeding from a pair of close bodies to others, so this answer is wrong as stated.
At least for sufficiently large $N$, the case of equal masses reduces to the case of unequal masses to give a positive probability of bounded motion.  This is because two bodies that are sufficiently close together relative to other masses can act as a body of twice the mass, since they will be energetically bound.  Thus if one can build a solar system out of unequal masses, one can build a stable solar system out of equal masses.
Edit: Since two masses are stable, there is a positive probability of stability for any $N$ by iteratively replacing one mass by two. 
