Numerical Simulation of the Percolation Threshold in Non-Overlapping Ellipsoid Composites: Toward Bottom-Up Approach for Carbon Based Electromagnetic Components Realization

: A Monte Carlo (MC) model for the calculation of the percolation threshold in the composite filled with ellipsoids of revolution is developed to simulate the real experimental situation of percolative composites in which functional additives do not penetrate each other. The important advantage is that the MC model can be easily applied to multi-components composites, e.g., containing graphene nanoplatelets, carbon black and carbon nanotubes, by means of utilising the ellipsoids of different aspect ratios with the filling fraction corresponding to concentrations of each type of inclusion. The developed model could be used in a pre-experimental step for producing effective close-to percolation and percolated nanocomposites for various electromagnetic applications to avoid time and resources consuming the “sort-out” experimental phase of composition optimization, and could be utilized as the first step of the bottom-up material approach to touch the macroscopic platform for antennas/circuit realization.


Introduction
The modeling of nanoantennas and nanocircuits depending on the material concept, i.e., a nanosized object as a working element of nano-device or a composite made of nano-inclusions as a macroscopic platform for antennas/circuit realization, could be provided in "top-down" and "bottom-up" approaches. The first one accesses the device realization by computation of microscopic electromagnetic response of nanoparticle (e.g., carbon nanotube) [1][2][3], whereas the second one is expected to describe the collective effects of nanoparticles electromagnetics in the final composition caused by individual nanoparticles antennas' properties through different homogenization and/or averaging techniques [4][5][6][7].
This communication is one of the second sort. It is motivated by an exceptional interest in nanocarbon based composites close or slightly above percolation concentration as the material support for design light and ultrathin effective electric components, optical and optoelectronic devices, constituent elements for electromagnetic compatibility realization [6,8] Moreover, polymer nanocomposites close to or slightly above the electrical percolation threshold are of tremendous interest as a material-science basis for conformal structurally integrated microwave-to THz antennas [1,9].
The percolation threshold for multi-facial irregular media is a problem that was studied for a long time (see [10,11] and Refs. therein). In modern material science, the main focus is the percolation phenomena in composite with randomly distributed carbon nanoinclusions like graphite nanoplatelets, carbon nanotubes, carbon black, onion-like carbon, etc. [12]. There are many approaches to simulate the nanocarbon composite system. For the case of nanotubes or fibers, the most obvious is to simulate it as a caped [13,14] or non-caped [15,16] cylinder. However, this approach is not universal because it does not allow to simulate all possible geometries of nanoinclusions.
All mentioned objects may be roughly modeled as ellipsoids of revolution with a different axial ratio. Many papers contribute to the modeling of the percolation threshold of overlapping ellipsoids [17][18][19]. The method of overlapping inclusions is faster; it describes the transport properties in porous media well [20], as well as the social phenomena [21], etc. However, it is not accurate in the case of nanocarbon composites, because these particles are robust and cannot overlap. Very few papers contribute to the modeling of the percolation threshold for the suspension of the non-overlapping ellipsoids [22,23].
We report the Monte Carlo model for the calculation of the percolation of the ellipsoids distributed inside a cubic unit cell. The dependence of the percolation threshold on the ellipsoids' distribution, unit cell size and the aspect ratio for the case of oblate ellipsoids will be discussed.

Modelling
The model generates the system of non-overlapping ellipsoids of revolution located inside the cubic unit cell. The dimensions UC = nd, where d is the maximal diameter of the ellipsoid and n > 1.

Positioning the Ellipsoid in 3D Space
The ellipsoid with semiaxes b i , i = 1, 3 is a set of points X, that satisfy the condition: where X 0 is vector of ellipsoid center, Q 0 is positively defined matrix. Any Q can be converted to diagonal matrix A as follows: Here A is diagonal matrix with elements a ii and R(φ, θ, Φ) is Euler's rotation matrix. The rotation matrix can be generally evaluated as a series of rotations: Let b i be the semiaxes of the ellipsoid then, elements of A may be computed as a ii = 1/b 2 i . If b 1 = b 2 (for the ellipsoid of revolution), the rotation R z (Φ) turns to unitary. So, the total number of the coordinates is 5: K = (X 0 , φ, θ). Using the described method, we can easily construct the ellipsoid of revolution with known orientation and position as a combination of the diagonal matrix A and 5-component vector K.

Distance between Ellipsoids
Minimal distance ∆D min between two ellipsoids (A, K) and (A , K ), is a solution of the minimisation problem: where S and S are the surfaces of ellipsoids. The equation 4 cannot be solved analyticaly, but there is a variety of numerical methods. All methods may be divided into geometrical [24,25] and algebraical [26,27].
Here we will use the geometrical approach developed by Lin and Han [25]. The idea of the method is to construct two balls completely inside each ellipsoid (see Figure 1).
can be used.

Composite Generation Procedure
The procedure of the composite creation is follows. On i-th step, we generate random K as: Next, two conditions should be satisfied: i i-th the ellipsoid should not intersect the walls of the unit cell. ii i-th ellipsoid should not penetrate into any ellipsoid of the already existing system of i−1.
For the real nanocarbon composite (ii), condition is more strict. Taking into account the Van der Waals separation, the minimal separation between two ellipsoids here is considered as 0.34 nm. If both conditions are satisfied then we store i-th ellipsoid, if not-we abandon it and make a new attempt for i-th (see Figure 2).

Percolation Computation
Once the system is computed, well-known Dijkstra's [28] algorithm may be applied for the calculation of the percolation path. The tunnelling distance of 2 nm was used as a connection criterion. Next, for Dijkstra protocol, we have to select an initial and final point (edge) of the graph. Obviously, these points are the ellipsoids located near i-th border of the unit cell and the border opposite to i-th correspondingly. Periodic boundary conditions are utilised. So the distance between the ellipsoid near the i-th boarder, and the opposite one shifted along the chosen direction on the unit cell, should be less than the predefined separation of 2 nm. If this condition is satisfied, the ellipsoids are considered as the initial and final points of the graph.
Finally, the boolean vector of the percolation may be computed. We consider the system percolated if it is percolated in at least one direction.
Generation of the composite of N(p) ellipsoids, where p is the volume fraction. A result of this step is an array of ellipsoid coordinates K. The number of strings is N(p).

2.
Checking if the system is percolated.
If it is, we then terminate, if not, we go back to (1) with p := p + δp and construct new composite. The implementation of the algorithm is written with Fortran (see Figure 3).

Two-Phase system
For the ellipsoid of revolution b 1 = b 2 , then the aspect ratio can be introduced as ρ = b 1 /b 3 , ρ ≥ 1. In this paper we consider b 3 = 2.38 nm, so the total thickness corresponds to 14 graphite interlayer distances. The formation of the percolation paths depends on the concentration of filler, and also the distribution of the inclusions. The distribution effect as studied for the system filled with ellipsoids of ρ = 5. We made 10,000 observations for n = 3, 7000 for n = 4 and 3000 for n = 5. The percolation probability function is presented in Figure 4.  The empiric probability functions are governed by the Weibull probability function, W(λ, k) = 1 − e −(x/λ) k [32]. The parameters of equation are follows: λ = 0.115, k = 8.797 for n = 3; λ = 0.113, k = 12.092 for n = 4, and λ = 0.113, k = 14.90 for n = 5. The Weibull function describes "time-to-failure" and it is typical for percolation processes [33][34][35]. All the probability functions for different cell size are closely collapsed, but the diffusion is slightly decreased from n = 3 to n = 5, that corresponds to the rise of k. This means that the percolation threshold is independent of the UC. The diffusion difference may be explained by the finite volume fraction of the single ellipsoid inside the unit cell. On the other hand, the calculation time increases drastically with UC rise, so the UC size of n = 3 is optimal for the calculation of the composites with ellipsoids of higher aspect ratio. Figure 5 shows the dependence of the percolation threshold on the ellipsoids aspect ratio. As expected, the percolation threshold decreases with the aspect ratio rise. The resulting dependence is in good agreement with experimental data [36][37][38] and theoretical modeling [17,39].

Three-phase system (Hybrids)
The calculation of the percolation of hybrid composite requires the modification of composite generation procedure Section 2.3. The generation is organised "one in time". We introduce the parameter prob-the probability of the GNP appearance, then 1-prob is the probability of the CNT appearance. So, upon the generation procedure, each i-th particle has the probability to be "born" as GNP or CNT. After, the percolation was computed using the standard protocol. For the hybrid composite, we consider p c as a total volume fraction of the inclusions required for the percolation. Obviously, for the hybrids total concentration may be presented as a sum of partial concentrations as p tot c = p GNP c + p CNT c . The dependence of the p tot c (prob) is presented on Figure 6. For the case of prob = 1 and 0, the data was verified with the model for 1 type of inclusions. For the calculations two types of ellipsoids were used: The excluded volume theory states that the percolation for the hybrid composites [40] occurs if the following condition is satisfied: where p GNP c and p CNT c are partial concentrations of the GNP and CNT in hybrid composite, and p tot c (1) and p tot c (1) are the critical concentrations for composites with GNP (prob = 1) and CNT (prob = 0). However, a lot of experimental works evident the percolation even if the condition Equation (6) is not satisfied [41][42][43]. In our case, the exact number of GNP and CNT is known for each configuration, so we can easily calculate the EVT-parameter. The result is presented in Figure 7. The mean value of EVT is below 1 for all studied configurations. The dependence demonstrates minimum at prob = 0.7-0.8. We can preliminarily conclude that there is some optimum combination of CNT and GNP in hybrid composites.

Conclusions
This study is the next computational physics step in the series of our papers [44][45][46] to approach electromagnetic components such as antennas, filters, polarizers, shields by using conductive/dielectric composite materials based on carbon additives. We present a Monte Carlo model for the calculation of the percolation threshold in the irregular system filled with ellipsoids of revolution. We demonstrate that the percolation concentration is independent of the unit cell size. The percolation concentration dependence on the aspect ratio of nanoinclusions is in the good agreement with that previously reported [17,22,23,39].
The reported model can be applied for nanocarbon containing composites as it supports the real experimental situation when percolated additives do not penetrate each other. The important advantage of the present model is that it can be easily applied for the modelling of multi-components composites, e.g., containing graphene nanoplatelets, carbon black and carbon nanotubes, by means of utilising the ellipsoids of different aspect ratio with the filling fraction corresponding to concentrations of each type of inclusion. This tested in simple case of oneand two-component composite (two and three-phase systems respectively) MC model could be used in a pre-experimental step for producing effective close-to percolation and percolated nanocomposites for various mechanical, thermal and electromagnetic applications to avoid time-and resources consuming "sort-out" experimental phase, that should lead to the design of optimal composition of many functional components providing e.g., the lowest overall percolative concentration at the lowest content of the most expensive functional additive, retaining at the same time high performance and functionality.
Author Contributions: A.P. and P.L. designed the computational experiments, A.P. developed the M.C. model, J.M. and G.S. tested the model vs. real experimental situation, A.P. and P.K. conceived the task as well as methodolgy and wrote the manuscript.
Funding: This research was funded by H2020 Marie Skłodowska-Curie Actions, grant number 734164 Graphene 3D, and Tomsk State University Competitiveness Improvement Program.