A Priori Analysis of a Higher Order Nonlinear Elasticity Model for an Atomistic Chain with Periodic Boundary Condition

TL;DR

This paper rigorously derives a higher order nonlinear elasticity model for one-dimensional crystals from atomistic descriptions, achieving fourth order accuracy and demonstrating convergence through numerical experiments.

Contribution

It introduces a novel higher order continuum model with fourth order accuracy derived from atomistic physics, extending beyond the traditional second order Cauchy-Born model.

Findings

Higher order model achieves fourth order accuracy

Numerical experiments confirm theoretical convergence

Model applicable to defect dynamics and nanotubes

Abstract

Nonlinear elastic models are widely used to describe the elastic response of crystalline solids, for example, the well-known Cauchy-Born model. While the Cauchy-Born model only depends on the strain, effects of higher order strain gradients are significant and higher order continuum models are preferred, in various applications such as defect dynamics and modeling of carbon nanotubes. In this paper, we rigorously derive a higher order nonlinear elasticity model for crystals from its atomistic description in one dimension. We show that, compared to the second order accuracy of the Cauchy-Born model, the higher order continuum model in this paper is of fourth oder accuracy with respect to the interatomic spacing in the thermal dynamic limit. In addition, we discuss the key issues for the derivation of higher order continuum models in more general cases. The theoretical convergence results are demonstrated by numerical experiments.

Equations

$$\Vert \nabla u^{\mathrm{a}}-\nabla u^{\mathrm{hoc}}{\Vert }_{L^2}\le C{\epsilon }^4,$$

$$⟨\delta \mathcal{F}(u),v⟩:=\underset{t\to 0}{\lim }{t}^{-1}(\mathcal{F}(u+tv)-\mathcal{F}(u))$$

$$\text{and}⟨{\delta }^2\mathcal{F}(u)v,w⟩:=\underset{t\to 0}{\lim }{t}^{-1}(\delta \mathcal{F}(u+tw)-\delta \mathcal{F}(u),v).$$

$$\mathcal{U}:=\left\{v:\Lambda \to \mathbb{R}\right\}.$$

$$D_{\rho }v(\xi ):=v(\xi +\rho )-v(\xi ),\text{for}\xi \in \Lambda ,\rho \in \mathcal{R},$$

$$\Pi v(\xi )$$

$$\nabla \Pi v(\xi )$$

$${\nabla }^2\Pi v(\xi )$$

$${\nabla }^3\Pi v(\xi )$$

$${\nabla }^4\Pi v(\xi )$$

$$\sum _{\xi \in \Lambda }(a+b\xi +c{\xi }^2+d{\xi }^3)\zeta (x-\xi )=a+bx+c{x}^2+d{x}^3,\forall a,b,c,d\in \mathbb{R}.$$

$$\hat{v}(x):=\sum _{\xi \in \Lambda }v(\xi )\zeta (x-\xi ),\text{for}v\in \mathcal{U}.$$

$$\tilde{v}(x):=(\zeta \ast \hat{v})(x)={\int }_{\mathbb{R}}\zeta (x-y)\hat{v}(y)dy.$$

$$D_{\rho }\tilde{v}(\xi )$$

$$={\int }_{\mathbb{R}}{\int }_0^1\zeta (\xi +t\rho -x)dt\cdot {\nabla }_{\rho }\hat{v}(x)dx$$

$$=:{\int }_{\mathbb{R}}{\chi }_{\xi ,\rho }(x)\cdot {\nabla }_{\rho }\hat{v}(x)dx.$$

$$||\nabla \hat{v}|{|}_{L^2}\lesssim ||\nabla \Pi v|{|}_{L^2}\lesssim ||\nabla \hat{v}|{|}_{L^2},\forall v\in \mathcal{U}.$$

$$||\nabla \hat{v}|{|}_{L^2}\lesssim ||\nabla \tilde{v}|{|}_{L^2}\lesssim ||\nabla \hat{v}|{|}_{L^2},\forall v\in \mathcal{U}.$$

$$\mathscr{U}^{{\rm{per}}}=\big{\{}u\in\mathscr{U}:u_{\xi+2N}=u_{\xi},{\textstyle\sum_{\xi\in\Omega^{\Lambda}}}u_{\xi}=0\big{\}}.$$

$$\mathscr{Y}^{{\rm{per}}}=\big{\{}y\in\mathbb{R}^{\mathbb{Z}}:y_{\xi}=F\xi+u_{\xi},u\in\mathscr{U}^{{\rm{per}}}\big{\}},$$

$${\mathcal{E}}_{\mathrm{a}}(y)=\sum _{\xi \in {\Omega }^{\Lambda }}\sum _{\rho \in \mathcal{R}}\varphi (D_{\rho }y(\xi )),$$

$${\mathcal{E}}_{\mathrm{a}}(u)=\sum _{\xi \in {\Omega }^{\Lambda }}\sum _{\rho \in \mathcal{R}}{\varphi }_{\rho }(D_{\rho }u(\xi )),$$

$$\displaystyle\mathscr{U}^{1,2}:=\{u\in\mathscr{U}^{{\rm{per}}}\big{|}||\nabla\hat{u}||_{L^{2}(\Omega)}<+\infty\}.$$

$$\displaystyle\mathscr{K}:=\{u\in\mathscr{U}^{1,2}\big{|}|D_{\rho}u(\xi)|\leq\kappa,\forall\xi\in\Omega^{\Lambda},\rho\in\mathscr{R}\},$$

$$M^{(j,s)}:=\sum _{\rho \in \mathcal{R}}{m}^{(j,s)}(\rho )<\infty ,$$

$$m^{(j,s)}(\rho ):={\rho }^{j+s}\underset{g}{\sup }|{\varphi }_{\rho }^{(j)}(g)|,$$

$${\mathcal{E}}_{\mathrm{a}}(u)=\sum _{\xi \in {\Omega }^{\Lambda }}\sum _{\rho \in \mathcal{R}}{\varphi }_{\rho }(D_{\rho }u(\xi )),$$

$$D_{\rho }u(\xi )=$$

$$=$$

$$\displaystyle-\Big{[}u(\xi^{\prime})-\frac{\rho}{2}\nabla u(\xi^{\prime})+\frac{1}{2}\nabla^{2}u(\xi^{\prime})(\frac{\rho}{2})^{2}-\frac{1}{6}\nabla^{3}u(\xi^{\prime})(\frac{\rho}{2})^{3}+...\Big{]}$$

$$\approx$$

$${\mathcal{E}}_{\mathrm{a}}(u)\approx \sum _{\xi \in {\Omega }^{\Lambda }}\sum _{\rho \in \mathcal{R}}{\varphi }_{\rho }(\rho \nabla u({\xi }^{\prime })+\frac{{\rho }^3}{24}{\nabla }^{3}u({\xi }^{\prime })).$$

$${\mathcal{E}}_{\mathrm{hoc}}(u)={\int }_{\Omega }\sum _{\rho \in \mathcal{R}}{\varphi }_{\rho }(\rho \nabla u+\frac{{\rho }^3}{24}{\nabla }^{3}u)dx,$$

$${\mathcal{E}}_{\mathrm{cb}}(u):={\int }_{\Omega }\sum _{\rho \in \mathcal{R}}{\varphi }_{\rho }(\rho \nabla u)dx.$$

$${\mathcal{E}}_{\mathrm{ill}}(u)={\int }_{\Omega }\frac{1}{2}(\nabla u{)}^2-\frac{1}{24}({\nabla }^{2}u{)}^{2}dx,$$

Full text

1. Introduction

Nonlinear elasticity models are widely used to describe the elastic response of crystalline materials. The Cauchy-Born model is probably the most well known nonlinear elasticity model which is consistent with the atomistic theory of crystals, and it is second order accurate with respect to the atomistic model under certain technical assumptions [5, 4, 8, 9, 20]. The Cauchy-Born energy density only depends on the strain, and can be interpreted as ’the stored energy per unit volume under a macroscopically homogeneous deformation equals the energy per unit volume in the corresponding homogeneous crystal’ [20].

The Cauchy-Born model is sufficiently accurate if the strain gradient is small. In various situations, nonlinear elastic models of higher order strain gradients are preferred. For example, the higher order strain gradients have significant impact for the defect zone [16], for curved crystalline sheets such as carbon nanotube [10], and for the wave propagation in crystals [3, 25].

In the mechanics literature, higher order continuum models were first derived in [28], where two different approaches were used to derive the continuum macro model from the discrete micro model. Later on, it was found that higher order continuum models can efficiently capture the inhomogeneous deformation of the underlying crystal [27] and the curvature effect of carbon nanotubes [26]. In [25], the higher order gradient model links the atomistic model and nonlocal models such as peridynamics, and is able to capture the correct dispersive behavior of crystals. Recently, [16] proposed a multiscale crystal defect dynamics (MCDD) model by adopting different higher order Cauchy-Born models (up to four) to construct atomistic informed constitutive relations for various defect process zones, and developed a hierachical strain gradient based finite element formulation.

On the contrary, only few works have been devoted to the mathematical analysis of higher order continuum models. In [4], it was proved that a higher order continuum model indeed has better accuracy in terms of energy. It was shown in [3] that the higher order continuum models in [28] might be ill-posed, and may lead to an uncontrolled behavior of the solution. The so-called ”inner expansion”, which is based on the formal Taylor expansion of the deformation gradient at some carefully chosen expansion points, was proposed to derive continuum models from the atomistic models (with pair interactions) and minimize the remainder terms of the energy. While a well-posed higher order continuum approximation was first developed in [3], a rigorous mathematical analysis was not included.

The main purpose of the current work is to derive a higher order nonlinear elasticity model from the atomistic model (with pair potential) in one dimension and present a rigorous a priori analysis of the obtained higher order model. The derivation of the model essentially follows the techniques of the ’inner expansion’ introduced in [3] which leads to a well-posed higher order continuum model.

The major contribution of the current work is that, to the best knowledge of the authors, it for the first time provides a rigorous analysis and error estimate for the energy minimizer of the higher order continuum model and numerically demonstrate the fourth order accuracy for such a model. To be precise, we will show that the approximation error, which is also known as the modeling error, is of $O(\varepsilon^{4})$. Namely, if we scale the system by in interatomic spacing $\varepsilon$, we have

[TABLE]

where $u^{\rm{a}}$ and $u^{\rm{hoc}}$ are the solutions to the atomistic model and the higher order continuum model, respectively, which will be defined in Section 3.3, and the constant $C$ depends only on some higher order partial derivatives of the interatomic potential $\phi$ and on the regularity of $u^{\rm{a}}$. We essentially extend the analytical framework in [20] to the higher order continuum model in one dimension with pair interactions, which include the analysis of the modeling error, stability and convergence estimates. In addition, we point out the possibilities and challenges to extend those results to the physically more relevant cases of multibody interactions and of higher dimensions in Section 8.

1.1. Outline

The paper is organized as follows.

§ 2 is a preliminary section with some interpolation results for lattice functions. Those results are the extensions of similar results in [18] to higher order interpolations, and they will be used extensively and play a key role in the forthcoming analysis.

In § 3, we set up the atomistic model and its continuum approximations. In particular, we derive the general formulation of the higher order continuum approximation. We carry out the analysis of the modeling error in § 4, and the analysis of the stability in § 5, respectively, for the higher order continuum model.

Our main result Theorem 6.2 in § 6 states that, the higher order continuum model (depending on the strain $\nabla u$ and the second order derivative of the strain $\nabla^{3}u$) which we derived from the original atomistic model has fourth order accuracy. Numerical experiments in § 7 complement and justify the theoretical analysis. We make concluding remarks and point out some promising directions for future work in § 8.

1.2. Summary of notations

We denote the directional derivative in the direction $\rho$ by $\nabla_{\rho}f(x)=\rho\cdot\nabla f(x)$. Let the symbol $\langle\cdot,\cdot\rangle$ denote the duality pairing, the first and second variations of $\mathcal{F}$ at $u$ are defined as

[TABLE]

[TABLE]

We use the convention that ’$\lesssim$’ stands for ’$\leq C$’, where $C$ is a generic constant that does not depend on the strain and its higher order derivatives.

2. Preliminary Results for the Interpolation of Lattice Functions

We define the reference lattice as $\Lambda:=\mathbb{Z}$. The space of lattice functions is given by

[TABLE]

The atomistic model is defined over the lattice functions, while the continuum model is defined over continuous functions. We will introduce certain interpolations to bridge lattice functions and continuous functions on the real line, which is the adaptation of the results in [18] to higher order interpolations.

For a lattice function $v\in\mathscr{U}$, we define the finite difference,

[TABLE]

where the finite set $\mathscr{R}\subset\{1,2,...\}$ represents the interaction range of the atomistic model. It is easy to observe that $D_{-\rho}v(\xi)=-D_{\rho}v(\xi-\rho)$.

For simplicity, we only consider periodic boundary condition in this paper, and we limit our analysis in the periodic domain $\Omega=[-N,N]$ for a fixed $N\in\mathbb{N}$. We also denote $\Omega\cap\Lambda$ by $\Omega^{\Lambda}$ .

In the following sections, we introduce three different types of interpolation for a lattice function $v\in\mathscr{U}$.

2.1. Interpolation based on finite differences

Since we are primarily interested in the a priori analysis of the higher order continuum model, we need smooth interpolates of the lattice functions which include the solutions to the atomistic model. A natural measure of the local smoothness of a lattice function $v\in\mathscr{U}$ would be the higher order finite differences, which are, however, cumbersome and of little use to our analysis. Therefore we define a smooth interpolation operator $\Pi:\mathscr{U}\rightarrow C^{4}$ such that $\Pi v\in W^{5,2}$ [6, Theorem 5.2] which is an Hermitian interpolation of degree 9, based on fourth order approximations of derivatives up to fourth order. One example of such interpolation could be defined as follows: for a lattice function $v\in\mathscr{U}$, we let

[TABLE]

where $\xi\in\Lambda$. Such interpolations will be used in the analysis of the modeling error since it satisfies both the requirement of smoothness and certain equivalence with the other two types of interpolations (c.f. Proposition 2.4 and Proposition 2.4).

2.2. Interpolation based on nodal basis functions

Let $\zeta\in W^{3,\infty}$ and $\zeta(\cdot-\xi)$ is a basis function associated with the lattice site $\xi$. We assume that $\zeta$ has a compact support and $\int_{\mathbb{R}}\zeta dx=1$. We also assume that the discrete convolution with the basis function $\zeta$ preserves cubic functions, namely

[TABLE]

One possible construction of such $\zeta$ can be the cubic spline basis function [11, Section 3.2].

We then define the standard interpolation as follows,

[TABLE]

2.3. Interpolation based on convolution

The third interpolation can be constructed by the convolution of the nodal basis interpolation $\hat{v}$ with $\zeta$ (c.f. [20]):

[TABLE]

We note that $\tilde{v}$ is in fact a quasi-interpolant of the lattice function $v$ since $\tilde{v}|_{\Lambda}\neq v|_{\Lambda}$ in general . The purpose of introducing $\tilde{v}$ is to construct the atomistic stress tensor, which will be defined in Section 4. The quasi-interpolation $\tilde{v}$ leads to the so-called localization formula [20]

[TABLE]

With the help of (2.3) we are able to rewrite the finite differences of test functions in an integral form.

2.4. Properties of the interpolation functions

We show the regularity and stability (in $W^{1,2}$ seminorm) of the convolution based interpolant $\tilde{v}$ and the finite difference based interplant $\nabla\Pi v$ in the following two propositions. The $W^{1,2}$ stability with respect to the nodal basis interpolant $\hat{v}$ shows that those three interpolations are essentially equivalent. The first proposition is similar to Proposition 3.1 of [14] and we follow the same lines of proof, while the second proposition can be found in [18].

**Proposition 2.1. ** * Let $v\in\mathscr{U}$ and $\Pi:\mathscr{U}\rightarrow C^{4}$ be a smooth interpolation operator satisfying (2.3). We then have *

[TABLE]

Proof.

The first inequality follows from Hölder’s inequality, the observation that $\int_{\xi}^{\xi+1}\nabla\Pi v(x)dx=v(\xi+1)-v(\xi)$, and Theorem 2 in [18]. The second inequality holds by Lemma 5.4 of [13]. ∎

*Proposition 2.2. ** * Let $v\in\mathscr{U}$. The nodal basis interpolation $\hat{v}$ of $v$ belongs to $W_{loc}^{3,\infty}$ and the convolution based interpolation $\tilde{v}$ belongs to $W_{loc}^{5,\infty}$. Moreover, the interpolants $\hat{v}$ and $\tilde{v}$ have the following norm equivalence

[TABLE]

Proof.

The regularity of $\hat{v}$ and $\tilde{v}$ follows from [18, Lemma 1] and the stability of $\nabla\tilde{v}$ with respect to $\nabla\hat{v}$ follows from [18, Theorem 2] since $W^{3,\infty}\subset W^{1,\infty}$. ∎

3. Atomistic Model and its Continuum Approximations

3.1. The atomistic model

In this section, we introduce the atomistic model as the ground truth description of the atomistic system. We impose a periodic boundary condition on the atomistic system, to avoid unnecessary technical difficulty which may prevent us from observing the correct convergence rate. For example, if a Dirichlet boundary condition is enforced, a suitable boundary layer of ”ghost atoms” should be added in order to guarantee the higher order convergence rate. For boundary value problems for Cauchy-Born model, please refer to [7].

Fix $N\in\mathbb{N}$, and define the space of $2N$-periodic mean zero displacements as

[TABLE]

The set of admissible deformations is given by

[TABLE]

where $F>0$ is a macroscopic deformation gradient.

The atomistic energy (per period) at a deformation $y\in\mathscr{Y}^{{\rm{per}}}$ is defined by

[TABLE]

where $\phi\in C^{3}((0,+\infty];{\mathbb{R}})$ (for example, a Lennard-Jones or Morse type potential). Using the relationship described in (3.2) and with a slight abuse of notation, the energy can be rewritten in the form of the displacement $u\in\mathscr{U}^{{\rm{per}}}$

[TABLE]

where $\phi_{\rho}(r)=\phi(r+F\rho)$ is the Lennard-Jones or Morse potential under macroscopic deformation gradient $F$. We will use (3.4) as the energy functional of our atomistic model throughout the paper.

We should also equip the space $\mathscr{U}^{{\rm{per}}}$ with $L^{2}$-norm using the interpolation of lattice function to obtain

[TABLE]

We equip the space $\mathscr{U}^{1,2}$ with the norm $||u||_{\mathscr{U}^{1,2}}:=||\nabla\hat{u}||_{L^{2}(\Omega)},\forall u\in\mathscr{U}^{1,2}$.

On the other hand, to exclude arbitrarily large deformations which are not covered by our results, we place an $L^{\infty}$-bound on the displacement gradient and define

[TABLE]

where $\kappa>0$ is a fixed constant. One reasonable choice of $\kappa$ is $\frac{1}{4}F$ which could be ensured through conditions on the external force.

It is straightforward to see that $\mathcal{E}_{\rm{a}}$ is well defined in $\mathscr{U}^{1,2}\cap\mathscr{K}$ (c.f. [20, Theorem 1]), which will be the solution space of the atomistic variational problem that will be introduced in Section 3.3.

Finally, we assume the decay hypothesis of the derivatives of the interaction potential $\phi_{\rho}$ [20, 14], which is a crucial ingredient in our analysis in Section 4. For $\rho\in\mathscr{R},1\leq j\leq k$, we require

[TABLE]

where

[TABLE]

where $\phi_{\rho}^{(j)}$ denotes the $j$th derivative of $\phi_{\rho}$. This will ensure that $\mathcal{E}_{\rm{a}}$ is $k$ times Fréchet differentiable. We note that in the current work $k\leq 6$ and the interaction range $\mathscr{R}=\{1,2,...,r_{{\rm{cut}}}\}$ is finite.

**Remark 3.1. ** We note that the assumption (3.7) does not hold at [math], for example, for Lennard-Jones or Morse potential. However for all practical purposes, we are concerned with configurations not far from reference configuration, and the atoms will not get accumulated. Therefore, it is reasonable to make the assumption (3.7). ∎

3.2. The continuum approximations

We introduce the continuum approximations of the atomistic model in this section. There are a number of approaches to obtain such approximations [28, 4, 3, 9]. We adopt the inner expansion technique in [3], which can easily satisfy the energy consistency and leads to a well-posed (the precise meaning of well-posedness will be made clear in Remark 3.2) higher order continuum model which depends on the first order and third order derivatives of $u$.

To introduce the continuum approximation, we begin with the atomistic model

[TABLE]

which is written as a sum of the site energy. Though we only require $u\in\mathscr{U}^{1,2}$ in the atomistic model, by Section 2, we can replace $u$ by its proper smooth interpolation, for example, the finite difference based interpolation $\Pi u$, without changing the atomistic energy. When no confusion occurs, we identify the discrete lattice function $u$ with its smooth interpolation in the following derivation.

After taking the Taylor expansion of the site energy $D_{\rho}u(\xi)$ at the midpoints $\xi^{\prime}:=\frac{\xi+\xi+\rho}{2}$ of the bonds $(\xi,\xi+\rho)$ and truncating at order three, we have

[TABLE]

The atomistic model can then be approximated as

[TABLE]

An approximation step similar to the Riemann sum leads to the following higher order continuum (HOC) approximation

[TABLE]

where $\phi_{\rho}$ is the Lennard-Jones or Morse potential under macroscopic deformation gradient $F$ introduced in Section 3.1 and $\nabla u=\nabla u(x)$ is the gradient of $u$ with respect to $x$.

The well-known Cauchy-Born approximation can be obtained by preserving only the first order term in (3.10)

[TABLE]

We note that the higher order energy functional defined in (3.10) depends both on $\nabla u$ and $\nabla^{3}u$, whereas the Cauchy-Born energy functional (3.11) depends only on $\nabla u$.

**Remark 3.2. ** In [28] the authors derived two higher order continuum models which contain both $\nabla u$ and $\nabla^{2}u$ terms from the atomistic model. However, it was discovered in [3, Section 4] that these two higher order continuum models are ill-posed, which led to an uncontrolled behavior of the solution. For example, if we take the harmonic potential $\phi(r)=\frac{1}{2}(r-1)^{2}$ as the atomistic potential, one of the (ill-posed) higher order continuum models is

[TABLE]

and the corresponding Euler-Lagrangian equation of (3.12) is

[TABLE]

We observe from (3.12) that the energy is not positive definite, and the differential operator on the left hand side of (3.13) is not elliptic. (3.12) is ill-posed in this sense, which also means the energy (3.12) is not stable in the sense of (3.25).

It is possible to derive a well-posed higher order continuum model (3.10) through the inner expansion technique in [3]. For example, for the above harmonic potential case, we can obtain

[TABLE]

and the corresponding Euler-Lagrange equation of (3.14) is

[TABLE]

which is actually a special case of the general (nonlinear) Euler-Lagrange equation (4.1). (3.14) is well-posed, and thus is stable. ∎

3.3. The Variational problems

In this section, we define the variational problem for both atomistic model and higher order continuum models.

3.3.1. The Atomistic problem

We apply an external force $f$ to the atomistic system, in order to generate a nontrivial solution to the atomistic model. Following previous literature [20, 9, 19, 21], the external force of the atomistic model is modeled as a dead load so that the work of the external force is given by

[TABLE]

where $u\in\mathscr{U}^{1,2}$ is a displacement and $f|_{\Lambda}\in\mathscr{U}^{1,2}$. The atomistic problem is: find a (local) minimizer $u^{\rm{a}}$ such that

[TABLE]

If $u^{{\rm{a}}}$ is a solution to (3.17), then it satisfies the first-order condition

[TABLE]

We call a solution $u^{{\rm{a}}}$ of (3.17) (strongly) stable if there exists $c_{0}>0$ such that

[TABLE]

3.3.2. The Higher order continuum problem

To define the variational problem with respect to $\mathcal{E}_{\rm{hoc}}$, we first introduce the following space

[TABLE]

In order to apply the inverse function theorem (Lemma 6.2) to obtain the error estimate $\|\nabla u^{\rm{a}}-\nabla u^{\rm{hoc}}\|_{L^{2}}$, we equip the space $\mathcal{U}^{1,2}$ with the $W^{1,2}$ norm

[TABLE]

We denote $W^{-1,2}$ as the standard topological dual of $W^{1,2}$.

We define the $L^{2}$ inner product for $u,v\in\mathcal{U}^{1,2}$,

[TABLE]

Similar to $\mathscr{K}$, we shall assume that all displacement gradients satisfy a uniform bound. To that end we define

[TABLE]

where $\kappa$ is the same constant as in the definition of $\mathscr{K}$. Notice that $\mathcal{K}$ is an open set in $\mathcal{U}^{1,2}$, and $\mathcal{E}_{\rm{hoc}}$ is well defined on $\mathcal{K}$.

For the higher order continuum model, assume that the external force $f\in\mathcal{U}^{1,2}$, we seek the solution for the following variational problem:

[TABLE]

The solution $u^{{\rm{hoc}}}$ to (3.23) satisfies the first-order condition

[TABLE]

We call the solution $u^{{\rm{hoc}}}$ of (3.23) (strongly) stable, if there exists a positive number $\gamma_{0}$ such that

[TABLE]

4. Modeling Error Analysis

In this section, we give a rigorous analysis of the modeling error. We first introduce the atomistic stress tensor $S^{{\rm{a}}}(u;x)$ and the stress of the higher continuum model $S^{{\rm{hoc}}}(u;x)$ in Section 4.1. Then we derive the pointwise error estimate in stress $R(u;x)=S^{{\rm{a}}}(u;x)-S^{{\rm{hoc}}}(u;x)$ in Section 4.2. Finally, we present the fourth-order consistency estimate of the higher order continuum model (3.10) in Section 4.3.

4.1. Atomistic and continuum stresses

The first variation of the atomistic energy functional $\mathcal{E}_{\rm{a}}$ (3.4) at $u\in\mathscr{U}^{1,2}$, is given by

[TABLE]

We replace the test function $v$ by its convolution based quasi-interpolation $\tilde{v}$ and apply the localization formula (2.3), it follows that

[TABLE]

where

[TABLE]

is defined to be the atomistic stress tensor. The last two identities in (4.2) hold because of the periodic boundary condition and $\chi_{\xi,\rho}(x)=0$ when $|\xi-x|>2r_{{\rm{cut}}}$.

The first variation of the higher order continuum energy functional $\mathcal{E}_{{\rm{hoc}}}$ defined in (3.10) is given by

[TABLE]

Integration by parts and the periodic boundary condition of the test function lead to

[TABLE]

where

[TABLE]

is defined to be the stress of the higher order continuum model (3.10).

By the integration by parts to (4.5) again, we obtain

[TABLE]

where

[TABLE]

$W^{{\rm{hoc}}}(u)=0$ is the Euler-Lagrange equation of the higher order continuum model, which is a sixth order nonlinear elliptic equation.

We now define the error in stress as $R(u;x):=S^{{\rm{a}}}(u;x)-S^{{\rm{hoc}}}(u;x)$. In the remaining part of this section, we will give the pointwise estimate of $R(u;x)$ and show the fourth order consistency of the higher order continuum model (3.10).

4.2. Pointwise estimate of the error in stress

In this section we prove the pointwise estimate of $R(u;x)$, the error in stress. We first introduce a useful lemma which is a direct extension of [20, Lemma 11].

*Lemma 4.1. ** * Let $x,\rho\in\Omega$, $k=0,1,2,3$, and $\chi_{\xi,\rho}(x)$ is defined by (2.3). We have

[TABLE]

Proof.

This result relies on the assumption that it is true on a shifted grid: if $v:\mathbb{R}\rightarrow\mathbb{R}$ is a polynomial whose order is less than $k$, where $k=0,1,2,3$, then for any $z,x\in\mathbb{R}$ we have

[TABLE]

To prove the result, let $s\in[0,1]$ be fixed, then

[TABLE]

where we substituted $\eta=-(\xi-x)$ and employing (4.10) with $v(x)\equiv(-x)^{k}$, we obtain

[TABLE]

By the definition of $\chi_{\xi,\rho}(x)$ in (2.3) and by integrating w.r.t. $s$, we have

[TABLE]

It is trivial to see that $\sum_{\xi\in\Lambda}\chi_{\xi,\rho}(x)=1$ if we let $k=0$. ∎

The pointwise estimate of the error in stress is given by the following lemma.

*Lemma 4.2. ** * Let $u\in W^{5,\infty}\cap\mathcal{K}$, and $x\in\Omega$, then

[TABLE]

*where $C$ depends on $M^{(j,4)},j=2,...,5$, defined in Section 3.1, and $v_{x}:=B_{2r_{cut}+1}(x)$ is the neighbourhood of some $x\in\mathbb{R}$ and $r_{cut}=\max_{r\in\mathscr{R}}|r|$. *

Proof.

In order to keep the notation concise, we first define

[TABLE]

By a direct Taylor expansion of (4.1) and using the fact $\sum_{\xi\in\Lambda}\chi_{\xi,\rho}(x)=1$, we can rewrite the stress of the higher order continuum model as

[TABLE]

where $\mu_{j}\in{\rm{conv}}\{\nabla_{\rho}u,\nabla_{\rho}u+\frac{\rho^{3}}{24}\nabla^{3}u\},j=1,2,3,4$.

We then turn our attention to the atomistic stress tensor in (4.3) where

[TABLE]

Since $\zeta$ has a compact support, we have $\chi_{\xi,\rho}(x)=0$ for all $\xi\in\Lambda$ with $|\xi-x|>2|\rho|$. We thus can apply Taylor expansion to the term $\phi^{\prime}_{\rho}(D_{\rho}u(\xi))$ at $x$. We begin by expanding $D_{\rho}u(\xi)$ for $\rho\in\mathscr{R}$ in $v_{x}$, in the neighbourhood of $x$, so that

[TABLE]

The fact that $D_{\rho}u(\xi)-\nabla_{\rho}u=O(\tau_{2})$ allows us to expand $\phi^{\prime}_{\rho}(D_{\rho}u(\xi))$ as

[TABLE]

where $\mu_{5}\in{\rm{conv}}\{\nabla_{\rho}u,D_{\rho}u(\xi)\}$.

Combining (4.14), (4.13), (4.15) and (4.2) and after some algebraic manipulation, we have

[TABLE]

Lemma 4.2 leads to the following identities, which result in the elimination of $(a)-(f)$ terms in (4.17).

[TABLE]

We now only need to estimate the remaining terms in (4.17). By the definition of $\chi_{\xi,\rho}$, we have the estimate $\sum_{\xi\in\Lambda}\chi_{\xi,\rho}h_{\xi}\leq\max_{\xi\in\Lambda,\chi_{\xi,\rho}\neq 0}h_{\xi}$, where $h_{\xi}$ is an arbitrary function with respect to $\xi$. Combined with the boundedness of the derivatives of the interaction potential $\phi_{\rho}$ assumed in (3.7), it is easy to show that

[TABLE]

and this finishes the proof. ∎

**Remark 4.3. ** Notice that the lower order terms ($\leq 4$) in (4.17) vanish for different reasons. Terms $(a)$, $(d)$, $(e)$, $(f)$ come from the atomistic model only. Terms $(b)$, $(c)$ come from both atomistic model and higher order continuum model, and the cancellation of those terms depends on the choice of expansion points in (3.8). For example, if we choose the lattice points instead of the mid points as expansion points in (3.8), we can get the following higher order continuum energy,

[TABLE]

which is only first order accurate. ∎

**Remark 4.4. **

In fact, the higher order continuum model which is originally derived in [28],

[TABLE]

also has fourth order estimate in stress. As a matter of fact, (3.10) and (4.20) differs only a null-Lagrangian for second order term. However, this model is ill-posed and thus is not stable. See Remark 3.2 for more details.

We give a simplified error analysis for the higher order continuum model with 6th order accuracy in Appendix A, which truncate the terms of order 5 onwards in (3.8).

From those observations, we conjecture that: we need to include higher order gradient up to $2k+1$th order to obtain a ”well-posed” higher order continuum model of order $2k+2$.

∎

Construction of higher order continuum model for more physical relevant cases of multi body iterations and/or higher dimensions will be discussed in Section 8. In those cases, a similar but more involved formulation of stress differences as (4.17) will serve as the key of developing and analyzing higher order continuum models.

4.3. Fourth-order consistency of the higher order continuum model

Lemma 4.2 gives us the upper bound of $|R(u;x)|$, we now convert this pointwise estimate into a global estimate. The main idea is to use the inverse estimates to obtain $L^{2}$ type bounds from the $L^{\infty}$ bounds. It is easy to show that

[TABLE]

where the interpolation operator $\Pi:\mathscr{U}^{1,2}\rightarrow C^{4}$ is defined in (2.3).

*Theorem 4.5. ** * Let $u\in\mathscr{U}^{1,2}\cap\mathscr{K}$ and $\tilde{v}=\zeta*\hat{v}$. The interpolation operator $\Pi:\mathscr{U}^{1,2}\rightarrow C^{4}$ is defined in Section 2.1 by (2.3). Then the model error is bounded by

[TABLE]

*where $C$ depends on $M^{(j,4)},j=2,...,5$ which are defined in Section 3.1. *

Proof.

From the definition of $\Pi$ by (2.3), we observe that $\Pi u|_{\Omega^{\Lambda}}=u|_{\Omega^{\Lambda}}$. It follows from the localization formula (2.3) that,

[TABLE]

An application of Cauchy-Schwarz inequality yields

[TABLE]

Using the inverse estimates (4.21), we obtain

[TABLE]

where $v_{x}$ is a compact support of $x$ defined in Lemma 4.2. Integrating (4.24) over $\Omega$, we have

[TABLE]

which yields the stated result. ∎

5. Stability

In this section, we present the stability estimate of the higher order continuum model (3.10), which extends the stability results for Cauchy-Born model in [12, Theorem 3.1]. The atomistic model (3.3) and its solution space $\mathscr{U}^{1,2}$ can be denoted as $\mathcal{E}_{{\rm{a}}}^{N}(u)$ and $\mathscr{U}^{1,2}_{N}$ since they actually depend on the computational domain $\Omega^{\Lambda}=[-N,N]\cap\Lambda$. For a fixed $N\in\mathbb{N}$, given potential $\phi_{\rho}(r)=\phi(r+F\rho)$ defined in Section 3.1, we call the homogeneous deformation $y=Fx$ is stable in the finite atomistic model if

[TABLE]

We require a stronger definition of the stability in the infinite atomistic model:

[TABLE]

Also, the homogeneous deformation is stable for the Cauchy-Born model (3.11) if

[TABLE]

and the homogeneous deformation is stable for the higher order continuum model (3.10) if

[TABLE]

The following lemma states the stability of the higher order continuum model (3.10) at the homogeneous deformation, namely: the stability of the atomistic model implies that of the higher order continuum model, and the stability of the higher order continuum model (3.10) is ”in between” the atomistic model and the Cauchy-Born model.

**Lemma 5.1. ** * If the deformation gradient $F$ introduced in Section 3.1 is positive, then $\Lambda_{{\rm{a}}}\leq\Lambda_{{\rm{hoc}}}\leq\Lambda_{{\rm{cb}}}$. *

Proof.

For the first inequality, we extend the proof given in [12, Section 3.1] to higher order continuum model. The energy $\mathcal{E}^{N}_{{\rm{a}}}$ and $\mathcal{E}_{{\rm{hoc}}}$ can be expanded up to second order for an arbitrary small $t>0$ and $\hat{v}\in C^{3}(\Omega)$,

[TABLE]

where

[TABLE]

By (3.10) and (3.4), we have $\mathcal{E}^{N}_{{\rm{a}}}(0)=\mathcal{E}_{{\rm{hoc}}}(0)$ and $\lim_{N\rightarrow\infty}\mathcal{E}^{N}_{{\rm{a}}}(0+t\hat{v})=\mathcal{E}_{{\rm{hoc}}}(0+t\hat{v})$. Hence we have that

[TABLE]

letting $t\rightarrow 0$, we have $\lim_{N\rightarrow\infty}\langle\delta^{2}\mathcal{E}^{N}_{{\rm{a}}}(0)\hat{v},\hat{v}\rangle=\langle\delta^{2}\mathcal{E}_{{\rm{hoc}}}(0)\hat{v},\hat{v}\rangle$. According to the definition of $\Lambda_{{\rm{hoc}}}$, there exists $\hat{v}_{\delta}\in\mathcal{U}^{1,2}$ such that $||\nabla\hat{v}_{\delta}||_{L^{2}(\Omega)}=1$ and $\langle\delta^{2}\mathcal{E}_{{\rm{hoc}}}(0)\hat{v}_{\delta},\hat{v}_{\delta}\rangle\leq\Lambda_{{\rm{hoc}}}+\delta$. We thus obtain

[TABLE]

The first inequality then follows since $\delta$ can be chosen arbitrarily small.

For the second inequality, the stability constants $\Lambda^{N}_{{\rm{a}}}$ and $\Lambda_{{\rm{cb}}}$ have explicit characterizations in [12, Section 3.2] using Fourier transform

[TABLE]

where $h(0)=\frac{\partial^{2}\mathcal{E}^{N}_{{\rm{a}}}(0)}{\partial u(0)^{2}}$ is positive and $\SS:=\{a\in\mathbb{R}:|a|=1\}$. For the higher order continuum model (3.10), we have

[TABLE]

We observe that $\varphi_{{\rm{hoc}}}$ is actually the truncated Taylor expansion of $\varphi_{{\rm{a}}}$ up to order 3, while $\varphi_{{\rm{cb}}}$ only preserves the first order term, which indicates the stated result $\Lambda_{{\rm{a}}}\leq\Lambda_{{\rm{hoc}}}\leq\Lambda_{{\rm{cb}}}$ if we assume $N$ is sufficiently large.

∎

For the stability of the higher order continuum model at small deformations, we have the following Theorem.

**Theorem 5.2. ** If the condition in Lemma 5 and (5.2) are satisfied, we then have

[TABLE]

Proof.

We note that $u\in\mathcal{U}^{1,2}\cap\mathcal{K}$ can be taken as a perturbation of the reference configuration. Combining the higher order continuum model (3.10) and the Lipschitz continuity of the potential $\phi_{\rho}$, we have

[TABLE]

We finish the proof by choosing $\kappa\leq\Lambda_{{\rm{hoc}}}/(2M^{(3,0)})$ and applying Lemma 5. ∎

6. A priori Error Estimates

In this section, we present the main result of the a priori error estimate, Theorem 6.2, which essentially shows that the minimizer of the higher order continuum model (3.10) has fourth-order accuracy.

6.1. Consistency error for the external work

We first present the following lemma which shows that the approximation error of the external energy is of fourth order.

*Lemma 6.1. ** * (Consistency error for the external work) Suppose $f\in\mathcal{U}^{1,2}$. Let $\langle f,\hat{v}\rangle_{\Omega^{\Lambda}}$ and $\langle f,\hat{v}\rangle_{\Omega}$ be defined in (3.16) and (3.21) respectively. We have the following estimate:

[TABLE]

Proof.

By the definition of $\tilde{v}(x)$ by (2.6), we have

[TABLE]

By the mean-zero condition for $f$ and the property that $\int\zeta dx=1$ , it is easy to show that $\int_{\Omega}g(x)dx=0$. Hence for an arbitrary constant $c_{\Omega}\in\mathbb{R}$, an application of Cauchy-Schwarz inequality yields that

[TABLE]

Choosing $c_{\Omega}=\frac{1}{|\Omega|}\int_{\Omega}\hat{v}dx$ and applying $Poincar\acute{e}$ inequality, we obtain the estimate that

[TABLE]

According to the standard Bramble-Hilbert lemma in [6, Theorem 6.4], we can estimate the $L^{2}$-norm of $g$ by

[TABLE]

which can also be obtained by a direct extension of [18, Lemma 13]. Combination of (6.4) and (6.5) leads to the required result. ∎

6.2. A priori error estimate

We first state the well-known inverse function theorem [17, Lemma 2.2].

Lemma 6.2. ** (Inverse function theorem) Let $\mathcal{A},\mathcal{B}$ be Banach spaces, $\mathcal{O}$ an open set of $\mathcal{A}$, and let $\mathcal{F}:\mathcal{O}\rightarrow\mathcal{B}$ be Fréchet differentiable with Lipschitz-continuous derivative $\delta\mathcal{F}$:

[TABLE]

where $M$ is a Lipschitz constant. Let $X\in\mathcal{O}$ and suppose also that there exists $\eta,\sigma>0$ such that

[TABLE]

[TABLE]

*Then there exists a locally unique $Y\in\mathcal{A}$ such that $\mathcal{F}(Y)=0$ and $||Y-X||_{\mathcal{A}}\leq 2\eta\sigma$.

The following theorem shows that for a stable and sufficiently small deformation, the solution of the higher order continuum model (3.10) is a good approximation to the solution of the atomistic model.

*Theorem 6.3. ** * (A priori error estimate) Suppose (5.2) is satisfied and $u^{\rm{a}}$ is a strongly stable atomistic solution of (3.17). $\Pi$ is defined in (2.3) and $f\in\mathcal{U}^{1,2}$. If we assume that $\eta_{1},\eta_{2}$ are sufficiently small such that $||\nabla^{j}\Pi u^{{\rm{a}}}||_{L^{2}(\Omega)}\leq\eta_{1}$, $j=2,3,4,5$, and $||\nabla^{4}f||_{L^{2}(\Omega)}\leq\eta_{2}$, there exists a stable solution $u^{{\rm{hoc}}}$ of problem (3.23) in $W^{1,2}$ such that

[TABLE]

*where $C$ depends on $\Lambda_{{\rm{hoc}}}$, $M^{(j,4)}$, $j=2,...,5$. *

Proof.

Following the framework of the a priori error estimates in [13, 14, 21], we divide the proof into three steps. Recalling the definition of the space $\mathcal{U}^{1,2}$ in Section 3.3, we apply Lemma 6.2 with $\mathcal{A}:=\mathcal{U}^{1,2},\mathcal{O}:=\mathcal{K},\mathcal{B}:=W^{-1,2},X:=\Pi u^{\rm{a}}$. We define the operator $\mathcal{F}:\mathcal{O}\rightarrow\mathcal{B}$ by

[TABLE]

$\delta\mathcal{F}$ is Lipschitz continuous due to the Lipschitz continuity of the potential $\phi_{\rho}$. We also note that $X\in\mathcal{O}$ since $u^{\rm{a}}\in\mathscr{U}^{1,2}\cap\mathscr{K}$ and $\eta_{1}$ is chosen to be sufficiently small.

Step 1: Stability.     In Section 5, we have already shown the stability of the higher order continuum model, that is,

[TABLE]

where $\Lambda_{{\rm{hoc}}}$ is positive. Hence, we have

[TABLE]

Step 2: Consistency.     We shall also require $\mathcal{F}(X)$ to be consistent. In fact, if we put the interpolation of the atomistic solution $u^{{\rm{a}}}$ in $\mathcal{E}_{{\rm{hoc}}}$, we have

[TABLE]

where we decompose the consistency into two parts: $T_{1}$ is the modeling error and $T_{2}$ is the consistency error of the external force. Applying Theorem 4.3 and Lemma 6.1, we have

[TABLE]

where $C$ depends on $M^{(j,4)}$, $j=2,...,5$.

Step 3: Inverse function theorem.     Combing the stability result in Step 1 and the consistency result in Step 2 and applying Lemma 6.2, we obtain the existence of the solution of the higher order continuum model $u^{{\rm{hoc}}}$ in $W^{1,2}$ and the error estimate

[TABLE]

which can be guaranteed if we choose $\eta_{1},\eta_{2}$ to be sufficiently small. ∎

**Remark 6.4 (Scaling). **

Up to now, we have taken the unit interatomic spacing in the reference lattice. In order to illustrate the order of accuracy, we scale the interatomic spacing by $\varepsilon$, that is, $X:=\varepsilon x$, $U:=\varepsilon u$ and $F:=\varepsilon^{-1}f$. Reversing the scaling, we have $u^{\rm{a}}(x):=\varepsilon^{-1}U^{\rm{a}}(\varepsilon x)$ of the atomistic problem (3.18) and the external force $f(x):=\varepsilon F(\varepsilon x)$. It can be easily shown that $||\nabla^{j}\Pi u^{{\rm{a}}}||_{L^{2}(\Omega)}=\varepsilon^{j-1-\frac{1}{2}}||\nabla^{j}\Pi_{\varepsilon}U^{{\rm{a}}}||_{L^{2}(\Omega)}$, $j=1,2,3,4,5$, where $\Pi_{\varepsilon}$ is the same interpolation as $\Pi$ under $\varepsilon$ scale. We also scale the estimate (6.5) by $||\nabla^{j}f||_{L^{2}(\Omega)}=\varepsilon^{j+1-\frac{1}{2}}||\nabla^{j}F||_{L^{2}(\Omega)}$, $j=0,1,2,3,4$. Hence, Theorem 6.2 essentially shows that our higher order continuum model is of fourth order accuracy, namely,

[TABLE]

∎

**Remark 6.5 (higher regularity). ** In Theorem 6.2, we only prove the existence of the minimizer in $W^{1,2}$, although the approximation space $\mathcal{U}^{1,2}\cap\mathcal{K}$ is more restrictive. It is possible to prove higher regularity ($W^{5,2}$) given ellipticity of the Euler-Lagrange equation, similar as the $W^{3,2}$ regularity of the Cauchy-Born solution as in [20, Proof of Theorem 3]. In this paper, we are mainly concerned with the $W^{1,2}$ accuracy for the solution of the higher order continuum model.

∎

Similar to [13, Proposition 3.2], we also have the estimate for the error in energy $|\mathcal{E}_{{\rm{a}}}(\Pi u^{{\rm{a}}})-\mathcal{E}_{{\rm{hoc}}}(u^{{\rm{hoc}}})|$.

*Theorem 6.6. ** * (Energy estimate) Under the conditions of Theorem 6.2, we have

[TABLE]

*where $C$ depends on $M^{(j,3)}$, $j=2,...,5$. *

7. Numerical Experiments

We present two numerical experiments to illustrate the analytical results of this paper. We include up to second nearest neighbor interactions in our energy functional such that

[TABLE]

and

[TABLE]

We set the computational domain to be $\Omega=[-1,1]$. In the reference configuration, there are $2N+1$ equally distributed atoms in $\Omega$, hence the scaling parameter is $\varepsilon:=1/N$. We choose

[TABLE]

as the external force so that a nonlinear but small enough displacement (or equivalently deformation) is generated. It is easy to see that $f\in\mathcal{U}^{1,2}$. We will carry out numerical experiments for both the harmonic potential \displaystyle\phi(r)=\frac{1}{2}\big{(}\frac{r}{\varepsilon}-1\big{)}^{2} and the Leonard-Jones potential \displaystyle\phi(r)=\big{(}\frac{r}{\varepsilon}\big{)}^{-12}-2\big{(}\frac{r}{\varepsilon}\big{)}^{-6}.

We use $C^{3}$ finite element to solve the variational form (4.4) of higher order continuum model. We denote the positions of the atoms to be $\{x_{i}\}^{2N+1}_{i=1}$, and let the nodes of the finite elements coincide with the atoms so that $\Omega$ is partitioned by $\mathcal{T}:=\{T_{i}\}_{i=1}^{2N+1}$ where $T_{i}:=[x_{i-1},x_{i}]$. The finite element solution space is then defined by

[TABLE]

where $\mathbb{P}_{5}$ is the quintic polynomial function space. We search the approximate local minimizer of $\mathcal{E}_{\rm{hoc}}(u)$ defined by (7.2) in $\mathcal{U}^{1,2}_{\varepsilon}$ using BFGS algorithm.

We use Guass-Legendre quadrature of order 5 to approximate the integral $\int_{\Omega}fudx$ which is the external work in the higher order continuum model to make the quadrature error negligible compared to the consistency error given in (6.1).

We use the following protocol to quantify the modeling error:

1. (1)

Let $\varepsilon_{1}=2^{-3}$, $\varepsilon_{2}=2^{-4}$, … , $\varepsilon_{8}=2^{-10}$ be the interatomic spacing which also define the reference lattice and the finite element mesh. 2. (2)

Compute the atomistic solution $u^{{\rm{a}},\varepsilon_{i}}$ and the higher order continuum solution $u^{{\rm{hoc}},\varepsilon_{i}}$ on different lattices (or corresponding meshes). 3. (3)

Compute the error $\|\nabla Iu^{{\rm{a}},\varepsilon_{i}}-\nabla u^{{\rm{hoc}},\varepsilon_{i}}\|_{L^{2}}$, where $I$ represents a smooth interpolation operator such that the interpolation error is negligible compared to the consistency error.

**Remark 7.1. ** The numerical error for the higher order continuum model is fifth order given the approximation space $\mathcal{U}^{1,2}_{\varepsilon}$ and the higher regularity of $u^{{\rm{hoc}}}$ (see Remark 6.2), we have,

[TABLE]

∎

7.1. The harmonic potential

We first give a numerical justification of our main result for the harmonic potential \displaystyle\phi(r)=\frac{1}{2}\big{(}\frac{r}{\varepsilon}-1\big{)}^{2}. Fixing $\varepsilon=2^{-3}$, we compute the solutions of the atomistic model (3.4), Cauchy-Born model (3.11), and higher order continuum model (3.10).

In Figure 1, we observe that the solution of the higher order continuum model ($\rhd$) is closer to the solution of the atomistic model ($\circ$) than that of the Cauchy-Born model ($\bullet$).

We then compute and plot the error $\|\nabla Iu^{{\rm{a}},\varepsilon}-\nabla u^{{\rm{hoc}},\varepsilon}\|_{L^{2}}$ in Figure 2, notice that the modeling error $\|\nabla\Pi u^{{\rm{a}}}-\nabla u^{{\rm{hoc}}}\|_{L^{2}}$ is the dominant part in $\|\nabla Iu^{{\rm{a}},\varepsilon}-\nabla u^{{\rm{hoc}},\varepsilon}\|_{L^{2}}$. Figure 2 clearly shows the fourth order accuracy of the higher order continuum model, compared with second order accuracy of the Cauchy-Born model. We need to mention here the importance of the proper choice of the interpolation operator $I$ for the atomistic solution $u^{{\rm{a}},\varepsilon}$.

We use different interpolation operators $I$ in Figure 2, interpolation operator $\Pi$ and quartic splines interpolation can preserve the 4th order accuracy, while cubic spline interpolation gives suboptimal results (3rd order accuracy).

7.2. Lennard-Jones Potential

Our second numerical example is for Lennard-Jones potential \displaystyle\phi(r)=\big{(}\frac{r}{\varepsilon}\big{)}^{-12}-2\big{(}\frac{r}{\varepsilon}\big{)}^{-6} with the interpolation $Iu^{{\rm{a}},\varepsilon}$ being the quartic spline. Figure 3 shows that the order of the modeling error is not affected by the nonlinearity of the potential.

We also plot the error in energy in Figure 4. We see that the error in energy is of fourth order for the higher order continuum model compared with the second order for the Cauchy-Born model, which is consistent with Theorem 6.2.

8. Conclusion and future work

In this paper, we derive a higher order nonlinear elasticity model from the atomistic model in one dimension, and present a rigorous a priori error analysis for this higher order continuum model. Using the techniques developed in [17, 20, 13], we prove that the modeling error of our higher order continuum model is fourth order, compared with the second order accuracy of the well-known Cauchy Born model. Numerical experiments are carried out to verify our theoretical results.

This work opens up several interesting research directions:

The first direction is the extension of the current work to general multibody interactions. We note that the inner expansion technique in [3] does not apply in this case, since it requires that the site potential $\phi_{\bar{\xi}}({y(\xi)}_{\xi\in\Lambda})$ for the site $\bar{\xi}$ can be written as

[TABLE]

where $\psi$ is a function and $a_{\xi}$ are constants. In fact, the essential problem here, is how to determine the ”optimal” expansion points in the Taylor expansion of the energy functional. While it is natural to use the midpoints of the bonds as the expansion points for pair interactions, multibody interactions involve a number of bonds and thus have cross terms in the Taylor expansion. This may lead to various possible formulations of the higher order continuum model corresponding to different choices of expansion points. We need to choose the expansion points ”optimally” such that the cancellation in (4.17) can be achieved at the highest possible order.

The second direction is the construction and analysis of the higher order continuum model in higher dimensions. Such extension seems to be straightforward following the framework proposed in the current work. However, we note that Lemma 4.2, which serves as the key to the cancellation of the lower order terms in (4.17), in general does not hold for $k>1$ in higher dimensions. To be more precise, in higher dimensions, the tensor product such as $\rho\otimes\rho$ will appear in the stress difference $R(u;x)$. Therefore an alternative identity for the cross terms should be sought after in higher dimensions to guarantee the accuracy.

The third direction is the development of atomistic/continuum coupling with the higher order continuum model. The key to design ”optimal” coupling method is to balance the modeling error with coarsening and truncation errors, within the analytical framework in [15, 29, 13, 24, 23]. The reduction of modeling error by the higher order models can facilitate the construction of coupling method with quasi-optimal convergence rate. In particular, for complex lattices, the modeling error of Cauchy Born is only first order due to lack of symmetry and becomes the bottleneck for the coupling method. We expect higher order models can be used to alleviate this problem. Once the coupling model is developed, the study of adaptivity should be under way where [1, 2, 22, 30, 32, 31] should provide good references.

Acknowledgements

Funding

YW and LZ were supported by NSFC grant 11871339, 11861131004, 11571314. HW was supported by National Science Foundation China grant No.11501389, 11471214.

Appendix A. Derivation and Justification of the higher order continuum model with sixth order accuracy

In this section, we present the derivation and the theoretical justification of the higher order continuum model which contains $\nabla u$, $\nabla^{3}u$ and $\nabla^{5}u$ with 6th order consistency. The stability and convergence of this model can be verified similar as 4th order model (3.10).

Following the derivation introduced in Section (3.2), we extend (3.8) by truncating the terms whose orders are higher than 5, we have

[TABLE]

The corresponding higher order continuum model reads

[TABLE]

The energy functional (.2) depends on $\nabla u$, $\nabla^{3}u$ and $\nabla^{5}u$, and we will show it has 6th order consistency. We introduce the space by imposing periodic boundary condition and mean zero condition on it:

[TABLE]

The first variation of the higher order continuum energy functional $\mathcal{E}^{*}_{{\rm{hoc}}}$ is given by

[TABLE]

Using the integration by parts and the periodic boundary condition of the test function, we obtain

[TABLE]

where

[TABLE]

Following the analysis in Section 4, we give a pointwise sixth-order consistency estimate of the stress error $R(u;x)=S^{{\rm{a}}}(u;x)-S^{{\rm{hoc}}}_{*}(u;x)$. The key step is the cancellation of the terms whose order is less than six in (Appendix A. Derivation and Justification of the higher order continuum model with sixth order accuracy), which can be achieved by the extension of Lemma 4.2 to $k=5$. If we choose the basis function $\zeta$ as the quintic spline basis function, then it preserves the fifth order polynomials. Following the proof of Lemma 4.2, we have

[TABLE]

The atomistic stress is the same as (4.3):

[TABLE]

We then expand $D_{\rho}u(\xi)$ at $x$ for $\rho\in\mathscr{R}$ in $v_{x}$, which is the neighbourhood of $x$.

[TABLE]

The fact that $D_{\rho}u(\xi)-\nabla_{\rho}u=O(\tau_{2})$ allows us to expand $\phi^{\prime}_{\rho}(D_{\rho}u(\xi))$:

[TABLE]

where $\mu_{7}\in{\rm{conv}}\{\nabla_{\rho}u,D_{\rho}u(\xi)\}$.

Combining (.8), (Appendix A. Derivation and Justification of the higher order continuum model with sixth order accuracy), (Appendix A. Derivation and Justification of the higher order continuum model with sixth order accuracy), (Appendix A. Derivation and Justification of the higher order continuum model with sixth order accuracy), after some algebraic manipulations, we can obtain the stress error $R(u;x)$ in (Appendix A. Derivation and Justification of the higher order continuum model with sixth order accuracy). With a slight abuse of notation, for the convenience and clearness of the expression, we divide the $R(u;x)$ into several parts which indicates the different orders if we rescale the displacement $u$ by $\varepsilon$. Using the identities (.7), all the terms in $R(u;x)$ whose order are lower than six cancel out.

[TABLE]

Hence we obtain the pointwise sixth order consistency estimate of the stress error $R(u;x)$. Following the analysis in Section (4.3), Section (5) and Section (6), we can similarly prove the higher order continuum model (.2) is of sixth order accuracy.