Strain-induced stacking transition in bilayer graphene

Strain, both naturally occurring and deliberately engineered, can have a considerable effect on the structural and electronic properties of 2D and layered materials. Uniaxial or biaxial heterostrain modifies the stacking arrangement of bilayer graphene (BLG) which subsequently influences the electronic structure of the bilayer. Here, we use Density Functional Theory (DFT) calculations to investigate the interplay between an external applied heterostrain and the resulting stacking in BLG. We determine how a strain applied to one layer is transferred to a second, 'free' layer and at what critical strain the ground-state AB-stacking is disrupted. To overcome limitations introduced by periodic boundary conditions, we consider an approximate system consisting of an infinite graphene sheet and an armchair graphene nanoribbon (AGNR). We find that above a critical strain of ~1%, it is energetically favourable for the free layer to be unstrained, indicating a transition between uniform AB-stacking and non-uniform mixed stacking. This is in agreement with a simple model estimate based on the individual energy contributions of strain and stacking effects. Our findings suggest that small levels of strain provide a platform to reversibly engineer stacking order and Moir\'e features in bilayers, providing a viable alternative to twistronics to engineer topological and exotic physical phenomena in such systems.

E-mail: nuala.caffrey@ucd.ie Abstract. Strain, both naturally occurring and deliberately engineered, can have a considerable effect on the structural and electronic properties of 2D and layered materials. Uniaxial or biaxial heterostrain modifies the stacking arrangement of bilayer graphene (BLG) which subsequently influences the electronic structure of the bilayer. Here, we use Density Functional Theory (DFT) calculations to investigate the interplay between an external applied heterostrain and the resulting stacking in BLG. We determine how a strain applied to one layer is transferred to a second, 'free' layer and at what critical strain the groundstate AB-stacking is disrupted. To overcome limitations introduced by periodic boundary conditions, we consider an approximate system consisting of an infinite graphene sheet and an armchair graphene nanoribbon (AGNR). We find that above a critical strain of ∼ 1%, it is energetically favourable for the free layer to be unstrained, indicating a transition between uniform AB-stacking and non-uniform mixed stacking. This is in agreement with a simple model estimate based on the individual energy contributions of strain and stacking effects. Our findings suggest that small levels of strain provide a platform to reversibly engineer stacking order and Moiré features in bilayers, providing a viable alternative to twistronics to engineer topological and exotic physical phenomena in such systems.
2D materials can further be combined to create heterostructures which can have different properties to their component layers due to interlayer interactions [6,7]. For example, the 2D elastic moduli of bilayer heterostructures, such as graphene/MoS 2 and MoS 2 /WS 2 , are smaller than the sum of the moduli of the individual layers [8].
Interlayer interactions can be tuned, for example, by changing the relative stacking of the layers. This allows a wide range of different behaviours to be observed, even in structures with multiple layers of the same material.
Two monolayer graphene layers (MLGs) can be stacked to form bilayer graphene (BLG) [9,10]. In the groundstate AB-stacking, half of the carbon atoms in each layer are directly above the centre of a hexagon on the other layer and the other half are directly on top of another carbon atom [11,12] (figure 1(a)). AAstacking has also been observed in which every carbon atom of one layer is directly above a carbon atom on the second layer ( figure 1(b)). AB-stacked BLG has parabolic electronic bands, whereas the bands of AA-stacked BLG remain linear, as in MLG [9,13]. As a result, AB-and AA-stacked BLG behave very differently under the application of interlayer bias, which opens a band gap proportional to the bias for AB-stacking [14], while the AA-system remains semimetallic [15]. Between the AB-and AA-stacking limits lie a range of different stacking possibilities. The electronic and topological properties of these systems varies as a function of relative shift between the layers [16,17]. For example, the transport properties of BLG depend sensitively on its stacking, with a large change in the transmission predicted when one layer is shifted relative to the other [16]. Thus, manipulating the stacking of a bilayer is a powerful tool to tune its electronic and transport behaviour. Twisted BLG (TBLG) harnesses this to great effect [18][19][20][21].
Both uniform and non-uniform strain can tune the electronic, transport and optical properties of both MLGs and BLGs [22][23][24][25][26][27][28][29][30]. Strain can arise naturally when graphene is placed on a substrate due to lattice mismatch between the two systems [31,32]. Strain can also be intentionally created and controlled in graphene by using a flexible substrate [33,34]. If BLG is placed on a flexible substrate, strain will be applied to the layer which is in direct contact with the substrate [35]. It is reasonable to expect that if the applied strain is sufficiently small, it will be entirely transferred to the second graphene layer, i.e., both layers will experience the same strain. However, for larger applied strains the second layer can exhibit a different strain profile [35][36][37]. A heterostrain modifies the stacking between layers in a manner similar to a relative twist angle, creating a Moiré superlattice [38][39][40]. Theoretical works predict that heterostrain can be used to open and tune an electronic energy gap in BLG [38], or to induce a transition from a direct to an indirect band gap in the presence of a bias field [41].
In this work, we investigate the interplay between uniaxial heterostrain and stacking effects in BLG. The system consists of a bottom layer which is uniaxially strained along the zigzag (ZZ) direction and a 'free' top layer as shown in figure 1(c). For small amounts of strain applied to the bottom layer, the free layer is expected to strain by the same amount, as the energetic cost of straining is small compared to that of breaking uniform AB-stacking. However, at a certain critical value of applied strain, the energetic cost of maintaining strain in the free layer exceeds the cost of breaking AB-stacking. At this point, strain is released in the free layer as the energy benefit in doing so is greater than the energy penalty to be paid for a less-than-ideal stacking configuration. At this critical strain, a transition between a uniform ABstacking order and a non-uniform stacking will occur, together with the formation of Moiré superlattices, as in TBLG [20,42]. Near this critical point, the change in stacking order can have significant implications on electronic properties [16,17,38].

Computational Details
Density functional theory calculations were performed using vasp-5.4.1 [43][44][45][46][47][48]. The Perdew-Burke-Ernzerhof (PBE) [49] parametrization of the generalized gradient approximation (GGA) was employed. Van der Waals (vdW) interactions are included using the D2 semi-empirical method of Grimme [50]. The plane wave basis set was converged using an 950 eV energy cutoff. A 13×21×1 mesh was used to determine the total energies of MLG and BLG. All structures were optimized until the residual forces were less than 0.01 eV/Å. A vacuum layer of at least 11.5 Å was included in the direction normal to MLG or BLG to ensure no spurious interactions between repeating slabs. The GGA calculated lattice constants of MLG and BLG were both found to be 2.47 Å, in good agreement with the experimental value of 2.46 Å. The interlayer distances in AB-and AA-stacked BLG were found to be 3.37 Å and 3.50 Å, respectively, in good agreement with previous studies [51][52][53].
To overcome restrictions caused by the use of periodic boundary conditions in the graphene plane, the heterostrained BLG system was approximated using a hydrogen-terminated armchair graphene nanoribbon (AGNR) adsorbed on top of a MLG. The dangling C bonds at the AGNR edges were passivated with H atoms. Heterostrain was then introduced by straining the MLG along the ZZ direction. A 21×5×1 k-point mesh was used to converge the total energy of the heterostrained structures. A distance of at least 12.5 Å was maintained between periodic replicas of adsorbed AGNRs to ensure that they do not interact. The optimal AGNR width was determined by comparing the stacking-dependent AGNR binding energy (E B ) to the stacking-dependent binding energy of BLG. In this case, the in-plane positions of the two central atoms of the AGNR were held fixed at the chosen stacking and all other AGNR atoms were allowed to relax both inplane and out-of-plane. To determine the lowest energy stacking configuration for different heterostrains, the carbon-hydrogen bonds and the interlayer distance between the MLG and the AGNR were relaxed while the in-plane positions of all other atoms were held fixed.

Results
To determine the critical strain for which a transition occurs from a uniformly strained BLG with ABstacking to a heterostrained system with disrupted stacking, we compare the energetics of two limiting cases: when the free layer either adopts the same strain as that applied to the bottom layer, or it remains completely unstrained. Possible intermediate scenarios, where the free layer adopts a non-zero strain different to that in the bottom layer or displays a non-uniform strain distribution, are not considered in this work due to their high computational cost. Similarly, in order to maintain periodicity for the DFT calculations, we neglect the role of contraction in the direction perpendicular to the applied strain. This is equivalent to setting the Poisson ratio ν = 0. We discuss how these approximations can be relaxed, and the expected consequences, in Section 4.

Simple Model
There are two principal energy costs, due to strain (∆E strain ) and stacking (∆E stack ), which determine the behaviour of the free layer when uniaxial strain is applied to the bottom layer. To get a rough estimate of where the transition between a strained and unstrained free layer occurs, we can compare the expected energy costs of straining the free layer in isolation, and of breaking AB-stacking in an unstrained BLG system.
We first consider ∆E strain = E −E unstrained , the energy cost associated with straining the free layer away from its relaxed structure to match the strain applied to the bottom layer. The energy cost of straining a graphene layer increases with the amount of strain considered, as shown in figure 2(a) for MLG with strain applied along the ZZ direction. ∆E strain displays almost identical behaviour in AA-and AB-stacked BLG, and for strains along the AC direction in all three systems, with a maximum variation of only 2.8 meV per atom. Since ∆E strain is not significantly affected by the strain direction, or the nature of the stacking in BLG systems, the curve in figure 2(a) should also be an excellent approximation to the energetic cost of straining a single layer in BLG in the absence of stacking effects.
If the amount of strain is different in the two layers, then the system is no longer able to maintain energetically-favourable AB-stacking, and instead must display a modulation of the stacking order with an associated energy cost ∆E stack = E − E AB . Although the modulation wavelength depends on the strain mismatch, the energy cost per atom is roughly constant, as a similar range of stackings will occur for any mismatch. Therefore, ∆E stack should not depend sensitively on the strain applied to the bottom layer, and we approximate it by considering different stacking configurations in unstrained BLG. Figure 2(b) shows the energetic cost of rigidly shifting one graphene layer over the other along ZZ or AC directions, starting from either an initial AB-or AA-stacking. As the layers are shifted, the in-plane positions of atoms are held fixed to maintain the desired stacking configuration, but the interlayer distance is allowed to relax. Any shift away from AB-stacking results in a positive ∆E stack , confirming that this is the preferred configuration. ∆E stack is maximum for AA-stacking with a value of 5.84 meV/atom, in excellent agreement with Ref. [16]. The stacking modulation arising from heterostrain contains a combination of various stackings, and the corresponding ∆E stack can be approximated as an appropriately weighted average of the values appearing in figure 2(b).
Comparing figure 2(a) and 2(b) allows us to understand the interplay between strain and stacking in heterostrained BLG. For small strains applied to the bottom layer, ∆E strain ∆E stack , and it is energetically favourable for the free layer adopt the same strain. However, as the strain applied to the bottom layer is increased, the cost of uniformly straining the free layer eventually balances the cost of breaking AB-stacking. A transition occurs above this critical strain, releasing the strain in the free layer and introducing a modulation of the stacking order. The maximum possible value of critical strain is restricted by the finite range of ∆E stack : the maximum possible cost of breaking AB-stacking (i.e., ∆E AA stack ∼ 5.84 meV/atom), corresponds to a strain of 1.4%. In reality, the mix of different stackings that occur in a hetero-strained BLG will give 0 < ∆E stack < ∆E AA stack . For an even distribution of stackings between AB and AA, we can estimate ∆E stack ∼ 1 2 ∆E AA stack , corresponding to a critical strain of 0.97%. Uneven stacking distributions can occur if the considered strain excludes certain stackings, or if a non-uniform strain distribution is allowed in the free layer. However, even accounting for a significant reduction in ∆E stack due to these effects does not dramatically change the expected critical strain. For example, assuming ∆E stack ∼ 1 4 ∆E AA stack still gives a critical strain of 0.66%. The results of this simple model strongly suggest that the critical strain is near 1%, and that uniform AB-stacking will be broken when larger strains are applied to the bottom layer.

AGNR on strained MLG
DFT calculations of heterostrained bilayers will now be used to test the prediction of this simple model that the critical strain occurs near 1%. However, periodic boundary conditions enforce a commensurability condition when dealing with two infinite graphene sheets.
Although neglecting contraction in the transverse direction (ν = 0) simplifies matters considerably, only certain values of strain would be achievable [38]. Furthermore, very large supercells would be required to investigate the relevant strain range: 1% strain requires 101 cells of the free layer and 100 cells of the strained layer. Getting sufficient resolution to determine the critical strain quickly becomes computationally prohibitive.
To overcome this constraint, we instead investigate the interplay of strain and stacking in a system where different regions of the free layer are modelled by finitewidth AGNRs. This is shown schematically on the left side of figure 3, for an exaggerated strain of 15% along the ZZ (vertical) direction applied to the bottom layer. A one-dimensional Moiré pattern is evident in the full BLG with regions of AB-stacking (lower dashed box) separated by other stacking types. Due to the ZZ strain direction and ν = 0 , no AA-stacking occurs and the furthest stacking from AB is that in the upper dashed box, which we denote 'Shifted' and corresponds to rigidly shifting one layer of unstrained BLG by half a graphene lattice constant in the ZZ direction. We aim to determine the energetics of the complete system (left hand side) by modelling different portions of it by a finite-width ribbon adsorbed onto an infinite bottom layer (right hand side). The AGNR can be rigidly shifted over the continuous bottom layer to approximate the different stackings that occur in a heterostrained bilayer system. As the free layer is no longer continuous, but now consists of a periodic array of AGNRs, we can consider different strains in each layer using a constant-size supercell. We consider hydrogen-passivated AGNRs to circumvent features including unpassivated bonds or localised edge states. These may occur in certain AGNRs but are not expected in extended bilayer systems [54][55][56].
The choice of AGNR width is determined by that which best approximates the stacking-dependent binding energy of BLG at a reasonable computational cost. To meaningfully compare the binding energy of BLG and the n-AGNR/MLG systems we normalise E B by the number of carbon atoms in the top-layer, N C : (dashed line). Agreement between the BLG and AGNR/MLG systems improves in general with the increase of n, but non-uniformly due to different behaviour of ribbons with widths n = 3q, 3q + 1 and 3q + 2, where q = 1, 2, 3, . . . . Similar trends have been noted, for example, for the band gap of AGNRs [57]. As we are interested in how stacking changes the energetics, and not the absolute magnitude of the E B , the 6-AGNR is deemed sufficiently wide for our purposes.
To estimate the critical strain in heterostrained BLG using the 6-AGNR/MLG system, we consider the energy difference, ∆E X , between strained and unstrained AGNR layers, for different stackings (X = AB, shifted). We emphasise that it is the unstrained case which gives rise to broken stacking, due to the strain applied to the bottom layer, whereas the strained case restores AB-stacking by matching the strain in the both layers. The energy difference between strained and unstrained AGNR layers is then given by: where we note that the energy of the final state E AB, strain is the same in each case, as AB-stacking has been restored. E X,unstrain corresponds to the case that the MLG is strained and the AGNR is unstrained, with the stacking, X, set by fixing the positions of the central carbon atoms of the AGNR. The in-plane positions of the AGNR carbon atoms are held fixed as determined by the strain and stacking, while the hydrogen-carbon bonds and the interlayer distance are allowed to relax. The full heterostrained system is considered as an average of the two stacking extremes, AB and Shifted (more general cases will be discussed later): (3) Negative values of ∆E av imply that the free layer prefers to be strained so that the bilayer system remains AB-stacked. Positive values indicate that the free layer prefers to be unstrained, and the system adopts a non-uniform stacking profile. Figure 5 shows ∆E and ∆E av , as a function of strain applied along the ZZ direction. The curve for the ABstacked AGNR/MLG system shows that it prefers to be unstrained, i.e., that it is energetically favourable for the AGNR to break perfect AB-stacking, instead of maintaining a strain of between 0.8% -1.3%. This is not surprising, as the stacking mismatch is not too significant in this part of the modulated structure (c.f. figure 3). However, in the shifted region, the stacking deviates furthest from AB when the free layer is unstrained, and the associated ∆E curve shows that a strained, uniform AB-stacking is preferred. The overall preference of the system is a competition between these different regions. For a uniform distribution of stackings in the modulated structure, the averaged case shown by the grey curve in Figure 5 indicates that a transition occurs at a critical strain of ∼ 1.2%. This is in very good agreement with the estimate given by the simple model.

Discussion
Modulated stacking profiles in BLG systems lead to a wide range of new properties, particularly when combined with an interlayer bias.
Such a bias will open a band gap in regions with particular stackings, leading to a complex distribution of gapped and conducting regions which follow the underlying stacking pattern. This has been widely studied in twisted bilayers where, for example, networks of 1D topological channels have been found between gapped AB-and BA-stacked regions of the Moiré pattern [58][59][60]. The localisation of electronic states in different regions of the lattice, together with the resulting strong interaction effects, are connected with the formation of correlated insulating states and unconventional superconductivity in twisted systems [61].
Heterostrained, untwisted BLG, as considered here, could potentially host a similar range of phenomena. The schematic structure in figure 3 shows the formation of a 1D Moiré pattern with different stacking profiles, which would also create a spatially-varying band gap landscape in the presence of an interlayer bias. To maintain periodicity for our calculations, we neglected a Poisson contraction perpendicular to the applied strain. Including such a contraction would lead to a 2D modulation of the stacking pattern and a more complete analogy with twisted systems. We note that setting ν = 0 allows a wider range of stackings, when the free layer is unstrained, than for the simple 1D modulation considered in figures 3 and 5. In particular, more energetically unfavourable stackings, such as AA, are now possible, which may slightly increase the critical strain due to the increased energetic cost of breaking uniform stacking.
Finally, the estimate of the critical strain in this work is based on the free layer being either unstrained, or uniformly adopting the same strain as the bottom layer. We have not considered, for example, an intermediate value of strain in the free layer. We expect such cases to be less energetically favourable than the unstrained free layer. This is because such cases will have to pay both strain and stackingrelated energy costs, and the latter is expected to be largely strain independent, as discussed in Section 3.1. Therefore the minimum energy heterostrained system should be that which minimises ∆E strain , namely that with an unstrained free layer. Due to the periodicity constraints of our calculations, we also not been able to explicitly consider the role of non-uniform strain in the free layer. Such non-uniform strain profiles could arise due to lattice relaxations, as occurs in TBLG [62][63][64], where it serves to reduce the area of the AA-stacking region and maximize the area of the lower energy ABstacking region [65].

Conclusions
Our results strongly indicate the existence of a critical strain of ∼ 1% which, when applied to one layer of BLG, can be used to tune the stacking profile of the system. Below the critical strain, it is energetically favourable to transfer the strain to the second layer in order to maintain a uniform AB-stacking configuration. Above the critical value, the cost of maintaining the strain in the second layer is too high and the system prefers to release it and adopt a non-uniform stacking profile. This finding is supported by DFT calculations which consider the energetic costs of strain and stacking independently in infinite, periodic systems, and by further calculations which consider both contributions simultaneously when one layer is represented by finite AGNRs.