{'fraction_non_alphanumeric': 0.039720977302298684, 'fraction_numerical': 0.01951713170449617, 'mean_word_length': 4.439158639669746, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Electric field-tunable layer polarization in graphene/boron nitride twisted quadrilayer superlattices

The recently observed unconventional ferroelectricity in AB bilayer graphene sandwiched by hexagonal Boron Nitride (hBN) presents a new platform to manipulate correlated phases in multilayered van der Waals heterostructures. We present a low-energy continuum model for AB bilayer graphene encapsulated by the top and bottom layers of either hBN or graphene, with two independent twist angles. For the graphene/hBN heterostructures, we show that twist angle asymmetry leads to a layer polarization of the valence and conduction bands. We also show that an out-of-plane displacement field not only tunes the layer polarization but also flattens the low-energy bands. We extend the model to show that the electronic structures of quadrilayer graphene heterostructure consisting of AB bilayer graphene encapsulated by the top and bottom graphene layers can similarly be tuned by an external electric field.

Moiré of moiré patterns are a consequence of the interference between two different bilayer moiré patterns, and they have length scales that are generally orders of magnitude longer than the first-order moiré length. As a result of the higher-order interference pattern, such a system lacks a periodic approximation even in the continuum limit, and consequently has no approximate Brillouin zone in reciprocal space [24,25,28,29]. The lack of periodicity and the large system size present significant challenges to the usual approaches of developing an accurate theoretical model. The electronic band structure of these moiré of moiré systems has been calculated using a momentum-space low-energy continuum model for twisted trilayer graphene [22,25] and twisted hexagonal Boron Nitride (hBN)/monolayer graphene/hBN heterostructures [28,29]. This approach is both computationally efficient and accurate for the low-energy features of interest.
By encapsulating monolayer or Bernal (AB) bilayer graphene in hBN, a moiré of moiré pattern can be created from the interference between two bilayer moiré superlattices: one from the twist angle, and the other from the ∼1.7% graphene-hBN lattice constant mismatch [1,[28][29][30][31][32][33][34][35][36][37]. In particular, unconventional ferroelectricity has been experimentally observed in heterostructures consisting of an AB bilayer graphene sandwiched between top and bottom hBN layers [1], in contrast to the conventional ferroelectricity in other 2D vdW heterostructures due to a net structural polarization [38][39][40][41][42][43][44][45][46]. In the graphene/hBN heterostructures, hysteretic loops in the top-gate/back-gate phase-space have been observed. The top-gate voltages corresponding to resistance peaks stay constant for a wide range of back-gate voltages, which is referred to as Layer Specific Anomalous Screening (LSAS). However, the ability to confirm the electronic origin of the observed ferroelectricity has been hindered by the lack of an accurate microscopic theoretical model. In this work we aim to provide a singleparticle low-energy continuum model for this quadrilayer graphene/hBN heterostructures and provide a theoretical foundation to understand the observed ferroelectricity in this system. We find that the graphene/hBN heterostructure exhibits a moiré-induced layer polarization that sensitively depends on the twist angle. The layer polarization can be tuned by the application of an external displacement field. This moiré-induced layer polarization helps explain the observed unconventional ferroelectricity and LSAS. In addition, we extend the model to study a similar quadrilayer system, AB bilayer graphene encapsulated by top and bottom monolayer graphene, whose electronic structure is also electric field-tunable.
The paper is organized as follows. In Section II, we discuss the system geometry and review the low-energy continuum model for moiré of moiré systems and generalize to twisted quadrilayer graphene/hBN heterostructures. In Section III, we present the main results and show the electric field-tunable layer polarization of the valence and conduction bands. In Section IV, we extend our low-energy continuum to study the electronic structure of quadrilayer graphene heterostructures. We summarize our findings and make connections to experimental observations in Section V.

A. Geometry
The real space geometry of the system is shown in Fig. 1. The top layer, layer 1 (L1), and the bottom layer, layer 4 (L4), are twisted with respect to the sandwiched AB bilayer graphene by θ 12 and θ 34 respectively, with θ 12 in the clockwise direction and θ 34 in the counterclockwise direction. The monolayer lattice constants are defined as the column vectors of the following matrix: where i = Gr/hBN, and a Gr = 2.4768 Å and a hBN = 2.5189 Å are the graphene and hBN lattice constants respectively. Defining the counter-clockwise rotation ma- the lattice vectors of the 4 layers are given as The monolayer reciprocal lattice vectors are given by the column vectors of G = 2πA −T for = 1, 2, 3, 4. In terms of G , the bilayer moiré Brillouin zones are spanned by the column vectors of , where = 1, 3 for the top or bottom Gr/hBN moiré interface respectively. The moiré supercell in real space between layers l and l + 1 is A , +1 = 1/(2π)G −T , +1 . The norm of the column vectors are the moiré length due to the twist and lattice mismatch, which is given by where δ = a hBN /a Gr = 1.017 is the lattice constant mismatch ratio between graphene and hBN [47]. Each moiré pattern of the two interfaces, L1/L2 and L3/L4, exhibits a single coherent moiré length (see the interference patterns between L1, L2 and L3, L4 in Fig. 1(a)). These two moiré supercells have different length scales and are rotated relative to each other, and their interference pattern forms the more complex higher-order moiré of moiré patterns ( Fig. 1(a)). Figure 1(b) shows a perspective view of the system. We assume that the top and bottom hBN layers are rotated by 180 • such that in the absence of the twist, the system has a mirror symmetry along the z-direction. With a twist angle, however, the stacking configuration varies in space and we expect all possible stackings to occur in the system, and as a result, the inversion symmetry is broken. We checked that when the top and bottom hBN layers are not rotated, the reported results have no qualitative difference, except for the slight asymmetry under the exchange of θ 12 and θ 34 . ical quadrilayer system, qualitatively similar to the graphene/hBN system. Note that due to broken sublattice symmetry, the Brillouin zone corners of the hBN monolayers (top and bottom layers) are X and Y and they are nonequivalent, unlike the K and K points in graphene (middle layers). There are two moiré Brillouin zones that correspond to the top and bottom graphene-hBN interfaces (dark red and blue hexagons in Fig. 2). These two moiré Brillouin zones are incommensurate, and thus the system has no overall Brillouin zone. In the band structures of graphene/hBN heterostructures in the rest of the paper, we plot the energy eigenstates along the high-symmetry line in the moiré Brillouin zone that corresponds to the L1/L2 interface (red dashed lines in Fig. 2).

B. Low-energy continuum model
In previous work we have derived a low-energy momentum space continuum model for the electronic structure of twisted trilayer graphene with two independent twist angles [25]. We extend the model to study the quadrilayer heterostructures. We note that we do not consider the effect of in-plane relaxation in this work, which forms domain walls in the small-angle limit [24,28,48]. Domain wall formation may affect the polarization of the system, especially the spatial distribution of the ferroelectric domains, and warrants future studies. The Hamiltonian can be formally written as the following 4 × 4 matrix: where q is a low-energy momentum around which the Hamiltonian is centered (hereafter referred to as the center site), the diagonal blocks are intralayer terms of the monolayers, and the off-diagonal blocks, H ij are the terms that describe the interlayer interactions. Note that the Hamiltonian describes a single valley because the inter-valley separation is much higher in energy than the relevant low-energy degrees of freedom. When θ 12 and θ 34 are not equal, the two valleys are nonequivalent because they are not related by a reciprocal lattice vector. However, we checked that the qualitative behaviors, including the layer polarization and band flatness, of the two valleys are the same and thus we focus on the K valley in the rest of the work.
For the intralayer terms, we take the Dirac Hamiltonian as the monolayer graphene Hamiltonian where σ x and σ y are Pauli matrices, and v F = 0.8 × 10 6 cm/s is the Fermi velocity obtained from DFT [49]. We adopt the hBN Hamiltonian from Moon and Koshino [50]: where the onsite energies V B = 3.34 eV and V N = −1.40 eV capture the sublattice asymmetry of hBN, with the assumption that all the graphene onsite energies are set to V C = 0 eV. This approximation is justified when the twist angle is small because the hBN potential is far from the low-energy features of interest near the graphene Dirac point.
The interlayer terms between the graphene -hBN interface, H 12 and H 34 , can be derived in a similar way as twisted trilayer graphene by performing a Fourier transform of the real space tight-binding Hamiltonian and taking a low-energy limit [25]. They are given by where q ( ) is the low-energy momentum space degree of freedom in layer , expanded around K for graphene and Y for hBN monolayers, and the symbols {α, β} label the sublattice degrees of freedom. We set the first scattering Colors of the band structure represent the projection of the wavefunction weights onto the center site along the high symmetry line for each layer. The color scale is indicated in the panel on the bottom left side of (a), with yellow, red, blue, green corresponding to layers 1-4 respectively. Transparency is inversely proportional to the magnitude of the projected wavefunction weights. The high symmetry line is indicated by the red dashed lines in Fig. 2.
The other two scattering directions are then generated through rotations q , +1 where the rotation R(θ) matrix is defined in Eq. (2). The tunneling terms T n that corresponds to each q , +1 n are given as follows, wherez is the complex conjugate ofz, φ = exp(2πi/3),φ = exp(−2πi/3), and u 0 = 0.152 eV [50]. Equation (7) dictates how momentum degrees of freedom between layers and + 1 are coupled, that is, it prescribes the non-zero matrix elements of H , +1 . Unlike in bilayer graphene, these degrees of freedom no longer form a simple momentum space lattice due to the interference between the two bilayer moiré patterns. Namely, these two moiré patterns arise from the lattice-mismatched (and possibly twisted) interfaces between the top hBN and AB bilayer graphene, and between the bottom hBN and AB bilayer graphene. Unless both interfacial twist angles are identical, these two moiré lengths are different, and the two moiré cells are rotated from each other. As a result, the interference pattern can be dominated by higher-order harmonics and is generally incommensurate even in the continuum limit [24,25,51].
The off-diagonal term H 23 represents the interlayer coupling between AB-stacked bilayer graphene, which is approximated to be parabolic: where ν = ±1 denotes the valley degree of freedom, the parameter γ 1 represents the band splitting and v 3 describes trigonal warping [52]. We take γ 1 = 0.34 eV and v 3 = 0.051 × 10 6 m/s [50]. When a vertical displacement field is applied, each layer has a different external potential energy, which modifies the intralayer terms. With a displacement field D, the total potential difference across the 4 layers is V = 3Dd, where d = 3.35 Å is the interlayer spacing. The potential energy of each layer is then Φ 1 = −eV /2, Φ 2 = −eV /6, Φ 3 = eV /6, Φ 4 = eV /2, and the intralayer term for each layer has an additional term that is Φ 1 where 1 is the 2 × 2 identity matrix. In this way, the positive electric field direction is upward-pointing from L4 to L1.

A. Electronic structures at zero electric field
In the absence of an external electric field, graphene/hBN heterostructures exhibit an intrinsic layer polarization that is dependent on the twist angles. Figure 3 shows the band structure and DOS of the graphene/hBN heterostructures for three different sets of twist angles. The DOS is projected onto the center site q along the high symmetry line. Without a twist angle, Fig. 3(a)-(b) shows that neither the valence band nor the conduction band exhibits layer polarization. The low-energy bands are parabolic and the system is gapless, similar to those of AB bilayer graphene. We compare Fig. 3(a) with the band structure of AB bilayer graphene with aligned hBN on one side in Ap-pendix B, in which case the low-energy bands near the charge neutrality point (CNP) exhibit some layer polarization from the hBN layer. In Fig. 3(a), note that the low-energy valence and conduction bands cross near the K point. This is analogous to trigonal warping in AB bilayer graphene: when the tunneling between non-dimer sites is included, the parabolic bands are split into four Dirac cones with linear dispersion [52,53]. Note that in AB bilayer graphene, trigonal warping is very close to the CNP and the saddle point is estimated to be on the order of 1 meV [54]. Here, while we do not include the non-dimer tunneling in the AB bilayer graphene Hamiltonian, the coupling between graphene and hBN serves the same role as remote hoppings and causes trigonal warping. When L4 is twisted while L1 is untwisted (Fig. 3(c)-(f)), bands near the CNP become clearly layer polarized: the valence band is L3/L4 polarized whereas the conduction band is L1/L2 polarized. When the two twist angles are equal in magnitude (Fig. 3(g)-(h)), the layer polarization again disappears and the gap closes because the existence of a gap and layer polarization requires the broken inversion symmetry between the top and bottom interfaces. Figure 4 shows the DOS near the CNP as a function of θ 34 while θ 12 is kept fixed at 1 • (note that here L1 and L4 are twisted in the opposite directions). Except for the θ 12 = θ 34 = 1 • case, Fig. 4(a) shows that the DOS is layer polarized. When the inversion symmetry is broken, the system is gapped, and the gap size increases away from θ 12 = θ 34 and reaches an asymptotic value of ∼ 16 meV when θ 34 3 • . Moreover, the layer polarization switches as θ 34 crosses 1 • : the valence band is L1/L2 polarized when θ 34 < 1 • and L3/L4 polarized when θ 34 > 1 • , suggesting that the valence band is polarized by the interface with a larger twist angle. Figure 5 shows the gap size between the valence and conduction bands, ∆, at the CNP in the (θ 12 , θ 34 ) phase space. When θ 12 < 0 • , the top and bottom layers are twisted in the same direction. Note that ∆ is symmetric around θ 12 = ±θ 34 , θ 12 = 0, and θ 34 = 0 (the slight asymmetry is due to numerical noise). When θ 12 = ±θ 34 , the gap disappears (∆ = 0) because of inversion symmetry. The system is gapped when this inversion symmetry is broken and when at least one of the twist angles is small (|θ 12 | 1 • or |θ 12 | 1 • ). When both twist angles are greater than 1 • , the gap size ∆ is drastically reduced. This is because when both twist angles are large, the top and bottom hBN layers are essentially decoupled from the AB bilayer graphene layer, and AB bilayer graphene itself has parabolic bands and is gapless. When one of the twist angles is small and the two twist angles are not equal, as the other twist angle increases, the gap size saturates to ∼ 16 meV as shown in Fig. 4.

B. Wavefunction localization
To obtain the layer-projected real space wavefunction distribution at a position r, we perform an inverse Fourier transform by summing over the wavefunction weights that correspond to each momentum degree of freedom q ( ) on layer l: where n is the band index, q is the center site of the Hamiltonian, and ψ n,q,α (q ( ) ) is the wavefunction weight that corresponds to momentum q ( ) and sublattice α. We sum over a total number of 5840 momentum degrees of freedom in the basis (values of q ( ) ). Figure 6 shows the K-point (q = K) valence and conduction band wavefunction distribution in real space for (θ 12 , θ 34 ) = (0 • , 1 • ) and (θ 12 , θ 34 ) = (0 • , 5 • ). In both cases, valence band wavefunctions localize at the moiré supercell that corresponds to θ 34 whereas the conduction band wavefunctions localize at the 0 • supercell from the lattice mismatch. This is consistent with the band layer polarization in Fig. 3 (c) and (e) -since the valence bands are L3/L4-polarized, the wavefunctions are also expected to have the length scale of the L3/L4 supercell, and vice versa for the conduction bands. We could consider the wavefunction distribution as the summation of localized electrons at the moiré scale and a uniform background. The scale of the localization depends on the interface while the uniform background can be attributed to band dispersion.

C. Electric field dependence
The layer polarization of the conduction and valence bands can be tuned by an out-of-plane displacement field. Figure 7 shows the band structure for (θ 12 , θ 23 ) = (−0.1 • , 15 • ) for four different values of the applied electric field. In absence of an external electric field (Fig. 7(c)), as shown in Section III A, the valence band is polarized to the interface with the larger twist angle, which is L3/L4. The polarization gives rise to a positive or upward-pointing intrinsic electric field, from the bottom L4 to the top L1. When a negative electric field is applied, there is a critical value of which the applied electric field cancels out the intrinsic polarization. In this case, the critical electric field is V c = −0.04 eV/nm (Fig. 7 (b)). In general, the value of the critical electric field depends on the single-particle gap size ∆ (Fig. 5), and the larger ∆ is, the larger the V c is. As the negative electric field increases, the polarization reverses (Fig. 7(a)) and the valence band becomes L1/L2 polarized while the conduction band is L3/L4 polarized, dominated by the direction of the external electric field. As the magnitude of the electric field becomes large, for both positive and negative electric fields, the valence and conduction bands flatten significantly and the gap increases ( Fig. 7(a)(d)). Figure 8 shows the DOS of (θ 12 , θ 34 ) = (−0.1 • , 15 • ) as a function of displacement fields. The projected DOS onto the L1/L2 and L3/L4 interfaces in Fig. 8 clearly shows the polarization reversal at the critical electric field −0.04 eV/nm. Both the gap size and the DOS maximum smoothly increase as the applied electric field increases. Moreover, we note that the DOS has a particle-hole asymmetry at large electric fields. For example, at a large positive electric field (say 0.25 eV/nm), the valence band peak clearly has a larger magnitude than the conduction band.
In addition to the band structures and DOS, the exter- nal electric field also changes the wavefunction localization. Figure 9 compares the valence and conduction band wavefunctions with positive and negative applied electric fields. When V < V c (Fig. 9(a)-(b)), the valence bands have the length scale of the |θ 12 | = 0.1 • moiré supercell, while the conduction band wavefunctions localize at the much smaller 15 • moiré supercell. This agrees with the layer polarization reversal by the negative electric field as shown in Figs. 7 and 8. When V > V c , Fig. 9(c)-(d) shows that the wavefunctions localize at the opposite length scales as in (b)-(d), which is also the same scales as the zero-electric field case because the applied electric field has the same direction as the intrinsic electric field, and thus there is no polarization reversal. Note that we do not include the dielectric permittivity of the hBN and graphene when modeling the effect of the electric field. Including the dielectric effect changes the critical electric field where the gap closes but is not expected to affect the qualitative behaviors discussed in this section including the gap closing, polarization reversal, and band flattening.

IV. QUADRILAYER GRAPHENE HETEROSTRUCTURES
We could employ the low-energy continuum model to study quadrilayer graphene heterostructures, with AB bilayer graphene encapsulated by top and bottom graphene layers. The intralayer terms for monolayer hBN, H hBN in Eq. (4) are replaced by rotated Dirac Hamiltonian, H Gr , which represents a slight modification to Eq. (5) and is defined as follows where σ θ x = σ x cos θ −σ y sin θ , σ θ y = σ x sin θ +σ y cos θ are rotated Pauli matrices with θ 1 = θ 12 , θ 2 = 0, θ 3 = −θ 34 . The interlayer interactions between the top and middle layers, H 12 and H 34 , have the same form as Eq. (8), but the three scattering vectors q , +1 j for j = 1, 2, 3, = 1, 3 are given by q , +1 , and q , +1 In addition, the scattering matrices T j that correspond to each q , +1 j are defined in the same way, where ω 0 = 0.07 eV and ω 1 = 0.11 eV/nm are the interlayer hopping parameters between AA and AB stackings respectively and the difference between ω 0 and ω 1 accounts for the out-of-plane relaxation [55][56][57]. Figure 10 shows the electronic structure of the quadrilayer graphene heterostructures with (θ 12 , θ 34 ) = (1.6 • , 3.0 • ). Unlike the hBN/graphene heterostructures, the quadrilayer graphene system exhibits no significant layer polarization or band gap in absence of an external electric field. The application of the electric field introduces layer polarization to the bands near the CNP, and the induced layer polarization agrees with the applied electric field direction ( Fig. 10(a)-(f)). In addition, for certain values of applied electric fields, (for this twist angle combination when V = 0.05 eV/nm), the low-energy bands near the CNP become significantly flatter and the symmetry between the valence and conduction is also broken. However, the system remains gapless under the applied electric field.

V. DISCUSSION AND CONCLUSION
In the experiment by Zheng et al. [1], the significant feature associated with the ferroelectric switching is an LSAS, namely, under certain conditions, the top (bottom) layer can completely screen the electric field from the top (bottom) while the electric field can still penetrate the bottom (top) layer. Such an asymmetric screening scenario points to vastly different electronic structures associated with different layers. Our calculations show that the twist angle difference between the top BN/graphene interface and the bottom BN/graphene interface could lead to layer-asymmetric electronic structures, which provides the starting point to understanding such unusual LSAS behaviors. In particular, if one interface (top) has a longer moiré wavelength and the other interface (bottom) has a shorter moiré interface (Figs. 8  and 9), under certain conditions the top layer will feature electronic structures with a large DOS and localized wavefunctions but not the bottom layer. This can support a scenario where the top layer may feature a stronger correlation effect and unique screening properties.
Based on the information from the band structure, we could provide a rough estimate of the intrinsic polarization as follows. If we fill the conduction band with one electron per the moiré cell that corresponds to the 0 • interface, according to, for instance, Fig. 3(e) and the weight distribution on different layers, effectively all the weights are concentrated at the L1/L2 interface (99.7%) and we can treat the electron to be purely localized at the L1/L2 interface. As a result, the dipole moment can be calculated as p = 3ed, where e is the electron charge and d = 3.35 A is the interlayer distance. The electric polarization is then P = p/(3A m d), where A m is the moiré superlattice area and A m = 1.69 × 10 −12 cm 2 for the 0 • interface. Using these, we estimate P = 0.095 C/cm −2 . Despite being a crude estimate that ignores the dielectric environment, this is on the same order of magnitude with the experimentally measured polarization by Zheng et al. [1], which is 0.12 C/cm −2 .
In addition to the hBN/graphene heterostructure, we showed that the low-energy band structure of the quadrilayer graphene heterostructures can be tuned by external electric fields in a similar way. While at the zero electric field the low-energy moiré bands do not exhibit moiré induced intrinsic polarization, an external electric field not only alters the band flatness but also breaks the symmetry between the valence and conduction bands. If the observed ferroelectricity in hBN-sandwiched AB bilayer graphene has an electronic origin, we expect the quadrilayer graphene system to exhibit similar hysteretic behavior in transport measurements. To properly normalize the DOS, we first calculate the DOS of the intralayer Hamiltonian only, which reduces to four independent monolayers. Near the charge-neutrality point, the DOS per eV per nm 2 is given by [58] where the prefactor includes a factor 4 from spin and valley degeneracies and D is system-dependent. For hBN/graphene heterostructures, D = 2 because only two graphene monolayers contribute to the low-energy DOS and the monolayer hBN is estimated to be a constant potential that is high in energy. For quadrilayer graphene heterostructures, D = 4 to account for the contribution from all 4 layers. We then obtain a normalization constant by fixing the prefactor to the expected slope given in Eq. (A1) and use the same constant for the DOS of the full Hamiltonian. The density of states (DOS) is obtained by discretizing the moiré Brillouin zone between L1 and L2 by 61 × 61. We adapt the Gaussian smearing full-widthhalf-maximum σ based on θ 12 because the moiré Brillouin zone area being sampled increases as the twist angle increases. For graphene/hBN heterostructures, θ 12 < 0.8 • , σ = 1.7 meV, θ 12 < 2 • , σ = 2.3 meV, θ 12 2 • , σ = 2.7 meV. For quadrilayer graphene heterostructure with θ 12 = 1.6 • , θ 34 = 3 • , we choose σ = 1.2 meV. Each point q on this BZ grid is then expanded into a truncated momentum-space Hamiltonian of dimension 2028 × 2028, which is equivalent to truncating the basis to the second shell of the monolayer lattice vectors before taking the low-energy expansion (namely, we constrain the magnitude of k ( ) = q ( ) + K K /Y ). We confirmed that the features of interest have converged with respect to the momentum-basis cutoff radius.  Figure 11 shows the electronic band structure of AB bilayer graphene on top of a monolayer hBN layer. The band structure agrees with the results in Moon and Koshino [50]. Unlike in the hBN-sandwiched AB bilayer graphene without a twist (Fig. 3(a)), this trilayer system is gapped because of the broken inversion symmetry due to the hBN monolayer. The low-energy bands have hBN characteristics because of the net potential from the hBN -we approximate the hBN potential to be constant, V B and V N , leading to a net constant potential of ∼ 2 eV and thus polarizes the low-energy bands. In contrast, the net hBN potential is canceled out in the hBN-sandwiched graphene heterostructures. Therefore, in the quadrilayer graphene/hBN heterostructures, in absence of a twist angle, we expect the system to be layer unpolarized, and the twist angle introduces moiré-induced layer polarization.