Close Banner |

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

Nature Communications volume 13, Article number: 3103 (2022 ) Cite this article

Two-dimensional (2D) bilayers, twisted to particular angles to display electronic flat bands, are being extensively explored for physics of strongly correlated 2D systems. However, the similar rich physics of one-dimensional (1D) strongly correlated systems remains elusive as it is largely inaccessible by twists. Here, a distinctive way to create 1D flat bands is proposed, by either stamping or growing a 2D monolayer on a non-Euclidean topography-patterned surface. Using boron nitride (hBN) as an example, our analysis employing elastic plate theory, density-functional and coarse-grained tight-binding method reveals that hBN’s bi-periodic sinusoidal deformation creates pseudo- electric and magnetic fields with unexpected spatial dependence. A combination of these fields leads to anisotropic confinement and 1D flat bands. Moreover, changing the periodic undulations can tune the bandwidth, to drive the system to different strongly correlated regimes such as density waves, Luttinger liquid, and Mott insulator. The 1D nature of these states differs from those obtained in twisted materials and can be exploited to study the exciting physics of 1D quantum systems.

Recently, twisted bilayer graphene (TBG) at magic angles1,2 and other van der Waals (vdW) heterostructures at small twist3,4 have garnered great attention as material platforms for realizing 2D correlated physics with an unprecedented level of control. Several interesting electronic phases have been observed in these systems, such as correlated insulator2,4, superconductivity1,5, non-trivial electronic topology6, and magnetism7,8. Physically, these emergent phases can be attributed to the existence of the Bloch flat bands9,10, where the kinetic energy scale is quenched, and the role of electronic interactions is enhanced. The flat bands originate from the perturbation of the electronic structure by the long-wavelength superlattice (moiré) period, which suppresses the group velocity in TBG at magic angles11,12 and creates electronic confinement in other vdW heterostructures13,14,15. The moiré periods arise from either lattice mismatch or rotational misalignment between the layers with fine-tuning of the twist, posing challenges16 in fabrication, variability between devices, and scalability. Moreover, examining the rich unexplored physics of 1D strongly correlated systems, likewise, is largely inaccessible by twists.

Below we describe creating flat bands through an alternative route, not requiring a twist angle. The strategy involves either growing17 or stamping18 a 2D material on a topographically patterned substrate with non-zero Gaussian curvature, that is non-Euclidean surface, Fig. 1a. To conform to such surfaces, a planar 2D crystal must deform, so the undulated topography imparts strain. Strain is known to perturb the crystal Hamiltonian through a deformation potential19, to a magnitude proportional to the strain, which in turn is determined by the topography of the surface. A periodic strain modulation will create a confining potential, which–if strong enough–can localize electrons and result in modulated super-lattice band whose bandwidth depends on the surface’s specific geometry. Hence, patterned surfaces with specific topography can, in general, create and fine-tune flat bands in any 2D semiconductor material. To the best of our knowledge, creating either 2D or 1D flat bands in monolayer semiconductors by undulation has not been discussed before.

a A 2D material stamped on a topography with bi-periodic sinusoidal height-modulation, causing in-plain strain. b Relaxed h-BN geometry and strain field ε ≡ uii at undulation aspect ratio A = 0.079. The left panels in c, d map the pseudo-magnetic (BP) and pseudo-electric potential (φp) fields, respectively; the right panels show the values of the respective fields along x = const, marked by a white dashed line on the left.

We illustrate this idea by theoretically investigating the electronic properties of monolayer hexagonal boron nitride (hBN) deformed by a bi-periodic sinusoidally modulated topography. Interestingly, strained hBN attains both pseudo-electric field EP (by virtue of deformation potential) and also (having a honeycomb lattice, like graphene, with two inequivalent basis atoms) pseudo magnetic field BP20. We find that for bi-sinusoidal deformation, EP and BP have very different spatial dependence. A combination of these leads to anisotropic confinement and creates one-dimensional (1D) flat bands, whose bandwidth can be varied by the surface topography. The 1D nature of these states can be exploited to probe the exciting physics of one-dimensional quantum systems, which have been predicted to exhibit interesting effects21 such as Luttinger liquid behavior, charge and spin density waves, Peierls instability, and deviation from Fermi-liquid theory. The origin and nature of these 1D states are different from the 2D flat bands observed in TBG and twisted vdW heterostructures, opening an exciting realm of exploring many-body effects in 1D quantum systems in a clean and controllable manner.

A 2D material conformed to a curved non-Euclidean surface undergoes a locally in-plain strain, which can be evaluated at a continuum level (for relatively smooth topography) by solving the second Föppl–von Kármán (FvK) equation22, Δ2χ = −Y(fxxfyy − fxy2). Here χ, Y, and f, are the Airy stress function, Young’s modulus, and the surface shape function, respectively. We consider bi-periodic sinusoidally modulated off-plane shape f(x,y) = hsinαx·sinαy, with α = 2π/L, akin to egg-cart, L and h are the undulation period, and amplitude, respectively. The displacement (ux, uy) and imparted strain (uij) fields were solved analytically from the FvK equation for the sinusoidal surface (for details see Methods section and Supplementary Note 1). The maximum tension or compression for this geometry depends on the aspect ratio A = h/L, as ε = |uii,max| = π2A2(1 − ν)/2, where ν = 0.31 is the Poisson’s ratio for hBN and uii = (uxx + uyy)/2. A relaxed atomic structure of a sinusoidal hBN can be further constructed from the solved displacements ux and uy, as in Fig. 1b, for A = 0.079, where color shows uii (for details see Supplementary Note 3). As expected, the hilltops and valley-bottoms are stretched while the saddle areas are compressed, to an amplitude ε~2.12%. We note here that, generally, the material strain depends on both the surface shape f(x,y) and boundary conditions (in case of growth, also on the chemical potential μhBN, controlled by the growth conditions). Here we allow full relaxation to minimize elastic energy while accommodating the substrate topography, that is no forces at (remote) layer’s perimeter and no friction to the substrate. This corresponds to “stamping” the 2D material onto the frictionless matrix, when the layer contracts laterally, with non-negligible displacements ux, uy (see Supplementary Note 1).

Analogous to graphene23,24,25,26, the strain in hBN generates pseudo- electric and magnetic fields20,23, significantly perturbing the crystal Hamiltonian. The low energy effective Hamiltonian in strained hBN in the vicinity of the K points is given by

where vF = 3|t|a/2, |t| is the nearest-neighbor (NN) hopping amplitude, and a is the interatomic distance. σ(τ)=(τσx,σy,σz) are defined in terms of the three Pauli matrices, and τ = +1 (−1) for K (K’). k = (kx,ky,Δ), where ħvFΔ is the difference in sublattice potential between B and N atoms, and kx,y is the electron crystal momentum measured relative to K or K’. Ap is the pseudo-vector potential caused by shear, Ap = (β0/√2a)[(uxx − uyy)/2, − uxy], where β0 = (a/t)∂t/∂a = −3.3. φp is the pseudo-electric potential (PEP) arising due to the hydrostatic component of strain, φp = −guii, where g ≈ 3.66 V27. Accordingly, these potentials generate pseudo-electric field (PEF) Ep = −∇ φp and pseudo-magnetic field (PMF) Bp = Bpz, where Bp = (ħ/e)(∂xAp,y − ∂yAp,x) and z the unit z-vector. One can already recognize that these additional pseudo fields in the Hamiltonian, arising due to strain, act as a perturbing confinement potential.

The strain fields obtained for sinusoidal surfaces allow us to derive the analytical expressions for pseudo electric and magnetic fields (for details see Supplementary Note 2).

Figure 1c, d shows the PMF and PEP for A = 0.079, and L = 6.35 nm. The spatial dependence for both fields is different and surprisingly, PMF depends only on y (Fig. 1c). It is known28 that periodic magnetic fields can lead to confinement and create localized electronic states. Similarly, we expect that for sinusoidally modulated hBN, a combination of both PEP and PMF will create flat bands. Sections of these fields along y-direction, at x = const are plotted in the right panels of Fig. 1c, d. The periodic PMF has an amplitude of Bp,max~420T, which corresponds to confinement energy of ~2μBBp,max = 49 meV, while the periodic PEP corresponds to confinement energy of ~76 meV. We will show that different spatial dependence of PMF and PEP leads to anisotropic confinement and results in the interesting electronic nature of the flat bands.

We next calculate the electronic bands of our topographically-strained hBN, using density-functional based tight-binding (DFTB) theory with a local orbital basis29. DFTB has been successfully applied to study various forms of hBN15,30, for which DFT calculations are intractable (see Supplementary Note 4 for details).

Monolayer hBN honeycomb lattice is akin to graphene, yet the different basis atoms break the sub-lattice A-B symmetry, and an energy gap opens, making hBN an insulator. The undeformed monolayer hBN shows a band gap of ~3.55 eV. Figure 2a shows the band structure under bi-sinusoidal strain, with A = 0.079, ε = 2.12%, and L = 6.35 nm. The Brillouin zone is defined based on the shape function. One can see additional bands appearing in the gap, looking like defect states which might arise due to electronic confinement. We find that the bandwidth (W) of the lowest unoccupied states (shown in red) is W = 39 meV, which is very small, and it is a flat band; in comparison, the effective W of pristine BN bands corresponding to nearest neighbor hopping t~2.16 eV20 is W~4t = 8.6 eV, which is much larger than the W of the modulated flat bands. These flat bands are well separated by >100 meV from the other states at higher energies. Interestingly, the bands are dispersive along Г-X and almost non-dispersive along R-X, which corresponds to kx and ky directions, respectively. This makes these flat bands one-dimensional and very different from those seen in TBG and other twisted vdW heterostructures. The band decomposed charge density |Ψnk|2 in Fig. 2b corresponds to one of the eigenstates at Г point. The electronic states are delocalized along the x- but are completely localized along the y-direction, confirming these flat bands’ 1D nature (see Supplementary Fig. 6 for the charge density corresponding to the whole flat band, which is similar to Fig. 2b).

a Calculated electronic band structure of sinusoidally deformed hBN, A = 0.079. The marked red defect-like flat bands are due to electronic confinement caused by deformation. n = 1, 2 correspond to the first and second eigenstate of the confinement, respectively. The n = 1 and 2 bands are dispersive and connected along the Г-X direction, while they are flat along the R-X direction. The splitting along the R-X direction of both n = 1, 2 states is due to interactions along the x-direction. Yellow-shaded areas mark the bands of pristine hBN with gap in between. b The band decomposed charge density corresponding to one of the flat bands at Г point | Ψnk | 2 shows the one-dimensional electronic nature of the bands. c An enlarged view of the flat bands, red in (a), plotted along X-Г-Y. The solid red line is the fit with a 8-band tight-binding model.

The flat bands are composed of 8 electronic states and are localized mainly on the four extremes of the sinusoidal modulation in Fig. 2b. Figure 2c shows the enlarged view of the flat bands plotted along X-Г-Y. We find that the flat bands dispersion and charge modulation can be described by a simple “coarse-grained” 8-band tight-binding (TB) Hamiltonian,

where tx and ty are hopping amplitudes along x-, and y-direction, respectively (for details see Supplementary Note 5). The 8 bands arise from the states localized on the 2 maxima and 2 minima, and each of them being doubly occupied. The maxima and minima act as artificial “quantum dots”. Fitting this TB model to DFTB results gives |tx| = 8.9 meV (~W/4) and |ty| = 0.9 meV. The ratio of hopping along x-, and y-direction is |ty|/|tx| = 0.1, again a manifestation of one-dimensionality of the electronic states. This is quite surprising at first, because the strain pattern appears to be isotropic along the x and y-direction (Fig. 1b).

To gain microscopic insights into the reasons behind these flat bands’ 1D nature, we calculated the electrostatic potential along x (Fig. 3a) or y (Fig. 3b) while averaged along the other two perpendicular directions. The sharp features are due to the approximations used to evaluate the diverging potential near atomic sites (for details see Supplementary Note 4). The potential rapid oscillations are due to periodic atomic sites, but a long-range modulation can also be seen. Gaussian averaging extracts the long-range modulation (Vconf., red lines in Fig. 3a, b), which is very different along x- and y-directions: eVconf. along y-direction has a depth of ~500 meV, larger than mere ~9 meV along x. This anisotropic confinement is expected because of the different spatial dependence of PEP and PMF, Fig. 1c, d. The smaller Vconf. along x-direction is mainly due to contributions from PEP only, while the larger Vconf. along the y-direction is contributed by both PEP and PMF. This signifies that the long-range potential modulations are due to pseudo-electric and magnetic fields, providing the anisotropic confinement needed to maintain one-dimensional flat bands. Additionally, the larger confinement energy along y-direction results in lower hopping amplitude ty in our coarse-grained model.

Electrostatic potential along (a) x, and (b) y and averaged along the other two perpendicular directions. The solid gray lines correspond to the calculated averaged potential, while the thin red lines correspond to the long-range potential modulation extracted by Gaussian averaging. The + and − mark tensile and compressive regions.

We find that the width of these flat bands can be tuned by changing either the aspect ratio of the topography, to alter the strain, or the period L of the undulation. Figure 4 shows the variation of W as a function of L, and ε. W goes below 10 meV for L > 8 nm. Moreover, W is found to depend exponentially on both L, and ε, a good fit W ∝ exp(−0.56*L–0.34*ε) is obtained (blue surface in Fig. 4); this is expected since W effectively corresponds to the hopping amplitude between the coarse-grained sites, defined by hybridization/overlap of wave function between them. Since the strength of hybridization decreases exponentially with distance L, W is found to show the same dependence. Additionally, the exponential dependence of W on ε can also be understood from the well-known dependence of hopping energy (t) with strain20, t(a) = t0exp(−|β0 | (a/a0−1)), is the hopping amplitude at the effective bond length a, and (a/a0−1) is its strain. We point out that flat bands appear even at strain as low as ~1.5%, and these bands are well separated by >70 meV from the other unoccupied bands above, which makes them accessible to experiments. Moreover, the quasi-1D nature of these flat bands remains robust against local strain imperfections (for details see Supplementary Note 2d) as well as small misorientation between hBN and the substrate (for details see Supplementary Note 2c).

The flat band width W versus topography period L and strain ε. The blue surface shows the fitted expression W ∝ exp(−0.56*L–0.34*ε), while the empty and filled circles are computed values. The red surface shows the estimate of the on-site Coulomb energy U, which depends on period L. In the (L, ε) parameter-plane, the outside of a shaded gray is the U > W region where physics of strongly correlated 1D phases can be realized.

Interesting strongly correlated physics in one-dimension31 is expected when the ratio of on-site Coulomb interaction U (responsible for electronic correlation) and the hopping amplitude t is large, U/4t = U/W > 1. Since U depends inversely on length L, U∝ 1/L13, and t∝ exp(-|β0 | L)20, the condition U/W > 1 should be easily achievable for reasonable L. We estimate U for our systems as U = e2/2πκL13, where L is the length of periodic modulation, and κ = 4.7332 is the effective dielectric constant of hBN. The red surface in Fig. 4 shows values of U. U/W can be enhanced by either increasing the aspect ratio (strain) and/or increasing the periodic length (Fig. 4). The region in (L, ε) with U < W is shown by the gray shaded area in Fig. 4. To achieve strongly correlated phases, where U > W, topography with L and ε values lying outside of the gray area in Fig. 4 must be chosen. Depending on the ratio of U/W and band filling, one can expect31 different phases such as Mott insulator (MI), Luttinger liquid, bond ordered wave (BOW), and band insulator (BI). E.g., at small band filling, gradually changing U/W from 0 to 2 can change the electronic phases in order BI → BOW → MI. Accordingly, one should expect that changing the periodic topography will provide a unique control to drive the system to different strongly correlated regimes exhibiting interesting physics. To realize the strongly correlated physics, the flat bands (Fig. 2a) must be partially filled, perhaps by electrostatic doping, as routinely done for 2D materials, including twisted bilayer TMDs4, and graphene2. Long-range ordered quantum phases in 1D tend to get destroyed due to thermal fluctuations at finite temperature21, hence, isolated 1D systems are not ideal to realize 1D physics. Importantly, our predicted system with parallel 1D states (resulting from the quasi 1D flat bands in a 2D material, Fig. 2b) will suppress such fluctuations, due to the finite coupling between them21, and will strengthen the physical effects in 1D, which will help to achieve interesting physics in 1D at finite temperatures.

To realize our predictions in experiments, hBN may be overlaid or stamped (or possibly grown directly) on patterned substrate, for instance SiO2, having a band gap of ~8.9 eV33, larger than hBN---so that there are no unwanted hybridization between hBNs electronic states and the substrate. Fabricating bi-periodic sinusoidal modulation on silicon is challenging but has already been attempted34. Delamination from the substrate, if the desired strain level is high may be a concern; by comparing the surface pressure due to substrate with the maximal adhesion forces of hBN to SiO2 with γad~13 meV/Å2 (see Supplementary Note 6), we estimate that a strain up to ε~2.75% at A~0.09 (see Supplementary Note 3b), should be sustainable in experiments.

In summary, we have shown that deforming a 2D semiconducting monolayer with a particular topography having non-zero Gaussian curvature, can be used as a unique and straightforward way to create flat bands and drive the system into different strongly correlated electronic regimes. Topographical modulation can be created by electron-beam lithography and does not require accurate fine-tuning of the twist angle and overcomes several challenges of twisted systems. For hBN as an example, we show that bi-periodic sinusoidal modulation generates pseudo-electric and magnetic fields, creating anisotropic electronic confinement and one-dimensional flat bands. Our proposed way to create 1D flat bands should be applicable to a variety of 2D systems like hBN with massive Dirac fermions. These flat bands are different from those observed with twisted bilayer graphene and other vdW heterostructures. These bands’ one-dimensional nature will pave the route to study the exciting physics of strongly correlated 1D systems, thereby going beyond what’s achievable with twisted materials. In addition to creating flat bands, substrate engineering can be used to realize intriguing electronic behaviors such as recently demonstrated electron optics and the valley Hall effect in undulated graphene35.

For a 2D material over a gently varying topography defined as f(x,y), its structural relaxation, on a continuum level, is governed by the Föppl-von Kármán (FvK) equation22:

where subscripts denote partial derivatives. Here χ(x, y) is the Airy stress function and Y the 2D Young’s modulus. In this work we consider a sinusoidal topography f(x,y) = hsinαx.sinβy, where h defines the height, and α = 2π/Lx, β = 2π/Ly defines the lateral periodicity Lx and Ly. Plugging f(x,y) into Eq. (4) we will find

And, integrating twice, we obtain the Airy function as

Integration constants are set to zero to ensure the lowest elastic energy. With the Airy function, the components of the strain tensor can be obtained as uij = (1/Y)(εikεjl - νδikδjl)∂k∂lχ:

With the strain fields solved from above, we can easily obtain the pseudo-fields as

φp = − g(uii/2) = g((1−ν)/16)h2(β2cos2αx + α2cos2βy)

Bp = (ħ/e)(∂xAp,y − ∂yAp,x) = −(β0ħ/√2ae)((1+ν)/8)h2α2β sin2αy

The pseudo-fields derived in Eq. 2 can be obtained by setting α = β.

The relaxed 2D materials crystal geometry on curved surfaces was constructed using the displacement fields of the respective strain fields. This requires (1) determining the amount of material in the periodic box and (2) find the correct displacement for each atom. From the definition of the strain uij = (1/2)(∂iuj + ∂jui + ∂if∂jf) we can integrate the components uxx = ∂xux + (fx)2/2 and uyy = ∂yuy + (fy)2/2 and find the displacement fields ux, uy

which connects the deformed coordinates (x, y) to the reference coordinates (X, Y) = (x − ux, y − uy). The first linear term in Eq. 10 determines the overall amount of lateral contraction, and the second term is an oscillating term causing periodic stretching/compression patterns. Hence we can write (Lx-Lx0)/Lx = −(h2α2/8) and (Ly-Ly0)/Ly = −(h2β2/8), with Lx and Ly the periodicity of the sinusoid and Lx0, Ly0 the periodicity of the original, reference flake. Plugging in the above definitions α = 2π/Lx, β = 2π/Ly we have (taking i = x or y) Li2 − Li0Li + (h2π2/2) = 0

The atomic structure can therefore be constructed as the following:

Create a hBN sample with dimensions Lx0, Ly0

Choose desired height h for sinusoid

Solve for the optimal Lx, Ly, or equivalently α = 2π/Lx, β = 2π/Ly for the sinusoid

Displace each atom from (X, Y, 0) → (X + ux, Y + uy, f(X + ux, Y + uy)) with ux, uy according to Eq. 10.

The electronic structure of the pristine and sinusoidally modulated boron nitride was calculated using the density functional based tight-binding approach implemented in DFTB + using atomic orbital basis. The self consistent charge calculation was performed using a varying k-grid of 6 × 6 × 1 − 12 × 12 × 1 depending on the size of the unit cell. The maximum angular momentum chosen for B and N atoms was p (l = 1). The pairwise B-B, N-N, and B-N Slater koster files (parameterization data) were obtained from the DFTB+29 repository, matsci. The electrostatic potential was estimated by taking the Mulliken-point charges and superposing the corresponding 1/r potentials as implemented in the DFTB + code. In the code, the 1/r potential is modified to remove the r = 0 divergence, and instead plots 1/√r2 + ε2, where ε = 10−4.

A 8 band tight-binding model was developed to fit the flat bands shown by red color in Fig. 2a. The 8 bands arise from the states localized on the 2 crests and 2 troughs, and each of them being doubly occupied. The crests and troughs act as an artificial “quantum dot”. The TB model can be described by the following Hamiltonian, H = ΣxtxcϮx,ycx+2,y + ΣytycϮx,ycx,y+1, where tx and ty are hopping amplitudes along x-, and y-direction, respectively. The matrix elements of the 8 × 8 Hamiltonian for the TB model are:

H(1,1) = H(2,2) = H(3,3) = H(4,4) = H(5,5) = H(6,6) = H(7,7) = H(8,8) = 0

H(1,3) = −tx*M*e(ikxa) − tx*conj(M)*e;(−ikxa) H(3,1) = conj(H(1,3))

H(2,4) = −tx*M*e(ikxa) + tx*conj(M)*e;(−ikxa) H(4,2) = conj(H(2,4))

H(5,7) = −tx*M*e(ikxa) − tx*conj(M)*e;(−ikxa) H(7,5) = conj(H(5,7))

H(6,8) = −tx*M*e(ikxa) + tx*conj(M)*e;(−ikxa) H(8,6) = conj(H(6,8))

H(1,5) = −ty*e(ikyb) − ty*e;(−ikyb) H(5,1) = conj(H(1,5))

H(2,6) = −ty*e(ikyb) − ty*e;(−ikyb) H(6,2) = conj(H(2,6))

H(3,7) = −ty*e(ikyb) − ty*e;(−ikyb) H(7,3) = conj(H(3,7))

H(4,8) = −ty*e(ikyb) − ty*e;(−ikyb) H(8,4) = conj(H(4,8))

M = e;(−iπ/4)a = Lx/2;b = Ly/2

A 2 × 2 unit cell of hBN and 7 layers of SiO2 along [001] was taken to create the hBN|SiO2 slab geometry. The SiO2 slab’s surface was passivated with hydrogens. The geometry was fully relaxed using first-principles density functional theory (DFT) implemented in VASP36. Ion-electron interactions were represented by all-electron projector augmented wave potentials. The generalized gradient approximation (GGA) parameterized by Perdew-Burke-Ernzerhof (PBE)37 was used to account for the electronic exchange and correlation. A plane wave basis with a kinetic energy cut-off of 500 eV was used for wave functions expansion and a Monkhorst-Pack grid of 12 × 12 × 1 k-points was used to sample the Brillouin Zone (BZ). A vacuum of 20 Å was used along the direction perpendicular to the slab to reduce the interaction between the periodic images. The structure was relaxed until the Hellmann-Feynman forces on the atoms were <0.01 eV/Å. The DFT-D2 method of Grimme was used to include van der Waals interaction.

The authors declare that the data supporting the findings of this study are available within the paper and its Supplementary Information files. An archive with relaxed BN structures at different aspect ratios is deposited on the Zenodo database under accession code relaxed_structures https://doi.org/10.5281/zenodo.6523272.

Cao, Y. et al. Unconventional superconductivity in magic-angle graphene superlattices. Nature 556, 43–50 (2018).