# Shock-driven nucleation and self-organization of dislocations in the dynamical Peierls model

Yves-Patrick Pellegrini, Marc Josien

► To cite this version:

Yves-Patrick Pellegrini, Marc Josien. Shock-driven nucleation and self-organization of dislocations in the dynamical Peierls model. Physical Review B, 2023, 108 (054309), [15 p.]. 10.1103/Phys-RevB.108.054309. cea-04179210

HAL Id: cea-04179210

https://cea.hal.science/cea-04179210

Submitted on 9 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Shock-driven nucleation and self-organization of dislocations in the dynamical Peierls model

Y.-P. Pellegrini¹,²,\* and M. Josien³,†

¹CEA, DAM, DIF, F-91297 Arpajon, France.

²Université Paris-Saclay, LMCE, 91680 Bruyères-le-Châtel, France

³CEA Cadarache, F-13108 St-Paul-lez-Durance, France.

(Dated: August 8, 2023)

Dynamic nucleation of dislocations caused by a stress front ('shock') of amplitude $\sigma_a$ moving with speed $V$ is investigated by solving numerically the dynamic Peierls equation with an efficient method. Speed $V$ and amplitude $\sigma_a$ are considered as independent variables, with $V$ possibly exceeding the longitudinal wavespeed $c_L$. Various reactions between dislocations take place such as scattering, dislocation-pair nucleation, annihilation, and crossing. Pairs of edge dislocation are always nucleated with speed $v \gtrsim c_L$ (and likewise for screws with $c_L$ replaced by $c_S$, the shear wavespeed). The plastic wave exhibits self-organization, forming distinct 'bulk' and 'front' zones. Nucleations occur either within the bulk or at the zone interface, depending on the value of $V$. The front zone accumulates dislocations that are expelled from the bulk or from the interface. In each zone, dislocation speeds and densities are measured as functions of simulation parameters. The densities exhibit a scaling behavior with stress, given by $((\sigma_a/\sigma_{th})^2 - 1)^\beta$, where $\sigma_{th}$ represents the nucleation threshold and $0 < \beta < 1$.

To appear in Physical Review B

# I. INTRODUCTION

Shock loading of metals [1, 2] involves a high rate of dislocation [3] generation [4, 5]. Realistically accounting for it remains a challenge for dislocation-based multiscale models of shock-loaded solids aimed at very high strain rates [5–9], in spite of a vast amount of data from atomistic simulations [10–16]. Shock-induced plasticity has also been explored using two-dimensional, non-supersonic, elastodynamic dislocation simulations with retarded interactions and separate nucleation criterion and mobility law [17]. At high strain rates, line multiplication from traditional sources [18] via the Orowan mechanism is no more the dominant generation mechanism [19], and nucleation must be considered instead [17]. Whereas standard thermally-activated homogeneous nucleation [20] is usually invoked below the nucleation stress threshold, no consensus exists for larger driving stresses. However newly-nucleated dislocations can trigger further nucleation events [21, 22].

Most dislocation-based models rely on the assumption that the plastic strain rate is proportional to the density and speed of mobile dislocations by Orowan's equation [23, 24]. Thus, speed is also a much debated issue. In uniaxial compression, the Hugoniot elastic limit [2, 25] where plasticity comes into play is first attained for shear on planes inclined $45^\circ$ from the compression axis. Dislocations accompanying a compression wave of (longitudinal) speed $c_L$ should then move with supersonic speed $v = \sqrt{2} c_L$ [26]. Early models [27] postulated layers of intersonic or supersonic edge dislocations [26, 28] at the shock front to accommodate shock-induced compression, with residual dislocation density in the bulk. To dispense with supersonic dislocations, alternative proposals [26, 29] considered dislocation loops homogeneously and repeatedly nucleated

behind the shock front by the deviatoric stress set up by uniaxial strain, moving thereafter only short distances at sub-sonic speeds. Scarce information is available concerning supersonic dislocation speeds in shocks, and in absence of direct evidence the experimental status of supersonic dislocations in metals remains unsettled [29]. Still, atomistic simulations have revealed that intersonic or supersonic dislocations are quite possibly involved [30–32], although current dislocation-based density evolution models disregard this possibility.

Leaving aside thermal effects [25, 33, 34], which are presumably irrelevant [26] at high stress except for changes in elastic constants, wavespeeds, and phonon drag, the present work investigates nucleation of perfect dislocations in a one-dimensional (1D) idealized model for straight dislocations on a single slip plane, with driving stress exceeding the nucleation threshold, and variable 'shock' speed. To this aim we use the (elasto-) Dynamic Peierls Equation (DPE) [35, 36], which we solve numerically [37]. This equation for the plastic slip rests on minimal assumptions within the small-strain formulation of continuum mechanics, assuming fixed elastic constants. Its properties with regard to dislocation nucleation and collective behavior are unexplored so far. By construction, it accounts for stress-wave emission and scattering [38] by/on dislocations in elastodynamic interaction. Dislocations are free to glide (i.e., all mobile) with no obstacles, nor Peierls stress. Supersonic motion is allowed. Once nucleated as pairs of opposite signs [26], they interact via retarded elastodynamic interactions. Dislocation trajectories as well as quantities usually considered in constitutive models, such as the plastic strain rate $\dot{\varepsilon}^P$, dislocation density $\rho$, speed $v$ and stress, are monitored. Our approach is motivated by simplicity and a direct connection between this continuum description and Orowan's equation. More accurate analogous lattice-based dislocation models [39] might be less suitable to this purpose.

The DPE is reviewed and cast in Sec. II in a form suitable to numerical solution (see Appendix B for the method). Sec. III is devoted to simulations. After the setup is explained,

<footnote>\* yves-patrick.pellegrini@cea.fr</footnote>

<footnote>† marc.josien@cea.fr</footnote>
the plastic structure arising from dislocation trajectories is discussed, with mean dislocation speeds and speed probability densities measured. Nucleation is addressed next, and mean dislocation densities are computed over a wide range of ‘shock’ speeds and driving stresses. A scaling in stress of the density is proposed on the basis of a supersonic-train solution derived in Appendix C. Concluding remarks close the paper (Sec. IV). Appendices are devoted to the numerical method. Our numerical results privilege the edge dislocation, as its behavior is more involved than the screw’s. Additional results for screws are reported in the Supplemental Material [40].

## II. DYNAMICAL PEIERLS EQUATION

### A. Overview and former results

The DPE extends the classical Peierls model of a dislocation [3, 41] to elastodynamics [42–45]. It reduces to the Weertman equation [46–48] for steady motion under uniform applied stess $\sigma_a$. In Ref. 36 – hereafter referred to as (I), and for a sinusoidal (Frenkel) pull-back force, the elastodynamic solution of the DPE was approached via a collective-variable approximation (CVA), using a single-dislocation ansatz depending on two collective variables – dislocation position and width, for which evolution equations were obtained. The CVA is such that its stable steady-velocity states at large times match those of the DPE for one dislocation, which are given by the single-dislocation solution of the Weertman equation. The resulting mobility law [36, 47] is recalled in Fig. 1, where $\sigma_{th} = \mu b/(2\pi d) \sim \mu/10$, is the theoretical shear stress (lattice shear strength) [3, 23, 49] defined in terms of the shear modulus $\mu$, the Burgers vector $b$ and the interplane distance $d$ [3]. Hereafter, $c_S$ and $c_L$ denote the shear and longitudinal wave speeds of the medium. The screw dislocation has only one subsonic branch that saturates at $v = c_S$. The edge dislocation has two stable branches: one subsonic ($|v| < c_S$), the other one intersonic ($c_S < |v| < c_L$). The subsonic branch saturates at the Rayleigh velocity $v = c_R \approx 0.93 c_S$. No single-dislocation mobility law exists for $\sigma_a > \sigma_{th}$ in this model, and the speed is undetermined if $\sigma_a = \sigma_{th}$. The inclusion of a gradient term provides a means to overcome this limitation [47]. However, our study demonstrates that under high-stress ‘shock’ conditions the DPE yields multi-dislocation supersonic solutions without needing any such modification.

Atomistic calculations for an edge dislocation [30] suggested a dependence of branch selection on the way the external loading is applied. This was further investigated in (I), leading to the following predictions: a) when $\sigma_a$ is suddenly applied and thereafter maintained (single-step loading), a critical stress $\sigma_c$ exists that separates regimes with subsonic or intersonic terminal dislocation speeds. This transition is represented at $\sigma_c \approx 0.41 \sigma_{th}$ in Fig. 1. Thus the model accounts for the jumps from one velocity branch of the mobility law to the other as stress increases; b) the terminal speed is dynamically selected and depends on the history of loading. Two-step loading, where a first stress level is imposed and subsequently modified after some time, grants access to asymptotic states

<figure><img src="image_2.png" /><figcaption>Figure 1. (Color online) Single-dislocation steady-state stable branches of the mobility law for screw and edge dislocation in the PND model, for $c_L/c_S = 2$ and $\alpha = 0.05$ (see below for a definition of this quantity). Edge case: the solid-line part depicts states accessible under single-step loading; the dashed part depicts states accessible under double-step loading [36]. The jump between subsonic and intersonic branches occurs at $\sigma_a = \sigma_c$ (see text).</figcaption></figure>

forbidden under single-step loading [36]. In anisotropic materials more branches come into play [31, 50–53].

### B. The equation

The DPE for the plastic slip function $\eta(x, t)$ of a system of straight dislocations moving along the x-axis reads

$$ -\frac{\mu}{\pi} \int dx' \int_{-\infty}^{t} dt' K_a(x - x', t - t') \frac{\partial \eta}{\partial x}(x', t') \\ - \kappa_a^{\alpha} \frac{\mu}{2c_S} \frac{\partial \eta}{\partial t}(x, t) + \sigma_a(x, t) = f(\eta(x, t) + \eta_e(x, t)). \quad (1) $$

The left-hand side of this equation represents the local stress $\sigma(x, t)$ on the slip plane. The quantity $\sigma_a(x, t)$ is an externally applied stress, which depends on both position and time. The local stress is balanced by the pull-back force $f(\eta)$, which derives from the generalized stacking-fault $\gamma$-potential of the lattice [54]. We use here the Frenkel force $f(\eta) = \sigma_{th} \sin(2\pi\eta/b)$, which implies perfect dislocations, but our numerical method applies as well to any other force. This force favors locking of $\eta$ on integer multiples of $b$. Dislocations lie where a transition between two successive multiples take place. The theoretical shear stress $\sigma_{th}$ is its maximal value, which identifies with the nucleation stress [10, 26, 55] (see Sec. III B 4 below). Equation (1) implies that the local stress in the system is always limited by $\sigma_{th}$, that is, $|\sigma(x, t)| \le \sigma_{th}$, whatever $\sigma_a$. In (1) the function $\eta_e$ represents the elastic part of the slip, i.e., the dislocation-free solution of the equation $f = \sigma_a$ if $|\sigma_a| \le \sigma_{th}$. It increases with $\sigma_a$ up to its saturation value $b/4$ for $\sigma_a = \sigma_{th}$, and by convention remains fixed to that value if $\sigma_a > \sigma_{th}$ [36]. Kernel $K_a(x, t)$ accounts for retarded elastodynamic self-interactions on the slip plane of the field $\eta(x, t)$. It depends on the dislocation character $a = s$ (screw), $g$ (glide edge), $c$ (climb edge). The local instantaneous term with $\partial\eta/\partial t$ is proportional to the di-
dimensionless constant [56]

$$ \kappa_a^\alpha = \kappa_a(1 + \alpha). \qquad (2) $$

Whereas $\kappa_a$ embodies radiative losses associated with $K_a$, the contribution $\alpha\kappa_a$ with $\alpha > 0$ accounts phenomenologically for phonon drag. The inclusion of the phonon-drag term ensures a well-defined mobility law in the subsonic range depicted in Fig. 1, while also providing the required linear mobility behavior at low speeds [36, 46, 47]. Expressions of $K_a(x, t)$ and $\kappa_a$ are reported in Appendix A. It is worth noting that the elastic waves propagated by $K_a$ undergo no damping. Indeed, bulk elasto-viscosity, or anelasticity, is insignificant in pure pristine metal crystals [57]. Additionally, the potential inclusion of inertia in the pull-back force [42], while unexplored at present and lacking documentation, is intentionally disregarded in order to maintain tractability of equations in the steady state. The emphasis in the current work is placed on simplicity to enable improved control and facilitate interpretations based on existing references [36, 47].

Evolution starts at $t = 0$. Prior to this instant, $\sigma_a(x, t) = 0$, and the system can, e.g., be assumed to contain one dislocation at rest, of shape $\eta_0$. The initial state $\eta_0(x) \equiv \eta(x, 0)$ is a solution of the Peierls-Nabarro (PN) equation

$$ -\frac{\mu K_a}{\pi} \text{p.v.} \int \frac{dx'}{x-x'} \frac{\partial \eta_0}{\partial x}(x') = f(\eta_0(x)), \qquad (3) $$

where p.v. stands for Cauchy's principal value, $K_s = 1/2$ or $K_{g,c} = 1/[2(1 - \nu)]$, and $\nu$ is Poisson's ratio. Indeed the DPE equation reduces to the latter in statics [35].

The dynamic solution is defined relatively to this state as

$$ \eta(x, t) = \eta_0(x) + \delta\eta(x, t), \qquad (4) $$

where $\delta\eta(x, t) \equiv 0$ for $t \le 0$. For $t > 0$ the equation to be solved for $\delta\eta(x, t)$ is, with same dependency on $(x, t)$ as above

$$ -\frac{\mu}{\pi} \int dx' \int_{-\infty}^{t} dt' K_a \frac{\partial \delta\eta}{\partial x} - \kappa_a^\alpha \frac{\mu}{2c_s} \frac{\partial \delta\eta}{\partial t} = F, \qquad (5a) $$

where $F(x, t)$ is

$$ F \equiv f(\eta_0 + \delta\eta + \eta_e) - f(\eta_0) - \sigma_a. \qquad (5b) $$

In absence of initial dislocation, one simply takes $\eta_0(x) \equiv 0$.

Throughout the paper, the longitudinal wavespeed is $c_L = 2c_s$ (which corresponds to a Poisson ratio $\nu = 1/3$) and the drag coefficient is $\alpha = 0.01$ [36]. Results will be given in dimensionless form, with time in units of $\tau_0 = d/c_s$, positions in units of $d$, speeds in units of $c_s$, and stresses in units of $\sigma_{th}$ so that the values of $\mu, b$, and $d$ need not be precised.

The numerical method of solution is sketched in Appendix B. Prior to going further, the numerical solution of the DPE for one pre-existing single dislocation has been cross-checked against the collective-variable approximation with good overall agreement. In particular, predictions of the CVA regarding the subsonic-intersonic transition for an edge dislocation [36]

are fully confirmed (see Supplemental Material [40]).

### C. Local form of Orowan's equation

The 'microscopic' dislocation density $\rho$ and plastic strain rate $\dot{\varepsilon}^P$ are

$$ \rho = -\frac{1}{db} \frac{\partial \eta}{\partial x}, \qquad \dot{\varepsilon}^P = \frac{1}{d} \frac{d\eta}{dt}, \qquad (6) $$

where $d$ is the interplane distance [3]. If we neglect core-shape changes and let $\eta(x, t) = \eta(x - \xi(t))$ represent a dislocation train moving with velocity $v = \dot{\xi}$, a local instance of Orowan's equation [23] is retrieved from (6) as

$$ \dot{\varepsilon}^P = b\rho v. \qquad (7) $$

However, accounting for dynamical core-shape changes via an intrinsic time-dependence of $\eta$ yields the more general relationship

$$ \dot{\varepsilon}^P = b\rho v + b\rho \frac{\partial \eta}{\partial t}. \qquad (8) $$

Upon ensemble and time averaging, this expression should produce a correction to (7) of a type different from the one, of density-rate origin, proposed by Armstrong and Zerilli [4, 5] to account for nucleation. Indeed, Orowan's equation (7) was derived assuming a constant dislocation density [23]. Nevertheless, this equation can always (formally) be enforced at the macroscopic level by defining the mean speed as

$$ \bar{v} := \bar{\varepsilon}^P / (b\bar{p}), \qquad (9) $$

where the overbar encompasses all relevant averages [58]. With Eqs. (8) and (9), the mean speed may turn out very different from that resulting from the mobility law of individual dislocations.

### III. 'SHOCK' SIMULATIONS

#### A. Set-up

<figure><img src="image_3.png" /><figcaption>Figure 2. (Color online) Time-dependent applied stress profile $\sigma_a(x, t)$ expanding with speed $V$, for $\sigma_a = 1.05$.</figcaption></figure>

Highly-simplified 'shock'-like conditions are simulated by uniformly loading with stress $\sigma_a > \sigma_{th}$ a region expanding
from the center of the system with expansion speed $V$ (hereafter referred to as ‘shock speed’); see Fig. 2(a). Hereafter $\sigma_a$ stands as the plateau value of the function $\sigma_a(x, t)$. The stress varies spatially from 0 to $\sigma_a$ over a small width constant over time. Thus, for simplicity, our ‘shocks’ are merely moving stress steps. For investigation purposes, we assume $\sigma_a$ and $V$ to be independent parameters. In shock physics, however, these parameters are connected via the pressure/shock speed relation determined by the Rayleigh line on the Hugoniot [2, 25]. The parametric domain explored is $1.05 \le \sigma_a/\sigma_{th} \le 2.0$ (11 values) and $0.1 \le V/c_s \le 3.0$ (31 values), using $c_L = 2c_s$. If $V$ is moderate with respect to wavespeeds, the loading operates over an expanding region as does a nanoindentation process in pressure [18, 59–61]; if instead $V > c_L$ for edges (resp. $c_s$ for screws), it mimicks a shock of speed exceeding the upper wavespeed relevant to the dislocation character considered. In shock physics the regime $V > c_L$ where the shock overruns elastic waves is often referred to as *overdriven* [14, 62, 63]. We emphasize that physical shocks possess inherent nonlinear characteristics that are beyond the scope of our ‘shock’ model, as further discussed in the concluding Sec. IV. For the time being, we will omit the quotation marks around the term ‘shock’, although it is crucial to remain aware of the simplifications made with regard to shock physics.

Simulations have been carried out on systems of size $L = (160\pi)d$ for screw dislocations, and $L = (320\pi)d$ for edge dislocations, discretized into 4096 spatial points, with periodic boundary conditions (PBCs), and inner time step $\delta t = 3 \cdot 10^{-2}\tau_0$ further subdivided by the algorithm when necessary. Derivatives $\partial\eta/\partial x$ and $\partial\eta/\partial t$ are estimated by finite differences using the inner time step. For each run, the various components of the DPE (self-stress, applied stress, viscous term) are output for 1000 equidistant times during the run, which results in an effective time step $\Delta t \gg \delta t$. These data are post-processed to extract dislocation positions and velocities. At each time step, the DPE is satisfied numerically with maximal absolute error of order $10^{-6}$ [37]. Due to PBCs, $\sigma_a(x, t)$ becomes uniform and steady once boundaries of the simulation box have been reached. Initial conditions consist of zero or one pre-existing dislocation at rest.

Typical evolutions of $\eta(x, t)$ are represented in Fig. 3. Hills at large times on both sides of the figures are due to collisions with reentrant dislocations. The asymmetry of the single-dislocation initial state in (a) is enhanced at large times through nucleation events.

## B. Results

### 1. Trajectories and overall system dynamics

Accounting for spatial discretization over $x$ of $\eta(x, t)$, individual dislocation trajectories [64] $\xi_i(t)$ ($i$ is a dislocation index) are identified by linear interpolation between discretization points as values $\xi_i$ of $\xi$ for which $|\eta(\xi, t)/b|$ is half-integer. Figure 4 displays trajectories for $\sigma_a = 1.2\sigma_{th}$, and increasing values of $V$, for an initial state with one dislocation

<figure><img src="image_4.png" /><figcaption>Figure 3. Edge dislocations. Evolution of $\eta(x, t)$ for $\sigma_a = 1.2\sigma_{th}$ (system size $L = 160\pi$). (a) one initial dislocation at rest; (b) void initial state. ‘Shock’ speeds $V$ as indicated.</figcaption></figure>

at rest (the smaller the slope of the trajectory, the larger the dislocation speed).

The plots show that the dynamics starts off by successive nucleations of pairs near the center. This transient step is soon followed by nucleation by bursts in the bulk, with strong apparent randomness in positions and instants. It is emphasized that through the function $\eta(x, t)$, the Peierls model naturally embodies the variability of the local Burgers vector at nucleation [65]. Nucleations occur where the applied stress is maximal at the beginning of simulations, and afterwards preferentially near to pre-existing moving dislocations in reactions such as in Fig. 5 below, as was previously reported in a Frenkel-Kontorova model [21]. Although the distance is much variable, this feature suggests a heterogeneous nucleation process, possibly involving nucleation avalanches [22] (but the system is too small to conclude on the latter point). Note, however, that the terminology *heterogeneous nucleation* [18, 66, 67] is usually reserved to nucleation near to pre-existing defects, other than previously-nucleated *moving* dislocations.

The leading dislocations cross the box’s boundaries via PBCs at time $t = t_{BC}$. Afterwards, the system is expected to enter a later steady regime (in a statistical sense) that has not been explored. In order to eliminate initial transients and the effect of PBCs, mean quantities will be extracted hereafter from data restricted to the time window $t \in T_{av} = [(2/3)t_{BC}, t_{BC}]$. A strong hypothesis of the forthcoming measurements of mean dislocation density and speeds is that within this window, the system achieves a state reasonably close to some steady-state. This might however be in-
<figure><img src="image_5.png" /><figcaption>Figure 4. (Color online) Edge dislocation trajectories for $\max_x \sigma_a(x, t) = 1.2 \sigma_{th}$, and various ‘shock’ speeds $V$ as indicated, with one dislocation in the initial state. Negative (resp., positive) dislocations are in blue (resp., orange).</figcaption></figure>

valid due to finite-size and time effects. This potential limitation could be overcome by considering larger system sizes and simulation times.

## 2. Reactions

The various reactions observed in Fig. 4 are categorized in Fig. 5. The quasi-elastic reaction #1 is observed at small $V$. Reactions #2 to #6 are the most frequent ones. Contrary to

the frequent isolated nucleation #2, annihilations are only observed in close connection with previous crossing events, #7. However, the annihilation-scattering process at the center of #8 (a time-reversed version of #5) has not been observed individually. The apparent ‘pre-acceleration’ of #5 is explained by shape changes in the initial stages of pair nucleation.

Dynamic collisions and crossings of dislocations of opposite signs and speeds occur via pair annihilation upon collision, and subsequent re-nucleation [68]. After annihilation, the available energy locally transferred to the elastic displacement
<figure><img src="image_6.png" /><figcaption>Figure 5. Dynamical reactions observed in calculations (blow-ups on trajectories such as in Fig. 4, with same orientation and color coding).</figcaption></figure>

ment field flows away with elastic wavespeed from the region of collision. In absence of applied stress the crossing can be completed by subsequent re-nucleation only if the dislocations possess enough initial momentum to re-create the pair before energy has escaped [68]. Here, with $\sigma_a > \sigma_{th}$, re-nucleation just after collision is almost always guaranteed unless some perturbation reduces the local stress level for long enough a time to prevent it. This composite sequence (typically, #3) may look quasi-instantaneous as in #4, without appreciable differences between the ‘in’ and ‘out’ dislocations speeds (here, $v \approx c_L$ for both).

### 3. Dislocation speeds

<figure><img src="image_7.png" /><figcaption>Figure 6. Speed $v_{lead}$ of leading dislocations versus ‘shock’ speed $V$ (method M1), for various $\sigma_a$. Inset in b): eye-estimated stress-dependent shock-speed threshold $V_c(\sigma_a)$ for transition from subsonic to intersonic branch, versus stress $\sigma_a$, with fit $V_c/c_s \approx 1.70 - 0.62(\sigma_a/\sigma_{th})$.</figcaption></figure>

Once positions $\xi_i(t)$ have been determined as explained above, individual dislocation speeds $v_i(t)$ are estimated either as $\Delta\xi_i(t)/\Delta t$ (method M1), or by Eqs. (6) and (7) as $v_i(t) \approx -\partial_t\eta(\xi_i(t), t)/\partial_x\eta(\xi_i(t), t)$ (method M2). Method M1 requires following the trajectories by continuity to perform the numerical time differentiation, which is problematic in presence of reactions such as above, except for isolated dislocations. Method M2 is free of the latter drawback, but can

be inaccurate due to the unavoidable time variations of the dislocation cores embodied in its defining formula; see Eq. (8).

Of special importance is the speed $v_{lead}$ of the leading (i.e., outermost) dislocations evoked in Sec. III B 1 above. Speed $v_{lead}$, estimated with method M1 (the trajectory being unambiguous here), is displayed in Fig. 6 versus $V$ for screw and edge dislocations, for a range of values of $\sigma_a/\sigma_{th}$. Speeds have been time-averaged over the interval $T_{av}$. The outlier in plot b) is an artifact. The plot shows that $v_{lead} \ge V$ whatever $V$, with $v_{lead} = V$ for $V \ge c_s$ (resp. $c_L$) for screw (resp. edge) dislocations. Therefore, dislocations are always present ahead of or on the ‘shock’ front. The observed low-stress dependence of the data for $v_{lead} < c_s$ is attributed to the dislocations not having reached their asymptotic subsonic velocity within the simulation window, due to slow relaxation [36]. Disregarding this, the data suggest that $v_{lead}$ is quasi-independent of $\sigma_a$, except for remarkable transitions of the edge dislocation from a subsonic branch of asymptotic states (with speed close to the Rayleigh velocity $c_R$) to an upper intersonic stable branch. This transition takes place at a stress-dependent shock-speed threshold $V_c(\sigma_a)$ that decreases linearly as $\sigma_a$ increases (see inset), which reveals a two-branch structure of $v_{lead}(V)$ akin to that of the stress/velocity curve of an edge dislocation (Fig. 1). Both are most certainly related. In fact, the simulations reported in Fig. 4 suggest that a two-step loading process (see Sec. II A) takes place for $V < c_L$, as the leading edge dislocation is nucleated within the loaded zone with local stress $\sigma \approx \sigma_{th}$, but is afterwards expelled out of the latter into the unloaded zone ahead where it only experiences the weaker long-range stress of the train of dislocations behind.

Remarkably, although the present model admits no steady-state mobility law of a single edge dislocation for $\sigma_a > \sigma_{th}$ and only allows for speeds restricted to the domain $v < c_L$, dynamically-nucleated dislocations can nonetheless reach steady speeds $v > c_L$. Such speeds are not determined by $\sigma_a$ but instead by the ‘shock’ speed. So to say, leading dislocations are ‘pushed’ forward by the dislocations continuously nucleated behind to accomodate shock advance.

<figure><img src="image_8.png" /><figcaption>Figure 7. Definition of the *Bulk* and *Front* averaging zones in $(x, t)$ representation for an edge dislocation and for: (a) ‘shock’ speed $V < c_L$; (b) shock-speed $V > c_L$. In both cases the stressed zone is $|x| < Vt$. For a screw dislocation, $c_L$ should be replaced by $c_s$.</figcaption></figure>

These differences of nature between the regimes $V \le c_L$ prompt us to introduce the notions of *bulk* and *front* zones, as depicted in Figs. 7 (a) and (b), of time-dependent widths
$\Delta_{\text{bulk}} = \min(V, c_L)t$ and $\Delta_{\text{front}} = \Delta - \Delta_{\text{bulk}}$, where $\Delta(t) = v_{\text{lead}}t$ is the total extension of the process zone for $x > 0$. The case $V > c_L$ figures an overdriven shock [14].

<figure><img src="image_9.png" /><figcaption>Figure 8. (Color online) Screw dislocations. Mean dislocation speed $\bar{v}$, Eq. (10), vs. shock speed $V$ for various $\sigma_a$ (data sets in different colors) obtained from averages over a) the front zone and b) the bulk zone.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. (Color online) Edge dislocations. Mean dislocation speed $\bar{v}$, Eq. (10), vs. shock speed $V$ for various $\sigma_a$ (data sets in different colors) obtained from averages over a) the front zone and b) the bulk zone.</figcaption></figure>

The steady-state speed of the ‘bulk’ dislocations after completion of the nucleation event (but not necessarily at larger times, due to reactions) can be understood from the following non-rigorous but intuitive argument. The absolute value of the local stress, $|\sigma(x, t)|$ is less than $\sigma_{\text{th}}$, but approaches that value in the majority of the bulk region, as Fig. 11 below demonstrates. Then, the driving stress experienced by each dislocation of a newly nucleated pair is presumably close to $\sigma_{\text{th}}$, too. It can thus be substituted to $\sigma_a$ in the stess-velocity law of Fig. 1. Since for edges it is obviously larger than the branch-selection threshold $\sigma_c$, the velocity branch selected is the upper one, with speed $c_L$, which is the maximal admissible one at $\sigma_a = \sigma_{\text{th}}$. Likewise, the mobility law of a screw dislocation is single-branched, with speed $c_S$ at $\sigma_a = \sigma_{\text{th}}$, which determines the nucleation speed of screws. The argument is made complete by considering that supersonic dislocations within each pair cannot interact via wave propagation, and can therefore be considered independent. Note that steady glide at the shear wavespeed of screw dislocations in dipole pairs was observed in atomistic simulations [69].

The dislocation speed $\bar{v}$ averaged separately over the front and bulk zones is represented versus $V$ in Fig. 8 (resp., 9) for

screw (resp., edge) dislocations. It is computed as

$$ \bar{v}_z = \left\langle \frac{1}{N_z(t)} \sum_{i \in z} |v_i(t)| \right\rangle_{t \in T_{\text{av}}} \quad (10) $$

where the $v_i$ are estimated by method M2, the zone label $z$ is bulk or front, and $N_z(t)$ is the number of dislocations in the zone. In subplots (a) front dislocations move with mean speed increasing with shock speed $V$, dominated by a leading-order proportionality to the latter. Figure 9(a) display weak overall dependence with respect to $\sigma_a$, which is absent for screw dislocations in Fig. 8(a), and might be a finite-size artifact. By contrast, the mean speed in the bulk in Fig. 9(b) is markedly stress-dependent. This has two different causes: for $V < c_S$, the lowering of the mean speed at high stress is presumably due to the contribution of bulk dislocations nucleated close to $v = c_L$ that become subsonic after having ‘bounced’ on the zone boundary. This effect, specific to edge dislocations, becomes majoritary in the range $c_R < V < c_S$ whatever $\sigma_a$, and is absent for screw dislocations in Fig. 8(b). On the other hand, the high nucleation rate for $V > c_L$ increases the number of transient speeds $v \gg c_L$ sampled close to nucleation. It is present for both edge and screw dislocations.

Typical multimodal speed distributions for $\sigma_a = 1.2\sigma_{\text{th}}$, normalized to unity, collected over the whole simulation duration and over both regions, are displayed in Fig. 10 for edge dislocations. A peak close to $v = c_L = 2.0c_S$ is always present, which corresponds to the ‘bulk’ dislocations at their nucleation speed. By contrast, no large-amplitude peak is associated with the shock speed $V$ until $V > c_L$. Moreover, the system becomes populated with an appreciable amount of intersonic dislocations ($c_S < v < c_L$) as soon as $V$ exceeds the transition threshold $V_c(\sigma_a = 1.2\sigma_{\text{th}}) \approx 0.96c_S$ (see inset of Fig. 6). These intersonic dislocations form a well-defined peak if $c_S < V < c_L$. Quite strikingly, the system stands free of subsonic dislocations ($v < c_S$) if $V > 1.3c_S$: then, plastic deformation essentially consist of intersonic/supersonic dislocations. Finally, even for $V < c_L$ distributions display a small but appreciable amount of speeds $v \gg c_L$. Indeed, newly-nucleated dislocations, until fully formed, induce huge apparent measured transient speeds, as can be seen from trajectories in Fig. 4.

### 4. Nucleation and local stress

As recalled in Sec. II B the athermal homogeneous nucleation threshold of a dipole identifies with the theoretical shear stress. Nucleation relieves the local excess stress in-between the two components of the nucleated dipole [13, 70]. As these components move away from one another (which would correspond to loop expansion in higher dimensions) this screening becomes less effective, until at some separation distance a new pair pops up to reinstate locally a condition $|\sigma|/\sigma_{\text{th}} \lesssim 1$.

Figure 11 displays, for subsonic $V$, a typical profile of $\sigma(x, t)$ at some intermediate time, together with the corresponding nucleated dislocation density (only half the system is shown). The local density $\rho$ being defined by Eq. (6)$_1$, the
<figure><img src="image_11.png" /><figcaption>Figure 10. (Color online) Edge dislocations. Speed distributions $P(v)$ for $\sigma_a/\sigma_{th} = 1.2$, and shock speeds $V$ (dashed, red) as indicated</figcaption></figure>

insets focused on the bulk/front boundary bewteen show that the stress and the dislocation cores are well resolved in our calculations. They also fluctuate, and stress heterogeneity [17] is different ahead of and behind the shock front. The main plot evidences a forward stress tail. It arises in the form of a wave emitted by the first nucleated dislocation, and is afterwards reinforced in a cumulative way as more dislocations are nucleated behind. Its own front (unnoticeable on the figure) propagates with speed $c_L$. At the time represented, it is quasi-complete and its maximal amplitude does not change much thereafter. Similar pulses were observed in atomistic simulations [14]. The strong shape variations of dislocation cores under the action of the surrounding wave system explains the difficulty of accurately estimating dislocation velocities with method M2 in the previous Section.

Although the exact triggering events for nucleation could not be pinpointed, obvious candidates are fluctuations in the bath of elastic waves emitted and scattered by dislocations. The role of waves in 'randomly' (i.e., pseudo-randomly) triggering nucleations is hinted at by considering cases $V > c_L$. Then the system entirely consists of supersonic dislocations with $v > c_L$ (see Fig. 9). Elastic waves from these supersonic dislocations are radiated as Mach cones out of the slip

<figure><img src="image_12.png" /><figcaption>Figure 11. Edge dislocations. Main plot: at $t/\tau_0 = 39.93$, local stress $\sigma(x, t)$ for $\sigma_a/\sigma_{th} = 1.5$ and $V = 0.30 c_S$ (black), superimposed with simultaneous $\sigma_a(x, t)$ (orange) and dislocation density $\rho(x, t)$ (green) – divided by 50 for better representation. Top inset: blow-up of same dislocation density for $0 \le x/d \le 20 \times 2\pi$. Bottom inset: blow-up of same $\sigma(x, t)$ within this range.</figcaption></figure>

plane without any possibility of disturbing $\eta$. This makes interactions bewteen edge dislocations essentially local, so that a decrease in the number of 'random' nucleation events is expected. This is precisely what Fig. 4(i) illustrates at the highest shock speed $V = 3 c_S$ where a regular array of supersonic dislocations [27] is produced, to be investigated further below and in Appendix C. Thus, for supersonic dislocations at very high shock speeds, nucleation turns out fully kinematic [19].

### 5. Dislocation spacings and densities

Dislocation densities can be understood in terms of dislocation mean spacings. The overdriven case $V \ge c_L$ such as in Figs. 4(g,h,i) is approximately addressed remarking that the front zone consists of a regularly-spaced array of supersonic dislocations of same sign moving with speed $v \approx V$. Based on the Weertman equation, a model of dislocation spacing $\Delta x(v, \tau)$ in the array is derived in Appendix C, leading to the dislocation spacing in the front zone

$$ \Delta x_f(V, \tau) = (2\pi) d B_\alpha(V) / \sqrt{\tau^2 - 1} \quad (V \ge c_L), \quad (11) $$

where $B_\alpha(v)$ is the drag coefficient at supersonic speed $v$ and $\tau = \sigma_a/\sigma_{th}$ is the reduced stress.

Figure 4 suggests that for $V \gg c_L$ dislocation pairs are almost exclusively nucleated at the bulk/front boundary, with one dislocation going into the bulk zone and the other one into the front zone. Under this assumption, the nucleation rate (number of dislocations nucleated per unit time) for each zone is the relative speed $\Delta v$ at which dislocations move away from the boundary, divided by the mean dislocation spacing. This speed being $\Delta v = 2c_L$ in the bulk zone, and $\Delta v = V - c_L$ in the front zone, equating nucleation rates $\Delta v/\Delta x$ gives the following estimate for the mean spacing for dislocations of same sign in the bulk zone as

$$ \Delta x_b(V, \tau) \approx \frac{2c_L}{V - c_L} \Delta x_f(V, \tau) \quad (V \gg c_L). \quad (12) $$
<figure><img src="image_13.png" /><figcaption>Figure 12. For $V \ge c_L$, dislocation spacings vs. $V$, for reduced stresses $\tau = 1.05, 1.1, 1.2, \dots, 2.0$ from top (blue) to bottom (red), in the front zone (a) and the bulk zone (b, in log. scale). Markers: measurements. Dashed lines: Eqs. (11) in (a) and (12) in (b).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13. (Color online) Edge dislocations. (a) Mean density $\bar{\rho}$ in the front zone vs. $V$ for various stresses and (b) data collapse with rescaling exponent $\beta = 1/2$.</figcaption></figure>

In Fig. 12 measurements of spacings between dislocations of same sign, averaged over over a few outermost dislocations in the front zone, and over a few dislocations in the bulk, are compared to Eqs. (11) and (12). Figure (a) shows excellent overall agreement with (11), in spite of deviations at small $V$ for the lowest $\sigma_a$s, which indicate that the steady-state has not been reached there. In (b) by constrast, Eq. (12) gives a correct order of magnitude only when $V \gg c_L$, due to the crude nature of our argument (data relative to parameter values for which a relevant spacing could not be determined due to excessive reaction-induced perturbations have been omitted in the plots). Indeed, an ordered ‘steady state’ in the bulk for $V > c_L$ consists of two arrays of dislocation of opposite signs steadily moving in opposite directions, which periodically annihilate one another (e.g., Fig. 4 (i)). Such a pattern is intrinsically dynamic and escapes the Weertman equation used to derive (11).

Using Eq. (6), the mean (unsigned) dislocation density $\bar{\rho}$ is estimated via measurements in bulk and front zones as

$$ \bar{\rho}_z = \left\langle \int_z dx \left| \rho(x, t) \right| / \Delta_z(t) \right\rangle_{T_{av}} \quad (13) $$

where the zone label $z$ is front or bulk. The integral is evaluated as a Riemann sum. A further average is then made over time within the useful time window.

For $V > c_L$, Eqs. (C6) and (C8) in Appendix C suggest a steady-state mean density in the front zone

$$ \bar{\rho}_f(V, \tau) = (b\Delta x_f)^{-1} \propto \sqrt{\tau^2 - 1}. \quad (14) $$

This estimate does not compare quantitatively well with the data and is not considered further. This may be due to various causes among which dislocation core shapes varying much strongly than in the theoretical model (see Fig. 9).

However, our data are consistent with a generalized scaling factor in stress of the form $S_\beta(\tau) = (\tau^2 - 1)^\beta$ with $0 < \beta < 1$ an empirical complexity-related zone-dependent rescaling exponent. Figure 13(a) displays measured densities in the front zone for edge dislocations. Apart from deviations attributed to finite-size effects, the rescaling in (b) is consistent with an exponent $\beta = 1/2$ over the whole range of $V$. The

transitions for $V < c_S$ in Fig. (b) should be brought together with those in Fig. 6. Mean density data in the bulk zone are

<figure><img src="image_15.png" /><figcaption>Figure 14. (Color online) Edge dislocations. (a) Mean density $\bar{\rho}$ in the bulk zone vs. $V$ for various stresses and (b) data collapse with different rescaling exponents $\beta = 3/4$ in region A, and $\beta = 1/2$ in region B.</figcaption></figure>

reported in Figure 14(a), while in (b) exponents were needed for $S(\tau)$ in regions A ($V < c_S$) and B ($V > c_S$) as indicated in the caption. See Supplemental Material [40] for data similar to Figs. 13 and 14 for screw dislocations.

## IV. CONCLUDING DISCUSSION

Within a 1D model, we investigated high-stress multiple nucleations of dislocation pairs on a single slip plane under (elasto)dynamic conditions by means of an original method. Although the framework is free of side complexities such as Peierls stress or forest hardening [71], the outcome turns out surprisingly rich.

With regard to dislocation-based constitutive models, the study suggests that the notion of a local dislocation-density nucleation rate (such as usually considered in models) may not straightforwardly apply to an advancing shock. Indeed, in absence of forest interactions, dislocations are continuously and kinematically nucleated so as to maintain a quasi-steady density of geometrically-necessary dislocations [72] that screens out the (macroscopic) applied stress. This is much alike the findings of Zhakhovsky et al. in atomistic simulations [14],
and of Gurrutxaga-Lerma et al. in subsonic discrete dislocation elastodynamic simulations [17]. Of relevance as well is the observation that whatever this applied stress, the inner stress that drives dislocations and determines their speed via the mobility law should never exceed the nucleation threshold. A distinction between the applied stress and this local stress should therefore presumably be included in models.

Applying the model to a highly idealized 'shock' situation, we found the system to spontaneously separate into two 'bulk' and 'front' zones. Whereas Meyers [26], opposing Smith's views [27], considered dislocation generation rather than motion responsible for the accommodation of lattice deviatoric strains, the present work suggests that both generation and motion are equally important in shocks: here dislocations are generated only either in the bulk zone or at the interface between the bulk and front zone, at supersonic speed, and can be expelled into the nucleation-free front zone. In particular, supersonic 'shock' speeds cause supersonic dislocations to transit from the bulk to the front zone in this model.

An accurate theory is presently lacking to account for the densities reported in Figs. 13 and 14 over the whole domain of parameters. However, dislocation densities were found to scale as $((\sigma_a/\sigma_{th})^2 - 1)^\beta$, where $\sigma_a$ is the applied stress, $\sigma_{th}$ is the theoretical shear stress (nucleation stress), and $0 < \beta < 1$. Exponent $\beta$ does not depend on the shock speed $V$ in the front zone, where it is of order $\beta \sim 0.5$. In the bulk zone, it differs depending on whether $V < c_S$ ($\beta \sim 0.75$ for edges, and $\sim 0.66$ for screws) or $V > c_S$ ($\beta \sim 0.5$ for both characters). These values should however be reexamined on larger systems for better statistics. Moreover, they may depend on the $\gamma$-potential, since the base form $[(\sigma_a/\sigma_{th})^2 - 1]^{1/2}$ (Appendix C) is very much related to the sine Frenkel form. A scaling of the type $[(\sigma_a/\sigma_{th})^{\beta_1} - 1]^{\beta_2}$ with two variable exponents might however hold more generally. Because the PND equation includes no non-trivial equation of state (EoS), in spite of some similarities the present two-zone picture differs from the one discovered by Zhakhovsky et al. in the overdriven regime [14]. One major difference is that our front zone is continuously expanding, but this can be due to 1D steric effects and may not occur in higher dimensions.

Our model's approximations may affect the results in at least three aspects. First, our framework assumes linear elasticity, except for the pull-back force. However, real shocks involve nonlinear elasticity due to the non-trivial EoS. This would result in pressure-dependent elastic constants and higher wave speeds behind the shock front compared to ahead of it. However, in compressed metals, Poisson's ratio remains approximately constant in the absence of phase transformation. Consequently, the wavespeed ratio $\gamma = c_L/c_S$ in the kernels presented in Appendix A would remain unchanged by compression [2]. Hence, uniformly compressed stress steps should yield similar outcomes within the loaded region (encompassing both the "bulk" and "front" zones) with same dislocation reactions. On the other hand, if the wavespeeds in the loaded "bulk" zone are higher than those in the unloaded "front" zone for the case of $V < c_L$, the self-organization of

the front might change.

Second, the underlying physics of the DPE may not be complete. Introducing local inertia in the pull-back force, as discussed in Section II B, would likely alter the dynamics and introduce additional delays in nucleation and reactions.

Third, in our 1D simulations, all dislocations turn out mobile. Overcoming this limitation would require expanding the analysis to higher dimensions to account for junctions and other reactions [17, 73]. These three issues require detailed investigations that extend beyond the scope of this paper.

Nevertheless, in spite of these possible limitations, the present study embodies full inertia of radiative origin [36, 74], drag, wave propagation and retarded interactions on the slip plane, nucleation and annihilations, and supersonic motion, within one self-contained field equation without the need for extra inputs. Therefore, it might serve as a one-dimensional, reference point for linear-elastic, discrete elastodynamic simulations, as well as a tool to evaluate future nucleation-related equations for dislocation density evolution.

## ACKNOWLEDGMENTS

Y.P.P. thanks his colleague C. Denoual for discussions and a remark about reactions.

## Appendix A: Kernels and their transforms

The various kernels needed are as follows (the climb-edge [75] case is reported for completeness only, as it is not used in our simulations). Kernels $K_a(x, t)$ [35, 74] vanish for $t < 0$ and are given below for $t > 0$ together with their associated constants $\kappa_a$. Their spatial Fourier transforms (FT) $K_a(k, t)$, of wavemode $k$, are expressed in terms of functions $C_\alpha(\tau)$, which also vanish for $t < 0$, as

$$K_a(k, t) = -i \frac{\pi}{2} c_S k C_a(c_S k t). \quad (A1)$$

Bearing in mind the usual correspondence between dislocation characters and fracture modes, the $C_a(\tau)$ are reproduced hereafter from Refs. [76] in terms of the wavespeed ratio $\gamma = c_L/c_S > 1$, and of the integral

$$W(\tau) = 1 - \int_{0}^{\tau} \frac{dq}{q} J_{1}(q), \quad (A2)$$

which admits a closed-form expression [77]. Finally, the numerical method of solution (see Appendix B) needs the Laplace transforms (LT)

$$C_a(s) = \int_{0}^{+\infty} d\tau C_a(\tau) e^{-s\tau}, \quad (A3)$$

which we give here in a slightly reorganized form with respect to [76] to prepare for Eq. (A5) below.
Thus for $t > 0$,

$$ \text{screw :} \qquad K_s(x, t) = \frac{1}{2c_s} \frac{x}{t^2} (c_s^2 t^2 - x^2)_+^{-1/2} \quad \text{with} \quad \kappa_s = 1, \qquad \text{(A4a)} $$

$$ C_s(\tau) = \frac{J_1(\tau)}{\tau}, \qquad \text{(A4b)} $$

$$ C_s(s) = \sqrt{1 + s^2} - s, \qquad \text{(A4c)} $$

$$ \text{glide edge :} \qquad K_g(x, t) = \frac{2c_s^2}{x^3} \left[ \frac{1}{c_L} (2c_L^2 t^2 - x^2)(c_L^2 t^2 - x^2)_+^{-1/2} - \frac{1}{c_s} (2c_s^2 t^2 - x^2)(c_s^2 t^2 - x^2)_+^{-1/2} \right] \\ \qquad \qquad + \frac{x}{2c_s t^2} (c_s^2 t^2 - x^2)_+^{-1/2} + \frac{c_s}{2} x \operatorname{Pf}(c_s^2 t^2 - x^2)_+^{-3/2} \quad \text{with} \quad \kappa_g = 1, \qquad \text{(A4d)} $$

$$ C_g(\tau) = \frac{J_1(\tau)}{\tau} + 4\tau [W(\gamma\tau) - W(\tau)] - \frac{4}{\gamma} J_0(\gamma\tau) + 3J_0(\tau), \qquad \text{(A4e)} $$

$$ C_g(s) = -\frac{4}{s^2} \left[ \sqrt{1 + (s/\gamma)^2} - \frac{(1 + s^2/2)^2}{\sqrt{1 + s^2}} \right] - s, \qquad \text{(A4f)} $$

$$ \text{climb edge :} \qquad K_c(x, t) = -\frac{2c_s^2}{x^3} \left[ \frac{1}{c_L} (2c_L^2 t^2 - x^2)(c_L^2 t^2 - x^2)_+^{-1/2} - \frac{1}{c_s} (2c_s^2 t^2 - x^2)(c_s^2 t^2 - x^2)_+^{-1/2} \right] \\ \qquad \qquad + \frac{c_L}{2c_s^2} \frac{x}{t^2} (c_L^2 t^2 - x^2)_+^{-1/2} + \frac{c_s^2}{2c_L} \left( \frac{c_L^2}{c_s^2} - 2 \right)^2 x \operatorname{Pf}(c_L^2 t^2 - x^2)_+^{-3/2} \quad \text{with} \quad \kappa_c = \gamma, \qquad \text{(A4g)} $$

$$ C_c(\tau) = \gamma^3 \frac{J_1(\gamma\tau)}{\gamma\tau} + 4\tau [W(\tau) - W(\gamma\tau)] + (4\gamma - \gamma^3)J_0(\gamma\tau) - 4J_0(\tau), \qquad \text{(A4h)} $$

$$ C_c(s) = -\frac{4}{s^2} \left[ \sqrt{1 + s^2} - \frac{(1 + s^2/2)^2}{\sqrt{1 + (s/\gamma)^2}} \right] - \gamma s, \qquad \text{(A4i)} $$

where Pf denotes Hadamard's finite part and $(x)_+^\alpha = x^\alpha$ if $x > 0$ and $0$ otherwise [78]. Equations (A4a) and (A4d) are, respectively, equations (7a) and (7b) of [74], written differently. Equation (A4g) derives from Eqs. (43a) and (43b) of [35] and the transformation explained in [74]. Letting $s = iv/c_s$, both terms within brackets in (A4f) and (A4i) are proportional to the Rayleigh function. Their common root $v^* = c_R$ is the Rayleigh wavespeed of free-surface waves [79, p. 18], equal to $c_R \approx 0.9325 c_s$ if $c_L = 2 c_s$ as in the present work.

Quite remarkably, expressions (A4c), (A4f), and (A4i) show Laplace transforms $C_a(s)$ to be directly connected to the so-called prelogarithmic steady-state Lagrangian functions [74, 80] $L_a(v)$ for each character $a$, with $v$ the dislocation velocity, via the relationship $(w_0 = \mu b^2/(4\pi))$ is the characteristic line energy density)

$$ C_a(s) = -\frac{1}{w_0} L(i c_s s) - \kappa_a s. \qquad \text{(A5)} $$

This connection went previously unnoticed to the best of our knowledge. Moreover, the constant $\kappa_a$ can be obtained in all cases from $L(v)$ as

$$ \kappa_a = -\frac{1}{w_0} \lim_{s \to \infty} \frac{1}{s} L(i c_s s), \qquad \text{(A6)} $$

which is the same limit as Eq. (66) of [74], formulated otherwise.

## Appendix B: Numerical method

The main difficulty of (5a) lies in the integrodifferential convolution operator, which is non-local in time and space. Nevertheless, this operator is diagonalized by the Fourier transform in space, followed by a Laplace transform in time. Our method crucially relies on this property. Although used in a different way (notably, with respect to the time convolution) those transforms are central as well in the approach of Rice and co-workers in the elastodynamic crack problem, which moreover involves the same elastodynamic kernels [76]; see Appendix A.

In the rest of the Section, the substitution $c_s t \to t$ is implied to simplify notations, so that time has dimension of space. Using Eq. (A1) the DPE reads in Fourier form

$$ \kappa_a^\alpha \partial_t \delta\eta(k, t) + k^2 \int_{-\infty}^{t} dt' C_a(|k|(t-t')) \delta\eta(k, t') = -(2/\mu)F(k, t). \qquad \text{(B1)} $$

The kernels involved are stiff in space and time, because of their proportionality to $k^2$, and because the kernel $C_a(|k|(t-t'))$ oscillates faster and faster as $k \to \infty$. The simulation interval of the $x$-axis is discretized with step $\Delta x$ and the continuous FT is approximated by the Discrete FT using FFT routines. For simplicity of exposition, the discrete character of the Fourier modes is left implicit hereafter.
In statics, the fast Fourier transform (FFT) technique has been employed repeatedly to solve the PN equation [48, 81]. Since here all functions are smooth, Gibbs oscillations are not a problem and no filtering or modification of the wavemodes is required. However, use of the DFT makes $\delta\eta$ periodic in space. No attempt has been made to mitigate this drawback, which could partially be overcome using zero-padding [37, 82]; see also [76].

Then, time integration requires first re-expressing (B1) as

$$\delta\eta(k, t) = -\frac{2}{\mu} \int_{0}^{t} dt' R_{a}(|k|(t - t')) F(k, t'), \quad (B2)$$

where the resolvent $R_a$ appears as the Fourier transform of the fundamental solution $g_a$ — more precisely $g_a(k, t) = -(2/\mu)R_a(|k|t)$. The latter is the solution of the instance of (5a) associated with a point source $F(x, t) \to F^{\text{point}}(x, t) = \delta(t)\delta(x)$. The advantage of (B2) over (B1) is that the stability properties of the convolution operator, which are hidden in the oscillations of $C_a$, appear more explicitly as a long-range decay of $R_a$ [37] (this has beneficial numerical consequences.)

Then, for time integration, we make use of Lubich and Schädle's so-called fast and oblivious method [83], where the adjective oblivious actually indicates that the information storage of the past times is 'parsimonious'. This method proved the most efficient one in terms of accuracy, speed, and memory requirements, among the wide variety of strategies thoroughly tested and compared for the present problem by the second author (M.J.) [37]. The main ideas of the method, which rests on a numerical inversion of the Laplace transform, are as follows.

Eq. (B2) is solved by a LT as

$$\delta\eta(k, s) = -\frac{2}{\mu} \frac{1}{|k|} \mathcal{R}_a(s/|k|) \mathcal{F}(k, s), \quad (B3)$$

where $\mathcal{F}(k, s)$ is the LT of $F(k, t)$, and where

$$\mathcal{R}_a(s) = \frac{1}{\kappa_a^\alpha s + C_a(s)} \quad (B4)$$

is the LT of the resolvent

$$R_a(\tau) = \int_{\Gamma_B} \frac{ds}{2i\pi} \mathcal{R}_a(s) e^{s\tau}, \quad (B5)$$

where $\Gamma_B$ is a Bromwich inversion contour. To evaluate (B2), the continuous integral in the Laplace inversion of $\mathcal{R}_a$ must be done numerically in discretized form. Knowledge of the poles and cuts in the complex s-plane of $\mathcal{R}_a$ is necessary. (These functions can be expressed so that their cuts lie horizontally parallel to the real axis.) Following Lubich, Schädle and coworkers, the Bromwich contour $\Gamma_B$ is deformed into a set of non-overlapping Talbot contour(s) $\Gamma_l$ with $l$ an integer index, suitably parametrized to enclose from the right the poles and branch cuts. Being confined within horizontal strips of the complex plane, Talbot contours are well-suited to this inversion problem.

Integral (B5) is then approximated by a Riemann sum of $d$

exponential terms, as

$$R_a(|k|t) \approx \sum_{j=1}^{d} \frac{w_j}{2i\pi} \mathcal{R}_a(s_j) e^{s_j |k|t}, \quad (B6)$$

where $\{s_j\}_{j=1,...,d}$ is a $|k|t$-dependent discrete sample of points on the set of contours, and the $w_j$ are parametrization-dependent weights for trapezoidal integration. We employed the discretization parametrization of Ref. [83], which enjoys spectral approximation error. Typically, the number $d$ is not too large —of order 100-500 for accuracy better than $10^{-8}$ on $R_a$; see Ref. [37, Fig. 8.7]. However, exponentials in (B6) explode for $t$ large when $\Re(s_j) > 0$. To ensure stability quantities $\Re(s_j)|k|t$ are kept under control by making contours tighter around singularities as time grows, whenever $|k|t$ falls into specific intervals of values.

Substituting (B6) into (B2) gives

$$\delta\eta(k, t) \approx \sum_{j=0}^{d} w_j u_j(k, t), \quad (B7a)$$

with for $j = 1, ..., d$,

$$u_j(k, t) \equiv -\frac{2}{\mu} \frac{\mathcal{R}_a(s_j)}{2i\pi} \int_0^t d\tau e^{s_j |k|(t-\tau)} F(k, \tau). \quad (B7b)$$

Further differentiating (B7b) with respect to time then yields evolution equations for the $u_j(k, t)$ as

$$\frac{du_j}{dt}(k, t) = s_j |k| u_j(k, t) - \frac{2}{\mu} \frac{\mathcal{R}_a(s_j)}{2i\pi} F(k, t). \quad (B8)$$

Thus, the integral equation (B1) has been replaced by a simpler system of $d$ first-order evolution equations. This approach avoids keeping in memory the $F(k, \tau)$ for $0 < \tau < t$, in favor of advancing the $u_j(k, t)$ in time for each $s_j$. A hierarchical handling of the past information on the contour-related intervals of $|k|t$ confers the algorithm its excellent scaling properties in time. We refer to the original publications [83] and to Ref. [48] for details on the (delicate) procedure for selecting the suitable contour at each time step. In practice, Eqs. (B8) are discretized in time by means of the so-called Radau II A Runge-Kutta method of order O($\Delta t^5$), with a fixed time step $\Delta t$ common to all modes $k$. Its stability [84] makes this method particularly advantageous to address the stiffness of Eq. (B8). If $T$ is the number of time steps, and $M$ is the number of points in the spatial grid the algorithm has a complexity of O($dM \log M \times T \log T$) and requires a memory of O($dM \log T$), to be compared with O($M^2 T^2$) and O($MT$), respectively, for a naive evaluation of the time convolution in (B1).

To summarize: starting from initial conditions $\eta_0(x)$, and $u_j(k, 0) = 0$ for each discrete wavemode $k$, and given: the expressions of $\mathcal{R}_a(s)$; sets of discretization points $s_j$ on suitable Talbot contours; and a model for $\sigma_a(x, t)$, the algorithm consists in looping over the following tasks for each time step $t_k$: (i) compute the elastic slip $\eta_e(x, t_k)$ associated with $\sigma_a$;
(ii) evaluate $F(x, t_k)$ in direct space as (5b); (iii) going to the Fourier domain, compute $F(k, t_k)$ by FFT; (iv) evolve for each discrete wavemode $k$ the $u_j(k, t)$ to the next time step $t_{k+1}$ by means of (B8); (v) reconstruct $\delta\eta(k, t_{k+1})$ from the $u_j(k, t_{k=1})$ using (B7a), and deduce $\delta\eta(x, t_{k+1})$ by inverse FFT. At $t_k$ the plastic slip is $\eta(x, t_k) = \eta_0(x) + \delta\eta(x, t_k)$.

## Appendix C: Supersonic dislocation array

The following applies to the front region, for $V > c_L$. Once divided by $\sigma_{th}$, and upon setting $\tau = \sigma_a/\sigma_{th}$, Weertman's equation [36, 47] reduces, in the supersonic regime $|v| > c_L$ for edges and $|v| > c_S$ for screws, to

$$2\pi \frac{d}{b} B_{\alpha}(v) \frac{\partial \eta}{\partial x}(x) + \tau = \sin \left( \frac{2\pi}{b} (\eta(x) + \eta_e) \right) \quad (C1)$$

in the co-moving frame of speed $v$. Here, $B_\alpha(v)$ is the radiative drag coefficient augmented by its empirical phonon-drag contribution of dimensionless viscosity parameter $\alpha > 0$, and $\eta_e$ is the background elastic slip [36]. Because of the supersonic regime, the integral term of the full Weertman equation [47] is absent from (C1). Letting $\beta_L = \sqrt{v^2/c_L^2 - 1}$, $\beta_S = \sqrt{v^2/c_S^2 - 1}$, and $s_v = \text{sign } v$, the supersonic radiative drag coefficient reads

$$B^{\text{edge}}(v) = 2s_v \left(\frac{c_S}{v}\right)^2 \left(\beta_L + \frac{[1 - v^2/(2c_S^2)]^2}{\beta_S}\right), \quad (\text{C2a})$$

$$B^{\text{screw}}(v) = \frac{s_v}{2} \beta_S, \quad (\text{C2b})$$

and $B_\alpha(v) = B(v) + \alpha v/c_S$.

In the strong-stress case $|\tau| > 1$, $\eta_e = (b/4)$ sign $\tau$ [36] and Eq. (C1) admits a multiple-dislocation solution akin to a Smith array of dislocations of speed $v$. Eshelby only briefly evoked its existence, describing it as “a linear function of $x - vt$ with a superimposed periodic ripple” [85]. Its mass-

ter analytical form was given by Movchan et al. [86], albeit in a different context that did not require the function $B_\alpha(v)$ (see the latter reference for a plot). It can be built from the following base solution of (C1) indexed by $k \in \mathbb{Z}$:

$$\eta_k^{\text{base}}(x) = -\frac{b}{4} \text{sign} \tau + \frac{b}{\pi} \arctan \left\{ \frac{1}{\tau} - \frac{\sqrt{\tau^2 - 1}}{\tau} \tan \left[ \pi \frac{x}{\Delta x(v, \tau)} - Q(\tau) \right] \right\} - k b, \quad (C3)$$

where the integration constant

$$Q(\tau) = \arctan \sqrt{(|\tau| - 1)/(|\tau| + 1)} \quad (C4)$$

has been adjusted so that $\eta_{k=0}^{\text{base}}(0) = 0$. The function $\eta_k^{\text{base}}$ is bounded, periodic, and discontinuous at points

$$x_k^{\text{disc}} = \Delta x(v, \tau) \left[ \left( k - \frac{1}{2} \right) + \frac{1}{\pi} Q(\tau) \right], \quad (C5)$$

but with continuous derivative. Its spatial period, which corresponds to the distance between dislocations in the array, is

$$\Delta x(v, \tau) = \frac{2\pi|B_{\alpha}(v)|d}{\sqrt{\tau^2 - 1}} \quad (C6)$$

The array solution itself stems from assembling pieces of $\eta_k^{\text{base}}$ in a continuous way, as

$$\eta^{\text{array}}(x) = \eta^{\text{base}}\left[\frac{x - x_0^{\text{disc}}}{\Delta x(v, \tau)}\right](x), \quad (C7)$$

where the value of $k$ in $\eta_k^{\text{base}}$ results from taking an integer part (denoted by $\ll \cdot \gg$), which involves $x_{k=0}^{\text{disc}}$. The mean dislocation density in the array is, as expected,

$$\bar{\rho} = \int_{0}^{\Delta x} \frac{dx}{\Delta x} \left| \frac{\partial \eta^{\text{array}}}{\partial x} \right| = \frac{1}{\Delta x(v, \tau)}. \quad (C8)$$

[1] M. Mogilevsky and P. Newman, Mechanisms of deformation under shock loading, Phys. Rep. 97, 357 (1983).

[2] G. I. Kanel, S. V. Razorenov, and V. E. Fortov, Shock-Wave Phenomena and the Properties of Condensed Matter, 4th ed., High-Pressure Shock Compression of Condensed Matter (Springer, New York, 2004); J. Forbes, Shock Wave Compression of Condensed Matter: A Primer, Shock Wave and High Pressure Phenomena (Springer, 2012).

[3] P. M. Anderson, J. P. Hirth, and J. Lothe, Theory of Dislocations, 3rd ed. (Cambridge University Press, Cambridge, 2017).

[4] R. W. Armstrong and F. J. Zerilli, High rate straining of tantalum and copper, J. Phys. D 43, 492002 (2010).

[5] R. Austin and D. McDowell, Parameterization of a rate-dependent model of shock-induced plasticity for copper, nickel, and aluminum, Int. J. Plast. 32–33, 134 (2012).

[6] N. R. Barton and M. Rhee, A multiscale strength model for tantalum over an extended range of strain rates, J. Appl. Phys.

114, 123507 (2013).

[7] R. Kositski and D. Mordehai, On the origin of the stress spike decay in the elastic precursor in shocked metals, J. Appl. Phys. 126, 085901 (2019); A dislocation-based dynamic strength model for tantalum across a large range of strain rates, J. Appl. Phys. 129, 165108 (2021).

[8] S. Yao, J. Yu, Y. Cui, X. Pei, Y. Yu, and Q. Wu, Revisiting the power law characteristics of the plastic shock front under shock loading, Phys. Rev. Lett. 126, 085503 (2021); S. Yao, X. Pei, J. Yu, and Q. Wu, Assessment of the time-dependent behavior of dislocation multiplication under shock loading, Int. J. Plast. 158, 103434 (2022).

[9] C. Ye, G. Liu, K. Chen, J. Liu, J. Hu, Y. Yu, Y. Mao, and Y. Shen, Unified crystal plasticity model for fcc metals: from quasistatic to shock loading, Phys. Rev. B 107, 024105 (2023).

[10] B. L. Holian and P. S. Lomdahl, Plasticity induced by shock waves in nonequilibrium molecular-dynamics simulations, Sci-
ence 280, 2085 (1998).

[11] T. C. Germann, B. L. Holian, P. S. Lomdahl, and R. Ravelo, Orientation dependence in molecular dynamics simulations of shocked single crystals, *Phys. Rev. Lett.* **84**, 5351 (2000); T. Germann, B. Holian, P. Lomdahl, D. Tanguy, M. Mareschal, and R. Ravelo, Dislocation structure behind a shock front in fcc perfect crystals: Atomistic simulation results, *Metall. Mater. Trans. A* **35**, 2609 (2004).

[12] D. Tanguy, M. Mareschal, P. S. Lomdahl, T. C. Germann, B. L. Holian, and R. Ravelo, Dislocation nucleation induced by a shock wave in a perfect crystal: Molecular dynamics simulations and elastic calculations, *Phys. Rev. B* **68**, 144111 (2003).

[13] M. Meyers, H. Jarmakani, E. Bringa, and B. Remington, Dislocations in shock compression and release, in *Dislocations in Solids*, Vol. 15, edited by J. Hirth and L. Kubin (Elsevier, 2009) Chap. 89, pp. 91–197.

[14] V. V. Zhakhovsky, M. M. Budzevich, N. A. Inogamov, I. I. Oleynik, and C. T. White, Two-zone elastic-plastic single shock waves in solids, *Phys. Rev. Lett.* **107**, 135502 (2011).

[15] M. M. Sichani and D. E. Spearot, A molecular dynamics study of dislocation density generation and plastic relaxation during shock of single crystal Cu, *J. Appl. Phys.* **120**, 045902 (2016).

[16] A. Bisht, A. Neogi, N. Mitra, G. Jagadeesh, and S. Suwas, Investigation of the elastically shock-compressed region and elastic-plastic shock transition in single-crystalline copper to understand the dislocation nucleation mechanism under shock compression, *Shock Waves* **29**, 913 (2019).

[17] B. Gurrutxaga-Lerma, D. S. Balint, D. Dini, D. E. Eakins, and A. P. Sutton, Attenuation of the dynamic yield point of shocked aluminum using elastodynamic simulations of dislocation dynamics, *Phys. Rev. Lett.* **114**, 174301 (2015); B. Gurrutxaga-Lerma, D. Balint, D. Dini, and A. Sutton, The mechanisms governing the activation of dislocation sources in aluminum at different strain rates, *J. Mech. Phys. Solids* **84**, 272 (2015); B. Gurrutxaga-Lerma, D. S. Balint, D. Dini, D. E. Eakins, and A. P. Sutton, The role of homogeneous nucleation in planar dynamic discrete dislocation plasticity, *ASME. J. Appl. Mech.* **82**, 071008 (2015).

[18] M. A. Shehadeh, E. M. Bringa, H. M. Zbib, J. M. McNaney, and B. A. Remington, Simulation of shock-induced plasticity including homogeneous and heterogeneous dislocation nucleations, *Appl. Phys. Lett.* **89**, 171918 (2006); H. Bei, Y. F. Gao, S. Shim, E. P. George, and G. M. Pharr, Strength differences arising from homogeneous versus heterogeneous dislocation nucleation, *Phys. Rev. B* **77**, 060103 (2008).

[19] J. Schiøtz, K. W. Jacobsen, and O. H. Nielsen, Kinematic generation of dislocations, *Phil. Mag. Lett.* **72**, 245 (1995); J. Schiøtz, T. Leffers, and B. Singh, Modeling of dislocation generation and interaction during high-speed deformation of metals, *Radiat. Eff. Defects Solids* **157**, 193 (2002).

[20] M. A. Shehadeh and H. M. Zbib, On the homogeneous nucleation and propagation of dislocations under shock compression, *Philos. Mag.* **96**, 2752 (2016).

[21] P. S. Lomdahl and D. J. Srolovitz, Dislocation generation in the two-dimensional Frenkel-Kontorova model at high stresses, *Phys. Rev. Lett.* **57**, 2702 (1986).

[22] R. Baggio, O. U. Salman, and L. Truskinovsky, Inelastic rotations and pseudoturbulent plastic avalanches in crystals, *Phys. Rev. E* **107**, 025004 (2023); R. Baggio, O. Salman, and L. Truskinovsky, Homogeneous nucleation of dislocations as a pattern formation phenomenon, *Eur. J. Mech. - A/Solids* **99**, 104897 (2023).

[23] E. Orowan, Problems of plastic gliding, *Proc. Phys. Soc.* **52**, 8 (1940).

[24] J. W. Taylor, Dislocation dynamics and dynamic yielding, *J. Appl. Phys.* **36**, 3146 (1965).

[25] L. Davison, *Fundamentals of Shock Wave Propagation In Solids*, High-Pressure Shock Compression of Condensed Matter (Springer, Berlin, 2008).

[26] M. A. Meyers, A mechanism for dislocation generation in shock-wave deformation, *Scr. Metall.* **12**, 21 (1978); J. Weertman, Moving dislocations in a shock front, in *Shock Waves and High-Strain-Rate Phenomena in Metals*, Elsevier Reference Monographs, edited by M. A. Meyers and L. E. Murr (Springer US, 1981) Chap. 29, pp. 469–486.

[27] C. S. Smith, Metallographic studies of metals after explosive shock, *Trans. Met. Soc. AIME* **212**, 574 (1958); E. Hornbogen, Shock-induced dislocations, *Acta Metall.* **10**, 978 (1962).

[28] J. Weertman, Uniformly moving transonic and supersonic dislocations, *J. Appl. Phys.* **38**, 5293 (1967); J. Weertman and P. Follansbee, Dislocation dynamics and plastic shock waves, *Mech. Mater.* **7**, 177 (1988).

[29] M.-A. Meyers, *Dynamic behavior of materials* (Wiley, 1994).

[30] P. Gumbsch and H. Gao, Dislocations faster than the speed of sound, *Science* **283**, 965 (1999); Z. Jin, H. Gao, and P. Gumbsch, Energy radiation and limiting speeds of fast moving edge dislocations in tungsten, *Phys. Rev. B* **77**, 094303 (2008).

[31] H. Tsuzuki, P. S. Branicio, and J. P. Rino, Accelerating dislocations to transonic and supersonic speeds in anisotropic metals, *Appl. Phys. Lett.* **92**, 191909 (2008).

[32] E. N. Hahn, S. Zhao, E. Bringa, and M. Meyers, Supersonic dislocation bursts in silicon, *Sci. Rep.* **6**, 26977 (2016).

[33] S. Ryu, K. Kang, and W. Cai, Entropic effect on the rate of dislocation nucleation, *Proc. Natl. Acad. Sci. U.S.A.* **108**, 5174 (2011).

[34] B. Gurrutxaga-Lerma, M. A. Shehadeh, D. S. Balint, D. Dini, L. Chen, and D. E. Eakins, The effect of temperature on the elastic precursor decay in shock loaded fcc aluminium and bcc iron, *Int. J. Plast.* **96**, 135 (2017).

[35] Y.-P. Pellegrini, Dynamic Peierls-Nabarro equations for elastically isotropic crystals, *Phys. Rev. B* **81**, 024101 (2010).

[36] Y.-P. Pellegrini, Equation of motion and subsonic-transonic transitions of rectilinear edge dislocations: A collective-variable approach, *Phys. Rev. B* **90**, 054120 (2014).

[37] M. Josien, *Etude mathématique et numérique de quelques modèles multi-échelles issus de la mécanique des matériaux*, Ph.D. thesis, École des Ponts (2018).

[38] A. Maurel, J.-F. Mercier, and F. Lund, Scattering of an elastic wave by a single dislocation, *J. Acoust. Soc. Am.* **115**, 2773 (2004).

[39] B. Gurrutxaga-Lerma and J. Verschueren, Generalized Kanzaki force field of extended defects in crystals with applications to the modeling of edge dislocations, *Phys. Rev. Mater.* **3**, 113801 (2019).

[40] Supplemental material at *URL will be inserted by publisher*. This resource contains Refs. [87–89] and Ref. Hirth et al. in [80].

[41] R. Peierls, The size of a dislocation, *Proc. Phys. Soc.* **52**, 34 (1940); F. R. N. Nabarro, Dislocations in a simple cubic lattice, *Proc. Phys. Soc.* **59**, 256 (1947).

[42] J. D. Eshelby, The equation of motion of a dislocation, *Phys. Rev.* **90**, 248 (1953).

[43] X. Markenscoff and S. Huang, Analysis for a screw dislocation accelerating through the shear-wave speed barrier, *J. Mech. Phys. Solids* **56**, 2225 (2008).

[44] M. Lazar and Y.-P. Pellegrini, Distributional and regularized radiation fields of non-uniformly moving straight dislocations, and elastodynamic Tamm problem, *J. Mech. Phys. Solids* **96**,
632 (2016).

[45] Y. Cui, G. Po, Y.-P. Pellegrini, M. Lazar, and N. Ghoniem, Computational 3-dimensional dislocation elastodynamics, *J. Mech. Phys. Solids.* **126**, 20 (2019).

[46] J. Weertman, Dislocations in uniform motion on slip or climb planes having periodic force laws, in *Mathematical Theory of Dislocations*, edited by T. Mura (ASME, New York, 1969) pp. 178-202; Stress dependence on the velocity of a dislocation moving on a viscously damped slip plane, in *Physics of Strength and Plasticity*, edited by A. Argon (The M.I.T. Press, Cambridge, MA., 1969) Chap. 7, pp. 75-83.

[47] P. Rosakis, Supersonic dislocation kinetics from an augmented Peierls model, *Phys. Rev. Lett.* **86**, 95 (2001).

[48] M. Josien, Y.-P. Pellegrini, F. Legoll, and C. Le Bris, Fourier-based numerical approximation of the Weertman equation for moving dislocations, *Int. J. Numer. Meth. Engrg.* **113**, 1827 (2018).

[49] S. Ogata, J. Li, N. Hirosaki, Y. Shibutani, and S. Yip, Ideal shear strain of metals and ceramics, *Phys. Rev. B* **70**, 104104 (2004).

[50] J. Marian and A. Caro, Moving dislocations in disordered alloys: Connecting continuum and discrete models with atomistic simulations, *Phys. Rev. B* **74**, 024113 (2006).

[51] D. L. Olmsted, L. G. Hector Jr, W. A. Curtin, and R. J. Clifton, Atomistic simulations of dislocation mobility in Al, Ni and Al/Mg alloys, *Model. Simul. Mater. Sci. Eng.* **13**, 371 (2005).

[52] N. Daphalapurkar, J. Wilkerson, T. Wright, and K. Ramesh, Kinetics of a fast moving twin boundary in nickel, *Acta Mater.* **68**, 82 (2014).

[53] E. Oren, E. Yahel, and G. Makov, Dislocation kinematics: a molecular dynamics study in Cu, *Model. Simul. Mat. Sci. Eng.* **25**, 025002 (2017).

[54] J. W. Christian and V. Vítek, Dislocations and stacking faults, *Rep. Prog. Phys.* **33**, 307 (1970).

[55] J. R. Rice, Dislocation nucleation from a crack tip: An analysis based on the peierls concept, *J. Mech. Phys. Solids* **40**, 239 (1992); G. Xu, A. S. Argon, and M. Ortiz, Nucleation of dislocations from crack tips under mixed modes of loading: Implications for brittle against ductile behaviour of crystals, *Philos. Mag. A* **72**, 415 (1995).

[56] In Refs. 36 and 48, this coefficient was defined as $\kappa_a^\alpha = \kappa_a + \alpha$. The present modification only (formally) impacts the *climb* kernel. It rationalizes the pole structure of the kernels in the plane $(\alpha, \gamma)$, where the value $\alpha = 1$ is special; see Supplemental Material [40].

[57] J. D. Eshelby, Dislocations in visco-elastic materials, *Phil. Mag.* **6**, 953 (1961); A. S. Nowick and B. S. Berry, *Anelastic Relaxation in Crystalline Solids*, edited by A. Alper, J. L. Margrave, and A. S. Nowick, Materials Science Series (Academic Press, 1972); M. A. Meyers, A model for elastic precursor waves in the shock loading of polycrystalline metals, *Mater. Sci. Eng.* **30**, 99 (1977).

[58] A particular instance of the latter definition is implicit in the rectangular loop model [90].

[59] D. Lorentz, A. Zeckzer, U. Hilpert, P. Grau, H. Johansen, and H. S. Leipner, Pop-in effect as homogeneous nucleation of dislocations during nanoindentation, *Phys. Rev. B* **67**, 172101 (2003).

[60] J. K. Mason, A. C. Lund, and C. A. Schuh, Determining the activation energy and volume for the onset of plasticity during nanoindentation, *Phys. Rev. B* **73**, 054102 (2006).

[61] S. Mridha, M. Sadeghilaridjani, and S. Mukherjee, Activation volume and energy for dislocation nucleation in multi-principal element alloys, *Metals* **9**, 263 (2019).

[62] D. Wallace, Nature of the process of overdriven shocks in metals, *Phys. Rev. B* **24**, 5607 (1981).

[63] D. Preston, D. Tonks, and D. Wallace, Model of plastic deformation for extreme loading conditions, *J. Appl. Phys.* **93**, 211 (2003).

[64] J. H. Weiner and M. Pear, Breakdown in high-speed dislocation motion, *Phil. Mag.* **31**, 679 (1975).

[65] S. Aubry, K. Kang, S. Ryu, and W. Cai, Energy barrier for homogeneous dislocation nucleation: Comparing atomistic and continuum models, *Scr. Mater.* **64**, 1043 (2011).

[66] Y. M. Gupta, G. E. Duvall, and G. R. Fowles, Dislocation mechanisms for stress relaxation in shocked LiF, *J. Appl. Phys.* **46**, 532 (1975).

[67] J. Li, The mechanics and physics of defect nucleation, *MRS Bulletin* **32**, 151 (2007).

[68] L. Pillon and C. Denoual, Inertial and retardation effects for dislocation interactions, *Philos. Mag.* **89**, 127 (2009).

[69] S. Peng, Y. Wei, Z. Jin, and W. Yang, Supersonic screw dislocations gliding at the shear wave speed, *Phys. Rev. Lett.* **122**, 045501 (2019).

[70] E. M. Bringa, K. Rosolankova, R. E. Rudd, B. A. Remington, J. S. Wark, M. Duchaineau, D. H. Kalantar, J. Hawreliak, and J. Belak, Shock deformation of face-centred-cubic metals on subnanosecond timescales, *Nat. Mater.* **5**, 805 (2006).

[71] R. Madec, B. Devincre, and L. P. Kubin, From dislocation junctions to forest hardening, *Phys. Rev. Lett.* **89**, 255508 (2002).

[72] A. Arsenlis and D. Parks, Crystallographic aspects of geometrically-necessary and statistically-stored dislocation density, *Acta Mater.* **47**, 1597 (1999).

[73] D. Gómez-García, B. Devincre, and L. P. Kubin, Dislocation patterns and the similitude principle: 2.5d mesoscale simulations, *Phys. Rev. Lett.* **96**, 125503 (2006).

[74] Y.-P. Pellegrini, Screw and edge dislocations with time-dependent core width: From dynamical core equations to an equation of motion, *J. Mech. Phys. Solids* **60**, 227 (2012).

[75] J. Weertman, Stress and displacement fields of an edge dislocation that climbs with a uniform velocity, *J. Appl. Phys.* **38**, 2612 (1967).

[76] A. Cochard and J. R. Rice, A spectral method for numerical elastodynamic fracture analysis without spatial replication of the rupture event, *J. Mech. Phys. Solids* **45**, 1393 (1997); P. H. Geubelle and J. R. Rice, A spectral method for three-dimensional elastodynamic fracture problems, *J. Mech. Phys. Solids* **43**, 1791 (1995).

[77] Indeed, $\int_0^x dy J_1(y)/y = \frac{\pi x}{2} [J_1(x)\mathbf{H}_0(x) - J_0(x)\mathbf{H}_1(x)] + xJ_0(x) - J_1(x)$, where $\mathbf{H}_\nu$ is Struve's function. This follows from formula 6.511.6 in I.S. Gradshteyn and I.M. Ryzhik, *Table of Integrals, Series and Products*, 7th. ed. (Academic Press, Amsterdam, 2007) and the identity $J_0(x) = J_1(x)/x + J'_1(x)$.

[78] I. M. Gel'fand and G. E. Shilov, *Generalized Functions*, Vol. 1 (Academic Press, New York, 1964).

[79] K. Ravi-Chandar, *Dynamic Fracture* (Elsevier, 2004).

[80] R. J. Beltz, T. L. Davis, and K. Malén, Some unifying relations for moving dislocations, *phys. stat. sol. (b)* **26**, 621 (1968); J. P. Hirth, H. M. Zbib, and J. Lothe, Forces on high velocity dislocations, *Model. Simul. Mat. Sci. Eng.* **6**, 165 (1998); Y.-P. Pellegrini, Dynamic Peach-Koehler self-force, inertia, and radiation damping of a regularized dislocation, preprint 10.48550/arXiv.2005.12704 (2020), (unpub.).

[81] Y. Xiang, H. Wei, P. Ming, and W. E, A generalized Peierls-Nabarro model for curved dislocations and core structures of dislocation loops in Al and Cu, *Acta Mater.* **56**, 1447 (2008).

[82] W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, *Numerical recipes, 3rd ed.* (Cambridge University Press, Cambridge,
2007)

[83] C. Lubich and A. Schädle, Fast convolution for nonreflecting boundary conditions, *SIAM J. Sci. Comput.* **24**, 161 (2002); A. Schädle, M. López-Fernández, and C. Lubich, Fast and oblivious convolution quadrature, *SIAM J. Sci. Comput.* **28**, 421 (2006).

[84] E. Hairer and G. Wanner, *Solving ordinary differential equations. II. Stiff and differential-algebraic problems*, 2nd ed. (Springer, Berlin, 1996).

[85] J. D. Eshelby, Supersonic dislocations and dislocations in dispersive media, *Proc. Phys. Soc.* **69**, 1013 (1956).

[86] A. B. Movchan, R. Bullough, and J. R. Willis, Stability of a dislocation: discrete model, *Eur. J. Appl. Math.* **9**, 373 (1998).

[87] L. Pillon, C. Denoual, and Y.-P. Pellegrini, Equation of motion for dislocations with inertial effects, *Phys. Rev. B* **76**, 224105 (2007).

[88] B. Gurrutxaga-Lerma, J. Verschueren, A. P. Sutton, and D. Dini, The mechanics and physics of high-speed dislocations: a critical review, *Int. Mater. Rev.* **66**, 215 (2020).

[89] D. G. Duffy, *Transform Methods for Solving Partial Differential Equations*, 2nd ed. (Chapman and Hall/CRC, New York, 2004).

[90] W. G. Johnston and J. J. Gilman, Dislocation velocities, dislocation densities, and plastic flow in lithium fluoride crystals, *J. Appl. Phys.* **30**, 129 (1959).
# Supplemental Material: Shock-driven nucleation and self-organization of dislocations in the dynamical Peierls model

Y.-P. Pellegrini¹,² and M. Josien³

¹ CEA, DAM, DIF, F-91297 Arpajon, France

² Université Paris-Saclay, LMCE, 91680 Bruyères-le-Châtel, France

³ CEA Cadarache, F-13108 St-Paul-lez-Durance, France

## S I. MEAN-FIELD COLLECTIVE-VARIABLE APPROXIMATION VS. FULL NUMERICAL SOLUTION

This section presents, for one edge dislocation, comparisons aimed at cross-checking the numerical solution of the Dynamical Peierls Equation (DPE) against its mean-field collective-variable-approximation (CVA) of Ref. [36] in the main text (hereafter referred to as (I)).

<figure><img src="image_16.png" /></figure>

Figure S1. Edge dislocation speed $v(t)$ vs. time $t$, for: (a) single-step loading, with: (blue) loadings $\sigma_a < \sigma_c$; (orange) loadings $\sigma_a > \sigma_c$, where $\sigma_c \approx 0.401\sigma_{th}$; (b) double-step loading, black: CVA, orange: present full solution. Compare with Fig. 3 and Fig. 11c of (I).

Figure S1(a) displays the dislocation speed vs. time for various single-step loadings $0.048 \le \sigma_a/\sigma_{th} \le 0.952$ from the full numerical solution. Comparing with Fig. 3 of (I), the overall behavior is qualitatively similar, with same asymptotic states (a.s.) (see Fig. 1 in the main text).

Time-dependent inertia of radiation-reaction origin is less [s1,s2] than inertia approximated with velocity-dependent masses [s3]. Figure S1(a) further shows that the full numerical solution experiences more inertia and damping than the CVA, as illustrated by Figs. S1(a,b): unlike the CVA in (I) the full solution is fully overdamped, with extra delay in the response to loading, and in the transitions. All else being equal, the CVA underestimates radiative damping since it relies on a fixed master core shape, whereas the large number of degrees of freedom of the full solution offers a larger number of dissipation routes.

In Fig. S1(b), analogous to Fig. 11(c) of (I) for double-step loading, both solutions are superimposed, with reasonable match. The finding of (I) that the asymptotic state strongly depends on the applied loading is thus confirmed by the DPE.

As in (I), and for a phonon-drag coefficient $\alpha = 0.01$, a critical loading $\sigma_c \approx 0.401 \sigma_{th}$ separates subsonic from intersonic asymptotic states. Although slightly lower, this value is close to the CVA value $\sigma_c \approx 0.415 \sigma_{th}$ reported in (I). The same holds for larger values of $\alpha$; see Fig. S2. Therefore, in spite of slight differences, overall good agreement between the dynamics in the

<figure><img src="image_17.png" /></figure>

Figure S2. Critical stress $\sigma_c(\alpha)$ vs. phonon drag coefficient $\alpha$. Solid: CVA theory; orange markers: full numerical solution. Compare with Fig. 5 of (I).

full numerical solution of the DPE and in its CVA is found. Specifically, the CVA predictions regarding the subsonic-intersonic transition of edge dislocations in the DPE are confirmed.
## S II. DENSITY DATA FOR SCREW DISLOCATIONS

Density data for screw dislocations similar to Figs. 13 and 14 in the main text are as follows. Figure S3(a) displays mean

<figure><img src="image_18.png" /><figcaption>Figure S3. Edge dislocations. (a) Mean density $\bar{\rho}$ in the front zone vs. V for various stresses and (b) data collapse with rescaling exponent $\beta = 1/2$.</figcaption></figure>

densities in the front zone. The rescaling in (b) is consistent with an exponent $\beta = 1/2$ over the whole range of V. Mean density

<figure><img src="image_19.png" /><figcaption>Figure S4. Edge dislocations. (a) Mean density $\bar{\rho}$ in the bulk zone vs. V for various stresses and (b) data collapse with different rescaling exponents $\beta = 0.66$ in region A, and $\beta = 0.45$ in region B.</figcaption></figure>

data in the bulk zone are reported in Figure S4(a), while in (b) different exponents were needed for $S(\tau)$ in regions A ($V < c_s$) and B ($V > c_s$) as indicated in the caption. As in the edge case, ‘bulk’ data thus needs two different exponents for $V < c_s$ and $V > c_s$.

## S III. SINGULARITIES OF LAPLACE TRANSFORM OF THE RESOLVENT

This section addresses the issue of the singularities – poles of cuts – in the complex plane of the Laplace transform (B4) of the resolvent, namely,

$$ \mathcal{R}_a(s) = \frac{1}{\kappa_a^\alpha s + C_a(s)}. \qquad (S1) $$

The necessary assumption for the numerical Laplace inversion technique described in Appendix B of the main text is that these singularities all lie in the half-plane $\Re(s) \le 0$.

To this aim a particular determination must be attributed to square roots in expressions of $C_a(s)$. Consider for instance the screw-dislocation case where, from Eqs. (2), (A4c), and (B4) of the main text

$$ \mathcal{R}_s(s) = \frac{1}{\alpha s + \sqrt{1 + s^2}}. \qquad (S2) $$

The function $\sqrt{1 + s^2} = \sqrt{(s + i)(s - i)}$ has branch points $s = \pm i$. We fix its determination by taking the branch cuts of the square root on the semi-infinite lines $s = \pm i + \lambda$, with $\lambda < 0$. The determination of the argument is taken such that $s - i = \rho_1 e^{i\theta_1}$ and $s + i = \rho_2 e^{i\theta_2}$, with the angles $\theta_1, \theta_2$ in $]-\pi, \pi[$. Since $\sqrt{(s + i)(s - i)} = \sqrt{1 + s^2} = \sqrt{\rho_1\rho_2e^{i(\theta_1+\theta_2)/2}}$, this choice of
branch cuts amounts to computing $\sqrt{1+s^2}$ as $\sqrt{s+i\sqrt{s-i}}$, using the principal determination of the square root in the last expression [s4]. For the 'edge' cases [Eqs. (A4f) and (A4i) of the main text] additional branch points are $s = \pm i\gamma$, and one would similarly take $\sqrt{1+s^2/\gamma^2} \rightarrow \sqrt{s/\gamma + i\sqrt{s/\gamma - i}}$ (recall that $\gamma = c_L/c_S$). Thus, each wavespeed generates two cuts of opposite imaginary part parallel to the negative real axis. Plots of $|R_a(z)|$ using expressions of Appendix A and square roots interpreted thus are displayed in Fig. S5.

<figure><img src="image_20.png" /><figcaption>Figure S5. Singularities of the function $R_a(s)$, denoted here as $f(z)$, for $\gamma = 2$, and $\alpha = 0.5$ (for better display of the cuts). From left to right: screw, glide edge, climb edge. The four poles of the 'edge' cases lie close to the branch cuts.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure S6. Number of poles of $R_a(z)$ in parameter space $(\alpha, \gamma)$.</figcaption></figure>

Turning now to the poles $R_a(s)$, their numbers in $(\alpha, \gamma)$ parameter space are reported in Fig. S6 for completeness (recall however that the physical order of magnitude of $\alpha$ is $\alpha \simeq 0.01$, and that $\gamma$ is of order 2). For all three dislocation characters, the value $\alpha = 1$ is a special one, if $\alpha$ is introduced as in Eq. (2) of the main text. In region $\alpha < 1$, all poles occur in conjugate pairs of nonzero imaginary part, wheras in region $\alpha > 1$ one of the poles is real negative, which indicates overdamped relaxation. All poles are simple ones, of non-positive real part.

Kernel $R_s(s)$ has only one pole, namely $s = -1/\sqrt{\alpha^2 - 1}$ if $\alpha > 1$, and no pole otherwise. In the 'edge' cases the algebraic analytical expressions of $R_{c,g}(s)$ allow for a determination of the poles by solving numerically a 12th-degree polynomial equation, and by selecting the relevant solutions by substitution into the denominator of $R_{c,g}(s)$. All poles are located in-between the cuts $s \in (-\infty, -\pm i\gamma, \pm i\gamma]$.

References:

[s1] L. Pillon, C. Denoual, and Y.-P. Pellegrini, Equation of motion for dislocations with inertial effects, Phys. Rev. B **76**, 224105 (2007).

[s2] B. Gurrutxaga-Lerma, J. Verschueren, A.P. Sutton, and D. Dini, The mechanics and physics of high-speed dislocations: a critical review, Int. Mater. Rev. **66**, 215-255 (2020).

[s3] J. P. Hirth, H. M. Zbib, and J. Lothe, Forces on high-velocity dislocations, Model. Simul. Mat. Sci. Eng. **6**, 165 (1998).

[s4] D. G. Duffy, *Transform Methods for Solving Partial Differential Equations, 2nd ed.*, (Chapman and Hall/CRC, New York, 2004).