N Gorbushin

A Vainchtein

L Truskinovsky

Transition fronts and their universality classes

Steadily moving transition (switching) fronts, associated with local transformation, symmetry breaking, or collapse, are among the most important dynamic coherent structures. The nonlinear mechanical waves of this type play a major role in many modern applications involving the transmission of mechanical information in systems ranging from crystal lattices and metamaterials to macroscopic civil engineering structures. While many different classes of such dynamic fronts are known, the interrelation between them remains obscure. Here we consider a minimal prototypical mechanical system, the Fermi-Pasta-Ulam (FPU) chain with piecewise linear nonlinearity, and show that there are exactly three distinct classes of switching fronts, which differ fundamentally in how (and whether) they produce and transport oscillations. The fact that all three types of fronts could be obtained as explicit Wiener-Hopf solutions of the same discrete FPU problem allows one to identify the exact mathematical origin of the particular features of each class. To make the underlying Hamiltonian dynamics analytically transparent, we construct a minimal quasicontinuum approximation of the FPU model that captures all three classes of the fronts and reveals interrelation between them. This approximation is of higher order than conventional ones (KdV, Boussinesq) and involves mixed space-time derivatives. The proposed framework unifies previous attempts to classify the mechanical transition fronts as radiative, dispersive, topological, or compressive and categorizes them instead as irreducible types of dynamic lattice defects.

I. INTRODUCTION

Transition fronts in discrete systems continue to attract a lot of attention because they represent examples of farfrom-equilibrium collective phenomena that emerge from the underlying many-body interactions. Interpreted as highly nonlinear coherent dynamic structures, such fronts play an important role in the energy transmission from macro to microscales. They are observed in both integrable and nonintegrable Hamiltonian systems [1,2], can be topological or nontopological [3][4][5], spreading or compact [6], compressive or undercompressive (non-Lax) [7], stable or unstable [8]. Together with solitons and breathers, they play a crucial role as building blocks in complex nonlinear wave patterns that emerge generically in mechanical systems ranging from crystals [9][10][11] to nanomechanical structures [12][START_REF] Raney | Proc. Natl. Acad. Sci[END_REF][START_REF] Nadkarni | [END_REF][15].

The concept of transition fronts is equally relevant for the description of pattern formation [16] and transport properties in nonmechanical dynamical systems, including coupled waveguide arrays [START_REF] Ricketts | Electrical Solitons: Theory, Design, and Applications[END_REF][START_REF] Fleischer | [END_REF], quantum systems [19,20], Bose-Einstein condensates [START_REF] Kevrekidis | Emergent Nonlinear Phenomena in Bose-Einstein Condensates: Theory and Experiment[END_REF][START_REF] Mossman | [END_REF][23][24], electronic liquids [25], ultracold quantum gases [26,27], rarefied plasma [28], intense electron beams [29], liquid helium [30], and exciton polaritons [31]. In this paper, we focus on mechanical switching fronts due to the importance of their dynamics for the design of modern metamaterials [START_REF] Raney | Proc. Natl. Acad. Sci[END_REF][32][33][34]. The term "mechanical metamaterials" is used here to describe high-contrast (soft-hard) composite structures with complex architecture at mesoscale. Characteristically, the macroscopic properties of such structures are controlled more by the structural stability of the subelements than by their material properties [33,[35][36][37][38][39][40][41]. The use of additive manufacturing techniques opened a way to exploit various elastic instabilities embedded in the metamaterial response and to creatively guide them using applied deformation [START_REF] Raney | Proc. Natl. Acad. Sci[END_REF]42,[START_REF] Chen | Proc. Natl. Acad. Sci. USA[END_REF]. Dynamic effects targeted by various metamaterial architectures include mitigation of impact loadings, nondestructive detection of inhomogeneities, suppression or amplification of internal instabilities, transmission, guiding, and encryption of mechanical information, including the enabling of logic operations, dynamic unfolding of deployable structures, energy harvesting, and even activating soft robotics [34,[START_REF] Foehr | [END_REF][45][46][47][48][START_REF] Dorin | Active and Passive Smart Structures and Integrated Systems XIII (International Society for Optics and Photonics[END_REF][START_REF] Kidambi | [END_REF][51][52][53][54].

One of the most interesting nonlinear dynamic effects that qualifies metamaterials as mesoscopic analogs of ordered solid-state materials, such as ferroelectrics, ferromagnets, and ferroelastics, is their ability to support moving transition fronts (analogs of domain boundaries), which enable the system to perform dynamic switching between different equilibrium states [32,[55][56][START_REF] Jin | Proc. Natl. Acad. Sci[END_REF][START_REF] Zareei | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Kang | [END_REF][START_REF] Paulose | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Rafsanjani | Proc. Natl. Acad. Sci. USA[END_REF]. There is already a rich body of theoretical and experimental literature devoted to the study of such dynamic snapping/switching waves in mesoscopic mechanical systems [15,[START_REF] Jin | Proc. Natl. Acad. Sci[END_REF][START_REF] Slepyan | [END_REF][63][64]. The ability to propagate transition fronts in metamaterials opens new ways towards potential applications in shape morphing, reconfigurable devices, mechanical logic, and controlled energy absorption [START_REF] Chen | Proc. Natl. Acad. Sci. USA[END_REF][65][66][START_REF] Novelino | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Deng | [END_REF][69]. Analysis of low-dimensional model systems can serve as a guide for the structural design and optimization of the actual three-dimensional (3D) mechanical systems.

Despite the ubiquity of transition fronts in metamaterials, the relation between different classes of such mobile nonlinear dynamic structures remains obscure. In this paper, we consider a well-known prototypical system, namely the FIG. 1. Piecewise linear stress-strain relation σ (ε) defined in Eq. ( 3) and the corresponding energy density φ(ε) for different values of σ .

Here E 1 = 1, E 2 = 4, ε c = 1, σ 0 = 1 (solid), σ 0 = 3 (dotted), and σ 0 = 5 (dashed).

Fermi-Pasta-Ulam (FPU) model [START_REF] Fermi | Studies of the nonlinear problems[END_REF][START_REF] Gallavotti | The Fermi-Pasta-Ulam Problem: A Status Report[END_REF][START_REF] Berman | [END_REF][73], and present a unified description of the three main types of steady transition fronts in this one-dimensional lattice, which we identify as subkinks, shocks, and superkinks. Various realizations of these archetypes have been previously encountered in applications and treated as unrelated: subkinks as subsonic phase boundaries [74][75][76], shocks as classical supersonic shock waves [77,[START_REF] Truskinovsky | Shock Induced Transitions and Phase Structures in General Media[END_REF], and superkinks as supersonic activity waves [54,[START_REF] Gorbushin | [END_REF]. They were first treated as disconnected solutions of the FPU model in [START_REF] Slepyan | [END_REF][START_REF] Slepyan | Models and Phenomena in Fracture Mechanics[END_REF]. Some conceptual links between subkinks and shocks have been previously established in [77,[START_REF] Truskinovsky | Shock Induced Transitions and Phase Structures in General Media[END_REF], while superkinks remain a disconnected class of transition fronts [START_REF] Gorbushin | [END_REF][START_REF] Iooss | [END_REF][82][83].

A unified description of all these transition fronts can be obtained if we use the simplest choice of nonlinearity and assume that the FPU interactions are piecewise linear. In fact, such interactions were already considered in Ref. [START_REF] Fermi | Studies of the nonlinear problems[END_REF] and have since been employed for the description of various dynamic nonlinear phenomena, e.g., [START_REF] Slepyan | Models and Phenomena in Fracture Mechanics[END_REF][84][85][86][87][88].

More specifically, we consider the Hamiltonian dynamics of a mass-spring chain with mass displacements u n (t ) satisfying the infinite system of equations

ρh d 2 u n (t ) dt 2 = σ u n+1 -u n h -σ u n -u n-1 h . ( 1
)
Here h is the equilibrium distance between the masses m = ρh, where ρ is the mass density. In terms of the strain variables

ε n (t ) = u n+1 (t ) -u n (t ) h ,
the equations become

ρh 2 d 2 ε n (t ) dt 2 = σ (ε n+1 ) -2σ (ε n ) + σ (ε n-1 ). ( 2
)
The assumed piecewise linear macroscopic stress-strain relation can be written as

σ (ε) = E 1 ε, ε < ε c , E 2 ε -σ 0 , ε > ε c , ( 3
)
where ε c is the critical (switching) strain, and E 1 , E 2 are the elastic moduli in the two linear regimes. We assume that E 2 > E 1 , so that the two characteristic speeds c 1,2 = E 1,2 /ρ satisfy c 2 > c 1 . The corresponding piecewise quadratic elastic energy density φ(ε) = σ (ε)dε is continuous:

φ(ε) = E 1 2 ε 2 , ε < ε c , E 2 2 ε 2 -ε 2 c -σ 0 (ε -ε c ) + E 1 2 ε 2 c , ε > ε c .
Examples of stress and energy density functions are shown in Fig. 1.

Note that as the stress jump at the critical strain

σ = σ (ε c -0) -σ (ε c + 0) = σ 0 -(E 2 -E 1 )ε c
varies from positive to negative values, we obtain two fundamentally different types of constitutive behavior. Thus, the elastic energy density φ(ε) is nonconvex when σ > 0 and convex for σ 0. When σ > 0, the different branches of the stress-strain curve can be considered as different "phases" of the material, with the spinodal region [where φ(ε) is concave in a smoother setting] represented by the single point ε = ε c . When σ < 0, the stress jump at ε = ε c is just a representation of the hardening-type nonlinearity, which is again concentrated at a single point.

The advantage of the piecewise linear choice for the stressstrain relation is the possibility to construct the corresponding traveling wave solutions of the FPU problem explicitly using the Wiener-Hopf transform technique [START_REF] Slepyan | Models and Phenomena in Fracture Mechanics[END_REF]. While smoothening the constitutive response around the singular point ε c could make the model more realistic, sometimes even without sacrificing much of analytical transparency [89][90][91], stronger nonlinearity is needed to capture such important physical effects as thermalization of the radiated phonons [92][93][94]. However, such a generalization of the model, which will make its analytical treatment almost impossible without contributing much to the classification of the transition fronts, is outside the scope of this paper.

To make the structure of the underlying Hamiltonian dynamics clearly visible, we pose the problem of constructing the minimal quasicontinuum (QC) approximation of the FPU model capturing all three classes of the fronts. The term "quasicontinuum" is used here in the sense that it is a continuum approximation of the discrete system, which is, however, not scale-free and carries a memory about the lattice discreteness [95]. Our analysis shows that the desired approximation must be necessarily of higher order than the conventional ones (KdV, classical "good" or "bad" Boussinesq) and should involve mixed space-time derivatives. The obtained minimal QC model with desired properties can be viewed as a higher-order mixed derivatives (temporal) analog of the "good" Boussinesq approximation [96]. In contrast to the more conventional approach of adding spatially nonlocal terms to the elastic energy [START_REF] Kunin | Elastic Media with Microstructure I: Onedimensional Models[END_REF][START_REF] Charlotte | [END_REF], it introduces the higherorder derivatives into the inertial part of the model (into the kinetic energy), as advocated earlier in [99].

The proposed QC framework not only provides a transparent interpretation of the three types of transition fronts as heteroclinic trajectories of different kinds in the phase space, but it also helps to explain in physical terms why some kinks are radiative (dissipative) while others are not, why some shocks are dispersive while others are not, and why kinks are topological while shocks are not. The comparison with the exact solutions of the discrete problem shows that, on both qualitative and quantitative levels, the relation between different classes of transition fronts is captured adequately by the proposed minimal QC approximation.

It is important to mention that while the nonstationary (spreading) dispersive shock waves (DSWs) [100][START_REF] Chong | Coherent Structures in Granular Crystals[END_REF][START_REF] Benzoni-Gavage | [END_REF] are not the focus of our study, which aims to classify steadily moving transition fronts, we show numerically that the DSWs replace the steady transition fronts in a subdomain of the parameter space. The adequacy of the QC approximation is corroborated by the fact that the DSW stability subdomains in discrete and QC models nearly overlap.

On a theoretical side, our approach unifies the previous attempts to classify the mechanical transition fronts as radiative, dispersive, topological, or compressive and categorizes them instead in a unified framework as fundamentally distinct types of dynamic lattice defects. The obtained analytical solutions can also be used in applications as a guidance in the design of new metamaterials exploiting structural nonlinearity at the scale of the periodicity cell. For instance, our analysis points to a particular type of nonlinearity which should be used if the goal is the suppression of shock loading by channeling the largest amount of energy from macro to microscales. It also makes clear that a different type of nonlinearity must be engineered if the task is to transmit mechanical information with minimal losses. There is of course still a long way from our prototypical 1D designs to the construction of the corresponding 3D mesoscopic composite structures.

The rest of the paper is organized as follows. In Sec. II we formulate the classical continuum approximation of the discrete problem and identify irreducible classes of transition fronts. Then in Sec. III we introduce a nonclassical quasicontinuum approximation of the same discrete problem and construct explicit solutions of the corresponding dispersive traveling-wave problem describing all three distinct types of transition fronts. In particular, we discuss the issues of solution admissibility in the piecewise linear model and the effective energy dissipation in this Hamiltonian framework. We then we present the results of direct numerical simulations that suggest stability of the obtained traveling waves. In Sec. IV we construct an explicit traveling wave solution of the original discrete problem providing a unified description of all three types of fronts. We then present numerical simulations illustrating stability of the different types of transition fronts in various domains of the parameter space. In Sec. V we briefly mention potential applications of our results for the design of metamaterials. A summary of the results and concluding remarks can be found in Sec. VI. Some asymptotic results are presented in Appendix.

II. CONTINUUM MODEL

In our search of a unified description for the different types of transition fronts, it is natural to start with the classical continuum approximation of the original discrete model (1). It can be obtained by taking a formal limit h → 0 and replacing finite differences by the lowest-order derivatives. Following [103], we obtain the standard nonlinear wave equation, which can be represented as the first-order system

∂ε ∂t = ∂v ∂x , ρ ∂v ∂t = ∂ ∂x σ (ε). (4)
Here ε(x, t ) = u x and v(x, t ) = u t are the strain and particle velocity, respectively. The system (4) has discontinuous solutions, which must satisfy the classical Rankine-Hugoniot (RH) conditions

v + V ε = 0, ρV v + σ (ε) = 0, ( 5
)
where V is the velocity of the jump discontinuity. The notation f ≡ f +f -will be used throughout the paper to describe the jump between the limiting values to the right and to the left of a discontinuity.

By changing the parameter σ and varying independently the velocity of the jump discontinuity, we can obtain three fundamentally different types of steadily moving transition fronts shown schematically in Fig. 2. Each transition front connects a state ε = ε + in front with a state ε = ε -behind. Both of these states ε ± belong to the stress-strain curve which is piecewise linear, and to be nontrivial the transition front must connect the states on two sides of the singular point ε = ε c . The RH conditions state that the slope of the Rayleigh line connecting (ε + , σ (ε + )) and (ε -, σ (ε -)) is proportional to the square of the velocity V of the front:

σ (ε + ) -σ (ε -) = ρV 2 (ε + -ε -). ( 6
)
The three different types of transition fronts are defined by the relation between their velocity V and the characteristic velocities c 1 and c 2 , which can be determined by comparing the slopes of the Rayleigh line and the corresponding linear regimes of the stress-strain curve. In what follows, we will refer to them as subkinks [subsonic kinks,

V < c 1 < c 2 , panel (a) of Fig. 2], shocks [intersonic fronts, c 1 < V < c 2 , panel (b)],
and superkinks [supersonic kinks,

c 1 < c 2 < V , panel (c)].

A. Well-posedness

Note that there are five variables to be determined for each discontinuity: v ± , ε ± , and V . Two relations between these five unknowns are furnished by the RH conditions (5). Figure 2 shows qualitatively the fundamentally different relations of this type. Additional information can be obtained by solving the problem (4) using the method of characteristics. Due to the piecewise linear nature of the problem, two families of characteristics with velocities ±c 1,2 can be defined on both sides of the moving front.

Figure 3 shows the arrangement of such characteristics in space-time for all three types of transition fronts. When V < FIG. 2. Rayleigh lines connecting (ε + , σ (ε + )) and (ε -, σ (ε -)) with the slope ρV 2 satisfying Eq. ( 6) for three distinct types of traveling wave solutions: (a) subsonic kinks,

V < c 1 < c 2 ; (b) shocks, c 1 < V < c 2 ; (c) supersonic kinks, c 1 < c 2 < V . The driving force is G = S 2 -S 1 ,
where S 1 (blue) and S 2 (pink) are the areas cut by the Rayleigh line from the stress-strain curve. c 1 (subkinks) or V > c 2 (superkinks), there are two incoming characteristics at the front, which reduces the number of unknowns to one, and therefore an additional condition is needed to find the remaining parameter, for instance V . If c 1 < V < c 2 (shocks), there are three incoming characteristics, which means that all five parameters can be determined without any additional conditions. In this sense, kinks are undercompressive (non-Lax) while shocks are compressive [START_REF] Lefloch | Hyperbolic Systems of Conservation Laws: The Theory of Classical and Nonclassical Shock Waves[END_REF]. The necessity of an additional "kinetic relation" on discontinuous transition fronts was first pointed out in [74,[START_REF] Truskinovsky | Shock Induced Transitions and Phase Structures in General Media[END_REF][START_REF] Truskinovskii | [END_REF]; see also [START_REF] Abeyaratne | Evolution of Phase transitions: A Continuum Theory[END_REF]. The difference between subkinks and superkinks, which both require an additional condition closing the problem, is not apparent in this purely continuum setting.

B. Dissipation rate

While the system of continuum equations ( 4) is conservative, it known that the corresponding discontinuous solutions may be dissipative. One way to supply the missing closure relations for subkinks and superkinks is to specify the dissipation rate at the moving transformation front.

For all three classes of fronts, the energy dissipation on the discontinuity can be written as the product [START_REF] Truskinovskii | [END_REF]

R = GV 0, (7)
FIG. 3. Characteristics η ± (c 1,2 ± V )t = const of the continuum problem in the moving frame with η = x -V t in phase 1 (blue) and phase 2 (red): (a) subkinks,

V < c 1 ; (b) shocks, c 1 < V < c 2 ; (c) superkinks, V > c 2 . Here η = x -V t.
where V is the velocity of the front, and G is the conjugate generalized (or driving) force, which is also known as the energy release rate. After appropriate symmetrization [74], it takes the form

G = φ(ε) -{σ (ε)} ε , ( 8
)
where we introduced a notation for the averaging over the jump

{ f } = ( f + + f -)/2.
The quasistatic notion of a driving force on a moving discontinuity dates back to Eshelby [START_REF] Eshelby | Inelastic Behavior of Solids[END_REF][START_REF] Knowles | [END_REF][109]. A recent application of this notion in inertial dynamics can be found, e.g., in [75].

In our piecewise linear continuum model, the driving force G can be computed explicitly. We obtain

G = E 2 -E 1 2 ε 2 c -ε + ε -+ σ 0 2 (ε + + ε --2ε c ). ( 9
)
In terms of the diagrams in Fig. 2, one can show that G can be represented as the difference between the two colored areas between the Raleigh line and the stress-strain curve: G = S 2 -S 1 . Given that V > 0, the area S 1 (blue) corresponds to the energy rate received on the jump while the area S 2 (red) describes the rate of energy loss. To ensure the overall dissipative nature of the jump encapsulated by the inequality (7), it is necessary that S 2 S 1 .

Note that according to Fig. 2, in the case of subkinks the energy is received at the frontal part and released at the back part of the transition front. Inside shocks the energy can only be released. For superkinks the energy is released in the frontal part and received in the back part.

C. Inner structure of the fronts

As we have seen, in the continuum model the transition region is infinitely localized in space (jump discontinuity). However, the different arrangements shown in Fig. 2 suggest that it may be of interest to reconstruct the energetic structure of each of the archetypal front in the configurational space of strains varying from ε + to ε -. The idea is that the energy transfers implied by the relative size of the areas S 1 and S 2 shown in Fig. 2 are accomplished by some microscopic dispersive mechanisms that are overlooked by the continuum approximation.

(ε, ε + ): (a) subkinks, V < c 1 ; (b) shocks, c 1 < V < c 2 ; (c) superkinks, V > c 2 .
For instance, in the case of subkinks, the continuously emerging energy in the frontal part of the transition region must be somehow transported from the back of the front where it is released. Such transport can be accomplished by the emitted subcontinuum (lattice) waves whose group velocity is larger than their phase velocity (which is equal to V ). In the case of superkinks, the energy released in the frontal part is at least partially reacquired in the back part, and for this the system can use lattice waves whose group velocity is smaller than the phase velocity. To support all three types of the fronts, the dispersion must be sufficiently complex, which is of course the case for the original discrete model.

To support this intuitive picture, it is instructive to introduce the notion of the local energy variation inside the strain interval connecting the limiting states ε + and ε -. Since the actual trajectory in the stress-strain space is not yet known, we can consider energy variation along the Rayleigh line, which ensures the conservation of the macroscopic mass and momentum. The corresponding auxiliary function was introduced in [START_REF] Truskinovsky | Dynamics of Crystal Surfaces and Interfaces[END_REF], and in our notation it takes the form

G(ε, ε + ) = φ(ε) -φ(ε + ) -(ε -ε + ) (ε), where (ε) = σ (ε + ) + ρV 2 2 (ε -ε + )
is the average of σ (ε + ) and the stress taken along the Rayleigh line at ε. One can show that the limiting states ε + and ε - correspond to the extrema of the potential G with respect to ε. Note also that the reference energy is chosen in such a way that

G(ε + , ε + ) = 0,
which means that the energy level assigned to the state ahead of the jump ε = ε + is zero. On the other hand, the overall dissipative (or nondissipative) nature of each type of front is reflected by the fact that at the final state ε = ε -we have

G(ε + , ε -) = -G 0.
In this way, the implied energy landscape describes the energy variation inside the moving front independently of its type. However, it is important to remember that the function G does not describe the actual variation of the energy inside the moving front, as we still do not refer to any particular dispersive mechanisms operating inside the transition zone.

The behavior of G as a function of ε for all three types of transition fronts is shown schematically in Fig. 4. As expected, the ensuing energy landscapes for different universality classes are qualitatively different. Thus, for subkinks, in addition to dissipation, which is expressed by the fact that the minimum at ε -is lower than the minimum at ε + , there is also an energy barrier in between that needs to be overcome. Crossing this barrier requires energy to be continuously transmitted by dispersion from the downstream, where it is continuously released. For shocks, there is no barrier, and the continuously released energy must be fully removed, with none of it being reabsorbed. Finally, for superkinks there is no dissipation (as will be confirmed later). However, in this case there is an antibarrier, and energy transmission by dispersion is still necessary, but now from upstream to downstream. Note also that since the barriers exist in the case of kinks and not shocks, the former can be considered as topological "lattice defects," while the latter remain nontopological.

III. QUASICONTINUUM MODEL

The scale-free approximation we used to obtain the continuum model does not reveal the fate of the energy dissipated on the localized transition front and does not explain which additional macroscopic jump condition must be chosen in the case of subkinks and superkinks. To answer these questions, we must solve the discrete problem. The qualitative information can also be obtained from a quasicontinuum (QC) approximation with sufficiently rich dispersion to adequately mimic the subcontinuum energy transport [95,96].

In this section, we show that the minimal QC approximation of the FPU model capturing all of the dynamic regimes of interest can be constructed following the general approach proposed in [99]. The idea is to focus on temporal dispersion and introduce internal scales into the expression of kinetic energy while keeping the elastic energy as in the scale-free theory. The focus on kinetic energy dates back to the theory of rotational inertia of beams by Rayleigh [START_REF] Rayleigh | The Theory of Sound[END_REF], with subsequent generalizations for other dispersive problems [START_REF] Benjamin | [END_REF]113]. While in the context of discrete lattices the QC theories of this type have been considered before [114][115][116][117][118], we show below that even the minimal QC approach, describing all three universality classes, includes new elements.

A. Main equations

To construct the QC approximation, we set x = nh, and we introduce the variables ε(x, t ) = ε n (t ) and σ (x, t ) = σ (ε n (t )), viewed as functions of continuous space and time. We can then rewrite the infinite system (2) as a single advance-delay partial differential equation, which after the spatial Fourier transform takes the form

ρh 2 ∂ 2 ε ∂t 2 = 4 sin 2 kh 2 σ , ( 10
)
where f (k, t ) = ∞ -∞ f (x, t ) exp (ikx) dx is the Fourier trans- form of f (x, t ).
To simplify the problem and develop the corresponding long-wave asymptotic expansion, we assume that kh 1. To adequately describe the temporal dispersion [99], we use the (2,4) Padé approximation of sin 2 (kh/2) in kh, which affects the kinetic energy, while preserving the classical continuum form of the elastic energy. This yields

4 sin 2 (kh/2) ≈ (kh) 2 1 + a 1 (kh) 2 + a 2 (kh) 4 , ( 11
)
where a 1 = 1/12 and a 2 = 1/240. The need to retain two subcontinuum terms in this approximation is dictated by the requirement that the resulting QC model is both comprehensive and minimal, as will be explained below. We remark that in addition to yielding bounded dispersion relations for the two linear regimes, the expansion in Eq. ( 11) is accurate up to O(k 8 ) near k = 0, and thus provides a long-wave approximation of the discrete Laplacian operator. Other choices of Padé approximations are discussed in [96,[119][120][121]. Substituting the expansion in Eq. ( 11) into Eq. ( 10) and mapping it back into physical space, we obtain, after integration,

ρ 1 -a 1 ∂ 2 ∂x 2 + a 2 ∂ ∂x 4 ∂ 2 u ∂t 2 = ∂σ ∂x , ( 12
)
where u(x, t ) is the displacement field defined by the relation u x = ε; here we also used the scaling x = x/h but dropped the tildes to simplify the expressions. The single partial differential equation ( 12) represents the desired QC approximation of the infinite FPU system (1) of ordinary differential equations.

To reveal the structure of the augmented kinetic energy term, we now derive Eq. ( 12) from the Hamiltonian action principle. We start with the sufficiently general action functional of the form

A = L(u ,i , u ,i j , u ,i jk ) dq 1 dq 2 , ( 13
)
where L is a Lagrangian density, q 1 = x is the spatial coordinate, q 2 = t denotes time, and the subscripts after the comma indicate partial derivatives with respect to q 1 and q 2 . The integration in Eq. ( 13) is over the two-dimensional space-time domain representing the evolving body between the time instants t = t 0 and t = t 1 . The deformation history is described by the function u(q a ), a = 1, 2. Given the structure of the action functional, we can write the Euler-Lagrange equations in the form [START_REF] Gelfand | Calculus of Variations[END_REF] (∂L/∂u ,i -

(∂L/∂u ,i j ) , j + (∂L/∂u ,i jk ) , jk ) ,i = 0. ( 14
)
Here and in what follows the summation over repeated indices is implied. To obtain Eq. ( 12) from Eq. ( 14), we need to specify the Lagrangian. It is not difficult to see that the desired equation will be obtained if we consider the Lagrangian in the form

L = (ρ/2) u 2 t + a 1 u 2 tx + a 2 u 2 txx -φ(u x ). ( 15
)
Here the density of the elastic energy φ(u x ) is the same as in the classical continuum theory, while two subcontinuum terms with mixed derivatives appear in the expression of the kinetic energy. While the "microkinetic" term a 1 u 2 tx is now standard (see, e.g., [116,[START_REF] Theil | [END_REF]), to our knowledge the next term in the expansion, a 2 u 2 txx , has not been used constructively before. The advantage of using the variational principle is that it allows one to derive not only the governing equation but also the corresponding jump conditions. This is relevant because despite regularization provided by the high derivative terms in the energy, our piecewise linear QC theory is still nonsmooth at the transition point u x = ε c . The corresponding generalization of the RH jump conditions, compatible with our higher-order QC theory, emerges as a natural consequence of extremality of the action functional. Indeed, if the spacetime domain contains a surface of discontinuity, the standard Euler-Lagrange equations must be supplemented by the additional necessary conditions of extremality on . In our case, the surface is characterized by the continuity condition u x = ε c , so we must also require that u = 0. While the particle trajectories are differentiable on , some higher-order derivatives of the displacement field may be discontinuous. We interpret the constraints on such singular surfaces imposed by the action principle as the dispersive Rankine-Hugoniot (DRH) jump conditions.

Using the standard manipulations detailed, for example, in [START_REF] Gelfand | Calculus of Variations[END_REF], we obtain ∂L/∂u ,i -(∂L/∂u ,i j ) , j + (∂L/∂u ,i jk ) , jk n i = 0, ( 16)

∂L/∂u ,i j -(∂L/∂u ,i jk ) ,k n i n j = 0. ( 17
) ∂L/∂u ,i jk n i n j n k = 0. ( 18
)
Here n a is the unit vector normal to facing the + direction; the spatial (n 1 ) and the temporal (n 2 ) components of such normal are related through n 2 = -n 1 V , where V is the velocity of the discontinuity. The necessary conditions ( 16), [START_REF] Ricketts | Electrical Solitons: Theory, Design, and Applications[END_REF], and ( 18) of extremality must be supplemented by the kinematic compatibility conditions

u ,i = μn i ,
where μ is a scalar. Eliminating μ, we obtain an auxiliary jump relation

u t + V u x = 0, (19)
which represents the balance of mass across the discontinuity. In our special case, the three DRH conditions ( 16), [START_REF] Ricketts | Electrical Solitons: Theory, Design, and Applications[END_REF], and (18) reduce to

ρV u t -a 1 u txx + a 2 u txxxx + σ (u x ) = 0, ( 20
)
a 1 u tx -a 2 u txxx = 0, ( 21
)
a 2 u txx = 0. ( 22
)
To satisfy all these conditions, we assume that u t = 0 and u tx = 0. Then u x = 0, while the two conditions ( 21) and ( 22) reduce to u txxx = 0 and u txx = 0, respectively. The condition (20) reduces to

ρa 2 V u txxxx + σ (u x ) = 0. ( 23
)
The derived jump conditions guarantee that the physical phenomena in the bulk and on the discontinuity surface are exactly the same.

B. Dimensionless formulation

In what follows, we use dimensionless variables

Ṽ = V c 1 , σ = σ E 1 , σ0 = σ 0 E 1 ,
with tildes dropped to simplify notation. The system is controlled by the dimensionless parameters ε c and

γ = E 2 E 1 > 1.
For the analysis presented below, it is convenient to work with the following equation obtained by differentiating the dimensionless version of Eq. ( 12) with respect to x:

1 -a 1 ∂ 2 ∂x 2 + a 2 ∂ ∂x 4 ∂ 2 ε ∂t 2 = ∂ 2 σ ∂x 2 . ( 24
)

C. Traveling waves

To find steadily moving transition fronts, we seek solutions of Eq. ( 24) in the form of traveling waves:

ε(x, t ) = ε(η), η = x -V t. ( 25
)
We place the front separating two linear regimes at η = 0 and thus require that the following consistency condition is satisfied:

ε(0) = ε c . ( 26
)
Moreover, we consider the solutions admissible only if they satisfy the inequalities

ε(η) > ε c for η < 0, ε(η) < ε c for η > 0. ( 27
)
Since our solutions can be expected to contain phonon radiation at ±∞, we formulate the boundary conditions in the form

ε(η) → ε ± as η → ±∞, (28)
with constant limits ε ± constrained by the standard RH condition ( 6) with stress-strain law given by Eq. ( 3), which in the dimensionless formulation becomes

ε -= (V 2 -1)ε + -σ 0 V 2 -γ 2 . ( 29
)
The angular brackets in Eq. ( 28) denote the average over the largest period of the short-wave oscillations representing phonon radiation; more generally,

ε(η) = lim τ →∞ 1 τ η+τ η ε(ζ )dζ .
The admissibility conditions in Eq. ( 27) require that ε + < ε c and ε -> ε c . Physically, this means that at ε = ε c the moving transition front performs the switching from one branch of the piecewise linear stress-strain curve to another. Substituting Eq. ( 25) into Eq. ( 24), integrating twice, and taking into account the boundary conditions (28), we obtain the ordinary differential equation

V 2 1 -a 1 d 2 dη 2 + a 2 d 4 dη 4 ε(η) = σ (η) + (V 2 -1)ε + , ( 30
)
where

σ (η) = ε(η)H (η) + [γ 2 ε(η) -σ 0 ]H (-η), (31)
and H (η) is the Heaviside function. We also need to apply the following jump conditions at η = 0:

ε = dε/dη = 0, ( 32
)
d 3 ε/dη 3 = 0, d 2 ε/dη 2 = 0. ( 33
)
It is straightforward to check that the condition (23), which takes the form σ (η)a 2 V 2 d 4 ε/dη 4 = 0, is satisfied automatically.

D. Mechanical radiation

Since Eq. ( 30) is piecewise linear, it can be solved explicitly. The analytical solution in each of the two linear regimes can be written as a combination of linear waves whose frequencies and wave numbers satisfy the characteristic equations

ω 2 ± (k) -V 2 k 2 = 0, (34)
where ω + (k) and ω -(k) are the dispersion relations defined by

ω 2 + (k) k 2 = ω 2 -(k) (γ k) 2 = 1 1 + a 1 k 2 + a 2 k 4 (35)
and shown in the insets of Fig. 5. The double root of (34) at k = 0 is responsible for a linear term in the solution, and in view of the assumption of boundedness of the solution, it contributes only constants in each domain of linearity. Due to the even symmetry of the functions ω ± (k), the four nonzero roots of (34), which we denote by k ± j , j = 1, 2, 3, 4, must

satisfy k ± 3 = -k ± 1 and k ± 4 = -k ± 2 .
Therefore, it suffices to seek nonzero roots with Imk > 0 and Rek > 0, where Rek and Imk are real and imaginary parts of k, respectively. The structure of the roots for three different types of fronts is shown in Fig. 5.

Of principal importance for the description of phonon radiation produced by the moving front are the nonzero real roots of (34). The corresponding points of intersection of ω ± (k) and V k are marked in the insets of Fig. 5. When V < 1 (subkinks), a symmetric pair of such roots ±k ± exists for each domain of linearity: when 1 < V < γ (shocks), only the roots ±k -remain, while in the case of superkinks V > γ , there are no nonzero real roots at all. Since each nonzero real root describes energy radiation to and from the moving front, the superkinks can potentially receive but cannot dissipate energy in the form of radiated waves.

To exclude the energy flux from infinity (antidissipation, which can sometimes be interpreted as an AC-type driving [124]), we must impose the radiation conditions disqualifying some of the waves associated with the real roots. In our case these conditions, comparing the velocity of the energy propagation (group velocity) with the velocity of the front, take the form [START_REF] Slepyan | Models and Phenomena in Fracture Mechanics[END_REF]87]

ω + (k) > V, ω -(k) < V. ( 36
)
Since the functions ω ± (k) are known, these conditions are explicit. They leave only one real root related component of the solution in the case of subkinks and shocks.

E. General solution

We observe that the whole configuration of the roots of the characteristic equations (real and complex) changes depending on the values of V . The nonzero roots are given by ±k ± 1,2 , where

k + 1,2 = √ 2 -5 ∓ 5(12 -7V 2 ) V , k - 1,2 = √ 2 -5 ∓ 5(12γ 2 -7V 2 ) V . ( 37
)
More specifically, for the state ahead of the moving front, we have the following three regimes:

k + 1 = ip, k + 2 = s, V < 1, k + 1,2 = ip 1,2 , 1 < V < V * , k + 1,2 = ∓id + f , V > V * . ( 38
)
For the state behind the front, we have the same three regimes but in different V ranges:

k - 1 = iq, k - 2 = r, V < γ , k - 1,2 = iq 1,2 , γ < V < V * * , k - 1,2 = ∓ig + w, V > V * * . ( 39
)
Explicit expressions for the real and positive functions p(V ),

s(V ), p 1,2 (V ), d (V ), f (V ), q(V ), r(V ), q 1,2 (V ), g(V )
, and w(V ) can be extracted from (37). The critical values

V * = 12/7, V * * = γ 12/7 > V *
are the artifacts of the QC approximation, and, as we show below, do not have any fundamental meaning.

Applying the radiation conditions (36) and the boundary conditions (28), we can write the general solutions corresponding to all three types of transition fronts. In particular, in the case of subkinks (V < 1), the solution takes the form

ε(η) = ε -+ B 1 e qη + B 2 cos(rη) + B 3 sin(rη), η < 0, ε + + A 1 e -pη , η > 0. ( 40
)
One can see that for subkinks there is one unknown coefficient on the + side and three on theside. All of them can be found from the consistency, continuity, RH, and DRH conditions. Indeed, the consistency condition (26) and the first of the continuity conditions in Eq. ( 32) yield in this case the relations

ε + + A 1 = ε c = ε -+ B 1 + B 2 . ( 41
)
This allows us to eliminate ε ± . Using the RH condition (29), the second continuity condition in Eq. (32), and the DRH conditions (33), we then obtain the system of linear equations for the coefficients in Eq. ( 40):

-C 0 A 1 + B 1 + B 2 = b, pA 1 + qB 1 + rB 3 = 0, p 2 A 1 -q 2 B 1 + r 2 B 2 = 0, p 3 A 1 + q 3 B 1 -r 3 B 3 = 0,
where

C 0 = V 2 -1 V 2 -γ 2 , b = (1 -C 0 )ε c + σ 0 V 2 -γ 2 .
The system yields explicit expressions for the four unknown coefficients A 1 , B 1 , B 2 , and B 3 as functions of V that are not provided here to simplify the exposition. The expressions for ε ± (V ) are then found from Eq. ( 41).

For shocks and superkinks (V > 1), the structure of the roots in Eqs. ( 38) and ( 39) changes depending on the value of V relative to the thresholds V * and V * * . To account for this, it is convenient to introduce the shortcuts

λ 1,2 = -p 1,2 , 1 < V < V * , -d ± i f , V > V * and μ 1,2 = q 1,2 , γ <V < V * * , g ± iw, V > V * * .
Then for shocks (1 < V < γ ) we have

ε(η) = ε -+ B 1 e qη + B 2 cos(rη) + B 3 sin(rη), η < 0, ε + + A 1 e λ 1 η + A 2 e λ 2 η , η > 0, (42)
with two unknown coefficients on the + side and three on theside (for V > V * the two nonconstant terms at η > 0 are complex conjugate since λ 2 = λ1 and A 2 = Ā1 ). The conditions ( 26), ( 29), (32), and (33

) yield ε + + A 1 + A 2 = ε c = ε -+ B 1 + B 2
and the following linear system for the coefficients in Eq. ( 42):

-C 0 (A 1 + A 2 ) + B 1 + B 2 = b λ 1 A 1 + λ 2 A 2 -qB 1 -rB 3 = 0, λ 2 1 A 1 + λ 2 2 A 2 -q 2 B 1 + r 2 B 2 = 0, λ 3 1 A 1 + λ 3 2 A 2 -q 3 B 1 + r 3 B 3 = 0. ( 43
)
This system of four equations does not allow one to find all five unknown coefficients A 1 , A 2 , B 1 , B 2 , and B 3 as functions of V . In other words, the structure of shocks is not fully determined internally, which in turn means that ε ± cannot be determined as functions of V . All parameters are fully defined in this case only if we provide one additional external condition, for example ε + = 0, which implies

A 1 + A 2 = ε c .
In the range V > γ (superkinks), the solution reads

ε(η) = ε -+ B 1 e μ 1 η + B 2 e μ 2 η , η < 0, ε + + A 1 e λ 1 η + A 2 e λ 2 η , η > 0, (44)
with

λ 2 = λ1 , A 2 = Ā1 for V > V * and μ 2 = μ1 , B 2 = B1 for V > V * * .
In this case, there are two unknown coefficients on each side of the front, so the solution is again fully specified by conditions ( 26), ( 29), (32), and (33), which yield the linear system

-C 0 (A 1 + A 2 ) + B 1 + B 2 = b, λ 1 A 1 + λ 2 A 2 -μ 1 B 1 -μ 2 B 2 = 0, λ 2 1 A 1 + λ 2 2 A 2 -μ 2 1 B 1 -μ 2 2 B 2 = 0, λ 3 1 A 1 + λ 3 2 A 2 -μ 3 1 B 1 -μ 3 2 B 2 = 0 for the four unknown coefficients A 1 , A 2 , B 1 , B 2 that
can be found as explicit functions of V , as well as the relations ε

+ + A 1 + A 2 = ε c = ε -+ B 1 + B 2 ,
which allows one to find the two remaining functions ε ± (V ).

To summarize, after using the conditions ( 36), ( 28), ( 26) and the first condition in Eq. ( 32), we are left in the range V < 1 (subkinks) with one unknown coefficient on the + side and three on theside (a single exponential boundary layer and a radiated wave). All of them can be found from the four conditions: the second condition in Eqs. ( 32), (29), and (33).

When 1 < V < γ (shocks) we are left with two coefficients on the + side and three on theside (a radiated wave and a single exponential boundary layer) and only four conditions. This leaves one of the constants in the corresponding linear system (43) undetermined. Finally, in the range V > γ (superkinks) there are two coefficients on each side, so the solution is again fully specified by the four conditions.

Once the strain field is determined in each regime, particle velocity is found from v(η) = -V ε(η).

F. Discussion

Now that the mathematical structure of traveling wave solutions is well understood, we provide a physical interpretation of the results that furnishes a somewhat more intuitive explanation of the fundamental differences between the three types of transition fronts.

Observe first that in all three cases, the traveling wave solutions describing the transition fronts can be written in the same general form

ε(η) = ε ± + ± (η) + ± (η), η ≷ 0. ( 45
)
Here the functions ± (η) depend on the real roots of the characteristic equation and describe the radiative part of the solution. The functions ± (η) depend on the nonreal complex roots and describe the exponentially localized boundary layers on both sides of the moving fronts. The constant terms in (45) are due to the double root at the origin; the strains ε ± correspond to the averaged states at η → ±∞ and satisfy the classical RH condition (29). We now consider in more detail the radiative component of the solution ± (η). We have seen that to exclude the energy flux from infinity (radiation condition), we need to set (in all three cases) that radiation is absent ahead of the front, so that

+ (η) = 0.
Moreover, while all three solutions obtained above in Eqs. (40), (42), and (44) have the form (45), the nontrivial radiation component (behind the moving front) exists only for subkinks and shocks and can be written as

-(η) = 2α -cos (rη + β -), (46)
with α -, β -expressed in terms of B 2 and B 3 in Eqs. (40) and (42). Thus, both subkinks and shocks radiate (dissipate) energy. In contrast, the superkinks are completely free from radiation (dissipation), since in this case we also have -(η) = 0. We now turn to the boundary layer terms ± (η). For subkinks they involve a single decaying exponential term on each side of the front [ + (η) = A 1 e -pη , -(η) = B 1 e qη ]; see Eq. (40). For shocks, there is a single exponential decay behind the front [ -(η) = B 1 e qη ], while ahead of it the decay is double exponential [ 42)]. For superkinks, there is a similar transition from double exponential to oscillatory decay ahead of the front at V = V * if γ < V * and behind it at V = V * * [see Eq. ( 44)]. As the analysis of the discrete problem presented below shows, both double exponential and oscillatory decays are artifacts of the chosen QC approximation.

+ (η) = A 1 e -p 1 η + A 2 e -p 2 η ] when 1 < V < V * and oscillatory ( + (η) = 2e -dη [Re(A 1 ) cos( f η) - Im(A 1 ) sin( f η)]) when V * < V < γ [see Eq. (
As we have seen, for both types of kinks all parameters of the traveling wave, and in particular the limiting states ε ± , are fully determined by the front velocity V . This means that the kinetic relations G = G(V ), whose absence in the classical continuum description produced the fundamental ill-posedness of the problem, are now fixed through the recovery of the internal structure of the kinks. In other words, such fronts are autonomous in the sense that their kinetics is controlled by the microscopic dispersion. For instance, if the state in front of the moving kink ε + is known, then both the state behind, ε -, and the velocity of the front V are determined.

In contrast, in the case of shocks, the knowledge of V is not sufficient to determine both ε ± , and one of the limiting strains remains as a free parameter. As a result, no particular kinetic relation in the form G = G(V ) emerges from the reconstruction of the internal structure of the transition front. In other words, in the case of shocks, the knowledge of the state ahead is not sufficient for complete specification of the remaining parameters and for fixing the internal structure of the transition. This means, for instance, that in addition to the state ahead of the front ε + , another piece of information has to be prescribed by the external (non-traveling-wave) solution in order to make the front velocity V known.

G. Characteristics

The obtained QC picture is in full agreement with what we have learned by studying the classical continuum approximation in Sec. II. There we found that kinks are dissimilar from shocks primarily due to the difference in the number of incoming characteristics shown in Fig. 3.

In particular, Fig. 3 shows that for both types of kinks two characteristics are bringing information to the front. Since in our analysis of the internal structure of the transition fronts we eliminated particle velocities v(η), we may always assume that this information concerns the limiting values v ± . Therefore, we can conclude that in the case of kinks, no information about one of the limiting strains ε ± is arriving from outside. Thus, to fix the unknown limiting strain and to ultimately specify the front velocity V , the system must rely exclusively on the internal dispersive machinery. The analysis of the QC approximation shows that such machinery is indeed in place delivering all of the unknown quantities.

In contrast, in the case of shocks, the classical continuum model tells us that the three characteristics are coming from outside. Therefore, the system can use one additional piece of external information to fix the limiting strains ε ± and to specify the front velocity V . In this case, the internal dispersive structure of the front does not have an autonomy and simply adjusts to the conditions imposed from the outside. Consistently, this is exactly what our study of the dispersive QC model has shown: for shocks the internal traveling wave solution is (one-parameter) underdetermined, and to make the global problem well-posed a single additional piece of information is needed. Such information is then naturally provided by the additional incoming characteristic that does exist in the case of shocks.

H. Dynamical system

Since all three types of transition fronts represent traveling wave solutions of the fourth-order ordinary differential equation (30), it is of interest to examine them from the point of view of the theory of dynamical systems. In this perspective, they emerge as fundamentally different types of heteroclinic trajectories connecting various types of attractors in the four-dimensional phase space. The nature of such attractors depends on the structure of the roots of the characteristic equations, which control the asymptotic behavior of the heteroclinic trajectories as η → ±∞. The knowledge of these asymptotics is sufficient to distinguish between the different universality classes of the transition fronts.

For example, in the case V < 1 (subkinks) the transition fronts correspond to heteroclinic trajectories of the type center-saddle to center-saddle. Such transitions are nongeneric and are possible due to the sufficiently high dimensionality of our dynamical system. More specifically, they are captured by our QC approximation because the latter includes the minimal number of the higher-order dispersive corrections to the classical continuum model, which makes the corresponding phase space four-dimensional. At η = -∞, the heteroclinic trajectory describing subkinks unwinds as the center-related separatrix. The corresponding two-dimensional center effectively describes the radiation behind the moving subkink, while the saddle-related component of the asymptotics describes the exponential boundary layer. At η = +∞ this trajectory ends as a saddle-related separatrix, which describes the exponential boundary layer ahead of the moving front.

Similar considerations can be applied to shock and superkink trajectories. For simplicity, we assume in what follows that γ < √ 12/7 and V < √ 12/7. This eliminates the oscillatory decay for shocks and superkinks, which, as we have discussed, is an artifact of the QC approximation. In the range 1 < V < γ (shocks), the corresponding heteroclinic orbits are of the type center-saddle to saddle-saddle. Such transitions are generic. At η = -∞, the heteroclinic trajectory unwinds as a center-related separatrix describing radiation behind the front. The center-related part of the asymptotics describes the exponential boundary layer. At η = +∞ the trajectory ends as a saddle-related separatrix describing the exponential decay ahead of the front. Finally, for V > γ (superkinks) the corresponding orbit is of saddle-saddle to saddle-saddle type. Such transitions are nongeneric. In this case, the heteroclinic trajectory starts as a saddle-related separatrix describing the exponential decay behind the front and ends as a saddle-related separatrix describing the boundary layer ahead of the front.

We have thus confirmed that the physical nature of all three types of the transition fronts described by the general Eq. ( 12) is fully consistent with the asymptotic behavior of the heteroclinic trajectories at η → ±∞. The fact that the latter is controlled by the structure of the roots of the characteristic equations characterizing the corresponding attractors goes beyond the adopted piecewise linear approximation of the stress-strain relation. Thus, even without such an assumption, the subkinks can be expected to correspond to nongeneric transition fronts that are described by center-saddle to center-saddle trajectories and generate their own kinetic relations. Such transitions, however, would be possible only if sufficiently higher-order dispersion is included into the model. Similarly, even in a smoother model shocks correspond (under our assumptions) to the heteroclinic orbits that are generic saddle-saddle to center-saddle trajectories, and therefore do not generate any specific kinetic relations. Finally, under the same assumptions superkinks are nongeneric transitions described by the saddle-saddle to saddle-saddle heteroclinic orbits. The fact that all possible types of sufficiently low-dimensional nondissipative attractors are accounted for suggests that the proposed classification of the transition fronts is exhaustive.

I. Dissipation rate

In the dispersively regularized setting, the jump discontinuities of strain and velocity that are present in the classical continuum theory are replaced by the extended transition zones. In addition, the energy released on such jumps in the continuum theory no longer disappears locally. Instead it is channeled by nonlinearity from long to short waves and radiated away from the moving front in the form of lattice waves. In the piecewise linear theory, it is transported by such waves to infinity. In other words, despite the absence of explicit damping, the effective dissipation takes place due to the energy escape through phonon radiation.

The developed QC model allows one to trace all these processes in full detail. In particular, one can compute explicitly the thermodynamic driving force G for all three types of transition fronts and determine the corresponding rate of energy dissipation R = GV 0. Based on the analysis of the corresponding modes of radiation, one can see that R is strictly positive for subkinks and shocks but equals zero for superkinks.

More specifically, depending on the structure of the real roots of the characteristic equations, the transition front may or may not emit elastic waves. In general, we have R = R + + R -, where

R + = k∈N + E + (k) [ω + (k) -V ] = G + V, R -= k∈N - E -(k) [V -ω -(k)] = G -V, ( 47
)
and G + and G -are the cumulative energy fluxes associated with emitted elastic waves ahead and behind the front, respectively. Here N ± = {k : Imk = 0, Rek > 0, ω ± (k) = V k, ω ± (k) ≷ V } is the set of positive real roots of the characteristic equation for the corresponding linear regime that satisfy the radiation conditions (36), and E ± (k) are the energy densities associated with the corresponding modes, averaged over the corresponding time period T = 2π/ω ± (k), with f = T -1 T 0 f dt. The energy is transported away from the front with relative velocities ω ± (k) -V [START_REF] Brillouin | Wave Propagation in Periodic Structures: Electric Filters and Crystal Lattices[END_REF].

From the structure of the exact solutions of the QC model, one can see that the set N + is empty for all transition fronts. Thus, independently of the front type there is no radiation of phonons ahead of the front, and G + = 0. In the superkink regime, N -is also empty, and therefore G -= 0 as well, yielding R = 0. In the case of subkinks and shocks there is a single emitted lattice wave mode with wave number r > 0 propagating in the region η < 0, so that N -= {r}. The associated energy with the density

E -(r) = V 2 2 ( 2 -+ a 1 ( -) 2 + a 2 ( -) 2 ) + γ 2 2 2 -,
averaged over the period 2π/ω -(r), is transported backwards relative to the moving front with the relative velocity ω -(r) -V [START_REF] Brillouin | Wave Propagation in Periodic Structures: Electric Filters and Crystal Lattices[END_REF]. This yields the driving force G = G -+ G + given by

G = G -= 2γ 2 (α -) 2 ω 2 -(r) 1 - ω -(r) V > 0,
where we recall that α -is half of the amplitude of the radiation contribution to the solution defined in Eq. ( 46) and can be obtained from Eqs. ( 40) and ( 42) for subkinks and shocks, respectively. The difference is that for subkinks the function G(V ) is fixed, while for shocks we obtain a one-parametric family of such functions.

J. Admissibility

We recall that the explicit expressions for the general solution of the piecewise linear problem are invalid if the admissibility conditions (27) are violated. Therefore, the inequalities ε(η) > ε c for η < 0 and ε(η) < ε c for η > 0 must be checked a posteriori, which means that some of the formally constructed solutions may have to be discarded [84,[START_REF] Marder | [END_REF].

The analysis of the global behavior of the obtained strain fields shows that all subkinks with V < 1 and all superkinks with V > γ are automatically admissible. In both of these cases, the transition fronts can be represented in the space of parameters ε + and σ by one-dimensional manifolds because the velocity of the front is determined uniquely by the corresponding kinetic relation. In the case of shocks, which can be either admissible or inadmissible, the velocity V is not determined internally. Therefore, shocks occupy a two-dimensional (2D) domain in the (ε + , σ ) plane. This domain is further divided into two subdomains: at sufficiently large values of σ shocks are admissible, while those located below a certain threshold are inadmissible. The inadmissible shocks show the repeated crossing of the ε c threshold by the oscillatory tail behind the moving front.

The admissibility diagram in the (V, σ ) plane is shown in Fig. 6, where we fixed ε + = 0. The insets illustrate the analytical solutions describing different types of transition fronts. The 2D domain of shocks on this diagram is bounded on two sides by the condition 1 < V < γ and from below by the dotted line below which ε -< ε c . One can see that only the shock solutions in the pink (upper) region above the threshold values σ * (V ) marked by a solid black curve are admissible, while the ones in the blue (lower) region are inadmissible. This is illustrated in the corresponding inset by the multiple FIG. 6. Admissibility sets of solutions of the QC problem. In the blue region we have ε(η) ε c for some intervals of η < ε c , and its dotted lower boundary marks the threshold ε -= ε c . The insets show examples of the strains ε(η), with horizontal lines marking ε = ε c . Here γ 2 = 1.5, ε c = 1, and we set ε + = 0. crossings of ε c (the dash-dotted horizontal line) by the strain profile ε(η).

To understand which solutions replace shocks in the "forbidden" region, we need to resort to simulations. Using direct numerical simulations of Eq. ( 12) for a sufficiently broad set of initial data, we can also numerically test the stability of the admissible transition fronts.

K. Numerical simulations

We solve Eq. ( 24) in the finite domain x ∈ (0, 200) with the Riemann-type initial data

ε(x, 0) = ε l , x < 100, 0, x 100, ∂ε ∂t (x, 0) = 0.
using the implicit fourth-order conservative finite-difference method developed in [127]. The first and second spatial derivatives of strain are set to 0 at the boundaries. The emergence of particular transition fronts, as an outcome of the breakdown of the initial state, will then depend on the choice of the parameters σ and ε l .

The results are summarized in Fig. 7, which shows time snapshots near the end of four different simulations. In each simulation we have chosen a particular set of parameters ε l and σ to reach one of the four structurally dissimilar regimes shown in Fig. 6.

While in all presented snapshots we observe complex breakdown patterns, most of their elements correspond to linear dispersive pulses with their characteristic overshoots. To identify genuinely nonlinear substructures, one needs to look for the patterns magnified in the insets in Fig. 6. Thus, the inset in Fig. 7(a) shows an admissible subkink moving to the right. The comparison of the internal structure of such numerically generated wave profile with the corresponding analytical solution shows perfect agreement, which confirms that the transformation fronts of this type can indeed serve as dynamical attractors. Similarly, the inset in Fig. 7(d) shows an admissible superkink moving to the right, which also matches the analytical waveform and points towards stability of the corresponding traveling wave solution. An admissible shock is shown in the inset of Fig. 7(b), and we again see that the analytical profile is reproduced faithfully and conclude that such transition fronts can be stable. The remaining panel (c) of Fig. 7 corresponds to parameter values that target inadmissible shocks. Not surprisingly, we do not observe a traveling wave profile in this case. Instead, the nonlinear structure that we see is reminiscent of a nonsteady dispersive shock wave (DSW).

Our broader numerical experiments strongly suggest that, in the whole domain of nonadmissibility, shock traveling waves are replaced by DSWs. This result, obtained so far only in the QC setting, will be confirmed below by a similar analysis of the original discrete problem. We recall that DSWs have been extensively studied using various other QC approximations of the FPU system (see, for example, [100,[START_REF] Benzoni-Gavage | [END_REF][128][129][130]). We conclude that in our regime diagram shown in Fig. 6, the domain of inadmissible shocks should be interpreted as a domain of stability of DSW-type nonsteady (spreading) transition fronts. The absence of steadily moving shock fronts in the FPU model with convex energy density ( σ 0 in our problem) is well known. It has been previously linked to the low dimensionality (lack of transversal radiation) and the absence of irreversibility (purely elastic constitutive modeling), which is a ubiquitous feature of the real crystals [131][132][133]. Here by allowing regimes with σ > 0 we acquire a limited parametric domain where stable stationary shocks exist. One can argue that the nonconvexity, which allows the one-dimensional system to accommodate large-amplitude lattice waves transmitting radiated energy away from the moving front, is the way to bring multivaluedness into the constitutive response, which ul-timately imitates the inherent multistability of the plastic response.

To summarize, the analysis of the dispersively regularized QC model allowed us to clarify the ambiguities left by the classical continuum description. In such an essentially microscopic model, all three classes of transition fronts acquired their natural raison d'être, with the numerical simulations providing confirmation of the exhaustiveness for the proposed classification. It is rather remarkable that such a task could be accomplished using a relatively simple QC approximation of the original discrete problem. Note, however, that the chosen approximation was not of the lowest order, and to capture the complete picture we had to introduce two internal time scales and modify the kinetic rather than elastic energy. As we show in the next section, the obtained description is fully adequate when compared to the discrete model.

IV. DISCRETE MODEL

We now analyze the dimensionless version of the original FPU problem (2), which takes the form

d 2 ε n (t ) dt 2 = σ (ε n+1 ) -2σ (ε n ) + σ (ε n-1 ), ( 48
)
with bilinear interactions σ (ε) = ε at ε < ε c , and σ (ε) = γ 2 ε -σ 0 at ε > ε c . The dispersion relations in each linear regime are defined by

ω 2 + (k) = ω 2 -(k)/γ 2 = 4 sin 2 (k/2) ( 49
)
and are much more intricate than in the QC model due to the presence of lattice resonances and the richness of the spectrum of available lattice-scale waves. Therefore, the analysis of the discrete problem can potentially challenge the description of the energy radiation provided by the QC model. To find the corresponding traveling wave solutions ε n (t ) = ε(η), η = n -V t of the discrete problem (48), we need to solve the advance-delay equation

V 2 d 2 ε dη 2 = σ (η + 1) + σ (η -1) -2σ (η), (50)
where the function σ (η) = σ (ε(η)) is given by Eq. ( 31). We will use the Fourier transform technique to solve Eq. ( 50) subject to the consistency condition (26), the boundary conditions (28), and the radiation conditions (36).

It is convenient to represent the transformed function in the form

ε(k) = ∞ -∞ ε(η)e ikη dη = ε+ (k) + ε-(k), where ε± (k) = ∞ -∞ ε(η)H (±η)e ikη dη
are analytic in Imk ≷ 0. The Fourier transform of ( 50) then yields where we introduced the parameter

M + ε+ + M - ε-= M --M + ik ε * , (51)
ε * = σ 0 γ 2 -1 ( 52
)
and the characteristic functions

M ± (k) = ω 2 ± (k) + (0 + ikV ) 2 . ( 53
)
Here 0 ± ikV = lim s→0+ (s ± ikV ), and we use the causality principle [START_REF] Slepyan | Models and Phenomena in Fracture Mechanics[END_REF] to handle the zero at the origin. A comparison of the characteristic functions (53) with their QC analogs in the whole complex plane shows that while the discrete dispersion relations ( 49) are more complex than their QC counterparts (35), the QC approximation captures the long-wave behavior adequately. More precisely, as shown in Fig. 8, the QC model gives an excellent approximation of the real and purely imaginary roots of Eq. ( 53) that have sufficiently small magnitude.

In general, it captures the four nonzero roots of each characteristic function that are closest to k = 0 qualitatively well but may represent purely imaginary roots by complex quadruples and vice versa.

A. Characteristic roots

Similar to the QC model, the solution of the discrete problem can be written in terms of elementary waveforms associated with the roots of the characteristic functions (53). In what follows, we consider the generic case when

V is nonresonant [V = ω + (k) and V = ω -(k) for any real k]. We can then define the sets Z = Z + r ∪ Z - r ∪ Z + c ∪ Z - c and P = P + r ∪ P - r ∪ P + c ∪ P - c containing nonzero roots of the charac- teristic equations M ± (k) = 0. Here Z ± r = {z : M + (z) = 0, z = 0, Imz = 0, ω + (z) ≷ V }, P ± r = {p : M -(p) = 0, p = 0, Imp = 0, ω -(p) ≷ V }, Z ± c = {z : M + (z) = 0, Imz ≶ 0}, P ± c = {p : M -(p) = 0, Imp ≶ 0}. ( 54
)
The structure of the roots of Eq. ( 53) is illustrated in Fig. 9, which can be compared the corresponding root structure for the QC model shown in Fig. 5 (see also Fig. 8, which compares the structure of real and purely imaginary in that case, the even symmetry of characteristic function implies that the roots are symmetric about the origin, and it suffices to consider the region Rek 0 and Imk 0.

Of particular importance are the sets of nonzero real roots

Z + r ∪ Z - r [roots of M + (k)] and P + r ∪ P - r [roots of M -(k)]
. As we will see, some of these roots correspond to radiated lattice waves. When these sets are nonempty for given nonresonant V , they contain an odd number of positive real roots, equal to 2l + 1 and 2m + 1, respectively. We arrange these roots in ascending order: z j < z j+1 , j = 1, . . . , 2l, and p j < p j+1 , j = 1, . . . , 2m.

We observe that in the case of superkinks (V > γ ) both functions M ± (k) have no nonzero real roots, as shown in Figs. 8(a) and 9(c), and hence there are no radiated waves in this case (no dissipation). For shocks (1 < V < γ ) only M -(k) has such roots [see Figs. 8(a) and 9(b)]. More specifically, we have m = 0 (i.e., one positive real root) for the values of velocity V above the first resonance velocity V 1 , which equals ω -(k) for some real k and satisfies the condition V 1 > 1 for large enough γ > 1. We then have m = 1 (three positive real roots) for the values of V between the first and second resonance velocities, where the second resonance velocity is defined accordingly, and so on. Finally, for subkinks (V < 1) each of the characteristic equations has at least one positive real root [see Figs. 8(a) and 9(a)], with l and m each increasing by 1 when the corresponding resonance velocity is crossed.

In addition to real roots, there are infinite sets of complex roots

Z + c ∪ Z - c [roots of M + (k)] and P + c ∪ P - c [roots of M -(k)
] with a nonzero imaginary part that can be seen in Fig. 9. These roots bifurcate from the maxima of the real-root curves shown in Fig. 8(a). This includes purely imaginary roots that bifurcate from the sonic maxima at k = 0 and are shown in Fig. 8(b). The nonreal roots define the structure of the boundary layers on both sides of the moving front.

B. Characteristics revisited

To make a connection with the classical continuum theory, we recall that the configuration of the real roots z j and p j around the origin k = 0 is intimately related to the structure of the characteristics in the continuum approximation. Therefore, by studying these roots one can expect to reconstruct the main subdivision of the transformation fronts into the three universality classes.

More precisely, we can exploit the fact that in the longwavelength limit, the discrete problem can be replaced by a single nonlinear wave equation. Indeed, in the limit k → 0, s → 0+ we can approximate the linear operators in Eq. ( 53) by

M + (k) = ω 2 + (k) + (s + ikV ) 2 ≈ g + (k, s) ≡ [(1 + V )(-ik) -s][(1 -V )(-ik) + s], M -(k) = ω 2 -(k) + (s + ikV ) 2 ≈ g -(k, s) ≡ [(γ + V )(-ik) -s][(γ -V )(-ik) + s]. (55)
Observe also that using the convective coordinate η = x -V t, we can rewrite the system (4) as a pair of linear wave equations for ε(η, t ) in each of the two domains of linearity:

(1 + V ) ∂ε ∂η - ∂ε ∂t (1 -V ) ∂ε ∂η + ∂ε ∂t = 0, η > 0, (γ + V ) ∂ε ∂η - ∂ε ∂t (γ -V ) ∂ε ∂η + ∂ε ∂t = 0, η < 0. (56)
Applying Fourier transform in η and Laplace transform in t turns Eq. ( 56) into the equations g ± (k, s) = 0, where the functions g ± (k, s) are defined in Eq. ( 55). Since the characteristics of Eq. ( 56) are defined by the equations η ± (1 ± V )t = const at η > 0 and η ± (γ ± V )t = const at η < 0, the location of the roots of the functions g ± (k, 0) is directly linked to the configuration of the characteristics relative to the line η = const. The configuration of the roots of the equations g ± (k, 0) = 0 is shown schematically in Fig. 10 separately for each class of the transition fronts. One FIG. 10. Schematic presentation of the roots g + (k, 0) (blue triangles) and g

-(k, 0) (red circles): (a) V < 1, (b) 1 < V < γ , (c) V > γ .
can see that in the range V < 1 (subkinks) the purely imaginary roots are located in two different complex half-planes for both g + (k, 0) = 0 and g -(k, 0) = 0. This is equivalent to the fact that there is one incoming and one outgoing characteristic on both sides of the line x -V t = const. Both roots of the equation g + (k, 0) = 0 end up in the upper complex half-plane in the range 1 < V < γ (shocks), producing two incoming characteristics on the right side of the line x -V t = const, while there is still one incoming and one outgoing characteristic on the left side. Finally, in the range V > γ (superkinks), the remaining roots of g -(k, 0) = 0 also shift into the upper complex half-plane, which produces two outgoing characteristics behind the moving front. One can see that the location of the roots in Fig. 10 is in full agreement with the propagation direction of the macroscopic perturbations with respect to the moving front for each of our universality classes, as shown in Fig. 3.

C. Solution of the discrete problem

We observe that ε± (k) can be written as

ε± (k) = ε ± 0 ∓ ik + χ ± (k), (57)
where the first term accounts for the boundary conditions (28), and the second term satisfies lim k→±i0 χ ± (k) = 0, so that lim η→±∞ ε(η) = lim k→±i0 ε± (k) = ε ± . To find χ ± (k), we use the Wiener-Hopf technique [START_REF] Slepyan | [END_REF]77,86,88]. To this end, we factorize the main linear operator

L(k) = M + (k) M -(k) = ω 2 + (k) + (0 + ikV ) 2 ω 2 -(k) + (0 + ikV ) 2 (58)
of the problem, which means representing it in the form

L(k) = L + (k)L -(k), (59)
where the superscripts ± identify functions that are regular (have no zeros or singularities) in Imk ≷ 0, respectively. Such factorization allows us to rewrite (51) as

L + (k) -ε * -ik χ + (k) + ε + 0 -ik = 1 L -(k) ik χ -(k) + ε - 0 + ik -ε * . ( 60
)
This representation ensures that the right-hand side is regular in the lower half-plane, while the left-hand side is regular in the upper half-plane, so that both can be analytically continued to the whole plane after we move the zeros and singularities on the real axis into the corresponding half-planes.

Using the infinite product theorem [START_REF] Noble | Methods Based on the Wiener-Hopf Technique for the Solution of Partial Differential Equations[END_REF], we can represent L ± (k) as follows [START_REF] Slepyan | [END_REF]:

L ± (k) = l ± (k)L 0 ± (k). (61)
Here the terms l ± (k) depend on nonzero real roots of the characteristic equations, while the terms L 0 ± (k) are defined by the remaining nonreal (complex) roots.

More specifically, we have

L ± 0 (k) = 1 -V 2 γ 2 -V 2 z∈Z ± c 1 -k z p∈P ± c 1 -k p , ( 62
)
where the products are over the sets Z ± c and P ± c of nonreal roots defined in Eq. (54). Note that the zeros and poles of L + 0 (k) (the set Z + c ∪ P + c ) are all located in Imk < 0, and the zeros and poles of L - 0 (k) (the set Z - c ∪ P - c ) are all in Imk > 0. Similarly, the functions l ± (k) can be expressed in terms of the nonzero real roots of the corresponding characteristic equations belonging to the sets Z ± r and P ± r in Eq. ( 54). These roots are placed into the "+" sets (which contribute to the solution at η > 0) if the associated group velocities ω (k) exceed the phase velocity V , and into the "-" sets (contributing to the solution at η < 0) if ω (k) < V . This ensures that the solution satisfies the radiation condition (36), and the radiated waves carry energy away from the front. Recalling the structure of the real roots discussed in Sec. IV A, we observe that for subkinks (V < 1) this implies that the roots ±z 2 j , j = 1, . . . , l, of M + (k) in Z + r and the roots ±p 2 j , j = 1, . . . , m, of M -(k) in P + r contribute to l + (k), while the remaining roots ±z 2 j-1 , j = 1, . . . , l + 1, of M + (k) in Z - r and ±p 2 j-1 , j = 1, . . . , m + 1, of M -(k) in P - r contribute to l -(k). We thus obtain

l + (k) = l j=1 1 + (0-ik) 2 z 2 2 j m j=1 1 + (0-ik) 2 p 2 2 j , l -(k) = l+1 j=1 1 + (0+ik) 2 z 2 2 j-1 m+1 j=1 1 + (0+ik) 2 p 2 2 j-1 (63)
for subkinks. When l = 0 or m = 0, the corresponding products equal unity. Here we combined symmetric pairs ±r of real roots using

1 - k r ± i0 1 - k -r ± i0 = 1 + (0 ± ik) 2 r 2 ,
where the notation r ± i0 underscores the fact that the real roots are effectively shifted into the half-planes Imk ≷ 0. In particular, the zeros and poles of l + (k) (the set Z + r ∪ P + r ) are moved into Imk < 0, while the zeros and poles of l -(k) (the set Z - r ∪ P - r ) are shifted into Imk > 0. In the case of shocks (1 < V < γ ), the sets Z ± c are empty, and we have

l + (k) = 1 m j=1 1 + (0-ik) 2 p 2 2 j , l -(k) = 1 m+1 j=1 1 + (0+ik) 2 p 2 2 j-1 . ( 64
)
Finally, in the superkink regime, both characteristic functions have no nonzero real roots, and thus we have

l ± (k) = 1. ( 65
)
We now consider the asymptotic behavior of the functions L ± (k). Note first that Eqs. ( 61)- (65) imply that

L ± (k) ∼ 1 -V 2 γ 2 -V 2 , k → ±i0, (66)
where we take the principal branch of the square root, which becomes purely imaginary when 1 < V < γ . As shown in Appendix, the asymptotic behavior at infinity is given by

L ± (k) ∼ R ∓1 , k → ±i∞, V < 1 or V > γ , L ± (k) ∼ R ∓1 k ±1 , k → ±i∞, 1 < V < γ ( 67
)
where R is given by

R = + z - p - z + p = l j=1 z 2 j m+1 j=1 p 2 j-1 l+1 j=1 z 2 j-1 m j=1 p 2 j , V < 1 ( 68
)
for subkinks,

R = - p + p = m+1 j=1 p 2 j-1 m j=1 p 2 j , 1 < V < γ ( 69
)
for shocks, while for superkinks the absence of radiation implies

R = 1, V > γ .
Following the standard Wiener-Hopf procedure [START_REF] Noble | Methods Based on the Wiener-Hopf Technique for the Solution of Partial Differential Equations[END_REF], we perform the analytic continuation of both sides of Eq. ( 60) to the entire complex plane and apply the Liouville theorem. Noting that the asymptotic estimates in Eq. ( 67) imply that both sides of Eq. ( 60) can be continued to a function that is at most linear in k, we obtain

L + (k) -ε * -ik χ + (k) + ε + 0 -ik = 1 L -(k) -ε * + ik χ -(k) + ε - 0 + ik = ψ 0 + ψ 1 k. ( 70
)
Here the constants ψ 0 and ψ 1 depend on the velocity regime due to the different asymptotic behavior in Eq. ( 67) for kinks and shocks. Taking the limit k → ±i0 in Eq. ( 70) and using the asymptotics Eq. ( 66), we obtain

1 -V 2 γ 2 -V 2 (ε + -ε * ) = γ 2 -V 2 1 -V 2 (ε --ε * ) = ψ 0 . ( 71
)
These relations hold for all velocities. Recalling Eq. ( 52), one can see that the first equality in Eq. ( 71) implies that the RH condition (29) automatically holds for ε ± .

Observe now that by Eq. ( 67), both sides of the first equality in Eq. ( 70) are constant at infinity when either V < 1 or V > γ . Therefore, we must set ψ 1 = 0 in these velocity ranges. For subkinks (V < 1) and superkinks (V > γ ), taking the limits of the two sides of the first equality in Eq. ( 70) as k → i∞ and k → -i∞, respectively, equating them to ψ 0 , and applying the consistency condition (26), which implies lim s→∞ [sε ± (±is)] = ε(0±) = ε c , then yields

ψ 0 = ε c -ε * R , ( 72
)
where we recall Eq. ( 57). Here R is defined in Eq. ( 68) for subkinks and R = 1 for superkinks. Equations ( 71) and [START_REF] Berman | [END_REF] then imply that in these regimes, the limiting states ε ± are fully determined by the velocity V via

ε ± = ε * + ε c -ε * R 1 -V 2 γ 2 -V 2 ∓1/2 . ( 73
)
Shocks (1 < V < γ ) correspond to the generic case when both constants ψ 0 and ψ 1 in Eq. ( 70) are nonzero. In this case, the zero-limit equation [START_REF] Gallavotti | The Fermi-Pasta-Ulam Problem: A Status Report[END_REF], which still holds, and the limits k → ±i∞ yield

ψ 0 = i γ 2 -V 2 V 2 -1 (ε --ε * ), ψ 1 = ε c -ε * R , ( 74
)
where R is defined in Eq. ( 69). Note, however, that although, as noted above, the RH condition ( 29) is automatically satisfied for all three types of fronts, in the case of shocks the limiting states ε ± are not uniquely determined by V , i.e., there is no condition that is equivalent to Eq. ( 73) that we have for subkinks and superkinks. Therefore, in the case of shocks one of the limiting strains remains a free parameter, which agrees with the conclusions we reached while considering the problem in both continuum and QC frameworks.

The solutions of the two equations in Eq. ( 70) thus take the form

χ ± (k) = ε * -ε ± 0 ∓ ik + ψ 0 + ψ 1 k 0 ∓ ik [L ± (k)] ∓1 . ( 75
)
Here ψ 0 is given by Eq. ( 72) and ψ 1 = 0 in the case of both kinks and subkinks. Instead, in the case of shocks ψ 0 and ψ 1 are given by Eq. ( 74). This yields the strains in the physical space given by

ε(η) = ε ± + 1 2π ∞ -∞ χ ± (k)e -ikη dk, η ≷ 0, (76)
where the integrals are computed by closing the contour of integration in Imk ≶ 0 for η ≷ 0 and applying the residue theorem. Here we recall that all real zeros and singularities have been effectively shifted off the real axis into the corresponding half-planes. As in the QC case, the solution can then be expressed in the general form (45). Recall that this form includes localized [ ± (η)] and radiative [ ± (η)] components.

The localized components ± (η) are given by exponentially decaying functions arranged in the infinite sums

+ (η) = z∈Z + c ω 2 -(z) -(zV ) 2 2z 2 V [ω + (z) -V ] L -(z)(ψ 0 + ψ 1 z)e -izη , -(η) = p∈P - c ω 2 + (p) -(pV ) 2 2p 2 V [ω -(p) -V ] (ψ 0 + ψ 1 p) L + (p) e -ipη . ( 77
)
The summation is over the sets of complex roots P - c [the poles of L -(k) in Imk > 0] and Z + c [the poles of 1/L + (k) in Imk < 0] defined in (54). To compute the residues, we used Eq. ( 58) and the identities 1/L + (k) = L -(k)/L(k) and L -(k) = L(k)/L + (k) that follow from Eq. [START_REF] Kang | [END_REF].

The radiative components ± (η) in Eq. ( 45) describe the lattice waves taking the energy from the moving front to infinity. For subkinks (V < 1), we have

-(η) = 2 m+1 j=1 α - j cos (p 2 j-1 η + β - j ), + (η) = 2 l j=1 α + j cos (z 2 j η + β + j ), (78)
where the second sum is zero when l = 0. For shocks (1 < V < γ ), there is no radiation ahead of the front, so + (η) ≡ 0, while -has the same form as above. The real coefficients α ± j and β ± j can be obtained from the polar representation

α + j e -iβ + j = L -(z 2 j )[ω 2 -(z 2 j ) -(z 2 j V ) 2 ] 2z 2 2 j V [ω + (z 2 j ) -V ] (ψ 0 + ψ 1 z 2 j ), α - j e -iβ - j = - [ω 2 + (p 2 j-1 ) -(p 2 j-1 V ) 2 ] 2p 2 2 j-1 V [V -ω -(p 2 j-1 )]L + (p 2 j-1 ) × (ψ 0 + ψ 1 p 2 j-1 )
with the corresponding values of ψ 0 and ψ 1 . Only the second equation is relevant for shocks since + = 0 in that case.

Here we used Eqs. (A4) and (A5) obtained in Appendix. Finally, for superkinks (V > γ ) there is no radiation either ahead or behind the propagating front, and so in this case

-(η) = + (η) ≡ 0.
In addition to strains, we can also explicitly compute the particle velocities v(η). To this end, we need to solve the equation v(η + 1)v(η) = -V ε (η), where ε(η) is given by Eqs. ( 45), (77), and [START_REF] Truskinovsky | Shock Induced Transitions and Phase Structures in General Media[END_REF]. Using Fourier transform, we obtain

v(η) = v ± + ± (η) + ϒ ± (η), η ≷ 1/2,
where v +v -= -V (ε + -ε -) coincides with the first RH condition in Eq. ( 5) for the continuum problem, and since one of v ± is arbitrary by Galilean invariance, we may set v ± = -V ε ± . Here we can explicitly identify the exponentially decaying terms

ϒ + = - z∈Z + c ω -(z) -(zV ) 2 4z sin z 2 [ω + (z) -V ] L -(z) × (ψ 0 + ψ 1 z)e -iz(η-1/2) , ϒ -= - p∈P - c ω + (p) -(pV ) 2 4p sin p 2 [ω -(p) -V ] (ψ 0 + ψ 1 p) L + (p) e -ip(η-1/2)
and the oscillatory terms ± (η) describing radiation. For subkinks (V < 1), we have

+ (η) = - l j=1 α + j z 2 j V sin z 2 j 2 cos (z 2 j (η -1/2) + β + j ), -(η) = - m+1 j=1 α - j p 2 j-1 V sin p 2 j-1 2 cos (p 2 j-1 (η -1/2) + β - j ), (79)
where the second sum is zero when l = 0. For shocks (1 < V < γ ), the function -(η) has the same form, while + (η) ≡ 0. For superkinks, -(η) = + (η) ≡ 0.

D. Dissipation rate

The knowledge of the exact solution of the discrete problem gives us access to the energy (phonon) radiation from the moving fronts to infinity. As we have already mentioned, since the radiated energy is lost by the front, the associated rate of the energy transport to infinity by lattice waves can be interpreted as the rate of dissipation.

Following the procedure we used for the QC model, we again consider the cumulative energy fluxes G + and G - emitted ahead and behind the front. Recalling Eq. ( 47), we find that dissipation rates R ± = G ± V on both sides are zero for superkinks, which involve no phonon radiation, and thus G + = G -= 0 in this case. For subkinks (V < 1) we obtain

R + = l j=1 E + (z 2 j ) [ω + (z 2 j ) -V ], R -= m+1 j=1 E -(p 2 j-1 ) [V -ω -(p 2 j-1 )],
where R + = 0 when l = 0, and

E + (z 2 j ) = v 2 j /2 + ε 2 j /2 and E -(p 2 j-1 ) = v 2 j /2 + γ 2 ε 2 j /2
are energy densities carried by individual lattice waves with (real and positive) wave numbers z 2 j ∈ Z + r and p 2 j-1 ∈ P - r , respectively, and the averaging is over the corresponding time periods. Using the expressions for strains ε j in Eq. ( 78) and particle velocities v j in Eq. ( 79) of the emitted waves with the corresponding wave numbers, we obtain

G + = 2 l j=1 (α + j ) 2 ω 2 + (z 2 j ) ω + (z 2 j ) V -1 , G -= 2γ 2 m+1 j=1 (α - j ) 2 ω 2 -(p 2 j-1 ) 1 - ω -(p 2 j-1 ) V , (80)
where G + = 0 when l = 0. For shocks (1 < V < γ ), G - has the same form, and G + = 0. This yields explicit expressions for the driving force G = G + + G -in different velocity regimes. Alternatively, we can compute the driving force from the macroscopic area-difference formula (9) (with E 1 = 1 and E 2 = γ 2 in the dimensionless formulation). Using Eq. ( 73) for the kink regimes, Eq. ( 29) for shocks, and recalling Eq. ( 52), we obtain

G = ⎧ ⎪ ⎨ ⎪ ⎩ γ 2 -1 2 1 -1 R 2 (ε c -ε * ) 2 , V < 1, γ 2 -1 2 [(ε c -ε * ) 2 + V 2 -1 γ 2 -V 2 (ε + -ε * ) 2 ], V ∈ (1, γ ), 0, V > γ .
For subkinks and superkinks, this yields the kinetic relations G = G(V ) [recall that R depends on V via Eq. ( 68) in the subkink regime], which complement the classical RH conditions, while for shocks the driving force remains dependent on the choice of ε + , which, as we recall, is a free parameter in this case. We have verified that these "macroscopic" expressions for G are equivalent to the ones obtained by computing the energy fluxes directly.

E. Admissibility

As in the case of the QC approximation, one still needs to check which of the obtained solutions are admissible, i.e., satisfy Eq. ( 27). In Fig. 11 we show the admissibility diagram for the discrete problem, which is a direct analog of the similar diagram for the QC model presented in Fig. 6. As in that case, admissible subkink and shock solutions in the discrete problem feature a single radiation mode propagating behind the front, where the wave number r is a positive root of the characteristic equation ω -(r) = V r, while + ≡ 0. In the superkink case, ± ≡ 0. In the case of shocks, one of the limiting states remains a free parameter, which agrees with both continuum and QC approximations. One can see that for V < 1, sufficiently fast subkinks are admissible. For V > γ , all superkinks satisfy the assumed inequalities. In the interval 1 < V < γ , the TW solutions describing shock waves are admissible inside the pink domain. In the blue domain, such TW solutions are not admissible and are replaced by the DSWs, as we will discuss in the next subsection.

We conclude that the main features of the QC regime diagram Fig. 6 are preserved in the full discrete model. Thus, both types of kinks, represented in Fig. 11 by one-dimensional manifolds, are admissible (for sufficiently large V in the case of subkinks). Shocks are again not defined uniquely for a given σ and are admissible for sufficiently large values of σ . The two diagrams differ significantly only at small V < 1, where the QC model, as expected, does not capture the complex resonant behavior of the (typically inadmissible) slow discrete subkinks.

Our comparison suggests that outside the regimes of particularly slow subkinks, all three types of transition fronts are adequately described by only a few roots of the characteristic equation capturing long (but not infinitely long) lattice waves. This implies that carefully designed QC theories with only a few parameters (describing the crucial mesoscopic scales) can be successful in capturing such a fundamental nonlinear dynamic effect as radiative friction. It also points to the paramount importance of the QC reproduction of the relevant mesoscopic time scales, in addition to the more conventional task of modeling the internal length scales. In other words, the task of the adequate dispersive approximation of the kinetic energy may be at least as crucial as the task of the satisfactory representation of the nonlocal elastic energy.

F. Numerical simulations

To test the stability of the obtained analytical solutions, we conducted a series of numerical simulations, in which, starting with Riemann initial data, we traced the emergence of the nonlinear transition fronts propagating at constant velocity. More specifically, we solved numerically the system (2) (rescaled so that ρ = 1 and h = 1) with N = 1000 springs and discontinuous initial conditions of the form

ε n (0) = ε l , n < 500, 0, n 500, dε n dt (0) = 0
and free boundary conditions. We used the Dormand-Prince algorithm (ode45 in MATLAB), and the duration of simulations was such that the boundaries did not affect the front dynamics.

In each simulation, we varied ε l and σ while keeping all other parameters fixed. As in the case of the QC model, we identified four generic types of traveling fronts which all emerged and stabilized by the numerical time t = 500. The results of the simulations are summarized in Fig. 12. They confirm the possibility of stable propagation of all three types of transition waves. Similar to the QC model, the transition fronts are accompanied by linear dispersive waves. In particular, in all cases such a wave appears behind the front and moves away from it with velocity -γ . In the case of a subkink shown in Fig. 12(a), there is also a linear dispersive wave propagating ahead of the transition front with velocity 1. In the superkink case [Fig. 12(d)], there are two linear dispersive waves moving behind the front with velocities γ and -γ .

Our results suggest stability of all three regimessubkinks, shocks, and superkinks-inside the corresponding admissible domains of the (V, σ ) plane. Recall that subkinks are admissible when V < 1 is sufficiently large. An example of a subkink propagation is shown in Fig. 12(a). We found that superkinks can only appear when V > 1 and σ < -ε c (γ 2 -1) < 0. An example is shown in Fig. 12(d). Recall also that shocks are only admissible when 1 < V < γ and σ is above a certain threshold, as shown in Fig. 11. An example of an admissible shock propagation is shown in Fig. 12(b). Inside the domain of inadmissible shocks, we expectedly do not find steady transition fronts but find instead the spreading transition profiles of DSW type [Fig. 12(c)], similar to the corresponding prediction of the QC model. We reiterate that the DSWs are mentioned here only for completeness. A detailed study of such nonsteady regimes is outside the scope of this paper, in part because these solutions are well documented in the literature. They appear here naturally as stable replacements for the inadmissible traveling waves.

V. APPLICATIONS IN METAMATERIAL DESIGN

The importance of metamaterials is due to their ability to exploit post-instability structural responses. Effectively, metamaterials utilize internal changes in the subelements, which imitate molecular phase transitions at supermolecular scales. The success of the metamaterial paradigm is due to the fact that artificial "metamolecules" with desired properties can be manufactured at the relevant scales.

The localized transition fronts studied in this paper can be viewed as elementary bites of mechanical information that can be generated, delivered, and erased in periodic lattice metamaterials. Due to the presence of stress-sensitive repeating structural units, such metamaterials can manipulate mechanical information using advantageously the dispersion of elastic waves. By carefully tailoring relationships between characteristic dimensions, one can design metamaterials combining the effects of strong dispersion with various forms of energy nonconvexity. One of the main challenges in the design of metamaterial structures is to ensure that the switching takes place at a predefined level of stress and that the particular switching waves are generated when the task is, for instance, to enhance actuation or perform energy harvesting.

In view of these and other potential applications, the prototypical FPU model studied in this work can serve as a proof of concept showing the broad variety of the functionally distinct switching regimes, which can be controlled by the deliberate parameter tuning. Even though the actual 3D metamaterials with the desired properties would still have to be designed, the results obtained in this paper already provide a specific guidance regarding, for instance, which metamaterial should be chosen to ensure a supersonic, dissipation-free communication of mechanical information, as opposed to a design favoring subsonic switching, which ensures a heavily dissipative response.

VI. CONCLUSIONS

The goal of this paper was to reveal the interrelations between structurally different steadily moving transition (switching) fronts in the classical FPU model. Our main result is the demonstration that this non-integrable Hamiltonian model supports three types of such fronts that can be classified as subsonic (subkinks), intersonic (shocks), and supersonic (superkinks).

To obtain analytical results, we limited our analysis to piecewise linear elastic responses. In this case, exact solutions of the discrete model for each class of fronts can be presented in the form of infinite series. Within this setting, we have shown that the proposed classification is exhaustive. The common framework considered in this work allows us to describe all three types of switching waves in a unified way and associate them with particular classes of elastic responses.

While the constructed explicit solutions of the discrete problem are sufficient to corroborate these qualitative claims, the origin of the difference between the three types of fronts remains relatively opaque in the FPU setting dealing with an infinite system of nonlinear ordinary differential equations. To achieve conceptual transparency, we constructed a QC approximation of the FPU problem. An excellent agreement with the behavior of the discrete model was obtained using a long-wave (infrared) approximation utilizing only two internal scales. We stress that the successful coarse-grained theory relies on the approximation of kinetic energy, in contrast to more conventional asymptotic approaches such as the KdV model and its higher-order analogs.

A detailed comparison of the exact solutions for the QC theory and the discrete problem showed that the chosen approximation adequately describes the complex interrelation between all three types of the transition fronts. This means that the whole complexity of the dispersive structure of the FPU model was not really necessary for the description of the main features of these special solutions. In other words, the dispersive properties of all three different classes of fronts can be satisfactorily captured using a simple QC model.

Our analysis also reveals that the obtained macroscopically dissipative subkink and shock front profiles cannot be adequately described by the continuum nonlinear wave equation, as may be suggested by a naive homogenization. Instead, they should be interpreted as microscopic descriptions of Whitham shocks connecting oscillatory and constant states [136,137]. Such generalized (dispersive) shocks usually correspond to heteroclinic traveling waves of a dispersive model connecting standard critical points with periodic orbits. To capture such connections in a PDE format, we had to use a higher-order QC model.

To fully understand the different structure of the three types of transition fronts, we have drawn upon a broad variety of physical and mathematical considerations, including characteristics, barriers, topological transitions, undercompressive nature, critical manifolds, and kinetic relations, which all point to the existence of exactly three universality classes of transition fronts. In this sense, the obtained perspective can be viewed as unifying not only for the description of switching waves but also for different analytical approaches to the analysis of nonlinear dispersive systems.

Several important issues have been naturally left for future studies. The traveling wave description of the switching waves is clearly incomplete when it comes to transient effects such as interaction with obstacles and multiple collisions. The approach to such problems proposed for special cases in [138] can also be generalized and applied in our more general framework. The present work does not address thermal effects, which may become relevant for metamaterial with a submicron scale mimicking cytoskeleton or extracellular environment. For these purposes, the approach proposed in [139] can be generalized here as well. Another issue that we have not addressed in this work concerns different modes of manipulation and control of transition fronts from a distance using DC-and AC-type dynamic loading, which is of particular interest for metamaterial applications. The successful use of such control was recently demonstrated for semilinear discrete systems in [124].

FIG. 4 .

 4 FIG. 4. Different behavior of the dissipation functionG(ε, ε + ): (a) subkinks, V < c 1 ; (b) shocks, c 1 < V < c 2 ; (c) superkinks, V > c 2 .

FIG. 5 .

 5 FIG. 5. The characteristic roots in the quasicontinuum model when (a) V < 1, (b) 1 < V < γ , and (c) V > γ . Due to symmetry, only the roots with Imk 0 and Rek 0 are shown. Insets show the dispersion relations and real roots as intersections with the line V k.

FIG. 7 .

 7 FIG. 7. Different regimes of front propagation in the QC model with the parameters γ 2 = 1.5 and ε c = 1 at t = 50: (a) subkink (ε l = 6, σ = 2.5); (b) conventional shock (ε l = 10, σ = 2.5); (c) dispersive shock (ε l = 10, σ = 0); (d) superkink (ε l = 6, σ = -1.5).

FIG. 8 .

 8 FIG.8. Magnitudes of (a) the real roots k = ±r and (b) the imaginary roots k = ±iq in the discrete (solid curves) and QC (dashed curves) forω 2 + (k)k 2 V 2 = 0 (blue curves) and ω 2 -(k)k 2 V 2 =0 (red curves). Black solid lines mark the sonic limits separating the velocity domains of different transition fronts. Complex roots with nonzero real and imaginary parts bifurcate at the velocities V * and V * * marked by dash-dotted lines in the QC model and from the nonsonic maxima of the real root curve in (a) in the discrete model. Here γ = 2.

FIG. 9 .

 9 FIG. 9. Distribution of the roots of M ± (k) in Eq. (53) for the discrete model when (a) V < 1, (b) 1 < V < γ , and (c) V > γ . Due to symmetry, only the roots with Rek 0 and Imk 0 are shown. Insets show the dispersion relations and real roots as intersections with the line V k.

FIG. 11 .

 11 FIG. 11. Admissibility sets of solutions of the discrete problem. In the blue region we observe ε(η) ε c when η < ε c , and the dashed lower boundary of the region marks the threshold ε -= ε c . The insets show examples of the strains ε(η). Here γ 2 = 1.5, ε c = 1, and we set ε + = 0.

FIG. 12 .

 12 FIG. 12. Different regimes of front propagation in FPU chain under Riemann-type initial conditions with different left strain ε l and σ : (a) subkink (ε l = 5, σ = 2.5); (b) conventional shock (ε l = 25, σ = 2.5); (c) dispersive shock (ε l = 25, σ = 0); (d) superkink (ε l = 5, σ = -1.5). Here γ 2 = 1.5, ε c = 1, and t = 300.

ACKNOWLEDGMENTS

The authors thank G. Mishuris for helpful discussions. The work of A.V. was supported by the NSF Grant No. DMS-1808956. L.T. and N.G. acknowledge the support of the French Agence Nationale de la Recherche under Grant No. ANR-17-CE08-0047-02.

APPENDIX: SOME ASYMPTOTIC RESULTS

To obtain the asymptotic behavior at infinity, we follow Ref. [START_REF] Slepyan | [END_REF] and observe that for subkinks (V < 1) we have

where

and

where

Here the expressions under the logarithms in the Cauchy-type factorization integrals are set up in such a way that they tend to 1 as k → ±i∞, while the logarithms remain real along the entire integration path [START_REF] Slepyan | [END_REF]. These asymptotic expressions imply that in the subkink regime

where R is given by Eq. ( 68). For shocks (1 < V < γ ) Eq. (A2) still holds, but due to the absence of nonzero real roots of M + (k) in this regime, Eq. (A1) is replaced by [START_REF] Slepyan | [END_REF]

so that

which together with

where R is given by Eq. ( 69). Finally, for superkinks (V > γ ), both characteristic functions have no nonzero real roots, and thus H ± z ∼ iV as in Eq. (A3) and H ± p ∼ iV in the limit k → ±i∞. Together with (65), this implies L ± ∼ 1 as k → ±i∞ in this velocity regime. Combining these results, we obtain Eq. [START_REF] Novelino | Proc. Natl. Acad. Sci. USA[END_REF].

Recalling Eq. ( 58) and Eq. ( 59), one can also show that near the real singularities,

and

, k → p 2 j-1 , (A5) with similar expressions for the negative real singular points.