Aitor Garcia-Ruiz,Vladimir Enaldiev,Andrew McEllistrim,Vladimir I. Fal’ko
Mixed-Stacking Few-Layer Graphene as an Elemental Weak Ferroelectric Material
Ferroelectricity,graphene,rhombohedral graphite,twin boundary,twistronics,screening
Ferroelectricity ( Valasek J. Phys. Rev. 1921 , 17 , 475 ), a spontaneous formation of electric polarization, is a solid state phenomenon, usually, associated with ionic compounds or complex materials. Here we show that, atypically for elemental solids, few-layer graphenes can host an equilibrium out-of-plane electric polarization, switchable by sliding the constituent graphene sheets. The systems hosting such effect include mixed-stacking tetralayers and thicker (5–9 layers) rhombohedral graphitic films with a twin boundary in the middle of a flake. The predicted electric polarization would also appear in marginally (small-angle) twisted few-layer flakes, where lattice reconstruction would give rise to networks of mesoscale domains with alternating value and sign of out-of-plane polarization.
Mixed-Stacking Few-Layer Graphene as an Elemental Weak Ferroelectric Material

Ferroelectricity (ValasekJ.Phys. Rev.1921, 17, 475), a spontaneous formation of electric polarization, is a solid state phenomenon, usually, associated with ionic compounds or complex materials. Here we show that, atypically for elemental solids, few-layer graphenes can host an equilibrium out-of-plane electric polarization, switchable by sliding the constituent graphene sheets. The systems hosting such effect include mixed-stacking tetralayers and thicker (5–9 layers) rhombohedral graphitic films with a twin boundary in the middle of a flake. The predicted electric polarization would also appear in marginally (small-angle) twisted few-layer flakes, where lattice reconstruction would give rise to networks of mesoscale domains with alternating value and sign of out-of-plane polarization.

Ferroelectricity, a spontaneous electric polarization in the absence of an external electric field, is a phenonmenon observed and thoroughly investigated in a broad range of solids. Caused by charge transfer between constituent atoms in a unit cell, up to now, ferroelectric polarization was observed only in chemically complex compounds, such as Rochelle salts, Pb[ZrxTi1–x]O3, BaTiO3, etc. Here, we show that there is one exception from this common rule, namely, several structural allotropes of multilayer graphene. The above statement is based on a theoretical study of various multilayer graphene structures with different stacking orders which, as an elemental material, is nonpolar: its intralayer bonding is dominantly covalent, whereas the interlayer adhesion has a van der Waals nature. There are two commonly studied multilayer graphene systems: thin films of Bernal graphite and rhombohedral (ABC) graphite. Those two possess inversion symmetry, or z → −z, respectively, which prohibit spontaneous ferroelectric polarization. However, graphitic films with mixed Bernal and rhombohedral stackings lack such symmetries, removing constraints on the formation of an out-of-plane electric dipole. One example of such a structure is a thin film of rhombohedral graphite with a twin boundary inside it, which consists of two ABC parts with n and m layers (n ≠ m), rotated by 180° with respect to each other, held together by an ABA-trilayer (see Figure 1). In the following section, we study nABAm, with n ≠ m, and demonstrate that they exhibit a weak spontaneous out-of-plane electric polarization, Pz, at zero doping. We determine the roles of parameters in the full multiparameter Slonczewski–Weiss–McClure (SWMcC) Hamiltonian that are related to the asymmetries that permit a finite Pz. Importantly, we find that self-consistently implemented screening of internal electric field not only strongly reduces Pz magnitude but also changes its sign for some structures, which is a result of different interlayer charge redistributions, caused individually by each symmetry-breaking term in the bare SWMcC Hamiltonian. For tetralayers, the known values and signs of SWMcC parameters are such that the largest contributions they produced individually mostly cancel each other, making the result sensitive to the precise choice of parameters. At the same time, for marginally twisted tetralayers, where lattice relaxation leads to the formation of equilibrium stacking domains (where polarization is small), the domain walls host seeds of different stacking arrangements for which such cancellation does not occur leading to an order of magnitude larger local Pz values. Below, we use the full Slonczewski–Weiss–McClure (SWMcC) model of graphite films discussed in ref (15). This model accounts for couplings sketched in Figure 1, which include both the closest and next-neighbor hoppings, and it is implemented in the framework of a hybrid k·p-tight-binding model, using a Hamiltonian, , specified in the Supporting Information. The diagonalization of gives the dispersions, εβ(p), and wave functions, ψβ,pαL, of bands β in the multilayer, which we use to compute both on-layer electron densities, nL, and electric polarization:Here e is the electron charge, d ≈ 3.35 Å the interlayer distance, and “L” the layer index, and αL = AL, BL sublattice indices. When computing the electron densities, nL, we take account of the screening the potential they generate on each individual layer, aswhere εz ≈ 2.6 is an effective out-of-plane dielectric permittivity of graphene stacks determined by the polarizability of carbon atoms. Among all SWMcC couplings, next-neighbor hoppings γ2, γ5, and v4, together with the energy difference between dimer and non-dimer sites, Δ′, are most important for a nonzero value of the ferroelectric polarization, Pz, to form across the structure. This is because such couplings break a hidden electron–hole symmetry,a ( is a unitary matrix equal to direct product of n+3+m-rank identity and third Pauli matrices, respectively), characteristic of the “reduced” models of graphene limited to the closest neighbor couplings only. The on-layer potentials in eq 2 induced by the screening can also break all symmetries, which will be important for understanding the result of self-consistent analysis. For a detailed quantitative self-consistent analysis of spontaneous electric polarization Pz, we implement the following steps. First, we compute the electric polarization induced by one of the four symmetry-breaking terms in the SWMcC Hamiltonian for each of the considered structures, both with and without self-consistent implementation of screening. Then, we check that the cumulative effect of all the terms in the full SWMcC model can be approximated as a sum of the individual symmetry-breaking contributions asHere, are dimensionless factors for which values for the ABCB tetralayer are listed in Table 1. Accounting for the charge redistribution self-consistently is an important part of the presented calculations. Taking into account redistribution of charges produced by screening of electric fields attributed to polarization not only reduces the value of Pz by 1 order of magnitude, but also can change the sign of some of the individual symmetry-breaking contributions, when compared with the unscreened case. We investigate numerically this effect by comparing the values for in ABCB graphenes, with and without the self-consistent implementation of screening in Table 1. In the absence of screening (ϵi = 0), we obtain a value of Pz ≈ 0.5 e/μm, whereas implementing screening changes this value to Pz = −0.07 e/μm. This is because the two largest contributions to the ferroelectric polarization, coming from and , have opposite signs and nearly cancel each other for the choice of SWMcC based on ref (16), leaving γ2 to define the value of Pz, so that the overall result changes sign after the implementation of screening. We suggest that such a sensitivity of the computed values of Pz to the input SWMcC parameters may be used to further narrow down their choices by comparing the computed Pz with the experimentally measured polarization of tetralayers. We find that screening is equally important for the analysis of polarization in thicker films, though the above-mentioned cancellation does not occur for all thicknesses. In Figure 1, we show the computed values of Pz in various nABAm structures with m > n (the mirror-symmetric mABAn configurations have the same magnitude of Pz but opposite sign). While the above-described twinned rhombohedral graphitic films may appear naturally among flakes exfoliated from bulk material, one can try to build weakly ferroelectric multilayers by a twistronic assembly of thinner flakes. A small misalignment of crystalline axes of the assembled graphene flakes, unavoidable in the mechanical transfer, leads to a long period variation of stacking, known as moiré structure, which upon lattice relaxation (promoting energetically preferred local stacking order) results in a network of domains with AB and BA interlayer preferences at the twisted interface. Here, we consider structures with graphene monolayers small-angle twisted with respect to Bernal and rhombohedral trilayer graphene (3ABA+1 and 3ABC+1, respectively) and a twisted double-bilayer (2AB+2BA). At the twisted interface, the moiré superlattice undergoes lattice reconstruction promoting formation of triangular domain arrays with local AB and BA stackings, shown in Figure 2, which are separated by a network of domain walls (DWs). We model lattice reconstruction in the structures following the approach of refs (24 and 25) to find local deformations of lattices in top/bottom flakes, ut/b, that minimize (via domain formation) a sum of elastic and interlayer adhesion energies (see Supporting Information). Then, we find the ferroelectric polarization of small-angle twisted structures using local stacking approximation. That is, we calculate the polarization based on a Hamiltonian shown in Methods, describing the aligned tetralayer films with an in-plane offset r0 = θẑ × r + ut – ub between constituent parts. As a reference, we set r0 = 0 for ABBA-stacking in 2AB+2BA, ABAA-stacking in 3ABA+1, and ABCC-stacking in 3ABC+1 areas. Having calculated Pz with a self-consistent analysis of screening and on a sufficiently dense grid of interlayer offsets r0, we interpolate it in the form of a harmonic series,where Λ labels reciprocal lattice vectors of graphene in the ascending order of the lengths (see Table 2). Note that mirror reflection of AB and BA bilayers in 2AB+2BA films translates to symmetry properties of Pz(r0). In particular, for a reference frame with x-axis along zigzag and y-axis along armchair, Pz(x0, y0) = Pz(−x0, y0), as the zigzag axis is perpendicular to the vertical mirror symmetry plane conserving Pz, and Pz(x0, y0) = −Pz(x0, −y0), as mirror reflection in the armchair axis is equivalent to exchange of AB ↔ BA stackings in each bilayer and consequently leads to inversion of Pz. This requires BΛ = 0 in eq 4 for 2AB+2BA films, whereas the other Fourier coefficients are listed in Table 2. In Figure 3, we illustrate the resulting real space distributions of the ferroelectric polarization in a representative region of moiré superlattice for each of the considered tetralayer films. The mirror reflection in armchair axes prescribes zero polarization for DWs in 2AB+2BA films, shown in Figure 2, oriented along armchair crystallographic axes and separating domains with opposite polarization, Pz = ±0.07 e/μm. This would correspond to the transfer of ∼1011 cm–2 electron density between the tetralayer film surfaces. Note that close to domain corners polarization reverses sign and its magnitude become four times higher than those in the middle of domains. DW intercrossings possess opposite polarization compared to the inner part of the DW. In contrast, supercells of 3ABA+1 (3ABC+1) twisted tetralayers contain one inversion symmetric domain, ABAB (ABCA) (see Figure 2), where ferroelectric polarization is forbidden by symmetry and another, ABAC (ABCB) domain (as in 2AB+2BA films), which carries an out-of-plane polarization. We note that in 3ABA+1 (3ABC+1) tetralayers DWs and their intercrossings host even an order of magnitude stronger local vertical charge transfer than in the inner areas of the domains. The spontaneous out-of-plane polarization that we predict here for few-layer rhombohedral graphene structures is a feature induced by an asymmetrically placed twin boundary which breaks inversion and mirror symmetry in the system, and we trace it to the effect of individual terms in the full Sloczewski–Weiss–McClure model of graphite involved with such asymmetries and responsible for electron–hole symmetry breaking of the single-particle spectra. The reported calculations show that the size and orientation of Pz are critically affected by the intrinsic screening inside the graphitic films, leading to the opposite signs of Pz in tetralayers analyzed with and without screening and much smaller values, about 1 order of magnitude smaller than those experimentally reported in WTe2, MoS2, or hBN. We also note that in marginally twisted multilayers, the largest local polarization appears to be at domain walls and the near the nodes of domain wall networks. Overall, such sensitivity to the choice of the values of input parameters in the tight-binding model of graphite, number of layers in the film, and variation across the moiré pattern in marginally twisted structures may be used to experimentally refine the parametrization of the Sloczewski–Weiss–McClure model.