{'fraction_non_alphanumeric': 0.05237946880598177, 'fraction_numerical': 0.015811200249240594, 'mean_word_length': 4.301641375038711, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An asymptotically compatible probabilistic collocation method for randomly heterogeneous nonlocal problems

In this paper we present an asymptotically compatible meshfree method for solving nonlocal equations with random coefficients, describing diffusion in heterogeneous media. In particular, the random diffusivity coefficient is described by a finite-dimensional random variable or a truncated combination of random variables with the Karhunen-Lo\`{e}ve decomposition, then a probabilistic collocation method (PCM) with sparse grids is employed to sample the stochastic process. On each sample, the deterministic nonlocal diffusion problem is discretized with an optimization-based meshfree quadrature rule. We present rigorous analysis for the proposed scheme and demonstrate convergence for a number of benchmark problems, showing that it sustains the asymptotic compatibility spatially and achieves an algebraic or sub-exponential convergence rate in the random coefficients space as the number of collocation points grows. Finally, to validate the applicability of this approach we consider a randomly heterogeneous nonlocal problem with a given spatial correlation structure, demonstrating that the proposed PCM approach achieves substantial speed-up compared to conventional Monte Carlo simulations.


Introduction
Since the last decade, there has been a great interest in using nonlocal integro-differential equations to describe physical systems, due to their natural ability to describe physical phenomena at small scales and their reduced regularity requirements which lead to greater flexibility . These nonlocal models are defined in terms of a lengthscale δ, referred to as a horizon, which denotes the extent of nonlocal interaction.
The nonlocal viewpoint allows a natural description of processes requiring reduced regularity in the relevant solution, such as the peridynamics model for fracture mechanics [2,24,25]. An important feature of such models is that they revert back to corresponding classical partial differential equation (PDE) models as the horizon δ → 0. When refining the spatial discretization characterized by grid size h such that h → 0, discretization methods which preserve correct local limits are termed asymptotically compatible (AC) schemes [26], and there has been significant work in recent years toward establishing such discretizations, see, e.g., [26][27][28][29][30][31][32][33][34][35][36]. Broadly, strategies either involve adopting traditional weak form via finite element shape functions and carefully performing geometric calculations to integrate over relevant horizon/element subdomains, or adopting a strong form meshfree discretization where particles are associated with abstract measure [37].
The former is more amenable to mathematical analysis due to a better variational setting, while the latter is simple to implement and generally with smaller computational cost [38,39]. In this work we pursue the asymptotically compatible meshfree approach.
One of the limitations of the current state-of-art works is that most of them consider a homogenized nonlocal model, which may not work well when the material is heterogeneous and its microsctructure plays a critical role. In a recent study [40], Zhao et al. found that a fully homogenized peridynamic model fails to capture certain correct fracture modes/patterns in reinforced concrete. Therefore, they have proposed a stochastic bond-based peridynamic model where the material property is described as random fields. The type of each bond connecting material points x and y was modeled by a random variable, and the discrete probability distribution of this random variable depends on the volume fraction of aggregate and cement on x and y. With this model, fracture patterns and the order in which various cracks develop match experimental observations. Their findings indicate the importance of considering the spatial variability of material properties in nonlocal models, especially when the physical parameters describing spatially varying properties of heterogeneous materials cannot be accurately characterized in all details.
In the present paper, we consider a stochastic nonlocal diffusion equation, where the heterogeneous material property is modeled by a random field. The solution of this stochastic equation describes the probability density function (PDF) of the state variable, e.g., the concentration. Differing from [40] which studied the solution pattern on each individual realization rather than the solution statistics, we focus on the numerical estimates of the first two statistical moments, i.e., the mean and (co)variance. The mean provides an unbiased estimate of the variables and the variance quantifies the uncertainty associated with this estimate.
The Monte Carlo (MC) method and its variations [41][42][43][44][45] are usually used to solve stochastic moment equations and often considered a reliable numerical tool [46]. In the MC method, a large number (denoted by K) of random realizations (samples) are generated for the prescribed random inputs and repetitive deterministic solvers are employed for each sample. The results are then statistically analyzed based on all K samples to calculate leading moments of variables of interest. However, its slow O(K −1/2 ) convergence rate hinders the application of MC method on relatively large scale problems, since one has to solve the differential equation for every sample. Moreover, comparing with local (classical) PDE models, numerically solving nonlocal equations is often more expensive due to its relative lack of sparsity. Therefore, the need for efficient and accurate stochastic numerical methods is even more pressing in the nonlocal setting.
To achieve a faster convergence rate, several stochastic numerical methods were developed for stochastic local (classical) PDE models, including probabilistic Galerkin methods (PGMs) [47][48][49][50][51][52][53], probablistic collocation methods (PCMs) [54][55][56][57][58], reduced basis methods [59][60][61][62][63][64], etc. Among these methods, the probabilistic collocation method with sparse grids inherits the ease of implementation in the MC methods since only solutions at sample points are needed. At the same time, it also reduces the required number of sample points to achieve a given numerical accuracy, especially on problems with small random dimensions and sufficient solution smoothness in the parameter space. Therefore, in this work we will employ the probabilistic collocation methods with sparse grids. Comparing with the attentions received by stochastic PDE problems, numerical studies of nonlocal problems in the uncertainty quantification setting remain limited. In [64], reduced-basis methods are developed for constructing surrogates of the solution of a parameterized nonlocal diffusion problem with random input data in a finite element method framework. However, to the authors' best knowledge, there exists no work on studying the solution smoothness and the theoretical limiting behavior of stochastic nonlocal problems when δ → 0, while these studies are crucial to the design of accurate and asymptotically compatible stochastic numerical schemes. Moreover, the application and rigorous error estimates of meshfree discretization method also remain limited for stochastic nonlocal problems.
The major contribution of the present work is to propose a complete workflow of asymptotically compatible stochastic numerical methods and rigorous mathematical analysis for randomly heterogeneous nonlocal diffusion problems. In particular, we propose to employ a meshfree method with optimization-based quadrature rule [33,37] for the discretization in the physical space, and a PCM with sparse grids for the discretization in the parameter space. By proving that the solution of the nonlocal equation with diffusivity coefficient described as a finite-dimensional random field is analytic in the input random variables, we show that the sparse grid PCM achieves at least algebraic convergence with the increase of sample points. Moreover, given sufficiently large level of sparse grid formulation, the sparse grid PCM converges sub-exponentially.
To characterize the convergence in physical space and the asymptotic compatibility, we for the first time provide analysis for the stochastic nonlocal diffusion problem with random coefficients, showing that its solution converges to the local solution when δ → 0. Based on these mathematical analysis, we provide error estimates for the optimization-based quadrature rule with both fixed horizon δ and also δ going to an asymptotic limit. Lastly, we develop a complete workflow to solve for randomly heterogeneous nonlocal problems, by representing the heterogeneous material coefficient as a random field with given spatial correlation structure and approximating the coefficient by a truncated combination of random variables using the Karhunen-Loève expansion. This work provides a road map to add uncertainty quantification functionality onto pre-existing asymptotically compatible code for deterministic nonlocal problems in a non-intrusive manner, which achieves algebraic or sub-exponential convergence in the solution mean and variance while sustaining the spatial asymptotic compatibility to the correct local limit.
We remark that the paper is organized to establish the rigorous mathematical underpinnings of the approach in the first half, while the second half focuses on a numerical verification and more engineeringoriented exploration of its application. The paper is organized as follows, with all major notations listed in Table 1. We recall first the relevant results in nonlocal calculus and provide mathematical analysis for the deterministic and stochastic nonlocal diffusion problems in Section 2. After establishing the continuous limits of the stochastic nonlocal problem, we next pursue a consistent discretization. In Section 3, we propose our numerical approach for stochastic nonlocal problems by employing the sparse grid PCM and an optimizationbased meshfree quadrature rule in the physical space and establish rigorous error estimates. In particular, we show that the proposed approach achieves algebraic or sub-exponential convergence in the parametric space, and address the convergence rates to the nonlocal and local limits respectively. The theoretical error estimates are verified on a number of one-dimensional and two-dimensional problems with analytic solutions for the local and nonlocal limits in Section 4. In Section 5, we further extend the proposed formulation to handle a more engineering-oriented problem, where the random diffusivity coefficient is modeled by a random field with a given spatial correlation structure. Section 6 summarizes our findings and discusses future research.

The Deterministic and Parametric Nonlocal Diffusion Problem
In this section we introduce the major notations and definitions will be used throughout this paper.
We begin with Table 1 and in Section 2.1 we introduce the deterministic nonlocal diffusion problem while Section 2.2 is dedicated to the stochastic nonlocal diffusion problem. Moreover, we provide novel theoretical analysis for the stochastic nonlocal diffusion problem, namely, its compatibility with the classical companion and analytic regularity.

Nonlocal Calculus and Deterministic Nonlocal Diffusion Problem
In this section, we review the the governing equations of deterministic nonlocal diffusion models which provide the foundation for the stochastic nonlocal problems of interest. Given that Ω ⊂ R d , d ∈ Z + , is a bounded Lipschitz domain, we consider the nonlocal elliptic equation in Ω. To do so, we first introduce the relevant nonlocal calculus. Let α(x, y) : and for a scalar function u(x) : As shown in [24], the adjoint operator of D with respect to the L 2 inner product is D * = −G. We then consider a nonlocal diffusion problem where every point x ∈ Ω is interacting with a neighborhood of points, and their interaction is described by a symmetric kernel function γ(x, y) := α(x, y) · α(x, y) and a two-point scalar function A(x, y) representing the nonlocal diffusion strength. Without loss of generality, we assume in this paper that A is symmetric in its two arguments, i.e., A(x, y) = A(y, x). A nonlocal diffusion operator on a scalar function u : R d → R is then given by In this paper we further assume that the interacting kernel function γ is radial and compactly supported on a Euclidean ball surrounding x, i.e., B δ (x) := {y ∈ R d : |y − x| < δ}: where γ 1 is a nonnegative and nonincreasing function with s-th order singularity, satisfying supp(γ 1 ) ⊂ B 1 (0) and´B 1(0) γ 1 (|z|)|z| 2 dz = d.
The i-th row j-th column element of matrix M .
Symmetric kernel function. s Order of singularity in kernel γ.
a : Ω × Γ → R Local random diffusion strength function. r and R Lower and upper bounds of A and a. S δ (Ω) and T δ Nonlocal energy space and the corresponding bilinear form. ρ : Γ → R + Probability density of the random variable. χ h = {x i } and h (Quasi-uniform) grid set and the grid size for spatial discretization. V h Space of functions to be integrated exactly in the spatial discretization. ω j,i Quadrature weight for x j to generate integral in B δ (x i ).
Prescribed nodes for the Lagrange interpolation in the random space Γ . µ k Corresponding quadrature weight for ξ k in the random space. ζ Sparseness parameter for the Smolyak sparse grid formulation.
Local solution of (2.15). Covariance kernel for the variances between material points x and y. The above kernel assumptions have implications on the boundary conditions that are prescribed on a collar of thickness δ outside the domain Ω, that we denote as and refer to as nonlocal boundary. Therefore, the static nonlocal diffusion problem in a deterministic parameters setting is given as: where u D is the given Dirichlet-type boundary datum in the nonlocal trace space [65]. Without loss of generality, for the analysis, we consider homogeneous Dirichlet boundary conditions u D (x) = 0, and the proposed method is applied to inhomogeneous Dirichlet-type problems in numerical tests of Section 4. Note that although the proposed model can be applied to other boundary conditions, e.g., the Neumann-type boundary conditions in [34,35], here we focus on the Dirichlet-type nonlocal constraint problem for simplicity.
We assume that the diffusion coefficient A ∈ L ∞ ((Ω ∪ Ω δ ) 2 ). To make the nonlocal diffusion problem (2.2) uniformly elliptic, we make the uniform boundedness assumption on the diffusion coefficient A(x, y): where r and R are two positive constants. (2.2) is then associated with the nonlocal energy semi-norm where the nonlocal energy space is defined as Moreover, we define the bilinear form T δ : With the boundedness assumption proposed in (2.3) we have 0 < r|v| . Notice that to solve for the weak solution of (2.2), we find u(x) ∈ S δ (Ω) such that As proved in [66], the bilinear form T δ holds uniform ellipticity as well as the nonlocal Poincaré inequality, and therefore the weak formulation for the deterministic nonlocal equation ( generic constants C and δ 0 > 0 such that for all 0 < δ < δ 0 , the nonlocal energy semi-norm | · | S δ (Ω) satisfies the nonlocal Poincaré inequality As a result of the Poincaré inequality, the semi-norm | · | S δ (Ω) defined is also a norm on S δ (Ω). For the rest of the paper, we will use · S δ (Ω) := | · | S δ (Ω) to denote the norm on S δ (Ω). Combining Lemma 2.1 with the properties of A in (2.3), we can see that the bilinear form T δ is bounded and coercive: Therefore, by the Lax-Milgram theorem, there exists a unique solution u δ ∈ S δ (Ω) for the deterministic . ·, · denotes the duality pairing between (S δ (Ω)) * and S δ (Ω), and f, v = (f, v) L 2 (Ω) when f ∈ L 2 (Ω) ⊂ (S δ (Ω)) * .
Next, we consider the compatibility of the nonlocal diffusion and the classical companion. To properly define the local limit of (2.2) as δ → 0, we need to make the following continuity assumption on the diffusion coefficient A(x, y): A(x, y) ∈ C((Ω ∪ Ω δ ) 2 ) and a(x) := A(x, x).
(2.8) Therefore, the nonlocal diffusion operator L δ [u] has a companion of the classical diffusion operator −∇ · (a(x)∇(u(x))), and (2.2) can be seen as a nonlocal analogue to the local diffusion equation with Dirichlettype boundary condition: From (2.3), we have 0 < r ≤ a(x) ≤ R for any x ∈ Ω, so (2.9) has a unique and bounded solution in When we consider u D = 0, the solution space is denoted by H 1 0 (Ω) := {u ∈ H 1 (Ω), u| ∂Ω = 0}. Notice that we take the minimal regularity assumptions on A(x, y) and a(x). Denoting the solution of local problem (2.9) as u 0 (x) and the solution of the nonlocal one (2.2) with a given horizon size δ as u δ (x), we now show the convergence theorem with the minimal regularity assumptions. Theorem 2.2. Assume that γ(x, y) satisfies the conditions in (2.1), and A(x, y) satisfies the conditions in (2.3). Let f ∈ (S δ (Ω)) * , the dual space of S δ (Ω), then we have In addition, if f (S δ (Ω)) * is uniformly bounded for all δ ∈ (0, δ 0 ) and A(x, y) satisfies (2.8), then the nonlocal and local diffusion problems are compatitble as δ → 0: Proof. We first show the proof of (2.10). Since u δ is a solution to the nonlocal problem, we have . Therefore, we have (2.10).
The proof of the second part involves two steps. In the first step, we assume that Then it is easy to see by Taylor expansion that for any Note that the above equalities are obtained since´B (y − x) i k (k) dy = 0 thanks to the symmetry of the kernel γ δ , where (y − x) (k) denotes the k-th component of (y − x) and i k ≥ 0 is the power on that component. Then we argue that the convergence is also in L 2 (Ω) since |D(AGv)(x)| is uniformly bounded for v ∈ C ∞ 0 (Ω) and δ ∈ (0, δ 0 ). Notice that from the assumption on f (S δ (Ω)) * , we have u δ S δ (Ω) being uniformly bounded for all δ ∈ (0, δ 0 ). Then using similar arguments in [26], we can show u δ − u 0 L 2 (Ω) → 0 as δ → 0. For the general case that A ∈ C((Ω ∪ Ω δ ) 2 ), we will use the mollification technique. Take standard We denote the solution to (2.2) associated with coefficient A to be u δ, . Then we can use the first step to conclude that u δ, − u 0, for any > 0. Therefore, we only need to show For the first equation in (2.11), we first notice that uniformly continuous on (Ω ∪ Ω δ ) 2 . Now since u δ, and u δ are solutions to (2.2) with different coefficients and the same right-hand side, we have where we have used u δ, S δ (Ω) ≤ C from (2.10). Therefore, we have as → 0 and the convergence in L 2 is then implied from the Poincaré inequality. The proof for the second equation in (2.11) is similar by noticing that a − a C(Ω) → 0 as → 0.

Parametric Nonlocal Diffusion Problem
We now consider the case in which the coefficient A is provided by a random field A(x, y, ω), where ω ∈ Ω p and Ω p is the sample space of a probability space (Ω p , F, P). Here, F is the σ-algebra of subsets of Ω p and P is the probability measure. In practice, this random field is usually represented in a "truncated" form using a limited number of random variables (see an example in Section 5.1). Thus, it can be rewritten as A(x, y, ξ), where ξ = (ξ (1) , ξ (2) , . . . , ξ (N ) ), N is a positive integer which denotes the dimension of the parametric space, and ξ (i) are random variables. In practice, we often assume they are independent and identically distributed (i.i.d.) random variables. Under this setting, we consider A(x, y, ξ) : is the space of ξ and it is typically called random space or parametric space. Without loss of generality, here 1], and the random variable ξ ∈ Γ has a probability density ρ : Γ → R + . We are interested in solving the family of nonlocal elliptic equations given by (2.12) For each ξ ∈ Γ , we assume that A(x, y, ξ) is symmetric in its first two variables and A(x, y, ξ) ∈ L ∞ ((Ω ∪ Ω δ ) 2 ). We also assume the uniform ellipticity of the nonlocal problems, i.e., Therefore the Lax-Milgram theorem ensures the well-posedness of nonlocal problem for each ξ ∈ Γ . In addition, in order to consider the limit δ → 0, we need to assume that for each ξ ∈ Γ , (2.14) Then we have the corresponding family of local elliptic equation For each given parameter ξ ∈ Γ , we denote the solution to the nonlocal equation (2.12) by u δ (x, ξ) and the solution to the corresponding local equation (2.15) by u 0 (x, ξ). A corollary of Theorem 2.2 is that Corollary 2.3. Assume that A(x, y, ξ) satisfies (2.13) and (2.14), then we have Proof. For any ξ ∈ Γ , and, we know from Theorem 2.2 that u δ (·, ξ) S δ (Ω) ≤ C for all δ ∈ (0, δ 0 ) and ξ ∈ Γ and δ ∈ (0, δ 0 ). By invoking the dominated convergence theorem, we have To discuss the regularity of solutions with respect to the parameter space, we need to assume the existence of a holomorphic extension of A(x, y, ξ).

Analytic regularity
In order for the function u δ (x, ξ) to be uniformly recovered by polynomial expansions in the parameter space, we will show the analyticity of the solution u δ with respect to the parameterization. By Assumption 2.4, we can extend the definition of A(x, y, ξ) to A(x, y,ξ) whereξ belongs to the complex domain Next, we need the complex uniform ellipticity assumption, namely there exists r c , R c ∈ R such that ( 2.16) forξ ∈Γ . Therefore the nonlocal problem with the complex coefficient A(x, y,ξ) is well-posed and the corresponding solution is denoted by u δ (x,ξ). Here we also remark that for the analyticity of solutions to hold, non-affine coefficients may also be used as long as ξ → A(x, y, ξ) possesses an analytic extension to the complex domains and the complex uniform ellipticity condition is satisfied, see related discussion for the local elliptic equations in [67]. We also define the set So it is clear that the set A containsΓ .
For the discussions in this Section, we need the complex function spaces. We let S δ (Ω, C) be the space of complex valued functions with norm ) and γ δ is still the real valued kernel function. (S δ (Ω, C)) * is the dual space of S δ (Ω, C) equipped with the induced norm. The main result of the this Section is the following theorem on the analyticity of solutions. To prove the theorem, we first need a stability lemma.
Lemma 2.6. Suppose u δ andũ δ are two solutions of (2.2) with the same right-hand side f and different Proof. Subtracting the variational formulations for u δ andũ δ , we find that From the proof of the stability lemma, it is easy to see that similar stability estimate holds when A and Proof of Theorem 2.5. For simplicity of notations, sometimes we only write the explicit dependence of functions on the parameterξ in the discussion, although they may also depend on the spatial variables.
First, it is easy to see that A is an open neighborhood ofΓ by the continuity of the mapsξ → ReA(ξ) andξ → |A(ξ)|. Now we will show that for allξ ∈ A, the functionξ → u δ (ξ) admits a complex derivative Since the mapsξ → ReA(ξ) andξ → |A(ξ)| are continuous, we have the boundedness of ReA(ξ +ĥe i ) and |A(ξ +ĥe i )| for sufficiently smallĥ. Therefore w δ h is well-defined for sufficiently smallĥ. Then for all v ∈ S δ (Ω; C), Therefore w δ h is the unique solution to the variational problem in which the last inequality comes from the stability estimate. Therefore, w δ h converges in S δ (Ω; C) to w δ 0 , which is the solution to¨(

Spatial: Optimization-Based Meshfree Quadrature Rules
In this section we introduce a strong form of particle discretizations of the deterministic nonlocal diffusion problem introduced in Section 2.1. This approach is based upon the optimization-based quadrature rule developed in [33,37]. Denoting the numerical solution of (2.2) as u δ h , two types of convergence are desired in the developed numerical scheme: Algorithm 1 Overall algorithm for the stochastic nonlocal diffusion problem (2.12) 1: Offline Stage: 1a) For spatial discretization, determine a meshfree grid set ∈ Ω ∪ Ω δ and calculate the optimization-based quadrature weights {ω j,i } j for each x i ∈ χ h by solving (3.4).
1b) For probabilistic collocation method, determine the (sparse) collocation points set in the parametric space Θ N = {ξ k } K k=1 ∈ Γ and the corresponding quadrature weights µ k following (3.19). 2: Online Stage: Solving for (2.12): for the deterministic problem corresponding to the k-th sample, by solving u = Q −1 f . 3: Postprocessing Stage: Generate statistical moments of the random solution following (3.14)-(3.15).
The first type of convergence indicates that the numerical discretization method is consistent with the nonlocal problem, while the second type shows that the nonlocal numerical solution preserves the correct local limit, or equivalently, the numerical scheme is asymptotically compatible. To maintain an easily scalable implementation, in asymptotic compatibility studies we assume δ to be chosen such that the ratio δ h is bound by a constant as δ → 0, restricting ourselves to the "δ-convergence" scenario [69] 1 . In this setting, one obtains banded stiffness matrices allowing scalable implementations. Typically in the literature a scheme is termed as asymptotically compatible (AC) if it recovers the classical solution for any δ, h → 0. Here we abuse the definition slightly and only require the δ-convergence. We will provide truncation error estimates for the quadrature error convergence rates to the nonlocal analytical solution and the local limit, respectively.

Mathematical Formulation and Implementation
Discretizing the whole interaction region Ω ∪ Ω δ by a collection of points we aim to solve for the solution u (i) ≈ u(x i ) on all x i ∈ χ h . Although the method can be applied to more general grids, for analysis we require χ h to be a uniform Cartesian grid: Here h is the spatial grid size. For the deterministic nonlocal diffusion model (2.2) we pursue a discretization through the following one point quadrature rule at χ h [70]: where we specify {ω j,i } as a to-be-determined collection of quadrature weights admitting interpretation as a measure associated with each collocation point x i . Note that although we only solve for u (i) on grid points in χ h , we will denote the numerical solution u δ h of the above nonlocal problem as the piecewise constant interpolation of u (i) .
We use an optimization-based approach to define these weights extending previous work [33,37], constructed to ensure consistency guarantees. Specifically, we seek quadrature weights for integrals supported on balls of the form where we include the subscript i in {ω j,i } to denote that we seek a different family of quadrature weights for different subdomains B δ (x i ). Denoting P m (R d ) as the space of m-th order polynomials, we obtain these weights from the following optimization problem where V h = S γ δ ,x := q = p(y − x)γ δ (|x − y|) | p(x) ∈ P 3 (R d ) denotes the space of functions which should be integrated exactly. Note that when p(x) = const, we have q(x, y) = γ δ (|x − y|). Therefore, reproducing this function requires the kernel function γ δ to be integrable, or equivalently, s < d. W (r) is a radially symmetric positive weight function supported in B δ (0). Following the discussions in [28], we take W (r) = γ δ (r).
As shown in [33], this particular choice of reproducing space provides the minimal reproducing set to achieve the optimal O(δ 2 ) asymptotic convergence rate in nonlocal problems with homogeneous diffusion coefficient.
Moreover, we note that provided the quadrature points are unisolvent over the desired reproducing space, problem [28]. For certain choices of V h , such as m-th order polynomials, unisolvency holds under the following assumptions: the domain Ω satisfies a cone condition, the pointset χ h ∩ B δ (x i ) is quasi-uniform, and h/δ is sufficiently small [71].
For each x i ∈ χ h ∩ Ω, we denote the set of its interacting neighbor points as {x j } j=1,··· ,Mi = χ h ∩ B δ (x i )\{x i } and the corresponding quadrature weights in a size M i vector, denoted as ω = [ω 1,i , · · · , ω Mi,i ] T .
For problems where the reproducing constraints are redundant, [HW −1 H ] −1 may be replaced by the pseudo-inverse.

Stability and Convergence Analysis
To provide the stability proof of our method, we first show that the quadrature weights are all positive:

7)
where C is a constant independent of h.
Proof. For a fixed δ, we denote the weights as ω δ j,i . Due to the scaling of the nonlocal kernel in (2.1), ω δ j,i has the following scaling ω δ j,i = δ d ω 1 j,i , and therefore it suffices to prove the estimate for δ = 1. With the symmetry property of γ 1 and the uniform Cartesian grid, we note that g (β) = 0 if one of the components of β is an odd number, and m (β,θ) = 0 if one of the components of β + θ is an odd number. Moreover, ifβ is a reordering of β, and m (β,θ) = m (β,θ) ifβ +θ is a reordering of β + θ. Denote b = [HW −1 H ] −1 g, we then have b (β) = 0 unless |β| l1 = 0 or |β| l1 = 2 with only one nonzero entry β k = 2 for some k ∈ {1, · · · , d}. With the symmetry properties of g (β) and m (β,θ) , we notice that b (β) = b (β) ifβ is a reordering of β. Therefore, we may denote b (0,··· ,0) := b 0 , b (2,0,··· ,0) = b (0,2,··· ,0) = · · · = b (0,··· ,0,2) := b 2 , and b 0 , b 2 satisfy: Notice that s < d − 1 and compare the quadrature rule with the Riemann sum estimates for integrals, we where B √ dh (0) is a sphere covering the origin point and the singularity on it. Note here we don't have an associated weight ω i,i , and the errors of integrals in the (hyper)cube containing the origin point with size h d were estimated separately. B √ dh is then chosen as a sphere to contain this (hyper)cube. We obtain From the above lemma and [71], we can see that there exists a constant C independent of δ, such that for all h < Cδ unisolvency holds for the optimization problem (3.4) and we have ω j,i > 0. In the following we will denote this constant as C pos and generally require h < C pos δ. Note that the symmetry property in Lemma 3.1 also yields ω j,i = ω i,j . Consider A satisfying the conditions in (2.3) and (2.8) and γ δ satisfying the conditions in Lemma 3.1, we obtain the discrete maximum principle for L δ h and the well-poseness property of the discretized nonlocal diffusion problem (3.2): Proof. We first prove the maximum principle for L δ h v h ≥ 0. Assume that max Apply the same argument to x j in Ω and to their neighbors, one may obtain u h = const on χ h and therefore max The discrete minimum principle can be similarly proved.
With the discrete maximum principle, the discretized nonlocal diffusion problem (3.2) is therefore wellposed, i.e., there exists a unique solution to the discretized nonlocal diffusion problem (3.2). We now consider the accuracy of the quadrature rule. In the following, we first present the truncation error estimate of the meshfree discretization for the nonlocal diffusion problem: Proof. We note that On the other hand, let G(x j ) denote the (hyper)cubic of size h d centered at x j , we have the error estimate of the Riemann sum formulation as: where |∇ y a(x, y)| denotes the maximum component of the first order partial derivatives with respect to y.

Therefore, substituting (3.7) into (3.2) yields
Note here the constant C is independent of h but may depends on δ.
To prove the asymptotic compatibility, we need the truncation error estimate to the local limit. In [33, Theorem 2.1], the authors have shown that for a sufficiently smooth u, the quadrature weights obtained from (3.4) with the choice of V h = S γ δ ,x provides an O(δ 2 ) pointwise error bound for the integral approximation of 2´B δ (xi) γ δ (|x i − y|)(u(y) − u(x i )) dy when the ratio h/δ is fixed. We can easily extend this error estimate to a nonlocal diffusion operator with heterogeneous diffusion coefficient: Consider a kernel γ δ satisfying the conditions in Lemma 3.1, A(x, y) and u(y) are C 4 with respect to y, and fixed ratio h/δ < C pos . Quadrature weights obtained from (3.4) with the choice of V h = S γ δ ,x satisfy the following pointwise error estimate, with C > 0 independent of δ: Proof. Note that for a given x we consider the estimate of A(x, y)u(x) ∈ C 4 and A(x, y)u(y) ∈ C 4 , the following truncation estimate is obtained immediately from [33, Theorem 2.1]: And the proof is finished by taking γ δ (|x − y|) = D0 δ d+2−s |xi−y| s as stated in (2.1).
From the derivation of Theorem 2.2, one can see that . We therefore obtain the following truncation error estimate to the local limit: With the discrete maximum principle and the above truncation estimate results, we finally get the main results on the stability in maximum norm: Assume that a ∈ C ∞ (Ω) and Ω ∪ Ω δ ∈ C 1 , then there exist generic constants C and δ 0 > 0 such that when δ < δ 0 and h < C pos δ, solution to the discretized nonlocal diffusion problem (3.2) satisfies: Here C is independent of both h and δ.
Proof. We first construct a barrier function ψ(x). With the properties of a, there exists a C ∞ solutionψ for the following classical diffusion problem [72]: where ∂(Ω ∪ Ω δ ) denotes the exterior boundary of Ω ∪ Ω δ . We then define ψ(x) :=ψ(x) − min which together with (3.9) finishes the proof.
With the stability property of L δ h and the truncation estimates in Lemmas 3.3-3.5 we proceed to prove the two types of convergence results in (3.1). We first consider the case with fixed δ and vanishing h. In particular, we investigate the convergence of numerical solution to nonlocal solution by combining Lemma

with the stability property:
Theorem 3.7 (Convergence to Deterministic Nonlocal Solution). Assume that the conditions in Lemma 3.3 and Lemma 3.6 are satisfied and u δ ∈ C 1 (Ω), then there exists a δ 0 > 0 such that for a fixed δ satisfying 0 < δ < δ 0 and h < C pos δ, the following convergence property holds for the numerical solution of (3.2): where C is a generic constant independent of h.
Proof. Apply the stability theorem to u δ h − u δ , we immediately obtain Next, we show the AC property of the meshfree method, when both δ and h vanish with a fixed ratio δ/h. In particular, we investigate the convergence rate of numerical solution to local limit by combining Lemma 3.4 and Corollary 3.5 with the stability property of L δ h : Theorem 3.8 (Asymptotic Compatibility in Deterministic Nonlocal Problems). Assume that the conditions in Corollary 3.5 and Lemma 3.6 are satisfied, and u δ , u 0 ∈ C 4 (Ω), then there exists a δ 0 > 0 such that for any 0 < δ < δ 0 and fixed ratio h/δ < C pos , the following error estimate holds for u δ h : Moreover, the meshfree scheme (3.2) is asymptotically compatible, i.e., where C is a generic constant independent of δ.
Proof. Apply the stability theorem to u δ h − u δ and u δ h − u 0 , we immediately obtain: Remark 3.9. When taking a flat kernel with s = 0, from Theorems 3.7 and 3.8, we note that the optimal convergence rate to the nonlocal solution is O(h) when h → 0, while the optimal convergence to the local limit is O(δ 2 ) when h, δ → 0.

Stochastic: Probabilistic Collocation Method with Sparse Grids
To solve the stochastic problem introduced in Section 2.2, we employ the probabilistic collocation method (PCM) in the parametric space for its high resolution and ease of implementation by sampling at discrete points in random space [54,73,74]. Consider the stochastic equation (2.12), PCM can be seen as a Lagrange interpolation in the random space. In particular, let Θ N = {ξ k } K k=1 ⊂ Γ be a set of prescribed nodes such that the Lagrange interpolation in the random space Γ is poised in an interpolation space Γ I , where N is the dimension of the parametric space. Then any function v : Γ → R can be approximated using the Lagrange interpolation polynomial: where J k (ξ) is the Lagrange polynomial satisfying J k (ξ) ∈ Γ I and J k (ξ j ) = δ kj . Denotingû(x, ξ) := K k=1 u(x, ξ k )J k (ξ), the collocation procedure to solve the stochastic nonlocal equation is where R is the residual of (2.12). With the property of Lagrange interpolation, we obtain for k = 1, · · · , K. Note that (3.13) is equivalent to solving K deterministic nonlocal diffusion problems, where the deterministic meshfree solver discussed in Section 3.1 can be readily applied. Therefore, the PCM approach can be implemented in an embarrassingly parallel way and the total computational cost is the product of the number of collocation points times the cost of the deterministic problem.
With the numerical solution of (3.13) on all collocation points ξ i , the statistical moments of the random solution can be evaluated: and so on. Here ρ is the PDF of random variable ξ. To further approximate the integral for above polynomials, we employ the quadrature rule approximation by choosing the set Θ N as quadrature point set: where {µ k } K k=1 is the set of corresponding quadrature weights. There are mainly two different strategies for the selection of collocation point sets: the tensor products of 1D collocation point sets and a sparse grid strategy for high dimensionality. In the tensor product strategy, one first construct a 1D interpolation for each dimension in the random space. For the i-th dimension, we take (i) numbers of nodal points Θ (i) , a 1D interpolation for a smooth function v on the i-th dimension then writes: where J i k (ξ (i) ) is the 1D Lagrange polynomial. Then for the multivariate case v : R N → R, the tensor product formula is: Notice here (3.17) requires K = Π N i=1 (i) numbers of collocation points in total, which grows quickly when N gets large. Therefore, the tensor product strategy may be employed for problems with a small number of random dimensions, but its required number of collocation points K generally grows exponentially as N increases and makes the simulation non-feasible (see, e.g., [58]). Hence for problems with a relatively large random dimension, we employ the sparse grid strategy. In particular, we employ the sparse grid constructed by the Smolyak algorithm [75], which is a linear combination of tensor product formulas: Here ζ is the sparseness parameter, = ( (1) , · · · , (N ) ) ∈ N N , | | l1 = N i=1 (i) , and (i) represents the number of collocation points in random dimension i. To compute (3.18), only evaluations on the sparse grids are needed: (

3.19)
As shown in [76,77], (3.18) is exact for p(ξ) ∈ P ζ−N (R N ) (all polynomials of degree less than ζ − N ) and the total number of nodes K ∼ 2N ζ−N (ζ−N )! . Therefore, we may see that the sparse grid formulation typically requires a much smaller number of collocation points K than the full tensor product set and we will refer η = ζ − N as the "level" of the Smolyak formulation.

An Asymptotically Compatible Meshfree PCM
Let δ be the horizon size, h be the grid size in space, and K = K(η, N ) be the total number of collocation points we use in the parameter space. We then denote the numerical solution to (2.12) by u δ h,K . As h → 0 and K → ∞, we expect the numerical solution to converge to the exact solution u δ for (2.12) in L 2 (Ω) ⊗ L 2 ρ (Γ ). In addition, our numerical method is asymptotically compatible, i.e., u δ h,K → u 0 as h → 0, K → ∞, and δ → 0, where u 0 is the exact solution for (2.15). The error estimate can be facilitated by introducing an intermediate function u δ K , the semi-discrete solution to (2.12). Then we split the errors into For each ξ ∈ Γ , The estimates of u δ (ξ k ) − u δ h (ξ k ) L 2 (Ω) for each ξ k are followed by the estimates in Section 3. So we have ρ (Γ ) converges to zero follows the estimates in Section 3. Moreover, we have u 0 (ξ) − u δ (ξ) L 2 (Ω) ≤ Cδ 2 for any ξ as long as u 0 ∈ C 4 (Ω). Therefore we only need to estimate u δ − u δ K . From the analytic regularity, ξ → u δ (ξ) is a S δ (Ω)-valued map that admits an analytic extension to an open neighborhood A ⊂ C N of Γ . Moreover, we know that maxξ ∈A u δ (ξ) S δ (Ω) is uniformly bounded and independent of δ. Therefore, the estimate of u δ − u δ K is followed exactly from [55]. Here we present the following result, whose proof can be found in [55, Theorems 3.10-3.11]. Lemma 3.10. There exists C 1 > 0 and β 1 > 0 depending on N and the analytic region A such that Moreover, when η > N log (2) , there exists C 2 > 0, C 3 > 0 and β 2 > 0 depending on N and the analytic region A, and β 3 > 0 depending only on N such that Equation (3.20) shows at least algebraic convergence with respect to the number of collocation points K, while equation (3.21) shows the subexponential convergence when the level of Smolyak formulation η > N log (2) . In the numerical experiments below, we will choose large enough η so as to observe the subexponential convergence. Based on Lemma 3.10 and the above discussions, we have the following convergence theorems.
Proof. Let (u δ − u δ K ) + (u δ K − u δ h,K ), then the theorem is a combination of Theorem 3.7 and Lemma 3.10.

Numerical Verification of Convergences
In this section, we numerically verify the proposed approach by investigating the two types of convergences: the consistency to nonlocal solutions and then the asymptotic compatibility to local companions 2 .
In particular, we study the L 2 errors for the first two statistical moments, the mean and standard deviation (std). Let u δ h,K represent the numerical solution with spatial grid size h in meshfree methods and K samples in sparse grid PCM, u δ represents the analytical nonlocal solution and u 0 stands for the analytical local limit, in Section 4.1 we investigate the convergence of numerical solutions to the nonlocal analytical solution with vanishing h and increasing sample numbers K by calculating: where the mean and standard deviation of u δ h,K are numerically evaluated with the quadrature rule approximation in (3.14)- (3.15). Similarly, in Section 4.2 we investigate the convergence of numerical solutions to the local limit as K increases and δ, h goes to 0 simultaneously under the δ-convergence limit. In particular we calculate:

Consistency with the Stochastic Nonlocal Solution
In this section we study the consistency of the numerical solution to the nonlocal analytical solution, on both 1D and 2D physical domains.
2 Notice that although we assumed Γ = N i=1 [−1, 1] ⊂ R N in Sections 2-3, in numerical tests we investigate and numerically verify our analysis on more general cases.  Three types of distributions are considered for ξ (i) , i = 1, · · · , 5: the uniform distribution ξ  Figures 1(b) and 1(c), we employ a fixed grid size h = 1/1280 and demonstrate the convergence of solution error with increasing sparse grid level η = 1, · · · , 5 in the parametric space. When η ≥ 3, the numerical error of mean for uniform distribution reaches a convergence plateau because of the spatial discretization error. Noting that we have η < N/ log(2) ≈ 7.5 in this case, algebraic convergence of the sparse grid PCM is verified for all three types of distributions. Therefore, the O(h) spatial convergence and the algebraic convergence of the sparse grid PCM together verify the estimates in Theorem 3.11.  With fixed δ = 0.525 and nonlocal diffusion coefficient A(x, y, ξ) = (2+ξ)[2+cos(A(x (1) +y (1) )) cos(B(x (2) + y (2) ))]/δ 4 , we consider the manufactured nonlocal analytical solution (2) , ξ) = cos(Ax (1) ) sin(Bx (2) )/(2 + ξ). we consider a fixed grid size 1/32 and study the convergence of solution error with increasing sample numbers in PCM. In particular, we take sparse grid levels η ∈ {1, · · · , 5}. When η ≥ 4, the numerical errors of mean for all distribution types reach a convergence plateau because of the spatial discretization error. When η > N/ log(2) ≈ 1.5, sub-exponential convergence is observed for all four types of considered distributions before reaching this plateau, which again verifies Theorem 3.11.

Asymptotic Compatibility (AC) to the Stochastic Local Limit
In this section we investigate the asymptotic compatibility (AC) of the proposed approach by studying the convergence of its numerical solution to the corresponding local limit when δ, h → 0. In particular, we focus on the δ-convergence limit and fix the ratio between δ and h. Following the conventions in [37,78,79], in all tests we consider the nonlocal diffusion coefficient as the harmonic mean of the local diffusion coefficient: A(x, y, ξ) = 2/(a −1 (x, ξ) + a −1 (y, ξ)).
Three types of distributions are considered for ξ (i) , i = 1, · · · , 5: the uniform distribution ξ proposed approach is therefore AC and the observed convergence rate is consistent with our AC analysis in Theorem 3.12. In Figures 3(b) and 3(c) we fix h = 1/5000 and δ = 3.8h, and show the convergence of solution error with increasing sparse grid level η = 1, · · · , 5 in the parametric space. Noting that we have η < N/ log(2) ≈ 7.5 in this case, algebraic convergence of the sparse grid PCM is verified for all types of distributions, which is also consistent with the analysis in Theorem 3.12.
Homogeneous Dirichlet-type boundary condition is applied on boundary ∂Ω. Since there is no analytical expression for the local limit, u 0 (x, ξ) is generated numerically based on a spectral method solver. The mean and standard deviation of the analytical local solution is calculated with the Monte Carlo (MC) method with 100, 000 samples.
Numerical results are provided in Figure 5. With 100 collocation points in PCM and fixed δ/h = 2.8, we (a) Convergence with δ, h → 0 in the physical space.
(b) Convergence with increasing sample numbers. investigate the spatial convergence of numerical solution with decreasing grid size from h = 1/4 to 1/32 in Figure 5(a). Second order convergence is observed, which is consistent with our analysis. In Figure 5(b) we employ fixed grid size h = 1/32, δ = 2.8h and show the convergence of numerical solutions obtained from MC simulations and the PCM results with respect to the increase of sample points. Here the sparse grid levels are taken as η = 1, · · · , 6. When η > N/ log(2) ≈ 3, sub-exponential convergence is observed for PCM, and the results also indicate that PCM can achieve a similar accuracy with much smaller number of sample points than MC.

Stochastic Representation and Karhunen-Loève Expansion
We use the Karhunen-Loève (K-L) expansion to represent the random field a(x, ω). In general, consider a square-integrable stochastic random field F (x, ω) defined on Ω × Ω p , where Ω is a subset of R d (d is the dimension) and Ω p is the sample space of a probability space (Ω p , F, P). If F (x, ω) has a constant mean and a continuous covariance function (also called kernel function) Ξ(x, y), then F (x, ω) can be represented by the following Karhunen-Loève (KL) expansion: where F 0 is the constant mean, (λ i , φ i ) are eigenpairs (i.e., eigenvalue and corresponding eigenfunction) of the kernel function Ξ, and ξ (i) are independent random variables with zero mean and unit variance. In practice, the summation is truncated up to N terms for computational purpose, where N is taken such that Notice that instead of 0.9, other numbers like 0.85 and 0.95 are also widely used. The notation of F (x, ω) is then replaced with F (x, ξ), where ξ = (ξ (1) , ξ (2) , . . . , ξ (N ) ). Here, we use the KL expansion to represent the random field a(x, ξ), and we consider the nonlocal diffusion coefficient as the harmonic mean of the local diffusion coefficient: A(x, y, ξ) = 2/(a −1 (x, ξ) + a −1 (y, ξ)) in our nonlocal model.  In this section, we consider the proposed approach on a randomly heterogeneous nonlocal problem with a given spatial correlation structure. We consider a 2D domain Ω = [−1, 1] × [−1, 1] with u 0 = 0 on ∂Ω. The local diffusivity coefficient a(x, ω) is modeled as a random field with constant mean a 0 = 4 and a continuous covariance function:
In this computational example, N (1) = N (2) = 2 is required to achieve the truncation error criteria (5.1). We assume ξ (i,j) , i, j = 1, 2, to be Gaussian random variables: ξ (i,j) ∼ N (0, 1), then investigate the convergence of u δ h,K to the local limit. Note that in this problem there exists no analytical expression for the local solution, we therefore generate u 0 (x, ξ) numerically based on a spectral method solver. Since the MC method suffers from a very slow convergence and requires more than 10 7 samples to achieve a desired precision, the mean and standard deviation of the local limit is calculated with the sparse grid PCM method with 2881 samples (sparse grid level η = 10).
Numerical results are provided in Figure 6. In Figure 6(a) we demonstrate the convergence in the physical space with uniform grids h ∈ {1/4, 1/8, 1/16, 1/32, 1/64} and fixed δ/h = 2.8, using 441 sample points (sparse grid level η = 5). Second order convergence is observed, which verifies the analysis in Theorem 3.12. To demonstrate the convergence in the parametric space with increasing sample numbers, in Figure 6(b) we employ a fixed grid size h = 1/64, δ = 2.8h, and show the numerical error to the local limit with sparse grid levels η = 1, · · · , 5 in PCM. The results again indicate that comparing with MC, the sparse grid PCM achieves a better accuracy with far smaller number of samples. An algebraic convergence is observed which is consistent with Theorem 3.12.

Summary and Discussion
Due to the limitation of computational resources and experimental resolutions, in many applications the physical parameters describing continuous properties of heterogeneous materials cannot be accurately characterized in all details. This problem becomes more acute in the nonlocal setting, due to its relatively lack of sparsity and correspondingly larger computational cost.
In this work we aim to consider the spatial variability of material properties in nonlocal models by solving randomly heterogeneous nonlocal diffusion problems. In particular, we have proposed an asymptotically compatible stochastic numerical method for randomly heterogeneous nonlocal diffusion problems, and provided rigorous mathematical analysis and error estimates. For spatial discretization, a meshfree discretization method with optimization-based quadrature rule is employed, which presents an up to O(h) consistency error to the nonlocal solution and an O(δ 2 ) convergence to the local limit. On the random parametric space, a probabilistic collocation method (PCM) with sparse grids is employed to sample the stochastic process. Since the fast convergence of the sparse grid PCM approach relies on smoothness of the solution in the random parametric space, we have proved that the nonlocal solution is analytic in the input random variable and therefore guarantees an at least algebraic convergence with increasing sample numbers. This work has for the first time provided a rigorous and comprehensive mathematical framework to add uncertainty quantification analysis onto pre-existing deterministic codes for nonlocal problems with guaranteed convergence in both physical and parametric spaces.
Future directions of this research will include the development and analysis of a generalization of this approach to the nonlocal mechanics problems, such as the peridynamics. A further extension of the meshfree discretization approach will also be pursued, to achieve a higher order convergence rate.