On Effective Conductivity on Zd Lattice Leonid G. Fel and Konstantin M. Khanin † ‡ 3 (Dedicated to D. Ruelle and Ya. G. Sinai on occasion of their 65th birthday) 0 0 2 n School of Physics and Astronomy of Exact Sciences a † J 7 Tel Aviv University, Ramat Aviv 69978, Israel 3 v and 8 4 Isaac Newton Institute for Mathematical Sciences, 0 ‡ 1 1 University of Cambridge, 1 0 / 20 Clarkson Road, Cambridge CB3 OEH, UK h p - Heriot-Watt University, Edinburgh h t a Landau Institute, Moscow m : v February 7, 2008 i X r a Abstract We study the effective conductivity σ for a random wire problem on the d- e dimensional cubic lattice Zd, d 2 in the case when random conductivities on bonds ≥ are independent identically distributed random variables. We give exact expressions for the expansion of the effective conductivity in terms of the moments of the disorder parameter up to the 5th order. In the 2D case using the duality symmetry we also derive the 6th order expansion. We compare our results with the Bruggeman approx- imation and show that in the 2D case it coincides with the exact solution up to the terms of 4th order but deviates from it for the higher order terms. Key words: effective conductivity, Bruggeman’s equation 1 1 Introduction. The problem of conductivity of the random composite medium and the equivalent problem of diffusion in a symmetric (self-adjoint) random environment has been a subject of intensive study for the last 25 years. It is virtually impossible to give a full reference list and we just mention few papers where the mathematical aspects of the theory were considered for the first time: [12],[15], [16], [1]. In the mathematical literature this problem usually is quoted as the problem of homogenization for the second order elliptic differential operators with random coefficients. Roughly speaking the main result can be formulated in the following way: there exists a non-random effective conductivity tensor or effective diffusion matrix such that the asymptotic properties of the system are the same as for a homogeneous system governed by the effective parameters. The subject is a very active research area till now with a vast number of papers publishing every year. However there are very few results related to the problem of calculation of effective conductivity and diffusion matrix. In addition to the trivial one-dimensional case such results are known only in the self-dual situation in dimension two (Keller-Dykhne duality) and in the case of two-component systems where the analytic continuation method is used to express the effective conductivity as an analytic function of theratio of the conductivities of two components (see [2], [14], [8], [4], [5]). Inthis paper we discuss a very general rigorous method in the lattice case which was developed in [1]. The method is based on a convergent power series expansion for the effective parameters and can be applied for arbitrary probability distribution of random conductivities. However, the combinatorics of this expansion is rather complicated. That is a reason why it was not used for concrete calculations in the past. The present paper has two main goals. First of all we demonstrate the constructive potential of the method in [1] and give exact formulae for the first 5 orders of the expan- sion for the effective conductivity in arbitrary dimension. In the 2D case we also calculate the 6th order terms. We then use our exact results to study the quality of the classical Bruggeman approximation. We show that in the 2D case the Bruggeman approximation is extremely accurate and coincides with the exact answer up to the terms of the 4th order. We assume everywhere that the random conductivities (jump rates) are independent iden- tically distributed random variables. Although we consider only the case of Zd lattice we strongly believe that the method can be generalized for other types of lattices and even for 2 the continuous situation. Yakov Sinai was a teacher of one of us and it is our pleasure to dedicate this paper to his 65th birthday. In fact one of the motivations for this paper was to illuminate the method developed together with Yakov Grigorevich and to demonstrate its effective power. 2 Effective conductivity on Zd Lattice. 2.1 Exact expansion for effective conductivity. We consider effective conductivity for a random wire problem on the d-dimensional cubic lattice Zd, d 2. Throughout the paper we assume that bond conductivities σ are indepen- ≥ dent identically distributed positive random variables. We are not making any assumptions on a probability distribution of σ which can be either discrete or continuous. As we have mentioned above the calculation of the effective conductivity is equivalent to the calculation of the effective diffusion matrix for the continuous time random walk in random environ- ment. In this case random conductivities should be understood as jump rates through the corresponding bond. We shall use the formula for the effective diffusion matrix M which e was obtained in [1]. This formula is given by a convergent series where the role of small parameter is played by a deviation of a random variable σ from its average value σ . Since h i we consider transitions only along the bonds of Zd lattice with i.i.d. transition rates σ, the effective diffusion matrix is a scalar matrix: M = 2σ I, where effective diffusion coefficient e e (or effective conductivity) σ can be expressed in terms of a convergent power series. We e first introduce the necessary notations. A path γ = (z ,α ),(z ,α ),...,(z ,α ) is a finite sequence of pairs (z,α) where z is a 1 1 2 2 k k { } point of lattice Zd and α = 1,2,...,d corresponds to one of the d possible directions. Notice thatz ,z arenot necessarily neighbours onthe lattice. The sum of two pathsγ = γ +γ is i i+1 1 2 simply the ordered union of two sequences where the pairs of the second path follow the pairs of the first one. With each pair (z,α) we associate a random variable σ (z) = σ(z,z +e ), α α where e is a unit vector in the direction α and σ(z,z + e ) is the random transition rate α α (conductivity) along the bond (z,z +e ). Denote by u (z) = σα(z)−hσi and define for each α α hσi 3 path γ = (z ,α ),(z ,α ),...,(z ,α ) the moment 1 1 2 2 k k { } k γ = u (z ) . (1) h i h αi i i i=1 Y A convergent expansion below for the effective conductivity is expressed through the mo- ments of a random variable u. We shall also need the following cumulant of a path γ: k m E(γ) = ( 1)m−1 γ , (2) j − h i mX=1 γ1+·X··+γm=γYj=1 where summation in (2) is taken over all possible partitions of the path γ into a sum of paths γ . Finally we define a kernel Γ (z): j αβ 1 1 sinπλ sinπλ cos2π((λ,z) 1λ + 1λ ) d Γ (z) = ... α β − 2 α 2 β dλ , (3) αβ − d sin2πλ γ Z0 Z0 γ=1 γ γ=1 Y P where λ = (λ ,...,λ ). Notice that Γ (0) = 1 and Γ (z) = Γ ( z). We can now write 1 d αα −d αβ βα − the following exact formula for σ : e ∞ σ = σ 1+ A(k) , (4) e h i ! k=2 X where k−1 A(k) = E(γ) Γ (z z ). (5) αiαi+1 i+1 − i γ={(z1,α1),.X..,(zk,αk)}∈G1(k) Yi=1 (k) Here is the set of all possible paths γ = (z ,α ),...,(z ,α ) such that z = 0 and G1 { 1 1 k k } 1 α = α = 1. It has been proven in [1] that the infinite sum in (5) is absolutely convergent. 1 d That is due to the fact that for the paths γ which might lead to divergence of A(k) one has E(γ) = 0. It was also shown that the expansion in (4) is absolutely convergent and gives an exact value of σ provided u u < 1/2. The last condition is technical and probably can e 0 | | ≤ be improved. In the following proposition we rewrite (4), (5) in a slightly different way. Proposition 1 ([1]). Assume that there exists a constant u < 1 such that u u with probability 1. Then for 0 2 | | ≤ 0 any dimension d [k] ∞ 2 σ = σ 1+ a(d) us1 ... usm , (6) e h i s1,...,smh i h i  Xk=2 mX=1 s1s1+,·.X·.·.+,ssmm≥=2k    4 (d) where the constants a depend only on dimension d and [ ] denotes the integer part. s1,...,sm · Moreover, for any n 1 the following estimate holds ≥ [k] n 2 (2u )n+1 σ σ 1+ a(d) us1 ... usm 0 . (7) (cid:12)(cid:12) e −h i s1,...,smh i h i(cid:12)(cid:12) ≤ 1 2u0 (cid:12)(cid:12)  Xk=2 mX=1 ss11+,·.X·.·.+,ssmm≥=2k (cid:12)(cid:12) − (cid:12)  (cid:12) (cid:12) (cid:12) Note that the series in (6) is absolutely convergent. (cid:12) (cid:12) 2.2 The 4th order expansion. It is easy to see that only those paths for which each pair (z,α) is present at least twice give nonzero contribution to (5). This immediately implies that u2 u3 A(2) = u2 Γ (0) = h i , A(3) = u3 Γ2 (0) = h i. (8) h i 11 − d h i 11 d2 Hence the 3-rd order approximation to σ is given by e u2 u3 σ(3) = σ 1 h i + h i . (9) e h i − d d2 (cid:18) (cid:19) In the 4-th order the combinatorics is slightly more complicated. Indeed, nonzero contribu- tions correspond to the paths γ(4) = (0,1),(0,1),(0,1),(0,1) , γ1 (4) = (0,1),(z,1),(z,1),(0,1) ,z = 0 , { } 1,z { } 6 γ (4) = (0,1),(z,α),(z,α),(0,1) ,α= 1 , γ2 (4) = (0,1),(z,1),(0,1),(z,1) ,z = 0 . α,z { } 6 1,z { } 6 Another possible type of paths γ3 (4) = (0,1),(0,1),(z,1),(z,1) ,z = 0 gives zero contri- 1,z { } 6 bution since E(γ3 (4)) = 0. Easy calculation gives 1,z A(4) = ( u4 u2 2)Γ3 (0) + u2 2Γ (0) Γ2 (z) Γ2 (0) h i−h i 11 h i 11 11 − 11 " !# (cid:20) (cid:21) zX∈Zd d + u2 2 Γ3 (z) Γ3 (0) + u2 2Γ (0) Γ2 (z) . (10) h i 11 − 11 h i αα 1α " !# " # zX∈Zd Xα=2 zX∈Zd Notice that 1 1 sin2πλ sin2πλ d Γ2 (z) = ... β α dλ . (11) βα ( d sin2πλ )2 γ z∈Zd Z0 Z0 γ=1 γ γ=1 X Y Hence P d Γ2 (z) = 1. (12) βα β,α=1z∈Zd X X 5 Since d Γ2 (z) (13) βα α=1z∈Zd X X does not depend on β we get d 1 Γ2 (z) = . (14) 1α d α=1z∈Zd X X Using (10, 14) we obtain 1 d 2 A(4) = u4 − u2 2 + u2 2 Γ3 (z). (15) −d3h i− d3 h i h i 11 z6=0 X The third term in (15) vanishes in the 2D case. Indeed, if z = (x,y) we have Γ (x,y) = 11 Γ (y,x). Obviously Γ (y,x) + Γ (y,x) = 0 if (y,x) = (0,0). Hence, for nonzero (x,y) 22 11 22 6 we have Γ (y,x) = Γ (x,y) which immediately implies Γ3 (z) = 0. As a result we 11 − 11 z6=0 11 obtain the 4-th order approximation for d = 2: P 1 1 1 σ(4) = σ 1 u2 + u3 u4 . (16) e h i − 2h i 4h i− 8h i (cid:18) (cid:19) We next demonstrate that for d 3 ≥ 1 Γ3 (z) = (17) 11 6 −d3 zX∈Zd which implies Γ3 (z) = 0. (18) 11 6 z6=0 X Denote H(d) = d3 Γ3 (z). Using simple Fourier analysis we have − z∈Zd 11 P 1 1 d d H(d) = ... H(λ,µ) dλ dµ , (19) γ γ Z0 Z0 γ=1 γ=1 Y Y where sin2(π(λ +µ )) sin2(πλ ) sin2(πµ ) 1 1 1 1 H(λ,µ) = . (20) 1 d sin2(π(λ +µ )) 1 d sin2(πλ ) 1 d sin2(πµ ) d γ=1 γ γ d γ=1 γ d γ=1 γ As we have explainedPabove the symmetry in tPhe 2D case givesP Γ3 (z) = 0 which is z6=0 11 equivalent to H(2) = 1. We conjecture that H(d) is a strictly decrPeasing function of d. The conjecture implies that Γ3 (z) > 0 for all d 3. Although the conjecture above was z6=0 11 ≥ not proven rigorously wePhave checked it numerically for 3 d 5: ≤ ≤ H(3) = 0.923, H(4) = 0.874, H(5) = 0.846. (21) 6 Finally, we get the following 4-th order approximation in an arbitrary dimension: 1 1 1 d+H(d) 3 σ(4) = σ 1 u2 + u3 u4 − u2 2 . (22) e h i − dh i d2h i− d3h i− d3 h i (cid:18) (cid:19) 2.3 The 5th order expansion. We proceed with the 5-thorder calculations. The following paths give nonzero contributions: γ(5) = (0,1),(0,1),(0,1),(0,1),(0,1) , γ1 (5) = (0,1),(z,α),(z,α),(z,α),(0,1) { } α,z { } γ2 (5) = (0,1),(0,1),(z,α),(z,α),(0,1) , γ3 (5) = (0,1),(z,α),(0,1),(z,α),(0,1) α,z { } α,z { } γ4 (5) = (0,1),(z,α),(z,α),(0,1),(0,1) , γ˜1 (5) = (0,1),(z,1),(0,1),(z,1),(z,1) α,z { } 1,z { } γ˜2 (5) = (0,1),(z,1),(z,1),(0,1),(z,1) , γ˜3 (5) = (0,1),(z,1),(0,1),(0,1),(z,1) 1,z { } 1,z { } γ˜4 (5) = (0,1),(0,1),(z,1),(0,1),(z,1) . 1,z { } Notice that in the case α = 1 the summation in the paths γs (5),γ˜s (5),1 s 4 is 1,z 1,z ≤ ≤ performed over all z = 0. Using (5) we get 6 1 A(5) = u5 +K (d) u2 u3 , (23) d4h i 5 h ih i where d 3 6 4 K (d) = Γ2 (z)+ Γ4 (z) Γ3 (z). (24) 5 d2 1α 1α − d4 − d 11 ! α=1 z∈Zd z∈Zd z6=0 X X X X This together with (14) gives d 3(d 2) 4 K (d) = − + Γ4 (z) Γ3 (z). (25) 5 d4 1α − d 11 Xα=1zX∈Zd Xz6=0 In the 2D case both the first and the last term in (25) vanish and K (2) = Γ4 (z)+ Γ4 (z) = I +I , (26) 5 11 12 1 2 z∈Z2 z∈Z2 X X where 1 1 I = h2(λ ,λ )dλ dλ , 1 1 1 2 1 2 Z0 Z0 1 1 sin2π(λ µ )sin2πµ dµ dµ 1 1 1 1 2 h (λ ,λ ) = − (27) 1 1 2 Z0 Z0 sin2π(λ µ )+sin2π(λ µ ) sin2πµ +sin2πµ 1 1 2 2 1 2 − − (cid:18) (cid:19)(cid:18) (cid:19) 7 and 1 1 I = h2(λ ,λ )dλ dλ , 2 2 1 2 1 2 Z0 Z0 1 1 sinπ(λ µ )sinπ(λ µ )sinπµ sinπµ dµ dµ 1 1 2 2 1 2 1 2 h (λ ,λ ) = − − .(28) 2 1 2 Z0 Z0 sin2π(λ µ )+sin2π(λ µ ) sin2πµ +sin2πµ 1 1 2 2 1 2 − − (cid:18) (cid:19)(cid:18) (cid:19) The values of I ,I were found numerically: I = 0.06391,I = 0.00439. As a result we get 1 2 1 2 in the 2D case the following 5-th order expansion: 1 1 1 1 σ(5) = σ 1 u2 + u3 u4 + u5 +I u2 u3 , (29) e h i − 2h i 4h i− 8h i 16h i h ih i (cid:18) (cid:19) where I = I +I = 0.0683. 1 2 In the general case d 3 we have ≥ d d Γ4 (z) = Γ4 (z)+ Γ4 (z) = I (d)+(d 1)I (d), (30) 1α 11 1α 1 − 2 α=1z∈Zd z∈Zd α=2z∈Zd X X X X X where 1 1 d I (d) = ... h2(λ) dλ , 1 1 γ Z0 Z0 γ=1 Y 1 1 sin2π(λ µ )sin2πµ d dµ h (λ) = ... 1 − 1 1 γ=1 γ (31) 1 Z0 Z0 d sin2(π(λ µ )) Qd sin2(πµ ) γ=1 γ − γ γ=1 γ (cid:18) (cid:19)(cid:18) (cid:19) P P and 1 1 d I (d) = ... h2(λ) dλ , 2 2 γ Z0 Z0 γ=1 Y 1 1 sinπ(λ µ )sinπ(λ µ )sinπµ sinπµ d dµ h (λ) = ... 1 − 1 2 − 2 1 2 γ=1 γ . (32) 2 Z0 Z0 d sin2(π(λ µ )) d sin2(πQµ ) γ=1 γ − γ γ=1 γ (cid:18) (cid:19)(cid:18) (cid:19) P P Collecting all the terms we get 1 1 1 d+H(d) 3 σ(5) = σ 1 u2 + u3 u4 − u2 2 (33) e h i − dh i d2h i− d3h i− d3 h i 1 3d+d4I(d)+4H(d) 10 + u5 + − u2 u3 , d4h i d4 h ih i ! where I(d) = I (d)+(d 1)I (d) and H(d) is given by (19), (20). 1 2 − 8 2.4 Keller-Dykhne duality and the 6th order expansion in the 2D case. Althoughitispossibleinprincipletocalculateanexpansionofanarbitraryordertheproblem becomes more and more cumbersome for higher order terms. However in the 2D case one can significantly simplify calculations using the duality symmetry which was discovered by Keller ([9]) and Dykhne ([6]). Consider duality transformation 1 σ . (34) → σ Denote by σ , σ−1 the probability distributions for positive random variables σ and σ−1 { } { } respectively. Then duality symmetry which holds only in the 2D case implies that σ ( σ−1 ) = σ−1( σ ). (35) e { } e { } Although both Keller and Dykhne considered only the continuous systems the symmetry (35) can be extended to the case of discrete lattice systems which we study in this paper (see [11]). The duality symmetry immediately implies that in the self-dual case, i.e. when the probability distributions σ and σ−1 coincide, the effective conductivity σ = 1. e { } { } It also gives an exact answer in the case which we call almost self-dual. We say that the probabilitydistributionforarandomvariableσ isalmost self-dualwithrespect totheduality transformation(34)ifthereexistsapositiveconstantσ suchthattheprobabilitydistribution 0 for σ σ is exactly self-dual, i.e. 0 σ σ = (σ σ)−1 . (36) 0 0 { } { } Since σ is a homogeneous function of the first order and σ ( σ σ ) = 1, it follows that e e 0 { } in the almost self-dual situation σ ( σ ) = σ−1. Notice that in the two-component case e { } 0 with equipartition, i.e. when σ takes values σ and σ with probabilities 1 the probability 1 2 2 distribution for σ is almost self-dual with σ = (√σ σ )−1. Hence, 0 1 2 σ = σ−1 = √σ σ . (37) e 0 1 2 This well-known result by Keller and Dykhne provides one of the very few exact solutions for the effective conductivity. We next show that the duality symmetry alone gives a lot of relations on the coefficients of the expansion (6). In fact we shall be able to recover the 6th order expansion using only 9 the 5th order and the symmetry. Consider the case when σ takes three values: 1 ǫ with − probability p, 1 αǫ with probability p and 1 with probability 1 2p. Correspondingly a − − random variable σ−1 takes values 1 and 1 with probabilities p and 1 with probability 1−ǫ 1−αǫ 1 2p. We shall use the formula (6) in order to calculate σ ( σ )σ ( σ−1 ) and check the e e − { } { } duality identity (35) subsequently in the 2nd, 4th, 6th and 8th orders of the power series expansion in ǫ. This inductive procedure allows to find all the relations on the coefficients (2) a . We performed calculations using the Maple symbolic package. In the 2nd order one s1,...,sm immediately gets a(2) = 1. The 4th order calculations give two relations: 2 −2 3 3 1 3 (2) (2) (2) (2) a = a , a = a . (38) 2,2 2 3 − 8 4 4 − 2 3 The 6th order expansion provides four more relations: 7 3 15 1 1 a(2) = a(2) + a(2) , a(2) = + (a(2))2 2a(2) a(2), 2,2,2 2 3 2 2,3 − 16 3,3 2 2 3 − 3 − 2,3 11 3 5 5 5 1 (2) (2) (2) (2) (2) (2) (2) a = 6a a + a , a = a a . (39) 2,4 8 − 3 − 2 2,3 2 5 6 2 3 − 2 5 − 2 Using (29) we have 1 1 (2) (2) (2) a = , a = , a = I = 0.0683 (40) 3 4 5 16 2,3 which immediately gives a(2) = 0, a(2) = 1 and 2,2 4 −8 3 1 1 1 3 1 (2) (2) (2) (2) a = I , a = I , a = I , a = . (41) 2,2,2 2 − 16 3,3 32 − 2,4 32 − 2 6 −32 As a result we obtain the 6th order expansion in the 2D case: 1 1 1 1 σ(6) = σ 1 u2 + u3 u4 + u5 +I u2 u3 e h i − 2h i 4h i− 8h i 16h i h ih i 1 3 1 1 u6 I u2 u4 I u3 2 (42) −32h i− 2 − 32 h ih i− − 32 h i (cid:18) (cid:19) (cid:18) (cid:19) 3 1 + I u2 3 . 2 − 16 h i ! (cid:18) (cid:19) 3 The Bruggeman Approximation. 3.1 Bruggeman’s equation. The Effective Medium Approximation (EMA) was invented by Bruggeman [3], and has re- mained one of the most popular approximations used for calculations of the linear bulk 10