Anomalous phonon dispersion near yielding in athermal crystals

Abstract

Vibrational properties of ordered athermal solids near yielding remain poorly understood. We show that yielding in a sheared crystal is governed not by a single localized instability but by directionally extended multimode softening that forms a cross-shaped low-frequency region in wave number space. Near yielding, the acoustic dispersion $ω\sim k$ is replaced by $ω\sim k^2$ along the soft direction, and the vibrational density of states crosses over from Debye to non-Debye scaling, with a diverging length scale. We analytically derive these scaling laws.

Figures (5)

• Figure 1: Schematic of a two-dimensional crystalline system composed of $N_x \times N_y$ Hertzian particles of diameter $d$ arranged on a triangular lattice. The initial center-to-center distance between neighboring particles is $(1-h)d$, where $h$ denotes the overlap parameter. Periodic boundary conditions are imposed in both the $x$ and $y$ directions. A quasistatic shear strain $\gamma$ is applied, and yielding occurs at a characteristic strain $\gamma_c$, defined as the point where the smallest eigenvalue of the Hessian first vanishes. In the theoretical analysis, a representative particle labeled $0$ and its six nearest neighbors labeled $1$–$6$ are explicitly considered.
• Figure 2: (a) Minimum Hessian eigenvalue $\min_{\bm{k}\ne \bm{0}}\omega_{-}^{2}(\bm{k})$ [Eq. \ref{['eqmain:eigenvalue']}], evaluated numerically, as a function of shear strain $\gamma$ for a system with $N_x=N_y=6$ and $h=0.05$. The minimum eigenvalue decreases monotonically and defines a yield strain, $\gamma_c$, where it vanishes (red dotted line, also shown in (b)). Inset: Logarithmic plot of $\omega_{-}^{2}$ versus the distance to yielding, $\Delta\gamma=\gamma_c-\gamma$, demonstrating the linear scaling $\omega_{-}^{2}\sim\Delta\gamma$ near yielding. (b) Shear stress $\sigma_{xy}$ [Eq. \ref{['eqmain:stress']}], evaluated numerically (red curve), as a function of $\gamma$ for the same system, in quantitative agreement with DEM simulation results (gray symbols). The stress drop in the DEM data occurs at the numerically predicted $\gamma_c$ (dotted line), validating the numerical framework. Inset: Numerical eigenmode [Eq. \ref{['eqmain:eigenvector']}] associated with the minimum eigenvalue at yielding, showing a plane-wave character (red arrows).
• Figure 3: (a) Two-dimensional dispersion relations $\omega_{-}(k_x,k_y)$ at $\gamma=0$ and near yielding ($\Delta\gamma=10^{-4},10^{-8}$). Color indicates $\omega_{-}$. At $\gamma=0$, low-frequency modes are concentrated near $k=0$. Near yielding, a cross-shaped low-frequency region appears, indicating directional soft modes. (b) Angular dependence of $\omega_{-}/k$ at several $k$ for various $\Delta\gamma$. Near yielding ($\Delta\gamma\to 0$), a minimum develops around $\theta/\pi\simeq0.5$. The critical angle $\theta_c$ is defined as the angle at which the minimum occurs at yielding. (c) Dispersion relation along the direction $\theta_c$, $\omega_{-}$ vs $k$, for various $\Delta\gamma$. Far from yielding, $\omega_{-}\sim k$; near yielding, the dispersion crosses over to $\omega_{-}\sim k^2$. (d) Vibrational density of states for various $\Delta\gamma$, computed from all $\omega_{\pm}(\bm{k})$. Far from yielding, $D(\omega)\sim\omega$ (Debye scaling); near yielding, $D(\omega)\sim\omega^{1/2}$. Panels (a)--(c) are obtained numerically from Eq. \ref{['eqmain:eigenvalue']} in the continuum limit ($N_{x}, N_{y}\to\infty$), while (d) uses a large system $N_{x}=N_{y}=4000$. The scalings $\omega_{-}\sim k^2$ and $D(\omega)\sim\omega^{1/2}$ are analytically derived, including prefactors, in the Appendix (Eqs. \ref{['eqmain:scaled_omega']} and \ref{['eqmain:analytical_vdos']}).
• Figure 4: Scaled dispersion relations near yielding for various $\Delta\gamma$ and $h$, plotted as $\omega_{-}\Delta\gamma^{-1}h^{11/4}$ versus $k\Delta\gamma^{-1/2}h^{2}$, showing data collapse. Symbols denote numerical results obtained from Eq. \ref{['eqmain:eigenvalue']}. The data exhibit linear ($\omega_{-}\sim k$) and quadratic ($\omega_{-}\sim k^2$) regimes with a crossover at $k\Delta\gamma^{-1/2}h^{2}\sim 1$. The black curve shows the analytical result [Eq. \ref{['eqmain:scaled_omega']}], in quantitative agreement with the numerical data. Inset: Unscaled dispersion relations for different $\Delta\gamma$ and $h$ before rescaling.
• Figure 5: Scaled dispersion relations near yielding for a two-dimensional crystal with WCA interactions for various $\Delta\gamma$, with initial nearest-neighbor distance $2^{1/6} - 0.1$, plotted as $\omega_{-}\Delta\gamma^{-1}$ versus $k\Delta\gamma^{-1/2}$, showing data collapse with linear ($\omega_{-}\sim k$) and quadratic ($\omega_{-}\sim k^2$) regimes and a crossover at $k\Delta\gamma^{-1/2}\sim 1$.