Percolation of rigid fractal carbon black aggregates

Carbon particles from pyrolysis have long been used as fillers in elastomer composites, for example, to improve the mechanical properties of tires and gaskets^1,2 and to lend polymers electrical conductivity in materials for electrostatic interference shielding or as positive temperature coefficient materials.^3,4 In response to a need for robust elastic sensor and actuator materials, the coupling of mechanical and elastic properties in elastomer–carbon particle composites has moved into the focus of research.^5–9 In all these types of applications, the filler has a significant impact on the composites’ electrical and mechanical properties if the filler forms a connected network throughout the material, i.e., if the concentration of the filler is close to or above the percolation threshold.^10,11 A thorough understanding of the percolation process of filler materials is therefore of immediate practical relevance.

The formation of percolating networks from rod-shaped particles, such as carbon nanotubes, and plate-like particles, such as graphene, has been studied extensively.^12,13 The effect of the filler particles’ geometry and simple interaction potentials on the percolating network has been studied theoretically and by means of computer simulation for spherical, rod-like, and plate-like fillers as well as their mixtures.^14–32 To our knowledge, however, the percolation behavior of fractal structures, such as carbon black (CB), has not yet been studied. Experimental work on the effect of CB fillers on mechanical^3,33,34 and electrical^35–37 properties of elastomer composites focused on the type of CB and its concentration but not on the role of the fractal shape in percolation. Systematic experimental studies of CB percolation are complicated by the small size of the quasi-spherical primary particles (typically below 100 nm) and the wide size distribution of the fractal aggregates that they form. Figure 1(a) shows an idealized cluster that represents a CB aggregate as it was observed by transmission electron microscopy (TEM) [Fig. 1(b)]. Primary particles are relatively uniform and quasi-spherical, but they form aggregates with wide size distributions as shown in Fig. 1(c). The aggregates form weak bonds and arrange into agglomerates. Experimental techniques alone cannot provide a holistic image of CB percolation at present. We therefore constructed a simulation model that is based on experimental data obtained from a standard type of branched acetylene CB and used it to study percolation. We finally compared the theoretical results to our experimental observations.

In the following, we present a detailed analysis of percolation in suspensions of fractal aggregates including polydispersity as well as interactions. Our main conclusions are as follows: (1) The percolation threshold depends quasi-universally on the averaged radius of gyration of the aggregates, i.e., the percolation threshold can be inferred from a single macroscopic number characterizing the aggregate ensemble, and (2) fractal particles have a lower percolation threshold than rods; hence, contrary to intuition, CB can be more efficient than, e.g., carbon nanotubes as a filler material.

The structure of CB is commonly described at three different levels. The primary CB particles are densely coalesced spheroids consisting of turbostratic graphene layers. The diameter of the primary spheres depends strongly on the production technique and ranges from 5 to 100 nm. Within a specific grade of CB, the primary particle size typically does not vary strongly. TEM images of the CB we used for our experiments (see details below) show a narrow distribution of primary particle diameters at roughly 40 nm.

At the first level of aggregation, the primary particles form primary aggregates, which are the smallest dispersible units, i.e., they consist of strongly bound primary particles, which can only be separated by fracture. The primary aggregates are self similar. The size of an aggregate is related to its number of primary particles through a power law. The exponent governing the aggregate growth is called mass fractal dimension d_m. For the majority of industrially common CB, the mass fractal dimension varies between values of 2.2 and 2.8.³⁸ 

At the second level, the primary aggregates form secondary aggregates (we call those agglomerates henceforth to prevent confusion) that are clusters of primary aggregates connected through weaker van der Waals (vdW) forces. As the surface of the primary aggregates can be extremely large, the sparse agglomerates form a bonded unit themselves, which is the backbone that enhances the mechanical strength of polymeric composites through filler materials. The agglomerate growth is statistically described by a power law with a smaller exponent d_s, which has been reported as roughly 1.8 for a variety of commonly used CB.^38,39 Remarkably, this value is a good match to the fractal dimension of agglomerates generated by cluster–cluster aggregation (CCA). Previous studies have thus used networks formed by CCA clusters on a lattice for a quantitative comparison to the conductivity of CB–rubber composites close to the percolation threshold.^40,41 Figure 2 schematically illustrates the different hierarchic layers of of carbon black.

It is often difficult to experimentally distinguish agglomerates from aggregates once they have formed. However, it is useful to treat them as separate structural entities because materials are formed by first dispersing CB as much as possible (such that ideally, only aggregates remain connected) and to then add them to a matrix, for example, a liquid elastomer precursor. Such liquids are imperfect solvents for CB, which leads to a certain level of attraction between the aggregates and slow agglomeration that affects material properties.

As the primary particles are approximately spherical, we model them as spheres. To take into account the strong bonds inherent to the primary structure, we model an aggregate as a rigid as well as impenetrable body. For simplicity, we assign a universal diameter d to our primary particles, which is a decent approximation to the almost monodisperse primary particle size distribution of many carbon blacks, including the one we use for experimental comparison. In order to create suitable primary aggregates, a formation scheme is required, which reproduces the appropriate power law statistics. Diffusion-limited aggregation (DLA) has been studied extensively in theory and simulation, and it is known to produce aggregates characterized by a fractal dimension of roughly 2.5. This value is consistent with the fractal dimension of 2.57 that we calculated from gas porosimetry experiments on our CB using Brunauer–Emmett–Teller (BET) theory as shown in Fig. 3.

We create aggregates by 3D continuum diffusion-limited aggregation of spheres that model primary particles. During the aggregation process, single primary particles diffuse freely, one-by-one, until they hit a pre-existing structure starting with a fixed primary particle as seed at the origin. When two primary particles overlap, they become rigidly linked to form an aggregate. We quantify the size of an aggregate in terms of its radius of gyration

where N denotes the number of primary particles per aggregate, r_S is the center of mass, and r_i is the position of primary particle i. We compute the average r_g for aggregates of a fixed number of primary particles N and obtain the mass fractal dimension of the DLA aggregation mechanism by means of a power law fit of the form

for large aggregates (see Fig. 4). In accordance with the literature (e.g., Ref. 42), we find a fractal dimension of roughly 2.5 for aggregates of more than 100 primary particles.

We take ensembles of these aggregates and carry out Monte Carlo simulations, sampling their positions and orientations within a box with periodic boundary conditions. Each simulation was performed with at least 1000 aggregates corresponding to up to several million primary particles. Two aggregates are considered “connected,” i.e., part of a conductive cluster, if the distance between at least one pair of primary particles of the corresponding aggregates is smaller than a threshold distance d + ϵ. The connectivity shell size ϵ roughly models the tunneling length of electrons between aggregates so that we can relate the percolating cluster to a conductive network. For a given aggregate ensemble, we vary the size of the simulation box and observe whether a system spanning connected component, i.e., a percolating cluster, exists. Averaging over a couple of thousand independent realizations of the system, we determine an ensemble average of the percolation probability. We define the percolation threshold p_c as the volume fraction for which a sigmoidal fit to the percolation probability returns 0.5. Naturally, periodic boundary conditions are not a perfect substitute of an infinite system so that finite size effects have to be discussed (see the  Appendix A).

In this paper, we analyze how the percolation threshold depends on the aggregate ensemble Γ, the connectivity shell size ϵ, and small added interactions between the aggregates.

Using the procedure described above, we systematically probed the influence of the aggregate composition on the percolation threshold. Ultimately, as the industrial application of composite materials often requires thrifty supplement of the filler material, a small critical filler fraction and high conductivity of the percolating network are generally desirable.

There are two different factors that impact the location of the percolation threshold in the simulation of aggregates with purely steric interaction—(A) the size distribution of aggregates and (B) the connectivity shell size ϵ. In the following, we will analyze the influences of these parameters separately and draw conclusions toward optimization.

We examine mono-disperse, bi-disperse, uniform, normal, and log-normal distributions of the number of primaries per aggregate N. The induced distribution of ⟨r_g⟩ is the composition of the specific distribution with the power law describing the aggregate growth [Eq. (2)]. Therefore, all distributions are inherently skewed. Experimentally, the distribution of Stokes diameters in realistic carbon black is often modeled as a log-normal distribution, which matches our experimental observations (see below). In order to compare the experiment to the simulations, we need to assume a relationship between the measured Stokes diameter and the radius of gyration of the associated aggregate. As both measures describe dimensions of effective spheres, we assume proportionality. For a first set of simulations, we choose the connectivity shell diameter as ϵ = 0.2d, i.e., 20% of the primary particle diameter. This value is essentially arbitrary and chosen as small as possible while still allowing to comfortably equilibrate samples of large aggregates. Naturally, as ϵ becomes arbitrarily large, the structure of an individual aggregate becomes irrelevant, while for arbitrarily small connectivity shells, we observe jamming rather than percolation. Thus, the range of reasonable connectivity shell sizes is limited.

In order to compare different distributions, we define for a given ensemble of aggregates $Γ=(γi)i∈N$ the weighted average radius of gyration $⟨rg⟩Γ$ as

where n_γ is the number of primary particles within aggregate γ. We finally determine the percolation threshold by Monte Carlo Simulation of a fixed aggregate ensemble in boxes of varying sizes. Figure 5 illustrates the findings.

The simulation results strongly suggest that the functional form of the aggregate size distribution has no influence whatsoever on the percolation threshold. Instead, the percolation threshold seems to depend exclusively on the mass weighted average radius of gyration of an aggregate ensemble. The percolation threshold decays monotonically as a function of $⟨rg⟩Γ$⁠. Thus, large aggregates induce low percolation thresholds.

However, the range of accessible radii of gyration is too small to confirm an asymptotic power law relationship. Ensembles of aggregates with $⟨rg⟩Γ>10$ are increasingly hard to equilibrate due to the intricate individual structure of each aggregate. To make sure that the quasi-universality prevails for different connectivity shell sizes, we analyzed the same ensembles for a connectivity shell size of ϵ = d/2. The result is again a unique relationship at an overall lower level (see Fig. 5).

The quasi-universal dependence of the percolation threshold on specific moments of size distributions has been observed before in polydisperse ensembles of rod-like particles²⁶ as well as fracture networks.⁴³ In Ref. 26, spherocylinders with different length distributions were examined and percolation thresholds for all distributions collapsed onto a single function of the average aspect ratio ⟨ζ⟩. To explain this behavior theoretically, a variation of the second virial approximation of continuum percolation theory has been studied, which predicted a quasi-universal relationship though quantitatively inaccurate. The second virial approximation assumes that the connected components grow tree-like until a percolating network emerges, i.e., finite clusters do not contain loops.⁴⁴ This has been found to be a decent approximation for systems of rods with large aspect ratios as in those systems the largest finite clusters remain tree-like up to the critical density. Indeed, as the percolation threshold decreases monotonically with the aspect ratio, we eventually reach arbitrarily low densities so that thermodynamic correlations become entirely negligible. In this limit, truncating a low density expansion at second order is always justified, i.e., the second virial expansion is expected to be asymptotically accurate. This applies to DLA aggregates as well as rods. However, rods require aspect ratios ζ > 100 before the second virial approximation becomes adequate and rods with a smaller aspect ratio have consistently higher percolation thresholds than DLA aggregates of the same radius of gyration (see Fig. 6). This is a strong indication that carbon black is more efficient in forming conductive nano-composites than, for instance, carbon nanotubes (on top of the obvious financial aspect). Yet, in their natural form, carbon black aggregates are limited in size. The major advantage of rod-like particles is therefore the accessibility of very large radii of gyration eventually compensating their minor efficiency.

Yet, the low density regime bears the disadvantage that the percolating cluster will be extremely scarce as well. The corresponding conductive network will heavily rely on specific nodal links essential to charge transport, i.e., there will be one or at most a few “red links,” which carry the entire current. Those links can easily be severed through external influences like shear deformation or temperature variations. Therefore, scarce networks allow for highly sensitive sensors and hence attracted a lot of attention. Previous studies have experimentally combined conductivity measurements and rheological measurements on molten CB composites, demonstrating that sufficiently strong shear destroys the networks formed by carbon fibers or carbon nanotubes.^47–47 In contrast to that, the conductivity of CB composites can even increase by shearing. Moreover, the conductivity of CB based composites was found to be rather insensitive to temperature fluctuations.⁴⁸ All these observations are consistent with the network structures induced by the different filler morphologies. The applicability of the second virial approximation is accompanied by a lack of resilience against breaking contacts. In contrast to that, the percolating network of DLA clusters features a more complex backbone with several paths that carry almost equal currents (see Fig. 7). It is hence more robust against external influences and thus more suited to technological applications that require long lifetimes.

We return to the quasi-universal behavior of the percolation threshold and note that, heuristically, the authors of Ref. 26 observed that a function of the form

with the aspect ratio ζ and constants A, B, and C excellently fit the simulation results asymptotically, resulting in the $1ζ$-dependence predicted by the second virial approximation. However, a function of that form cannot adequately describe the percolation threshold of DLA clusters. Yet, for a fixed connectivity shell size, a function of the form

allows for a remarkably accurate fit, which is displayed in Fig. 5. Importantly, this fit suggests an asymptotic $⟨rg⟩−1$ decay of the percolation threshold matching the behavior of the spherocylinders. However, simulations of much larger DLA clusters would be required to confidently resolve this regime.⁴⁹ 

Larger connectivity shells reduce the percolation threshold. Yet, the influence of the connectivity shell diminishes when the aggregates grow large. DLA aggregates are too dense in their core to heavily intertwine so that a larger portion of volume corresponding to conductive contact is effectively blocked. The accessible and connected volume varies less rapidly with ϵ for large aggregates. Thus, ⟨r_g⟩ and ϵ are coupled to each other, and the strength of this coupling is connected to the aggregation method. Figure 6 illustrates the simulation results. In order to characterize the coupling, we consider the crossover length

and, inspired by Eq. (4), fit our simulation data to

The parameters A = 0.62, B = 2.63, and C = −0.89 lead to remarkably good agreement for all simulations with r_g > 1. Thus, representing our simulation results for all distributions and connectivity ranges in terms of ξ causes all percolation thresholds to coincide on a common curve (see Fig. 8). This is true regardless of the parameters A and C so that the crossover length is defined on the basis of a one parameter fit to roughly 100 different aggregate ensembles.

The model above only considers steric repulsion. Common attractive interactions may additionally induce agglomeration of the filler particles. Flandin et al. identified an interplay between short-ranged vdW attraction and long-ranged Coulomb repulsion as a dominating factor in the conductivity of CB–epoxy composites.⁵⁰ Other works reported a crossover from a diffusion-limited to a reaction-limited aggregation mechanism as the filler fraction was increased.³⁹ Attractive interactions may hence significantly alter the volume fraction of the percolation transition. Structural analysis of the agglomerates can provide insight into the underlying aggregation mechanism.

Industrial CB–elastomer composites often contain multiple components that are processed at different temperatures and undergo multiple chemical reactions. It is useful to reduce the number of components in order to study the possible formation of agglomerates systematically. To this end, we prepared a simplified composite by mixing our standard CB from above with a silicon elastomer [polydimethylsiloxane (PDMS), SYLGARD^TM 184] that can be cross-linked at moderate temperatures through a hydrosilylation process that does not affect CB chemistry. The mixture was prepared by vigorous stirring of the dry CB powder with the liquid PDMS precursor mixture, which yielded a highly viscous black paste.

Model composites were prepared using different volume fractions of CB while keeping all other processing parameters unchanged. The resulting materials were structurally analyzed using ultra-small angle x-ray scattering (USAXS) at a synchrotron (DESY in Hamburg, Germany) in order to examine the fractal structure of the composite material. Scattering from a fractal results in a characteristic power-law decay of the scattering intensity with the wave vector. A comprehensive overview can be found in Ref. 51. We have to distinguish two different types of fractals, mass fractals and surface fractals. Mass fractals describe how the mass of a structure scales with its size. For instance, the exponent we determined in Fig. 4 corresponds to a mass fractal as by construction the mass of our DLA aggregates grows linear with the number of primary particles. The scattering intensity I of a mass fractal on length scales smaller than the size R_M of the fractal object is governed by a power law

with the mass fractal dimension d_m. The mass fractal dimension is naturally smaller than the dimension of the embedding space, i.e., d_m < 3. In contrast to that, surface fractals describe how the surface of a fractal object grows with its size. The corresponding scattering intensity I_S scales like

with the surface fractal dimension d_s and the embedding dimension d. Importantly, d − 1 < d_s < d, which implies that scattering from a surface fractal in a three-dimensional space always corresponds to a power law decay with an exponent smaller than −3. We multiply our scattering intensities I by q³ because scattering from all hierarchical levels exhibits power-law-like behavior with exponents close to −3. Our representation stresses the differences between those exponents, thus emphasizing structural features. The relevant surface fractal in our system is the roughness of the primary particles. We can decide whether the primary source of scattering on a certain length scale is the mass fractal of the aggregate or the surface fractal of the primary particles by analyzing the slope of the scaled scattering intensity Iq³. Figure 9 shows the scattering of the composites for different CB loads.

We can clearly discern three hierarchical levels due to the local extrema in the scaled scattering intensity. The regime of negative slope corresponds to the primary particle surface fractal, which dominates the scattering signal on length scales between roughly 5 and 25 nm. This is consistent with the diameter of the primary spheres 40 nm that we reported above from TEM images. The mass fractal that emerges on smaller length scales corresponds to the internal structure of the primary particles, i.e., turbostratic sheets of carbon. A power law fit to that regime yields an exponent of roughly 2.5.

In both of the aforementioned regimes, no qualitative change in the scattering signal with the CB load is observed. This is expected as the primary particles are supposed to be indispersible. The same would be expected for the aggregates; however, the arrangement of aggregates changes. Due to the distinct polydispersity of the aggregates, the small q-regime in Fig. 9 is a superposition of the aggregate and the agglomerate level. As the carbon content increases, the ascent in this regime becomes steeper, indicating a lower fractal dimension and thus the formation of agglomerates. A quantitative analysis of the aggregate fractal dimension is unfeasible due to strong aggregate polydispersity (see also Fig. 1), which blurs the borderlines between hierarchic levels. Therefore, owing to the limited information from USAXS, we cannot assign definite geometries to the formed agglomerates. However, the experimental results are consistent with short-ranged vdW attraction between the CB aggregates that slowly diffuse in the PDMS matrix and form agglomerates. Large aggregate sizes and the high viscosity of typical CB formulations slow down diffusion. It is therefore reasonable to model the agglomeration process by adding a short-ranged attractive interaction between the aggregates to the simulation.

We apply a square-well potential of depth ΔE and range [d, d + ϵ] between primary particles in order to render conductive contacts between aggregates energetically beneficial. This supports the formation of agglomerates from existing aggregates. We analyze the impact of the strength of the interaction onto ensembles of constant primary particle number—the results are displayed in Fig. 10. We find that scaled versions of the previous fit, Eq. (7), accurately describe the percolation threshold for all interaction strengths, i.e., all lines in Fig. 10 are parallel. This observation supports the idea that the effective aggregate size is enhanced by a supplementary attraction. We conclude that if agglomeration occurs, the percolation threshold drops in general. In principle, the location of the percolation threshold depends strongly on the size distribution of the aggregates.

However, the homogeneous agglomeration enhancing potential we used only shifts the percolation threshold by a factor that does not depend on the higher order moments of the size distribution. Under realistic conditions, the agglomeration mechanism and force might depend on the aggregate size, which adds a level of complexity that exceeds the scope of this study.

The theoretical predictions made above suggest that larger CB agglomerates reduce the percolation threshold. This result is unaltered by a short range attractive interaction and holds even for broad aggregate size distributions.

It would be cumbersome to experimentally evaluate the prediction using different CB types. Industrial CB types not only differ in average agglomerate size but may also have different surface chemistry and electrical conductivity of the primary particles. Therefore, we devised experiments based on our standard CB type and changed the average aggregate size by ball milling, assuming that the primary particles were left unchanged. Thus, CB powders were ground in a planetary ball mill with the addition of isopropanol and loose grind balls. In contrast to high-energy ball milling, our approach limits the energy input and protects the primary particles because milling ceases when the viscosity reaches a certain level.

To ensure that the structures formed can be related to the experiment, we determined the hydrodynamic diameter distribution by analytical centrifugation after dispersing the CB in toluene. The diameter distribution of unground CB was monomodal with most aggregate sizes between 60 and 800 nm, i.e., 1.5–20 primary diameters. The measured distribution roughly followed a log-normal distribution as reported for CB aggregates in the literature⁵²—a corresponding fit is depicted in Fig. 1. Table I illustrates the decrease in the modal particle size during ball milling as a function of the mass of milling balls, i.e., milling intensity. The modal particle size consistently decreased with the intensity of ball milling and the overall distribution remained monomodal. We confirmed the effect of milling using additional USAXS measurements, the results of which are displayed in Fig. 11. Scattering from milled samples led to similar results in the q-range pertaining to primary particles compared to the unmilled sample (cf. Fig. 9) but changed for smaller q values. The left turning point shifts from 0.044 to 0.089 nm⁻¹, indicating a change in the aggregates and agglomerates. As aggregates become smaller, they contest the q-range previously dominated by primary particles. Scattering intensity from aggregates generally shifts to higher q values interfering with the unaltered primary particle scattering while reducing the intensity for small q’s. This observation is hence consistent with the analytical centrifugation (AC) measurements. Yet, aggregate polydispersity still prevents a comprehensive quantitative analysis of the fractals corresponding to the aggregate and agglomerate levels, respectively.

The conductivity of PDMS-CB composites as a function of the CB volume fraction is shown in Fig. 12. The results show that the composite’s conductivity increased with CB volume fraction for each carbon black set. Renormalization theory for continuum percolation predicts the conductivity to scale as a power law on approach of the percolation transition from the percolating regime

where η is the CB volume fraction, η_c is the volume fraction at the percolation threshold, t is the conductivity exponent, σ is the composite conductivity, and σ_∞ is a prefactor characterizing the asymptotic compound conductivity at η = 1. In experiments, exponents different from the known theoretical value for the rudimentary continuum percolation problem are frequently observed (see, for instance, Ref. 41). Real chemically interacting systems contain additional layers of complexity so that it is expected to leave the universality class of the purely geometric percolation problem. Therefore, the conductivity exponent is typically regarded as a fitting parameter. We measured the conductivity as a function of CB load for three different milling protocols and fitted the results to the scaling law (10). The resultant percolation thresholds are displayed in Table I. Reduction in the aggregate size through ball milling thus reduced the conductivity of the sample at constant CB volume fraction and shifts the percolation threshold to larger volume fractions. The values for the conductivity exponent should be interpreted with care because the number of fitted points and the range they cover are simply insufficient for a reliable quantitative analysis. However, the conductivity data for milling protocols “No” and “A” can be adequately fitted to a power law with specified conductivity exponent t ≈ 2.0, while this is not possible for protocol “B.” This means that networks formed by unmilled CB and mildly milled CB are consistent with the universal conductivity exponent, but the intensely milled sample is not. In previous studies,⁵³ a wider microscopic distribution of tunneling conductances could be linked to an anomalous increase in the conductivity exponent. We are able to observe the same phenomenon in simulations. We can adjust the distribution of tunneling conductances by varying the connectivity shell size without changing the aggregate ensemble.  Appendix B contains details on the conductivity computation in simulations and examples of conductivity calculations for the same aggregate ensemble showing both universal and non-universal conductivity exponents. The simulation results underline that non-universal exponents can be the consequence of microscopic reordering. The experimentally measured non-universal exponent thus indicates such a reordering. On top of that, the scattering data show the most pronounced differences in the regime of small wave vectors, i.e., large length scales. Therefore, we suggest that the destruction of large aggregates allows smaller structures to disperse more freely, resulting in a macroscopically more disperse arrangement inducing the non-universal conductivity exponent. Yet, ultimately, we cannot entirely exclude the possibility of the primary CB being impacted by ball milling as well. The combination of analytical centrifugation results, scattering experiments, and electrical measurements thus experimentally confirms the trend observed in the simulations: large CB aggregates significantly increase the electrical conductivity and reduce the percolation threshold. Thus, in order to construct highly conductive composite materials with minimal addition of the filler material, the CB should be prepared accordingly. Moreover, as previous experimental studies already suggested,⁵⁴ a low fractal dimension of aggregates is desirable because it allows us to grow larger aggregates with less material.

We have studied the percolation of fractal aggregates by means of Monte Carlo simulation and experiments. The fractal aggregates were modeled based on experimental data from analytical centrifugation, porosimetry, and x-ray scattering of a specific type of carbon black that is commonly used as a conductive filler. The physical aggregates have broad size distributions that depend on their processing history. Our simulations show that their percolation thresholds coincide if they are expressed as a function of the ensemble averaged radius of gyration of the aggregates. Aggregates with larger average sizes reduced the percolation threshold; higher order moments of the size distribution did not affect the percolation threshold. A similar observation has been made in a study on fracture networks⁴³ in which percolation of two-dimensional polygons in three-dimensional space was simulated. We tested this prediction experimentally by decreasing the aggregates’ average sizes through ball milling. This decrease increased the percolation threshold in a way that was consistent with the predictions.

We found that the connectivity range (i.e., the distance over which charges can be transported between aggregates) can also be incorporated in a quasi-universal scaling law as well as the interaction strength between the aggregates. Furthermore, we noted that the backbones of percolating networks formed by fractal aggregates consist of many conducting pathways, i.e., there is not just one “red link” that carries the entire current. Hence, in contrast to networks formed by rod-like particles, we expect networks of fractal aggregates to be resilient against the fracture of individual bonds. In addition, the fractal aggregates also have lower percolation thresholds than rod-like particles, which suggests that they are interesting materials as fillers in conductive composites.

Synchrotron ultra-small angle x-ray scattering experiments indicated that the carbon black forms hierarchical structures when added to a silicone polymer matrix. We modeled such agglomeration processes by adding short-ranged attractive interactions to the aggregates and found that in an ideal case, short-ranged attraction reduces the percolation threshold uniformly. The formation of agglomerates in the actual material is affected by the polymer matrix and process details such as the mixing technique, too. Notwithstanding such limitations, we find that our prediction that smaller average aggregate sizes generally result in larger percolation thresholds holds experimentally, too, even when agglomeration was detected. This shows that systematic studies comparing experiment and simulation can be done even for a filler system as complex as carbon black.

Materials: Carbon black {carbon black, acetylene, 100% compressed, 99.9+%, bulk density 170–230 g/l, Alfa Aesar [Chemical Abstracts Service Registry Number (CAS): 1333-86-4]}, polydimethylsiloxane (SYLGARD^TM 184 Silicone Elastomer Kit), and toluene [toluene anhydrous, 99.8%, Sigma-Aldrich (CAS: 108-88-3)] were used as received without further purification.

Ball milling of carbon blacks: We performed the ball milling experiments to investigate the relationship between filler size and composites’ conductivity. For all experiments, 1.5 g carbon black powders were placed in a 50 ml zirconia grinding jar with 30 ml isopropanol, and ball milling was performed at 400 rpm for 2 h, with reversal of direction each 10 min. 60 and 120 g milling balls (1 nm in diameter) were used, respectively, to reach two different size distributions. The ball-milled carbon blacks were washed with isopropanol and centrifuged and then dried at 100 °C overnight. The dried carbon blacks were hand-milled again and collected for further use.

Size distribution characterization: Carbon black diameters were analyzed by analytical centrifugation (AC) using a LUMiSizer 651 (LUM GmbH, Berlin, Germany). The carbon black/toluene suspensions of 0.0014 wt. % were prepared using 30 min ultrasound in an ultrasonic bath (Elmasonic S 100H). An ultralow concentration was used here, considering the transmittance and the stability of the suspension. The temperature was set to 20 °C for all tests. For a typical AC test, the blue wavelength of 470 nm was used. 2 mm optical path length polyamide (PA) cells filled with 430 μl suspensions were used for all experiments. Centrifugation tests were performed immediately after the cells were filled. A volume-weighted particle size distribution analysis was performed after the measurement.

Composite preparation: The composites were prepared by carbon black fillers and PDMS in different volume ratios in a Speedmixer (DAC 600.2 VAC-P Vacuum Mixer System). First, carbon blacks and PDMS base were mixed at 2350 rpm for 3 min under vacuum. Second, the curing agent was added to the PDMS-CB mixture with 1/10 of the PDMS mass, and the mixture was then remixed under the same condition. The mixture was cured at 100 °C for 2 h. After cooling down, the conductivity of the samples was measured by a four-point probe method triggered by a Keithley 2450 source meter at room temperature. The CB density was determined as 1.7 g/cm³ using a helium Pycnometer AccuPyc 1330. The density of PDMS was assumed to be identical to that of the liquid 0.965 g/cm³. All volume fractions were calculated using these values.^55,56

Ultra-small angle x-ray scattering: The USAXS experiment in transmission geometry was carried out at the P03/MiNaXS beamline at the synchrotron source PETRA III at Deutsches Elektronen-Synchrotron (DESY), Hamburg, Germany.⁵⁷ The experiments were performed during two campaigns: (1) x-ray beam with the energy 12.85 keV (wavelength λ = 0.0965 nm, Δλ/λ = 10⁻⁴) and (2) 11.8 keV with a wavelength of λ = 0.105 nm. In both cases, the x-ray beam was focused on the detector by beryllium compound refractive lenses (CRL) with the Sample–Detector Distance SDD = 9760 and 9348 mm, respectively. Transmission was calculated by a diode inside the beamstop. To increase resolution, the beamstop was deliberately shifted away from the transmitted x-ray direct beam.⁵⁸ This allowed for reaching ultra-small scattering angles. A PILATUS 300 K detector (Dectris Ltd., Switzerland) with the pixel size of 172 × 172 μm² and dimensions of 487 × 619 square pixels (83.8 × 106.5 mm²) was used as an area detector. Azimuthally averaged radial distribution of intensity was performed after the analysis of USAXS data using the DPDAK software package.⁵⁹ 

Other characterizations: The TEM image was recorded using a JEM 2010 (JEOL, Germany), working at 200 kV. The BET measurement was carried out with Quadrasorb EVO, measured at 77.35 K. The sample was filled in a sample cell and put in the degasser. Under vacuum, it was heated up to 130 °C and degassed for 10 h. After cooling, the sample weight was determined, the filler rod was inserted into the sample cell, and then the sample cell was set in the measuring instrument.

We acknowledge funding from the German Research Foundation under Project No. 404913146. L.G.G. acknowledges the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. 949785 ELECTROFLUID).

There are no conflicts to declare.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

A standard Metropolis Monte Carlo algorithm was used sampling rotational as well as translational degrees of freedoms of the aggregates. As long as aggregates only interact hard, every allowed configuration, i.e., configuration without overlap of hard structures, is an equilibrium state. Due to the built-in randomness of the aggregates, however, there is no straightforward systematic way to set up the system in an allowed state. Therefore, we setup an aggregate ensemble by positioning the center of mass of each aggregate on a vertex of a cubic lattice with random orientation. Instead of the hard potential, we apply a potential that penalizes overlap between different aggregates based on the extent of mutual permeation. Then, we slowly reduce the system temperature until all overlaps are eliminated so that we can switch to hard interaction. To make sure that the system evolves away from the initial allowed state, we observe the diffusive behavior of the mean square displacement. For systems with a supplemental potential, we apply the same technique to obtain an allowed state that is subsequently evolved until the system energy fluctuates around a stable minimum. We then apply the same sanity checks to check for ergodicity.

One primary sphere is placed at the origin as the aggregate seed before other primary particles randomly diffuse one by one until they hit a pre-existing structure. Those singular primary particles are randomly displaced with a maximum step length of 0.1d so that the maximum overlap of primaries in the final aggregate is controlled. The formed aggregate is stationary, i.e., its rotational diffusion is neglected. This simplification physically encodes a time-scale separation between the movement of the formed aggregate and primary particle but is primarily chosen to keep the model as simple and reproducible as possible. During the aggregation process, the box size changes dynamically with the size of the formed aggregate. After a primary sphere is attached to the designated aggregate, its center of mass is shifted to the origin of the box. Adding three primary diameters to the maximum norm of a constituent primary particle of the aggregate becomes half the new box length. The next mobile particle commences diffusion from a surface of a sphere with the diameter of the new box length. Naturally, the primary spheres forming the final aggregate permeate each other so that their combined volume is less than the sum of two sphere volumina. The actual hard volume of each aggregate can be computed analytically as all overlap volumes can be decomposed into pairwise sphere–sphere intersections. For simplicity, however, we instead obtain an accurate approximation to the hard volume by brute force Monte Carlo integration. We vary the volume fraction exclusively via the box length so that the hard particle content never changes. For each volume fraction simulated, we probe the box at 10⁷ independent spots and calculate the corresponding estimate for the hard volume. Finally, we take the average over all box lengths. It is this hard volume that we report as reference for our simulation results.

Using periodic boundary conditions, a cluster of aggregates is percolating if a member of that cluster is connected to one of its mirror images. In turn, two aggregates are connected if the connectivity shells of any of their corresponding primary particles overlap. In order to determine the percolation threshold, we measure the percolation probability Θ(η) as a function of the volume fraction η, i.e., for a given volume fraction, we take 10 000 decorrelated snapshots of the system and probe for the existence of a percolating cluster. The ratio of percolating snapshots is the simulational estimate for the percolation probability. In an infinite system, Θ(η) would be a step function with the step demarcating the percolation threshold p_c. Periodic boundary conditions, however, induce finite size effects, which give Θ(η) a sigmoidal shape. Following a common protocol,⁶⁰ we fit the percolation probability to a hyperbolic tangent of the form

where a is the parameter characterizing the width of the percolation probability.

The simulation results presented above were obtained using cubic boxes of finite size with periodic boundary conditions. As the correlation length diverges in proximity to a critical point, finite size effects have to be taken into account. Thus, this section discusses how the observed percolation thresholds depend on the length L of the encompassing box.

Since we always determined the percolation threshold of an aggregate ensemble by varying the box size, we need to simulate aggregate ensembles of different sizes to systematically probe for finite size effects.

Figure 13 displays the percolation probability in the vicinity of the percolation threshold for equivalent aggregate ensembles of varying sizes. Clearly, even for the largest system sampled corresponding to L > 140d, the percolation probability still differs significantly from the step function expecting an infinite system. As expected, the smaller the system, the broader the transition. The width of the percolation probability can be used to determine the critical exponent ν governing the divergence of the correlation length. The fitting parameter a from the sigmoidal fits to the percolation probability like a⁻¹ ∝ L^−1/ν. Treating ν as a fit parameter, we recover the expected universal exponent of three-dimensional percolation ν ≈ 0.8774.⁶¹ Recall that we chose the volume fraction at percolation probability 0.5 as percolation threshold. Figure 13 shows that this value substantially drifts to smaller values as the system size is increased. Percolation theory quantifies this drift as

with p_eff being the observed percolation threshold, p_∞ being the percolation threshold of the infinite system, and ν being the critical exponent. Having previously verified that ν takes the expected value, we can determine the pre-factor of the power law Eq. (A2) as well as p_∞ by a linear fit to the simulation results for p_eff. Note that the values for L that we use in that fit are the box lengths at percolation probability 0.5. As the box lengths are varied to adjust the volume fraction while the number of aggregates N_agg is kept constant, one might consider N_agg the more natural scaling variable. As N_agg ∝ L³ at constant volume fraction, we can equivalently use $peff−p∞∝Nagg−1/(3ν)$ to obtain the same estimate for p_∞.

Figure 14 illustrates that the percolation threshold as initially evaluated overestimates the percolation threshold of the infinite system by roughly 2.5%, which is less than the marker size in the corresponding figures in the main text. The pre-factor according to the fit above is roughly 0.31; however, this pre-factor may vary strongly with the aggregate ensemble. A rigorous finite-size scaling for every ensemble we analyzed would require disproportionate effort. Nevertheless, we repeated the procedure for aggregates with 50 and 200 primary particles, respectively, and deviation due to finite size effects is almost the same although the pre-factor grows slightly with the number of primary particles. For the remaining ensembles, the width of the percolation probability curve as shown in Fig. 13 is a good indicator of the strength of finite size effects. For all percolation thresholds reported above, the system size was chosen so that the range between percolation probabilities 0.2 and 0.8, respectively, was smaller than 10% of the volume fraction at probability 0.5. Thus, we expect the percolation thresholds displayed in the main text to be collectively overestimated by roughly 2.5% without altering the relationships found and the conclusions drawn above.

We analyze the simulated aggregate networks by assuming aggregates to be perfectly conductive rendering tunneling contacts the only source of resistance. The tunneling conductance is set to unity for hard contact between aggregates and decays exponentially to 0.01 for a surface separation of 2ϵ, i.e., when the connectivity shells just barely still overlap. Beyond a surface separation of 2ϵ, there is no tunneling contact. We can now vary the microscopic conductance distribution by adjusting ϵ. We pick two opposite faces of the simulation box for which we suppress periodic boundary conditions. Instead, we consider each aggregate permeating one of the faces as connected to a plate in a parallel-plate capacitor setup. Finally, we apply a unit current between the plates and solve the ensuing random resistor network for the resistance between the parallel plates. We find that simulations with ϵ = 0.2 are consistent with the universal exponent t ≈ 2.0 irrespective of the aggregate size distribution. Conversely, simulations with ϵ = 0.5 are not consistent with the universal exponent but can be adequately fitted to a power law with a larger exponent. Figures 15 and 16 show how the conductivity changes with the volume fraction. Notice that each data point corresponds to an average over 400 networks generated from decorrelated equilibrium configurations. The error bars mark the standard deviation of each set of conductivities. Non-percolating networks are treated as non-conductive.