Strong Effect of Weak Charging in Suspensions of Anisotropic Colloids

Suspensions of hard colloidal particles frequently serve as model systems in studies on fundamental aspects of phase transitions. But often colloidal particles that are considered as ``hard'' are in fact weakly charged. If the colloids are spherical, weak charging has a only a weak effect on the structural properties of the suspension, which can be easily corrected for. However, this does not hold for anisotropic particles. We introduce a model for the interaction potential between charged ellipsoids of revolution (spheroids) based on the Derjaguin approximation of Debye--H\"uckel Theory and present a computer simulation study on aspects of the system's structural properties and phase behaviour. In line with previous experimental observations, we find that even a weak surface charge has a strong impact on the correlation functions. A likewise strong impact is seen on the phase behaviour, in particular, we find stable cubatic order in suspensions of oblate ellipsoids.


Introduction
Colloids are widely used as models to study basic questions of statistical mechanics. In particular, "hard" particles that only interact by volume exclusion have been studied intensively since the 1950s 1 . Hard particle systems are appealing because their phase behaviour is of purely entropic origin and they can easily be treated by computer simulation. For example, in the context of liquid crystals, studies of hard ellipsoids, spherocylinders and platelets have provided valuable insight into the basic phase transition mechanisms [2][3][4][5] .
In one of the first computer simulation studies of a phase diagram of hard, oblate particles, Veerman and Frenkel 5 observed that the particles arranged parallely in stacks which in turn formed a suprastructure of perpendicular orientations. They named this phase "cubatic". The existence of the cubatic phase has since been under heated debate, and recently several simulation studies [6][7][8] showed that the cubatic is always metastable with respect to either the isotropic or the columnar phase for various round hard platelet models (i.e. platelets with circular cross section). In contrast, the cubatic phase is stable for square plates 9 . Experimentally cubatic order has recently been detected in dispersions of hexagonal, charged plate-like particles 10 . We will show below that for oblate ellipsoids cubatic order becomes stabilized as soon as there is a small surface charge.
Over the past years, experimental methods to synthesize and characterize suspensions of ellipsoidal colloids have been advanced and theoretical predictions have been tested experimentally [11][12][13][14][15][16] . In 2011, Cohen et al. 17 measured the structural properties of a PMMA ellipsoid system and showed significant differences when comparing their data to theoretical predictions for hard ellipsoids given by Percus-Yevick theory 18 and simulations 19 . In the following, we will test our theoretical treatment of weakly charged ellipsoids against the results of this experimental study. We will show that weak charging, as it is often present in PMMA-colloid suspensions, changes the pair correlations such that they match those observed experimentally and alters the phase diagram considerably.
The paper is organized as follows. We first derive the interaction potential. Then we present simulation results on the positional correlations in the system and compare them to the experimental results of ref. 17 . Finally we present a scan through the phase diagram for changing surface charge density and show that the cubatic phase is stabilized.

Derivation of the interaction potential
The interaction of weakly charged colloids in an electrolyte suspension can be treated in an adiabatic fashion: one assumes that co-and counterions instantaneously readjust upon a change in the colloidal positions, giving rise to an effective interaction between the colloids, possibly of multi-body nature. If the Debye-Hückel screening length is smaller than the extensions of the colloid, then this effective potential can be well approximated by a sum of two-body terms. For the interaction between two colloids, we use the Debye-Hückel approximation (linearized Poisson-Boltzmann (PB) theory). For two spheroids (ellipsoids with one rotational symmetry axis), the effective potential depends on four variables (the center-to-center distance and three angles) such that an explicit tabulation of the PB solutions, let alone the determination of PB solutions "on the fly" in a simulation code appears forbidding 20,21 . Here, we resort to the venerable Derjaguin approximation which has been often used to calculate effective colloid-wall or colloid-colloid interactions in the literature but we are not aware of its practical use in further simulation or theoretical studies of concentrated solutions involving anisotropic particles. The Derjaguin approximation rests upon the following argument: Suppose the free energy of the interaction between two planar walls is known, and its density will be denoted by f (h) where h is the distance between the walls. The interaction potential between two convex bodies (of the same type as the walls) can be approximated by just integrating over geometrically opposing area elements dA (at distance h) where the free energy of interaction between the area elements is given by the wall free energy f (h)dA. The mathematical elaboration of this approxi-mation is given in App. A and results in the following expression for the free energy of interaction F(H 0 ) between two convex bodies with the minimal distance H 0 between their surfaces Here, the product εε is given by It involves the principal curvatures ε i , ε i of the surface of body i = 1, 2 in the planes tangential to the distance vector between the bodies, and also the angle ω between the coordinate systems in the tangential planes with coordinate axes given by the directions of the principal curvatures. For these geometric definitions, see Fig. 1.

Charged ellipsoids
We will apply these ideas to the interaction between hard, charged ellipsoids (spheroids) with main axes a and b where b is the main axis in the plane perpendicular to the rotational symmetry axis. The aspect ratio is given by t = a/b. In Debye-Hückel approximation, the electrostatic potential ψ fulfills where κ −1 is the Debye-Hückel screening length. For a charged wall with charge density σ , the solution is with the wall contact potential φ 0 = σ /(ε s κ) (ε s is the dielectric constant of the solvent). We approximate the solution for two charged walls at distance h by and the pressure (force density per unit area between the plates)f 2w is obtained most easily by evaluating the stress tensor at the midplane z = h/2 which has there only a contribution from the ion osmotic pressure: The free energy density f is then found through integration Using this, the Derjaguin free energy (1) becomes The Derjaguin free energy decays exponentially with H 0 , as expected for the Debye-Hückel approximation. This decay is fast enough that the approximation is accurate enough for practical purposes, see App. B for a discussion. Note that the anisotropy in the free energy has two sources: H 0 depends on the different orientations as well as the curvature term 1/ √ εε . The latter one has a strong influence on the interaction of oblate ellipsoids (see below). We write the Derjaguin free energy as Note that V (H 0 ) is dimensionless in the last equation. The prefactor (with dimension of energy) in Eq. (9), V 0 = 2σ 2 b/(ε s κ 2 ), contains the charge density σ as a parameter. It is advantageous to introduce the dimensionless charge densityσ = σ eβ /(κε s ) (where e is the elementary charge and β = 1/(k B T ) is the inverse temperature) as well as the Bjerrum length of the solvent, λ B = β e 2 /(4πε s ). With these definitions the prefactor becomes Below, in the comparison with the experiments of ref. 17 , we will treatσ as a fitting parameter. Regarding the interpretation of this value one should keep in mind that it is an effective or renormalized charge density. Forσ 1, the bare and the renormalized charge density are approximately the same, whereas in the limit of very large bare charge densities the renormalized charge density approaches a constant,σ → 4 22 . Thus the approach is consistent only forσ < 4.

Numerical implementation
We computed numerically the dimensionless, exponentiated free energy exp(−V (H 0 )) (Eq. 9) on a 4-dimensional grid with axes characterizing the relative configurational state of two ellipsoids (center-to-center distance and three angles). For each configuration, the minimal distance H 0 was determined by a conjugate gradient routine and the radii and directions of principal curvature were determined through the first and second fundamental form of the ellipsoid surfaces at the two points O and O (whose distance is H 0 , see Fig. 1). In the Monte Carlo simulations (see below), the such tabulated free energy was used as the acting potential between pairs of ellipsoids, together with linear interpolation to determine the potential at off-grid values of the variables characterizing the relative configurational state. We also tested a further approximation to the Derjaguin free energy in which the curvature term 1/ √ εε is replaced by a constant, the average radius of the ellipsoid. Then V (H 0 ) only depends on the minimal distance H 0 , which can be well approximated by an extension of the Perram-Wertheim routine 23 frequently used for checking overlap of hard ellipsoids (see App. C). In this way, the potential can be determined "on the fly", and it works reasonably well for aspect ratios 0.8 t 2. Note, however, that the short-range anisotropy of the potential increases rapidly with the aspect ratio becoming small. For oblate ellipsoids (t < 1, disk-like particles), the ratio between the potential at contact in side-side configuration (flat sides of the disk touching) and in edge-edge configuration (rims of the disk are touching) is (t 2 + 1)/(2t 3 ) and thus scales for small t as 1/t 3 . Therefore, this further approximation to the Derjaguin free energy is not applicable to flat oblates.

Simulation Results
We have carried out Monte Carlo simulations at constant temperature T, constant number of particles N and volume V with periodic boundary conditions, and computed equilibrium structural properties of the system as a function of the effective surface charge density and the packing fraction. The particle number N ranged from 3000 to 3840, T was set to 300K.
First, we discuss the structure of the isotropic phase. This part of our work has been inspired by recent experimental measurements of the radial, orientation-averaged pair correlation function g(r) in suspensions of prolate ellipsoids (aspect ratio t = 1.6, a = 3.2 µm, b = 2.0 µm) and oblate ellipsoids (t ≈ 0.25, a ≈ 0.96 µm b ≈ 3.8 µm, with a considerable experimental uncertainty on the polydispersity and thus on the value of t) 17 .
The structural correlations that are presented in ref. 17 are much stronger than one would expect for a system of hard ellipsoids, and this was taken as an indication that on the theory side, the correlations in suspensions of ellipsoids are not sufficiently well understood. However, the experimental suspension was additionally stabilized by a surfactant which introduced a small amount of charge on the particles. In a later study 12 , the authors investigated the influence of charge on g(r) for the prolate particles by simulation and found it non-negligible: with a small charge density of σ ≈ 9 e/µm 2 (where e = 1.6 · 10 −19 C is the elementary charge), distributed on particles modelled by an assembly of three cut spheres to approximate the shape of the ellipsoids, the experimental g(r) could be reproduced. No corresponding results for the oblate particles have been reported, though.
In order to model the experiment of ref. 17 , we set λ B ≈ 22 nm appropriate for a solvent with an average dielectric constant of ε s = 2.5 and κ −1 = 0.3 µm for the Debye-Hückel screening length. We only vary the effective surface charge density to reproduce the experimental data. Fig. 2(top) shows the radial positional distribution function g(r) of ellipsoids (circles) for an aspect ratio t = 1.6 and a packing fraction φ = 0.31. The simulation data perfectly match the experimental data (triangles, from Fig. 2 in ref. 17 ). The corresponding dimensionless charge density is given byσ = 0.83, i. e. the effective charge density is σ = 10 e/µm 2 . This value is reasonable for the experimental system used in ref. 17 , it is in good agreement with the deduced effective charge on the colloids of ref. 12 , and it justifies in retrospect the assumption of the Debye-Hückel approximation used to derive the interaction potential. (Note that the modelling of the electrostatic particle interactions in ref. 12 approximately corresponds to the Derjaguin free energy with curvature neglected (see Sec. 1.2 above) which works well for the moderate aspect ratio of 1.6 but not for particles with larger curvatures. This rationalizes the good agreement between our results and those of ref. 12 for σ and g(r).) The second data set discussed in ref. 17 has been measured in a more dilute suspension of prolate ellipsoids at a packing fraction φ = 0.26. Using the same system parameters for λ B , κ −1 and even σ as for φ = 0.31, we again obtain very good agreement of the radial distribution functions, see Fig. 2(bottom). The data set for hard ellipsoids from ref. 17 is presented as well for comparison.
The last radial distribution function that is presented in ref. 17 , was measured in a suspension of oblate ellipsoids of an aspect ratio of "t ≈ 0.25" (with a larger polydispersity than in the prolate case). As explained in Sec. 1.2, the curvature around the rim of the particles in this case is important for the electrostatic interactions, hence the approach of ref. 12 could not be applied here. Fig. 3 shows our simulation results for oblate ellipsoids. We set again the same value for λ B , κ −1 and σ . The experimental and theoretical data agree reasonably well given the uncertainty of the aspect ratio of the experimental system.  To conclude this section, we validated the Derjaguin approximation for the electrostatic interaction of charged ellipsoids. A small amount of surface charge has a strong influence on the pair correlations, due to the small dielectric constant (large Bjerrum length) of the solvent. The short-range anisotropy of the electrostatic interaction is especially important for oblate ellipsoids.

Impact of the surface charges on the nematic phase
We now consider oblate ellipsoids of aspect ratio t = 0.25 at a packing fraction of φ = 0.48. In suspensions of hard ellipsoids the nematic phase is located at packing fractions φ > 0.4 2 . We study the effect of increasing surface charge density, ranging from σ = 0.00002 e/µm 2 to σ = 2.0 e/µm 2 , on the structural properties of the liquid. Note that these surface charge densities are even lower than the value that was needed to reproduce the experimental findings of ref. 17 . The Bjerrum length λ B and the Debye-Hückel screening length κ −1 are not modified with respect to the previous section.
As the initial configuration of a first set of simulation runs we used an fcc crystal in which the ellipsoids were oriented in parallel. We let the system relax until its energy had reached a stable value. In the case of perfect charge screening (i.e. almost hard ellipsoids) the system relaxed into the expected nematic phase, see Fig. 4. A similar degree of nematic ordering formed for a surface charge density of σ = 0.02 e/µm 2 . In contrast we observe qualitatively different behaviour for surface charge densities σ = 0.2 e/µm 2 and σ = 2.0 e/µm 2 , see snapshots in Fig. 5 and Fig. 6. Note that these configurations evolved from an initial fcc configuration with parallel orientation of the ellipsoids. Fig. 7 shows g(r) for different surface charge densities σ at φ = 0.48. In the nematic phase, g(r) has a cusp at a distance r that corresponds to in-plane rim-rim configurations. In contrast, at σ = 0.2 e/µm 2 there is pronounced positional order with peak positions at multiples of the length of the small axis a plus a small distance to account for the electrostatic repulsion. Fig. 8 shows the orientational distribution function g 2 (r) = 1 g(r) where u i is the unit vector along the axis of particle i and the average is over all pairs of particles and the canonical ensemble. At small σ g 2 (r) decays smoothly to a nonvanishing value at large distances, which is characteristic for nematic ordering. At σ = 0.2 e/µm 2 and above there is parallel order at short distances and random orientation at large distances, i.e g 2 (r) decays to zero. The first "perpendicular peak" with g 2 < 0 appears at the distance that corresponds to a configuration in which the rim of one ellipsoid points to the pole of the other. This peak is superposed with the second layer of parallel stacking in g(r). We conclude that stacks of ellipsoids are arranged perpendicular to each other to form cubatic order.
To test whether this phase is metastable, we then initialized simulations in the nematic phase, the columnar phase and a perfect long-range cubatic phase. All runs equilibrated into the phase discussed above, thus we conclude that it is the most stable.

Conclusion
We have added surface charges to hard ellipsoids and treated them numerically using the Derjaguin approximation. With the Derjaguin approximation, the short-range anisotropy of the electrostatic interactions is captured quantitatively correctly as long as the screening length is smaller than the extensions of the particles. Even very small surface charges, as they often are present in experiments on "hard" particles have a strong effect on the structure of the suspension.
We showed for charged oblate ellipsoids that a phase of perpendicularly oriented short stacks (a cubatic) is thermodynamically more stable than the nematic phase, which is stable for uncharged ellipsoids at the same packing fraction. The cubatic phase is also more stable than the crystal and the columnar phase. It should be accessible to experiments on suspensions of PMMA ellipsoids.
The Derjaguin approximation 24 involves the following steps: Suppose the free energy of the interaction between two planar walls is known, and its density will be denoted by f (h) where h is the distance between the walls. For calculating the interaction potential between two convex bodies (of the same type as the walls) one determines the minimal distance H 0 between the two surfaces and the tangential planes (see Fig. 1). The free energy of interaction between the two bodies is approximated by where the integral runs over (one of) the tangential planes and z 1 , z 2 are the distances of the point on surface i described by z i (x, y) to their respective tangential plane (see Fig. 1). If f (h) quickly decays to zero, it is safe to integrate over the whole plane. This integral is greatly simplified if we approximate the surfaces by quadratic forms around the points O and O respectively: where ε 1 , ε 1 are the principal curvatures of surface 1 at point O and x 1 , y 1 are coordinates in the tangential plane in the direction of the principal curvatures. Likewise for surface 2. The directions of the principal curvatures of surface 1 and 2 do not agree but include an angle ω, thus: The distance z 1 + z 2 becomes This distance is a quadratic form. We may perform a rotation to another coordinate system x, y where the off-diagonal matrix elements become zero, i.e.
Using that result, the free energy becomes Introducing new coordinates r, φ via x = r cos φ / √ ε and y = r sin φ / √ ε : where the second line (which is Eq. (1)) follows from the substitution h = H 0 + r 2 /2. The force (in direction of OO ) between the two bodies is just given by The product εε is the determinant of the matrix in Eq. (20) and must be equal to the determinant of the matrix in Eq. (16). Thus we find:

B Validity of the Derjaguin approximation
We can estimate the validity of the Derjaguin approximation by considering the example of two interacting, charged spheres with radius r 0 and charge density σ at center distance d for which we can compare the "exact" Debye-Hückel result with the corresponding Derjaguin approximated result. In Debye-Hückel approximation, the potential of a single sphere is given by and the charge Q eff is determined through the boundary condition ∂ Φ/∂ r| r=r 0 = −σ /ε s , giving In superposition approximation (as before), the interaction free energy of two such spheres is given by where H 0 = d − 2r 0 is the minimal surface-to-surface distance. The Derjaguin approximated free energy follows from Eq. (8) using εε = 4/r 2 0 : It is the limit κr 0 1, H 0 r 0 of Eq. (28). The ratio is given by From this equation one sees that the Derjaguin approximation produces a large relative error for H 0 2r 0 . Appropriate for the examples studied in Sec. 2, we set r 0 = 1 µm and the screening length κ −1 = 300 nm. Thus, κr 0 ≈ 1/3, and the Derjaguin approximation overestimates the free energy by a factor 3.4 at H 0 = 2r 0 = 2 µm. However, at that distance the potential has dropped by a factor exp(−2κr 0 ) ≈ 0.0013 compared to its value at contact, so the free energy itself at that distance is small and likewise the absolute error the Derjaguin approximation produces is small. The effect is presumably negligible in our simulations which are sensitive to the short-range behaviour of the effective free energy between the particles.

C Perram-Wertheim approximation for the minimal distance
Perram and Wertheim 23 developed a criterion to check for overlap of two ellipsoids. The algorithm can also be used to compute an approximate minimal distance. We summarize this approach as follows: Suppose we have two ellipsoids with half axes a 1 , a 2 , a 3 and corresponding axis orientation vectors (of unit length) in a lab-fixed coordinate system u 1 , u 2 , u 3 (ellipsoid 1) and v 1 , v 2 , v 3 (ellipsoid 2). One defines two matrices: Let r a and r b denote the center positions of ellipsoid 1 and 2, respectively. We define quadratic forms and its minimum, depending on λ ∈ [0, 1]: For λ = 0, r(0) = r b and for λ = 1, r(1) = r a . Thus as λ varies from 0 to 1, then r(λ ) moves from the center of ellipsoid 2 to the center of ellipsoid 1. If the two ellipsoids do not overlap, then there exists a particular λ for sure for which r(λ ) is outside both ellipsoids and F(r(λ ), λ ) > 1 there. If the ellipsoids overlap, then F(r, λ ) < 1 in the overlap region and thus for each λ the minimal point r(λ ) can not lie outside both ellipsoids since there F(r, λ ) > 1. Therefore a useful overlap criterion is formulated with introducing s = max λ ∈(0,1) F(r(λ ), λ ) This criterion is convenient to use due to the explicit form for F(r(λ ), λ ) which Perram and Wertheim provide 23 : The maximization needed in Eq. (37) has to be done numerically, though. The value of s can also be used to calculate an approximative minimal distance through The interpretation, according to Paramonov and Yaliraki 25 , is as follows. First, the geometrical meaning of the λ maximization in Eq. (37) becomes clear by looking at The term ∇F is zero by virtue of the definition of F(r(λ ), λ ) in Eq. (36) and thus the derivative above is zero when s = F A (r(λ max ), λ max ) = F B (r(λ max ), λ max ). This, however, describes the tangential contact between ellipsoids with scaled half-axes √ sa 1 , √ sa 2 , √ sa 3 . (r(λ max ) is the tangential contact point since ∇F A ||∇F B there.) Then the points s a , s b defined by (likewise for a → b) lie on the surface of ellipsoid 1 and 2, respectively. We see that s b − s a is parallel to the center distance vector r b − r a and that d PW = |s b − s a |. Thus d PW is the minimal directional distance between the ellipsoids (i.e. minimal distance between two points on the surfaces of ellipsoid 1 and 2 in the direction of the center distance vector. This approximation is not exact but as presented in Fig. 9 the error is small and hence the approach using the simple determination of the Perram-Wertheim distance justified for the calculation.
The relation is not one to one because of the different relative orientations but there is agreement for moderate aspect ratios. Fig. 9 For given configurations of two prolates of aspect ratio t = 1.6, we plot the minimal distance relative to the Perram Wertheim distance.