Relative Stability of Bernal and Rhombohedral Stackings in Trilayer Graphene under Distortions

Stackings in graphene have a pivotal role in properties to be discussed in the future, as seen in the recently found superconductivity of twisted bilayer graphene. Beyond bilayer graphene, the stacking order of multilayer graphene can be rhombohedral, which shows flat bands near the Fermi level that are associated with interesting phenomena, such as tunable conducting surface states expected to exhibit spontaneous quantum Hall effect, surface superconductivity, and even topological order. However, the difficulty in exploring rhombohedral graphenes is that in experiments, the alternative, hexagonal stacking is the most commonly found geometry and has been considered the most stable configuration for many years. Here we reexamine this stability issue in line with current ongoing studies in various laboratories. We conducted a detailed investigation of the relative stability of trilayer graphene stackings and showed how delicate this subject is. These few-layer graphenes appear to have two basic stackings with similar energies. The rhombohedral and Bernal stackings are selected using not only compressions but anisotropic in-plane distortions. Furthermore, switching between stable stackings is more clearly induced by deformations such as shear and breaking of the symmetries between graphene sublattices, which can be accessed during selective synthesis approaches. We seek a guide on how to better control -- by preserving and changing -- the stackings in multilayer graphene samples.


Introduction
Graphene remains a key two-dimensional material to still successfully reveal intriguing properties when stacked into just a few layers. As an example, bilayer graphene in twisted form has been found to exhibit superconductivity properties [1]. Large graphene samples are obtained experimentally by mechanical exfoliation of graphite, so they are supposed to favor the Bernal stacking, which for graphite has been considered more stable than the rhombohedral stacking [2]. However, this stability is far from being fully settled.
Firstly, recent calculations showed that graphite is energetically more stable in the rhombohedral stacking, so it is suggested that crystal growing conditions under high temperatures and pressures in a geological scale can stand behind graphite in the Bernal stacking [3]. Secondly, experimental and synthesizing techniques in 2D materials have been so improved over the years that graphene is being investigated nowadays with a defined number of layers [4,5]. Research works that focused on the same exfoliated sample of trilayer graphene (TLG) simultaneously found regions with the Bernal and rhombohedral stacking [6,7,8,9,10]. Rhombohedral graphenes are today being synthesized with a revival in a research race to show unique properties mostly due to their flat bands [11]. The stability between the Bernal and rhombohedral stacking in few layer graphenes therefore stands as a crucial question to be discussed.
Techniques are currently being developed on how to obtain, preserve or change stacking in trilayer graphene (TLG). The phase change between rhombohedral (ABC) and Bernal (ABA) stacking can be induced by external driving forces. The ABC stacking was experimentally found when graphene layers were grown on SiC substrates, and it was stable up to a temperature of around 1,200 C • when it reverts to the ABA stacking [12,13]. The transition from ABA to ABC stacking can be caused by a scanning tunneling microscope (STM) tip on highly ordered pyrolytic graphite (HOPG) samples due to the small barriers involved (∼ 1.0 mev/atom), as recently reported [7,14]. The transition between ABC and ABA stacking was also achieved using triazine decoration, which causes a large energy difference to induce the stacking transition [15]. The stacking in few-layer graphenes was stabilized mechanically by strain, using electric fields and even doping [16,17,18]. In order to fabricate high-quality encapsulated ABC trilayer -hBN devices used for quantum transport investigations, it has been shown that often the ABC trilayer converts to ABA during van-der-Waals transfer, even though ABA regions had been carefully removed prior to stamping [19]. The rhombohedral to Bernal transition has been theoretically studied by modeling the continuous sliding deformation of a side layer and by parameterizing results for two layers, calculations that were performed within density functional theory (DFT) using local density approximation (LDA) [20,21]. Although the van-der-Waals energies per atom are small, the total bonding energies are sizeable when the sample areas are considered, total energies that are large but needed to keep a certain phase locked in a metastable stacking after being synthesized. These findings suggest that not only the relative stability between these TLG stackings should be studied, but that their dependence on the lattice deformations deserves further study.
In this work, we address the stability between rhombohedral and Bernal stacking in trilayer graphene by calculations performed within density functional theory (DFT) using explicitly van der Waals functionals. We study different deformations that induced the energetic exchange between these stackings by looking at the relative stability. We first show that as the TLG ground state, the rhombohedral stacking is slightly more stable than the Bernal one. Then we perform several deformations in the quest to induce the transition from rhombohedral stacking to the Bernal one. In our simulations we find that the Bernal stacking is energetically more stable than the rhombohedral one under some distortions based on the planar expansion of the lattice parameters and their perpendicular compression. Furthermore, we show that under distortions such as shear and the atomic shifts in graphene sublattices, the structural anisotropy helps to stabilize the ABA stackings in certain directions. These mechanisms point out that the stacking transition is achieved by breaking the sublattices symmetries of layer-layer interactions involved in each stacking. These findings indicate that different practical realizations are possible, such as depositing samples in sandwiches, substrates, and molecule decoration, where stacking transitions may take place. Figure 1(a) shows the relaxed unit cells of trilayer graphene with the Bernal stacking on the left and the rhombohedral stacking on the right. Each stacking has a relaxed intralayer atomic distance of a c−c = 1.43Å, and the interlayer distance of d = 3.55Å in agreement with experiments [22]. The relative stability between Bernal and rhombohedral stackings is next analyzed by looking at the difference between their total energies. After geometry relaxations, the energy differences per atom are calculated as ∆E = (E B − E r )/N , where E B and E r are the total energies of Bernal and rhombohedral stacking respectively, and N = 6 is being the number of atoms in the unit cells. For pristine TLG our former results show a total energy difference of 0.079 meV/atom between both stackings. We have performed a detailed analysis of converging parameters to the accuracy required here, as shown in supplementary data. It seems theoretically established that the Bernal stacking is more stable than rhombohedral stacking in the literature, but this statement is based on inconsistent results [3]. In fact, we find that in TLG the rhombohedral stacking is slightly more favorable than the Bernal one. The energy differences although small are in agreement with enthalpy experiments between hexagonal and rhombohedral graphene [18,23]. Kinetic conditions such as temperatures and pressures seem to be the key in obtaining the Bernal stacking. All atoms have to shift the entire layer in order to pass from ABA to ABC stacking, and the other way around. We propose that these processes have to be considered in order to compare with the room temperature energy taking into account the total contact area. This means that there is a minimum area in the flakes as to obtain a stacking stable at room temperature. Depending on the value of stacking energy times area, there is a blocking temperature above which the separate layers can slide. In other words there is also a minimum area at a given temperature so that the flake stacking remain stable to be observed in experiments. An implication is that stackings can be different depending on the synthesis methods using exfoliation and CVD growth. We compare the electronic band structures of the Bernal and rhombohedral stackings along the ΓKM path near the Fermi energy, as shown in Figure  1(a). At the K point the Bernal stacking bands are displaying superposed a massive and massless electronic behavior, and the rhombohedral electronic structure shows massive electrons only; these two types of dispersion bands for each stacking are in agreement with Ref. [24]. The Bernal stacking can be seen as a superposition of single layer graphene bands and Bernal bilayer graphene bands; however, the graphenelike bands are slightly shifted up indicating the source of some intrinsic electric field between the layers to be added to external fields [25]. For the rhombohedral stacking there is a parabolic behavior with a flat region around the K point, source of many interesting properties, such as tunable conducting surface states [11] expected to exhibit spontaneous quantum Hall effect [26], surface superconductivity [27], and even topological order [28]. The massless band in the Bernal stacking now has acquired mass in the rhombohedral geometry. The conduction and valence bands are degenerated and slightly displaced from the K point in the KM direction. At lower energies, the rhombohedral stacking presents split bands what justifies its stability found in the calculations. A simple image can allow us to explain this stability. The Bernal stacking has directly connected three atoms between layers so that there are states at zero energy. However, the rhombohedral case has two pairs of atoms between layers so the corresponding electronic states are split from the zero energy, a fact that could be somehow linked to its higher stability.

In-plane deformations
We now study different in-plane deformations of the lattice vectors. Each deformation is characterized by a stretch factor λ that operates in the lattice vectors on different ways. First, we consider a homogeneous (or isotropic) deformation of lattice vectors described by a 1,2 = λ(a xî +a yĵ ). In this sense, having λ=1 (or strain δ = 0%) shows the previous relaxed lattice vectors, while a value of λ = 1.05 indicates an expansion strain of 5% (λ = 0.95 is a compression of 5%). This causes a proportional deformation of the lattice vectors keeping the interlayer distance constant. Figure 1(b) shows the total energy differences per atom between Bernal and rhombohedral stacking for the in-plane deformations explained above. Positive values of the energy differences (colored in red) indicate that the rhombohedral stacking is more stable than the Bernal one. On the contrary, negative values for the energy differences (colored in blue) refer to the case of the Bernal stacking being more stable. We show that although under homogeneous deformation the Bernal stacking becomes more favorable on expansion, the energy difference decreases for positive strains till it changes sign for ∼ 1%. At this point the Bernal stacking becomes more favourable. Next, we consider anisotropic deformations with uniaxial strain without and with constant area, described by the lattice vectors a 1,2 = λa xî + a yĵ ( a 1,2 = a xx + λa yŷ ) and a 1,2 = λa xî + a y /(λ)ĵ alonĝ x (ŷ) direction, respectively. The energy difference induced by the anisotropic deformation considering a constant area is slightly decreasing and remains almost constant in the whole range of considered positive or negative strains in the order of 5 %, as shown in Figure 1 (c). Furthermore, this figure includes the energy differences with the anisotropic uniaxial deformations, and we find that the Bernal stacking becomes more stable during stretching. The values of energy in the strain range ±5% show smaller response compared with those of homogeneous deformation, 2 versus 5 tenths of meV. The rhombohedral stacking with anisotropic strains is more stable for strain values till δ ∼ 2.2%, to be compared with the Bernal stacking with isotropic strains that becomes more stable under small strains δ ∼ 1.1% . Figure 2 summarizes the trends found between the Bernal and rhombohedral energy differences versus the in-plane deformations in TLG. There is an asymmetry in strains that induces that the most stable rhombohedral stacking switches into the Bernal stacking. The strain values when this transition is taking place are within the experimental range of reversible strains attained for samples in substrates, already found for a single layer of graphene [29]. In fact, strain values of beyond 1% could be artificially applied to graphene heterostructures [30] and strain up to 0.3% has been seen in graphene encapsulated in hBN [31]. Thus detailed care must then be taken in experiments when depositing exfoliated graphenes on substrates, and when using patterned contacts [10,32], so that strain is induced by top contacts while edge contacts can help in reducing strain [33]. Summary of energy differences with respect to in-plane deformations δx and δy. Note the anisotropy of the rhombohedral-Bernal stability versus in-plane deformations.

Out of plane deformations
We further study the relative stability order induced by the strain ǫ in the interlayer distance by stretching and compression between layers, as shown in Figure. 3. Although the energies of both stackings are almost degenerated at zero strain, ∆E under stretching shows that the rhombohedral continues being more stable even for an interlayer distance up to d = 4.26Å (ǫ ≤ 20%). For a compression at d = 3.11Å (ǫ = −2.5%), the ∆E values show that the Bernal stacking becomes more favorable, in agreement with another current calculation using empirically fitted vdW Grimme functional [34]. Small compressions induced by processes during transfer and protective layers will be enough to end into the ABA stacking. While we are not aware of direct measurements of compressions during h-BN encapsulation or subsequent cleaning by squeezing out residues [35], one can anticipate that compression of a trilayer can take place during these processing steps.
The next question is to determine how under a perpendicular pressure also the atomic in-plane layers are modified. It is noteworthy that under constant volume a perpendicular pressure implies that the layers are expanding on the plane. For instance assuming a constant volume, a value of ǫ = −3.5% corresponds to a homogenous expansion of δ = 1.88%. The critical δ value reinforces that the transition to Bernal stacking comes coupled to the in-plane deformations, as shown in Figure 1(b) and (c) above. These results indicate that the TLG stackings have to be analyzed in the sample once they are covered and isolated in the operating device whenever is possible to be certain that the rhombohedral stacking has not changed into Bernal.

Shear strain deformations
To assess the role of symmetry breaking we analyze the shear strain deformations. We consider the shear deformations along the perpendicular layer direction using the azimuthal φ and polar θ angles, as shown in the inset of Figure 4 and the movie included in supplementary data. The question is what θ values must be chosen before entering into different stackings from the rhombohedral and Bernal ones, i.e. how close the structures with shear can be seen as deformations of the original ABA and ABC stackings. The shear structures remain close to the pristine ABC and ABA stackings for the values with θ < 15 • , as shown in supplementary data, where shears beyond this critical angle show layers that can hardly be seen anymore as small deformations of the Bernal and rhombohedral stackings. Figure 4 shows energy differences between Bernal and rhombohedral stackings for the shears with θ = 10 • , 5 • changing φ angles. The shear for θ = 10 • also provokes that the interlayer distance d decreases from the value in pristine lattices so that ǫ = −1.5%; however, this ǫ value is still well below the critical distance to get the ABA stacking with perpendicular pressure. So that when obtaining Bernal stackings is not just correlated with the distance but better with the shear angle. Note that for θ = 0 • the shear is displacing the layers in the in-plane with respect to the perfect stacking as cones, shown in the upper panel, so that the overlap of the carbon p z orbitals between layers change related not only to distances but also to angles. In fact the ǫ value for θ = 5 • (ǫ ∼ 0.3%) just decreases slightly the interlayer distance, being larger than the critical distance to obtain ABA stacking under pressure. On the one hand, the energies of Bernal stacking exhibit a 60 • periodicity in φ, while increasing θ makes the energy profile just wider. That means that these deformations destabilize Bernal stacking minimally up to φ = 30 • . On the other hand, the energies for rhombohedral stacking does not show the same periodicity along φ angle, because their periodicity (120 • in φ) is twice than one found for Bernal stacking. Thus we find a region where the rhombohedral stacking is more destabilized than the Bernal one.
In Figure 4 we show the energy difference per atom between the two stackings under shear deformations. The curve for θ = 5 • shows the ∆E in the range of ∼ ±0.1meV /atom decreasing as φ values increase. The case of θ = 10 • where the ∆E are in the range of ∼ ±6 tenths of meV/atom much above the hundreths of meV/atom found for previous studied deformations; in fact they become an order of magnitude larger. For values of φ less than 30 • the rhombohedral stacking still remains as more favorable, and for the values of φ > 30 • , we find that the Bernal stacking is more stable than the rhombohedral one. Taking into account the symmetry, we find that the Bernal stacking becomes favorable for three fold directions of φ = 60 • as shown in the inset included in Figure 4. Interestingly, the ABA stackings for θ = 5 • are observed even for very small interlayer distance, values with very small ǫ. Even for a test case with ǫ = 0, i.e with no change in interlayer distance, the energies remain in the order of tenths of meV. In general, considering that those ∆E values are larger than previous considered deformation cases, the shear deformation of the stackings promote the Bernal stacking in certain directions, so that the area of stacking regions can even be reduced by a factor ten to be stable using shear deformations. These findings are in line with previous experimental results in which shear is induced by contacts [10], by other encapsulating layered graphenes within BN ‡ [16], and by explicitly putting shear to the samples [36].

Perpendicular sublattice distortions
Up to now all the investigated deformations affected both graphene sublattices within the same layer in the same way. We have then shown that shear introduced mostly anisotropy in plane. For the sake of completeness we now asses the role of deformations out of plane by looking at the anisotropy induced by the sublattices sites. We study the displacement up and down for each of the sublattices in both stackings. We characterize how much the same perpendicular shift affects the relative stability of each stacking. Here, we apply a perpendicular displacement to A or B sublattice type nodes, shifts that are breaking the sublattice layer symmetry. We shift an atom of a particular sublattice in the direction perpendicular to the layer plane by less than the layer compression to get ABA stacking, and compare how much the total ‡ Along zigzag and armchair direction the ABC and ABA stackings were favored, respectively energy of each stacking is changed. First we note that in each stacking, atoms of both side layers have the same number of interlayer neighbors. The energy response for atom shifts in side layers is very small, it does not exceed 0.01 meV. Thus we focus on the middle layer analyzing the coordination numbers of its atoms. Figure 5 (a) shows the Bernal and rhombohedral stacking schemes indicating the nearest neighbours in next layers. The vertical lines are linking the middle layer atoms to the first neighbors belonging to the side layers and shaded prims with numbers mark the second interlayer neighbors. In both stackings the A and B type nodes have in total 2 first neighbors and 18 second neighbors belonging to the side layers. However, the number of neighbors for each sublattice calculated separately is different. We can expect different responses when breaking the symmetry by shifting A or B atoms. In the case of Bernal stacking the numbers of interlayer neighbors for a sublattice site are the same with top and bottom layer. Thus the up or down shifts will be symmetric. In the rhombohedral stacking the numbers of neighbors for each sublattice site are "antisymmetric" i.e. shifting the A atom up is equivalent to shift the B type atom down and vice versa. For the rhombohedral geometry the number of cases to be discussed is reduced to two: where subscripts refer to the layer number, and arrows indicate the perpendicular displacement. Figure 5(b) shows the comparison of energy response difference defined as E shif t − E 0 , where for each stacking in the two types of displacements E 0 is the ground state energy, and E shif t is the total energy when the atoms are shifted. We shift the atoms by 1% of interlayer distance, i.e. 0.0335Å. In all the considered cases the total energy in each stacking increases by ≈ 0.5 meV. The energy response due to each sublattice shift remains nearly the same for the ABA stacking, and for the ABC stacking it changes more depending on the shift direction. For ABC stacking, moving the A atom up -equivalent to shift the B atom down -causes different energy response in the two directions. In the case of rhombohedral stacking, the energy increase is higher by 0.07 meV. As a result the Bernal stacking becomes more stable. The situation is just opposite when moving the A atom down (equivalent to shift the B atom up) -in this case the Bernal stacking is more destabilized, and the rhombohedral geometry remains favourable. The perpendicular sublattice shift studied here is a lattice deformation that causes an energy difference between stackings about 0.07 meV. This finding points out that breaking the sublattice symmetry is another way to stabilize the Bernal and rhombohedral stackings. These energy differences due to breaking of sublattice symmetry can be explained by having the TLG stackings interpreted again in terms of trimers and dimers between layers.
The total energy can be written in two parts as a sum of eigenvalues and a sum of interatomic potentials with other atoms. Looking at the total contribution of interatomic potentials the second neighbours atoms in the nearest layers are held responsible for these energy differences. This energy contribution adds to the difference in the sum of eigenvalues, already commented, in which the levels for the case of Bernal stacking split from the zero energy when the atoms are shifted up or down. These calculations raise further intringuing experiments regarding how the sublattice symmetry of stackings could be considerably broken by substrates and adatoms, or even better by having molecules and nanoparticles conforming these graphene heterostructures [37,38].

Conclusions
The aim of the present paper was to determine when the graphene stacking changes from rhombohedral to Bernal due to small lattice deformations. Under compression the rhombohedral stacking suffers a phase transformation to Bernal with strains ǫ < 2.5%. For even smaller compressions, shear deformations and shifts breaking graphene sublattice symmetries induce an anisotropy that stabilizes the Bernal stacking. These findings provide insights into the role of substrate-associated strains when graphene layers are integrated into devices, so that the stacking order and consequently their ultimate electronic properties are modified. This work would be of interest in relevant technological areas such as patterning contacts, and encapsulating graphene flakes between other materials. Our results are then asking for further experiments looking at the role of directional shears along and the adsorption of adatoms and molecules on graphene stackings. (project no. IT-1246-19), and the European Commission NRG-STORAGE project (project no. GA 870114). This research was conducted in the scope of the Transnational Common Laboratory (LTC) "Aquitaine-Euskadi Network in Green Concrete and Cement-based Materials."

Methods
We performed our study within density functional theory (DFT) using the ab-initio simulation package VASP [39,40,41]. This method based on plane waves is applied using a well converged kinetic energy cut-off of 700 Ry. In the energy range that we are interested, the mesh grid on the reciprocal space is key to get consistent results (see supplementary data), a 288×288×1 centered at the Γ point is used. For the dispersive interactions between layers, we perform our calculations using the Van der Waals functional vdW-DF2 [42]. The electronic convergence was performed with 10 −7 eV. We relax the cell shape and keep the volume constant, and we fix the ions in the xy plane while allowing them to move in the z direction. Tests performed with other functionals, as well as detailed convergence test in k-points, and electronic and ionic relaxations are included in supplementary data. Note that a fine k-mesh is required to be applied with the number of sampling k points in the x and y directions being a multiple of 3 to explicitly include the points K, K´and a region around them, regions in the reciprocal space where the low energy physics of multilayer graphene occurs.