Modeling of monolayer charge-stabilized colloidal crystals with static hexagonal crystal lattice

The mathematical model of monolayer colloidal crystals of charged hard spheres in liquid electrolyte is proposed. The particles in the monolayer are arranged into the two-dimensional hexagonal crystal lattice. The model enables finding elastic constants of the crystals from the stress-strain dependencies. The model is based on the nonlinear Poisson-Boltzmann differential equation. The Poisson-Boltzmann equation is solved numerically by the finite element method for any spatial configuration. The model has five geometrical and electrical parameters. The model is used to study the crystal with particles comparable in size with the Debye length of the electrolyte. The first- and second-order elastic constants are found for a broad range of densities. The model crystal turns out to be stable relative to small uniform stretching and shearing. It is also demonstrated that the Cauchy relation is not fulfilled in the crystal. This means that the pair effective interaction of any kind is not sufficient to proper model the elasticity of colloids within the one-component approach.


Introduction
Charge-stabilized colloidal systems are stable suspensions of electrically charged colloidal particles immersed in a liquid electrolyte [1]. In the colloidal crystals, particles are spatially ordered in such a way that they form a crystalline structure. These systems are of great interest as a basis for the photonic crystals and self-organizing structures [2]. Colloidal crystals are also interesting as model systems in the study of ordinary molecular crystals, as well as disordered colloidal systems. Monolayer colloidal systems contain one layer of particles usually near the substrate or at the interface. In the absence of displacement of particles in a direction perpendicular to the layer, such structures can be considered as two-dimensional.
Colloidal crystals of like-charged particles are stable only in the presence of an external pressure. Elastic properties of such systems have some specificity compared to conventional media [3,4]. In particular, the series expansion of the stress with respect to strain includes a constant term corresponding to the stress in the absence of deformation.
The nonlinear Poisson-Boltzmann equation proved to be a good basis for theoretical study of the charge-stabilized colloids. It does not refer to any prescribed effective interaction between colloidal particles. In the present paper, we describe the numerical method based on the Poisson-Boltzmann equation to study elastic properties of the monolayer colloidal crystals. We then use this method to find elastic constants of the model crystals in the broad range of the crystals' densities.

Description of the Model
The model monolayer charge-stabilized colloidal crystal is shown schematically in figure 1. It is a system of identical absolutely rigid spherical particles of radius R located between two parallel plates. The distance from the edge of the particles to the plates on both sides is the same and equal to H . The particles and plates are charged with a constant surface charge density σ and  respectively. The centers of colloidal particles are located at the nodes of a two-dimensional hexagonal crystal lattice with the lattice parameter a . The space between the plates and the particles is filled with electrolyte. The static lattice approximation is used, so that the thermal motion of the colloidal particles is neglected. in the text. The electrostatic potential φ in the electrolyte region obeys the nonlinear differential Poisson-Boltzmann equation [1]. We consider the case of a binary symmetric univalent electrolyte, or 1:1 electrolyte, which has two species with valences 1 1 z  and where 0 ε is the electric constant, ε is the relative dielectric permittivity of the electrolyte, (2) In what follows, the case of (infinitely) large dielectric permittivity of the electrolyte in comparison with the solid phase is considered, that is an appropriate limit case for aqua-based electrolytes. Then, the electric potential in the electrolyte region is independent of the potential inside the particles and plates and is determined only by the magnitude of the charge on their surfaces. This leads to the following Neumann boundary conditions on the surface of the particles: φσ and on the surface of the plates: where n is the unit normal vector directed into the interior of the electrolyte. To find the electric potential, only one unit cell is needed due to the spatial periodicity of the crystal. The domain in the initial configuration is based on the Wigner-Seitz cell of one of the particles. This cell is a prism with a hexagonal base and height   2 HR  . The lateral size of the domain is characterized by the lattice parameter a shown in the figure 1. The space of the colloidal particle in the center of the cell is subtracted from the domain. Only the lower half of the cell is sufficient due to the horizontal mirror symmetry plane passing through the centers of the particles. When a deformation is imposed, the domain is obtained by the corresponding transformation of the initial region.
The homogeneous Neumann boundary condition holds on the plane of mirror symmetry: On the side faces of the domain, the periodic boundary conditions hold for the potential: and for the normal component of the potential gradient: Here m , , are the vectors of primitive translations that separate these opposite faces. The electric potential in any instantaneous crystal configuration, both initial and deformed, is determined by the solution of the boundary value problem for equation (2) in the domain mentioned above with the boundary conditions (3)-(7) on its boundaries.

Numerical Experiment
The elastic constants of the crystal are found from the stress-strain relations obtained in the numerical experiments. We consider only lateral deformations so that the colloidal particles are displaced in directions parallel to the plates with the centers remaining in the middle plane. Then, the system can be treated as a two-dimensional one. The symmetry analysis shows that such a system has three independent second-order elastic constants (elastic moduli), where ε is the strain parameter. In this case, the stress-strain relations are the following: 2 TB ε  (12) for shearing [3]. Here dots denote terms of the orders higher than the first. Three expressions (10)-(12) are sufficient to find all elastic constants of the first and second order. The stress tensor ij T in a two-dimensional crystal is expressed in terms of the usual osmotic stress tensor ij S (in three-dimensional space) as follows: where   2 HR  is the height of the cell and ij δ is the Kronecker delta symbol. Stress b S is the isotropic background stress, which is present in the very dilute system, when the particles are separated by sufficiently large distances and do not interact with each other.
The osmotic stress ij S is calculated using the fundamental stress tensor S is also calculated by expression (14) but for an "empty" system, that is, a system with charged plates but without the particles.
The tensor ij  for equation (2)  T ε are obtained in two numerical experiments described above for every fixed set of the model's parameters σ ,  , R , H and a . Then, the polynomial approximation of these relations is carried out by the standard least square method. The coefficients of the approximation give the elastic constants in accordance with formulas (10)-(12). The numerical values of the strain parameter ε vary in the range 0.01 0.01  for stretching-compression and 0 0.01  for shearing in step of 0.001. At each step of deformation, the new domain is constructed and the corresponding boundary value problem is solved numerically. The solution is carried out by the finite element method using gradient irregular meshes of tetrahedral elements. The calculations were partly supported by the Supercomputer Centre of Lomonosov Moscow State University [5].

Results and Discussion
Elastic constants were calculated for the crystal with the following model parameters: Bp  , tends monotonically to zero together with the grows of the lattice parameter. Behavior of the second-order elastic constants is also monotonic. While the lattice parameter decreases, the curves rise almost exponentially initially. Then, the rate of growth increases for smaller values of the lattice parameter. Positivity of all the three second order elastic constants means that the crystal is stable relative to the uniform stretch and shear strains within the present model. Another set of elastic constants, C-constants, was also calculated. These constants are the coefficients in the series expansion of the second Piola-Kirchhoff stress tensor with respect to the components of the Lagrange strain tensor [6]. B-constants and C-constants are related to each other by the following expressions:  Instead of the theory based on the PB equation, the simpler one-component model [1] can be used to describe the colloidal systems. Within the one-component model, the complex interaction in the system is reduced to the effective interaction between the particles only. The effective potential is often chosen to be the central pair one. Having the C-constants, we can verify, whether the elastic properties of the particular crystal can be described by central pair effective potential of any kind. For pair effective interaction the Cauchy relation CC definitely deviates from unity even for large lattice parameters where the pure pair interaction could be expected. This means that the effective interaction between colloids in the crystal is not purely pairwise, and three-and many-body effective interactions can play a role in the system.

Conclusion
The model of monolayer charge-stabilized colloidal crystals with two-dimensional monatomic hexagonal crystal lattice and constant surface charge density on the particles is proposed. The model is based on the Poisson-Boltzmann nonlinear differential equation. Numerical procedure for obtaining elastic constants from the stress-strain relations is proposed and applied to the crystals with different values of the lattice parameter. It was found that the crystals are stable relative to small uniform strains in the whole range of the lattice parameters studied. The breakdown of the Cauchy relation shows than any purely pair effective potential is not sufficient to proper simulate the elasticity of the colloidal crystals within the one-component model.