# Learning Homogenization For Elliptic Operators

Kaushik Bhattacharya  California Institute of Technology, Pasadena, United
States  Nikola Kovachki  NVIDIA, Santa Clara, United States  Aakila Rajan
California Institute of Technology, Pasadena, United States  Andrew Stuart
California Institute of Technology, Pasadena, United States  Margaret Trautner
California Institute of Technology, Pasadena, United States

###### Abstract

Multiscale partial differential equations (PDEs) arise in various
applications, and several schemes have been developed to solve them
efficiently. Homogenization theory is a powerful methodology that eliminates
the small-scale dependence, resulting in simplified equations that are
computationally tractable. In the field of continuum mechanics, homogenization
is crucial for deriving constitutive laws that incorporate microscale physics
in order to formulate balance laws for the macroscopic quantities of interest.
However, obtaining homogenized constitutive laws is often challenging as they
do not in general have an analytic form and can exhibit phenomena not present
on the microscale. In response, data-driven learning of the constitutive law
has been proposed as appropriate for this task. However, a major challenge in
data-driven learning approaches for this problem has remained unexplored: the
impact of discontinuities and corner interfaces in the underlying material.
These discontinuities in the coefficients affect the smoothness of the
solutions of the underlying equations. Given the prevalence of discontinuous
materials in continuum mechanics applications, it is important to address the
challenge of learning in this context; in particular to develop underpinning
theory to establish the reliability of data-driven methods in this scientific
domain. The paper addresses this unexplored challenge by investigating the
learnability of homogenized constitutive laws for elliptic operators in the
presence of such complexities. Approximation theory is presented, and
numerical experiments are performed which validate the theory for the solution
operator defined by the cell-problem arising in homogenization for elliptic
PDEs.

##  1 Introduction

Homogenization theory is a well-established methodology that aims to eliminate
fast-scale dependence in partial differential equations (PDEs) to obtain
homogenized PDEs which produce a good approximate solution of the problem with
fast-scales while being more computationally tractable. In continuum
mechanics, this methodology is of great practical importance as the
constitutive laws derived from physical principles are governed by material
behavior at small scales, but the quantities of interest are often relevant on
larger scales. These homogenized constitutive laws often do not have a closed
analytic form and may have new features not present in the microscale laws.
Consequently, there has been a recent surge of interest in employing data-
driven methods to learn homogenized constitutive laws.

The goal of this paper is to study the learnability of homogenized
constitutive laws in the context of one of the canonical model problems of
homogenization, namely the divergence form elliptic PDE. One significant
challenge in applications of homogenization in material science arises from
the presence of discontinuities and corner interfaces in the underlying
material. This leads to a lack of smoothness in the coefficients and solutions
of the associated equations, a phenomenon extensively studied in numerical
methods for PDEs. Addressing this challenge in the context of learning remains
largely unexplored and is the focus of our work. We develop underlying theory
and provide accompanying numerical studies to address learnability in this
context.

In Subsection 1.1 we establish the mathematical framework and notation for the
problem of interest, state the three main contributions of the paper and
overview the contents of each section of the paper. In Subsection 1.2 we
provide a detailed literature review. Subsection 1.3 states the stability
estimates that are key for the approximation theory developed in the paper and
discusses the remainder of the paper in the context of these estimates.

###  1.1 Problem Formulation

We consider the following linear multiscale elliptic equation on a bounded
domain $\Omega\subset\mathbb{R}^{d}$:

$\displaystyle-\nabla\cdot\left(A^{\epsilon}\nabla u^{\epsilon}\right)$
$\displaystyle=f\quad x\in\Omega,$ (1.1a) $\displaystyle u^{\epsilon}$
$\displaystyle=0\quad x\in\partial\Omega.$ (1.1b)

Here $A^{\epsilon}(x)=A\left(\frac{x}{\epsilon}\right)$ for $A(\cdot)$ which
is $1$-periodic and positive definite: $A:\mathbb{T}^{d}\to\mathbb{R}^{d\times
d}_{{\rm sym},\succ 0}$. Our focus is on linking this multiscale problem to
the homogenized form of equation (1.1), which is

$\displaystyle-\nabla\cdot\left(\overline{A}\nabla u\right)$
$\displaystyle=f\quad x\in\Omega,$ (1.2a) $\displaystyle u$
$\displaystyle=0\quad x\in\partial\Omega,$ (1.2b)

where $\overline{A}$ is given by

$\overline{A}=\int_{\mathbb{T}^{d}}\left(A(y)+A(y)\nabla\chi(y)^{T}\right)\;\>\mathsf{d}y,$
(1.3)

and $\chi:\mathbb{T}^{d}\to\mathbb{R}^{d}$ solves the cell problem

$-\nabla_{y}\cdot(\nabla_{y}\chi A)=\nabla_{y}\cdot A,\quad\chi\;\text{is
}1\text{-periodic}.$ (1.4)

For $0<\epsilon\ll 1$, the solution $u^{\epsilon}$ of (1.1) is approximated by
the solution $u$ of (1.2) and the error converges to zero as $\epsilon\to 0$
in various topologies [5, 10, 49].

We assume $f\in L^{2}(\Omega;\mathbb{R})$ in equation (1.1) is independent of
the microscale variable $\epsilon$, and that

$\|A\|_{L^{\infty}}:=\sup_{y\in\mathbb{T}^{d}}|A(y)|_{F}<\infty$

where $|\cdot|_{F}$ is the Frobenius norm. Hence $A\in
L^{\infty}\left(\mathbb{T}^{d};\mathbb{R}^{d\times d}\right)$ and
$A^{\epsilon}\in L^{\infty}\left(\Omega;\mathbb{R}^{d\times d}\right).$
Similarly, for $A\in L^{2}(\mathbb{T}^{d};\mathbb{R}^{d\times d})$, we define

$\|A\|_{L^{2}}^{2}:=\int_{\mathbb{T}^{d}}|A(y)|_{F}^{2}\>\mathsf{d}y=\sum_{j=1}^{d}\|a_{j}\|_{L^{2}}^{2}$

where $a_{j}$ denotes the $j$-th column (or row) of $A$. We also, for given
$\beta\geq\alpha>0$, define the following subset of 1-periodic, positive-
definite matrix fields in $L^{\infty}\left(\mathbb{T}^{d};\mathbb{R}^{d\times
d}\right)$ by

$\mathsf{PD}_{\alpha,\beta}=\\{A\in
L^{\infty}(\mathbb{T}^{d};\mathbb{R}^{d\times
d}):\;\forall(y,\xi)\in\mathbb{T}^{d}\times\mathbb{R}^{d},\;\alpha|\xi|^{2}\leq\langle\xi,A(y)\xi\rangle\leq\beta|\xi|^{2}\\}.$

Finally, we often work in the Sobolev space $H^{1}$ restricted to spatially
mean-zero periodic functions, denoted

$\dot{H}^{1}:=\\{u\in W^{1,2}(\mathbb{T}^{d})\;|\;u\text{ is
}1\text{-periodic},\int_{\mathbb{T}^{d}}u\;\>\mathsf{d}y=0\\};$

the norm on this space is defined by

$\|g\|_{\dot{H}^{1}}:=\|\nabla g\|_{L^{2}}.$ (1.5)

Numerically solving (1.1) is far more computationally expensive than solving
the homogenized equation (1.2), motivating the wish to find the homogenized
coefficient $\overline{A}$ defining equation (1.2). The difficult part of
obtaining the equation (1.2) is solving the cell problem (1.4). Although
explicit solutions exist in the one-dimensional setting for piecewise constant
$A$ [9] and in the two-dimensional setting where $A$ is a layered material
[49], in general a closed form solution is not available and the cell problem
must be solved numerically. Note that in general the right hand side
$\nabla_{y}\cdot A$ of the cell-problem can only be defined weakly for
$A\notin C^{1}(\mathbb{T}^{d},\mathbb{R}^{d\times d})$, a commonly occuring
situation in applications such as those arising from porous medium flow, or to
vector-valued generalizations of the setting here to elasticity, rendering the
numerical solution non-trivial. For this reason, it is potentially valuable to
approximate the solution map

$G:\;A\mapsto\chi,$ (1.6)

defined by the cell problem, using a map defined by a neural operator. More
generally it is foundational to the broader program of learning homogenized
constitutive models from data to thoroughly study this issue for the
divergence form elliptic equation as the insights gained will be important for
understanding the learning of more complex parameterized homogenized models,
such as those arising in nonlinear elasticity, viscoleasticity and plasticity.

The full map from $A$ to the homogenized tensor $\overline{A}$ is expressed by
$A\mapsto\chi\mapsto\overline{A}$, and one could instead learn the map

$F:A\mapsto\overline{A}.$ (1.7)

Since the map $\chi\mapsto\overline{A}$ is simply a linear integration of
$\nabla\chi$, we focus on the approximation of $A\mapsto\chi$ and state
equivalent results for the map $A\mapsto\overline{A}$ that emerge as
consequences of the approximation of $\chi$. In this paper we make the
following contributions:

  1. 1.

We state and prove universal approximation theorems for the maps $G$ defined
by (1.4) and (1.6) and $F$ defined by (1.3), (1.4) and (1.7), in various
topologies and for a pair of different neural operator architectures.

  2. 2.

We provide explicit examples of microstructures which satisfy the hypotheses
of our theorems.

  3. 3.

We provide numerical experiments to demonstrate the ability of a neural
operators to approximate the solution map on three different classes of
material parameters $A$, including one class which goes beyond the confines of
the theory.

In Subsection 1.2 we provide an overview of the literature, followed in
Subsection 1.3 by a discussion of stability estimates for (1.3), with respect
to variations in $A$; these are at the heart of the analysis of universal
approximation. The main body of the text then commences with Section 2 which
characterizes the microstructures of interest to us in the context of
continuum mechanics. Section 3 states universal approximation theorems for
$G(\cdot)$ and $F(\cdot)$, using the Fourier neural operator and a more
general neural operator, respectively. In Section 4 we give numerical
experiments illustrating the approximation of map $G$ defined by (1.6) on
microstructures of interest in continuum mechanics. Details of the stability
estimates, the proofs of universal approximation theorems and properties of
the microstructures are given in Appendices A, B and C respectively.

###  1.2 Literature Review

Homogenization aims to derive macroscopic equations that describe the
effective properties and behavior of solutions to problems at larger scales
given a system that exhibits behaviour at (possibly multiple) smaller scales.
Although it is developed for the various cases of random, statistically
stationary, and periodic small scale structures, we work here entirely in the
periodic setting. The underlying assumption of periodic homogenization theory
is that the coefficient is periodic in the small scale variable, and that the
scale separation is large compared to the macroscopic scales of interest.
Convergence of the solution to the multiscale problem to the homogenized
solution is well-studied; see [2, 14]. We refer to the texts [5, 10, 49] for
more comprehensive citations to the literature. Homogenization has found
extensive application in the setting of continuum mechanics [23] where, for
many multiscale materials, the scale-separation assumption is natural. In this
work, we are motivated in part by learning constitutive models for solid
materials, where crystalline microstructure renders the material parameters
discontinuous, and may include corner interfaces. This difficulty has been
explored extensively in the context of numerical methods for PDEs,
particularly with adaptive finite element methods [28, 11, 46, 48].

There is a significant body of work on the approximation theory associated
with parametrically dependent solutions of PDEs, including viewing these
solution as a map between the function space of the parameter and the function
space of the solution, especially for problems possessing holomorphic
regularity [15, 16, 13]. This work could potentially be used to study the cell
problem for homogenization that is our focus here. However, there has been
recent interest in taking a data-driven approach to solving PDEs via machine
learning because of its flexibility and ease of implementation. A particular
approach to learning solutions to PDEs is operator learning, a machine
learning methodology where the map to be learned is viewed as an operator
acting between infinite-dimensional function spaces rather than between
finite-dimensional spaces [8, 34, 41, 45, 31]. Determining whether, and then
when, operator learning models have advantages over classical numerical
methods in solving PDEs remains an active area of research [4]. However, in
the setting of learning homogenized constitutive laws, a few characteristics
make operator learning a promising option. First, machine learning has been
groundbreaking in application settings with no clear underlying equations,
such as computer vision and language models [26, 12]. In constitutive
modeling, though the microscale constitutive laws are known, the homogenized
equations are generally unknown and can incorporate dependencies that are not
present on the microscale, such as history dependence, anisotropy, and slip-
stick behavior [50, 7]. Thus, constitutive models lie in a partially equation-
free setting where data-driven methods could be useful. Second, machine
learned models as surrogates for expensive computation can be valuable when
the cost of producing data and training the model can be amortized over many
forward uses of the trained model. Since the same materials are often used for
fabrication over long time periods, this can be a setting where the up-front
cost of data production and model training is justified.

Other work has already begun to explore the use of data-driven methods for
constitutive modeling; a general review of the problem and its challenges, in
the context of constitutive modeling of composite materials, may be found in
[39]. Several works use the popular approach of physics-informed machine
learning to approach the problem [21, 53, 42, 24]. In [3], physical
constraints are enforced on the network architecture while learning nonlinear
elastic constitutive laws. In [35], the model is given access to additional
problem-specific physical knowledge. Similarly, the work of [54] predicts the
Cholesky factor of the tangent stiffness matrix from which the stress may be
calculated; this method enforces certain physical criteria. The paper [29]
studies approximation error and uncertainty quantification for this learning
problem. In [25], a derivative-free approach is taken to learning homogenized
solutions, but this method assumes regularity of the material coefficient. In
[36], the potential of operator learning methodology to model constitutive
laws with history dependence, such as those that arise in crystal plasticity,
is illustrated. Finally, a number of further works demonstrate empirically the
potential of learning constitutive models, including [44, 40, 55, 38].

However, the underlying theory behind operator learning for constitutive
models lags behind its empirical application. In [9], approximation theories
are developed to justify the use of a recurrent Markovian architecture that
performs well in application settings with history dependence. This
architecture is further explored in [37] with more complex microstructures.
Universal approximation results are a first step in developing theory for
learning because they guarantee that there exists an $\epsilon$-approximate
operator within the operator approximation class, which is consistent with an
assumed true model underlying the data [17, 33, 31, 30]. In addition to
universal approximation, further insight may be gained by seeking to quantify
the data or model size required to obtain a given level of accuracy; the
papers [33, 30, 43] also contain work in this direction, as do the papers [27,
47], which build on the analysis developed in [15, 16, 13] referred to above.
In our work we leverage two existing universal approximation theorems for
neural operators, one from [31] for general neural operators (NOs) and one
from [30] for Fourier neural operators (FNOs), a particular practically useful
architecture from within the NO class. We take two different approaches to
proving approximation theorems based on separate PDE solution stability
results in pursuit of a more robust understanding of the learning problem.
Since the state of the field is in its infancy, it is valuable to have
different approaches to these analysis problems. Finally, we perform numerical
experiments on various microstructures to understand the practical effects of
non-smooth PDE coefficients in learning solutions. We highlight the fact that
in this paper we do not tackle issues related to the non-convex optimization
problem at the heart of training neural networks; we simply use state of the
art stochastic gradient descent for training, noting that theory explaining
its excellent empirical behaviour is lacking.

Throughout this paper we focus on equation (1.1), which describes the
constitutive law of linear elasticity. Though it is a linear elliptic
equation, we echo the sentiment of Blanc and Le Bris [10] with their warning
“do not underestimate the difficulty of equation (1.1).” There are many
effects to be understood in this setting, and resolving learning challenges is
a key step towards understanding similar questions for the learning of
parametric dependence in more complex homogenized constitutive laws where
machine-learning may prove particularly useful.

###  1.3 Stability Estimates

At the heart of universal approximation theorems is stability of the solution
map (1.6); in particular continuity of the map for certain classes of $A$. In
this subsection, we present three key stability results that are used to prove
the approximation theorems in Section 3. The proofs of the following stability
estimates may all be found in Appendix A.

A first strike at the stability of the solution map (1.6) is a modification of
the classic $L^{\infty}/H^{1}$ Lipschitz continuity result for dependence of
the solution of elliptic PDEs on the coefficient; here generalization is
necessary because the coefficient also appears on the right-hand side of the
equation defining $G(\cdot):$

######  Proposition 1.1.

Consider the cell problem defined by equation (1.4). The following hold:

  1. 1.

If $A\in\mathsf{PD}_{\alpha,\beta}$, then (1.4) has a unique solution
$\chi\in\dot{H}^{1}(\mathbb{T}^{d};\mathbb{R}^{d})$ and

$\|\chi\|_{\dot{H}^{1}(\mathbb{T}^{d};\mathbb{R}^{d})}\leq\frac{\sqrt{d}\beta}{\alpha}.$

  2. 2.

For $\chi^{(1)}$ and $\chi^{(2)}$ solutions to the cell problem in equation
(1.4) associated with coefficients
$A^{(1)},A^{(2)}\in\mathsf{PD}_{\alpha,\beta}$, respectively, it follows that

$\|\chi^{(2)}-\chi^{(1)}\|_{\dot{H}^{1}(\mathbb{T}^{d};\mathbb{R}^{d})}\leq\frac{\sqrt{d}}{\alpha}\left(1+\frac{\beta}{\alpha}\right)\|A^{(1)}-A^{(2)}\|_{L^{\infty}(\mathbb{T}^{d};\mathbb{R}^{d\times
d})}.$ (1.8)

However, this perturbation result is insufficient for approximation theory
because the space $L^{\infty}$ is not separable and approximation is hardly
possible in such spaces [19, Chapter 9]. While one may define the problem on a
separable subspace of $L^{\infty}$, see Lemma A.1, such spaces are not
particularly useful in applications to micromechanics. This is because many of
the models for realistic microstructures include functions which can be
discontinuous on an uncountable set of points in the domain and such functions
cannot be contained in any separable subspace of $L^{\infty}$; see Lemma A.2.
Instead, we require continuity from $L^{q}$ to $\dot{H}^{1}$ for some
$q\in[2,\infty)$. To this end, we provide two additional stability results.
The first stability result gives continuity, but not Lipschitz continuity,
from $L^{2}$ to $\dot{H}^{1}$. The second stability result gives Lipschitz
continuity from $L^{q}$ to $\dot{H}^{1}$ some $q\in(2,\infty)$.

######  Proposition 1.2.

Endow $\mathsf{PD}_{\alpha,\beta}$ with the
$L^{2}(\mathbb{T}^{d};\mathbb{R}^{d\times d})$ induced topology and let
$K\subset\mathsf{PD}_{\alpha,\beta}$ be a closed set. Define the mapping
$G:K\to\dot{H}^{1}(\mathbb{T}^{d};\mathbb{R}^{d})$ by $A\mapsto\chi$ as given
by (1.4). Then there exists a bounded continuous mapping $\mathcal{G}\in
C(L^{2}(\mathbb{T}^{d};\;\mathbb{R}^{d\times
d});\dot{H}^{1}(\mathbb{T}^{d};\mathbb{R}^{d}))$ such that
$\mathcal{G}(A)=G(A)$ for any $A\in K$.

The preceding $L^{2}$ continuity proposition is used to prove the
approximation results for the FNO in Theorems 3.3 and 3.4. For the second
approximation theorem, we need the following proposition on Lipschitz
continuity from $L^{q}$ to $\dot{H}^{1}$.

######  Proposition 1.3.

If $A\in\mathsf{PD}_{\alpha,\beta}$ in (1.4), there exists $q_{0}$ satisfying
$2<q_{0}<\infty$ such that for all $q$ satisfying $q_{0}<q\leq\infty$, the
solution map $A\mapsto\chi$ of (1.4) is Lipschitz-continuous as a map from
$L^{q}(\mathbb{T}^{d};\mathbb{R}^{d\times d})$ to
$\dot{H}^{1}(\mathbb{T}^{d};\mathbb{R}^{d})$.

######  Remark 1.4.

Explicit upper bounds for $q_{0}$ in Proposition 1.3 exist and are discussed
in Remark A.14.

Proposition 1.3 is used to prove the approximation theorem for the NO in
Theorems 3.5 and 3.6.

##  2 Microstructures

The main application area of this work is constitutive modeling. In this
section we describe various classes of microstructures that our theory covers.
In particular, we describe four classes of microstructures in two dimensions:

  1. 1.

Smooth microstructures generated via truncated, rescaled log-normal random
fields.

  2. 2.

Discontinuous microstructures with smooth interfaces generated by Lipschitz
star-shaped inclusions.

  3. 3.

Discontinuous microstructures with square inclusions.

  4. 4.

Voronoi crystal microstructures.

Visualizations of examples of these microstructures may be found in Figure 1.
Proofs that the non-smooth classes of microstructures satisfy the assumptions
of the theorems in Section 3 may be found in Appendix C.

Figure 1: Microstructure Examples

#### Smooth Microstructures

The smooth microstructures are generated by exponentiating a rescaled Gaussian
random field. $A$ is symmetric and coercive everywhere in the domain with a
bounded eigenvalue ratio. Furthermore, the smooth function $A$ and its
derivatives are Lipschitz. Our theory is developed specifically to analyze
non-smooth microstructures, so this example is used mainly as a point of
comparison.

#### Star-Shaped Inclusions

For the star-shaped inclusion microstructure, $A$ is taken to be constant
inside and outside the star-shaped boundary. The boundary function is smooth
and Lipschitz in each of its derivatives. $A$ is positive and coercive in both
regions with a bounded eigenvalue ratio. This microstructure introduces
discontinuities, but the boundary remains smooth.

#### Square Inclusions

The square inclusion microstructures are composed of two materials; one
constant $A$ inside the square inclusion, and another constant $A$ outside the
square inclusion. Since we assume periodicity, without loss of generality the
square inclusion is centered. The size of the square inclusion within the cell
is varied between samples as are the constant values of $A$. This
microstructure builds on the complexity of the star inclusion microstructure
by adding corners to the inclusion boundary.

#### Voronoi Interfaces

The Voronoi crystal microstructures are generated by assuming a random Voronoi
tessellation and letting $A$ be piecewise-constant taking a single value on
each Voronoi cell. The number of cells, values of $A$ on the cells, and
locations of the cell centers may all be varied. This is the most complex
microstructure among our examples and is a primary motivation for this work as
Voronoi tessellations are a common model for crystal structure in materials.

##  3 Universal Approximation Results

In this section we state the four approximation theorems for learning solution
operators to the cell problem. Theorems 3.3 and 3.5 concern learning the map
$A\to\chi$ in equation (1.4), and Theorems 3.4 and 3.6 concern learning the
map $A\to\overline{A}$ described by the combination of equations (1.4) and
(1.3). Theorems 3.4 and 3.3 are specific to learning a Fourier neural operator
(FNO), which is a subclass of the general neural operator described by
Theorems 3.6 and 3.5. The proofs of the theorems in this section may be found
in Appendix B.

###  3.1 Definitions of Neural Operators

First, we define a general neural operator (NO) and the Fourier neural
operator (FNO), which is a subclass of the former. The definitions are largely
taken from [31], and we refer to this work for a more in-depth understanding
of these operators. In this work, we restrict the domain to the torus.

######  Definition 3.1 (General Neural Operator).

Let $\mathcal{A}$ and $\mathcal{U}$ be two Banach spaces of real vector-valued
functions over domain $\mathbb{T}^{d}$. Assume input functions
$a\in\mathcal{A}$ are $\mathbb{R}^{d_{a}}$-valued while the output functions
$u\in\mathcal{U}$ are $\mathbb{R}^{d_{u}}$-valued. The neural operator
architecture $\mathcal{G}_{\theta}:\mathcal{A}\to\mathcal{U}$ is

$\displaystyle\mathcal{G}_{\theta}$
$\displaystyle=\mathcal{Q}\circ\mathsf{L}_{T-1}\circ\dots\circ\mathsf{L}_{0}\circ\mathcal{P},$
$\displaystyle v_{t+1}$
$\displaystyle=\mathsf{L}_{t}v_{t}=\sigma_{t}(W_{t}v_{t}+\mathcal{K}_{t}v_{t}+b_{t}),\quad
t=0,1,\dots,T-1$

with $v_{0}=\mathcal{P}(a)$, $u=\mathcal{Q}(v_{T})$ and
$\mathcal{G}_{\theta}(v_{0})=u$. Here,
$\mathcal{P}:\mathbb{R}^{d_{a}}\to\mathbb{R}^{d_{v_{0}}}$ is a local lifting
map, $\mathcal{Q}:\mathbb{R}^{d_{v_{T}}}\to\mathbb{R}^{d_{u}}$ is a local
projection map and the $\sigma_{t}$ are fixed nonlinear activation functions
acting locally as maps $\mathbb{R}^{d_{v_{t+1}}}\to\mathbb{R}^{d_{v_{t+1}}}$
in each layer (with all of $\mathcal{P}$, $\mathcal{Q}$ and the $\sigma_{t}$
viewed as operators acting pointwise, or pointwise almost everywhere, over the
domain $\mathbb{T}^{d}$), $W_{t}\in\mathbb{R}^{d_{v_{t+1}}\times d_{v_{t}}}$
are matrices,
$\mathcal{K}_{t}:\\{v_{t}:\mathbb{T}^{d}\to\mathbb{R}^{d_{v_{t}}}\\}\to\\{v_{t+1}:\mathbb{T}^{d}\to\mathbb{R}^{d_{v_{t+1}}}\\}$
are integral kernel operators and
$b_{t}:\mathbb{T}^{d}\to\mathbb{R}^{d_{v_{t+1}}}$ are bias functions. For any
$m\in\mathbb{N}_{0}$, the activation functions $\sigma_{t}$ are restricted to
the set of continuous $\mathbb{R}\to\mathbb{R}$ maps which make real-valued,
feed-forward neural networks dense in $C^{m}(\mathbb{R}^{d})$ on compact sets
for any fixed network depth. We note that all globally Lipschitz, non-
polynomial, $C^{m}(\mathbb{R})$ functions belong to this class. The integral
kernel operators $\mathcal{K}_{t}$ are defined as

$(\mathcal{K}_{t}v_{t})(x)=\int_{\mathbb{T}^{d}}\kappa_{t}(x,y)v_{t}(y)\>dy$

with standard multi-layered perceptrons (MLP)
$\kappa_{t}:\mathbb{T}^{d}\times\mathbb{T}^{d}\to\mathbb{R}^{d_{v_{t+1}}\times
d_{v_{t}}}$. We denote by $\theta$ the collection of parameters that specify
$\mathcal{G}_{\theta}$, which include the weights $W_{t}$, biases $b_{t}$,
parameters of the kernels $\kappa_{t}$, and the parameters describing the
lifting and projection maps $\mathcal{P}$ and $\mathcal{Q}$ (usually also
MLPs).

The FNO is a subclass of the NO.

######  Definition 3.2 (Fourier Neural Operator).

The FNO inherits the structure and definition of the NO in Definition 3.1,
together with some specific design choices. We fix $d_{v_{t}}=d_{v}$ for all
$t$, where $d_{v}$ is referred to as the number of channels, or model width,
of the FNO. We fix $\sigma_{t}=\sigma$ to be a globally Lipschitz, non-
polynomial, $C^{\infty}$ function.111In this work in all numerical experiments
we use the GeLU activation function as in [34]. Finally, the kernel operators
$\mathcal{K}_{t}$ are parameterized in the Fourier domain in the following
manner. Let

$\psi_{k}(x)=\text{e}^{2\pi i\langle k,x\rangle},\quad
x\in\mathbb{T}^{d},\;k\in\mathbb{Z}^{d},$

denote the Fourier basis for $L^{2}(\mathbb{T}^{d};\mathbb{C})$ where
$i=\sqrt{-1}$ is the imaginary unit. Then, for each $t$, kernel operator
$\mathcal{K}_{t}$ is parameterized by

$(\mathcal{K}_{t}v_{t})_{l}(x)=\sum_{\begin{subarray}{c}k\in\mathbb{Z}^{d}\\\
|k|\leq
k_{\max}\end{subarray}}\left(\sum_{j=1}^{d_{v}}P_{lj}^{k}\langle(v_{t})_{j},\psi_{k}\rangle_{L^{2}(\mathbb{T}^{d};\mathbb{C})}\right)\psi_{k}(x).$

Here, $l=1,\dots,d_{v_{l+1}}$ and each $P^{k}\in\mathbb{C}^{d_{v_{l+1}}\times
d_{v_{l}}}$ constitute the learnable parameters of the integral operator.

From the definition for the FNO, we note that parameterizing the kernels in
the Fourier domain allows for efficient computation using the FFT. We refer to
[31, 34] for additional details.

###  3.2 Main Theorems

The first two theorems guarantee the existence of an FNO approximating the
maps $A\to\chi$ and $A\to\overline{A}$ and are based on the stability estimate
for continuity from $L^{2}\to\dot{H}^{1}$ obtained in Proposition 1.2.

######  Theorem 3.3.

Let $K\subset\mathsf{PD}_{\alpha,\beta}$ and define the mapping
$G:K\to\dot{H}^{1}(\mathbb{T}^{d};\mathbb{R}^{d})$ by $A\mapsto\chi$ as given
by (1.4). Then, for any $\epsilon>0$ and $K$ compact in
$L^{2}(\mathbb{T}^{d};\mathbb{R}^{d\times d})$, there exists a FNO
$\Psi:K\to\dot{H}^{1}(\mathbb{T}^{d};\mathbb{R}^{d})$ such that

$\sup_{A\in K}\|G(A)-\Psi(A)\|_{\dot{H}^{1}}<\epsilon.$

######  Theorem 3.4.

Let $K\subset\mathsf{PD}_{\alpha,\beta}$ and define the mapping
$F:K\to\mathbb{R}^{d\times d}$ by $A\mapsto\bar{A}$ as given by (1.3). Then,
for any $\epsilon>0$ and $K$ compact in
$L^{2}(\mathbb{T}^{d};\mathbb{R}^{d\times d})$ there exists an FNO
$\Phi:K\to\mathbb{R}^{d\times d}$ such that

$\sup_{A\in K}|F(A)-\Phi(A)|_{F}<\epsilon.$

The remaining two theorems guarantee the existence of a general neural
operator approximating the maps $A\to\chi$ and $A\to\overline{A}$. Although
the FNO is a subclass of general neural operators, from a theoretical
perspective, the stability estimates used to obtain the results for general
neural operators are more concrete. Indeed, the following two theorems are
based on the stability estimate of Proposition 1.3, which obtains Lipschitz
continuity from $L^{q}\to\dot{H}^{1}$.

######  Theorem 3.5.

Let $K\subset\mathsf{PD}_{\alpha,\beta}$ and define the mapping
$G:K\to\dot{H}^{1}(\mathbb{T}^{d};\mathbb{R}^{d})$ by $A\mapsto\chi$ as given
by (1.4). Let $q_{0}$ be as in Proposition 1.3. Then, for any $q$ satisfying
$q_{0}<q<\infty$ and for any $K$ compact in
$L^{q}(\mathbb{T}^{d};\mathbb{R}^{d\times d})$, it holds that for any
$\epsilon>0$, there exists a neural operator
$\Psi:K\to\dot{H}^{1}(\mathbb{T}^{d};\mathbb{R}^{d})$ such that

$\sup_{A\in K}\|G(A)-\Psi(A)\|_{\dot{H}^{1}}<\epsilon.$

######  Theorem 3.6.

Let $K\subset\mathsf{PD}_{\alpha,\beta}$ and define the mapping
$F:K\to\mathbb{R}^{d\times d}$ by $A\mapsto\bar{A}$ as given by (1.3). Let
$q_{0}$ be as in Proposition 1.3. Then, for any $q$ satisfying
$q_{0}<q<\infty$ and for any $K$ compact in
$L^{q}(\mathbb{T}^{d};\mathbb{R}^{d\times d})$, it holds that for any
$\epsilon>0$, there exists a neural operator $\Phi:K\to\mathbb{R}^{d\times d}$
such that

$\sup_{A\in K}|F(A)-\Phi(A)|_{F}<\epsilon.$

######  Remark 3.7.

Although the statements of Theorems 3.5 and 3.6 seem almost identical to those
of Theorems 3.3 and 3.4, the proofs are quite different as they rely on
stability estimates obtained through entirely different methods. Additionally,
they depend on different universal approximation theorems. The FNO results use
the universal approximation theorem for FNOs in [30], while the NO results use
the universal approximation theorem for NOs in [31].

The above approximation results can also be formulated to hold, on average,
over any probability measure with a finite second moment that is supported on
$\mathsf{PD}_{\alpha,\beta}$. In particular, if we let $\mu$ be such a
probability measure then there exists an FNO or a neural operator $\Psi$ such
that

$\mathbb{E}_{A\sim\mu}\|G(A)-\Psi(A)\|_{\dot{H}^{1}}<\epsilon.$

This follows by simple exchanging the appropriate results from [30] or [31] in
the respective proofs. We do not carry out the full details here. While this
allows approximation over the non-compact set $\mathsf{PD}_{\alpha,\beta}$,
the error can only be controlled on average instead of uniformly. Such results
might find applications in situations where the microstructures are modeled
stochastically.

##  4 Numerical Experiments

In this section we perform numerical experiments to illustrate the fact that
it is possible to find good operator approximations of the homogenization maps
(1.6), (1.7) in practice. We focus on use of the FNO and note that, whilst
Theorems 3.3 and 3.4 assert the existence of desirable operator
approximations, they are not constructive and they do not come equipped with
error estimates. We find approximations using standard empirical loss
minimization techniques and quantify the complexity with respect to data and
parametric size of the approximations.

To verify that Theorems 3.3 and 3.4 apply, we have to show that the subset of
coefficient functions employed are compact. Lemma C.1 applies to show
compactness for the square inclusion and Voronoi crystal interfaces of
Examples 3 and 4. Lemma C.2 applies to show compactness for the star-shaped
inclusions of Example 2. The smooth microstructure example serves as a
comparison case for examining the impact of discontinuous coefficients on the
learning accuracy.

The experiments are all conducted using an FNO with a fixed number $T=4$ of
hidden layers. The two remaining parameters to vary are the channel width
$d_{v}$ and the number of Fourier modes $k_{max}$. We make the following
observations based on the numerical experiments:

  1. 1.

The effective $\overline{A}$ tensors computed from the model predicted
solutions exhibit very low error.

  2. 2.

The error in the learned $\chi$ is significantly higher along discontinuous
material boundaries and corner interfaces, as expected. However, the FNO
operator approximation is able to approximate the solution with reasonable
relative error even for the most complex case of a set of input functions with
varying Voronoi geometry and varying microstructural properties within
domains.

  3. 3.

In comparison with the smooth microstructure case, learning the map for the
Voronoi microstructure requires substantially more data to avoid training a
model which plateaus at a poor level of accuracy.

  4. 4.

When compared with the smooth microstructure case, error for the Voronoi
microstructure decreases less rapidly with respect to increasing model width,
but shows more favourable response with respect to increasing the number of
Fourier modes.

We first describe implementation details of each of the microstructures in
Subsection 4.1. Then we show outcomes of the numerical experiments in
Subsection 4.2, which we discuss in Subsection 4.3.

###  4.1 Microstructure Implementation

#### Smooth Microstructures

The smooth microstructures are generated by exponentiating a rescaled
approximation of a Gaussian random field. We consider material tensor $A(x)$
given by

$A(x)=\begin{bmatrix}\lambda_{1}&0\\\ 0&\lambda_{2}\end{bmatrix}.$ (4.1)

The random field used to generate the diagonal entries of $A(x)$ is
constructed as follows:

$\displaystyle\widehat{\lambda}_{i}(x)$
$\displaystyle=\sum_{k_{1},k_{2}=1}^{4}\xi^{(1)}_{k_{1},k_{2}}\sin\left(2\pi
k_{1}x_{1}\right)\cos(2\pi k_{2}x_{2})+\xi^{(2)}_{k_{1},k_{2}}\cos\left(2\pi
k_{1}x_{1}\right)\sin\left(2\pi k_{2}x_{2}\right),$ (4.2)
$\displaystyle\lambda_{i}(x)$
$\displaystyle=\exp\left(\frac{\widehat{\lambda}_{i}(x)}{\max_{x^{\prime}\in[0,1]^{2}}|\widehat{\lambda}_{i}(x^{\prime})|}\right),$
(4.3)

where $\xi^{(j)}_{k_{1},k_{2}}$ are i.i.d. normal Gaussian random variables.
Due to the rescaling and exponentiation, the ratio of the eigenvalues of
$A(x)$ is at most $e^{2}$.

#### Star-Shaped Inclusions

The star-shaped inclusions are generated by defining a random Lipschitz polar
boundary function as

$\displaystyle r(\theta)$
$\displaystyle=\alpha+\beta\sum_{k=1}^{5}\xi_{k}\sin(k\theta)$ (4.4)

where $\xi_{k}$ are i.i.d. normal random variables, and $\alpha$ and $\beta$
are constants. Then $A(x)$ is constant inside and outside the boundary. We
randomly sample eigenvalues for $A$ on each domain via
$\lambda_{i}=0.1+\zeta_{i}$ where $\zeta_{i}$ are uniform random variables on
$[0,1]$. Three components of the eigenvectors are i.i.d. normal random
variables, and the fourth component enforces orthogonality to guarantee
symmetry of $A$. In this manner, $A$ is symmetric and coercive and has a
bounded eigenvalue ratio.

#### Square Inclusions

The radius of the square is randomly generated via

$r=\alpha+\beta\zeta$

where $\zeta$ is a uniform random variable on $[0,1]$ and $\alpha$ and $\beta$
are positive constants. The values of $A$ on each of the constant domains are
chosen in the same manner as in the star-shaped inclusion case.

#### Voronoi Interfaces

The Voronoi crystal microstructure has $A=aI$ where $a$ is constant on each
Voronoi cell and is chosen uniformly at random in the same manner as for the
star inclusions. Voronoi tessellations are a common model for crystal
structure in materials. In one Voronoi example, we we fix the geometry for all
data, and in a second Voronoi example we vary the geometry by randomly
sampling five cell centers from a uniform distribution on the unit square.

###  4.2 Results

Since the cell problem is separable into two problems for each component of
$\chi$, we learn $\chi_{1}$ and $\chi_{2}$ separately; in this paragraph,
therefore, $j\in\\{1,2\\}.$ Each FNO model was trained using the empirical
estimate of the mean squared $H^{1}$ norm:

$\text{Loss}(\theta)=\frac{1}{N}\sum_{n=1}^{N}\Bigl{(}\|\chi^{(n)}_{j}-\widehat{\chi}^{(n)}_{j}\|^{2}_{L^{2}}+\|\nabla\chi_{j}^{(n)}-\nabla\widehat{\chi}^{(n)}_{j}\|^{2}_{L^{2}}\Bigr{)}$
(4.5)

where $n$ is the sample index, $\chi_{j}$ is the true solution, and
$\widehat{\chi}_{j}$ is the FNO approximation of the solution, parameterized
by $\theta.$ In the analysis, we examine several different measures of error,
including the following relative $H^{1}$ and relative $W^{1,10}$ errors.

$\text{Relative $H^{1}$ Error
(RHE)}=\frac{1}{N}\sum_{n=1}^{N}\left(\frac{\|\chi_{j}^{(n)}-\widehat{\chi}_{j}^{(n)}\|^{2}_{L^{2}}+\|\nabla\chi_{j}^{(n)}-\nabla\widehat{\chi}_{j}^{(n)}\|^{2}_{L^{2}}}{\|\chi_{j}^{(n)}\|^{2}_{L^{2}}+\|\nabla\chi_{j}^{(n)}\|^{2}_{L^{2}}}\right)^{1/2}$
(4.6a) $\text{Relative $W^{1,10}$ Error
(RWE)}=\frac{1}{N}\sum_{n=1}^{N}\left(\frac{\|\chi_{j}^{(n)}-\widehat{\chi}_{j}^{(n)}\|^{10}_{L^{10}}+\|\nabla\chi_{j}^{(n)}-\nabla\widehat{\chi}_{j}^{(n)}\|^{10}_{L^{10}}}{\|\chi_{j}^{(n)}\|^{10}_{L^{10}}+\|\nabla\chi_{j}^{(n)}\|^{10}_{L^{10}}}\right)^{1/10}$
(4.6b)

The error visualized in Figures 2 and 3 is spatially-varying and computed via

$\text{Scaled $L^{2}$ Error}(x)=\frac{1}{\|g\|_{L^{2}}}|g(x)-\widehat{g}(x)|$
(4.7)

where $g$ is $\chi^{(n)}_{j}$ and $\nabla\chi^{(n)}_{j}$ for a particular
sample $n$.

Finally, we also look at error in $\overline{A}$, which we scale by the
difference between the arithmetic and harmonic mean of $A$. Any effective
$\overline{A}$ should fall between the arithmetic and harmonic means of $A$.
Thus, this error is given by

$\text{Relative $\overline{A}$ Error (RAE)
}=\frac{1}{N}\sum_{n=1}^{N}\frac{\|\overline{A}^{(n)}-\widehat{\overline{A}}^{(n)}\|_{F}}{a^{(n)}_{m}-a^{(n)}_{h}}$
(4.8)

where the arithmetic mean $a_{m}$ and harmonic mean $a_{h}$ are given by

$\displaystyle a_{m}$
$\displaystyle=\left\|\int_{\mathbb{T}^{2}}A(x)\>\mathsf{d}x\right\|_{F}$
(4.9) $\displaystyle a_{h}$
$\displaystyle=\left\|\left(\int_{\mathbb{T}^{2}}A^{-1}(x)\>\mathsf{d}x\right)^{-1}\right\|_{F}.$
(4.10)

We trained six different models whose details may be viewed in Table 1. Each
of these models was trained on $9500$ data samples generated using a FEM
solver on a triangular mesh with the solution interpolated to a $128\times
128$ grid. The models were each trained with an FNO architecture for $300$
epochs with the same architecture of $24$ Fourier modes, a model width of
$48$, and $4$ layers. Table 1 also includes a quantification of the error.
Visualizations of the median-error test samples for each example may be viewed
in Figures 2 and 3. These figures are similar in form, and may be compared to,
those in [18].

We also investigated the effects of the number of training data and the model
size on the error for the smooth and Voronoi microstructures described as
Examples 1 and 5 in Table 1. A plot of error versus training data may be found
in Figure 5, and plots of error versus the number of Fourier modes for fixed
total model size, as measured by (model width) $\times$ (number Fourier
modes), may be found in Figure 5. Figure 5 is similar to, and may be examined
along with, those in [18]. Figure 5 addresses the question of how to optimally
distribute computational budget through different parameterizations to achieve
minimum error at given cost as measured by number of parameters; it should be
compared to similar experiments in [32].

Ex. | Structure | Form of $A$ | Mean RHE | Med. RHE | Mean RWE | Med. RWE | Med. RAE
---|---|---|---|---|---|---|---
1 | Smooth | $\begin{bmatrix}a_{11}&0\\\ 0&a_{22}\end{bmatrix}$ | 0.0053 | 0.0051 | 0.0074 | 0.0067 | 0.00058
2 | Star | $\begin{bmatrix}a_{11}&a_{12}\\\ a_{12}&a_{22}\end{bmatrix}$ | 0.0397 | 0.0350 | 0.1490 | 0.1415 | 0.00228
3 | Square | $\begin{bmatrix}a_{11}&a_{12}\\\ a_{12}&a_{22}\end{bmatrix}$ | 0.0940 | 0.0806 | 0.2294 | 0.2093 | 0.00893
4 | Voronoi | $aI$ | 0.0045 | 0.0041 | 0.0072 | 0.0067 | 0.00027
5 | Voronoi | $aI$ | 0.0557 | 0.0504 | 0.1987 | 0.1828 | 0.00151
6 | Voronoi | $aI$ | 0.0501 | 0.0450 | 0.1864 | 0.1686 | 0.00155
Table 1: Errors and key information about each numerical example. The
expressions for the RHE (Relative $H^{1}$ Error), RWE (Relative $W^{1,10}$
Error) and RAE (Relative $\overline{A}$ Error) may be found in equations (4.6)
and (4.8). For each error measurement, the mean and median error over the
$500$ test samples is shown as well as the median error in the effective
tensor $\overline{A}$. The data for Examples 1-5 were generated using an FEM
mesh of 90804 elements, and the data for Example 6 were generated using a mesh
of 551062 elements. All solutions were then interpolated to a uniform
$128\times 128$ grid. Voronoi Example $4$ had fixed geometry in the
microstructure data whereas the Examples $5,6$ had varying geometry in the
microstructure data.

Smooth microstructure (Example 1)

TS

Star inclusion microstructure (Example 2)

Square inclusion microstructure (Example 3)

Voronoi microstructure (Example 4)

TS

Figure 2: Visualization of the median error samples for the FNO Approximations
of Examples 1-4 in Table 1. For each example, Subplot (a) is the true solution
component $\chi_{1}$ for the median-error sample. Subplot (b) is the scaled
absolute error in $\chi_{1}$ given by (4.7). Subplot (c) is the first
component of the true solution gradient $\nabla\chi_{1}$. Finally, Subplot (d)
is the scaled absolute error in the first component of the gradient
$\nabla\chi_{1}$.

TS

Figure 3: FNO approximation of $\chi_{1}$ and $\nabla\chi_{1}$ for varying-
geometry Voronoi microstructure for two different levels of data accuracy
corresponding to Examples 5 and 6 in Table 1.

TS

###  4.3 Discussion

As can be seen from the data in Table 1, the microstructures exhibiting
discontinuities lead to higher model error than the smooth microstructure. The
visualizations of the median-error test samples in Figure 2 give some
intuition; error is an order of magnitude higher along discontinuous
boundaries; this is most apparent in the gradient. The true solution gradient
often takes its most extreme values along the discontinuities, and the RWE
gives an indication of how well the model captures the most extreme values in
the solution. Unsurprisingly, this error is much higher than the RHE, but we
note that it is confined to a small area of the domain along discontinuous
boundaries and corner interfaces.

Voronoi Examples 5 and 6 in Table 1 are identical except that the numerical
solution for the two data sets is obtained using solvers with different levels
of accuracy. For Example $5$, the accuracy is the same as for Examples 1-4,
but for 6 the solver uses $551062$ elements instead of $90804$. Both solutions
are then interpolated to a $128\times 128$ grid to train the model. A common
belief among machine learning practitioners is that the error from the model
approximation is independent from the error of the numerical solution used to
obtain the data. Since the model only “sees” the data, the reasoning goes that
model will obtain similar train and test error relative to the data regardless
of the numerical accuracy of the data relative to the true solution. With this
experiment, we see that this does not hold for the Voronoi microstructure;
using a more accurate solver results in lower model error in both the RHE and
RWE measurements. These observations lead to a hypothesis that the errors
produced via a numerical solver are harder to learn than the true solution
behavior, but we leave the formalization and validation of this hypothesis to
future work.

We also examine the effect of the number of training data samples and the FNO
size on model accuracy for the smooth and Voronoi microstructures
corresponding to Examples 1 and 4 in Table 1. In Figure 5 we see the effect of
increasing the number of training data. The Voronoi microstructure example
benefits greatly from additional data, but the smooth microstructure example
saturates. We note that this is in contrast to the behavior during training
over $300$ epochs; the test error for the smooth microstructure continues to
decrease over the entire training periodic, but the test error for the Voronoi
microstructure plateaus by around 100-150 epochs. The model size also presents
a qualitatively different effect on error for the smooth and Voronoi
microstructures. In Figure 5, we see the tradeoff between the number of
Fourier modes and the model width for approximately constant model size,
measured as the product of the width and number of modes. The Voronoi example
benefits from additional Fourier modes, whereas the smooth example worsens. On
the other hand, the smooth model benefits more from an increase in model
width. We refer to [18, 32] for in-depth numerical studies of errors, choice
of hyperparameters, and parameter distributions for FNO; here we highlight
only the qualitative differences between the model behavior for different
microstructures.

Finally, we compare the error in the effective $\overline{A}$ defined in
(1.3). This error is scaled by a difference between the Frobenius norms of the
arithmetic and harmonic means of the true $A$ because the Frobenius norm of
the true $\overline{A}$ should fall within that range. For this reason, in the
case where the arithmetic and harmonic means are very close, it is not
valuable to learn the true $\overline{A}$. The varying-geometry Voronoi
microstructure examples on average have about 100 times greater difference
between the means than the star and square microstructure examples. Thus, it
is valuable to note that the median relative $\overline{A}$ error shown in
Table 1 is lower for the Voronoi examples than for the star and square
inclusion examples because in the Voronoi setting, the arithmetic and harmonic
means are poor approximations of $\overline{A}$. This characteristic of the
Voronoi microstructure further underscores the value of learning in this
setting.

Figure 4: Relative $H^{1}$ error versus the number of training data used for
the smooth and Voronoi Examples 1 and 5. Model size was fixed at 12 modes and
32 width.

Figure 5: Relative $H^{1}$ error versus model size for the smooth and Voronoi
Examples 1 and 5. The number of Fourier modes in each direction and the model
width were varied. Each line indicates a constant product of
modes$\times$width. Training data size was fixed at 9500 samples.

##  5 Conclusions

In this work, we establish approximation theory for learning the solution
operator arising from the elliptic homogenization cell problem (1.4) viewed as
a mapping from the coefficient to the solution; the theory allows for
discontinuous coefficients. We also perform numerical experiments that
validate the theory, explore qualitative differences between various
microstructures, and quantify error/cost trade-offs in the approximation. We
provide two different proof approaches that rely on different stability
results for the underlying solutions. These stability results, when combined
with existing universal approximation results for neural operators, result in
rigorous approximation theory for learning in this problem setting. On the
numerical side, we provide examples of various microstructures that satisfy
the conditions of the approximation theory. We observe that model error is
dominated by error along discontinuous and corner interfaces, and that
discontinuous microstructures give rise to qualitatively different learning
behavior. Finally, we remark that the learned effective properties are highly
accurate, especially in the case of the Voronoi microstructure that we regard
as the most complex. Since discontinuous microstructures arise naturally in
solid mechanics, understanding learning behavior in this context is an
important prerequisite for using machine learning for applications. In this
area and others, numerous questions remain which address the rigor necessary
for use of machine learning in scientific applications.

We have confined our studies to one of the canonical model problems of
homogenization theory, the divergence-form elliptic setting, with periodic
microstructure, to obtain deeper understanding of the learning constitutive
laws. One interesting potential extension of this work is the setting in which
the material coefficient $A$ is not periodic but random with respect to the
microstructure. Another is where it is only locally periodic and has
dependence on the macroscale variable as well; thus
$A^{\epsilon}=A(x,\frac{x}{\epsilon})$. In this case, the form of the cell
problem (1.4) and homogenized coefficient (1.3) remain the same, but $A$ and
$\chi$ both have parametric dependence on $x$. The approximation theory, and
the empirical learning problem, would grow in complexity in comparison to what
is developed here, but the resulting methodology could be useful and
foundational for understanding more complex constitutive models in which the
force balance equation couples to other variables. Indeed, the need for
efficient learning of constitutive models is particularly pressing in complex
settings such as crystal plasticity. We anticipate that the potential use of
machine learning, to determine parametric dependence of constitutive models
defined by homogenization, will be for these more complex problems. The work
described in this paper provides an underpinning conceptual approach,
foundational analysis and set of numerical experiments that serve to underpin
more applied work in this field.

#### Acknowledgements

The work of KB, NK, AR and AMS was sponsored by the Army Research Laboratory
and was accomplished under Cooperative Agreement W911NF-12-2-0022. AMS is also
grateful to the DoD for support through a Vannevar Bush Faculty Fellowship.
The work of MT was funded by the Department of Energy Computational Science
Graduate Fellowship under award DE-SC002111.

## References

  * [1] R. A. Adams and J. J. Fournier, Sobolev Spaces, Elsevier, 2003.
  * [2] G. Allaire, Homogenization and two-scale convergence, SIAM Journal on Mathematical Analysis, 23 (1992), pp. 1482–1518.
  * [3] F. As’ ad, P. Avery, and C. Farhat, A mechanics-informed artificial neural network approach in data-driven constitutive modeling, International Journal for Numerical Methods in Engineering, 123 (2022), pp. 2738–2759.
  * [4] P. Batlle, M. Darcy, B. Hosseini, and H. Owhadi, Kernel methods are competitive for operator learning, arXiv preprint arXiv:2304.13202, (2023).
  * [5] A. Bensoussan, J.-L. Lions, and G. Papanicolaou, Asymptotic analysis for periodic structures, vol. 374, American Mathematical Soc., 2011.
  * [6] J. Bergh and J. Löfström, Interpolation Spaces: An Introduction, vol. 223, Springer Science & Business Media, 2012.
  * [7] K. Bhattacharya, Phase boundary propagation in a heterogeneous body, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 455 (1999), pp. 757–766.
  * [8] K. Bhattacharya, B. Hosseini, N. B. Kovachki, and A. M. Stuart, Model reduction and neural networks for parametric pdes, The SMAI journal of computational mathematics, 7 (2021), pp. 121–157.
  * [9] K. Bhattacharya, B. Liu, A. Stuart, and M. Trautner, Learning Markovian homogenized models in viscoelasticity, Multiscale Modeling & Simulation, 21 (2023), pp. 641–679.
  * [10] X. Blanc and C. Le Bris, Homogenization Theory for Multiscale Problems: An Introduction, vol. 21, Springer Nature, 2023.
  * [11] A. Bonito, R. A. DeVore, and R. H. Nochetto, Adaptive finite element methods for elliptic problems with discontinuous coefficients, SIAM Journal on Numerical Analysis, 51 (2013), pp. 3106–3134.
  * [12] T. Brown, B. Mann, N. Ryder, M. Subbiah, J. D. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell, et al., Language models are few-shot learners, Advances in neural information processing systems, 33 (2020), pp. 1877–1901.
  * [13] A. Chkifa, A. Cohen, R. DeVore, and C. Schwab, Sparse adaptive taylor approximation algorithms for parametric and stochastic elliptic pdes, ESAIM: Mathematical Modelling and Numerical Analysis, 47 (2013), pp. 253–280.
  * [14] D. Cioranescu and P. Donato, An Introduction To Homogenization, vol. 17, Oxford university press Oxford, 1999.
  * [15] A. Cohen, R. DeVore, and C. Schwab, Convergence rates of best n-term galerkin approximations for a class of elliptic spdes, Foundations of Computational Mathematics, 10 (2010), pp. 615–646.
  * [16] A. Cohen, R. Devore, and C. Schwab, Analytic regularity and polynomial approximation of parametric and stochastic elliptic pde’s, Analysis and Applications, 9 (2011), pp. 11–47.
  * [17] G. Cybenko, Approximation by superpositions of a sigmoidal function, Mathematics of Control, Signals and Systems, 2 (1989), pp. 303–314.
  * [18] M. V. de Hoop, Z. Huang, E. Qian, and A. M. Stuart, The cost-accuracy trade-off in operator learning with neural networks, Journal of Machine Learning, 1 (2022), pp. 299–341.
  * [19] R. DeVore and G. Lorentz, Constructive Approximation, Grundlehren der mathematischen Wissenschaften, Springer Berlin Heidelberg, 1993.
  * [20] J. Dugundji, An extension of Tietze’s theorem., Pacific Journal of Mathematics, 1 (1951), pp. 353–367.
  * [21] J. N. Fuhg and N. Bouklas, On physics-informed data-driven isotropic and anisotropic constitutive models through probabilistic machine learning and space-filling sampling, Computer Methods in Applied Mechanics and Engineering, 394 (2022).
  * [22] J.-L. Guermond, The lbb condition in fractional sobolev spaces and applications, IMA journal of numerical analysis, 29 (2009), pp. 790–805.
  * [23] M. E. Gurtin, An Introduction To Continuum Mechanics, Academic press, 1982.
  * [24] E. Haghighat, S. Abouali, and R. Vaziri, Constitutive model characterization and discovery using physics-informed deep learning, Engineering Applications of Artificial Intelligence, 120 (2023), p. 105828.
  * [25] J. Han and Y. Lee, A neural network approach for homogenization of multiscale problems, Multiscale Modeling & Simulation, 21 (2023), pp. 716–734.
  * [26] K. He, X. Zhang, S. Ren, and J. Sun, Deep residual learning for image recognition, in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 770–778.
  * [27] L. Herrmann, C. Schwab, and J. Zech, Neural and gpc operator surrogates: construction and expression rate bounds, arXiv preprint arXiv:2207.04950, (2022).
  * [28] T. Y. Hou and X.-H. Wu, A multiscale finite element method for elliptic problems in composite materials and porous media, Journal of computational physics, 134 (1997), pp. 169–189.
  * [29] D. Z. Huang, K. Xu, C. Farhat, and E. Darve, Learning constitutive relations from indirect observations using deep neural networks, Journal of Computational Physics, 416 (2020), p. 109491.
  * [30] N. Kovachki, S. Lanthaler, and S. Mishra, On universal approximation and error bounds for Fourier Neural Operators, The Journal of Machine Learning Research, 22 (2021), pp. 1–76.
  * [31] N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, and A. Anandkumar, Neural operator: Learning maps between function spaces with applications to pdes, Journal of Machine Learning Research, 24 (2023), pp. 1–97.
  * [32] S. Lanthaler, Z. Li, and A. M. Stuart, The nonlocal neural operator: universal approximation, arXiv preprint arXiv:2304.13221, (2023).
  * [33] S. Lanthaler, S. Mishra, and G. E. Karniadakis, Error estimates for deeponets: A deep learning framework in infinite dimensions, Transactions of Mathematics and Its Applications, 6 (2022).
  * [34] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, Fourier neural operator for parametric partial differential equations, International Conference on Learning Representations, (2021).
  * [35] K. Linka, M. Hillgärtner, K. P. Abdolazizi, R. C. Aydin, M. Itskov, and C. J. Cyron, Constitutive artificial neural networks: A fast and general approach to predictive data-driven constitutive modeling by deep learning, Journal of Computational Physics, 429 (2021), p. 110010.
  * [36] B. Liu, N. Kovachki, Z. Li, K. Azizzadenesheli, A. Anandkumar, A. M. Stuart, and K. Bhattacharya, A learning-based multiscale method and its application to inelastic impact problems, Journal of the Mechanics and Physics of Solids, 158 (2022).
  * [37] B. Liu, E. Ocegueda, M. Trautner, A. M. Stuart, and K. Bhattacharya, Learning macroscopic internal variables and history dependence from microscopic models, Journal of the Mechanics and Physics of Solids, (2023).
  * [38] X. Liu, F. Tao, and W. Yu, A neural network enhanced system for learning nonlinear constitutive law and failure initiation criterion of composites using indirectly measurable data, Composite Structures, 252 (2020), p. 112658.
  * [39] X. Liu, S. Tian, F. Tao, and W. Yu, A review of artificial neural networks in the constitutive modeling of composite materials, Composites Part B: Engineering, 224 (2021), p. 109152.
  * [40] H. J. Logarzo, G. Capuano, and J. J. Rimoli, Smart constitutive laws: Inelastic homogenization through machine learning, Computer methods in applied mechanics and engineering, 373 (2021), p. 113482.
  * [41] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis, Learning nonlinear operators via deeponet based on the universal approximation theorem of operators, Nature machine intelligence, 3 (2021), pp. 218–229.
  * [42] P. Ma, P. Y. Chen, B. Deng, J. B. Tenenbaum, T. Du, C. Gan, and W. Matusik, Learning neural constitutive laws from motion observations for generalizable pde dynamics, arXiv preprint arXiv:2304.14369, (2023).
  * [43] C. Marcati, J. A. Opschoor, P. C. Petersen, and C. Schwab, Exponential relu neural network approximation rates for point and edge singularities, Foundations of Computational Mathematics, (2022), pp. 1–85.
  * [44] M. Mozaffar, R. Bostanabad, W. Chen, K. Ehmann, J. Cao, and M. Bessa, Deep learning predicts path-dependent plasticity, Proceedings of the National Academy of Sciences, 116 (2019), pp. 26414–26420.
  * [45] N. H. Nelsen and A. M. Stuart, The random feature model for input-output maps between banach spaces, SIAM Journal on Scientific Computing, 43 (2021), pp. A3212–A3243.
  * [46] R. H. Nochetto, K. G. Siebert, and A. Veeser, Theory of adaptive finite element methods: an introduction, in Multiscale, Nonlinear and Adaptive Approximation: Dedicated to Wolfgang Dahmen on the Occasion of his 60th Birthday, Springer, 2009, pp. 409–542.
  * [47] J. A. Opschoor, C. Schwab, and J. Zech, Exponential relu dnn expression of holomorphic maps in high dimension, Constructive Approximation, 55 (2022), pp. 537–582.
  * [48] H. Owhadi and L. Zhang, Numerical homogenization of the acoustic wave equations with a continuum of scales, Computer Methods in Applied Mechanics and Engineering, 198 (2008), pp. 397–406.
  * [49] G. Pavliotis and A. Stuart, Multiscale Methods: Averaging and Homogenization, Springer Science & Business Media, 2008.
  * [50] R. Phillips, Crystals, defects and microstructures: modeling across scales, Cambridge University Press, 2001.
  * [51] H.-J. Schmeisser and H. Triebel, Topics in Fourier Analysis and Function Spaces, vol. 1, Wiley, 1987.
  * [52] E. M. Stein and G. Weiss, Introduction to Fourier analysis on Euclidean spaces, vol. 1, Princeton University Press, 1971.
  * [53] A. M. Tartakovsky, C. O. Marrero, P. Perdikaris, G. D. Tartakovsky, and D. Barajas-Solano, Physics-informed deep neural networks for learning parameters and constitutive relationships in subsurface flow problems, Water Resources Research, 56 (2020), p. e2019WR026731.
  * [54] K. Xu, D. Z. Huang, and E. Darve, Learning constitutive relations using symmetric positive definite neural networks, Journal of Computational Physics, 428 (2021), p. 110072.
  * [55] H. You, Y. Yu, S. Silling, and M. D’Elia, Data-driven learning of nonlocal models: from high-fidelity simulations to constitutive laws, CEUR-WS, 2964 (2021).

## Appendices

##  Appendix A Proofs of Stability Estimates

In this section, we prove the stability estimates stated in Section 1.3. The
following lemma is a modification of the standard estimate for parametric
dependence of elliptic equations on their coefficient. We include it here for
completeness.

See 1.1

###### Proof.

For existence and uniqueness of the solution to the cell problem using Lax-
Milgram, we refer to the texts [10, 49]; we simply derive the bounds and
stability estimate. First, note that (1.4) decouples, in particular,

$-\nabla\cdot(A\nabla\chi_{\ell})=\nabla\cdot Ae_{\ell},\qquad
y\in\mathbb{T}^{d}$ (A.1)

for $l=1,\dots,d$ where $e_{\ell}$ is the $\ell$-th standard basis vector of
$\mathbb{R}^{d}$ and each $\chi_{\ell}\in\dot{H}^{1}(\mathbb{T}^{1})$.
Multiplying by $\chi_{\ell}$ and integrating by parts shows

$\displaystyle\alpha\|\nabla\chi_{\ell}\|_{L^{2}}^{2}$
$\displaystyle\leq\int_{\mathbb{T}^{d}}\langle
A\nabla\chi_{\ell},\nabla\chi_{\ell}\rangle\>\mathsf{d}y$
$\displaystyle=-\int_{\mathbb{T}^{d}}\langle
Ae_{\ell},\nabla\chi_{\ell}\rangle\>\mathsf{d}y$
$\displaystyle\leq\int_{\mathbb{T}^{d}}|Ae_{\ell}||\nabla\chi_{\ell}|\>\mathsf{d}y$
$\displaystyle\leq\left(\int_{\mathbb{T}^{d}}|Ae_{\ell}|^{2}\>\mathsf{d}y\right)^{\frac{1}{2}}\left(\int_{\mathbb{T}^{d}}|\nabla\chi_{\ell}|^{2}\>\mathsf{d}y\right)^{\frac{1}{2}}$
$\displaystyle\leq\|A\|_{L^{\infty}}\|\nabla\chi_{\ell}\|_{L^{2}}.$

Therefore

$\|\nabla\chi\|_{L^{2}}^{2}=\sum_{l=1}^{d}\|\nabla\chi_{\ell}\|_{L^{2}}^{2}\leq\frac{d\|A\|_{L^{\infty}}^{2}}{\alpha^{2}}\leq\frac{d\beta^{2}}{\alpha^{2}},$

which implies the first result.

To prove the second result, we denote the right hand side of A.1 by
$f^{(i)}_{\ell}=\nabla\cdot A^{(i)}e_{\ell}$ in what follows. For any
$v\in\dot{H}^{1}(\mathbb{T}^{d};\mathbb{R})$, we have that

$\displaystyle-\int_{\Omega}\nabla\cdot(A^{(1)}\nabla\chi^{(1)}_{\ell})v\;\>\mathsf{d}x$
$\displaystyle=\int_{\Omega}f^{(1)}_{\ell}v\;\>\mathsf{d}x$
$\displaystyle-\int_{\partial\Omega}vA^{(1)}\nabla\chi^{(1)}_{\ell}\cdot\widehat{n}\;\>\mathsf{d}x+\int_{\Omega}\nabla
v\cdot A^{(1)}\nabla\chi^{(1)}_{\ell}\;\>\mathsf{d}x$
$\displaystyle=\int_{\Omega}f^{(1)}_{\ell}v\;\>\mathsf{d}x.$

Since $v$, $A^{(1)}$, and the solution $\chi^{(1)}_{\ell}$ are all periodic on
$\mathbb{T}^{d}$, the first term is $0$. Combining with the equation for
$\chi^{(2)}_{\ell}$, we get

$\displaystyle\int_{\Omega}\nabla
v\cdot\left(A^{(1)}-A^{(2)}\right)\nabla\chi^{(1)}_{\ell}\;\>\mathsf{d}x$
$\displaystyle=\int_{\Omega}(f^{(1)}_{\ell}-f^{(2)}_{\ell})v+\nabla
v\cdot\left(A^{(2)}\left(\nabla\chi^{(2)}_{\ell}-\nabla\chi^{(1)}_{\ell}\right)\right)\>\mathsf{d}x.$

Setting $v=\chi^{(2)}_{\ell}-\chi^{(1)}_{\ell}$, we have

$\displaystyle\int_{\Omega}\left(\nabla\chi^{(2)}_{\ell}-\nabla\chi^{(1)}_{\ell}\right)\cdot\left(\left(A^{(1)}-A^{(2)}\right)\nabla\chi^{(1)}_{\ell}\right)\;\>\mathsf{d}x$
$\displaystyle=\int_{\Omega}(f^{(1)}_{\ell}-f^{(2)}_{\ell})\left(\chi^{(2)}_{\ell}-\chi^{(1)}_{\ell}\right)\;\>\mathsf{d}x,$
$\displaystyle+\int_{\Omega}\left(\nabla\chi^{(2)}_{\ell}-\nabla\chi^{(1)}_{\ell}\right)\cdot\left(A^{(2)}\left(\nabla\chi^{(2)}_{\ell}-\nabla\chi^{(1)}_{\ell}\right)\right)\;\>\mathsf{d}x,$
$\displaystyle\alpha\|\nabla\chi^{(2)}_{\ell}-\nabla\chi^{(1)}_{\ell}\|^{2}\leq\|A^{(1)}-A^{(2)}\|_{L^{\infty}}\|\nabla\chi^{(1)}_{\ell}\|_{L^{2}}\|\nabla\chi^{(2)}_{\ell}$
$\displaystyle-\nabla\chi^{(1)}_{\ell}\|_{L^{2}}+\|f^{(1)}_{\ell}-f^{(2)}_{\ell}\|_{\dot{H}^{-1}}\|\nabla\chi^{(2)}_{\ell}-\nabla\chi^{(1)}_{\ell}\|_{L^{2}},$
$\|\chi^{(2)}_{\ell}-\chi^{(1)}_{\ell}\|_{\dot{H}^{1}}\leq\frac{1}{\alpha}\left(\|A^{(1)}-A^{(2)}\|_{L^{\infty}}\|\nabla\chi^{(1)}_{\ell}\|_{L^{2}}+\|f^{(1)}_{\ell}-f^{(2)}_{\ell}\|_{\dot{H}^{-1}}\right).$
(A.2)

Evaluating,

$\displaystyle\|f^{(1)}_{\ell}-f^{(2)}_{\ell}\|_{\dot{H}^{-1}}$
$\displaystyle=\|\nabla\cdot A^{(1)}e_{\ell}-\nabla\cdot
A^{(2)}e_{\ell}\|_{\dot{H}^{-1}},$ (A.3)
$\displaystyle=\sup_{\|\xi\|_{\dot{H}^{1}}=1}\int_{\Omega}\xi\nabla\cdot(A^{(1)}-A^{(2)})e_{\ell}\;\>\mathsf{d}x,$
(A.4)
$\displaystyle\leq\sup_{\|\xi\|_{\dot{H}^{1}}=1}\|(A^{(1)}-A^{(2)})e_{\ell}\|_{L^{2}}\|\nabla\xi\|_{L^{2}},$
(A.5)
$\displaystyle\leq\|A^{(1)}-A^{(2)}\|_{L^{2}}\leq\|A^{(1)}-A^{(2)}\|_{L^{\infty}}.$
(A.6)

since our domain is $\mathbb{T}^{d}$. Combining this with (A.2) and the bound
of $\|\nabla\chi_{\ell}\|_{L^{2}}\leq\frac{\beta}{\alpha}$ obtained in the
first part of this proposition, we have

$\|\chi^{(2)}_{\ell}-\chi^{(1)}_{\ell}\|_{\dot{H}^{1}}\leq\frac{1}{\alpha}\left(1+\frac{\beta}{\alpha}\right)\|A^{(1)}-A^{(2)}\|_{L^{\infty}}.$
(A.7)

Returning to $d$ dimensions yields the result. ∎

The following result shows that the mapping $A\mapsto\bar{A}$ is continuous on
separable subspaces of $L^{\infty}(\mathbb{T}^{d};\mathbb{R}^{d\times d})$.

######  Lemma A.1.

Let $\mathcal{A}\subset L^{\infty}(\mathbb{T}^{d};\mathbb{R}^{d\times d})$ be
a separable subspace and $K\subset\mathcal{A}\cap\mathsf{PD}_{\alpha,\beta}$ a
closed set. Define the mapping $G:K\to\mathbb{R}^{d\times d}$ by
$A\mapsto\bar{A}$ as given by (1.3). Then there exists a continuous mapping
$\mathcal{G}\in C(\mathcal{A};\mathbb{R}^{d\times d})$ such that
$\mathcal{G}(A)=G(A)$ for any $A\in K$.

###### Proof.

Let $A^{(1)},A^{(2)}\in K$ then, by Proposition 1.1,

$\displaystyle\big{|}G\big{(}A^{(1)}\big{)}-G\big{(}A^{(2)}\big{)}\big{|}_{F}$
$\displaystyle\leq\int_{\mathbb{T}^{d}}|A^{(1)}-A^{(2)}|_{F}\big{(}1+|\nabla\chi^{(1)}|_{F}\big{)}\>\mathsf{d}y+\int_{\mathbb{T}^{d}}|A^{(2)}|_{F}|\nabla\chi^{(1)}-\nabla\chi^{(2)}|_{F}\>\mathsf{d}y$
$\displaystyle\leq\|A^{(1)}-A^{(2)}\|_{L^{\infty}}\big{(}1+\|\nabla\chi^{(1)}\|_{L^{2}}\big{)}+\|A^{(2)}\|_{L^{\infty}}\|\nabla\chi^{(1)}-\nabla\chi^{(2)}\|_{L^{2}}$
$\displaystyle\leq\left(1+\frac{\sqrt{d}}{\alpha}\left(\|A^{(1)}\|_{L^{\infty}}+\|A^{(2)}\|_{L^{\infty}}\left(\frac{\min\big{(}\|A^{(1)}\|_{L^{\infty}},\|A^{(2)}\|_{L^{\infty}}\big{)}}{\alpha}+1\right)\right)\right)$
$\displaystyle\;\;\;\;\cdot\|A^{(1)}-A^{(2)}\|_{L^{\infty}}$

hence $G\in C(K;\mathbb{R}^{d\times d})$. Applying the Tietze extension
theorem [20] to $G$ implies the existence of $\mathcal{G}$. ∎

The following lemma shows that, unfortunately, separable subspaces of
$L^{\infty}(\mathbb{T}^{d};\mathbb{R}^{d\times d})$ are not very useful.
Indeed, in the desired area of application of continuum mechanics, we ought to
be able to place a boundary of material discontinuity anywhere in the domain.
The following result shows that that would not be possible for a subset of
$\mathsf{PD}_{\alpha,\beta}$ which lies only in a separable subspace of
$L^{\infty}(\mathbb{T}^{d};\mathbb{R}^{d\times d})$.

######  Lemma A.2.

For any $t\in[0,1]$ define $c_{t}:[0,1]\to\mathbb{R}$ by

$c_{t}(x)=\begin{cases}1,&x\leq t\\\
0,&x>t\end{cases},\qquad\forall\>x\in[0,1].$

Define $E=\\{c_{t}:t\in[0,1]\\}\subset L^{\infty}([0,1])$. There exists no
separable subspace $\mathcal{A}\subset L^{\infty}([0,1])$ such that
$E\subseteq\mathcal{A}$.

###### Proof.

Suppose otherwise. Since $(\mathcal{A},\|\cdot\|_{L^{\infty}})$ is a separable
metric space, $(E,\|\cdot\|_{L^{\infty}})$ must be separable since
$E\subseteq\mathcal{A}$; this is a contradiction since
$(E,\|\cdot\|_{L^{\infty}})$ is not separable. To see this, let
$\\{c_{t_{j}}\\}_{j=1}^{\infty}$ be an arbitrary countable susbset of $E$.
Then for any $t\not\in\\{t_{j}\\}_{j=1}^{\infty}$, we have,

$\inf_{\\{t_{j}\\}_{j=1}^{\infty}}\|c_{t}-c_{t_{j}}\|_{L^{\infty}}=1.$

Hence no countable subset can be dense. ∎

Instead of working on a compact subset of a separable subspace of
$L^{\infty}(\mathbb{T}^{d};\mathbb{R}^{d\times d})$, we may instead try to
find a suitable probability measure which contains the discontinous functions
of interest. The following remarks makes clear why such an approch would still
be problematic for the purposes of approximation.

######  Remark A.3 (Gaussian Threshholding).

Let $\mu$ be a Gaussian measure on $L^{2}([0,1])$. Define

$T(x)=\begin{cases}1,&x\geq 0\\\ 0,&x<0\end{cases},\qquad\forall\>x\in[0,1]$

and consider the corresponding Nemytskii operator $N_{T}:L^{2}([0,1])\to
L^{\infty}([0,1])$. Then, working with the definitions in Lemma A.2, it is
easy to see that $E\subset\text{supp }{N_{T}}^{\sharp}\mu$. Therefore there
exists no separable subspace of $L^{\infty}([0,1])$ which contains $\text{supp
}{N_{T}}^{\sharp}\mu$.

We therefore abandon $L^{\infty}$ and instead show continuity and Lipschitz
continuity for some $L^{q}$ with $q<\infty$ to $\dot{H}^{1}$. The following
lemma is a general result for convergence of sequences in metric spaces which
is used in a more specific context in the next lemma.

######  Lemma A.4.

Let $(M,d)$ be a metric space and $(a_{n})\subset M$ a sequence. If every
subsequence $(a_{n_{k}})\subset(a_{n})$ contains a subsequence
$(a_{n_{k_{l}}})\subset(a_{n_{k}})$ such that $(a_{n_{k_{l}}})\to a\in M$ then
$(a_{n})\to a$.

###### Proof.

Suppose otherwise. Then, there exists some $\epsilon>0$ such that, for every
$N\in\mathbb{Z}^{+}$, there exists some $n=n(N)>N$ such that

$d(a_{n},a)\geq\epsilon.$

Then we can construct a subsequence $(a_{n_{j}})\subset(a_{n})$ such that
$d(a_{n_{j}},a)\geq\epsilon\forall n_{j}$. Therefore $a_{n_{j}}$ does not have
a subsequence converging to $a$, which is a contradiction. ∎

The following lemma proves existence of a limit in $L^{2}(D;\mathbb{R}^{d})$
of a sequence in $L^{\infty}(D;\mathbb{R}^{d\times d})$.

######  Lemma A.5.

Let $D\subseteq\mathbb{R}^{d}$ be an open set and $(A_{n})\subset
L^{\infty}(D;\mathbb{R}^{d\times d})$ a sequence satisfying the following.

  1. 1.

$A_{n}\in\mathsf{PD}_{\alpha,\beta}$ for all $n$,

  2. 2.

There exists $A\in L^{\infty}(D;\mathbb{R}^{d\times d})$ such that $(A_{n})\to
A\text{ in }L^{2}(D;\mathbb{R}^{d\times d}).$

Then, for any $g\in L^{2}(D;\mathbb{R}^{d})$, we have that $(A_{n}g)\to Ag$ in
$L^{2}(D;\mathbb{R}^{d})$.

###### Proof.

We have

$\|A_{n}g\|_{L^{2}}\leq\beta\|g\|_{L^{2}}$

hence $(A_{n}g)\subset L^{2}(D;\mathbb{R}^{d})$ and, similarly, by finite-
dimensional norm equivalence, there is a constant $C_{1}>0$ such that

$\|Ag\|_{L^{2}}\leq C_{1}\|A\|_{L^{\infty}}\|g\|_{L^{2}}$

hence $Ag\in L^{2}(D;\mathbb{R}^{d})$. Again, by finite-dimensional norm
equivalence, we that there exists a constant $C_{2}>0$ such that, for
$j\in\\{1,\dots,d\\}$ and almost every $y\in D$, we have

$(A_{n}g)_{j}(y)^{2}\leq|A_{n}^{(j)}(y)|^{2}|g(y)|^{2}\leq
C_{2}\beta^{2}|g(y)|^{2}$

where $A_{n}^{(j)}(y)$ denotes the $j$-th row of $A_{n}^{(j)}(y)$. In
particular,

$|(A_{n}g)_{j}(y)|\leq\sqrt{C_{2}}\beta|g(y)|.$

Let $(A_{n_{k}})\subset(A_{n})$ be an arbitrary subsequence. Since $(A_{n})\to
A$, we have that $(A_{n_{k}})\to A$ in $L^{2}(D;\mathbb{R}^{d\times d})$.
Therefore, there exists a subsequence $(A_{n_{k_{l}}})\subset(A_{n_{k}})$ such
that $A_{n_{k_{l}}}(y)\to A(y)$ for almost every $y\in D$. Therefore
$A_{n_{k_{l}}}(y)g(y)\to A(y)g(y)$ for almost every $y\in D$. Since $|g|\in
L^{2}(\mathbb{R}^{d})$, we have, by the dominated convergence theorem, that
$(A_{n_{k_{l}}}g)_{j}\to(Ag)_{j}$ in $L^{2}(D)$ for every
$j\in\\{1,\dots,d\\}$. Therefore $(A_{n_{k_{l}}}g)\to Ag$ in
$L^{2}(D;\mathbb{R}^{d})$. Since the subsequence $(A_{n_{k}})$ was arbitrary,
Lemma A.4 implies the result.

∎

Finally, we may prove Proposition 1.2.

See 1.2

###### Proof.

Consider the PDE

$-\nabla\cdot(A\nabla u)=\nabla\cdot Ae,\qquad y\in\mathbb{T}^{d}$ (A.8)

where $e$ is some standard basis vector of $\mathbb{R}^{d}$. Let
$(A_{n})\subset K$ be a sequence such that $(A_{n})\to A\in K$ in
$L^{2}(\mathbb{T}^{d};\mathbb{R}^{d\times d})$. Denote by
$u_{n}\in\dot{H}^{1}(\mathbb{T}^{d})$ the solution to (A.8) corresponding to
each $A_{n}$ and by $u\in\dot{H}^{1}(\mathbb{T}^{d})$ the solution
corresponding to the limiting $A$. A similar calculation as in the proof of
Proposition 1.1 shows

$\displaystyle\alpha\|u_{n}-u\|_{\dot{H}^{1}}^{2}$
$\displaystyle\leq\int_{\mathbb{T}^{d}}\langle(A-A_{n})(\nabla u+e),\nabla
u_{n}-\nabla u\rangle\>\mathsf{d}y$
$\displaystyle\leq\|u_{n}-u\|_{\dot{H}^{1}}\|(A_{n}-A)(\nabla u+e)\|_{L^{2}}.$

Since $\nabla u+e\in L^{2}(\mathbb{T}^{d};\mathbb{R}^{d})$, by Lemma A.5,
$\big{(}A_{n}(\nabla u+e)\big{)}\to A(\nabla u+e)$ in
$L^{2}(\mathbb{T}^{d};\mathbb{R}^{d})$ hence $(u_{n})\to u$ in
$\dot{H}^{1}(\mathbb{T}^{d})$. In particular, the mapping $A\mapsto u$ defined
by (A.8) is continuous. Since the problem (1.4) decouples as shown by (A.1),
we have that each component mapping $G_{l}:K\to\dot{H}^{1}(\mathbb{T}^{d})$
defined by $A\mapsto\chi_{\ell}$ is continuous thus $G$ is continuous.
Applying the Tietze extension theorem [20] to $G$ implies the existence of
$\mathcal{G}$. ∎

The following is a straightforward consequence of Proposition 1.2 that
establishes continuity of the map $A\mapsto\overline{A}$ defined in (1.3) as
well.

######  Lemma A.6.

Endow $\mathsf{PD}_{\alpha,\beta}$ with the
$L^{2}(\mathbb{T}^{d};\mathbb{R}^{d\times d})$ induced topology and let
$K\subset\mathsf{PD}_{\alpha,\beta}$ be a closed set. Define the mapping
$G:K\to\mathbb{R}^{d\times d}$ by $A\mapsto\bar{A}$ as given by (1.3). Then
there exists a bounded continuous mapping $\mathcal{G}\in
C\bigl{(}L^{2}(\mathbb{T}^{d};\mathbb{R}^{d\times d});\mathbb{R}^{d\times
d}\bigr{)}$ such that $\mathcal{G}(A)=G(A)$ for any $A\in K$.

###### Proof.

Since $\nabla:\dot{H}^{1}(\mathbb{T}^{d};\mathbb{R}^{d})\to
L^{2}(\mathbb{T}^{d};\mathbb{R}^{d\times d})$ is a bounded operator, Lemma 1.2
implies that the mapping $A\mapsto A+A\nabla\chi^{T}$ is continuous as
compositions, sums, and products of continuous functions are continuous. Now
let $A\in\mathsf{PD}_{\alpha,\beta}$ then $A\in
L^{1}(\mathbb{T}^{d};\mathbb{R}^{d\times d})$ since $A\in
L^{\infty}(\mathbb{T}^{d};\mathbb{R}^{d\times d})$. Thus

$\left|\int_{\mathbb{T}^{d}}A\>\mathsf{d}y\right|_{F}\leq\int_{\mathbb{T}^{d}}|A|_{F}\>\mathsf{d}y\leq\|A\|_{L^{2}}$

by Hölder’s inequality and the fact that
$\int_{\mathbb{T}^{d}}\>\mathsf{d}y=1$. Hence $G\in C(K;\mathbb{R}^{d\times
d})$ as a composition of continuous maps. Again applying the Tietze extension
theorem [20] to $G$ implies the existence of $\mathcal{G}$. ∎

To prove Proposition 1.3, we need to establish Lipschitz continuity. We first
establish the following result, which is similar to the one proved in [11] in
Theorem 2.1. We show it again here both for completeness and because we
specify to the case of the cell problem (1.4) with periodic boundary
conditions rather than the system (1.1) with Dirichlet boundary conditions.

######  Lemma A.7.

If $A^{(1)},A^{(2)}\in\mathsf{PD}_{\alpha,\beta}$ and let
$\chi^{(1)},\chi^{(2)}$ be the corresponding solutions to (1.4). Then

$\|\chi^{(1)}-\chi^{(2)}\|_{\dot{H}^{1}}\leq\frac{\sqrt{d}}{\alpha}\left(\|A^{(2)}-A^{(1)}\|_{L^{2}}+\|\nabla\chi^{(2)}\|_{L^{p}}\|A^{(2)}-A^{(1)}\|_{L^{q}}\right)$
(A.9)

for $p\geq 2$ and $q=\frac{2p}{p-2}$.

###### Proof.

As in the proof of Proposition 1.1, we denote $f^{(i)}=\nabla\cdot A^{(i)}$
for $i\in\\{1,2\\}$ for simplicity of notation and to be easily comparable to
the proof of Theorem 2.1 in [11]. Since both sides of the cell problem
equation (1.4) depend on $A^{(i)}$, we introduce $\widetilde{\chi}$ as the
solution of

$\displaystyle-\nabla\cdot\left(\nabla\widetilde{\chi}A^{(2)}\right)$
$\displaystyle=\nabla\cdot
A^{(1)},\quad\widetilde{\chi}\in\dot{H}^{1}(\Omega)$ (A.10)

as an intermediate function. We obtain bounds using $\widetilde{\chi}$ and
apply the triangle inequality to

$\|(\chi^{(1)}-\widetilde{\chi})+(\widetilde{\chi}-\chi^{(2)})\|_{\dot{H}^{1}}$

to obtain a bound on $\|\chi^{(1)}-\chi^{(2)}\|_{\dot{H}^{1}}$. From the naïve
perturbation bound in (A.2) we have

$\|\widetilde{\chi}_{\ell}-\chi^{(2)}_{\ell}\|_{\dot{H}^{1}}\leq\frac{1}{\alpha}\|f^{(1)}-f^{(2)}\|_{\dot{H}^{-1}},$

so we are left to bound
$\|\chi^{(1)}_{\ell}-\widetilde{\chi}_{\ell}\|_{\dot{H}^{1}}$. We note that

$\displaystyle\nabla\cdot\left(A^{(2)}\nabla\widetilde{\chi}_{\ell}\right)$
$\displaystyle=\nabla\cdot\left(A^{(1)}\nabla\chi^{(1)}_{\ell}\right)$
$\displaystyle\int_{\Omega}A^{(2)}\nabla\widetilde{\chi}_{\ell}\cdot\nabla
v\;\>\mathsf{d}x$
$\displaystyle=\int_{\Omega}A^{(1)}\nabla\chi^{(1)}_{\ell}\cdot\nabla
v\>\mathsf{d}x\quad\forall v\in\dot{H}^{1}(\Omega)$

Letting $v=\chi^{(1)}_{\ell}-\widetilde{\chi}_{\ell}$,

$\displaystyle\int_{\Omega}A^{(2)}\nabla\widetilde{\chi}_{\ell}\cdot\left(\nabla\chi^{(1)}_{\ell}-\nabla\widetilde{\chi}_{\ell}\right)\;\>\mathsf{d}x$
$\displaystyle=\int_{\Omega}A^{(1)}\nabla\chi^{(1)}_{\ell}\cdot\left(\nabla\chi^{(1)}_{\ell}-\nabla\widetilde{\chi}_{\ell}\right)\;\>\mathsf{d}x$
$\displaystyle\int_{\Omega}A^{(2)}\left(\nabla\widetilde{\chi}_{\ell}-\nabla\chi^{(1)}_{\ell}\right)\cdot\left(\nabla\widetilde{\chi}_{\ell}-\nabla\chi^{(1)}_{\ell}\right)\;\>\mathsf{d}x$
$\displaystyle=\int_{\Omega}\left(A^{(2)}-A^{(1)}\right)\nabla\chi^{(1)}_{\ell}\cdot\left(\nabla\chi^{(1)}_{\ell}-\nabla\widetilde{\chi}_{\ell}\right)\;\>\mathsf{d}x$
$\displaystyle\alpha\|\widetilde{\chi}_{\ell}-\chi^{(1)}_{\ell}\|_{\dot{H}^{1}}$
$\displaystyle\leq\|(A^{(2)}-A^{(1)})(\nabla\chi^{(1)}_{\ell})\|_{L^{2}}$

Applying Hölder for $L^{2}$, we get

$\|\widetilde{\chi}_{\ell}-\chi^{(1)}_{\ell}\|_{\dot{H}^{1}}\leq\frac{1}{\alpha}\|\nabla\chi^{(1)}_{\ell}\|_{L^{p}}\|A^{(2)}-A^{(1)}\|_{L^{q}}$
(A.11)

for $q=\frac{2p}{p-2}$ where $p\in[2,\infty]$. Putting the two parts together,
we have that

$\displaystyle\|\chi^{(2)}_{\ell}-\chi^{(1)}_{\ell}\|_{\dot{H}^{1}}$
$\displaystyle\leq\frac{1}{\alpha}\|\nabla\cdot A^{(2)}e_{\ell}-\nabla\cdot
A^{(1)}e_{\ell}\|_{\dot{H}^{-1}}+\frac{1}{\alpha}\|\nabla\chi^{(1)}_{\ell}\|_{L^{p}}\|A^{(2)}-A^{(1)}\|_{L^{q}}$
$\displaystyle\leq\frac{1}{\alpha}\|A^{(2)}-A^{(1)}\|_{L^{2}}+\frac{1}{\alpha}\|\nabla\chi^{(1)}_{\ell}\|_{L^{p}}\|A^{(2)}-A^{(1)}\|_{L^{q}}$

Combining bounds for all $d$ dimensions yields the result. ∎

######  Remark A.8.

Since $L^{q}(\Omega)\hookrightarrow L^{2}(\Omega)$ for bounded
$\Omega\subset\mathbb{R}^{d}$ and $q\geq 2$, we could also write the bound of
Lemma A.7 as

$\|\chi^{(2)}_{\ell}-\chi^{(1)}_{\ell}\|_{\dot{H}^{1}}\leq\frac{1}{\alpha}\left(C+\|\nabla\chi^{(1)}_{\ell}\|_{L^{p}}\right)\|A^{(2)}-A^{(1)}\|_{L^{q}}$

for some $C$ dependent only on $q$ and $\Omega$.

The result of Lemma A.7 is unhelpful if $\|\nabla\chi\|_{L^{p}}$ is unbounded,
as it is for the case of $p=2$ with sets of $A$ containing discontinuous
corner interfaces appearing in the microstructure examples of square
inclusions and Voronoi crystals as described in Section 2.222To see this,
consider local solutions $u:\mathbb{R}^{2}\to\mathbb{R}^{2}$ of
$\nabla\cdot(a\nabla u)=0$ where $a=a(\theta)$ is
$a(\theta)=\begin{cases}a_{-},\quad 0<\theta<\theta^{*}\\\
a_{+},\quad\theta^{*}<\theta<2\pi\end{cases}$ with the ansatz
$u(r,\theta)=h(r)g(\theta)$. The solution component $h(r)$ takes the form
$h(r)=r^{b},\;b>0$, whose gradient is singular for $0<b<1$. Solving the
associated eigenvalue problem for $b$ gives a singular component for
$\theta^{*}\neq\pi$. In this setting, it is not possible for Lemma A.7 to
result in Lipschitz continuity as a map from $L^{2}$ to $\dot{H}^{1}$.
Instead, we seek to bound $\|\nabla\chi\|_{L^{p}}$ for $p$ satisfying
$2<p<\infty$.

Before continuing, we establish a bound on the gradient of the solution to the
Poisson equation on the torus. This follows the strategy of [11] for the
Dirichlet problem. In order to avoid extra factors of $2\pi$ in all formulae
we assume $\mathbb{T}^{d}=[0,2\pi]^{d}$ with opposite faces identified. As we
work on the torus, it is useful to first set up notation for the function
spaces of interest. Let

$\mathcal{D}(\mathbb{T}^{d})=C^{\infty}_{c}(\mathbb{T}^{d})=C^{\infty}(\mathbb{T}^{d})$

be the space of test functions where the last equality follows from
compactness of the torus. Functions can be either $\mathbb{R}$ or $\mathbb{C}$
valued hence we do not explicitly specify the range. We equip
$\mathcal{D}(\mathbb{T}^{d})$ with a locally convex topology generated by an
appropriate family of semi-norms, see, for example, [51, Section 3.2.1]. Any
function $g\in\mathcal{D}(\mathbb{T}^{d})$ can be represented by its Fourier
series

$g(x)=\sum_{k\in\mathbb{Z}^{d}}\widehat{g}(k)\text{e}^{ix\cdot k}$

where $\widehat{g}$ denotes the Fourier transform of $g$ and convergence of
the r.h.s. sum is with respect to the topology of
$\mathcal{D}(\mathbb{T}^{d})$. It holds that
$\widehat{g}\in\mathcal{S}(\mathbb{Z}^{d})$, the Schwartz space of rapidly
decreasing functions on the integer lattice, so we have

$|\widehat{g}(k)|\leq c_{m}(1+|k|)^{-m},\quad m=0,1,\dots$

for some constants $c_{m}$. We may then define the topological (continuous)
dual space of $\mathcal{D}(\mathbb{T}^{d})$, the space of distributions,
denoted $\mathcal{D}^{\prime}(\mathbb{T}^{d})$, which can be described as
follows: the condition that $f\in\mathcal{D}^{\prime}(\mathbb{T}^{d})$ is
characterized by the property

$|\widehat{f}(k)|\leq b_{m}(1+|k|)^{m},\quad m=0,1,\dots$

for some constants $b_{m}$. We take the weak${\text{-}}^{*}$ topology on
$\mathcal{D}^{\prime}(\mathbb{T}^{d})$ and generally use the prime notation
for any such dual space. For any $-\infty<s<\infty$, we define the fractional
Laplacian as

$(-\Delta)^{s}f=\sum_{k\in\mathbb{Z}^{d}\setminus\\{0\\}}|k|^{2s}\widehat{f}(k)\text{e}^{ik\cdot
x}$ (A.12)

where the r.h.s. sum converges in the topology of
$\mathcal{D}^{\prime}(\mathbb{T}^{d})$. It is easy to see that
$(-\Delta)^{s}:\mathcal{D}^{\prime}(\mathbb{T}^{d})\to\mathcal{D}^{\prime}(\mathbb{T}^{d})$
is continuous. Furthermore, for any $j\in\\{1,\dots,d\\}$, we define the
family of operators
$\tilde{R}_{j}:\mathcal{D}^{\prime}(\mathbb{T}^{d})\to\mathcal{D}^{\prime}(\mathbb{T}^{d})$,
defining periodic Riesz transforms, by

$\tilde{R}_{j}f=\sum_{k\in\mathbb{Z}^{d}}-\frac{ik_{j}}{|k|}\hat{f}(k)\text{e}^{ik\cdot
x}$ (A.13)

where we identify $\frac{k_{j}}{|k|}|_{k=0}=\lim_{|k|\to
0}\frac{k_{j}}{|k|}=0$. Again, we stress that convergence of the r.h.s. sum is
in the topology of $\mathcal{D}^{\prime}(\mathbb{T}^{d})$. Lastly, we denote
by $\mathcal{S}(\mathbb{R}^{d})$ and $\mathcal{S}^{\prime}(\mathbb{R}^{d})$
the Schwartz space and the space of tempered distributions on $\mathbb{R}^{d}$
respectively; see, for example, [52, Chapter 1] for the precise definitions.

The following lemma establishes boundedness of the periodic Riesz transform on
$L^{p}(\mathbb{T}^{d})$. It is essential in proving boundedness of the
gradient to the solution of the Poisson equation on the torus. The result is
essentially proven in [52]. We include it here, in our specific torus setting,
giving the full argument for completeness.

######  Lemma A.9.

There exists a constant $c=c(d,p)>0$ such that, for any $j\in\\{1,\dots,d\\}$
and any $f\in L^{p}(\mathbb{T}^{d})$ for some $2\leq p<\infty$, we have

$\|\tilde{R}_{j}f\|_{L^{p}(\mathbb{T}^{d})}\leq
c\|f\|_{L^{p}(\mathbb{T}^{d})}.$

###### Proof.

Let $g\in L^{2}(\mathbb{R}^{d})\cap L^{p}(\mathbb{R}^{d})$ for some
$1<p<\infty$. For any $j\in\\{1,\dots,d\\}$, define the family of operators
$R_{j}$ by

$(R_{j}g)(x)=\lim_{\delta^{-1},\epsilon\to
0^{+}}\int_{\delta\geq|t|\geq\epsilon}g(x-t)K_{j}(t)\>dt,$

where

$K_{j}(t)=\frac{\Gamma\big{(}(d+1)/2\big{)}t_{j}}{\pi^{(d+1)/2}|t|^{d+1}}$

and $\Gamma$ denotes the Euler-Gamma function. By [52, Chapter 4, Theorem
4.5], $K_{j}\in\mathcal{S}^{\prime}(\mathbb{R}^{d})$ and its Fourier transform
satisfies

$\hat{K}_{j}(t)=-\frac{it_{j}}{|t|}$

where $i=\sqrt{-1}$ denotes the imaginary unit. Therefore, for any
$\phi\in\mathcal{S}(\mathbb{R}^{d})$, we have

$(K_{j}*\phi)^{\widehat{}}\>(t)=-\frac{it_{j}}{|t|}\hat{\phi}(t)$

where $*$ denotes convolution, see, for example, [52, Chapter 1, Theorem
3.18]. Since $g\in L^{2}(\mathbb{R}^{d})$, we therefore find that, by [52,
Chapter 6, Theorem 2.6],

$(R_{j}g)^{\widehat{}}\>(x)=-\frac{ix_{j}}{|x|}\hat{g}(x)$ (A.14)

for Lebesgue almost every $x\in\mathbb{R}^{d}$. The result [52, Chapter 6,
Theorem 2.6] further shows that there exists a constant $c=c(d,p)>0$ such that

$\|R_{j}g\|_{L^{p}(\mathbb{R}^{d})}\leq c\|g\|_{L^{p}(\mathbb{R}^{d})}.$

We note from (A.14) and the definition (A.13) that $\tilde{R}_{j}$ may be
viewed as $R_{j}$ with the restriction of the Fourier multiplier
$-\frac{ix_{j}}{|x|}$ to the lattice $\mathbb{Z}^{d}$. We can therefore use
the transference theory of [52] to establish boundedness of $\tilde{R}_{j}$
from the boundedness of $R_{j}$. In particular, note that the mapping
$x\mapsto-\frac{ix_{j}}{|x|}$ is continuous at all $x\in\mathbb{R}^{d}$ except
$x=0$. However, by symmetry, we have that, for all $\epsilon>0$

$\int_{|x|\leq\epsilon}-\frac{ix_{j}}{|x|}\>dx=0.$

Therefore we can apply [52, Chapter 7, Theorem 3.8, Corollary 3.16] to
conclude that, since $R_{j}$ is bounded from $L^{p}(\mathbb{R}^{d})$ to
$L^{p}(\mathbb{R}^{d})$, $\tilde{R}_{j}$ is bounded from
$L^{p}(\mathbb{T}^{d})$ to $L^{p}(\mathbb{T}^{d})$ with

$\|\tilde{R}_{j}\|_{L^{p}(\mathbb{T}^{d})\to
L^{p}(\mathbb{T}^{d})}\leq\|R_{j}\|_{L^{p}(\mathbb{R}^{d})\to
L^{p}(\mathbb{R}^{d})}.$

This implies the desired result. ∎

We define the Bessel potential spaces by

$L^{s,p}(\mathbb{T}^{d})=\\{u\in\mathcal{D}^{\prime}(\mathbb{T}^{d})\;|\;\|u\|_{L^{s,p}(\mathbb{T}^{d})}\coloneqq\|(I-\Delta)^{s/2}u\|_{L^{p}(\mathbb{T}^{d})}<\infty\\}$

for any $-\infty<s<\infty$ and $1<p<\infty$. We also define the homogeneous
version of these spaces, sometimes called the Riesz potential spaces, by

$\dot{L}^{s,p}(\mathbb{T}^{d})=\\{u\in\mathcal{D}^{\prime}(\mathbb{T}^{d})\;|\;\|u\|_{\dot{L}^{s,p}(\mathbb{T}^{d})}\coloneqq\|(-\Delta)^{s/2}u\|_{L^{p}(\mathbb{T}^{d})}<\infty,\;\int_{\mathbb{T}^{d}}u\>\mathsf{d}y=0\\}.$

It is clear that $\dot{L}^{s,p}(\mathbb{T}^{d})\subset
L^{s,p}(\mathbb{T}^{d})$ is closed subspace. We then have the following result
for the Poisson equation.

######  Lemma A.10.

For each $f\in L^{s,p}(\mathbb{T}^{d})$, for $-\infty<s<\infty$ and $2\leq
p<\infty$, the solution $u$ of the equation

$-\Delta u=f,\quad u\text{ 1-periodic, }\int_{\mathbb{T}^{d}}u\>\mathsf{d}x=0$
(A.15)

satisfies

$\|\nabla u\|_{\dot{L}^{s+1,p}(\mathbb{T}^{d})}\leq
K\|f\|_{\dot{L}^{s,p}(\mathbb{T}^{d})}$ (A.16)

for some finite $K>0$ depending only on $p$ and $d$.

###### Proof.

From the definitions (A.12) and (A.13), it is easy to see that the Riesz
transform can be written as

$\tilde{R}_{j}=-\partial_{x_{j}}(-\Delta)^{-1/2}$

in the sense of distributions. Consider now equation (A.15) with $f\in
L^{s,p}(\mathbb{T}^{d})$ for $2\leq p<\infty$. We have that

$\displaystyle\|\partial_{x_{j}}u\|_{\dot{L}^{s+1,p}(\mathbb{T}^{d})}$
$\displaystyle=\|\partial_{x_{j}}(-\Delta)^{-1}f\|_{\dot{L}^{s+1,p}(\mathbb{T}^{d})}$
$\displaystyle=\|\partial_{x_{j}}(-\Delta)^{-1/2}(-\Delta)^{s/2}f\|_{L^{p}(\mathbb{T}^{d})}$
$\displaystyle=\|\tilde{R}_{j}(-\Delta)^{s/2}f\|_{L^{p}(\mathbb{T}^{d})}.$

It is clear that

$\|(-\Delta)^{s/2}f\|_{L^{p}(\mathbb{T}^{d})}=\|f\|_{\dot{L}^{s,p}(\mathbb{T}^{d})}<\infty$

hence $(-\Delta)^{s/2}f\in L^{p}(\mathbb{T}^{d})$. We can thus apply Lemma A.9
to find a constant $c=c(d,p)>0$ such that

$\|\partial_{x_{j}}u\|_{\dot{L}^{s+1,p}(\mathbb{T}^{d})}\leq
c\|(-\Delta)^{s/2}f\|_{L^{p}(\mathbb{T}^{d})}=c\|f\|_{\dot{L}^{s,p}(\mathbb{T}^{d})}.$

The result follows by finite-dimensional norm equivalence.

∎

Next we define the homogeneous Sobolev spaces on the torus as

$\dot{W}^{k,p}(\mathbb{T}^{d})=\\{u\in W^{k,p}(\mathbb{T}^{d})\;|\;u\text{ is
}2\pi\text{-periodic },\;\int_{\mathbb{T}^{d}}u\>\mathsf{d}y=0\\}$ (A.17)

for $k=0,1,\dots,$ and $1\leq p\leq\infty$ with the standard norm on
$W^{k,p}$, see, for example [1].

######  Remark A.11.

By [51, Section 3.5.4], we have that, for any $k=0,1,\dots$ and $1<p<\infty$,

$L^{k,p}(\mathbb{T}^{d})=W^{k,p}(\mathbb{T}^{d}),\qquad\dot{L}^{k,p}(\mathbb{T}^{d})=\dot{W}^{k,p}(\mathbb{T}^{d}).$

Furthermore, by [51, Section 3.5.6],

$\displaystyle W^{-k,p^{\prime}}(\mathbb{T}^{d})$
$\displaystyle=\big{(}W^{k,p}(\mathbb{T}^{d})\big{)}^{\prime}=\big{(}L^{k,p}(\mathbb{T}^{d})\big{)}^{\prime}=L^{-k,p^{\prime}}(\mathbb{T}^{d}),$
$\displaystyle\dot{W}^{-k,p^{\prime}}(\mathbb{T}^{d})$
$\displaystyle=\big{(}\dot{W}^{k,p}(\mathbb{T}^{d})\big{)}^{\prime}=\big{(}\dot{L}^{k,p}(\mathbb{T}^{d})\big{)}^{\prime}=\dot{L}^{-k,p^{\prime}}(\mathbb{T}^{d})$

where $p^{\prime}$ is the Hölder conjugate of $p$ i.e. $1/p+1/p^{\prime}=1$.

In the following, we use the notation

$[K_{0},K_{1}]_{\theta,q}$ (A.18)

to denote the real interpolation between two Banach spaces continuously
embedded in the same Hausdorff topological space, as described in [1]. We also
need Lemma A1 from [22], which we have copied below as Lemma A.12 to ease
readability. Although this lemma was written only for $q=2$, the result still
holds for our $q>2$ with a very similar proof.

######  Lemma A.12.

Let $E_{1}\subset E_{0}$ be two Banach spaces with $E_{1}$ continuously
embedded in $E_{0}$. Let $T:E_{j}\to E_{j}$ be a bounded operator with closed
range and assume that $T$ is a projection, $j\in\\{0,1\\}$. Denote by $K_{0}$
and $K_{1}$ the ranges of $T|_{E_{0}}$ and $T|_{E_{1}}$ respectively. Then the
following two spaces coincide with equivalent norms:

$[K_{0},K_{1}]_{\theta,q}=[E_{0},E_{1}]_{\theta,q}\cap K_{0}\quad\forall
s\in(0,1).$

We now state the result for the bound on $\|\nabla\chi\|_{L^{p}}$.

######  Lemma A.13.

Let $\chi$ solve (1.4) for $A\in\mathsf{PD}_{\alpha,\beta}$. Then

$\|\nabla\chi\|_{L^{p}}\leq\frac{1}{\beta}\frac{K^{\eta^{*}(p)}}{1-K^{\eta^{*}(p)}\left(1-\frac{\alpha}{\beta}\right)}$
(A.19)

for $2<p<p^{*}\left(\frac{\alpha}{\beta}\right)$ where

$p^{*}(t):=\text{max}\left\\{p\;|\;K^{-\eta(p)}\geq 1-t,\;2<p<Q\right\\}$
(A.20)

for $\eta(p)=\frac{1/2-1/p}{1/2-1/Q}$and $K=K(d,Q)$ is the constant in Lemma
A.10, for any choice of $Q>2$.

###### Proof.

The operator $T=-\Delta$ is invertible from $\dot{H}^{-1}$ to $\dot{H}^{1}$,
and the inverse $T^{-1}$ is bounded with norm $1$ since the Poisson equation
with periodic boundary conditions has a unique solution in $\dot{H}^{1}$ for
$f\in H^{-1}$ with bound $\|u\|_{\dot{H}^{1}}\leq\|f\|_{\dot{H}^{-1}}$. From
Lemma A.10 it is also bounded with norm $K$ from $\dot{W}^{-1,Q}$ to
$\dot{W}^{1,Q}$ for any $Q>2$. By the real method of interpolation [1], we
have that

$W^{1,p}=\left[H^{1},W^{1,Q}\right]_{\eta(p),p}$ (A.21)

using the notation of [1] where $\eta(p)=\frac{1/2-1/p}{1/2-1/Q}$. From the
duality theorem (Theorem 3.7.1. of [6]), we have that

$\left[H^{-1},W^{-1,Q}\right]_{\eta(p),p}=\left(\left[H^{1},W^{1,Q^{\prime}}\right]_{\eta(p),p^{\prime}}\right)^{\prime}$
(A.22)

From real interpolation, the right hand side equals
$(W^{1,p^{\prime}})^{\prime}=W^{-1,p}$ in our notation. Therefore, we have the
necessary dual statement that parallels (A.21):

$W^{-1,p}=\left[H^{-1},W^{-1,Q}\right]_{\eta(p),p}.$ (A.23)

These interpolation spaces are not yet restricted to the space of functions on
$\Omega$ with periodic boundary conditions. Using the projection onto the
space of continuous functions on $\Omega$ with periodic boundary conditions as
$T$ and noticing that $W^{1,Q}\hookrightarrow H^{1}$, we apply Lemma A.12 with
$K_{0}=\dot{H}^{1}$ and have

$\dot{W}^{1,s}=[\dot{H}^{1},\dot{W}^{1,Q}]_{\eta(p),p}.$ (A.24)

The duality theorem still holds in this setting, so we also have

$\dot{W}^{-1,p}=\left[\dot{H}^{-1},\dot{W}^{-1,Q}\right]_{\eta(p),p}.$ (A.25)

Using the exact interpolation theorem, Theorem 7.23 of [1], $T^{-1}$ is also a
bounded map from $\dot{W}^{-1,p}$ to $\dot{W}^{1,p}$ with norm $K^{\eta(p)}$:

$\|T^{-1}f\|_{\dot{W}^{1,p}}\leq K^{\eta(p)}\|f\|_{W^{-1,p}}.$ (A.26)

The remainder of the proof is identical to that of the proof of Proposition
$1$ in [11], but we state it here in our notation for completeness. Define
$S$: $\dot{W}^{1,p}\to W^{-1,p}$ as the operator
$Su=-\nabla\cdot\left(\frac{1}{\beta}A\nabla u\right)$. Let $V$ be the
perturbation operator $V:=T-S$. Since $A\in\mathsf{PD}_{\alpha,\beta}$, we
have $\|S\|\leq 1$ and $\|V\|\leq 1-\frac{\alpha}{\beta}$. Therefore, as a
mapping from $\dot{W}^{1,p}$ to $\dot{W}^{-1,p}$,

$\|T^{-1}V\|\leq\|T^{-1}\|\|V\|\leq
K^{\eta(p)}\left(1-\frac{\alpha}{\beta}\right)$ (A.27)

Since $T$ is invertible, $S=T(I-T^{-1}V)$ is invertible provided
$K^{\eta(p)}\left(1-\frac{\alpha}{\beta}\right)<1$. Moreover, as a mapping
from $W^{-1,p}$ to $W^{1,p}$,

$\|S^{-1}\|\leq\|(I-T^{-1}V)^{-1}\|\|T^{-1}\|\leq\frac{K^{\eta(p)}}{1-K^{\eta(p)}\left(1-\frac{\alpha}{\beta}\right)}.$
(A.28)

Therefore,

$\|\nabla
u\|_{L^{p}}=\|u\|_{\dot{W}^{1,p}}\leq\frac{1}{\beta}\frac{K^{\eta(p)}}{1-K^{\eta(p)}\left(1-\frac{\alpha}{\beta}\right)}$
(A.29)

for provided $K^{\eta(p)}\left(1-\frac{\alpha}{\beta}\right)<1$. The bound and
specified range of $p$ follow. ∎

Finally, we may prove Proposition 1.3

See 1.3

###### Proof.

Lemma A.13 guarantees a $p_{0}>2$ such that $\|\nabla\chi^{(2)}\|_{L^{p}}$ in
Lemma A.7 is bounded above by a constant for $2<p<p_{0}$. Then Lemma A.7 gives
Lipschitz continuity of the solution map from
$L^{q}(\mathbb{T}^{d})\mapsto\dot{H}^{1}(\mathbb{T}^{d})$ for $q$ satisfying
$q_{0}<q<\infty$ for some $q_{0}>2$. ∎

######  Remark A.14.

From the results of Lemma A.13 and Lemma A.7, we have that we can take
$q_{0}=\frac{2p_{0}}{p_{0}-2}$ where

$p_{0}=\max\\{p\;|\;K^{-1+2/p}\geq 1-t,\;2<p<\infty\\}.$

Therefore, bounds on $p_{0}$ may be inherited from bounds on $K$ that appears
in Lemma A.10.

We can leverage the result of Proposition (1.3) to also show continuity in the
map $A\mapsto\overline{A}$ from $L^{q}$ to $\dot{H}^{1}$ in equation (1.3).

######  Lemma A.15.

If $A\in\mathsf{PD}_{\alpha,\beta}$ in (1.4), there exists $q_{0}<\infty$ such
that for all $q$ satisfying $q_{0}<q\leq\infty$, the solution map
$A\mapsto\overline{A}$ of equations (1.3) and 1.4 is Lipschitz-continuous as a
map from $L^{q}(\mathbb{T}^{d})$ to $\mathbb{R}^{d\times d}$.

###### Proof.

Recall equation (1.3):

$\overline{A}=\int_{\mathbb{T}^{d}}\left(A(y)+A(y)\nabla\chi(y)^{T}\right)\>\mathsf{d}y$

For two different coefficient functions
$A^{(1)},A^{(2)}\in\mathsf{PD}_{\alpha,\beta}$, $\chi^{(1)},\chi^{(2)}$ the
associated solutions of (1.4), and $\overline{A}^{(1)},\overline{A}^{(2)}$ the
associated homogenized coefficients of (1.3), we may write

$\displaystyle\overline{A}^{(1)}-\overline{A}^{(2)}$
$\displaystyle=\int_{\mathbb{T}^{d}}A^{(1)}-A^{(2)}+A^{(1)}\nabla(\chi^{(1)})^{T}-A^{(2)}\nabla(\chi^{(2)})^{T}\>\mathsf{d}y$
(A.30) $\displaystyle\left|\overline{A}^{(1)}-\overline{A}^{(2)}\right|_{F}$
$\displaystyle\leq\left|\int_{\mathbb{T}^{d}}A^{(1)}-A^{(2)}\>\mathsf{d}y\right|_{F}+\left|\int_{\mathbb{T}^{d}}A^{(1)}\nabla\chi^{(1)}-A^{(1)}\nabla\chi^{(2)}+A^{(1)}\nabla\chi^{(2)}-A^{(2)}\nabla\chi^{(2)}\>\mathsf{d}y\right|_{F}$
(A.31)
$\displaystyle\leq\|A^{(1)}-A^{(2)}\|_{L^{2}}+\left|\int_{\mathbb{T}^{d}}A^{(1)}\left(\nabla\chi^{(1)}-\nabla\chi^{(2)}\right)+\nabla\chi^{(2)}\left(A^{(1)}-A^{(2)}\right)\right|_{F}$
(A.32)
$\displaystyle\leq\|A^{(1)}-A^{(2)}\|_{L^{2}}+\|\chi^{(1)}-\chi^{(2)}\|_{\dot{H}^{1}}\|A^{(1)}\|_{L^{2}}+\|\nabla\chi^{(2)}\|_{L^{2}}\|A^{(1)}-A^{(2)}\|_{L^{2}}$
(A.33)
$\displaystyle\leq\left(1+\frac{\sqrt{d}\beta}{\alpha}\right)\|A^{(1)}-A^{(2)}\|_{L^{2}}+C\beta\|A^{(1)}-A^{(2)}\|_{L^{q}}$
(A.34)

where $C$ is the Lipschitz constant from Proposition 1.3. The embedding
$L^{q}\hookrightarrow L^{2}$ gives Lipschitz continuity from $L^{q}$ to
$\mathbb{R}^{d\times d}$ of the map. ∎

##  Appendix B Proofs of Approximation Theorems

In this section we prove the approximation theorems stated in Section 3.

See 3.3

###### Proof.

By Proposition 1.2, there exists a continuous map $\mathcal{G}\in
C(L^{2}(\mathbb{T}^{d};\mathbb{R}^{d\times
d});\dot{H}^{1}(\mathbb{T}^{d};\mathbb{R}^{d}))$ such that
$\mathcal{G}(A)=G(A)$ for any $A\in K$. By [30, Theorem 2.5], there exists a
FNO $\Psi:L^{2}(\mathbb{T}^{d};\mathbb{R}^{d\times
d})\to\dot{H}^{1}(\mathbb{T}^{d};\mathbb{R}^{d})$ such that

$\sup_{A\in K}\|\mathcal{G}(A)-\Psi(A)\|_{\dot{H}^{1}}<\epsilon.$

Therefore

$\sup_{A\in K}\|G(A)-\Psi(A)\|_{\dot{H}^{1}}=\sup_{A\in
K}\|\mathcal{G}(A)-\Psi(A)\|_{\dot{H}^{1}}<\epsilon$

as desired. ∎

See 3.4

###### Proof.

The result follows as in Theorem 3.3 by applying Lemma A.6 instead of
Proposition 1.2. ∎

See 3.5

###### Proof.

By Proposition 1.3, there exists a Lipschitz-continuous map $G\in
C(L^{q}(\mathbb{T}^{d};\mathbb{R}^{d\times
d});\dot{H}^{1}(\mathbb{T}^{d};\mathbb{R}^{d}))$. By [31, Theorem 11], there
exists a NO $\Psi:L^{q}(\mathbb{T}^{d};\mathbb{R}^{d\times
d})\to\dot{H}^{1}(\mathbb{T}^{d};\mathbb{R}^{d})$ such that

$\sup_{A\in K}\|G(A)-\Psi(A)\|_{\dot{H}^{1}}<\epsilon.$ (B.1)

∎

See 3.6

###### Proof.

The result follows as in Theorem 3.3 by applying Lemma A.15 instead of
Proposition 1.3. ∎

##  Appendix C Proofs for Microstructure Examples

######  Lemma C.1.

For $d=2$, let $\mathcal{A}\subset\mathsf{PD}_{\alpha,\beta}$ be a set of
piecewise-constant functions on $\mathbb{T}^{d}$ whose level sets are
nontrivial for at most $M$ constants, and each level set consists of the union
of at most $M$ convex polygons. The set $\mathcal{A}$ is compact in
$L^{2}(\Omega;\mathbb{R}^{d\times d})$.

###### Proof.

For arbitrary $\epsilon>0$, we identify a finite $\epsilon$-net, denoted
$\mathcal{N}_{\epsilon}$, for $\mathcal{A}$. Partition the unit square by a
uniform grid of $2^{2\ell}$ boxes for integer $\ell>1$, and denote by
$B_{\ell}$ this set of boxes. Let $\mathcal{N}_{\epsilon}$ be the set of
piecewise-constant functions that are constant over each box and take values
only in the set $\\{z:\;z=\alpha+\eta
k,k\in\\{0,1,\dots,\lceil{\frac{\beta-\alpha}{\eta}\rceil}\\}\\}$, for
$\eta<\epsilon$.

Since each function takes at most $M$ values on at most $M$ convex polygons,
we have naive upper bounds of $M^{2}$ total polygons partitioned from one
another by at most $M^{4}$ line segments, not including the domain boundaries.
To see this, note that each $n$-sided convex polygon bordered only by other
convex polygons must border at least $n$ other polygons, and thus the number
of sides of each polygon not including the boundaries of $\Omega$ is upper
bounded by the total number of polygons. For any function $A\in\mathcal{A}$,
we can find a function $A^{\prime}\in\mathcal{N}_{\epsilon}$ such that, on
boxes with no intersections by line segments of $A$,
$\|A-A^{\prime}\|_{L^{\infty}}\leq\frac{\eta}{2}$, and on boxes with line
intersections, $\|A-A^{\prime}\|_{L^{\infty}}<|\alpha-\beta|$. Since the
latter error is not controlled, we must bound the number of boxes that may be
intersected by line segments of $A$. Each line segment may pass through at
most $2^{\ell+1}$ boxes, so a set, denoted $B_{i}$, of at most
$M^{4}2^{\ell+1}$ boxes may have errors of $|\alpha-\beta|$. Therefore,

$\displaystyle\|A-A^{\prime}\|^{2}_{L^{2}}$ $\displaystyle\leq\sum_{b\in
B_{\ell}}\|A-A^{\prime}\|^{2}_{L^{\infty}(b)}\left(\frac{1}{2^{2\ell}}\right)$
$\displaystyle\leq\sum_{b\in
B_{i}}|\alpha-\beta|^{2}\left(\frac{1}{2^{2\ell}}\right)+\sum_{b\in B\setminus
B_{i}}\frac{\eta}{2}\left(\frac{1}{2^{2\ell}}\right)$
$\displaystyle\leq\frac{M^{4}2^{\ell+1}}{2^{2\ell}}|\alpha-\beta|^{2}+\frac{\eta}{2}.$

A choice of
$\ell>\log_{2}\left(\frac{M^{4}|\alpha-\beta|^{2}}{\epsilon}\right)+2$ and
$\eta<\epsilon$ gives the result.

∎

######  Lemma C.2.

Let $\mathcal{A}\subset\mathsf{PD}_{\alpha,\beta}$ be the set of
$\gamma$-Lipschitz star-shaped inclusions, defined as functions that takes one
constant value inside a domain $B$, whose boundary is defined by a
$\gamma$-Lipschitz polar function $r=g(\theta)$ contained in the unit square,
and another value outside the domain $B$. The set $\mathcal{A}$ is compact.

###### Proof.

For arbitrary $\epsilon>0$, we identify a finite $\epsilon$-net, denoted
$\mathcal{N}_{\epsilon}$, for $\mathcal{A}$. Partition the unit square by a
uniform grid of $2^{2\ell}$ boxes for integer $\ell>1$, and denote by
$B_{\ell}$ this set of boxes. Let $\mathcal{N}_{\epsilon}$ be the set of
piecewise-constant functions that are constant over each box and take values
only in the set $\\{z:\;z=\alpha+\eta
k,k\in\\{0,1,\dots,\lceil{\frac{\beta-\alpha}{\eta}\rceil}\\}\\}$, for
$\eta<\epsilon$.

We need to bound the number of grid boxes of height and width
$h=\frac{1}{2^{\ell}}$ that can be intercepted by the parametric curve
$C=(g(\theta),\theta)$ in the unit square. First observe that if the curve is
partitioned into segments of length $h$, each segment can intersect at most
$4$ grid boxes. This can be seen by noticing that, although a curve can pass
through four grid boxes that meet at a corner with arbitrarily small arc
length, to reach another one, the curve would have to cross one of the four
boxes entirely, which would entail an arc length longer than $h$. Now note
that the length of the curve is bounded above by

$\displaystyle L$
$\displaystyle=\int_{C}\;\>\mathsf{d}s\leq\int_{0}^{2\pi}\sqrt{r^{2}+\left(\frac{dr}{\>\mathsf{d}\theta}\right)^{2}}\;\>\mathsf{d}\theta\leq\int_{0}^{2\pi}\sqrt{2+\gamma^{2}}\;\>\mathsf{d}\theta\leq
2\pi\sqrt{2+\gamma^{2}}.$

Therefore, the line intercepts at most
$2^{\ell+1}\left(4\pi\sqrt{2+\gamma^{2}}\right)$ boxes. Now, for any function
$A\in\mathcal{A}$, we can find a function
$A^{\prime}\in\mathcal{N}_{\epsilon}$ such that, on boxes with no
intersections by the curve of $A$,
$\|A-A^{\prime}\|_{L^{\infty}}\leq\frac{\eta}{2}$, and on boxes with
intersections, $\|A-A^{\prime}\|_{L^{\infty}}<|\beta-\alpha|$. Therefore, we
have

$\displaystyle\|A-A^{\prime}\|_{L^{2}}^{2}$ $\displaystyle\leq\sum_{b\in
B_{\ell}}\|A-A^{\prime}\|_{L^{\infty}}^{2}\left(\frac{1}{2^{2\ell}}\right)$
$\displaystyle\leq\sum_{b\in
B_{i}}|\alpha-\beta|^{2}\left(\frac{1}{2^{2^{\ell}}}\right)+\sum_{b\in
B\setminus B_{i}}\frac{\eta}{2}\left(\frac{1}{2^{2\ell}}\right)$
$\displaystyle\leq\frac{2^{\ell+1}\left(4\pi\sqrt{2+\gamma^{2}}\right)}{2^{2\ell}}|\alpha-\beta|^{2}+\frac{\eta}{2}.$

Picking
$\ell>\log_{2}\left(\frac{\left(4\pi\sqrt{2+\gamma^{2}}\right)|\alpha-\beta|^{2}}{\epsilon}\right)$
and $\eta<\epsilon$ gives the result. ∎