Analysis of nonlinear elastic aspects of precursor attenuation in shock-compressed metallic crystals

Unsteady characteristics of shock waves in metals, for example elastic precursor decay, have often eluded a complete model description. Historic continuum elastic-plastic theories tend to require excessive initial dislocation densities in order to match experimental observations. Studies incorporating superposition of linear elastodynamic solutions for dislocations, either in analytical solutions or in discrete numerical simulations, omit nonlinear elastic effects and only consider effects of defects immediately at the elastic shock front. Prior analytical treatments of hydrodynamic attenuation consider nonlinearity manifesting as effects of stress gradients or particle velocity gradients immediately behind the front. The present analysis seeks to augment the predictive precursor decay equation from linear elastodynamics to account for elastic nonlinearity and dislocation nucleation in the wake of the precursor shock. A complete solution for the precursor magnitude at a material point is shown to require consideration of the prior history of dislocation generation and the entire flow field behind the elastic shock up to that material point. However, introduction of reasonable simplifying assumptions and basic models for dislocation generation and glide resistance enable derivation of a mathematically tractable relation for precursor decay. This relation is a nonlinear first-order ordinary differential equation that, in non-dimensional form, contains only one scalar parameter controlling the rate of dislocation generation behind the shock. Model predictions that include nonlinear effects are shown to provide a better match to experimental data for three metals. Effects of nonlinearity are shown to emerge early, but not immediately, in the shock attenuation process, and these effects increase in prominence with increasing impact stress.


Introduction
A standard means of quantifying the mechanical response of metallic single and polycrystals subjected to large compressive stress at high rates of deformation invokes planar shock loading, e.g., via explosively generated waves or flyer-plate impact. Such studies have served as the ubiquitous shock physics experiments for quantifying equations-of-state and dynamic strengths of solids for over six decades [1,2]. Besides metals, the focus of the present article, plate impact methods are also widely used to study the high-pressure, high-rate response of polymers, ceramics, and glass [3][4][5], each of which demonstrates particular behaviors with regard to elastic nonlinearity, structural changes, and possible inelasticity.
In metals, with increasing impact stress, shocks can usually be labeled as elastic, elastic-plastic, or overdriven. Elastic shocks in ductile solids occur only for small to moderate impact stress, by definition stresses not exceeding the Hugoniot Elastic Limit (HEL). Elastic-plastic shocks occur for moderate stresses (above the HEL), wherein local stress components (e.g., resolved shear stresses) are sufficient to enable dislocation motion and generation, but at stresses lower than those corresponding to overdriven shocks. A shock is said to be overdriven if its velocity exceeds that of any elastic precursor, which in turn leads an overtaking of the precursor by the plastic wave and formation of a single-wave structure. A shock that is overdriven is often labeled as a strong shock, while a plastic shock whose velocity is lower than that of its precursor is labeled as a weak shock [6]. An elastic shock in a metal is necessarily also weak since its stress is even lower than that of a weak plastic shock in the same material. A shock is said to be steady if its profile or structure does not change with time, whereas a shock is unsteady if its structure evolves. Viewed together, the elastic and plastic waveforms comprising a weak shock are unsteady since the distance between the precursor front and the plastic wave increases continuously with time, tending to spread the overall structure.
In planar elastic-plastic shock loading of materials of adequate symmetry (i.e., pure mode directions), a pure two-wave structure exists, with a longitudinal plastic shock trailing a faster moving longitudinal elastic precursor. A pure mode direction [7,8] enables propagation of one purely longitudinal mode, with its polarization vector parallel to the propagation direction, and two transverse modes, with polarization normal to the propagation direction. For planar impact loading along other directions different from pure mode directions, quasi-longitudinal and quasi-transverse elastic and plastic waves appear.
In the elastic-plastic regime, contributions of shear strength and mean stress may both be significant, depending on the material. The traditional overdriven regime [9] is attained when, as a result of stiffening of the material behind the front due to nonlinear elastic and entropic effects, the plastic wave outpaces any elastic precursor. Relative contributions of dislocation generation/multiplication to movement of pre-existing dislocations are thought to become more important as impact stresses increase, and homogeneous nucleation of dislocations becomes more important as shocks increase in strength from weak to overdriven [10][11][12]. On the other hand, the relative magnitude of deviatoric stress responsible for dislocation activity to mean stress (i.e., hydrostatic pressure) decreases in the plastic region with increasing shock stress, leading to adequacy of the hydrodynamic approximation (i.e., disregard of shear stress) for analysis of many features of sufficiently strong shocks.
Precursor decay can be thought of as a decrease in a transient HEL of the material from the initial/impact stress to an often steady minimum value. The HEL referred to in the context of shock regimes discussed above is this final steady value. In the weak shock regime (elastic-plastic shocks), a number of theoretical and numerical studies [13][14][15][16] have sought to describe the transient decay of elastic precursors observed experimentally in metals. If dislocation density is assumed constant, use of Orowanʼs equation in continuum elastic-plastic analysis requires dislocation densities several orders of magnitude larger than those measured in the bulk crystal [13]. Incorporation of larger transient dislocation density via multiplication laws allows for closer agreement of numerical results with experimental precursor wave profiles [14]. Importance of consideration of elastic nonlinearity, especially effects of stress on effective shear stiffness, for correctly predicting precursor behavior have been noted [14,15]. The importance of nonlinear elastic effects on weak shock wave structure was first quantified in experiments on rock materials [17], and further examined in the context of ductile nonmetallic crystals (specifically, lithium fluoride [18][19][20][21]).
Analysis incorporating linear dislocation elastodynamics [16], using superposed isotropic solutions for effects of motion of existing edge dislocations on the particle velocity field [22], demonstrated that the aforementioned continuum analyses may be inaccurate by a factor of ≈2 because they do not take into account differences between dislocations moving towards versus away from the wave front. In that work [16], which considered straight edge dislocations, it was speculated that large dislocation densities or fluxes at free surfaces may be the main source of discrepancy between observations and model representations of precursor decay. Linear elastodynamic analysis considering dislocation loops [23] produced similar conclusions. Linear elastodynamics is defined in this context as the treatment of transient waves generated by and emanated from the nucleation and acceleration of dislocations, respectively, in a linear elastic medium.
Recent work [24][25][26][27][28] has extended the elastodynamic analyses of Markenscoff and Clifton [16,22] to consider injection and subsequent non-uniform motion of edge dislocations. This numerical work improves existing (quasi-static) discrete dislocation plasticity simulations wherein spurious dislocations can appear ahead of a shock wave front because stress fields are transmitted instantaneously, violating the causality principle. Application of the method to a study of precursor decay in aluminum [26] demonstrated reasonable agreement of stress relaxation in shock fronts at various strain rates, considering only effects of nucleation of new dislocations (i.e., no pre-existing defects) in a linear elastic medium. It was noted how cumulative effects of shielding dislocations exceed those of anti-shielding dislocations, and defect morphologies reminiscent of Smith-Hornbogen interfaces [29][30][31][32] were predicted by the simulations. Significant extrapolation of numerical results to compare with longer-term decay rates measured experimentally was required, however. Furthermore, all of these aforementioned discrete dislocation simulations, as well as the prior analytical treatments upon which they are based, only consider effects of defect nucleation and/or motion immediately at the elastic shock front since linear elasticity is invoked for dislocations' stress fields and corresponding wave speeds. Any defect activity occurring at a finite distance behind the front does not affect the front since the longitudinal linear elastic wave speed is lower than the shock velocity, or at most equal to it in a linear elastic approximation of the shock itself.
The present work seeks an improved understanding of precursor decay via consideration of elastic coefficients and sound velocities for the medium behind the elastic wave front that have been updated to account for stiffening due to the elastic compressive wave [33,34]. First, in section 2, key equations describing precursor decay used in prior analyses are reviewed. Then in section 3, shock wave and release velocities are computed via anisotropic Eulerian elasticity [11,[35][36][37] over stresses ranging from the (final minimum) HEL to that of a plastic wave at the overdriven limit wherein shock stress attains its imposed maximum for a two-wave structure. Single crystals of three metals of vastly differing impedance and crystal structure-copper (Cu), magnesium (Mg), and tungsten (W)-are compared for shocks in pure mode directions. Next, in section 4, consideration of transient effects of dislocation generation and motion at points behind the precursor shock on its magnitude leads to derivation of a new expression for precursor decay accounting for nonlinearity. Physically reasonable simplifying assumptions and dimensional analysis are used to reduce the expression to a nonlinear first-order ordinary differential equation (ODE). This equation contains a single scalar parameter controlling dislocation generation effects that dominate the nonlinear contribution to decay of the elastic wave front. The treatment of dislocation mechanics invoked in the derivation essentially superposes elastic fields from those dislocations generated over relevant time-and space-intervals that would affect the front. Such fields are interpreted in the context of small elastic deformations from a potentially finitely deformed state behind the shock. In section 5, solutions to this ODE are analyzed and compared to experimental data. Conclusions follow in section 6.
No attempt is made here to explicitly compute precursor decay in a fully nonlinear setting for dislocation dynamics. Nonlinear elastic dislocation solutions are available only for special geometries and static problems [38], and standard superposition required for consideration of effects of multiple dislocations cannot be used in an exact nonlinear setting.

Elastic precursor decay 2.1. Continuum plasticity analysis
Elastic precursor decay is first considered from the viewpoint of nonlinear elastic-plastic continuum theory [17,39]. Let P denote the axial stress, positive in sign under compressive strain corresponding to planar shock loading. The original analysis of Ahrens and Duvall [17] invokes the following constitutive equation for a stressrelaxing solid in the region of material immediately behind a planar elastic shock front propagating in the xdirection: Here (·) t d d denotes the material time derivative (i.e., ∂ (·)/∂t at fixed Lagrangian coordinate X), a is the Eulerian sound speed for isentropic longitudinal elastic compression waves, ρ denotes the mass density of the material at time t, and with G an elastic shear modulus and D P the isochoric Eulerian plastic velocity gradient, i.e., a plastic strain rate referred to spatial coordinates. The sign convention used here and in [17] is > D 0 xx P for compression in the xdirection. The usual continuum balance laws for conservation of momentum and mass are applied to the smoothly deforming region behind the shock, while the usual Rankine-Hugoniot jump conditions [2] are applied across the surface of singularity at the elastic shock front. The rate of change of P corresponding to precursor decay at the shock front is obtained by differentiation along the shock path: is the Eulerian elastic shock velocity, here equal to the Lagrangian shock velocity X t D D since the material ahead of the shock is quiescent and undeformed [39].
Upon use of (2.1), (2.2), (2.3) and the aforementioned balance laws and jump conditions, the following relationship is derived as in [17]: H H the Eulerian Hugoniot sound speed. Particle velocity in the shocked state is υ, and the initial mass density of the solid, ahead of the precursor shock, is ρ 0 . If the release isentrope is close to the Hugoniot, then a≈c H , and if the tangent and secant to the Hugoniot are similar (i.e., relatively weak shocks), then a≈U. Invoking these two approximations in [17], (2.4) reduces to The first term in the numerator on the right side of (2.5) accounts for changes in particle velocity behind the wave front, i.e., the so-called 'hydrodynamic attenuation' from transient effects including rarefactions that overtake the front from the rear. As shown in figure 1, if the stress gradient ∂P/∂X immediately behind the front is positive in sign, corresponding to particle deceleration ( u < t d d 0) then this term hastens precursor decay. Conversely, a negative stress gradient corresponding to particle acceleration immediately behind the precursor front would inhibit its decay. The former situation (∂P/∂X>0) seems much more prevalent in data records from experiments [17][18][19][20]40] and nonlinear continuum simulations [14,15]. Physically, the drop in stress behind the precursor spike (figure 1(a)) has been attributed in ductile crystals to stress relaxation caused by multiplication of mobile dislocations [14,19,40,41].
The second term in the numerator of (2.5) accounts for 'intrinsic attenuation' resulting from stress relaxation mechanisms such as dislocation activity or fractures occurring immediately at the front. Equation (2.4) with c H ≈a=ρ 0 c l /ρ can be rearranged to [39]  with tangent shear modulus G depending on the deformation or stress state immediately behind the elastic shock front. Strain rate component D P xx is equal in magnitude to the axial rate of plastic deformation gradient since the deformation itself at this location is elastic and uniaxial [39]. Experimental measurements of the decay of peak stress of the elastic shock with propagation time in conjunction with measurements of the rate of change of stress at a material particle behind the shock can be used with a modified form of (2.5) or (2.6) to calculate intrinsic attenuation function F directly from time-resolved plate impact data [18,19].
In the linear approximation, which omits hydrodynamic attenuation (i.e., = = U c constant l and = G constant), (2.6) reduces to Values of F calculated using either expression (2.6) or (2.7) can be compared with those resulting from Orowanʼs equation: Here ρ D is the density of dislocations at the shock front gliding with velocity υ D , b is the relevant Burgerʼs vector component, and α is a constant related to slip system geometry. This equation can be obtained directly by substituting . The linear form in (2.7) was applied to interpret precursor decay data from shocked polycrystals [42] and single crystals [13]; a nonlinear form analogous to (2.6) was applied to evaluate precursor decay data from single crystals [19][20][21]. Alternative derivations of equations identical to or very similar to (2.4)-(2.7) can be found in [14,15,18,20,21,39].

Discrete dislocation analysis
In lieu of use of a continuum plasticity constitutive equation such as (2.1) or Orowanʼs equation (2.8), alternative studies of precursor decay have considered explicit effects of stress relaxation due to transient stress fields of individual dislocations, albeit often in highly idealized geometric configurations. Invoking isotropic linear elastodynamics, the following precursor decay equation was derived in [16] by superposing changes in particle velocity fields induced by pre-existing mobile edge dislocations of constant density ρ D =ρ D0 : a q r u u q = -- with θ the orientation of the slip plane relative to the shock direction x. Equation (2.9) accounts for the influence of dislocations of opposite signs (via the algebraic signs of υ D and θ), with those moving towards the front accounting for a larger fraction of precursor decay than those moving away from the front. The maximum magnitude of the denominator of (2.9) is around 1 2 , leading to a maximum difference between (2.9) and that of (2.7) and (2.8) of about a factor of two. Given a constitutive model for the stress dependence of dislocation velocity υ D , values of initial dislocation density, b, and θ, and linear elastic properties G=G 0 , c l =c l0 , the differential equation (2.9) [or even more simply (2.7) with (2.8)] can be integrated numerically to predict precursor decay, i.e., stress P at the front versus time t or propagation distance X. The recent numerical study of [26] invokes dynamic discrete dislocation mechanics and analytical expressions from [16,22,24,43] to compute effects of injected and subsequently moving and interacting dislocations, with potentially complex morphologies, on precursor decay.
Both [16] and [26] omit the possible contributions of hydrodynamic attenuation and effects of dislocation nucleation, motion, and interactions occurring at a finite distance behind the precursor front since all stress waves are assumed to propagate at the linear elastic sound speeds of the undeformed medium. The longitudinal speed c l0 is equal to U in the linear approximation invoked in such studies, so trailing rarefaction waves never reach the shock front and thus cannot contribute to its decay. Subsequent analysis in the present paper removes such restrictions associated with equal longitudinal wave velocities, thus permitting consideration of effects of rarefaction sources behind the shock front omitted in the linear setting.

Wave speed analysis
3.1. Nonlinear elastic precursor Anisotropic, finite strain Eulerian thermoelastic theory advanced in [35] is invoked here for analysis of elastic shocks in crystals. This theory has been demonstrated [36,37] to be more accurate, upon inclusion of elastic constants of the same order, than Lagrangian theory [8,34,44,45] for modeling the high pressure response of ductile metals under uniaxial shock loading, as well as under hydrostatic loading, whereby the Eulerian theory reduces to the pressure-volume equation-of-state of Murnaghan and Birch [46].
Consider an elastic shock of strength P (longitudinal Cauchy stress, positive in compression) traveling in the X=X 1 direction in a crystal at instantaneous Lagrangian velocity U. Particle velocity, internal energy density per unit reference volume, and compression ratio in the shocked state immediately behind the surface of singularity are denoted by υ, e, and J=V/V 0 =ρ 0 /ρ, respectively. Mass densities ahead of and immediately behind the front are ρ 0 and ρ. Material ahead of the shock is undeformed and stress-free. Denote the Eulerian uniaxial strain measure, negative in compression, by Let C 11 and C 111 be the usual second-and third-order Lagrangian elastic stiffness constants, Γ the axial Gr ü neisen parameter, and T 0 the ambient temperature ahead of the shock. The crystal internal energy potential under uniaxial strain is, to order three in strain and order one in entropy density η, The analytical solution to the Rankine-Hugoniot equations derived in [35] for stress, particle velocity, and shock velocity is Properties for three metals considered in the calculations are given in table 1. Elastic shock propagation is prescribed along certain pure mode directions: along a cube axis for Cu and W and an a-axis for Mg.
The following definitions of dimensionless shock stress are invoked to cover the domain of shocked states wherein a two-wave structure exists: As introduced in the first of (3.6),P is the dimensionless shock stress associated with true shock stress P in the elastic precursor wave, P OD is the imposed shock stress at the overdriven limit, and P H is the final (minimum) HEL stress after any precursor decay. The functionP sets the precursor stress to a value between the HEL at = P 0 to a maximum of = P 1 for an overdriven shock wherein the elastic and plastic shock wave speeds just become equal. The second of (3.6) defines a dimensionless stress for the precursor, with a value of = P 0 at the final HEL. Over the duration of a precursor decay event, the value ofP would decrease with time from = = -P P P P 1 H 0 0 , with P 0 the imposed/initial impact stress, to = P 0. Note that P H and P OD are treated as constants for a given material, with the latter determined via use of constitutive relations for the plastic wave to be presented in section 3.2. Definitions (3.6) are clarified in figure 2(a) for a hypothetical elastic-plastic solid with a rather large HEL.
Given an input volume ratio J for the elastic shock state, equations (3.1), (3.3), and (3.5) are sequentially invoked to obtain the strain D, shock stress P, and shock velocity U for this value of J. As shown in figure 2(b), the shock velocity always exceeds the linear elastic wave speed r = c C l0 1 1 0 , with the margin of increase becoming larger with increasing stress (i.e., with increasingP). Even at the minimum value of precursor stress at = P 0, U>c l0 for all three materials. The ratio U/c l0 attains a maximum on the order of 1.5 for Mg and W for precursor strengths at the overdriven limit. The maximum for Cu is significantly lower, on the order of 1.1 at = P 1 (figure 2(b)).

Plastic shock and material response along the principal Hugoniot
The plastic wave velocity U P is used to obtain the overdriven threshold stress P OD that limits the maximum imposed stress considered in calculations. The constitutive model used for the plastic shock at impact stresses exceeding P H is the standard linear shock velocity versus particle velocity relation [2,39]: This equation is only used in the present work to model the plastic wave, specifically to set the domain of possible precursor stresses to . It is not used to compute the behavior of the elastic precursor, which instead will be addressed later using a more detailed dislocation mechanics-based approach. Material constants are c 0 and s; polycrystalline values are implemented (table 1) since these are readily available and since the axial stress-strain response at high pressure is dominated by compressibility rather than anisotropy or shear.
Hugoniot curves are constructed as follows. The material response at stresses up to P H is computed incrementally using the nonlinear elastic theory of section 3.1, with the threshold particle velocity υ H , volume ratio J H , and elastic shock velocity U H recorded when P H is reached. The particle velocity is increased above this threshold incrementally, and (3.7) is invoked to obtain the corresponding Lagrangian plastic wave velocity U P . The Rankine-Hugoniot equations for conservation of mass and momentum are solved for the volume J and stress P in the plastic regime, where the state ahead of the plastic shock corresponds to the HEL: Incrementation continues until U P =U H . Subsequently, in the overdriven regime,  U U P , and (3.8) becomes J=1−υ/U and P=ρ 0 U υ. Hugoniot stresses normalized by initial isentropic bulk modulus B 0 are reported in figure 2(c).
The ratio of elastic shock velocity at the HEL, U H to the plastic wave speed U P is shown in figure 2(d), where hereP is computed via the first of (3.6) with P the stress in the plastic wave rather than the precursor. The ratio of elastic to plastic shock velocities decreases with increasing stress. Unsteady states associated with precursor decay (i.e., elastic shock stresses exceeding P H ) do not fall on the Hugoniot curves of figures 2(a) and (c) or enter ratios in figure 2(d).

Elastic wave speeds behind the shock front
Cartesian coordinates of material particles in (current, pre-stressed, unstressed) states. The current state corresponds to the motion of the body at time t. The unstressed state corresponds to a natural minimum energy configuration of the body, assumed homogeneous. The pre-stressed state corresponds to a condition of static equilibrium and homogeneous deformation, and will later serve as an idealization of the region just behind the precursor front at a fixed instant in time. The identical idealization, though unstated, was invoked in prior work [18,19,21] wherein third-order Lagrangian elasticity was used to compute c l . It is most reasonable for steady shocks wherein particle velocity behind the precursor is constant, since deformation gradients and equations of motion are invariant upon superposition of a constant background velocity field (i.e., rigid body translation) in the pre-stressed state. Deformations from the pre-stressed to current state are assumed isentropic.
Functional forms of coordinate changes among states are Define displacement from pre-stressed state to current state and particle velocity from pre-stressed to current state as Material time differentiation occurs equivalently at fixed X I or X I 0 since both sets of coordinates are static in time. Mass densities in the three states are (r r r , , 0 ), related as usual by Jacobian determinants of deformation gradients. Letting σ ij denote symmetric Cauchy stress components, in the absence of body forces the local balance of linear momentum with respect to the current state is To first order in displacement gradients ∂u i /∂X I , the linear momentum balance can be written with respect to the intermediate state as [33] d where, with s IJ the symmetric pre-stress referred to X K and e(D IJ , η) the internal energy per unit initial volume, Transformation formulae derived in [35] have been applied to convert the Lagrangian elastic coefficients used in [33,34] to the Eulerian representation. Coefficients B IJKL do not, in general, satisfy full Voigt symmetries, nor do they possess the initial material symmetries of isentropic second-order elastic constants

Longitudinal waves
Now consider uniaxial elastic deformation corresponding to a region immediately behind the precursor front, with respective volume ratio and strain J and D related via (3.1). In this case, noting that s = -P 11 is the shock stress from (3.3), the longitudinal wave propagation coefficient obtained from (3.13) and use of the third-order internal energy function in (3.2) is Here B 11 is the equivalent expression for B 1111 in Voigt notation. The longitudinal wave velocity referred to the initial, undeformed configuration is In the present context, terms in (3.13) become r r  and . The factor of 1/J transforms the wave velocity a relative to the compressed state at density ρ to c l relative to the natural state at density ρ 0 [18,21].
The ratio of longitudinal wave speed c l computed via (3.16) to shock velocity U computed via (3.5) is shown in figure 3(a) for precursor shock strengths ranging from the HEL at = P 0 to the overdriven threshold at = P 1, for Cu, Mg, and W single crystals. This ratio increases monotonically with increasing P orP, and even at = P 0 always exceeds unity. Specifically, at = « = P P P 0 H , c l /U=1.002 for Cu, c l /U=1.004 for Mg, and c l /U=1.010 for W. The ratio at higher stresses is significantly smaller for Cu than Mg and W. Shown in figure 3(b) is c l divided by the shock velocity at the HEL, U H , for back stresses less than or equal to the HEL stress P H . This situation corresponds to a scenario such as depicted in figure 1(a), wherein the region behind the precursor could conceivably be of lower stress than the precursor front featuring a spike. Recall from (3.6) that = -P P P 1 H . Here,P corresponds to the dimensionless back stress behind a precursor spike, where a value of zero would mean no spike or dip. Results shown in figure 3(b) demonstrate that the longitudinal wave speed c l would still exceed the precursor shock speed U H for  -P 0.5, i.e., for back stresses in excess of about P H /2.
Collective results in figures 3(a) and (b) suggest that longitudinal waves emanated from sources behind the precursor front should eventually reach the front so long as the stresses behind the front exceed about half the HEL stress. The time required for such sources to reach and interact with the front would depend on the time and position of their origin and any changes in the compressive state behind the front during their transit. Summarizing, c l Uc l0 , with equalities holding in the linear elastic limit.

Shear waves
In anisotropic media, shear waves travel at different velocities depending on their directions of polarization and propagation. However, shear wave velocities for directions of interest in the present work, in single crystals of cubic and hexagonal symmetries, can be computed explicitly for pre-stress conditions corresponding to finite uniaxial elastic strain using second-and third-order isentropic elastic constants.
First consider propagation, in the longitudinal X 1 -direction, of a shear wave polarized in a direction orthogonal to X 1 . For a cubic crystal, propagation is along  An incremental or tangent elastic shear modulus G, which can be parameterized versus shock stress = (¯) P P F 11 for uniaxial elastic shock loading, is introduced along with its value at the ambient state as 3 . 1 9 The ratio of shear wave velocity c s from (3.18) to the elastic shock velocity U of (3.5) is shown in figure 3(c). Shear waves clearly lag the shock over longitudinal stresses ranging from the HEL at = P 0 to the overdriven threshold at = P 1. Ratios of c s /U are largest for Mg at high stresses but always remain less than unity. Similar conclusions are obtained from computed ratios of the Lagrangian shear wave speed obtained from G of (3.19) to the elastic shock speed (not shown). These results suggest that shear waves emanated from sources behind the precursor front, in contrast to longitudinal waves, should neither reach nor explicitly interact with the front.

Nonlinear precursor decay
Direct application of existing relations such as (2.5) or (2.6) that account for nonlinear effects to predict precursor decay is challenging. Stress gradients or material accelerations behind the precursor front must be obtained from sophisticated numerical schemes [14,15] that invoke a constitutive relation between D xx P and stress or other state variables, to apply these equations.
The objective of the present section is derivation of a differential equation for precursor decay -accounting for both the intrinsic attenuation from pre-existing dislocations at the front and for rarefaction signals due to dislocations generated behind the front-that can be solved via basic numerical integration. The starting point for the derivation is (2.6), reproduced from equation (10.55) in [39], that can be rearranged as Each group of terms on the right side of (4.1) will be represented in terms of quantities obtained from fundamental dislocation mechanics. Specifically, the intrinsic term will be represented via the solution for mobile edge dislocations at the precursor front derived in [16], suitably modified to account for the stress dependence of the shear modulus via (3.19) and the nonlinear elastic contribution in the denominator with c l /U calculated via the method described in section 3.3.1. The contribution from trailing rarefaction sources will be attributed to dislocation generation, noting that the stress gradient is thought to result from dislocation multiplication behind the front [14,19,40,41]. Recall that the linear decay equation (2.7) used in [13] is obtained from (4.1) by setting c l =U, reiterating that the contribution from trailing rarefaction sources vanishes identically in the linear model.

Problem geometry
Dislocation generation and motion for shock loading is analyzed for plane strain geometry in configurations considered in [25,26], which account for three possible slip systems oriented at (0°,±θ) from the direction of shock propagation, as shown in figure 4(a). Coordinate frames (x, z) and (X, Z) are parallel and share the same origin. Values of θ for FCC and BCC crystals derived by Rice [53] are 35.25°and 54.7°, respectively, corresponding to á ñ{¯} 110 111 slip in FCC and á ñ{¯} 111 110 and á ñ{¯} 111 121 slip in BCC. Also considered here is prismatic slip in hexagonal crystals, á ñ {¯} 1120 1100 , whereby θ=60°. Planar shock loading is directed along θ=0°, so that, from symmetry considerations, only two systems at ±θ are subjected to nonzero resolved shear stress and can contribute to plastic deformation ( figure 4(a)). The plane strain geometry is representative of dislocation loops in the Smith-Hornbogen interface [29,30,32], with only edge components of the loops moving inside the shock front and contributing to plastic flow [25].
Isotropic polycrystals could be addressed by selection of two effective slip planes oriented at θ=45°c orresponding to planes of maximum resolved shear stress [16,39] and use of isotropic elastic constants. Extension of the model to textured polycrystals is not as straightforward but could be attempted if effective anisotropic polycrystalline elastic constants and preferred directions for plastic shear can be estimated. Shown in figure 4(b) is the elastic precursor front centered at point X for time instant t. The impact face corresponding to the initiation of wave propagation is located at X=x=0. The width of the shock in the transverse (Z) direction is W=constant. Mobile dislocations in the wave front of infinitesimal thickness = X U t d d , as well as some dislocations behind it, contribute to stress relaxation at the wave front. The variable propagation speed of the shock is U[X(t)], while each rarefaction source located at coordinate ξ contributes a longitudinal waveform moving at time-and position-dependent velocity c l (ξ, t). A differential element of material containing sources is labeled x d . Following results presented in section 3.3.3 and figure 3(c), shear waves emanated from dislocations at the front and behind it propagate too slowly to reach or affect the front. This assumption (i.e., null shear wave contributions) is also used elsewhere in the linear theory for dislocations existing [16] or homogeneously generated [26] at the precursor front.

Intrinsic attenuation
Considered now is the first term on the right side of (4.1). Generalizing the linear anisotropic representation of relaxation function F in [13] to account for elastic nonlinearity results in The plane strain geometry of figure 4 is invoked, as is plastic incompressibility = -D D zz P xx P . Recall that by the present sign convention, D xx P and P are taken as positive for compression. Voigt notation applies for B IJ . The anisotropic stress-state dependent shear modulus G(P) is first defined in (3.19). Modulus G in (4.2) replaces the isotropic pressure-dependent version used in [39]. The plastic velocity gradient (i.e., Eulerian plastic strain rate) at the shock front contributing to F is assumed to consist solely of contributions of mobile edge dislocations located at the front [14][15][16]; generation of new dislocations is assumed to require a finite time such that their subsequent motion only affects the second (trailing rarefaction) term in (4.1). For the plane strain geometry and two slip systems shown in figure 4(a), following [16], Recall that ρ D0 is the total, unpartitioned initial mobile dislocation density encountered by the precursor front, and υ D is the dislocation velocity taken here and henceforth as positive in sign. In the final expression, equal numbers of forward and backward moving dislocations are assumed, where the former correspond to the multiplier q -( ) 1 2 cos and contribute more strongly to precursor decay. An explicit partitioning of the plastic strain rate into contributions from lattice displacement incurred by dislocation generation and inelastic slip from dislocation glide has been invoked elsewhere [10,54,55]. Such a partitioning is not invoked here; rather, for simplicity, the standard Orowan-type prescription inherent in (4.3) is used.
A simple constitutive model is used to relate dislocation velocity to shock stress P. As proposed in [39], a linear viscous drag model of the following form is implemented: Here, υ D0 is a nominal velocity that will be related to a drag coefficient below, τ is the resolved shear stress for the relevant slip plane and direction, and τ H is the resolved shear stress at the (final, fully decayed) HEL. The ratio of shear stresses τ/τ H is exactly equal to the ratio of shock stresses P/P H only for linear elasticity, which is reasonably accurate for the small elastic strains encountered at the final minimum HEL. Denoting the dislocation drag coefficient byB, The viscous drag coefficient is computed using the following relation from Kocks et al [56] that requires no tunable parameters: Temperature rise contributing to phonon drag is ignored in the constant temperature (T=T 0 ) approximation, a reasonable assumption for the elastic precursor but perhaps less so for the plastic wave. Boltzmannʼs constant is k B and the atomic volume in the undeformed state is Ω 0 . Since (4.4) omits relativistic effects that would prohibit supersonic dislocation glide, in calculations the velocity υ D is set to the minimum of that computed from (4.4) and 0.9c s0 , i.e., dislocation velocities are restricted to a maximum of 90% of the linear elastic shear wave speed. Computed values of υ D0 are listed in table 1 for Cu, W, and Mg.
Combining (4.2), (4.3), and (4.4), the intrinsic contribution F becomes the following function of normalized shock stress at the precursor frontP and slip plane orientation θ: Since θ is constant for a given crystal structure, it will be dropped from the list of arguments hereafter.

Hydrodynamic attenuation
The second term on the right side of (4.1) requires more in-depth consideration. The stress gradient immediately behind the front, ∂P/∂X, evidently cannot be readily computed without enacting a full-field simulation of the elastic-plastic wave structure in time and space, e.g., via finite difference methods (e.g., with artificial viscosity) or minimally via a semi-linear method of characteristics [14,15]. The approach taken in the present work involves substitution of this hydrodynamic attenuation term with an alternative expression related explicitly to the dislocation generation rate behind the precursor front. A direct computation of the stress gradient is avoided in the present approach, but the end result of an increased rate of precursor decay due to rarefaction sources behind the front is physically relevant and qualitatively, if not quantitatively, similar.

Wave interactions
Relaxation waves generated by dislocations located at the front considered in section 4.2 and in [16] are assumed to be immediately swept into the precursor and subsequently travel at the instantaneous elastic shock velocity U. These longitudinal waves simply decrement the precursor magnitude in a additive manner, without reflection or transmission. In contrast, trailing rarefaction waves emanated some distance behind that ultimately reach the shock front at a higher velocity (c l >U ) should interact differently. The following interaction scheme is proposed for those longitudinal wave signals arising from generated mobile dislocations that reach the precursor front: • The shock acts as a moving boundary separating regions of impedance ρ c l (behind) and ρ 0 c l0 (ahead), leading to a fraction f t of the arriving waveform to be transmitted and a fraction f r to be reflected [57,58].
• The transmitted fraction f t of a waveform passes through the shock and immediately decelerates since it enters the region of lower impedance. This fraction is then swept into the precursor front in the same way as signals generated by dislocations at the front of density ρ D0 that cause intrinsic attenuation.
• In accordance with the point above, an effective density change of mobile dislocations δρ D (t) is defined that would lead to the same decay rate at the precursor front at time t if these dislocations were to have existed at the front at time t. In other words, the product f t δρ D contributes to precursor decay completely analogously to ρ D0 . However, those dislocations entering δρ D are not physically located at the front at time t, but rather are distributed over the region behind it, at points ξ<X in figure 4(b).
The fractions of transmitted and reflected signals are estimated here using the linear elastic relations [59] z z z z z z z z z r z r = + = - Since ζ 0 <ζ, f r <0 and the reflected signal experiences a sign change relative to the incident wave: here the reflected part of the incident pulse is compressive. Though such reflected pulses may influence the material state some distance behind the shock, they do not affect the precursor at time t and will not be considered further. In later calculations, most of the signal is transmitted rather than reflected due to the relatively small impedance mismatch across the elastic shock front, e.g., f t  0.9 and  | | f 0.1 r . Though of relatively low magnitude in such cases, reflected pulses could alter the stress state responsible for dislocation generation behind the precursor front and in the plastic waveform, and would also interact with oncoming release waves from other dislocations. Such alterations and interactions are beyond the scope of the present analytical treatment. Interactions among simple nonlinear waveforms can be addressed graphically using the method of characteristics [39,58], but the complex waveforms of present interest associated with dislocation dynamics appear to require more sophisticated numerical methods of analysis.
Combining the results of (4.7) derived in section 4.2 with the above-stated assumptions for the form of the hydrodynamic term, (4.1) can be rewritten as In the denominator of the prefactor, the ratio c U l 2 2 -where here c l is the limiting value of isentropic sound speed just behind the precursor front-can be obtained as a function of shock stressP using (3.5) and (3.16), with characteristic results given already in figure 3(a). Notice that both intrinsic and hydrodynamic contributions from dislocations (terms in parentheses on the far right) are affected by elastic nonlinearity via the stress dependent shear modulus G and the effect of c l >U, both of which tend to increase the decay rate. The factor f t can likewise be obtained as a function ofP via (3.16) and (4.8). The only remaining unknown term on the right side of (4.9) is δρ D , the cumulative contribution of generated trailing dislocations to the effective plastic strain rate at the precursor wave front.
The crux of the present model is that the effects of stress signals produced by the virtual density of dislocations generated in the wake of the shock, δρ D , are essentially equivalent in (4.9) to those intrinsically arising from sources located at the front. Their transmitted fraction f t δρ D interacts with the precursor shock in the same manner as those pre-existing dislocations set in motion immediately at the front. It will be defined precisely later how δρ D accounts only for those signals that reach the front in a given time increment. Linear elastic concepts are necessarily invoked in (4.8), and likewise implicitly for superposition of wavelets from individual dislocations comprising the defect density field. Analysis of interactions among finite amplitude, nonlinear elastic waves from transient discrete dislocations or their mobile densities would require an extensive numerical treatment beyond the scope of this work. Physical validation of the current treatment would perhaps benefit most from comparison with results from molecular dynamics (MD) simulations that should be able to realistically address transient, fully nonlinear dislocation stress fields and their interactions with stress or shock waves. However, MD simulations that focus on the desired physical interactions do not seem to exist in the current literature, and length and time scales achievable are far smaller than those commensurate with standard precursor decay experiments.

Dislocation generation
Referring to figure 4(b), let the dislocation generation rate per unit area at material point ξ<X be denoted by r ẋ (ˆ) t , D . Argument Î [ ] t t 0, covers the time history of precursor decay up to current time t corresponding to the Lagrangian location of the precursor front X(t). A constitutive equation for ṙ D will be introduced later. Consider a differential element of material x d centered at ξ with this generation rate. The absolute number of dislocations generated per unit time in this differential area is with W the transverse width of the shocked domain as defined in figure 4(b). The earliest start time t s at which generation is possible is the time at which the shock front crosses point ξ. Rarefaction signals emanated from point ξ could reach the shock front by the time the shock is located at point X only for those dislocations generated over a window of time [t s , t f ], where t f is the latest possible time that signal emanated from (ξ, t f ) could reach (X, t) by time t>t f . For example, waveforms from dislocations generated farther behind the precursor front will require more time for their signals to reach point (X, t) than those from dislocations generated closer to it, assuming for the moment that the propagation speed c l >U is uniform among signals. Signals emanated from ξ at > t t f do not have time to reach the front at or before the frontʼs arrival at point X at time t.
Since the precursor stress decreases with increasing time and the separation distance between the precursor and plastic wave increases with increasing time, the effective wave speedc l over the path from ξ to X also generally decreases with increasing time. This suggests that a signal emanated from ξ at a time >t t 2 1 would never catch a signal emanated from the same point ξ at an earlier timet 1 .
Now consider a small time window  + t t t d . During this window, the location of the precursor front will advance as  + = + X X X X U t d d. A rarefaction signal trailing closely behind the front will advance as For the latter signal to just reach the precursor front at the end of the time window, the starting point of the rarefaction signal must be = -- The time increment dt required for the rarefaction wavelet moving at velocity c l to cover the difference -X X is where c l and U are the longitudinal wave speed just behind the front and precursor shock velocity at time t.
Combining this result with (4.10), the increment in effective dislocation content acting at the precursor front at (X, t) due to the generation of dislocations in the element x d located at ξ is Taking a larger value of dt than that prescribed by (4.11) would incorrectly lead to contributions of signals from dislocations that trail too far behind the front to be able to reach it during the window d . These signals instead can only affect the precursor front at later times. Time increment dt must be 'infinitesimal' since a finite time interval would lead to 'piling up' of relaxation contributions that already affected precursor decay at earlier times. As  c U l , no catching up with the precursor front is possible and The temporal argument of the dislocation density rate r ẋ (ˆ) t , D is defined as follows: Here,c l is the effective velocity of the rarefaction signal over its transit from ξ to the shock front at X. Recall that t f is the latest time instant at which a signal generated at ξ would just reach the front at (X, t). Since signals emanated earlier in time from the location ξ are presumed to not be overtaken by those emanated from the same location later, wavelets from generation at ξ at times earlier than t f will have already interacted with the precursor front. The total contribution of relaxation sources produced by all material elements behind the front is obtained by integrating (4.12) over the shock path: Here, W is a constant and U and c l are values at the front, independent of ξ, that can be extracted from the integral as shown. On the other hand, ṙ D generally depends on ξ. The effective density change of dislocations in the shock front is obtained by dividing (4.14) by the differential area element swept by the front = W X WU t d d: where now all quantities are evaluated at time t and point X as the elastic shock front is approached from the rear. A constitutive model is needed for the source generation rate. Sophisticated expressions accounting for effects of mobile versus immobile dislocations, dislocation interactions, saturation, and recovery have been proposed and implemented elsewhere [10][11][12]60]. Accompanying the linear drag model used for dislocation velocity, a basic breeding model [60] is invoked for defect generation: The breeding coefficient m is a material parameter with dimensions of length, and M has the same dimensions as dislocation density. At the moving elastic precursor front, ρ D =ρ D0 (a fixed density of sources), and the dislocation velocity is given by (4.4). Thus, at the precursor front, where t 0 is a constant with dimensions of time. According to (4.21), the dislocation density rate simply increases linearly with increasing stress above the HEL. Inserting (4.21) into (4. 19), and then subsequently using the result in (4.9) gives the following final expression for precursor decay: This is a nonlinear first-order ODE with all terms on the right side depending potentially only on P and/or t. Subsequent derivations and results invoke (4.22) in non-dimensional form.

Dimensionless form
The following dimensionless quantities are defined, noting thatP is repeated from (3.6):

4.25
This differential equation will be solved subseqently in section 5 for precursor decay in Cu, Mg, and W single crystals. The initial condition is where P 0 =P(t=0) is the impact stress, i.e., the axial Hugoniot stress in the plastic wave, and P H is recalled as the final equilibrium HEL after precursor decay is complete. All parameters and functions entering the right side of (4.25) have been defined or calculated from the crystal structure, (non)linear elasticity, and basic dislocation mechanics, with the exception ofM . The latter parameter, which will be prescribed later, controls the dislocation generation rate and the resulting contribution of hydrodynamic attenuation to precursor decay. Solutions to (4.25) are obtained by direct numerical integration. Furthermore, onceP is known at timet , a cumulative dislocation density ρ D at any dimensionless timet 1 is obtained by integration of (4.21): Accordingly, the precursor stress in dimensional form decreases exponentially from P 0 to P H with increasing time t.  D0 is initially viewed as adjustable, recalling that m is the breeding coefficient for first-order dislocation kinetics [60].

Results
Results are presented in figure 5 for four values ofM , noting that = M 0 suppresses hydrodynamic attenuation from dislocation generation in the present nonlinear model of (4.25). The linear result excludes hydrodynamic attenuation from dislocation generation, the stress-state dependence of the shear modulus, and the possibility of c l >U. The difference between the linear result and the nonlinear result with = M 0 is very small in Cu and W but visibly measurable in Mg. Incorporation of dislocation generation via > M 0 is demonstrated to strongly affect precursor decay according to the present theory.
Model predictions are compared with experimental data in figure 6, where stress and propagation distance are of their true physical dimensions. The position of the shock front is obtained from the decay time t via use of the first of (4.16); the constant shock velocity approximation is not needed. Only the nonlinear result with = · M 5 10 5 is shown in each case, as it provides the closest match to experiment. This finding is reassuring, if not somewhat remarkable, that the same value of this parameter can be used to fit all three data sets. As is evident in figure 6, the linear model grossly under-predicts the magnitude of precursor decay when the measured values of initial dislocation density ρ D0 are used (table 1). This finding is in agreement with early results of [13,18] where similar, but not identical, models of dislocation velocity accounting for drag also produced much too low rates of decay for the same initial defect densities.
Consider the chosen value of = · M 5 10 5 in more detail. For the initial dislocation densities used here, this choice corresponds to = » --· · M m b 5 10 5 10 m 15 16 2 . Values of M tabulated by Gilman [60] for various metals are M≈10 13 −10 16 m −2 . Thus the value ofM used here, though near the upper end of the range reported in [60], appears physically admissible. Furthermore, even if M were reduced by two orders of magnitude, the effects of dislocation generation and corresponding hydrodynamic attenuation on precursor decay would remain substantial, as is evident from results for = · M 5 10 3 in figure 5, though lower values of breeding coefficient would not enable as close agreement with experiment.
Results corresponding to these experimental conditions address relatively large shock stresses, over 20 or 30 times the final equilibrium HEL for Cu and Mg in particular. Lower impact stresses, corresponding to only 20% in excess of the HEL, are considered in model predictions shown in figure 7. Differences between the nonlinear and linear results, though smaller than those observed at higher stresses in figure 5, are still substantial.  . Precursor stress versus propagation distance compared with experimental data for Cu [13], Mg [40], and W [13]. For the nonlinear results, = · M 5 10 5 .
Differences are smallest for Cu as a result of its relatively low ratio of c l /U compared to that for Mg and W ( figure 3(a)), which in turn leads to a relatively low value of b entering (4.25). The exponential decay of the linear solutions is more evident now than in figure 5 due to the longer domain of timet chosen here.
Remarks on experimental uncertainty and sparsity of test data are in order. Uncertainty in HEL data on Mg from [40] is reported as ±0.015 GPa, on the order of ±5% of the values of P corresponding to data points in figure 6(b). Corresponding error bars are not shown since these would be smaller than the data points themselves in the figure and would not be visible. Uncertainties in data for Cu and W are not reported in [13], but contemporary experimental data for LiF [19,21] suggest uncertainties in the measured HEL ranging from 1% to 10% and uncertainty in the impact stress on the order of 2%. Because the experimental data are sparse, little constraint on model predictions is enabled, and the apparent agreement of the nonlinear model with the data confirms only that the calculations are physically reasonable but does not provide strong validation of the theory. However, sparse data for these and other crystal orientations and materials were also used in [13] to invalidate linear precursor decay models that assume constant dislocation densities. Similar conclusions were drawn from analysis of more extensive data on LiF in [18]; mechanisms responsible for precursor decay in LiF were deemed not well understood [20].
Contributions of nonlinearity to precursor decay are now closely examined. Calculated effects ofP on functions f t of (4.8) and ā and b of (4.23) are listed in table 2 for each material. With increasing precursor stress, f t decreases but both ā and b = -· ( ) f U c 1 t l increase. The latter increases contribute to the relatively greater effects of nonlinearity on precursor decay at large stresses evident upon comparison of figure 5 with figure 7. Recall from (4.25) that the hydrodynamic contribution from dislocations is represented via the product b ·¯·P t. The factor b accounts for the approach and interaction of trailing rarefaction signals with the precursor front. The factorP arises from the dependence of the generation rate of dislocations, the sources of such signals, on stress. Finally, the factort accounts for the integrated contributions of all such sources over the domain ξ ä [0, X), where X is the Lagrangian distance traversed by the shock at timet . During the attenuation process,P decreases rapidly andt increases linearly. Thus, the maximal net contribution of hydrodynamic attenuation in the present model is not initial since a finite time is needed for rarefaction signals to generate and propagate to the wave front. Nor is it towards the end of the attenuation process since thenP is very small and the generation rate of new dislocations behind the front is low.  Prior experimental studies [17,18], where hydrodynamic attenuation was computed from stress gradients via relations such as (2.5) or (2.6), found nonlinear effects to be important but not dominant, for example responsible for the order of 10%-40% of the total decay rate, with intrinsic attenuation accounting for the remainder. However, as noted in [16,21], the measurement of stress gradients in such studies suffers from extreme imprecision: the partitioning of hydrodynamic and intrinsic attenuation is so uncertain that even the sign of F may be in error for X>0.5mm. In experiments on LiF [18], the relative contribution of hydrodynamic to intrinsic decay was found to increase with propagation distance X. In contrast, the relative contribution of nonlinear elastic effects was found to dominate the precursor decay rate at early times in numerical simulations of Cu [14]: a decay rate of 50 times higher in nonlinear simulations than in simulations with a constant sound speed was reported. Significant effects of stress dependence of both compressibility and shear stiffness on precursor decay were also stated in [15]. Numerical results [14,15] suffer from sources of imprecision inherent in usual continuum plasticity simulations, for example those associated with viscoplastic constitutive models and numerical integration methods.
Cumulative dislocation densities computed via (4.27) are reported in figure 8 for impact stressesP 0 corresponding to the experimental conditions examined in figure 5 and the lower stress regime examined in figure 7. For the former higher stress conditions, the cumulative dislocation content over the distance traversed by the precursor saturates at values on the order of 10 3 ∼10 4 times the initial dislocation density ρ D0 , corresponding to dislocation densities on the order of 10 13 ∼10 15 m −2 . These values are comparable in magnitude to defect densities observed in recovered experimental samples of single crystals of Cu [61] and in dislocation dynamics simulations of single crystals of aluminum [26].

Discussion
With regard to precursor decay, the main findings from the present analysis are summarized as follows: • Effects of elastic nonlinearity are mildly important when dislocation generation behind the wave front is omitted; • Hydrodynamic attenuation associated with the combination of elastic nonlinearity and generation of rarefaction sources (dislocations) behind the wave front appears to be very important; • A single value of the breeding coefficient, within the upper limit of experimental validity, can be used to enable close agreement between model and experiment for all three metals (Cu, Mg, and W); • Hydrodynamic attenuation can dominate intrinsic attenuation as the initial impact stress increases, and it remains significant even at impact stresses not greatly exceeding the final HEL.
It is conceded that a number of simplifying assumptions have been invoked to arrive at an ODE for nonlinear precursor decay. As also assumed in [18,19,21], the state behind the shock has been idealized as homogeneous and non-inertial at each increment in time during the decay process in order to compute instantaneous wave velocity c l for rarefaction signals via third-order isentropic elasticity. Basic constitutive equations, essentially linear in stress, have been implemented for dislocation velocity and dislocation generation. The non-reflected fraction of rarefaction signals that catch the precursor front has been assumed to affect the front in a simple way similar to intrinsic attenuation from pre-existing dislocations. Furthermore, transient conditions at or immediately behind the precursor front have been used to estimate rates of generation at locations further behind the front, presumably tending to under-predict generation for regions in the plastic rise. Nevertheless, cumulative dislocation densities computed according to the present model appear physically reasonable ( figure 8). Dislocation interactions, including immobilization and annihilation, for example, have not been addressed. As shown in appendix, a reduction in hydrodynamic attenuation would be expected to occur if rarefaction waves were themselves subject to decay in magnitude while transiting from their source points to the shock front. It is recommended that future numerical continuum studies re-examine effects of hydrodynamic attenuation, modernizing the prior numerical work reported in [14,15], for example, with state-of-the-art models for viscoplastic crystal response and associated dislocation kinetics [10][11][12]62]. In particular, the finite difference model reported in [11], which invokes third-order Eulerian elastic constants, has been used to predict precursor decay in textured [111] Al crystals [63] within a factor of ≈2 accuracy. A one-to-one comparison with the present analytical treatment is not possible since different materials, slip geometries, viscoplastic flow rules, and dislocation evolution laws are invoked in [11], along with artificial viscosity that would be expected to smooth any precursor spike. Exploratory calculations [unpublished] with various choices of nonlinear elasticity parameters (Lagrangian and Eulerian models, with and without higher order elastic constants of orders three and four) have demonstrated different wave profiles than third-order Eulerian results reported in [11]; however, such differences do not enable isolation of hydrodynamic effects since the model in [11] is geometrically nonlinear, with a variable sound speed that would induce hydrodynamic attenuation even if all higher-order elastic constants are zeroed. A more thorough study of effects of features of this and relatedly sophisticated models [10][11][12]62], outside the scope of the present paper, on precursor decay appears needed.
The primary weakness of the present approach is the unknown accuracy of its major assumptions, in particular its treatment of how rarefaction waves from trailing dislocations are generated, combined, propagated, and interact with the precursor front. The validity of these assumptions can only be verified by comparison with other methods, experimental or theoretical-computational, that directly interrogate the physics involved. Experimental methods with the requisite fidelity at such small (nano) length and time scales appear to be presently unavailable, as do discrete dislocation dynamics computations that account for nonlinear elasticity essential to capture hydrodynamic attenuation. Discrete lattice dynamics [64] or modern molecular dynamics [65,66] can be used to consider elastic nonlinearity and/or complex dislocation kinetics at the nanoscale. For example, metastable elastic zones of finite length, essentially resulting from nonlinear elastic behavior and pulses emitted from trailing dislocations, were observed in simulated overdriven shocks in aluminum [65]. However, classical MD simulations cannot be conducted on sample sizes large enough for direct comparison to experiments at the mm scale. To address all of these issues with 'first-principles' methods, implementation of atomistic scale-bridging techniques and tools for interrogating discrete nonlinear stresswave interactions in the context of shock wave propagation would seem essential.

Conclusions
A theoretical investigation of effects of elastic nonlinearity on precursor decay has been reported. Decay rates depend on intrinsic attenuation, due to pre-existing dislocations at the wave front, and on hydrodynamic attenuation, due to rarefaction signals from sources behind the front that travel at a higher velocity. Nonlinear anisotropic elasticity with elastic constants up to order three has been invoked in Eulerian form to estimate velocities of rarefaction signals, which here always exceed the shock velocity, for three metals of interest: Cu, Mg, and W. An equation for precursor decay accounting for intrinsic attenuation and hydrodynamic attenuation from elastic nonlinearity with dislocation generation has been newly derived. Introduction of basic constitutive models for dislocation velocity and generation permits reduction of this relation to a nonlinear ordinary differential equation. The hydrodynamic term from dislocations is controlled by a single scalar parameter, the breeding coefficient for dislocation generation. Solutions for Cu, Mg, and W have demonstrated how nonlinearity can dominate the overall precursor decay rate for physically bounded dislocation generation rates behind the shock; such solutions also show reasonable agreement with experimental data. Future numerical studies should invoke models of higher fidelity to further investigate nonlinear effects and confirm or refute their significance.

Appendix. Decaying rarefaction signals
The model proposed in section 4.3 of this paper treats those rarefaction signals that catch the precursor front, and do not reflect off of it, as similar in form and magnitude to those produced by dislocations at the front. In this context, referring to figure 4(b), generated dislocations along a line { ξ=constant, 0ZW } collectively produce a 1-D rarefaction wavelet traveling in the X-direction at velocity c l . The strength and form of this wavelet are assumed unchanged as it propagates from ξ to the front at X. However, if individual dislocations identical to those reported in figure 5 and likewise adequately match experimental data in figure 10. Such cases correspond to the limiting result of long t c or large R c noted above.
Ast c decreases, the decay of trailing rarefaction signals becomes more rapid, and hydrodynamic attenuation is impeded. Trends are most similar for Cu and Mg; the same values oft c more strongly impede hydrodynamic attenuation in W. The difference between curves for =t 10 c 2 and =t 10 c 3 becomes apparent for solution times t 10 3 , whereby the distance between early sources and the precursor front becomes large enough that source attenuation becomes important. Mathematically, such importance corresponds to non-negligible higher-order terms in the series expansion of (A.6).
Prescription of =t 10 c 2 corresponds to decay distances R c of approximately 13 mm, 3.1 mm, and 1.9 mm respectively in Cu, Mg, and W. Referring to figure 10, these distances are on the order of the propagation distance X traveled by the precursor shock over its period of prominent attenuation. Thus, the present findings show that any decay of rarefaction signals that might occur during their transit to the precursor front would reduce their contribution to precursor attenuation. On the other hand, especially for Cu and Mg, hydrodynamic attenuation is still significant relative to intrinsic attenuation (e.g., compare with the linear results) even under a reduction in the decay time by three orders of magnitude to =t 10 c 5 , which would correspond to decay distances R c of only 13 μm in Cu and 3.1 μm in Mg. A quantitative estimate of the true extent of any decay of rarefaction signals from dislocations nucleated in single crystals with dimensions of mm could perhaps be estimated using dislocation elasto-dynamics simulations, presumably omitting nonlinear elastic effects.