Justification of Leading Order Quasicontinuum Approximations of Strongly Nonlinear Lattices

We consider the leading order quasicontinuum limits of a one-dimensional granular medium governed by the Hertz contact law under precompression. The approximate model which is derived in this limit is justified by establishing asymptotic bounds for the error with the help of energy estimates. The continuum model predicts the development of shock waves, which are also studied in the full system with the aid of numerical simulations. We also show that existing results concerning the Nonlinear Schrodinger (NLS) and Korteweg de-Vries (KdV) approximation of FPU models apply directly to a precompressed granular medium in the weakly nonlinear regime.

1. Introduction. We consider a one-dimensional granular medium which is governed by the Hertz law. Denote q n the relative displacement from equilibrium of the n-th particle. Then the renormalized equations of motion describing the q n 's have the form [12], where W ′ (u) = [δ 0 + u] p + , [u] + = u Θ(u) = max(0, u) (2) where Θ is the Heaviside function and δ 0 is the static load (precompression) applied to the chain at all times. There is a double nonlinearity stemming from this system due to the lack of tensile strength (resulting in an asymmetric potential function W ) and the nonlinear coupling. For spherical particles we have p = 3/2, which corresponds to the classical Hertz contact law.
With the long wave ansatz u n (t) = A(X, T ), X = εn, T = εt where X, T, A(X, T ) ∈ R and ε ≪ 1 is a small perturbation parameter, one can derive the leading order continuum model approximation to the discrete dynamics The main result of the present paper is the following approximation theorem. be a solution of (4). Then for C 1 > 0 sufficiently small there exists C, ε 0 > 0 such that for all ε ∈ (0, ε 0 ) there are solutions (u n (t)) n∈Z of (3) satisfying  In the weakly nonlinear case, i.e. if the precompression is much greater than the amplitude of the solution, then the Nonlinear Schrödinger (NLS) and Korteweg de-Vries (KdV) equations can be derived as continuum models for solutions that have amplitude of order O(ε) and O(ε 2 ) respectively, see Section 8. For solutions of order O(1) (i.e. solutions with amplitude that are on the same order as the precompression) model (4) is relevant. In order to establish local existence and uniqueness of the continuum model and to handle the nonlinearity in Fourier space, we will expand the nonlinearity as a series. Since the amplitude is O(1) we cannot truncate this series (which is done when deriving the NLS and KdV equations). Therefore, we require that the precompression is of the same order as the amplitude (but not greater). It is unclear if this is only a technical assumption or if it is really necessary for having the approximation property. It will be the subject of future research to explore whether or not the assumptions of Theorem 1.1 can be relaxed, and if so, how the continuum model (4) needs to be modified.
If one does not ignore higher order terms in the derivation (see next section), then one arrives at the following O(1) continuum model for the strain in the purely nonlinear case (if δ 0 = 0) which is the equation derived in [1] (where the small parameter ε is formally set to unity). From a physical standpoint, the reason for doing so is understandable, as the higher order term in the model affords the existence of exact localized traveling wave solutions (which do not exist in the continuum model (4)). An alternative approach towards such a (higher order) quasi-continuum model is the one spearheaded in the earlier work of [12], where one does a similar formal Taylor-expansion based calculation at the displacement level and then differentiates with respect to x to derive the effective long-wavelength equation for the strains r = q x . Remarkably, it should be pointed out, as indicated in [1], that these two procedures (reverting to strain variables and Taylor expanding to go to long wavelengths) do not commute, a feature which poses a mathematical challenge in its own right (about their respective validity). It should be noted that from a visual inspection the solitary wave profiles predicted by (5) (and the corresponding ones of [12]) compare quite well with the "numerically exact" traveling wave solution of the granular crystal model (3) [1]. Nevertheless, without an approximation theorem, it is not clear to what level/extent these kinds of approximations can be used. For example, for the initial value problem (and not just for very special solutions, such as the exact traveling wave ones), there are several nontrivial concerns about models such as those of [1,12] in connection to their local uniqueness and existence properties and the fact the solutions will depend on the small parameter ǫ (hence it is unclear what, if any, influence the O(ǫ) terms have on the dynamics on the O(1) time scale w.r.t. T ). Thus, in terms of providing rigorous error estimates for the initial value problem, model (4) is more appropriate, at least as a starting point or a leading order approximation. Although Eq. (4) fails to describe exact traveling solitary waves (and hence the proper, controllable generalization to higher order so as to capture this important trait is, from a rigorous perspective, an open problem), the model does well in other regards, such as in the description of shock wave formation. This, and other properties of (4) are described in Sections 4 and 7.
Notation. Throughout this paper many possibly different constants are denoted with the same symbol C if they can be chosen independently of the small perturbation parameter 0 < ε ≪ 1.

2.
Fourier transform as fundamental tool. The Fourier transform is the major tool in the proof of the approximation result. In this section we recall some basic facts and establish notation conventions.
Since the solutions u of (3) live on Z we need the Fourier transform on Z leading to periodic functions in Fourier space. Fourier transform on Z: System (3) can be transferred into Fourier space bŷ The inverse of F is given by For every s ≥ 0 the Fourier transform F is continuous from The inverse Fourier transform F −1 is continuous from H s per into ℓ 2 s and from L 1 |û(k, t)|dk and u(·, t) ℓ ∞ = sup n∈Z |u n (t)|.
Fourier transform on R: Beside the Fourier transform on Z we need the Fourier transform on the real line in order to handle the continuum model which lives on the real line. We set The inverse is given by For every s ≥ 0 the Fourier transform F is continuous from The Fourier transform F and its inverse F −1 are continuous from H s to L 2 s and vice versa. Moreover, they are continuous from L 1

3.
Derivation of the continuum model. We assume the existence of a δ 0 > 0 such that inf n∈Z u n (t) ≥ −δ 0 for all t ≥ 0. This is equivalent to requiring that the beads remain in contact for all times. The existence of such solutions will be justified below. With this assumption we may ignore the Heaviside function in the potential, i.e. we have W ′ (u) = (δ 0 + u) p . Thus, we may expand the nonlinearity in a Taylor series, With this expansion we find Taking the Fourier-transform of the right hand side of (10) yields where u * ℓ denotes the (ℓ−1)-times convolution of u with itself, ω(k) 2 = 2(1−cos(k)) and u(k, t) = u(k + 2π, t). In Fourier space, the system (3) is therefore given by We should note here in passing that the use of Fourier space techniques in the context of granular systems was pioneered in the work of [3], where it was used to develop an understanding of the decay properties of traveling waves in the absence of precompression, as well as to offer an efficient numerical tool for computing them. Among the recent ramifications of this approach are the proof of the existence of such bell-shaped traveling waves without [20] and with [19] precompression. The long wave limit ansatz u(n, t) = A(εn, εt) is given in Fourier space by with A : R → C a function decaying to zero for |k| → ∞. Inserting this ansatz into (11), rescaling the integrals, taking formally the limit . Ignoring the higher order terms and taking the inverse Fourier transform of this expression yields our continuum model (4). Note that keeping the next term in the expansion of ω 2 would yield Eq. (5) which, as mentioned above, is not covered by the proof presented herein. Before we turn these formal calculations into rigorous arguments we consider properties of the continuum model (4) itself. 4. Local existence and uniqueness of the continuum equation. We will need a certain regularity of the solutions of the continuum model (4), therefore, we prove the following existence and uniqueness result for the limit equation.
Proof. We know A satisfies . As before, we expand the nonlinearity in a series with real-valued coefficients b ℓ as defined before. The series is convergent for |A| < δ 0 . Applying ∂ −1 X to (4) then multiplying the resulting equation with ∂ −1 X ∂ T A and then integrating w.r.t. X yields Proceeding similarly for the derivatives yields By partial integration we obtain where G s , G s only contain terms with at most s+1 spatial and temporal derivatives. Therefore if A H s+1 and ∂ T A H s are smaller than half of the radius of convergence of the involved series. Note that all series have the same radius of convergence, namely δ 0 . Since a multiple of E s+1 = E 0 + . . . + E s+1 is an upper bound of the squared H s+1 -norm we obtain an estimate Hence by Gronwall's inequality we can guarantee that A H s+1 stays in between half the radius of convergence for all t ∈ [0, T 0 ] if T 0 > 0 and C 1 > 0 are chosen sufficiently small. Since for x ∈ R the sup-norm can be estimated by the H 1 -norm, the second inequality of (13) follows too. Since the previous a priori estimates guarantee that we have a quasilinear system in the sense of [7], the local existence and uniqueness of solutions follows.

5.
Estimates for the residual. For the proof of the approximation result we need a way to measure how ansatz (12) fails to satisfy (3), i.e. we will need estimates for the residual.
It turns out to be advantageous to work in Fourier space, i.e., to work with (11) instead of (3). The error and so where stands for the residual terms, i.e., for the terms which we neglected within the leading order quasicontinuum approximation. In order to solve (14) with boundary conditions R(k, t) = R(k + 2π, t) we have to modify A(k, t) which decays to zero for |k| → ∞. We multiply A(k, t) with a cut-off function χ [−π/2,π/2] . The segment from [−π, π] is then extended periodically with a period 2π to the entire real axis. Call the outcome A # (k, t). In exactly the same way we modify −k 2 in the residual, which is then denoted with −k 2 # .
We will need an estimate of the difference A(k, t) − A # (k, t) and the error ε β R = u − A # which satisfies Thus, in order to bound the error R, we will need estimates for the residual Res # (A # ). Since L 2 is closed under convolution on [−π, π] it turns out to be sufficient to make the estimates in L 2 . We have (4). Then there exist ε 0 , C > 0 such that for all ε ∈ (0, ε 0 ) we have for ε → 0, (note the index of G starts at ℓ = 2). The loss of ε −1/2 comes from the scaling properties of the L 2 -norm. Using that Applying the same argument as above with ( We close this section with an estimate for the difference A(k, t) − A # (k, t).

Proof. We have
For the difference we obtain uniformly in n. The loss of ε −1/2 again comes from the scaling properties of the L 2 -norm. 6. The error estimates. It remains to bound the solutions of (15). In accordance with Lemma 5.1 we choose β = 3/2. We proceed as in Section 4 using energy estimates. From Section 5 we know Define the energy

Making use of the fact
we can compute Note that the autonomous linear terms have canceled, explaining why the sum begins at ℓ = 2. Recall that the application of ω −1 to the residual terms is well defined, see Lemma 5.1. We show below (see e.g. (15)) that the application of ω −1 to ∂ tR andR is also well defined. Using the Plancherel's identity allows us to rewrite where A n (resp. R n ) is the inverse discrete Fourier transform of A # (resp. R) evaluated at n. This motivates the definition of a modified energy which by construction satisfies, Let A(·, T ) ℓ 2 ≤ C 2 and ∂ t A(·, T ) ℓ 2 ≤ C 3 ε and R(·, t) ℓ 2 ≤ C E with C E defined below in (15). Then The first series is convergent for C 1 (and hence C 2 ) sufficiently small. By changing indices, factorizing and ignoring a finite number of terms in the second series results in the expression C 5 (C E ) ≈ δ p 0 (1 + (C 2 + ε β C E )/(δ 0 )) p and thus the series is convergent if C 1 and ε are sufficiently small. Due to the Plancherel identity the energy E 1 is an upper bound for the squared L 2 -norm for C 1 > 0 sufficiently small. Then since ω is bounded and since L 2 ⊂ L 1 on bounded domains (as a consequence L 2 is closed under convolution), we find with constants C j independent of ε. For the inequality above we used the fact that which is a consequence of the Cauchy-Schwarz inequality, Combining this estimate with the one from Lemma 5.2 gives the assertion of Theorem 1.1, completing our proof.
7. Shock Formation in Granular Media. We now briefly discuss some analytical and numerical observations that stem from the leading order quasicontinuum approximation developed in Theorem 1.1. The continuum model can be written as a system of conservation laws These conservation laws have the form of a so-called p-system [10,18]. Defining In both panels it can be seen that the pulse separates into two counter-propagating waves. A zoom of the wave before and after the point of shock development is shown in Fig. 2.
which has the eigenvalues λ ± (U ) = ± p(δ 0 + A) p−1 . Thus, solutions of (4) will consist of two counter-propagating waves traveling with velocity ± p(δ 0 + A) p−1 . Since the wave speed will depend on the amplitude of the solution, bell-shaped initial data will deform and steepen and a shock wave will form in finite time. Discrete shock-like structures (which we will simply call shock waves) have been studied in granular media, e.g. in homogeneous and periodic chains with p = 3/2 and δ 0 = 0 [5,11]. In those works, however, the shock wave is generated by applying a velocity to a single bead [5] or by imparting velocity to the end of the chain continuously [11]. In this paper, the mechanism for the development of the shock wave is fundamentally different. It manifests from an arbitrary non-monotonicallyincreasing initial strain profile under precompression and given a sufficiently long time to develop. A natural question is if the shock wave formation predicted by the continuum model (4) is also present in full system (3). Theorem 1.1 no longer applies in this case, as the the shock wave violates the required smoothness condition. Nonetheless, we carry out a numerical simulation to address the relevant question with a smooth initial condition. Figure 1 shows the development of a shock wave in the discrete model (3)  shape parameters. The wave propagation closely follows the theoretical expectation on the basis of the leading order quasicontinuum approximation. Moreover, the predicted velocity ± p(δ 0 + A) p−1 proved to be very accurate (for example, there was less than a %0.01 relative error in the case shown in Fig. 1). We would like to make a direct comparison of the solutions of the continuum model and the discrete model. However, one has to be careful when using numerical approximations of the continuum model. For example, if one uses a finite difference approximation for a spatial discretization of (4), then we arrive at a model identical to the granular model (3). One could use other numerical schemes, such as those based on adding artificial dispersion [10], but we will proceed in an alternative way to predict the development of a shock wave. Using the velocity relationship ± p(δ 0 + A) p−1 , we can construct the profile of the continuum model for an arbitrary time, starting with the left or right wave (once they are separated) as an initial profile, see Fig. 2. A numerical computation is used to separate the profiles (i.e. we simulate (4) with a finite difference method until separation but before the development of any shock wave, thus avoiding any issues with smoothness). The continuum model predicts a shock wave for any time past the point of non-single-valuedness, as shown in right panel Fig. 2.
In FPU lattices, it is well known that dispersive shocks can develop, in which microscopic oscillations spread out in space and time [6]. From Fig. 2 one clearly sees near the point of wave breaking the development of such oscillations, which are absent in the continuum model [6]. Thus, it would be relevant to extend works like [6] in order to better understand shock waves in the granular crystal model (3). 8. Justification of the KdV and NLS approximation. In this section we would like to contrast the previous result with approximation results for the KdV and NLS approximation. The major difference lies in the ratio between the amplitude and the precompression. For the KdV and NLS approximation this ratio is O(ε 2 ) resp. O(ε) where 0 < ε ≪ 1 is the small perturbation parameter, whereas for the quasicontinuum approximation the ratio is O(1), i.e., of a comparable order.
Since Eq. (10) is exactly of the form of the FPU systems considered in [2,15] the approximation theorems of these papers also apply here. The subsequent solutions u n will be O(ε 2 ), resp. O(ε) and hence will always live in the ball of the convergence with radius δ 0 if the perturbation parameter 0 < ε ≪ 1 is chosen to be sufficiently small.
Even spatially periodic arrangements can be considered here, namely with a n = a n+N for a fixed N . We formulate the relevant approximation theorems in the homogeneous case N = 1. For the KdV approximation we have, with suitable chosen coefficients ν 1 , ν 2 ∈ R. Then there exist ε 0 > 0, C > 0 such that for all ε ∈ (0, ε 0 ) we have solutions (u n ) n∈Z of (18) with Proof. The proof follows trivially from [2, Theorem 3.1] or [16]. We note there is no gap-opening in the small amplitude limit (i.e. the Heaviside function in the potential will play no role).
The theorem can be generalized easily to an approximation theorem for two decoupled KdV equations describing counter-propagating waves, cf. [16].

9.
Conclusions and Future Challenges. The main result of this paper derives in a rigorous way (and with controllable corrections) the leading order quasicontinuum approximation for long wavelength solutions in the granular crystal model (3), in accordance with the formal derivation of Nesterenko [12]. As a technical assumption, we required the presence of a precompression factor (while the original Nesterenko model has been developed also in the case of the so-called "sonic vacuum" i.e., without precompression). One obvious avenue of future research is to investigate a proof without this assumption, which, however, would demand a fundamentally different technique than the one presented herein.
On the other hand, perhaps an even more important aspect of investigation concerns the well-posedness theory of the full Nesterenko model [12] or of the variant developed by Ahnert and Pikovsky in [1]. One important consequence of keeping only the leading order terms in the continuum model (as done herein) is the inability to capture the exact solitary wave solutions (which are known to exist in granular crystals [4,3]) and are at the core of experimental observations in such systems [12,17,8]. The methods of this paper cannot be directly applied to that case, although we should note that we suspect that these higher order long wavelength models suffer (especially so in the case of precompression) from the type of pathologies that were identified by Rosenau and led him to devise appropriate regularizations [13,14]; see also the more recent discussion of [9]. It would be especially relevant to consider such regularizations of the higher order long wavelength models both from a rigorous, as well as from a numerical perspective.
Another important consequence of keeping only first order terms in the continuum model is the prediction of the development of shock waves for a suitable (yet broad) class of initial data. Although the main theorem of this paper does not apply to the case of shock waves, due to smoothness considerations, numerical simulations indicate the steepening of relevant initial data towards a shock structure and suggest that this is indeed an issue worthy of further exploration, with an aim towards transferring these results to the discrete model. Indeed, it is known in FPU lattices that this procedure fails [6], due to the existence of high frequency oscillations (resulting from so-called dispersion shock waves). Thus, a different continuum model (than the one derived herein) will most likely be needed to fully characterize the emerging dispersive shock wave case. These topics are currently under consideration and will be reported in future publications.