# Lattice dynamics from a continuum viewpoint

M. Charlotte, Lev Truskinovsky

► To cite this version:

M. Charlotte, Lev Truskinovsky. Lattice dynamics from a continuum viewpoint. Journal of the Mechanics and Physics of Solids, 2012, 60 (8), pp.1508-1544. 10.1016/j.jmps.2012.03.004. hal-01853071

HAL Id: hal-01853071

https://hal.science/hal-01853071

Submitted on 2 Aug 2018

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
**To cite this version:** Charlotte, Miguel and Truskinovsky, Lev *Lattice dynamics from a continuum viewpoint*. (2012) Journal of the Mechanics and Physics of Solids, vol. 60 (n° 8). pp. 1508-1544. ISSN 0022-5096
# Lattice dynamics from a continuum viewpoint

Miguel Charlotte<sup>a,*</sup>, Lev Truskinovsky<sup>b</sup>

<sup>a</sup> Université de Toulouse, ISAE, INSA, UPS, EMAC, Institut Clément ADER (ICA), BP 54032, 31055 Toulouse, France

<sup>b</sup> Laboratoire de Mécanique des Solides, Ecole Polytechnique, 91128 Palaiseau, France

## ABSTRACT

We develop a novel continuum description of lattice dynamics for a particle chain with nearest neighbor interactions. Our continuum model interpolates discrete solutions and allows one to deal adequately with singular and impact loadings. The resulting theory is local in space and nonlocal in time. It is characterized by a nontrivial hereditary structure which blends inertial and elastic forces. The proposed methodology can be used to build continuum approximations for lattice dynamics which incorporate successively finer scales. It can also serve as a seamless interface between discrete and continuum elasticity theories.

Keywords:

Lattice dynamics

Quasi-continuum modeling

Nonlocal elasticity

Discrete to continuum

Multiscale modeling

Non-Newtonian inertia

Meta-materials

## 0. Introduction

The two most widely used theories of elastic deformation in crystalline materials are the atomic scale molecular dynamics (MD) (Born and Huang, 1954) and the classical continuum (CC) thermo-elasticity (Landau et al., 1986). Continuum elasticity is a scale free phenomenological theory representing a formal coarse grained approximation of molecular dynamics. Recently, considerable efforts have been made to obtain the CC theory as a rigorous mathematical limit of MD, however, the progress is still modest outside the well understood case of linear elasto-statics (Hays et al., 1994; Curtin and Miller, 2003; Park et al., 2005; E and Li, 2005; Charlotte and Truskinovsky, 2008; Marian et al., 2010; Tadmor and Miller, 2011). Despite the remaining difficulties, the classical theory shows validity over a remarkably broad range of scales, from kilometers to micrometers, and allows one to use powerful analytical tools of differential calculus to obtain relatively simple description of the phenomena involving large number of particles participating in complex nonlinear interactions.

The continuum elasto-dynamics fails at the atomic scale where it misses micro-instabilities (Charlotte and Truskinovsky, 2002; Friesecke and Theil, 2002), phonon gaps (Brillouin and Parodi, 1956; Vasseur et al., 2001), binary oscillations (Deift and McLaughlin, 1998; Vasiliev et al., 2010), dispersive overshoot (Slepyan, 1972), dispersive energy transport (Mielke, 2006; Maciá, 2004), radiative damping (Atkinson and Cabrera, 1965; Truskinovsky and Vainchtein, 2005; Carpio and Bonilla, 2005) and lattice trapping (Slepyan, 1981, 2002; Truskinovsky and Vainchtein, 2005), to mention just a few of the physical effects invisible in the CC model. Some general questions requiring a micro-scale analysis concern the dissipative closure relations and the nature of the size effect. Lattice or, more generally, micro-structural or

<footnote>* Corresponding author. Tel.: +33 5 6133 8957; fax: +33 5 6133 9095.</footnote>

<footnote>E-mail address: miguel.charlotte@isae.fr (M. Charlotte).</footnote>
anti-continuum effects are becoming dominant as one tries to miniaturize engineering devices and optimize them through
the fine material selection. It is hardly an exaggeration to say that recent technological interest has been mostly focussed
on the scales where the CC theory fails (Massobrio et al., 2010).

To capture the micro-scale effects one can directly apply the lattice MD theory which is, however, analytically opaque
and can be viewed only as a numerical tool. More importantly, adequate MD simulations of crystalline materials are in
most cases still too time consuming to be handled by modern computers. It is therefore an important challenge to develop
a hybrid theory of dynamic crystal lattices which combines the advantages of the continuum description with the ability to
reproduce some micro-structural features, for instance the nontrivial dispersion relations (Collins, 1981; DiVincenzo,
1986; Kevrekidis et al., 2002; Arndt and Griebel, 2005; Tang et al., 2006; Askes et al., 2008; Andrianov et al., 2010). Among
many attempts to construct such quasi-continuum (QC) theories (continuum theories incorporating internal length
scales)<fnref n="1" /> one can naturally distinguish two major approaches.

In the first group of theories, one extends the classical elasticity by introducing one or several nonclassical continuum-
based ingredients, such as internal variables (Mindlin, 1964; Kröner, 1968; Kunin, 1982), higher gradients (Toupin, 1962;
Krumshansl, 1965; Mindlin, 1965; Blanc et al., 2002; Sunyk and Steinmann, 2003), nonlocal interactions (Eringen, 1972;
Kunin, 1982; Rogula, 1982; Silling, 2000), etc., which are then considered operative everywhere in the domain of interest.
As one example, we can mention a static multi-field theory (Charlotte and Truskinovsky, 2008) allowing one to reduce the
discrete problem to a set of coupled second order partial differential equations. Another example is a dynamic multi-field
theory (Il’ushina, 1969; Kunin, 1982; Vasiliev et al., 2010) introducing internal variables with nontrivial inertia. The
problem with this general approach is that the new continuum structures are accounted for not only where they are
essential, but also where the simpler CC theory could be used as well.

In the second, mostly numerical approach, one applies continuum thermo-elasticity everywhere outside some small
domains (usually around singularities) and uses discrete MD inside these domains (Ortiz et al., 1996; E and Huang, 2002;
E et al., 2006; Dobson et al., 2007; Lin, 2007; Braides and Truskinovsy, 2008). Despite notable successes in recent years the
existing methods of this type still generate errors which cannot be fully controlled (Curtin and Miller, 2003; Liu et al.,
2004; Park et al., 2005; E and Li, 2005; Tadmor and Miller, 2011). More importantly, the full microscopic resolution is often
excessive even in the small domains adjacent to lattice defects.

A synthesis of the two approaches can be achieved through the systematic development of a set of nested QC theories
that, on one side, resolve a succession of increasingly fine length and time scales and, on the other side, can be matched
with each other and with the CC theory. If the finest scale of the QC sequence represents an exact interpolation of the
discrete theory, such a gradual continualization can ensure a seamless link between the classical thermo-elasticity and the
MD theory. In this way one can formulate both limiting theories in the same framework and resolve the fine scale QC
equations on numerical meshes that are coarser than atomic scales.

As a step in this direction, we consider in this paper a prototypical setting and pose the problem of interpolating a
discrete Hamiltonian particle dynamics by a QC theory with hereditary (history dependent) elasticity and generalized
inertia. The advantage of working with a drastically simplified problem (a linear particle chain with nearest neighbor
interactions) is that we know the dynamic behavior of the discrete model in full detail. However, by omitting nonlinearity
we exclude from the picture such important effects as thermalization and heat propagation. Also, by insisting on exact
interpolation and avoiding coarse graining we exclude dissipation and focus exclusively on dispersion and nonlocality. We
therefore remain in the conservative QC framework where no information is getting lost in the process of transition from
discrete to continuum.

We begin by demonstrating that in the case of localized loadings the classical elasto-dynamics fails to reproduce several
basic features of the discrete dynamics and creates artificial singularities. To remedy these shortcomings in a QC
framework one needs to find a continuum representation of the discrete elasticity operator and a consistent continuum
interpretation of point forces and initial conditions.

The most well known approach to the interpolation problem is based on the idea of projecting the dynamics onto the
first Brillouin zone (e.g. Kunin, 1982; Rogula, 1982; Eringen, 1972; Rosenau, 2003). In this way one can achieve an
isomorphism between discrete and QC theories at the expense of making the elasticity operator nonlocal in space. The
main limitation of this approach is its inability to deal with point forces and, more generally, with forces whose Fourier
spectra are not band-limited to the first Brillouin zone. Despite this inconsistency, spatially nonlocal phenomenological
theories with localized forces have been widely used in numerical simulations and are often presented as a form of a
coarse grained MD (e.g. Seleson et al., 2009).

To deal with localized and impact loadings, we propose in this paper a fundamentally different way of building a QC
interpolation of discrete dynamics. Our approach was first developed in statics (Charlotte and Truskinovsky, 2008) where
we could handle interactions of arbitrary range that are known to generate an infinite sequence of internal length scales.
The main idea is to drop the requirement of maximum smoothness of the interpolation in order to preserve the spatially
local character of the resulting QC theory. The discreteness at the continuum level is represented in this approach by a set

<footnote n="1">¹ The term *quasi-continuum* was probably first introduced by Kunin in the 1960s who used it to describe a maximally smooth continuum interpolation of a discrete model (see references to earlier papers in Kunin, 1982). The term was later used in Tadmor et al. (1996) as a general reference to a wide variety of mixed discrete-continuum numerical codes.</footnote>
of PDEs incorporating increasingly fine microcomponents of the displacement field. This allows one to deal with classical sets of boundary conditions and make the QC theory fully compatible with continuum elasticity. The nonlocal model is recovered in this framework if the microscopic components of the displacement field are eliminated by means of nonlocal relations involving the associated Green's functions.

In dynamics, which is the main focus of the present paper, we intentionally limit our description to a chain with nearest neighbors (NN) interactions where the nontrivial internal length scales disappear. As a result, the spatial discreteness in the NN theory can be adequately represented by a single displacement field. By using the Mittag-Leffler fractional decomposition of Green's function in the Fourier-Laplace space we construct a continuum description of the NN chain in terms of a single field obtaining a theory which is local in space and nonlocal in time. The salient feature of the resulting temporally nonlocal (TN) dynamic theory is that both inertial and elastic forces are nonlocal in time and are actually hardly distinguishable which has been anticipated in Willis (1981, 1997) and Milton and Willis (2007). The main advantage of such a theory is that it does not place any restrictions on the support of the applied loading and can handle point impact loadings as naturally as the original discrete theory. It is also important that the proposed general interpolation scheme can be used as the basis for constructing a sequence of successive approximations which are local in both space and time. Expectedly, the memory/hereditary structure of the approximate inertial and elastic constitutive relations becomes increasingly more complex as the order of approximation increases and finer time scales are taken into consideration.

The paper is organized as follows. In Section 1 we introduce our prototypical discrete linear system and solve explicitly the generic initial value problem. In particular, we consider a benchmark point impact test and discuss the peculiar oscillatory structure of the ensuing discrete solution. We then show that the CC theory, formulated in Section 2, is unable to reproduce some fundamental microscopic features of the energy distribution associated with this discrete solution. In Section 3 we discuss the general principles of constructing QC interpolating models. We identify two equally important problems: the approximation of the dynamic operator in the Fourier-Laplace space and choosing the regularized representation for point forces. In Section 4 we review the previously known interpolation schemes leading to spatially nonlocal theories and show that they are unable to deal with localized loadings. The new TN theory is developed in Section 5 where we show that it reproduces adequately the discrete solution of a generic point impact problem which is the most difficult test for such a theory. In Section 6 we discuss several ways of constructing spatially and temporally local approximations of the exact TN theory characterized by progressively shorter memory, all the way to CC elasticity. The last section contains our conclusions.

# 1. Discrete model

Consider an unbounded particle chain shown in Fig. 1. Suppose that the mass of a single particle is $ρa$ where a is the lattice length scale and $ρ$ is the macroscopic mass density. For simplicity we assume that each particle interacts only with its nearest neighbors through linear elastic springs with elasticity constant $αa$ where $α$ is the elastic modulus (NN chain). The time dependent displacements of particles $u(t) = (u_k(t))_{k \in \mathbb{Z}}$ are assumed to be continuous and particle collisions are excluded. We demand that the external forces $af(t) = a(f_k(t))_{k \in \mathbb{Z}}$ have appropriated growth properties in space and time. We also allow for the impact loadings defined by their impulses $a \int_{-\infty}^{t} f(\hat{t}) d\hat{t}$.

The ensuing mechanical model is fully characterized by the elastic energy, the kinetic energy and the external work:

$$ \mathcal{E}_e(\mathbf{u},t) = a \sum_{k \in \mathbb{Z}} \frac{\alpha}{2} \left[ \frac{u_k(t) - u_{k-1}(t)}{a} \right]^2, \qquad (1a) $$

$$ \mathcal{E}_k(\mathbf{u},t) = a \sum_{k \in \mathbb{Z}} \frac{\rho}{2} [D_t u_k(t)]^2, \qquad (1b) $$

$$ \mathcal{P}_r(\mathbf{f}, \mathbf{u}, t) = a \sum_{k \in \mathbb{Z}} \int_{-\infty}^{t} f_k(\hat{t}) D_t u_k(\hat{t}) d\hat{t}, \qquad (1c) $$

where the operator $D_t$ denotes appropriately generalized partial time derivative.

To make our discussion more transparent, we mostly focus on a problem with zero initial conditions

$$ \mathbf{u}(t) = 0 \quad \text{for } t \le 0 \qquad (2) $$

and only briefly discuss the general initial value problem with the given $(\mathbf{u}(0^{-}), D_t \mathbf{u}(0^{-}))$. We also assume that $f(t) = 0$ for $t \le 0$. To illustrate the difference between the discrete and the classical continuum (CC) theories we consider three tests: a step load, an impact load and an initial displacement in a single point.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Mass spring chain with nearest neighbor (NN) interactions. Masses are subjected to time dependent external forces $f_k(t)$.</figcaption></figure>
The dynamics of the NN chain is governed by the following classical system of equations (Brillouin and Parodi, 1956; Maradudin et al., 1963)

$$ \rho D_t^2 u_k = \frac{\alpha}{a^2} [u_{k+1} + u_{k-1} - 2u_k] + f_k. \quad (3) $$

Here

$$ \omega_* = \frac{2}{a} \sqrt{\frac{\alpha}{\rho}} = \frac{2c}{a} \quad (4) $$

is the characteristic frequency providing the time scale and $c = \sqrt{\alpha/\rho}$ is the macroscopic characteristic velocity.

To solve a generic problem for the linear system (3) one can use continuous Laplace (L) transform in time and a discrete Fourier (DF) transform in space. The L transform is defined as follows:

$$ \tilde{y}(t) \xrightarrow{L} \tilde{y}(\omega) = \int_{0}^{\infty} \tilde{y}(t) e^{-i\omega t} dt, \quad (5) $$

where $\Im(m(\omega)) < -\tilde{\omega}_b \le 0$ and the positive real number $\tilde{\omega}_b$ must be sufficiently large to ensure the existence of the integral. For the DF transform we write

$$ y_k \xrightarrow{\text{DF}} y(\lambda) = a \sum_{k \in \mathbb{Z}} y_k e^{-ik\lambda a}, \quad (6) $$

where $\lambda \in \mathbb{K} = [-\pi/a, \pi/a]$ and $\mathbb{K}$ represents the fundamental Brillouin's interval in the lattice reciprocal space.

By applying to (3) first the L transform and then the DF transform we obtain

$$ \Phi(\lambda, \omega)u(\lambda, \omega) = f(\lambda, \omega), \quad (7) $$

where the complex function

$$ \Phi(\lambda, \omega) = \rho[\omega_*^2 \sin^2(\lambda a/2) - \omega^2] \quad (8) $$

fully characterizes elastodynamic properties of the system.

The solution of the problem (7) defined for $t \ge 0$ is given by the following inversion formulas:

$$ u_k(t) = \frac{1}{2\pi} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} u_k(\omega) e^{i\omega t} d\omega, \quad \forall \omega_b > \tilde{\omega}_b \ge 0, \quad (9a) $$

where

$$ u_k(\omega) = \sum_{p \in \mathbb{Z}} \frac{1}{2\pi} \int_{-\pi/a}^{\pi/a} \frac{af_p(\omega)}{\Phi(\lambda, \omega)} e^{ik(p-\lambda)a} d\lambda. \quad (9b) $$

Notice that here and in what follows we imply that $u_k(t) \equiv H(t)u_k(t)$ where

$$ H(t) \stackrel{\text{def}}{=} \begin{cases} 0 & \text{if } t < 0, \\ 1 & \text{if } t \ge 0 \end{cases} $$

is Heaviside's step function.

To compute the integrals (9a) and (9b) we must first analyze the characteristic equation

$$ \Phi(\lambda, \omega) = 0. \quad (10) $$

If we partition the $\lambda$-complex domain $a^{-1}\mathbb{C}$ into simply connected strips $\cup_{q \in \mathbb{Z}} \mathbb{B}_q$, with

$$ \mathbb{B}_q = \{ \lambda = \lambda_r + i\lambda_i \in \mathbb{C}; |\lambda_r a - 2q\pi| \le \pi, (\lambda_r a, \lambda_i a) \in \mathbb{R}^2 \}, \quad (11) $$

we can solve the characteristic equation (10) in each of these strips with respect to $\omega$. The roots form an infinite sequence $\{\lambda_k(\omega)\}_{k \in \mathbb{Z}}$ implicitly given by

$$ e^{i\lambda_k a} = 1 - \frac{2\omega^2}{\omega_*^2} \pm 2i \frac{\omega}{\omega_*} \sqrt{1 - \frac{\omega^2}{\omega_*^2}}. $$

It will be convenient to arrange the different branches of the function $\{\lambda_k(\omega)\}_{k \in \mathbb{Z}}$ into doublets

$$ \lambda_{2q}(\omega) = \lambda_0(\omega) + 2q\pi/a \quad \text{and} \quad \lambda_{2q+1}(\omega) = -\lambda_{2q}(\omega), \quad (12a) $$

where we defined

$$ \lambda_0(\omega) = \frac{-i}{a} \log \left( 1 - \frac{2\omega^2}{\omega_*^2} - 2i \frac{\omega}{\omega_*} \sqrt{1 - \frac{\omega^2}{\omega_*^2}} \right) = -\frac{2}{a} \arcsin \left( \frac{\omega}{\omega_*} \right). \quad (12b) $$

The explicit relation (12b) is well-defined everywhere in $\omega_*\mathbb{C}$ outside the branch-cuts

$$ c^{\text{def}} = -\infty, -\omega_*] \cup [\omega_*, \infty[. $$
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** The dispersion relation $\lambda = \lambda_0(\omega)$ with $\omega = \omega_r + i0^\pm$. The arrows indicate the paths corresponding to increasing $\omega_r \in \omega_*\mathbb{R}$. Four branches at $\omega_r < -\omega_*$ merge into two branches at $|\omega_r| \le \omega_*$ which then split again into four branches at $\omega_r > \omega_*$. The doublet $(\lambda_{2q}, \lambda_{1-2q})$ lies inside the strip $\mathbb{B}_q$ in (11).</figcaption></figure>

Then for each $q \in \mathbb{Z}$, the doublet $(\lambda_{2q}, \lambda_{1-2q})$ lays inside the strip $\mathbb{B}_q$ in (11). The first strip $\mathbb{B}_0$ contains the fundamental Brillouin's interval $K = [-\pi/a, \pi/a]$. Inside this strip the characteristic equation has two complex roots: $\lambda_0$ and $\lambda_1 = -\lambda_0$. One period of the dispersion relation $\lambda(\omega)$ with real frequencies $\omega = \omega_r \in \omega_*\mathbb{R}$ is shown in Fig. 2.

We can now invert the discrete Fourier transform according to Eq. (9b) to obtain

$$u_k(\omega) = \sum_{p \in \mathbb{Z}} \frac{ia^2 e^{i|k-p|\lambda_0 a}}{2\alpha \sin(\lambda_0 a)} f_p(\omega),$$

where $f_p(\omega)$ is the Laplace image of the loading. Given the initial conditions (2) the inversion of the continuous Laplace transform gives

$$u_k(t) = \sum_{p \in \mathbb{Z}} \int_0^t G((k-p)a, t-\hat{t}) \frac{f_p(\hat{t})}{\rho} d\hat{t}, \qquad (13)$$

where we introduced the continuum interpolation of the discrete Green's function<fnref n="2" />

$$G(s,t) = \frac{1}{2\pi} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} \frac{i\rho a^2 e^{i|s|\lambda_0 + \omega t}}{2\alpha \sin(\lambda_0 a)} ds. \qquad (14)$$

The expression (14) can also be rewritten as

$$G(s,t) = \frac{H(t)}{2\pi\omega_*} \left\{ \begin{aligned} & \pi + 2 \int_0^{\omega_*} \frac{\sin(\omega_r t - \frac{|s|}{\ell_-})}{\omega_r \sqrt{1 - \frac{\omega_r^2}{\omega_*^2}}} d\omega_r - 2 \int_{\omega_*}^{+\infty} \frac{e^{-|s|/\ell_+} \cos(\omega_r t - \frac{|s|}{\pi/a})}{\omega_r \sqrt{\omega_r^2 - 1}} d\omega_r, \\ & = \frac{aH(t)}{\pi} \left\{ \int_0^{\pi/a} \frac{\sin(\omega_* t \sin(\lambda_r a/2)) \cos(\lambda_r S)}{\omega_* \sin(\lambda_r a/2)} d\lambda_r - \sin(|s|/\pi/a) \int_0^{+\infty} \frac{e^{-|s|/\ell_+} \sin(\omega_* t \cosh(\lambda_i a/2))}{\omega_* \cosh(\lambda_i a/2)} d\lambda_i \right\}. \end{aligned} \right. \qquad (15)$$

<footnote n="2">² More precisely, the discrete (causal) Green's function is $[\rho a]^{-1}G(s,t)$ as the loading is defined as $af$.</footnote>
where we defined two characteristic length scales $\ell_{\pm}(\omega)$

$$ \frac{a}{\ell_{-}(\omega)} = \cos^{-1}\left(1 - 2\frac{\omega^2}{\omega_{*}^2}\right) \in [0, \pi] \quad \text{for } \omega \in [-\omega_{*}, \omega_{*}], \qquad (16a) $$

$$ \frac{a}{\ell_{+}(\omega)} = \ln\left[\frac{2\omega^2}{\omega_{*}^2} - 1 + \frac{2|\omega|}{\omega_{*}}\sqrt{\frac{\omega^2}{\omega_{*}^2} - 1}\right] \geq 0 \quad \text{for } \omega \in \omega_{*}\mathbb{R}\setminus[-\omega_{*}, \omega_{*}]. \qquad (16b) $$

The discrete Green's function associated with Cauchy's problem (3), (2) solves the following equation:

$$ \rho D_t^2 G(ka,t) - \frac{\alpha}{a^2} [G(ka+a,t) + G(ka-a,t) - 2G(ka,t)] = \rho \delta_{k,0} \delta_+(t), \qquad (17) $$

and the solution of this equation can be written explicitly as

$$ G(ka,t) = \frac{2H(t)}{\pi\omega_{*}} \int_{0}^{\omega_{*}} \frac{\sin(\omega_{r}t)\cos\left(\frac{ka}{\ell_{r}}\right)}{\omega_{r}\sqrt{1-\frac{\omega_{r}^{2}}{\omega_{*}^{2}}}} d\omega_{r} = H(t) \int_{0}^{t} J_{2k}(\omega_{*}\hat{t}) d\hat{t}. \qquad (18) $$

Here $\delta_{k,p}$ denotes Kronecker's symbol, $\delta_{+}(t) = D_t H(t)$ is the "causal" Dirac's delta function (Schwartz, 1966) and

$$ J_n(\omega_*) = \frac{a}{\pi} \int_0^{\pi/a} \cos(\omega_* t \sin(\lambda_* a) - \lambda_* na) d\lambda_* \qquad (19) $$

is Bessel's function of the first type with integral order $n \in \mathbb{Z}$ (Abramowitz and Stegun, 1972). To study velocity distribution we would also need the time derivative of $G(ka,t)$ given by

$$ D_t G(ka,t) = H(t) J_{2k}(\omega_* t). $$

The behavior of the function $G(s,t) \ge 0$ is illustrated in Fig. 3(a). The asymptotic analysis based on the steepest-descend method shows that if either $\omega_* t \to \infty$ (with $|s|/a\omega_* t$ fixed) or $|s|/a \to \infty$ (with $a\omega_* t/|s|$ fixed), the function $G(s,t)$ approaches the CC limit shown in Fig. 3(b)

$$ G_c(s,t) = \omega_*^{-1} H(t-|s|/c). \qquad (20) $$

The time dependence of the functions $G(ka,t)$ and $D_t G(ka,t)$ at discrete points $s = ka \in a\mathbb{Z}$ is illustrated in Fig. 4.

**Impact test:** When applied loads change sufficiently slow in space and time one can expect that the discrete nature of the chain does not play an important role and one can safely use the CC elasto-dynamics. The problems arise when the loading is localized at a length scale comparable to $a$ or at a time scale comparable to $\omega_*^{-1}$. To illustrate the mechanism of failure of the CC theory it is instructive to consider a singular discrete loading imitating point impact:

$$ af_k(t) = \bar{p} \delta_{k,0} \delta_+(t-t_0). \qquad (21) $$

The parameter $\bar{p}$ prescribes the impulse of force applied at time $t = t_0$.

The knowledge of the response of the system to impact loading is of fundamental importance in both experimental and theoretical characterization of linear discrete systems because it allows one to fully identify the corresponding Green's functions. Indeed, the spectral representation of (21) can be written as $f(\lambda, \omega) = \bar{p}$ and if we solve the corresponding

<figure><img src="image_4.png" /><figcaption>Fig. 3. Green's functions $G(s,t)$ for the discrete model (a) and $G_c(s,t)$ for the CC model (b).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. (a) Lattice Green's function $G(ka,t)$ and (b) its time derivative $D_t G(ka,t)$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Particle motions generated by a localized impact loading. (a) Discrete model. (b) CC model. Trajectories are colored according to the magnitude of the ratio $2\rho cu_k(t)/\bar{p} \equiv \omega_* G(ka,t)$ in (a) and $2\rho cu_c(ka,t)/\bar{p} \equiv \omega_* G_c(ka,t) \equiv H(t-t_0-|s|/c)$ (see Section 2) in (b). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

discrete problem with zero initial data we obtain

$$u_k(t) = G(ka, t-t_0) \frac{\bar{p}}{\rho a}. \qquad (22)$$

To illustrate this result we show in Fig. 5(a) the particle trajectories at $t > 0$. One can distinguish two macro-domains. One, in a shape of a cone, shows oscillations, while the complementary domain is essentially oscillations free. The oscillations spread from the transition zone into the cone shaped domain where they eventually decay as the particles approach their classical macroscopic positions (see Fig. 6(a) and Section 2)

$$u_c(s,t) = H(t-t_0-|s|/c) \frac{\bar{p}}{2\rho c}. \qquad (23)$$

It is instructive to trace the evolution of the microscopic boundaries of the oscillation-generating fronts. To this end we can either trace the points $|s| = \ell_{G_+}(t-t_0)$ where the oscillations-induced overshoot reaches its maximal value or the points $|s| = \ell_{G_-}(t-t_0)$ where the displacements first reach their long time asymptotic limit $\bar{p}/(\rho\omega_*a)$ given by the CC theory (see Section 2). These two functions are illustrated in Fig. 6(b); notice that both fronts are always behind the classical Stokes' rays $|s| = c(t-t_0)$. The increasing separation of the two discrete fronts is an indication of the dispersion-induced spreading (delocalization) of the strain singularities.
<figure><img src="image_7.png" /><figcaption>Fig. 6. (a) Comparison between G(s,t) (in black) and its classical analog G_c(s,t) (in green). (b) The characteristic trajectories |s| = ℓ_G-(t-t_0), |s| = ℓ_G+(t-t_0) and the classical wave-front |s| = c(t-t_0) (with t_0 = 0). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

To better understand the energy distribution in the discrete solution one can compute the work of the loading device (1c) with the force distribution (21)

$$P_r(\mathbf{f}, \mathbf{u}, t) = \bar{p}H(t-t_0)D_t u_0(t_0^+) = H(t-t_0) \frac{\bar{p}^2}{\rho a}. \qquad (24)$$

Notice that at (s,t) = (0,0) the displacement field is continuous while the velocity field suffers a discontinuity. The instantaneous work (24) of the impact loading (21) is transformed into micro-oscillations. By a slight abuse of terminology one can say that the impact results in propagation of two “hot” rays; to render the thermodynamic language fully adequate we need sufficiently strong nonlinearity ensuring the equipartition of energy in the sense of Fermi-Pasta-Ulam system (e.g. Efendiev and Truskinovsky, 2010 and reference therein).

As another illustration, consider a spatially localized step loading in the form

$$af_k(t) = \omega_* \bar{p} \delta_{k,0} H(t-t_0). \qquad (25)$$

In this case the solution of the dynamic problem can be written as

$$u_k(t) = \frac{\omega_* \bar{p}}{\rho} [A(s,t) - A(s,t_0)], \qquad (26)$$

where we introduced the anti-derivative of Green's function from (14) which can be written in the form

$$A(s,t) = \int_{0}^{t} G(s, \hat{t}) \, d\hat{t} \equiv \frac{1}{2\pi\omega_{*}} \int_{-i\omega_{b}-\infty}^{-i\omega_{b}+\infty} \frac{(1-e^{i\omega t})e^{i|s|\lambda_{0}}}{\omega^{2}\sqrt{1-\frac{\omega_{*}^{2}}{\omega_{*}^{2}}}} \, d\omega. \qquad (27)$$

This expression can be further simplified as

$$A(s,t) = \frac{4H(t)}{\pi\omega_*} \left\{ \int_0^{\omega_*} \frac{\sin^2\left(\frac{\omega_*t}{2}\right) \cos\left(\frac{s}{\bar{p}}\right)}{\omega_r^2 \sqrt{1 - \frac{\omega_*^2}{\omega_*^2}}} d\omega_r - \sin(|s|\pi/a) \int_{\omega_*}^{+\infty} \frac{e^{-|s|/\ell_*} \sin^2\left(\frac{\omega_r t}{2}\right)}{\omega_r^2 \sqrt{\frac{\omega_*^2}{\omega_*^2} - 1}} d\omega_r \right\} \qquad (28a)$$

$$= \frac{2aH(t)}{\pi} \left\{ \int_{0}^{\pi/a} \frac{\sin^2(\omega_*t \sin(\lambda_r a/2)/2) \cos(\lambda_r s)}{\omega_*^2 \sin^2(\lambda_r a/2)} d\lambda_r - \sin(|s|\pi/a) \int_{0}^{+\infty} \frac{e^{-|s|/\ell_*} \sin^2(\omega_*t \cosh(\lambda_r a/2)/2)}{\omega_*^2 \cosh^2(\lambda_r a/2)} d\lambda_r \right\}. \qquad (28b)$$

In Fig. 7 we show the resulting displacement field and compare it with the limiting CC solution

$$u(s,t) \sim u_c(s,t) = \frac{\omega_* \bar{p}}{\rho} A_c(s,t-t_0), \qquad (29a)$$
<figure><img src="image_8.png" /><figcaption>Fig. 7. The function $A = \int_0^t G(s, \hat{t}) \, d\hat{t}$ (a) and the comparison with its CC analog $A_c = \int_0^t G_c(s, \hat{t}) \, d\hat{t}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Particle trajectories resulting from a displacement of a single particle. (a) Discrete theory with initial data: $u_k(0) = \bar{u}\delta_{k,0}$ at $k=0$. (b) CC theory with initial data: $u_c(s,0) = \kappa_2(s/a)\bar{u}$ with a triangular hat-shape function $\kappa_2$ (see (69) and Fig. 14). In each figure, the trajectories are colored according to the displacement ratio $u_k(t)/\bar{u}$ in (a) and $u_c(ka,t)/\bar{u}$ in (b). The particle trajectories plotted in Fig. 8b correspond to the displacement field plotted in Fig.13. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

where

$$A_c(s,t) = \int_0^t G_c(s, \hat{t}) \, d\hat{t} \equiv \left[ t - \frac{|s|}{c} \right] G_c(s,t) \equiv \frac{ct - |s|}{\omega_* c} H\left( t - \frac{|s|}{c} \right). \qquad (29b)$$

One can see that despite the less singular loading in this case, the deviation from the classical theory is considerable around the time of the impact and remains finite in the long time limit near the propagating deformation fronts (Figs. 7 and 8).

**Nonzero initial conditions:** In order to complete the analysis of the discrete model, it is necessary to consider the problem with nonzero initial conditions. Suppose that at $t = 0^-$ we know the displacements $\mathbf{u}(0^-)$ and the velocities $D_t \mathbf{u}(0^-)$. We also assume that $\mathbf{u}(0^+) = \mathbf{u}(0^-)$ as required by our assumption of continuity. Following Schwartz (1966, 1983), we reinterpret the nonzero initial conditions in terms of localized inertial forces

$$\rho a \mathcal{D}_t^2 \mathbf{u}(t) = \rho a [D_t \mathbf{u}(0^-) \delta_+(t) + \mathbf{u}(0^-) D_t \delta_+(t)]. \qquad (30)$$

The generalized equations incorporating these initial conditions can be rewritten in the form

$$\rho D_t^2 [H u_k] = \frac{\alpha}{a^2} H[u_{k+1} + u_{k-1} - 2u_k] + f_k + \rho \mathcal{D}_t^2 u_k. \qquad (31)$$

The additional singular loading (30) is required to move the chain instantaneously from a rest state to the configuration with the given initial displacement and velocity (Roddier, 1971). Notice that the inertial loading $\rho a \mathcal{D}_t^2 \mathbf{u}(t)$ in (30) involves
a dipole term $D_t \delta_+(t)$ whose impulse is equal to zero so that

$$ \int_{-\infty}^{t} \rho a \mathfrak{T}_t^2 \mathbf{u}(\hat{t}) d\hat{t} = H(t) \rho a D_t \mathbf{u}(0^{-}) $$

but which nevertheless contributes to the instantaneous external work

$$ \mathcal{P}_r(\rho \mathfrak{T}_t^2 \mathbf{u}, \mathbf{u}, t) = H(t) \rho a \sum_{k \in \mathbb{Z}} [D_t u_k(0^{-}) D_t u_k(0^{+}) - u_k(0^{-}) D_t^2 u_k(0^{+})]. \quad (32) $$

In spectral space the generalized dynamic problem (31) can be written as

$$ \Phi(\lambda, \omega)u(\lambda, \omega) = f(\lambda, \omega) + \rho \mathfrak{T}_t^2 u(\lambda, \omega), \quad (33a) $$

where

$$ \mathfrak{T}_t^2 u(\lambda, \omega) = i\omega u(\lambda, 0) + \dot{u}(\lambda, 0^{-}) \quad (33b) $$

and $u(\lambda, 0)$ and $\dot{u}(\lambda, 0^{-})$ denote the DF transforms of the initial displacements and velocities.

In real space the solution of the resulting Cauchy's problem can be written as (Kunin, 1982; Maradudin et al., 1963; Wolf, 1979)

$$ u_k(t) = \sum_{p \in \mathbb{Z}} \left[ \int_0^t G(ka-pa,t-\hat{t}) \frac{f_p(\hat{t})}{\rho} d\hat{t} + u_p(0) D_t \hat{G}(ka-pa,t) + \dot{u}_p(0^{-}) \hat{G}(ka-pa,t) \right], \quad (34) $$

where we introduced

$$ \hat{G}(s,t) = G(s,t) - G(s,-t) = \operatorname{sign}(t)G(s,|t|). \quad (35) $$

One can see that if in a system without body forces a single particle is given an initial velocity, the resulting motion is the same as in the case of localized impact loading.

Suppose for instance that we prescribe the initial velocity $D_t u_k(0^{-}) = \bar{v} \delta_{k,0}$, which means that $\mathfrak{T}_t^2 u_k(t) = D_t u_k(0^{-}) \delta_+(t)$. The ensuing displacement field can be written as $u_k(t) = \bar{v} \hat{G}(ka,t)$. The computation of the instantaneous work gives $\mathcal{P}_r(\rho \mathfrak{T}_t^2 \mathbf{u}, \mathbf{u}, t) = \rho a \bar{v}^2 H(t)$ which could be also obtained from (21) and (24) by the change of variables $\bar{p} = \rho a \bar{v}$ and letting $t_0 = 0$.

If, instead, a single particle is given an initial displacement, it first oscillates but then comes back to its reference position. More precisely, consider the problem with $\dot{u}_k(0^{-}) = 0$ and $u_k(0^{-}) = \bar{u} \delta_{k0}$. The solution takes the form

$$ u_k(t) = \bar{u} D_t \hat{G}(ka,t), \quad (36) $$

where the function $D_t \hat{G}(s,t)$ is illustrated in Fig. 9(a).

The function

$$ \kappa_0(s/a) = D_t G(s, 0^+) \quad (37) $$

is shown in Fig. 9(b). The corresponding particle trajectories are shown in Fig. 8(a); the displacement field is also illustrated in Fig. 12. The instantaneous work (32) can be computed explicitly $\mathcal{P}_r(\rho \mathfrak{T}_t^2 \mathbf{u}, \mathbf{u}, t) = (2\rho c^2 \bar{u}^2 / a)H(t)$.

An interesting prediction of the discrete theory for a particle chain with massless springs is that localized initial disturbances (kinematic or dynamic) propagate with infinite velocity, which means that every mass point in the infinite

<figure><img src="image_10.png" /><figcaption>Fig. 9. The function $D_t G(s,t)$ (a) its zero time limit $\kappa_0(\eta) = D_t G(\eta a, 0^+)$ (b).</figcaption></figure>
lattice feels the disturbance instantaneously (Wolf, 1979). Indeed, consider the initial data and the impact loading of the form

$$u_k(0) = \bar{u}\delta_{k,p}, \quad D_t u_k(0^-) = \bar{v}\delta_{k,p}, \quad af_k(t) = \bar{p}\delta_+(t)\delta_{k,p}. \qquad (38)$$

Then the general solution (34) can be written as

$$u_k(t) = \bar{u}D_t \hat{G}(ka-pa,t) + \left(\frac{\bar{p}}{\rho a} H(t) + \bar{v}\right) \hat{G}(ka-pa,t). \qquad (39)$$

By using (18) one can show that

$$\begin{aligned} u_k(0) &= \bar{u}\delta_{k,p}, \quad D_t u_k(0^+) = (\bar{p}/\rho a) + \bar{v}\delta_{k,p}, \\ D_t^2 u_k(0^+) &= (\omega_*^2/4)\bar{u}[\delta_{k,p+1} + \delta_{k,p-1} - 2\delta_{k,p}], \quad \text{etc.} \end{aligned}$$

One can see that the higher the order of the derivative, the larger is the number of particles affected by the impact. In this sense the concept of a sharp deformation front which is common in continuum models can be viewed in the discrete theory only as an approximation.

## 2. Classical continuum theory

The most well-known continuum approximation of the discrete dynamics is provided by the CC elasto-dynamics which replaces the lattice model by its formal long wave and low frequency approximation. In the derivation of this theory one assumes $|\lambda a| \ll 1$ and makes a crucial assumption that the kinematic field $u(s,t)$ is sufficiently smooth which naturally places restrictions on the loading field $f(s,t)$. The main outcome of the formal long wave expansion is the replacement of the discrete spectral operator $\Phi(\lambda, \omega)$ by the first term of its Taylor expansion around $\lambda = 0$

$$\Phi_c^\dagger(\lambda, \omega) = \alpha\lambda^2 - \rho\omega^2 = \alpha[\lambda^2 - \lambda_0^{*2}(\omega)], \qquad (40a)$$

where

$$\lambda_0^*(\omega) = -\omega \sqrt{\frac{\bar{p}}{\alpha}} = -\frac{\omega}{c} = -\frac{2\omega}{\omega_* a}. \qquad (40b)$$

For the limiting model one can again define the elastic energy, the kinetic energy and the external work:

$$\varepsilon_c^e(u_c,t) = \frac{\alpha}{2} \int_{-\infty}^{\infty} [D_s u_c(s,t)]^2 ds \qquad (41a)$$

$$\varepsilon_k^c(u_c,t) = \frac{\rho}{2} \int_{-\infty}^{\infty} [D_t u_c(s,t)]^2 ds \qquad (41b)$$

$$\mathcal{P}_t^c(f, u_c, t) = \int_{-\infty}^{t} \int_{-\infty}^{\infty} f(s, \tilde{t}) D_t u_c(s, \tilde{t}) ds \, d\tilde{t}. \qquad (41c)$$

The dynamic equation in the CC theory is a well known wave equation

$$\rho D_t^2 u_c(s,t) - \alpha D_s^2 u_c(s,t) = f(s,t). \qquad (42)$$

Here the operator $D_s$ represents appropriately generalized partial space derivative.

To solve the continuum problem (42) we need to introduce the Integral Fourier (IF) transform and its inverse

$$y(s) \xrightarrow{\text{IF}} y^\dagger(\lambda) = \int_{-\infty}^{\infty} y(s)e^{-i\lambda s} d\lambda, \qquad (43a)$$

$$y^\dagger(\lambda) \rightarrow y(s) = \frac{1}{2\pi} \int_{-\infty}^{\infty} y^\dagger(\lambda)e^{i\lambda s} ds. \qquad (43b)$$

By applying first the L transform and then the IF transform, we obtain the following algebraic equation:

$$\Phi_c^\dagger(\lambda, \omega) u_c^\dagger(\lambda, \omega) = f^\dagger(\lambda, \omega). \qquad (44)$$

After inverting the IF transform we obtain

$$u_c(s, \omega) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \frac{f^\dagger(\lambda, \omega)}{\Phi_c^\dagger(\lambda, \omega)} e^{is\lambda} d\lambda. \qquad (45)$$

Finally, by inverting the L transform we obtain the real space solution

$$u_c(s,t) = \frac{1}{2\pi} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} u_c(s, \omega) e^{i\omega t} d\omega = \int_0^t \int_{-\infty}^{\infty} G_c(s-\hat{s}, t-\hat{t}) \frac{f(\hat{s}, \hat{t})}{\rho a} d\hat{s} \, d\hat{t}. \qquad (46)$$
Recall that here we again imply that $H(t)u_c(s,t) = u_c(s,t)$ which means in particular that $u_c(s,0^-) = 0$ and $D_t u_c(s,0^-) = 0$. In (46) we introduced

$$G_c(s,t) = \frac{i\rho a}{4\pi\tau} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} \frac{e^{i[(\omega t + s)/\tau_0^*]}}{\tau_0^*(\omega)} d\omega = \frac{a}{2c} H(t-s/c), \quad (47)$$

which is Green's function of the classical continuum elasticity describing D'Alembert's solution (see Fig. 3(b)). It solves Cauchy's problem for

$$\rho D_c^2 G_c(s,t) - \alpha D_s^2 G_c(s,t) = \rho a \delta(s) \delta_+(t) \quad (48)$$

with the initial conditions $G_c(s,t) = 0$ for $t \le 0$. Here $\delta(s) = D_s^2|s|/2$ is the “acausal” Dirac's delta function (Schwartz, 1966).

In addition to replacing the discrete elasticity operator by its long-wave classical analog, the continuum elastodynamics description must propose a continuum representation $f(s,t)$ of the discrete loading $f(t) = \{f_k(t)\}_{k \in \mathbb{Z}}$. While it is clear that at any instant $t$ the resultant force must be preserved

$$\int_{-\infty}^{\infty} f(s,t) ds = a \sum_{k \in \mathbb{Z}} f_k(t), \quad (49)$$

the detailed continuum representation of the localized discrete forces is non-unique. In what follows we consider two most natural candidates and demonstrate that the predictions of the classical model depend on the chosen mode of representation of the loading.

**Impact test: singular loading:** The most straightforward choice for the continuum distribution representing $f(t)$ is the singular loading of the form

$$f(s,t) = a \sum_{k \in \mathbb{Z}} f_k(t) \delta(s-ka). \quad (50)$$

In this case the discrete impact loading (21) should be interpreted in the classical continuum theory as

$$f(s,t) = \sum_{k \in \mathbb{Z}} \bar{p} \delta_{k,0} \delta_{+}(t-t_0) \delta(s-ka). \quad (51)$$

The main problem with the representation (50) is that it violates the long wave character of the Taylor approximation involved in the construction of classical elastodynamics. This crucial inconsistency is often ignored because the dynamic problem (42) can still be formally solved.

Indeed, if we assume that the generalized function $f(s,t)$ in (50) takes the form (51) we obtain

$$u_c(s,t) = G_c(s,t-t_0) \frac{\bar{p}}{\rho a} = H(t-t_0-|s|/c) \frac{\bar{p}}{2\rho c}. \quad (52)$$

We have already illustrated the function $u_c(s,t)$ in Fig. 5(b) and our Fig. 6(a) provides further comparison of the discrete and continuum solutions. One can see that no oscillations occur in the classical elasto-dynamic solution. Outside the cone-shaped domain bounded by Stokes's rays $|s| = c(t-t_0)$ the displacements are equal to zero while inside this domain the displacements are constant.

Notice that, in contrast to the discrete case, the continuous displacement has a singularity at $(s,t) = (0,t_0)$. This singularity produces two discontinuities traveling in opposite directions with velocities $c$ and $-c$. As a consequence, the velocity field $D_t u_c$ becomes sufficiently singular for the integral (41c) to diverge at each $t \ge t_0$.

Indeed, the impact loading (51) can be viewed as the limit of a regularized sequence (Schwartz, 1966, 1983)

$$f^c(s,t) = \sum_{k \in \mathbb{Z}} a f_k^c(t) \delta(s-ka), \quad (53)$$

where

$$a f_k^c(t) = \bar{p} \delta_{k,0} \frac{H(t-t_0)-H(t-t_0-\epsilon)}{\epsilon},$$

and $\epsilon > 0$ is a small time scale. An additional regularization in space does not change the result as shown in Appendix A.

By solving the CC problem with the regularized loading (53) we obtain

$$u_c^c(s,t) = \frac{\bar{p}}{2\rho c} \left\{ \frac{c(t-t_0)-|s|}{c\epsilon} \left[ H(t-t_0-|s|/c) - H(t-t_0-\epsilon-|s|/c) \right] + H(t-t_0-\epsilon-|s|/c) \right\}. \quad (54)$$

At finite $\epsilon$ the displacement field is continuous, however, the function $u_c^c$ becomes singular as $\epsilon \to 0^+$, approaching the limit field $u_c$ given by (52). In contrast, the regularized velocity field is discontinuous

$$\dot{u}_c^c(s,t) = \frac{\bar{p}}{2\rho c\epsilon} [H(t-t_0-|s|/c)-H(t-t_0-\epsilon-|s|/c)], \quad (55)$$
while the work $\mathcal{P}_r^c(f^c, u_c^c, t)$ (41c) is finite

$$ \mathcal{P}_r^c(f^c, u_c^c, t) = \int_0^t \frac{\bar{p}^2}{2\rho c \epsilon^2} [H(\hat{t}-t_0)-H(\hat{t}-t_0-\epsilon)] d\hat{t} = H(t-t_0)\bar{p} \frac{u_c^c(0,t)-u_c^c(0,t_0)}{\epsilon}. $$

If we now take into account that $u_c^c(0,t_0) = 0$, the first term of the asymptotic expansion $\epsilon \to 0^+$ can be written as

$$ \mathcal{P}_r^c(f^c, u_c^c, t) \sim H(t-t_0) \frac{\bar{p} u_c(0, t_0^+)}{\epsilon} = \frac{1}{\epsilon \omega_*} \mathcal{P}_r(\mathbf{f}, \mathbf{u}, t). $$

An immediate consequence of this formula is that in classical elasto-dynamics the work of the loading device required to generate the displacement field (52) is infinite. Instead, in the discrete theory, which the CC theory is expected to approximate, the corresponding work was finite and equal to $\mathcal{P}_r(\mathbf{f}, \mathbf{u}, t)$ (see (24)). The problem is that in CC theory we apply a finite force impulse to an effective particle with zero mass and, as a result, it starts moving with an infinite velocity. In the discrete theory the mass of the particle receiving the force impulse is finite which leads to finite force increment and finite work of the loading device.

It is of interest to trace the evolution of the infinite energy supplied by the loading device in the continuum problem. Since none of the two macro-domains (unperturbed and translated) carries any energy, all the energy must be concentrated in the infinitely thin transition zones $|s| = c(t-t_0)$. The work of the loading device is therefore localized in the singularities of the displacement field propagating in both directions with the characteristic speed $c$. This energy localization is another spurious artifact of the long wave approximation because the latter is incompatible with the singular structure of impact loading.

**Impact test: regularized loading:** To correct the inconsistency between the long wave character of the elasticity operator and the short wave structure of the loading, one can try to modify the continuum representation of the localized forces. The interpolation condition can still be satisfied if we restrict the Fourier image of the discrete system of forces to the first Brillouin zone. More precisely, we can use the following construction (Krumshansl, 1965; Eringen, 1972; Kunin, 1982; Rogula, 1982):

$$ f^{\dagger}(\lambda, \omega) = 1_{\mathbb{B}_0}(\lambda)f(\lambda, \omega), \qquad (56) $$

where $1_{\mathbb{B}_0}$ is the characteristic function of the domain $\mathbb{B}_0$. In the physical space we can write

$$ f(s,t) = \sum_{k \in \mathbb{Z}} f_k(t) \kappa_1(s/a-k), \qquad (57) $$

where the weight function is now regular (see Fig. 10)

$$ \kappa_1(\eta) = \mathrm{sinc}(\eta) = \frac{\sin(\pi\eta)}{\pi\eta}. \qquad (58) $$

One can show that the continuum distribution of forces (57) provides the smoothest (Whittaker-Shannon-Kotel'nikov or WSK) interpolation of the discrete system of loads $f(t)$. Observe that while in the case of the singular loads (50) the Fourier spectrum is completely de-localized, in the case of the regularized forces (57) it is fully contained in the strip $\mathbb{B}_0$. In the physical space this results in delocalization of the loading at the scale $a$ and brings an internal length scale into the continuum theory.

To solve the continuum problem (42) with the source term given by (57) it is convenient to apply first the **IF** transform and then the **L** transform. This again brings an algebraic problem (44)

$$ \Phi_c^*(\lambda, \omega) u_c^*(\lambda, \omega) = f^*(\lambda, \omega). $$

<figure><img src="image_11.png" /><figcaption>Fig. 10. The function $\kappa_1(\eta) = \mathrm{sinc}(\eta)$.</figcaption></figure>
After inverting the L transform we obtain

$$u_c^\dagger(\lambda, t) = \frac{1}{2\pi} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} \frac{f^\dagger(\lambda, \omega)}{\Phi_c^\dagger(\lambda, \omega)} e^{i\omega t} d\omega, \qquad (59)$$

which can be rewritten as

$$u_c^\dagger(\lambda, t) \equiv a_1 \mathbb{E}_0(\lambda) \sum_{p \in \mathbb{Z}} H(t) \int_0^t \frac{2 \sin(\omega_* (t-\hat{t}) \lambda a/2) f_p(\hat{t}) e^{-i\lambda p a}}{\omega_* \lambda a} d\hat{t}.$$

Next, we invert the IF transform to obtain

$$u_c(s,t) = \sum_{p \in \mathbb{Z}} \int_0^t G_c^*(s-pa,t-\hat{t}) \frac{f_p(\hat{t})}{\rho} d\hat{t}, \qquad (60)$$

where

$$G_c^*(s,t) = H(t) \int_0^{\pi/a} \frac{2 \sin(\omega_* t \lambda a/2) \cos(\lambda s)}{\omega_* \lambda \pi} d\lambda \equiv H(t) \int_0^t \frac{\kappa_1 \left( \frac{s+c\hat{t}}{a} \right) + \kappa_1 \left( \frac{s-c\hat{t}}{a} \right)}{2} d\hat{t}. \qquad (61)$$

In Fig. 11 we compare the kernels $G_c(ka,t)$ and $G_c^*(ka,t)$. One can see that the regularization of the force allows one to recover the lattice scale oscillatory structure of the discrete solution without any additional improvements at the level of the operator which remains scale-free.

We can now consider the impact test with the loading (21) to obtain

$$u_c(s,t) = G_c^*(s,t-t_0) \frac{\bar{p}}{\rho a}. \qquad (62)$$

The computation of the work gives

$$P_r^c(f, u_c, t) = H(t-t_0) \bar{p} D_t u_c(0, t_0^+) = H(t-t_0) \frac{\bar{p}^2}{\rho a} D_t G_c^*(0, 0^+) = P_r(f, u, t). \qquad (63)$$

As we see, now the work of the loading device is finite and is exactly equal to the analogous term in the discrete theory.

Despite this success in matching the energy, the point-wise agreement between the discrete solution (22) and its continuum analog (62) is at most qualitative. In particular, while continuum loading interpolates the discrete one, the continuum displacement field does not. This is related to the fact that the discrete elasticity operator $\Phi(\lambda, \omega)$ is represented in classical elastodynamics rather crudely.

<figure><img src="image_12.png" /><figcaption>Fig. 11. Comparison of the discrete Green's function $G(ka,t)$, the CC Green's function $G_c(ka,t)$ and the pseudo-classical approximation $G_c^*(ka,t)$ for $(k,t) \in \mathbb{N} \times \omega_*^{-1}\mathbb{R}^+$.</figcaption></figure>
**Nonzero initial conditions:** Suppose now that in the CC theory we pose the problem with initial data $u_c(s,0^{-})$ and $\dot{u}_c(s,0^{-})$. The corresponding spectral problem can be written as

$$ \Phi_c^t(\lambda, \omega)u_c(\lambda, \omega) = f(\lambda, \omega) + i\omega\rho u_c(\lambda, 0) + \rho\dot{u}_c(\lambda, 0^{-}), \quad (64) $$

where we introduced the **IF** transforms of the initial data. In real space we obtain D'Alembert's solution (e.g. Wolf, 1979)

$$ u_c(s,t) = H(t) \int_0^t \int_{-\infty}^{\infty} G_c(s-\hat{s}, t-\hat{t}) \frac{f(\hat{s}, \hat{t})}{\rho a} d\hat{s} d\hat{t} + \frac{u_c(s+ct, 0)+u_c(s-ct, 0)}{2} + \frac{1}{2c} \int_{s-ct}^{s+ct} \dot{u}_c(\hat{s}, 0) d\hat{s}, \quad (65) $$

which is continuous as long as the loading $f$ does not involve any singular impact component.

By following the pattern developed in the discrete theory, we can again interpret the nonzero initial conditions in the CC setting as an instantaneous distribution of inertial forces

$$ \rho a \mathfrak{D}_t^2 u_c(s,t) = \rho a \{ D_t u_c(s,0^{-}) \delta_+(t) + u_c(s,0^{-}) D_t \delta_+(t) \}. \quad (66) $$

We then need to solve the problem (Schwartz, 1966, 1983)

$$ \rho D_t^2 [Hu_c] - \alpha H D_s^2 u_c = f + \rho \mathfrak{D}_t^2 u_c. \quad (67) $$

The correspondence between the localized initial data, $u(0)$ and $D_t u(0^{-})$, and the continuously distributed initial conditions $u_c(s,0^{-})$ and $D_t u_c(s,0^{-})$ can be assumed in the following general form:

$$ u_c(s, 0^{-}) = \sum_{k \in \mathbb{Z}} \kappa(s/a-k), \quad D_t u_c(s, 0^{-}) = \sum_{k \in \mathbb{Z}} D_t u_k(0^{-}) \kappa(s/a-k), \quad (68) $$

where we introduced the weight functions $\kappa(\eta)$ to be specified later; notice that the weight functions used in the displacement field $u_c(s,0^{-})$ may be different from those used in the velocity field $D_t u_c(s,0^{-})$.

In the example with an impact load we have seen that the equality between the discrete and the continuum instantaneous work terms can be achieved when the force distribution is represented by the WSK relations with $\kappa = \kappa_1(\eta)$. A direct inspection of the exact solution of the discrete problem seemingly suggests another weight function $\kappa_0(\eta)$ defined in (37) and depicted in Fig. 9(b). In addition to these specially designed weight functions $\kappa_0$ and $\kappa_1$ one can also use more formal regularization of the point loads as it is done in the standard FE codes. For instance, it is customary to consider polynomial shape functions, in particular, the simplest representative of this class

$$ \kappa_2(\eta) = \frac{|\eta+1| + |\eta-1| - 2|\eta|}{2} \quad (69) $$

is often used in the numerical quasi-continuum codes (Tadmor et al., 1996; Miller and Tadmor, 2002).

To illustrate different choices, consider the CC analog of the problem with initial displacement of a single mass point which means that $u_k(0) = \bar{u}\delta_{kp}$ and $D_t u_k(0) = 0$. Then

$$ u_c(s, 0) = \sum_{k \in \mathbb{Z}} \kappa(s/a-k)u_k(0) = \kappa(s/a)\bar{u} $$

and the resulting motion takes the form

$$ u_c(s,t) = \frac{\kappa((s+ct)/a) + \kappa((s-ct)/a)}{2} \bar{u}. \quad (70) $$

In Fig. 12 we show the displacement field obtained by using the weight function $\kappa = \kappa_0$ and compare it with the interpolation of the corresponding exact solution of the discrete problem $u(s,t) = \bar{u}D_t \hat{G}(s,t)$ (36). Our Fig. 13 shows the dynamics in the CC theory resulting from the initial point displacement regularized as $u_c(s,0) = \kappa_2(s/a)\bar{u}$.

While various choices of weight functions ensure the continuity of the field $u_c(s,t)$, they do not necessarily guarantee the interpolation property $u_c(ka,t) = u_k(t)$. Also neither the local and global impulse balances

$$ \rho a \mathfrak{D}_t^2 u_c(ka, t) = \rho a \mathfrak{D}_t^2 u_k(t), \quad (71) $$

$$ \int_{-\infty}^{\infty} \rho \mathfrak{D}_t^2 u_c(s,t) ds = \sum_{k \in \mathbb{Z}} \rho a \mathfrak{D}_t^2 u_k(t), \quad (72) $$

nor the work balance

$$ \mathcal{P}_r^c(\rho \mathfrak{D}_t^2 u_c, u_c, t) = \mathcal{P}_r(\rho \mathfrak{D}_t^2 \mathbf{u}, \mathbf{u}, t), \quad (73) $$

where

$$ \mathcal{P}_r^c(\rho \mathfrak{D}_t^2 u_c, u_c, t) = H(t) \int_{-\infty}^{\infty} \rho(D_t u_c(s, 0^{-}) D_t u_c(s, 0^{+}) - u_c(s, 0^{-}) D_t^2 u_c(s, 0^{+})) ds $$

are automatic. The local impulse equivalence (71) is satisfied if the interpolation condition $\kappa(k) = \delta_{k,0}$ is satisfied as it is the case for $\kappa_1$ and $\kappa_0$. The global impulse equivalence (72) requires that the weight functions satisfy the interpolation condition $\int_{-\infty}^{\infty} \kappa(s) ds = a$ as in the case of $\kappa_1$. However, in the case when $u_k(0) = \bar{u}\delta_{kp}$ and $D_t u_k(0) = 0$, the weight function $\kappa_1$ does not satisfy (73) because $\mathcal{P}_r^c(\rho \mathfrak{D}_t^2 u_c, u_c, t) = (\pi^2/6)\mathcal{P}_r(\rho \mathfrak{D}_t^2 \mathbf{u}, \mathbf{u}, t)$. The weight function $\kappa_2$ satisfies (71) and (72);
<figure><img src="image_13.png" /><figcaption>Fig. 12. Displacement fields resulting from different continuum approximations of the discrete initial data $u_k(0) = \bar{u}\delta_{k,0}$ and $\dot{u}_k = 0$. Discrete displacement field $u(s,t) = \bar{u}D_tG(s,t)$ corresponding to the initial displacement $u(s,0) = \kappa_0(s/a)\bar{u} = D_tG(s,0^+)\bar{u}$ (in black) and the continuum solution $u_c(s,t) = \bar{u}[\kappa_0((s+ct)/a) + \kappa_0((s-ct)/a)]/2$ corresponding to the same initial displacement $u_c(s,0) = \kappa_0(s/a)\bar{u}$ (in green). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. The solution of CC equations corresponding to initial regularized displacement field $u_c(s,0) = \kappa_2(s/a)\bar{u}$ and zero initial velocity.</figcaption></figure>

however, it satisfies (73) only in the present case of prescribed displacement

$$P_r^c(\rho \mathfrak{D}_t^2 u_c, u_c, t) = \frac{2\rho c^2 \bar{u}^2}{a} H(t) = P_r(\rho \mathfrak{D}_t^2 \mathbf{u}, \mathbf{u}, t).$$

The comparison of Fig. 13 and Fig. 12 shows that in the case of $\kappa_2(k)$ the agreement with the discrete theory is rather poor. The inadequacy of the popular weight function $\kappa_2(k)$ becomes even more visible if we consider the problem with localized initial velocity $D_t u_k(0^-) = \bar{v}\delta_{k,0}$ and approximate $D_t u_c(s,0^-) = \bar{v}\kappa_2(s/a)$. Then $u_c(s,t) = \bar{v}\hat{G}_c^{**}(s,t-t_0)$ where $\hat{G}_c^{**}(s,t) = \text{sign}(t)G_c^{**}(s,|t|)$ and

$$G_c^{**}(s,t) = H(t) \int_0^t \frac{\kappa_2((s+c\tau)/a) + \kappa_2((s-c\tau)/a)}{2} d\tau.$$
The computation of the work gives $P_r(\rho \mathfrak{T}_t^2 u_c, u_c, t) = \frac{2}{3} P_r(\rho \mathfrak{T}_t^2 \mathbf{u}, \mathbf{u}, t)$, which means that one-third of the initial energy is lost in this approximation. Similar results can be obtained for the $\kappa_2$-regularized impact loading $f(s,t) = f_k(t)\kappa_2(s/a)$. In this case the displacement field is $u(s,t) = (\bar{p}/\rho a)\hat{G}_c^*(s,t-t_0)$ and the work is again $P_r^c(f, u_c, t) = \frac{2}{3}P_r(\mathbf{f}, \mathbf{u}, t)$.

In the following sections we consider a set of QC models providing more complete force and energy equivalences between discrete and continuum theories.

### 3. Continuum interpolation

To achieve micro-scale consistency between continuum and discrete theories it is natural to look for a QC interpolation of lattice dynamics. The task is to replace the discrete sequence of functions, $u_k(t)$, representing a response of the system to loads, $f_k(t)$, by a continuous function of two variables $u(s,t)$ satisfying

$$u(ka,t) = u_k(t) \qquad (74)$$

and solving a continuum problem with distributed loads $f(s,t)$ such that

$$f(ka,t) = f_k(t). \qquad (75)$$

From the spectral viewpoint, the crucial step is to replace the exact equation (7) in the Fourier-Laplace space, which corresponds to the problem with trivial initial data, by an interpolating equation

$$\Phi^{\dagger}(\lambda, \omega) u^{\dagger}(\lambda, \omega) = f^{\dagger}(\lambda, \omega). \qquad (76)$$

Here the new functions$^3$ $\Phi^\dagger(\lambda, \omega)$ and $f^\dagger(\lambda, \omega)$ replace the periodic in $\lambda$ functions $\Phi(\lambda, \omega)$ and $f(\lambda, \omega)$. After the function $u^\dagger(\lambda, \omega)$ is recovered as a solution of the algebraic problem (76), the interpolating continuum displacement field $u(s,t)$ in real space can be obtained from the following inversion formulas:

$$u(s,t) = \frac{1}{2\pi} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} u^\dagger(s, \omega) e^{i\omega t} d\omega, \qquad (77a)$$

where

$$u^{\dagger}(s, \omega) = \frac{1}{2\pi} \int_{-\infty}^{\infty} u^{\dagger}(\lambda, \omega) e^{i\lambda s} d\lambda = \frac{1}{2\pi} \int_{-\infty}^{\infty} \frac{f^{\dagger}(\lambda, \omega) e^{i\lambda s}}{\Phi^{\dagger}(\lambda, \omega)} d\lambda. \qquad (77b)$$

In some cases it may be more convenient to first invert the L transform and then the IF transform. Whatever is the choice, the function $u(s,t)$ must be a solution of the following nonlocal (in both space and time) equation:

$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \mathcal{L} \left( \frac{s-\tilde{s}}{a}, t-\tilde{t} \right) u(\tilde{s}, \tilde{t}) \, d\tilde{s} \, d\tilde{t} = f(s,t), \qquad (78a)$$

where

$$\mathcal{L}(\eta, t) = \frac{1}{4\pi^2} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} \int_{-\infty}^{\infty} \Phi^{\dagger}(\lambda, \omega) e^{i(\omega t + \eta \lambda a)} d\lambda \, d\omega. \qquad (78b)$$

While looking at the force balance equation (78a) one is tempted to distinguish between the inertial and the elastic forces and rewrite (78a) in the form

$$I_1[u] - I_2[u] = f, \qquad (79a)$$

where

$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \mathcal{L} \left( \frac{s-\tilde{s}}{a}, t-\tilde{t} \right) u(\tilde{s}, \tilde{t}) \, d\tilde{s} \, d\tilde{t} = I_1[u](s,t) - I_2[u](s,t). \qquad (79b)$$

Although below we present several cases when the introduction of the functionals $I_1[u]$ and $I_2[u]$ is natural, in the general theory the inertia-elasticity splitting appear to be rather arbitrary.

To ensure that the continuous field $u(s,t)$ (with $s \in a\mathbb{R}$) interpolates the discrete field $u_k(t)$ (with $k \in \mathbb{Z}$) one needs to check that (77b) and (9b) agree at $s = ka \in a\mathbb{Z}$ for all $t \in \omega_*^{-1}\mathbb{R}$. The general interpolation condition (74) may be rewritten as

$$\int_{-i\omega_b - \infty}^{-i\omega_b + \infty} \int_{-\infty}^{\infty} \frac{f^{\dagger}(\lambda, \omega) e^{i(\lambda ka + \omega t)}}{\Phi^{\dagger}(\lambda, \omega)} d\lambda \, d\omega \equiv \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} \int_{-\pi/a}^{\pi/a} \frac{f(\lambda, \omega) e^{i(\lambda ka + \omega t)}}{\Phi(\lambda, \omega)} d\lambda \, d\omega, \qquad (80)$$

<footnote>$^3$ Starting from this section, we denote by $y^\dagger(\lambda, \cdot)$ the IF transform of any continuum field $y(s, \cdot)$ when it is necessary to distinguish it from the DF transform $y(\lambda, \cdot)$ of a corresponding discrete field $\mathbf{y}(\cdot) \stackrel{\text{def}}{=} \{y_k(\cdot)\}_{k \in \mathbb{Z}}$.</footnote>
where the integration order can also be reversed. We recall that

$$f(\lambda, \omega) = a \sum_{k \in \mathbb{Z}} f_k(\omega) e^{-ika\lambda} \quad \text{with} \quad f_k(\omega) = \int_0^\infty f_k(t) e^{-i\omega t} dt$$

and view Eq. (80) as a restriction on both $f^\dagger(\lambda, \omega)$ and $\Phi^\dagger(\lambda, \omega)$.

As we have already seen in the previous section, an important step is to specify the interpolation method for the collection of point forces $f_k(t)$ (and for the discrete initial data, if necessary). The two natural choices for the body forces have been already mentioned in the previous section: the singular loading

$$f(s,t) = \sum_{k \in \mathbb{Z}} f_k(t) \delta(s/a - k), \qquad (81)$$

and the WSK regularized loading

$$f(s,t) = \sum_{k \in \mathbb{Z}} f_k(t) \mathrm{sinc}(s/a - k). \qquad (82)$$

It is clear that both representations ensure the equivalence of the resultants (49).

Now, since the force distribution $f_k(t)$ may be arbitrary, the interpolation property can be reformulated as a restriction on the elasticity operator $\Phi^\dagger(\lambda, \omega)$ only. For instance, in the case of singular loading (81) we obtain

$$\int_{-\infty}^{\infty} \frac{e^{i\lambda ka}}{\Phi^\dagger(\lambda, \omega)} d\lambda = \int_{-\pi/a}^{\pi/a} \frac{e^{i\lambda ka}}{\Phi(\lambda, \omega)} d\lambda. \qquad (83)$$

In the case of regularized loading (82) we have instead

$$\int_{-\pi/a}^{\pi/a} \frac{e^{i\lambda ka}}{\Phi^\dagger(\lambda, \omega)} d\lambda = \int_{-\pi/a}^{\pi/a} \frac{e^{i\lambda ka}}{\Phi(\lambda, \omega)} d\lambda. \qquad (84)$$

One can view Eqs. (83) and (84) as a set of (weak) restrictions on the inverse elasticity operator $[\Phi^\dagger(\lambda, \omega)]^{-1}$ (Green's function).

Below we present in two separate sections the analysis of two types of QC interpolating theories. One of them, which leads to spatially nonlocal and temporally local problem, requires nonlocal boundary conditions and can deal adequately only with regularized loadings (82). The other one operates with spatially local and temporally nonlocal equations, requires only conventional boundary conditions and is compatible with both regularized and singular loadings (81). Both types of theories ensure the energetic equivalence of discrete and continuum models by preserving the value of the external work

$$\int_0^t \int_{-\infty}^{\infty} f(s, \tilde{t}) D_t u(s, \tilde{t}) ds d\tilde{t} = a \sum_{k \in \mathbb{Z}} \int_0^t f_k(\hat{t}) D_t u_k(\hat{t}) d\hat{t}. \qquad (85)$$

## 4. Spatially nonlocal QC models

Assume first that the point forces appearing in the lattice model are interpolated in the continuum framework by the WSK interpolants (82). This case can be handled by two known models which differ only in details. Throughout most of this section we assume that $u(s,t) = 0$ for $t \le 0$, at the end of this section we briefly discuss a nontrivial initial value problem.

**Model with unlimited Fourier spectrum:** The most straightforward way to satisfy condition (84) is to assume (Eringen, 1972; Kunin, 1982; Rogula, 1982)

$$\Phi^\dagger(\lambda, \omega) = \Phi(\lambda, \omega). \qquad (86)$$

This strong equivalence assumption yields the following dynamical equation:

$$\rho D_t^2 u(s,t) - \alpha \frac{u(s+a,t) + u(s-a,t) - 2u(s,t)}{a^2} = f(s,t), \qquad (87)$$

where the separation into inertial and elastic terms mentioned in (79b) is straightforward

$$I_1[u](s,t) = \rho D_t^2 u(s,t), \qquad (88a)$$

$$I_2[u](s,t) = \alpha \frac{u(s+a,t) + u(s-a,t) - 2u(s,t)}{a^2}. \qquad (88b)$$

The corresponding nonlocal problem (78) has the kernel

$$\mathcal{L}(\eta, t) = \frac{\rho}{a} \delta(\eta) D_t^2 \delta_+(t) - \frac{\alpha}{a^3} \delta_+(t) D_\eta^2 \kappa_2(\eta), \qquad (89)$$
where $\kappa_2(\eta)$ was introduced in (69) and

$$D_{\eta}^{2}\kappa_{2}(\eta) = \delta(\eta+1) + \delta(\eta-1) - 2\delta(\eta).$$

The above relation allows one to rewrite (88b) in the form

$$I_2[u](s,t) = \frac{\alpha}{a} \int_{-\infty}^{\infty} \kappa_2((s-\tilde{s})/a) D_s^2 u(\tilde{s}, t) d\tilde{s}.$$

In this framework the elastic energy, the kinetic energy and the external work take the form (Eringen, 1972; Rosenau, 2003; Silling, 2000; Seleson et al., 2009)

$$\varepsilon_e^*(u,t) = \frac{\alpha}{2} \int_{-\infty}^{\infty} \left[ \frac{u(s,t) - u(s-a,t)}{a} \right]^2 ds, \qquad (90a)$$

$$\varepsilon_k^*(u,t) = \frac{\rho}{2} \int_{-\infty}^{\infty} [D_t u(s,t)]^2 ds, \qquad (90b)$$

$$\mathcal{P}_r^*(f, u, t) = \int_0^t \int_{-\infty}^{\infty} f(s, \tilde{s}) D_t u(s, \tilde{s}) ds \, d\tilde{s}. \qquad (90c)$$

For general continuum fields ($u(s,t), f(s,t)$) interpolating the discrete fields ($u(t), f(t)$) we have

$$\mathcal{P}_r^*(f, u, t) \notin \mathcal{P}_r(\mathbf{f}, u, t), \quad \varepsilon_e^*(u, t) \ge \varepsilon_e(\mathbf{u}, t) \quad \text{and} \quad \varepsilon_k^*(u, t) \ge \varepsilon_k(\mathbf{u}, t).$$

These inequalities turn into equalities for the WKS class of loadings and displacement fields whose support is confined to $\mathbb{K}$.

We observe that while the kinetic energy and the external work are both local, the elastic energy in this model is (strongly) nonlocal in space. This type of models have been recently widely used for numerical simulation of solid materials with discontinuous displacements (Silling, 2000; Seleson et al., 2009).

Interestingly, spatial nonlocality can be “moved” from elastic to kinetic energy if one uses instead of the field $u(s,t)$ a different field $\tilde{U}(s,t)$ with (Rosenau, 2003; Truskinovsky and Vainchtein, 2006)

$$u(s,t) = \int_{-\infty}^{\infty} \mathcal{M}(s-\tilde{s}) \tilde{U}(\tilde{s}, t) d\tilde{s}. \qquad (91)$$

The nonlocal operator $\mathcal{M}$ may be defined in different ways; for instance, one can adopt the following implicit characterization:

$$aD_s \tilde{U}(s,t) = u(s+a/2,t) - u(s-a/2,t). \qquad (92)$$

At the spectral level this corresponds to the following explicit assumption:

$$\mathcal{M}^{\dagger}(\lambda) = \frac{\lambda a}{2 \sin(\lambda a/2)}. \qquad (93)$$

The resulting spatially local equations for the new variable $\tilde{U}(s,t)$ are derived in Appendix B.

**Model with $\mathbb{B}_0$-band limited Fourier spectrum:** In the class of loadings and displacement fields whose spectral support is localized in $\mathbb{K}$, it is natural to replace (86) by (Krumshansl, 1965; Kunin, 1982; Rogula, 1982)

$$\Phi^{\dagger}(\lambda, \omega) = 1_{\mathbb{B}_0}(\lambda) \Phi(\lambda, \omega). \qquad (94)$$

This is another strong equivalence assumption which establishes a formal one-to-one mapping between the QC problem defined in the WSK interpolation space and the original discrete problem (Kunin, 1982). The associated nonlocal problem in the physical space (see (78)) is characterized by the kernel

$$\mathcal{L}(\eta, t) = \frac{\rho}{a} \mathrm{sinc}(\eta) D_t^2 \delta_+(t) - \frac{\alpha}{a^3} \delta_+(t) D_\eta^2 \kappa_3(\eta), \qquad (95)$$

where

$$\kappa_3(\eta) = \int_{-\infty}^{\infty} \kappa_1(\eta - \hat{\eta}) \kappa_2(\hat{\eta}) d\hat{\eta} = (\eta + 1) \int_{0}^{\eta+1} \mathrm{sinc}(\eta) d\eta + (\eta - 1) \int_{0}^{\eta-1} \mathrm{sinc}(\eta) d\eta - 2\eta \int_{0}^{\eta} \mathrm{sinc}(\eta) d\eta - \frac{4}{\pi^2} \cos(\eta\pi). \qquad (96)$$

The second derivative of the function $\kappa_3(\eta)$ can be written as

$$D_{\eta}^{2}\kappa_{3}(\eta) = -4\frac{\eta^{2}-1/2}{\eta^{2}-1}\mathrm{sinc}(\eta).$$
For this model, the inertial contribution to the force balance $I_1[u]$ in (79b) remains as in (88a), while the elastic contribution takes the form

$$I_2[u](s,t) = \frac{\alpha}{a} \int_{-\infty}^{\infty} \kappa_3((s-\check{s})/a) D_s^2 u(\check{s}, t) \, d\check{s}. \qquad (97a)$$

For comparison with the classical continuum approach we mention that similarly regularized CC model with

$$\Phi^{\dagger}(\lambda, \omega) = 1_{\mathbb{B}_0}(\lambda) \Phi_c^{\dagger}(\lambda, \omega)$$

is equivalent in the physical space to a nonlocal problem of the type (42) with the kernel

$$\mathcal{L}_c^*(\eta, t) = \frac{\rho}{a} \mathrm{sinc}(\eta) D_t^2 \delta_+(t) - \frac{\alpha}{a^3} \delta_+(t) D_\eta^2 \kappa_1(\eta), \qquad (98)$$

where the second derivative of $\kappa_1(\eta)$ in (58) is given by

$$D_\eta^2 \kappa_1(\eta) = \frac{(2-\pi^2\eta^2)\mathrm{sinc}(\eta) - 2\cos(\pi\eta)}{\eta^2}.$$

In such a hybrid model the inertial functional $I_1[u](s,t)$ in (79b) remains like in (88a), while the elastic contribution to the force balance takes the form

$$I_2[u](s,t) = \frac{\alpha}{a} \int_{-\infty}^{\infty} \kappa_1((s-\check{s})/a) D_s^2 u(\check{s}, t) \, d\check{s}. \qquad (99a)$$

Our Fig. 14 presents a comparison of the three elastic kernels $\kappa_1$, $\kappa_2$, and $\kappa_3$.

The main limitation of the QC theory (94) is the restriction on the space of shape functions whose spectral support must be contained in $\mathbb{K}$ and which must be regular. In particular, the functions describing localized impact are automatically excluded. The model based on (86) (but not on (94)) can be in principle used in a broader space; however, in the case of singular loading and nonzero initial conditions, both models (86) and (94) produce artifacts (see Appendix C).

To illustrate the inner working of this class of models, consider the response of the system to a general regularized discrete loading in the form (82). In this case the two models presented in this section are equivalent. The solution of (76) can be easily found if we first invert the L transform and then the IF transform (in the previous sections the inversion was done in an opposite order). We obtain

$$u^{\dagger}(\lambda, t) = 1_{\mathbb{K}}(\lambda) \sum_{p \in \mathbb{Z}} \int_{0}^{t} G_{\mathbb{K}}^{\dagger}(\lambda, t-\hat{t}) \frac{f_{p}(\hat{t}) e^{-i \lambda p a}}{\rho} d\hat{t} \qquad (100)$$

$$u(s,t) = \sum_{p \in \mathbb{Z}} \int_{0}^{t} G_{\mathbb{K}}(s-pa, t-\hat{t}) \frac{f_{p}(\hat{t})}{\rho} d\hat{t}. \qquad (101)$$

Here we introduced the “pseudo” Green's function

$$G_{\mathbb{K}}(s,t) = \frac{1}{2\pi} \int_{-\pi/a}^{\pi/a} G_{\mathbb{K}}^{\dagger}(\lambda,t) e^{i\lambda s} d\lambda \quad \text{where} \quad G_{\mathbb{K}}^{\dagger}(\lambda,t) = \frac{a}{2\pi} \int_{-i\omega_b-\infty}^{-i\omega_b+\infty} \frac{e^{i\lambda\omega t}}{\omega_*^2 \sin^2(\lambda a/2) - \omega^2} d\omega.$$

<figure><img src="image_15.png" /><figcaption>Fig. 14. Comparison of the three weight functions $\kappa_1$ in (58), $\kappa_2$ in (69), and $\kappa_3$ in (96).</figcaption></figure>
A straightforward computation gives

$$G_K^\dagger(\lambda, t) = H(t) \frac{a \sin(\omega_* t \sin(\lambda a/2))}{\omega_* \sin(\lambda a/2)}$$

and

$$G_K(s,t) = H(t) \int_0^{\pi/a} \frac{a \sin(\omega_* t \sin(\lambda a/2)) \cos(\lambda s)}{\omega_* \sin(\lambda a/2) \pi} d\lambda = \sum_{q \in \mathbb{Z}} \text{sinc}(s/a-q) G(qa,t). \quad (102)$$

An important consequence of the last formula is that $D_t G_K(s, 0^+) = \text{sinc}(s/a) = \kappa_1(s/a)$ (see Fig. 10) which is different from the function $D_t G(s, 0^+) = \kappa_0(s/a)$ shown in Fig. 9.

**Impact test:** We are now in a position to present the solution of a localized impact problem under the assumption that it has been WSK-regularized. By direct substitution of (21) into (101), we obtain

$$u(s,t) = \sum_{q \in \mathbb{Z}} \text{sinc}(s/a-q) G_K(qa, t-t_0) \frac{\bar{p}}{\rho a} = G_K(s,t-t_0) \frac{\bar{p}}{\rho a}. \quad (103)$$

The solution (103) is illustrated in Fig. 15. Notice that $G_K(s,t) \ge 0$ at the lattice points $s \in a\mathbb{Z}$, which agrees with the fact that $G_K(ka,t) = G(ka,t) \ge 0$. In Fig. 15(b) we show the nonphysical oscillations originating from the replacement of the discrete model by our spatially nonlocal QC interpolation.

**Nonzero initial conditions:** Consider now a nontrivial initial value problem in the model with band limited spectra. According to WKS interpolation hypothesis (Kunin, 1982) the continuous initial conditions must be written as

$$u(s,0) = \sum_{k \in \mathbb{Z}} u_k(0) \text{sinc}(s/a-k) = \sum_{k \in \mathbb{Z}} u_k(0) \kappa_1(s/a-k), \quad (104a)$$

$$\dot{u}(s,0^{-}) = \sum_{k \in \mathbb{Z}} \dot{u}_k(0^{-}) \text{sinc}(s/a-k) = \sum_{k \in \mathbb{Z}} \dot{u}_k(0^{-}) \kappa_1(s/a-k). \quad (104b)$$

As before, we can represent the initial conditions as singular inertial forces. Then the equation of motion (79a) takes the form

$$I_1[u] - I_2[u] = f + \rho \mathfrak{D}_t^2 u, \quad (105)$$

where we recall that in this class of models $I_1[u] = D_t^2[Hu]$ and

$$\rho a \mathfrak{D}_t^2 u(s,t) = \rho a \sum_{k \in \mathbb{Z}} \mathfrak{D}_t^2 u_k(t) \text{sinc}(s/a-k). \quad (106)$$

The solution of this initial value problem can be written as

$$u(s,t) = \sum_{k \in \mathbb{Z}} \left[ \int_0^t G_K(s-ka, t-\hat{t}) \frac{f_k(\hat{t})}{\rho} d\hat{t} + u_k(0) D_t \hat{G}_K(s-ka, t) + \dot{u}_k(0^{-}) \hat{G}_K(s-ka, t) \right], \quad (107)$$

where

$$\hat{G}_K(s,t) = G_K(s,t) - G_K(s, -t) = \text{sign}(t) G_K(s, |t|). \quad (107)$$

<figure><img src="image_16.png" /><figcaption>Fig. 15. (a) The pseudo-Green's function $G_K(s,t)$ for the QC model. (b) The difference $G_K(s,t)-G(s,t)$.</figcaption></figure>
Due to the known properties of the WKS interpolation, this way of representing the discrete initial data ensures the continuity of $u(s,t)$, guarantees pointwise equivalences between continuum and discrete external loads $f(s,t)$ and $f(t)$, and preserves the work as required by (85); such an equivalence also holds for the inertial loads $\rho a \mathbb{T}_i^2 u(s,t)$ and $\rho a \mathbb{T}_i^2 u(t)$.

The main problem with the two interpolation schemes described in this section is that the corresponding models are spatially nonlocal which raises the issue of nontrivial boundary conditions and requires consideration of extended interaction layers replacing sharp interfaces between the bodies (Kunin, 1982). In the next section we discuss how the spatial nonlocality can be systematically circumvented.

## 5. New temporally nonlocal (TN) model

In the case of generic loadings whose spectrum $f(\lambda, t)$ is not confined to the interval $\mathbb{K}$, the continuum interpolation proposed in the previous section fails (see Appendix C) and one has to look for the alternatives. In this section we present a dynamic extension of the static theory developed in Charlotte and Truskinovsky (2008) which is free of the shortcomings discussed above. Since we are no longer dealing with regularized loading, our construction will be based on the interpolation condition (83).

We begin by rewriting the spectral Green's function of the discrete model $[\Phi(\lambda, \omega)]^{-1}$ as a Mittag-Leffler's fractional decomposition by the poles in the complex plane $\lambda$ at fixed frequency $\omega$

$$[\Phi(\lambda, \omega)]^{-1} = 2 \sum_{q \in \mathbb{Z}} \frac{\lambda_{2q}(\omega)}{\Phi'(\lambda_{2q}(\omega), \omega)} \frac{1}{\lambda^2 - \lambda_{2q}^2(\omega)}, \quad (109)$$

where $\Phi'(\lambda, \omega) = \partial\Phi(\lambda, \omega)/\partial\lambda$. Following Charlotte and Truskinovsky (2008) we can now truncate the expansion (109) by including only the roots of the characteristic equation located inside the $\mathbb{B}_0$-strip. We then obtain a new quasi-continuum model with elasticity operator

$$[\Phi^{\dagger}(\lambda, \omega)]^{-1} = \frac{2\lambda_0}{\Phi'(\lambda_0, \omega)} \frac{1}{\lambda^2 - \lambda_0^2} = \frac{a}{\alpha \sin(\lambda_0 a)} \frac{1}{\lambda^2 - \lambda_0^2}. \quad (110)$$

For each value of $\omega \in \omega_* \mathbb{C} \backslash \mathbb{C}$, the spectral function $[\Phi^\dagger(\lambda, \omega)]^{-1}$ has the same residues in the domain $\mathbb{B}_0$ as the spectral function of the original discrete model. This ensures that the interpolating condition (83) is satisfied. Notice that the elasticity operator of the new theory

$$\Phi^{\dagger}(\lambda, \omega) = \frac{\alpha[\lambda^2 - \lambda_0^2]\sin(\lambda_0 a)}{\lambda_0 a} \quad (111)$$

is no longer periodic in $\lambda$ and that the resulting model is spatially local. Since the Mittag-Leffler decomposition could be performed only in one complex dimension $\lambda$, the remaining dependence on $\omega$ is nontrivial which means strong temporal (time) nonlocality.

While the TN operator (111) appears to be very different from the original spectral operator $\Phi(\lambda, \omega)$ of the discrete theory (8), which is local in time and nonlocal in space, both models converge to the CC model $\Phi_c^\dagger(\lambda, \omega)$ in the limit of low frequencies ($|\omega| \ll \omega_*$) and long wavelengths ($|\lambda|^{-1} \gg a$). More precisely, if $(\omega_*^{-1}, a) \to (0, 0)$ while the product $\omega_* a = 2c$ remains finite and $|\lambda_*^2 a| = 2|\omega/\omega_*| \ll 1$ then $\lambda_0 \to \lambda_0^*$ and $\Phi^\dagger(\lambda, \omega) \to \Phi_c^\dagger(\lambda, \omega)$. We also recall that $\Phi(\lambda_0, \omega) = 0$, $\Phi^\dagger(\lambda_0, \omega) = 0$ and $\Phi_c^\dagger(\lambda_0, \omega) = 0$ for all $\omega \in \omega_* \mathbb{C} \backslash \mathbb{C}$.

In the physical space the dynamic equations of the TN theory take the form

$$-\frac{\alpha}{2\pi} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} \frac{[D_s^2 + \lambda_0^2]\sin(\lambda_0 a)}{\lambda_0 a} u(s, \omega) e^{i\omega t} d\omega = f(s, t), \quad (112)$$

where $f(s,t)$ is defined as in (50). The inertial and elastic terms can be separated as in (79a) with

$$I_1[u](s,t) = \frac{\alpha}{\pi} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} \frac{\lambda_0 \omega}{a\omega_*} \sqrt{1 - \frac{\omega^2}{\omega_*^2}} u(s, \omega) e^{i\omega t} d\omega, \quad (113a)$$

$$I_2[u](s,t) = -\frac{\alpha}{\pi} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} \frac{\omega}{a\lambda_0 \omega_*} \sqrt{1 - \frac{\omega^2}{\omega_*^2}} D_s^2 u(s, \omega) e^{i\omega t} d\omega. \quad (113b)$$

Notice that, the integration in (113) is performed over the measure generated by the phonon frequency distribution function

$$v(\omega) = - \frac{a}{\pi} \frac{d\lambda_0}{d\omega}(\omega) = \frac{2}{\pi \omega_* \sqrt{1 - \frac{\omega_*^2}{\omega^2}}}$$

which is normalized according to $\int_0^{\omega_*} v(\omega) d\omega = 1$ (Kunin, 1982).
To compute the integrals in (113) we need to first rewrite their integrands as a combination of terms satisfying the conditions of Jordan's integration lemma. We obtain

$$I_1[u](s,t) = \frac{1}{2\pi} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} [i\omega]^2 \{i\omega \gamma_1(\omega)\} \left\{ \left[1 - \frac{\omega^2}{\omega_*^2}\right] u(s, \omega) \right\} e^{i\omega t} d\omega \quad (114a)$$

$$I_2[u](s,t) = \frac{1}{2\pi} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} [i\omega \gamma_2(\omega)] \left\{ \left[ 1 - \frac{\omega^2}{\omega_*^2} \right] D_s^2 u(s, \omega) \right\} e^{i\omega t} d\omega, \quad (114b)$$

where

$$\gamma_1(\omega) \stackrel{\text{def}}{=} \frac{2i\alpha\lambda_0 a}{a^2 \omega_* \omega^2 \sqrt{1 - \frac{\omega^2}{\omega_*^2}}} = \frac{i\rho\omega_*\lambda_0 a}{2\omega^2 \sqrt{1 - \frac{\omega^2}{\omega_*^2}}}, \quad (115a)$$

$$\gamma_2(\omega) \stackrel{\text{def}}{=} \frac{-2\alpha}{i\omega_*\lambda_0 a \sqrt{1 - \frac{\omega^2}{\omega_*^2}}} = \frac{2i\alpha}{\omega_*\lambda_0 a \sqrt{1 - \frac{\omega^2}{\omega_*^2}}}. \quad (115b)$$

By mapping these two kernels into the real space we obtain

$$\tilde{\gamma}_p(t) = \frac{1}{2\pi} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} \gamma_p(\omega)e^{i\omega t} d\omega = \tilde{\gamma}_p(\omega_*)t \quad \text{for } p=1,2. \quad (116)$$

The function $\gamma_1(t) = \tilde{\gamma}_1(\omega_*)t$ can be interpreted as a mass density kernel while the function $\gamma_2(t) = \tilde{\gamma}_2(\omega_*)t$ can be viewed as an elastic kernel. The integrals in (116) can be written as

$$\gamma_1(t) = H(t)\rho \left\{ 1 - \int_{\omega_*}^{\infty} \frac{\omega_* \cos(\omega_r t)}{\omega_r^2 \sqrt{\frac{\omega_*^2}{\omega_r^2} - 1}} d\omega_r \right\}, \quad (117a)$$

$$\gamma_2(t) = H(t)\alpha \left\{ 1 - \int_{\omega_*}^{\infty} \frac{4\cos(\omega_r t) d\omega_r}{\omega_* \left[ \frac{a^2}{\lambda_r^2} + \pi^2 \right] \sqrt{\frac{\omega_r^2}{\omega_*^2} - 1}} \right\}. \quad (117b)$$

By using a different path in Laplace's space we can obtain another representation which is easier to handle

$$\gamma_1(t) = H(t) \frac{\rho a^2}{2\pi} \left\{ \int_0^{\pi/a} \frac{\lambda_r \sin(\sin(\lambda_r a/2)\omega_*t)}{\sin^2(\lambda_r a/2)} d\lambda_r - \int_0^{\infty} \frac{\lambda_r \sin(\cosh(\lambda_r a/2)\omega_*t)}{\cosh^2(\lambda_r a/2)} d\lambda_r \right\}, \quad (118a)$$

$$\gamma_2(t) = H(t) \frac{2\alpha}{\pi} \left\{ \int_0^{\pi/a} \frac{\sin(\sin(\lambda_r a/2)\omega_*t)}{\lambda_r} d\lambda_r + \int_0^{\infty} \frac{\lambda_r \sin(\cosh(\lambda_r a/2)\omega_*t)}{\lambda_r^2 + \frac{\pi^2}{a^2}} d\lambda_r \right\}. \quad (118b)$$

<figure><img src="image_17.png" /><figcaption>Fig. 16. The behavior of the functions $1-\gamma_p(t)/\gamma_p(\infty) \equiv 1-\tilde{\gamma}_p(\omega_*)t/\tilde{\gamma}_p(\infty)$ for $p=1,2$.</figcaption></figure>
Our Fig. 16 shows the transient behavior of the effective mass density $\tilde{Y}_1$ and the effective elastic modulus $\tilde{Y}_2$ as they approach their long-time asymptotes

$$ \tilde{Y}_1(+\infty) = \rho \quad \text{and} \quad \tilde{Y}_2(+\infty) = \alpha. \qquad (119) $$

Observe that both kernels $Y_p(t)$ are positive. On can also obtain from (117) and (118) the following useful relations:

$$ Y_1(t) = 0, \quad Y_2(t) = 0, \quad D_t Y_1(t) = 0, \quad D_t^2 Y_1(t) = 0 \quad \text{for } t \le 0. \qquad (120) $$

To express the integrals in (113) more transparently, we assume temporarily that the following initial condition is satisfied:

$$ u(s,t) = 0 \quad \text{for } t \le 0. \qquad (121) $$

We also exclude impact loads, so that $u(s,t^+) = u(s,t^-)$ and $D_t u(s,t^+) = D_t u(s,t^-)$. Then, we can write for $t > 0$

$$ I_1[u](s,t) = \left[ 1 + \frac{1}{\omega_*^2} D_t^2 \right] \int_0^t u(s,t-\hat{t}) D_t^3 Y_1(\hat{t}) \, d\hat{t}, \qquad (122a) $$

$$ I_2[u](s,t) = \left[ 1 + \frac{1}{\omega_*^2} D_t^2 \right] D_s^2 \int_0^t u(s,t-\hat{t}) D_t Y_2(\hat{t}) \, d\hat{t}. \qquad (122b) $$

Substituting expressions (122) in the dynamic equation (79a) we obtain

$$ \left[ 1 + \frac{1}{\omega_*^2} D_t^2 \right] \int_0^t \left[ u(s,t-\hat{t}) D_t^3 Y_1(\hat{t}) - D_s^2 u(s,t-\hat{t}) D_t Y_2(\hat{t}) \right] d\hat{t} = f(s,t). \qquad (123) $$

A unique solution of Eq. (123) that satisfies the initial conditions (121) can be written in the form

$$ u(s,t) = \int_{0}^{t} \int_{-\infty}^{\infty} G(s-\hat{s}, t-\hat{t}) \frac{f(\hat{s}, \hat{t})}{\rho a} d\hat{s} d\hat{t} \qquad (124) $$

where Green's function $G(s,t)$ is given by (15). The appearance in this formula of the continuous analog of the discrete causal Green's function (18) means that the TN theory provides the most natural continuous interpolation of the solution for the original discrete problem (13).

**Impact test:** An important property of the TN model (123) is that it can adequately handle singular loadings even though the formal derivation of Eq. (123) seems to have excluded the impact loads. Indeed, consider a regular sequence of loadings in (124) converging to classical impact in the limit $(\beta, \epsilon) \to (0^+, 0^+)$

$$ f^{\beta, \epsilon}(s,t) = \frac{H(t)-H(t-\epsilon)}{\epsilon} \frac{\bar{p}}{\beta a} \kappa_i(s/\beta a) \to \bar{p} \delta_+(t) \delta(s) \quad \text{with } i=1,2. $$

The corresponding displacement field in (124) can be shown to have a limit

$$ u^{\beta, \epsilon}(s,t) = \int_{0}^{t} \int_{-\infty}^{\infty} G(s-\hat{s}, t-\hat{t}) \frac{f^{\beta, \epsilon}(\hat{s}, \hat{t})}{\rho a} d\hat{s} d\hat{t} \to \frac{\bar{p}}{\rho a} G(s,t). $$

In the distribution theory sense, the function $G(s,t)$ in (15) represents the fundamental solution of Cauchy's initial value problem

$$ I_1[G](s,t) - I_2[G](s,t) = \rho a \delta_+(t) \delta(s), \qquad (125a) $$

$$ G(s,t) = 0 \quad \text{for } t \le 0. \qquad (125b) $$

In particular, the particle displacement field in the benchmark test (51) can be written in the form

$$ u(s,t) = G(s,t-t_0) \frac{\bar{p}}{\rho a}. \qquad (126) $$

This continuous solution interpolates the discrete solution (22) and therefore the TN theory preserves all the details of the discrete solution.

If instead of spatially singular loads (81), we apply spatially regularized loads (82), the response in (126) remains relatively close to the solution of the discrete problem. The exact solution in this case can be written as in (101) with the pseudo-Green function kernel $G_K$ replaced by

$$ G^*(s,t) = \int_{-\infty}^{\infty} G(s-\eta a, t) \kappa_1(\eta) \, d\eta. \qquad (127) $$

A quantitative comparison of the kernels $G$, $G_K$, $G_c$, $G_c^*$ and $G^*$ computed along the classical Stokes' rays is shown in Fig. 17. One can see that, as expected, the discrete kernel $G(ka, t)$ (discrete solution) is interpolated by both $G_K(s,t)$ (solution of spatially nonlocal QC theory) and $G(s,t)$ (solution of the TN theory) while the functions $G_c$ (classical solution for singular impact problem) and $G_c^*$ (classical solution for regularized impact problem) clearly stand aside. In contrast to the classical elasto-dynamic solution $G_c$, the pseudo-kernel $G_c^*$ exhibits oscillations but still remains rather far from the discrete solution. Another pseudo-kernel $G^*$ (solution of the TN model for regularized impact) while being clearly different from $G_K$
<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Wave-front values of G, G<sub>K</sub>, G*, G<sub>c</sub> and G<sub>c</sub>* (respectively given in (15), (102), (127), (47) and (61)) along the CC wave-front set |s| = c(t-t<sub>0</sub>). Red circles corresponds to the discrete points ω<sub>a</sub>t = 2 s/a ∈ Z where G(ct,t) = G<sub>K</sub>(ct,t). After a transient phase the functions G<sub>K</sub>(s,t), G(s,t) and G*(s,t) converge to the same asymptotic limit (0.3 × ω<sub>a</sub><sup>-1</sup>). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

near the impact source, converges relatively fast to it; this shows the robustness of the temporally nonlocal (TN) model even when the singular loading is approximated by the regularized loading. Most important, however, is that the response of the TN theory (110) to singular loading agrees with the discrete theory without introducing any spurious oscillations. This is achieved due to a nontrivial memory structure in this continuum model which imposes additional smoothness on the time variations of the displacement field.

**Nonzero initial conditions:** As we have seen, the TN model is characterized by complex internal inertia and nontrivial rate dependent internal elasticity. However, the model remains spatially local and requires only two boundary conditions. It is also helpful that exactly as the original discrete model, the TN model can be viewed as a theory requiring only two initial conditions.

Indeed, consider a general discrete problem with initial conditions **u**(0) and D<sub>t</sub>**u**(0<sup>-</sup>). In contrast to the CC model and the continuous interpolations studied in the previous section, in TN model we do not need to choose the continuum interpolation of the initial conditions **u**(s,0<sup>-</sup>) and D<sub>t</sub>**u**(s,0<sup>-</sup>) and can instead use directly the lattice values

$$u(ka, 0^{\pm}) = u_k(0^{-}), \quad D_t u(ka, 0^{-}) = D_t u_k(0^{-}). \qquad (128)$$

We can again write the dynamic equation (123) while taking into account the inertial forces

$$\left[ 1 + \frac{1}{\omega_*^2} D_t^2 \right] \int_0^t [u(s,t-\hat{t}) D_t^3 \gamma_1(\hat{t}) - D_s^2 u(s,t-\hat{t}) D_t \gamma_2(\hat{t})] d\hat{t} = f(s,t) + \rho \mathfrak{T}_t^2 u(s,t). \qquad (129)$$

Here we introduced

$$\rho a \mathfrak{T}_t^2 u(s,t) = \rho a \sum_{k \in \mathbb{Z}} \mathfrak{T}_t^2 u_k(t) \delta(s/a-k). \qquad (130)$$

The spectral form of the equation (129) is

$$f^{\dagger}(\lambda, \omega) + \rho \mathfrak{T}_t^2 u^{\dagger}(\lambda, \omega) = \Phi^{\dagger}(\lambda, \omega) \tilde{u}^{\dagger}(\lambda, \omega), \qquad (131a)$$

where

$$\rho \mathfrak{T}_t^2 u^{\dagger}(\lambda, \omega) = \rho [D_t u(\lambda, 0^{-}) + i\omega u(\lambda, 0)]. \qquad (131b)$$

The solution of (131a) for the singular loads (81) reads

$$u(s,t) = \sum_{p \in \mathbb{Z}} \left[ \int_{0}^{t} G(s-pa,t-\hat{t}) \frac{f_{p}(\hat{t})}{\rho} d\hat{t} + D_{t} \hat{G}(s-pa,t) u_{p}(0) + \hat{G}(s-pa,t) \dot{u}_{p}(0^{-}) \right], \qquad (132)$$
<figure><img src="image_19.png" /><figcaption>Fig. 18. TN and CC displacement fields resulting from smooth continuum representations of the discrete initial data: $u_k(0) = \bar{u} \delta_{k,0}$ and $\dot{u}_k(0) = 0$. The TN solution $u = D_t \hat{G}^{**}$ $\bar{u}$ (in black) and the CC solution $u_c(s,t) = \bar{u}[\kappa((s+ct)/a) + \kappa((s-ct)/a)]/2$ (in green) for the cases $\kappa(\eta) = \kappa_1(\eta)$ in (a) and $\kappa(\eta) = \kappa_2(\eta)$ in (b). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

where the functions $G$ and $\hat{G}$ are defined in (15) and in (35), respectively. One can check that the continuum representation of the initial conditions takes exactly the same form as in the original discrete theory (34)

$$u(s,0) = \sum_{p \in \mathbb{Z}} D_t G(s-pa, 0^+) u_p(t=0) = \sum_{p \in \mathbb{Z}} \kappa_0(s/a-p) u_p(0), \qquad (133)$$

$$D_t u(s, 0^{-}) = \sum_{p \in \mathbb{Z}} D_t G(s - pa, 0^{+}) \dot{u}_p(0^{-}) = \sum_{p \in \mathbb{Z}} \kappa_0(s/a - p) \dot{u}_p(0^{-}). \qquad (134)$$

The correspondence between fundamental solutions ensures an equivalence between the work in the discrete model $\mathcal{P}_t(\mathbf{f}, \mathbf{u}, t)$ and in the TN model $\mathcal{P}_t^c(\mathbf{f}, \mathbf{u}, t)$.

To highlight the role of singular inertia loads $\rho a \mathfrak{D}_t^2 u$ in (130), one can try to regularize them by assuming that

$$\rho a \mathfrak{D}_t^2 u(s,t) = \rho a \sum_{k \in \mathbb{Z}} \mathfrak{D}_t^2 u_k(t) \kappa(s/a-k), \qquad (135)$$

where we use regular spatial weight function $\kappa(\eta)$ instead of the singular Dirac function $\delta(\eta)$. The solution of the corresponding TN problem (129) can be written as

$$u(s,t) = \sum_{p \in \mathbb{Z}} \left[ \int_{0}^{t} G(s-pa,t-\hat{t}) \frac{f_p(\hat{t})}{\rho} d\hat{t} + D_t \hat{G}^{**} (s-pa,t) u_p(0) + \hat{G}^{**} (s-pa,t) \dot{u}_p(0^{-}) \right], \qquad (136)$$

where the pseudo-kernel

$$\hat{G}^{**}(s,t) = G^{**}(s,t) - G^{**}(s, -t) = \operatorname{sign}(t) G^{**}(s, |t|) \qquad (137)$$

is based on function $G^{**}(s,t)$ given by

$$G^{**}(s,t) = \int_{-\infty}^{\infty} G(s-\eta a, t) \kappa(\eta) d\eta. \qquad (138)$$

In Fig. 18 we show the displacement field (136) for two different continuum representations of the discrete initial data $u_k(0) = \bar{u} \delta_{k,0}$ and $\dot{u}_k(0) = 0$. One representation is based on the weight function $\kappa_1$ (see Fig. 18(a)) and another one on $\kappa_2$ (see Fig. 18(b)). Both CC solutions shown in Fig. 18 satisfy the condition $u_c(ka, 0) = u_k(0)$, and conditions (71) and (72), however, none of them satisfies the kinematic equivalence condition (74). The TN solution $u = D_t \hat{G}^{**} \bar{u}$ computed with the regularized representation of the initial conditions similarly does not satisfy any of the above conditions, in particular, we lose the kinematic equivalence (74). However, the TN solution (136) can still be shown to approach the interpolating solution with singular initial data (132) at sufficiently large times $t > 0$ (see Fig. 12).

# 6. Alternative formulations and approximations

The integro-differential equation of the TN theory (123) is rather complicated and it is natural to try to simplify it. Below we propose several ways of rewriting this equation which makes its structure more transparent. We again assume
for simplicity that the initial conditions are trivial (121). The external load $f(s,t)$ with support in $[t_0, \infty[$ can be arbitrary, however, in this section we exclude point impact.

Under these conditions the general solution of the discrete problem can be interpolated by the function (124)

$$u(s,t) = \int_{t_0}^{t} \int_{-\infty}^{+\infty} D_t A(s-\hat{s}, t-\hat{t}) \frac{f(\hat{s}, \hat{t})}{\rho a} d\hat{s} d\hat{t}, \qquad (139)$$

which we have presented in terms of the anti-derivative $A(s,t)$ from (28). As a benchmark test we consider here a spatially localized step loading (25)

$$f(s,t) = a \sum_{k \in \mathbb{Z}} f_k(t) \delta(s-ka) \quad \text{with} \quad af_k(t) = \omega_* \bar{p} \delta_{k,0} H(t-t_0), \qquad (140)$$

which produces the displacement field

$$u(s,t) = \frac{\omega_* \bar{p}}{\rho} [A(s,t) - A(s,t_0)]. \qquad (141)$$

In this section we compare this exact solution with solutions generated by several modifications of the general TN theory.

First we notice that by taking into account (120) and the assumed initial conditions, we can rewrite (123) in the form

$$f(s,t) = D_t^3 \int_{t_0}^{t} \gamma_1(t-\hat{t}) \hat{u}_0(s, \hat{t}) d\hat{t} - D_s^2 D_t \int_{t_0}^{t} \gamma_2(t-\hat{t}) \hat{u}_0(s, \hat{t}) d\hat{t} - D_s^2 \int_{t_0}^{t} D_t \gamma_1(t-\hat{t}) \hat{u}_0(s, \hat{t}) d\hat{t}, \qquad (142)$$

where we introduced a new displacement field

$$\hat{u}_0(s,t) = \left[ 1 + \frac{1}{\omega_*^2 D_t^2} \right] u(s,t). \qquad (143)$$

The nonlocal relation between $\hat{u}_0(s,t)$ and $u(s,t)$ can be explicitly inverted

$$u(s,t) = H(t-t_0) \int_{t_0}^{t} \frac{\sin(\omega_*(t-\hat{t}))}{\omega_*} \hat{u}_0(s, \hat{t}) d\hat{t}. \qquad (144a)$$

The solution $\hat{u}_0(s,t)$ can be written as

$$\hat{u}_0(s,t) = \int_{t_0}^{t} \int_{-\infty}^{+\infty} D_t A_0(s-\hat{s}, t-\hat{t}) \frac{f(\hat{s}, \hat{t})}{\rho a} d\hat{s} d\hat{t}, \qquad (144b)$$

where we defined the new kernel

$$\begin{aligned} A_0(s,t) &= \frac{1}{2\pi\omega_*} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} \frac{\sqrt{1-\frac{\omega_*^2}{\omega_*^2}(1-e^{i\omega_*t})e^{i\lambda_0|s|}}}{\omega_*^2} d\omega \\ &= \frac{aH(t)}{\pi} \left\{ \int_0^{\pi/a} \frac{2 \sin^2\left(\frac{\omega_* \sin(\lambda_r a/2)}{2}\right) \cos(\lambda_r s)}{\omega_*^2 \tan^2(\lambda_r a/2)} d\lambda_r \right. \\ & \qquad \left. + \sin(|s|\pi/a) \int_0^{+\infty} \frac{2e^{-\lambda_i|s|}}{\omega_*^2} \tanh^2(\lambda_r a/2) \sin^2\left(\frac{\omega_* t \cosh(\lambda_i a/2)}{2}\right) d\lambda_i \right\}. \end{aligned} \qquad (145a)$$

In the case of step loading the solution $\hat{u}_0$ can be written in the form (141), modulo replacement of $A$ by $A_0$. As one can see in Fig. 19, this solution captures the main macroscopic phenomenon, however, the recovery of the original displacement field $u$ brings additional lattice scale oscillations generated by the nonlocal operator (144a).

An interesting feature of the new representation is that both momentum and stress have now similar relation to the displacement field. Indeed we can rewrite (142) as

$$D_t p - D_s \sigma = f. \qquad (146)$$

Here the stress $\sigma$ and the linear momentum $p$ are defined by the classical constitutive relations $p = \rho v^*$ and $\sigma = \alpha \varepsilon^*$, where the generalized strain

$$\varepsilon^*(s,t) = D_s \int_{t_0}^{t} \tilde{\gamma}_2^{-1}(+\infty) D_t \gamma_2(t-\hat{t}) \hat{u}_0(s, \hat{t}) d\hat{t} \qquad (147)$$

and the generalized velocity

$$v^*(s,t) = D_t \int_{t_0}^{t} \tilde{\gamma}_1^{-1}(+\infty) D_t \gamma_1(t-\hat{t}) \hat{u}_0(s, \hat{t}) d\hat{t} \qquad (148)$$

are represented as strikingly similar nonlocal functions of the displacement field. Such a ‘blending’ of stress–strain relations and momentum–velocity relations is typical for the homogenized description of dynamic media with micro-structure (Willis, 1981, 1997; Milton and Willis, 2007).
<figure><img src="image_20.png" /><figcaption>Fig. 19. Comparison of the functions A(s,t) and A₀(s,t); in (a), (b)) the function A_c(s,t) is shown in green. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

It is also of interest that the dynamic equations (142) and (146) can be reformulated as a classical wave equation if we introduce two different (but dependent) displacement fields

$$ \tilde{Y}_1(\infty)D_t^2 \hat{u}_1(s,t) - \tilde{Y}_2(\infty)D_s^2 \hat{u}_2(s,t) = f(s,t). \quad (149) $$

The fields {$\hat{u}_p(s,t)$}$_{p=1,2}$ are defined by the following relations:

$$ \hat{u}_p(s,t) = \int_{t_0}^{t} \tilde{Y}_p^{-1}(\infty) D_t \hat{Y}_p(t-\hat{t}) \hat{u}_0(s, \hat{t}) d\hat{t}, \quad (150) $$

which implies a nonlocal relation between them.

In addition to symmetric formulations of the TN model where time nonlocality is present in both inertial and elastic terms, it is instructive to consider formulations where temporal nonlocality is concentrated either only in the elastic or only in the inertial term. For instance, in terms of the field $\hat{u}_1$ introduced above, the TN dynamic problem can be rewritten as a problem with conventional inertia and time nonlocal elasticity

$$ f(s,t) = \rho D_t^2 \hat{u}_1(s,t) - \frac{\alpha}{\omega_*^2} D_t \int_{t_0}^{t} \Theta_1(t-\hat{t}) D_s^2 \hat{u}_1(s, \hat{t}) d\hat{t} = \rho D_t^2 \hat{u}_1(s,t) - \frac{\alpha}{\omega_*^2} D_t^2 D_s^2 \int_{t_0}^{t} \hat{u}_1(s,t-\hat{t}) D_t \Theta_1(\hat{t}) d\hat{t}, \quad (151) $$

where

$$ \Theta_1(t) = -\frac{1}{2i\pi} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} \frac{\omega_*^2 e^{i\omega t}}{\omega c^2 \lambda_0^2(\omega)} d\omega = \mathrm{H}(t) \left[ \frac{\omega_*^2 t^2}{2} + \frac{1}{3} - 16 \int_{\omega_*}^{\infty} \frac{a}{\ell^+} \frac{\cos(\omega_r t)}{\omega_r \left[ \frac{a^2}{\ell^+} + \pi^2 \right]^2} d\omega_r \right]. \quad (152) $$
<figure><img src="image_21.png" /><figcaption>Fig. 20. The transient behaviors of the functions $\Theta_1^a - \Theta_1$ where $\Theta_1^a(t) = (\omega_*t)^2/2 + 1/3$ (a) and $1 - \Theta_2(t)$ (b).</figcaption></figure>

The behavior of the kernel $\Theta_1(t)$ is shown in Fig. 20(a); one can show that $\Theta_1(0) = 0$, and that $\Theta_1(t) \sim \Theta_1^a(t) = (\omega_*t)^2/2 + 1/3$ at large $\omega_*t$.

After the field $\hat{u}_1$ is known the original displacement field $u$ can be reconstructed explicitly

$$u(s,t) = \frac{1}{\alpha} \int_{t_0}^{t} D_t \gamma_2(t-\hat{t}) \hat{u}_1(s, \hat{t}) d\hat{t}. \qquad (153)$$

Next, suppose that $\hat{u}_1$ solves the trivial initial value problem for (151) which means that

$$\hat{u}_1(s,t) = \int_{t_0}^{t} \int_{-\infty}^{+\infty} D_t A_1(s-\hat{s}, t-\hat{t}) \frac{f(\hat{s}, \hat{t})}{\rho a} d\hat{s} d\hat{t} \qquad (154)$$

with

$$A_1(s,t) = \int_0^t \frac{1}{2\pi} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} \frac{i\lambda_0 a e^{i(\lambda_0 |s| + \omega\hat{t})}}{2\omega^2} d\omega d\hat{t} = \frac{aH(t)}{\pi} \left\{ \begin{aligned} & \int_0^{\pi/a} \lambda_r a \frac{\sin(\omega_*t \sin(\lambda_r a/2)) \sin(\lambda_r |s|)}{2\omega_*^2 \sin^2(\lambda_r a/2) \tan(\lambda_r a/2)} d\lambda_r \\ & - \pi \sin(|s|/\pi/a) \int_0^{+\infty} \frac{e^{-\lambda_i s} \sin(\omega_*t \cosh(\lambda_i a/2))}{2\omega_*^2 \cosh^2(\lambda_i a/2)} \tanh(\lambda_i a/2) d\lambda_i \\ & - \cos(s\pi/a) \int_0^{+\infty} \frac{e^{-\lambda_i s} \sin(\omega_*t \cosh(\lambda_i a/2))}{2\omega_*^2 \cosh^2(\lambda_i a/2)} \lambda_i a \tanh(\lambda_i a/2) d\lambda_i \end{aligned} \right\}. \qquad (155a)$$

In the case of step loading (140), the solution $\hat{u}_1$ of (151) can be again expressed as in (141) with $A$ replaced by $A_1$. This solution is illustrated in Fig. 21(a) and (c). Notice that already $\hat{u}_1$ captures the main macroscopic phenomenon; however, it still misses lattice scale oscillations ahead of the deformation fronts. The recovery of the microscopic solution is performed by using the nonlocal operator (153) which restores the missing lattice scale details (compare Fig. 19(c) with Fig. 21(c)). Observe also that the field $\hat{u}_1$ is more oscillatory than the field $\hat{u}_0$.

Similarly, the TN dynamic problem can be reformulated as a problem with nontrivial inertia and standard CC elasticity in terms of the field $\hat{u}_2$

$$f(s,t) = \rho D_t \int_{t_0}^{t} \Theta_2(t-\hat{t}) D_t^2 \hat{u}_2(s, \hat{t}) d\hat{t} - \alpha D_s^2 \hat{u}_2(s, \hat{t}) = \rho D_t^2 \int_{t_0}^{t} \hat{u}_2(s, t-\hat{t}) D_t \Theta_2(\hat{t}) d\hat{t} - \alpha D_s^2 \hat{u}_2(s, \hat{t}). \qquad (156)$$

Here we introduced another kernel

$$\Theta_2(t) = \frac{c^2}{2\pi i} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} \frac{\lambda_0^2(\omega) e^{i\omega t}}{\omega^3} d\omega = H(t) \left[ 1 - \omega_*^2 \int_{\omega_*}^{\infty} \frac{a}{\ell_+} \frac{\cos(\omega_r t)}{\omega_r^3} d\omega_r \right]. \qquad (157)$$

This function is illustrated in Fig. 20(b); one can see that $\Theta_2(0) = D_t \Theta_2(0) = 0$, and $\Theta_2(\infty) \sim 1$ in the long-time limit.

Once again the original field $u(s,t)$ can be obtained explicitly as

$$u(s,t) = \frac{1}{\rho} \int_{t_0}^{t} D_t \gamma_1(t-\hat{t}) \hat{u}_2(s, \hat{t}) d\hat{t}, \qquad (158)$$
<figure><img src="image_22.png" /><figcaption>Fig. 21. Comparison at fixed instants of $A_1$ (in black) in (a), and $A_2$ (in black) in (b), versus $A_c$ (in green). The differences $A_c(s,t)-A_1(s,t)$ in (c) and $A_c(s,t)-A_2(s,t)$ in (d) at fixed instants. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

where the field $\hat{u}_2(s,t)$ solves the dynamic problem (156)

$$ \hat{u}_2(s,t) = \int_{t_0}^{t} \int_{-\infty}^{+\infty} D_t A_2(s-\hat{s}, t-\hat{t}) \frac{f(\hat{s}, \hat{t})}{\rho a} d\hat{s} d\hat{t}. \qquad (159) $$

Here

$$ A_2(s,t) = \frac{\rho a^2}{4\pi\alpha} \int_{-i\omega_b - \infty}^{-i\omega_b + \infty} \frac{(e^{i\omega t} - 1)e^{i\lambda_0 |s|}}{\omega\lambda_0 a} d\omega = \frac{aH(t)}{\pi} \left\{ \begin{aligned} & \pi t - \int_0^{\pi/a} \frac{2 \sin(\omega_* t \sin(\lambda_r a/2)) \sin(\lambda_r |s|)}{\lambda_r a \omega_*^2 \tan(\lambda_r a/2)} d\lambda_r \\ & - \pi \sin(|s|\pi/a) \int_0^{+\infty} \frac{2e^{-\lambda_i s} \sin(\omega_* t \cosh(\lambda_i a/2))}{\omega_*^2[\lambda_i^2 a^2 + \pi^2]} \tanh(\lambda_i a/2) d\lambda_i \\ & + \cos(s\pi/a) \int_0^{+\infty} \frac{2e^{-\lambda_i s} \sin(\omega_* t \cosh(\lambda_i a/2))}{\omega_*^2[\lambda_i^2 a^2 + \pi^2]} \lambda_i a \tanh(\lambda_i a/2) d\lambda_i \end{aligned} \right\}. \qquad (160a) $$

and in the case of step loading (140), the solution $\hat{u}_2$ of (156) can be expressed as in (141), by replacing $A$ by $A_2$ (see Fig. 21(b) and (d)). Again, the solution $\hat{u}_2$ captures the macroscopic features of dynamics but misses the detailed structure of lattice oscillations. Those can be recovered by applying the operator (158) (compare Fig. 19(c) with Fig. 21(d)). Although the fields $\hat{u}_2$ and $\hat{u}_1$ are qualitatively similar in the sense that they are both more oscillatory than the field $\hat{u}_0(s,t)$, they remain quantitatively different as one can see from Fig. 22. The model based on $\hat{u}_2$ with all temporal nonlocality contained in the inertial term (see (156)) is particularly promising because it can be easily generalized to the case of nonlinear elasticity.
<figure><img src="image_23.png" /><figcaption>Fig. 22. The difference between the kernels $A_1$ and $A_2$. (a) Large scale structure. (b) Small scale structure at different instants.</figcaption></figure>

We now briefly describe one approach to constructing a sequence of approximate models capturing successively shorter time scales. Starting with the classical form of dynamic equations (146) we observe that

$$ \sigma(s,t) = \tilde{Y}_2(\infty) D_s \hat{u}_2(s,t) = \int_{t_0}^{t} D_t Y_2(t-\hat{t}) \left[ 1 + \frac{1}{\omega_*^2} D_t^2 \right] D_s u(s, \hat{t}) \, d\hat{t}, \quad (161a) $$

$$ \mathbf{p}(s,t) = \tilde{Y}_1(\infty) D_t \hat{u}_1(s,t) = D_t \int_{t_0}^{t} D_t Y_1(t-\hat{t}) \left[ 1 + \frac{1}{\omega_*^2} D_t^2 \right] u(s, \hat{t}) \, d\hat{t}. \quad (161b) $$

In the spectral space the relations (161a) and (161b) take the form

$$ \sigma(s, \omega) = \alpha \frac{\sin(\lambda_0 a)}{\lambda_0 a} \varepsilon(s, \omega), \quad (162a) $$

$$ \mathbf{p}(s, \omega) = \rho \omega_*^2 \lambda_0 a \frac{\sin(\lambda_0 a)}{4\omega_*^2} v(s, \omega), \quad (162b) $$

where we introduced the classical velocity field $v = D_t u$ and the classical strain field $\varepsilon = D_s u$. This spectral representation of the nonlocal constitutive relations suggests several ways of constructing approximate QC models. For instance, by using a combination of Taylor and Padé expansions we can generate a series of local constitutive models of the form

$$ \left[ 1 + \sum_{n=1}^{N_\alpha} \frac{(-1)^n}{\omega_*^{2n}} Q_n D_t^{2n} \right] \sigma = \alpha \left[ 1 + \sum_{n=1}^{N_\alpha} \frac{(-1)^n}{\omega_*^{2n}} R_n D_t^{2n} \right] \varepsilon, \quad (163a) $$

$$ \left[ 1 + \sum_{n=1}^{N_\nu} \frac{(-1)^n}{\omega_*^{2n}} S_n D_t^{2n} \right] \mathbf{p} = \rho \left[ 1 + \sum_{n=1}^{N_\nu} \frac{(-1)^n}{\omega_*^{2n}} T_n D_t^{2n} \right] \mathbf{v}, \quad (163b) $$

where the coefficients $(Q_n, R_n, S_n, T_n)$ can be found explicitly. As a basis of such an expansion one can use one of the following representations:

$$ \varepsilon(s, \omega) = \frac{1}{\alpha c_{\text{gr}}(\omega) c_{\text{ph}}(\omega)} c^2 \sigma(s, \omega) \quad \text{and} \quad \mathbf{p}(s, \omega) = \frac{c_{\text{gr}}(\omega)}{c_{\text{ph}}(\omega)} \rho \mathbf{v}(s, \omega), \quad (164a) $$

$$ \frac{c}{c_{\text{ph}}(\omega)} \sigma(s, \omega) = \frac{c_{\text{gr}}(\omega)}{c} \alpha \varepsilon(s, \omega) \quad \text{and} \quad \frac{c}{c_{\text{gr}}(\omega)} \mathbf{p}(s, \omega) = \frac{c}{c_{\text{ph}}(\omega)} \rho \mathbf{v}(s, \omega). \quad (164b) $$

Here we introduced the group and phase velocities

$$ c_{\text{gr}}(\omega) = - \frac{d\omega}{d\lambda_0}(\lambda_0(\omega)) = \frac{\omega_* a}{2} \cos(\lambda_0 a / 2), \quad (165a) $$

$$ c_{\text{ph}}(\omega) = - \frac{\omega}{\lambda_0} = \frac{\omega_* \sin(\lambda_0 a / 2)}{\lambda_0}. \quad (165b) $$
In the low frequency regime $|\omega| < \omega_*$ one can use the Taylor expansions

$$ \frac{c}{c_{ph}(\omega)} = \frac{1}{\Gamma(1/2)} \sum_{n \in \mathbb{N}} \frac{\Gamma(n+1/2)}{\Gamma(n+1)(2n+1)} \left(\frac{\omega}{\omega_*}\right)^{2n}, $$

$$ \left[ \frac{c}{c_{gr}(\omega)} \right]^p = \frac{1}{\Gamma(p/2)} \sum_{n \in \mathbb{N}} \frac{\Gamma(n+p/2)}{\Gamma(n+1)} \left( \frac{\omega}{\omega_*} \right)^{2n}, \quad p = \pm 1, $$

where $\Gamma$ is Euler's function (Abramowitz and Stegun, 1972). If we identify the coefficients in (163) according to (164a) or (164b), we obtain

$$ Q_n = \frac{1}{[\Gamma(1/2)]^2} \sum_{k=0}^{n} \frac{1}{2k+1} \frac{\Gamma(k+1/2) \Gamma(n-k+1/2)}{\Gamma(k+1) \Gamma(n-k+1)}, $$

$$ T_n = \frac{1}{\Gamma(1/2)\Gamma(-1/2)} \sum_{k=0}^{n} \frac{1}{2k+1} \frac{\Gamma(k+1/2) \Gamma(n-k-1/2)}{\Gamma(k+1) \Gamma(n-k+1)}, $$

$$ R_n = S_n = 0. $$

for the case (164a) and

$$ Q_n = T_n = \frac{S_n}{2n+1} = \frac{1}{2n+1} \frac{\Gamma(n+1/2)}{\Gamma(n+1)\Gamma(1/2)}, \quad R_n = \frac{1}{\Gamma(-1/2)} \frac{\Gamma(n-1/2)}{\Gamma(n+1)}. $$

for the case (164b). In order to decide which approximation is appropriate one needs to specify the nature of the physical problem. Several interesting attempts in this direction can be already found in the literature (e.g. Mindlin, 1964; Polyzos and Fotiadis, 2012).

## 7. Discussion

In this paper we made an attempt to extend the framework of continuum elasto-dynamics into the domain where the behavior of a crystal lattice is strongly discrete. This means that we consider problems where both spatial and temporal external scales, describing either geometry or loading, are comparable to their internal counterparts dictated by the lattice itself. In order to focus on nonlocality and dispersion we excluded the effects of nonlinearity. In this way we avoided the localization, on one side, and the energy tunneling from long to short waves with subsequent thermalization, on the other side. We also minimized the geometric and structural complexity of the problem by focussing on the simplest 1D chain with NN interactions.

The main goal of our study was to construct a QC model whose solutions provide an *exact interpolation* of the solution of the fully inertial discrete model. Since the interpolation problem has typically multiple solutions one is confronted with the choice that can be guided by physical, analytical or computational considerations.

The previous interpolation related efforts, mainly associated with the names of Kunin (1982) and Rogula (1982), were aimed at the models with classical inertia and strong spatial nonlocality. Such models have been recently rediscovered by the computational community under the name of peridynamic theory (e.g Silling, 2000) and are widely used as a coarse grained form of molecular dynamics. However, the analytical complexity of this approach and the necessity to deal with long range spatial interactions extending beyond the actual boundaries of the bodies, makes the detailed mathematical analysis elusive. Therefore our main goal was to develop a spatially local QC theory which is nevertheless fully faithful to the original discrete model and still allows for temporally local approximations.

In our previous work (Charlotte and Truskinovsky, 2008) devoted to statics we have shown that a *spatially nonlocal QC model* of a simple lattice with arbitrary long interactions can always be replaced by a series of *spatially local models* incorporating auxiliary variables describing micro-corrections to the macroscopic displacement field. In this paper we extended this framework to dynamics. Assuming that the static part of the model is understood we focussed on the simplest inertial effects. To this end we have drastically simplified the spatial structure of the lattice model by considering a simple particle chain without nontrivial internal length scales.

Our main result is that lattice scale elasto-dynamics can be adequately represented at the continuum level by a temporal memory structure. More precisely, we proposed a QC model with nontrivial inertia and hereditary elasticity (temporally nonlocal or TN model) which accomplishes an *exact interpolation* of our prototypical discrete model. The unusual non-Newtonian inertia can be understood in the framework of a meta-material paradigm (Milton and Willis, 2007) which presumes that a “visible” continuum particle carries a variety of internal degrees of freedom representing locally non-affine (non-Cauchy-Born) dynamic responses. The memory structure of our elasticity can be linked to the “structural attenuation” (Brillouin and Parodi, 1956) involving energy redistributions between microscopic degrees of freedom. In this sense both the nonclassical inertia and the hereditary/memory structure of elasticity are the means of adequate continuum representation of high frequency vibrations of the lattice without introducing temperature (see also Park and Liu, 2004; Park et al., 2005; Tang et al., 2006; Karpov et al., 2007; Tadmor and Miller, 2011). The extension of the
proposed model towards including interactions beyond nearest neighbors is straightforward in view of the results obtained in Charlotte and Truskinovsky (2008) and will require adding appropriate internal variables equipped with their own temporal nonlocality.

Despite the unquestionable complexity of the derived TN model in the time domain, it remains *local in space* which facilitates spatial matching of the QC based computations with the MD based computations. In particular, the proposed QC scheme can be used as an interface between discrete and CC representations of lattice vibrations (Farrell et al., 2007). The development of such an interface is crucial for building effective hybrid discrete-continuous codes aimed at solving a variety of multi-scale problems involving nucleation and interaction of crystal defects, as well as addressing other challenges raised by nanotechnologies (Fish, 2006).

## Acknowledgments

L.T. was supported by the EU contract MRTN-CT-2004-505226. The authors would like to thank the anonymous reviewers for their valuable suggestions.

## Appendix A. Spatial regularization in the impact test

Consider the following delta sequence of loading functions which are regular both in space and time:

$$f^{\beta, \epsilon}(s,t) = \sum_{k \in \mathbb{Z}} \frac{f_k^{\epsilon}(t)}{\beta} \text{sinc}((s-ka)/\beta a). \qquad (A.1)$$

Here $a f_k^{\epsilon}(t)$ is given by (53). The sequence (A.1) preserves the resultant force $\int_{-\infty}^{\infty} f^{\beta, \epsilon}(s,t) ds = a \sum_{k \in \mathbb{Z}} f_k^{\epsilon}(t)$ and in the limit $\beta \to 0^+$ and $\epsilon \omega_* \to 0^+$ we recover the localized impact loading (21). The solution of the continuum problem (42) with the source term given by (A.1) can be written as

$$u_c^{\beta, \epsilon}(s,t) = \sum_{p \in \mathbb{Z}} \int_0^t G_c^*(s-pa)/\beta(t-\hat{t})/\beta) \frac{f_p^{\epsilon}(\hat{t})}{\rho} d\hat{t} = \sum_{k \in \mathbb{Z}} u_c^{\beta, \epsilon}(k\beta a, t) \text{sinc}(s/\beta a - k). \qquad (A.2)$$

According to (61), the function $G_c^*(s/\beta, t/\beta)$ converges to $G_c(s,t)$ in (47) when $\beta \to 0$. If we now compute the work (41c) for the regularized impact solution, we obtain

$$\begin{aligned} \mathcal{P}_r(f^{\beta, \epsilon}, u_c^{\beta, \epsilon}) &= \sum_{k \in \mathbb{Z}} \int_0^t \int_{-\infty}^{\infty} \frac{f^{\beta, \epsilon}(s, \hat{t})}{\beta} \hat{u}_c^{\beta, \epsilon}(s, \hat{t}) ds d\hat{t} \\ &= \frac{\bar{p}}{\epsilon} \{[H(t-t_0)-H(t-t_0-\epsilon)][u_c^{\beta, \epsilon}(0,t)-u_c^{\beta, \epsilon}(0,t_0)]+H(t-t_0-\epsilon)[u_c^{\beta, \epsilon}(0,t_0+\epsilon)-u_c^{\beta, \epsilon}(0,t_0)]\}. \end{aligned}$$

Suppose that $\beta > 0$ is fixed and consider the limit $\epsilon \omega_* \to 0^+$. Then the loading and the resulting displacement field remain spatially regular as we see from

$$f^{\beta}(s,t) = \lim_{\omega_* \epsilon \to 0} f^{\beta, \epsilon}(s,t) = \sum_{k \in \mathbb{Z}} \frac{f_k(t)}{\beta} \text{sinc}((s-ka)/\beta a) \qquad (A.3)$$

and

$$u_c^{\beta}(s,t) \stackrel{\text{def}}{=} \lim_{\omega_* \epsilon \to 0} u_c^{\beta, \epsilon}(s,t) = G_c^*(s/\beta, (t-t_0)/\beta) \frac{\bar{p}}{\rho a}. \qquad (A.4)$$

If now $\beta \to 0^+$ we recover the CC solution for the impact loading

$$u_c(s,t) = G_c(s,t-t_0) \frac{\bar{p}}{\rho a}.$$

The computation of the work in the limit $\epsilon \to 0^+$ gives

$$\mathcal{P}_r^c(f, u_c^\beta, t) = H(t-t_0) \bar{p} D_t u_c^\beta(0, t_0^+) = H(t-t_0) \frac{\bar{p}^2}{\beta \rho a} = \frac{1}{\beta} \mathcal{P}_r(f, u, t). \qquad (A.5)$$

One can see that in the limit $\beta \to 0^+$ this expression diverges. We also notice that the value $\beta = 1$ (WSK regularized load) ensures a perfect match of the discrete work in the CC theory.

## Appendix B. An alternative representation of the spatially nonlocal QC model

Given that (76) can be rewritten as

$$[\Phi(\lambda, 0) - \rho\omega^2]u^+(\lambda, \omega) = f^\dagger(\lambda, \omega) \quad \text{with} \quad \Phi(\lambda, 0) = \frac{4\alpha}{a^2} \sin^2(\lambda a/2) = \frac{\alpha\lambda^2}{M^2(\lambda)}$$
and that (91) provides $u^{\dagger}(\lambda, \omega) = M^{\dagger}(\lambda)\tilde{U}^{\dagger}(\lambda, \omega)$, it is clear that the new variable $\tilde{U}$ must satisfy the following spectral equation:

$$[\alpha\lambda^2 - \rho\omega^2 R^{\dagger}(\lambda)]\tilde{U}^{\dagger}(\lambda, \omega) = M^{\dagger}(\lambda)f^{\dagger}(\lambda, \omega) \quad \text{where } R^{\dagger}(\lambda) = M^{\dagger 2}(\lambda) = \frac{\alpha\lambda^2}{\Phi(\lambda, 0)}. \qquad (B.1)$$

In the physical space this leads to the following dynamic equation:

$$\int_{-\infty}^{\infty} \rho R(s-\tilde{s}) D_t^2 \tilde{U}(\tilde{s}, t) d\tilde{s} - D_s^2 \tilde{U}(s,t) = \int_{-\infty}^{\infty} M(s-\tilde{s}) f(\tilde{s}, t) d\tilde{s}. \qquad (B.2)$$

The next step is to compute explicitly the kernels $M(s)$ and $R(s)$ by using the inverse **IF** transform (43a) in the sense of generalized tempered functions (Schwartz, 1966, 1983; Bremermann, 1965) and the following Mittag-Leffler fractional decompositions:

$$M^{\dagger}(\lambda) = 1 + \lambda^2 \sum_{q \ge 1} \frac{(-1)^q}{\lambda_{2q}(0)} \left[ \frac{1}{\lambda - \lambda_{2q}(0)} - \frac{1}{\lambda + \lambda_{2q}(0)} \right], \qquad (B.3)$$

$$R^{\dagger}(\lambda) = 1 + \lambda^2 \sum_{q \ge 1} \frac{2}{\lambda_{2q}^2(0)} + 2\lambda^3 \sum_{q \ge 1} \frac{1}{\lambda_{2q}(0)} \left( \frac{1}{[\lambda - \lambda_{2q}(0)]^2} - \frac{1}{[\lambda + \lambda_{2q}(0)]^2} \right) - \lambda^4 \sum_{q \ge 1} \frac{1}{\lambda_{2q}^2(0)} \left( \frac{1}{[\lambda - \lambda_{2q}(0)]^2} + \frac{1}{[\lambda + \lambda_{2q}(0)]^2} \right), \qquad (B.4)$$

where according to (12) we must put $\lambda_{2q}(0) = 2q\pi a^{-1}$. The explicit representations for $M(s)$ and $R(s)$ can be written either as

$$M(s) = D_s^2 \left\{ \frac{|s|}{2} + |s| \sum_{q \ge 1} (-1)^q \text{sinc}(2qs/a) \right\}, \qquad (B.5)$$

$$R(s) = D_s^2 \left\{ \frac{|s|}{2} + 2D_s \sum_{q \ge 1} |s| \text{sinc}(2qs/a) - \frac{1}{2} D_s^2 \sum_{q \ge 1} |s|^3 \text{sinc}^2(qs/a) \right\} \qquad (B.6)$$

or as (see also Higgins, 1996; Marks, 1991)

$$M(s) = D_s^2 \left[ \frac{|s|}{2} + \sum_{p \in \mathbb{Z}} \frac{|s| [pa-s]}{2s} H\left(1 - \frac{|s|}{a} - p\right) \right], \qquad (B.7)$$

$$R(s) = D_s^2 \left\{ \frac{|s|}{2} + D_s \sum_{p \in \mathbb{Z}} \frac{|s|[(2p+1)a-2s]}{2} H\left(1 - \frac{|2s|}{a} - 2p - 1\right) - \frac{1}{4} D_s^2 \sum_{p \in \mathbb{Z}} |s|[s-pa|(p+1)a-s|H\left(1 - \frac{|2s|}{a} - 2p - 1\right)| \right\}. \qquad (B.8)$$

For comparison, we mention that in the classical (long-wave) continuum approximation with $u \equiv \tilde{U}$ we preserve only the first terms in these expansions

$$M(s) = R(s) = D_s^2 \frac{|s|}{2} = \delta(s) \Leftrightarrow R(\lambda) = 1 = M(\lambda).$$

Now, by using the variable $\tilde{U}$ we can rewrite the elastic energy, the kinetic energy and the external work as

$$\tilde{\varepsilon}_e^*(\tilde{U}, t) = \frac{\alpha}{2} \int_{-\infty}^{\infty} [D_s \tilde{U}(s, t)]^2 ds = \varepsilon_e^c(\tilde{U}, t), \qquad (B.9a)$$

$$\tilde{\varepsilon}_k^*(\tilde{U}, t) = \frac{\rho}{2} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} R(s-\tilde{s}, t) D_t \tilde{U}(s, t) D_t \tilde{U}(\tilde{s}, t) d\tilde{s} ds, \qquad (B.9b)$$

$$\tilde{P}_r^*(f, \tilde{U}, t) = \int_0^t \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} M(s-\tilde{s}) f(\tilde{s}, \tilde{t}) D_t \tilde{U}(s, \tilde{t}) d\tilde{s} ds d\tilde{t} \equiv P_r^c(f, u, t). \qquad (B.9c)$$

The comparison with energies introduced in (90) shows that

$$\tilde{P}_r^*(f, \tilde{U}, t) = P_r^*(f, u, t), \quad \tilde{\varepsilon}_e^*(\tilde{U}, t) = \varepsilon_e^*(u, t) \quad \text{and} \quad \tilde{\varepsilon}_k^*(\tilde{U}, t) = \varepsilon_k^*(u, t).$$

The solution of the spectral problem (B.1) with the spatially smooth loading force (82) can be written as

$$\tilde{U}(s,t) = \sum_{p \in \mathbb{Z}} \int_{0}^{t} D_{t}A_{K}(s-pa,t-\hat{t}) \frac{f_{p}(\hat{t})}{\rho} d\hat{t} \qquad (B.10)$$
<figure><img src="image_24.png" /><figcaption>Fig. 23. Comparison of the kernels $A_K$ (in black) and $A_c$ (in green) (a). The difference $A_c(s,t)-A_K(s,t)$ (b). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

with

$$A_K(s,t) = \frac{4H(t)}{\pi\omega_*^2} \int_0^{\pi/a} \frac{\sin^2(\omega_*t \sin(\lambda_r a/2)/2) \cos(\lambda_r s)}{\lambda_r \sin(\lambda_r a/2)} d\lambda_r. \qquad (B.11)$$

In particular, for a spatially smooth, time-step, loading combining (25) with (82) this solution reads

$$\tilde{U}(s,t) = \frac{\omega_* \bar{p}}{\rho} [A_K(s,t) - A_K(s,t_0)]. \qquad (B.12)$$

Fig. 23 illustrates this solution and compares it with the corresponding solution of the CC theory. Interestingly, the new representation captures lattice oscillations ahead of the propagating front.

Despite the obvious appeal of having a local elastic energy in the new variables the complexity of the kinetic energy kernel makes it hardly practical outside some particularly simple dynamic situations (e.g. traveling waves as in Truskinovsky and Vainchtein, 2006). In contrast, the alternative formulation in term of the original field $u(s,t)$ is quite easy to implement numerically and different nonlinear extensions of the model (87) have recently gained a lot of attention (e.g. Silling, 2000; Seleson et al., 2009).

## Appendix C. Artifacts associated with singular forces in spatially nonlocal QC model

In order to show that the general approach introduced in Section 4 provides at most an approximation in the case of localized loadings, it is sufficient to consider a static problem.

Suppose that the self-equilibrated force distribution $\tilde{f}(s)$ is such that its IF-spectral support $\tilde{f}^{\dagger}(\lambda)$ on $a^{-1}\mathbb{R}$ is not limited to $\mathbb{K}$. If we use the approximation (86), then the equilibrium displacement field must solve the following set of difference equations:

$$\frac{\alpha}{a^2} [\tilde{u}(s+a) + \tilde{u}(s-a) - 2\tilde{u}(s)] = -\tilde{f}(s). \qquad (C.1)$$

The solution $\tilde{u}(s)$ can be written in the form

$$\tilde{u}(s) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \frac{\tilde{f}^{\dagger}(\lambda)}{\Phi(\lambda, 0)} e^{is\lambda} d\lambda + \tilde{u}^{h}(s), \qquad (C.2)$$

where according to (4) and (8)

$$\Phi(\lambda, 0) = \rho \omega_*^2 \sin^2(\lambda a/2) = \frac{4\alpha}{a^2} \sin^2(\lambda a/2).$$

Here $\tilde{u}^h(s) = \sum_{q \in \mathbb{Z}} [A_q + B_q s] e^{is\lambda_{2q}(0)}$ is the general solution of the homogeneous equation and $\lambda_{2q}(0) = 2q\pi/a = -\lambda_{2q+1}(0)$. The unknown coefficients $(A_q, B_q)$ and $(A_{-q}, B_{-q})$ are complex conjugate; if the spectrum of $\tilde{u}^h(s)$ is inside $\mathbb{K} \subset \mathbb{B}_0$, only the two real coefficients $(A_0, B_0)$ are different from zero.
By using the Mittag-Leffler fractional decomposition (B.4) for $[\Phi(\lambda,0)]^{-1} = \mathcal{R}(\lambda)/\alpha\lambda^2$, we can rewrite (C.2) as

$$ \tilde{u}(s) - \tilde{u}^h(s) = \sum_{p=0}^{2} a^{p+1} \int_{-\infty}^{\infty} S_p \left(\frac{s-s^*}{a}\right) D_s^p \tilde{f}(s) ds^*, \quad (C.3) $$

where the functions $S_p(\eta)$ for $\eta \in \mathbb{R}$ are defined as follows:

$$ S_0(\eta) = -\frac{|\eta|}{2\alpha}, \quad (C.4a) $$

$$ S_1(\eta) = -\frac{2|\eta|\eta}{\alpha} \sum_{q \ge 1} \text{sinc}(2q\eta) \stackrel{\text{a.e.}}{=} -\frac{|\eta|}{2\alpha} \sum_{p \in \mathbb{Z}} [2p+1-2\eta]\text{H}(1-|2\eta-2p-1|), \quad (C.4b) $$

$$ S_2(\eta) = \frac{|\eta|^3}{2\alpha} \sum_{q \ge 1} \text{sinc}^2(q\eta) = \frac{|\eta|}{4\alpha} \sum_{p \in \mathbb{Z}} [\eta - p][p + 1 - \eta]\text{H}(1 - |2\eta - 2p - 1|). \quad (C.4c) $$

We can now compare $\tilde{u}(s)$ in (C.2) with the solution of the discrete static problem $\tilde{u} = \{\tilde{u}_k\}_{k \in \mathbb{Z}}$, which satisfies the equation

$$ \frac{\alpha}{a^2} [\tilde{u}_{k+1} + \tilde{u}_{k-1} - 2\tilde{u}_k] = -\tilde{f}_k, \quad (C.5) $$

and can be written as (Charlotte and Truskinovsky, 2008)

$$ \tilde{u}_k = a \sum_{p \in \mathbb{Z}} \tilde{G}_c(ka - pa) \tilde{f}_p + \tilde{u}_c^h(ka), \quad (C.6a) $$

where

$$ \tilde{G}_c(s) = a S_0(s/a) = -\frac{|s|}{2\alpha}. \quad (C.6b) $$

Here $\tilde{u}_c^h(s) \stackrel{\text{def}}{=} A_c + B_c s$, with real constants ($A_c, B_c$), is the general solution of the homogeneous equation (C.5).

One can see that for the special solutions of the inhomogeneous problem, $\tilde{u}(s)$ in (C.3) and $u$ in (C.6a), to agree at the lattice points $a\mathbb{Z}$, we need to place special restrictions on $\tilde{f}(s)$ and $\tilde{f} = (\tilde{f}_k)_{k \in \mathbb{Z}}$, namely that the spectral support of $\tilde{f}(s)$ on $a\mathbb{R}$ is limited to $\mathbb{K}$ and that it represents a maximally smooth interpolation of the discrete system $\tilde{f}$

$$ \tilde{f}(s) = \sum_{k \in \mathbb{Z}} \tilde{f}_k \text{sinc}(s/a-k). \quad (C.7) $$

In this case (C.3) can be written as

$$ \tilde{u}(s) - \tilde{u}^h(s) = a \sum_{p \in \mathbb{Z}} \tilde{G}_K(s - pa) \tilde{f}_p, \quad (C.8a) $$

where

$$ \tilde{G}_K(s) \stackrel{\text{def}}{=} \frac{1}{2\pi} \int_{-\pi/a}^{\pi/a} \frac{e^{is\lambda}}{\Phi(\lambda, 0)} d\lambda \equiv \tilde{G}_c(s) + \frac{\sin(\pi|s|/a)}{\pi} \int_0^{\infty} \frac{a^2 e^{-|s|^\lambda}}{4\alpha \cosh^2(\lambda a/2)} d\lambda. \quad (C.8b) $$

Now it is not difficult to check that $\tilde{u}(ka) - \tilde{u}^h(ka) = \tilde{u}_k - \tilde{u}_c^h(ka)$ for $k \in \mathbb{Z}$. In the opposite case when the loads are singular $\tilde{f}(s) = \sum_{k \in \mathbb{Z}} \tilde{f}_k \delta(s/a-k)$ the function (C.3) obviously disagrees with the discrete solution.

## References

* Abramowitz, M., Stegun, I.A. (Eds.), 1972. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover, New York.

* Andrianov, I.V., Awrejcewicz, J., Weichert, D., 2010. Improved continuous models for discrete media. Math. Probl. Eng., 35.

* Arndt, M., Griebel, M., 2005. Derivation of higher order gradient continuum models from atomistic models for crystalline solids. Multiscale Model. Simul. 4, 531–562.

* Askes, H., Metrikine, A., Pichugin, A.V., Bennett, T., 2008. Four simplified gradient elasticity models for the simulation of dispersive wave propagation. Philos. Mag. 88 (28–29), 3415–3443.

* Atkinson, W., Cabrera, N., 1965. Motion of a Frenkel–Kontorova dislocation in a one-dimensional crystal. Phys. Rev. A 138 (3), A763–A766.

* Blanc, X., Le Bris, C., Lions, P.-L., 2002. From molecular models to continuum mechanics. Arch. Ration. Mech. Anal. 164, 341–381.

* Born, M., Huang, K., 1954. Dynamical Theory of Crystal Lattices. Oxford University Press, Oxford.

* Bremermann, H., 1965. Distributions, Complex Variables and Fourier Transforms. Addison-Wesley Publishing Company, Inc.

* Braides, A., Truskinovsky, L., 2008. Asymptotic expansions by $\Gamma$ convergence. Continuum Mech. Thermodyn. 20, 21–62.

* Brillouin, L., Parodi, M., 1956. Propagation des ondes dans les milieux periodiques. Masson, Dunod, Paris.

* Carpio, A., Bonilla, L.L., 2005. Discrete models of dislocations and their motion in cubic crystals. Phys. Rev. B 71, 134105.

* Charlotte, M., Truskinovsky, L., 2002. Linear chain with a hyper-pre-stress. J. Mech. Phys. Solids 50, 217–251.

* Charlotte, M., Truskinovsky, L., 2008. Towards multi-scale continuum elasticity theory. Continuum Mech. Thermodyn. 20, 133–161.

* Collins, M.A., 1981. A quasi-continuum approximation for solitons in an atomic chain. Chem. Phys. Lett. 77, 342–347.

* Curtin, W.A., Miller, R.E., 2003. Atomistic/continuum coupling in computational materials science. Model. Simul. Mater. Sci. Eng. 11, R33–R68.

* Deift, P., McLaughlin, T.-R., 1998. A Continuum Limit of the Toda Lattice. Memoirs of the American Mathematical Society, vol. 131/624. American Mathematical Society.

* DiVincenzo, D.P., 1986. Dispersive corrections to continuum elastic theory in cubic crystals. Phys. Rev. B 34, 5450–5465.
Dobson, M., Elliot, R.S., Luskin, M., Tadmor, E.B., 2007. A multilattice quasicontinuum for phase transforming materials: cascading Cauchy-Born kinematics. J. Comput. Aided Mater. Des. 14 (suppl. 1), 219-237.

Efendiev, Y.R., Truskinovsky, L., 2010. Thermalization of a driven bi-stable FPU chain. Continuum Mech. Thermodyn. 22, 679-698.

Eringen, A.C., 1972. Linear theory of nonlocal elasticity and dispersion of plane waves. Int. J. Eng. Sci. 10, 425-435.

Farrell, D.E., Karpov, E.G., Liu, W.K., 2007. Algorithms for bridging scale method parameters. Comput. Mech. 40, 965-978.

Fish, J., 2006. Bridging the scales in nano engineering and science. J. Nanoparticle Res. 8, 577-594.

Friesecke, G., Theil, F., 2002. Validity and failure of the Cauchy-Born rule in a 2D mass-spring lattice. J. Nonlinear Sci. 12, 445-478.

Hays, M.H., Levermore, C.D., Miller, P.D., 1994. Macroscopic lattice dynamics. Physica D 79 (1), 1-15.

Higgins, J.R., 1996. Sampling Theory in Fourier and Signal Analysis. Foundations. Clarendon Press, Oxford.

E, W., Huang, Z., 2002. A dynamic atomistic continuum method for the simulation of crystalline materials. J. Comput. Phys. 182, 234-261.

Il'ushina, E.A., 1969. On a model of continuous medium, taking into account the microstructure. Prikl. Matem. Mekhan. 33 (5), 917-923.

Karpov, E.G., Park, H.S., Liu, W.K., 2007. A phonon heat bath approach for the atomistic and multiscale simulation of solids. Int. J. Numer. Methods Eng. 70, 351-378.

Kevrekidis, P.G., Kevrekidis, I.G., Bishop, A.R., Titi, E.S., 2002. Continuum approach to discreteness. Phys. Rev. E 65, 046613.

Krumshansl, J.A. 1965. Generalized continuum field representation for lattice vibrations. In: Wallis, R.F. (Ed.), Lattice Dynamics. Pergamon, London, pp. 627-634.

Kröner, E. (Ed.), 1968. Mechanics of generalized continua. Proceedings of the IUTAM-Symposium. Springer, Berlin.

Kunin, I., 1982. Elastic Media with Microstructure, v.I (One Dimensional Models). Springer, Berlin.

Landau, L.D., Pitaevskii, L.P., Lifshitz, E.M., Kosevich, A.M., 1986. Theory of Elasticity. Butterworth-Heinemann.

E, W., Li, X., 2005. Multiscale modeling of crystalline solids. In: Handbook of Materials Modeling. Springer, pp. 1491-1506.

Lin, P., 2007. Convergence analysis of a quasi-continuum approximation for a two-dimensional material without defects. SIAM J. Numer. Anal. 45, 313-332.

Liu, W.K., Karpov, E.G., Zhang, S., Park, H.S., 2004. An introduction to computational nanomechanics and materials. Comput. Methods Appl. Mech. Eng. 193, 1529-1578.

E, W., Lu, J., Yang, J., 2006. Uniform accuracy of the quasicontinuum method. Phys. Rev. B 74, 214115.

Maciá, F., 2004. Wigner measures in the discrete setting: high-frequency analysis of sampling and reconstruction operators. SIAM J. Math. Anal. 36 (2), 347-383.

Maradudin, A.A., Montroll, E.W., Weiss, G.H., 1963. Theory of Lattice Dynamics in the Harmonic Approximation. Academic Press, New York.

Marks II, R.J., 1991. Introduction to Shannon Sampling and Interpolation Theory. Springer-Verlag, New York.

Massobrio, C., Bulou, H., Goyhenex, C. (Eds.), 2010. Atomic-Scale Modeling of Nanosystems and Nanostructured Materials. Lecture Notes in Physics, vol. 795. Springer.

Marian, J., Venturini, G., Hansen, B.L., Knap, J., Ortiz, M., Campbell, G.H., 2010. Finite-temperature extension of the quasicontinuum method using Langevin dynamics: entropy losses and analysis of errors. Modell. Simul. Mater. Sci. Eng. 18, 015003.

Mielke, A., 2006. Macroscopic behavior of microscopic oscillations in harmonic lattices via Wigner-Husimi transforms. Arch. Ration. Mech. Anal. 181 (3), 401-448.

Miller, R.E., Tadmor, E.B., 2002. The quasi-continuum method: overview, applications and current directions. J. Comput. Aided Mater. Des. 9, 203-239.

Milton, G.W., Willis, J.R., 2007. On modifications of Newton's second law and linear continuum elastodynamics. Proc. R. Soc. A 463, 855-880.

Mindlin, R., 1964. Micro-structure in Linear Elasticity. Arch. Ration. Mech. Anal. 10, 51-77.

Mindlin, R., 1965. Second gradient of strain and surface tension in linear elasticity. Int. J. Solids Struct. 1, 417-438.

Ortiz, M., Philips, R., Tadmor, E.B., 1996. Quasicontinuum analysis of defects in solids. Philos. Mag. A 73, 1529-1563.

Park, H.S., Liu, W.K., 2004. An introduction and tutorial on multi-scale analysis in solids. Comput. Methods Appl. Mech. Eng. 193, 1733-1772.

Park, H.S., Karpov, E.G., Liu, W.K., Klein, P.A., 2005. The bridging scale for two-dimensional atomistic/continuum coupling. Philos. Mag. 85 (1), 79-113.

Polyzos, D., Fotiadis, D.I., 2012. Derivation of Mindlin's first and second strain gradient elastic theory via simple lattice and continuum models. Int. J. Solids Struct. 49, 470-480.

Roddier, F., 1971. Distributions et transformation de Fourier. Ediscience, Paris, pp. 57-64 (Chapter 4).

Rogula, D., 1982. Introduction to nonlocal theory of material media. In: Nonlocal Theory of Elastic Media. CISM Courses. Springer, Berlin, vol. 268, pp. 123-222.

Rosenau, P., 2003. Hamiltonian dynamics of dense chains and lattices: or how to correct the continuum. Phys. Lett. A 311, 39-52.

Schwartz, L., 1966. Théorie des distributions. Hermann, Paris.

Schwartz, L., 1983. Méthodes mathématiques pour les sciences physiques. Hermann, Paris.

Slepyan, L.I., 1972. Non-stationary Elastic Waves. Sudostroenie, Leningrad (in Russian).

Silling, S.A., 2000. Reformulation of elasticity theory for discontinuities and long-range forces. J. Mech. Phys. Solids 48, 175-209.

Seleson, P., Parks, M.L., Ginzburger, M., Lehoucq, R.B., 2009. Peridynamics as upscaling of molecular dynamics. Multiscale Model. Simul. 8 (1), 204-227.

Slepyan, 1981. Dynamics of a crack in a lattice. Sov. Phys.—Dokl. 26 (5), 538-540.

Slepyan, L.I., 2002. Models and Phenomena in Fracture Mechanics. Springer-Verlag.

Sunyk, R., Steinmann, P., 2003. On higher gradients in continuum-atomistic modelling. Int. J. Solids Struct. 40, 6877-6896.

Tadmor, E.B., Ortiz, M., Philips, R., 1996. Quasicontinuum analysis of defects in solids. Philos. Mag. A 73 (6), 1529-1563.

Tadmor, E.B., Miller, R.E., 2011. Modeling Materials: Continuum, Atomistic and Multiscale Techniques. Cambridge University Press.

Tang, S., Hou, T.Y., Liu, W.K., 2006. A mathematical framework of the bridging scale method. Int. J. Numer. Methods Eng. 65, 1688-1713.

Toupin, R.A., 1962. Elastic materials with couple stresses. Arch. Ration. Mech. Anal. 11, 385-414.

Truskinovsky, L., Vainchtein, A., 2003. Peierls-Nabarro landscape for martensitic phase transitions. Phys. Rev. B 67, 172103.

Truskinovsky, L., Vainchtein, A., 2005. Kinetics of martensitic phase transitions: lattice model. SIAM J. Appl. Math. 66, 533-553.

Truskinovsky, L., Vainchtein, A., 2006. Quasicontinuum models of dynamic phase transitions. Continuum Mech. Thermodyn. 18, 1-21.

Vasiliev, A.A., Dmitriev, S.V., Bishop, A.R., Miroshnichenko, A.E., 2010. Multi-field approach in mechanics of structural solids. Int. J. Solids Struct. 47 (3-4), 510-525.

Vasseur, J.O., Deymier, P.A., Chenni, B., Djafari-Rouhani, B., Dobrzynski, L., Prevost, D., 2001. Experimental and theoretical evidence for the existence of absolute acoustic band gaps in two-dimensional solid phononic crystals. Phys. Rev. Lett. 86 (14), 3012-3015.

Willis, J.R., 1981. Variational and related methods for the overall properties of composites. Adv. Appl. Mech. 21, 1-78.

Willis, J.R., 1997. Dynamics of composites. In: Suquet, P. (Ed.), Continuum Micromechanics: CISM Courses and Lectures, vol. 377, Springer-Verlag-Wien, New York, pp. 265-290.

Wolf, K.B., 1979. Integral Transforms in Science and Engineering. Plenum Press, New York.