Snap-through transition of graphene membranes for memcapacitor applications: A combined study using MD, DFT and elasticity theory

Using computational and theoretical approaches, we investigate the snap-through transition of buckled graphene membranes. Our main interest is related to the possibility of using the buckled membrane as a plate of capacitor with memory (memcapacitor). For this purpose, we performed molecular-dynamics (MD) simulations and elasticity theory calculations of the up-to-down and down-to-up snap-through transitions for membranes of several sizes. We have obtained expressions for the threshold switching forces for both up-to-down and down-to-up transitions. Moreover, the up-to-down threshold switching force was calculated using the density functional theory (DFT). Our DFT results are in general agreement with MD and analytical theory findings. Our systematic approach can be used for the description of other structures, including nanomechanical and biological ones, experiencing the snap-through transition.


Introduction
Memcapacitors 1 are an emerging type of circuit elements with memory whose instantaneous response depends on the internal state and input signal. Such devices are prospective candidates for applications in information storage and processing 2, 3 technologies as their states can be manipulated by the applied voltages or charges and can store information for long intervals of time. Several possible realizations of memcapacitors were suggested by using micro-electro-mechanical systems 4 , ionic transport 5 , electronic effects 6 , superconducting qubits 7 , etc. 8 Generally, voltage-controlled memcapacitive systems (memcapacitors) are described by 1 where q(t) is the charge on the capacitor at time t, V (t) is the applied voltage, C is the memcapacitance (memory capacitance), x is a set of n internal state variables, and f is a continuous n-dimensional vector function. In some cases, it is more convenient to consider charge-controlled memcapacitors 1 such that the charge instead of voltage is considered as input.
Among several possible realizations of memcapacitors, the membrane-based memcapacitors 4 are of significant interest as their geometry makes them intrinsically suitable for non-volatile storage of binary information. Indeed, the buckled membrane used as the top capacitor plate (see Fig. 1 for schematics) has two stable buckled states corresponding to two distinct values of capacitance. It was suggested 4 that the switching between these states can be performed using the attractive interaction of oppositely charged capacitor plates. Moreover, it was demonstrated theoretically that simple circuits of membrane memcapacitors offer an in-memory computing functionality 3 .
In this work, we consider a possible realization of membrane-based memcapacitor 4 employing a single-or multi-layer graphene membrane 9-13 as its bistable plate (see Fig. 1). Our aim is to understand the basic physical processes and parameters underlying the snap-through transition of such membrane including details of the membrane dynamics, threshold forces, etc. For this purpose, we perform a combined study using MD, DFT and elasticity theory focusing on a single-layer graphene membrane with clamped boundary conditions. This choice of boundary conditions is justified by the typically strong adhesion of graphene to substrates. Our results extend our prior DFT investigation 14 of the up-to-down snap-through transition of graphene membrane with hinged boundary conditions. The combination of computational/theoretical methods adds breadth and depth to our analysis. Using MD simulations, we were able to understand main features of the membrane dynamics in the presence of an external force and after the force removal. This understanding has helped us to develop analytical models that resulted in compact algebraic expressions for the threshold switching forces. DFT calculations were used to validate MD results for the up-to-down transition.
This paper is organized as follows. In Sec. "Molecular Dynamics Simulations" we investigate the snap-through transition of graphene membranes using molecular dynamics simulations. In particular, MD simulations of the up-to-down and down-to-up transitions are reported in Subsec. "Up-to-down transition" and "Down-to-up transition", respectively, while MD simulation details can be found in Supplementary Information (SI) Sec. "MD Simulation details". The standard elasticity theory is applied to the membrane switching in Subsec. "Buckling and snapping-through within the theory of elasticity". A phenomenological analytical model of the snap-through transition is presented in Subsec. "Phenomenological elasticity theory" and in SI Sec. "Phenomenological Elasticity Theory". Our DFT calculations are summarized in Sec. "Density Functional Theory". The results obtained within different approaches as well as their implications are discussed in Sec. "Discussion".

Molecular Dynamics Simulations
MD simulations are a well established modeling tool frequently employed in studies of nanoscale carbon-based materials  . We used classical MD simulations to investigate the dynamics of snap-through transition of buckled graphene membranes. Zigzag graphene nanoribbons (membranes) of two lengths were considered: the nanoribbon A, L = 54Å (22 rings), and nanoribbon B, L = 103Å (42 rings). Both nanoribbons were of the same width (w = 41Å). In order to implement the clamped boundary conditions, the first two lines of carbon atoms at shorter sides were kept fixed. The buckling was realized by changing  the distance between the fixed sides from L to d < L. The same external force was applied to each atom in the downward direction to simulate the attractive interaction between the plates. See SI Sec. "MD Simulation details" for the details of MD simulations.

Up-to-down transition
In order to simulate the up-to-down transition, the double-clamped graphene nanoribbon buckled upwards was subjected to the force in the downward direction. The final state of membrane was found using MD simulations as described in SI Sec. "MD Simulation details". Figure 2 shows the final position of a central atom of membrane versus the applied force for several values of d/L and two membrane sizes. According to Fig. 2, at fixed L, the up-to-down threshold switching force (the minimal force required for the up-to-down transition) is larger for smaller values of d/L. Moreover, at fixed d/L, the threshold switching force is smaller for longer membranes. Additionally, some curves in Fig. 2 exhibit a noisy threshold (such as d/L = 0.8 and d/L = 0.98 in (a)), which can be related to thermal fluctuations of membranes. All these observations are intuitively reasonable. Depending on the force magnitude, the membrane switching occurs either through the symmetric or non-symmetric membrane profile (see Fig. 3). Figure 3 was obtained using overdamped simulations of membrane dynamics (additional results of these simulations can be found in SI Sec. "MD Simulation details"). The non-symmetric profile is associated with a smaller energy barrier 14 and involved in switching by smaller forces. Larger applied forces result in the switching through the symmetric profile. According to our observations and previous work 14 , in all cases, the membrane profile is symmetric at short times. If the force magnitude is sufficient to overcome the energy barrier for the symmetric profile, then, typically, the switching takes place through the symmetric path. Otherwise, a symmetry breaking occurs leading to the switching through the lower energy barrier associated with the non-symmetric membrane profile. Thermal fluctuations help the symmetry breaking.

Down-to-up transition
It was suggested in Ref. 4 that the memcapacitor membrane can be set into the buckled upwards state 0 (see Fig. 1) by also using the electrostatic attraction between the capacitor plates. When the pulled-down membrane is suddenly released, there are conditions such that the membrane overcomes the potential barrier and ends up in the buckled upwards state. As the kinetic energy plays an important role in the down-to-up transition, this process can not be analyzed using the overdamped dynamics. In order to simulate the down-to-up transition, every atom of a double-clamped buckled membrane was subjected to a constant force in -z direction. After an equilibration period, the forces were removed, and the system was simulated for a time interval sufficient to reach the steady state. The final position of a central atom of membrane is presented in Fig. 4 as a function of the applied force for several values of d/L. Figure 4 shows that the threshold forces for the down-to-up transition are larger than those for the up-to-down transition (see Fig. 2). Moreover, the threshold forces for the down-to-up transition are larger for the shorter membrane A compared to the longer membrane B. It is interesting that the final state of membrane oscillates as a function of the applied force (see d/L = 0.98 curves in Fig. 4). Qualitatively, having a high kinetic energy, the membrane moves up and down while its energy dissipates. Finally, it gets trapped in one of two stable buckled states.
Overall, our MD results indicate the possibility of writing one bit of information into the state of membrane memcapacitor. In this interpretation, logic 1 corresponds to the membrane buckled downwards, and logic 0 -to the membrane buckled upwards. In addition to 0 → 1 and 1 → 0 transitions presented in Figs. 2 and 4, we have verified the occurrence of 1 → 1 and 0 → 0 transitions at the same values of forces required for 0 → 1 and 1 → 0 transitions, respectively. Therefore, logic 0 or 1 can be written into the memcapacitor just by selecting a suitable value of the applied force.
In our work, we considered relatively small nanoribbons, which are mechanically rigid. This is an important requirement for the memcapacitor application as the mechanical nonvolatile information storage is not possible with flexible graphene. From some previous studies 37 it is known that local structural corrugations (ripples) 38 disappear in the transition from flexible to rigid graphene. In agreement with this previous work 37 , our molecular dynamics simulations have shown the absence of ripples (spontaneous height fluctuations) in nanoribbons of reported sizes and their existence in larger nanoribbons. The latter, however, are not suitable for the memcapacitor application because of their flexibility. To summarize, while our simulation tools provide a means for ripples detection, we did not observe these in the reported structures that are mechanically rigid.

Theory of Elasticity
Even though the graphene has a thickness of one atom, the graphene membranes are quantitatively good described by the theory of elasticity [39][40][41][42] . This allows us, on one hand, to obtain analytical expressions for the buckled membranes [43][44][45] and, on the other hand, perform a comparison with MD simulations.
There is a number of publications dealing with buckled beams and membranes under the transverse load (see, for example, Refs. [46][47][48][49]). Such systems are frequently described in the framework of Euler-Bernoulli theory, which, however, leads to complex analytically unsolvable equations. Bubnov-Galerkin decomposition is a one of the best approaches to find the approximate analytical solution of these equations, although it still requires a complex phase-diagram analysis. In particular, using such analysis, the authors of Refs. [46][47][48][49] investigated an electrostatically loaded double-clamped membrane above a rigid flat electrode and derived cumbersome conditions for the symmetric snap-through transition, symmetry breaking, existence of bifurcation and pull-in instability.
In this Section, we derive compact but sufficiently precise expressions for the snap-through switching forces based on the theory of elasticity.

Description of buckled membranes within the theory of elasticity
Consider a 2D membrane within the theory of elasticity 40,41,50 . The total potential energy of membrane is defined by the deflection ζ (along the normal direction z) as Here, ∆ stands for the 2D Laplacian, T = T 0 + δ T is the total tensile force, T 0 is the force applied by the support and δ T is the bending tension resulted from the extension, S 0 = wL, and F is the external force density. The compression of membrane corresponds to T 0 < 0. Given the potential energy, one may also be interested in the dynamics of membrane, which is defined by the equation where µ is the 2D mass density.

Eigenmodes and buckling
The spatial harmonics (eigenmodes) of clamped membranes can be written as where the numbers b n are the solutions of the equation cosh b n cos b n = 1. One can find that b n ≈ 3π/2 + nπ (in particular, b 0 = 4.73 is close to 3π/2 = 4.65). The functions ζ n (x) are orthonormal. Consider a buckled membrane shown in Fig. 1. Its potential energy (Eq. (3)) calculated for the fundamental n = 0 mode where , and f = πFdw/6. Here, T c = 4π 2 D/d 2 is the critical tension. At α > 0, that is at |T 0 | > T c , the potential is the double-well potential with minima at symmetric deflections, ζ c = ± α/2β , and the potential barrier δU = α 2 /4β . (We note that a less realistic case of hinged boundary conditions would lead to simpler expressions for eigenmodes, ∼ sin (πnx/d), compared to Eq. (6), making the hinged boundary conditions preferable for some calculations. However, this case results in quantitatively different values. For example, the critical tension differs four times.) We note that the nanoribbon length is given by

5/22
Expanding this expression to the second order, one obtains the interrelation between L, d, and the maximal deflection ζ c : On the other hand, the maximal deflection of buckled membrane can be inferred from the minima of the potential energy, Eq. (7): Equations (9) and (10) link the tensile force T 0 to the parameter d/L, which can be considered as a measure of deformation of buckled membrane. These equations can be used to express the tensile force T 0 through d/L.

Snap-through transition
Here we present a convenient method to describe the dynamics of membrane subjected to an external force and find the minimal force causing its snap-through transition. The main idea is to expand ζ (x,t) in harmonics ζ n (x) (given by Eq. (6)) with amplitudes q n (t) as and limit the sum to few first terms. We found that in order to describe the symmetric and non-symmetric transitions, it is sufficient to consider n = 0, 2 and n = 0, 1 terms in Eq. (11), respectively. The calculation consists in the following. First, we substitute the expansion (11) in Eq. (5). Second, the resulting equation is multiplied by ζ m (x) and integrated taking into account the orthogonality of harmonics. Here we emphasize that for the harmonics (6), the integrals are readily calculated, and instead of the integro-differential equation one obtains a system of differential equations for the functions q n (t). Inserting the obtained functions q n (t) back into Eq. (11) gives us the description of membrane dynamics. In the following, the results of these calculations (performed using the time normalized by the membrane characteristic frequency ω c = (1/L 2 ) D/µ and t c = 1/ω c ) 51 are analyzed for membrane B.
Consider the results presented in Fig. 5. Fig. 5(a) demonstrates the case of fully symmetric switching, when the dynamics can be described by n = 0 and 2 harmonics. Next, we introduce a small asymmetry into the problem via a small non-symmetric modification of the force (of the order of 0.1%). This results in the non-symmetric dynamics of Fig. 5(b), which can be described in terms of n = 0 and 1 harmonics. We note that the force needed for the non-symmetric snap-through transition is smaller than that for the symmetric one, which will also be analyzed below in more detail. Finally, in Fig. 5(c) we illustrate a combination of the above regimes found using a symmetric force, when a tiny fluctuation in the numerical solution changes the symmetric dynamics to the non-symmetric one. Note that this case is close to the one discussed in Ref. 14. Figure 6(a) depicts the dependence of the final position of membrane center, ζ c (t f ), on the applied force F. At a certain value of force, F = F ↓ , the membrane switches from the buckled upwards state to the buckled downwards state. Our calculations show that the threshold force is proportional to the initial central-point deflection ζ c and is different for the symmetric (s) and non-symmetric (ns) dynamics. Namely, this force, being multiplied by the membrane area wL, reads Since the non-symmetric transition is energetically favorable, we expect that this value of force is the one to be used in the device design/experiment analysis. We note that estimations based on Eq. (12) are in a good agreement with results found by other methods as we discuss later and illustrate in Fig. 8. Figure 6(b) presents the time-dependence of harmonics' amplitudes q n (t) in the up-to-down transition. Importantly, the dynamics is well described by n = 0, 1, and 2 harmonics, while n = 3 and 4 harmonics amplitudes are negligibly small. Our calculations (similar to the consideration of the up-to-down force above) demonstrate that the minimal force needed for the down-to-up transition is proportional to ζ 2 c as See also SI Sec. "Down-to-up transition", where the down-to-up transition is illustrated graphically.

Phenomenological elasticity theory
In this Section, in order to avoid all complications associated with the Bubnov-Galerkin decomposition, we elaborate a different approach describing the up-to-down and down-to-up transitions in a compact analytical form.

Ansatz
Following our observation of the existence of the symmetric and non-symmetric membrane profiles in the snap-through transition (see Fig. 3), we introduce two simplest polynomial functions to describe such symmetric (s) and non-symmetric (ns) shapes of membrane: Here θ i (s) is the angle that the membrane makes with the horizontal, s is the internal coordinate that changes between −1/2 and 1/2, A i and c i are coefficients, and i = {s, ns}. Clearly, Eqs. (16) and (17) describe the double-clamped membrane as θ i (±1/2) = 0.

Up-to-down transition
Our goal here is to estimate the minimal force required for the snap-through transition. In this subsection, we assume that the snap-through transition is induced by a slowly increasing force such that the system always stays in the potential energy minimum. At zero applied force, there are two minima of the potential energy corresponding to the two stable states of membrane. The applied force modifies the potential energy landscape such that, starting at a certain force, the potential energy has a single minimum. Here, this value of force is found and considered as an estimation for the threshold switching force. In the presence of a constant force F in the downward direction, the membrane potential energy is given by where U b is the bending energy, U str is the stretching energy, U ext is the potential energy due to the applied force. Then, using U ext = Fz cm , where z cm is z coordinate of the center of mass, and neglecting U str term not important for the up-to-down transition (see, e.g., Fig. S1 that indicates the insignificance of U str in buckled membrane), we get the following equation for the potential energy extrema: where i = {0, 1} corresponds to {s, ns}, respectively. Equation (19) allows finding c 0 and c 1 for a given strength of the applied force.
The minimal force needed for the snap-through transition corresponds to the maximum of F(z cm ) and, for the transition through the symmetric shape, is given by taking place at c 0 = 0.3589 (the corresponding membrane profile is presented in Fig. S4). At this value of c 0 , the bending energy is Performing the same calculations for the non-symmetric shape, one can find that the force needed to support the equilibrium non-symmetric shape decreases with the progress of switching. A zero force is reached at c 1 = −c 2 = 1/(2 √ 5) that corresponds to the maximum of U b,ns = 90Dw(L − d)/L 2 . A possible (rough) estimation for the threshold switching force can be obtained taking the limit c 1 → ∞ leading to and U b,ns = 42Dw(L − d)/L 2 .
In fact, a realistic scenario of the switching through the non-symmetric shape can be imaged as follows. We start with a symmetric membrane at F = 0 and slowly increase the force. The symmetry breaking occurs at a certain value of force, say, F ↓ . The dynamics of switching is a complex process significantly relying on thermal fluctuations. As the switching dynamics can

8/22
not be reached within the framework of our model, for estimation purposes, we assume that F ↓ corresponds to the point where the bending energies and applied forces for the symmetric and non-symmetric shapes are the same. One can find that both conditions are satisfied at c 0 = 0.4683 and c 1 = 0.6948, so that U b = 47.8Dw(L − d)/L 2 and Consequently, at F < F ↓ , the membrane keeps the symmetric shape and switches to the non-symmetric one as soon as F reaches F ↓ . As there is no barrier involved in the non-symmetric switching, no further increase in the applied force is required to complete the snap-through transition through the non-symmetric shape.

Down-to-up transition
As the stretching energy U str (see Eq. (18)) plays an important role in the down-to-up transition, it needs to be taken into account. For our purposes, it is sufficient to approximate U str as where ∆L is the change in the membrane length. In the approximation of small elongation, ∆L L, one can assume that both U b and the shape of the stretched membrane are not significantly modified compared to F = 0 case. In this situation, at the threshold of transition, the stretching energy (Eq. (23)) is equal to the difference between the maximal and minimal bending energies (the energy conservation condition). Using the energies given below Eq. (S8) we find Moreover, at equilibrium Under the condition of small ∆L L, the center of mass position of membrane (buckled downwards) can be expressed as .
Using Eqs. (24) and (27), Eq. (25) can be rewritten as We have verified that the threshold switching force for the down-to-up transition given by Eq. (28) is in agreement with our typical MD simulations results. A very good agreement was obtained with MD simulations performed with very small energy dissipations.

Density Functional Theory
A DFT investigation of the up-to-down transition was performed following the previously reported study of the up-to-down switching of buckled graphene membrane 14 . In the present calculations, we focused on the effects of boundary conditions and buckling strength on the up-to-down threshold switching force. DFT computational details can be found in SI Sec. "DFT Calculations". We investigated a buckled zigzag graphene nanoribbon specified in SI. Two types of boundary conditions were examined: 1) hinged boundary conditions (configurations C1 with d/L = 0.96, Fig. 7(a), C2 with d/L = 0.78, and C3 with d/L = 0.63), and 2) clamped boundary conditions (configuration D1 with d/L = 0.96, Fig. 7(b) and (c)). The initial geometries of nanoribbons corresponded to the buckled upwards state. In order to simulate the up-to-down transition, we used the approach developed in Ref. [14]. Specifically, we performed a series of DFT calculations with the central chain of membrane atoms (marked red in Figure 7. Geometries of membranes C1 (hinged boundary conditions, (a)) and D1 (clamped boundary, (b) and (c)) in the process of the up-to-down switching. The switching occurs through non-symmetric paths (a) and (b), and initially symmetric path (c). These geometries were found using DFT optimization of geometries obtained utilizing a progressive displacement of central atoms. Fig. 7) gradually displaced in −z direction from its position in the buckled upwards state. The membrane geometry found at the preceding deformation step was used to build its subsequent modification. In contrast to Ref. [14], in the present calculations only z coordinates of the central chain of atoms were constrained.
The deformation energies of buckled graphene are reported in SI. Our force estimation shows that the threshold switching force for the up-to-down transition is about 3.8 pN/atom for C1, 11.1 pN/atom for C2, and 16.4 pN/atom for C3. Moreover, the non-symmetric switching force for D1 is 3.7 pN/atom and the symmetric one is about 11.5 pN/atom. As the non-symmetric switching forces for D1 and C1 are close to each other, we use the results for C1-C3 for the order-of-magnitude comparison with MD and elasticity theory results.

Four methods to describe the snap-through transitions
As the key finding, Figure 8 presents the dependence of the up-to-down threshold switching force (calculated with different approaches) on the membrane deformation parameter d/L. In particular, this figure demonstrates that the threshold forces obtained by a number of different methods are in good agreement with each other.
We emphasize that: (i) MD is widely used approach to describe the dynamics of nanoscale systems. In our work, MD clearly shows the possibility of 0 → 1 and 1 → 0 transitions of graphene membrane. We have also verified the occurrence of 1 → 1 and 0 → 0 transitions at the same values of forces required for 0 → 1 and 1 → 0 transitions, respectively.
(ii) In the framework of elasticity theory, the membrane is described by means of an integro-differential equation for the deflection ζ (x,t). Expanding the function ζ (x,t) into membrane's harmonics allows reducing the integro-differential equation to a system of differential equations. This approach was utilized in Sec. Buckling and snapping-through within the theory of elasticity. There it was shown that the up-to-down threshold switching force depends linearly on the initial central-point deflection, F ↓ ∝ ζ c , while the down-to-up snap-through force at low dissipation depends quadratically on the initial deflection, F ↑ ∝ ζ 2 c . Besides, we have demonstrated that with a high accuracy, the relevant membrane's dynamics can be described (and visualized) just by two harmonics.
(iii) The phenomenological approach based on the theory of elasticity has allowed us to obtain analytical expressions for the up-to-down and down-to-up switching forces for buckled graphene membrane. The forces for 0 → 1 and 1 → 0 transitions are in a good agreement with MD simulations results.
(iv) DFT is a non-standard method for studying mechanical properties of membranes. Similarly to other methods, DFT has shown that the non-symmetric switching pathway is the preferable one. This method also provides a consistent estimation for the up-to-down threshold switching force.

Implications for memcapacitor design
From the application point of view, the membrane memcapacitor should have a strong 'memory content', namely, the device characteristics in its two logic states should be sufficiently different so as to provide a significant influence on other elements of electronic circuits 52 . For this purpose, it is desirable that both the distance between the fixed edge of membrane and second electrode (h in Fig. 1) and the maximal deflection of membrane (z s (0)) are chosen (much) smaller than the membrane length. Moreover, the gap between electrodes should be larger than the maximal membrane deflection, h > z s (0) (see Eq. (S2) for the definition of z i (s)). For given capacitances C 0 and C 1 of the states 0 and 1 (logic 1 corresponds to the membrane buckled downwards, and logic 0 -to the membrane buckled upwards) such that C 1 /C 0 < 3.01, one can find that Next, we need an expression for the force applied to the membrane expressed through the input variable of memcapacitor such as the applied voltage V (see Eqs. (1) and (2)). Since the energy of the plate can be written as where ε 0 is the vacuum permittivity and A = d · w is the plate area, the force is given by The overall scheme of memcapacitor switching is schematically presented in Fig. 9. Let us consider first the up-to-down transition (path (a-c,f) in Fig. 9). As mentioned in SI Sec. "Phenomenological Elasticity Theory" , the smallest threshold force for this transition is associated with the switching through the non-symmetric state and corresponds to c * 1 = 0.6948. In this way, the buckled upwards membrane (structure (a) in Fig. 9) evolves into the buckled downwards state according to (b) and (c) in Fig. 9. After the voltage removal, the membrane remains in the buckled downwards state ((e) in Fig. 9). In order to estimate the minimal (threshold) voltage required for this transition, Eq. (22) is rewritten accounting for Eq. (31): Consequently, where z cm (c * 1 ) = 0.3203 L(L − d). Eq. (33) provides an easy way to find the threshold voltage necessary for the up-to-down transition of membrane.
In the down-to-up transition, the initial buckled downwards membrane ((e) in Fig. 9) is first stretched by the applied force ((c) in Fig. 9). When the force is removed, the membrane state becomes (d). If the potential energy of (d) is sufficient to overcome the potential barrier height symmetrically (f), then the membrane may end up in the buckled upwards state (a). These general processes are partially described by Eqs. (24), (25). Because of the membrane elongation in (c,d), a refined value of z cm should be used in relevant calculations.
The threshold voltage can be estimated based on the energy conservation (Eq. (24)) and force equilibrium (Eq. (25)) conditions accounting for z * cm from Eqs. (26) and (27). One can find that and An additional aspect of membrane dynamics -an estimation of thermally-induced switching time -is considered in SI Sec. "Stability of buckled membrane".

Concluding Remarks
We have investigated the snap-through transition of a buckled graphene membrane using a variety of computational and theoretical tools. The main results of this paper are the expressions for the threshold switching forces (Eq. (12,22) for the up-to-down transition and Eq. (14,28) for the down-to-up one) and corresponding voltages (Eqs. (33) and (35)). Our analytical results are supported by the results of numerical simulations, MD and DFT calculations. We expect that our findings will find applications in the design, fabrication and analysis of membrane-based memory devices.
In graphene, the interactions between the carbon atoms were described using the standard 2-body spring bond, 3-body angular bond (including the Urey-Bradley term), 4-body torsion angle and Lennard-Jones potential energy terms 2 . The interaction constants were optimized in order to fit the experimentally observed properties of graphene. The known values were selected for the equilibrium angles. The Lennard-Jones coefficients were chosen to match the AB stacking distance and energy of graphite 3 . A global optimization was performed over the remaining parameters to match the in-plane stiffness (E 2D = 342 N/m), bending rigidity (D = 1.6 eV) and equilibrium bond length (a = 1.421Å) of graphene. We performed a series of test calculations that have demonstrated that the in-plane stiffness and bending rigidity of a selected graphene nanoribbon are in perfect agreement with the listed above values (see Fig. S10). A 1 fs time step was used and the system temperature is kept at room temperature with a Langevin damping parameter of γ = 0.2 ps −1 in the equations of motion. The van der Waals interactions were gradually cut off starting at 10Å from the atom until reaching zero interaction 12Å away. In each simulation, the buckled membrane was used as an initial condition. Typically, the system was equilibrated for 100 ps in the presence of the applied force and its subsequent dynamics (after the force removal) was typically simulated for 100 ps (these times were changed to 500 ps in the case of overdamped simulations). The final state was determined according to z-coordinate of an atom in the cental part of the membrane.
We also performed a series of overdamped simulations of membrane dynamics. This type of simulations was made using zero temperature and strong damping that significantly suppress the kinetic energy of the membrane. As a result, we have obtained more regular final states of the membrane (see Fig. S11) and slightly different values of the threshold switching force, which, however, can be more easily interpreted as now the kinetic energy can be disregarded. Comparing Fig. 2(b) and Fig. S11 we note that the overdamped calculations provide somewhat larger/comparable estimations for the threshold switching force.

Down-to-up transition
A study of the down-to-up transition based on the standard elasticity theory is presented in Fig. S12. Fig. S12(a) shows the final position of membrane depending on the applied force magnitude. In these calculations, the force was applied for a finite interval of time and removed at ω c t = 0.6. It was observed that the final state of membrane depends on few modeling parameters (the applied force, damping rate, etc.). The possibility of the down-to-up transition is clearly seen in Fig. S12(a).
The time-dependence of harmonics amplitudes is exemplified in Fig. S12(b) for a particular value of force. We note that in Fig. S12(b) the amplitudes q n are normalized to q 0 at t = 0. Fig. S12(b) shows that the harmonics amplitudes are modified by the applied force (e.g., it is clearly seen that |q 0 (t)/q 0 (0)| > 1 at ω c t < 0.6). Importantly, the amplitudes of higher harmonics are small and can be neglected. Time-dependence of the first few harmonics amplitudes q n (t).

Phenomenological Elasticity Theory
The Cartesian coordinates, which give the position of a membrane element defined by the internal coordinate s, are related to the angle θ i (see Eqs. (16) and (17)) as According to the above equations, x i (−1/2) = 0 and ζ i (−1/2) = 0. Another pair of boundary conditions is given by In the limit of small deflections, the boundary conditions, Eq. (S38), can be resolved for the coefficients A i and c i in Eqs. (16) and (17): A ns = ± c 1 240 Therefore, the value of a single parameter c 0 (or c 1 ) defines completely the symmetric (or non-symmetric) shape of membrane.
In the case of the symmetric shape of membrane (Eq. (16) Fig. S13. Importantly, the snap-through transition of membrane takes place as c 0 changes from 0.6954 to 0.2775. We note that Eq. (19) can be used to find F s needed to 'support' the membrane profile defined by a given value of c 0 . The force as a function of z cm is presented in Fig. S14 for several values of d/L.

Dynamic regime
Here we assume that the membrane transition is induced by an abruptly applied force (similar to the case of our MD simulations). Under the assumption of energy conservation, the threshold switching force can be found from the condition that the initial energy equals the potential energy barrier height (that depends on the applied force). Technically, the calculation can be performed as follows. The value of parameter c i , c * i (F), corresponding to the maximum of U = U b +U ext (in the presence of F, see Eq. (18)) can be obtained from dU/dc i = 0. Next, c * i (F) is substituted into where c 0 i corresponds to the potential energy minimum at F = 0 (note that U b (c 0 i ) is the energy corresponding to k 1 given below Eq. (S43)). Finally, the threshold switching force is found from Eq. (S44).
If one assumes that, at t = 0, the membrane has the non-symmetric profile, then one can get that the dynamic threshold switching force is the same as F ↓ ns given by Eq. (21). This result, however, should only be used as an upper estimate. Indeed, initially, the membrane has the symmetric profile and determining the exact point of switching from the symmetric to non-symmetric shape is a challenge. Unfortunately, the above mentioned equations can not be solved analytically for the transition through the symmetric profile. Figure S15 shows the numerically found potential energy as a function of position of the center of mass for several representative values of the applied force. This plot demonstrates that the potential barrier height becomes lower than the initial energy when the applied force becomes larger than the threshold switching force F sw .
Various threshold switching forces for the up-to-down transitions are summarized in Fig. 8 together with the results of our MD and DFT calculations. This plot shows that F ↓ provides the best phenomenological approach estimation for our MD results. The deviation of MD points from F ↓ curve at smaller d/L can be related to the approximation of small deflections used in our analytical model. Additionally, MD simulations involve energy relaxation channels not included into the simplified model given by Eq. (S44). Overall, however, we are quite satisfied with the results of our phenomenological model.

Computational details
Our DFT calculations were carried out using the Jaguar quantum-chemistry 5 package within the electron density functional approach with the use of a Becke three-parameter hybrid functional 6 , a Lee-Yung-Parr correlation functional 7 (the B3LYP method), and a 6-31G basis set of atomic orbitals. Optimization of structures was performed by an analytical method up to a gradient of the displacement of atomic positions of 10 −4 au. The minimum on the potential energy surface was identified by the absence of imaginary values in the matrix of second derivatives. As a test calculation, we have identified that the stretching of an armchair-edged graphene nanoribbon (L = 20.9Å, w = 14.7Å) yields the in-plane stiffness E DFT 2D = 338 N/m, which is in good agreement with the experimental value E 2D = 340 N/m. 8 We investigated a zigzag graphene nanoribbon of the length of L = 24.6Å (10 rings) and of the width of w = 11.4Å. The dangling bonds at nanoribbon boundary were terminated by hydrogen atoms. The hinged boundary conditions were realized by fixing the armchair edges at all calculation steps. The nanoribbon D1 with clamped boundary conditions was constructed from the flattest configuration C1 (with hinged boundary conditions) by an elongation of each armchair side by two rings kept fixed during calculations.

Up-to-down transition
First of all, we note that the optimization of unloaded configurations C1, C2, and C3 gives the bending rigidity D DFT = 1.75 eV, which is in good agreement with the experimental value D = 1.6 eV 9 .
As reported in Ref. [4], while, the non-symmetric path of switching is preferable due to the lower energy barrier, at small deflections from the unloaded state, the nanoribbon tends to be symmetric. A detailed investigation of the up-to-down transition of D1 shows that at small deformations, the difference between the energies of the symmetric and assymetric shapes is small (for instance, the initial displacement of membrane by 0.5Å requires about 0.02 eV and 0.05 eV for the symmetric and non-symmetric paths, respectively, see Fig. S16(a)). However, the energy difference between the paths increases significantly with deformation. When the positions of central atoms are close to the level of fixed edges, the deformation energy has a maximum, which was used to estimate the up-to-down threshold switching force 4 . We note that the symmetric path in Fig. S16(a) was generated by an increase of the z-coordinate of central atoms starting at z ∼ 0. The deformation energy of buckled graphene membrane with hinged boundary conditions is plotted in Fig. S16(b). In the case of hinged boundary conditions, only the non-symmetric path has been identified in the framework of our DFT approach. It was found that the transition energy barrier is smaller for larger values of d/L.
Overall, the up-to-down threshold switching force estimated from DFT calculations agrees with both MD simulation and elasticity theory results (see Fig. 8). This plot presents the threshold switching force in the units of Dw/L 2 for C1 and C2 configurations. We note that the DFT value for the threshold force for C3 is about 235.3Dw/L 2 . As d/L of C3 is significantly less than 1, the DFT threshold switching force for C3 is not the right quantity to compare with results of the linear elasticity theory.

Stability of buckled membrane
The information storage time in the state of buckled membrane is limited by the finite temperature which leads to thermallyactivated switching between two stable membrane states at long times. In the classical approach 10 , the rate constant k of 20/22 Figure S16. DFT deformation energy versus the position of central atoms in the process of the up-to-down switching. Based on DFT approach, we were able to identify (a) both symmetric and non-symmetric paths for the case of clamped boundary conditions, and (b) only the non-symmetric path for the case of hinged boundary conditions. thermally-activated switching can be approximately described by Arrhenius equation where A is the pre-exponential factor (according to some more advanced theories 11,12 , it describes the energy distribution in the system) and E b is the potential barrier height. In the case of buckled membrane, E b corresponds to the energy of the intermediate asymmetric state of membrane during its switching (the saddle point state) with respect to the energy of its stable states.
Next we make an estimation for the thermally-activated switching time of d/L = 0.94 membrane of type B (a detailed study of the membrane stability is beyond the scope of this paper). For this purpose we first performed a series of MD simulations of the dynamics of d/L = 0.98 membrane without applied force. These simulations have indicated that d/L = 0.98 membrane stays in its initial state for at least 10 ns without switching (this is, of course, an underestimate of its storage time). Then, assuming that the change of the pre-exponential factor in Eq. (S45) with d/L compared to the corresponding change in the exponential factor is not significant, and using E 0.98 = 0.60 eV, E 0.94 = 1.8 eV, and T = 300 K one finds k −1 > 40, 000 years. Clearly, such stability is sufficient for practical applications.