Collective dynamics and shattering of disturbed two-dimensional Lennard-Jones crystals

We first consider the analytically tractable cases of 6- and 7-particle hexagonal configurations as shown in Fig. 1, and analyze the normal modes in the perturbation regime, where the interaction between neighboring particles is approximated by the harmonic potential. These two kinds of elementary hexagonal configurations are denoted as $H_{6}$ and $H_{7}$. We analytically derive for all of the normal modes and their frequencies from the eigenvalue equations goldstein2011classical :

where $V_{ij}$ and $T_{ij}$ are the matrices related to the potential and kinetic energies. $\omega_{k}$ and $a_{jk}$ are the frequency and eigenvector of the $k$-mode. The complete results and derivation details are presented in SI. In all of the normal modes, the total momentum and angular moment are conserved, and the center of mass is invariant in time.

In the extension of the system from the elementary hexagonal configurations to a large triangular lattice of hexagonal shape, the analysis of the normal modes in the $H_{6}$ and $H_{7}$ configurations inspires a series of inquiries. First of all, among all of the normal modes, it is found that only the radial modes shown in Fig. 1 possesses the $C_{6}$ symmetry. Will this symmetry be preserved in the large system? In general, will the velocity field as a whole respect the symmetry of the boundary? Furthermore, for the radial mode in the $H_{7}$ case, it is noticed that the velocity of the central particle is zero. Does it imply the appearance of a singularity point in the velocity field in the continuum limit? If yes, will the singularity structure in the velocity field manifest itself as a point or other dynamical morphologies? It is also noticed that the frequency of the radial mode in the $H_{6}$ configuration is smaller than that in the $H_{7}$ configuration due to the vacancy at the center of the $H_{6}$ system; in fact, the radial mode has the lowest frequency among all the normal modes in the $H_{6}$ configuration (see SI). According to Maxwell’s rule regarding the stability of a mechanical structure, removing the central particle in the $H_{7}$ configuration leads to floppy modes in the resulting $H_{6}$ configuration, and thus fundamentally changes the dynamical state of the system maxwell1864calculation ; vitelli2012topological . It is therefore of interest to examine the vacancy-related dynamics in a large lattice.

To address the questions proposed in the preceding subsection, we first resort to high-precision numerical integration of the equations of motion to construct high-quality particle trajectories with well conserved total energy, momentum and angular momentum.

In a disturbed triangular lattice of hexagonal shape, we observe that individual motions of the particles are quickly (in comparison with the characteristic time $\tau_{0}$) organized to form coherent dynamical states. A typical instantaneous velocity field for $b=0.01\%$ (measured in the unit of the lattice spacing $r_{m}$) is presented in Fig. 2(a) for a system of 1951 particles. We see that the entire velocity field exhibits $C_{6}$ symmetry that is consistent with the crystal boundary. Remarkably, the velocity field maintains the $C_{6}$ symmetry in the evolution of the ever-changing pattern. We also investigate the case of rectangular boundaries, and find that the boundary symmetry also dictates the symmetry of the velocity field. Typical instantaneous velocity fields of rectangular crystals are presented in SI.

A salient common feature of the velocity field driven by random perturbation is that the velocity vectors exhibit distinct behaviors in different regions, as shown in Fig. 2(a). A zoomed-in plot of a sixth of the entire $\vec{v}$-field is shown in the second panel. In some region, the motions of the particles follow the average orientation of the neighboring velocity vectors. In the remaining region, the orientations of adjacent velocity vectors seem uncorrelated. To quantify the distinct behaviors of the velocity vectors, we introduce a dynamical order parameter $\alpha_{i}$ defined at particle $i$:

$\langle\vec{v}_{j}\rangle$ is the average velocity of the neighboring particles of particle $i$; the neighbors of a particle are determined by the standard Delaunay triangulation nelson2002defects . $\langle\vec{v}_{j}\rangle=\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}\vec{v}_{j}$, where $n_{i}$ is the coordination number of particle $i$. When $\vec{v}_{i}$ follows the orientation of $\langle\vec{v}_{j}\rangle$, $\alpha_{i}=1$, and the cluster of the particles tends to move coherently. The order parameter $\alpha_{i}$ thus characterizes the local alignment of velocity vectors regardless of their magnitude.

In the third panel in Fig. 2(a), we show the instantaneous $\alpha$-field with $C_{6}$ symmetry. In the brighter region, the value of $\alpha$ is larger. The $\alpha$-field clearly shows that the velocity field is sharply divided into ordered (bright) and disordered (dark) regions. Here, we shall note that the degree of order in the velocity field is specified by the value of the order parameter $\alpha_{i}$ as defined in eqn (3); the system is still in the crystalline state under mild disturbance. Extensive data analysis of both hexagonal and rectangular systems shows that the coexistence of ordered and disordered regions is a common feature of the velocity field in the perturbation regime. Our proposal of the dynamical order parameter $\alpha_{i}$ is inspired by the seminal work of active matter physics, where the rule of specifying particle velocity by the average velocity of surrounding particles has been employed to create coherent collective dynamics vicsek1995novel .

The origin of the disordered region in the velocity field is closely related to the key observation that the boundary particles tend to move perpendicular to the boundary. Such a boundary dynamical state leads to a nonzero winding number in the interior velocity vector field aminov2000geometry . Consequently, singularities are inevitable in the interior velocity field. The disordered region in the velocity field may be regarded as an area full of singularities in the continuum limit. In contrast, by fixing the positions of the boundary particles, our simulations show that the entire velocity field becomes disordered everywhere. Note that the disordered region in the velocity field of the disturbed L-J lattice is distinct from the node structure arising in standing waves.

With the increase of the disturbance amplitude, we observe the dynamical transition of the velocity field from the coexistence of ordered and disordered regions to completely disordered state. In Fig.2(b), we show that when the value of $b$ increases to $5\%$ of the lattice spacing, both the velocity vectors and the $\alpha$-field become disordered. We characterize this dynamical transition in the plot of $\langle\alpha_{i}\rangle$ versus $b$ in Fig.2(c). $\langle\alpha_{i}\rangle$ is the spatiotemporal averaging of the order parameter $\alpha_{i}$ over tens of instantaneous velocity fields during the time interval of $100\tau_{0}$. From Fig.2(c), we see that the velocity field becomes completely disordered ($\langle\alpha_{i}\rangle=0$) when the amplitude of disturbance exceeds about $2\%$ of the lattice spacing. Dynamical transition is also reflected in the kinetic energy curves: the originally periodic energy curve becomes irregular at short time scale when dynamical transition occurs (more information about the energy curves is presented in SI).

Here, we shall emphasize that the crystalline order is still well preserved in the dynamical transition of the velocity field in preceding discussions. The lattice becomes slightly wiggly with the increasing disturbance strength. The resulting variation of the bond-orientational order of the lattice can be quantified by the order parameter $\psi_{6}(\vec{x}_{i})$ Bruinsma1981 :

where $\theta_{j}$ is the angle of the bond connecting the particle $i$ and its neighbor $j$ with respect to some chosen reference line. $n_{i}$ is the coordination number of the particle $i$. From Fig.2(d), we see the slight decline of the spatiotemporally averaged $\psi_{6}$ with the increase of $b$, indicating that the crystalline order is well preserved up to $b\approx 6\%$.

We proceed to explore the regime of stronger disturbance. At $b\approx 7\%$, we observe the transient existence of a few quadrupoles in the crystal. The bond-orientational order is still preserved. A schematic plot of a quadrupole is shown in the right panel in Fig. 3(h). The red and blue dots represent five- and seven-fold disclinations, which are elementary topological defects in two-dimensional triangular lattices. An $n$-fold disclination in triangular lattice is a vertex whose coordination number $n$ is deviated from six. The coordination number of a particle can be uniquely determined by the standard Delaunay triangulation nelson2002defects . One may define the topological charge of an $n$-fold disclination as $q_{i}=(6-n)\pi/3$. Topological charges of the same sign repel and those of opposite signs attract, which resembles the behavior of electric charges nelson2002defects . A pair of five- and seven-fold disclinations form a dislocation. And a pair of dislocations form a quadrupole as shown in Fig. 3(h). While the emergence of an isolated disclination requires a global reconfiguration of the particle array, a quadrupole can be excited by a slight local compression of the particles indicated by the red arrows in Fig. 3(h). Energetically, it required much less energy to excite a quadrupole than either a dislocation or a disclination Chaikin1995 .

According to the elegant KTHNY theory, the disruption of two-dimensional crystals under thermal agitation (melting) involves the consecutive solid-hexatic and hexatic-liquid transitions RevModPhys.60.161 ; nelson2002defects . The intermediate hexatic phase is characterized by the proliferation of dislocations which still preserves the bond-orientational order; the dislocation-unbinding transition leads to the isotropic liquid phase.

Here, for our L-J crystal system of finite size under the stress-free boundary condition, we discover the featured vacancy-driven disruption mode that is distinct from the dislocation-unbinding mechanism in the melting of two-dimensional crystals. Typical particle configurations for $b=8\%$ and $10\%$ are presented in Figs. 3(a)-3(d) and 3(e)-3(g), respectively. Simulations show that the disruption of the crystalline order is initiated from the center of the system. For clarity, only the central circular regions of the hexagonal crystals are shown in Fig. 3. For the case of $b=8\%$, we observe the unbinding of the quadrupole into a pair of dislocations [see Fig. 3(b)], and the appearance of void regions (bubbles) extending several lattice spacings [see Fig. 3(c)]. Line tension emerges once a bubble is formed, and it tends to inhibit the growing of the bubble. At the disturbance amplitude of $b=8\%$, the line-tension effect wins, and we observe the closure of the bubbles, resulting in different kinds of vacancies [see Fig. 3(d); also see SI] yao2017topological .

By increasing the disturbance amplitude to $b=10\%$, we find that the line tension fails to suppress the growing of the bubbles, and the entire crystal is shattered into several pieces of crystallites, as shown in Fig. 3(e)-3(g). It is noticed that the critical value of the disturbance amplitude ($b=10\%$) for the mechanical disruption of the two-dimensional L-J lattice agrees with both the conventional Lindemann melting criterion for three-dimensional solid Lindeman1910 and the modified Lindemann’s criterion for two-dimensional solid zheng1998lindemann ; khrapak2020lindemann . The conventional Lindemann’s criterion states that melting of three-dimensional solid occurs when the square root of the particle mean-squared displacement from the equilibrium position reaches a threshold value of about one tenth of the interparticle distance Lindeman1910 . This criterion is generalized to two-dimensional solid by measuring the displacements of particles in local coordinate zheng1998lindemann ; khrapak2020lindemann . Here, we shall emphasize that in our work the L-J crystal is two-dimensional, and it is treated as a purely classical mechanical system whose evolution conforms to deterministic Hamiltonian dynamics. The Lindemann’s criterion is applicable to the melting of solid as a thermodynamic system zheng1998lindemann ; khrapak2020lindemann .

In contrast, under periodic boundary condition, simulations of L-J systems show the existence of a metastable state between the solid and liquid phases chen1995melting , and the featured defect fraction as predicted by the KTHNY theory wierschem2011simulation . Furthermore, for the case of a hexagonal L-J crystal with fixed boundary, our simulations show that the constraint of fixed area suppresses the development of bubble structures, and the disruption process resembles the melting scenario of two-dimensional crystals according to the KTHNY theory (more information is presented in SI). To conclude, the disruption of the two-dimensional L-J crystal occurs in the form of the vacancy-driven shattering under the stress-free boundary condition, which is in contrast to the dislocation-unbinding mechanism (the scenario of the KTHNY theory) at constant density as implemented by either the periodic boundary condition or the fixed boundary condition.

Fluctuation of the bubble structure in its size inspires us to ask how a void region is healed under line tension. To address this question, we create an $n$-point linear vacancy at the center of a hexagonal L-J crystal simply by removing $n$ particles. Typical instantaneous velocity fields prior to the closure of the $n$-point vacancy are shown in Figs. 4(a) and 4(b). The particles near the vacancy move coherently to shrink the width $d_{\textrm{vac}}$ of the linear vacancy as indicated by the light green bar. Here, the introduction of the linear vacancy breaks the $C_{6}$ symmetry of the system. The symmetry of the resulting velocity field changes accordingly. We track the temporally-varying vacancy width $d_{\textrm{vac}}$ at varying vacancy length $\ell_{\textrm{vac}}$ (without imposing any disturbance), and the results are presented in Fig. 4(c). It turns out that the short- and long-vacancies exhibit distinct dynamical behaviors. The short vacancy ($\ell_{\textrm{vac}}=5a_{0}$) alternately opens and closes. In contrast, the long vacancy ($\ell_{\textrm{vac}}=21a_{0}$) is quickly healed, leaving out a pair of dislocations at its two ends with opposite Burgers vectors. Note that our previous molecular dynamics simulations have shown that infinitely large two-dimensional L-J crystal exhibits similar behaviors yao2014dynamics2 . For the intermediate case of $\ell_{\textrm{vac}}=11a_{0}$, our numerical solution captures the transient open state of the vacancy for a duration of $\Delta t\approx 11$ [see the blue curve in Fig. 4(c)]. By imposing perturbation, the $d_{vac}(t)$ curves at varying vacancy size for $b=1\%$ are almost identical with those in Fig. 4(c).