# A variance reduction strategy for numerical random homogenization based on the equivalent inclusion method

Sébastien Brisard, Michael Bertin, Frederic Legoll

<figure><img src="image_1.png" /></figure>

## To cite this version:

Sébastien Brisard, Michael Bertin, Frederic Legoll. A variance reduction strategy for numerical random homogenization based on the equivalent inclusion method. Computer Methods in Applied Mechanics and Engineering, 2023, 417, pp.116389. 10.1016/j.cma.2023.116389. hal-04058543

HAL Id: hal-04058543

https://hal.science/hal-04058543v1

Submitted on 1 Oct 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A variance reduction strategy for numerical random homogenization based on the equivalent inclusion method

Sébastien Brisard<sup>a,*</sup>, Michaël Bertin<sup>a,b,c</sup>, Frédéric Legoll<sup>a,c</sup>

<sup>a</sup>Navier, Ecole des Ponts, Univ Gustave Eiffel, CNRS, Marne-la-Vallée, France

<sup>b</sup>CERMICS, Ecole des Ponts, Marne-la-Vallée, France

<sup>c</sup>MATHERIALS project-team, Inria, Paris, France

## Abstract

Using the equivalent inclusion method (a method strongly related to the Hashin-Shtrikman variational principle) as a surrogate model, we propose a variance reduction strategy for the numerical homogenization of random composites made of "inclusions" (or rather inhomogeneities) embedded in a homogeneous matrix. The efficiency of this strategy is demonstrated within the framework of two-dimensional, linear conductivity. Significant computational gains vs full-field simulations are obtained even for high contrast values. We also show that our strategy allows to investigate the influence of parameters of the microstructure on the macroscopic response. A specific attention is given to the computational cost of the surrogate model. In particular, an inexpensive approximation of the so-called influence tensors (that are used to compute the surrogate model) is proposed.

**Keywords:** Homogenization, Full-field simulations, Conductivity, Variance reduction, Surrogate model

This is the accepted version of the following article: "A variance reduction strategy for numerical random homogenization based on the equivalent inclusion method", which was sent to the publisher on 27 July 2023 and was subsequently published at https://doi.org/10.1016/j.cma.2023.116389.

This manuscript version is made available under the CC-BY-NC-ND 4.0 license. See Elsevier Sharing Policy.

## 1. Introduction

Most natural and manufactured materials are heterogeneous, and the prediction of their effective properties must account for this small-scale heterogeneity by means of an appropriate upscaling strategy. When the classical mean-field/effective-field techniques (Hashin-Shtrikman bounds, Mori-Tanaka bounds, Self-Consistent scheme, see e.g. the review articles [1, 2]) fail to deliver a prediction that is accurate enough, expensive full-field numerical simulations must be performed. These simulations compute the solution to the so-called *corrector problem* within an “elementary cell” (to be defined below), under specific local constitutive laws. The elementary cell is subjected to an elementary macroscopic loading that ensures that the Hill–Mandel lemma holds.

The extent of the elementary cell depends on the geometric structure of the heterogeneities. For *periodic homogenization*, the elementary cell is well-defined and corresponds to the unit-cell [3, Chap. 1], [4, 5, 6, 7]. For *random homogenization* (the topic of the present article), the elementary cell should be a realization of the random microstructure over the whole space [8, 5], which is of course not practical. The computation is therefore carried out on a large, but finite-size realization of the random material, which we call here a *statistical volume element* (SVE), following Ostoja-Starzewski [9]. The solution to the corrector problem on the SVE delivers the *apparent properties*. In the limit of asymptotically large SVEs, the apparent properties converge almost surely to the true *effective (homogenized) properties*.

<footnote>*Corresponding author.</footnote>

<footnote>Email addresses: sebastien.brisard@univ-eiffel.fr (Sébastien Brisard), frederic.legoll@enpc.fr (Frédéric Legoll)</footnote>
The apparent properties thus defined appear as a size-dependent random variable. The above-discussed workflow introduces two types of errors: (i) a deterministic bias and (ii) random errors. The first type of error relates to the finite size of the SVE: the average of the apparent properties over all possible microstructures in the SVE is different from the exact effective, homogenized properties. The second type of error relates to the realization itself: the microstructure in the SVE being random, its associated apparent properties are random. In this article, because this is overwhelmingly the case in the numerical practice, we consider corrector problems in the SVE complemented by *periodic boundary conditions*. This choice is motivated by considerations about the decrease of the systematic error (i.e. the deterministic bias) when the SVE size increases. It has been numerically observed (see e.g. [10]) and theoretically demonstrated in some cases (see [11] and the works of A. Gloria, F. Otto and their collaborators initiated in [12]) that the statistical error decays with a slower rate (with respect to the SVE size) than the systematic error. For large SVEs, the statistical error is therefore dominating the systematic error. In what follows, we focus on reducing the statistical error, and describe an approach to more efficiently compute the *average apparent properties*.

Estimating the average apparent properties naturally requires the Monte-Carlo method, which is known to converge slowly. Indeed, the central limit theorem delivers an estimate of the confidence interval that scales as $M^{-1/2}$, where M is the number of realizations (the number of full-field simulations in the present case). In other words, dividing the amplitude of the confidence interval by 10 requires to multiply the number of realizations by 100. This can lead to costly simulations.

Observing that the amplitude of the confidence interval is also proportional to the standard deviation of the random variable under consideration, variance reduction techniques aim at refining the estimate of the statistical error by controlling this standard deviation rather than the sample size. Many types of such approaches have been proposed in the literature (antithetic variables, control variate, importance sampling, stratification based approaches, to name but a few), in various contexts where the Monte-Carlo method is used (we refer to [13] for a review of these approaches). We mention here that we have already adapted some of these methods to the present context of numerical random homogenization, for linear and for some nonlinear problems, in [14, 15, 16, 17, 18, 19] (see also the review articles [20, 21]).

In this article, we focus on the class of control variate approaches, which are based on the use of a so-called *surrogate model*. This new random variable should be strongly correlated to the random variable of interest and its expectation should be known exactly (alternatively, evaluation of the surrogate model should be inexpensive enough to allow for a very accurate Monte-Carlo estimation of its expectation). The Monte-Carlo method is then run on a linear combination of the random variable of interest with the surrogate model, the coefficients of the linear combination being chosen such that its standard deviation is smaller than that of the random variable of interest.

Within the framework of linear electric conductivity, we present here a control-variate-based variance reduction strategy for the determination of the effective conductivity of an assembly of circular inclusions embedded in a homogeneous matrix. We do not consider size-effects, which (as argued previously) are a separate issue. We therefore work on *fixed-size* SVEs and focus on the determination of the *average apparent conductivity*.

We use the equivalent inclusion method (EIM) to produce inexpensive estimates of the apparent conductivity. The EIM was first introduced within the framework of linear elasticity by Moschovidis and Mura [22]. It relies on the formulation of the corrector problem as an integral equation known as the Lippmann-Schwinger equation [23, 24, 25]. As suggested by Moschovidis and Mura [22], the EIM seeks an approximate solution to the Lippmann-Schwinger equation that is constant in each inclusion. Since their seminal work, various strategies (mostly relying on collocation) have been proposed to derive EIM equations. Following previous work [26], we use here a Galerkin approach in a periodic setting.

The EIM estimates of the apparent conductivity then serve as a surrogate model which is used in conjunction with full-field simulations performed using the finite element method (FEM). We show that this surrogate model is indeed effective at reducing the amplitude of the confidence interval, even in the difficult cases of very large contrast between the matrix and the inclusions. In order for the surrogate model to also be cost-effective, an inexpensive approximation of the so-called influence tensors is proposed. This approximation sheds new light on the recent work by Zecevic and Lebensohn [27], where a similar approximation is considered.

We underline here the strong mechanical content of the EIM, which is a surrogate model well-tailored to our application, namely the homogenization of random assemblies of inclusions embedded in a homogeneous medium. We also emphasize the strong link between EIM, the Hashin-Shtrikman variational principle and the Hashin-Shtrikman
bounds. This link is further discussed in Sec. 3.2.2.

The article is organized as follows. Variance reduction techniques based on the use of a surrogate model (a.k.a. control variate techniques) are briefly reviewed in Sec. 2. Then, Section 3 provides an overview of the equivalent inclusion method, in particular its variational form. Our variance reduction strategy is outlined in Sec. 4 where some specificities of the EIM are underlined. The efficiency of the proposed strategy is demonstrated by numerical experiments in a two-dimensional setting in Sec. 5. We collect some conclusions in Sec. 6.

Note that all data necessary to reproduce the results presented in Sec. 5 and Appendix F are available at the Recherche Data Gouv open data repository [28, 29].

## 2. Variance reduction through a surrogate model

In the present section, we give a brief overview of the control variate variance reduction technique. Our goal is to compute a numerical estimate of the ensemble average $E(X)$ of a scalar random variable $X(\omega)$. To do so, the Monte-Carlo approach consists in drawing $M$ independent realizations of $X(\omega)$, denoted $X_1(\omega), \dots, X_M(\omega)$, and evaluating the empirical average

$$ \mu_M(X) = \frac{1}{M} \sum_{m=1}^{M} X_m(\omega). $$

The law of large numbers says that the empirical average $\mu_M(X)$ converges, when $M \to +\infty$, to the exact ensemble average $E(X)$. The empirical average $\mu_M(X)$ is unbiased, in the sense that, for any $M$, its expectation is equal to $E(X)$. The only error between $\mu_M(X)$ and $E(X)$ is therefore a statistical error, which is estimated by the central limit theorem:

$$ \lim_{M \to +\infty} \mathrm{Pr} \left[ \frac{1}{M} \sum_{m=1}^{M} X_m(\omega) - E(X) \le \alpha \frac{\sigma_X}{M^{1/2}} \right] = \Phi(\alpha), $$

where $\Phi$ is the standard normal cumulative distribution function and $\sigma_X^2$ is the variance of $X$ ($\sigma_X$ is the so-called standard deviation of $X$). The statistical error thus scales as $\alpha \sigma_X M^{-1/2}$. For example, for the 99% confidence interval, $\alpha = 2.6$: with a probability of 99%, the expectation $E(X)$ belongs to the interval

$$ \left[ \frac{1}{M} \sum_{m=1}^{M} X_m(\omega) - 2.6 \frac{\sigma_X}{M^{1/2}} ; \frac{1}{M} \sum_{m=1}^{M} X_m(\omega) + 2.6 \frac{\sigma_X}{M^{1/2}} \right]. $$

As expected, increasing the number $M$ of samples reduces the statistical error, albeit at the slow rate $M^{-1/2}$.

As pointed out above, several variance reduction strategies have been proposed in the literature with the aim to provide a more efficient approach than the Monte-Carlo method described above. The control variate approach [see e.g. 13, p. 277] is one of them. It is based on the observation that a reduction of the statistical error can be obtained by operating on the variance of $X$ itself. More precisely, let us consider a second random variable $Y(\omega)$, denoted below the surrogate model. A crucial assumption is that the ensemble average $E(Y)$ is known. We then consider the new random variable

$$ Z(\omega) = X(\omega) - \xi [Y(\omega) - E(Y)], $$

where $\xi$ is a free (deterministic) parameter. Since $Z$ and $X$ have the same ensemble average, we are left with evaluating the ensemble average of $Z$. To do so, we again compute an empirical average

$$ \mu_M(Z) = \frac{1}{M} \sum_{m=1}^{M} Z_m(\omega), $$

where the realizations $Z_m(\omega)$ are given by $Z_m(\omega) = X_m(\omega) - \xi [Y_m(\omega) - E(Y)]$. Using that $E(Z) = E(X)$, the central limit theorem yields

$$ \lim_{M \to +\infty} \mathrm{Pr} \left[ \frac{1}{M} \sum_{m=1}^{M} Z_m(\omega) - E(X) \le \alpha \frac{\sigma_Z}{M^{1/2}} \right] = \Phi(\alpha). $$
Note that the accuracy is now controlled by the standard deviation $\sigma_Z$ of the random variable Z. The variance of Z is given by

$$ \mathrm{var}(Z) = \mathrm{var}(X) - 2\xi \mathrm{cov}(X, Y) + \xi^2 \mathrm{var}(Y), $$

where $\mathrm{cov}(X, Y) = \mathbb{E}[(X - \mathbb{E}(X))(Y - \mathbb{E}(Y))]$ denotes the covariance of X and Y. Therefore, any choice of $\xi$ that makes $\xi^2 \mathrm{var}(Y) - 2\xi \mathrm{cov}(X, Y) \le 0$ effectively reduces the statistical error on the numerical estimation of $\mathbb{E}(X)$. The optimal choice $\xi^*$ of $\xi$ is the one which minimizes $\mathrm{var}(Z)$:

$$ \xi^* = \frac{\mathrm{cov}(X, Y)}{\mathrm{var}(Y)} \quad \text{and} \quad \mathrm{var}(Z^*) = [1 - \mathrm{corr}(X, Y)^2] \mathrm{var}(X), \qquad (1) $$

where $Z^*$ is the optimum value of Z and $\mathrm{corr}(X, Y) = \mathrm{cov}(X, Y)[\mathrm{var}(X)\mathrm{var}(Y)]^{-1/2}$ denotes the correlation of X and Y (by construction, $|\mathrm{corr}(X, Y)| \le 1$). The above equation shows that, if X and Y are strongly correlated, then the variance of $Z^*$ (hence the statistical error on the numerical estimate of its ensemble average) is much smaller than the variance of X.

In order to evaluate the ensemble average of X, the Monte-Carlo strategy is therefore applied to Z (with reduced statistical error). The quantity of interest $\mathbb{E}(X)$ is then obtained from the identity $\mathbb{E}(X) = \mathbb{E}(Z)$.

In practice, the quantities $\mathrm{cov}(X, Y)$ and $\mathrm{var}(Y)$, which are needed to find the optimal value of $\xi$ (see (1)), are not known. They are estimated numerically as follows:

$$ \mathrm{cov}(X, Y) \approx \frac{1}{M} \sum_{m=1}^{M} X_m(\omega) Y_m(\omega) - \frac{1}{M^2} \left[ \sum_{m=1}^{M} X_m(\omega) \right] \left[ \sum_{m=1}^{M} Y_m(\omega) \right] \quad \text{and} \quad \mathrm{var}(Y) \approx \frac{1}{M} \sum_{m=1}^{M} [Y_m(\omega)]^2 - \frac{1}{M^2} \left[ \sum_{m=1}^{M} Y_m(\omega) \right]^2. $$

In practice, the random variable Y is often cheap to evaluate in comparison to the evaluation of X (actual computational costs are provided in Sec. 5 when we consider our example of interest). The cost to compute the empirical average $\frac{1}{M} \sum_{m=1}^{M} X_m(\omega)$ is therefore identical to the cost to compute $\frac{1}{M} \sum_{m=1}^{M} Z_m(\omega)$. The latter estimate is a more accurate approximation of $\mathbb{E}(X)$ than the former: after variance reduction, and at equal computational cost, the statistical error on $\mathbb{E}(X)$ is indeed reduced by the factor

$$ \frac{\sigma_X}{\sigma_{Z^*}} = \sqrt{\frac{1}{1 - \mathrm{corr}(X, Y)^2}}, $$

defined as the efficiency index (higher is better).

**Remark 1.** In the case when the expectation of Y is not known, it is still possible to use the above control variate approach by replacing $\mathbb{E}(Y)$ by an accurate estimate of the form $\frac{1}{M} \sum_{m=1}^{M} Y_m(\omega)$, where the number M of realizations that are considered is sufficiently large such that the additional statistical error introduced by this approximation is much smaller than the statistical error between $\mu_M(X)$ and $\mathbb{E}(X)$. We typically take $M \gg M$. In the case when Y is inexpensive to evaluate, the additional cost is negligible.

### 3. Overview of the equivalent inclusion method for homogenization

In this section, we first present the electric conductivity problem of interest in this article, as well as the associated homogenization results (see Sec. 3.1). We next turn to the equivalent inclusion method (in Sec. 3.2), which is a strategy to approximate the apparent conductivity that we are going to use here as a surrogate model, as explained in Sec. 4.

#### 3.1. Effective and apparent conductivity

The homogenization problems considered in the present article are expressed within the framework of electric conductivity. We recall here this setting as well as the notion of effective and apparent conductivities.
We consider a random, heterogeneous material with local conductivity $\sigma_{\varepsilon}(\mathbf{x}, \omega)$ (second-order, symmetric, positive definite tensor) where $\mathbf{x} \in \mathbb{R}^d$ denotes the observation point, while $\omega$ refers to the realization. The subscript $\varepsilon$ highlights the fact that the characteristic scale $\varepsilon$ of variations of the conductivity is small. We consider the classical random ergodic homogenization setting, in which $\sigma_{\varepsilon}(\mathbf{x}, \omega) = \sigma(\mathbf{x}/\varepsilon, \omega)$ for a fixed stationary tensor $\sigma$, that we assume to be (uniformly in $\mathbf{x}$ and $\omega$) bounded from below and from above. Each realization of $\sigma(\cdot, \omega)$ fills the whole space $\mathbb{R}^d$.

The heterogeneous electric conductivity problem is given by

$$ \mathrm{div} \, \mathbf{J}_{\varepsilon} = f, \quad \mathbf{J}_{\varepsilon} = -\sigma_{\varepsilon} \, \mathrm{grad} \, \Phi_{\varepsilon}, $$

in some computational domain $\mathcal{D}$, where $f$ is the load (which is assumed independent of $\varepsilon$ and $\omega$), $\Phi_{\varepsilon}$ is the electric potential, $-\mathrm{grad} \, \Phi_{\varepsilon}$ is the electric field and $\mathbf{J}_{\varepsilon}$ is the electric current. The above problem is complemented with some boundary conditions that make it well-posed (typically Dirichlet, Neumann or periodic boundary conditions). It is then well-known (see e.g. the textbooks [8, 5, 30], the lectures [31, Part 7], the seminal works [32, 33] and also the review article [34] and the extensive bibliography contained therein) that the corresponding homogenized problem reads

$$ \mathrm{div} \, \mathbf{J}^{\mathrm{eff}} = f, \quad \mathbf{J}^{\mathrm{eff}} = -\sigma^{\mathrm{eff}} \, \mathrm{grad} \, \Phi^{\mathrm{eff}}, \qquad (2) $$

where $\sigma^{\mathrm{eff}}$ is the effective (here constant and deterministic) electric conductivity. This effective conductivity is defined in terms of the so-called corrector functions $\Phi$, which are solutions to the following problem: for any fixed constant macroscopic electric field $\bar{\mathbf{E}} \in \mathbb{R}^d$, find $\Phi$ such that

$$ \mathrm{div} \, \mathbf{J} = 0, \quad \mathbf{J} = \sigma \mathbf{E}, \quad \mathbf{E} = \bar{\mathbf{E}} - \mathrm{grad} \, \Phi \quad \text{in } \mathbb{R}^d, \quad \mathrm{grad} \, \Phi \text{ is stationary and } \mathbf{E}(\mathrm{grad} \, \Phi) = 0. \qquad (3) $$

In practice, solving (3) is not possible, since the problem is posed on the whole space $\mathbb{R}^d$. For practical computations, we therefore consider restrictions of this problem (and thus of the realizations of the conductivity $\sigma$) to the finite-size domain $Q = (-L/2, L/2)^d$ for some large $L$. Following Ostoja-Starzewski [9], we will call these restrictions statistical volume elements.

Owing to the periodic boundary conditions to be applied here, it will be convenient to define homologous points on opposite sides of the boundary $\partial Q$ as two points $\mathbf{x}^{\pm} \in \partial Q$ such that there exists $i \in \{1, \dots, d\}$ such that $\mathbf{x}_j^+ = \mathbf{x}_j^- + \delta_{ij} L$ for all $1 \le j \le d$ (see Fig. 1).

<figure><img src="image_3.png" /><figcaption>Figure 1: $(\mathbf{x}_1^+, \mathbf{x}_1^+)$ and $(\mathbf{x}_2^-, \mathbf{x}_2^+)$ are pairs of homologous points on the boundary of the two-dimensional domain $Q$.</figcaption></figure>

The apparent conductivity $\sigma^{\mathrm{app}}(\omega, L)$ of the SVE is defined through the solution to the following boundary-value problem on $Q$, denoted as the corrector problem, the localization problem or the auxiliary problem in the literature:

**Problem 1 (Corrector problem).** For any fixed constant macroscopic electric field $\bar{\mathbf{E}} \in \mathbb{R}^d$, find the $Q$-periodic fluctuations of the electric potential $\Phi^{\mathrm{per}}$ such that

$$ \mathrm{div} \, \mathbf{J} = 0 \quad \text{in } Q, \qquad (4) $$

$$ \mathbf{J} = \sigma \mathbf{E} \quad \text{in } Q, \qquad (5) $$

$$ \mathbf{E} = \bar{\mathbf{E}} - \mathrm{grad} \, \Phi^{\mathrm{per}} \quad \text{in } Q, \qquad (6) $$

$$ \mathbf{J}^{+} \cdot \mathbf{n}^{+} = -\mathbf{J}^{-} \cdot \mathbf{n}^{-} \quad \text{on } \partial Q. \qquad (7) $$
Equation (7) expresses that the flux $\mathbf{J} \cdot \mathbf{n}$ must be **anti-periodic** at the boundary of $Q$: $\mathbf{x}^-$ and $\mathbf{x}^+$ being two homologous points on opposite faces of $\partial Q$, $\mathbf{n}^\pm$ denotes the outer unit normal to $Q$ at point $\mathbf{x}^\pm$.

**Remark 2.** The solution $\Phi^{\text{per}}$ to Problem 1 is naturally searched in $H_\text{per}^1(Q)$. We observe that $\mathbf{J} \in H(\text{div}, Q)$, and its normal trace $\mathbf{J} \cdot \mathbf{n}$ (which appears in (7)) is thus well defined (see Appendix A).

Owing to its linearity, the solution to Problem 1 depends linearly on the sole parameter, namely the constant macroscopic electric field $\bar{\mathbf{E}} \in \mathbb{R}^d$. In particular, the apparent conductivity $\sigma^{\text{app}}(\omega, L)$ is defined as the linear mapping between the macroscopic electric field $\bar{\mathbf{E}}$ and the volume average of the electric current $\langle \mathbf{J} \rangle$:

$$ \langle \mathbf{J} \rangle = \sigma^{\text{app}}(\omega, L) \bar{\mathbf{E}}, \qquad (8) $$

where angle brackets refer to volume averages over $Q$. The above relation also writes

$$ \frac{1}{|Q|} \int_Q \sigma(\bar{\mathbf{E}} - \mathbf{grad} \, \Phi^{\text{per}}) = \sigma^{\text{app}}(\omega, L) \bar{\mathbf{E}}. $$

It should be observed that (6) and the periodicity of $\Phi^{\text{per}}$ lead to $\langle \mathbf{E} \rangle = \bar{\mathbf{E}}$ and the equality (8) therefore also reads

$$ \langle \mathbf{J} \rangle = \sigma^{\text{app}}(\omega, L) \langle \mathbf{E} \rangle. $$

We emphasize that the apparent conductivity depends on the realization $\omega$ (which introduces a random error with respect to $\sigma^{\text{eff}}$) and on the size $L$ of the SVE (which introduces a bias since, in general, $\mathbb{E}[\sigma^{\text{app}}(\omega, L)] \neq \sigma^{\text{eff}}$; see e.g. [10]).

For statistically homogeneous, ergodic systems, the effective conductivity introduced in the homogenized problem (2) (which, we recall, is defined in terms of the correctors $\Phi$ solution to (3)) is the limit of the above apparent conductivity when $L$ grows to infinity [35, 36] (we also refer to [37, 38] for related studies):

$$ \sigma^{\text{eff}} = \lim_{L \to +\infty} \sigma^{\text{app}}(\omega, L) \quad \text{almost surely in } \omega. $$

We recall that the realization $\sigma(\cdot, \omega)$ fills the whole space: the apparent conductivity is therefore defined for all $L$, and the above limit is meaningful. It turns out (this is the manifestation of the law of large numbers) that this limit does not depend on the realization $\omega$ under consideration. However, in order to obtain a better approximation of $\sigma^{\text{eff}}$, the above formula is often applied as

$$ \sigma^{\text{eff}} = \lim_{L \to +\infty} \mathbb{E}[\sigma^{\text{app}}(\omega, L)]. $$

The above expectation is evaluated as an empirical average. In practice, the typical homogenization workflow is as follows:

1. generate $M$ realizations of the microstructure $\sigma_1(\cdot, \omega), \dots, \sigma_M(\cdot, \omega)$;

2. for a fixed value of $L$, solve Problem 1 for the $M$ realizations $\sigma_m(\cdot, \omega)$;

3. compute the apparent conductivities $\sigma_m^{\text{app}}(\omega, L)$ for all $1 \le m \le M$;

4. estimate the average apparent conductivity

$$ \mathbb{E}[\sigma^{\text{app}}(\omega, L)] \approx \frac{1}{M} \sum_{m=1}^{M} \sigma_{m}^{\text{app}}(\omega, L); $$

5. possibly repeat for larger values of $L$ (typically, $2L, 4L, \dots$), until convergence (according to a specified practical criterion) of the apparent conductivity is observed. The “converged” value is the estimate of the effective conductivity.

The above workflow entails three types of numerical errors. In Step 2, discretization errors can typically be controlled by the coarseness of the mesh. In Step 4, statistical errors are controlled through the number $M$ of realizations (see below), while finite-size errors can be reduced by increasing the size $L$ of the SVE. The present article is devoted to statistical errors: it will therefore be assumed that the discretization error is negligible and that the size $L$ of the SVE is fixed. The fixed domain $Q$ occupied by the SVE will be referred to as the unit-cell, and we will drop any references to its size $L$: we write $\sigma^{\text{app}}(\omega)$ rather than $\sigma^{\text{app}}(\omega, L)$ in the remainder of this article.
## 3.2. Overview of the equivalent inclusion method

The equivalent inclusion method (EIM) was first introduced by Moschovidis and Mura [22] to compute approximations of the effective response of an assembly of inhomogeneities embedded in a homogeneous matrix. In the present section, we give a brief overview of the variational form of this method, proposed in [26].

Note that, following Eshelby [39] as well as Moschovidis and Mura [22], we define an *inclusion* as a bounded domain with same conductivity as the surrounding matrix, subjected to a non-zero electric polarization (to be defined below) and an *inhomogeneity* as a bounded domain with homogeneous conductivity that does not match that of the surrounding matrix.

### 3.2.1. The Lippmann–Schwinger equation

The starting point of the EIM is the Lippmann–Schwinger equation, which is a reformulation of the corrector problem (4)–(7) [23, 24, 25]. This equation requires a so-called *reference material* with homogeneous conductivity $\sigma_0$ (symmetric, positive definite second-order tensor). The reference material being selected, we define the following auxiliary problem:

**Problem 2.** The electric polarization P being prescribed (as a Q-periodic vector field with square integrable components), find the Q-periodic electric potential $\Phi^{\text{per}}$ such that

$$ \operatorname{div} \mathbf{J} = 0 \quad \text{in } \mathcal{Q}, $$

$$ \mathbf{J} = \sigma_0 \mathbf{E} + \mathbf{P} \quad \text{in } \mathcal{Q}, $$

$$ \mathbf{E} = -\operatorname{grad} \Phi^{\text{per}} \quad \text{in } \mathcal{Q}, $$

$$ \mathbf{J}^{+} \cdot \mathbf{n}^{+} = -\mathbf{J}^{-} \cdot \mathbf{n}^{-} \quad \text{on } \partial\mathcal{Q}. $$

For any $\mathbf{P} \in (L^2(\mathcal{Q}))^d$, we observe that the solution $\mathbf{J}$ to Problem 2 belongs to $H(\operatorname{div}, \mathcal{Q})$, and its normal trace (which appears in the last equation above) is thus well defined (see Appendix A). Note that we do not request $\mathbf{P}$ to have a well-defined normal trace.

Problem 2 is linear with respect to the parameter $\mathbf{P}$. The Green operator $\Gamma_0$ is defined as the linear operator that maps the vector field $\mathbf{P}$ onto $-\mathbf{E}$, where $\mathbf{E}$ denotes the unique solution to Problem 2:

$$ \Gamma_0(\mathbf{P}) = -\mathbf{E} = \operatorname{grad} \Phi^{\text{per}}. \qquad (9) $$

In the particular case of isotropic homogeneous reference materials, the conductivity tensor $\sigma_0$ is proportional to the second-order identity tensor $\mathbf{I}$, that is $\sigma_0 = \sigma_0 \mathbf{I}$, and the Green operator $\Gamma_0$ has a closed form expression in Fourier space [5]:

$$ \Gamma_0(\mathbf{P})(\mathbf{x}) = \sum_{\mathbf{n} \in \mathbb{Z}^d} \hat{\Gamma}_0(\mathbf{k}_\mathbf{n}) \tilde{\mathbf{P}}_\mathbf{n} \exp(i \mathbf{k}_\mathbf{n} \cdot \mathbf{x}), \qquad (10) $$

where $\mathbf{k}_\mathbf{n} = 2\pi \mathbf{n}/L$, $\tilde{\mathbf{P}}_\mathbf{n}$ is the n-th Fourier coefficient of $\mathbf{P}$ (our convention on Fourier transforms are recalled in Appendix B), and where the second-order tensor $\hat{\Gamma}_0(\mathbf{k})$ is given, for any $\mathbf{k} \in \mathbb{R}^d$, by

$$ \hat{\Gamma}_0(\mathbf{k}) = \frac{\mathbf{k} \otimes \mathbf{k}}{\sigma_0 |\mathbf{k}|^2} \quad \text{for } \mathbf{k} \neq \mathbf{0} \quad \text{and} \quad \hat{\Gamma}_0(\mathbf{0}) = \mathbf{0}. \qquad (11) $$

The Green operator $\Gamma_0$ being thus defined by (9), it can be shown that Problem 1 is equivalent to the Lippmann–Schwinger equation

$$ \mathbf{E} + \Gamma_0[(\sigma - \sigma_0)\mathbf{E}] = \bar{\mathbf{E}}, \qquad (12) $$

where the unknown is the electric field $\mathbf{E} \in (L^2(\mathcal{Q}))^d$. In turn, introducing the polarization $\mathbf{P} = (\sigma - \sigma_0)\mathbf{E}$, the above equation can be transformed as

$$ (\sigma - \sigma_0)^{-1} \mathbf{P} + \Gamma_0(\mathbf{P}) = \bar{\mathbf{E}}, \qquad (13) $$

and the apparent conductivity defined by (8) is obtained from the average polarization by

$$ \sigma^{\text{app}} \bar{\mathbf{E}} = \langle \mathbf{J} \rangle = \langle \sigma \mathbf{E} \rangle = \langle \sigma_0 \mathbf{E} + \mathbf{P} \rangle = \sigma_0 \langle \mathbf{E} \rangle + \langle \mathbf{P} \rangle = \sigma_0 \bar{\mathbf{E}} + \langle \mathbf{P} \rangle. \qquad (14) $$

Solving (13) is equivalent to solving Problem 1. In the following, we turn to discretization approaches for (13).
### 3.2.2. The equivalent inclusion method

The equivalent inclusion method of Moschovidis and Mura [22] delivers an approximate solution to the Lippmann–Schwinger equation (13) for materials made of well-defined individual phases.

We consider that the heterogeneous material under investigation is made of $N + 1$ phases $\Omega_1, \dots, \Omega_{N+1} \subset \mathbb{R}^d$ that form a partition of $\mathbb{R}^d$. “Phase” should be understood here in a loose sense: each $\Omega_\alpha$ has uniform conductivity $\sigma_\alpha$. As an example, we refer to the microstructure shown in Fig. 3 below, which contains $N$ inhomogeneities in $Q$. Each inhomogeneity (along with its periodic images in $\mathbb{R}^d$) is considered to be a phase. Phase $N + 1$ will be particularized, and referred to as the “matrix”.

We introduce the indicator function $\chi_\alpha: \mathbb{R}^d \to \{0, 1\}$ of $\Omega_\alpha$ and the volume fraction of phase $\alpha$ in the unit-cell $Q$, $f_\alpha = \langle \chi_\alpha \rangle = |\Omega_\alpha \cap Q|/|Q|$.

In the most simple form of the EIM, a piecewise constant approximate solution of (13) is constructed and the reference material coincides with the matrix: $\sigma_0 = \sigma_{N+1}$. The polarization $\mathbf{P}$ therefore vanishes over $\Omega_{N+1}$ (since $\mathbf{P}_{N+1} = (\sigma_{N+1} - \sigma_0)\mathbf{E}$). Based on this approximation space, discretization of the Lippmann–Schwinger equation initially relied on collocation [22]. More recently, a Galerkin discretization of the variational form of this equation was proposed [26]. As briefly outlined in Appendix C, Equation (13) reduces after discretization to the following linear system with unknowns $\mathbf{P}_1, \dots, \mathbf{P}_N \in \mathbb{R}^d$ (the constant polarizations in each phase; recall that $\mathbf{P}_{N+1} = \mathbf{0}$ in the matrix phase): for any $1 \le \alpha \le N$,

$$ (\sigma_{\alpha} - \sigma_{0})^{-1} \mathbf{P}_{\alpha} + \sum_{\beta=1}^{N} \Gamma_{\alpha\beta} \mathbf{P}_{\beta} = \bar{\mathbf{E}}. \qquad (15) $$

We are thus left with solving a linear system with $Nd$ scalar unknowns. In (15), the so-called second-order influence tensors $\Gamma_{\alpha\beta}$ are defined by how they operate on (constant) vectors: for any $\mathbf{P} \in \mathbb{R}^d$,

$$ \Gamma_{\alpha\beta} \mathbf{P} = |\Omega_{\alpha} \cap Q|^{-1} \int_{\Omega_{\alpha} \cap Q} \Gamma_0(\chi_{\beta} \mathbf{P})(\mathbf{x}) \, dx. \qquad (16) $$

The quantity $-\Gamma_{\alpha\beta}\mathbf{P}$ is the average electric field (in a periodic setting) induced in phase $\alpha$ by a polarization which is constant and equal to $\mathbf{P}$ in the phase $\beta$ and vanishes everywhere else. Note that, from the point of view of (16), $\Omega_\alpha$ and $\Omega_\beta$ should be seen as inclusions embedded in a homogeneous matrix with conductivity $\sigma_0$ (as discussed at the beginning of Sec. 3.2).

The EIM approximation of the apparent conductivity is deduced from the solution to (15) by the following identity, which directly follows from (14):

$$ \sigma^{\text{EIM}} \bar{\mathbf{E}} = \sigma_0 \bar{\mathbf{E}} + \sum_{\alpha=1}^{N} f_{\alpha} \mathbf{P}_{\alpha}. \qquad (17) $$

As usual in computational homogenization, the full determination of the conductivity tensor $\sigma^{\text{EIM}}$ requires (15) to be solved and (17) to be evaluated for $d$ linearly independent values of the constant macroscopic electric field $\bar{\mathbf{E}} \in \mathbb{R}^d$. Since we consider a random setting and a bounded SVE, the EIM approximation $\sigma^{\text{EIM}}$ of course depends on $\omega$.

**Remark 3** (Assembly of the linear system). The linear system $Ax = b$ corresponding to (15) is assembled by blocks:

$$ A_{(\alpha-1)d+i, (\beta-1)d+j} = \mathbf{e}_i \cdot [\delta_{\alpha\beta} (\sigma_{\alpha} - \sigma_0)^{-1} + \Gamma_{\alpha\beta}] \mathbf{e}_j \quad \text{and} \quad b_{(\alpha-1)d+i} = \bar{\mathbf{E}} \cdot \mathbf{e}_i, $$

where $1 \le \alpha, \beta \le N$ and $1 \le i, j \le d$. The solution $x$ to the linear system $Ax = b$ delivers the polarizations in each inclusion: $\mathbf{P}_\alpha \cdot \mathbf{e}_i = x_{(\alpha-1)d+i}$.

As pointed out in the introduction, there is a strong link between the EIM method and the Hashin-Shtrikman variational principle [40] and the Hashin-Shtrikman bounds [41]. The reformulation of the corrector problem 1 as the Lippmann–Schwinger equation (12), which is next recast as the equation (13) on the polarization $\mathbf{P}$, is at the basis of the Hashin-Shtrikman variational principle. The Hashin-Shtrikman bounds are obtained by taking $\sigma_0 = \sigma_{N+1}$ (the reference material coincides with the matrix) and searching an approximate solution to (13) in the form of a
polarization **P** which is constant in each inclusion, the constant being *the same* for all the inclusions and all the realizations. Here, we also take $\sigma_0 = \sigma_{N+1}$ but we search for an approximate solution to (13) in the form of a polarization **P** which is constant in each inclusion, the constant being a priori *different* for each inclusion and for each realization.

### 3.2.3. Fourier expansion of the influence tensors

To compute the above approximation $\sigma^{\text{EIM}}$, we need to compute beforehand the influence tensors $\Gamma_{\alpha\beta}$. We under-line in this section that this computation is actually challenging. However, in the present context where we use the EIM method only as a surrogate model (and not as an approximation of the apparent conductivity $\sigma^{\text{app}}$), it turns out, as explained in Sec. 4.2, that we can circumvent this difficult and adopt an alternative strategy.

The standard approach to compute $\Gamma_{\alpha\beta}$ is as follows. The following expression of the influence tensors is derived in Appendix D:

$$f_{\alpha} \Gamma_{\alpha\beta} = \sum_{\mathbf{n} \in \mathbb{Z}^d} \tilde{\chi}_{\alpha,n}^* \tilde{\chi}_{\beta,n} \widehat{\Gamma}_0(\mathbf{k}_n), \qquad (18)$$

where $\tilde{\chi}_{\alpha,n}^*$ is the complex conjugate of the Fourier coefficient $\tilde{\chi}_{\alpha,n}$ of the indicator function $\chi_{\alpha}$ of the phase $\alpha$,

$$\tilde{\chi}_{\alpha,n} = \frac{1}{|Q|} \int_Q \chi_{\alpha}(\mathbf{x}) \exp(-i \mathbf{k}_n \cdot \mathbf{x}) d\mathbf{x}. \qquad (19)$$

In the particular case when $d=2$ and when the phase $\alpha$ occupies a disk $\Omega_{\alpha}$ of center $\mathbf{x}_{\alpha}$ and of radius $a_{\alpha}$, we have, as shown in Appendix D, that

$$\tilde{\chi}_{\alpha,n} = f_{\alpha} F(k_{n} a_{\alpha}) e^{-i k_{n} \cdot x_{\alpha}} \quad \text{with} \quad F(\xi) = 2 J_{1}(\xi) / \xi, \qquad (20)$$

where $J_1$ denotes the Bessel function of the first kind and $k_n = |k_n|$. Inserting (20) into (18), we find, with $\mathbf{r}_{\alpha\beta} = \mathbf{x}_{\beta} - \mathbf{x}_{\alpha}$, that

$$\Gamma_{\alpha\beta} = f_{\beta} \sum_{\mathbf{n} \in \mathbb{Z}^2} F(k_{\mathbf{n}} a_{\alpha}) F(k_{\mathbf{n}} a_{\beta}) e^{-i \mathbf{k}_{\mathbf{n}} \cdot \mathbf{r}_{\alpha\beta}} \widehat{\Gamma}_0(\mathbf{k}_{\mathbf{n}}).$$

In dimension $d=2$, and when the phases occupy a regular shape, the general term of the series (18) behaves as $n^{-2}$ [42]. In other words, convergence is very slow, which makes the numerical evaluation of the influence tensors $\Gamma_{\alpha\beta}$ difficult. Several acceleration techniques have been proposed for an accurate evaluation of these coefficients [43]. We do not detail them here since, as pointed out above, in our particular context, we can adopt an alternative strategy, as described in Sec.4.2.

## 4. EIM-based variance reduction

We have shown in Sec. 3 how the EIM can deliver an approximation of the apparent conductivity. This estimate, given by (17), is generally much cheaper to compute than a full-field estimate, since it involves $Nd$ scalar unkowns only (where $N$ is the number of phases when ignoring the matrix phase). Hoping for a good correlation between the true apparent conductivity $\sigma^{\text{app}}(\omega)$ and its EIM estimate $\sigma^{\text{EIM}}(\omega)$, we now discuss how to use the EIM estimate as a surrogate model for variance reduction.

### 4.1. *Outline of the strategy*

Given a multi-phase material, our goal is to estimate the ensemble average $\mathbb{E}[\sigma^{\text{app}}(\omega)]$ (which is an average over all realizations of the microstructure). We recall that we are not concerned with finite-size effects in this article (the size $L$ of the SVE is therefore fixed).

As discussed in Sec. 3.1, the ensemble average $\mathbb{E}[\sigma^{\text{app}}(\omega)]$ is in practice approximated by an empirical average over $M$ realizations of the apparent conductivities $\sigma_m^{\text{app}}(\omega)$, $1 \le m \le M$, themselves evaluated by means of a full-field simulation (typically, using the finite element method – FEM):

$$\mathbb{E}(\sigma^{\text{app}}) \approx \frac{1}{M} \sum_{m=1}^{M} \sigma_{m}^{\text{FEM}}(\omega),$$
where $\sigma_m^{\text{FEM}}(\omega)$ denotes the finite element numerical approximation of the apparent conductivity of the m-th realization.

In order to increase the reliability of this estimate, we apply the variance reduction strategy introduced in Sec. 2 with the EIM estimate of the apparent conductivity as surrogate model. More precisely, in order to reduce the variance on the full-field estimate of one component of the apparent conductivity, we use the EIM estimate of the same component as a surrogate model.

**Remark 4** (Ensemble average of the surrogate model). Note that, for any component $1 \le i, j \le d$ of the conductivity tensor, the ensemble average $\mathbb{E}(\sigma_{ij}^{\text{EIM}})$ is not known. In the vein of Remark 1, a Monte-Carlo estimate is therefore used here, with a number $M$ of realizations that is significantly larger than the number $M$ of full-field simulations. Owing to the comparatively low cost of the equivalent inclusion method, these extra numerical computations do not jeopardize the overall efficiency of the proposed strategy.

**Remark 5.** It is in theory possible to use a tensor version of the variance reduction method described in Sec. 2, where any component of $\sigma^{\text{FEM}}$ is cross-correlated with all components of $\sigma^{\text{EIM}}$. We have implemented this approach and observed that it led to negligible improvements. The present article is therefore devoted to the case where only the component $\sigma_{ij}^{\text{EIM}}$ is used to improve the computation on $\mathbb{E}(\sigma_{ij}^{\text{FEM}})$.

To summarize the above discussion, the proposed workflow is outlined as follows (as recalled above, the size $L$ of the SVE is fixed):

1. Generate $M \gg M$ realizations of the microstructure.

2. Perform FEM simulations on the first $M$ microstructures. This step delivers the apparent conductivities $\sigma_1^{\text{FEM}}(\omega)$, ..., $\sigma_M^{\text{FEM}}(\omega)$.

3. Perform EIM simulations on all $M$ microstructures. For each realization, this step requires to assemble and solve the linear system (15); the EIM approximations $\sigma_1^{\text{EIM}}(\omega)$, ..., $\sigma_M^{\text{EIM}}(\omega)$ are obtained by (17).

4. Use all the results of Step 3 to compute an accurate estimate of $\mathbb{E}(\sigma^{\text{EIM}})$, in the vein of Remarks 1 and 4.

5. For any $1 \le i, j \le d$, use the results from Steps 2, 3 and 4 to apply the variance reduction technique described in Sec. 2 to $X = \sigma_{ij}^{\text{FEM}}$ and $Y = \sigma_{ij}^{\text{EIM}}$.

Of course, for this strategy to be viable, it is essential that Step 3 be performed in negligible time. The linear system to be solved is small. However, its assembly is complex, because of the slow convergence of the Fourier series (18) that defines the influence tensors. This problem is addressed in the next section.

## 4.2. Efficient evaluation of the influence tensors

The crucial observation is that we do not need an exact evaluation of the influence tensors $\Gamma_{\alpha\beta}$. Indeed, as a surrogate model, $\sigma^{\text{EIM}}$ is not required to be an accurate estimate of the true conductivity $\sigma^{\text{app}}$. The only requirement is that both random variables $\sigma^{\text{FEM}}$ and $\sigma^{\text{EIM}}$ be correlated. In other words, alterations of the EIM linear system are possible, as long as such alterations do not pollute the correlations between FEM and EIM simulations.

A natural idea is to replace the infinite series (18) with the truncated sum

$$f_{\alpha} \Gamma_{\alpha \beta} \approx \sum_{-P \le n_1, \dots, n_d \le P} \tilde{\chi}_{\alpha, n}^{*} \tilde{\chi}_{\beta, n} \widehat{\Gamma}_0(\mathbf{k}_n),$$

where $P$ is sufficiently small to keep the cost of this $d$-dimensional sum to a minimum (typically, $P \le 3$). Our numerical experiments (not reproduced here) show that such an approach does not lead to satisfactory variance reduction, which probably means that the neglected high-frequency terms have an important correlation content.

A different approach was suggested in [27]. Within the framework of “FFT-based homogenization methods” [44, 45], Zecevic and Lebensohn considered series similar to (18), $\chi_\alpha$ and $\chi_\beta$ being the indicator functions of two “voxels”. Using the Poisson formula, they first transformed the series in Fourier space into a series in the real-space, that they then truncate. We propose the same transformation based on the Poisson formula. However, unlike in [27], we will keep a minimum number of terms in the truncated sum. Our goal in this section is to provide an insight into the physical meaning of the transformation introduced in [27] and also propose a correction to this formula.
<figure><img src="image_4.png" /><figcaption>Figure 2: Unit-cell $\mathcal{Q}$ (solid blue line) and reference domain $\Omega_{\alpha}^0$ (black dashed line). The domain $\Omega_{\alpha}$ occupied by the phase $\alpha$ is the union of the principal image (centered at $\mathbf{x}_{\alpha}$, dark grey; by definition, $\mathbf{x}_{\alpha} \in \mathcal{Q}$; note that, in the case represented on the figure, the principal image $\mathbf{x}_{\alpha} + \Omega_{\alpha}^0$ is not a subset of $\mathcal{Q}$) and of its periodic images (light gray).</figcaption></figure>

The microstructure being periodic, it is interesting to define, for any phase $1 \le \alpha \le N$, the reference domain $\Omega_{\alpha}^0$, such that the whole phase $\Omega_{\alpha}$ is covered by periodic copies of $\mathbf{x}_{\alpha} + \Omega_{\alpha}^0$, $\mathbf{x}_{\alpha}$ being a fixed position in the unit-cell $\mathcal{Q}$. For example, in the case of circular inclusions, $\Omega_{\alpha}^0$ is the disk centered at the origin (with radius $a_{\alpha}$) and $\mathbf{x}_{\alpha}$ is the center of the sole periodic image the center of which lies in the unit-cell $\mathcal{Q}$ (see Fig. 2). From the above definitions, we have ('+' referring to translation)

$$ \Omega_{\alpha} = \bigcup_{\mathbf{n} \in \mathbb{Z}^d} [(\mathbf{x}_{\alpha} + \mathbf{n} \mathcal{L}) + \Omega_{\alpha}^0]. $$

It is recalled (see (16)) that, up to a sign change, $\Gamma_{\alpha\beta}$ measures the average electric field induced in phase $\Omega_{\alpha}$ by a constant polarization of phase $\beta$. Owing to the periodicity of the microstructure, averaging over phase $\Omega_{\alpha}$ is equivalent to averaging over the principal image $\mathbf{x}_{\alpha} + \Omega_{\alpha}^0$. Besides, $\Omega_{\beta}$ is made of the *principal image* $\mathbf{x}_{\beta} + \Omega_{\beta}^0$ and all its periodic images. The tensor $\Gamma_{\alpha\beta}$ therefore measures the average electric field induced in the fundamental image $\mathbf{x}_{\alpha} + \Omega_{\alpha}^0$ by a constant polarization of all periodic images of $\mathbf{x}_{\beta} + \Omega_{\beta}^0$. In short, $\Gamma_{\alpha\beta}$ measures the interaction between $\mathbf{x}_{\alpha} + \Omega_{\alpha}^0$ and all periodic images of $\mathbf{x}_{\beta} + \Omega_{\beta}^0$.

Since the conductivity problem under consideration is linear, these interactions are *additive*, and we can focus on the interaction between $\mathbf{x}_{\alpha} + \Omega_{\alpha}^0$ and one periodic image $\Omega_{\beta}^1 = \mathbf{x}_{\beta} + \mathbf{n}_1 \mathcal{L} + \Omega_{\beta}^0$ (for a given $\mathbf{n}_1 \in \mathbb{Z}^d$) of $\Omega_{\beta}^0$. By invariance upon translation, we focus on the interaction between $\Omega_{\alpha}^0$ and $\Omega_{\beta}^1 = \mathbf{r} + \Omega_{\beta}^0$, where $\mathbf{r} = \mathbf{x}_{\beta} - \mathbf{x}_{\alpha} + \mathbf{n}_1 \mathcal{L}$.

Intuition suggests that we are no longer working with periodic boundary conditions. Rather, $\Omega_{\alpha}^0$ and $\Omega_{\beta}^1$ are embedded in the infinite space $\mathbb{R}^d$, and boundary conditions are applied at infinity (to ensure square integrability of the electric field). We then consider the following problem, for a constant vector $\mathbf{P}_{\beta} \in \mathbb{R}^d$:

$$ \mathrm{div} \, \mathbf{J} = 0 \quad \text{in} \, \mathbb{R}^d, \qquad (21) $$

$$ \mathbf{J} = \sigma_0 \mathbf{E} + \mathbf{P}_{\beta} \quad \text{in} \quad \Omega_{\beta}^{1}, \qquad (22) $$

$$ \mathbf{J} = \sigma_0 \mathbf{E} \quad \text{in} \, \mathbb{R}^d \setminus \Omega_{\beta}^{1}, \qquad (23) $$

$$ \mathbf{E} = - \mathrm{grad} \, \Phi \quad \text{in} \, \mathbb{R}^d, \qquad (24) $$

where $\Phi \in L_{\mathrm{loc}}^2(\mathbb{R}^d)$ and $\mathbf{E} \in (L^2(\mathbb{R}^d))^d$, and we are interested in the average electric field over $\Omega_{\alpha}^0$. In other words, in an infinite medium with conductivity $\sigma_0$, we apply the polarization $\mathbf{P}_{\beta}$ to the domain $\Omega_{\beta}^1 = \mathbf{r} + \Omega_{\beta}^0$, and measure the
average induced electric field over $\Omega_{\alpha}^0$. By definition, the influence tensor $\Gamma_{\alpha\beta}^{\infty}(\mathbf{r})$ maps the polarization $\mathbf{P}_{\beta}$ onto the opposite of the average electric field:

$$ \Gamma_{\alpha\beta}^{\infty}(\mathbf{r}) \mathbf{P}_{\beta} = -|\Omega_{\alpha}^{0}|^{-1} \int_{\Omega_{\alpha}^{0}} \mathbf{E}, \qquad (25) $$

where $\mathbf{E}$ is the solution to (21)–(24). The influence tensor $\Gamma_{\alpha\beta}^{\infty}(\mathbf{r})$ of course depends on the geometry of the reference domains $\Omega_{\alpha}^0$ and $\Omega_{\beta}^0$, and on the vector $\mathbf{r} = \mathbf{x}_{\beta} - \mathbf{x}_{\alpha} + \mathbf{n}_1 L$.

As an example, in Appendix E, we consider the two-dimensional case, for isotropic homogeneous reference materials ($\sigma_0 = \sigma_0 \mathbf{I}$) and when the phases are made of circular inclusions. We denote $a_{\beta}$ the radius of the inclusions of phase $\beta$. In this case, we show (see (E.8)) that

$$ \Gamma_{\alpha\beta}^{\infty}(\mathbf{r}) = \frac{a_{\beta}^2}{2\sigma_0 |\mathbf{r}|^2} \left( \mathbf{I} - 2 \frac{\mathbf{r} \otimes \mathbf{r}}{|\mathbf{r}|^2} \right), \qquad (26) $$

where $\mathbf{r}$ is then the vector going from the center of $\Omega_{\alpha}^0$ to the center of $\Omega_{\beta}^1$.

Note that the only difference between $\Gamma_{\alpha\beta}$ and $\Gamma_{\alpha\beta}^{\infty}$ lies in the boundary conditions (Problem 2 is complemented with periodic boundary conditions, while Problem (21)–(24) is complemented with the “boundary condition” $\mathbf{E} \in (L^2(\mathbb{R}^d))^d$).

The above discussion leads to propose the following formula:

$$ \Gamma_{\alpha\beta} \stackrel{?}{=} \sum_{\mathbf{n} \in \mathbb{Z}^d} \Gamma_{\alpha\beta}^{\infty}(\mathbf{x}_{\beta} - \mathbf{x}_{\alpha} - \mathbf{n} L), \qquad (27) $$

which expresses the fact that $\Gamma_{\alpha\beta}$ measures the (additive) interactions between $\mathbf{x}_{\alpha} + \Omega_{\alpha}^0$ and all periodic images of $\mathbf{x}_{\beta} + \Omega_{\beta}^0$ (as measured by $\Gamma_{\alpha\beta}^{\infty}$).

This is, in essence, the formula initially proposed in [27] (based on the Poisson summation formula). Note that a similar formula was also proposed by Molinari and El Mouden [46]. However, Zecevic and Lebensohn [27] observed that this sum is not absolutely convergent. It is therefore at best conditionally convergent and the summation scheme must be defined precisely. In [27], the following scheme was adopted:

$$ \Gamma_{\alpha\beta} \stackrel{?}{=} \lim_{P \to +\infty} \sum_{-P \le n_1, \dots, n_d \le P} \Gamma_{\alpha\beta}^{\infty}(\mathbf{x}_{\beta} - \mathbf{x}_{\alpha} - \mathbf{n} L). \qquad (28) $$

While appealing because of its reasonable physical meaning, the above formula cannot be correct. Consider indeed the case of isotropic reference materials, in dimension $d=2$ and for circular inclusions. It is shown in Appendix D (see (D.1)) that, for any $\alpha \neq \beta$,

$$ \mathrm{tr} \Gamma_{\alpha\beta} = -\sigma_0^{-1} f_{\beta}, $$

while we deduce from (26) that, for any $\alpha \neq \beta$,

$$ \mathrm{tr} \Gamma_{\alpha\beta}^{\infty} = 0. $$

Equation (28) is inconsistent with the two above relations. Instead, we (heuristically) postulate the following correction (for the two-dimensional case of disks and isotropic reference materials):

$$ \Gamma_{\alpha\beta} = \lim_{P \to +\infty} \sum_{-P \le n_1, \dots, n_d \le P} \Gamma_{\alpha\beta}^{\infty}(\mathbf{x}_{\beta} - \mathbf{x}_{\alpha} - \mathbf{n} L) - \frac{f_{\beta}}{d \sigma_0} \mathbf{I}, \qquad (29) $$

where $\mathbf{I}$ is the identity tensor, which ensures that the traces of both sides coincide by sharing equally between the diagonal coefficients of $\Gamma_{\alpha\beta}$ the error on the trace. We have not been able to give a mathematical proof of (29). However, a numerical validation in dimension $d=2$ is proposed in Appendix F. Besides, it is shown below that this correction leads to the FEM and EIM results being satisfactorily correlated.

In practice, for the applications presented below, we use a truncated form of (29):

$$ \Gamma_{\alpha\beta} \approx \sum_{-P \le n_1, \dots, n_d \le P} \Gamma_{\alpha\beta}^{\infty}(\mathbf{x}_{\beta} - \mathbf{x}_{\alpha} - \mathbf{n} L) - \frac{f_{\beta}}{2\sigma_0} \mathbf{I}, $$
with a low value of P (typically P = 1, 2, 3). Unless otherwise stated, we set P = 2 in Sec. 5. Note also that we have considered the approximation

$$ \Gamma_{\alpha\beta} \approx \sum_{-P \le n_1, \dots, n_d \le P} \Gamma_{\alpha\beta}^{\infty}(\mathbf{x}_{\beta} - \mathbf{x}_{\alpha} - \mathbf{n} \cdot \mathbf{L}), $$

which leads to poor results (not reproduced here) in terms of variance reduction efficiency.

**Remark 6.** The additional term introduced in (29) is certainly **not universal**. Indeed, this term being isotropic, it is questionable that it should apply to situations where inclusions are not circular and/or the unit-cell is not a square. We leave the study of such general settings to future works.

**Remark 7.** Zecevic and Lebensohn [27] showed that (27) delivers good results for the numerical applications they considered. At this stage, the reason for this is unclear to us.

The workflow outlined in Sec. 4 above is applied in the next section to random assemblies of circular inclusions.

## 5. Numerical examples

Note: the following numerical examples can be reproduced from the dataset at the Recherche Data Gouv open data repository [28].

In this section, we consider two applications of the EIM-based variance reduction technique to 2D random assemblies of N = 32 circular inhomogeneities. The size L of the square simulation box and the radius a of the disks are adjusted so as to ensure a prescribed value of the fraction f of inhomogeneities.

A standard Monte-Carlo algorithm as described by Allen and Tildesley [47, Chap. 4] is implemented and used to generate all microstructures below. Starting from a centered square orthorombic lattice<fnref n="1" />, the simulation is run for 100 000 cycles (during one cycle, a trial move is proposed for each particle). The random amplitude of the trial moves is adjusted every 50 cycles in order to ensure that the acceptance ratio is about 0.3: if the acceptance ratio is too large (resp. too small), the amplitude of the trial moves is scaled by 1.05 (resp. 0.95). The microstructures thus generated are subsequently meshed (mesh element size: h ≈ a/8), as shown in Fig. 3.

The (isotropic) conductivity of the matrix is σ₀. The conductivity of the inhomogeneities is χσ₀, where χ denotes the contrast. We use the theoretical developments of Sec. 3 and 4, where each circular inhomogeneity in Q (along with its periodic images in ℝ<sup>d</sup>) represents a phase. The number of realizations considered for the FEM simulations is M = 100. The number of realizations considered for the EIM simulations is M = 10 000. On the laptop used to perform the computations (equipped with an Intel(R) Core(TM) i9-9880H at 2.30 GHz), solving one full-field problem takes on average 1.5 second, while solving one EIM model is 1700 times less expensive and takes on average 0.88 millisecond.

**Remark 8** (On the periodicity of the microstructures). The microstructures considered until now are **not necessarily periodic**. Rather, the SVE is in fact the **restriction** to the finite-size unit-cell Q of an infinite realization of the random heterogeneous material under consideration. As a consequence, the volume fraction within the unit-cell Q is in general a random variable that may fluctuate around its expectation.

Previous studies [18, 17] have shown that these fluctuations are indeed the leading sources of statistical errors and can be used as control variates (or stratification criteria) to efficiently reduce the variance. Our goal in the present article is to investigate **further** sources of statistical errors, related to the spatial distribution of the phases (rather than their mere volume fractions). In order to focus on these other sources, we fix the volume fractions in the microstructures that we generate, thus killing the fluctuations of the apparent conductivity that are induced by the fluctuations of volume fractions. This is most easily achieved by **enforcing Q-periodicity of the microstructure**, which we will systematically do in the remainder of this article. In other words, all microstructures considered in the sequel are geometrically periodic.

**Remark 9.** We again note that we are not concerned with size-effects in this work. However, it is generally accepted that the assemblies considered here (N = 32 inhomogeneities) are large enough to deliver an estimate of the effective properties with limited bias. The dispersion remains large enough for variance reduction techniques to be worthwhile [48, 10].

<footnote n="1">¹https://en.wikipedia.org/wiki/Orthorhombic_crystal_system#2D_Bravais_lattices, last retrieved 2022-07-27.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 3: A realization of the microstructure considered in Sec. 5.</figcaption></figure>

### 5.1. Effect of elastic contrast

In this section, we select a high value for the fraction of inclusions, $f = 50\%$ (this corresponds to $L/(2a) \approx 7.08$) and investigate how the value of the contrast $\chi$ affects the efficiency of the variance reduction. To do so, we run simulations for $\chi \in \{10^{-6}, 10^{-5}, \dots, 10^{-1}, 10, \dots, 10^5, 10^6\}$.

Figure 4 displays, for the $M = 100$ full-field computations, the EIM estimate of $\sigma_{11}$ vs. the FEM estimate of the same quantity. More precisely, $M = 100$ realisations of the relative error (with respect to $\mathbb{E}(\sigma_{11}^{\text{FEM}})$)

$$ \varepsilon^*(\omega) = \frac{\sigma_{11}^*(\omega) - \mathbb{E}(\sigma_{11}^{\text{FEM}})}{\mathbb{E}(\sigma_{11}^{\text{FEM}})}, \qquad (30) $$

where $\star$ stands for “FEM” or “EIM”, are displayed in Fig. 4. We first observe that the point clouds on that figure are rather far from the first bisector. In other words, the EIM provides a *biased* estimate of the conductivity. The rather large bias can be attributed to the poor discretization of the true polarization over each inhomogeneity (here, the polarization is locally simply approximated by a constant; polynomial approximations do improve the quality of the estimate, at the expense of a higher computational cost [26]). We also observe that the two point clouds roughly form two lines: FEM and EIM values are *strongly correlated* and the proposed variance reduction method is therefore expected to perform efficiently.

Figure 5 shows the FEM and EIM estimates of the reduced conductivity $\sigma_{11}/\sigma_0$ vs. the contrast $\chi$ (more precisely, we plot the empirical mean over $M$ – resp. $M$ – realisations as well as the 99% confidence interval for the FEM estimate – resp. the EIM estimate). As expected, $\sigma_{11}$ increases when $\chi$ increases. It is also observed that the error bars are much smaller for the EIM estimates than for the FEM estimates: this is simply due to the fact that the former are computed on the basis of a much larger number of realisations ($M = 10000 \gg M = 100$) than the latter. As already observed, the EIM does not deliver a good estimate of the conductivity (the green and purple curves of Fig. 5 are significantly different). Figures 5 and 6 also display the numerical estimates of the reduced conductivity after application of the proposed variance reduction strategy. Clearly, the resulting error bars (on the blue curve) are reduced significantly in comparison to the original error bars (on the purple curve), leading to a much more accurate estimation of $\mathbb{E}(\sigma_{11}^{\text{FEM}})$ at constant computational cost.

More quantitatively, Figure 7 shows the efficiency index (defined at the very end of Sec. 2) of the variance reduction technique vs. the contrast $\chi$, for $P = 1, 2$ and $3$. Even for highly contrasted microstructures, the standard deviation is
<figure><img src="image_6.png" /><figcaption>Figure 4: Plot of ε<sup>EIM</sup>(ω) as a function of ε<sup>FEM</sup>(ω) (these two quantities are defined by (30)) for M = 100 realisations (we consider two values of the contrast: χ = 10<sup>-6</sup> and χ = 10<sup>6</sup>).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5: FEM and EIM estimates of the reduced conductivity σ<sub>11</sub>/σ<sub>0</sub>, for various values of the contrast χ (note that the vertical axis is not the same on the left plot and on the right plot). For the FEM curve (resp. EIM curve), we show the empirical mean over M (resp. M) realisations as well as the 99 % confidence interval (displayed as error bars). The values obtained after variance reduction (following the workflow detailed in Section 4.1) are also displayed on the plot. A close-up of these curves is shown in Fig. 6.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: Same as in Fig. 5, with dilated vertical scales (the EIM results do not appear on the scale of the plot).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Efficiency index of the proposed variance reduction technique, as a function of contrast $\chi$ and truncation parameter $P$.</figcaption></figure>

significantly reduced by a factor of about 2.3 to 2.4. This demonstrates the value of the proposed technique. For the case considered here, the efficiency index seems to have reached a plateau for any value of $\chi \le 10^{-4}$ (resp. any value of $\chi \ge 10^4$), which is an encouraging evidence that the efficiency may remain of the same order of magnitude for even smaller contrasts than $10^{-6}$ (resp. larger contrasts than $10^6$).

We observe that the efficiency index increases with increasing truncation parameter $P$. However, improvements are barely noticeable for $P = 3$, and our recommended value is $P = 2$. We also observe that, for contrasts smaller than 1, the efficiency increases when $\chi$ becomes closer to 1 (and similarly for $\chi > 1$). In the limit $\chi \to 1$, the variance (denoted $\sigma_Y^2$ and $\sigma_X^2$, following the notation of Sec. 2) of the EIM and the FEM conductivities converges to 0 (when $\chi = 1$, the materials is simply homogeneous), the variance $\sigma_Z^2$ of the reduced variable $Z$ converges to 0 as well, and numerical experiments show that the efficiency index $\sigma_X/\sigma_Z$ actually increases.

**Remark 10.** When $\chi = 10^6$ and $P = 2$, the statistical error on $E(\sigma_{11}^{\text{EIM}}/\sigma_0)$ and on $E(\sigma_{11}^{\text{FEM}}/\sigma_0)$ (corresponding to a 99% confidence interval) are about 0.0024 and 0.0514, respectively. The efficiency index is about 2.333, and the optimum value of the variable $\xi$ is $\xi^* = 1.8926$. We make here a more precise computation of the efficiency index, taking into account the uncertainty on the expectation of $\sigma_{11}^{\text{EIM}}$. After variance reduction, and since the random variable $\sigma_{11}^{\text{FEM}} - \xi^* \sigma_{11}^{\text{EIM}}$ is independent from the random variables used to estimate the expectation of $\sigma_{11}^{\text{EIM}}$, the total
<figure><img src="image_10.png" /><figcaption>Figure 8: Examples of microstructures considered in Sec. 5.2. Left: g/a = 0; Center: g/a = 0.121; Right: g/a = 0.236 (we recall that a is the radius of the inhomogeneities and that g is the minimum gap between inhomogeneities).</figcaption></figure>

<table><caption>Table 1: Parameters for the Monte-Carlo simulations of Sec. 5.2, with the following symbols: f' and f are the initial and final inhomogeneities fraction; a' and a are the initial and final radius; g is the minimum gap between disks.</caption><thead><tr><th>f'</th><th>f</th><th>a/a'</th><th>g/a</th></tr></thead><tbody><tr><td>0.40</td><td>0.40</td><td>1.000</td><td>0.000</td></tr><tr><td>0.45</td><td>0.40</td><td>0.943</td><td>0.121</td></tr><tr><td>0.50</td><td>0.40</td><td>0.894</td><td>0.236</td></tr></tbody></table>

error on $\sigma_{11}^{app}/\sigma_0$ is

$$ \sqrt{\left(\frac{0.0514}{2.333}\right)^2 + (1.8926 \times 0.0024)^2} \approx 0.0225, $$

leading to an **effective** efficiency index of $0.0514/0.0225 \approx 2.285$, just slightly smaller than 2.333. Should need be, increasing $M$ would improve this effective efficiency index.

## 5.2. Application to the discrimination of statistically different microstructures

In the present application, we consider three families of assemblies of disks. For all three families, the total number of inhomogeneities is $N = 32$. They cover $f = 40\%$ of the unit-cell $Q$. As shown in Tab. 1, the three microstructures differ by the minimum allowed gap between two inhomogeneities, ranging from 0 (the inhomogeneities may touch) to 23.6% of their radius (see Fig. 8).

The three families of inclusions are generated with the same Monte-Carlo code as previously, starting from an initial fraction $f' \ge f$ of inhomogeneities (we consider $N = 32$ inhomogeneities with a slightly larger radius $a'$). After equilibration, the initial radius $a'$ of the inclusions is scaled down to the final value $a$, so as to reach the desired volume fraction $f$. The minimum allowed gap $g$ is defined by the relation $2a' = 2a + g$. It corresponds to the gap between two inhomogeneities of radius $a$ that would have been tangent one with each other when their radius was equal to $a'$.

The contrast is set to $\chi = 10^6$, and the apparent conductivities of the microstructures are estimated with $M = 10$ finite element simulations ($M$ is kept very low on purpose in this application). The results are shown in Fig. 9, where it is observed that the error bars are too large to allow for the discrimination of the three types of microstructures. The EIM-based variance reduction strategy is next applied (the expectation of the surrogate model is again estimated using $M = 10000$ realisations). As can be seen in Fig. 9, discrimination is possible after variance reduction: the apparent conductivity does depend on the gap, and actually decreases when the gap increases.

**Remark 11.** In the present article, no attempt has been made to optimize the number $M$ of realizations used to evaluate the ensemble average of the surrogate model. The large value ($M = 10000$) selected here leads to a small statistical error on the evaluation of $E(\sigma^{EIM})$. We have checked that, although the error bars on the quantity of interest get larger, reducing $M$ down to $M = 100$ still allows to separate the three families of microstructures (see Table 2).
<figure><img src="image_11.png" /><figcaption>Figure 9: Estimates of the (1, 1) component of the apparent conductivity for the microstructures considered in Sec. 5.2. Without variance reduction (“w/o”), the error bars are so large that the three microstructures are statistically undistinguishable. With the variance reduction strategy (“w”), the error bars are sufficiently small so that it becomes apparent that the confidence intervals do not overlap.</figcaption></figure>

<table><caption>Table 2: Bounds of the 99% confidence interval on the (1, 1) component of the apparent conductivity for the three microstructures considered in Sec. 5.2, in the case when $M = 100$. The values shown in this table account for the statistical errors on the surrogate model.</caption><thead><tr><th>f'</th><th>Lower bound</th><th>Upper bound</th></tr></thead><tbody><tr><td>0.40</td><td>2.563</td><td>2.701</td></tr><tr><td>0.45</td><td>2.467</td><td>2.525</td></tr><tr><td>0.50</td><td>2.425</td><td>2.426</td></tr></tbody></table>

Note also that $M$ should be chosen sufficiently large so that the statistical error on the evaluation of $\mathbb{E}(\sigma^{\text{EIM}})$ does not become the dominant error.

# 6. Conclusion

Within the framework of linear conductivity, we have proposed a variance reduction strategy based on the use of a surrogate model, here the equivalent inclusion method, for the estimation of the macroscopic (apparent) properties. This strategy has been implemented in a two-dimensional setting, for assemblies of disks. Even for high values of the matrix-inhomogeneity conductivity contrast (namely $10^6$ and $10^{-6}$), the proposed strategy performs very satisfactorily: its efficiency index (which quantifies the accuracy gain at constant computational cost) remains above 2.2. This is quite remarkable, since the equivalent inclusion method is known to deliver poor estimates of the apparent properties (especially in the regime of high contrast and/or high volume fraction). Still, the EIM method seems to be able to capture some salient features of the microstructures that lead to a high correlation with the true apparent conductivity.

There are several perspectives to the present work. First, from a methodological point of view, the framework presented here readily extends to linear elasticity: we refer to [26] for a presentation of the equivalent inclusion method in that setting. In fact, the proposed technique can be applied to any physical setting (including coupled problems) for which a Green operator can be defined (see [49] for the case of piezoelectric composites). In the same spirit, our approach also carries over, from a methodological standpoint, to 3D computations.

Second, owing to the similarities between the equivalent inclusion method and the *self-consistent clustering* method [50], a similar variance reduction strategy could be proposed in situations where the microstructure is not an assembly of elementary geometrical objects or does not behave linearly. A full-field approach based on the fast Fourier transform [44, 45] must then be adopted. Of course, this increases the cost of the reduced model (one evaluation of the surrogate model being roughly equivalent to a few iterations of the full-field simulation), but the strategy might still be competitive. Likewise, there are deep connexions between the equivalent inclusion method and the
recent works by Bonnet, Huang, Monchiet, Nguyen and To [43, 51, 52], that should be reassessed from the point of view of variance reduction.

Finally, a heuristic correction has been proposed to a recent formula by Zecevic and Lebensohn [27]. This formula (actually its corrected version) is used here to efficiently compute the influence tensors when periodic boundary conditions are considered. We have numerically investigated the validity of this corrected formula in the case of circular inclusions embedded in a square unit-cell. Extending the correction to other settings and producing a mathematical proof of the corrected formula is another exciting perspective.

### A. $H(\text{div}, Q)$ space

We recall here some results about the $H(\text{div}, Q)$ space. For any bounded domain $Q$ of $\mathbb{R}^d$, we recall that

$$H(\text{div}, Q) := \{v \in (L^2(Q))^d, \text{div } v \in L^2(Q)\},$$

and that $H(\text{div}, Q)$ is a Hilbert space. The space $(C^\infty(\bar{Q}))^d$ is dense in $H(\text{div}, Q)$. Consider the normal trace application

$$\gamma: \begin{cases} (C^\infty(\bar{Q}))^d & \to C^0(\partial Q) \\ v & \mapsto (v \cdot n)|_{\partial Q} \end{cases},$$

where $n$ denotes the unit exterior normal vector to $\partial Q$. The application $\gamma$ can be uniquely extended (as a continuous application) on the whole space $H(\text{div}, Q)$, which in particular implies that the normal trace $v \cdot n$ of any field $v \in H(\text{div}, Q)$ is well-defined.

### B. Fourier series and Fourier transforms

We define here our conventions for the Fourier transform of a periodic function. We consider a $Q$-periodic tensor field $T$, with components that are square-integrable over $Q$ (we recall that $Q = (-L/2, L/2)^d$). For any $n \in \mathbb{Z}^d$, its Fourier coefficient $\tilde{T}_n$ is defined by

$$\tilde{T}_n = \frac{1}{|Q|} \int_Q T(x) e^{-i k_n \cdot x} dx, \quad \text{with} \quad k_n = \frac{2\pi n}{L}.$$

The synthesis formula reads $T(x) = \sum_{n \in \mathbb{Z}^d} \tilde{T}_n e^{i k_n \cdot x}$, where the right-hand side is a converging series in $(L^2(Q))^{d \times d}$.

### C. Variational formulation of the equivalent inclusion method

Multiplying (13) with an arbitrary test function $Q$ (a vector-valued field with square integrable components) and averaging over the unit-cell delivers the weak form of the Lippmann-Schwinger equation (13), that is written as Problem 3 below:

**Problem 3.** Find $P \in (L^2(Q))^d$ such that, for any $Q \in (L^2(Q))^d$,

$$\langle Q \cdot (\sigma - \sigma_0)^{-1} P \rangle + \langle Q \cdot \Gamma_0(P) \rangle = \bar{E} \cdot \langle Q \rangle. \quad (C.1)$$

Note that the bilinear form in the above variational equation is symmetric (the Green operator $\Gamma_0$ owes its symmetry to the symmetry of the conductivity tensor $\sigma_0$). The well-posedness of (C.1) can be studied directly (without using the equivalence of the Lippmann-Schwinger equation (13) with Problem 1). Such a study (which has been been performed in, e.g., [53, 54]) opens the way to establishing the well-posedness of Galerkin discretizations of (C.1).

The variational form of the EIM results from the Galerkin discretization of (C.1) with phase-wise constant trial and test fields:

$$\mathbf{P}(\mathbf{x}) = \sum_{\alpha=1}^{N} \chi_{\alpha}(\mathbf{x}) \mathbf{P}_{\alpha} \quad \text{and} \quad \mathbf{Q}(\mathbf{x}) = \sum_{\alpha=1}^{N} \chi_{\alpha}(\mathbf{x}) \mathbf{Q}_{\alpha}, \quad (C.2)$$
where $\mathbf{P}_\alpha, \mathbf{Q}_\alpha \in \mathbb{R}^d$ are constant vectors (we recall that $\chi_\alpha$ is the indicator function of the phase $\alpha$). Inserting (C.2) in (C.1) and using the definition (16) of the influence tensors $\Gamma_{\alpha\beta}$ yields the discrete variational Problem 4 below.

**Problem 4.** Find $\mathbf{P}_1, \dots, \mathbf{P}_N \in \mathbb{R}^d$ such that, for all $\mathbf{Q}_1, \dots, \mathbf{Q}_N \in \mathbb{R}^d$,

$$ \sum_{\alpha=1}^{N} f_{\alpha} \mathbf{Q}_{\alpha} \cdot (\sigma_{\alpha} - \sigma_{0})^{-1} \mathbf{P}_{\alpha} + \sum_{\alpha, \beta=1}^{N} f_{\alpha} \mathbf{Q}_{\alpha} \cdot \Gamma_{\alpha\beta} \mathbf{P}_{\beta} = \sum_{\alpha=1}^{N} f_{\alpha} \bar{\mathbf{E}} \cdot \mathbf{Q}_{\alpha}. $$

Since the vectors $\mathbf{Q}_1, \dots, \mathbf{Q}_N$ are arbitrary, Problem 4 results in the linear system (15).

### D. On the influence tensor (periodic BCs)

We derive here the expression (18) of the influence tensor of two inhomogeneities in a periodic setting. Recalling that the influence tensor $\Gamma_{\alpha\beta}$ is defined by (16), we first note, owing to $\Gamma_0$ being self-adjoint, that $f_\alpha \Gamma_{\alpha\beta} = f_\beta \Gamma_{\beta\alpha}$, where we recall that $f_\alpha$ denotes the volume fraction of inhomogeneity $\alpha$.

When the phase $\beta$ is subjected to the constant polarization $\mathbf{P}_\beta$, we have $\mathbf{P}(\mathbf{x}) = \chi_\beta(\mathbf{x}) \mathbf{P}_\beta$ at any point $\mathbf{x} \in \mathcal{Q}$. In view of the expression (10) of the Green operator $\Gamma_0$ in the Fourier space, we can express the induced electric field $\mathbf{E} = -\Gamma_0(\chi_\beta \mathbf{P}_\beta)$ at any point as

$$ \mathbf{E}(\mathbf{x}) = - \sum_{\mathbf{n} \in \mathbb{Z}^d} \tilde{\chi}_{\beta, \mathbf{n}} \exp(i \mathbf{k}_{\mathbf{n}} \cdot \mathbf{x}) \hat{\Gamma}_0(\mathbf{k}_{\mathbf{n}}) \mathbf{P}_{\beta}, $$

where the Fourier coefficients $\tilde{\chi}_{\beta, \mathbf{n}}$ of the indicator function $\chi_\beta$ are defined by (19).

The above expression is then averaged over the phase $\alpha$. In view of (16), we have that, for any $\mathbf{P}_\beta$,

$$ \Gamma_{\alpha\beta} \mathbf{P}_{\beta} = -|\Omega_{\alpha} \cap \mathcal{Q}|^{-1} \int_{\mathcal{Q}} \chi_{\alpha}(\mathbf{x}) \mathbf{E}(\mathbf{x}) \, d\mathbf{x}, $$

from which we deduce, observing that $|\Omega_\alpha \cap \mathcal{Q}| = f_\alpha |\mathcal{Q}|$, that

$$ f_{\alpha} \Gamma_{\alpha\beta} = |\mathcal{Q}|^{-1} \int_{\mathcal{Q}} \sum_{\mathbf{n} \in \mathbb{Z}^d} \chi_{\alpha}(\mathbf{x}) \tilde{\chi}_{\beta, \mathbf{n}} \hat{\Gamma}_0(\mathbf{k}_{\mathbf{n}}) \exp(i \mathbf{k}_{\mathbf{n}} \cdot \mathbf{x}) \, d\mathbf{x}. $$

Recognizing the expression (19) of $\tilde{\chi}_{\alpha, \mathbf{n}}$, Equation (18) is eventually retrieved.

We next derive a closed form expression of the trace of the influence tensor $\Gamma_{\alpha\beta}$, in the case of isotropic reference materials. From the general expression (18), we infer

$$ f_{\alpha} \operatorname{tr} \Gamma_{\alpha\beta} = \sum_{\mathbf{n} \in \mathbb{Z}^d} \tilde{\chi}_{\alpha, \mathbf{n}}^{*} \tilde{\chi}_{\beta, \mathbf{n}} \operatorname{tr} [\hat{\Gamma}_0(\mathbf{k}_{\mathbf{n}})]. $$

In view of (11), we have $\operatorname{tr}[\hat{\Gamma}_0(\mathbf{0})] = 0$ and $\operatorname{tr}[\hat{\Gamma}_0(\mathbf{k})] = \sigma_0^{-1}$ for any $\mathbf{k} \neq \mathbf{0}$. We therefore deduce that

$$ \sigma_0 f_{\alpha} \operatorname{tr} \Gamma_{\alpha\beta} = \sum_{\mathbf{n} \in \mathbb{Z}^d} \tilde{\chi}_{\alpha, \mathbf{n}}^{*} \tilde{\chi}_{\beta, \mathbf{n}} - \tilde{\chi}_{\alpha, \mathbf{0}}^{*} \tilde{\chi}_{\beta, \mathbf{0}}. $$

Parseval's theorem is applied to the first term, and we notice that $\tilde{\chi}_{\alpha, \mathbf{0}} = \langle \chi_\alpha \rangle = f_\alpha$ in the second term. We thus obtain

$$ \sigma_0 f_{\alpha} \operatorname{tr} \Gamma_{\alpha\beta} = \langle \chi_{\alpha} \chi_{\beta} \rangle - f_{\alpha} f_{\beta}. $$

If $\alpha \neq \beta$, the domains $\Omega_\alpha$ and $\Omega_\beta$ do not overlap, and $\langle \chi_\alpha \chi_\beta \rangle = 0$. If $\alpha = \beta$, $\chi_\alpha \chi_\beta = \chi_\alpha$ and $\langle \chi_\alpha \chi_\beta \rangle = f_\alpha$. We eventually obtain

$$ \operatorname{tr} \Gamma_{\alpha\beta} = \begin{cases} -f_{\beta}/\sigma_0 & \text{when } \alpha \neq \beta, \\ (1-f_{\alpha})/\sigma_0 & \text{when } \alpha = \beta. \end{cases} \qquad (D.1) $$
<figure><img src="image_12.png" /><figcaption>Figure E.10: Geometric setting for the derivation of $\Gamma_{12}^{\infty}$.</figcaption></figure>

To close this section, we consider the particular case when $d = 2$ and when the phase $\alpha$ occupies a disk $\Omega_\alpha$, centered at $\mathbf{x}_\alpha$ and of radius $a_\alpha$. We then have

$$|Q| \tilde{\chi}_{\alpha,n} = \int_{|\mathbf{x}-\mathbf{x}_\alpha| \le a_\alpha} \exp(-i \mathbf{k}_n \cdot \mathbf{x}) \, d\mathbf{x} = \int_{|\mathbf{r}| \le a_\alpha} \exp(-i \mathbf{k}_n \cdot (\mathbf{r} + \mathbf{x}_\alpha)) \, d\mathbf{r} = e^{-i \mathbf{k}_n \cdot \mathbf{x}_\alpha} \int_{|\mathbf{r}| \le a_\alpha} \exp(-i \mathbf{k}_n \cdot \mathbf{r}) \, d\mathbf{r}.$$

Recalling the notation $k_n = |k_n|$, we deduce that

$$|Q| \tilde{\chi}_{\alpha,n} = e^{-i \mathbf{k}_n \cdot \mathbf{x}_\alpha} \int_0^{a_\alpha} \int_0^{2\pi} \exp(-i \mathbf{k}_n r \cos \theta) r d\theta dr = 2\pi e^{-i \mathbf{k}_n \cdot \mathbf{x}_\alpha} \int_0^{a_\alpha} J_0(\mathbf{k}_n r) r dr \\ = \frac{2\pi}{k_n^2} e^{-i \mathbf{k}_n \cdot \mathbf{x}_\alpha} \int_0^{k_n a_\alpha} J_0(\xi) \xi d\xi = \frac{2\pi a_\alpha}{k_n} e^{-i \mathbf{k}_n \cdot \mathbf{x}_\alpha} J_1(\mathbf{k}_n a_\alpha),$$

where, in the above derivation, $J_0$ and $J_1$ are the Bessel functions of the first kind defined by

$$J_0(r) = \frac{1}{2\pi} \int_{-\pi}^{\pi} \exp(-i r \sin \theta) d\theta, \quad \text{and} \quad J_1(r) = \frac{1}{2\pi} \int_{-\pi}^{\pi} \exp[i(\theta - r \sin \theta)] d\theta.$$

Note that, in the above derivation, we have used Eqs. 10.9.1 and 10.6.6 from the *Digital Library of Mathematical Functions* [55]. We have thus retrieved the expression (20).

## E. On the influence tensor (BCs at infinity)

In the present appendix, we restrict ourselves to the two-dimensional setting and we derive the closed form expression of the influence tensor $\Gamma_{12}^{\infty}$ (defined in Sec. 4.2, see in particular (25)) of two circular inclusions $\Omega_1$ and $\Omega_2$ with radii $a_1$ and $a_2$, centered at $\mathbf{x}_1$ and $\mathbf{x}_2 = \mathbf{x}_1 + \mathbf{r}_{12}$ (see Fig. E.10). We also derive the expression of the self-influence tensor $\Gamma_{11}^{\infty}$. The two inclusions $\Omega_1$ and $\Omega_2$ are embedded in a homogeneous, isotropic, infinite matrix with conductivity $\sigma_0 = \sigma_0 \mathbf{I}$. We also assume that the two inclusions do not overlap: $|\mathbf{r}_{12}| \ge a_1 + a_2$.

### E.1. Analytical solution

Inclusion $\Omega_2$ is subjected to the (constant) polarization $\mathbf{P}$. We seek the (average) induced electric field in $\Omega_1$. The problem at hand is therefore (see (21)–(24))

$$\operatorname{div} \mathbf{J} = 0 \quad \text{in } \mathbb{R}^2, \qquad (E.1)$$

$$\mathbf{J} = \sigma_0 \mathbf{E} + \mathbf{P} \quad \text{in } \Omega_2, \qquad (E.2)$$

$$\mathbf{J} = \sigma_0 \mathbf{E} \quad \text{in } \mathbb{R}^2 \setminus \Omega_2, \qquad (E.3)$$

$$\mathbf{E} = -\operatorname{grad} \Phi \quad \text{in } \mathbb{R}^2. \qquad (E.4)$$
The boundary conditions are encoded in the fact that the scalar field $\Phi$ (which belongs to $L_{\text{loc}}^2(\mathbb{R}^2)$) is such that $\mathbf{E}$ (and therefore $\mathbf{J}$) has components that are square integrable over $\mathbb{R}^2$: $\mathbf{E} \in (L^2(\mathbb{R}^2))^2$. Introducing the center-to-center distance $r_{12} = |\mathbf{r}_{12}|$ and the unit vectors $\mathbf{e}_1 = \mathbf{r}_{12}/r_{12}$ and $\mathbf{e}_2$ such that $\angle(\mathbf{e}_1, \mathbf{e}_2) = \pi/2$ (see Fig. E.10), the solution to problem (E.1)–(E.4) is classically found in cylindrical coordinates $(r, \theta)$ with pole $\mathbf{x}_2$ and polar axis $\mathbf{e}_1$:

$$ \mathbf{x} = \mathbf{x}_2 + r \cos \theta \mathbf{e}_1 + r \sin \theta \mathbf{e}_2. $$

The (unique up to the addition of a constant) electric potential is seeked in the form

$$ \Phi(\mathbf{x}) = f(r) \cos \theta + g(r) \sin \theta, $$

and it can be shown that $f$ and $g$ both satisfy the same differential equation,

$$ r^2 f'' + r f' - f = 0 \quad \text{and} \quad r^2 g'' + r g' - g = 0, $$

for $r \in (0, a_2)$ and for $r \in (a_2, +\infty)$. This implies that, on each of the two intervals $(0, a_2)$ and $(a_2, +\infty)$, $f$ and $g$ are linear combinations of $r$ and $1/r$. Since $\Phi$ should be square integrable in the neighborhood of the origin, and since $\nabla\Phi$ should belong to $(L^2(\mathbb{R}^2 \setminus \Omega_2))^2$, we obtain that

$$ f(r) = C_{f,\text{int}} \frac{r}{a_2} \quad \text{on } (0, a_2), \qquad f(r) = C_{f,\text{ext}} \frac{a_2}{r} \quad \text{on } (a_2, +\infty), $$

and likewise for $g$. The integration constants $C_{f,\text{int}}$, $C_{f,\text{ext}}$, $C_{g,\text{int}}$ and $C_{g,\text{ext}}$ are found from the continuity of the electric potential (which implies that $C_{f,\text{int}} = C_{f,\text{ext}}$ and $C_{g,\text{int}} = C_{g,\text{ext}}$) and of the radial flux, leading to

$$ \Phi(\mathbf{x}) = \begin{cases} \frac{1}{2\sigma_0} \mathbf{P} \cdot (\mathbf{x} - \mathbf{x}_2) & \text{if } |\mathbf{x} - \mathbf{x}_2| \le a_2, \\ \frac{a_2^2}{2\sigma_0} \frac{\mathbf{P} \cdot (\mathbf{x} - \mathbf{x}_2)}{|\mathbf{x} - \mathbf{x}_2|^2} & \text{otherwise.} \end{cases} \qquad (\text{E.5}) $$

This expression is next used to evaluate the average electric field over the two inclusions.

### E.2. Evaluation of the influence tensor $\mathbf{\Gamma}_{12}^{\infty}$

By definition, $-\mathbf{\Gamma}_{12}^{\infty} \mathbf{P}$ is the average of the electric field $\mathbf{E} = -\mathbf{grad} \Phi$ over the inclusion $\Omega_1$. Since the inclusions do not overlap, the electric potential in inclusion $\Omega_1$ reads (see (E.5))

$$ \Phi(\mathbf{x}) = \frac{a_2^2}{2\sigma_0} \frac{\mathbf{P} \cdot (\mathbf{x} - \mathbf{x}_2)}{|\mathbf{x} - \mathbf{x}_2|^2}. \qquad (\text{E.6}) $$

We use Green's formula to transform the average over $\Omega_1$ into an integral over the boundary of the inclusion:

$$ \pi a_1^2 \mathbf{\Gamma}_{12}^{\infty} \mathbf{P} = \int_{\Omega_1} \mathbf{grad} \Phi = \int_{\partial\Omega_1} \Phi(\mathbf{x}') \mathbf{n}_1(\mathbf{x}') \, \mathrm{d}\mathbf{s}, \qquad (\text{E.7}) $$

where $\mathbf{n}_1$ denotes the outer unit normal to $\partial\Omega_1$ and $\mathbf{s}$ is the arc-length measured on the boundary $\partial\Omega_1$. Inserting (E.6) in (E.7), we find

$$ \mathbf{\Gamma}_{12}^{\infty} = \frac{1}{2\pi\sigma_0} \frac{a_2^2}{a_1^2} \int_{\mathbf{x}' \in \partial\Omega_1} \mathbf{n}_1(\mathbf{x}') \otimes \frac{\mathbf{x}' - \mathbf{x}_2}{|\mathbf{x}' - \mathbf{x}_2|^2} \, \mathrm{d}\mathbf{s}. $$

The above integral is evaluated in polar coordinates with pole $\mathbf{x}_1$ and polar axis $\mathbf{e}_1$:

$$ \mathbf{x}' = \mathbf{x}_1 + a_1(\cos \varphi \mathbf{e}_1 + \sin \varphi \mathbf{e}_2), $$

with $s = a_1 \varphi$ and $\mathbf{n}_1(\mathbf{x}') = \cos \varphi \mathbf{e}_1 + \sin \varphi \mathbf{e}_2$. Then

$$ \mathbf{x}' - \mathbf{x}_2 = a_1(\cos \varphi \mathbf{e}_1 + \sin \varphi \mathbf{e}_2) - r_{12} \mathbf{e}_1 $$

and, upon integration (see Appendix G for details), we obtain

$$ \mathbf{\Gamma}_{12}^{\infty} = \frac{1}{2\pi\sigma_0} \frac{a_2^2}{a_1} \int_0^{2\pi} [\cos \varphi \mathbf{e}_1 + \sin \varphi \mathbf{e}_2] \otimes \frac{a_1(\cos \varphi \mathbf{e}_1 + \sin \varphi \mathbf{e}_2) - r_{12} \mathbf{e}_1}{|a_1(\cos \varphi \mathbf{e}_1 + \sin \varphi \mathbf{e}_2) - r_{12} \mathbf{e}_1|^2} \mathrm{d}\varphi = \frac{1}{2\sigma_0} \frac{a_2^2}{r_{12}^2} (\mathbf{e}_2 \otimes \mathbf{e}_2 - \mathbf{e}_1 \otimes \mathbf{e}_1), \quad (\text{E.8}) $$

which delivers the intrinsic expression (26).
<table><caption>Table F.3: Parameters used for the numerical validation of (29) (see Appendix F).</caption><thead><tr><th>Symbol</th><th>Value</th><th>Comment</th></tr></thead><tbody><tr><td><i>L</i></td><td>1.0</td><td>Size of the square unit-cell</td></tr><tr><td><i>σ</i><sub>0</sub></td><td>1.0</td><td>Uniform conductivity</td></tr><tr><td><b>x</b><sub>1</sub></td><td>(0.5<i>L</i>, 0.5<i>L</i>)</td><td>Center of disk 1</td></tr><tr><td><i>a</i><sub>1</sub></td><td>0.35<i>L</i></td><td>Radius of disk 1</td></tr><tr><td><b>x</b><sub>2</sub></td><td>(0.1<i>L</i>, 0.06<i>L</i>)</td><td>Center of disk 2</td></tr><tr><td><i>a</i><sub>2</sub></td><td>0.2<i>L</i></td><td>Radius of disk 2</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>Figure F.11: The configuration considered for the numerical validation of (29) (see Appendix F). The mesh-size is $h = L/16$ in the case shown on the figure.</figcaption></figure>

### E.3. Evaluation of the self-influence tensor $\Gamma_{22}^{\infty}$

Inside the inclusion $\Omega_2$, we infer from (E.5) that the electric field is constant and reads $\mathbf{E} = -\frac{1}{2\sigma_0} \mathbf{P}$, which implies that $\Gamma_{22}^{\infty} = \frac{1}{2\sigma_0} \mathbf{I}$.

## F. Numerical validation of formula (29)

Note: the following numerical validation can be reproduced from the dataset at the Recherche Data Gouv open data repository [29].

In the relation (29), we have proposed a correction to the expression (27) initially introduced by Zecevic and Lebensohn [27]. The validity of this correction has not been proved mathematically. However, we provide in this appendix a numerical validation (in dimension $d = 2$) of this formula for the set of parameters displayed in Table F.3. The resulting configuration is represented in Fig. F.11.

The left-hand side of (29) (“true” periodic influence tensor) is evaluated by means of a finite element simulation (see Fig. F.11) with increasing element fineness (from $L/4$ to $L/512$).
<figure><img src="image_14.png" /><figcaption>Figure F.12: FEM approximation of the (1, 1) coefficient of the influence tensor $\Gamma_{12}$, as a function of the mesh size $h$. To improve visibility, we plot $100\Gamma_{12,11}$ rather than $\Gamma_{12,11}$.</figcaption></figure>

The variational problems to be solved are: find $\Phi_i$, $i = 1, 2$, such that, for all $\Psi$,

$$ \int_Q \mathbf{grad} \Phi_i \cdot \mathbf{grad} \Psi = \int_{\Omega_2 \cap Q} \partial_i \Psi, $$

where $\Phi_1$, $\Phi_2$ and $\Psi$ are $Q$-periodic. The coefficients of $\Gamma_{12}$ are then computed as

$$ \Gamma_{12,ji} = \mathbf{e}_j \cdot \Gamma_{12} \mathbf{e}_i = \frac{1}{\sigma_0 |\Omega_1|} \int_{\Omega_1} \partial_j \Phi_i. $$

The values of $\Gamma_{12,11}$ are plotted in Fig. F.12 versus the element-size $h$. Convergence as $h \to 0$ seems to be observed. This is confirmed by Fig. F.13, which represents the relative error on $\Gamma_{12,11}$, defined by

$$ \left| \frac{\Gamma_{12,11}(h) - \Gamma_{12,11}(h_{\min})}{\Gamma_{12,11}(h_{\min})} \right|, $$

where the reference value is taken as the FE result for the finest mesh size considered here, that is $h_{\min} = L/512$. The figure clearly shows the expected $h^2$ convergence rate.

We now turn to the right-hand side of (29) (Poisson summation). This right-hand side is estimated for various values of $P$ (ranging from 2 to 4096), using (26). The obtained values of the (1, 1) coefficient of the influence tensor are plotted in Fig. F.14. Again, convergence seems to be observed as $P \to +\infty$. This is confirmed by Fig. F.15, which represents the relative error on $\Gamma_{12,11}$, now defined by

$$ \left| \frac{\Gamma_{12,11}(P) - \Gamma_{12,11}(P_{\max})}{\Gamma_{12,11}(P_{\max})} \right|, $$

where the reference value is taken as the truncated sum for the largest value of $P$ considered here, that is $P_{\max} = 4096$. The truncated sum seems to converge as $P^{-2}$.

The FEM (resp. Poisson) value of the complete $2 \times 2$ matrix $\Gamma_{12}$ for the smallest grid-size $h$ (resp. for the largest truncation parameter $P$) is

$$ \Gamma_{12} = \begin{bmatrix} -0.0567582 & -0.0113464 \\ -0.0113464 & -0.0689035 \end{bmatrix} \quad (\text{FEM}), \quad \text{and} \quad \Gamma_{12} = \begin{bmatrix} -0.0567591 & -0.0113466 \\ -0.0113466 & -0.0689045 \end{bmatrix} \quad (\text{Poisson}). $$

Both results are in excellent agreement, which provides a first validation of the proposed corrected formula (29).
<figure><img src="image_15.png" /><figcaption>Figure F.13: Relative error on the (1, 1) coefficient of the influence tensor $\Gamma_{12}$ (when computed by FEM), as a function of the mesh size $h$. The straight line represents the convergence rate $h^2$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure F.14: Estimates of the (1, 1) coefficient of the influence tensor $\Gamma_{12}$ based on the Poisson summation formula, as a function of the truncation parameter $P$. As in Fig. F.12, to improve visibility, we plot $100\Gamma_{12,11}$ rather than $\Gamma_{12,11}$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure F.15: Relative error on the (1, 1) coefficient of the influence tensor $\Gamma_{12}$ (when computed by the Poisson summation formula), as a function of the truncation parameter $P$. The straight line represents the convergence rate $P^{-2}$.</figcaption></figure>
## G. Details on formula (E.8)

We show in this appendix how to obtain the last equality in (E.8). It is first observed that

$$|a_1(\cos \varphi \mathbf{e}_1 + \sin \varphi \mathbf{e}_2) - r_{12} \mathbf{e}_1|^2 = (a_1 \cos \varphi - r_{12})^2 + a_1^2 \sin^2 \varphi = (a_1^2 + r_{12}^2) \left(1 - \frac{2a_1 r_{12}}{a_1^2 + r_{12}^2} \cos \varphi\right), \quad (G.1)$$

with $r_{12} \ge a_1$ since the two disks do not overlap. We introduce the angle $\psi \in [0, \pi/2]$ defined by $\cos \frac{\psi}{2} = \frac{r_{12}}{\sqrt{a_1^2 + r_{12}^2}}$ and $\sin \frac{\psi}{2} = \frac{a_1}{\sqrt{a_1^2 + r_{12}^2}}$. Eq. (G.1) then reads $|a_1(\cos \varphi \mathbf{e}_1 + \sin \varphi \mathbf{e}_2) - r_{12} \mathbf{e}_1|^2 = (a_1^2 + r_{12}^2)(1 - \sin \psi \cos \varphi)$, which delivers, in view of (E.8), the following expression of the influence tensor:

$$\mathbf{\Gamma}_{12}^{\infty} = \frac{1}{2\pi\sigma_0} \frac{a_1^2}{a_1(a_1^2 + r_{12}^2)} \int_0^{2\pi} [\cos \varphi \mathbf{e}_1 + \sin \varphi \mathbf{e}_2] \otimes \frac{a_1(\cos \varphi \mathbf{e}_1 + \sin \varphi \mathbf{e}_2) - r_{12} \mathbf{e}_1}{1 - \sin \psi \cos \varphi} d\varphi.$$

The influence tensor therefore appears as a linear combination of integrals of the form

$$I_{mn}(\psi) = \int_{0}^{2\pi} \frac{\sin^m \varphi \cos^n \varphi}{1 - \sin \psi \cos \varphi} d\varphi. \quad (G.2)$$

The integral (G.2) is split in four integrals over the intervals $(0, \pi/2)$, $(\pi/2, \pi)$, $(\pi, 3\pi/2)$ and $(3\pi/2, 2\pi)$. Using the change of variable $\varphi \mapsto \pi + \varphi$ in the second and third intervals and $\varphi \mapsto 2\pi + \varphi$ in the fourth interval, we obtain

$$I_{mn}(\psi) = \int_{-\pi/2}^{\pi/2} \frac{\sin^m \varphi \cos^n \varphi}{1 - \sin \psi \cos \varphi} d\varphi + (-1)^{m+n} \int_{-\pi/2}^{\pi/2} \frac{\sin^m \varphi \cos^n \varphi}{1 + \sin \psi \cos \varphi} d\varphi.$$

If $m$ is odd, the integral vanishes since the integrand is an odd function. We thus conclude that $I_{11}(\psi) = I_{10}(\psi) = 0$, which yields

$$\mathbf{e}_1 \cdot \mathbf{\Gamma}_{12}^{\infty} \mathbf{e}_2 = \mathbf{e}_2 \cdot \mathbf{\Gamma}_{12}^{\infty} \mathbf{e}_1 = 0. \quad (G.3)$$

To proceed, we introduce the new integration variable $t = \tan \frac{\varphi}{2}$ and write that

$$I_{mn}(\psi) = \frac{2^{m+1}}{1 + \sin \psi} \int_{-1}^{1} \frac{t^m (1 - t^2)^n}{(1 + t^2)^{m+n} (\alpha^2 + t^2)} dt + \frac{2^{m+1} (-1)^{m+n}}{1 - \sin \psi} \int_{-1}^{1} \frac{t^m (1 - t^2)^n}{(1 + t^2)^{m+n} (\beta^2 + t^2)} dt,$$

where we have introduced $\alpha = \sqrt{\frac{1 - \sin \psi}{1 + \sin \psi}}$ and $\beta = \sqrt{\frac{1 + \sin \psi}{1 - \sin \psi}}$. Note that $\alpha\beta = 1$, $\tan(\pi/4 - \psi/2) = \frac{1 - \tan \psi/2}{1 + \tan \psi/2} = \alpha$ and $\tan(\pi/4 + \psi/2) = \frac{1 + \tan \psi/2}{1 - \tan \psi/2} = \beta$, which can be recast as $\arctan \alpha = \frac{\pi}{4} - \frac{\psi}{2}$ and $\arctan \beta = \frac{\pi}{4} + \frac{\psi}{2}$.

### G.1. The case $m = 0, n = 0$

The term $I_{00}(\psi)$ does not appear in the expansion of $\mathbf{\Gamma}_{12}^{\infty}$, but its computation is useful. In view of the identity $\int_{-1}^{1} \frac{1}{\alpha^2 + t^2} dt = \frac{2}{\alpha}(\frac{\pi}{4} + \frac{\psi}{2})$, we obtain that $I_{00}(\psi) = \frac{2\pi}{\cos \psi}$.

### G.2. The case $m = 0, n = 1$

It is easy to see that

$$I_{01}(\psi) = \frac{2}{1 + \sin \psi} \int_{-1}^{1} \left( \frac{c_1(\alpha)}{1 + t^2} + \frac{c_2(\alpha)}{\alpha^2 + t^2} \right) dt - \frac{2}{1 - \sin \psi} \int_{-1}^{1} \left( \frac{c_1(\beta)}{1 + t^2} + \frac{c_2(\beta)}{\beta^2 + t^2} \right) dt,$$

with $c_1(\alpha) = \frac{2}{\alpha^2 - 1} = -\frac{1 + \sin \psi}{\sin \psi}$ and $c_2(\alpha) = -1 - c_1(\alpha) = \frac{1}{\sin \psi}$ on the one hand, and $c_1(\beta) = \frac{2}{\beta^2 - 1} = \frac{1 - \sin \psi}{\sin \psi}$ and $c_2(\beta) = -1 - c_1(\beta) = -\frac{1}{\sin \psi}$ on the other hand. We then deduce that $I_{01}(\psi) = -\frac{2\pi}{\sin \psi} + \frac{2\pi}{\sin \psi \cos \psi}$.
### G.3. Trace of the influence tensor

Using the fact that $I_{02}(\psi) = I_{00}(\psi) - I_{20}(\psi)$, we see that $\text{tr} \Gamma_{12}^{\infty} = \frac{1}{2\pi\sigma_0} \frac{a_2^2}{a_1(a_1^2 + r_{12}^2)} (a_1 I_{00}(\psi) - r_{12} I_{01}(\psi))$. Using the above results, we obtain

$$ \frac{a_1}{r_{12}} I_{00}(\psi) - I_{01}(\psi) = \frac{a_1}{r_{12}} \frac{2\pi}{\cos \psi} + \frac{2\pi}{\sin \psi} - \frac{2\pi}{\sin \psi \cos \psi} = \tan \frac{\psi}{2} \frac{2\pi}{\cos \psi} + \frac{2\pi}{\sin \psi} - \frac{2\pi}{\sin \psi \cos \psi} = 0, $$

where we have used the generic relation $\tan \frac{\psi}{2} = \frac{\sin \psi}{1 + \cos \psi}$ to conclude. This directly implies that $\text{tr} \Gamma_{12}^{\infty} = 0$.

### G.4. The case $m = 2, n = 0$

We write that

$$ I_{20}(\psi) = \frac{8}{1 + \sin \psi} \int_{-1}^{1} \left( \frac{c_3(\alpha)}{1 + t^2} + \frac{c_4(\alpha)}{(1 + t^2)^2} + \frac{c_5(\alpha)}{\alpha^2 + t^2} \right) dt + \frac{8}{1 - \sin \psi} \int_{-1}^{1} \left( \frac{c_3(\beta)}{1 + t^2} + \frac{c_4(\beta)}{(1 + t^2)^2} + \frac{c_5(\beta)}{\beta^2 + t^2} \right) dt, $$

with $c_3(\alpha) = \frac{\alpha^2}{(1 - \alpha^2)^2} = \frac{\cos^2 \psi}{4 \sin^2 \psi}$, $c_4(\alpha) = \frac{1}{1 - \alpha^2} = \frac{1 + \sin \psi}{2 \sin \psi}$ and $c_5(\alpha) = -\frac{\cos^2 \psi}{4 \sin^2 \psi}$ on the one hand, and $c_3(\beta) = \frac{\beta^2}{(1 - \beta^2)^2} = \frac{\cos^2 \psi}{4 \sin^2 \psi}$, $c_4(\beta) = \frac{1}{1 - \beta^2} = \frac{\sin \psi - 1}{2 \sin \psi}$ and $c_5(\beta) = -\frac{\cos^2 \psi}{4 \sin^2 \psi}$ on the other hand. We then deduce that $I_{20}(\psi) = \frac{2\pi}{\sin^2 \psi} - \frac{2\pi \cos \psi}{\sin^2 \psi}$.

### G.5. Diagonal terms of the influence tensor

We are now in position to identify the diagonal terms of the influence tensor. Since its trace vanishes, we have

$$ -\mathbf{e}_1 \cdot \Gamma_{12}^{\infty} \mathbf{e}_1 = \mathbf{e}_2 \cdot \Gamma_{12}^{\infty} \mathbf{e}_2 = \frac{1}{2\pi\sigma_0} \frac{a_2^2}{a_1(a_1^2 + r_{12}^2)} a_1 I_{20}(\psi) = \frac{1}{2\sigma_0} \frac{a_2^2}{a_1^2 + r_{12}^2} \frac{1}{\cos^2 \psi / 2} = \frac{1}{2\sigma_0} \frac{a_2^2}{r_{12}^2}, $$

where the last identity stems from the definition of $\psi$. Collecting the above relation with (G.3), we obtain (E.8).

## Acknowledgements

This work has benefited from the French government grant ANR-11-LABX-022-01 (Labex MMCD, Multi-Scale Modelling & Experimentation of Materials for Sustainable Construction) managed by ANR within the frame of the national program Investments for the Future. The authors also thank the anonymous referees for their constructive comments.

## Software

The meshes required by the full-field simulations have been generated with the Gmsh software [56], using constructive geometry primitives provided by the OpenCascade geometry kernel<fnref n="2" />. The FreeFem++ software [57] has been used for all FEM simulations. All other computations have been performed with the Julia programming language [58].

<footnote n="2">²https://dev.opencascade.org/, last retrieved 2022-07-27.</footnote>
References

[1] A. Zaoui, Continuum Micromechanics: Survey, Journal of Engineering Mechanics 128 (8) (2002) 808-816.
doi:10.1061/(ASCE)0733-9399(2002)128:8(808).

[2] P. Suquet, Effective Properties of Nonlinear Composites, in: P. Suquet (Ed.), Continuum Micromechanics, International Centre for Mechanical Sciences, Springer Vienna, Vienna, 1997, pp. 197–264. doi:10.1007/978-3-7091-2662-2_4.

[3] G. Allaire, Shape Optimization by the Homogenization Method, Springer, 2002. doi:10.1007/978-1-4684-9286-6.

[4] A. Bensoussan, J.-L. Lions, G. Papanicolaou, Asymptotic Analysis for Periodic Structures, American Mathematical Society, 1978.

[5] G. W. Milton, The Theory of Composites, no. 6 in Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, Cambridge; New York, 2002. doi:10.1017/CBO9780511613357.

[6] E. Sanchez-Palencia, Non-Homogeneous Media and Vibration Theory, no. 127 in Lecture Notes in Physics, Springer-Verlag, Berlin Heidelberg, 1980. doi:10.1007/3-540-10000-8.

[7] J. Sanchez-Hubert, E. Sanchez-Palencia, Introduction aux méthodes asymptotiques et à l'homogénéisation. Application à la mécanique des milieux continus, Masson, 1992.

[8] V. V. Jikov, S. M. Kozlov, O. A. Oleinik, Homogenization of Differential Operators and Integral Functionals, Springer-Verlag, Berlin Heidelberg, 1994. doi:10.1007/978-3-642-84659-5.

[9] M. Ostoja-Starzewski, Material spatial randomness: From statistical to representative volume element, Probabilistic Engineering Mechanics 21 (2) (2006) 112–132. doi:10.1016/j.probengmech.2005.07.007.

[10] T. Kanit, S. Forest, I. Galliet, V. Mounoury, D. Jeulin, Determination of the size of the representative volume element for random composites: Statistical and numerical approach, International Journal of Solids and Structures 40 (13–14) (2003) 3647–3679. doi:10.1016/S0020-7683(03)00143-4.

[11] J. Nolen, Normal approximation for a random elliptic equation, Probab. Theory Related Fields 159 (3–4) (2014) 661–700. doi:10.1007/s00440-013-0517-9.

[12] A. Gloria, F. Otto, An optimal variance estimate in stochastic homogenization of discrete elliptic equations, Ann. Prob. 39 (2011) 779–856. doi:10.1214/10-AOP571.

[13] G. S. Fishman, Monte Carlo: concepts, algorithms, and applications, Springer, 1996.

[14] X. Blanc, R. Costaouec, C. Le Bris, F. Legoll, Variance reduction in stochastic homogenization using antithetic variables, Markov Process. Relat. Fields 18 (2012) 31–66, (preliminary version available at http://cermics.enpc.fr/~legoll/blanc_costaouec_lebris_legoll_MPRF_2012.pdf).

[15] R. Costaouec, C. Le Bris, F. Legoll, Variance reduction in stochastic homogenization: proof of concept, using antithetic variables, Bol. Soc. Esp. Mat. Apl. 50 (2010) 9–27. doi:10.1007/BF03322539.

[16] X. Blanc, R. Costaouec, C. Le Bris, F. Legoll, Variance reduction in stochastic homogenization: the technique of antithetic variables, in: B. Engquist, O. Runborg, R. Tsai (Eds.), Numerical Analysis of Multiscale Computations, Vol. 82 of Lect. Notes Comput. Sci. Eng., Springer, 2012, pp. 47–70.

[17] C. Le Bris, F. Legoll, W. Minvielle, Special Quasirandom Structures: a selection approach for stochastic homogenization, Monte Carlo Methods and Applications 22 (1) (2016) 25–54. doi:10.1515/mcma-2016-0101.

[18] F. Legoll, W. Minvielle, A Control Variate Approach Based on a Defect-Type Theory for Variance Reduction in Stochastic Homogenization, Multiscale Modeling & Simulation 13 (2) (2015) 519–550. doi:10.1137/140980120.
[19] F. Legoll, W. Minvielle, Variance reduction using antithetic variables for a nonlinear convex stochastic homogenization problem, Discrete and Continuous Dynamical Systems – Series S 8 (1) (2015) 1–27. doi: 10.3934/dcdss.2015.8.1.

[20] X. Blanc, C. Le Bris, F. Legoll, Some variance reduction methods for numerical stochastic homogenization, Phil. Trans. R. Soc. A 374 (2066) (2016) 20150168. doi:10.1098/rsta.2015.0168.

[21] C. Le Bris, F. Legoll, Examples of computational approaches for elliptic, possibly multiscale PDEs with random inputs, J. Comput. Phys. 328 (2017) 455–473. doi:10.1016/j.jcp.2016.10.027.

[22] Z. A. Moschovidis, T. Mura, Two-Ellipsoidal Inhomogeneities by the Equivalent Inclusion Method, Journal of Applied Mechanics 42 (4) (1975) 847–852. doi:10.1115/1.3423718.

[23] J. Korringa, Theory of elastic constants of heterogeneous media, Journal of Mathematical Physics 14 (4) (1973) 509–513. doi:doi:10.1063/1.1666346.

[24] R. Zeller, P. H. Dederichs, Elastic Constants of Polycrystals, Physica Status Solidi (B) 55 (2) (1973) 831–842. doi:10.1002/pssb.2220550241.

[25] E. Kröner, On the Physics and Mathematics of Self-Stresses, in: J. L. Zeman, F. Ziegler (Eds.), Topics in Applied Continuum Mechanics, Springer Verlag, Wien, 1974, pp. 22–38.

[26] S. Brisard, L. Dormieux, K. Sab, A variational form of the equivalent inclusion method for numerical homogenization, International Journal of Solids and Structures 51 (3–4) (2014) 716–728. doi:10.1016/j.ijsolstr.2013.10.037.

[27] M. Zecevic, R. A. Lebensohn, Approximation of periodic Green’s operator in real space using numerical integration and its use in fast Fourier transform-based micromechanical models, International Journal for Numerical Methods in Engineering 122 (24) (2021) 7536–7552. doi:10.1002/nme.6841.

[28] S. Brisard, F. Legoll, M. Bertin, A variance reduction strategy for numerical random homogenization based on the equivalent inclusion method — Supporting data (2023). doi:10.57745/QRTYI6.

[29] S. Brisard, F. Legoll, M. Bertin, Evaluation of the influence tensor of two disks in a periodic setting through the Poisson summation formula (2023). doi:10.57745/EKZGIL.

[30] X. Blanc, C. Le Bris, Homogenization Theory for Multiscale Problems, no. 21 in Modeling, Simulation and Applications, Springer, 2023.

[31] S. Brisard, M. Schneider, F. Willot, Introduction to FFT-based numerical methods for the homogenization of random materials, (HAL preprint, available at https://hal.science/hal-03631887) (2022).

[32] R. Hill, The Elastic Behaviour of a Crystalline Aggregate, Proceedings of the Physical Society. Section A 65 (5) (1952) 349. doi:10.1088/0370-1298/65/5/307.

[33] R. Hill, Elastic properties of reinforced solids: Some theoretical principles, Journal of the Mechanics and Physics of Solids 11 (5) (1963) 357–372. doi:10.1016/0022-5096(63)90036-X.

[34] A. Anantharaman, R. Costaouec, C. Le Bris, F. Legoll, F. Thomines, Introduction to numerical stochastic homogenization and the related computational challenges: some recent developments, in: W. Bao, Q. Du (Eds.), Multiscale modeling and analysis for materials simulation, Vol. 22 of Lect. Notes Series, Institute for Mathematical Sciences, National University of Singapore, World Sci. Publ., Hackensack, NJ, 2011, pp. 197–272. doi:10.1142/9789814360906_0004.

[35] K. Sab, On the homogenization and the simulation of random materials, European Journal of Mechanics – A/Solids 11 (5) (1992) 585–607.
[36] A. Bourgeat, A. Piatnitski, Approximation of effective coefficients in stochastic homogenization, Ann. I. H. Poincaré – PR 40 (2) (2004) 153–165. doi:10.1016/j.anihpb.2003.07.003.

[37] C. Huet, Application of variational concepts to size effects in elastic heterogeneous bodies, Journal of the Mechanics and Physics of Solids 38 (6) (1990) 813–841. doi:10.1016/0022-5096(90)90041-2.

[38] K. Sab, B. Nedjar, Periodization of random media and representative volume element size for linear composites, Comptes Rendus Mécanique 333 (2) (2005) 187–195. doi:10.1016/j.crme.2004.10.003.

[39] J. D. Eshelby, The Determination of the Elastic Field of an Ellipsoidal Inclusion, and Related Problems, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 241 (1226) (1957) 376–396. doi:10.1098/rspa.1957.0133.

[40] Z. Hashin, S. Shtrikman, On some variational principles in anisotropic and nonhomogeneous elasticity, Journal of the Mechanics and Physics of Solids 10 (4) (1962) 335–342. doi:10.1016/0022-5096(62)90004-2.

[41] Z. Hashin, S. Shtrikman, A variational approach to the theory of the elastic behaviour of polycrystals, Journal of the Mechanics and Physics of Solids 10 (4) (1962) 343–352. doi:10.1016/0022-5096(62)90005-4.

[42] P. Debye, H. R. Anderson, H. Brumberger, Scattering by an Inhomogeneous Solid. II. The Correlation Function and Its Application, Journal of Applied Physics 28 (6) (1957) 679–683. doi:10.1063/1.1722830.

[43] Q.-D. To, G. Bonnet, D.-H. Hoang, Explicit effective elasticity tensors of two-phase periodic composites with spherical or ellipsoidal inclusions, International Journal of Solids and Structures 94–95 (2016) 100–111. doi:10.1016/j.ijsolstr.2016.05.005.

[44] H. Moulinec, P. Suquet, A fast numerical method for computing the linear and nonlinear mechanical properties of composites, Comptes Rendus de l'Académie des Sciences. Série II, Mécanique, physique, chimie, astronomie 318 (11) (1994) 1417–1423.

[45] H. Moulinec, P. Suquet, A numerical method for computing the overall response of nonlinear composites with complex microstructure, Computer Methods in Applied Mechanics and Engineering 157 (1–2) (1998) 69–94. doi:10.1016/S0045-7825(97)00218-1.

[46] A. Molinari, M. El Mouden, The problem of elastic inclusions at finite concentration, International Journal of Solids and Structures 33 (20–22) (1996) 3131–3150. doi:10.1016/0020-7683(95)00275-8.

[47] M. P. Allen, D. J. Tildesley, Computer Simulation of Liquids, Clarendon press, Oxford, 1994.

[48] A. A. Gusev, Representative volume element size for elastic composites: A numerical study, Journal of the Mechanics and Physics of Solids 45 (9) (1997) 1449–1459. doi:10.1016/S0022-5096(97)00016-1.

[49] R. Brenner, J. Bravo-Castillero, D. M. Léon, Investigation of the effective response of 2-1-2 piezoelectric composites, Procedia IUTAM 3 (2012) 292–300. doi:10.1016/j.piutam.2012.03.018.

[50] M. Schneider, On the mathematical foundations of the self-consistent clustering analysis for non-linear materials at small strains, Computer Methods in Applied Mechanics and Engineering 354 (2019) 783–801. doi:10.1016/j.cma.2019.06.003.

[51] M.-T. Nguyen, V. Monchiet, G. Bonnet, Q.-D. To, Conductivity estimates of spherical-particle suspensions based on triplet structure factors, Physical Review E 93 (2) (2016) 022105. doi:10.1103/PhysRevE.93.022105.

[52] Q.-D. To, M.-T. Nguyen, G. Bonnet, V. Monchiet, V.-T. To, Overall elastic properties of composites from optimal strong contrast expansion, International Journal of Solids and Structures 120 (2017) 245–256. doi:10.1016/j.ijsolstr.2017.05.006.
[53] S. Brisard, L. Dormieux, Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFT-based numerical method for the homogenization of composites, Computer Methods in Applied Mechanics and Engineering 217–220 (2012) 197–212. doi:10.1016/j.cma.2012.01.003.

[54] M. Schneider, Convergence of FFT-based homogenization for strongly heterogeneous media, Math. Meth. Appl. Sci. 38 (2015) 2761–2778. doi:10.1002/mma.3259.

[55] F. W. J. Olver, A. B. O. Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, M. A. McClain, NIST digital library of mathematical functions – Release 1.1.5, http://dlmf.nist.gov/ (Mar. 2022).

[56] C. Geuzaine, J.-F. Remacle, Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering 79 (11) (2009) 1309–1331. doi:10.1002/nme.2579.

[57] F. Hecht, New development in FreeFem++, Journal of Numerical Mathematics 20 (3-4) (2012) 251–265. doi:10.1515/jnum-2012-0013.

[58] J. Bezanson, A. Edelman, S. Karpinski, V. Shah, Julia: A Fresh Approach to Numerical Computing, SIAM Review 59 (1) (2017) 65–98. doi:10.1137/141000671.