We first discuss, within perturbation theory, the scalings of the three contributions to the depletion interaction in Eq.(32). The first order term influences the depth of the attractive (primary) minimum, which is of order , and this is the same scaling as for depletion by small spheres whose diameter is equal to the length of the rods L, Ref.. For a given , however, much larger attractions are possible with rods since the maximum is of order (corresponding to ) which far exceeds that attainable for spheres, of order . The second order contribution in (32), which provides a repulsive barrier, arises directly from the excluded volume interaction among rods, and accordingly is smaller than the first by one power of . The resulting barrier height scales as in the perturbative regime. Finally, the third order correction multiplies a term reduced by one further power of , on which basis we expect the previous estimate of the barrier height to become unreliable for of order unity. In practice, the form of and are such that the range of validity of the second/third order treatment is roughly , depending on the separation h as can be seen in Figures 3 and 4. These show the full numerical results compared with the perturbation treatment, for the cases .
At low concentrations, the numerical results and perturbation treatment agree well, which is expected since the expansion leading to Eq.(16) is then justified. Put differently, the second and especially the third order results provide very powerful checks of our numerical solution; for intermediate concentrations ( ) the effect of the third order correction is quantitatively reproduced by the full results. This makes us confident that our computer programme indeed provides a correct solution of the full self-consistent integral equation. Another check is provided by the limit for which the Onsager pressure is recovered, as mentioned above.
At concentrations above the full results show the second order perturbation theory to somewhat overestimate the main depletion attraction. A more surprising and important feature of the full calculation is the relative smallness of the repulsive barrier now found in the depletion potential for of order unity. For example, in a system of size ratios R/L=10,L/D=20 (reasonable parameters for colloidal sphere-rod mixtures), it was argued in Ref. that a barrier of could be expected at and one of for . In contrast to these predictions, which were based on second order perturbation theory, the full numerical solution presented here shows that the maximum repulsive barrier height in the depletion potential is of order , attained at a reduced density of , giving for the aforementioned parameters a maximum barrier height of . This will now have very little effect on either phase equilibria or kinetic stabilization, although the repulsive features of the potential might still be measurable between parallel plates in a force machine.
As Fig.2 shows, the third order correction has an attractive part largely coinciding with the repulsive part of . This overlap seems to be mainly responsible for the reduction (and perhaps ultimately the disappearance) of the repulsive depletion potential at high densities. Although in our earlier predictions we were aware of the dangers of extrapolating the second order perturbation results beyond their valid density range, it seems very surprising that the discrepancy should already be so large for of order unity. We can find no deep-seated reason for this except perhaps that in going from the force between plates to the interaction energy between particles there are two integrations. This means that the magnitude (once) and the range (twice) of the repulsive force enter multiplicatively; at each factor is perturbatively overestimated by roughly a factor two which leads to an order of magnitude reduction in barrier height. Even this does not explain the unusual smallness of the dimensionless prefactor in the expression for the maximum barrier height ( ).
We now consider an important limitation to our treatment, arising from the finite aspect ratio of the rods. For practical situations the aspect ratio L/D is never more than 100, so the Onsager limit ( ), on which Eq.(14) is based, may not be fully satisfied. A conservative estimate of this error can be obtained from examining the virial series for rods in bulk. Onsager showed that in the bulk the virial expansion of pressure would terminate at the second order in for . For finite aspect ratio rods, the virial coefficients have been computed numerically by Frenkel , and we find that at and L/D=100 this pressure deviation from the Onsager limit is of the same order of magnitude as the maximum repulsive depletion force (which arises for ). This deviation increases with increasing concentration and with decreasing aspect ratio. Accordingly we cannot rule out the possibility that for aspect ratio of (say) 20, as used in the numerical examples discussed above, a significant repulsive barrier could once again be restored by these higher order virial effects. This issue is probably best addressed by computer simulation.
A second limitation is our assumption that the density distribution function possess the uniaxial symmetry imparted by the parallel plate geometry (or the geometry of two spheres in solution). This would of course be untrue if the rods spontaneously broke the uniaxial symmetry to form a biaxial system. This excluded-volume driven transition could occur near that occurring in bulk at the I-N (Isotropic-Nematic) transition ( ); however the role of the plates should help to suppress such a transition, because the restricted geometry partially aligns the rods thereby alleviating the excluded-volume interaction. Accordingly for hard plates (or spheres) without specific surface treatments we do not expect the uniaxial symmetry to be broken way before the bulk I-N transition is reached. Indeed the stability analysis in Ref. suggest that the uniaxial symmetry should be preserved for . There should be some adjustment, however, for rods whose aspect ratio is not large. In the case of uniaxial symmetry being broken, a grid of higher dimension will be needed in the numerical solution, demanding significantly more computing power in solving the Eq.(14).