Stress Wave Propagation in Two-dimensional Buckyball Lattice

Orderly arrayed granular crystals exhibit extraordinary capability to tune stress wave propagation. Granular system of higher dimension renders many more stress wave patterns, showing its great potential for physical and engineering applications. At nanoscale, one-dimensionally arranged buckyball (C60) system has shown the ability to support solitary wave. In this paper, stress wave behaviors of two-dimensional buckyball (C60) lattice are investigated based on square close packing and hexagonal close packing. We show that the square close packed system supports highly directional Nesterenko solitary waves along initially excited chains and hexagonal close packed system tends to distribute the impulse and dissipates impact exponentially. Results of numerical calculations based on a two-dimensional nonlinear spring model are in a good agreement with the results of molecular dynamics simulations. This work enhances the understanding of wave properties and allows manipulations of nanoscale lattice and novel design of shock mitigation and nanoscale energy harvesting devices.

well as an analogue to Nesterenko solitary wave. Substituting buckyballs with carbon nanotubes (resembling cylindrical or tubal particles), we have demonstrated that 1D single-walled carbon nanotube (SWNT) system serves as a highly effective and reusable energy absorber 62 . Further, we studied a C 60 -SWNT hybrid system, having achieved a quantitative tuning of nanoscale solitary waves with good precision 63 .
In this study, we investigate stress wave behaviors in 2D buckyball system via molecular dynamics (MD) simulation, as an extension of our previous work on 1D nanoscale lattice. As coherent macroscale studies suggest, 2D constructions render more types of wave patterns, suggesting a very promising potential for stress wave-related applications. The presentation of our work is structured as follows. In Sec. 2, we first introduce the setup of investigated system, including the square close packing (scp) and the hexagonal close packing (hcp) system. In Sec. 3, we put forward a theoretical model to describe 2D problem. In Sec. 4, we present the results of MD simulations and numerical calculations based on the 2D NS model of the two packing modes. Finally, in Sec. 5, we present some concluding remarks and possibilities for future investigations.

System descriptions
In this study, we investigate two basic 2D uniform packing modes, i.e. scp and hcp, whose coordination numbers 64 (or contact numbers) are 4 and 6 respectively and hcp is the densest packing mode possible (with a packing density of ~0.9609), as is illustrated in Fig. 1. In this paper, the size of investigated systems is chosen as 10-order. The shortest distance between two C 60 molecules equals to the equilibrium spacing d 0 = 10.05Ǻ 65 , corresponding to the close packing modes. A C 60 molecule on one edge away from the corners is selected as an impactor to generate a stress wave with an initial velocity of 2000 m/s, a moderate impacting velocity for studying nanoscale impact dynamics. The impacting direction is determined by an impacting angle θ, defined as the angle between the velocity and the y axis (see Fig. 1). The MD investigation of stress wave behaviors of the system is carried out based on the open-source program LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) 66 and visualized using molecular visualization program VMD (Visual Molecular Dynamics) 67 . To minimize the disturbance of thermal vibration, the initial temperature is set as T = 10 K and the systems are run for equilibrium in the canonical ensemble (NVT) for 3 ps, which is a sufficient time duration in this study (see Figure S1). Major investigations of wave behaviors are in the micro-canonical ensemble (NVE). More detailed simulation descriptions are provided in the Supplementary Information.

Theoretical modeling
The dynamics of the 2D system is modeled by a 2D NS model, as a dimensional extension of 1D NS model. In this study, the values of both stiffness parameter k and nonlinearity index n are identical to the 1D NS model, whose theoretical background is briefly reviewed in the Supplementary Information. The Hamiltonian of the whole C 60 system can be given as follows where V(r) is the potential function of two interacting C 60 molecule with respect to distance r; u i is the displacement vector from the initial position of i th C 60 molecule; e ij is the unit vector pointing from i th C 60 molecule to j th C 60 molecule.
In numerical calculations of scp and hcp arrangements all interactions of a C 60 molecule with 8 surrounding molecules are taken into consideration (see the regions surrounded by dashed lines in Fig. 1) without any symmetrical or collision mode-related simplifications. It is advantageous over previous numerical methods with respect to calculation precision [46][47][48][49]55 . The 2D equation of motion of i th C 60 can be given as The numerical calculations are performed using a forth order Runge-Kutta algorithm to integrate the nonlinear ordinary differential equations 68 . As an initial condition, the displacement vector u i of each C 60 equals to zero vector and the velocity vector u i is also zero vector except for the impactor, whose initial velocity is nonzero. The numerical calculation lasts 3000 fs each time.

Results
We choose two impacting angles to study the impact responses of nanoscale lattice, i.e. θ = 0° and θ = 30°. The magnitude of particle velocity is extracted to characterize stress wave propagation. The investigation of wave behaviors in scp and hcp configurations is organized as follows: (a) particle velocity magnitude v M distributions at given instants; (b) particle velocity amplitude v Amp distributions; (c) comparison between MD simulation results and numerical calculation results, i.e. solutions of Eq. (2); (d) Discussions on the unique wave properties in scp and hcp configurations.
Square close packing. Particle velocity magnitude distributions at given instants and particle velocity amplitude distributions of scp configuration under two different impacting angles are presented in Fig. 2.
As is illustrated in Fig. 2, stress wave in scp system is highly directional where the major portion of wave energy transmits through initially excited chains. In fact, waves traveling in initially excited chains are solitary waves, whose shape and wave speed remain constant as they travel (see Figure S2). Later we will show that these solitary waves are in good agreement with Nesterenko's well-known 1D solitary wave in many aspects. Because of the highly directionality of scp packing modes and the dispersionless characteristics of translational solitary waves, there is potential for scp granular configuration to be applied to stress wave guidance and energy focusing etc.
The comparison between MD simulations and numerical calculations of scp system is shown in Fig. 3, where the time histories of particle velocity magnitude and acceleration on certain direction of 3 representative C 60 s are compared.
According to the comparison in Fig. 3, the results of MD simulations and numerical calculations are in excellent agreement. To confirm that the prediction accuracy of NS model is not sensitive to the position of impactor C 60 molecule, another comparison is made where the impactor is chosen as the C 60 molecule at the bottom left, as is shown in Figure S3.
A further exploration to the properties of excited solitary waves in 2D scp system is as follows. A set of MD simulations are conducted where the impacting angles range from 0 deg to 75 deg, resulting in a series of solitary waves of various force amplitudes (presented in Table S1) propagating along the C 60 chain in the dashed square in Fig. 4. With the theory of 1D NS model (see the Supplementary Information), wave speed can by predicted by inputting wave amplitude and the predicted values match well with simulation results (see Fig. 4). The scaling relation between particle velocity amplitude and wave speed of Nesterenko solitary wave is   Fig. 4, demonstrating that the investigated system satisfies Eq. (3) very well. In this strongly nonlinear system, the width of the wave (characteristic spatial length) L n , according to Nesterenko's theory, can be given as where a is the distance between adjacent particles (in our case a = d 0 ). For buckyball system, Eq. (4) outputs L n ≈ 3a, consistent with the simulation results. It is smaller than 5a in Hertzian chains with very small initial compression due to higher nonlinearity. Above analysis further confirms that this solitary wave obeys Nesterenko's theory.
Hexagonal close packing. Particle velocity magnitude distributions at given instants and particle velocity amplitude distributions of hcp configuration under two different impacting angles are presented in Fig. 5. As can be seen, although stress wave propagation in hcp configuration tends to be directional, the amplitude decays as the wave travels and the input energy tends to be distributed. The reason is that in hcp system adjacent row and column will be affected by moving particles regardless of their velocity directions and thus energy is distributed to two or more surrounding particles. This pattern will continue spreading to more and more particles and energy or particle velocity amplitude is gradually decreasing. Therefore, hcp configuration of C 60 molecules can serve as a nanoscale protective device.
The comparison between MD simulations and numerical calculations of scp system is shown in Fig. 6. Although the agreement between the results of MD simulations and numerical calculations is generally good, it can be observed that there appears a little discrepancy between simulations and model predictions for the velocity and acceleration curves of some selected C 60 molecules and it is not surprising due to the atoms vibrating about the equilibrium positions instead of being absolutely static after stress wave passes through. Therefore, the position or time of a C 60 molecule when it interacts with other surrounding molecules is actually random. In addition, the parameters of NS model are obtained in 1D system where two degrees of freedom are eliminated, which may not be utterly suitable for 2D situation in this study. Due to error accumulation, the precision of model prediction decreases as the length of the force chain increases, and therefore, precisely predicting the velocity or acceleration of C 60 molecule with a complex stress wave path can be fairly difficult. This also explains why the predictions on scp configurations are usually more accurate than hcp configuration on the condition that traveling distances of the stress waves are similar.
To further investigate the wave spreading behavior of nanoscale hcp configuration, a 30° observation line is defined (see the right panel of Fig. 7), along which stress wave propagation is studied. For a θ = 0° impact, the time histories of particle velocity magnitude of C 60 molecules on the observation line is extracted and amplitudes of each molecule are marked, as is shown in the left panel of Fig. 7. Previous research on high dimensional-structured granular system tends to model the spreading of wave energy partition as exponentially decay [55][56][57] . In this study, we fit the particle velocity amplitude in the same manner and the fitting results are satisfactory (see the dashed line in Fig. 7). Moreover, the pulse duration t p on each C 60 molecule is subplotted in Fig. 7, showing that as the amplitude decreases, stress wave stays longer on a single C 60 molecule.

Concluding Remarks
At macroscale, orderly arranged granular materials have demonstrated its power to tune stress waves. At nanoscale, however, stress wave tuning and manipulations based on discrete granular configurations is a rather unexplored topic. In this work, we studied the 2D problem of nanoscale stress wave propagation as an extension of 1D propagation in buckyball system via MD simulation. We show that scp supports highly directional Nesterenko solitary waves along initially excited chain while hcp attenuates impact exponentially by continuously spreading wave energy to adjacent particles. We are able to describe the wave behaviors of both configurations in a good agreement with a 2D NS model. This work further validates the NS model as a feasible simplification of complicated van der Waals interaction in modeling the stress wave behaviors in nanoscale lattice of higher dimension and suggests the possibilities of 2D nanoscale system acting as energy harvesting, guiding or shock mitigating devices. More types of nanoscale granular configurations and their potential applications will be covered in our future work.