Snap-through transition of buckled graphene membranes for memcapacitor applications

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.

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".
In this paper, the following notations are used: q -the charge on capacitor (see Eq. (1)) V -the applied voltage (see Eqs (1) and (2)) C -the (memory) capacitance (see Eq. (1)) d -the distance between fixed sides of membrane (see Fig. 1) h -the distance between the bottom plate and the level of fixed sides (see Fig. 1) L -the membrane length w -the membrane width D = 16 eV -the bending rigidity of graphene E 2D = 340 N/m -the 2D Young's module ζ -the deflection of membrane (see Eq. (3)) ζ c -the maximum deflection of membrane (see Eq. (7)) θ i (s) -the angle that the membrane makes with the horizontal (see Eqs (16) and (17)), s -the internal coordinate that changes between −1/2 and 1/2 (see Eqs (16) and (17)) A i and c i -coefficients (see Eqs (16) and (17)) z cm -the center of mass position (see Eq. (19)) U b , U str , U ext -the bending, stretching and external potential contributions to the potential energy of membrane (see Eq. (18)) F ↓ -the up-to-down threshold switching force (see Eqs (12), (20), (21) and (22)) F ↑ -the down-to-up threshold switching force (see Eqs (14) and (28)) ε 0 -the vacuum permittivity

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 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.
Buckling and snapping-through within 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, 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 the potential is the double-well potential with minima at symmetric deflections, ζ α β = ± /2 c , 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, nx d sin( / ) π ∼ , 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 Expanding this expression to the second order, one obtains the interrelation between L, d, and the maximal deflection ζ c : c 8 On the other hand, the maximal deflection of buckled membrane can be inferred from the minima of the potential energy, Eq. (7): 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 n n n 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 ω μ = L D (1/ ) / c 2 and t c = 1/ω c ) 51 are analyzed for membrane B.
Consider the results presented in Fig. 5. Figure 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 Sec. "Discussion" below. 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 ζ c 2 as 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: 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 U b,s = 108Dw(L − d)/L 2 . 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 c 1/(2 5 ) 1 2 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 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 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 2 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

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,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 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.

Discussion
Four methods to describe the snap-through transitions. As the key finding, Fig. 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 Subsec. "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 c 2 . 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 1 0 1 0 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 ext cm cm where ε 0 is the vacuum permittivity and A = d · w is the plate area, the force is given by cm cm 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 = . 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. Consequently, . 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) and (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   . Schematic representation of transition pathways in the switching of membrane memcapacitor (red arrows represent the applied electrostatic force, bold blue arrows represent the movement of membrane). The buckled upwards membrane (structure (a)) evolves into the buckled downwards state (c) through a nonsymmetric transition state (b). Note that in the presence of force, (c) is also a transition state. After a slow voltage removal (dashed green arrow), the membrane remains in the buckled downwards state (e) (the up-to-down transition). In case of fast force removal, the membrane does not have enough time to relax its tension (d) and start moving into the upward state (a) through a symmetric transition state (f) (the down-to-up transition).