{'fraction_non_alphanumeric': 0.03048806524616764, 'fraction_numerical': 0.010162688415389214, 'mean_word_length': 4.367866147146459, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Ground state of a polydisperse electrorheological solid: Beyond the dipole approximation

The ground state of an electrorheological (ER) fluid has been studied based on our recently proposed dipole-induced dipole (DID) model. We obtained an analytic expression of the interaction between chains of particles which are of the same or different dielectric constants. The effects of dielectric constants on the structure formation in monodisperse and polydisperse electrorheological fluids are studied in a wide range of dielectric contrasts between the particles and the base fluid. Our results showed that the established body-centered tetragonal ground state in monodisperse ER fluids may become unstable due to a polydispersity in the particle dielectric constants. While our results agree with that of the fully multipole theory, the DID model is much simpler, which offers a basis for computer simulations in polydisperse ER fluids.


I. INTRODUCTION
The study of structure formation of electrorheological (ER) fluids has attracted increasing interest in recent years for its fundamental and technological values. Upon the application of an external electric field, the suspended particles in an ER fluid aggregate into chains and then columns parallel to the field, and drastically change the rheology of the suspension [1]. Tao and his coworker first suggested the existence of microcrystalline structures inside the columns and identified its ground state to be a body-centered tetragonal (bct) lattice [2]. These authors developed an analytic method based on a point-dipole (PD) approximation, i.e., treating the dielectric spheres as point dipoles interacting with one another. The possible ground state is the configuration that minimizes the dipole interaction energy and consequently the total Coulomb energy. The idea was soon confirmed by computer simulations and experiments [3,4]. The discovery of this property is not only helpful to understand further the ER response to the external field, but also offers a new technique to form mesocrystals with unique photonic properties [5,6].
After many efforts to reveal the details of structure formation in an ER fluid, it has been known that what structure will be formed in such a system actually depends on many factors, such as volume fractions [7], the size distribution of particles [8,9], the field frequency [10], etc. Combining the electrorheological and magnetorheological effects [5,11] or applying a rotating electric field [12] was also found to cause a structure transition from the bct to face-centered cubic (fcc) lattices [13].
The influence of the dielectric constant is another interesting topic since dielectric mismatch between the particles and the fluid is widely accepted as the main reason for the ER phenomenon. Davis found that the bct, fcc and hexagonal closed packed (hcp) lattices degenerate when the particle permittivity (ǫ p ) is much larger than that of the base fluid (ǫ f ) [14]. Clercx and Bossis [15] developed a fully multipolar treatment and compared the bct and fcc lattices for some values of ǫ p /ǫ f . Their results are in agreement with Davis. Lukkarinen and Kaski studied both the free energies of many types of lattices and dynamical effects in an ER fluid containing particles with the dielectric constant being greater and/or smaller than that of the fluid [16,17]. However, all of the previous work has either limited to the extreme case (ǫ p ≫ ǫ f ) or been too complicated to be adopted in further studies. And the understanding of the way dielectric constant affects structure formation is far from being sufficient.
Apart from theoretical work, computer simulation is another effective method to gain insight into ER effects. The point-dipole approximation is routinely adopted in simulation [18] because of its simplicity. Since many-body and multipolar interactions are neglected in this approximation, the predicted strength of ER effects is of an order lower than the experimental results . On the other hand, the accurate theoretical models are usually too complex to use in dynamic simulation of ER fluids. Hence, a model which is easy to use but beyond the point-dipole results is needed in computer simulations.
In this paper, we use a dipole-induced-dipole (DID) model [19] proposed by one of the authors to study the ground state of ER fluids in a wide range of the permittivity ratios: 0 < ǫ p /ǫ f < ∞. This model accounts for the multipolar interaction but is significantly simpler compared with existing multipolar theories [15,20,21], therefore it can serve as a candidate for computer simulation instead of the traditional point-dipole approximation. Some computation has been carried out such as the calculation of the interacting force between particles of different sizes and various dielectric constants and the simulation of the athermal aggregation of particles in ER fluids, both in uniaxial and rotating fields [19,22,23]. The purpose of this paper is to use this model to deal with structure formation in ER fluids and to offer some theoretical prediction as instructions to further dynamic simulations of ER fluids. It is found the ground state of ER fluids may vary with the dielectric constant, and the critical dielectric contrast is estimated. Furthermore, we obtain the interaction between two chains containing particles of various dielectric constants and apply it to polydisperse ER fluids where the particles have the same size but different permittivities. The effects of polydispersity on structure formation are investigated.

II. INTERACTION BETWEEN TWO PARTICLES
We start with briefly reviewing the DID model and then applying it to deal with the interaction energy between two particles. The DID model is defined from a multiple image method [19] which is based on a generalization of the image method to dielectric spheres. First consider a simple situation in which a point dipole p is placed at a distance r from the center of a sphere. The orientation of the dipole is perpendicular to the line joining the dipole and the center of the sphere. If the sphere is conducting, the induced image dipole is exactly given by p ′ = −p(a/r) 3 and at a distance r ′ = a 2 /r from the center. Generalizing this result to a sphere of dielectric constant ǫ p placed in a host medium ǫ f reads p ′ = −βp(a/r) 3 with β as the dipolar factor β = ǫp−ǫ f ǫp+2ǫ f . If the orientation of the point dipole is parallel to the axis, then p ′ = 2βp(a/r) 3 . In the limit β → 1, the above results reduce to the conducting sphere case.
Then consider a pair of dielectric spheres A and B, of the same radius a but different dielectric constants ǫ p1 and ǫ p2 , separated by a distance r, in a base fluid of a dielectric constant ǫ f . Upon the application of an electric field E 0 , the induced dipole moments in the individual spheres are respectively given by with the dipole factors β = ǫ p2 +2ǫ f . The initial dipole moment p a0 induces an image dipole p b1 in sphere B, while p b1 induces another image dipole p a2 in sphere A. As a result an infinite series of dipoles p a0 , p a1 , p a2 , ... are formed inside sphere A and the total dipole moments for transverse and longitudinal fields are respectively given by The parameter α in Eq. (1) and (2) satisfies Similar expressions for the total dipole moments inside sphere B, namely p bL and p bT , can be obtained by interchanging β and β ′ . We should remark that the present generalization is only approximate because there is no image method for a dielectric sphere. An integral equation approach was proposed to examine the validity of the multiple image method in Ref. [19]. It was shown that Eq. (2) and (3) can produce good results at high dielectric contrast β → 1 as expected, and the model keeps reasonable even in the low contrast case [19,22].
How many terms are retained in the analytic multiple image result determines different models. With only the n = 1 term in the series of Eq. (2) and (3), it reduces to the traditional point-dipole approximation, while the DID model is defined by retaining the first two terms (n = 1 and n = 2). Taking part of the multipole effects into account, the DID model is generally better than the PD results, especially in polydisperse systems, and can be chosen as an effective approximation to calculate multipolar interaction [19,22,23].
The incorporation of multipole effects leads to one of the obvious differences between the PD and DID model: the symmetry of negative and positive βs will be broken in the latter, which is illustrated in Fig. 1. This asymmetry of the dipolar factor will also be exhibited in the interaction energy between particles, as shown below, and finally affect the structure formation in ER fluids. It implies that different structures may form when the dielectric constant of particles is smaller than that of the host fluid.
Now we begin to apply the DID model to deal with the interaction energy between two particles. The electrostatic energy of the two particles upon the application of an electric field E 0 at an arbitrary angle θ (Fig. 2a) is given by [24] Substituting the dipole moment expressions into Eq. (5) and letting λ = β ′ /β, p 0 = p a0 , we can get the energy using the DID approximation as where the dimensionless parameters have been defined: The first term W f in Eq. (6) is the energy of the individual particles by the applied field, which is independent of the relative position of the particles, while the remaining terms correspond to the interaction energy which determines the formation of ground states [2]. Note that the term W 0 is, if λ = 1, exactly the result of Tao et al.'s work [2] , in which they dealt with the monodisperse case in the use of a PD approximation and concluded that a bct lattice will form as the ground state. So what we are interested in is the last two terms which describe the DID correction to the total interparticle energy. Adding the hard-core repulsion which keeps the particles from coalescing (r ′ > 2), W 1 and W 2 can be rewritten as Here we have taken being the coordinates of the particles and z axis is chosen as along the field direction). Obviously, the DID correction cannot be ignored when β is large, while the term W 1 being proportional to β instead of |β| shows the asymmetry of β in the DID model.

III. MONODISPERSE SYSTEMS
Now we study a monodisperse ER fluid composed of particles with a dielectric constant ǫ p , which can be larger (β > 0) or smaller (β < 0) than that of the fluid ǫ f . It has been mentioned that the interacting particles in a ER fluid will first form chains between electrode plates, then the chains aggregate into columns containing microstructures. Combined with the images, the chains can be treated as infinite if the distance between two plates is large enough [2]. The interaction energy per particle is divided into two parts: one is from the self-energy of an infinite chain, i.e., the interaction energy between the particles belonging to the same chain; the other from interaction between different chains [2].
A. u s , the self-energy Consider an infinite chain containing particles of radius a and dielectric constant ǫ p at r j = 2ajẑ (j = 0, ±1, ±2, ...) (Fig. 2b). The self-energy of the chain is given by W i refers to the interaction energy between two particles. Substituting the DID terms in the particle interaction, W 1 and W 2 , into Eq. (10), we get where ζ(n) is the Zeta function defined as ζ(n) = ∞ n=1 1 j n . Again, the DID correction is comparable to the PD result [2], −0.300514, when β is large. The total self-energy u s in a DID model is the sum of u s 1 , u s 2 and that from PD assumption which represents a correction to the established PD results.

B. u i , the interaction energy
The interchain energy between two parallel infinite chains, with vertical shift z and separated by distance ρ, is given by 1 2 u i (ρ, z), where u i (ρ, z) is the interaction between one dielectric particle at r p = ρ + zẑ near the infinite chain in which particles locate at r j = 2ajẑ (j = 0, ±1, ±2, ...) (Fig. 2c), Hence we have the DID correction of u i and Following a Fourier series technique proposed by Tao et al. [2] we expand u i 1 and u i 2 into with ω ≡ sπ/a.
Γ(x) and K i (x) in the above equations are the gamma function and i-th order modified Bessel function, respectively. The sums above can be easily evaluated numerically, and the expression of u i (ρ, z) is given by C. Possible ground states The total interaction energy per particle in a certain configuration is [11] where ′ k denotes the summation over all chains labeled k except the one containing the considered particle. Since the self-energy of a chain is independent of the structure, what really affects the energy difference between various lattices is u i (ρ, z). Fig. 3 shows the dependence of u i on the shift z for different values of ρ when β > 0. The results from the PD approximation are also plotted for comparison. The interaction may be either attractive or repulsive depending on the shift z, and the range of z in which two chains attract each other enlarges quickly when β increases. This implies a tendency to form a more close packed structure than the bct lattice with increasing dielectric contract. The most possible candidate is the fcc lattice, which is nearest to the bct one in energy [2]. An estimation including the nearest and the next-nearest neighboring chains can give the energy gap between the two lattices which decreases when β increases from 0 to 1. For the limit case β = 1 (ǫ p ≫ ǫ f ), ∆u ≃ 6.6× 10 −4 and the two phases can be regarded as degenerate, in agreement with the conclusion of previous work [14,15]. Reversed situation happens when β < 0 (Fig. 4 ). Now it's the repulsive region that enlarges and the system tends to form a looser structure as β approaches the negative limit − 1 2 (ǫ p ≪ ǫ f ). The loosest structure is certainly the one composed of separate chains. Calculating the difference between the bct lattice and separate chains, we get When β = −0.323 a transition from the bct lattice to separate chains happens and the latter becomes more and more stable while the absolute value of β increases. Finally we want to make a comparison with Clercx et al.'s elegant work based on a multipole-expansion theory [15]. They expanded the potentials inside and outside the spheres in terms of solid spherical harmonics and gained the values of multipole moments Q i lm by solving a set of linear equations determined by the boundary conditions at the grain surfaces. The upper limit of l, denoted by L, determines how many multipole effects are considered: L = 1 refers to the simple dipole approximation and L = ∞ to the exact calculation where all multipole effects are included. This multipole expansion theory and our DID model are based on different pictures and there exists no direct and exact equivalent relationship between them. The calculation of the interparticle energy W i using Eq. (6) and Clercx et al.'s model shows that the DID results reflects, although not exactly, some characteristics of the third order(L = 3) multipole-expansion theory, i.e., an octupole effect. This is particularly obvious when θ = π/2 even in the touching-particle case (See Fig. 5), where the point-dipole approximation is known to err considerably. We expect that the DID results may be better than the first approximation of Clercx et al.'s model since the former takes higher multipole effects into consideration. Unfortunately the available data related with the energy ground state are limited in Ref. [15] and we cannot give a thorough comparison. However, from Table 1 in which the ratios of the total induced dipole per particle between the bcc and the fcc lattice calculated with our model are compared with the available results of Clercx et al., we can still see that the DID results is closer to the exact solutions than the dipole approximation in a wide range of dielectric mismatch.

IV. POLYDISPERSE SYSTEMS
So far we have discussed the interaction energy and ground states in an ER fluid containing the same particles. Now we begin to investigate a more "natural" polydisperse system in which the particle dielectric constant ǫ p , and consequently the dipolar factor β, has a probability distribution. Let the average interaction energy between a pair of particles in this system denoted by W * i . Assuming the dipolar factors β and β ′ of the two particles are independent random variables which have the same distribution we can get W * i by taking an average of Eq.(6) over β and β ′ as where W i is the interaction energy, including the PD results and two DID terms, between identical particles with the average dipolar factor β 0 , while W δ is the polydispersity correction with σ and β 0 as the standard deviation and mean of the distribution of β. On the basis of the expression of interparticle energy, we can easily obtain the average self-energy of an infinite chain u s * and the interaction energy between two chains u i * , in a polydisperse system where particle permittivities are distributed randomly, as where the energy increment caused by polydispersity u s δ and u i δ are given by In the above equations, δ is defined as σ/β 0 to describe the degree of polydispersity, u s , u i are the total self-energy and interchain energy for a monodisperse system in which all particles have the same dipolar factor β 0 , while u s 1 , u s 2 and u i 1 , u i 2 are the corresponding DID corrections. They can be calculated using Eq. (10) - (22).
As a consequence of the change of the chain energies caused by polydispersity, the energy gaps between different lattices also vary and hence affect structure formation in ER fluids. Fig. 6 shows both ∆u − β 0 and ∆u δ − β 0 with δ = 0.1 to exhibit the effect of polydispersity on ground states. ∆u is defined as u bct − u f cc when β 0 > 0 and u bct − u sc when β 0 < 0 (Here "sc" means separated chains) for a system composed of particles β 0 , while ∆u δ is that caused by polydispersity. It is shown that ∆u δ keeps opposite to ∆u when β 0 is larger than −0.323 and not very near 1.0. In this range the bct lattice has been proved to be the stable ground state in a monodisperse system, and the existence of polydispersity tends to destroy its stableness. However, the magnitude of the ratio ∆u δ /∆u is only of the order 10 −2 , so we can conclude that, in a polydisperse system where the dielectric mismatch between particles and the fluid is not very large, similar structures will form but the ground state is not so stable as that in a monodisperse system.
All the results discussed above are based on the assumption that the dielectric constants of particles are distributed at random in the configuration. But the question is: Is the distribution really random, especially when the dielectric constants of particles differ greatly from each other? In the extreme case, will the phenomenon "phase separation" happen, i.e., identical particles gather together upon the application of the electric field? In order to investigate this problem, we consider a two-component fluid containing particles with β 1 and β 2 . The fractions of particles are supposed to be p 1 and p 2 , respectively. First assume a totally random configuration will form in this system. The average self-energy u s mix and interchain energy u i mix in such a system can be calculated by means of the method described according to the binary distribution Then consider another configuration in which only identical particles aggregate into "uniform" chains. Two types of chains are formed in this system: one, called chain A, is composed of particles β 1 , and another, chain B, of β 2 . The average self-energy u s unif is determined by and the interaction energy between two uniform chains satisfies with λ is the ratio of βs of the two chains and other parameters are defined as previously.
We have compared the energies in these two configurations and plotted the results in Fig. 7. The upper graph in Fig. 7 shows the dependence of ∆u s ≡ u s mix − u s unif on β 1 and β 2 , and the lower one shows ∆u i ≡ u i mix − u i unif . Note that here u i unif is the interaction between chain A and chain B. The results when we select other chain configuration such as chains AA or BB are similar. The fraction of the particles with β 1 is 0.5, and the relative coordinates of the interacting chains are chosen as ρ = √ 3a, z = a. It is shown that u s mis is always larger than u s unif , while ∆u i is usually negative for most values of β 1 and β 2 but much smaller than ∆u s in magnitude. Since the self-energy is usually dominant in the total interaction energy per particle, chains of identical particles may be more stable than those containing different types of particles, particularly when β 1 and β 2 approach the positive limit 1.0 and negative limit -0.5 simultaneously.
The interaction form in polydisperse fluids is more complex than what we have considered because it is sensitive to the micro structure formed in the ER solid and the number of possible configurations in a polydisperse system is much larger than that in a monodisperse one. The energy gap between the random and phase separation configurations may be narrowed because of the contribution of interaction between chains. But we can still expect that the ground state of such a two-component system may contain quite a few of these "uniform" chains instead of a totally random configuration. And the simplicity of our DID model will also make it possible to carry out computer simulations in polydisperse systems and study the dielectric effects in a more detailed way [23].

V. CONCLUSIONS
In conclusion, we have presented a DID model to study the dielectric effects on structure formation in ER fluids. Based on the DID expression of the interacting energy of two particles, we have corrected the PD results of the self-energy and the interaction energy of chains, including multipole effects partially. Series expressions of these energies are obtained and used to calculate the effects of dielectric mismatch between particles and the fluid. Both monodisperse and polydisperse cases are discussed and some interesting results are obtained as follows.
(1) In monodisperse systems where particles are negatively polarized, there may exist a phase transition form the bct lattice to the configuration of separate chains when β < −0.323; (2) Polydispersity of particle dielectric constants in ER fluids will cause the bct ground state not so stable as that in monodisperse systems when the dielectric mismatch between the particle and the fluid is not very large; (3) When the dielectric constants of particles differ much with each other, identical particles tend to aggregate into uniform chains.
The ER effects in polydisperse system will be studied more deeply in our future simulation work based on the DID model.    TABLES   TABLE I. The comparison between the DID model and the full multipole theory proposed by Clercx and Bossis . α is defined as ǫ p /ǫ f , where ǫ p and ǫ f are the dielectric constants of the particles and the fluid, respectively. ∆ DID is the polarization ratio between the bct and fcc lattices using the present model. ∆ 1 and ∆ ∞ are from the first-order approximation and the exact results of the full multipole theory, respectively.