# From discrete to continuum: a generalized Euler–Maclaurin framework for scale effects in nanomechanics

Gabriele La Valle, Christian Soize

► To cite this version:

Gabriele La Valle, Christian Soize. From discrete to continuum: a generalized Euler–Maclaurin framework for scale effects in nanomechanics. Journal of Elasticity, 2025, 157 (3), pp.67. 10.1007/s10659-025-10158-0. hal-05225295

HAL Id: hal-05225295

https://univ-eiffel.hal.science/hal-05225295v1

Submitted on 8 Dec 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# From Discrete to Continuum: A Generalized Euler-Maclaurin Framework for Scale Effects in Nanomechanics

G. La Valle and C. Soize

Université Gustave Eiffel, MSME UMR 8208, Marne-la-Vallée, France

**Abstract** This paper introduces a novel approach for constructing, at the same scale, a continuum model equivalent to a given nanoscale discrete system, effectively capturing scale effects. Starting from a general formulation of *m*-particle interaction potentials for discrete particle systems, we propose the use of the Euler-Maclaurin (E-M) summation formula to construct the equivalent continuum model at the same scale. The proposed theory is developed for arbitrary 3D domains. The resulting novel continuum model captures scale effects in both statics and dynamics through additional edge, surface, and volume integrals, which are analytically obtained and driven by the nanostructure. For an arbitrary domain, the proposed approach provides a pathway for its integration into a computational framework. Since nanoscale and microscale systems are inevitably affected by uncertainties, the geometric and constitutive parameters must be modeled as random fields and their identification must be conducted within a statistical framework. Consequently, we present a discussion on this identification in a probabilistic framework.

**Keywords** Equivalent continuum, Scale effects, Discrete modeling, Euler-Maclaurin summation formula, Nanomechanics

## 1 Introduction

(i) *Motivation of the paper: connecting discrete to continuum models in nanomechanics*. Physical and biological phenomena are modeled using four principal theories: quantum mechanics and nanomechanics, which describe materials as discrete systems, and micromechanics and macromechanics, which model materials as continuum media. In recent years, technological advances have driven the development of nanomaterials, enabling the conception of mechanical metamaterials with unusual behaviors. Examples include smart and cognitive materials, nanoelectromechanical systems (NEMS), and microelectromechanical systems (MEMS) [1]. These fields offer significant opportunities for advancement and exploration in science and engineering. The progress in nanotechnology has also enabled the development of innovations such as protective coatings, energy-absorbing polymer-based nanomaterials for ballistic protection, self-healing materials capable of autonomously repairing cracks, etc.

To fully exploit the potential of nanomaterials, we must deepen our understanding of their behavior and interactions at the molecular level, as well as at the nanoscale and microscale. Despite remarkable advancements in computational performance, simulations of nanomaterials within the framework of nanomechanics, modeled as discrete systems, often remain prohibitive and challenging, especially as the number of degrees of freedom increase. These challenges have driven the development of micromechanics, which involves continuum models that can be seen as generalizations of classical theories, such as the theory of elasticity. These models reflect an inherent tendency toward a rapprochement between mechanics and physics. The development of such models is not only intrinsically valuable for advancing knowledge but also holds significant potential for multidisciplinary applications. Within the framework of micromechanics, two main approaches dominate.

(ii) *State of the art: homogenization and generalized continua*. For some, the term micromechanics is inseparable from the concept of *homogenization* [2–6], which allows a heterogeneous material modeled at microscale or mesoscale to be replaced by an equivalent homogeneous material at the macroscale, enabling the estimation of its effective properties (assuming scale separation). For example, when dealing with microstructures defined by distinct phases, the modeling can be performed at the microscale. These properties are obtained based on the knowledge of its constituents and the spatial distribution of its phases. If the microstructure cannot be
described at the microscale because of complexity (such as biological live tissues), then a stochastic modeling of the microstructure can be performed at the mesoscale, yielding apparent properties. Within this framework, it is essential to recall the *stochastic* homogenization, which has led to several works [7–15] and several analyses aimed at estimating a Representative Volume Element (RVE) [16–19].

For others, micromechanics evokes the *generalized continua models* pioneered by G. Piola [20, 21] and E. and F. Cosserat [22], which remained unnoticed for a long time. On the one hand, G. Piola was a pioneer in the theory of higher-order continua, extending classical elasticity to include higher-order gradients of displacement, whereas classical elasticity considers only first-order gradients. Higher-order continua have recently found applications in the modeling, development, and design of architectured metamaterials [23–27]. A second-gradient continuum model, which includes Piola peridynamics and Eringen nonlocal elasticity, has been proposed in [28]. This model has also been integrated into a stochastic framework [29] and adapted for an identification procedure [30]. On the other hand, E. and P. Cosserat pioneered the so-called Cosserat continua, which introduce the kinematic field of micro-rotations alongside the displacement field. Consequently, each material point is characterized by six degrees of freedom. In the literature, Cosserat continua are also referred to as oriented media, asymmetric continua, multipolar media, atomic lattices, couple-stress theories, etc. Essential contributions to the development of Cosserat continua were made by researchers such as Aero, Eringen, Green, Grioli, Günther, Koiter, Mindlin, Toupin, among others, though the list is far from exhaustive. A comprehensive overview of the subject can be found in [31–33]. A further development involves micromorphic continua, where the kinematic variables include both the displacement field and the set of micro-deformations, rather than just the micro-rotations as in Cosserat continua. In this regard, we refer to [34–36].

Finally, following Kröner work [37], we find a further classification of generalized continua that distinguishes between the so-called by Kröner *differential theories* and *integral theories* for modeling short- and long-range interactions, such as Van der Waals cohesive forces. Let us consider two points in a discrete system, with the free energy described in terms of pairwise interactions and the associated displacement fields. In his differential theory, a continuous displacement field is introduced, which can be expanded as a Taylor series around a reference point. This approach transforms the discrete description into a local representation, where the free energy is represented as a volume integral involving the continuous displacement field and its spatial derivatives. Within the framework of his integral theory, the Taylor expansion is discarded, and the discrete free energy is converted into a formulation involving integrals of the displacement field and its derivatives at the points of interest. This is achieved using the Euler-Maclaurin (E-M) summation formula, incorporating surface integrals to account for surface scale effects. The E-M summation formula was used by Kröner to establish a nonlocal continuum theory, which was later abandoned in the literature and subsequently replaced by the Eringen nonlocal models [38–42]. Kröner work and the use of the E-M summation formula, has been underestimated in the continuum mechanics literature and remains largely unknown. Actually, it could be an useful approach to transform a discrete model into an *equivalent continuum at the same scale* (or *quasi-continuum*).

(iii) *State of the art: from discrete to quasi-continuum models.* The transformation of discrete models into equivalent continuum ones at the same scale has been addressed in the literature. Notable approaches include the Piola Ansatz [43] and the Cauchy-Born rule [44]. A comprehensive treatment of the topic can be found in [45], which explains many of the key theoretical ideas behind multiscale modeling. Among the different approaches, we also recall the bridging scale methods, which involve a two-scale decomposition: the coarse scale is simulated using continuum methods, while the fine scale is modeled using discrete approaches [46–48]. Within a stochastic framework, [49] provides an additional contribution by presenting a methodology, inspired by meshless methods [50–55], of stochastic homogenization formulated on the domain of a microstructure that is not a RVE, implying the absence of scale separation.

(iv) *Novelties of the paper.* This paper discusses, for the first time, the Euler-Maclaurin (E-M) summation formula as an efficient tool to obtain an equivalent continuum model at the same scale, starting from *m*-particle interaction potential of discrete particle systems, and connecting differential and integral theories (as called by Kröner). It is important to emphasize that the proposed model is not intended as a substitute for molecular dynamics or quasi-continuum approaches. Rather, it aims to provide an analytically tractable and computationally efficient alternative in cases where discrete simulations become prohibitive and volume, surface, and edge size effects must be taken into account.

(v) *Organization of the paper.* In Section 2, recognizing the fundamental importance of the mechanics of discrete particle systems, we first establish a general form of *m*-particle interaction potentials that satisfy the objectivity principle and then consider the specific case that also satisfies invariance under point symmetry. In Section 3, we use Taylor series expansions to express nonlocal discrete potentials by means of a local representation. In Section 4, we discuss the further transformation to continuum models through the application of E-M summation formula. However, as E-M summation formula for functions on 3D domains was currently available only for discrete systems on regular grids in rectangular parallelepiped domains, we propose an approximation for arbitrary domains. This approximation, which is also a novel contribution, is based on coordinate transformations inspired by the finite element method. A notable novel consequence of
the proposed approach is the introduction of additional integrals, including novel edge, surface, and volume integrals, enabling the capture of scale effects in statics and dynamics arising from particle systems. Moreover, for an arbitrary domain, the proposed approach provides a pathway for its integration into a computational framework. Finally, in Section 5, we briefly discuss the model parameters identification within a probabilistic framework.

## 2 Interaction potential in particle systems

This section introduces a general form for the m-interaction potential for discrete particle systems. This general form is inspired by [56], which addresses the simplest case of interaction forces in systems characterized by pairwise interactions. While [56] considers pairwise interactions based on an axiomatic framework where the concept of force is primitive, here we extend this approach to m-particle interactions (beyond pairwise) within an axiomatic framework where energy is the primitive concept. Such m-particle interactions are fundamental for complex systems [57].

### 2.1 Definition of the particle system

Let $V$ be a bounded open domain of $\mathbb{R}^3$ and $\bar{V} = V \cup \partial V$ its closure, $\partial V$ is the boundary of $V$. If $\partial V$ also contains edges, we denote as $\partial\partial V$ the union of the edges and vertices (intersection of edges). For instance, if $V$ was a cube, $\partial V$ would consist of six planar faces (with edges and vertices), and $\partial\partial V$ would be the set of edges and vertices of the cube. While, if $V$ was a sphere, $\partial\partial V = \emptyset$, because the surface of a sphere (its boundary $\partial V$) has no boundary itself. We consider a system of $n$ particles in $\bar{V}$. Consequently, it is assumed that the number of particles $n$ is the sum $n = n_V + n_S + n_L$, where $n_V$ is the number of particles in volume $V$, $n_S$ on surface $\partial V \setminus \partial\partial V$, and $n_L$ on edges $\partial\partial V$. For instance, for the cube, there are $n_V$ particles inside the cube, $n_S$ inside the faces (do not include the edges and vertices), and $n_L$ on the edges and vertices. These particles are assumed to be located at positions $\mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_2}, \dots, \mathbf{x}^{\alpha_n}$ in the undeformed configuration and at positions $\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}, \dots, \mathbf{r}^{\alpha_n}$ (not necessarily in $\bar{V}$) in the deformed configuration. Moreover, it is assumed that domain $\bar{V}$ can be decomposed into a union of general curved hexahedron subdomains $\bar{V}^e = V^e \cup \partial V^e$, with $\cap_{e \in \mathcal{E}} V^e = \emptyset$ and such that $\bar{V} = \cup_{e \in \mathcal{E}} \bar{V}^e$, where $\mathcal{E} = \{1, 2, \dots, n_e\}$ denotes the set of indices of subdomains $\bar{V}^e$. Since any point in $\mathbb{R}^3$ is of zero measure for the surface measure and the volume measure in $\mathbb{R}^3$, we can always consider that there are no particles on $\cup_{e \in \mathcal{E}} V^e \setminus \partial V$, which means that all the particles in $V$ are inside $\cup_{e \in \mathcal{E}} V^e$. In this paper, we use the following convention. The symbol $\mathcal{A}$ represents the set of positive indices associated with each particle. For a fixed value of $m = 2, \dots, n$, the symbol $\mathcal{A}_m$ refers to the set of positive indices corresponding to all m-tuples formed within the system. The set $\mathcal{A}_1$ is defined as the intersection between these two, $\mathcal{A}_1 = \mathcal{A} \cap \mathcal{A}_m$.

**Hypothesis 1 (Related to the construction of an equivalent continuum)** Within the framework of the construction of an equivalent continuum model, it is assumed that $n_V \gg n_{eV}$, $n_S \gg n_{eS}$, and $n_L \gg n_{eL}$, where $n_{eV}$ is the number of subdomains $V^e$ (volumes), $n_{eS}$ is the number of subdomain surfaces that intersect $\partial V$, and $n_{eL}$ is the number of subdomain edges that intersect $\partial\partial V$.

### 2.2 Interaction potential satisfying objectivity principle

We consider interaction potential $\phi^{\text{def}}$ whose value in the deformed configuration is written as

$$ \phi^{\text{def}}(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}, \dots, \mathbf{r}^{\alpha_n}) = \sum_{(\alpha_1, \alpha_2) \in \mathcal{A}_2} \phi_2(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}) + \dots + \sum_{(\alpha_1, \alpha_2, \dots, \alpha_n) \in \mathcal{A}_n} \phi_n(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}, \dots, \mathbf{r}^{\alpha_n}), \quad (1) $$

where $\mathbf{r}^{\alpha_i} = \mathbf{r}(\mathbf{x}^{\alpha_i})$ denotes the position of the particle at position $\mathbf{x}^{\alpha_i}$ in the undeformed configuration, and $\phi_m$ denotes the m-interaction potential. The first term describes pairwise interactions, the second accounts for three-particle interactions, and so on. Each term reflects contributions involving progressively higher-order particle interactions. Proposition 1 provides a general form of the m-interaction potential $\phi_m$ consistent with the objectivity principle. In the following, we assume that the distance between two any particles (represented by material points) is not equal to zero for undeformed and deformed configurations, which implies that $\mathbf{r}(\mathbf{x}) \neq \mathbf{r}(\bar{\mathbf{x}})$ for all $\mathbf{x} \neq \bar{\mathbf{x}}$.

**Proposition 1** The m-interaction potential $\phi_m(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}, \dots, \mathbf{r}^{\alpha_m})$ of the system of particles satisfies the objectivity principle if
- for $m = 2$,

$$ \phi_2(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}) = h_2 (\|\mathbf{v}^{\alpha_2}\|^2), \qquad (2) $$

- for $m = 3$,

$$ \phi_3(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}, \mathbf{r}^{\alpha_3}) = h_3 (\|\mathbf{v}^{\alpha_2}\|^2, \|\mathbf{v}^{\alpha_3}\|^2, \langle\mathbf{v}^{\alpha_2}, \mathbf{v}^{\alpha_3}\rangle), \qquad (3) $$

- for $m \ge 4$,

$$ \phi_m(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}, \dots, \mathbf{r}^{\alpha_m}) = h_m (\{\|\mathbf{v}^{\alpha_i}\|^2, \langle\mathbf{v}^{\alpha_i}, \mathbf{v}^{\alpha_j}\rangle, \det([\mathbf{v}^{\alpha_i}, \mathbf{v}^{\alpha_j}, \mathbf{v}^{\alpha_k}])\}_{i \in \mathcal{I}_1, (i,j) \in \mathcal{I}_2, (i,j,k) \in \mathcal{I}_3}), \qquad (4) $$

$$ \text{where } \mathcal{I}_1 = \{2 \le i \le m\}, \mathcal{I}_2 = \{(i,j), 2 \le i < j \le m\}, \text{ and } \mathcal{I}_3 = \{(i,j,k), 2 \le i < j < k \le m\}. $$

Here, $\mathbf{v}^{\alpha_j} = \mathbf{r}^{\alpha_j} - \mathbf{r}^{\alpha_1}$ is the relative position vector from particle $\alpha_1$ to particle $\alpha_j$ in the deformed configuration, and $h_m$ is a function that satisfies the objectivity principle, i.e., it remains invariant under a change in observer frame.

*Proof.* For $m = 2$, the proof is known [56]. Let $m \ge 3$. If $\phi_m$ satisfies the objectivity principle, then for any special orthogonal matrix $[R] \in SO(3)$ ($[R]^T[R] = [I]$, $\det[R] = 1$) and vector $\mathbf{b} \in \mathbb{R}^3$, we must have

$$ \phi_m([R] \mathbf{r}^{\alpha_1} + \mathbf{b}, [R] \mathbf{r}^{\alpha_2} + \mathbf{b}, \dots, [R] \mathbf{r}^{\alpha_m} + \mathbf{b}) = \phi_m(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}, \dots, \mathbf{r}^{\alpha_m}), \qquad (5) $$

Without loss of generality, we select one reference vector $\mathbf{r}^{\alpha_1}$ and express the remaining vectors relative to $\mathbf{r}^{\alpha_1}$,

$$ \mathbf{v}^{\alpha_j} = \mathbf{r}^{\alpha_j} - \mathbf{r}^{\alpha_1}, \quad \text{for } j = 2, 3, \dots, m. \qquad (6) $$

Function $\phi_m(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}, \dots, \mathbf{r}^{\alpha_m})$ can now be expressed as

$$ \phi_m(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}, \dots, \mathbf{r}^{\alpha_m}) = \tilde{\phi}_m(\mathbf{v}^{\alpha_2}, \mathbf{v}^{\alpha_3}, \dots, \mathbf{v}^{\alpha_m}). \qquad (7) $$

The function $\tilde{\phi}_m$ is translation-invariant and must now remains invariant under rotations in SO(3). If $\tilde{\phi}_m$ depends on invariants under rotation of $\mathbf{v}^{\alpha_2}, \mathbf{v}^{\alpha_3}, \dots, \mathbf{v}^{\alpha_m}$, then it satisfies the objectivity principle. The considered invariants are

(i) **scalar products:** pairwise scalar product $\mathbf{v}^{\alpha_i} \cdot \mathbf{v}^{\alpha_j}$, which captures magnitudes and angle between vectors $\mathbf{v}^{\alpha_i}$ and $\mathbf{v}^{\alpha_j}$,

$$ \|\mathbf{v}^{\alpha_i}\|^2 = \langle \mathbf{v}^{\alpha_i}, \mathbf{v}^{\alpha_i} \rangle, \quad \langle \mathbf{v}^{\alpha_i}, \mathbf{v}^{\alpha_j} \rangle, \qquad (8) $$

(ii) **determinants:** determinant of 3 × 3 matrix formed by three vectors $\mathbf{v}^{\alpha_i}, \mathbf{v}^{\alpha_j}, \mathbf{v}^{\alpha_k}$ in $\mathbb{R}^3$, which is related to the volume defined by the three vectors,

$$ \det([\mathbf{v}^{\alpha_i}, \mathbf{v}^{\alpha_j}, \mathbf{v}^{\alpha_k}]). \qquad (9) $$

Let us prove that this determinant is invariant under rotation. Consider a rotation of the coordinate system represented by $[R] \in SO(3)$. Vectors $\mathbf{v}^{\alpha_i}, \mathbf{v}^{\alpha_j}, \mathbf{v}^{\alpha_k}$ are transformed as

$$ \mathbf{v}^{\alpha_i'} = [R] \mathbf{v}^{\alpha_i}, \quad \mathbf{v}^{\alpha_j'} = [R] \mathbf{v}^{\alpha_j}, \quad \mathbf{v}^{\alpha_k'} = [R] \mathbf{v}^{\alpha_k}. \qquad (10) $$

The determinant of the matrix formed by the rotated vectors is written as

$$ \det ([\mathbf{v}^{\alpha_i'} \mathbf{v}^{\alpha_j'} \mathbf{v}^{\alpha_k'}]) = \det ([ [R] \mathbf{v}^{\alpha_i} [R] \mathbf{v}^{\alpha_j} [R] \mathbf{v}^{\alpha_k} ] ). \qquad (11) $$

Since,

$$ [[R] \mathbf{v}^{\alpha_i} [R] \mathbf{v}^{\alpha_j} [R] \mathbf{v}^{\alpha_k}] = [R] [\mathbf{v}^{\alpha_i} \mathbf{v}^{\alpha_j} \mathbf{v}^{\alpha_k}], \qquad (12) $$

using the property of determinants for matrix products, we have

$$ \det ([[R] \mathbf{v}^{\alpha_i} [R] \mathbf{v}^{\alpha_j} [R] \mathbf{v}^{\alpha_k}]) = \det([R]) \det ([\mathbf{v}^{\alpha_i} \mathbf{v}^{\alpha_j} \mathbf{v}^{\alpha_k}]). \qquad (13) $$

Since $\det([R]) = 1$, we obtain

$$ \det ([\mathbf{v}^{\alpha_i'} \mathbf{v}^{\alpha_j'} \mathbf{v}^{\alpha_k'}]) = \det ([\mathbf{v}^{\alpha_i} \mathbf{v}^{\alpha_j} \mathbf{v}^{\alpha_k}]), \qquad (14) $$

which finishes the proof.

Thus, for $m = 2, m = 3$, and $m \ge 4$, function $\phi_m$ can be expressed as an arbitrary function $h_m$ of these invariants, as shown in Eqs. (2), (3), and (4). $\square$

*Remark 1*
(i) It should be noted that det([v<sup>αi</sup> v<sup>αj</sup> v<sup>αk</sup>]) is invariant under rotation but not under point symmetry (also called inversion symmetry), because det([-v<sup>αi</sup> - v<sup>αj</sup> - v<sup>αk</sup>]) = -det([v<sup>αi</sup> v<sup>αj</sup> v<sup>αk</sup>]).

(ii) Among systems that are not invariant under point symmetry, there are particle systems with m-interaction potential, with m ≥ 4, having chiral properties. From now on, we will consider systems that are invariant under point symmetry.

**Proposition 2** The m-interaction potential $\phi_m(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}, \dots, \mathbf{r}^{\alpha_m})$ of the system of particles satisfies the objectivity principle and is invariant under point symmetry if

- for $m=2$,

$$ \phi_2(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}) = h_2 (\|\mathbf{v}^{\alpha_2}\|^2), \qquad (15) $$

- for $m \ge 3$,

$$ \phi_m(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}, \dots, \mathbf{r}^{\alpha_m}) = h_m (\{\|\mathbf{v}^{\alpha_i}\|^2, \langle\mathbf{v}^{\alpha_i}, \mathbf{v}^{\alpha_j}\rangle\}_{i \in \mathcal{I}_1, (i,j) \in \mathcal{I}_2}), \qquad (16) $$

$$ \text{where } \mathcal{I}_1 = \{2 \le i \le m\} \text{ and } \mathcal{I}_2 = \{(i, j), 2 \le i < j \le m\}. $$

Here, $\mathbf{v}^{\alpha_j} = \mathbf{r}^{\alpha_j} - \mathbf{r}^{\alpha_1}$ is the relative position vector from particle $\alpha_1$ to particle $\alpha_j$ in the deformed configuration, and $h_m$ is a function that satisfies the objectivity principle and is invariant under point symmetry.

*Proof.* This result follows directly from Proposition 1. □

### 3 Transforming nonlocal interaction potentials into local representations

We consider systems invariant under point symmetry, i.e., systems in which the simultaneous point symmetry of all three coordinates transforms the system into an equivalent one. In the previous section, we obtained a general form of an m-interaction potential $\phi_m$ that satisfies both the principle of objectivity and invariance under point symmetry. Potential $\phi_m$ is expressed as a function of the particle positions $\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}, \dots, \mathbf{r}^{\alpha_m}$ in the deformed configuration. For this reason, it is referred to as nonlocal interaction potential, as it explicitly depends on the overall spatial configuration of the system.

In the following, we use a Taylor expansion combined with an approach inspired by Piola [58], to show the possibility of representing the nonlocal interaction potential $\phi_m$ in a local representation. Specifically, this local representation depends explicitly on the radius vector $\mathbf{x} \mapsto \mathbf{r}(\mathbf{x})$ only at $\mathbf{x}^{\alpha_1}$, its spatial gradient $\mathbf{x} \mapsto \nabla_{\mathbf{x}}\mathbf{r}(\mathbf{x})$ only at $\mathbf{x}^{\alpha_1}$, and the initial positions $\mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_2}, \dots, \mathbf{x}^{\alpha_m}$.

For simplicity, we consider models where higher-order derivatives of $\mathbf{r}$ are negligible when compared to the first-order ones.

#### 3.1 The 2-interaction potential

The 2-interaction potential (or pairwise potential) $\phi_2$ describes the interaction between pairs in the particle system. For particles located at $\mathbf{r}^{\alpha_1}$ and $\mathbf{r}^{\alpha_2}$, in accordance with Section 2, we consider

$$ \phi_2(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}) = h_2 (\|\mathbf{r}^{\alpha_1} - \mathbf{r}^{\alpha_2}\|^2). \qquad (17) $$

When only the first-order derivatives of $\mathbf{r}$ can be used, the first-order Taylor expansion of $\mathbf{r}^{\alpha_2}$ around $\mathbf{r}^{\alpha_1}$ allows $\mathbf{r}^{\alpha_2}$ to be defined as

$$ \mathbf{r}^{\alpha_2} = \mathbf{r}^{\alpha_1} + [\mathbf{F}(\mathbf{x}^{\alpha_1})](\mathbf{x}^{\alpha_2} - \mathbf{x}^{\alpha_1}), \qquad (18) $$

where $[\mathbf{F}(\mathbf{x}^{\alpha_1})] = [\nabla_{\mathbf{x}}\mathbf{r}(\mathbf{x}^{\alpha_1})]$ is the gradient of $\mathbf{r}$ at $\mathbf{x}^{\alpha_1}$. Consequently, the difference between the squared distances in the deformed and undeformed configurations of particles $\alpha_1$ and $\alpha_2$ is given by

$$ \|\mathbf{r}^{\alpha_2} - \mathbf{r}^{\alpha_1}\|^2 - \|\mathbf{x}^{\alpha_2} - \mathbf{x}^{\alpha_1}\|^2 = 2 \langle [\mathbf{E}(\mathbf{x}^{\alpha_1})](\mathbf{x}^{\alpha_2} - \mathbf{x}^{\alpha_1}), (\mathbf{x}^{\alpha_2} - \mathbf{x}^{\alpha_1}) \rangle, \qquad (19) $$

where $[\mathbf{E}(\mathbf{x}^{\alpha_1})] = \frac{1}{2}([\mathbf{F}(\mathbf{x}^{\alpha_1})]^T[\mathbf{F}(\mathbf{x}^{\alpha_1})] - [I])$ is the Green-Saint Venant strain tensor at $\mathbf{x}^{\alpha_1}$. The second-order Taylor expansion of $\phi_2$ around $\|\mathbf{x}^{\alpha_2} - \mathbf{x}^{\alpha_1}\|^2$ yields

$$ \phi_2(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}) - \phi_2(\mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_2}) = a_{ij}(\mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_2})[\mathbf{E}(\mathbf{x}^{\alpha_1})]_{ij} + \frac{1}{2} a_{ijkh}(\mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_2})[\mathbf{E}(\mathbf{x}^{\alpha_1})]_{kh}[\mathbf{E}(\mathbf{x}^{\alpha_1})]_{ij}, \qquad (20) $$

where

$$ a_{ij}(\mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_2}) = 2 h_2' (\|\mathbf{x}^{\alpha_2} - \mathbf{x}^{\alpha_1}\|^2) (x_i^{\alpha_2} - x_i^{\alpha_1})(x_j^{\alpha_2} - x_j^{\alpha_1}), \qquad (21) $$

$$ a_{ijkh}(\mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_2}) = 4 h_2'' (\|\mathbf{x}^{\alpha_2} - \mathbf{x}^{\alpha_1}\|^2) (x_i^{\alpha_2} - x_i^{\alpha_1})(x_j^{\alpha_2} - x_j^{\alpha_1})(x_k^{\alpha_2} - x_k^{\alpha_1})(x_h^{\alpha_2} - x_h^{\alpha_1}). \qquad (22) $$
where $h'_2$ and $h''_2$ are the first and second derivatives of $y \mapsto h_2(y)$ with respect to $y$. We define the deformation potential $\varphi_2$, which serves as a local representation of $\phi_2$, as follows

$$ \begin{aligned} \varphi_2(\mathbf{x}^{\alpha_1}, [\mathrm{E}(\mathbf{x}^{\alpha_1})]) &= \sum_{\alpha_2 \mid (\alpha_1, \alpha_2) \in \mathcal{A}_2} \left( \phi_2(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}) - \phi_2(\mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_2}) \right) \\ &= c_{ij}(\mathbf{x}^{\alpha_1})[\mathrm{E}(\mathbf{x}^{\alpha_1})]_{ij} + \frac{1}{2} c_{ijkh}(\mathbf{x}^{\alpha_1})[\mathrm{E}(\mathbf{x}^{\alpha_1})]_{kh} [\mathrm{E}(\mathbf{x}^{\alpha_1})]_{ij}, \end{aligned} \quad (23) $$

where

$$ c_{ij}(\mathbf{x}^{\alpha_1}) = \sum_{\alpha_2 \mid (\alpha_1, \alpha_2) \in \mathcal{A}_2} a_{ij}(\mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_2}), \quad (24) $$

and

$$ c_{ijkh}(\mathbf{x}^{\alpha_1}) = \sum_{\alpha_2 \mid (\alpha_1, \alpha_2) \in \mathcal{A}_2} a_{ijkh}(\mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_2}). \quad (25) $$

*Remark 2* Note that $c_{ij}$ is symmetric, and $c_{ijkh}$ satisfies the usual symmetry properties.

## 3.2 The 3-interaction potential

Proposition 2 provides a very general form for $\phi_3$ satisfying the objectivity principle and invariance under point symmetry. For simplicity, we remove the dependency on the norms and we write $\phi_3$ as

$$ \phi_3(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}, \mathbf{r}^{\alpha_3}) = h_3(\langle \mathbf{r}^{\alpha_2} - \mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_3} - \mathbf{r}^{\alpha_1} \rangle). \quad (26) $$

The first-order Taylor expansion of $\mathbf{r}^{\alpha_2}$ and $\mathbf{r}^{\alpha_3}$ around $\mathbf{r}^{\alpha_1}$ yields

$$ \langle \mathbf{r}^{\alpha_2} - \mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_3} - \mathbf{r}^{\alpha_1} \rangle - \langle \mathbf{x}^{\alpha_2} - \mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_3} - \mathbf{x}^{\alpha_1} \rangle = 2 \langle [\mathrm{E}(\mathbf{x}^{\alpha_1})](\mathbf{x}^{\alpha_2} - \mathbf{x}^{\alpha_1}), (\mathbf{x}^{\alpha_3} - \mathbf{x}^{\alpha_1}) \rangle. \quad (27) $$

Thus, the second-order Taylor expansion of $\phi_3$ around $\langle \mathbf{x}^{\alpha_2} - \mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_3} - \mathbf{x}^{\alpha_1} \rangle$ gives

$$ \phi_3(\mathbf{r}^{\alpha_1}, \mathbf{r}^{\alpha_2}, \mathbf{r}^{\alpha_3}) - \phi_3(\mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_2}, \mathbf{x}^{\alpha_3}) = a_{ij}(\mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_2}, \mathbf{x}^{\alpha_3})[\mathrm{E}(\mathbf{x}^{\alpha_1})]_{ij} + \frac{1}{2} a_{ijkh}(\mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_2}, \mathbf{x}^{\alpha_3})[\mathrm{E}(\mathbf{x}^{\alpha_1})]_{ij} [\mathrm{E}(\mathbf{x}^{\alpha_1})]_{kh}, \quad (28) $$

where the coefficients are

$$ a_{ij}(\mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_2}, \mathbf{x}^{\alpha_3}) = 2 h'_3 (\langle \mathbf{x}^{\alpha_2} - \mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_3} - \mathbf{x}^{\alpha_1} \rangle) (x_i^{\alpha_3} - x_i^{\alpha_1}) (x_j^{\alpha_2} - x_j^{\alpha_1}), \quad (29) $$

and

$$ a_{ijkh}(\mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_2}, \mathbf{x}^{\alpha_3}) = 4 h_3'' (\langle \mathbf{x}^{\alpha_3} - \mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_2} - \mathbf{x}^{\alpha_1} \rangle) (x_i^{\alpha_3} - x_i^{\alpha_1}) (x_j^{\alpha_2} - x_j^{\alpha_1}) (x_k^{\alpha_3} - x_k^{\alpha_1}) (x_h^{\alpha_2} - x_h^{\alpha_1}). \quad (30) $$

Summing over $\alpha_2$ and $\alpha_3$, we define the deformation potential $\varphi_3$, which serves as a local representation of $\phi_3$, as follows

$$ \varphi_3(\mathbf{x}^{\alpha_1}, [\mathrm{E}(\mathbf{x}^{\alpha_1})]) = c_{ij}(\mathbf{x}^{\alpha_1})[\mathrm{E}(\mathbf{x}^{\alpha_1})]_{ij} + \frac{1}{2} c_{ijkh}(\mathbf{x}^{\alpha_1})[\mathrm{E}(\mathbf{x}^{\alpha_1})]_{ij} [\mathrm{E}(\mathbf{x}^{\alpha_1})]_{kh}, \quad (31) $$

where

$$ c_{ij}(\mathbf{x}) = \sum_{(\alpha_2, \alpha_3) \mid \{\alpha_1, \alpha_2, \alpha_3\} \in \mathcal{A}_3} a_{ij}(\mathbf{x}, \mathbf{x}^{\alpha_2}, \mathbf{x}^{\alpha_3}), \quad (32) $$

$$ c_{ijkh}(\mathbf{x}) = \sum_{(\alpha_2, \alpha_3) \mid \{\alpha_1, \alpha_2, \alpha_3\} \in \mathcal{A}_3} a_{ijkh}(\mathbf{x}, \mathbf{x}^{\alpha_2}, \mathbf{x}^{\alpha_3}). \quad (33) $$

*Remark 3* Note that, in general, $c_{ij}$ is not symmetric and $c_{ijkh}$ does not satisfy $c_{ijkh} = c_{jikh}$ and $c_{ijkh} = c_{ikhj}$. It satisfies only $c_{ijkh} = c_{khij}$.
### 3.3 Generalizing to the m-interaction potential

The same reasoning applies to the m-interaction potential $\phi_m$. If an m-interaction potential $\phi_m$ satisfies Proposition 2 (under objectivity principle and point symmetry invariance) and higher-order derivatives are negligible, the deformation potential $\phi_m$, which is a local representation of $\phi_m$, can be defined as

$$ \varphi_m(\mathbf{x}^{\alpha_1}, [\mathrm{E}(\mathbf{x}^{\alpha_1})]) = c_{ij}(\mathbf{x}^{\alpha_1})[\mathrm{E}(\mathbf{x}^{\alpha_1})]_{ij} + \frac{1}{2} c_{ijkh}(\mathbf{x}^{\alpha_1})[\mathrm{E}(\mathbf{x}^{\alpha_1})]_{kh}[\mathrm{E}(\mathbf{x}^{\alpha_1})]_{ij}, \quad (34) $$

where $c_{ijkh} = c_{khij}$. Tensors $c_{ij}$, $a_{ij}$, $c_{ijkh}$, and $a_{ijkh}$ are directly constructed from the m-interaction potential $\phi_m$ and the undeformed configuration of the discrete particle system, as detailed in Subsections 3.1 and 3.2.

**Remark 4** These results show that a discrete nonlocal model, using truncated Taylor expansions, can be rewritten by a local representation. Based on this local representation and on the use of the Euler-Maclaurin summation formula, we construct an equivalent continuum model while considering scale effects.

## 4 From the particle system to an equivalent continuum

In this section, we construct a continuum model equivalent to the particle system via the Euler-Maclaurin (E-M) summation formula. This approach is inspired by [37], where the transformation from discrete to continuum systems using E-M summation formula is briefly mentioned but such an approach has not been developed in the literature. The reasons why this approach fell out of use in the literature is unclear. One plausible explanation could be the difficulty in generalizing the E-M summation formula from 1D systems to 3D systems with arbitrary shapes, except for the simple parallelepiped case, where the extension is relatively straightforward.

In this section, we propose an approximate extension of the E-M summation formula for arbitrary systems using coordinate transformations inspired by finite element method [59,60]. This result enables the transformation from any discrete particle system, consisting of a sufficiently large number of particles, to an equivalent continuum representation capable of capturing scale effects. The proposed methodology is particularly useful when the nano- or microstructure is unknown a priori, conditions under which traditional techniques developed for architected materials are challenging to apply, particularly in the case of systems with random structures.

### 4.1 The Euler-Maclaurin summation formula

Under appropriate hypotheses, the E-M summation formula relates the sum of a real function $\xi \mapsto f(\xi)$ over the discrete set $\{0, 1, \dots, n\}$, corresponding to the integer points within the interval $[0, n]$, to the corresponding integral, with correction terms accounting for the discrepancy between the sum and the integral. According to [61], the E-M summation formula for a unit step size is expressed as

$$ \sum_{k=1}^{n-1} f_k = \int_{0}^{n} f(\xi) d\xi - \frac{f(0) + f(n)}{2} + \dots \quad (35) $$

Let $\xi \mapsto g(\xi)$ be a real function defined on an interval $[-1, 1]$ and satisfying appropriate hypotheses, where $\xi$ takes discrete values $\xi_k = -1 + kh$ with $h = 2/n$ and $n-1$ being the number of internal discretization points. Then, Eqs. (35) give

$$ \sum_{k=1}^{n-1} g(\xi_k) = \frac{1}{h} \int_{-1}^{1} g(\xi) d\xi - \frac{g(-1) + g(1)}{2} \dots \quad (36) $$

### 4.2 The Euler-Maclaurin summation formula for 3D summation in a cube

Let $\Omega =] -1, 1]^3$ be the open cube in $\mathbb{R}^3$ with boundary $\partial\Omega$. We denote as $\partial\partial\Omega$ the union of edges and vertices. Let $\xi \mapsto f(\xi)$ be the real function defined on $\bar{\Omega} = \Omega \cup \partial\Omega$. We define the subsets of integers $\mathcal{J}_n = \{1, \dots, n-1\}$. We now consider a regular cubic grid of domain $\bar{\Omega}$ with the same step size $h = 2/n$ and $n-1$ being the number of internal discretization points in each direction $\xi_1, \xi_2$, and $\xi_3$. The internal grid
points are $\xi^{ijk} = (-1 + i h, -1 + j h, -1 + k h)$ with $(i, j, k)$ belonging to $\mathcal{J}_n$. Under appropriate hypotheses for $f$, applying the E-M summation formula in Eq. (36) successively in each direction, yields

$$ \sum_{(i,j,k) \in \mathcal{J}_n} f(\xi^{ijk}) = \frac{1}{h^3} \int_{\Omega} f(\xi) d\xi - \frac{1}{2h^2} \int_{\partial\Omega} f(\xi) d\sigma \dots \quad (37) $$

This approach enables the transformation from a discrete to a continuum model, relying entirely on energetic considerations.

## 4.3 Construction of the equivalent continuum via coordinate transformation and statistical approach

Let $V$ be the open domain previously defined (see Section 2), which contains the particles located at positions $\mathbf{x}^{\alpha_1}, \mathbf{x}^{\alpha_2}, \dots, \mathbf{x}^{\alpha_n}$. It is assumed that domain $\bar{V}$ can be decomposed into a union of general curved hexahedron subdomains $\bar{V}^e = V^e \cup \partial V^e$, such that $\bar{V} = \bigcup_{e \in \mathcal{E}} \bar{V}^e$, where $\mathcal{E}$ denotes the set of indices associated with subdomains $\bar{V}^e$. We denote by $\partial \partial V^e$ the union of edges and vertices. Let $s^e$ be the coordinate transformation

$$ s^e : \xi \in \bar{\Omega} \mapsto \mathbf{x} = s^e(\xi) \in \bar{V}^e, \quad (38) $$

which maps the local coordinates $\xi$ of the reference cube $\bar{\Omega} = [-1, 1]^3$ to the Cartesian coordinates $\mathbf{x}$ of sub-domain $\bar{V}^e$. For any system of particles characterized by an $m$-interaction potential, as discussed in Section 3, we can introduce the deformation potential $\varphi_m(\mathbf{x}^{\alpha_1}, [\mathbf{E}(\mathbf{x}^{\alpha_1})])$, where $[\mathbf{E}(\mathbf{x})]$, representing the Green-Saint Venant strain tensor, can be expressed in terms of the displacement $\mathbf{u}(\mathbf{x}) = \mathbf{r}(\mathbf{x}) - \mathbf{x}$, which serves as the unknown of the mechanical problem. As a consequence, we can define the function $\psi_m$ such that

$$ \psi_m(\mathbf{x}, [\nabla_{\mathbf{x}} \mathbf{u}(\mathbf{x})]) = \varphi_m(\mathbf{x}, [\mathbf{E}(\mathbf{x})]), \quad (39) $$

which implies that

$$ \sum_{\alpha_1 \in \mathcal{A}_1} \psi_m(\mathbf{x}^{\alpha_1}, [\nabla_{\mathbf{x}} \mathbf{u}(\mathbf{x}^{\alpha_1})]) = \sum_{\alpha_1 \in \mathcal{A}_1} \varphi_m(\mathbf{x}^{\alpha_1}, [\mathbf{E}(\mathbf{x}^{\alpha_1})]). \quad (40) $$

By applying the coordinate transformation in Eq. (38) to the second term of Eq. (40), yields

$$ \sum_{\alpha_1 \in \mathcal{A}_1} \varphi_m(\mathbf{x}^{\alpha_1}, [\mathbf{E}(\mathbf{x}^{\alpha_1})]) = \sum_{e \in \mathcal{E}, \alpha_1^e \in \mathcal{A}_1^e} f_m^e(\xi_1^{\alpha_1^e}), \quad (41) $$

where, to simplify, the notation we introduce the function $\xi \mapsto f_m^e(\xi)$ such that

$$ f_m^e(\xi_1^{\alpha_1^e}) = \psi_m(s^e(\xi_1^{\alpha_1^e}), [\nabla_{\xi} \mathbf{u}(s^e(\xi_1^{\alpha_1^e}))][\mathbf{J}(\xi_1^{\alpha_1^e})]^{-1}). \quad (42) $$

In Eq. (41), $\mathcal{A}_1^e$ denotes the set of indices associated with the particles in $V^e$, $\xi_1^{\alpha_1^e}$ are the local coordinates corresponding to the Cartesian coordinates $\mathbf{x}_1^{\alpha_1^e}$ of particles inside the subdomain $V^e$, and $[\mathbf{J}(\xi)]_{ij} = \partial x_i / \partial \xi_j$ denotes the Jacobian matrix in local coordinates.

**Definition 1 (Step sizes definition related to the particle system)** We define the step sizes $h_{V^e}, h_{S^e}$, and $h_{L^e}$ for the volume, surface, and edge effects. These are written as

$$ h_{V^e} = \frac{2}{n_{V^e}^{1/3} + 1}, \quad h_{S^e} = \frac{2}{n_{S^e}^{1/2} + 1}, \quad h_{L^e} = \frac{2}{n_{L^e} + 1}, \quad (43) $$

where $n_{V^e}, n_{S^e}$, and $n_{L^e}$ are the number of particles within the subdomains $V^e, \partial V^e \cap (\partial V \setminus \partial \partial V)$, and $\partial \partial V^e \cap \partial \partial V$.

**Hypothesis 2 (Hypothesis on the spatial distribution of particles)** The E-M summation formula is based on a regular grid of points with constant step size in $\Omega$. On the contrary, the spatial distribution of the particles in $V^e$ and, consequently, in $\Omega$ are randomly distributed. We then assume that $\mathbf{x}_1^{\alpha_1^e}$ is the realization of a random variable $\mathbf{X}_1^{\alpha_1^e}$ and that $\xi_1^{\alpha_1^e}$ is the realization of the random variable $\Xi_1^{\alpha_1^e} = s^e(\mathbf{X}_1^{\alpha_1^e})$ in $\Omega$. Random variable $\Xi_1^{\alpha_1^e}$ is written as $\Xi_1^{\alpha_1^e} = \xi_1^{\alpha_1^e} + \mathbf{Z}_1^{\alpha_1^e}$, where $\xi_1^{\alpha_1^e}$ is the point of the regular grid with step size $h_{V^e}$ for each direction in $\Omega$ and $\mathbf{Z}_1^{\alpha_1^e}$ is a second-order real-valued random variable with probability distribution $P_{\mathbf{Z}_1^{\alpha_1^e}}(dz)$ on $\mathbb{R}^3$, which is centered and such that $\mathbb{E}\{(Z_1^{\alpha_1^e})^2\} = \mathcal{O}(h_{V^e}^2)$, which implies $\mathbb{E}\{Z_1^{\alpha_1^e} Z_j^{\alpha_1^e}\} = \mathcal{O}(h_{V^e}^2)$, where $\mathbb{E}$ denotes the mathematical expectation. Since $\bar{\Omega}$ is a cube, an adapted similar hypothesis is used for the surfaces and edges contributions.
Remark 5 For instance, if $P_{\mathbf{Z}^{\alpha_1^e}}(dz)$ was $\frac{1}{h_{Ve}^3}1_{[-h_{Ve}/2, h_{Ve}/2]^3}(z) dz$, then we would have $\mathbb{E}\{Z_i^{\alpha_1^e}\} = 0$ and $\mathbb{E}\{(Z_i^{\alpha_1^e})^2\} = h_{Ve}^2/12$.

**Hypothesis 3 (Hypothesis for constructing an equivalent continuum)** It is assumed that there exists a probability distribution $P_{\mathbf{Z}^{\alpha_1^e}}(dz)$ such that the total contribution of the interaction potential, expressed in Cartesian coordinates $\mathbf{x}^{\alpha_1}$, on volume $V$, surface $\partial V \setminus \partial \partial V$, edges and vertices $\partial \partial V$, is equal to the mathematical expectation of the same total contribution, expressed in the random local coordinates $\Xi^{\alpha_1^e}$, which, for each type of contribution, is written as

$$ \sum_{\alpha_1 \in \mathcal{A}_{1V}} \varphi_m(\mathbf{x}^{\alpha_1}, [\nabla_{\mathbf{x}} \mathbf{u}(\mathbf{x}^{\alpha_1})]) = \sum_{e \in \mathcal{E}, \alpha_1^e \in \mathcal{A}_{1V}^e} \mathbb{E}\{f_m^e(\Xi^{\alpha_1^e})\}, \quad (44) $$

$$ \sum_{\alpha_1 \in \mathcal{A}_{1\partial V \setminus \partial \partial V}} \varphi_m(\mathbf{x}^{\alpha_1}, [\nabla_{\mathbf{x}} \mathbf{u}(\mathbf{x}^{\alpha_1})]) = \sum_{e \in \mathcal{E}, \alpha_1^e \in \mathcal{A}_{1\partial V \setminus \partial \partial V}^e} \mathbb{E}\{f_m^e(\Xi^{\alpha_1^e})\}, \quad (45) $$

and

$$ \sum_{\alpha_1 \in \mathcal{A}_{1\partial \partial V}} \varphi_m(\mathbf{x}^{\alpha_1}, [\nabla_{\mathbf{x}} \mathbf{u}(\mathbf{x}^{\alpha_1})]) = \sum_{e \in \mathcal{E}, \alpha_1^e \in \mathcal{A}_{1\partial \partial V}^e} \mathbb{E}\{f_m^e(\Xi^{\alpha_1^e})\}, \quad (46) $$

where $\mathcal{A}_{1V}$, $\mathcal{A}_{1\partial V \setminus \partial \partial V}$, and $\mathcal{A}_{1\partial \partial V}$ denote the sets of indices $\alpha_1 \in \mathcal{A}_{1}$ such that $\mathbf{x}^{\alpha_1} \in V$, $\mathbf{x}^{\alpha_1} \in (\partial V \setminus \partial \partial V)$, and $\mathbf{x}^{\alpha_1} \in \partial \partial V$, respectively.

**Lemma 1 (Rewriting Hypothesis 3 in terms of local coordinates of a regular grid)** Under Hypotheses 2 and 3, Eqs. (44), (45), and (46) can be rewritten, at order $\mathcal{O}(h_{Ve}^2)$, $\mathcal{O}(h_{Se}^2)$, and $\mathcal{O}(h_{Le}^2)$, as

$$ \sum_{\alpha_1 \in \mathcal{A}_{1V}} \varphi_m(\mathbf{x}^{\alpha_1}, [\nabla_{\mathbf{x}} \mathbf{u}(\mathbf{x}^{\alpha_1})]) = \sum_{e \in \mathcal{E}, \alpha_1^e \in \mathcal{A}_{1V}^e} f_m^e(\underline{\xi}^{\alpha_1^e}). \quad (47) $$

$$ \sum_{\alpha_1 \in \mathcal{A}_{1\partial V \setminus \partial \partial V}} \varphi_m(\mathbf{x}^{\alpha_1}, [\nabla_{\mathbf{x}} \mathbf{u}(\mathbf{x}^{\alpha_1})]) = \sum_{e \in \mathcal{E}, \alpha_1^e \in \mathcal{A}_{1\partial V \setminus \partial \partial V}^e} f_m^e(\underline{\xi}^{\alpha_1^e}), \quad (48) $$

and

$$ \sum_{\alpha_1 \in \mathcal{A}_{1\partial \partial V}} \varphi_m(\mathbf{x}^{\alpha_1}, [\nabla_{\mathbf{x}} \mathbf{u}(\mathbf{x}^{\alpha_1})]) = \sum_{e \in \mathcal{E}, \alpha_1^e \in \mathcal{A}_{1\partial \partial V}^e} f_m^e(\underline{\xi}^{\alpha_1^e}), \quad (49) $$

*Proof of Lemma 1.* Let us prove Eq. (47). The Taylor expansion of $f_m^e(\Xi^{\alpha_1^e})$ is written as

$$ f_m^e(\Xi^{\alpha_1^e}) = f_m^e(\underline{\xi}^{\alpha_1^e}) + \frac{\partial f_m^e(\underline{\xi}^{\alpha_1^e})}{\partial \xi_i} Z_i^{\alpha_1} + \frac{1}{2} \frac{\partial^2 f_m^e(\underline{\xi}^{\alpha_1^e})}{\partial \xi_i \partial \xi_j} Z_i^{\alpha_1} Z_j^{\alpha_1} + \dots \quad (50) $$

Taking the mathematical expectation of the two member of Eq. (50) and since $Z^{\alpha_i}$ and $Z^{\alpha_j}$ are centered and independent for $i$ not equal to $j$, we obtain

$$ \mathbb{E}\{f_m^e(\Xi^{\alpha_1^e})\} = f_m^e(\underline{\xi}^{\alpha_1^e}) + \frac{1}{2} \frac{\partial^2 f_m^e(\underline{\xi}^{\alpha_1^e})}{\partial \xi_i^2} \mathbb{E}\{(Z_i^{\alpha_1^e})^2\} + \dots \quad (51) $$

Since $\mathbb{E}\{(Z_i^{\alpha_1^e})^2\} = \mathcal{O}(h_{Ve}^2)$ we have proven Eq. (47). The same reasoning can be performed to proof the formula for contributions related to Eqs. 48 and 49. $\square$

## 4.4 Formulation of the continuum model using the Euler-Maclaurin summation formula

The explicit formulation of the equivalent continuum is based on Eqs. (47), (48) and (49) resulting from Eqs. (44), (45), and (46).

**Hypothesis 4 (Negligible contribution of derivatives of $f_m^e$)** As in [37]), it is assumed that the rapidly varying part of $f_m^e$ is related to short-range effects. Therefore, it is expected the rapid convergence of the right member of Eq. (37), where $f$ is substituted by $f_m^e$. Consequently, the contributions related to the derivatives of $f_m^e$ are assumed to be negligible.
**Proposition 3 (Equivalent continuum model)** Under hypotheses 1 to 4, and assuming that $h_{V^e}$, $h_{S^e}$, and $h_{L^e}$ are $\ll 1$, the equivalent continuum of the particle system introduced in Section 2 is defined by the following contributions,

$$ \sum_{\alpha_1 \in A_{1V}} \varphi_m(\mathbf{x}^{\alpha_1}, [\nabla_{\mathbf{x}} \mathbf{u}(\mathbf{x}^{\alpha_1})]) = \sum_{e \in \mathcal{E}} \left( \frac{1}{h_{V^e}^3} \int_{\Omega} f_m^e(\xi) \, d\xi - \frac{1}{2h_{V^e}^2} \int_{\partial\Omega} f_m^e(\xi) \, d\sigma \right), \quad (52) $$

$$ \sum_{\alpha_1 \in A_{1\partial V}} \varphi_m(\mathbf{x}^{\alpha_1}, [\nabla_{\mathbf{x}} \mathbf{u}(\mathbf{x}^{\alpha_1})]) = \sum_{e \in \mathcal{E}} \left( \frac{1}{h_{S^e}^2} \int_{\partial\Omega} f_m^e(\xi) \mathbf{1}_{\partial V}(\mathbf{s}^e(\xi)) \, d\sigma - \frac{1}{2h_{S^e}^2} \int_{\partial\partial\Omega} f_m^e(\xi) \mathbf{1}_{\partial V}(\mathbf{s}^e(\xi)) \, d\ell \right), \quad (53) $$

and

$$ \sum_{\alpha_1 \in A_{1\partial\partial V}} \varphi_m(\mathbf{x}^{\alpha_1}, [\nabla_{\mathbf{x}} \mathbf{u}(\mathbf{x}^{\alpha_1})]) = \sum_{e \in \mathcal{E}} \left( \frac{1}{h_{L^e}} \int_{\partial\partial\Omega} f_m^e(\xi) \mathbf{1}_{\partial\partial V}(\mathbf{s}^e(\xi)) \, d\ell + \text{neglected term related to } \partial\partial\partial\Omega \right), \quad (54) $$

where $\mathbf{x} \mapsto \mathbf{1}_{\partial V}(\mathbf{x})$ $\mathbf{x} \mapsto \mathbf{1}_{\partial\partial V}(\mathbf{x})$ are the indicator functions of the sets $\partial V$ and $\partial\partial V$. Then, we obtain

$$ \sum_{\alpha_1 \in A_1} \varphi_m(\mathbf{x}^{\alpha_1}, [\nabla_{\mathbf{x}} \mathbf{u}(\mathbf{x}^{\alpha_1})]) = \sum_{e \in \mathcal{E}} \left( \frac{1}{h_{V^e}^3} \int_{\Omega} f_m^e(\xi) \, d\xi + \int_{\partial\Omega} f_m^e(\xi) \left( \frac{\mathbf{1}_{\partial V}(\mathbf{s}^e(\xi))}{h_{S^e}^2} - \frac{1}{2h_{V^e}^2} \right) \, d\sigma \right. \\ \left. + \int_{\partial\partial\Omega} f_m^e(\xi) \left( \frac{\mathbf{1}_{\partial\partial V}(\mathbf{s}^e(\xi))}{h_{L^e}} - \frac{1}{2h_{S^e}} \right) \, d\ell \right). \quad (55) $$

*Proof.* By applying the Euler-Maclaurin summation formula, defined by Eq. (37), to the right member of Eqs. (47), (48), and (49), Proposition 3 is proven. $\square$

**Remark 6** The validity of the Euler-Maclaurin-based approximation in Proposition 3 requires a sufficiently large number of particles within each geometric entity (volume, surface, edge). However, this condition does not imply the need for global scale separation, as in classical homogenization theory. Rather, it ensures that discrete sums can be accurately approximated by integrals with bounded error. This is a convergence requirement, not a global scale separation assumption. Importantly, the model remains formulated at the same scale as the discrete system and does not rely on periodicity, statistical homogeneity, or the existence of a representative volume element (RVE). As such, it remains applicable even in the presence of non-periodic or irregular microstructures.

## 4.5 Weighting functions defining the coordinate transformation

We use the classical finite element method [59,60] for defining the coordinate transformation $\mathbf{s}^e$ based on the isoparametric weighting functions, which allows for writing,

$$ \mathbf{s}^e(\xi) = [N(\xi)] \mathbf{x}^e, \quad \mathbf{u}(\mathbf{s}^e(\xi)) = [N(\xi)] \mathbf{u}^e, \quad (56) $$

where $\mathbf{x}^e$ and $\mathbf{u}^e$ are the vectors of coordinates and displacements of all the nodes of the finite element $e$, where $[N(\xi)]$ is the matrix of the weighting functions expressed in terms of local coordinates $\xi$.

### 4.5.1 Case of an equivalent continuum to account for scale effects corresponding to linear elasticity

We consider the case where the term $c_{ij}(\mathbf{x}^{\alpha_1})[E(\mathbf{x}^{\alpha_1})]$ does not exist in Eq. (34) and $c_{ijkh}(\mathbf{x}^{\alpha_1})$ satisfies the usual properties of symmetry and positive definiteness. Therefore, referring to Section 3 and adding the hypothesis of small deformations, we have

$$ \varphi_m(\mathbf{x}^{\alpha_1}, [E(\mathbf{x}^{\alpha_1})]) \approx \varphi_m(\mathbf{x}^{\alpha_1}, [\varepsilon(\mathbf{x}^{\alpha_1})]) = \frac{1}{2} c_{ijkh}(\mathbf{x}^{\alpha_1})[\varepsilon(\mathbf{x}^{\alpha_1})]_{kh}[\varepsilon(\mathbf{x}^{\alpha_1})]_{ij}, \quad (57) $$
where $[\varepsilon(\mathbf{x})] = \frac{1}{2}(\nabla\mathbf{u}(\mathbf{x})^T + \nabla\mathbf{u}(\mathbf{x}))$ is the linearized strain tensor. Let $\varepsilon(\mathbf{x})$ denote the vector representation of $[\varepsilon(\mathbf{x})]$ in Voigt notation. We can write

$$ \varepsilon(\mathbf{s}^e(\boldsymbol{\xi})) = [\mathbf{B}(\boldsymbol{\xi})] \mathbf{u}^e, \qquad (58) $$

in which $[\mathbf{B}(\boldsymbol{\xi})]$ contains the first-order derivatives of weighting functions (as in [59]). Therefore, in Voigt notation and considering Eq. (58), Eq. (57) can be rewritten as

$$ \varphi_m(\mathbf{x}^{\alpha_1^e}, [\boldsymbol{\varepsilon}(\mathbf{x}^{\alpha_1^e})]) \approx \hat{\varphi}_m(\mathbf{x}^{\alpha_1^e}, \boldsymbol{\varepsilon}(\mathbf{x}^{\alpha_1^e})) = \frac{1}{2} \mathbf{u}^{eT} [\mathbf{B}(\boldsymbol{\xi}^{\alpha_1^e})]^T [\mathbf{D}^e(\boldsymbol{\xi}^{\alpha_1^e})] [\mathbf{B}(\boldsymbol{\xi}^{\alpha_1^e})] \mathbf{u}^e, \qquad (59) $$

where $[\mathbf{D}^e(\boldsymbol{\xi})]$ is the Voigt representation of the constitutive tensor $c_{ijkh}([\mathbf{N}(\boldsymbol{\xi})]\mathbf{x}^e)$, expressed in local coordinates. By summing with respect to $\alpha_1^e$, we can rewrite Eqs. (55) as follows

$$ \sum_{\alpha_1^e \in A_1^e} \hat{\varphi}_m(\mathbf{x}^{\alpha_1^e}, \boldsymbol{\varepsilon}(\mathbf{x}^{\alpha_1^e})) = \frac{1}{2} \mathbf{u}^{eT} [\mathbf{K}^e] \mathbf{u}^e, \qquad (60) $$

and $[\mathbf{K}^e]$ is the stiffness matrix defined by

$$ [\mathbf{K}^e] = \frac{1}{h_{V^e}^3} \int_{\Omega} [\mathbf{B}(\boldsymbol{\xi})]^T [\mathbf{D}^e(\boldsymbol{\xi})] [\mathbf{B}(\boldsymbol{\xi})] d\boldsymbol{\xi} + \int_{\partial\Omega} [\mathbf{B}(\boldsymbol{\xi})]^T [\mathbf{D}^e(\boldsymbol{\xi})] [\mathbf{B}(\boldsymbol{\xi})] \left( \frac{1\partial V([\mathbf{N}(\boldsymbol{\xi})]\mathbf{x}^e)}{h_{S^e}^2} - \frac{1}{2h_{V^e}^2} \right) d\sigma \\ + \int_{\partial\partial\Omega} [\mathbf{B}(\boldsymbol{\xi})]^T [\mathbf{D}^e(\boldsymbol{\xi})] [\mathbf{B}(\boldsymbol{\xi})] \left( \frac{1\partial V([\mathbf{N}(\boldsymbol{\xi})]\mathbf{x}^e)}{h_{L^e}} - \frac{1\partial V([\mathbf{N}(\boldsymbol{\xi})]\mathbf{x}^e)}{2h_{S^e}} \right) d\ell. \qquad (61) $$

The positivity of $[\mathbf{K}^e]$ is guaranteed for $h_{V^e}, h_{S^e}$, and $h_{L^e}$ that are $\ll 1$.

### Remark 7

(i) The technique of isoparametric finite elements is well-established in the literature and supported by theorems that we will not show here, referring instead to [59,60]. It should be noted that in the finite element stiffness matrix $[\mathbf{K}^e]$ in Eq. (61) a surface integral appears to account for scale effects for an arbitrary particle system.

(ii) The stiffness matrix follows directly from the particle system description.

(iii) Note that matrix $[\mathbf{D}^e(\boldsymbol{\xi})]$ is related to the constitutive property of the particle system not per unit volume as usually considered in the common finite element method in continuum mechanics.

### 4.5.2 Case of an equivalent continuum to account for scale effects in dynamics

When a deformable body is in motion, the displacement field $\mathbf{u}(\mathbf{x}, t)$ becomes a function of both space and time. In the quasi-static case, discussed earlier, dynamic effects were neglected. Within a dynamic framework, the kinetic energy must be taken into account. For each particle at $\mathbf{x}$, let us consider the kinetic energy $\varphi^{\text{kin}}$ expressed as,

$$ \varphi^{\text{kin}}(\mathbf{x}, \dot{\mathbf{u}}(\mathbf{x}, t)) = \frac{1}{2} m(\mathbf{x}) \langle \dot{\mathbf{u}}(\mathbf{x}, t), \dot{\mathbf{u}}(\mathbf{x}, t) \rangle, \qquad (62) $$

then, using the coordinate transformation defined in Eq. (56), summing with respect to $\alpha_1^e$, and applying the E-M summation formula as in Subsection 4.5.1, we obtain

$$ \sum_{\alpha_1^e \in A_1^e} \varphi^{\text{kin}}(\mathbf{x}^{\alpha_1^e}, \dot{\mathbf{u}}(\mathbf{x}^{\alpha_1^e}, t)) = \frac{1}{2} \dot{\mathbf{u}}^{eT}(t) [\mathbf{M}^e] \dot{\mathbf{u}}^{e}(t), \qquad (63) $$

in which $[\mathbf{M}^e]$ is the mass matrix defined by

$$ [\mathbf{M}^e] = \frac{1}{h_{V^e}^3} \int_{\Omega} \mu(\boldsymbol{\xi}) [\mathbf{N}(\boldsymbol{\xi})]^T [\mathbf{N}(\boldsymbol{\xi})] d\boldsymbol{\xi} + \int_{\partial\Omega} \mu(\boldsymbol{\xi}) [\mathbf{N}(\boldsymbol{\xi})]^T [\mathbf{N}(\boldsymbol{\xi})] \left( \frac{1\partial V([\mathbf{N}(\boldsymbol{\xi})]\mathbf{x}^e)}{h_{S^e}^2} - \frac{1}{2h_{V^e}^2} \right) d\sigma \\ + \int_{\partial\partial\Omega} \mu(\boldsymbol{\xi}) [\mathbf{N}(\boldsymbol{\xi})]^T [\mathbf{N}(\boldsymbol{\xi})] \left( \frac{1\partial V([\mathbf{N}(\boldsymbol{\xi})]\mathbf{x}^e)}{h_{L^e}} - \frac{1\partial V([\mathbf{N}(\boldsymbol{\xi})]\mathbf{x}^e)}{2h_{S^e}} \right) d\ell, \qquad (64) $$

where $\mu(\boldsymbol{\xi}) = m([\mathbf{N}(\boldsymbol{\xi})]\mathbf{x}^e)$ is the mass (not per unit volume) of the particle at $\mathbf{x}$ in local coordinates. The positivity of mass matrix $[\mathbf{M}^e]$ is guaranteed for $h_{V^e}, h_{S^e}$, and $h_{L^e}$ that are $\ll 1$. Mass matrix $[\mathbf{M}^e]$ defined by Eq. (64), accounts for surface effects driven by the particle system.
# 5 Discussion on the model parameters identification in a probabilistic framework

The discrete models and, consequently, the equivalent continuum model built upon them are characterized by geometric and constitutive parameters. Since nanoscale and microscale systems are inevitably affected by uncertainties, these parameters must be modeled by random fields and their identification must be conducted within a statistical framework. A probabilistic model for these parameters is therefore introduced, with a probability measure that depends on hyperparameters that must also be identified. This identification process corresponds to a statistical inverse problem, which requires target values associated with a given set of observations. These observations cannot be defined for each particle individually, but instead correspond to time-dependent displacements at various points in the equivalent continuum.

## 5.1 Formulation of an optimization problem for identifying the hyperparameters of the probabilistic model

(i) *Parametric probabilistic formulation.* Let $\mathbf{k}$ denote the vector of parameters on which the equivalent continuum model depends. For real-world applications, and assuming constitutive equations without memory, this vector parameter $\mathbf{k}$ is modeled as a random field,

$$ \{\mathbf{K}(\mathbf{x}; \delta) : \mathbf{x} \in \bar{\mathbf{V}}\}, \qquad (65) $$

defined on a probability space $(\Theta, \mathcal{F}, \mathbb{P})$. The random field $\mathbf{K}$ depends on a hyperparameter $\delta \in C_{\delta} \subset \mathbb{R}^{n_{\delta}}$ that governs its prior probability measure. We consider a random set

$$ \{\mathbf{W}(\mathbf{x}_i, t_j) \mid i = 1, \dots, n_x; j = 1, \dots, n_t\} \qquad (66) $$

of observations of the random system, where $\mathbf{W}(\mathbf{x}_i, t_j)$ is an $\mathbb{R}^{n_{\mathbf{w}}}$-valued random variable. Let

$$ \mathcal{E}^{\text{targ}} = \{ \mathbf{w}^{\text{targ},\ell}(\mathbf{x}_i, t_j) \mid i = 1, \dots, n_x; j = 1, \dots, n_t; \ell = 1, \dots, n_{\text{targ}} \} \qquad (67) $$

denote the targets associated with the defined observations. The optimal value $\delta^{\text{opt}}$ of $\delta$ is obtained by solving the optimization problem:

$$ \delta^{\text{opt}} = \arg \min_{\delta \in C_{\delta}} \mathbb{E} \{ \mathcal{J}(\delta; \mathcal{E}^{\text{targ}}) \}, \qquad (68) $$

where $\mathcal{J}(\delta; \mathcal{E}^{\text{targ}})$ is the random cost function.

(ii) *Nonparametric probabilistic formulation.* A nonparametric identification problem can also be performed as briefly presented in the next three subsections.

## 5.2 Brief comments concerning the classical methods to solve statistical inverse problems

In stochastic nonlinear computational models, the identification of parameters and hyperparameters can be performed using the following classical methods, within parametric and nonparametric frameworks. The Bayesian approach provides a systematic methodology to combine prior knowledge with target data [62–64]. When direct model evaluations are computationally prohibitive, statistical surrogate model can be constructed. For instance, the Gaussian Process regression offers a surrogate modeling technique that enables exploration of posterior distributions [65,66]. In parallel, information-theoretic and nonparametric probabilistic approaches [67] have provided robust frameworks for representing model-form uncertainty without assuming specific parametric distributions, particularly through the use of maximum entropy principles and random matrix theory [68,69], and extend to nonlinear stochastic dynamical systems [70]. These methods have proven effective for high-dimensional problems involving nonlinear dynamics and non-Gaussian behavior. Sequential inference in time-dependent systems is often handled [71].

## 5.3 Physics-Informed Neural Networks for nonlinear stochastic partial differential equations

Physics-Informed Neural Networks (PINNs) have emerged as a tool for solving partial differential equations by incorporating physical laws into the training of deep neural networks.

While effective for deterministic systems, extending PINNs to nonlinear stochastic PDEs (SPDEs) remains challenging due to issues of computational complexity, closure of statistical moments, and the inadequacy of generative models for realistic high-dimensional physics. In the deterministic setting, Physics-Informed Neural Networks (PINNs) [72,73] provide a framework for solving dynamical systems governed by ordinary or partial
differential equations. The key idea is to approximate the unknown solution using a deep neural network. Rather than relying solely on empirical data, PINNs incorporate the governing differential equations directly into the learning process by constructing a loss function that combines a *data loss*, ensuring consistency with available observations, and a *physics loss*, penalizing residuals of the differential equation. The total loss is typically a weighted sum of these components. By minimizing this objective, the neural network learns a solution that satisfies both the physical model and the available data.

In the probabilistic framework, extending PINNs to nonlinear stochastic partial differential equations (SPDEs) poses significant challenges. Several strategies have been proposed in recent literature. (i) Moment-based PINNs attempt to learn statistical moments of the stochastic solution [74]. However, for nonlinear SPDEs, the evolution of moments generates an infinite hierarchy of coupled equations. Closing this hierarchy requires strong assumptions (e.g., moment truncation, Gaussian closure), which limits generality and accuracy. (ii) Monte Carlo PINNs, where a separate neural network is trained for each stochastic realization [75], are conceptually straightforward but incur exponential computational costs in high-dimensional random spaces. Moreover, treating random network weights is nontrivial within standard architectures. (iii) Generative PINNs, including those based on variational autoencoders (VAE-PINNs) or generative adversarial networks (GAN-PINNs) [76, 77], aim to model the full probability distribution of the solution field. However, such methods present difficulties in complex physical systems. While PINNs demonstrate significant potential for solving deterministic nonlinear systems, their application to nonlinear SPDEs remains an open research challenge. Existing approaches face limitations in terms of scalability, accuracy, and robustness in realistic, high-dimensional, and physically complex settings.

## 5.4 Probabilistic Learning on Manifolds (PLoM) constrained by nonlinear stochastic partial differential equations for small datasets

In contrast to neural-network-based methods such as Physics-Informed Neural Networks (PINNs), which typically require large datasets and complex training pipelines, a complementary approach has been developed for small training datasets. The Probabilistic Learning on Manifolds (PLoM) method [78, 79], extended to incorporate constraints from nonlinear stochastic partial differential equations (PDEs) [80], constructs a statistical surrogate model that remains consistent with the governing physical laws. Rather than approximating the PDE solution directly via neural architectures, the PLoM approach builds a probabilistic representation of the solution manifold from a small training dataset. New realizations are generated by preserving the structure of the estimated probability measure on the intrinsic manifold, which is learned using diffusion maps. In its constrained form, PLoM further ensures that the generated samples satisfy the stochastic PDE in a mean-square sense, providing a physics-informed statistical generator suitable for complex and high-dimensional systems. This approach is particularly well-suited to scenarios where data is scarce but simulations are computationally expensive, and it offers a viable alternative to deep learning methods for uncertainty quantification and stochastic modeling in nonlinear systems.

## 6 Conclusions

In this paper, we present a novel general methodology for transforming a discrete particle system occupying an arbitrary 3D nanoscale domain, characterized by an *m*-interaction potential, into an equivalent continuum model at the same scale, which captures scale effects. We obtain a general expression for the nonlocal *m*-interaction potential that satisfies the objectivity principle and point-symmetry invariance. We introduce a local representation of this nonlocal *m*-interaction potential using Taylor expansions. Finally, we construct the equivalent continuum model at the same scale by applying the Euler-Maclaurin (E-M) summation formula. As part of this work, we propose a generalization for an arbitrary domain based on coordinate transformations. For such an arbitrary domain, the proposed approach provides a pathway for integration into a computational framework. While the proposed model is not intended as a replacement for molecular dynamics or quasi-continuum approaches, it offers a computationally efficient and analytically grounded alternative for systems affected by volume, surface, edges size effects in statics and dynamics. Since nanoscale and microscale systems are inevitably affected by uncertainties, the geometric and constitutive parameters should be modeled as random fields and the procedure for their identification has been discussed within a statistical framework. Future work will focus on the numerical implementation of the equivalent continuum model to investigate its capability in capturing size effects by comparison with discrete simulations. Moreover, the development and integration of the stochastic formulation will allow for a systematic quantification of uncertainties in geometric and constitutive parameters, enabling a probabilistic validation of the model predictions and enhancing its robustness and applicability in nanoscale mechanics.
## Declarations

**Ethical approval:** not applicable. **Consent for publication:** the authors give their consent for publication. **Availability of data and materials:** not applicable. **Competing interests:** the authors declare that they have no interests of a financial or personal nature that might be perceived to influence the results reported in this paper. **Funding:** not applicable. **Authors' contributions:** G.L.V. developed the theory, wrote, and reviewed the manuscript, C.S. developed the theory, wrote, and reviewed the manuscript.

## References

1. P. Jiao, J. Mueller, J. R. Raney, X. Zheng, and A. H. Alavi, "Mechanical metamaterials and beyond," *Nature communications*, vol. 14, no. 1, p. 6004, 2023. DOI: https://doi.org/10.1038/s41467-023-41679-8.

2. J. Fish, "Bridging the scales in nano engineering and science," *Journal of Nanoparticle Research*, vol. 8, pp. 577-594, 2006. DOI: https://doi.org/10.1007/s11051-006-9090-9.

3. C. Soize, "Non-Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators," *Computer Methods in Applied Mechanics and Engineering*, vol. 195, no. 1, pp. 26-64, 2006. DOI: https://doi.org/10.1016/j.cma.2004.12.014.

4. J. Fish, *Practical multiscaling*. John Wiley & Sons, 2013.

5. A. Malyarenko and M. Ostoja-Starzewski, "Statistically isotropic tensor random fields: correlation structures," *Mathematics and Mechanics of Complex Systems*, vol. 2, no. 2, pp. 209-231, 2014. DOI: dx.doi.org/10.2140/memocs.2014.2.209.

6. G. J. Fish, G. J. Wagner, and S. Keten, "Mesoscopic and multiscale modelling in materials," *Nature materials*, vol. 20, no. 6, pp. 774-786, 2021. DOI: https://doi.org/10.1038/s41563-020-00913-0.

7. G. C. Papanicolaou, "Boundary value problems with rapidly oscillating random coefficients," in *Colloquia Math. Soc., Janos Bolyai*, vol. 27, pp. 853-873, 1979.

8. S. Torquato and G. Stell, "Microstructure of two-phase random media. V. The n-point matrix probability functions for impenetrable spheres," *The Journal of chemical physics*, vol. 82, no. 2, pp. 980-987, 1985. DOI: https://doi.org/10.1063/1.448475.

9. G. Nguetseng, "A general convergence result for a functional related to the theory of homogenization," *SIAM Journal on Mathematical Analysis*, vol. 20, no. 3, pp. 608-623, 1989. DOI: https://doi.org/10.1137/0520043.

10. G. Allaire, "Homogenization and two-scale convergence," *SIAM Journal on Mathematical Analysis*, vol. 23, no. 6, pp. 1482-1518, 1992. DOI: https://doi.org/10.1137/0523084.

11. K. Sab, "On the homogenization and the simulation of random materials," *Eur. J. Mech., A/Solids*, no. 5, pp. 585-607, 1992.

12. K. T. Andrews and S. Wright, "Stochastic homogenization of elliptic boundary-value problems with $L^p$-data," *Asymp-totic Analysis*, vol. 17, no. 3, pp. 165-184, 1998. DOI: https://doi.org/10.3233/ASY-1998-30.

13. M. Ostoja-Starzewski, "Random field models of heterogeneous materials," *International Journal of Solids and Structures*, vol. 35, no. 19, pp. 2429-2455, 1998. DOI: https://doi.org/10.1016/S0020-7683(97)00144-3.

14. V. V. Jikov, S. M. Kozlov, and O. A. Oleinik, *Homogenization of differential operators and integral functionals*. Springer Science & Business Media, 2012. DOI: https://doi.org/10.1007/978-3-642-84659-5.

15. C. Soize, "Stochastic elliptic operators defined by non-gaussian random fields with uncertain spectrum," *Theory of Probability and Mathematical Statistics*, vol. 105, pp. 113-136, 2021. DOI: https://doi.org/10.1090/tpms/1159.

16. M. Ostoja-Starzewski, "Material spatial randomness: From statistical to representative volume element," *Probabilistic engineering mechanics*, vol. 21, no. 2, pp. 112-132, 2006. DOI: https://doi.org/10.1016/j probengmech.2005.07.007.

17. M. Ostoja-Starzewski, X. Du, Z. Khisaeva, and W. Li, "Comparisons of the size of the representative volume element in elastic, plastic, thermoelastic, and permeable random microstructures," *International Journal for Multiscale Computational Engineering*, vol. 5, no. 2, 2007. DOI: https://doi.org/10.1615/IntJMultCompEng.v5.i2.10.

18. C. Soize, "Tensor-valued random fields for meso-scale stochastic model of anisotropic elastic microstructure and probabilistic analysis of representative volume element size," *Probabilistic Engineering Mechanics*, vol. 23, no. 2-3, pp. 307-323, 2008. DOI: https://doi.org/10.1016/j.probengmech.2007.12.019.

19. T. Zhang, F. Pled, and C. Desceliers, "Robust multiscale identification of apparent elastic properties at mesoscale for random heterogeneous materials with multiscale field measurements," *Materials*, vol. 13, no. 12, p. 2826, 2020. DOI: https://doi.org/10.3390/ma13122826.

20. F. Dell'Isola, G. Maier, U. Perego, U. Andreaus, R. Esposito, and S. Forest, *The complete works of Gabrio Piola: volume I: commented English translation*, vol. 38. Springer, 2014. DOI: https://doi.org/10.1007/978-3-319-00263-7.

21. F. Dell'Isola, U. Andreaus, A. Cazzani, R. Esposito, L. Placidi, U. Perego, G. Maier, P. Seppecher, et al., *The Complete Works of Gabrio Piola: Volume II-Commented English Translation*. Springer International Publishing AG, 2019. DOI: https://doi.org/10.1007/978-3-319-70692-4.

22. E. Cosserat and F. Cosserat, *Theorie des Corps Deformables*. Hermann, Paris, 1970.

23. J.-J. Alibert, P. Seppecher, and F. Dell'Isola, "Truss modular beams with deformation energy depending on higher displacement gradients," *Mathematics and Mechanics of Solids*, vol. 8, no. 1, pp. 51-73, 2003. DOI: https://doi.org/10.1177/1081286503008001658.

24. H. Abdoul-Anziz and P. Seppecher, "Strain gradient and generalized continua obtained by homogenizing frame lattices," *Mathematics and mechanics of complex systems*, vol. 6, no. 3, pp. 213-250, 2018. DOI: dx.doi.org/10.2140/memocs.2018.6.213.

25. S. R. Eugster and E. Barchiesi, "A second gradient continuum formulation for bi-pantographic fabrics," *PAMM*, vol. 21, no. 1, p. e202100192, 2021. DOI: https://doi.org/10.1002/pamm.202100192.

26. A. Ciallella, F. D'Annibale, D. Del Vescovo, and I. Giorgio, "Deformation patterns in a second-gradient lattice annular plate composed of "spira mirabilis" fibers," *Continuum Mechanics and Thermodynamics*, vol. 35, no. 4, pp. 1561-1580, 2023. DOI: https://doi.org/10.1007/s00161-022-01169-6.

27. A. Ciallella and D. J. Steigmann, "Unusual deformation patterns in a second-gradient cylindrical lattice shell: Numerical experiments," *Mathematics and Mechanics of Solids*, vol. 28, no. 1, pp. 141-153, 2023. DOI: https://doi.org/10.1177/10812865221101820.
28. G. La Valle and C. Soize, "A higher-order nonlocal elasticity continuum model for deterministic and stochastic particle-based materials," *Zeitschrift für angewandte Mathematik und Physik*, vol. 75, no. 2, p. 49, 2024. DOI: https://doi.org/10.1007/s00033-024-02196-w.

29. G. La Valle and C. Soize, "Stochastic second-gradient continuum theory for particle-based materials: part II," *Zeitschrift für angewandte Mathematik und Physik*, vol. 75, no. 3, p. 93, 2024. DOI: https://doi.org/10.1007/s00033-024-02232-9.

30. G. La Valle and C. Soize, "Identifying second-gradient continuum models in particle-based materials with pairwise interactions using acoustic tensor methodology," *Journal of Elasticity*, pp. 1-17, 2024. DOI: https://doi.org/10.1007/s10659-024-10067-8.

31. I. A. Kunin, *Elastic Media with Microstructure I: One-Dimensional Models*, vol. 26. Springer-Verlag, Berlin, Heidelberg, New York, Tokyo, 1982. https://doi.org/10.1007/978-3-642-81748-9.

32. I. A. Kunin, *Elastic Media with Microstructure II: Three-Dimensional Models*, vol. 44. Springer-Verlag, Berlin, Heidelberg, New York, Tokyo, 1983. https://doi.org/10.1007/978-3-642-81960-5.

33. V. A. Eremeyev, L. P. Lebedev, and H. Altenbach, *Foundations of Micropolar Mechanics*. Springer Science & Business Media, 2012. DOI: https://doi.org/10.1007/978-3-642-28353-6.

34. P. Germain, "The method of virtual power in continuum mechanics. Part 2: Microstructure," *SIAM Journal on Applied Mathematics*, vol. 25, no. 3, pp. 556-575, 1973. DOI: https://doi.org/10.1137/0125053.

35. S. Forest and K. Sab, "Finite-deformation second-order micromorphic theory and its relations to strain and stress gradient models," *Mathematics and Mechanics of Solids*, vol. 25, no. 7, pp. 1429-1449, 2020. DOI: https://doi.org/10.1177/1081286517720844.

36. S. Forest, "Continuum thermomechanics of nonlinear micromorphic, strain and stress gradient media," *Philosophical Transactions of the Royal Society A*, vol. 378, no. 2170, p. 20190169, 2020. DOI: https://doi.org/10.1098/rsta.2019.0169.

37. E. Kröner, "Elasticity theory of materials with long range cohesive forces," *International Journal of Solids and Structures*, vol. 3, no. 5, pp. 731-742, 1967. DOI: https://doi.org/10.1016/0020-7683(67)90049-2.

38. A. C. Eringen, "Linear theory of nonlocal elasticity and dispersion of plane waves," *International Journal of Engineering Science*, vol. 10, no. 5, pp. 425-435, 1972. DOI: https://doi.org/10.1016/0020-7225(72)90050-X.

39. A. C. Eringen and D. G. B. Edelen, "On nonlocal elasticity," *International Journal of Engineering Science*, vol. 10, no. 3, pp. 233-248, 1972. DOI: https://doi.org/10.1016/0020-7225(72)90039-0.

40. A. C. Eringen, "Vistas of nonlocal continuum physics," *International Journal of Engineering Science*, vol. 30, no. 10, pp. 1551-1565, 1992. DOI: https://doi.org/10.1016/0020-7225(92)90165-D.

41. A. C. Eringen, *Nonlocal Continuum Field Theories*. New York: Springer, 2002. DOI: https://doi.org/10.1007/b97697.

42. C. Polizzotto, P. Fuschi, and A. A. Pisano, "A strain-difference-based nonlocal elasticity model," *International Journal of Solids and Structures*, vol. 41, no. 9, pp. 2383-2401, 2004. DOI: https://doi.org/10.1016/j.ijsolstr.2003.12.013.

43. L. Placidi, E. Barchiesi, A. Misra, and D. Timofeev, "Micromechanics-based elasto-plastic-damage energy formulation for strain gradient solids with granular microstructure," *Continuum Mechanics and Thermodynamics*, vol. 33, no. 5, pp. 2213-2241, 2021. DOI: https://doi.org/10.1007/s00161-021-01023-1.

44. E. B. Tadmor, M. Ortiz, and R. Phillips, "Quasicontinuum analysis of defects in solids," *Philosophical magazine A*, vol. 73, no. 6, pp. 1529-1563, 1996. DOI: https://doi.org/10.1080/01418619608243000.

45. E. B. Tadmor and R. E. Miller, *Modeling materials: Continuum, Atomistic and Multiscale Techniques*. Cambridge University Press, 2011. DOI: https://doi.org/10.1017/CBO9781139003582.

46. G. J. Wagner and W. K. Liu, "Coupling of atomistic and continuum simulations using a bridging scale decomposition," *Journal of Computational Physics*, vol. 190, no. 1, pp. 249-274, 2003. DOI: https://doi.org/10.1016/S0021-9991(03)00273-0.

47. H. S. Park, E. G. Karpov, W. K. Liu, and P. A. Klein, "The bridging scale for two-dimensional atomistic/continuum coupling," *Philosophical magazine*, vol. 85, no. 1, pp. 79-113, 2005. DOI: https://doi.org/10.1080/14786430412331300163.

48. W. K. Liu, H. S. Park, D. Qian, E. G. Karpov, H. Kadowaki, and G. J. Wagner, "Bridging scale methods for nanomechanics and materials," *Computer Methods in Applied Mechanics and Engineering*, vol. 195, no. 13-16, pp. 1407-1421, 2006. DOI: https://doi.org/10.1016/j.cma.2005.05.042.

49. C. Soize, "Probabilistic learning inference of boundary value problem with uncertainties based on Kullback-Leibler divergence under implicit constraints," *Computer Methods in Applied Mechanics and Engineering*, vol. 395, p. 115078, 2022. DOI: https://doi.org/10.1016/j.cma.2022.115078.

50. B. Nayroles, G. Touzot, and P. Villon, "Generalizing the finite element method: Diffuse approximation and diffuse elements," *Computational mechanics*, vol. 10, no. 5, pp. 307-318, 1992. DOI: https://doi.org/10.1007/BF00364252.

51. T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, and P. Krysl, "Meshless methods: An overview and recent developments," *Computer methods in applied mechanics and engineering*, vol. 139, no. 1-4, pp. 3-47, 1996. DOI: https://doi.org/10.1016/S0045-7825(96)01078-X.

52. C. A. Duarte and J. T. Oden, "H-p clouds—an h-p meshless method," *Numerical Methods for Partial Differential Equations: An International Journal*, vol. 12, no. 6, pp. 673-705, 1996. DOI: https://doi.org/10.1002/(SICI)1098-2426(199611)12:6<673::AID-NUM3>3.0.CO;2-P.

53. P. Breitkopf, A. Rassineux, G. Touzot, and P. Villon, "Explicit form and efficient computation of MLS shape functions and their derivatives," *International Journal for Numerical Methods in Engineering*, vol. 48, no. 3, pp. 451-466, 2000. DOI: https://doi.org/10.1002/(SICI)1097-0207(20000530)48:3<451::AID-NME892>3.0.CO;2-1.

54. A. Rassineux, P. Villon, J.-M. Savignat, and O. Stab, "Surface remeshing by local hermite diffuse interpolation," *International Journal for numerical methods in Engineering*, vol. 49, no. 1-2, pp. 31-49, 2000. DOI: https://doi.org/10.1002/1097-0207(20000910/20)49:1/2<31::AID-NME921>3.0.CO;2-6.

55. X. Zhang, K. Z. Song, M. W. Lu, and X. Liu, "Meshless methods based on collocation with radial basis functions," *Computational mechanics*, vol. 26, pp. 333-343, 2000. DOI: https://doi.org/10.1007/s004600000181.

56. A. Romano, *Lezioni di meccanica razionale*. Liguori, 1977.

57. W. K. Liu, E. G. Karpov, and H. S. Park, *Nano Mechanics and Materials: Theory, Multiscale Methods and Applications*. John Wiley & Sons, 2006. DOI: https://doi.org/10.1002/0470034106.

58. F. dell'Isola, U. Andreaus, and L. Placidi, "At the origins and in the vanguard of peridynamics, non-local and higher-gradient continuum mechanics: An underestimated and still topical contribution of Gabrio Piola," *Mathematics and Mechanics of Solids*, vol. 20, no. 8, pp. 887-928, 2015. DOI: https://doi.org/10.1177/1081286513509811.

59. O. C. Zienkiewicz, *La méthode des éléments finis: appliquées à l'art de l'ingénieur*. Ediscience, Paris, France, 1973.

60. R. Dautray and J.-L. Lions, *Mathematical analysis and numerical methods for science and technology*. Springer Science & Business Media, 2012.
61. M. Abramowitz and I. A. Stegun, *Handbook of mathematical functions with formulas, graphs, and mathematical tables*, vol. 55. US Government printing office, 1948.

62. J. P. Kaipio and E. Somersalo, *Statistical and computational inverse problems*. Springer, 2005. DOI: https://doi.org/10.1007/b138659.

63. A. M. Stuart, "Inverse problems: a bayesian perspective," *Acta numerica*, vol. 19, pp. 451-559, 2010. DOI: https://doi.org/10.1017/S0962492910000061.

64. A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin, *Bayesian Data Analysis*. CRC Press, 3rd ed., 2013. DOI: https://doi.org/10.1201/b16018.

65. C. E. Rasmussen and C. K. I. Williams, *Gaussian Processes for Machine Learning*. The MIT Press, 11 2005. DOI: 10.7551/mitpress/3206.001.0001.

66. M. C. Kennedy and A. O'Hagan, "Bayesian calibration of computer models," *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, vol. 63, no. 3, pp. 425-464, 2001. DOI: https://doi.org/10.1111/1467-9868.00294.

67. C. Soize, *Uncertainty quantification*. Springer, 2017. DOI: https://doi.org/10.1007/978-3-319-54339-0.

68. C. Soize, "Construction of probability distributions in high dimension using the maximum entropy principle: Applications to stochastic processes, random fields and random matrices," *International Journal for Numerical Methods in Engineering*, vol. 76, no. 10, pp. 1583-1611, 2008. DOI: https://doi.org/10.1002/nme.2385.

69. C. Soize, "Random matrix models and nonparametric method for uncertainty quantification," in *Handbook of uncertainty quantification*, pp. 1-69, Springer, 2015. DOI: https://doi.org/10.1007/978-3-319-11259-6_5_1.

70. C. Soize and C. Farhat, "A nonparametric probabilistic approach for quantifying uncertainties in low-dimensional and high-dimensional nonlinear models," *International Journal for Numerical methods in engineering*, vol. 109, no. 6, pp. 837-888, 2017. DOI: https://doi.org/10.1002/nme.5312.

71. A. Doucet, N. De Freitas, N. J. Gordon, et al., *Sequential Monte Carlo methods in practice*. Springer, 2001. DOI: https://doi.org/10.1007/978-1-4757-3437-9.

72. M. Raissi, P. Perdikaris, and G. E. Karniadakis, "Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations," *Journal of Computational physics*, vol. 378, pp. 686-707, 2019. DOI: https://doi.org/10.1016/j.jcp.2018.10.045.

73. B. Yuan, C. S. Choo, L. Y. Yeo, Y. Wang, Z. Yang, Q. Guan, S. Suryasentana, J. Choo, H. Shen, M. Megia, et al., "Physics-informed machine learning in geotechnical engineering: a direction paper," *Geomechanics and Geoengineering*, pp. 1-32, 2025. DOI: https://doi.org/10.1080/17486025.2025.2502029.

74. Z. Chen, Y. Liu, and H. Sun, "Physics-informed learning of governing equations from scarce data," *Nature communications*, vol. 12, no. 1, p. 6136, 2021. DOI: https://doi.org/10.1038/s41467-021-26434-1.

75. L. Yang, X. Meng, and G. E. Karniadakis, "B-PINNs: Bayesian physics-informed neural networks for forward and inverse PDE problems with noisy data," *Journal of Computational Physics*, vol. 425, p. 109913, 2021. DOI: https://doi.org/10.1016/j.jcp.2020.109913.

76. T.-Y. Yang, J. Rosca, K. Narasimhan, and P. J. Ramadge, "Learning physics constrained dynamics using autoencoders," in *Advances in Neural Information Processing Systems* (S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, eds.), vol. 35, pp. 17157-17172, Curran Associates, Inc., 2022.

77. W. Zhong and H. Meidani, "Pi-vae: Physics-informed variational auto-encoder for stochastic differential equations," *Computer Methods in Applied Mechanics and Engineering*, vol. 403, p. 115664, 2023. DOI: https://doi.org/10.1016/j.cma.2022.115664.

78. C. Soize and R. Ghanem, "Data-driven probability concentration and sampling on manifold," *Journal of Computational Physics*, vol. 321, pp. 242-258, 2016. https://doi.org/10.1016/j.jcp.2016.05.044.

79. C. Soize and R. Ghanem, "Transient anisotropic kernel for probabilistic learning on manifolds," *Computer Methods in Applied Mechanics and Engineering*, vol. 432, p. 117453, 2024. DOI: https://doi.org/10.1016/j.cma.2024.117453.

80. C. Soize and R. Ghanem, "Probabilistic learning on manifolds constrained by nonlinear partial differential equations for small datasets," *Computer Methods in Applied Mechanics and Engineering*, vol. 380, p. 113777, 2021. DOI: https://doi.org/10.1016/j.cma.2021.113777.