FPU physics with nanomechanical graphene resonators: intrinsic relaxation and thermalization from flexural mode coupling
Daniel Midtvedt, Zenan Qi, Alexander Croy, Harold S. Park, Andreas Isacsson
FFPU physics with nanomechanical graphene resonators: intrinsic relaxation andthermalization from flexural mode coupling
Daniel Midtvedt, ∗ Alexander Croy, and Andreas Isacsson
Department of Applied Physics, Chalmers University of Technology, SE-412 96 G¨oteborg, Sweden
Zenan Qi and Harold S. Park
Department of Mechanical Engineering, Boston University, Boston, MA 02215, USA (Dated: September 26, 2018)Thermalization in nonlinear systems is a central concept in statistical mechanics and has beenextensively studied theoretically since the seminal work of Fermi, Pasta and Ulam (FPU). Usingmolecular dynamics and continuum modeling of a ring-down setup, we show that thermalizationdue to nonlinear mode coupling intrinsically limits the quality factor of nanomechanical graphenedrums and turns them into potential test beds for FPU physics. We find the thermalization rate Γto be independent of radius and scaling as Γ ∼ T ∗ /(cid:15) , where T ∗ and (cid:15) pre are effective resonatortemperature and prestrain. Advances in fabrication techniques enable produc-tion and characterization of one and two dimensionalnanoscale mechanical resonators [1–4]. In particular,carbon-based resonators are considered to be promisingfor many applications due to their low mass and highquality factors (Q-factors) [5–8]. It is also known thatthese systems display strongly nonlinear behavior [4, 9],which makes them interesting for investigations of non-linear dynamics. The nonlinearities lead to a couplingbetween the vibrational modes [10, 11]. This couplingallows for intermodal energy-transfer, which facilitatesthe redistribution of energy initially localized in a singlemode.In this respect, the mode-coupling provides a dissipa-tion channel for the fundamental mode (FM) dynamics.In contrast to other dissipation mechanisms previouslystudied in nanomechanical resonators [12–17], this is afundamental intrinsic mechanism and therefore consti-tutes a lower limit on the relaxation rate of the FM. Atfinite temperatures, the effect of the mode couplings willbe two-fold. First, they give rise to fluctuations in theresonator strain leading to dephasing or spectral broad-ening of the resonator [18]. Second, as we show in thisLetter, they allow for energy redistribution among themodes. To distinguish the two effects we consider a ring-down setup, where the total energy of the resonator isconserved and the evolution of the spectral distributionof energy is monitored. This allows us to access the dy-namics of the FM energy.The process of thermalization in a system of non-linearly coupled oscillators was originally considered byFermi, Pasta and Ulam (FPU) in their famous computerexperiment in 1955 [19, 20], and spawned an impressiveamount of research that eventually resulted in the devel-opment of chaos theory [21] and the discovery of solitons[22]. For the FPU problem, it is known that above acertain critical energy density, energy initially fed intothe FM is quickly redistributed among all modes, andthe system approaches a thermal state. This threshold is connected to the onset of widespread chaos in the modedynamics [20, 21] and the stability of localized modes (”q-breathers”) [23, 24]. In recent years, the consensus hasbeen reached that the main features observed in the FPUproblem are not specific to the original model Hamilto-nian [25, 26]. A natural question, which is still underdebate [27], is whether those features can be observed ina physical system. For this to be possible two require-ments need to be fulfilled: first, the nonlinearity has to besufficiently strong to allow for appreciable coupling be-tween resonator modes already at low energies. Second,the time-scale of energy dissipation to the environmentmust be long compared to that for thermalization dueto the mode coupling. We propose that nanomembraneresonators can be used to test the persistence of the FPUphenomena in the thermodynamic limit.While mode-coupling is present in any crystallinemembrane, we here focus on graphene due to the manyrealizations of membrane resonators using this material.To study the intrinsic loss mechanisms in such systems,we utilized classical molecular dynamics (MD) simula-tions to systematically investigate the free vibrations ofpristine circular graphene monolayers with varying ra-dius, pre-strain, temperature and excitation energy. TheMD simulations were performed using LAMMPS [28, 29].The covalent bonds between carbon atoms were modeledby the AIREBO potential [30–33].After an initial relaxation stage, the graphene mono-layer was strained and fixed at its edges. The systemwas then equilibrated at a specific temperature withinthe canonical ensemble (NVT) for 10 ps. Thereafter, themonolayer was actuated by assigning an initial velocityprofile in the out of plane direction corresponding to thefundamental mode shape of the resonator. After thatpoint, the system was allowed to vibrate freely in the mi-cro canonical (NVE, or energy conserving) ensemble for5000 ps with a time step of 1 fs. The entire simulationwas divided into 10 time windows of 500 ps each; thefirst 100 ps of each time window were used for further a r X i v : . [ c ond - m a t . m e s - h a ll ] A p r analysis. Specifically, total kinetic energy and velocityhistory of each atom were recorded to calculate the totalkinetic energy spectrum. We studied four temperatures( T = ∼
0, 50, 100 and 300 K), four radii ( R = 5, 7, 9 and11 nm) and four initial excitation energies ( v max = 2, 5,10 and 15 ˚A/ps). Additionally, the pre-strain (cid:15) pre of thestructures was varied.In Fig. 1, the time evolution of the kinetic energy spec-tral density of a circular graphene sheet of radius 5 nm, (cid:15) pre = 0 .
5% kept at a temperature of 300 K is reportedup to a time of 5 ns. Initially, the FM is excited by avelocity of 10 ˚A/ps.The spectral density of the kinetic energy in Fig. 1 cor-responds to the frequency distribution of the contributionto the kinetic energy from the out-of-plane motion. Thein-plane contribution is negligible in the frequency rangerelevant for flexural vibrations. The area under a peakgives the kinetic energy of that mode. The prominentpeak initially located around ∼
200 GHz corresponds tothe FM, and its energy decreases continuously during thetime evolution. Simultaneously, the frequency of the FMdecreases. As the FM frequency is energy-dependent,this suggests that energy is redistributed among normalmodes in the system.This redistribution is known to occur in the problemof coupled Duffing oscillators. In fact, in a continuummechanics (CM) approach [9, 11] the present system isalso described by a system of coupled Duffing oscillators.In dimensionless form, the equations of motion for themode amplitudes q n are written as (see supplement) ∂ τ q n + ˜ ω n q n + X ijk W ij ; kn q i q j q k = 0 . (1)The coupling matrix W ij ; kn depends only on the geom-etry of the resonator. For drum geometries it is, in con-trast to the FPU case [34], a dense matrix, with permuta-tion symmetry in the indices i ↔ j and k ↔ n [11]. Thedimensionless time and energy are related to the physicalunits through τ = √ (cid:15) pre c L R t, ˜ E = E π(cid:15) c L ρ G R , (2)where c L = p Y /ρ G is the longitudinal speed of soundin graphene, with Y ≈
350 N/m being the two dimen-sional Young’s modulus of graphene [35] and ρ G = 0 . / m the graphene mass density. In deriving (1)we have further scaled the radial coordinate r and thevertical displacement w according to ˜ r = r/R , ˜ w =(1 / w/ ( (cid:15) R ). The linear frequencies are given by thezeros ξ ,n of the zeroth order Bessel function, ˜ ω n = ξ ,n .The two models (CM and MD) complement each other.The MD is derived from an atomistic approach, butis computationally heavy and is restricted to relativelysmall systems. The CM equations are derived assuming
100 200 300 400 500 600Frequency (GHz)012345 T i m e ( n s ) E n e r gy S p ec t r a l D e n s it y ( l og a r it h m i c s ca l e ) FM Energy/Total Energy
Time (dim. less) -3 -2 -1 M od e E n e r g i e s ( d i m . l e ss ) FMn=1n=2n=3 ... 10
FIG. 1. Time evolution of the energy spectral density of acircular graphene sheet obtained by MD simulations for theout-of-plane motion. Parameters are given in the text. Notethe decay of the FM energy which is indicated by the filledarea under the curve, as well as the shift in the FM frequency.The dashed black line is the FM frequency as estimated fromthe CM model. The inset shows the evolution of the individ-ual mode energies obtained from the CM model in logarithmicscale. negligible bending rigidity and long wavelength deforma-tions, but require less computational power and furtherallow to predict scaling behaviors for various physicalparameters. Since the CM model is a long wavelengthapproximation, only the low-frequency modes can be ac-curately described within the model. Additionally, in theCM-simulations presented here only radially symmetricmodes are considered, as the interaction term in Eq. (1)conserves the radial symmetry of the fundamental mode.The dimensionless and parameter free form of the CMequations implies that the dynamics is completely deter-mined by the initial conditions. Further, strong mixingin phase space causes the distribution functions of themodes to decouple. The dynamics of the system is thendescribed by the total dimensionless energy ˜ E and the ra-tio of applied FM energy ˜ E to total energy, η = ˜ E / ˜ E .There will also be a dependence on the number of de-grees of freedom in the system, set by the number ofatoms. In the CM model, this is introduced artificiallythrough the number of modes N that are considered inthe simulations. This number is always much smallerthan the number of atoms in the physical system.The dimensionless energy can be written in terms of T ~0 K V =15 Å/ps V =10 Å/ps V =5 Å/ps V =2 Å/ps T =200 K T =100 K T =300 K ○□ ◇ ✳ Fundamental mode energy (dim. less) F und a m e n t a l m od e fr e q . ( GH z ) (a) + η (0)=1 ○ η (0)=0.98 ✕ η (0)=0.95 △ η (0)=0.85 -12000 -8000 -4000 0 4000 Time (units of periods) N o r m a li ze d e n e r gy η τ m τ tr (b) Effective temperature T * (K) Q =240 Q =2400 Q =24000 Q =240000 R e l a x a ti on r a t e ˜ -5 -4 -3 -2 MD - R =5 nmMD - R =7 nmMD - R =9 nmMD - R =11 nmCM - N =24CM - N =32CM(T) - N =32CM - N =48 (c) FIG. 2. (a) Fundamental mode frequency as a function of mode energy. Symbols correspond to MD results, while the fullline is the curve predicted from CM. The finite length of the time window used to calculate the frequency spectrum limits thefrequency resolution. This is represented by the size of the symbols. (b) Simulation of Eq. (1) for a fixed total energy, butwith varying initial FM energy. By shifting time, the curves align. The solid line is a fitted sigmoid function used to extractthe transition time τ tr . (c) Extracted rates Γ = τ − from MD simulations (filled symbols) and CM (open symbols) for fixedinitial value of η = 1 / T ∗ for fixed (cid:15) pre = 0 . T ∗ . the average energy per atom κ as˜ E = 12 (cid:15) κm c c , (3)where m c is the carbon atom mass. For a system of un-coupled harmonic oscillators in equilibrium, κ = k B T .We study a non equilibrium situation, but an effectivetemperature may be defined by considering the energynot residing in the FM, k B T ∗ ≡ ( E − E ) /N where N isthe number of degrees of freedom. The relation (3) showsthe correspondence between temperature and strain inthis system. The importance of the mode coupling is de-termined by the thermal fluctuations of the membrane.These may be enhanced either by increasing the temper-ature, or by decreasing the strain.The shift of the FM frequency in Fig. 1 can be repro-duced in the CM model. In Fig. 2(a), the energy depen-dence of the FM frequency is reported for the MD simu-lations (symbols), together with the predicted curve fromthe CM model. The frequency shift is a result of non-linearities, and the overall agreement indicates that thenonlinearities are well described within the CM model.Next we consider the dynamics of the FM energy. Theinset of Fig. 1 shows the temporal evolution of the in-dividual mode energies in the CM model, when initiallyall energy is fed into the fundamental mode. Note theappearance of an initial metastable state, with nearly allenergy localized in the fundamental mode. Among modeswith mode number >
3, the energy is always equiparti-tioned, indicating that the assumption of strong mixingis valid. These modes define an instantaneous effectivetemperature of the system, which due to the redistribu-tion increases in time. This effect is more pronounced inthe CM-simulations due to the mode number cut-off.In Fig. 2(b) the normalized FM energy η is monitoredas a function of time, by calculating the mode dynamics from Eq. (1) for a system of N = 32 modes for a fixedtotal energy ˜ E = 3 .
25. The degree of initial excitationover the thermal background, i.e. η , is varied in the sim-ulations. Shifting the time, so that η = 1 / t = 0, thecurves align. This indicates that the important param-eters are the total energy and η . Further, two separatetime-scales associated with the non-equilibrium dynam-ics of the system can be identified; initially, when mostenergy is in the FM, the system relaxes to a long-livedmetastable state with energy dependent life-time τ m . Inthis region, the system is sensitive to fluctuations and therelaxation is therefore strongly temperature dependent.Thereafter the system undergoes a sharp transition to anequipartition region during a time τ tr = Γ − .The MD simulations are best suited to investigatethe shortest of these timescales, the transition time τ tr .We performed MD simulations at various total energieswith initial data chosen such that the system was ini-tially in the transition stage. The velocity applied to theFM was chosen such that η ≈ / ∂ τ η = ξ , ∂ t ˜ E /ω ˜ E . This quantity may be con-sidered as an ”instantaneous Q-factor” of the resonator.The Q-factors reported here are evaluated at η ≈ / .
396 ˚A [30]).The results from the MD simulations (filled symbols)and CM model (open symbols) are compared in Fig. 2(c).There is no dependence of the dimensionless relaxationrate on the size of the drum. This is consistent with theobservation that the dimensionless energy in Eq. (3) doesnot contain any length scale. The transition rate is linearin energy for both models, i.e. Γ ∝ T ∗ /(cid:15) .Note that previous studies on CVD (Chemical VaporDeposition) grown resonators show a power law depen-dence of the Q-factor on radius [6], in contrast to what isreported here. However, for CVD grown graphene grainboundaries will introduce an additional length scale thatdetermines the mode coupling [32].The existence of an initial metastable state for initialconditions far from equilibrium is a well known featureof the FPU-problem [20, 25, 27], and occurs also in othernonlinear lattices [25] (see supplement). The metastablestate is strongly localized in mode space, and is relatedto so called q-breathers [23, 24]. The metastable state isobtained when the energy in the FM is much larger thanthe thermal energy. External losses and thermal noiseare expected to destroy this state, but traces of it maybe found by considering the relaxation as a function ofexcitation energy. Numerical evidence on the FPU-chainsuggests a stretched exponential dependence of the life-time with excitation energy, τ m h exp (cid:16) ˜ E − α (cid:17) [27, 36].We probe the metastable state by initially feeding allenergy into the FM, and monitoring the evolution of theFM energy. As the time-scales for the decay of this stateis very long compared to the FM oscillations, this stateis most readily investigated using the CM model. Forthe present model, we obtain an exponent α ≈ .
18 fromCM simulations with N = 32 and N = 40 modes (seeFig. 3). The life-time of the metastable state does notshow a dependence on the number of modes. Note thatthe exponent α is model dependent [37], e.g. for theFPU- β -problem α = 0 .
25 has been reported [38].The energy is strongly localized to the FM during themetastable state. The localization of the metastable statein mode space depends on the frequency spectrum. Forthe drum resonator, the frequencies can be approximatedby ω n, drum ≈ ξ , + nπ . In comparison, for small wavenumbers the frequencies of the FPU chain are given by ω n, FPU ≈ nπ [39]. Low frequencies of the FPU chain arealmost resonant, and so the formation of a q-breatherwill consist of a cascade of modes being excited. Thepresent system, being far from resonant, displays a stronglocalization of energy.The CM model allows us to make quantitative predic-tions for the relaxation rate due to energy redistributionin resonators of arbitrary size and tension. We take asan example a drum of arbitrary radius, (cid:15) pre = 0 .
2% and
Dimensionless energy ˜ E M e t a s t a b l e li f e t i m e τ m N=40N=32 --- α =0.18 α =0.25 FIG. 3. Dependence of the metastable lifetime τ m on sys-tem energy, for the CM model with 32 modes (triangles) and40 modes (squares). The lines corresponds to stretched ex-ponentials of the form τ m h exp (cid:16) ˜ E − α (cid:17) . The dashed linecorresponds to α = 0 .
18, the dotted corresponds to α = 0 . T = 4 K. The relaxation rate can be read off from Fig.2(c), giving a transition time of approximately 24000 os-cillations, decreasing to ∼
200 oscillations at 300 K. Forthis relaxation to be observable, the Q-factor arising fromexternal losses must exceed these values.In summary, the redistribution of energy in grapheneresonators due to nonlinear mode coupling has been in-vestigated numerically in a ring-down setup. This cou-pling is a limiting factor for the stability of excitationsof individual modes. At low temperatures, we find evi-dence for a state akin to the metastable state found inthe FPU problem. After a time increasing exponentiallywith inverse excitation energy, the system relaxes towardits equilibrium with an energy dependent rate.The rate of relaxation from a strong non-equilibriumsituation with all energy in the FM to a situation closeto equipartition depends only on the dimensionless en-ergy of the resonator, which in turn depends on tem-perature and strain through the ratio T ∗ /(cid:15) , but noton resonator size. The Q-factors obtained in this Letterare comparable in magnitude to experimentally observedQ-factors [4]. Therefore, it should be possible to experi-mentally verify the proposed mechanisms.Since the system is closed, the rates reported here con-stitute a lower bound on the dissipation of graphene res-onators. Many applications of nanomechanical systemsrely on high Q-factors [40, 41], and therefore an improvedunderstanding of the physical processes limiting this is offundamental interest. Our results demonstrate the pos-sibility of using graphene resonators as a test bed forFPU physics. By this approach long standing questionsin non-equilibrium statistical mechanics might eventuallybe within experimental reach.The authors acknowledge funding from the Euro-pean Union through FP7 project no. 246026 (RODIN)(DM, AI), the Swedish Research Council (DM, AC),the Knut & Alice Wallenberg foundation (DM, AI) andNSF CMMI-1036460 (ZQ,HSP). The authors acknowl-edge Prof. David Campbell for valuable input. ∗ [email protected][1] J. S. Bunch, A. M. van der Zande, S. S. Verbridge, I. W.Frank, D. M. Tanenbaum, J. M. Parpia, H. G. Craighead,and P. L. McEuen, Science , 490 (2007).[2] A. Eriksson, S. Lee, A. A. Sourab, A. Isacsson, R. Kau-nisto, J. M. Kinaret, and E. E. B. Campbell, Nano Lett. , 1224 (2008).[3] C. Chen, S. Rosenblatt, K. I. Bolotin, W. Kalb, P. Kim,I. Kymissis, H. L. Stormer, T. F. Heinz, and J. Hone,Nat. Nanotechnol. , 861 (2009).[4] A. Eichler, J. Moser, J. Chaste, M. Zdrojek, I. Wilson-Rae, and A. Bachtold, Nat. Nanotechnol. , 339 (2011).[5] Y. Oshidari, T. Hatakeyama, R. Kometani, S. Warisawa,and S. Ishihara, Appl. Phys. Express , 117201 (2012).[6] R. A. Barton, B. Ilic, A. M. van der Zande, W. S. Whit-ney, P. L. McEuen, J. M. Parpia, and H. G. Craighead,Nano Letters , 1232 (2011).[7] A. M. v. d. Zande, R. A. Barton, J. S. Alden, C. S.Ruiz-Vargas, W. S. Whitney, P. H. Q. Pham, J. Park,J. M. Parpia, H. G. Craighead, and P. L. McEuen, NanoLetters , 4869 (2010).[8] K. Eom, H. S. Park, D. S. Yoon, and T. Kwon, PhysicsReports , 115 (2011).[9] J. Atalaya, A. Isacsson, and J. M. Kinaret, Nano Lett. , 4196 (2008).[10] M. H. Matheny, L. G. Villanueva, R. B. Karabalin, J. E.Sader, and M. L. Roukes, Nano Letters , 1622 (2013).[11] A. M. Eriksson, D. Midtvedt, A. Croy, and A. Isacsson,Nanotechnology , 39570 (2013).[12] R. Lifshitz and M. L. Roukes, Phys. Rev. B , 5600(2000).[13] M. C. Cross and R. Lifshitz, Phys. Rev. B , 085324(2001).[14] I. Wilson-Rae, Phys. Rev. B , 245418 (2008).[15] L. G. Remus, M. P. Blencowe, and Y. Tanaka, Phys.Rev. B , 174103 (2009).[16] A. Croy, D. Midtvedt, A. Isacsson, and J. M. Kinaret,Phys. Rev. B , 235435 (2012). [17] M. Imboden and P. Mohanty, Physics Reports , (inpress).[18] A. W. Barnard, V. Sazonova, A. M. van der Zande, andP. L. McEuen, Proceedings of the National Academy ofSciences , 19093 (2012).[19] E. Fermi, J. R. Pasta, and S. Ulam, Studies of nonlinearproblems , Tech. Rep. (Los Alamos Scientific Laboratory,1955).[20] D. K. Campbell, P. Rosenau, and G. Zaslavsky, Chaos(2005).[21] F. Izrailev and B. V. Chirikov, Soviet Physics Doklady , 30 (1966).[22] N. Zabusky and M. D. Kruskal, Physical Review Letters , 240 (1965).[23] T. Penati and S. Flach, Chaos (2007).[24] S. Flach, M. V. Ivanchenko, and O. I. Kanakov, Phys.Rev. Lett. , 064102 (2005).[25] F. Fucito, Marchesoni, F. Marinari, G. Parisi, L. Peliti,and S. Ruffo, Le Journal de Physique (1982).[26] D. Bambusi and A. Ponno, in The Fermi-Pasta-Ulamproblem (Springer, 2008).[27] G. Benettin, A. Carati, L. Galgani, and A. Giorgilli, in
The Fermi-Pasta-Ulam problem (Springer, 2008).[28] Lammps, http://lammps.sandia.gov (2012).[29] S. J. Plimpton, Journal of Computational Physics ,1 (1995).[30] S. J. Stuart, A. B. Tutein, and J. A. Harrison, Journalof Chemical Physics , 6472 (2000).[31] H. Zhao, K. Min, and N. R. Aluru, Nano Letters , 3012(2009).[32] Z. Qi and H. S. Park, Nanoscale , 3460 (2012).[33] Z. Qi, F. Zhao, X. Zhou, Z. Sun, H. S. Park, and H. Wu,Nanotechnology , 265702 (2010).[34] A. J. Lichtenberg, R. Livi, and M. Pettini, in The Fermi-Pasta-Ulam problem , edited by G. Gallavotti (Springer,2008).[35] C. Lee, X. Wei, J. W. Kysar, and J. Hone, Science ,385 (2008).[36] G. Benettin, L. Galgani, and A. Giorgilli, Nature ,444 (1984).[37] M. Pettini and M. Landolfi, Phys. Rev. A , 768 (1990).[38] L. Berchialla, A. Giorgilli, and S. Paleari, Physics LettersA , 167 (2004).[39] G. P. Berman and F. M. Izrailev, Chaos (2005).[40] B. Ilic, Y. Yang, and H. Craighead, Applied PhysicsLetters , 2604 (2004).[41] M. D. LaHaye, O. Buu, B. Camarota, and K. C. Schwab,Science , 74 (2004). upplementary information for: FPU physics with nanomechanical grapheneresonators: intrinsic relaxation and thermalization from flexural mode coupling Daniel Midtvedt, ∗ Alexander Croy, and Andreas Isacsson
Department of Applied Physics, Chalmers University of Technology, SE-412 96 G¨oteborg, Sweden
Zenan Qi and Harold S. Park
Department of Mechanical Engineering, Boston University, Boston, MA 02215, USA
DERIVATION OF COUPLED MODE EQUATIONS
The dynamics of suspended graphene resonators are described within continuum mechanics by the F¨oppl-vonKarmann equations [1]. For the drum geometry (Fig. 1), they read [2] ρ G ¨ u r − (cid:2) ∂ r σ rr + r − ∂ φ σ rφ + r − ( σ rr − σ φφ ) (cid:3) = 0 , (1a) ρ G ¨ u φ − (cid:2) ∂ r σ rφ + 2 r − σ rφ + r − ∂ φ σ φφ (cid:3) = 0 , (1b) ρ G ¨ w − (cid:15) pre ( λ + 2 µ ) r − ∂ r ( r∂ r w ) − r − (cid:2) ∂ r ( rσ rr ∂ r w + σ rφ ∂ φ w ) + ∂ φ (cid:0) σ rφ ∂ r w + r − σ φφ ∂ φ w (cid:1)(cid:3) = 0 , (1c)where u r , u φ and w are the radial, angular and out-of-plane displacement fields, σ ij is the stress tensor, ρ G is themass density of graphene, (cid:15) pre is the pre-strain of the sheet and λ and µ are the Lam´e parameters of graphene. Toderive the equation of motion for the mode dynamics, the in-plane dynamics is first eliminated adiabatically by setting¨ u r = 0 and ¨ u φ = 0 and solve the remaining stationary equations for the in-plane fields. The equations are madedimensionless through the following scaled variables,˜ r = rR ; ˜ t = √ (cid:15) pre c L R t ; ˜ w = 12 w R (cid:15) pre ; ˜ E = E (cid:15) c ρ G R , (2)where c = ( λ + 2 µ ) /ρ G is the longitudinal speed of sound in graphene and R is the radius of the drum. In thefollowing, only radially symmetric deformations are considered.The dimensionless equations for the in- and out-of-plane components are, skipping the tildes over the scaled variablesfor convenience, (cid:20) u r,rr + 1 r u r,r − r u r (cid:21) − − ν r w ,r − ∂ r w ,r = 0 , (3a)¨ w − (cid:18) ∂ r + 1 r (cid:19) w ,r + w ,r ! − (cid:18) ∂ r + 1 r (cid:19) u r,r w ,r − ν r ∂ r u r w ,r = 0 , (3b)with boundary conditions w (1 , t ) = 0, u r (1 , t ) = 0. Subscript , r denotes radial differentiation and ν is the Poissonratio. The in-plane field is expanded as u r ( r, t ) = P k Q k ( t ) J ( ξ k r ) where ξ nk is the k:th zero of the n:th Besselfunction, J n ( ξ nk ) = 0. The amplitudes Q k ( t ) are given by Q k = J ( ξ k ) ξ k C k ( t ) , (4)where C k ( t ) ≡ Z drrJ ( ξ k r ) (cid:20) − ν r w ,r + 12 ∂ r w ,r (cid:21) . (5) ∗ [email protected] a r X i v : . [ c ond - m a t . m e s - h a ll ] A p r Inserting this result into the equation for w we obtain¨ w − (cid:18) ∂ r + 1 r (cid:19) w ,r + w ,r ! − X k J ( ξ k ) C k ξ k (cid:20)(cid:18) ∂ r + 1 r (cid:19) J ( ξ k r ) w ,r + νr ∂ r ( J ( ξ k r ) w ,r ) (cid:21) = 0 . (6)Next we expand w ( r, t ) = P n q n ( t ) J ( ξ n r ) ≡ P n q n ( t ) φ n ( r ) for brevity. The equation of motion then becomes2 J ( ξ n ) ¨ q n + 2 J ( ξ n ) ξ n q n + 12 X k,l,m M klmn q k q l q m − X k,l J ( ξ k ) C k ξ k K kln q l , (7)with M klmn = Z drrφ k,r φ l,r φ m,r φ n,r , (8a) K kln = − Z dr [ rJ ( ξ ,k r ) φ n,r φ l,r + νJ ( ξ ,k r ) φ n,r φ l,r ] . (8b)Inserting the expansion of w into the definition of C k we find C k = 12 X ij K kij q i q j . (9)Consequently, the equation of motion for the out-of-plane motion becomes¨ q n + ξ n q n + X i,j,m W ij ; mn q i q j q m = 0 , (10)with W ij ; mn = J ( ξ n ) M ijmn − J ( ξ n ) X k J ( ξ ,k ) ξ ,k K kim K kjn . (11)being the effective coupling matrix, which enters Eq. (1) in the main text. From the definition, it is clear that M ijmn is symmetric in all indices. The second term in the definition of the coupling matrix is symmetric under i ↔ j and m ↔ n , which is reflected in the notation W ij ; mn . An alternative derivation based on the Airy stress function isoutlined in [2]. MAGNITUDE OF COUPLING
Since we are interested in the dynamics of the fundamental mode, the structure of the matrix can be visualized byfixing the first index to i = 0, and loop over the remaining indices, treating ij and mn as superindices. From Fig. 2,it is clear that the strength of the coupling as well as the density of modes to which the fundamental mode couplesincreases with mode number. METASTABLE STATES IN FPU SYSTEMS
One of the striking features of the FPU problem is the existence of a long-lived metastable state far from equipar-tition. An initial understanding of the apparent lack of equipartition was based on the similarity of the FPU modelto the Korteweg-de Vries (KdV) equation [3]. The KdV equation is integrable, and the lack of equipartition wasattributed to this fact. However, the metastability persists only for energies below a certain threshold. Above thethreshold, the metastable state disappears and the system approaches equipartition on a relatively short time scale [4].The metastability is therefore not completely explained by the similarity of the FPU problem to the integrable KdVequation. A complementary approach to describe the metastable state was developed by Flach and coworkers [5],which focusses on the energy spectrum in mode space (q-space). For vanishing nonlinearity, any excitation consisting 𝑅 FIG. 1. Schematic image of the simulated structure.
Index m,n I nde x i , j log W ij ; mn FIG. 2. Visualization of the nonlinear mode coupling to the fundamental mode. The axes provide the mode number combina-tions. The coupling generally increases with mode number. of only one mode is a periodic orbit. For nonzero nonlinearity, the orbit will spread out to nearby modes, constitutingwhat is denoted a q-breather. These q-breathers are exponentially localized in mode space, and are obtained by treat-ing the FPU system as a perturbation to a system of uncoupled harmonic oscillators. Consequently, the q-breatherformalism applies to generic nonlinear lattices as long as a non-resonance condition is fulfilled and for sufficiently lowenergies. The exponential localization in mode space is present also in the system considered in this paper, see Fig.3. Here, the four lowest modes are initialized with equal energy. The localization length shows a clear dependence onthe excitation energy, and in the limit of large energies the spectrum is essentially flat, signifying equipartition.
Mode number M ode ene r g y ( no r m a li z ed ) − − − − E = 2 E = 10 E = 5 · FIG. 3. Mode energy distribution during the metastable state. The lowest four modes where initialized with equal energy. Theenergy is exponentially localized in mode space, with an increasing localization length with energy. Blue squares correspond tototal energy E = 5 · − , black circles to E = 10 − and green diamonds to E = 2. The mode energies are normalized to theenergy of the fundamental mode.[1] J. Atalaya, A. Isacsson, and J. M. Kinaret, Nano Lett. , 4196 (2008).[2] A. M. Eriksson, D. Midtvedt, A. Croy, and A. Isacsson, Nanotechnology , 39570 (2013).[3] N. Zabusky and M. D. Kruskal, Physical Review Letters , 240 (1965).[4] F. Izrailev and B. V. Chirikov, Soviet Physics Doklady , 30 (1966).[5] S. Flach, M. V. Ivanchenko, and O. I. Kanakov, Phys. Rev. Lett.95