Q:

How to quantify the mixing of two hard-spheres gases?

Suppose that there are two types of hard-spheres gases ($a$ and $b$) in a box. Suppose that their radius is much smaller than the box's characteristic size. These two gases tend to repel each other. You want to compute the Entropy of mixing of these gases. The following examples constitute two limiting situations which should make the problem clear.
Case 1:
If the inter-gas repulsion tends to zero, both gases will tend to spread over the entire box and will be perfectly mixed.
The tricky aspect is when the inter-gas repulsion is present. In this circumstance, the two gases will not completely mix. For example, there will be regions where one gas will form a cluster and where the other gas will be scarcely present. 
Case 2:
In the limit of large inter-gas repulsion, it could happen that gas $a$ occupies the left part of the box while gas $b$ occupies the right part of the box.
Clearly, since both gases are made up of hard spheres, the volumes that they occupy are always the same, no matter how much the inter-gas repulsion is. What changes, according to the inter-gas repulsion, is their degree of miscibility. Therefore the traditional formula for the Entropy of mixing don't work.
My target is to find an indicator which effectively measures the degree of mixing of these two gases. My idea is to start from the already enstablished concept of "Entropy of mixing" and modify it. Moreover, I was thinking of subdividing the total volume of the box in many sub-regions $j=1,2,\dots,M$. In each sub-region, I can count the number of $a-$ spheres, $a_j$ and the number of $b-$ spheres, $b_j$. (Let's not talk, for the moment, about the optimal size of these sub-regions).
The difficult thing is to come up with a function
$$
  S=S(a_1,\dots,a_M,b_1,\dots,b_M)
$$
capable of quantifying the degree of mixing in a consistent and reasonable way. 
For example, the following properties must be matched by my indicator:

$S$ must be zero when the inter-gas repulsion is zero, i.e. when the two gases are perfectly mixed, i.e. when $a_i=a_j$, and $b_i=b_j$ $\forall i,j$.
$S$ must be maximum (but finite) when the inter-gas repulsion tends to infinite. In this circumstance, the species-separation will be perfect, i.e. $a_i=0$ if $b_i>0$ and, viceversa, $b_i=0$ if $a_i>0$. This because, in this limit, the two gases tend to exclude each other.
$S$ should not be affected by the relative number of hard spheres in the two gases. E.g. If the inter-gas repulsion is zero, you can have perfect mixing, that means $S=0$, both if the two gases have the same number of hard spheres (in this case you have that $a_i=b_i$, $\forall i$ ) and if one gas constitutes an impurity with respect to the other (in this case you have that $a_i=C_a$, $\forall i$  and $b_i=C_b$ $\forall i$, but $C_b\ll C_a$).
If $S$ is computed as an average over the different sub-regions $j$ (meaning that in each sub-region $j$, you preliminarly compute $S_j$), this average should be weighted in a reasonable and consistent way. For example: $S_j$ has a bigger weight if the total number of particles in region $j$ is bigger. 

A:

I'll give an answer that summarizes our current understanding of the statistical thermodynamics of hard sphere binary mixtures. It probably isn't exactly what you were asking for, but hopefully will be (at least) useful background information.
Hard spheres and hard sphere mixtures have been extensively studied theoretically, notably by integral equation methods and "experimentally" by Monte Carlo simulation. One of the main integral equations, the Percus-Yevick approximation, can be solved exactly for hard spheres. This leads directly to the radial distribution function $g(r)$ and, indirectly, to the free energy and other thermodynamic properties. Typically you get the pressure $P$ as a function of density, and can integrate this to get the free energy. The famous Carnahan-Starling equation of state (pressure as a function of density) is fairly accurate over the whole fluid range. Specifically it gives
$$
\frac{F^{\text{ex}}}{Nk_BT} = -\frac{S^{\text{ex}}}{Nk_B} = 
\frac{\eta(4-3\eta)}{(1-\eta)^2},\qquad
\text{where}\qquad
\eta=\pi\rho d^3/6
$$
where $\rho$ is the number density and $d$ the sphere diameter.
Notice that this is the excess entropy: you simply add it to the ideal gas entropy to get the total entropy. There is no "excess internal energy" for hard spheres.
Computer simulations have been used to check the accuracy of this equation,
and have led to small improvements to the formula.
The story for binary mixtures of hard spheres is similar. One can solve the integral equations to get the various radial distribution functions $g_{aa}(r)$,
$g_{bb}(r)$ and $g_{ab}(r)$. These are already relevant to your question, because by integrating, say, $g_{aa}(r)$ you get the number of neighbours of type $a$ surrounding an atom of type $a$, within a desired distance $R$ (and similarly for $b$ around $a$ etc). Also, because they represent density-density correlations, I believe they can be used to calculate the fluctuations of the numbers of atoms of each type within a macroscopic volume, something like $\langle N_a^2\rangle - \langle N_a\rangle^2$, which is close
to what you want. The analogous expression for the equation of state is attributed to Boublik, J Chem Phys, 53, 47 (1970) and to Mansoori, Carnahan, Starling and Leland, J Chem Phys, 54, 1523 (1971) (see also here) and, of course,
there are slightly improved versions. From the pressure, you can get $F^{\text{ex}}/k_BT$, which is equal to $S^{\text{ex}}/k_B$, as a function of
overall density and composition. This is added to the ideal entropy term, which
of course includes the ideal entropy of mixing.
Concerning your case 2, the possibility of demixing into two phases (one rich in $a$ and the other rich in $b$) I believe that the current understanding is that this never happens within the fluid phase of additive binary hard sphere mixtures. (Additive simply means that $a$ interacts with $b$ with a hard sphere interaction of diameter $\frac{1}{2}(d_a+d_b)$). The excess entropy is never strong enough to overcome the ideal entropy of mixing (which, of course, always favours a single disordered phase). There was a burst of interest in this dating from 1991, when Biben and Hansen, Phys Rev Lett, 66, 2215 (1991) carried out more accurate integral equation calculations for the case where the hard-sphere diameters were very different from each other. This suggested that demixing might occur due to something called the depletion force. However, subsequent theories and simulations suggest that this doesn't happen. For an open-access paper giving a lot of background on this, see Lopez de Haro et al,  J Chem Phys, 138, 161104 (2013).

