We report Monte Carlo simulations of systems of polydisperse prolate and oblate ellipsoids using the critical path based tunneling-percolation model. For polydisperse prolate ellipsoids, the critical percolation volume fraction $\varphi c$ is shown to have a quasi-universal dependence on weight-averaged aspect ratio. For polydisperse oblate ellipsoids, $\varphi c$ is shown to have a quasi-universal dependence on the apparent aspect ratio, which is a function of up to fourth moment of the size distribution, as given by the generalized connectedness percolation theory. The functions are observed to approach the theoretical predictions for higher volume fractions and higher aspect ratios. The model predictions are compared with experimental data available on polydisperse multi-walled nanotubes (prolate ellipsoids) and graphene nanoplatelets (oblate ellipsoids) to estimate the tunneling lengthscale which is found to be well within the expected range.

## I. INTRODUCTION

Nanofillers with highly anisotropic shapes, such as carbon nanotubes (CNTs), graphene nanoplatelets (GNPs), carbon black (CB), and metallic nanowires,^{1} are used as inclusions in polymer matrix materials to generate nanocomposites with superior electrical, mechanical, and thermal properties. In this study, the focus is on the enhanced effective electrical conductivity. Electrically conductive polymer nanocomposites can be used in static dissipative, semi-conductive, conductive, or electromagnetic shielding applications. The increased conductivity is also desired in preventing the aircraft structures made of carbon based polymer composites from localized damage due to lightning strikes.^{2} The strong dependence of the electrical conductivity on interfiller separation has also been exploited by using nanocomposites for strain measurement applications.^{3–5} Given the wide range of applications, electrical conductivity of nanocomposites has been an intense topic of research over the past few decades as an in-depth understanding of the various factors governing the electrical behavior is a crucial step in the advancement of nanocomposites based technology.

Electrical conductivity of insulating polymers can be increased by several orders of magnitude by adding a small amount of highly conductive nanofillers. The sharp rise in the electrical conductivity of polymer nanocomposites results from formation of a percolating network of electrically interconnected nanofillers through electron tunneling.^{6,7} For a given volume fraction of nanofillers, the polymer matrix layer separating two fillers is sufficiently thin so as to allow electrons to tunnel across the barrier. The tunneling conductance between two nanoparticles can be expressed as^{8}

where *δ _{ij}* is the shortest distance between two particles

*i*and

*j*,

*ξ*is the tunneling lengthscale, and

*g*

_{0}is a constant proportional to the contact resistance.

Given the extreme contrast between the electrical conductivity of the fillers and the polymer matrix, the effective electrical conductivity of the nanocomposite can be modeled as a tunneling-percolation problem with conductive fillers distributed in a continuous non-conductive matrix. Modeling this behavior is a challenging task given the large number of variables (i.e., “all thinkable parameters” as concluded in Ref. 9) that play a key role in determining the percolation behavior of the system.

It was shown in Ref. 8 that effective electrical conductivity of nanocomposites can be obtained from the global tunneling network (GTN) of the interparticle tunneling conductors as defined in Eq. (1). Furthermore, as the tunneling conductance varies exponentially with the interparticle distance, the effective conductivity of the GTN can also be obtained using the critical path approximation.^{10} The major advantage of the critical path method is that the complexity of the problem of finding the effective conductivity of the network is significantly reduced. Instead of finding the effective conductivity by numerically decimating the GTN, one can estimate the conductivity geometrically using

where $\delta /2$ is the thickness of the soft penetrable shell surrounding the hard particles and *δ _{c}* is the smallest

*δ*for which the fillers at a given volume fraction form a percolating network. A model based on the critical path method,

^{8}as explained in detail in Section II, is used in this study.

Aspect ratio of nanofillers is identified as the primary factor that has a significant impact on the electrical conductivity of conductive polymer nanocomposites.^{1} The effect of aspect ratio on the conductivity has been extensively studied through experiments,^{11} modeling,^{8} and theoretical investigations.^{12,13} However, in these studies, it is often assumed that all the fillers are identical in size and shape to simplify the analysis.

Comparison of the simulations or theoretical predictions of the idealized models with experimental data is not straightforward. Many other factors such as the size distribution of fillers, changes in shape, size, and aspect ratio due to breaking of fillers during processing, and agglomeration of fillers lead to a system that differs significantly from the idealized model^{8} with uniformly dispersed identical spheroidal fillers.^{1} As a result, the idealized model may fail to represent the actual composite and deviation from experimental data can be observed. Therefore, understanding the effect of such parameters is highly important when comparing the results of the model with experiments, and more so, when using the model as a predictive tool.

Analytical studies on polydisperse rods, based on the mean-field assumption,^{9,13–15} have shown that the percolation threshold $\varphi c$ is extremely sensitive to polydispersity. It was found that $\varphi c$ is inversely proportional to the weight average of the length distribution $Lw=\u3008L2\u3009/\u3008L\u3009$, for penetrable and impenetrable systems of rods. Therefore, with increasing polydispersity, the $\varphi c$ decreases even if the mean length $\u3008L\u3009$ is fixed. The longer rods contribute strongly in establishing the percolating network, especially when the length difference between longer and shorter rods increases. The analytical models, although valid strictly in the mean-field sense and in the slender-rod limit (high aspect ratios), provide a powerful way to interpret the experimental and simulation results in a polydispersed system of elongated fillers like CNTs.

The effect of length polydispersity on a system of spherocylinders with a fixed diameter *D* is studied in Ref. 16 using Monte Carlo (MC) simulations. Spherocylinders are cylinders with diameter *D*, length *L*, and hemispherical caps of diameter *D* at the ends. For a polydispersed system of penetrable (soft-core) spherocylinders, $\varphi c$ is shown to be inversely proportional to *L _{w}*, as predicted theoretically in the slender-rod limit.

^{9,13–15}For a sufficiently large $\u3008L2\u3009/D$, a universal dependence of the form $\rho \delta c3=F(\u3008L2\u3009/\delta c)$ has been observed, where

*ρ*is the number density. Mutiso

*et al.*

^{17}have used a generalized form of the excluded volume model and percolation simulations of soft-core polydispersed cylinders to model the effect of finite

*L*/

*D*and size polydispersity on $\varphi c$. Their model is not limited by the slender-rod limit

^{14}and works for finite, but penetrable, rods as well. Chen

*et al.*

^{18}have presented a generalization of the excluded volume model

^{12}that can predict the percolation threshold for penetrable polydisperse spherocylinders. Recently, Meyer

*et al.*

^{19}studied the effect of polydispersity on percolation behavior of rod-like inclusions using theoretical and Monte Carlo simulation approaches. They modeled Inclusions as spherocylinders with a hard core and a penetrable shell and accounted for polydispersity in shell thickness of the fillers as well as length and diameter. Such analytical and simulation based studies have greatly enhanced our understanding of percolation behavior of penetrable and impenetrable polydisperse rodlike fillers.

Based on a generalized connectedness percolation theory, Otten and van der Shoot^{9} have theoretically predicted the percolation threshold for hard polydisperse platelike objects with soft penetrable shells. For monodisperse hard discs, $\varphi c$ is shown to be independent of the aspect ratio. However, their model predicts that $\varphi c$ of polydisperse hard discs depends on up to fourth moment of the size distribution, and so a strong influence of the diameter polydispersity is expected. However, the theoretical prediction is only valid in the limiting case of thin discs.

In this study, we analyze the effect of size polydispersity on percolation of randomly oriented prolate and oblate ellipsoids for a wide range of aspect ratios using MC simulations and compare our results with the theoretical predictions. The results for polydisperse prolate ellipsoids are shown to be in agreement with the theoretical predictions^{9} and polydisperse spherocylinders simulations.^{16} We show that the results for monodisperse and polydisperse oblate ellipsoids can be collapsed on a universal functional form when the aspect ratio is sufficiently large, but still much smaller than the thin plate limit. Therefore, the theoretical predictions based on connectedness percolation^{9} can be confidently used in practice where the aspect ratios of the fillers used is sufficiently large.

The paper is organized as follows. The tunneling-percolation model is presented in Section II. The MC results for polydisperse prolate and oblate ellipsoids representing CNTs and GNPs, respectively, are compared with the theoretical model in Section III, and with experimental results in Section IV, followed by conclusion in Section V.

## II. MODEL

The fillers are modeled as impenetrable spheroids with aspect ratio *α* that best captures the morphology of the filler. *α* is defined as the ratio of polar to equitorial semi-axes of the spheroid. For CNTs, the cylinder-like shape is obtained using $\alpha \u226b1$, while for GNPs the plate-like shape is obtained with $\alpha \u226a1$. The random distribution of fillers is generated using the random sequential addition (RSA) method.^{20} In this method, the desired fillers are added one by one such that the recently added particle does not collide with the previously added ones. The test of collision between two ellipsoids is implemented using a procedure proposed by Perram and Wertheim.^{21} The process is continued until the desired volume fraction is reached. Although simple, RSA suffers from a major drawback that for a given *α* there is a theoretical limit on the maximum volume fraction that can be attained by finite random trials.^{20} However, the maximum attainable volume fractions by RSA were found to be sufficient for the current study. It is also worth noting that this technique generates configurations that are not equilibrated, which can be problematic for configurations with high density of fillers. Nevertheless, the discrepancy is not significant at small volume fractions^{22} that are considered in the present work. MC relaxation technique is generally used to obtain equilibrated configurations. Also, percolation behavior of RSA generated monodisperse oblate ellipsoids with and without MC relaxations was found not to differ significantly in Ref. 8.

The position of the filler is chosen within a cube using three random numbers between 0 and *l*, where *l* is the edge length of the cube. Periodic boundary conditions are enforced along all directions. The random orientation of the polar axis of spheroids is generated with a uniformly distributed random variable on a unit sphere surface using $\theta =2\pi m$ and $\varphi =arccos(2n\u22121)$, where *m* and *n* are random numbers in $[0,1]$, while *θ* and $\varphi $ are the spherical coordinates giving the polar and azimuthal angles of a point on the unit-sphere surface, respectively.

In the simulations performed, the box size *l* is chosen such that it is at least 5 times the largest filler dimension and that the total number of fillers is at least 25 000. In this way, it is ensured that the finite-size effects does not significantly affect the percolation based results. The box is further divided into cubical sub-blocks of the size of the largest filler dimension. Such a subdivision narrows down the collision check only to the surrounding sub-blocks. The expensive ellipsoid collision check is bypassed if the distance between the centers of the two fillers is larger than the sum of their largest semi-axes, thus further improving the speed and efficiency of the collision check operation.

After generating the RSA configuration of fillers with a desired volume fraction $\varphi $, all fillers are coated with a soft penetrable shell of uniform thickness $\delta /2$. The objective is to find *δ _{c}* which is the minimum

*δ*for which a percolating cluster of the interconnected fillers spanning the entire domain is formed. Here, for simplicity, the underlying assumption is that the tunneling mechanism is isotropic, i.e., has no directional dependence.

The percolation check is performed using the tree-based union/find algorithm with path compression proposed by Newmann and Ziff.^{23} In this method, the hard core fillers with a penetrable shell of thickness $\delta /2$ are added one by one to the box. After adding each filler, the cluster information, stored using a tree datastructure, is updated based on the contacts established with the neighboring fillers through the penetrable soft shell. The objective is to check if for a given $\varphi c$ and *δ*, a percolating cluster spanning the cube in a given direction is formed. The details of the percolation check algorithm and its implementation can be found in Ref. 23.

Starting with an initial guess for *δ*, for a RSA generated filler configuration at a desired volume fraction $\varphi $, the percolation check is performed as described above. If the system is found to percolate, then *δ* is decreased, and if percolation is not achieved, then *δ* is increased. Therefore, *δ _{c}* is obtained to within a desired accuracy using a bounded root finding methods such as the bisection method.

## III. RESULTS AND DISCUSSION

### A. Polydisperse prolate ellipsoids

The nanofiller mixing processes such as sonication and shear mixing, which improve the filler dispersion, lead to significant changes in the size distribution of the CNTs Refs. 7 and 24. Experimental studies that have analyzed CNT size distribution have shown that Weibull distribution effectively captures the polydispersity.^{25,26} Such distributions with significantly higher moments have a notable effect on the percolation behavior, and hence the effective conductivity of the composite.

As mentioned before, the effect of polydispersity on rodlike fillers has been investigated using analytical methods^{9} and MC simulations.^{1,16,17} In this section, the results for polydisperse prolate ellipsoids are compared with the theoretical predictions in the slender-rod limit. The analysis is along the same lines as that of MC simulations of spherocylinders in Ref. 16, and we expect to see similar trends for prolate ellipsoids, especially at higher aspect ratios.

Using the liquid-state integral equation theory, the critical percolation volume fraction $\varphi c$ for polydisperse impenetrable rods with diameter *D* and length *L* is given as^{9}

Then,

The theoretical prediction, applicable in the slender-rod limit, i.e., $L/D\u226b1$, was used in Ref. 16 to analyze the percolation behavior of impenetrable spherocylinders. It was shown that, for different polydisperse systems with the same $\u3008L2\u3009/D$, the $\rho \delta 3$ values are essentially independent of the actual size distribution. Moreover, with increasing $\u3008L2\u3009/D$, the simulation data tend to converge towards a common functional form which approaches Eq. (4) for $\u3008L2\u3009/\delta c\u226b1$.

Similar arguments can be extended for monodisperse and polydisperse systems of prolate fillers studied here (Fig. 1). The semi-axes of the prolate ellipsoids are given as ${L/2,D/2,D/2}$. Prolate ellipsoids following a bimodal distribution $f(L)=p\delta (L\u2212L1)+(1\u2212p)\delta (L\u2212L2)$ are used to account for length polydispersity, where *p* is the number fraction. All prolate ellipsoids are assumed to have the same *D*, and all lengths are normalized with *D*. The aspect ratio *α* is given as *L* / *D*. Similar to spherocylinders,^{16} different size distributions with the same $\u3008L2\u3009/D$ collapse on a single curve as shown in Fig. 2. Moreover, for larger $\u3008L2\u3009/D$, the collapse of data suggests a common functional form that is independent of $\u3008L2\u3009/D$. Therefore, as expected, the quasi-universal behavior for a given $\u3008L2\u3009/D$ that was reported for spherocylinders^{16} is observed for polydisperse prolate ellipsoids as well. A similar analysis on polydisperse oblate ellipsoids is presented next.

### B. Polydisperse oblate ellipsoids

The inverse proportionality between $\varphi c$ and aspect ratio that is observed for prolate ellipsoids is absent in the case of oblate ellipsoids.^{9,27,28} After an initial drop in $\varphi c$ with increasing aspect ratio, $\varphi c$ saturates to a limiting value independent of the aspect ratios. This behavior of a system of disc-like fillers is suggested to be the outcome of interference of the percolation transition with the isotropic-nematic phase transition.^{28} However, $\varphi c$ is shown to be significantly altered in the presence of diameter and thickness polydispersity.^{9}

The analytical expression for $\varphi c$ for polydisperse hard discs based on connectedness percolation theory and second-virial approximation to the pair connectedness function was derived in Ref. 9 as

where *L* is the disc thickness and *D* is the disc diameter. Higher moments of the size-distribution appear in the expression, unlike in the slender-rod case, where only the first and second moments appear. Equation (5) is applicable when $\delta c\u226aD$ and when *L* and *δ _{c}* are of the same order of magnitude. Unlike for rodlike particles, the second-virial approximation is supposed to break down for platelike and spherical particles as the contribution of the higher order terms cannot be neglected.

^{9}However, for $\delta \u226aD$, it turns out to be a reasonable approximation for spheres, and hence is presumed to be applicable to platelike particles as well when deriving Eq. (5).

^{9}Therefore, Eq. (5) is expected to work at least qualitatively within the second-virial approximation. Moreover, $\varphi c$ predicted using Eq. (5) was shown to match well with experimental data

^{29}where the GNP size polydispersity information was obtained using dynamic light scattering.

For a polydisperse distribution of hard discs, $\varphi c=\rho \pi \u3008D2\u3009L/4$. Then, using Eq. (5)

where $1/A$ represents the expression in the square brackets of Eq. (5). Therefore, in the thin disk limit, scaling of the form Eq. (6) is expected.

Using the percolation-tunneling model of ellipsoidal fillers, we have simulated monodispersed and polydispersed configurations of oblate ellipsoids (Fig. 3). The semi-axes of oblate ellipsoids are represented as ${D/2,D/2,L/2}$, and the aspect ratio is given as $\alpha =L/D$. The results are presented by plotting $\rho \delta c3$ as a function of $A/\delta c$. The plots for monodisperse cases with different aspect ratios are shown in Fig. 4. With increasing aspect ratio, the curves are observed to converge on to a universal functional form $\rho \delta c3=F(A/\delta c)$. Equation (6) for the limiting case of thin hard discs is also presented for comparison. The simulated oblate ellipsoids approach the thin disc behavior when *α* is very small and $\delta c/D\u226a1$.

A bidisperse system of oblate ellipsoids described by a bimodal distribution $f(D)=p\delta (D\u2212D1)+(1\u2212p)\delta (D\u2212D2)$ is simulated, where $D1>D2$ and $p\u2208[0,1]$ is the number fraction of the larger ellipsoids. To simplify the analysis and to focus on polydispersity in the larger dimension, all oblate ellipsoids in the simulations have the same *L,* i.e., all GNPs are assumed to have the same thickness. Also, all dimensions are normalized with *L*.

The effect of polydispersity on the critical distance *δ _{c}* can be examined by using the reduction factor

*R*, which is the ratio of

*δ*to $\delta c0$, where $\delta c0$ is the critical distance for the monodisperse system of ellipsoidal fillers with semi-axes equal to the mean semi-axes of the polydisperse case, i.e., ${\u3008D\u3009/2,\u3008D\u3009/2,L/2}$. For the bimodal distribution, $\u3008D\u3009=pD1+(1\u2212p)D2$.

_{c}An analytical expression for the reduction factor *R*_{0} is obtained using Eq. (5) as

where

is the critical distance for monodisperse hard discs at a volume fraction $\varphi c$. Both *δ _{c}* and $\delta c0$ are for the same volume fraction $\varphi c$.

The reduction factors *R* obtained from simulating oblate ellipsoids following the bimodal distribution for $D2/L=5$ and $D1/D2=4$ are plotted in Fig. 5 as a function of *p*. The plots are generated for different total volume fraction $\varphi c$ of the fillers. The maximum $\varphi c$ used is limited by the RSA procedure used to generate the realizations. It can be observed that *R* obtained from the simulations becomes qualitatively similar to *R*_{0} as $\varphi c$ increases. Similar trends can be observed in Fig. 6 for $D2/L=5$ and $D1/D2=2$. While for small $\varphi c$, the condition necessary for the theoretical prediction *R*_{0} to hold, i.e., $\delta c\u226aD$ is not satisfied, and hence there is a large difference between *R* and *R*_{0}. With increasing $\varphi c,\u2009\delta c/D$ decreases and the theoretical result is approached by the simulated data. Apart from $\varphi c$, the aspect ratio $D2/L$ is also expected to affect *R*. The effect of $D2/L$ on *R* is shown in Fig. 7, where *R* is plotted as a function of *p* for $\varphi c=0.025,\u2009D1/D2=2$, and $D2/L={2,4,10,20}$. It is again observed that with increasing $D2/L$, *R* approaches the analytical prediction *R*_{0}.

The data for oblate ellipsoids following a bimodal distribution with $\alpha 1=L/D1=1/10,\u2009\alpha 2=L/D2=1/2$, and $p={0,0.5,1}$ obtained for a range of $\varphi c$ is shown in Fig. 8, where $\rho \delta c3$ is plotted as a function of $A/\delta c$. The data for monodisperse oblate ellipsoids with $D/L=A/(10\pi +12)/L$, where *A* is calculated for the bimodal distribution given above at *p* = 0.5 is also shown. For $\alpha 1=L/D1=1/10,\u2009\alpha 2=L/D2=1/2$, and *p* = 0.5, $D/L=7.435$. For a given $A/L$, the $\rho \delta c3$ against $A/\delta c$ curves collapse well and can be concluded to be independent of the specific size-distribution. The data for a wide range of monodisperse and bidisperse oblate ellipsoids simulated in this study are collectively shown in Fig. 9. For oblate ellipsoids with sufficiently large $A/L$, the curves can be observed to collapse well, suggesting a universal functional form $\rho \delta c3=F(A/\delta c)$. Also, as $A\u226b\delta c$ and for a system of hard discs, the scaling function approaches Eq. (6).

## IV. COMPARISON WITH EXPERIMENTS

Based on the observations above, we have analyzed some experimental data available in literature for CNT and GNP nanocomposites. As we are interested in the effect of polydispersity, we have chosen the experimental data from sources that have reported the actual size distribution of the nanofillers. We have used the data from Castillo *et al.*^{7} for multi-walled CNTs (MWCNTs) and from Tkalya *et al.*^{29} for GNPs.

The experimental data can be compared with the tunneling-percolation based theories in the following two ways:

The first approach is based on continuum percolation of soft-shell hard-core fillers (or cherry-pit model), where the soft shell thickness corresponds to a tunneling-induced lengthscale *λ*. The second approach uses the exponential dependence of tunneling conductance (Eq. (1)) on the filler separation, but does not impart a sharp percolation cutoff and involves the tunneling lengthscale *ξ*, which is expected to be proportional but not necessarily equal to the lengthscale *λ* of the first approach.

It is clear from Figs. 2 and 9 that predictions based on the percolation connectedness theory are accurate for higher aspect ratios that are generally used in experiments. Therefore, Eqs. (3) and (5) are used to compare the model predictions with experimental data in the following. The objective is to estimate *λ* and *ξ* using the size polydispersity information along with the $\sigma \u2212\varphi c$ data from the experiments. Consistent estimations of *λ* and *ξ* within the expected range^{8} would indicate the importance of the polydispersity based formulation.

### A. Polydisperse MWCNTs

Following the first approach, $\varphi c$ is obtained from Eq. (3) as $\varphi c=D2/2\lambda Lw$, where $\lambda /2=\delta c/2$ is tunneling soft-shell thickness. Here, *L _{w}* can be calculated from the MWCNT length distributions obtained experimentally. The polydispersity and conductivity data provided for MWCNTs from five different suppliers presented in Ref. 7 are used. We assume that all MWCNTs have the same

*D*that is provided by the manufacturers.

^{7}Therefore, the only unknown,

*λ*, can be estimated using the experimentally obtained $\varphi c$. The length distributions for the five cases are extracted from Fig. 1 of Ref. 7 and are named as A–E following the same order therein. The

*λ*values estimated using the experimentally obtained $\varphi c$ and

*L*are summarized in Table I. Interestingly, for all cases, the estimated

_{w}*λ*falls within 7–11 nm and is comparable to the MWCNT diameter

*D*.

Samples . | D (nm)
. | $L50$ (nm) . | $\u3008L\u3009$ (nm) . | L (nm)
. _{w} | $\varphi c$ (vol. %) . | λ (nm)
. |
---|---|---|---|---|---|---|

A | 7.8 | 735 | 974 | 1523 | 0.23 | 8.82 |

B | 10.5 | 770 | 1104 | 1905 | 0.42 | 6.91 |

C | 10 | 1341 | 1604 | 2393 | 0.19 | 10.87 |

D | 10.5 | 727 | 1027 | 1772 | 0.34 | 9.25 |

E | 6 | 332 | 382 | 484 | 0.34 | 10.82 |

Samples . | D (nm)
. | $L50$ (nm) . | $\u3008L\u3009$ (nm) . | L (nm)
. _{w} | $\varphi c$ (vol. %) . | λ (nm)
. |
---|---|---|---|---|---|---|

A | 7.8 | 735 | 974 | 1523 | 0.23 | 8.82 |

B | 10.5 | 770 | 1104 | 1905 | 0.42 | 6.91 |

C | 10 | 1341 | 1604 | 2393 | 0.19 | 10.87 |

D | 10.5 | 727 | 1027 | 1772 | 0.34 | 9.25 |

E | 6 | 332 | 382 | 484 | 0.34 | 10.82 |

The effect of polydispersity is apparent when we observe the large difference between $\u3008L\u3009$ and *L _{w}* for each case. In their study, Castillo

*et al.*

^{7}used the continuum percolation based formula $\varphi c=D/2L50$ for theoretical estimation of $\varphi c$, where $L50$ is the median of the MWCNT length distribution. However, such a form clearly over-predicts $\varphi c$ because the simple formula accounts neither for tunneling nor polydispersity (Fig. 10). A somewhat better match is obtained if $\u3008L\u3009$ is used instead of $L50$, i.e., $\varphi c=D/2\u3008L\u3009$ which can be obtained as a special case of Eq. (3) for a monodisperse CNT distribution with lengths $\u3008L\u3009$, diameter

*D*, and $\delta c=\lambda =D$. However, using Eq. (3) with

*L*and the $\lambda \u22489.5$ nm (obtained from the average of

_{w}*λ*for all cases in Table I) provides excellent comparison with the experimental data, as shown in Fig. 10.

Next, the second approach based on Eq. (2) is used to compare the model results with experimental data. For this purpose, we need conductivity *σ* as a function of $\varphi c$ which is provided in Fig. 2 of Ref. 7 for cases A–E. Using $\varphi c$ values, the corresponding *δ _{c}* is obtained using Eq. (3), i.e., $\delta c=D2/2\varphi cLw$. Then, using

*ξ* can be estimated from the slope of the line fitted to points plotted on a $ln\sigma \u2212\delta c$ scale. The line is fitted only through the points above the percolation threshold where the critical path assumption is valid due to minimal effect of finite matrix conductivity.^{8} The points corresponding to volume fractions below the percolation threshold were observed to deviate significantly from the predicted behavior. The estimated *ξ* values for cases A–E (Fig. 11), summarized in Table II, are found to fall within 1.5–4 nm, which is in agreement with the expected *ξ.*^{8} Using the polydispersity information along with Eq. (3) have sufficiently narrowed down the range for *ξ*. The estimated *ξ* can be used as an input to predict *σ* of MWCNT nanocomposite for a given $\varphi c$.

### B. Polydisperse GNPs

Tkalya *et al.*^{29} used dynamic light scattering technique to obtain diameter distribution of GNPs and have shown that theoretical predictions using Eq. (5) are in qualitative agreement with experimentally obtained $\varphi c$ for polydisperse GNPs. The theoretical prediction is obtained by ignoring the thickness polydispersity and focusing only on the diameter polydispersity of the GNPs. Graphene sheet thickness *L* and tunneling-induced soft-shell thickness $\lambda /2$ remain as fitting parameters because both are not known *a priori*. However, by fitting, they have found $\lambda /L\u22484.5$ which correctly captures the trends and qualitatively agrees with the experimental data. This is the first method of comparison as described above. Their results show that size polydispersity plays an important role in determining the percolation threshold.

Following the second method of comparison, Eq. (9) used along with Eq. (5), which relates *δ _{c}* with $\varphi c$ for polydisperse discs, gives

Therefore, plotting $ln\u2009\sigma $ as a function of $8\u3008D2\u3009/\varphi cA$, we can estimate $\xi /L$ from the slope of the line fitted to the experimental data (Fig. 12) for different samples from Ref. 29. Similar to the MWCNTs case above, only points above the percolation threshold are considered to obtain the fits.^{8} Points below the percolation threshold can be observed to deviate from the fitted line in Fig. (12) for some cases. The diameter distributions reported therein are used to calculate *A* and $\u3008D2\u3009$. The estimated $\xi /L$ are summarized in Table III. The GNP sheet thickness *L* depends on the degree of exfoliation among other factors and, although not exactly known, is expected to be between 0.3 and 2 nm.^{29} For the first three cases, $\xi /L\u223c1$, which results in an admissible estimate of *ξ* within 1–2 nm. For the last case, however, $\xi /L=0.6$, which results in a considerably smaller estimate of *ξ*. A similar issue with case *D* is reported in Ref. 29 where the experimentally obtained $\varphi c$ is much higher than the theoretical prediction from Eq. (5). For all other cases, the theoretical predictions are in excellent quantitative agreement with the experiments. Thus, by correctly accounting for the size polydispersity, accurate estimates of *ξ* confined to a narrower range are obtained. Using *ξ* values estimated above, Eq. (10) can be used to obtain conductivity of GNP nanocomposite at a given $\varphi c$.

Samples . | A . | A-LC . | A-HC . | B . |
---|---|---|---|---|

$\xi /L$ | 0.8 ± 0.2 | 0.8 ± 0.4 | 1.1 ± 0.1 | 0.6 ± 0.2 |

Samples . | A . | A-LC . | A-HC . | B . |
---|---|---|---|---|

$\xi /L$ | 0.8 ± 0.2 | 0.8 ± 0.4 | 1.1 ± 0.1 | 0.6 ± 0.2 |

Both approaches used herein provide reasonable estimates of *λ* and *ξ* within the expected range. However, the tunneling lengthscale *ξ* in the second approach is precisely defined through Eq. (1) compared to *λ* in the first approach which is somewhat vague. Another estimate of *λ* can be obtained from the crossover point between tunneling dominated and matrix dominated regimes using fits obtained using Eq. (9). Using this apparent percolation threshold would result in lower *λ* values.^{8} However, fits based on Eq. (9) in Figs. 11 and 12 are observed to consistently over predict the effective conductivity near and below the percolation threshold. Therefore, although this apparent percolation threshold based on Eq. (9) would give lower values of *λ*, it may not be accurate as Eq. (9) does not seem to perform well near the percolation transition region.

It should also be noted that although polydispersity is accounted for in our model, we have neglected the effect of particle agglomerations and segregation in this study. Incorporating such effects in the model would involve quantifying the effect of anisotropic spatial distribution of fillers on the $\delta c\u2212\varphi c$ relations. This would help in further refining the estimates of *λ* and *ξ*, and improve predictability of the models.

Theory, simulations, and experiments have all confirmed that the percolation threshold is very sensitive to the degree of polydispersity. Therefore, along with conductivity measurements, a detailed characterization of nanofiller size distribution and spatial distribution should be reported in experimental studies on conductive polymer nanocomposites.^{1,7,29}

## V. CONCLUSIONS

In summary, we have performed MC simulations of polydisperse prolate and oblate ellipsoids using the critical path based tunneling-percolation model. For polydisperse prolate ellipsoids, provided a sufficiently large $\u3008L2/D\u3009$, the quasi-universal behavior is obtained, as shown for polydisperse spherocylinders in Ref. 16. For polydisperse oblate ellipsoids, the data are observed to follow a quasi-universal functional form for a sufficiently large $A/L$. Moreover, the quasi-universal functions approach the second-virial approximation based theoretical predictions^{9} for rodlike and platelike fillers as the volume fraction increases.

As the connectedness percolation theory^{9} is in good agreement with MC results at higher aspect ratios, experimental data on polydisperse MWCNTs^{7} and GNPs^{29} are analyzed using theoretical predictions. The tunneling lengthscale *ξ* predicted using Eq. (2) along with the polydispersity and conductivity data from experiments is found to be well within the expected range. The estimated *ξ* can be used as an input to predict *σ* of MWCNT or GNP nanocomposites for a given $\varphi c$.

Apart from polydispersity, other factors such as clustering of nanofillers, curling of CNTs, crumpling of GNPs, etc., are also expected to affect the percolation behavior.^{13} The analytical models have provided some important insights regarding the impact of such properties on percolation threshold, but are applicable, in general, for the limiting cases. Further work through modeling is in progress to understand the effect of such factors on the percolation threshold and the effective electrical conductivity of nanocomposites.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge the use of the Taub cluster resources provided under the CSE program at University of Illinois. This work was financially supported by the NSF Center for Novel High Voltage/Temperature Materials and Structures [NSF I/UCRC (IIP-1362146)].