The origin of interparticle potential of electrorheological fluids

The particles of electrorheological fluids can be modelled as dielectric spheres (DS) immersed in a continuum dielectric. When an external field is applied, polarization charges are induced on the surfaces of the spheres and can be represented as point dipoles placed in the centres of the spheres. When the DSs are close to each other, the induced charge distributions are distorted by the electric field of the neighbouring DSs. This is the origin of the interaction potential between the DSs. The calculation of this energy is very time consuming, therefore, the DS model cannot be used in molecular simulations. In this paper, we show that the interaction between the point dipoles appropriately approximates the interaction of DSs. The polarizable point dipole model provides better results, but this model is not pair-wise additive, so it is not that practical in particle simulations.


Introduction
Electrorheological (ER) fluids are suspensions of fine non-conducting solid particles (up to 50 µm diameter) in an electrically insulating liquid. The dielectric permittivity of dispersed particles is usually higher than that of the carrier liquid [1]. The rheological properties of ER fluids are controllable by the application of an electric field [2]. The apparent viscosity of an ER fluid increases abruptly by the application of a strong electric field of the order of kilovolts per millimeter [silica particles (SiO 2 ) dispersed in silicone oil is a typical electrorheological fluid]. The electric field causes a reversible change in the viscosity. The increase of the apparent viscosity is caused by the chain and column formation of the grains carrying electric dipole moments induced by the external field [1,3]. Due to the electric-field-induced aggregation, the dielectric properties of the ER fluids are also changed [4].
Electrorheological particles are beyond the molecular scale, therefore, their modelling necessarily includes some coarse graining. Coarse-graining means that the many-atom system is modelled by averaging certain degrees of freedom into a response function. The atoms of the ER grains are not modelled explicitly, instead, their dielectric response is taken into account by a dielectric continuum characterized by a certain dielectric constant. The same is performed for the carrier liquid, but with a different dielectric constant.
Therefore, in this coarse-grained but realistic picture, the particle is modelled as a dielectric sphere (DS) immersed in a dielectric continuum. The two interacting particles carry three-dimensional surface charge distributions on their surfaces induced by the external electric field and the electric field exerted by the other particles. Computation of this induced charge distribution is a non-trivial and time consuming process. Therefore, this model is not feasible in computer simulations. Thus, a more simplified model is needed for simulations, possibly a model with pair-wise additive interactions. Such a model for the particles of ER fluids is a sphere carrying a point dipole in its center (we refer to this model as DD). The dipoles are induced by an external electric field E, all being aligned in its direction. The lowest energy configuration of two dipoles are the head-to-tail position where the two dipoles are aligned in the same direction along the same line. Therefore, as we have mentioned, the particles in electrorheological fluids form chains in the presence of an external field. This chain formation is responsible for the increase of viscosity when an electric field is applied.
Simulation study of chain formation of ER fluids is based on different models. Klingenberg et al. [5] used the interaction of single point dipoles with hard core repulsion while Bonnecaze and Brady [6] used a sophisticated polarization model for the dipole-dipole interaction. Chain formation was also found in fluids of particles carrying strong permanent dipoles [7,8]. The orientation of these dipoles, however, was not restricted. This distinction makes that case different from ER fluids, where the dipoles are induced dipoles with directions fixed by the external field.
The ideal point dipole is clearly an approximation to the charge distribution of the more realistic DS model. The DD model ignores the fact that the spheres are polarized not only by the external field but also by the other spheres. This effect can be taken into account by the polarizable dipole (PD) model, which is a step further but still an approximation to the DS model. The PD model is not pair-wise additive, but still feasible in simulations because we have to compute the potentials/fields only at the particle centers rather than on the whole particle surfaces.
The PD model was also used in our simulation study [9] of the correction to the Clausius-Mosotti equation describing the dielectric constant of non-polar fluids. In this work, the non-polar particles were also polarized by a uniform external field. The analogy with electrorheological fluids is unmistakable.
In this paper, we study the DS model and compare its energetics with that given by the DD and PD models by computing the interaction energy between two spheres using all the three models. We conclude that the PD model is an excellent, while the DD model is a reliable approximation to the DS model.

The dielectric sphere model
The dielectric constant inside the sphere, ǫ i , is different from that outside the sphere, ǫ w , (the subscript i refers to an ER particle species here). Then, a dielectric boundary is formed at the surface of the sphere. The external electric field induces a surface charge distribution on this boundary. Since there is no free charge inside the sphere (it is neutral), the net induced charge is zero according to Gauss's law. Thus, the separation of charge on the surface of the sphere corresponds to a dipole-like distribution. Using basic electrostatics in terms of Legendre polynomials [10,11], the dipole moment can be computed as where a i is the radius of the sphere and is the polarizability of the sphere. When the external field is uniform, this is an exact solution: the electric field outside the sphere is equal to the electric field of a point dipole in the center of the sphere.
These charges on the surface of the spheres are not free charges, but rather they are bound (induced) charges. This means that they do not get there from some external circuit, because they have always been there. Their appearance is due to polarization: an electric field separates positive and negative charges of the dielectric. When we compute the energy of a macroscopic dielectric system, the interaction between bound charges does not appear in the formulation. The total energy of the system is the work done against electric field as we bring in the free charges from infinity, namely the work needed to build up the charge distribution in a charge-up process. If we denote the distribution of free charges in the system by q(r), and that of the induced charges by h(r), then this work is computed as follows: 3)

43002-2
where Ψ(r) = Ψ q (r) + Ψ h (r) contains the potentials produced both by the free charges and by the induced (bound) charges Equations use Gaussian units. If the free charges are point charges, the integral in equation (2.5) becomes a sum. If the dielectric boundary is sharp, the induced charge distribution is a surface charge, so the integral becomes a surface integral. In an experiment, the external electric field is produced by surface charges σ 1 and σ 2 on the plates of a capacitor. In this case, the integral in equation (2.4) also becomes a surface integral.
In the above equations, the free chargefree charge interaction and the free chargeinduced charge interaction appear. The induced chargeinduced charge interaction is missing. If we write up the energy as the sum of the interactions between all charges (including free and induced charges), an additional term has to be added. This is the work involved in stretching the dielectric molecules, namely, the work necessary to polarize the dielectric. This work is equal in magnitude and opposite in sign to the induced chargeinduced charge interaction, so they cancel. That is why only the energy of the free charges appears in the equation for the total electrostatic energy of the system [see equation (2.3)].
This result is counter-intuitive: one might think that we have to take into account the interaction between all charges in the system. Moreover, it is counter-intuitive regarding the dipoles of the ER particles. In the DD approach, the interaction energy is computed from the interaction of the dipoles. The dipoles are interpretations of the induced charges. This stands a paradox: why do we take into account the direct interaction between the charge distributions of the particles in the DD approach and why do not we take it into account in the DS approach? A goal of this paper is to resolve the apparent conflict between the two approaches. We will shed some light on the mechanism of the interaction between the ER particles.
The induced charge is calculated by the Induced Charge Computation (ICC) method [12][13][14]. This is a boundary element method where the dielectric boundary surface is divided into surface elements. Poisson's equation is transformed into an integral equation where the unknown variable is the discretized induced charge treated as constant on a given surface element. These charges are included in a column vector h. This vector can be computed from a matrix-vector multiplication h = A −1 c, (2.6) where vector c contains the normal components of the electric field in the centers of the tiles and the matrix A depends on the geometry of the dielectric boundary. Filling and inverting the matrix is a very time consuming process. In our simulations for ion channels [14][15][16] we used the fact that the dielectric boundary at the surface of the protein does not change during the simulation. Thus, the matrix can be precalculated at the beginning of the simulation and it does not really contribute to simulation time. The matrix-vector multiplication is also a time consuming step, but the simulations are still feasible. In the case of an ER fluid, on the contrary, the particles are moving during the simulation, the geometry of the dielectric boundary is constantly changing, and the matrix should be filled and inverted in every simulation step. This is why simulation is technically impossible using the DS model.

The interaction energy between two dielectric spheres in an electric field
Let us consider two DSs at a distance r 12 from each other. The dielectric constant inside sphere 1 (S 1 ) is ǫ 1 and its radius is a 1 . Similarly, the dielectric constant inside sphere 2 (S 2 ) is ǫ 2 and its radius is a 2 . In our calculations, we use spheres of equal unit radii a 1 = a 2 = a = 1. The spheres are embedded in a dielectric ǫ w . A homogeneous external electric field E is exerted on the system with strength of unity E = 1, so the only free charges in the system are the electrode charges σ and −σ that raise this electric field. The potential of this field is Ψ q (r) = −E · r.

43002-3
The total energy of the system is W = W q +W h , where W q is the free chargefree charge interaction energy, while W h is the free chargeinduced charge interaction energy. The latter, by symmetry, can be computed as the energy of the induced charges in the potential field of the free charges: Since the external field is constant, E can be brought out from the integral and the energy becomes The integral can be divided into two integrals on the two spheres where p i is the dipole moment on the S i sphere. So the energy is where we defined the electric field as pointing in the direction of the z-axis E = E k and p i is the zcomponent of the dipole moment. Of course, the dipole moment depends on the mutual position r 12 = r 1 − r 2 (with respect to the electric field) of the particles: p i = p i (r 12 ).
The interparticle potential energy between the two spheres is defined as the difference between the total energies for a given mutual position and for the case when the spheres are infinitely far from each other:

The point dipole models
The interaction potential between two point dipoles is φ DD (r 12 , p 1 , p 2 ) = −3 (p 1 · r 12 )(p 2 · r 12 ) where p i is the point dipole moment in the center of sphere S i [given by equation (2.1), therefore, p i is also a result of polarization] and r 12 = |r 12 |. Since p 1 and p 2 are parallel to each other and to E, this potential can be written as follows: φ DD (r 12 , cos θ, p 1 , where θ is the angle between the vectors r 12 and p i . In the above potential, the dipole moments induced by the external fields on individual isolated dielectric spheres were used. When the spheres are close to each other, nevertheless, not only the external electric field, but also the electric field produced by the dipole of the other particle acts on this sphere. This effect can be taken into account with the polarizable point dipole model, where, in addition to the permanent dipole induced by the external field, an induced dipole appears Since the induced dipoles produce fields that, in turn, induce dipoles, equations (4.3) and (4.4) have to be solved iteratively [17]. After the induced dipoles are obtained, the interaction between p 1 and p ind 2 as well as the interaction between p 2 and p ind 1 are computed using equation (4.1). Then, we add this correction to φ DD (r 12 , p 1 , p 2 ) thus obtaining φ PD (r 12 , p 1 , p ind 1 , p 2 , p ind 2 ), where PD stands for "polarizable dipole".

Results and discussion
We consider two cases that we call parallel and antiparallel cases. In the parallel case, both spheres have dielectric constants smaller than that of the surrounding medium. Thus, the dipoles induced on the spheres point into the same direction. We use the value ǫ 1 = ǫ 2 = 2 in this study. The dielectric constant of the solvent is ǫ w = 4. In the antiparallel case, the dielectric constant in one sphere is larger than ǫ w , while the dielectric constant in the other sphere is smaller than ǫ w . Thus, the dipoles induced by E on the two spheres point into opposite directions. We use the values ǫ 1 = 6 and ǫ 2 = 2 in this study.  Figure 1 shows the interaction energies as a function of the distance between the two spheres (this distance will be denoted by r = r 12 henceforth) for various situations. Based on the dielectric constants of the spheres, we can consider the parallel (ǫ 1 = ǫ 2 = 2) and the antiparallel (ǫ 1 = 2 and ǫ 2 = 6) cases. Based on the mutual position of the spheres, we can consider the aligned (E ∥ r 12 ) and the alongside (E ⊥ r 12 )

43002-5
positions. As seen in figure 1, the aligned position is the low-energy position for the parallel case (the classical head-to-tail situation), while the alongside position is the low-energy position for the antiparallel case (negative, attractive interaction energies). The interactions are repulsive for the other cases (aligned antiparallel and alongside parallel). All these energies are larger in magnitude when the particles are closer to each other. The agreement between the DD and DS potentials is surprisingly good (dashed lines vs. symbols in figure 1). Some deviation occurs for small interparticle distances because the DD approximation is not satisfactory when the two charge distributions are close to each other. The mutual polarization of the spheres can be taken into account by the PD model (solid lines in figure 1). The agreement with the DS results is excellent.
These results imply that the DD potential is an appropriate model of ER fluids in computer simulations. The good agreement is, nonetheless, surprising and counter-intuitive. It is not obvious from the first glance that the potential in equation (3.5) agrees with the potential in equation (4.1). To shed light on this, let us consider that the dipole moment induced on a sphere S i is (5.1) where h 0,i (r) is the distribution on the isolated sphere and ∆h i (r) is the distortion resulting from the effect of the other particle. The energy is obtained by multiplying by − 1 2 E . The term containing h 0,i (r) is canceled by W h (∞), so the interaction energy is simply The distortion of the induced charge distribution on sphere S 1 , for example, is proportional to the induced charge on this sphere: more induced charge can be distorted more (∆h 1 ∝ h 1 ∝ p 1 ). It is also proportional to the induced charge on the other sphere, because the electric field of sphere 2 that polarizes sphere 1 linearly depends on p 2 . This electric field depends on the cube of the distance inversely (∝ r −3 ). Finally, it will depend on the mutual position of the spheres which means an angle-dependent factor. In summary, which is exactly the form of the dipole-dipole potential in equation (4.2) for the special case of parallel dipoles. The mutual polarization of the two spheres is illustrated in figure 2 for the aligned parallel case. The induced charges on the surfaces of the spheres are plotted as functions of the angle with the electric field for various interparticle distances. For large distances, the charge distribution is symmetrical and the dipole moments on both of the spheres are the values for isolated spheres α i = (2 − 4)/(2 + 2 × 4) = −0.2 (for E = 1). The zero value of θ corresponds to the direction of the field, so the induced charge is negative on the hemisphere in the direction of the field and positive on the opposite hemisphere. This corresponds to a dipole moment whose direction is opposite to that of the field in agreement with the negative value of the polarizability.
When the distance between the spheres decreases, the symmetrical charge distribution is distorted. The charge distribution of sphere S 2 pulls some extra negative charge on the tip of sphere S 1 that is in the closest proximity to sphere S 2 . This extra negative charge is taken from all over the surface of the sphere, so the opposite positive charge appears there but on a larger surface, and this deviation in surface density is hardly distinguishable on the scale of figure 2.
The distortion of the h i (r) surface charge is exactly the ∆h i (r) surface charge introduced in equation (5.1) and shown in figure 3 for the case considered in figure 2. The ∆h i (θ) function is multiplied by sin θ, so the total charge on a ring at angle θ is plotted. It is seen now that the integral of this charge distribution is zero.
The interaction energy of this charge with the external field [see equation (5.2)] produces almost the same energy as that computed from the DD interaction potential [see equation The negative value means that the dipole is directed opposite to the electric field. Therefore, the total energy in this case is 0.2. When the spheres are close to each other, the charge distributions are slightly distorted on them due to polarization. This distortion is larger where the two spheres are in the closest proximity to each other. This corresponds to a slightly different dipole moment and a slightly different energy. The difference between this energy and 0.2 gives the φ DS (r 12 ) interaction energy. This energy is very close to the φ DD interaction energy between the dipoles, which is quite surprising given the fact that the φ DS energy is the result of a distortion of the induced charges and thus, is the result of a change of the dipoles, while φ DD is computed from the dipole moments fixed at their values at infinite separation.
In the next step, we study the effect of the mutual angular position of the spheres for a given distance: we fix r /a = 2.5 and change the angle between r 12 and E from 0 to π/2. Figure 4 shows the results as obtained from various models. Similar conclusions can be drawn as in the case of figure 1 for the 43002-7 comparison of the various methods. For the parallel case, the angle θ = 0 corresponds to the minimum energy head-to-tail position. The angle θ = π/2 corresponds to the maximum energy position where the dipoles of the spheres are next to each other pointing to the same direction. This means that the chains will repulse each other if the particles are in the same planes. A shifted position of chains corresponds to a more stable configuration at high densities when the chains are forced to be close to each other.  In the antiparallel case, similar conclusions can be drawn except that the configurations for the minimum and the maximum energy are interchanged. Here the minimum energy position is when the two spheres are next to each other with dipoles directed in opposite directions (blue symbols and curves in figure 1 and θ = π/2 in figure 4). The maximum energy position is when the two dipoles are aligned on the same line, a "head-to-head" position (green symbols and curves in figure 1 and θ = 0 in figure 4). The induced charge is shown in figure 5. The profiles are different for the two spheres because the polarizabilities of the two spheres are different now; different both in sign and magnitude. The profiles for the right hand sphere (ǫ 2 = 2) are similar to those in figure 2, while the profiles for the left hand sphere (ǫ 1 = 6) have decreasing tendency as a function of θ. The polarizability of this sphere is positive and smaller in magnitude than the polarizability of the other sphere: |α 1 | = 0.143, while |α 2 | = 0.2.

Summary
We have presented calculations for the interaction potential between two electrorheological particles, which, in a more detailed description, can be modelled as DSs immersed in a continuum dielectric that has a dielectric constant different from that of the sphere. We have shown that the interaction energy originated from the distortion of the induced charge distribution on the surface of the sphere as an effect of the presence of the other sphere is well reproduced by point dipoles placed in the centers of the spheres. Surprisingly, even the DD model (where this dipole is fixed at the value of the isolated sphere) gives a reasonable description.
We conclude that the DD or the PD models are useful simplified representations of the DS model for application in computer simulations. The potential acting between ER particles is used to calculate the energies in Monte Carlo simulations. Forces used in molecular dynamics or Brownian dynamics simulations can be straightforwardly derived from the potentials.