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.


tonian takes the form
where we introduce the Sloczewski-Weiss-McClure (SWMcC) parameters 1-3 , v ≈ 1.02 · 10 6 m/s, Hamiltonian of tetralayer graphite with a marginally twisted interface Local stacking characterised by an interlayer offset, r 0 at the twsited interface. The following Hamiltonians generalise those derived in Ref. 5 . For 2AB+2BA, Hamiltonian reads where we introduce K-valleys momenta, K

Convergence analysis
In our analysis, we compute the electron density in each layer by brute-force diagonalisation of the Hamiltonian (1) for a sufficiently dense grid of points in reciprocal space within a circular area with cut-off wavevector, k c . For each nABAm film, as k c increases, the ferroelectric polarisation converges to a P z -value, displayed below in the left panel of Fig with deeper bands, contributed to P z in negative direction of z-axis. As shown Fig. 1 (a), it is necessary to sample an area of 5 · 10 −2 Å −2 (or a value for the cut-off energy of ∼ −0.13 eV) to reach convergence.
In contrast, the low-energy dispersion of 4ABA3 exhibits two pairs of flat bands [see inset in Figure 1: Value for the ferroelectric polarisation, P z , as a function of the radius of the circle centred at K ± in reciprocal space, k c taken to compute the electron density in each layer, for 4ABA (a) and 4ABA3 (b) tetralayer graphene. For the latter, it suffices to sample the area on which flat bands extend to obtained a converged value.

Self-consistent implementation of screening
Here, we extend the above-described SWMcC model to incorporate the effects of internal electric fields produced by the charge redistribution. These fields produce an energy offset in each layer, which can be captured in our model by adding a term ϵ i12 to the i-th diagonal block of the Hamiltonians in Eqs. (1), (2) and (3) using the formula 6 with d ≈ 3.35 Å being the interlayer distance in graphene stacks, ε z ≈ 2.6 6 the relative permittivity of graphene. We determine the Fermi level for neutral structures, and calculate the charge redistribution across the layers, taking into account only the eigenvectors associated to states below the Fermi level. Because the charge densities n i also depend implicitly on the on-layer potentials, 5 the expression in Eq. (4) needs to be solved self-consistently. In the first iteration, we set all onlayer energies to zero, ϵ (0) i = 0, which gives the value of polarisation for the unscreened system, P u z . In each iteration, we use Eq. (4) to compute a new set of on-layer energies,ε (n) i , which are used to update the value employed in the following iteration using where η determines the fraction of the value for the on-layer energies that is updated in each iteration, and was taken to be 4 · 10 −4 in our work. This linear-mixing algorithm is necessary in a broad range of self-consistent methods to avoid overshooting 7 .
In Figs. 2(a) and 2(b), we show the value of P z and electron densities in each layer as a function of the iteration for 1ABA and 4ABA3 films, highlighting importance of self-consistent calculations. Histograms in Fig. 2(c) show comparison of P z values for mABAn films with less than 8 layers computed with (right) and without (left) account of self-consistent screening. We also note that our results strongly depend on the input parameters, and suggest that measuring P z in these systems could be an experimental route to narrow down the possible values of the SWMcC parameters.

Analysis of parametric dependence
Equation (3) of the main text is based on the analysis shown in Fig. 3, where each panel represents dependence of P u z of ABCB tetralayer graphene on one of the SMWcC parameters. The blue dots are numerically computed values of P u z and the red lines are the linear interpolations, which we use to obtain the coefficients X 4,D,2,5 . We observe the linear dependence of P z on the symmetrybreaking parameters v 4 , ∆ ′ , γ 2 and γ 5 . The middle and rightmost panels of Fig. 3 (a), shows the dependence on γ 1 and v, when v 4 is kept as the only non-zero electron-hole symmetry-breaking parameter, which is proportional to γ 1 |γ 1 | × v −3 , respectively. Likewise, in the same panels of Figs. 3 (b), (c) and (d), we demonstrate that P u z as the power law |γ 1 | × v −2 when ∆ ′ , γ 2 and γ 5 Figure 2: Convergence analysis of P z and n i for (a) 1ABA and (b) 4ABA3 films. The latter evidences that, for larger films, electron charge primarily redistributes to the surface and twin boundary layers. (c) Electric polarisation of mABAn, P u z (P z ), displayed in the form of a histogram, using a model that neglects (accounts for) screening effects of electron-density redistribution. Figure 3: Numerical values (blue points) and fitting curves (red) computed to support the parametric dependence of the first, second, third and fourth terms in Eq. (2) in the main text are analysed in the triad of panels (a), (b), (c) and (d), respectively, for ABCB tetralayer graphene. In the leftmost panels, we demonstrate the linear dependence of each SWMcC parameter that induces ferroelec- Figure 4: Dependence of P u z (black) and P z (blue) as a function of each individual SWMcC parameter responsible for electron-hole asymmetry. Vertical red dotted lines lie on the values given in 4 , which we use in our work. Adding all the values where they intersect the black and blue interpolating lines amounts approximately the value por P u z and P z obtained using the full SWMcC model, and shown as the first and last values in Fig. 2(a).
are kept as the only non-zero parameters, respectively.
In Fig. 4, we extend analysis of P z when effect of self-consistent screening is taken into account. In particular, we compare the leftmost panels of Fig. 3 using a model with (blue dots) and without (black dots) the effects of screening, and demonstrate that, while X 4 , X D and X 5 reduce their values by a factor of ∼4, 1.4 and 2.2, respectively, the coefficient X 2 changes sign. As stated in the main text, this peculiar behaviour plays a prominent role in the change of sign of polarisation, being P u z ≈ 0.54 e/µm in the unscreened model and P z ≈ −0.07 e/µm in the screened model [see Fig. 2(a)]. It is also worth mentioning that adding the individual contributions to P z (or P u z ) from the four symmetry-breaking terms amounts approximately to the to computed using the full SWMcC model with (without) the account for screening effects.

Lattice relaxation in twisted tetralayer structures
To describe lattice reconstruction of moiré superlattice formed at the twisted interface we introduce  Table 1: Main non-zero Fourier series coefficients (in ) of displacement fields describing reconstruction in twistronic tetralayers 3ABA+1, 3ABC+1 and 2AB+2BA structures with θ = 0.05 • , which were used in Fig. 2 of the main manuscript.