Ising Model for Conductive Percolation in Electrically Conductive Adhesives by Considering Random Arrangement of Conducting Particles

Problem statement: Thermodynamically, particles in composites will ar range in a way such that the Helmholtz free energy is minimized. H owever, even a single structure has the lowest free energy, it should not ignore the probability of oth er structures having larger energies to occur, alth oug at small chances. Approach: All possible arrangements of particles in the comp sites, therefore, must be taken into account in the theory or simulation d evelopment. Results: The composite energy depends on the interaction between components in th e composites. To consider the effect of interactions on energy, in this study we used a sim ple Ising model incorporated with the BraggWilliams approximation. We used the model to predic t the average packing fraction and the percolation threshold in composites as well as othe r quantities related to percolation phenomenon. Conclusion/Recommendations: We found several predictions that have not been re port d by previous authors. This model can be important in the underst anding conductivity development in electrically conductive adhesive composites.


INTRODUCTION
The properties of composites of particles dispersed in continuum media have been reported by many authors several decades, covering numerous composites, such as conductive particles dispersed in insulating adhesives (electrically conductive adhesives) for microelectronic applications (Mikrajuddin et al., 2000;Yim et al., 2008;Lin and Chen, 2008;Li et al., 2008;Lin and Chiu, 2008;Morris and Lee, 2008;Zenner et al., 2008;Inoue et al., 2008;2009;Inoue and Suganuma, 2007;2009;Kim and Paik, 2008;Mundlein et al., 2002;Li et al., 1993;Tongxiang et al., 2008;Lee et al., 2005;Sander et al., 2002;, insulating particles dispersed in Ionic adhesive for use in solid batteries or fuel cells (Mikrajuddin et al., 2000;1999), colloidal systems (Lu et al., 2006;, pharmaceutical tablets (Stromme et al., 2003) and other composites such as fibers in concrete (Newman, 2002) and fly ash-concrete (Mohan et al., 2012). One interesting phenomenon exhibited by these composites is the occurrence of percolation threshold, a quantity which divides two strongly different states, such as conductive and insulating states (Mikrajuddin et al., 2000;1999), magnetic and non-magnetic states (Thorpe, 1978), crystal and amorphous states, Surprisingly, in recent progress, the percolation phenomena are also applied to other fields, in the past of which likely did not show any relation with material composites such as communication systems such as the internet and other Peer-To-Peer networks (Pastor-Satorras and Vespignani, 2007), dynamics of epidemic spreading (Anderson and May, 1992;Meyers, 2007), HIV infection to AIDS (Kamp and Bornholdt, 2002), social networks (Chen et al., 2007a;2007b).
Most theories and simulations developed so far confirmed the coordination number, γ, (number of particle nearest neighbors) affects the percolation threshold, vc . Generally, the percolation threshold decreases when the coordination number increases. Some authors reported that the percolation threshold depends on the coordination number as v c = D/γ(D-1) with D is the dimensionality (2 for two dimensions and 3 for three dimensions) (Chelidze, 1982). Using the modified effective medium approximation, we previously reported that the percolation threshold occurs at c v f / 2 / = γ , with f is the packing fraction (Mikrajuddin et al., 1999) and by adopting a theory for sol-gel development in polymerization for describing the connectedness of particles in composites, we have predicted the occurrence of the percolation threshold at vc = f/ (γ-1) (Mikrajuddin et al., 2001). In most theories and simulations, the arrangement of particles in the composite was taken as an initial axiom. For example, some authors treated the particles have been arranged in a simple cubic structure and other authors assumed the particles have been arranged in a body centered cubic structure, prior to theory or simulation development. Indeed, this assumption may raise some questions. For example, how can the particles in the composite arrange at a certain structure? This question has become more significant when applied to liquid-like (easy flow) composites, such as colloidal systems, where microscopically the particles may displace locally. Slight displacement in the particle position can create any chances for the occurrence of different arrangements at different times or different locations.
The assumption that the particles have been arranged in a specific three dimensional structure is equivalent to assuming that the particles have undergone a process like a self organization or externally induced organization. However, it is hardly to accept the particles in the composites have passed such organization history. Instead, the particles in the composite may develop arbitrary structures. This is likely one reason why a model to describe the inhomogeneous particle distribution in the composite and the effect of the eddy current have been proposed (Vinogradov et al., 2009). In addition, based on numerous simulations and observations reported so far, no clear evidence the particles were arranged in a certain structure when forming large clusters (Lu et al., 2006;Zaccarelli et al., 2008).
Furthermore, recent reports might strengthen questionability of the above axiom. For example a "non permanent bond" model in colloidal systems as simulated by Coniglio et al. (2004) and Candia et al. (2005) might suggest that the arrangement of particles is possibly not fixed. They therefore hypothesized that the bonds between colloidal particles have a lifetime, which is extremely high at low temperature and decreases as the temperature increases. The finite lifetime bonds were also introduced in the study of a lattice model for gelation phenomena (Zaccarelli et al., 2009).
Thermodynamically, the particles in the composite will arrange in a way such that the Helmholtz free energy is minimized. However, even a single structure has the lowest free energy, it must not ignore the possibility of other structures having larger energies to occur, although at small probabilities. All possible arrangements, therefore, must be taken into account in the theoretical formulation and the microscopically observed arrangement should be the statistical average of all the possible arrangements. Since each arrangement is attributed to a specific coordination number, this speculation leads to the existence of an average coordination number. Indeed, the concept of average coordination number is logically accepted when we consider some previous reports related to average bond connection between particles per nodes (Dhydkov, 2009;Blumenfeld et al., 2005) and the concept of non permanent bonds (Coniglio et al., 2004). Lagemaat et al. (2001) reported that in porous TiO 2 film, at a film porosity of 58%, about 10% of particles have two neighbors, 25% of the particles have four neighbors and less than 1% of the particles have eight neighbors and the average coordination number is 4.1. These reports clearly proved the presence of a non fixed coordination number in composites.
At the present work we will focus our consideration on the conductivity development in electrically conductive adhesive composites. The probability of the occurrence of a specific arrangement of particles depends on the configuration energy, which is dependent on the inter particle interaction, inter adhesive interaction and interaction between particle and adhesive. To consider the effect of interactions on the configuration energy, in this study we used a simple Ising model incorporated with the Bragg-Williams approximation. The developed model was then used to predict the average packing fraction and the percolation threshold in the electrically conductive adhesives as well as other quantities related to percolation phenomenon. With best our knowledge, no report has been published by any authors on this topic.

MATERIALS AND METHODS
We divide the composite into lattices. Each lattice point is occupied either by a particle or a matrix element. All possible lattice structures can occur in the composite, the probability of which depends on the configuration energy.
Assume the number of lattice points is N and the number of particles is N p . The number of lattice points occupied by matrix elements is a N m = N-N p . Suppose the interaction energy between particles is ε p , between matrix elements is ε m and between particle and matrix element is ε pm . In the composite, three kinds of contacts occur, i.e., particles-particle contact, matrix-matrix contact and particles-matrix contact. We assume the number of particle-particle contacts, matrix-matrix contacts and particle-matrix contacts as N pp , N mm and N pm , respectively. The energy for a certain configuration {N pp , N mm , N pm ,} can be written as: By a simple calculation (Huang, 1987), we can rewrite the configuration energy in Eq. 1 as: With, E 1 = ε p + ε m -2 εpm , E 2 = ε pm -ε m and E 3 = ε m /2.
The Bragg-Williams approximation assumes the relationship of N pp /(γN/2) ≈(N p /N) 2 to hold (Huang, 1987). By introducing a long-range order parameter, L, through the following relation Eq. 3: with, -1≤ L ≤ 1, Eq. 2 can be written as: or the energy per site is: Here L depends on γ. For the same volume fraction, different arrangements (different γ) lead to different L. We can then simply write the configuration energy in Eq. 4 as: The energy expressed in Eq. 5 will be used for estimating the configuration probability.
We assumed the particles are spherical in shape. For each arrangement of particles in the lattice, there is corresponding parameter of packing fraction, ƒ. If v is the volume fraction of particle we can show: This relation states that if all lattice points are occupied by particles (N p = N), the volume fraction of particles is exactly equal to the packing fraction (v = ƒ). Based on Eq. 3 and 7, the long range order parameter depends on the volume fraction as: For three dimension arrangements, we assumed that only diamond (γ = 4, ƒ d = 0.340), simple cubic (γ = 6, = ƒ SC 0.523), body-centered cubic (γ = 8, ƒ bcc = 0.680), facecentered cubic and hexagonal closed packing (γ= 12, ƒ ƒcc = = ƒ hcp 0.740) possibly occur.
We proposed here, if the volume fraction v is less than ƒ d all arrangements possibly occur. If ƒ d < v < ƒ sc , the diamond structure is absent since the particle content is higher than the maximum content allowed for diamond arrangement. If ƒ sc < v < ƒ bcc , the diamond and simple cubic arrangements are absent and if ƒ bcc < v < ƒ fcc = ƒ hcp , only the face-centered cubic and hexagonal close packing are allowed.
We now define the average coordination number of different volume fractions of particles as Eq. 8: where, Θ(x) is the Heaviside step function, satisfying Θ(x) = 1 if x>0 and Θ(x) = 0 if x≤0 ƒ γ is the packing fraction for arrangement with coordination number γ and β = -1/kT with k is the Boltzmann constant. For simplicity, we selected the energy reference such that ε m = 0. This selection leads to E 1 = ε p -2ε pm , E 2 = ε pm and E 3 = 0. With these assumptions, Eq. 6 can be rewritten as: with, ξ = ε pm /ε p and φ = Nε p /8kT. We used this model to estimate the effect of particle volume fraction on the conductivity of electrically conductive adhesive composites. Previously we have proposed a modified effective medium approximation to calculate such conductivity using the following Eq. 9 (Mikrajuddin et al., 1999;Abdullah et al., 2003) with P pp is the probability of particle-particle contact P pm , is the probability of particle matrix contact, P mm , is the probability of matrix-matrix contact, where P pp + P pm + P mm = 1, σ pp is the conductivity at contact point between particles σ pm , is the conductivity between contact point between particle and matrix σ mm , is the conductivity at contact point between matrix elements (the conductivity of adhesive) and σ e is the effective conductivity of the composite.

RESULTS
We did simulations for different values of ξ (fraction of particle-matrix interaction energy and particle-particle interaction energy) and different values of η (represents the strength of interaction energy). Figure 1 shows the effect of volume fraction of particles in the average coordination number when taking η = 2. The positive sign for η means that ε p is more positive than ε m , to state that replacement of matrix element with particles increases the configuration energy (the attraction force between particles is weaker than between matrix elements). The effect of η of the average coordination number is shown in Fig. 2. This figure shows the effect of filler volume fraction on the average coordination number when η = 0.1 for different ξ.
We also inspected the effect of negative values of η on the average coordination number. Figure 3 and 4 were obtained using η = -2 and η = -0.1, respectively. Figure 5 shows the conductivity as a function of volume fraction of particles for η = 2 and Fig. 6 for η = 0.1 at different ξ. Apparently, the conductivity percolations occur at volume fractions between 0.2 and 0.4. The exact location of the percolation threshold and how to find it will be explained later. The volume percolation decreases when ξ increases. For large η, the percolation thresholds for different ξ occurs at a wider range of volume fractions.

DISCUSSION
We identified in Fig. 1 a discontinuity in the average coordination number at some volume fractions. The existence of this kind of discontinuity has also been discussed by Zhang et al. (2008;2009). At high ξ, or ε pm >>ε p , the average coordination number quickly decreases to 4 (the coordination number of the diamond structure) when increasing the volume fraction of particles. A high ξ means two contacting particles are bound strongly. The binding energy between particle and matrix is more positive than the binding energy between particles, so that in order to minimize the composite free energy, the arrangement of particles in the composite must be in such a way that the contact area between particles and matrix is as small as possible (contribution of contact surface energy must be suppressed) and this is found when the arrangement is in small coordination number (which leads to small contact surface). The volume of particles per cell is ƒ and the contact area of a particle and matrix in a cell is s ∝ f 2/3 , to prove that reducing in the coordination number will reduce the packing fraction and therefore will reduce the contact area. As the interaction energy between particles is more positive than the interaction energy between particle and matrix (ξ< 1), the composite will increase the contact area between particles and matrix as much as possible to reduce the free energy. It is reached by increasing in the coordination number so that the number of particles in a lattice cell increases (leads to increasing the surface area). It results in the increase in the coordination number. As shown in Fig. 1, when ξ < 0.5, even the simple cubic arrangement does not occur. Only the body centered cubic (bcc) and face centered cubic (ƒcc) or hexagonal closest packing (hcp) arrangements can occur. The bcc arrangement occurs when v > 0.68 and the ƒcc and hcp occur only when v > 0.68. At low interaction energy between particles,η = 0.1, low coordination number is hardly to occur (Fig. 2). The diamond structure is only obtained for ξ> 5. For ξ < 0.75, even the simple cubic arrangement does not occur. The η small means that the interaction energy between particles approaches the interaction energy between matrix elements.
When the interaction energy between particles precisely equals to interaction energy between particles and matrix and between matrix elements we have ε p = ε pm = ε m = 0. If this condition occurs, we have Ω = 0 and the average coordination number becomes: The η small means that the interaction energy between particles approaches the interaction energy between matrix elements.
Different results were obtained when the energetic interaction between particles and between particle and matrix element are more negative that the interaction energy between matrix elements. In this situation, the attractive forces between particles and between particle and matrix element are stronger than between matrix elements. Many colloidal systems such as those investigated comprehensively by Lu et al. (2006; are a class of this composite. The dependence of the average coordination number on the volume fraction of η = -2.0 is shown in Fig. 3 and four η = -0.1 is shown in Fig. 4. Based on Fig. 3, only for very large we can obtain the low coordination number until volume fraction of 0.34 (packing fraction of diamond). For small ξ(< 0.75) we did not obtain structures with small coordination number, even a simple cubic structure.
The condition of η < 0 implies that dispersing of particles into the matrix reduces the configuration energy. The most stable structure will be achieved when the particle content is as much as possible and it is obtained at a high coordination number. The more drastic situation was obtained whenη = -0.1 (Fig. 4), where structures with low coordination numbers (diamond and simple cubic) never appear.
Based on data in Fig. 1-4 we can conclude that the average coordination number is resulted by competition between η and ξ. Large η (large interaction energy between particles) tends to produce the low coordination number and large ξ tends to produce large coordination number.
By approximating σ pm ≈ σ mm we can write Eq. 9 as: ( ) pp e pp pp pp e mm e mm e P 1 P ( / 2 1) 0 ( / 2 1) The solution for σ e in Eq. 10 can be obtained easily by simply solving the root of a quadratic equation, i.e.: Let us assume σ pp >> σ mm . This is the general condition satisfied by conducting particles dispersed in insulating matrix. If γP pp /2-1≠ 0 we can approximate: We then obtain the following approximation for the effective conductivities. If γP pp /2-1 < 0 we have: or σ e is the same order of magnitude with σ mm . On the contrary, when γP pp /2-1 > 0 we have: or σ e is the same order of magnitude with σ pp .
From this result we concluded that γP pp /2-1= 0 is the condition when the conductivity changes abruptly from that of the matrix to that of the particles, i.e., The condition of percolation. Thus, the percolation threshold is the volume fraction when γP pp /2-1= 0 or P pp = 2/γ is satisfied. Using Eq. 3 and the Bragg-Williams approximation, the condition for percolation threshold is: Since the arrangement of particles in the composite may have differed γ and different L, here the percolation threshold is assumed to occur when the averages of γ and L satisfy Eq. 12: where, 〈γ〉 is given by Eq. 8 and 〈L〉 satisfies

∑ ∑
Remember that both 〈γ〉 depend of volume fraction of particles. To determine the volume fraction when Eq. 15 holds, we plotted 2/ 〈γ〉 and [(1+〈L〉)/2] 2 as a function of v in the same curve as shown in Fig. 7. These two curves intersect at the point of percolation. Figure 8 shows the percolation thresholds as a function of ξ at different η. For the three simulation parameters, the percolation threshold occurred in the volume fraction range between 0.24-0.34.
It should be noted that no single percolation threshold has been reported by the authors, either those obtained numerically or experimentally. Different authors reported different percolation threshold. For example, based on simulation up to 10 5 particles in a composite, Rintoul and Torquato (1997) shown a percolation threshold at φ c = 0.2895±0.0005. Sancaktar and Liu (2003) showed a percolation threshold of 30% in a composite of emeraldine salt polyethylene powder dispersed in an adhesive polymer.
Wen and Chung observed a percolation threshold occurred at sand fractions between 24 and 30% vol. In a mixture of sand and cement paste (Wen and Chung, 2007). Stromme et al. (2003) reported the percolation volume of 0.238 and 0.249 for MMC equilibrated at 45 and 75% RH, respectively. Simulation using lattice models on a recursive square Husimi lattice, Corsi and Gujrati (2006) found a percolation threshold between 0.245-0.25.  showed the occurrence of the percolation threshold at 22% in a composite of synthetic graphite particles dispersed in an epoxy-polyurethane resin. Very low percolation thresholds are usually observed in elongated filler such as carbon nanotube. For example, the percolation thresholds between 0.017-0.019 have been observed in polystyrene-MWCNT nanocomposites (Kota et al., 2007). For all η, the percolation threshold decreases when ξ increases. It mentions that the interaction energy between particle and matrix element influences the percolation thresholds. The same conclusions have been reported by Miyasaka et al. (1982); Sumita et al. (1986) and Corsi and Gujrati (2006). The increase in ξ means that the interaction energy between particles and matrix element is more positive that between particles. In order to achieve lower Helmholtz free energy, the particles are more favorable to create contacts between particles rather than between particles and matrix element. Therefore, connections between particles take place more easily to result a lower in the percolation threshold.
It is interesting to see that at low ξ, the percolation threshold increases with η and at high ξ the percolation threshold decreases with η. All curves belong to different η seem to meet at a single point of v c = 0.267 and ξ = 0.606. This point states that, as long as the fraction of energy interactions between particles and between particle and matrix element is constant at ξ = 0.606, the percolation threshold is always v c = 0.267, irrespective to the strength of interaction between particles and between particle and matrix element.
At percolation threshold, Eq. 11 happens. Based on Eq. 3 and 7 we can write: c v 1 L 2 f + = Therefore, we obtain the average packing fraction in the composite as: c 2 v f 1 L = + Figure 9 is a plot of the average packing fraction as a function ξ of at different values of η. The average packing fraction decreases with increasing of ξ (when particle-particle interaction is more positive than matrix elements' interaction). The average packing fraction is also dependent on the interaction energy between particles, where the average packing fraction decreases with the interaction energy when ξ > 0.45 and decreases with the interaction energy when ξ < 0.45. Figure 9 mentions that the packing fraction is not a fixed parameter even in a specific composite. Its value may change when other parameter in the composite changes. The non constant in packing fraction was also reported by Kim et al. (2007) in a colloidal gel system. Based on information presented in Fig. 8 and 9 we can conclude that the average coordination number per particles is an increasing function of the packing fraction. The same conclusion came from Buhot (1999).

CONCLUSION
The Ising model combined with Bragg-Williams approximation has been used to estimate the percolation phenomena in electrically conductive adhesives. The main point proposed here was all packing arrangements of particles are possible to occur with probabilities that depend on configuration energies. We found the average coordination number and the average packing fraction depend on the volume fraction of filler. The average coordination number was also dependent on the strength of interaction between particle-particles, between particlematrix and between matrix-matrix. We also identified the discontinuity in average coordination number at several volume fractions.