Metastable states and energy flow pathway in square graphene resonators

Nonlinear interaction between flexural modes is critical to heat conductivity and mechanical vibration of two-dimensional materials such as graphene. Much effort has been devoted to understand the underlying mechanism. In this paper, we examine solely the out-of-plane flexural modes and identify their energy flow pathway during thermalization process. The key is the development of a universal scheme that numerically characterizes the strength of nonlinear interactions between normal modes. In particular, for our square graphene system, the modes are grouped into four classes by their distinct symmetries. The couplings are significantly larger within a class than between classes. As a result, the equations for the normal modes in the same class as the initially excited one can be approximated by driven harmonic oscillators, therefore, they get energy almost instantaneously. Because of the hierarchical organization of the mode coupling, the energy distribution among the modes will arrive at a stable profile, where most of the energy is localized on a few modes, leading to the formation of “natural package” and metastable states. The dynamics for modes in other symmetry classes follows a Mathieu type of equation, thus, interclass energy flow, when the initial excitation energy is small, starts typically when there is a mode that lies in the unstable region in the parameter space of Mathieu equation. Due to strong coupling of the modes inside the class, the whole class will get energy and be lifted up by the unstable mode. This characterizes the energy flow pathway of the system. These results bring fundamental understandings to the Fermi-Pasta-Ulam problem in two-dimensional systems with complex potentials, and reveal clearly the physical picture of dynamical interactions between the flexural modes, which will be crucial to the understanding of their abnormal contribution to heat conduction and nonlinear mechanical vibrations.

Figure 1

(a) Setup of the square graphene resonator. The sheet is approximately 6.9 nm by 7.4 nm. It contains 1992 carbon sites in total, the outmost two layers are fixed, and the 1748 inner sites are able to move. (b) The first four normal modes with different symmetries: symmetric-symmetric (SS), antisymmetric-symmetric (AS), SA, and AA along $x$ and $y$ axes. Light gray (yellow) and black (dark blue) indicate maximum and minimum values of the normal modes, respectively.

Figure 2

The coupling strength $|S_{ik}|$ between the initially excited mode $k$ and the other modes. The initially excited modes are 1, 592, 1165, and 1746 for (a)–(d), belonging to the four different symmetry classes, respectively. The corresponding specific energies [Eq. (2)] are $\epsilon =1.7524\times {10}^{-22}$ J, $2.1343\times {10}^{-22}$ J, $1.2134\times {10}^{-22}$ J, and $3.3978\times {10}^{-22}$ J. The system is 6.9 nm by 7.2 nm (1748 movable sites). The first half modes, say $i\le 874$, are in acoustical branch. The last half modes ($i>874$) are in optical branch. The four symmetry classes are marked by different symbols, as shown in (a).

Figure 3

(a)–(d) The nonlinear mode coupling $|S_{ik}|$ between the initially excited mode $k$ and the other modes. The initially excited mode $k$ is 1, 592, 1165, and 1746 for (a)–(d), which belong to four different symmetry classes, respectively. The symbols are the same as Fig. 2, but only the modes in the same symmetry class are plotted. The corresponding excitation specific energies are smaller than that in Fig. 2, i.e., $\epsilon =7.8919\times {10}^{-27}$, $2.1240\times {10}^{-26},1.8944\times {10}^{-26}$, and $1.5085\times {10}^{-26}$ J. (e)–(h) The plots of $\sqrt{E_i\left(t\right)/f_i\left(t\right)}$ versus $i$ except mode $k$ for the same initial excitations as (a)–(d), respectively, at $t=12$ ps.

Figure 4

Comparison of the energy spectra obtained from full molecular dynamics simulation [Eq. (3), panels (a) and (e)], the full mode coupling dynamics simulation [Eq. (7), panels (b) and (f)], the simplified version of the mode coupling dynamics simulation [Eq. (10), panels (c) and (g)], and the further simplified version of the mode coupling dynamics simulation [Eq. (13), panels (d) and (h)]. The excitation is on the first mode. The initial energy is $2.7242\times {10}^{-23}$ J. The upper panels are for $t=0.027$ ps, and the lower panels are for $t=23.995$ ps.

Figure 5

(a) The stable (white) and unstable (blue or gray) regions of Mathieu equation. Initially, mode 1 is excited with specific energy $\epsilon =1.1148\times {10}^{-24}$ J. The crosses mark the parameters $(Q_i,A_i)$ calculated from Eqs. (24) and (25) for the first four modes. (b) The time evolution of the harmonic energy for the first four modes. (c)–(f) The energy spectrum at different time instances as marked by blue arrows in (b), the corresponding time is 0.02, 4.18, 8.90, and 300 ns, respectively. Inset of (e): zoom-in of the energy spectrum for the lower modes.

Figure 6

The same plots as Fig. 5 except that the initial excitation energy is a little bit larger: $\epsilon =1.4428\times {10}^{-24}$ J. The time instances for the energy spectrum (c)–(f) are 0.03, 14.97, 167.90, and 300 ns, respectively. Inset of (f): zoom-in of the energy spectrum for the lower modes.