{'fraction_non_alphanumeric': 0.04606101491642306, 'fraction_numerical': 0.017558653774079905, 'mean_word_length': 4.194685096855836, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dielectric resonances of binary random networks

We investigate the AC conductivity of binary random impedance networks, with emphasis on its dependence on the ratio of the complex conductances of both phases. We propose an algorithm to determine the conductance of a finite network, in terms of its poles and of the associated residues. A numerical implementation of the algorithm, on the example of the square lattice, allows a detailed investigation of the resonant dielectric response of the binary model, including the density of resonances, the spectral function, the distribution of spacings between resonances. The distribution of local electric fields at resonance is found to be multifractal.

Abstract. We investigate the AC conductivity of binary random impedance networks, with emphasis on its dependence on the ratio h = σ 1 /σ 0 , with σ 0 and σ 1 being the complex conductances of both phases, occurring with respective probabilities p and 1 − p. We propose an algorithm to determine the rational h-dependence of the conductance of a finite network, in terms of its poles and of the associated residues. The poles, which lie on the negative real h-axis, are called resonances, since they show up as narrow resonances in the AC conductance of the RL − C model of a metal-dielectric composite with a high quality factor Q. This approach is an extension of a previous work devoted to the dielectric resonances of isolated finite clusters. A numerical implementation of the algorithm, on the example of the square lattice, allows a detailed investigation of the resonant dielectric response of the binary model, including the p-dependence of the density of resonances and the associated spectral function, the Lifshitz behaviour of these quantities near the endpoints of the spectrum of resonances, the distribution of spacings between neighbouring resonances, and the Q-dependence of the fraction of visible resonances in the RL−C model. The distribution of the local electric fields at resonance is found to be multifractal. This result is put in perspective with the giant surface-enhanced Raman scattering observed e.g. in semicontinuous metal films.

INTRODUCTION
Random networks of complex impedances are currently used to model electrical and optical properties of disordered inhomogeneous media. The most common situation is that of a binary composite medium, modeled by attributing a random conductance to each bond (x, y) of a lattice, according to the binary law σ x,y = σ 0 with probability p, σ 1 with probability q = 1 − p, (1.1) in correspondence with the bond percolation problem (see ref. [1] for a review).
The conductances (inverse impedances, or admittances) σ 0 and σ 1 of both phases take arbitrary frequency-dependent complex values. Hereafter we follow the notations of ref. [2]. The dimensionless complex ratio and the concentration p are the essential parameters of the model. As far as static (DC) properties are concerned, the limiting case h = 0 embraces the conductor-insulator mixture (σ 1 = 0) and the superconductor-conductor mixture (σ 0 = ∞). In both situations the conductivity exhibits power-law behaviour for p close to the critical concentration p c , corresponding to the percolation threshold.
More generally, the conductivity of the binary model obeys a scaling law in the critical regions defined by |h| ≪ 1 and |p − p c | ≪ 1, or equivalently |h| ≫ 1 and |q − p c | ≪ 1 [1]. Frequency-dependent (AC) properties of metal-dielectric composites, such as cermets or thin films, are often investigated by means of either the R−C and the RL−C models [3]. In both cases the low-frequency regime corresponds to |h| ≫ 1, and thus to critical behaviour when the metallic concentration q is near the percolation threshold.
The purpose of the present work is to shed some new light on the resonant behaviour of the binary model. The emphasis will be put on the analytic structure of the conductivity in the ratio h, or in the equivalent complex variable For any value of the concentration p, the conductivity has singularities for h real and negative, i.e., in the range 0 ≤ λ ≤ 1. A quantitative investigation of these singularities, and of related quantities, is the main goal of the present paper. Our aim is twofold. First, the analytic structure of the conductivity of the binary model is a classical subject, since the developments of the Bergman-Milton theory [4,5]. The key ingredient of this formalism is the spectral function H(p, x), which has only been the subject of a limited number of investigations so far [6,7]. The present approach provides a direct accurate numerical evaluation of the spectral function of binary random networks. Second, the singularities of the conductivity have a physical interpretation in terms of dielectric resonances in the RL − C model, and of the relaxation times in the transient response in the R − C model. The regime of dielectric resonances has been argued to provide a natural explanation for the anomalous fluctuations of the local electric field [8], which are responsible for giant surface-enhanced Raman scattering observed e.g. in semicontinuous metal films [9].
The dielectric resonances of isolated clusters have been investigated in ref. [2]. The situation considered there was a finite set of (metallic) bonds with conductance σ 1 , embedded in an infinite (dielectric) host lattice, whose bonds have a conductance σ 0 . It was shown there that the conductance of such a system is entirely characterised by a finite number of resonances. The positions λ a of the resonances, and the associated cross-sections γ a , are expressed in terms of the eigenvalues and eigenvectors of a finite matrix M, dictated by the geometry of the clusters. The present paper is an extension of the method of ref. [2] to an arbitrary binary network. We shall make use of an efficient algorithm, which allows for an exact determination of all the resonances of a finite sample. The setup of this paper is as follows. In section 2, we gather definitions and results on various features of the conductivity of the binary model. Section 3 is devoted to the presentation of the algorithm. Section 4 contains a variety of numerical results, concerning especially the spectral function of the Bergman-Milton theory and the spectral density of resonances, their Lifshitz behaviour near the endpoints of the spectrum (λ = 0 and λ = 1), the number of visible resonances of a finite sample and its dependence on the quality factor Q of the RL − C model, and the distribution of the local electric fields at resonance, which turns out to be multifractal. The bonds with a conductance σ 0 are called the P-bonds, which form the P-set. There are n P of them, among which n H P are horizontal and n V P are vertical. Similarly, the bonds with a conductance σ 1 are called the Q-bonds, which form the Q-set. There are n Q of them, among which n H Q are horizontal and n V Q are vertical. Let Y be the conductance (admittance) of the network, measured between the plane electrodes (bus bars) shown in Figure 1, and let Z = 1/Y be its impedance. These quantities are rational functions of the dimensionless complex variables h or λ. Let us anticipate that it is more convenient to use the variable λ. We have The prefactors of these expressions are the conductance and the impedance of the uniform network whose bonds have a conductance σ 0 , since λ = ∞ corresponds to σ 1 = σ 0 (no disorder). The λ a involved in the product expressions are the poles of the conductance and the zeros of the impedance, while the λ a are the zeros of the conductance and the poles of the impedance. The number n R of poles and zeros depends on the configuration of the random bonds. It is bounded by the number of sites of the network: 0 ≤ n R ≤ n S . Finally, the partial-fraction expansions involve residues α a and β a .
It follows from considerations on the dissipated power that Y /σ 0 is a Stieltjes function, namely its imaginary part has the sign of Im h, or equivalently of Im λ [4,5,1]. This property implies that the poles and zeros of Y and of Z alternate, according to We have λ 1 = 0 if, and only if, the Q-set is conducting, i.e., it contains at least one connected path between both electrodes. The conductance of this Q-set then reads Similarly, we have λ n R = 1 if, and only if, the P-set is not conducting.
The Stieltjes property also implies that the residues α a and β a are positive. These two sets of residues are different from each other in general. They obey, however, a remarkable sum rule, which can be proved by expanding the conductance Y around λ = ∞, to first order in 1/λ. Indeed, 1/λ = (σ 0 − σ 1 )/σ 0 is the dimensionless contrast between the conductances of both phases. As a consequence, for a given realisation of the random network, up to first order in 1/λ included, the conductance reads Y = (N/M )σ H , where is the average value of the conductances of the horizontal bonds of the network. By inserting this estimate into the second expressions of eqs. (2.1a, 1b), we obtain the relation In the case of isolated finite clusters [2], the poles and zeros of the conductance form very tight doublets. We have indeed The sum rule (2.5) then becomes where n H is the number of horizontal bonds in the clusters (perpendicular to the electrodes). This identity was not noticed in ref. [2].

R − C model
The R − C model, already mentioned in the Introduction, is defined as follows. The Q-set is a metallic phase, whose bonds consist of a pure resistance R, while the P-set is a dielectric phase, whose bonds consist of a perfect capacitance C. The complex conductances at frequency f = ω/(2π) thus read Along the lines of refs. [1,2], we introduce the microscopic relaxation time We thus have h → ∞ and λ → 0 at low frequency.
The poles of the conductance Y show up as relaxation times in the transient response of the model [1]. Consider indeed the R − C model on a finite network, submitted to a voltage V (t) = V 0 θ(t) applied between the electrodes, with θ(t) being the Heaviside stepfunction. The intensity I(t) across the network in the transient regime can be evaluated by means of the Fourier transformation. We thus obtain The second expression of eq. (2.1a) for the conductance Y , with λ given in eq. (2.10), yields the result where the relaxation times read If the metallic Q-phase is conducting, we have λ 1 = 0. The term a = 1 in eq. (2.12) yields the DC current through the network, I 0 = N V 0 α 1 /(M R), in agreement with the result (2.3).

RL − C model
In the RL − C model, already mentioned in the Introduction, the metallic bonds of the Q-set now consist of an inductance L in series with a weak resistance R, while the dielectric bonds of the P-set still consist of a perfect capacitance C. The conductances at frequency f = ω/(2π) now read (2.14) Along the lines of refs. [1][2][3], we introduce the microscopic resonance frequency 16) and the quality factor which is a dimensionless measure of the dissipation rate.
In the following, we shall mostly consider the case of a weak dissipation, corresponding to a large quality factor (Q ≫ 1). We have then We notice that the low-frequency regime (y → 0) again corresponds to h → ∞ and λ → 0.
Since the variables h and λ have a small (negative) imaginary part, proportional to 1/Q, the poles of the conductance Y and of the impedance Z of the network show up in the frequency dependence of these quantities as narrow resonances, respectively located at ω a = ω 0 y a and ω a = ω 0 y a , with Eq. (2.20a) shows that the real part Re Y (ω) of the conductance exhibits n R narrow resonances, at ω a = ω 0 y a . The resonance peaks have a Lorentzian shape, with a common absolute width ∆ω = ω 0 /(2Q). The maxima at resonance read while the area under each resonance peak is (

2.22)
A similar pattern of resonances can be observed on the real part Re Z(ω) of the impedance (2.20b), with maxima at resonance and areas under resonance peaks respectively given by (2.23)

Spectral function and density of resonances
The conductance Y of an infinitely large network drawn on any regular lattice becomes a self-averaging quantity. In other words, the conductivity is an intrinsic characteristic of the binary model, which depends on the lattice under consideration, on the concentration p, and on the complex variable h or λ. The conductivity of the binary model has been the subject of many investigations [1]. We gather below some definitions and properties which will be useful in the sequel.
The second expression of eq. (2.1a) yields the Bergman-Milton integral representation [4,5] Σ(p, λ) = σ 0 1 − is called the spectral function of the binary model. In more precise terms, it is the density of a positive measure, supported by the interval 0 ≤ x ≤ 1. Indeed, as a consequence of eq. (2.3), the spectral function has a singular component at x = 0, of the form with the notation to be introduced in eq. (2.38), whenever the Q-set is conducting, i.e., for 1 − p = q > p c .
The spectral function entirely determines the conductivity Σ(p, λ), by means of the integral representation (2.25). Conversely, it is given by the inverse formula In other words, the conductivity is analytic in the complex λ-plane cut along the interval 0 ≤ λ ≤ 1, and its discontinuity along this cut reads Disc Σ(p, λ) = 2iπσ 0 H(p, λ). (2.29) The spectral function directly yields the transient intensity I(t) of the R − C model on a very large network, namely To close up, we introduce the spectral density of resonances (2.31) and the total density of resonances representing the mean number of resonances per site. This quantity has a non-trivial dependence on the concentration p. It will be expressed in eq. (3.14) in terms of geometrical quantities of the bond percolation problem.

Homogeneity and duality
First, the conductivity of the binary model is invariant under the simultaneous interchange p ↔ q = 1 − p, σ 0 ↔ σ 1 . Indeed, let Y = σ 0 F (h) = σ 0 f (λ) be the conductance of the network G, shown in Figure 1. Then the network G ′ , obtained from G by interchanging all the bond conductances according to σ 0 ↔ σ 1 , has a conductance . For a large enough network, if G is typical of the concentration p, then G ′ is typical of the concentration 1−p. The conductivity, the spectral function, and the density of resonances therefore obey the following identities, on any regular lattice More interestingly, the square lattice is self-dual, i.e., invariant under the geometric transformation called duality. This concept has been introduced in physics by Kramers and Wannier [10], while its consequences on random resistor networks have been explored in a systematic way by Straley [11]. The dual G of a planar network G has bonds which cross those of G, the conductances of any pair of crossing bonds being inverse to each other. The duality property implies that the conductance of the whole network G reads Y = 1/Y = Z. With the same notations as above, we Because of its non-linearity, the duality identity (2.34) for the conductivity does not yield any identity involving the spectral function H(p, x) only.
The above identities can be combined with the homogeneity properties (2.33). We thus obtain that the percolation threshold is p c = 1/2, and that the conductivity and the spectral function right at this point read Furthermore, the spectral and total densities of resonances obey the relations

Critical behaviour and scaling
As mentioned in the Introduction, the conductivity exhibits scaling behaviour around the percolation threshold p = p c . The conductivity of the conductor-insulator mixture vanishes for p ≤ p c , while it reads Similarly, the conductivity of the superconductor-conductor mixture is infinite for p ≥ p c , while it reads Both amplitudes A(p) and B(p) have a power-law behaviour for p near p c : The conductivity of the binary mixture obeys a scaling law in the critical region defined by |h| ≪ 1 and |p − p c | ≪ 1, of the form [1] where the Φ ± are scaling functions of one complex variable, with ± referring to the sign of p − p c . The homogeneity relation (2.33a) allows to describe the vicinity of the other critical point, |h| ≫ 1 and |q − p c | ≪ 1.
As a consequence of eq. (2.28), the spectral function also obeys a scaling law of the form (2.41), namely (2.44) for 1 − x ≪ 1 and |p − p c | ≪ 1, and a similar law around the other critical point, i.e., for x ≪ 1 and |q − p c | ≪ 1. We shall come back to the scaling law (2.44) in section 2.9.
Consider now the binary model on a large but finite sample, of size M × N . In the critical region, its mean conductance obeys the finite-size scaling law where Ψ is a three-variable scaling function. As a consequence, right at the percolation threshold, the critical region extends over a range On the square lattice, the duality symmetry implies p c = 1/2, and A(p)B(1 − p) = 1, hence s = t and u = 1/2, in agreement with eq. (2.36). The common numerical value of these exponents is s/ν = t/ν = 0.9745 ± 0.0015 [12], with the exponent of the correlation length being exactly ν = 4/3, hence s = t = 1.300. Hence, for p = p c = 1/2, the same critical singularity simultaneously affects both endpoints of the spectrum, λ = 1 and λ = 0, corresponding respectively to h = 0 and h = ∞, again in agreement with eq. (2.36).

Sum rules
The representation (2.25) of the conductivity can be expanded as the following series in inverse powers of λ where the coefficients are the moments of the spectral function H(p, x).
As we have already noticed in section 2.1, the expansion variable 1/λ is the dimensionless contrast between the conductances of both phases. As a consequence, the expansion (2.47) can be viewed as a special case of the weak-disorder expansion of the conductivity of a random network with an arbitrary narrow distribution of bond conductances [13]. The general results to sixth order derived there can be transcribed in the present case of binary disorder on the square lattice. We thus obtain the following expressions for the first six moments, , with the notations of ref. [13], is the only non-trivial quantity occurring in the sixth-order expansion. We have thus derived explicit sum rules for the spectral function, which agree with the expressions given in ref. [7]. The number J, denoted there as 1 − a 4 5 , was estimated there from numerical data to be J = 0.9 ± 0.5.

Effective-medium approximation
The effective-medium approximation (EMA), introduced by Bruggeman in 1935 [14], is a self-consistent approximate scheme to evaluate the conductivity of random impedance networks [15,16,13], which is still being very widely used [17].
In the present case of the binary model on the square lattice, the EMA prediction for the conductivity is given by This EMA formula for the conductivity is analytic in the λ-plane cut along the interval λ min ≤ λ ≤ λ max , with The prediction for the associated spectral function [cf. eq. (2.28)] reads The EMA formula gives a very accurate approximation to the conductivity of the binary model in generic circumstances. For instance, the 1/λ-expansion of the EMA formula (2.51) gives expressions H EMA k (p) for the moments which only differ from the true results (2.49), starting with H 4 (p), by replacing J by J EMA = 1 [13]. The EMA scheme also respects the duality symmetry (2.34). The predictions (2.51), (2.53) therefore agree with the exact results (2.36) for the conductivity and the spectral function at p = p c = 1/2. The EMA also correctly predicts the equality s = t, but not the common value of these exponents (s = t = 1 instead of 1.300).
For a generic value of the concentration (p = p c ), the endpoints (2.52) are such that 0 < λ min < λ max < 1. The support of the EMA prediction (2.53) for the spectral function thus does not extend over the whole interval 0 ≤ λ ≤ 1. In the vicinity of the percolation threshold (p → p c = 1/2), we have λ min = 1 − λ max ≈ (p − p c ) 2 , in agreement with the estimate (2.43) of the crossover scale h ⋆ , with s + t = 2.

Lifshitz tails
It has been argued [6] that the spectral function H(p, x) of the true conductivity of the binary model extends over the whole allowed spectrum 0 ≤ λ ≤ 1, for any value of the concentration p, at variance with the EMA formula (2.53).
This argument has then been put in perspective [1] with Lifshitz singularities [18][19][20]. These singularities are caused by the presence of very large, and thus very improbable, ordered regions in a randomly disordered system. The original example considered by Lifshitz [18] is that of the phonon spectrum of a binary harmonic alloy, consisting of light atoms, with mass m and concentration p, and heavy atoms, with mass M > m and concentration q = 1 −p. Lifshitz has argued that the vicinity of the upper edge ω max of the phonon spectrum of the alloy is dominated by large ordered regions, almost spherical in shape, consisting only of light atoms. He thus showed that ω max coincides with the upper edge of the pure lattice consisting only of light atoms, and that the density of states of the alloy vanishes exponentially fast near ω max , as where c is a lattice-dependent constant, which can be evaluated exactly. Along this line of thought, it has been argued in ref. [1] that the spectral function of the binary model has an exponentially small Lifshitz tail, extending all the way to the endpoints λ = 0 and λ = 1 of the spectrum of resonances, of the form and similarly for x → 1. This expression was rather conjectural, as the determination of the relevant ordered regions was left as an open question, so that the prefactor C(p) was not predicted.
Hesselbo [21] then argued that the relevant ordered regions are hairpin configurations, as shown in Figure 2a. Let Y n be the transversal conductance of the hairpin consisting of n cells, measured between the point electrodes shown in Figure 2a, considered as an isolated network (not embedded in the square lattice). This quantity can be evaluated from the recursion relation Consider first a value of h not on the negative real axis, and set with Re µ > 0 and |Im µ| < π. The Möbius map involved in the recursion (2.56) has two fixed points, For a strong dielectric contrast (h → 0), the correlation length of currents along the hairpin diverges according to ξ = 1/(2µ) ≈ (2h) −1/2 . In this regime, the conductance of long hairpins scales as This scaling form of the conductance exhibits an infinite array of alternating zeros, lying at and poles, lying at with a ≥ 1. The first zero λ 1 , corresponding to a pole in the dual configuration, yields, according to Hesselbo, the Lifshitz behaviour of the spectral function, at least for a small enough concentration p.
It turns out that the formula (2.59) also gives the conductance of worm-like networks, such as the configuration shown in Figure 2b, where the square cells are put together in any random fashion, respecting the linear structure and the constraint of self-avoidance. The number of such worm-like configurations with n cells is of order exp(nS), with S being the associated configurational entropy. On the other hand, a hairpin with n cells occurs with a probability of order p 2n , at least for p small enough. Altogether, the first zero of eq. (2.61) is expected to show up with a probability weight of order (p 2 e S ) n . By eliminating the number n between the above estimates, we obtain the following analytical form for the Lifshitz tail of the density of resonances and of the spectral function and a similar formula for x → 1. This result, with an inverse-square-root behaviour in the spectral variable x, is characteristic of Lifshitz tails in one-dimensional systems [19]. This is due to the fact that the relevant structures are linear objects. Furthermore, the identities (2.33), (2.37) imply C(p) = C(1 − p). The above argument also leads to a prediction for the small-p behaviour of the amplitude C(p), namely These predictions will be compared to numerical data in section 4.
The Lifshitz tail of the spectral function manifests itself in the long-time tail of the transient intensity response of the R − C model. Indeed, in the regime of long times (t ≫ τ ), eq. (2.30) is dominated by the vicinity of x = 0, hence In the critical region, the prediction (2.63) for the Lifshitz tail is compatible with the scaling law (2.44) for the spectral function, provided the amplitude C(p) vanishes near the percolation threshold, according to The Lifshitz behaviour (2.63) is expected to hold only deep in the tails, for x (or 1 − x) much smaller than the crossover scale h ⋆ , defined in eq. (2.43).

ALGORITHM
We now turn to the presentation of our algorithm for evaluating the rational hdependence of the conductance of a finite binary network, such as that shown in Figure 1. This approach is an extension of the method of ref. [2]. It turns out that a very similar formalism was proposed by Straley [6] some twenty years ago, but this work has apparently not been noticed since then. The conductance will be determined in the second form of eq. (2.1a), namely we shall calculate first the poles λ a of the conductance, giving the positions of the resonances, and then the associated residues α a , giving the strengths of the resonances.

Generalities
Along the lines of ref. [2], our starting point is the Kirchhoff equations for V x , the electric potential at site x: There is one such equation per site x inside the network. The notation y(x) means that y is a neighbour of x, i.e., there exists a bond (x, y), and the sum possibly includes sites y belonging to either electrode. Eqs. (3.1) have to be complemented by the boundary conditions V y = 0 for the sites y on the left electrode, and V y = V 0 for the sites y on the right electrode.
We define the topological Laplace operator ∆ on the network as again with the convention that the sum possibly includes sites y belonging to either electrode, in which case the Dirichlet boundary condition V y = 0 is assumed. The operator ∆ can be written as the sum ∆ = ∆ P + ∆ Q of its components ∆ P on the P-set and ∆ Q on the Q-set, respectively defined as where y ∈ P(x) (respectively, y ∈ Q(x)) means that (x, y) is a P-bond (respectively, a Q-bond). We also introduce the quantities   5) or equivalently, in vector and matrix notation, This reduced form only involves the complex variable λ defined in eq. (1.3). The conductance of the network is given by where the total current I flowing into the network from the left electrode reads 8) or equivalently, in vector notation, (3.9)

Poles of the conductance
In analogy with ref. [2], the poles of the conductance are the non-trivial values λ a of λ for which the homogeneous Kirchhoff equations (3.1) with V 0 = 0 have a non-zero solution, namely This is a well-posed generalised eigenvalue problem, since ∆ and ∆ Q are two real symmetric matrices, of size n S × n S , and (−∆) is a positive definite matrix.
It turns out that the endpoints λ = 0 or λ = 1 are in general extensively degenerate eigenvalues of eq. (3.10). These eigenvalues do not correspond to resonances. Indeed, we know from section 2.1 that the conductance has no pole at λ = 1, while it has a simple pole at λ = 0 if, and only if, the Q-phase is conducting. Let us set ε Q = 1 in this situation, ε Q = 0 else. The number of resonances then reads where n 0 and n 1 denote the respective multiplicities of the endpoint eigenvalues λ = 0 and λ = 1. More precisely, n 0 is the number of zero-modes of the operator ∆ Q , i.e., the dimension of its kernel. Since ∆ Q has one zero-mode per cluster C of the Q-phase which is disconnected from the electrodes, we have In this formula, the sum runs over the clusters C of the Q-phase, s(C) denotes the number of sites of the cluster C, and the characteristic function χ(C) is unity if the cluster C overlaps with either of the electrodes, and zero else. The multiplicity n 1 of the eigenvalue λ = 1 can be expressed similarly, in terms of the clusters of the P-phase.
In the thermodynamic limit, one can derive from eq. (3.12) expressions for the fractions of eigenvalues which are condensed at λ = 0 and λ = 1. Indeed, the terms χ(C) are negligible, while the other terms can be expressed as functions of geometrical characteristics of the bond percolation problem [22], namely where n c (p) is the mean number of finite clusters per site, while P (p), the percolation probability, is the probability for any given bond to belong to the infinite cluster. The latter quantity is non-vanishing only for p > p c .
The density of resonances ρ R (p), defined in eq. (2.32), then reads This quantity is symmetric in the exchange p ↔ q, in agreement with eq. (2.37). For a small concentration p, the contributions to eq. (3.12) of all the P-clusters and Q-clusters consisting of up to four bonds can be enumerated by hand. We thus obtain At the percolation threshold, ρ R (p) takes its maximal value, which can be determined as follows. The percolation probability P (p c ) vanishes, while n c (p c ) is known exactly [23], hence (3.16)

Residues of the conductance
At any resonance corresponding to a non-trivial eigenvalue (λ a = 0 and 1) of eq. (3.10), the map of the electric potentials on the network is given by the components (X a ) x of the associated right eigenvector X a . Since eq. (3.10) is symmetric, the X a are simultaneously its left and its right eigenvectors: with the row vector X t being the transposed of the column vector X. The n S eigenvectors {X a } form a basis. They are orthogonal to each other with respect to the metric (−∆), namely X t a (−∆)X b = 0 for a = b. We normalise them as The squared norm of the eigenvector X a reads X t a (−∆)X a = (x,y) where is the local electric field on the bond (x, y).
The solution V to the inhomogeneous Kirchhoff equations (3.1), (3.6) can be obtained in terms of the eigenvectors X a . Indeed, let us expand V on the basis of the X a : (3.21) By inserting this expansion into eq. (3.6), and multiplying to the left by X t b , we readily obtain the amplitude c b (λ) in the form Finally, by inserting this result into eqs. (3.9), (3.7), we obtain the following formula for the residues α a ≥ 0 of the conductance at its non-trivial poles (λ a = 0 and 1) As recalled in the beginning of section 3.1, the conductance also exhibits a simple pole at λ = 0, if the Q-phase is conducting. The corresponding residue α 1 , which yields the conductance of the Q-phase by means of eq. (2.3), is not directly given by eq. (3.23). It can, however, be determined from the other ones (α 2 , · · · , α n R ), by using the sum rule (2.5).

NUMERICAL RESULTS
We have shown in section 3 that the full rational λ-dependence of the conductance of a finite binary network can be expressed in terms of the generalised eigenvalues and eigenvectors of eq. (3.10). We have implemented this algorithm numerically, using the IMSL routine EIGZS, in order to obtain numerical data concerning several quantities of interest, which will be discuss successively throughout this section. The CPU time for solving the spectral problem for each sample grows rapidly with the system size, proportionally to n 3 S , i.e., to N 6 for a square sample of size N × N .
The optimal use of our algorithm therefore consists in obtaining good statistics on samples of moderate sizes. We have commonly used sample sizes such as N = 16 or N = 20, and a statistical ensemble of several times 10 4 samples, having some 10 7 random bond conductances in total, a good enough statistics to obtain very accurate data. No observable systematic finite-size effects have been found, even in the critical regions (p → p c and λ → 0 or λ → 1). This observation is in agreement with the argument, given at the end of section 4.1, showing that the critical regions are indeed very small in the variable λ.

Density of resonances and spectral function
We have evaluated the spectral density of resonances ρ(p, x) and the spectral function H(p, x), by means of their respective definitions (2.31) and (2.26). The spectral function, which plays a central role in the Bergman-Milton formalism [4,5], had only been the subject of a limited amount of work. In ref. [7] it has been extracted from the imaginary part of the conductance, measured at a small but finite distance ǫ from the negative real h-axis. The present method yields a direct measurement of the spectral function, avoiding especially any contamination from the delta-function at x = 0 [see eq. (2.27)], which can be discarded in an exact way.  Figure 4, is symmetric under x ↔ 1−x. Both ρ(p, x) and xH(p, x) exhibit a rich structure, down to the scale of the resolution (each plot contains 100 bins). It will be demonstrated more clearly in section 4.2 that they are non-vanishing over the whole allowed spectrum (0 ≤ x ≤ 1).
For a small enough concentration, the most salient structures in ρ(p, x) and H(p, x) can be predicted from the analysis of the resonances of isolated finite clusters. For consistency with ref. [2], we shall consider the regime q = 1 − p → 0. To leading order in q, the relevant configuration consists of one isolated Q-bond embedded in a host lattice consisting of P-bonds. This one-bond cluster, shown as configuration A in Figure 5, has one single resonance, at λ = λ A1 = 1/2, yielding the observed leading peak in ρ(p, x) and xH(p, x). This one-bond cluster has two possible orientations, but only the horizontal case yields a non-vanishing residue [2], to leading order as q → 0. The spectral function, the density of resonances, and the mean number of resonances therefore behave as in agreement with the sum rules (2.49) and with eq. (3.15), to leading order in p. To second order in p, ρ(p, x) and H(p, x) consist of a countable infinity of discrete components (deltafunctions), corresponding to the resonances of configurations consisting of two bonds, in arbitrary relative position and orientation [2], and so on. The most salient subleading peaks have been marked in Figures 3 and 4 (dashed verticals) by some of the resonances of the two-bond and three-bond configurations shown in Figure 5. The configurations A, B, and E are self-dual, while (C,D) and (F,G) form dual pairs. The resonances of these clusters, determined exactly along the lines of ref. [2], are given in Table 1.
Right at the percolation threshold (p = p c = 1/2), the data for the spectral function ( Figure 4e) agree with the exact analytical result (2.36), which coincides with the EMA prediction. The accuracy of this agreement provides another check of the quality of the numerical simulations. The density of resonances remains a non-trivial function ρ(p c , x) at the percolation threshold (Figure 3e). Integrating the data of this Figure leads to the estimate ρ R (p c ) ≈ 0.80, again in good agreement with the exact result (3.16). It is worth noticing that the data shown in Figures 3e and 4e are practically not affected by the critical singularities at the endpoints of the spectrum. Indeed, the size of the critical region, given by eq. (2.46), can be estimated to be a few times 10 −3 , i.e., smaller than the width of the first or last bin.
As the concentration varies from p = 0 to p = p c = 1/2, the profiles of the density of resonances and of the spectral function deform in a progressive way. They get smoother and smoother, with their maxima moving in a continuous way. The dashed verticals are shown as guides for the eye on all plots of Figures 3 and 4, although they only label the most salient structures for a small enough concentration. The spectral function also gets progressively in better and better agreement with the EMA prediction.

Lifshitz tails
We have also investigated the behaviour of the density of resonances near the endpoints x = 0 and x = 1, in order to check the prediction (2.63) of the Lifshitz tail. We have chosen to investigate the density of resonances ρ(p, x), rather than the spectral function, because ρ(p, x) can be expected on general grounds to exhibit a clearer signal, in analogy with the one-dimensional situation [19]. Furthermore, the statistics can be doubled by using the symmetry relations (2.37), and each sample requires less CPU time, since the calculation of the eigenvectors X a is not required for ρ(p, x). against x −1/2 , for a concentration p = 0.1. The data correspond to over 103,000 samples of size N = 16 × 16, i.e., to 5 × 10 7 random bonds. The range of the plotted data corresponds to x ≤ λ min and λ max ≤ x ≤ 1, with λ min and λ max being the endpoints (2.52) of the EMA prediction for the spectral function. Indeed, the Lifshitz behaviour is expected to manifest itself mostly out of the "bulk" of the spectrum, the latter being conveniently defined as the support of the EMA formula [20]. A linear behaviour is clearly observed, confirming the analytical form (2.63) of the Lifshitz tail. A further qualitative confirmation of Hesselbo's argument on the Lifshitz behaviour is as follows. The data of Figure 6 exhibit oscillations around to the fitted straight line, and the top of each of the most prominent of these oscillations corresponds, with a good accuracy, to the lowest resonance of the hairpin structures shown in Figure 2a, embedded in the square lattice, with n = 1 to 4 cells. The case n = 1 corresponds to configuration G of Figure 5, with its lowest resonance λ G1 .
From a quantitative viewpoint, the slope of the fitted straight line in Figure 6 yields C(p = 0.1) ≈ 3. 19. The amplitude C(p) has been similarly measured for p = 0.05, 0.15, and 0.2. The results are plotted in Figure 7 against |ln p|. The data are nicely fit to the straight line C(p) ≈ 1.98|ln p| − 1.39, to be compared with the analytical prediction (2.64). The slope 1.98 is some 10% smaller than the analytical value π/ √ 2 = 2.2214. The intercept yields the estimate S ≈ 1.40 for the configurational entropy per cell, a significantly larger value than the entropy of self-avoiding walks on the square lattice, S SAW = ln µ SAW = 0.970 [24]. These observations suggest that other types of linear extended structures, besides the worm-like ones identified in the framework of Hesselbo's argument, may contribute to the Lifshitz behaviour of the conductivity.

Distribution of spacings between resonances
The distribution of spacings between successive energy levels has been extensively investigated in a variety of quantum systems, ranging from nuclei to billiards [25]. Generic spectra belong to three universality classes of level spacing distributions, according to their symmetry properties, in correspondence with the classical ensembles of real symmetric, Hermitian, and symplectic random matrices [26], respectively called GOE, GUE, and GSE.
We have investigated the distribution of spacings between successive resonances in the range where the spectral density of resonances is very flat, i.e., very close to being a constant, ρ = 0.658 (see Figure 3e). The range (4.3) can be considered as fully generic, although the concentration assumes its critical value p c , since the critical singularities only influence very small regions around the endpoints x = 0 or x = 1, as explained in section 4.1.
(4.4) Figure 8 shows a histogram plot of the distribution P (s) of the spacings s a , obtained from an ensemble of networks of size 20 × 20 having 10 7 random bonds in the range (4.3) in total. This distribution obeys by construction ∞ 0 P (s)ds = ∞ 0 sP (s)ds = 1. Since the generalised eigenvalue problem (3.10) involves two real symmetric matrices, ∆ and ∆ Q , the natural universality class to which the data for P (s) should be compared is the Gaussian orthogonal ensemble (GOE). The law P GOE (s) [26] is shown as a full line in Figure 8. The data share characteristic qualitative features of the GOE spacing distribution, including a linear repulsion at short spacings, i.e., P (s) ∼ s for s ≪ 1, and a fast fall-off at large spacings. There is, however, a small but significant quantitative difference between the data in the range (4.3) and the GOE prediction.

Number of visible resonances
We now turn to the number of visible resonances of the RL − C model on a finite network. A resonance is said to be visible if it corresponds to a true maximum in the frequency dependence of the real part of the admittance, as given by eq. (2.20a).
For a finite network of size N × N , the typical spacing between resonances scales as 1/N 2 , as already mentioned in section 4.3. Since two resonances of comparable strengths are resolved, i.e., separately visible, if their spacing is larger than the width 1/Q of each of them, we are led to propose the following finite-size scaling law for the fraction of visible resonances The scaling function F (X) is expected to grow monotonically from F (0) = 0, since a vanishing fraction of the resonances is visible if Q does not scale as N 2 , to F (∞) = 1, since all the resonances of a finite sample are eventually visible for a large enough quality factor. Figure 9 shows a plot of the fraction of visible resonances over the whole spectrum (0 ≤ λ ≤ 1), for p = p c = 1/2. The observed collapse of the data for the sizes 12 × 12 and 20 × 20 clearly demonstrates the validity of the scaling law (4.5). The full line shows a numerical fit of both series of data to the common analytical form 1/ 1 − F (X) 2 = 1 + XP 2 (X), with P 2 (X) being a quadratic polynomial. The quality of the fitted curve, meant as a guide to the eye, suggests a linear behaviour of the scaling function at small X, as well as a power-law convergence of the form The fraction of visible resonances depends a priori on how uniformly the resonances are distributed, and on the dispersion in the corresponding residues α a . Both features can depend quantitatively on the concentration p. The scaling function F (X) is therefore expected not to be universal, but rather to weakly depend on the concentration p and on the range of values of λ under consideration. Its main qualitative features, such as the power law (4.6), are however expected to be universal. The same remarks apply to the distribution of spacings between resonances, investigated in section 4.3.

Distribution of local electric fields
The algorithm presented in section 3 also gives access to the spatial structure of the resonances. Indeed, the eigenvector X a directly provides a map of the electric potentials at the resonance corresponding to the eigenvalue λ a . For each resonance, we define the local electric field on the bond (x, y) as (4.7) The electric fields are defined up to an overall multiplicative constant c. We choose this constant to be c = √ n B , so that the normalisation (3.19) of the eigenvectors implies (x,y) One could think of many ways of analysing the spatial structure of the electric field at resonance. We have chosen to investigate the distribution of the local field living on any given bond. More precisely, we have evaluated the successive moments of this distribution on square samples of size N × N , namely The normalisation (4.8) ensures the identities S 0 (N ) = S 2 (N ) = 1. Figure 10 shows a log-log plot of the first non-trivial moments, of index k = 1, 3, 4, 5, and 6, against the linear size N of the sample, from N = 6 to N = 24. Data are obtained in the range (4.3), with around 10 7 random bonds for each sample size. Power laws of the form are clearly apparent. This scaling behaviour, with a non-trivial dependence of the exponent x k on the index k, is a signature of multifractality [27,28]. Along the lines of the multifractal formalism, we introduce the generalised (Rényi) dimensions d k , such that The d k are expected to decrease from d 0 = 2, the dimension of the network, to d ∞ = 0.
In physical terms, multifractality implies that the patterns of resonant electric fields exhibit strong local fluctuations, rather similar to those observed in wavefunctions of the Anderson model, in the marginal two-dimensional situation [29]. In particular the resonances are neither localised nor extended, in the conventional sense of these expressions. Indeed, extended patterns of electric fields would correspond to d k = 2 for all k ≥ 0, while localised ones would have d k = 0 for k ≥ 1. We recall that a similar phenomenon of multifractality has been reported for the DC problem of the conductor-insulator mixture, right at the percolation threshold [30,31].
From a quantitative viewpoint, a more refined fit of the data to the power laws (4.10), including a relative correction in 1/N , yields more accurate estimates for the exponents x k and the associated dimensions d k , listed in Table 2. Figure 11 shows a plot of the d k against the index k. An approximate linear decay of the form with β ≈ 0.194, is observed over a fairly broad range (0 ≤ k ≤ 4). A similar linear behaviour has been predicted analytically for the two-dimensional Anderson model in the weak-disorder regime [29]. The linear law (4.12) corresponds to a log-normal (LN) distribution of the local fields. Indeed, let us set ℓ = ln |E|. The scaling law (4.10) then reads exp(kℓ) ∼ exp(Bk(k − 2)/2), with B = β ln N . Neglecting the k-dependence of the prefactor, this last expression corresponds to a log-normal distribution for |E|, i.e., to a Gaussian law for the logarithmic variable ℓ, namely (4.14) The actual probability density Π(ℓ) is shown in Figure 12, for networks of size 16 × 16 in the range (4.3). This very asymmetric distribution looks quite different from a lognormal law. In particular, it falls off as Π(ℓ) ∼ exp(−ℓ) as ℓ → −∞, yielding a finite value at |E| = 0 of the probability density P (|E|) = e ℓ Π(ℓ). Most of the dependence of the distribution on the sample size N takes place to the right of the plot, for large values of ℓ, close to the upper bound ℓ max = (ln n B )/2, where a very fast decay is observed.

DISCUSSION
We have investigated the AC conductivity of binary random networks of complex impedances, with emphasis on its analytic structure in the complex variable h or λ, and on the corresponding resonant behaviour. The present analysis is an extension to the general binary case of a previous work [2], devoted to the resonant response of a finite cluster, or set of clusters, embedded in an infinite homogeneous host lattice. Along the lines of refs. [1,2], the poles of the conductance are interpreted in terms of the resonances which show up in the AC conductivity of the RL − C model, and of the relaxation times in the transient response of the R − C model. We have proposed an efficient algorithm, which allows a determination of the rational λ-dependence of the conductance Y of a finite sample, in terms of its poles λ a and of the associated residues α a . A very similar formalism had been proposed long ago by Straley [6].
We want to underline again that the main advantage of the present approach is to give at once the analytic structure of the conductance in h or λ, for any finite sample. As far as a numerical investigation of the resonant response is concerned, this approach is therefore more suitable than the usual numerical methods, which can only yield the conductivity of the binary model for a fixed value of the ratio h, such as the transfer-matrix method [12], or the iterative algorithm based on the Y − ∆ transformation [32], used e.g. in ref. [7]. It is also worth noticing that our approach yields the full analytic structure of the conductance, including the static conductance of the Q-phase of any finite sample. Indeed the latter quantity is given by eq. (2.3), with the corresponding residue α 1 being given by the sum rule (2.5) in terms of all the other residues, corresponding to genuine resonances.
An extensive use of this algorithm, in the case of the square lattice, has allowed us to investigate in detail many aspects of the resonant dielectric response of the binary model. In general we have used 10 7 random bonds or more per measurement, a good enough statistics to obtain very accurate data. We have investigated the density of resonances ρ(p, x) and the spectral function H(p, x). This approach yields a better evaluation of the spectral function than the more direct method, consisting in measuring the imaginary part of the conductivity. The most salient structures have been labeled, at least for a small enough concentration p, by resonances of configurations made of one to three bonds, which can be determined exactly, along the lines of ref. [2]. The data for the spectral function have also been compared with the EMA prediction.
The Lifshitz behaviour of the density of resonances near the endpoints has been successfully characterised. A good qualitative agreement is found with the argument of Hesselbo [21], according to whom the analogues of the Lifshitz sphere are linear extended objects, such as hairpins. From a quantitative viewpoint, our data suggest that other classes of extended structures may also contribute to the Lifshitz tails, even for a small concentration p. The present situation is a lucky one, since numerical investigations of Lifshitz tails in more than one dimension are known to be a very difficult task in general, especially in the case of binary disorder [20]. The distribution of the spacings between neighbouring resonances has been found in qualitative agreement with the universal distribution of the GOE universality class of random matrices, although small but definite differences show up at a quantitative level. The number of visible resonances of the RL−C model on a finite sample of size N × N , as a function of the quality factor Q, has been shown to obey a finite-size scaling law (4.5), involving a scaling function F (X) of the variable X = Q/N 2 . Quantities such as the scaling function F (X), or the spacing distribution P (s), are expected not to be universal: quantitative features of these functions should rather weakly depend on the concentration p and on the range of resonances considered, the range (4.3) being meant as a generic example.
More generally, for the binary model on a d-dimensional lattice, the appropriate finitesize scaling variable describing dielectric resonances reads X = Q/N d . Indeed, the number of resonances on a sample of linear size N scales as its volume N d . This observation yields in particular a prediction for the divergence of the current correlation length ξ in the weak-dissipation regime (Q ≫ 1). By setting X ∼ 1 for N ∼ ξ, we get We thus recover a simple result due to Hesselbo [21], which has been corroborated by numerical simulations, yielding ν = 0.4 ± 0.1 in two dimensions [8].
Finally, we have investigated the distribution of the resonant electric field living on any given bond. The moments of this distribution obey power laws with non-trivial exponents x k , a characteristic feature of multifractality. The associated generalised dimensions d k are found to behave similarly to those observed in the Anderson model of localisation, in the marginal two-dimensional case [29]. Multifractality thus appears to be a quite generic feature of the resonant response of binary networks. In particular, this phenomenon is unrelated to the percolation transition. In analogy with the fraction of visible resonances or the spacing distribution P (s), the exponents x k and the dimensions d k are expected not to be fully universal, but to exhibit a weak dependence on the concentration p, and possibly on the range of resonances considered.
This multifractal picture provides a quantitative characterisation of local features of the fluctuations in electric fields at resonance observed previously [8]. These giant fluctuations have been argued to be responsible for surface-enhanced Raman scattering, this phenomenon being especially pronounced in strongly disordered semicontinuous films. There is a regime where the enhancement factor is predicted to be proportional to E 4 [9], whence the relevance of the dimension d 4 , in our notation. In a realistic system, with a small but finite dissipation rate 1/Q, these fluctuations are expected to be critical, i.e., to exhibit strong spatial correlations, on scales smaller than the current correlation length ξ, estimated in eq. (5.1). Since the algorithm used in this work gives direct access to the full map of electric fields at resonance, it could be used to investigate other aspects of dielectric resonances, including their spatial correlations, on which some information is already available [8]. Figure 1: Schema of the binary network under consideration, with M = N = 8. The conductance Y is measured between the plane electrodes (bus bars). P-bonds with conductance σ 0 are shown as solid lines, Q-bonds with conductance σ 1 are shown as dotted lines.   Table 1, of the configurations made of one to three bonds, shown in Figure 5.   Table  1, and used to label the most salient peaks in Figures 3 and 4, for p small enough.         Figure 5, and used to label the most salient peaks in Figures 3 and 4, for p small enough. Table 2: Exponents x k and associated generalised dimensions d k characterising the multifractal distribution of local electric fields at resonance, measured in the range (4.3). Table 1 configuration resonances A λ A1 = 1/2 = 0.500000 B λ B1 = 1/π = 0.318310 λ B2 = 1 − 1/π = 0.681690 C λ C1 = 1 − 2/π = 0.363380 λ C2 = 2/π = 0.636620