April 02, 2024

This paper proposes a novel collocation-type numerical stochastic homogenization method for prototypical stochastic homogenization problems with random coefficient fields of small correlation lengths. The presented method is based on a recently introduced localization technique that enforces a super-exponential decay of the basis functions relative to the underlying coarse mesh, resulting in considerable computational savings during the sampling phase. More generally, the collocation-type structure offers a particularly simple and computationally efficient construction in the stochastic setting with minimized communication between the patches where the basis functions of the method are computed. An error analysis that bridges numerical homogenization and the quantitative theory of stochastic homogenization is performed. In a series of numerical experiments, we study the effect of the correlation length and the discretization parameters on the approximation quality of the method.

**Keywords:** numerical homogenization, stochastic homogenization, super-localization, quantitative theory, error estimates, uncertainty

**AMS subject classifications:** 35R60, 65N12, 65N15, 65N35, 73B27

This paper presents a novel numerical stochastic homogenization method for the prototypical random diffusion problem \[-\mathrm{div}(\boldsymbol{A}\nabla \boldsymbol{u}) = f\] subject to homogeneous Dirichlet boundary conditions. Microscopic features of the problem are encapsulated in the random diffusion coefficient \(\boldsymbol{A}\). In this paper, we are particularly interested in coefficients with small correlation lengths.

For deterministic coefficients, numerical homogenization techniques have been studied extensively in the last decades. For recent monographs and reviews on this topic, we refer to [1]–[5]. The random case has not received a similar attention. However, there are several numerical approaches. Let us mention, for example, MsFEM-based numerical stochastic homogenization methods that assume so-called weakly random coefficients; see [6] for an overview. Furthermore, there are many approaches to approximate the effective coefficient in stochastic homogenization, such as [7]–[11]. Closely related to the present work is the numerical stochastic homogenization method proposed and analyzed in [12], [13]. This method is based on the Localized Orthogonal Decomposition (LOD) introduced in [14], [15]; see also [3], [4] for an overview. More precisely, using a reformulation of the LOD for deterministic problems based on a quasi-local discrete integral operator, as discussed in [16], one can derive an effective model of the problem at hand by taking the expectation. This effective model is deterministic, and its solution gives an accurate coarse scale approximation to the expected value of the solution. However, the method is strongly tied to linear finite elements on simplicial meshes (with piecewise constant gradients), which seems to be an artificial limitation and practically unfavorable, especially since structured deterministic and random diffusion coefficients are often based on Cartesian meshes.

In addition to this technical shortcoming, the localization technique underlying the LOD framework has recently been improved, evolving into the Super-Localized Orthogonal Decomposition (SLOD) introduced in [17] (see also [18]–[22]). While the LOD has exponentially decaying basis functions that lead to an exponentially decaying localization error with respect to the diameter of the basis supports relative to the underlying mesh, the localization error of the SLOD actually decays super-exponentially, cf. [17]. This results in smaller local fine-scale problems when computing the basis functions and increased sparsity for a given tolerance.

In this paper, we propose a computationally simple and efficient numerical stochastic homogenization method based on a special collocation-type formulation of the SLOD. This collocation-type formulation leads to a coarse stiffness-type matrix that can be assembled without any communication between the basis functions defined on the patches of the coarse mesh. This allows each patch to be considered separately for sampling, allowing for improved parallelization and a significant speed-up of the method’s assembly process. Furthermore, the favorable localization properties of the SLOD allow for a computationally efficient sampling procedure. In the case of a random diffusion coefficient with a small correlation length and under standard assumptions of quantitative stochastic homogenization, this paper provides an error estimate for the coarse scale approximation of the proposed method, where certain SLOD-specific quantities contribute in an a posteriori manner. The proof of this error estimate is based on the theory of quantitative stochastic homogenization; see, e.g., [23]–[26]. Classical LOD-techniques [3], [14], [15] are used to further evaluate these SLOD-specific quantities, and a worst-case a priori error analysis is conducted for one of them. Several numerical experiments are performed to quantitatively study the effect of the correlation length and other discretization parameters on the accuracy of the approximation.

This manuscript is structured as follows. First, in 2 we state the model problem in weak form. Then, in 3, we introduce the novel numerical stochastic homogenization method. An a posteriori error analysis of the method is performed in 4. 5 derives a worst-case estimate for the quantity appearing in the a posteriori error bound. Practical aspects of implementation are addressed in 6. Finally, 7 provides numerical experiments that underline the theoretical results of this paper.

We consider the model problem \[\label{e:modelstrong} \left\{ \begin{align} -\mathrm{div} (\boldsymbol{A}(\omega)(x) \nabla \boldsymbol{u}(\omega)(x) ) &= f(x),&x\in D\\ \boldsymbol{u}(\omega)(x)&=0,&x\in \partial D \end{align} \right\}\quad\text{for almost all }\omega\in\Omega,\tag{1}\] where \((\Omega,\mathcal{F},\mathbb{P})\) denotes the underlying probability space, \(f \in L^2(D)\) is a deterministic right-hand side and \(D\) is a \(d\)-dimensional bounded Lipschitz polytope with \(d\in\{1,2,3\}\). Without loss of generality, we assume that \(D\) is scaled to unit size. Suppose that \(\boldsymbol{A}\) is a \(\mathbb{R}^{d\times d }\)-valued pointwise symmetric Bochner measurable function, which is uniformly elliptic and bounded, i.e., there exist \(0<\alpha \leq \beta < \infty\) such that for almost all \(\omega \in \Omega\) \[\label{eq:coeffspd} \alpha |\xi|^2 \leq \langle \xi,\boldsymbol{A}(\omega)(x)\xi \rangle \leq \beta |\xi|^2\tag{2}\] holds for all \(\xi \in \mathbb{R}^d\) and almost all \(x \in \mathbb{R}^d\), where \(\langle\cdot,\cdot\rangle\) denotes the Euclidean inner product of \(\mathbb{R}^d\) and \(|\cdot|\) its induced norm.

The weak formulation of the model problem 1 seeks a \(H^1_0(D)\)-valued random field \(\boldsymbol{u}\) such that for almost all \(\omega\in\Omega\) it holds that \[\label{e:modelweak} \boldsymbol{a}_\omega(\boldsymbol{u}(\omega),v):={( \boldsymbol{A}(\omega)\nabla\boldsymbol{u}(\omega)\,,\,\nabla v )}_{L^2(D)}= {( f\,,\,v )}_{L^2(D)}\quad\text{for all }v\in H^1_0(D).\tag{3}\] Here, \({\left( \cdot\,,\,\cdot \right)}_{L^2(D)}\) denotes the inner product on \(L^2(D)\) or \((L^2(D))^d\).

Subsequently, we introduce a shorthand notation for norms and inner products of Bochner spaces. Let \(X\) be a Hilbert space equipped with the inner product \((\cdot,\cdot)_X\). In this case, the Bochner space \(L^2(\Omega;X)\), denoting the space of \(X\)-valued random fields with finite second moments, is also a Hilbert space with the inner product \[(\boldsymbol{v},\boldsymbol{w})_{L^2(\Omega;X)}\mathrel{\vcenter{:}}= {\mathbb{E}}\big[(\boldsymbol{v}(\omega),\boldsymbol{w}(\omega))_{X}\big].\] We write \(\|\cdot\|_{L^2(\Omega;X)}^2 \mathrel{\vcenter{:}}= (\cdot,\cdot)_{L^2(\Omega;X)}\) for the induced norm of this inner product.

Under the given assumptions, the bilinear from \(\boldsymbol{a}_\omega\) depends continuously on \(\boldsymbol{A}\) and, in particular, is measurable as a function of \(\omega\). Hence, the above problem can be reformulated in the Hilbert space \(L^2(\Omega;H^1_0(D))\). The Lax–Milgram theorem then proves its well-posedness, i.e., there exists a unique solution \({\boldsymbol{u}\in L^2(\Omega;H^1_0(D))}\) satisfying \[\label{eq:stabilitysol} \|\nabla\boldsymbol{u}\|_{L^2(\Omega;L^2(D))} \leq \alpha^{-1}C_{\mathrm F}\|f\|_{L^2(D)},\tag{4}\] where \({\| \cdot \|}_{L^2(D)}\) denotes the \(L^2(D)\)-norm and \(C_F\) is the Friedrichs constant of \(D\).

The construction of the novel stochastic homogenization method is based on ideas of the SLOD introduced in [17]. In the deterministic setting, the SLOD identifies an almost local basis of the space obtained by applying the solution operator to \(\mathbb{P}^0(\mathcal{T}_H)\), the space of piecewise constants with respect to some coarse mesh \(\mathcal{T}_H\) of the domain \(D\). This is achieved by identifying piecewise constant right-hand sides supported on patches of the mesh \(\mathcal{T}_H\) such that their responses under the corresponding localized solution operator have a minimal conormal derivative. These localized responses are then used as the basis functions of the SLOD. In the stochastic setting, an adaptation of this approach is required, which involves identifying deterministic local source terms such that the conormal derivative of the localized responses is small in expectation.

In the following we assume that the considered family of meshes \(\{\mathcal{T}_H\}_{H}\) is quasi-uniform and consists of meshes with closed, convex and shape-regular elements. The parameter \(H>0\) specifies the maximal element diameter of the mesh \(\mathcal{T}_H\). We denote by \(\Pi_H\colon L^2(D) \rightarrow \mathbb{P}^0({\mathcal{T}}_H)\) the \(L^2\)-orthogonal projection onto \(\mathbb{P}^0({\mathcal{T}}_H)\). Let us also give a precise definition of the concept of patches with respect to \(\mathcal{T}_H\). The first-order patch \({\mathsf N}(S) = {\mathsf N}^1(S)\) of a union of elements \(S \subset D\) is defined by \[{\mathsf N}^1(S) \mathrel{\vcenter{:}}= \bigcup \{T\in{\mathcal{T}}_H \,\colon\,T\cap S \neq \emptyset \}.\] For any \(\ell = 2,3,4,\dots\), the \(\ell\)-th order patch \({\mathsf N}^\ell(T)\) of \(T\) is then given recursively by \[\label{eq:patch} {\mathsf N}^\ell(T) \mathrel{\vcenter{:}}= {\mathsf N}^1({\mathsf N}^{\ell-1}(T)).\tag{5}\]

The following derivation of the basis functions considers a fixed element \(T\in {\mathcal{T}}_H\) and oversampling parameter \(\ell \in \mathbb{N}\), where we assume that the patch \(D_T \mathrel{\vcenter{:}}= {\mathsf N}^\ell(T)\) does not coincide with the whole domain. We denote the deterministic source term corresponding to \(T\) by \(g_T \in \mathbb{P}^0({\mathcal{T}}_{H,D_T})\) with the submesh \({\mathcal{T}}_{H,D_T} \mathrel{\vcenter{:}}= \{K \in {\mathcal{T}}_H \,\colon\,K \subset D_T\}\). The global response \(\boldsymbol{\varphi}_T\in L^2(\Omega; H^1_0(D))\) to \(g_T\) is defined for almost all \(\omega \in \Omega\) by \[\label{eq:basisglob} \boldsymbol{a}_\omega(\boldsymbol{\varphi}_T(\omega),v) = {( g_T\,,\,v )}_{L^2(D)}\quad\text{for all }v\in H^1_0(D).\tag{6}\] Its localized version \(\boldsymbol{\varphi}_T^\mathrm{loc}\in L^2(\Omega; H^1_0(D_T))\) is for almost all \(\omega \in \Omega\) defined by \[\label{eq:basisloc} \boldsymbol{a}_\omega(\boldsymbol{\varphi}_T^\mathrm{loc}(\omega),v) = {( g_T\,,\,v )}_{L^2(D_T)}\quad\text{for all }v\in H^1_0(D_T).\tag{7}\] From now on, the dependence of stochastic variables is only indicated by a bold symbol for better readability.

To define the conormal derivative of the localized basis function \(\boldsymbol{\varphi}_T^\mathrm{loc}\), we need to introduce some preliminaries on traces and extensions. We denote by \(H^1_\Gamma(D_T)\) the complete subspace of \(H^1(D_T)\) consisting of functions with trace zero at the boundary segment \(\Gamma \mathrel{\vcenter{:}}= \partial D_T \cap \partial D\). Furthermore, let \[\begin{align} {\operatorname{tr}}\colon H^1_\Gamma(D_T)\to X \mathrel{\vcenter{:}}= \operatorname {im} {\operatorname{tr}}\subset H^{1/2}(\partial D_T) \end{align}\] denote the classical trace operator restricted to \(H^1_\Gamma(D_T)\). As an extension operator, we henceforth consider the \(\boldsymbol{A}\)-harmonic extension operator \(\boldsymbol{\operatorname{tr}}^{-1}\colon L^2(\Omega;X) \to L^2(\Omega; H^1_\Gamma(D_T))\) defined as follows: For almost all \(\omega \in \Omega\) and for any given \(\boldsymbol{b}\in L^2(\Omega;X)\), we set \(({\operatorname{tr}}\boldsymbol{\operatorname{tr}}^{-1}\boldsymbol{b})(\omega) = \boldsymbol{b}(\omega)\) and demand that \[\begin{align} \label{eq:trinv} \boldsymbol{a}(\boldsymbol{\operatorname{tr}}^{-1}\boldsymbol{b}, v ) = 0\quad \text{for all } v\in H^1_0(D_T). \end{align}\tag{8}\] The space of locally \(\boldsymbol{A}\)-harmonic functions satisfying homogeneous Dirichlet boundary conditions on \(\Gamma\) can then be defined as \[\begin{align} \boldsymbol{Y}\mathrel{\vcenter{:}}= \boldsymbol{\operatorname{tr}}^{-1}L^2(\Omega;X) \subset L^2(\Omega; H^1_\Gamma(D_T)). \end{align}\] For more details on trace- and extension operators we refer, e.g., to [27].

Combining 7 8 yields for almost all \(\omega \in \Omega\) and all \(\boldsymbol{v}\in L^2(\Omega;H^1_0(D))\) the identity \[\boldsymbol{a}(\boldsymbol{\varphi}_T^\mathrm{loc}, \boldsymbol{v}) = \boldsymbol{a}(\boldsymbol{\varphi}_T^\mathrm{loc}, \boldsymbol{v}- \boldsymbol{\operatorname{tr}}^{-1}{\operatorname{tr}}\boldsymbol{v}) = (g_T, \boldsymbol{v}- \boldsymbol{\operatorname{tr}}^{-1}{\operatorname{tr}}\boldsymbol{v})_{L^2(D_T)},\] where we used that \((\boldsymbol{v}- \boldsymbol{\operatorname{tr}}^{-1}{\operatorname{tr}}\boldsymbol{v})|_{D_T} \in L^2(\Omega;H^1_0(D_T))\). With this identity, the definition of \(\boldsymbol{\varphi}_T\) in 6 , and \(\operatorname{supp} g_T \subset D_T\), it follows that \[\begin{align} \boldsymbol{a}(\boldsymbol{\varphi}_T- \boldsymbol{\varphi}_T^\mathrm{loc}, \boldsymbol{v}) = (g_T, \boldsymbol{v})_{L^2(D_T)} - \boldsymbol{a}(\boldsymbol{\varphi}_T^\mathrm{loc}, \boldsymbol{v}) = (g_T, \boldsymbol{\operatorname{tr}}^{-1}{\operatorname{tr}}\boldsymbol{v})_{L^2(D_T)}. \end{align}\] Taking the expectation, we obtain for any \(\boldsymbol{v}\in L^2(\Omega;H^1_0(D))\) that \[\begin{align} \label{eq:loc-error} {\mathbb{E}}\big[\boldsymbol{a}(\boldsymbol{\varphi}_T- \boldsymbol{\varphi}_T^\mathrm{loc}, \boldsymbol{v}) \big] = {\mathbb{E}}\big[(g_T, \boldsymbol{\operatorname{tr}}^{-1}{\operatorname{tr}}\boldsymbol{v})_{L^2(D_T)}\big] = \left(g_T, \Pi_H\; {\mathbb{E}}[ \boldsymbol{\operatorname{tr}}^{-1}{\operatorname{tr}}\boldsymbol{v}]\right)_{L^2(D_T)}. \end{align}\tag{9}\] As a consequence, the (almost) \(L^2\)-orthogonality of \(g_T\) to the space \({\mathbb{E}}[\boldsymbol{Y}] \subset H^1_\Gamma(D_T)\) leads to a small expected localization error for the basis functions \(\boldsymbol{\varphi}_T^\mathrm{loc}\).

Therefore, we can obtain an optimal choice of \(g_T\) by performing a singular value decomposition (SVD) of the compact operator \((\Pi_{H, D_T} \circ {\mathbb{E}})\vert_{\boldsymbol{Y}}: \boldsymbol{Y}\rightarrow \mathbb{P}^0 ({\mathcal{T}}_{H, D_T})\) restricted to the complete subspace \(\boldsymbol{Y}\). Note that the rank of \((\Pi_{H, D_T} \circ {\mathbb{E}})\vert_{\boldsymbol{Y}}\) is less than or equal to \(N \mathrel{\vcenter{:}}= \# {\mathcal{T}}_{H,D_T}\). Hence, the SVD is given by \[\label{eq:svd} (\Pi_{H, D_T} \circ {\mathbb{E}})\vert_{\boldsymbol{Y}} \; \boldsymbol{v}= \sum_{k=1}^{N} \sigma_k (\boldsymbol{v}, \boldsymbol{w}_k)_{L^2(\Omega; H^1(D_T))} g_k\tag{10}\] with singular values \(\sigma_1 \geq \dots \geq\sigma_N\geq 0\), \(L^2(\Omega; H^1(D_T))\)-orthonormal right singular vectors \(\boldsymbol{w}_1,\dots,\boldsymbol{w}_N\) and \(L^2(D_T)\)-orthonormal left singular vectors \(g_1,\dots, g_N\). The choice \(g_T=g_N\) as the left singular vector corresponding to the smallest singular value \(\sigma_N\) is optimal in the sense that \[g_N \in \underset{g\in\mathbb{P}^0({\mathcal{T}}_{H,D_T})\;:\; \Vert g \Vert_{L^2(D_T)}=1}{\text{argmin}} \quad \underset{\boldsymbol{v}\in \boldsymbol{Y}\; :\; \Vert \boldsymbol{v}\Vert_{L^2(\Omega;H^1(D_T))}=1}{\sup} (g,{\mathbb{E}}[\boldsymbol{v}])_{L^2(D_T)}.\] The corresponding smallest singular value \(\sigma_N\) is a measure of the (quasi-)orthogonality between \(g_T\) and \({\mathbb{E}}[\boldsymbol{Y}]\). We hence define \[\label{eq:32sigma} \sigma_T(H,\varepsilon,\ell) \mathrel{\vcenter{:}}= \sigma_N = \underset{\boldsymbol{v}\in\boldsymbol{Y}\; :\; \Vert \boldsymbol{v}\Vert_{L^2(\Omega;H^1(D_T))}=1}{\sup} (g_N,{\mathbb{E}}[\boldsymbol{v}])_{L^2(D_T)}.\tag{11}\] We emphasize that the practical implementation of the SVD in 10 is difficult due to the stochasticity involved; for a practical implementation based on sampling, see 6. For the error analysis in the following section, we introduce the quantity \[\label{eq:locerrind} \sigma\mathrel{\vcenter{:}}=\sigma(H,\varepsilon,\ell) \mathrel{\vcenter{:}}= \max_{T\in \mathcal{T}_H} \sigma_T(H,\varepsilon,\ell),\tag{12}\] which is an indicator for the overall localization error.

Given that, in expectation, \(\boldsymbol{\varphi}_T^\mathrm{loc}\) closely approximates the response of the global solution operator applied to \(g_T\), it is reasonable to define the approximation of a non-Galerkin, collocation-type numerical stochastic homogenization method by \[\label{eq:stochhommethod} \bar{u}_{H,\ell}\mathrel{\vcenter{:}}= \sum_{T\in{\mathcal{T}}_{H}} c_T \Pi_H\mathbb{E}[\boldsymbol{\varphi}_T^\mathrm{loc}],\tag{13}\] where \((c_T)_{T\in {\mathcal{T}}_{H}}\) are the coefficients of the expansion of \(\Pi_Hf\) in terms of the basis functions \(\{g_T \,\colon\,T \in \mathcal{T}_H\}\). An illustration of the deterministic basis functions \({\mathbb{E}}[\boldsymbol{\varphi}_T^\mathrm{loc}]\) can be found in 1. Error estimates for this method are derived in the following sections.

In this section, we perform an error analysis of the proposed stochastic homogenization method based on results from the theory of quantitative stochastic homogenization. This theory requires structural conditions on the randomness of the coefficient field \(\boldsymbol{A}\). For simplicity, the conditions are formulated for coefficient fields defined on \(\mathbb{R}^d\). Hence, the following assumptions implicitly assume that the coefficient field is defined on the full space \(\mathbb{R}^d\). A random field defined on the bounded domain \(D\) can be obtained by restriction.

*Assumption 1* (Stationarity and decorrelation). Assume that the random coefficient field \(\boldsymbol{A}\) is

*stationary*, i.e., the law of the shifted coefficient field \(\boldsymbol{A}(\omega)(\cdot+x)\) coincides with the law of \(\boldsymbol{A}(\omega)(\cdot)\) for all \(x\in \mathbb{R}^d\),*quantitatively decorrelated*on scales larger than \(\varepsilon\) in the sense of the spectral gap inequality with correlation length \(\varepsilon>0\), i.e., there exists a constant \(\rho>0\) such that for any Fréchet differentiable random variable \(F=F(\boldsymbol{A})\) the estimate \[\begin{align} \label{SpectralGap} \mathbb{E}\big[|F-\mathbb{E}[F]|^2\big] \leq \frac{\varepsilon^d}{\rho} \mathbb{E}\Bigg[\int_{\mathbb{R}^d} \bigg(\fint_{B_\varepsilon(x)} \bigg|\frac{\partial F}{\partial \boldsymbol{A}}(\tilde{x})\bigg| \,\mathrm{d}\tilde{x}\bigg)^2 \,\mathrm{d}x\Bigg] \end{align}\tag{14}\] holds.

We emphasize that these conditions are standard in the theory of quantitative stochastic homogenization; see, e.g., [26].

The error bound presented in this section is an a posteriori bound including the constant \(\sigma\) from 12 and the Riesz stability constant of the local source terms \(\left\{g_T\,\colon\,T\in\mathcal{T}_H\right\}\), which quantifies their linear independence. Both constants can be computed a posteriori as outlined in 6. Additionally, we provide a worst-case a priori upper bound on \(\sigma\) in 5. Note that in a practical implementation, the Riesz stability of the local source terms can be ensured as outlined in 6 or [17].

*Assumption 2* (Riesz stability). The set \(\left\{g_T\,\colon\,T\in\mathcal{T}_H\right\}\) is a Riesz basis of \(\mathbb{P}^0(\mathcal{T}_H)\), i.e., there exists \(C_{\mathrm{rb}}(H, \ell)>0\) such that for all possible choices of \((c_T)_{T\in\mathcal{T}_H}\) it holds \[C^{-1}_{\mathrm{rb}}(H,\ell) \sum_{T\in\mathcal{T}_H}
c^2_T \leq \bigg\| \sum_{T\in\mathcal{T}_H} c_T g_T \bigg\|^2_{L^2(D)},\] where \(C_\mathrm{rb}(H,\ell)\) depends polynomially on \(H^{-1}\) and \(\ell\).

To handle the stochasticity in the error analysis, we need to estimate the variance of the random variables \((\boldsymbol{\varphi}_T^\mathrm{loc}, \mathbb{1}_K)_{L^2(K)}\) for any \(T,K \in \mathcal{T}_H\), where \(\mathbb{1}_K\) denotes the indicator function of the element \(K\). To achieve this, we employ the spectral gap inequality 14 from 1. The following lemma provides a representation of the Fréchet derivative of \((\boldsymbol{\varphi}_T^\mathrm{loc}, \mathbb{1}_K)_{L^2(K)}\), a crucial element for this particular step.

*Lemma 1* (\(L^2\)-representation of Fréchet derivative). Let \(\boldsymbol{v}\in L^2(\Omega; H^1_0(D_T))\) for almost all \(\omega \in \Omega\) be
defined as the weak solution to \[\begin{align}
\label{def95v95help} \left\{ \begin{array}{rll} -\operatorname{div} (\boldsymbol{A}\nabla \boldsymbol{v})&=\mathbb{1}_K &\text{in } D_T,\\ \boldsymbol{v}&= 0 &\text{on } \partial D_T. \end{array}\right.
\end{align}\tag{15}\] The \(L^2\)-representation of the Fréchet derivative of \((\boldsymbol{\varphi}_T^\mathrm{loc}, \mathbb{1}_K)_{L^2(K)}\) is then given by \[\frac{\partial}{\partial \boldsymbol{A}} (\boldsymbol{\varphi}_T^\mathrm{loc}, \mathbb{1}_K)_{L^2(K)} = - \nabla \boldsymbol{\varphi}_T^\mathrm{loc}\otimes \nabla \boldsymbol{v},\] where \(\otimes\colon {\mathbb{R}}^d \times {\mathbb{R}}^d \rightarrow {\mathbb{R}}^{d\times d}\) denotes the outer product.

*Proof.* Let \(\omega \in \Omega\) be fixed. We rewrite the Fréchet derivative as \[\label{eq:prooffrechet} \frac{\partial}{\partial
\boldsymbol{A}} (\boldsymbol{\varphi}_T^\mathrm{loc}, \mathbb{1}_K)_{L^2(K)}(\delta \boldsymbol{A}) = \left(\frac{\partial \boldsymbol{\varphi}_T^\mathrm{loc}}{\partial \boldsymbol{A}}(\delta \boldsymbol{A}),\mathbb{1}_K\right)_{L^2(K)}= \int_{D_T}
\left(\boldsymbol{A}\nabla \boldsymbol{v}\right) \cdot \nabla \frac{\partial \boldsymbol{\varphi}_T^\mathrm{loc}}{\partial \boldsymbol{A}}(\delta \boldsymbol{A}) \,\mathrm{d}x,\tag{16}\] where we tested the weak formulation of 15 with \(\frac{\partial \boldsymbol{\varphi}_T^\mathrm{loc}}{\partial \boldsymbol{A}}(\delta \boldsymbol{A})(\omega)\in H^1_0(D_T)\). To further simplify the expression on the right-hand side,
we differentiate 7 with respect to \(\boldsymbol{A}\) using the product rule. This gives for any \(w \in H^1_0(D_T)\) that \[\begin{align} \int_{D_T} \big(\delta \boldsymbol{A}\nabla \boldsymbol{\varphi}_T^\mathrm{loc}\big)\cdot\nabla w \,\mathrm{d}x+ \int_{D_T} \boldsymbol{A}\frac{\partial \boldsymbol{\varphi}_T^\mathrm{loc}}{\partial
\boldsymbol{A}}(\delta \boldsymbol{A}) \cdot \nabla w \,\mathrm{d}x= 0.
\end{align}\] Using the test function \(w = \boldsymbol{v}(\omega) \in H^1_0(D_T)\) and combining the previous two identities yields the assertion. ◻

Another ingredient in the error analysis is the following regularity result for the localized basis functions. The result is needed to further estimate the term we get after applying the spectral gap inequality. The proof of this result relies on the condition that the patches take the form of \(d\)-dimensional bricks, cf. 5. This condition can be guaranteed, for example, by considering a brick-shaped domain equipped with a Cartesian mesh.

*Remark 1* (Tilde notation). In the following, we will write \(a \lesssim b\) or \(b\gtrsim a\) if it holds that \(a \leq C b\) or \(a \geq C b\), respectively, where \(C>0\) is a constant that may depend on the domain, the shape of the elements and the bounds \(\alpha,\beta\) of \(\boldsymbol{A}\), but is independent of the discretization parameters \(H\), \(\ell\) and the variations of \(\boldsymbol{A}\).

*Lemma 2* (\(L^4\)-regularity of localized basis functions). Let \(\boldsymbol{A}\) be a random coefficient field subject to 1. Then, assuming that the patches \(D_T\) take the form of bricks, the localized basis functions \(\boldsymbol{\varphi}_T^\mathrm{loc}\) satisfy that \[\int_{D_T} \mathbb{E}\bigg[\Big(\fint_{B_\varepsilon(x)} |\nabla \boldsymbol{\varphi}_T^\mathrm{loc}|^2 \,\mathrm{d}\tilde{x}\Big)^{2} \bigg]
\,\mathrm{d}x\lesssim \left(\ell H\right)^{4-d}.\]

*Proof.* In order to apply 5, we need to construct a function \(b_T\) such that the localized basis function
\(\boldsymbol{\varphi}_T^\mathrm{loc}\) is the weak solution to \[-\nabla\cdot (\boldsymbol{A}\nabla \boldsymbol{\varphi}_T^\mathrm{loc}) = \nabla \cdot b_T \quad \text{on} \;\; D_T\]
subject to homogeneous Dirichlet boundary conditions on \(\partial D_T\). To this end, one may choose \(b_T \mathrel{\vcenter{:}}= \nabla r\) for \(r\)
solving the Laplace problem \(\Delta r = g_T\) subject to homogeneous Dirichlet boundary conditions on \(\partial D_T\). With 5 we then obtain that \[\int_{D_T} \mathbb{E}\bigg[\Big(\fint_{B_\varepsilon(x)} |\nabla \boldsymbol{\varphi}_T^\mathrm{loc}|^2 \,\mathrm{d}\tilde{x}\Big)^{2} \bigg]
\,\mathrm{d}x\lesssim \vert D_T\vert^{1-4/q}\Big(\int_{D_T}\vert b_T\vert^q \,\mathrm{d}x\Big)^{4/q}\] for any \(4<q<\infty\). Using standard elliptic regularity on convex domains yields that \[\|\nabla b_T\|_{L^2(D_T)} = \|D^2 r \|_{L^2(D_T)}\lesssim \| g_T \|_{L^2(D_T)}=1,\] since \(g_T\) is \(L^2\)-normalized. Using Friedrichs’ inequality on \(D_T\) for \(r\in H^1_0(D_T)\), we similarly obtain that \[\|b_T\|_{L^2(D_T)}\lesssim \ell H\|g_T\|_{L^2(D_T)} = \ell H.\] Applying the Sobolev embedding (\(q=6\) is the critical exponent for \(d=3\)) and a scaling argument (the embedding constant scales with the diameter of \(D_T\)), we obtain that \[\int_{D_T} \vert b_T \vert^q \,\mathrm{d}x\lesssim (\ell H)^{d-qd/2}\|b_T\|_{L^2(D_T)}^q + (\ell H)^{d+q(2-d)/2}\|\nabla b_T \|_{L^2(D_T)}^q \lesssim (\ell H)^{d+q(2-d)/2}.\] Combining the previous inequalities and setting
\(q = 5\) gives the assertion. ◻

The following theorem encapsulates the main result of this work, giving an a posteriori error bound for the proposed numerical stochastic homogenization method.

*Theorem 2* (A posteriori error bound). Let \(\boldsymbol{A}\) be a random coefficient field subject to 1.
Then, if 2 is satisfied, the solution 13 of the proposed numerical stochastic homogenization method satisfies for any \(f \in L^2(D)\) that \[\|\boldsymbol{u}- \bar{u}_{H,\ell}\|_{L^2(\Omega,L^2(D))}\lesssim H + C_\mathrm{rb}^{1/2}(H,\ell) \ell^{d/2} \big(\sigma(H,\varepsilon,\ell) + \varepsilon^{d/2} \ell^2
H^{(4-d)/2}\big) \|f\|_{L^2(D)}\] with \(C_\mathrm{rb}\) from 2.

*Proof.* For the error analysis, we introduce the function \[\label{eq:dechelpfun} \boldsymbol{u}_{H,\ell}\mathrel{\vcenter{:}}= \sum_{T\in\mathcal{T}_H}c_T
\boldsymbol{\varphi}_T^\mathrm{loc},\tag{17}\] where \((c_T)_{T \in \mathcal{T}_H}\) are the coefficients of the representation of \(\Pi_Hf\) in terms of the local source terms
\(\{g_T\,\colon\,T \in \mathcal{T}_H\}\). Using the triangle inequality, we obtain that \[\begin{align} \begin{aligned} & \|\boldsymbol{u}- \bar{u}_{H,\ell}\|_{L^2(\Omega,L^2(D))} \\ &
\qquad \leq \|\boldsymbol{u}- \Pi_H\boldsymbol{u}\|_{L^2(\Omega,L^2(D))} + \|\Pi_H(\boldsymbol{u}- \boldsymbol{u}_{H,\ell})\|_{L^2(\Omega,L^2(D))} + \|\Pi_H\boldsymbol{u}_{H,\ell}- \bar{u}_{H,\ell}\|_{L^2(\Omega,L^2(D))} \\ &\qquad \eqqcolon
\Xi_1+\Xi_2+\Xi_3. \end{aligned}
\end{align}\] In the subsequent analysis, we will estimate the terms \(\Xi_1\), \(\Xi_2\), and \(\Xi_3\) separately. Prior to this, we mention the
following approximation result for \(\Pi_H\), the \(L^2\)-orthogonal projection onto \(\mathcal{T}_H\)-piecewise constants: It holds that
\[\Vert v- \Pi_H v \Vert_{L^2(T)} \leq \pi^{-1} H \Vert \nabla v \Vert_{L^2(T)}, \quad v\in H^1(T),\; T \in {\mathcal{T}}_H\label{L2-approx-properties};\tag{18}\]
see, e.g., [28], [29]. For the term \(\Xi_1\),
we obtain using the approximation result 18 and the stability estimate 4 that \[\Xi_1^2 = {\mathbb{E}}\big[\|\boldsymbol{u}-
\Pi_H\boldsymbol{u}\|_{L^2(D)}^2\big] \leq \pi^{-2} H^2 {\mathbb{E}}\big[\|\nabla \boldsymbol{u}\|_{L^2(D)}^2\big]\leq \pi^{-2}\alpha^{-2}C_\mathrm{F}^2 H^2 \|f\|_{L^2(D)}^2.\]

For estimating the term \(\Xi_2\), we first apply the \(L^2\)-stability of \(\Pi_H\) and Friedrichs’ inequality. Then, following the lines of the convergence proof of the SLOD in the deterministic setting, cf. [17], we obtain that \[\begin{align} \Xi_2 & \leq C_\mathrm{F}\|\nabla (\boldsymbol{u}- \boldsymbol{u}_{H,\ell})\|_{L^2(\Omega,L^2(D))} \lesssim (H + C_\mathrm{rb}^{1/2}(H,\ell) \ell^{d/2} \sigma(H,\varepsilon,\ell)) \|f\|_{L^2(D)}. \end{align}\]

In order to estimate term \(\Xi_3\), we recall definitions 17 13 and use the Cauchy–Schwarz inequality to obtain that \[\begin{align} \label{eq:xi953} \begin{aligned} \Xi_3^2 & = \sum_{T\in\mathcal{T}_H} c_T \;{\mathbb{E}}\Big[\big(\Pi_H\boldsymbol{\varphi}_T^\mathrm{loc}-\Pi_H{\mathbb{E}}[\boldsymbol{\varphi}_T^\mathrm{loc}], \Pi_H\boldsymbol{u}_{H,\ell}- \bar{u}_{H,\ell}\big)_{L^2(D_T)}\Big] \\ &\leq \sum_{T\in\mathcal{T}_H} |c_T|\|\Pi_H\boldsymbol{\varphi}_T^\mathrm{loc}-\Pi_H{\mathbb{E}}[\boldsymbol{\varphi}_T^\mathrm{loc}]\|_{L^2(\Omega,L^2(D_T))} \|\Pi_H\boldsymbol{u}_{H,\ell}- \bar{u}_{H,\ell}\|_{L^2(\Omega,L^2(D_T))}. \end{aligned} \end{align}\tag{19}\] Algebraic manipulations then yield for the first term of each summand on the right-hand side of the previous inequality that \[\begin{align} \label{eq:referencephi} \begin{aligned} &\|\Pi_H\boldsymbol{\varphi}_T^\mathrm{loc}- \Pi_H{\mathbb{E}}[\boldsymbol{\varphi}_T^\mathrm{loc}]\|_{L^2(\Omega,L^2(D_T))}^2\\ &\qquad = {\mathbb{E}}\bigg[\int_{D_T} \Big(\sum_{K \subset D_T} \big((\boldsymbol{\varphi}_T^\mathrm{loc},\mathbb{1}_K)_{L^2(K)} - {\mathbb{E}}[(\boldsymbol{\varphi}_T^\mathrm{loc}, \mathbb{1}_K)_{L^2(K)}]\big)\vert K\vert^{-1} \mathbb{1}_K \Big)^2\,\mathrm{d}x\bigg] \\ &\qquad = \sum_{K \subset D_T}\vert K\vert^{-1}\; {\mathbb{E}}\Big[\big((\boldsymbol{\varphi}_T^\mathrm{loc},\mathbb{1}_K)_{L^2(K)} - {\mathbb{E}}[(\boldsymbol{\varphi}_T^\mathrm{loc}, \mathbb{1}_K)_{L^2(K)}]\big)^2\Big]. \end{aligned} \end{align}\tag{20}\] Applying the spectral gap inequality 14 and using the \(L^2\)-representation of the Fréchet derivative from 1, we obtain that \[\begin{align} &{\mathbb{E}}\Big[\big((\boldsymbol{\varphi}_T^\mathrm{loc},\mathbb{1}_K)_{L^2(K)} - {\mathbb{E}}[(\boldsymbol{\varphi}_T^\mathrm{loc}, \mathbb{1}_K)_{L^2(K)}]\big)^2\Big] \\[.5ex] &\qquad \lesssim \varepsilon^d \; {\mathbb{E}}\bigg[\int_{D_T}\Big(\fint_{B_\varepsilon(x)} \big\vert \nabla\boldsymbol{\varphi}_T^\mathrm{loc}\otimes\nabla \boldsymbol{v}\big\vert \text{d}\tilde{x} \Big)^2\,\mathrm{d}x\bigg]\\ &\qquad \leq \varepsilon^d \bigg(\int_{D_T}{\mathbb{E}}\bigg[\Big(\fint_{B_\varepsilon(x)}\vert \nabla\boldsymbol{\varphi}_T^\mathrm{loc}\vert^2\text{d}\tilde{x}\Big)^2\bigg]\,\mathrm{d}x\bigg)^{1/2} \bigg(\int_{D_T}{\mathbb{E}}\bigg[\Big(\fint_{B_\varepsilon(x)}\vert \nabla \boldsymbol{v}\vert^2\text{d}\tilde{x}\Big)^2\bigg]\,\mathrm{d}x\bigg)^{1/2}, \end{align}\] where we used the Cauchy–Schwarz inequality. 2 can be employed to bound the first factor on the right-hand side of the preceding inequality. For estimating the second factor, we note that problem 15 for \(\boldsymbol{v}\in L^2(\Omega;H^1_0(D_T))\) has the same structure as problem 7 for the localized basis functions. Consequently, a result analogous to 2 also holds for \(\boldsymbol{v}\), leading to \[\begin{align} \int_{D_T}{\mathbb{E}}\bigg[\Big(\fint_{B_\varepsilon(x)}\vert \nabla \boldsymbol{v}\vert^2\text{d}\tilde{x}\Big)^2\bigg]\,\mathrm{d}x \lesssim (\ell H)^{4-d} \|\mathbb{1}_K\|_{L^2(D_T)}^4. \end{align}\] Inserting the estimates for \(\boldsymbol{\varphi}_T^\mathrm{loc}\) and \(\boldsymbol{v}\), we get that \[\begin{align} {\mathbb{E}}\Big[\big((\boldsymbol{\varphi}_T^\mathrm{loc},\mathbb{1}_K)_{L^2(K)} - {\mathbb{E}}[(\boldsymbol{\varphi}_T^\mathrm{loc}, \mathbb{1}_K)_{L^2(K)}]\big)^2\Big] \lesssim \varepsilon^d (\ell H)^{4-d} \|\mathbb{1}_K\|_{L^2(D_T)}^2 = \varepsilon^d (\ell H)^{4-d} \vert K\vert . \end{align}\] Using this, we continue to estimate 20 as follows \[\label{eq:estpihphiloc} \|\Pi_H\boldsymbol{\varphi}_T^\mathrm{loc}- \Pi_H{\mathbb{E}}[\boldsymbol{\varphi}_T^\mathrm{loc}]\|_{L^2(\Omega,L^2(D_T))}^2 \lesssim \varepsilon^d \ell^4 H^{4-d}.\tag{21}\] Inserting this estimate into 19 , applying the Cauchy-Schwarz inequality, recalling the finite overlap of the patches, and utilizing 2, we finally obtain for \(\Xi_3\) that \[\begin{align} \Xi_3^2 & \lesssim \varepsilon^{d/2} \ell^2 H^{(4-d)/2} \sqrt{\sum_{T\in\mathcal{T}_H} c_T^2} \sqrt{\sum_{T\in\mathcal{T}_H}\|\Pi_H\boldsymbol{u}_{H,\ell}- \bar{u}_{H,\ell}\|_{L^2(\Omega,L^2(D_T))}^2}\\ & \lesssim \varepsilon^{d/2} \ell^{2+d/2} H^{(4-d)/2} C_\mathrm{rb}^{1/2}(H, \ell) \,\Xi_3. \end{align}\] The assertion follows immediately after combining the estimates for \(\Xi_1\), \(\Xi_2\), and \(\Xi_3\). ◻

This section utilizes LOD theory to derive an upper bound for the quantity \(\sigma\) that appears in the error estimate from 2. We further estimate \(C_\mathrm{rb}\) for the choice of LOD basis functions made for the upper bound on \(\sigma\).

We first derive an upper bound for the localization error \(\sigma\) defined in 12 . The bound is based on the lowest-order LOD from [30]–[32], whose construction uses non-negative bubbles \(\{b_T \,\colon\,T \in \mathcal{T}_H\}\); see also [33]. The bubble function \(b_T\in H_0^1(T)\) is chosen such that \(\Pi_Hb_T = \mathbb{1}_T\) and \[\label{eq:bubble} \pi \Vert b_T \Vert_{L^2(T)} \leq H \Vert \nabla b_T \Vert_{L^2(T)} \leq C_\mathrm{b} \sqrt{\vert T \vert}\tag{22}\] holds, where \(C_\mathrm{b}>0\) is independent of \(H\). Recalling the abbreviation \(D_T = \mathsf N^\ell(T)\) for the \(\ell\)-th order patch around \(T\), cf. 5 , we introduce the space of fine-scale functions supported on \(D_T\) by \(\mathcal{W}_{T,\ell} \mathrel{\vcenter{:}}= \{w\in H_0^1(D_T) \,\colon\,\Pi_H\vert_{D_T} w = 0\}\). The LOD basis function corresponding to the element \(T\in\mathcal{T}_H\) is then defined by \[\label{eq:lodbasis} \boldsymbol{\varphi}_{T,\ell}^\mathrm{LOD}\mathrel{\vcenter{:}}= (1-\boldsymbol{\mathcal{C}}_{T,\ell})b_T \in L^2(\Omega; H^1_0(D_T)),\tag{23}\] where \(\boldsymbol{\mathcal{C}}_{T,\ell}b_T \in L^2(\Omega;\mathcal{W}_{T,\ell})\) denotes the fine-scale correction of the bubble \(b_T\), which is defined for almost all \(\omega \in \Omega\) by \[\label{eq:corrc} \boldsymbol{a}(\boldsymbol{\mathcal{C}}_{T,\ell}b_T, w) = \boldsymbol{a}(b_T,w) \quad \text{ for all } w\in\mathcal{W}_{T,\ell}.\tag{24}\]

The bound on \(\sigma\) in the following lemma is based on the observation that the LOD basis function \(\boldsymbol{\varphi}_{T,\ell}^\mathrm{LOD}\) possesses a \(\mathcal{T}_H\)-piecewise constant source term \[\boldsymbol{g}_{T,\ell}^\mathrm{LOD}\mathrel{\vcenter{:}}= - \mathrm{div} \boldsymbol{A}\nabla \boldsymbol{\varphi}_{T,\ell}^\mathrm{LOD}\in L^2(\Omega;\mathbb{P}^0(\mathcal{T}_{H,D_T}));\] see, e.g., [17].

*Lemma 3* (Upper bound on \(\sigma\)). Choosing an \(L^2\)-normalized version of \(g_T \mathrel{\vcenter{:}}=
{\mathbb{E}}[\boldsymbol{g}_{T,\ell}^\mathrm{LOD}]\) in 11 yields the upper bound \[\label{eq:estsigma} \sigma \lesssim \ell^2H^{-1} \exp(-C_\mathrm{d}
\ell) + \ell^4\left(\frac{\varepsilon}{H}\right)^{d/2}\tag{25}\] with \(C_\mathrm{d}>0\) independent of \(H\) and \(\ell\), provided that
\(\varepsilon\) satisfies the smallness assumption \[\label{assumption95eps1} \varepsilon^d \lesssim \ell^{-8}H^d.\tag{26}\]

*Proof.* For all \(\boldsymbol{v}\in \boldsymbol{Y}\subset L^2(\Omega; H^1_\Gamma(D_T))\) it holds that \({\mathbb{E}}[\boldsymbol{\operatorname{tr}}^{-1}{\operatorname{tr}}\boldsymbol{v}] = {\mathbb{E}}[\boldsymbol{v}]\). Hence, by inserting \(g_T =
{\mathbb{E}}[\boldsymbol{g}_{T,\ell}^\mathrm{LOD}]\) into 11 , we obtain that \[\sigma_T(H,\varepsilon,\ell) \leq \frac{1}{\| g_T \|_{L^2(D_T)}}\; \underset{\substack{\boldsymbol{v}\in
L^2(\Omega;H^1_\Gamma(D_T)) \\ \Vert \boldsymbol{v}\Vert_{L^2(\Omega;H^1(D_T))}=1}}{\sup} (g_T,{\mathbb{E}}[\boldsymbol{\operatorname{tr}}^{-1}{\operatorname{tr}}\boldsymbol{v}])_{L^2(D_T)}.\] Note that by dividing by the norm of \(g_T\), we account for the fact that \(g_T\) may not be normalized. We denote by \(\boldsymbol{{\mathcal{A}}}^{-1}_{T, \ell}\colon L^2(\Omega;L^2(D_T))\to
L^2(\Omega;H^1_0(D_T))\) the local solution operator defined on the patch \(D_T\), which satisfies the following stability estimate \[\label{eq:stabilitysolop} \|\nabla \boldsymbol{{\mathcal{A}}}^{-1}_{T, \ell}\boldsymbol{g}\|_{L^2(\Omega,L^2(D_T))} \lesssim \|\boldsymbol{g}\|_{L^2(\Omega,L^2(D_T))}.\tag{27}\] Therefore, we obtain for any \(\boldsymbol{v}\in L^2(\Omega; H^1_\Gamma(D_T))\) that \[\begin{align} (g_T,{\mathbb{E}}[\boldsymbol{\operatorname{tr}}^{-1}{\operatorname{tr}}\boldsymbol{v}])_{L^2(D_T)} &=
{\mathbb{E}}[(\boldsymbol{g}_{T,\ell}^\mathrm{LOD}, \boldsymbol{v})_{L^2(D_T)} - \boldsymbol{a}(\boldsymbol{{\mathcal{A}}}^{-1}_{T, \ell}\boldsymbol{g}_{T,\ell}^\mathrm{LOD},\boldsymbol{v})] \\&\qquad + {\mathbb{E}}[(g_T -
\boldsymbol{g}_{T,\ell}^\mathrm{LOD}, \boldsymbol{v})_{L^2(D_T)} - \boldsymbol{a}(\boldsymbol{{\mathcal{A}}}^{-1}_{T, \ell}(g_T-\boldsymbol{g}_{T,\ell}^\mathrm{LOD}),\boldsymbol{v})]\\ &\eqqcolon \Xi_1 + \Xi_2.
\end{align}\]

To estimate the term \(\Xi_1\), we apply the deterministic result [17] for any \(\omega \in \Omega\) and use the Cauchy–Schwarz inequality to get that \[\Xi_1\lesssim H^{-1} \exp(-C_\mathrm{d} \ell) \|\boldsymbol{v}\|_{L^2(\Omega,H^1(D_T))}\,\|\boldsymbol{g}_{T,\ell}^\mathrm{LOD}\|_{L^2(\Omega,L^2(D_T))},\] where \(C_\mathrm{d}>0\) is independent of \(H\) and \(\ell\). Using the estimate \[\label{eq:norm95glod} \|\boldsymbol{g}_{T,\ell}^\mathrm{LOD}\|_{L^2(\Omega,L^2(D_T))} \lesssim H^{d/2-2},\tag{28}\] which can be derived by taking the expectation of the corresponding deterministic identity from [17], yields that \[\begin{align} \Xi_1 \lesssim H^{d/2-3} \exp(-C_\mathrm{d} \ell)\, \|\boldsymbol{v}\|_{L^2(\Omega,H^1(D_T))}. \end{align}\]

For the term \(\Xi_2\), we obtain using 27 and the Cauchy–Schwarz inequality that \[\begin{align} \Xi_2 \lesssim \|g_T - \boldsymbol{g}_{T,\ell}^\mathrm{LOD}\|_{L^2(\Omega,L^2(D_T))} \|\boldsymbol{v}\|_{L^2(\Omega,H^1(D_T))}. \end{align}\] In order to estimate the first factor on the right-hand side, we proceed similarly as in the proof of 2 to obtain that \[\begin{align} \|g_T - \boldsymbol{g}_{T,\ell}^\mathrm{LOD}\|_{L^2(\Omega,L^2(D_T))}^2 = \sum_{K\subset D_T}\vert K\vert^{-1}\; {\mathbb{E}}\Big[\big((\boldsymbol{g}_{T,\ell}^\mathrm{LOD}, \mathbb{1}_K)_{L^2(K)} - {\mathbb{E}}\left[(\boldsymbol{g}_{T,\ell}^\mathrm{LOD}, \mathbb{1}_K)_{L^2(K)}\right]\big)^2\Big]. \end{align}\] Using the spectral gap inequality 14 , we obtain for each summand that \[\begin{align} \label{eq:rhsspectralgapie} \begin{aligned} &{\mathbb{E}}\Big[\big((\boldsymbol{g}_{T,\ell}^\mathrm{LOD},\mathbb{1}_K)_{L^2(K)} - {\mathbb{E}}[(\boldsymbol{g}_{T,\ell}^\mathrm{LOD}, \mathbb{1}_K)_{L^2(K)}]\big)^2\Big]\\* & \qquad \lesssim \varepsilon^d \; {\mathbb{E}}\bigg[\int_{{\mathbb{R}}^d} \Big(\fint_{B_\varepsilon(x)} \Big|\frac{\partial (\boldsymbol{g}_{T,\ell}^\mathrm{LOD}, \mathbb{1}_K)_{L^2(K)}}{\partial \boldsymbol{A}}(\tilde{x})\Big| \mathrm{d}\tilde{x} \Big)^2 \,\mathrm{d}x\bigg]. \end{aligned} \end{align}\tag{29}\] The \(L^2\)-representation of the Fréchet derivative of \((\boldsymbol{g}_{T,\ell}^\mathrm{LOD},\mathbb{1}_K)_{L^2(K)}\) is derived in 6. It consists of a sum of outer products of the gradients of combinations of \(b_T\), \(b_K\), \(\boldsymbol{\mathcal{C}}_{T,\ell}b_T\) and \(\boldsymbol{\mathcal{C}}_{T,\ell}b_K\). To estimate the summands involving bubble functions, we utilize the property 22 for all \(K \subset D_T\) and derive the estimate \[\begin{align} \label{eq:bubbleest} \int_{D_T} \Big(\fint_{B_\varepsilon(x)} |\nabla b_K|^2 \,\mathrm{d}\tilde{x}\Big)^2 \,\mathrm{d}x\lesssim H^{d-4}. \end{align}\tag{30}\] To proceed with the estimation of 29 , we need to estimate the four terms resulting from the summands of the Fréchet derivative, cf. 6. In the following, we present the estimate for the second term, noting that all other estimates follow analogously. By employing the regularity result from 7 and 30 , we obtain that \[\begin{align} &\mathbb{E} \bigg[\int_{D_T} \Big(\fint_{B_\varepsilon(x)}| \nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_T \otimes \nabla b_K| \,\mathrm{d}\tilde{x}\Big)^2 \,\mathrm{d}x\bigg]\\ & \quad \leq \bigg(\int_{D_T}\mathbb{E} \bigg[ \Big(\fint_{B_\varepsilon(x)} |\nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_T|^2 \,\mathrm{d}\tilde{x}\Big)^2 \bigg] \,\mathrm{d}x\bigg)^{1/2} \bigg(\int_{D_T} \Big(\fint_{B_\varepsilon(x)}|\nabla b_K |^2 \,\mathrm{d}\tilde{x}\Big)^2 \,\mathrm{d}x\bigg)^{1/2}\\ &\quad \lesssim \ell^{2-d/2}H^{d-4}, \end{align}\] where we used the Cauchy–Schwarz inequality. Note that all four terms can be majorized by \(\ell^{4-d}H^{d-4}\), which results from estimating the last summand. The combination of the previous estimates yields that \[\begin{align} \label{eq:32diff95g95glod} \|g_T - \boldsymbol{g}_{T,\ell}^\mathrm{LOD}\|_{L^2(\Omega,L^2(D_T))} \lesssim \bigg(\sum_{K\subset D_T }\vert K\vert^{-1}\; \varepsilon^d \ell^{4-d}H^{d-4}\bigg)^{1/2} \lesssim \varepsilon^{d/2}\ell^2H^{-2}. \end{align}\tag{31}\]

Using the estimate \[\|\boldsymbol{g}_{T,\ell}^\mathrm{LOD}\|_{L^2(\Omega,L^2(D_T))} \gtrsim \ell^{-2}H^{d/2-2},\] which can be derived by taking the expectation of the corresponding deterministic identity from [17], we can derive the following lower bound for the \(L^2\)-norm of \(g_T\) \[\begin{align} \label{eq:32norm95g} \begin{aligned} \| g_T \|_{L^2(D_T)}^2 &= \| g_T\|_{L^2(\Omega,L^2(D_T))}^2 \\ &\geq \frac{1}{2}\|\boldsymbol{g}_{T,\ell}^\mathrm{LOD}\|_{L^2(\Omega,L^2(D_T))}^2 - \|g_T-\boldsymbol{g}_{T,\ell}^\mathrm{LOD}\|_{L^2(\Omega,L^2(D_T))}^2 \\ & \gtrsim \frac{1}{2}\,\ell^{-4}H^{d-4} - \ell^4H^{-4} \varepsilon^{d} \gtrsim \ell^{-4}H^{d-4}. \end{aligned} \end{align}\tag{32}\] Here, we used the reverse triangle inequality, the weighted Young’s inequality for showing that for \(a,b\geq0\) it holds that \(|a-b|^2\geq \tfrac{a^2}{2}-b^2\), as well as the smallness assumption 26 . Finally, combining all estimates leads to \[\begin{align} \sigma_T \lesssim \frac{1}{\| g_T \|_{L^2(D_T)}}\big(H^{d/2-3} \exp(-C_\mathrm{d} \ell)+ \ell^2H^{-2} \varepsilon^{d/2}\big)\lesssim \ell^{2}H^{-1} \exp(-C \ell) + \ell^4\left(\frac{\varepsilon}{H}\right)^{d/2}. \end{align}\] The assertion follows directly when taking the maximum over all \(T \in \mathcal{T}_H\). ◻

Combining this a priori result for \(\sigma\) with 2 yields the error estimate given in the following corollary. The Riesz constant \(C_\mathrm{rb}\) can be computed a posteriori, cf. 6.

*Corollary 1* (Combined error bound). Suppose that the assumptions of 2 3
are fulfilled and that \(\ell\gtrsim |\log H|\) holds. Then, the solution 13 of the proposed numerical stochastic homogenization method satisfies, for any \(f
\in L^2(D)\), that \[\|\boldsymbol{u}- \bar{u}_{H,\ell}\|_{L^2(\Omega,L^2(D))}\lesssim H + C_\mathrm{rb}^{1/2}(H,\ell) \ell^{4+d/2} \left(\frac{\varepsilon}{H}\right)^{d/2}\|f\|_{L^2(D)}.\]

In a next step, we show that the local source terms corresponding to the LOD basis functions 23 are Riesz stable in the sense of 2.

*Lemma 4* (Riesz stability of LOD source terms). Suppose that \(\ell\) is chosen such that \(\ell \gtrsim |\log(H)|\) and that \(\varepsilon\)
satisfies the smallness assumption \[\label{assumption95eps2} \varepsilon^d \lesssim\ell^{-(8+d)}H^{4+d}.\tag{33}\] Then, for the local source terms \(g_T = {\mathbb{E}}[\boldsymbol{g}_{T,\ell}^\mathrm{LOD}]\) it holds for all \((c_T)_{T\in\mathcal{T}_H}\) that \[\label{eq:reisz}
H^4 \sum_{T\in\mathcal{T}_H} c_T^2 \lesssim \bigg\|\sum_{T\in\mathcal{T}_H} c_T \frac{g_T}{\|g_T\|_{L^2(D_T)}}\bigg\|_{L^2(D)}^2.\tag{34}\]

*Proof.* We begin the proof by noting that applying the weighted Young inequality twice gives the elementary estimate \(|a-b-c|^2 \geq \tfrac14|a|^2-|b|^2-|c|^2\) for any \(a,b,c\geq0\). Combining this with the inverse triangle inequality, we obtain that \[\begin{align} &\bigg\|\sum_{T\in\mathcal{T}_H} c_T \frac{g_T}{\|g_T\|_{L^2(D_T)}}\bigg\|_{L^2(D)}^2
=\bigg\|\sum_{T\in\mathcal{T}_H} c_T \frac{g_T}{\|g_T\|_{L^2(D_T)}}\bigg\|_{L^2(\Omega,L^2(D))}^2 \\ &\qquad \geq \frac{1}{4}\,\bigg\|\sum_{T\in\mathcal{T}_H} c_T
\frac{\boldsymbol{g}_{T,\ell}^\mathrm{LOD}}{\|\boldsymbol{g}_{T,\ell}^\mathrm{LOD}\|_{L^2(D_T)}}\bigg\|_{L^2(\Omega,L^2(D))}^2 - \bigg\|\sum_{T\in\mathcal{T}_H} c_T \frac{g_T -
\boldsymbol{g}_{T,\ell}^\mathrm{LOD}}{\|g_T\|_{L^2(D_T)}}\bigg\|_{L^2(\Omega,L^2(D))}^2\\* & \qquad\qquad - \bigg\|\sum_{T\in\mathcal{T}_H} c_T \bigg(\frac{\boldsymbol{g}_{T,\ell}^\mathrm{LOD}}{\|g_T\|_{L^2(D_T)}} -
\frac{\boldsymbol{g}_{T,\ell}^\mathrm{LOD}}{\|\boldsymbol{g}_{T,\ell}^\mathrm{LOD}\|_{L^2(D_T)}}\bigg)\bigg\|_{L^2(\Omega,L^2(D))}^2 \\* & \qquad \eqqcolon \frac{1}{4}\,\Xi_1 - \Xi_2 - \Xi_3.
\end{align}\] For estimating the term \(\Xi_1\) from below, we use the corresponding deterministic result from [17]
and take the expectation which yields that \[\Xi_1 \gtrsim H^4 \sum_{T\in\mathcal{T}_H} c_T^2.\] To estimate the term \(\Xi_2\) from above, we use the finite overlap of the patches \(D_T\) as well as estimates 31 32 to get that \[\begin{align} \Xi_2 &\lesssim \ell^{4+d} H^{4-d} \sum_{T\in\mathcal{T}_H} c_T^2
\, \|g_T - \boldsymbol{g}_{T,\ell}^\mathrm{LOD}\|_{L^2(\Omega,L^2(D_T))}^2 \lesssim \ell^{8+d} \varepsilon^dH^{-d} \sum_{T\in\mathcal{T}_H} c_T^2.
\end{align}\] The estimate for \(\Xi_3\) can be derived similarly using again the finite overlap of the patches \(D_T\), the reverse triangle inequality, 31 32 . We obtain that \[\begin{align} \Xi_3 &=\bigg\|\sum_{T\in\mathcal{T}_H} c_T
\frac{\boldsymbol{g}_{T,\ell}^\mathrm{LOD}(\|\boldsymbol{g}_{T,\ell}^\mathrm{LOD}\|_{L^2(D_T)}- \|g_T\|_{L^2(D_T)})}{\|g_T\|_{L^2(D_T)}\|\boldsymbol{g}_{T,\ell}^\mathrm{LOD}\|_{L^2(D_T)}} \bigg\|_{L^2(\Omega,L^2(D))}^2 \\ & \lesssim \ell^d
\sum_{T\in\mathcal{T}_H} c_T^2 \, {\mathbb{E}}\Bigg[ \frac{\|g_T - \boldsymbol{g}_{T,\ell}^\mathrm{LOD}\|_{L^2(D_T)}^2}{\|g_T\|_{L^2(D_T)}^2} \Bigg] \lesssim \ell^{8+d} \varepsilon^dH^{-d} \sum_{T\in\mathcal{T}_H} c_T^2.
\end{align}\] Combining the previous estimates and using the smallness assumption 33 yields the assertion. ◻

To effectively implement the proposed numerical stochastic homogenization method, it is crucial to employ an efficient sampling strategy for the space \(\boldsymbol{Y}\) and ensure that the local source terms \(\{g_T\,\colon\,T\in \mathcal{T}_H\}\) form a stable basis of \(\mathbb{P}^0(\mathcal{T}_H)\). These aspects will be addressed in the following two subsections.

We consider an arbitrary patch \(D_T\) and denote the number of coarse elements in this patch by \(N \mathrel{\vcenter{:}}= \#{\mathcal{T}}_{H, D_T}\). In a practical implementation, all local infinite-dimensional problems that appear in the derivation of the basis functions must be replaced by finite-dimensional counterparts. To obtain these finite-dimensional counterparts, we perform a discretization using the \(\mathcal{Q}^1\)-finite element method with respect to the fine mesh \(\mathcal{T}_{h,D_T}\) constructed by uniform refinements of \({\mathcal{T}}_{H,D_T}\). The number of elements of \(\mathcal{T}_{h,D_T}\) is denoted by \(n\).

To handle the stochasticity in the definition of \(\boldsymbol{Y}\), our implementation draws \(M\) samples of the random coefficient \(\boldsymbol{A}\) and, for each sample, closely follows the methodology outlined in [17] for the deterministic case. Specifically, we generate a matrix \(\mathbf{S}_i \in \mathbb{R}^{n\times m}\) for \(i=1,\dots,M\), whose columns represent the coordinate vectors of the discrete \(\boldsymbol{A}(\omega_i)\)-harmonic extensions of \(m\in\mathbb{N}\) samples of random boundary data on \(\partial D_T \backslash \partial D\). Then we compute the matrices \(\mathbf{P}_i\in\mathbb{R}^{N\times m}\) by applying the \(L^2\)-orthogonal projection onto the characteristic functions \(\lbrace \mathbb{1}_T\,\colon\,T\in\mathcal{T}_{H,D_T} \rbrace\) column by column to \(\mathbf{S}_i\). Finally, the SVD of the matrix \(\mathbf{X} \mathrel{\vcenter{:}}= [\mathbf{P}_1,\dots,\mathbf{P}_M]\) is computed, yielding coordinate vectors of potential right-hand sides \(g_T\). For details on the practical realization of this SVD, we refer to [17]. Finally, the localized deterministic basis functions are computed as empirical means, again using \(M\) samples of the random coefficient. In the numerical experiments performed in 7, the number of random boundary samples is set to \(m=3N\). For the number of random coefficient samples, we use \(M = 5000\).

Next, we discuss how the stability of the local source terms \(\{g_T \,\colon\,T \in \mathcal{T}_H\}\) can be ensured in a practical implementation. Our implementation achieves stability by an additional optimization step, similar to the one used in [20]. Given the singular values \(\sigma_1\geq\sigma_2\geq\dots \geq \sigma_N\geq 0\) of the matrix \(\mathbf{X}\) associated with the patch \(D_T\), we consider all indices \(1\leq i \leq N\) such that \[\begin{align} \label{eq:choicelocrhs} \frac{\sigma_i}{\sigma_1}\leq \max \Big\{\Big(\frac{\sigma_N}{\sigma_1}\Big)^{1/p},10^{-10}\Big\} \end{align}\tag{35}\] and denote the resulting set of indices by \(\mathcal{I}\). Each index in the set \(\mathcal{I}\) corresponds to a potential candidate for a local source term. For the choice \(p = 1\) only the smallest singular value is considered. Since our optimization problem is meaningful whenever multiple functions are considered, we restrict ourselves to the choices \(p > 1\).

Among these candidate functions, we choose the one that maximizes a weighted \(L^2(D_T)\)-norm under the unit mass constraint. The weighted \(L^2(D_T)\)-norm is defined using a piecewise constant weighting function that is zero in the central element \(T\) and grows polynomially as the distance from the center increases. This enforces a concentration of mass in the center of each patch, resulting in linearly independent local source terms \(\{g_T \,\colon\,T \in \mathcal{T}_H\}\) in practice. More specifically, we introduce the distance function \(\mathrm{dist}(T, K)\) between the elements \(T,K \in \mathcal{T}_H\) as \[\begin{align} \mathrm{dist}(T, K) \mathrel{\vcenter{:}}= H^{-1}|m_K - m_T | \in \mathbb{N}^d, \end{align}\] where \(m_T , m_K \in \mathbb{R}^d\) are the midpoints of the elements \(T\) and \(K\), respectively. The weighting function is then defined for each element \(K\in{\mathcal{T}}_{H, D_T}\) as \[\begin{align} w_T(K) \mathrel{\vcenter{:}}= \big\vert \mathrm{dist}(T,K)\big\vert_\infty^{r} \end{align}\] for a parameter \(r\geq1\), where \(|\cdot|_\infty\) denotes the infinity norm on \(\mathbb{R}^d\). 2 provides an illustration of this weighting function in two spatial dimensions. In our numerical experiments in 7, we use \(p=1.5\) and \(r = 6\).

*Remark 3* (Computation of \(C_\mathrm{rb}\)). Given the local source terms \(\{g_{T_i}\,\colon\,i=1,\dots,\# \mathcal{T}_H\}\), the Riesz stability constant \(C_\mathrm{rb}\) appearing in 1 equals the reciprocal of the smallest eigenvalue of the matrix \(G\in \mathbb{R}^{\#\mathcal{T}_H \times \#\mathcal{T}_H}\) with entries given by \(G_{ij}=(g_{T_i},g_{T_j})_{L^2(D)}\).

*Remark 4* (Uniform Cartesian meshes). Note that in the case of uniform Cartesian meshes, the computational complexity of the method can be significantly reduced when utilizing the stationarity of the coefficient \(\boldsymbol{A}\), cf. 1. In fact, only \(\mathcal{O} (\ell^d)\) reference patches need
to be considered for the computation of the basis functions and local source terms of the method. All other basis functions and local source terms can then be obtained by translation; see, e.g., [34].

The following numerical experiments are intended to demonstrate the effectiveness of the proposed numerical homogenization method. In our implementation, we consider uniform Cartesian meshes of the domain \(D = (0,1)^d\) with \(d\in\{1,2\}\). Note that from now on we use \(H\) to denote the side length of the elements instead of their diameter. For the solution of the local patch problems and the computation of the reference solution \(\boldsymbol{u}_h\) we employ the \(\mathcal{Q}_1\)-finite element method on the fine mesh \({\mathcal{T}}_h\) with \(h=2^{-10}\). We denote by \(\bar{u}_{H,h,\ell}\) the fully discrete numerical approximation to \({\mathbb{E}}[\boldsymbol{u}]\). In the following all expected values are replaced by appropriate empirical means.

The random coefficients \(\boldsymbol{A}\) that are considered in the following numerical experiments are piecewise constant with respect to the uniform Cartesian meshes \(\mathcal{T}_\varepsilon\) with mesh sizes \(\varepsilon \in \lbrace 2^{-5},2^{-6},2^{-7},2^{-8},2^{-9}\rbrace\). These coefficients take independent and identically distributed element values in the interval \([0.1,1]\). We further consider the sequence of coarse meshes \(\mathcal{T}_H\) with mesh sizes \(H\in \lbrace 2^{-3},2^{-4},2^{-5},2^{-6}\rbrace\). Note that we only consider coarse mesh sizes \(H>\varepsilon\) for which the coarse mesh does not resolve the minimal length scale of the random coefficient. We also exclude combinations of \(H\) and \(\ell\) for which a patch coincides with the whole domain \(D\). To calculate the reference solution, we employ \(M=5000\) samples, which is consistent with the number used for the local patch problems. The samples are obtained by a quasi-Monte Carlo sampling strategy in one spatial dimension and a Monte Carlo sampling strategy in two spatial dimensions.

We first examine the behavior of the localization error indicator \(\sigma\) as a function of the coarse mesh size \(H\) and the correlation length \(\varepsilon\). For this, we consider the case \(d = 2\) and utilize the sequences of coarse meshes and correlation lengths mentioned above. 3 visualizes the values of \(\sigma\) for a fixed correlation length \(\varepsilon\) and varying mesh sizes \(H\) (left) and for fixed \(H\) and varying \(\varepsilon\) (right). In both cases one observes a scaling like \(\tfrac{\varepsilon}{H}\), which numerically validates the upper bound for \(\sigma\) from 3 in the case \(d = 2\). Note that the stochastic errors dominate, and consequently, the first term in 25 , which decays exponentially in \(\ell\), is not visible. Plotting \(\sigma\) as a function of \(\ell\) would give a scaling like \(\ell^{-1/2}\).

Next we examine the behavior of the Riesz stability constant \(C_\mathrm{rb}\) of the local SLOD source terms as a function of \(H\). In 4 we observe that \(C_\mathrm{rb}\) scales like \(H^{-4}\), which is consistent with the results for the stochastically averaged LOD source terms proved in 4. Our numerical experiments indicate no dependency of the Riesz stability constant on \(\varepsilon\) or \(\ell\), which is also in line with the findings from 4.

To numerically verify the convergence of the proposed numerical stochastic homogenization method, we consider the source terms \[f_1(x)= 2\pi^2\sin(x),\quad f(x,y)= 2\pi^2\sin(x)\sin(y)\] in one and two spatial dimensions, respectively. 5 6 show the resulting relative \(L^2\)-errors computed using the reference solution \(\boldsymbol{u}_h\). For fixed \(H\) and varying \(\varepsilon\) we observe the rate \(\varepsilon^{d/2}\), which is in agreement with 1. When considering the converse case, we have to distinguish between one and two spatial dimensions. In one dimension, the expected negative power of \(H\) does not manifest itself, and in our numerical experiments the error remains relatively constant with respect to \(H\) (provided the coarse mesh is sufficiently coarse compared to \(\varepsilon\)). In the two-dimensional case, we observe a negative dependence on \(H\), which is much weaker than the \(H^{-2}\) predicted by 1. The error rather seems to scale like \(H^{-1/3}\).

The error analysis of the numerical homogenization method proposed in this paper is based on a regularity theory of quantitative stochastic homogenization. It is established within the framework of an equation on the entire space \({\mathbb{R}}^d\) as outlined in [35], thereby extending earlier results from [36], [37]. In this paper we require its subsequent adaptation to the Dirichlet problem on a bounded box. The results were previously used in [13]. Their proof, which is beyond the scope of this manuscript, is analogous to the full-space case and uses the boundary regularity theory of [38] as well as a classical regularity theory at edges and corners.

*Lemma 5* (Annealed large-scale \(L^p\) regularity). Let \(d \in \lbrace2,3\rbrace\), and let \(\boldsymbol{A}\) be a random coefficient field
subject to 2 and 1. Let \(Q\subset \mathbb{R}^d\) be a box, let \(\boldsymbol{b}\in L^2(\Omega; L^2(Q))\), and let \(\boldsymbol{u}\in L^2(\Omega;H^1_0(Q))\) be a solution to the linear elliptic PDE \[\label{e:rhs} \begin{align} \begin{aligned} -\nabla \cdot (\boldsymbol{A}\nabla \boldsymbol{u})&=\nabla \cdot \boldsymbol{b}\quad &&\text{on }Q,\\ \boldsymbol{u}&\equiv 0 &&\text{on }\partial Q. \end{aligned}
\end{align}\tag{36}\] Then for any \(2\leq p<\infty\) and any \(p<q<\infty\) there holds a regularity estimate of the form \[\begin{align}
\fint_Q \mathbb{E}\Bigg[\bigg(\fint_{B_\varepsilon(x)} |\nabla \boldsymbol{u}|^2 \,d\tilde{x}\bigg)^{p/2} \Bigg] \,\mathrm{d}x \leq C(\lambda,\Lambda,\rho,p,q) \bigg(\fint_Q \mathbb{E}\Bigg[\bigg(\fint_{B_\varepsilon(x)} |\boldsymbol{b}|^2
\,d\tilde{x}\bigg)^{q/2} \Bigg] \,\mathrm{d}x\bigg)^{p/q}.
\end{align}\]

In the following, we present two results used in the proof of 3. The first provides a \(L^2\)-representation of the Fréchet derivative, which is needed to apply the spectral gap inequality.

*Lemma 6* (Fréchet derivative of LOD right-hand sides). The \(L^2\)-representation of the Fréchet derivative of \((\boldsymbol{g}_{T,\ell}^\mathrm{LOD}, \mathbb{1}_K )_{L^2(K)}\)
is given by \[\frac{\partial}{\partial \boldsymbol{A}} (\boldsymbol{g}_{T,\ell}^\mathrm{LOD}, \mathbb{1}_K )_{L^2(K)} = \nabla b_T \otimes \nabla b_K - \nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_T \otimes \nabla b_K - \nabla b_T
\otimes \nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_K + \nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_T \otimes \nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_K.\]

*Proof.* Since \(\boldsymbol{g}_{T,\ell}^\mathrm{LOD}\) is piecewise constant and by the definition of \(\boldsymbol{\varphi}_{T,\ell}^\mathrm{LOD}\), we obtain that \[(\boldsymbol{g}_{T,\ell}^\mathrm{LOD}, \mathbb{1}_K)_{L^2(K)} = (\boldsymbol{g}_{T,\ell}^\mathrm{LOD}, b_K)_{L^2(K)}=\boldsymbol{a}(\boldsymbol{\varphi}_{T,\ell}^\mathrm{LOD}, b_K) =
\boldsymbol{a}((1-\boldsymbol{\mathcal{C}}_{T,\ell})b_T, b_K).\] Hence, the Fréchet derivative of \((\boldsymbol{g}_{T,\ell}^\mathrm{LOD}, \mathbb{1}_K )_{L^2(K)}\) equals \[\frac{\partial}{\partial \boldsymbol{A}} (\boldsymbol{g}_{T,\ell}^\mathrm{LOD}, \mathbb{1}_K )_{L^2(K)} (\delta \boldsymbol{A}) = \frac{\partial}{\partial \boldsymbol{A}} \boldsymbol{a}(b_T, b_K) (\delta \boldsymbol{A}) -
\frac{\partial}{\partial \boldsymbol{A}} \boldsymbol{a}(\boldsymbol{\mathcal{C}}_{T,\ell}b_T, b_K) (\delta \boldsymbol{A}).\] The first term is easily calculated, yielding \[\frac{\partial \boldsymbol{a}(b_T,
b_K)}{\partial \boldsymbol{A}} (\delta \boldsymbol{A}) = \int_{D_T} \delta \boldsymbol{A}\nabla b_T \cdot \nabla b_K \,\mathrm{d}x.\] For the second term, we obtain with the product rule that \[\frac{\partial\boldsymbol{a}(\boldsymbol{\mathcal{C}}_{T,\ell}b_T, b_K)}{\partial \boldsymbol{A}}(\delta \boldsymbol{A}) = \int_{D_T} \delta \boldsymbol{A}\nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_T \cdot \nabla b_K \,\mathrm{d}x+
\int_{D_T}\boldsymbol{A}\nabla \frac{\partial \boldsymbol{\mathcal{C}}_{T,\ell}b_T }{\partial \boldsymbol{A}}(\delta \boldsymbol{A})\cdot \nabla b_K \,\mathrm{d}x.\] Using 24 , the fact that \(\frac{\partial \boldsymbol{\mathcal{C}}_{T,\ell}b_T }{\partial \boldsymbol{A}}(\delta \boldsymbol{A}) \in \mathcal{W}_{T,\ell}\) and the symmetry of \(\boldsymbol{A}\), yields that
\[\int_{D_T}\boldsymbol{A}\frac{\partial \boldsymbol{\mathcal{C}}_{T,\ell}b_T }{\partial \boldsymbol{A}}(\delta \boldsymbol{A})\cdot \nabla b_K \,\mathrm{d}x=
\int_{D_T}\boldsymbol{A}\frac{\partial \boldsymbol{\mathcal{C}}_{T,\ell}b_T }{\partial \boldsymbol{A}}(\delta \boldsymbol{A})\cdot \nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_K \,\mathrm{d}x. \label{eq:defcorrc}\tag{37}\] Furthermore, by
differentiating 24 , we get for any \(w \in \mathcal{W}_{T,\ell}\) that \[\begin{align} \begin{aligned} \int_{D_T} \delta
\boldsymbol{A}\nabla b_T \cdot \nabla w \,\mathrm{d}x= \int_{D_T} \delta \boldsymbol{A}\nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_T \cdot \nabla w\,\mathrm{d}x+ \int_{D_T }\boldsymbol{A}\nabla \frac{\partial \boldsymbol{\mathcal{C}}_{T,\ell}b_T }{\partial
\boldsymbol{A}}(\delta \boldsymbol{A})\cdot \nabla w \,\mathrm{d}x. \end{aligned} \label{eq:corrcderivate}
\end{align}\tag{38}\] Using 37 and 38 for \(w = \boldsymbol{\mathcal{C}}_{T,\ell}b_K\), we obtain for the Fréchet derivative that \[\begin{align} \frac{\partial \boldsymbol{a}(\boldsymbol{\mathcal{C}}_{T,\ell}b_T, b_K)}{\partial {\boldsymbol{A}}}(\delta {\boldsymbol{A}}) &= \int_{D_T} \delta \boldsymbol{A}\nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_T \cdot
\nabla b_K \,\mathrm{d}x+ \int_{D_T}\delta \boldsymbol{A}\nabla b_T \cdot \nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_K \,\mathrm{d}x\\&\quad- \int_{D_T }\delta \boldsymbol{A}\nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_T \cdot \nabla
\boldsymbol{\mathcal{C}}_{T,\ell}b_K \,\mathrm{d}x.
\end{align}\] The \(L^2\)-representation of the Fréchet derivative of \(\boldsymbol{a}(\boldsymbol{\mathcal{C}}_{T,\ell}b_T, b_K)\) is therefore given by \[\begin{align} \frac{\partial}{\partial \boldsymbol{A}} \boldsymbol{a}(\boldsymbol{\mathcal{C}}_{T,\ell}b_T, b_K) = \nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_T \otimes \nabla b_K + \nabla b_T \otimes \nabla
\boldsymbol{\mathcal{C}}_{T,\ell}b_K - \nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_T \otimes \nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_K.
\end{align}\] The combination of the above results yields the assertion. ◻

The following result is needed to estimate the terms appearing after applying the spectral gap inequality in the proof of 3.

*Lemma 7* (\(L^4\)-regularity estimate for LOD correction operators). Let \(\boldsymbol{A}\) be a random coefficient field subject to 1. Then, the correction of the bubble functions \(\boldsymbol{\mathcal{C}}_{T,\ell}b_T\) satisfies the following regularity estimate \[\int_{D_T} \mathbb{E}\Bigg[\bigg(\fint_{B_\varepsilon(x)} |\nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_T |^2 \,\mathrm{d}\tilde{x}\bigg)^{2} \Bigg] \,\mathrm{d}x\lesssim \left(\frac{\ell}{H}\right)^{4-d}.\]

*Proof.* First, let \(\omega \in \Omega\) be arbitrary but fixed. In order to apply 5, we need to establish
the appropriate right-hand side, which results in the equation for \(\boldsymbol{\mathcal{C}}_{T,\ell}b_T\) taking the form as in 36 . Naturally, \(\boldsymbol{\mathcal{C}}_{T,\ell}b_T\) solves, together with the Lagrange multiplier \(\boldsymbol{p}_{T, \ell}\) the following saddle-point problem
\[\begin{align}
\label{eq:saddle-point} \begin{pmatrix} \boldsymbol{{\mathcal{A}}}_{T, \ell}& {\mathcal{B}}^T\\{\mathcal{B}}& 0 \end{pmatrix} \begin{pmatrix} \boldsymbol{\mathcal{C}}_{T,\ell}b_T\\\boldsymbol{p}_{T, \ell} \end{pmatrix} = \begin{pmatrix}
\boldsymbol{{\mathcal{A}}}_{T, \ell}\, b_T\\0 \end{pmatrix}
\end{align}\tag{39}\] with the patch-local operators \(\boldsymbol{{\mathcal{A}}}_{T, \ell}\colon H^1_0(D_T )\rightarrow H^{-1}(D_T), u\mapsto - \nabla \cdot (\boldsymbol{A}\nabla u)\), \({\mathcal{B}}\colon H^1_0(D_T )\rightarrow \mathbb{P}^0({\mathcal{T}}_{H, D_T}),\) \(v\mapsto \Pi_H\vert_{D_T} v\), and its transpose defined by \({\mathcal{B}}^T\colon
\mathbb{P}^0({\mathcal{T}}_{H, D_T})\rightarrow H^{-1}(D_T), p \mapsto \{v\in H_0^1(D_T)\mapsto \int_{D_T}p\,v\,\mathrm{d}x\}\).

It is a direct consequence that \(\boldsymbol{\mathcal{C}}_{T,\ell}b_T\) solves \[\nabla \cdot (\boldsymbol{A}\nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_T) = \nabla \cdot (\boldsymbol{A}\nabla b_T ) + {\mathcal{B}}^T\boldsymbol{p}_{T, \ell},\] which, for some \(\boldsymbol{q}_{T, \ell}\in L^2(D_T)\), can be rewritten as \[\nabla \cdot (\boldsymbol{A}\nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_T) = \nabla \cdot (\boldsymbol{A}\nabla b_T + \boldsymbol{q}_{T, \ell}).\] To see this, we set \(\boldsymbol{q}_{T, \ell}\mathrel{\vcenter{:}}= \nabla v\), where \(v\) solves \(\Delta v = \boldsymbol{p}_{T, \ell}\) with homogeneous Dirichlet boundary conditions in a ball of radius \(C\ell H\), where the constant \(C>0\) is chosen such that the ball contains \(D_T\).

Furthermore, the local LOD source terms satisfy \(\boldsymbol{g}_{T,\ell}^\mathrm{LOD}= \boldsymbol{p}_{T, \ell}\); see [17]. Hence, using 28 and following the proofs of [13] and 2 yields that \[\label{eq:estgradq} \int_{D_T}|\nabla \boldsymbol{q}_{T, \ell}|^2 \,\mathrm{d}x= \int_{D_T}|D^2 v|^2 \,\mathrm{d}x\lesssim \int_{D_T}|\boldsymbol{p}_{T, \ell}|^2 \,\mathrm{d}x\lesssim H^{d-4},\tag{40}\] as well as \(\|\boldsymbol{q}_{T, \ell}\|_{L^2(D_T)}\lesssim \ell H \|\boldsymbol{p}_{T, \ell}\|_{L^2(D_T)}\), leading to \(\int_{D_T}|\boldsymbol{q}_{T, \ell}|^q\,\mathrm{d}x\lesssim \ell^{d+q(2-d)/2}H^{d-q}\). Moreover, using \(\|\nabla b_T\|_{L^\infty} \approx H^{-1}\) we obtain that \(\int_{T}|\boldsymbol{A}\nabla b_T |^q\,\mathrm{d}x\lesssim H^{d-q}\). Therefore, applying 5 for \(Q = D_T\), \(p=4\) and \(\boldsymbol{b}= \boldsymbol{A}\nabla b_T - \boldsymbol{q}_{T,\ell}\) yields that \[\begin{align} &\int_{D_T} \mathbb{E}\Bigg[\bigg(\fint_{B_\varepsilon(x)} |\nabla \boldsymbol{\mathcal{C}}_{T,\ell}b_T |^2 \,\mathrm{d}\tilde{x}\bigg)^2 \Bigg] \,\mathrm{d}x\\ &\quad\lesssim |D_T|^{(q-4)/q} \bigg(\mathbb{E} \Bigg[\int_{D_T } |\boldsymbol{A}\nabla b_T |^{q}\,\mathrm{d}x+ \int_{D_T}|\boldsymbol{q}_{T, \ell}|^{q}\,\mathrm{d}x\Bigg]\bigg)^{4/q}\\ &\quad \lesssim \left(\frac{\ell}{H}\right)^{4-d}, \end{align}\] which is the assertion. ◻

[1]

E. Chung, Y. Efendiev, and T. Y. Hou. *Multiscale model reduction—multiscale finite element methods and their generalizations*, volume 212 of *Applied Mathematical
Sciences*. Springer, Cham, 2023.

[2]

X. Blanc and C. Le Bris. *Homogenization Theory for Multiscale Problems*. Springer, Cham, 2023.

[3]

R. Altmann, P. Henning, and D. Peterseim. Numerical homogenization beyond scale separation. *Acta Numer.*, 30:1–86, 2021.

[4]

A. Målqvist and D. Peterseim. *Numerical Homogenization by Localized Orthogonal Decomposition*. Society for Industrial and Applied Mathematics (SIAM), 2020.

[5]

H. Owhadi and C. Scovel. *Operator-adapted wavelets, fast solvers, and numerical homogenization*, volume 35 of *Cambridge Monographs on Applied and Computational
Mathematics*. Cambridge University Press, Cambridge, 2019.

[6]

A. Anantharaman, R. Costaouec, C. Le Bris, F. Legoll, and F. Thomines. Introduction to numerical stochastic homogenization and the related computational challenges: some recent
developments. In *Multiscale Modeling and Analysis for Materials Simulation*, pages 197–272. World Scientific, 2011.

[7]

A. Bourgeat and A. Piatnitski. Approximations of effective coefficients in stochastic homogenization. *Ann. Inst. Henri Poincaré Probab. Stat.*, 40(2):153–165,
2004.

[8]

É. Cancès, V. Ehrlacher, F. Legoll, and B. Stamm. An embedded corrector problem to approximate the homogenized coefficients of an elliptic equation. *C. r.,
Math.*, 353(9):801–806, 2015.

[9]

J.-C. Mourrat. Efficient methods for the estimation of homogenized coefficients. *Found. Comput. Math.*, 19(2):435–483, 2018.

[10]

J. Fischer. The choice of representative volumes in the approximation of effective properties of random materials. *Arch. Ration. Mech. Anal.*, 234(2):635–726, 2019.

[11]

V. Khoromskaia, B. N. Khoromskij, and F. Otto. Numerical study in stochastic homogenization for elliptic partial differential equations: Convergence rate in the size of representative
volume elements. *Numer. Linear Algebra Appl.*, 27(3), 2020.

[12]

D. Gallistl and D. Peterseim. Numerical stochastic homogenization by quasilocal effective diffusion tensors. *Commun. Math. Sci.*, 17(3):637–651, 2019.

[13]

J. Fischer, D. Gallistl, and D. Peterseim. A priori error analysis of a numerical stochastic homogenization method. *SIAM J. Numer. Anal.*, 59(2):660–674, 2021.

[14]

A. Målqvist and D. Peterseim. Localization of elliptic multiscale problems. *Math. Comp.*, 83(290):2583–2603, 2014.

[15]

P. Henning and D. Peterseim. Oversampling for the multiscale finite element method. *Multiscale Model. Simul.*, 11(4):1149–1175, 2013.

[16]

D. Gallistl and D. Peterseim. Computation of quasi-local effective diffusion tensors and connections to the mathematical theory of homogenization. *Multiscale Model. Simul.*,
15(4):1530–1552, 2017.

[17]

M. Hauck and D. Peterseim. Super-localization of elliptic multiscale problems. *Math. Comp.*, 92(342):981–1003, 2022.

[18]

P. Freese, M. Hauck, and D. Peterseim. Super-localized orthogonal decomposition for high-frequency Helmholtz problems. *Accepted for publication in SIAM J. Sci. Comput.
(arXiv preprint 2112.11368)*, 2024.

[19]

F. Bonizzoni, M. Hauck, and D. Peterseim. A reduced basis super-localized orthogonal decomposition for reaction–convection–diffusion problems. *J. Comput. Phys.*, 499:112698,
2024.

[20]

F. Bonizzoni, P. Freese, and D. Peterseim. Super-localized orthogonal decomposition for convection-dominated diffusion problems. *arXiv preprint 2206.01975*, 2022.

[21]

P. Freese, M. Hauck, T. Keil, and D. Peterseim. A super-localized generalized finite element method. *Numer. Math.*, 156(1):205–235, 2023.

[22]

M. Hauck and A. Målqvist. Super-localization of spatial network models. *Accepted for publication in Numer. Math. (arXiv preprint 2210.07860)*, 2024.

[23]

A. Gloria and F. Otto. . *Ann. of Prob.*, 39(3):779 – 856, 2011.

[24]

A. Gloria and F. Otto. An optimal error estimate in stochastic homogenization of discrete elliptic equations. *Ann. Appl. Probab.*, 22(1), 2012.

[25]

A. Gloria, S. Neukamm, and F. Otto. Quantification of ergodicity in stochastic homogenization: optimal bounds via spectral gap on glauber dynamics. *Invent. Math.*,
199(2):455–515, 2014.

[26]

A. Gloria, S. Neukamm, and F. Otto. A regularity theory for random elliptic operators. *Milan J. Math.*, 88(1):99–170, 2020.

[27]

J.-L. Lions and E. Magenes. *Non-homogeneous boundary value problems and applications. Vol. I*. Die Grundlehren der mathematischen Wissenschaften, Band
181. Springer-Verlag, New York-Heidelberg, 1972.

[28]

L. E. Payne and H. F. Weinberger. An optimal Poincaré inequality for convex domains. *Arch. Rational Mech. Anal.*, 5:286–292 (1960), 1960.

[29]

M. Bebendorf. A note on the Poincaré inequality for convex domains. *Z. Anal. Anwendungen*, 22(4):751–756, 2003.

[30]

R. Maier. A high-order approach to elliptic multiscale problems with general unstructured coefficients. *SIAM J. Numer. Anal.*, 59(2):1067–1089, 2021.

[31]

M. Hauck and D. Peterseim. Multi-resolution localized orthogonal decomposition for Helmholtz problems. *Multiscale Model. Simul.*, 20(2):657–684, 2022.

[32]

Z. Dong, M. Hauck, and R. Maier. An improved high-order method for elliptic multiscale problems. *SIAM J. Numer. Anal.*, 61(4):1918–1937, 2023.

[33]

M. Feischl and D. Peterseim. Sparse compression of expected solution operators. *SIAM J. Numer. Anal.*, 58(6):3144–3164, 2020.

[34]

D. Gallistl and D. Peterseim. Stable multiscale Petrov-Galerkin finite element method for high frequency acoustic scattering. *Comput. Methods Appl. Mech.
Eng.*, 295:1–17, 2015.

[35]

M. Duerinckx and F. Otto. Higher-order pathwise theory of fluctuations in stochastic homogenization. *Stoch. PDE: Anal. Comp.*, 8:625 – 692, 2020.

[36]

S. Armstrong and J.-P. Daniel. Calderón–Zygmund estimates for stochastic homogenization. *J. Funct. Anal.*, 270(1):312–329, 2016.

[37]

M. Duerinckx, A. Gloria, and F. Otto. The structure of fluctuations in stochastic homogenization. *Commun. Math. Phys.*, 377:259 – 306, 2020.

[38]

J. Fischer and C. Raithel. Liouville principles and a large-scale regularity theory for random elliptic operators on the half-space. *SIAM J. Math. Anal.*, 49(1):82–114,
2017.

The work of Moritz Hauck, Hannah Mohr, and Daniel Peterseim is part of a project that has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 865751 – RandomMultiScales). The work of Moritz Hauck is also supported by the Knut and Alice Wallenberg foundation postdoctoral program in mathematics for researchers from outside Sweden (Grant No. KAW 2022.0260).↩︎