Understanding the anomalous frequency responses of composite materials using very large random resistor-capacitor networks

In this paper large resistor-capacitor (RC) networks that consist of randomly distributed conductive and capacitive elements which are much larger than those previously explored are studied using an efficient algorithm. We investigate the emergent power-law scaling of the conductance and the percolation and saturation limits of the networks at the high and low frequency bounds in order to compare with a modification of the classical Effective Medium Approximation (EMA) that enables its extension to finite network sizes. It is shown that the new formula provides a simple analytical description of the network response that accurately predicts the effects of finite network size and composition and it agrees well with the new numerical calculations on large networks and is a significant improvement on earlier EMA formulae. Avenues for future improvement and explanation of the formula are highlighted. Finally, the statistical variation of network conductivity with network size is observed and explained. This work provides a deeper insight into the response of large resistor-capacitor networks to understand the AC electrical properties, size effects, composition effects and statistical variation of properties of a range of heterogeneous materials and composite systems.


Introduction
Composite materials that consist of conductive and dielectric phases are present in a wide variety of manufactured materials and in nature. Therefore an understanding of the behaviour of such materials in terms of composition effects, size effects and statistical variation of properties is critical to enable the design of composites with specific properties. An intriguing aspect of heterogeneous materials and composite systems is the anomalous power law frequency dependence of the AC conductivity; termed by Jonscher as the "universal dielectric response" (UDR) [1][2][3]. A variety of systems have been observed to follow the universal dielectric response; examples include ionically conducting ceramics in single crystal and polycrystalline form, glassy materials, mixtures of rock and water, polymer composites loaded with a variety of conductive fillers, and biological materials such as skin and nail.
A characteristic of the universal dielectric response is that at low frequencies the AC conductivity, Y (ω) is frequency independent [4], Y (0) which can be considered the DC conductivity, while at higher frequencies the AC cona e-mail: c.r.bowen@bath.ac.uk ductivity increases with increasing frequency and follows a power law, where ω is angular frequency (2πf ), A is a constant and 0 < n < 1.0. The relative permittivity (k) also follows a power-law decay such that, k ∼ ω n−1 . Figure 1a shows an example of the UDR from experimental data of a water saturated porous ceramic with 22%vol. porosity whereby measurements were made at regular intervals until the ceramic was dry [5]. The frequency independent region, Y (0), and frequency dependent regions can be clearly observed and Y (0) decreases as the ceramic dries and the fraction of conducting water decreases. A recent interpretation of the UDR is that it is a characteristic of the random resistive and capacitive paths formed by the microstructure of the material. These flowpaths can be modelled using 2D lattices of randomly positioned resistors and capacitors; for example in the water saturated ceramic in Figure 1a the resistors represent conductive water filled pores and the capacitors represent the regions insulating ceramic. Such models do not directly represent the 3D physical structure of the disordered binary mixture, but are rather a more abstract representation of the various flow-paths though the modelled Example of a 2D resistor-capacitor (R-C) network to represent a microstructure containing randomly distributed conductors and insulators. (c) Simulations of conductivity and capacitance of a relatively small 2D square RC network containing 512 randomly spaced components, 60% 1 kΩ resistors and 40% 1 nF capacitors. R −1 = ωC is satisfied at ω = 10 6 s −1 .
medium. The resistors have a fixed conductivity, strictly the complex admittance, and map to conductive flowpaths in the material, whereas the admittance of the capacitors varies according to the applied frequency and represent paths through the dielectric regions of the medium.
An example of an RC network is shown in Figure 1b which consists of a two-dimensional array of randomly positioned resistors and capacitors. The analysis of the properties of random resistor-capacitor networks has been presented by Almond and co-workers [4][5][6][7][8]; typically for networks that contain a much smaller number of resistors or capacitors (less than 10 4 components) than will be examined in this work. Good correlation of experimental data with the RC microstructural network approach has been reported using a water infiltrated porous ferroelectric ceramic which can be considered as a random microstructure consisting of conducting water filled pores that represent the resistive regions (R in Fig. 1b) and insulating ceramic that represent the capacitive regions (C in Fig. 1b) [9].
An example of a typical response of a RC network is presented in Figure 1c, which shows the frequency dependent conductivity and capacitance of a small network of 512 components where 60% of the components are 1 kΩ resistors and 40% are 1 nF capacitors; as determined using circuit simulation software. The conductivity and capacitance in Figure 1b reveals that the frequency response is qualitatively similar to the UDR; i.e. equations (1) and (2). The network response is best understood by considering the frequency dependent conductivity of the individual resistor and capacitor components in the network. The resistor conductance (R −1 ) is frequency independent while the capacitors capacitor conductance (iωC) rises linearly with frequency. This leads to the following observations [4,7,9]: (i) at low frequencies the conductance of the resistors is much larger than that of the capacitors which effectively act as open circuits, R −1 > ωC. The AC currents will flow through the resistors and if a percolation path of resistors exists across the network the current will flow preferably through these resistors; (ii) at high frequencies the conductance of the capacitors is much larger than that of the resistors, R −1 < ωC. The capacitors effectively act as short circuits and AC currents will flow through the capacitors. If a percolation path of capacitors exists across the network the AC current will flow preferably through these capacitors; (iii) at "intermediate" frequencies the AC currents flow through both resistors and capacitors, since R −1 ∼ ωC. In this frequency range the presence of percolated resistors or capacitors has little influence on network response since both resistors and capacitors have a similar conductance.
The power law response of the electrical networks, such as that observed in Figure 1c, occurs at intermediate frequencies where the resistor and capacitor conductances are similar (R −1 ∼ ωC). For the network in Figure 1b the condition R −1 = ωC is at ω = 10 6 s −1 , since R = 1 kΩ and C = 1 nF. Almond et al. have shown empirically that the network power law response [4,9] can be related to the component values by a simple logarithmic mixing rule where the network complex conductivity (Y * network ) is where p is the fraction of capacitors in the network. Truong et al. [10] have presented a possible explanation of the logarithmic mixing rule for random mixtures as a general case of the series and parallel mixing rules. The logarithmic mixing rule also predicts that the phase angle in radians, is The frequency response of resistor-capacitor networks has been employed to understand the frequency dependent conductivity and permittivity of a number of materials and composites at a range of scales of conductorinsulator heterogeneity. This includes zirconia ceramics [11], ferroelectric-conducting polymer composite systems [7], granular thin films [12], and water adsorbed in mesoporous materials [13]. The networks have also been used to explain the influence of conductivity on the high permittivity of lead halide perovskite solar cell materials [14], understand electromagnetic moisture sensing of materials [15] and dielectric loss in ferrite materials [16]. Two-dimensional arrays of resistors, capacitors and diodes have been used by Bychanok et al. to examine the electromagnetic properties of polymers filled with carbon nanotubes [17]. Mechanical networks that consist of elements of different stiffness have also been considered [18]. Previous analysis of relatively small RC networks [4][5][6][7]11,19] has shown that the AC conductivity exhibits a large distribution of magnitudes at low (R −1 > ωC) and high (R −1 < ωC) frequencies. These large distributions are due to the different possible combinations of percolating or non-percolating resistors or capacitors that form in the network. It is clear that whether or not percolation paths are formed is dependent on the proportion of resistors or capacitors in the network. The power law region at intermediate frequencies is far more sharply defined for range for networks of the same composition since the AC current is flowing though both components in this frequency range; it has been suggested [4] that the power law response in this frequency range is an emergent property of RC networks. To date simulations have been undertaken on relatively small random RC networks with a small number of random networks realizations. This has limited the possible insights on aspects of the influence of RC network size on the low and high frequency distribution of AC conductivity, the statistical variation of conductivity with network size and finite size effects on percolation. As an example, Almond et al. [8] examined networks where the number of components (N ) was approaching 10 000 and Hamou et al. provided a detailed study on large numbers of networks with up to 500 000 components [20]. The influence of network size is of interest to understand the influence of scale in real materials or composite systems; for example bulk materials can be effectively considered as infinitely large networks while composite thin films at the micrometer scale containing nanometer sized inclusions are effectively smaller finite electrical networks.
In this paper we explore the validity of theoretical approximations in various regimes (network size and composition) by comparing them to large ensembles of 2D square electrical networks at a very large range of sizes. In Section 2 we summarise the various analytical and numerical methods from the previous work of ourselves and others. The analytical expressions presented in Section 2.1 combine classical averaging approaches with spectral methods into a suggestive asymptotic formula. This formula enables extension of the Effective Medium Approximation (EMA) to finite network sizes, and avenues for future improvement and explanation of the formula are highlighted. The RC network models are solved numerically using an efficient algorithm based on the reduction method of Frank and Lobb [21], which is outlined in Section 2.2. The approach allows numerical calculations of the electrical properties of RC networks with network sizes and random realizations that are much larger than previously explored, approaching up to 2 × 10 6 components with 1024 random realization of each network type. The algorithm is used to calculate the equivalent properties of networks of a range of sizes and proportions of resistors and capacitors. This enables investigation of the validity of these approximation schemes to determine the network properties, as presented in Section 3. In Section 3.1 the typical responses of the RC networks in different regimes of composition (p) and size (N ) are demonstrated, before comparison with the asymptotic formula in Sections 3.2-3.4. Finally, in Section 3.5 some interesting observations of the statistical variation of network conductivity with network size at the different ranges of frequencies, R −1 < ωC, R −1 > ωC, and R −1 ∼ ωC are explored and explained.

Resistor-capacitor (RC) network calculation methodologies 2.1 Analytical approximation methods
For networks of randomly sited single conductive components, classical percolation theory gives good results [22]. However, in the case of binary mixtures where the components have a variable conductivity ratio (here due to the capacitive components having frequency-dependent admittance), other analytical approaches are required. These include: empirically derived mixing laws, such as equation (2), which are valid in the intermediate powerlaw region [4]; averaging approaches such as the Effective Medium Approximation [23], valid away from the critical percolation probability where network size (which represents granularity in a real material) does not have an effect; approaches based on the distribution of the polezero (generalised eigenvalue) spectra [8], which account for network size effects near the critical mixing ratio but only approximately capture the percolation and satura-tion limits away from criticality; as well as more speculative combinations thereof, as described next.

Effective medium approximation (EMA) for an infinite network
The purely empirical power-law mixing rules described earlier, equations (2) and (3), are limited to providing an estimate of the properties of RC networks at intermediate frequencies where R −1 ∼ ωC. However, the approach does not provide an indication of the AC conductivity at low or high frequencies where currents percolate through resistors or capacitors, respectively [4]. Here, the "classical" Effective Medium Approximation (EMA) formula is derived using an averaging method, which provides an approximation to the conductance as a function of ω for a homogeneous RC network. The approximation assumes perfect mixing, implying an infinitely large network, and therefore works only in cases where network size is unimportant, away from the critical mixing ratio of the two components. Empirical corrections have previously been made in cases of finite-size networks [24], but these are far from adequate, as will be demonstrated later (Sect. 3.3) when we examine a wide range of RC network sizes.
We initially consider resistors with a frequency independent conductance y 1 = 1/R and capacitors with a frequency dependent conductance y 2 = iωC that are in proportion 1−p and p, respectively. The EMA states that the effective medium conductance Y for a two-dimensional square lattice where the number of components N → ∞ is given by solutions of the quadratic equation, for which solutions can be obtained either numerically, or in the usual way. This formula can further be rearranged, by defining into the following symmetric form, which will be compared to the spectral formula in the next section, As stated above, the EMA assumes perfect mixing in an infinitely large network and it fails to capture any size, N , dependent scaling of the network response. We will see later that as the fractions of capacitors or resistors in the network approaches the percolation limit when |ε| is small, which for a two-dimensional square lattice p c = 0.5 the network response becomes more dependent on network size [8]. Therefore, the EMA is most appropriate when p is not too close to the network percolation limit. We can draw a number of conclusions from equation (6).
Firstly, if ε = 0 then we have that which is the well known Keller mixing law [25].
Secondly, if ε is positive, and |μ| is large, then there is an asymptotic solution of the form Thirdly if |μ| is small then Note that the ratio of these two solutions is given by 1 ε 2 . If ε is negative then the roles of y 1 and y 2 are reversed. It is of interest to observe that the first of these formulae becomes unbounded as we approach the critical value of ε = 0. This is because we have made the assumption that the networks we are considering have an infinite size.

Spectral formulae for a finite network
A complementary method, using a combination of statistical techniques and spectral analysis, has been used by Budd et al. [8] to derive analytical expressions that are valid for critical values of p c = 0.5 in the case of a finite network size. These methods accurately capture the effect of the network size on the conductivity frequency response. This approach works by considering the statistics of the poles and zeros of the rational function of μ, which describes the admittance frequency response of the network. For RC networks these poles and zeros all lie on the imaginary axis, and interlace with each other. Indeed, the poles and zeros correspond to the eigenvalues of the Laplacian matrix of the RC network formed by applying Kirchoff's laws to it, with extensions to include the boundary conditions which apply at the edges of a finite network. The method uses statistical arguments to show that these eigenvalues satisfy a log-Normal distribution and determines certain statistical regularities between their spacings, which apply in the case of ε = 0. By applying these statistical regularities to an evaluation of the admittance function, it is shown in [8] that it is possible to derive four formulae for the admittance function for each of the four combinations of percolating or non-percolating resistor or capacitor bonds, which is strictly valid in the critical case of p = p c = 0.5.
In the following treatment we illustrate only the case with percolating resistor components (formula 40 in [8]), which has a response resembling networks with p < 0.5, ε > 0, but the approach can be used for the other three cases starting from the results in [8].
Restating the working in [8] for the case considered with percolation of resistors, but not capacitors, if the network is a square lattice with a large number of components, N , then the admittance We can rearrange the above equation following [8], in a similar way to the classical EMA and using the same definitions, into the symmetric form, equation (8), This equation has a very similar form to equation (6), with the term of size 1/N acting as a perturbation of the Keller mixing law. It is shown in [8] that equation (8) provides a close agreement to the calculated percolation/saturation limits for a finite network, as well as the scaling behaviour of the conductivity with network size for intermediate frequencies. To summarise this: Firstly, as N → ∞ we recover the exact Keller mixing law with θ = 1. Secondly, if N is large, and |μ| is also large, then equation (8) has the asymptotic solution describing a percolation situation which takes the form, Thirdly, if N is large and |μ| is small, then equation (8) has a further asymptotic solution describing a percolation situation which takes the form In this case the ratio of the two results is 1/N . Fourthly, the mixing region with θ ≈ 1 is valid for the range in which, and hence the width of the mixing region depends critically on the network size.
We shall see a numerical validation of all of these results in the next section on our simulations for RC networks for a wide range of sizes.

Combining the two formulae to give a Modified EMA (MEMA) formula valid for finite networks
We now make the assumption that there is a more general formula for the scaled admittance, θ, which depends both on the deviation from criticality ε and the deviation δ = 1/N from an infinite network size. The derivation of such a general formula is likely to be problematic, however we can deduce its approximate form in the case of small ε and δ from equations (6) and (8). Indeed, both of these can be considered as perturbations of the Keller mixing law θ − 1 θ = 0, as either small values of ε or δ vary from zero. Therefore, using the multivariate Taylor expansion, we can combine these two formulae to give the following more general approximation [8] that is valid over a range of values for which ε is small and N is large: We call this new formula the modified effective medium approximation, MEMA. We note that it can easily be rearranged into a simple quadratic equation and can thus be readily solved to find θ for all values of N and ε. For very large N it reduces to the EMA approximation (Eq. (6)) (and the size of the network becomes unimportant), and for very small ε to the spectral approximation (Eq. (8)), where the size of the network becomes critical. There is thus a point where one approximation becomes more valid than the other. To see this more precisely we consider the case of large |μ| in which case we have the approximate formula, rearrange and ignore terms of size 1/μ to give (to leading order) Hence we have, on taking the positive root of this equation, Hence on rescaling we recover the previous two percolation limits. We draw the important conclusion that the behaviour of the network is approximately independent of N when the network is sufficiently large such that [8], otherwise we will see the effects of the finite network size. The results of solving the MEMA analytical approximation (Eq. (9)) will be compared to the numerical solution of the wide range of network sizes solved using the Frank and Lobb approach in Section 3.4.

Numerical Frank and Lobb method
In this section the method employed to evaluate the frequency dependent admittance and phase of the large RC networks is described. We consider a square network with S nodes on the horizontal side and therefore a total of N = 2S 2 components; for the sake of simplicity the diagrams in Figures 3-5 show only resistors but we have generalised the approach to consider networks containing both resistors and capacitors. An example network with S = 4 and thus N = 32 components is shown in Figure 2a. The total admittance Y (ω) of this network is obtained through the equivalent impedance Z eq (ω) calculated using  the Frank and Lobb [20] reduction scheme. The method of Frank and Lobb is based on transformation of star-delta and delta-star, as shown in Figure 2b, and the transformation is defined in both directions.
This method applies a sequence of delta-star transformations to the components (the resistors or capacitors) that form the bonds between the nodes of the lattice. An example of this transformation, which we call the propagator transformation, is illustrated in Figures 3a-3d. This can be decomposed into three separate parts: first a deltastar (Δ → Y ) transformation (Figs. 3a and 3b), then a redefinition of the lattice point, and finally a star-delta (Y → Δ) transformation (Figs. 3c and 3d). The middle step is included primarily for clarity. In transforming from Figures 3a to 3b, admittances 6, 7, and 8 are determined from 1, 2, and 3 by equation (11): Then, in the transformation from Figures 3c to 3d, admittances 9, 10, and 11 are determined from 4, 5 and 8 by equation (12): When some of the admittances are missing, the propagator usually becomes more simple. An example of this is shown in Figures 3e-3h, where one bond is missing [20]. The final result of this sequence of transformations, as in Figure 4, is to reduce the lattice to a single bond that has the same admittance (Fig. 4i) as the entire lattice in Figure 4 [21]. The benefits of the approach compared to other approaches, such as efficient sparse matrix techniques, as used in the simulation algorithms from SPICE (Simulation Program with Integrated Circuit Emphasis) [4] is that reducing the network to a single bond is computationally less expensive. Figure 5 compares the CPU time as a function of system size solved via Frank and Lobb approach with "Double" and "Quad" precision calculations run on 512 CPU cores in parallel. The times are for a single reduction of a matrix at a single frequency. The corresponding longer times for sparse matrix techniques are shown for comparison; the Sparse Matrix involves solving the Kirchoff matrices using sparse matrix routines in Python [26] and using Berkley SPICE circuit simulation via the "NGSPICE" open-source implementation [22]. The increasing solve time with network size for "NGSPICE" in Figure 5 makes solving large networks impractical, while the Frank and Lobb has a consistent solve time t ∼ N 1.5 and is much more rapid than Sparse Matrix or NGSPICE method. This computationally efficient approach enables this work to solve large numbers of networks that are much larger than previously explored.
The outputs of the numerical simulations of the Frank and Lobb RC networks will now be described. Network results will be compared with both the logarithmic mixing (Eqs. (2) and (3)) and MEMA (Eq. (9)) to demonstrate the extent to which the MEMA formula agrees with the numerical simulations both close to the percolation threshold (p ∼ p c ) and away from the percolation threshold as the size of the network is increased.   is p c = 0.5 [28], networks with the following proportions of capacitors were selected: In Figure 6a for the case p = 0.4 there are percolation paths of resistors and not capacitors and the low frequency conductivity is frequency independent; with slope = 0. At low frequencies, where R −1 ωC, the capacitors are considered open circuit and the AC currents flow through percolated resistors. This can also be seen in Figure 6b where the phase is zero at low frequencies since both current and voltage are in phase, as would be expected for a resistor. At intermediate frequencies, R −1 ∼ ωC and the power law region can be observed where the AC conductivity rises with a slope ∼0.4, as predicted from equation (2). This region spans approximately four orders of magnitude (ω ≈ 10 4 -10 8 s −1 ). The phase of the network in the intermediate frequency region agrees well with that predicted by equation (3) which is 0.4 × −90 • = −36 • ; the frequency range at which the phase is at this value is relatively narrow compared to other studies restricted to smaller networks (e.g. N = 512 [4]). At higher frequencies, where ωC R −1 the capacitors act as short circuits, however, there are insufficient capacitors to percolate the network so the conduction paths across the network all retain a number of resistors which impose an upper limit on the admittance and set its high frequency independence (Fig. 6a); again the phase approaches zero at these high frequencies (Fig. 6b).

Initial Frank and Lobb
For the case p = 0.6 in Figure 7a there is percolation of capacitors, but not of resistors. At low frequencies the admittance continues to fall as the frequency decreases (slope = 1) despite R −1 ωC since there are insufficient resistors to percolate the network. As a result all conduction paths across the network must employ one or more capacitors which sets the low frequency admittance and its frequency dependence. This can be seen in Figure 7b where the phase is -90 • since current leads voltage in a capacitor. At intermediate frequencies (ω ≈ 10 4 -10 8 s −1 ) the power law region can again be observed and rises with frequency with a slope ≈0.6, as predicted from equation (2). The phase in this frequency range also agrees well with -55 • predicted by equation (3), see Figure 7b. The agreement in phase is best at ω ≈ 10 6 where R −1 = ωC and this is true for all network compositions; see Figures 6b, 7b and 8b. At higher frequencies, ωC R −1 , the percolation of capacitors results in the conductivity rising linearly with frequency (slope = 1) as the AC currents flow through the percolated capacitors with the resistors acting as open circuits. The percolation of capacitors at high frequencies leads to a phase of -90 • (Fig. 7b).
For the case p = 0.5 in Figure 8 there is an equal probability of the percolation of capacitors or resistors in the network. In this case the power law region of the conductivity spans a much wider frequency range with a slope of 0.5 (ω ≈ 10 2 -10 10 s −1 ) and a phase of -45 • , as predicted by equations (2) and (3), respectively. The conductivity or phase at low and high frequency depends on the existence of either percolated capacitors or resistors, as indicated in Figure 8. It is of interest to note that the frequency range of the power law response is largest for this network composition, the phase is also constant at a wider range of frequencies (Fig. 8b) and we will see later that for two-dimensional RC networks with p = 0.5 the power law region is dependent on network size [8], as predicted by the MEMA in equation (8c).

Frank and Lobb results with network size and comparison with MEMA
The low and high frequency network response of the networks in Figures 6 to 8 are clearly influenced by whether or not there are percolation paths of resistors or capacitors within the network. Previous work has been limited to smaller networks [4,5,[7][8][9]18,19] and the distributions of conductance at high and low frequencies at typically much larger than those observed in Figures 6 to 8. The outputs of the numerical Frank and Lobb networks and the MEMA approximation will be compared to show how the form of the magnitude of admittance |Y (ω)| depends on the number of components, N , in networks of different size and proportion. We compare the MEMA predictions using equation (9) with numerical results for the Frank and Lobb RC networks described in the previous section.
For our first computation we take networks with percolated resistors with p = 0.4 (Fig. 9) and consider increasing numbers of components, N , from 32 768 to 2 097 152 (S from 128 to 1024). We compute |Y (ω)| using the  (3) and (4), prediction for the networks. Unit of |Y | is Siemens. MEMA from equation (9) and compare these values with the results of Frank and Lobb computations from 1024 random realizations of each network size. In Figure 9 we see that the power-law region is typically in the range (ω ≈ 10 4 -10 8 s −1 ) and the size of the power-law region does not change significantly with the wide range of network sizes examined for p = 0.4. It can also be seen that the MEMA as indicated by the solid green line, also indicates that the power law region is also independent of network size when p = 0.4 since |1 − 2p| is relatively large, see equation (10). The broad range of conductivities at low or high frequencies is due to different percolation paths in networks of the same size and composition and the distribution at high and low frequencies clearly decreases as the network size increases; for example compare Figures 9a and 9d. This is to be expected since as the networks size increases the percolation paths will become similar in form to those in infinitely large networks. Good agreement between the MEMA and Frank-Lobb approach is observed in the power-law region. For small networks the MEMA is within the broad distribution of network results at high and low frequencies (e.g. Fig. 9a when N is 32 768); and this was also the case for smaller networks [8]. However, as the network size increases and the distribution decreases the agreement is reduced with the MEMA overestimating conductivity at low frequency and under-estimating conductivity at high frequencies; see Figure 9d.
We now examine in Figure 10 computations at p = 0.6 that contain percolation of capacitors. We present the result of computing |Y (ω)| from a number of realizations of the network for different size of N up to 2 097 152; for the purposes of brevity only data at N = 8192 and N = 2 097 152 are shown. Again the power-law region is typically in the range (ω ≈ 10 4 -10 8 s −1 ) and, as observed for p = 0.4, the range does not change significantly with network size; again this is well predicted by the MEMA (green line in Fig. 10). The range of admittances at low or high frequencies is due to different capacitive percolation paths in networks at a particular size and composition decreases as the network size increases. Good agreement between the MEMA and Frank and Lobb approach is observed in the power-law region but again the agreement is reduced at large network sizes with the MEMA overestimating the admittance at low frequency and underestimating it at high frequencies; see Figure 10b.
Data for RC networks at the critical mixing ratio of p = p c = 0.5 and N from 2048 to 2 097 152 are shown along with a comparison with the MEMA is shown in Figure 11. We focus only on the high ω limit for ease of visualisation, given that the curves are symmetrical, and the results only include the saturating cases (non-percolated capacitors), in order to analyse and compare with the saturation limits of the MEMA. In this limit ε = 1−2p = 0 and the MEMA reverts to the rigorous spectrally derived formula, equation (8) [8]. At the percolation limit, there is clearly an increase in the frequency range of the power-law as the network size increases and this is well predicted by the MEMA; see solid lines in Figure 11.

Deviation of MEMA from Frank and Lobb numerical results for p = p c
The MEMA is shown to predict the numerical results well when p = p c = 0.5 (Fig. 11) but less well when p = p c (Figs. 9 and 10), where the 1/N term in equation (9) becomes insignificant and the MEMA reverts back to the classical EMA. Our investigation of the large networks indicates that a more accurate formula would retain some N dependence up to larger network sizes and the response of the MEMA and the numerical calculations are now compared in more detail for networks where p < p c . These networks have a resistive percolation path but no capacitive path, so that the conduction becomes constant at the limits of low and high frequency, as in Figure 9, allowing ease of comparison with the MEMA. This is in contrast to networks with p > p c where the presence of percolated capacitors leads to a frequency dependent conductivity at all frequencies, as in Figure 10. It can be seen from the comparison of the MEMA with the numerical Frank and Lobb calculations for p = 0.4 that the MEMA matches the power law in the emergent region and qualitatively the shape of the curve in the transition to the percolation and saturation limits. However, it fails to capture the precise conduction values of these limits.
Finite size networks were investigated by Jonckheere and Luck [24], who gave some empirical evidence that, rather than having percolation limits as given in equations (6b) and (6c), proportional to ε or, 1 ε in the limit of |1 − 2p| 1 that they are more closely approximated by expressions of the form |Y | is proportional to |1 − 2p| ±β . From numerical results, on relatively small networks of typical size N = 512, it was suggested that β ≈ 1.
3. An examination of this is shown in Figure 12, where the Frank and Lobb numerical calculations at the low conductivity limit, |Y | ω→0 (Y F L red points with error bars) are compared to those from the MEMA, plotted both unadjusted (Y M ) as well as the whole value raised to some power β; Y β M . The value of β = 1.3 [24] appears to be approximately accurate in the small N limit (points labelled Y 1.3 M ), but our ability to explore large networks indicates that as N becomes very large the exponent tends instead to β ≈ 1.11 (points labelled Y 1.11 M ) for networks at p = 0.4. This is an observation of interest that suggests that it should be the focus of analytical attention for future investigations as the reasons for this are currently unknown. We will see later that the value of β has some dependency on network composition and β ≈ 1.15 is more appropriate for a range of network compositions.

The dynamic range (|Y| ω→∞ /|Y| ω→0 ) for p = p c
For networks on the percolation threshold p = p c = 0.5 the maximum conductance at high frequency scales strongly with N , as can been seen in Figure 11. This data focusses only on the high ω limit for ease of visualisation; given that the curves are symmetrical. It can be seen that the curves given by the MEMA formula capture the upper limits of the distribution of numerical results, which can be understood heuristically in that realisations that are closest to the percolation limit will be those with the minimum number of conductive (or maximum capacitive) bonds needed to percolate, and any statistical deviation from this will result in smaller conductance values. We define here the dynamic range of the network,Ŷ , as:  We observe from this figure that there is a power-law mixing region, which evolves into a percolation region. The boundary for the mixing region increases (as predicted by MEMA) with the network size.
Following the results in equations (6b) and (6c) and also equations (8a) and (8b) we predict that the dynamic range is given by, if N is large, and N ε 2 1 (13a)  Note that this result can be obtained directly from equation (9a) by takingŶ = φ 2 so that equation (9a) allows us to see how the dynamic range depends on both the closeness to criticality and the network size. In Figure 13a,Ŷ is plotted as a function of N for a range of capacitor fractions with p < 0.5 (percolated resistor networks). If p = 0.5 then ε = 0 and, as predicted by equations (13a),Ŷ is expected to be directly proportional to N for all values of N . However, if p < 0.5 thenŶ is only directly proportional to N for smaller values of N, but becomes asymptotic to a finite valueŶ (p) as N → ∞. The transition between these two forms of behaviour is seen in Figure 13a and occurs when N = |1 − 2p| −2 , as indicated by equations (10) and (13a). This is best observed in Figure 13a for the cases close to percolation at p = 0.49 and p = 0.48 whereŶ is predicted by equation (10) to be proportional to network size up to N = 10 4 and N = 10 3.4 respectively; this confirms the results in [8] up to much larger values of N than previously explored. For the case of p = 0.4 the transition is estimated by equation (10) at N = 10 2 and thereforeŶ is independent of N in Figure 13a (and Fig. 9). This can be of interest in real systems where it is possible to determine the necessary size of the effective material/composite RC network to achieve properties that are independent of the number of components.
It can also be seen that, in the transition region, the distribution ofŶ from the random Frank and Lobb realisations is large for relatively small networks and we will see later that the standard deviation of the network conductivity falls with 1/N 0.5 (or 1/S). The size dependent term in the MEMA (Eq. (9)) falls with 1/N , therefore an important point to notice in Figure 13a is that this results in the overshoot ofŶ of the Frank and Lobb network simulations at intermediate values of N compared to its final asymptotic value at large N ; i.e. once the standard deviation shrinks with the network size.
There is clearly a discrepancy between the Frank and Lobb numerical data ofŶ and the solid MEMA lines with the data points in Figure 13. As discussed in Section 3.3, previous suggestions by Jonckheere and Luck for relatively small networks (N = 512) state that the EMA can be corrected by raising its predicted conductance to the power of 1.3 [24] seem to be valid only for low to intermediate values of N , but a better fit across the range of p values considered appears to be some intermediate value, in this case 1.15. This is confirmed by consideringŶ as a function of p as p → 0.5 in Figure 13b, showing its limiting value as N → ∞ (S is increased from 16 to 1024 in the figure). It has been shown in previous work [8] that the MEMA predictsŶ ∼ |1 − 2p| −2 with some accuracy, althougĥ Y ∼ 1/|1 − 2p| −2.6 was a better fit for networks; although the study was limited to N = 10 000. Our ability to explore larger networks in this work show thatŶ ∼ |1−2p| −2 andŶ ∼ 1/|1 − 2p| −2.6 under and over-predict, respectively, except for low N , where the statistical variation is a dominating factor, while a better fit across the range of p is forŶ ∼ |1 − 2p| −2.3 . However, it is clear that there remains a discrepancy between the MEMA, simply raised to any power, and the numerical results. It can be speculated that the power itself has some dependence on p, as might be expected by re-considering spectral derivations of the conductance, such as that in [8], or a p dependence in the right hand term of the MEMA.

Scaling of network variation with size
In Figure 14 we plot the standard deviation (σ) and the standard deviation multiplied by the network size (Sσ)  of the admittance with ω for different network sizes S, for p = 0.4. Both (σ) and (Sσ) are smallest at the condition R −1 = ωC since in this frequency range AC currents flow through both resistors and capacitors so that there is minimal variation between random networks of the same composition. As the frequency is increased or decreased relative to the frequency at R −1 = ωC the standard deviation increases linearly and eventually plateaus to a constant value at low and high frequencies where there is percolation, see Figure 14a. Interestingly the product (Sσ) for networks of different sizes fall on the same line in Figure 14b. The origin of this behaviour stems from the binary choice of components in the networks being taken from a Bernoulli probability distribution since they are created with each of the N = 2S 2 components used to fill the network lattice chosen randomly to be either a capacitor or resistor with probability p or (1 − p). This results in a normal distribution for the number of resistors and capacitors over the whole network. The mean number of the capacitors is therefore pN and the variance is p(1 − p)N in each bond. As would be expected from central limit theorem the variance of the mean is p(1 − p)/N (strong law of large numbers) and the standard deviation is p(1 − p)/N or p(1 − p) 2S 2 ; therefore Sσ is expected to be constant. This simple method of understanding network size on the distribution of conductivity, or other properties such as permittivity, is of interest to understand the influence of sample size or variability is properties of materials or composite systems.

Conclusions
In this paper we have investigated and compared models representing the conduction flow-paths though disor-dered random composite media in terms of new MEMA formula and RC networks that are much larger than previously explored. An efficient algorithm based on the Frank and Lobb reduction method has been used to study in detail the frequency dependent admittance of large twodimensional square lattice resistor-capacitor networks. The approach employed involves a sequence of transformations to reduce the lattice to a single bond that has the same admittance as the network and this has allowed networks to be studied with up to 2 × 10 6 resistor or capacitor components. The ability to assess such large networks has allowed us to study the statistical variation of the network properties with network size by analysis up to 1024 networks of a specific size, and component fraction. Results from the networks in terms of emergent powerlaw scaling of the conductance and the percolation and saturation limits at the high and low frequency bounds were compared with a Modified Effective Medium Approximation (MEMA) that extends the Effective Medium Approximation (EMA) to finite network sizes. It shows that new MEMA formula gives good results over a wide range of values of network compositions and sizes, and provides a simple analytical formula for the network response that accurately predicts the effects of (finite) network size and composition, while at the same time showing where it remains to be improved and explained. It is also seen that: (i) Good agreement between the power law mixing rules, MEMA and the Frank and Lobb numerical results is observed in the power-law region. The power law mixing rules described earlier are limited to providing an estimate of the properties of RC networks at intermediate frequencies where R −1 ≈ ωC. (ii) For small networks the MEMA is within the broad distribution of network results at high and low frequencies However, as the network size is increased to the large sizes explored in this work the distribution of conductance decreases significantly and the agreement is reduced As an example for p = 0.4 (percolated resistors) the MEMA overestimates conductivity at low frequency and under-estimates conductivity at high frequencies, so that the dynamic range is not predicted accurately. Our results over a large range of network sizes indicate the percolation admittances are more closely approximated by expressions of the form |Y | proportional to |1 − 2p| ±β , β = 1.11 when p = 0.4 while the previously reported value of β = 1.3 [24] appears to be approximately accurate in the small N limit. (iii) For networks on the percolation threshold p = p c = 0.5 the conductance given by the MEMA scales well with the network size, N, with the dynamic range, |Y | ω→∞ /|Y | ω→0 , being directly proportional to N for all values of N . If p < 0.5 thenŶ is only directly proportional to N for small values of N and p close to the percolation limit and becomes asymptotic to a finite value as the network size increases. A transition between these two forms at N = |1 − 2p| −2 is observed and confirms the results in [8], but up to much larger values of N than previously explored. (vi) Previous suggestions that the EMA could be corrected by raising its predicted conductance to the power of 1.3 [24] seem to be only valid for low to intermediate values of N , but our ability to explore very large values indicate a better fit appears to be some intermediate value, here 1.15, for p < 0.5 so that the dynamic range,Ŷ is proportional to |1 − 2p| −2.3 . (v) The distributions of network admittance values at high and low frequencies decrease as the network size increase and a scaling of the conductivity distribution was observed. The standard deviation of the distribution of admittance values was found to be inversely proportional to the network size and the product (Sσ) for networks of different sizes fall on the same line. The origin of this behaviour stems from the binary choice of components in the networks being taken, which results in a normal distribution for the number of resistors and capacitors over the whole network. The scaling condition is explained by the network having the variance of the mean being p(1 − p)/N (strong law of large numbers) and the standard deviation is therefore p(1 − p)/N .
The data from the very large two dimensional square resistor-capacitor networks provides a deeper insight into the responses of such networks and their relationships to the limits of analytical approximations, such as logarithmic mixing and our MEMA for prediction of network response. Methods of correction to improve analytical predications are proposed. While we have focussed on network conductance, the approach can also be easily used to explore other frequency dependent properties, such as capacitance. The data also elucidates the statistical nature of the networks as a function of frequency, composition and size. It should be highlighted that in many cases materials exhibit three-dimensional conduction where the percolation limits are different. For example there is a possibility to achieve percolation of both resistors and capacitors in the same 3D network which is a common characteristic of the universal dielectric response (UDR), such as in Figure 1a. This data on large RC networks all contributes to a deeper understanding of the frequency dependent properties of real heterogeneous materials and composite material systems to understand size effects, composition effects and statistical variation of properties.