Mechanical, electronic, optical, piezoelectric and ferroic properties of strained graphene and other strained monolayers and multilayers: an update

There are quantitative discrepancies between calculated electronic properties using DFT and TB methods without ab initio input. The π-electron TB Hamiltonian for strained graphene is given by [1]:

$$\begin{align} \boldsymbol{H} = -\sum_{\boldsymbol{r}^{^{\prime}},n} t_{\boldsymbol{r}^{^{\prime}},\boldsymbol{d}_{n}^{^{\prime}}\left(\boldsymbol{r}\right) }c_{\boldsymbol{r}^{^{\prime}}}^{\dagger} c_{\boldsymbol{r}^{^{\prime}}+\boldsymbol{d}_{n}^{^{\prime}}\left(\boldsymbol{r}\right) }+H.c., \end{align} \tag{ 4 }$$

where $\boldsymbol{r}^{^{\prime}}$ runs over all sites of the deformed honeycomb lattice and the hopping integral $t_{\boldsymbol{r}^{^{\prime}},\boldsymbol{d}_{n}^{^{\prime}}(\boldsymbol{r}) }$ varies due to changing interatomic distances. The operators $c_{\boldsymbol{r}^{^{\prime}}}^{\dagger}$ and $c_{\boldsymbol{r}^{^{\prime}}+\boldsymbol{d}_{n}^{^{\prime}}(\boldsymbol{r})}$ create an electron at $\boldsymbol{r}^{^{\prime}}$ and annihilate another at $\boldsymbol{r}^{^{\prime}}+\boldsymbol{d}_{n}^{^{\prime}}(\boldsymbol{r})$. Vectors $\boldsymbol{d}_{n}^{^{\prime}}(\boldsymbol{r})$ include multiple nearest neighbours and thus update equation (51) in [1], which included first-nearest neighbours only.

Considering up to third-nearest neighbours, the disagreement between DFT and a first-nearest neighbour TB Hamiltonian without ab initio corrections was tracked down to an angular dependence of the TB hopping parameters for second- and third-nearest-neighbours [40]. Figure 1 shows that the TB parameters up to third-nearest neighbours are modified depending on the angle of the tensile deformation. The radius of the blue and red circles in figure 1 is the deviation from the TB hopping parameter. Large corrections of the TB Hamiltonian are observed in second- and third-nearest neighbour interactions [40]; those corrections must be added to the Hamiltonian set up by Hasegawa et al [41] to describe graphene under tensile strain appropriately.

Zoom In Zoom Out Reset image size

In-plane strain leads to quadratic corrections on spin–orbit coupling [42]. Out of plane deformations, on the other hand, hybridise π-orbitals with $\sigma-$electrons and produce a first-order contribution in the spin-orbit interaction strength [43] with three contributions [42]:

$$\begin{align} \boldsymbol{H}_\mathrm{so} = \boldsymbol{H}_{A_1}+\boldsymbol{H}_{B_2}+\boldsymbol{H}_{G^{^{\prime}}}, \end{align} \tag{ 5 }$$

where the labels A₁, B₂, and Gʹ represent the irreducible representations of the group $C_6^{^{\prime\prime}}$, which results from considering a $\sqrt{3}\times \sqrt{3}$ supercell with six atoms. The non-primitive cell was employed to avoid dealing with degenerate states at two inequivalent Dirac points [42]. The corrections give rise to a Kane–Mele mass and to a Rashba-like coupling, which is present only when a mirror symmetry is broken [42]. Both coupling effects are weak (in the $1-15$ µeV range [42, 44]).

Low-energy models for non-interacting electrons in graphene with (in-plane) strain and (out-of-plane) flexural deformations [44–48] have found unusual applications, such as the holographic imaging of black holes [49], and other proposals for cosmological models [50].

Other studies relying on low-energy approximations dwell on the optical and electronic properties of graphene [51, 52], on curvature-induced quantum spin-Hall effects in Möbius strips [53], on the topological Hall effect in strained twisted graphene bilayers with enhanced interactions [54], and on the valley-dependent time evolution of coherent electron states in tilted anisotropic Dirac materials [55]. In addition, the optical properties of massive anisotropic tilted Dirac materials [56], valley polarisers and filters [57], the creation of complex magnetic fields [58] or of Landau levels in curved spaces [59], the propagation of pseudo-electromagnetic waves [60], and comparisons between twistronics and straintronics in twisted bilayers of graphene and TMDCs [61] have also appeared in the literature.

The low-energy continuum approximation takes a Dirac equation in two dimensions and a minimal coupling to effective pseudo-electromagnetic fields. The dynamics may include additional contributions caused by $\pi-\sigma$ band hybridisation (which is proportional to curvature) and/or by external electromagnetic fields. Nevertheless, this approach only works for small deformations; something most clearly seen when performing the low-energy approximation at each unit cell within the atomistic lattice [62].

Under large strain, higher-order contributions give rise to sizeable pseudomagnetic fields [63, 64]. Those corrections to the low-energy approximation can be captured through gradient expansions in real space. According to Roberts and Wiseman, corrections with nontrivial frame contribute at the same order as higher covariant derivative terms, and therefore ought to be included [64]. Comparisons against numerical diagonalisation of the empirical TB model seem to be in agreement with their observations [65]. An effective Hamiltonian for non-uniform strain taking the displacement of Dirac points into account looks as follows [66, 67]:

$$\begin{align} H_{ps} = &-i\hbar\boldsymbol{\sigma}\cdot\bar{\boldsymbol{v}}\left(\boldsymbol{r}\right)\cdot\nabla -v_{F}\boldsymbol{\sigma}\cdot\boldsymbol{A}_{s}\left(\boldsymbol{r}\right)\nonumber\\ & -\hbar v_{F}\boldsymbol{\sigma}\cdot\boldsymbol{\Gamma}_{s} \left(\boldsymbol{r}\right)+\sigma_0V\left(\boldsymbol{r}\right), \end{align} \tag{ 6 }$$

where the position-dependent Fermi velocity tensor $\bar{\boldsymbol{v}}(\boldsymbol{r})$ is given by:

$$\begin{align} \bar{\boldsymbol{v}}\left(\boldsymbol{r}\right) = v_{F}\left(\sigma_0+\bar{\boldsymbol{\epsilon}}\left(\boldsymbol{r}\right)-\beta\bar{\boldsymbol{\epsilon}}\left(\boldsymbol{r}\right)\right), \end{align} \tag{ 7 }$$

with v_F the Fermi velocity ($v_{F}/c \approx 1/300$, where c is the speed of light in vacuum), $\boldsymbol{\sigma} = (\eta \sigma_{x}, \sigma_{y})$ is the Pauli matrix vector and $\eta = \pm 1$ labels the valley, and σ₀ is the $2\times2$ identity matrix. In turn, the (deformation) scalar and (pseudomagnetic) vector potentials are given by:

$$\begin{align} V\left(\boldsymbol{r}\right) & = g\left(\varepsilon_{xx}+ \varepsilon_{yy}\right)~\text{and} \end{align} \tag{ 8 }$$

$$\begin{align} \boldsymbol{A}_s\left(\boldsymbol{r}\right) & = \left(A_{x},A_{y}\right) = \frac{\hbar \beta}{2 a}\left(\varepsilon_{xx}-\varepsilon_{yy},-2 \varepsilon_{xy}\right), \end{align} \tag{ 9 }$$

respectively. The dimensionless coefficient β ≈ 2.0–3.0 [68] measures the effect of the deformation on the first-nearest neighbour TB hopping parameter. The coupling g has been renormalised to about 4 eV [69].

The l-component of the corresponding complex vector field $\boldsymbol{\Gamma}_{s}(\boldsymbol{r})$ ($l = x,y,z$) is:

$$\begin{align} \Gamma_{s,l} = \sum_j\frac{i}{2 v_{F}}\partial_{j}{v}_{lj}\left(\boldsymbol{r}\right) = \sum_j\frac{i\left(1-\beta\right)}{2}\partial_{j}{\epsilon}_{lj}\left(\boldsymbol{r}\right). \end{align} \tag{ 10 }$$

Unlike A _s , $\boldsymbol{\Gamma}_{s}$ is purely imaginary and cannot be interpreted as a gauge field [67, 70]. Experimental signatures of the imaginary field remain elusive, but the term may become important at boundaries. (Boundaries are accounted for in effective-mass approaches of semiconductor heterojunctions by establishing consistency relationships that match atomically-resolved TB states to boundary envelope wave functions [71].) One is to remember that the low-energy approach for strain engineering is a semiclassical approach that mixes real and reciprocal space pictures, which breaks down at boundaries (periodic within neighbouring cells [62] or at edges) unless treated with care.

Indeed, a consistent treatment of boundaries within graphene strain engineering remains lacking. Nevertheless, multiple deformation fields are uniquely determined by boundaries and by sharp strain profiles: as it is the case in semiconductor heterojunctions, effective theories based on envelope wave functions call for supplemental boundary conditions of the form $\psi = M \psi$ to retain hermiticity, and for self-adjoint extensions to preserve currents [72–75]. M is a matrix containing microscopic details and symmetries, and ψ is the electron/hole wavefunction at the boundary. Here, $\boldsymbol{\Gamma}_{s}(\boldsymbol{r})$ restores hermiticity. The need for non-Hermitian corrections can also be ascertained by the breakdown of hermiticity as phasors do not add up to zero within individual unit cells under non-homogeneous strain [1, 62].

Debus, Mendoza and Herrmann matched the Hamiltonian in equation (6) and the more common version (denoted with the superindex c here) by means of the following transformations [76]:

$$\begin{align} v_{\mu} = \sqrt{g}e_{\mu}^{i} , \hspace{3mm} A_{\mu} = \sqrt{g}e_{\mu}^{i}A_{i}^{c}, \hspace{3mm} \Gamma = \sqrt{g}e_{\mu}^{i}\Gamma_{i}^{c}, \end{align} \tag{ 11 }$$

where

$$\begin{align} e_{\mu}^{i} = \frac{g_{ai}+\delta_{ai}\sqrt{g}}{\sqrt{Tr\left(g\right)+2\sqrt{g}}}, \end{align} \tag{ 12 }$$

$g_{\mu\nu}$ is the metric tensor associated with strain, δ_ai is the Kronecker delta, $\sqrt{g}$ the square root of the determinant of the metric tensor, and $Tr(g) = \sum_{\mu} g_{\mu\mu}$ is its trace.

The position-dependent Fermi velocity tensor components obtained from equation (6) were validated on graphene under uniaxial strain ε_xx , and used to tune the sensitivity factor of a graphene pressure sensor $GF = ((\Delta R)/R)/\varepsilon_{xx}$, where R was the resistance and $\Delta R$ the change of R due to strain [77]. Modifications of optical properties under uniaxial strain are detailed elsewhere [78, 79].

On samples with large out-of-plane deformations, $\sigma-\pi$ orbital hybridisation leads to corrections [47, 80]:

$$\begin{align} V_{\pi\sigma}\left(\boldsymbol{r}\right) = -g_1\left(\nabla ^{2} h\right)^{2} \end{align} \tag{ 13 }$$

that must be added to the scalar potential $V(\boldsymbol{r})$. Here, $g_1 = \alpha (3/4) a^{2}$ and α ≈ 9.23 eV.

Substrate deformations can also make the on-site energies vary. Such deformation potential [68] can be treated by decomposing the interaction into a smooth spatial effective potential $V_\mathrm{sub}(\boldsymbol{r})\sigma_0$ [47] and, if the substrate produces a bipartite symmetry breaking, an extra term $\Delta V(\boldsymbol{r})\sigma_z = (V_A(\boldsymbol{r})-V_B(\boldsymbol{r}))\sigma_z$ [44], which measures the difference of the electrostatic potential in the A and B sublattices [81].

The tearing and subsequent folding of the topmost graphitic layer by an STM tip is depicted in figure 2 [82]. There, the tip approaches graphite right before a step edge, moving across it afterwards (figure 2(a)) [82]. As a result, the tip (i) tears and (ii) folds the uppermost monolayer (figures 2(b)–(d)), creating nonuniform strain at the folded edge. The effects of the resulting effective gauge potential were seen by $\mathrm{d}I/\mathrm{d}V$ measurements taken at locations marked by circles in figures 2(e) and (f) [85]. The splitting of $\mathrm{d}I/\mathrm{d}V$ peaks (and their weight redistribution) in figure 2(g) are manifestations of localised states on folded graphene [82, 83].

Zoom In Zoom Out Reset image size

One-dimensional electronic bands were created along the longest fold in figure 2(d). There, the axial direction of the fold turned out to be parallel to graphene's armchair direction, and the experimental electronic density [82] was contrasted against a linear chain TB Hamiltonian model, where the strain appears as a modulation of hopping parameters and self-energies [1]. According to the model, the flat bands are topological in nature [86]. Folds also create curvature [87] that pushes π-electrons onto the outer side of the fold [80]. The creation of valley polarisation through folds will be discussed in subsequent pages.

Wu et al [84] (figure 3) revealed a Coulomb blockade, so that their fold acted as a quantum dot. Linear folds as the one in figure 3(a) [84, 88–90] were theoretically studied using continuous models with effective pseudomagnetic fields and TB models. Those folds feature bands originating from a mass term in the Dirac equation [91]. Calculations using the Dirac equation in a continuum and including strain were consistent with measured bound state energies, and predicted valley-polarised currents [84], opening an avenue for graphene folds as straintronic quantum wires [92, 93].

Zoom In Zoom Out Reset image size

Triangular periodic ripples in suspended graphene were found in another work. There, bending was limited to a region few unit cells wide [94]. Triangular ripples are unlike the sinusoidal ripples found in fabrics.

The effects of folds in the electronic structure are modeled with effective pseudomagnetic fields within a continuum picture [95]: one considers a pure deformation perpendicular to graphene's plane, i.e. $\boldsymbol{u}(\boldsymbol{r}) = 0$. From equations (3) and (9), one gets:

$$\begin{align} V\left(\boldsymbol{r}\right) & = \frac{g}{2}|\nabla h|^{2},~\text{and} \end{align} \tag{ 14 }$$

$$\begin{align} \boldsymbol{A}_s\left(\boldsymbol{r}\right) & = \left(A_{x},A_{y}\right) = \frac{\hbar \beta}{4 a}\left(\partial_{x}h\partial_{x}h-\partial_{y}h\partial_{y}h,-2\partial_{x}h\partial_{y}h\right). \end{align} \tag{ 15 }$$

By the definition of uniaxial fold, the deformation field is constant in one direction, say ($h(x,y) = h(y)$). If periodic, it is expressed by a Fourier expansion:

$$\begin{align} h\left(y\right) = \sum_{k = -k_c}^{k_c} a_k e^{iky}, \end{align} \tag{ 16 }$$

in which the Fourier coefficients obey $a_{-{k}} = a_{{k}}^{*}$ (because $h({y})$ is real), and k_c is a cutoff parameter that can be estimated from experiment [95]. This way, $\mathbf{A}_s(\mathbf{r})$ has only one non-zero component [95]:

$$\begin{align} A_x\left(y\right) = \frac{\hbar \beta}{4 a} \left[ \sum_{k} a_k ke^{iky} \right]^{2}. \end{align} \tag{ 17 }$$

Equation (17) implies that $\nabla \cdot \mathbf{A}_s(\mathbf{r}) = 0$, so that $\mathbf{A}_s(\mathbf{r})$ is in the Coulomb gauge. Thus it can be derived from a pseudoscalar magnetic potential:

$$\begin{align} A_i = \epsilon_{ij}\partial_j \Phi\left(\boldsymbol{r}\right), \end{align} \tag{ 18 }$$

with $i,j = x,y$ and _ij a 2D Levi-Civita tensor. Expressing $\Phi(y)$ as a Fourier series:

$$\begin{align} \Phi\left(y\right) = \Phi_0\left(y\right)+ \sum_{k \ne 0}e^{iky} \tilde{\Phi}\left(k\right), \end{align} \tag{ 19 }$$

with

$$\begin{align} \Phi_0\left(y\right) = \frac{\hbar \beta}{4 a}\left(\sum_k k^{2}|a_k|^{2}\right)y, \end{align} \tag{ 20 }$$

and

$$\begin{align} \tilde{\Phi}\left(k\right) = -i\frac{\hbar \beta}{4 a k} \left[ \sum_{k^{^{\prime}}}a_{k}a_{k^{^{\prime}}-k}^{*}k^{^{\prime}}\left(k^{^{\prime}}-k\right) \right]. \end{align} \tag{ 21 }$$

The associated pseudomagnetic field becomes:

$$\begin{align} \mathbf{B}_s = \mathbf{\nabla}^{2} \Phi \left(\mathbf{r}\right). \end{align} \tag{ 22 }$$

$\Phi_0(y)$ does not produce a pseudomagnetic field, but instead an Aharonov–Bohm effect manifested by a phase change along closed circuits enclosing folds. As $\mathbf{A}_s(\mathbf{r})$ is derived from $\Phi(\mathbf{r})$, zero-modes can always be found for any fold with an arbitrary cross-section (periodic, random, quasi-periodic, etc assuming a large unit cell for the periodic function h(y)), and being topologically protected, their existence is independent of $V(\mathbf{r})$. Using equation (6) for E = 0, two wave functions are obtained:

$$\begin{align} \psi_{\pm}\left(\mathbf{r}\right) = \left(\mathrm{const.}\right) \left(1\pm \sigma_z\right) \left( \begin{array}{c} e^{\Phi\left(y\right)} \\ e^{-\Phi\left(y\right)} \end{array} \right), \end{align} \tag{ 23 }$$

where σ_z is the Pauli z-matrix. Those topological modes appear as a consequence of the Atiyah-Singer theorem, and correspond to soliton solutions of the Dirac equation with a space-dependent mass [96].

This approach permits exploiting similarities with the integer quantum Hall under a random magnetic field [97]. In such field, the vector potential is given by a gaussian distribution with zero mean and variance $\Delta_{A}$. The average of the coefficients in equation (19) can be expressed as:

$$\begin{align} \langle \tilde{\Phi}\left(k\right)\tilde{\Phi}\left(k^{^{\prime}}\right)\rangle = \left(2\pi\right)^{2}\delta\left(k-k^{^{\prime}}\right)\frac{\Delta_A}{k^{2}}. \end{align} \tag{ 24 }$$

From this analogy with real magnetic fields it can be concluded that the wave-function is multifractal [97]. The multifractal localisation can be determined from the participation ratio moments $P_q(L)$ where L is the sample length [98]:

$$\begin{align} P_q\left(L\right) = \langle |\psi\left(\mathbf{r}\right)|^{2q}\rangle, \end{align} \tag{ 25 }$$

from where it follows that [97]:

$$\begin{align} P_q\left(L\right) \approx \frac{1}{L^{2+\tau\left(q\right)}}, \end{align} \tag{ 26 }$$

with

$$\begin{align} \tau\left(q\right) = 2\left(q-1\right)+\frac{\Delta_A}{\pi}q\left(1-q\right). \end{align} \tag{ 27 }$$

Such multifractal conductance fluctuations have been measured in high-mobility mesoscopic graphene devices under external magnetic fields [99]. In addition, edge modes arising in graphene under periodical strain-induced pseudomagnetic fields have been contrasted to those arising under real, periodic magnetic fields [86, 100].

It is now time to address the creation of sublattice charge imbalance, which can be visualised in local density of states (LDOS) calculations and $dI/dV$ STM measurements [62, 91, 101–103]. Strain-induced sublattice imbalances create an energy separation of the two pseudospin degrees of freedom due to the pseudomagnetic field B_s . The creation of charge imbalance can be seen by squaring graphene's Dirac Hamiltonian with the strain-induced pseudo-magnetic field B_s :

$$\begin{align} E^2\Psi^A& = v_F^2\left[\boldsymbol{\pi}^2+e\hbar B_{s}\right]\Psi^A, \nonumber\\[8pt] E^2\Psi^B& = v_F^2\left[\boldsymbol{\pi}^2-e\hbar B_{s}\right]\Psi^B, \end{align} \tag{ 28 }$$

where E is the electron's energy, v_F the Fermi velocity, $\Psi^{A}$ ($\Psi^{B}$) the wave function amplitude on the A (B) sublattice, and π the canonical momentum. The term proportional to B_s appears with opposite signs at each sublattice, shifting their energy in opposite directions and leading to a (Zeeman-like) pseudospin (sublattice) polarisation.

An STM setup was designed to induce and visualise sublattice imbalance [101]. Similar to [104], the tip lifted graphene from the substrate, creating a gaussian-shaped deformation that moved along with the tip (see figures 4(a) and (b)) and a pseudomagnetic field with 3-fold symmetry (figure 4(c)). $\mathrm{d}I/\mathrm{d}V$ measurements (figures 4(d)–(g)) performed by the same STM tip show increasing sublattice LDOS imbalance with increasing tip current. Analytical expressions for sublattice contrast [103] and molecular dynamic simulations showed good agreement with experiment. Gaussian deformations have been proposed as valley filters [84, 105–110], and as Aharonov–Bohm interferometers [111]. Strain fields similar to those created by gaussian-like deformations were observed in graphene nanobubbles [34].

Zoom In Zoom Out Reset image size

Producing valley-polarised states is the aim of valleytronics [113, 114]. However, there are important differences between graphene's valley, sublattice, and spin degrees of freedom: spins can be controlled by magnetic fields directly, and the sublattice can be polarised directly by strain as well. Unfortunately, controlling the valley polarisation of charge carriers in graphene is not that straightforward [46]. For example, while a real magnetic field B_R splits spin polarised states, a pseudo-magnetic field B_s does not split valley polarised states [46].

Some types of strain can be employed to control the valley pseudospin of graphene's carriers [112, 115–118]. Depicted in figure 5, one method involves the combination of strain-induced pseudo-magnetic fields B _s and real magnetic fields B _R [116, 117]. The approach works because B _s points in opposite directions at each valley to satisfy time-inversion symmetry, while the direction of B _R is the same for all carriers. This way, magnetic and pseudo-magnetic fields add up in one valley, and suppress one another at the other valley, leading to different energy locations of Landau levels for each valley $K_\pm$:

$$\begin{align} E_{n\pm}\sim\sqrt{n|B_R\pm B_{s}|}. \end{align} \tag{ 29 }$$

Thus, energy levels of the $K_+$ valley are pushed up, while those of the $K_-$ valley are pushed down.

Zoom In Zoom Out Reset image size

This effect was observed in rippled graphene grown on rhodium foil [112, 119–121]: atomically-resolved STM images along the ripple reveal compressive strains along graphene's zigzag direction (figure 5(a)). A TB Hamiltonian accounting for deformations along the ripple was used to approximate the resulting pseudo-magnetic field B _s , which was non-uniform along the ridge of the ripple (dropping to zero and reversing its sign) as shown in figure 5(b).

Spectra recorded at different positions along the ripple for $\boldsymbol{B}_R = 0$ T (figure 5(c)) show pseudo-magnetic Landau levels only at points where $\boldsymbol{B}_{s}\neq 0$. The ($n\neq0)$ pseudo-magnetic Landau levels further split in two when a magnetic field is introduced. This splitting is only seen at points along the ripple where both B _s and B _R are nonzero, in agreement with equation (29).

The valley energy levels invert their energy as the pseudomagnetic field ${B}_{s}$ changes sign in equation (29): beyond their valley polarisation, charge carriers also exhibit valley inversion along the ripple because B _s reverses sign along the y direction (figure 5(b)). This is shown at the top panel of figure 5(d): the energy levels of valley $K_+$ ($n_+ = 1,2,3$) lie lower in energy than those of the $K_-$ valley ($n_-$) when B _s is positive. B _s is negative at the lower panel of figure 5(d) and the levels split again, but now the order is inverted and the LLs of the $K_+$ valley appear at higher energy than those of the $K_-$ valley. The non-strained region where the B _s reverses its direction is not necessary to create valley polarisation, but it is critical for valley inversion.

There is a compromise between the magnitude and the spatial extent of magnetic fields at the nanoscale: fields created by large magnets tend to be homogeneous across the samples' dimensions, while inhomogeneous fields produced by micromagnets or by superconducting gates are rather weak. However, strain makes it possible to produce pseudogauge fields $\mathbf{B}_s(\mathbf{r})$ that are strong and vary at nanometer scales.

For instance, strong and inhomogeneous pseudo-magnetic fields were created in graphene deposited on Cu foils in [122] (figure 6): those foils exhibit flat terraces separated by vertical steps up to ${\sim} 35$ nm in height. Graphene drapes over the steps when it is grown, becoming pulled by contact forces at the terraces. This gives rise to a periodically modulated strain profile in the (vertical) draped region, and to a graphene superlattice with pseudomagnetic fields of alternating signs modelled in figure 7(a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The nanometer-scale modulation of the pseudomagnetic field affects the electronic spectrum in radical ways. Indeed, the spectrum does not exhibit Landau quantisation of Dirac fermions of the form $E_n\propto \sqrt{n}$, but to a quantisation appearing linear in energy ($E_n\propto n$). Although a similar type of quantisation can occur due to confinement, such explanation is discarded because the gap between peaks is independent of the ripples' wavelength [122]. Instead, calculations with a TB model indicate that this quantisation is a consequence of a pseudomagnetic field oscillating with a period L_B of the same order of magnitude than the magnetic length ($l_B = \sqrt{\hbar/eB_s}$).

To reproduce the LDOS spectra, Banerjee et al considered both the pseudomagnetic and deformation potential fields created by the periodic strain [122]. Figures 7(c) and (d) contrast the calculated LDOS from a TB Hamiltonian considering the strained graphene system, and those based on a Dirac equation with a periodic pseudomagnetic field and electric potential of the form:

$$\begin{align} \boldsymbol{B}_s = B_\mathrm{max}\sin\left(4\pi y/\lambda\right)\mathbf{e}_z, \end{align} \tag{ 30 }$$

$$\begin{align} V = V_\mathrm{max}\cos\left(4\pi y/\lambda\right), \end{align} \tag{ 31 }$$

respectively. Both types of calculations showed good agreement with experimental measurements [122]. Strain patterns with alternating zones of up/down pseudomagnetic fields can create valley-Hall edge states, which could be observable at room temperature due to the pseudo-magnetic fields ($B_s\sim100$ T) that can be experimentally created: graphene ripples might provide a platform for valley-dependent electron transport.

Methods to induce pre-designed strains are highly desirable, and a promising avenue is to support graphene on patterned substrates [28, 123–126] which create a network of wrinkles (figure 8). Depending on the height-to-separation ratio of those pillars, sections of graphene might rest on the substrate, or the entire monolayer may be suspended like a canopy.

Zoom In Zoom Out Reset image size

Estimations of strain with a STM are usually derived from $\mathrm{d}I/\mathrm{d}V$ spectra. Nevertheless, the magnifying effect of moiré patterns allowed to characterise the strain in graphene directly in [28]. The setup consisted of a gold nanopillar array over a hBN flake, which produces a moiré due to the mismatch with graphene (see figure 9). The pillar-induced strain on graphene modifies the moiré [127], making it possible to quantify the strain in graphene (whose maximum value turned out to be about $4.5\%$) with standard STM resolution. The measured strain was in agreement with estimations made from the energy separation of pseudo-Landau levels from the $\mathrm{d}I/\mathrm{d}V$ spectra. Such levels of strain would have been difficult to detect directly without the additional moiré pattern, which allowed for a 20-fold increase in atomistic precision.

Zoom In Zoom Out Reset image size

Similar rectangular nanostructure arrays were used to strain graphene in [126]. There, wrinkles–quasi-1D 'channels' with a near-uniform pseudo magnetic field–developed on graphene when it was placed between two nanostructures having a critical separation. Calculations revealed that parallel arrays of those channels can be used to realise valley polarisation. Other works have studied strained graphene over SiO₂ nanospheres [124], or over black phosphorus substrates [128].

As discussed in section 3.1.5, a non-zero sublattice imbalance (this is, an in-plane electric polarisation P) can be created by strain that breaks the inversion symmetry of the graphene lattice. Sublattice charge imbalance has also been verified through optical second-harmonics: the schematics on figure 10(a) show graphene on nanopillars again (see figures 8 and 9), and this sample is irradiated with a pump laser at 1035 nm. Figure 10(b) is a measure of strain through Raman spectra: according to [129], the sample is non-uniformly strained up to ∼2.4% by the nanopillar array. The temperature-dependent peak at half wavelength (1035/2 nm) displayed in figure 10(c) confirms the creation of a second harmonic on non-uniformly-strained graphene. That a pseudomagnetic field can be obtained from within atomistic displacements was stated in [102], and the fact that the polarisation differs in between lattices under anisotropic strain can be traced down to the fact that the effect on the diagonal entries of the Dirac Hamiltonian can be written as [81]:

$$\begin{align} \begin{pmatrix} E_{s,A} & 0\\ 0 & E_{s,B} \end{pmatrix} = \bar{E}\sigma_0+\Delta_E\sigma_z, \end{align} \tag{ 32 }$$

where $\bar{E} = (E_{s,A}+E_{s,B})/2$ and $\Delta_E = (E_{s,A}-E_{s,B})/2$.

Zoom In Zoom Out Reset image size

Second-harmonic generation has also been experimentally employed to sense strain in MoS₂ monolayers [130], and it is a powerful and promising technique suitable for strain characterisation of arbitrary 2D materials. More will be said on structures lacking a centre of inversion, on the existence or creation of electric dipoles, and on second harmonic generation in forthcoming pages.

The absence of superconductivity (SC) in intrinsic graphene (with its Fermi energy located at the zero-density Dirac point) is expected from the Bardeen–Cooper–Schrieffer (BCS) theory: in the limit of small DOS at the Fermi level $\rho(E_F)$, that theory indicates a critical temperature depending exponentially on the coupling strength λ, $T_c\sim e^{-1/(|\lambda|\rho(E_F))}$ [131], and a large charge doping would be required to produce a T_c of a few Kelvin.

Because flat bands enhance electronic interactions, strained graphene can be of interest for studying many-body phases arising due to a flat electronic spectrum. Flattening of the electronic bands provides an alternative to increase the DOS and hence to induce SC, because the BCS theory leads to a linear dependence between T_c and the coupling strength λ, $T_c\sim \lambda$ [132–134] for large $\rho(E_F)$.

Now, although large magnetic fields can induce flat bands (Landau levels), they also break time reversal symmetry and suppress (singlet) SC order. This is due to the Zeeman effect, which causes alignment of the electron spins and breaks Cooper pairs. Strain-induced flat bands, on the other hand, do not break time reversal symmetry and open the possibility to engineer SC [133–139]. Furthermore, one can choose the strain field's symmetry, period, and strength separately. (This situation is unlike that observed in twisted bilayer graphene (TBG), where the corresponding parameters are mostly determined by the moiré pattern. Note that flat bands can also be created in few-layer rhombohedral graphene [140].)

Strain-induced flat bands in graphene are due to wavefunctions localised at domain walls where the strain-induced potential changes sign [134, 141]. This permits mapping the system into a continuum Jackiw–Rebbi model with pseudo-Landau levels [141]. In fact, an exact Su–Schrieffer–Hegger (SSH) model is obtained when the momentum is set parallel to the direction of the strain modulation [134]. Interestingly, the inclusion of the electron–electron interaction leads to magnetic domains and a pairing effect [141]. A model for superconductivity in graphene with strain-induced flat bands was studied in [139], and it is discussed next. Strain was added to the low-energy dispersion of graphene via periodic (1D or 2D cosine) pseudo-vector potentials A _s of period d and strain strength u₀ (see equation (9)):

$$\begin{align} \boldsymbol{A}_{\cos}^{1D}\left(x,y\right)& = \frac{\hbar u_0}{d}\left(0,\cos \frac{2\pi x}{d}\right),~{\text{or}} \end{align} \tag{ 33 }$$

$$\begin{align} \boldsymbol{A}_{\cos}^{2D}\left(x,y\right)& = \frac{\hbar u_0}{d}\left(\cos \frac{2\pi y}{d},\cos \frac{2\pi x}{d}\right), \end{align} \tag{ 34 }$$

which can be realised by the following in-plane displacement fields:

$$\begin{align} \boldsymbol{u}_{\cos}^{1D}\left(x,y\right)& =- \frac{u_0 a}{\beta \pi} \left(0,\sin\frac{2\pi x}{d}\right),~{\text{or}} \end{align} \tag{ 35 }$$

$$\begin{align} \boldsymbol{u}_{\cos}^{2D}\left(x,y\right)& =- \frac{u_0 a}{\beta \pi} \left(0,\sin\frac{2\pi y}{d}+\sin\frac{2\pi x}{d}\right), \end{align} \tag{ 36 }$$

(see figures 11(a) and (b)) with β the Grüneisen parameter for graphene. The vector potentials can be obtained from these displacement fields by following equations (3) and (9), and the notation in [139] was modified for consistency with the one used here. The resulting pseudomagnetic field $B_{s} = |\nabla \times \mathbf{A}_s|$ is shown in figures 11(c) and (d). (The pseudomagnetic field strength ($B_{s} \approx$ 100 T) and period ($d\approx$ 14 nm) visualised in experiments [142] would lead to  u₀≈ 5 for the displacement fields defined here.)

Zoom In Zoom Out Reset image size

The SC state was modeled by the Bogoliubov–de Gennes mean-field theory [143] once strain was introduced, assuming an intervalley local attractive interaction leading to a s-wave pairing Δ order parameter (figures 11(e) and (f)). Critical values are studied as a function of the interaction strength. The low-energy states of strained graphene exhibit a sublattice polarisation that follows the pseudo-magnetic field and, accordingly, the SC order parameter $\Delta_{A/B}$ is sublattice-dependent. $\Delta_{A/B}$ always peaks at the largest value of the pseudomagnetic field. This situation can be contrasted against TBG, where Δ is localised around AA stacking regions for the first-magic angle.

The electronic energy dispersion of graphene with strain fields and without them is shown in figure 12. The case with $\boldsymbol{A}_{\cos}^{1D}$ is shown both in the mixed and reduced zone schemes for easier comparison. Substantial band flattening occurs at low energies, similar to those observed in TBG.

Zoom In Zoom Out Reset image size

The SC order parameter Δ exceeds the flat-band bandwidth for high enough values of the strain strength u₀ and the interaction strength λ, and its dependence is linear in λ [132–134]. At lower values of u₀ and λ, the dispersive behaviour of the bands becomes important and the order parameter is exponentially suppressed.

While the buckled graphene studied by Mao et al [142] is a promising platform for correlated phenomena, it is not in the flat band regime (where the linear dependence $\Delta\sim\lambda$ occurs) and it would require increasing the strain strength u₀ by a factor of 4 in order to lead to a SC transition at a temperature on the order of 1 K. Moreover, for this strain strength, Peltonen and Heikkilä [139] estimated a transition temperature of ≈11 K for a period d ≈ 8 nm. The calculation is done using the same λ estimated for TBG [144]. SC could still be possible for lower strain fields outside of the flat band regime, but a T_c well under 1 K would be expected. Regardless, the possibility of having tunable SC and order parameters via strain is attractive, and the generality of the model might be of relevance in the study of SC in other graphene systems like TBG–where intrinsic strain is always present (more on this later).

Correlations were introduced in a similar model of periodically strained graphene via a spatially-dependent antiferromagnetic coupling [138], showing that a chiral, d-wave SC can be stabilised via strain as well.

It was just argued that flat bands are a necessary ingredient to create superconductivity in graphene. Additional approaches to achieve flat electron bands are arising because of that.

Stiff membranes can undergo a buckling transition when exceeding a critical value of in-plane compressive strain (because out-of-plane distortions reduce the elastic energy of the membrane). Those distortions lead to periodic strain patterns, which in turn create periodic pseudomagnetic fields that (as seen in the previous section) can give rise to flat electronic bands [142].

Triangular buckling patterns consisting of alternating crests and troughs have been observed in graphene deposited over NbSe₂ [142]. Those superstructures had a height modulation of 0.17 nm and a period of about $a_b = 14.4$ nm. A schematic of buckled graphene is shown in figure 13(a). In the crest regions, the $\mathrm{d}I/\mathrm{d}V$ spectra exhibited peaks that were identified as the pseudo-Landau levels induced by a pseudomagnetic field B_s ($E_n\sim \text{sgn}(n)\sqrt{|n|B_{s}}$ for $n = 0,\pm1,\pm2,\ldots$), confirming the existence of a pseudo-magnetic field with a magnitude of about $B_{s} = 108$ T. Furthermore, sublattice polarisation, a hallmark of the wavefunction of the n = 0 pseudo-Landau level, was visualised through STM measurements at crests and troughs [142]. The sublattice polarisation (A or B) seen in the STM measurements changes when moving from a crest to a trough, indicating a change in the sign of B_s between these two regions (there was no sublattice polarisation at the interface between crests and troughs, indicating that $B_{s} = 0$ there).

Zoom In Zoom Out Reset image size

The $\mathrm{d}I/\mathrm{d}V$ spectra and the pseudo-Landau levels were described by a TB model of graphene in the presence of a periodic pseudo-magnetic field with triangular symmetry [142]:

$$\begin{align} B_{s}\left(x,y\right) = B\left[\cos\left(\boldsymbol{b}_1\cdot\boldsymbol{r}\right)+\cos\left(\boldsymbol{b}_2\cdot\boldsymbol{r}\right)+\cos\left(\boldsymbol{b}_3\cdot\boldsymbol{r}\right)\right], \end{align} \tag{ 37 }$$

where B is the amplitude of the pseudo-magnetic field, a_b the superlattice period, $\boldsymbol{b}_1 = \frac{2\pi}{a_b}(1,-\frac{1}{\sqrt{3}})$, $\boldsymbol{b}_2 = \frac{2\pi}{a_b}(0,\frac{2}{\sqrt{3}})$ and $\boldsymbol{b}_3 = \boldsymbol{b}_1+\boldsymbol{b}_2$. The choice of coordinates is consistent with C₆ symmetry and with the fact that the total pseudo-magnetic flux has to vanish (since time-reversal symmetry is not broken by strain). When the sample is doped to partially fill the n = 0 pseudo-Landau level, a pseudo-gap feature splits the peak in two, indicating the appearance of a correlation-induced state. Moreover, a spectral weight redistribution between the two peaks is observed when doping away from charge neutrality. Similar signatures have been observed upon partially filling the flat band of magic-angle TBG.

Milovanović et al studied the flat bands in periodically buckled graphene [146] considering three different geometries: a triangular pseudo-magnetic field mode (experimentally observed in [142]), a hexagonal buckling mode (which leads to a pseudo-magnetic field similar to the ones obtained for gaussian bumps and bubbles), and a herringbone buckling mode corresponding to an undulating 1D strain pattern. The herringbone buckling has the lowest energy of those three configurations at large strain.

The amplitudes and periods of pseudo-magnetic fields required to access correlated phases were ascertained by comparing the bandwidth $\Delta E$ to the characteristic energy scale for the interactions, E_I . The critical pseudo-magnetic field takes this form when employing the triangular pseudo-magnetic field mode [146]:

$$\begin{align} B^\mathrm{int}_n = c_n \Phi_0/S, \end{align} \tag{ 38 }$$

at which $\Delta E = E_I$ for the nth band. Φ₀ is the magnetic flux quantum and S is the area of the unit cell. In the triangular pseudo-magnetic field mode, $c_1 = 6.5$ for the lowest band. For the next band $c_2 = 4.4$. Because higher bands have smaller bandwidths and are also gapped, they might also be relevant for studies of correlated phenomena.

Bands are not as flat in the hexagonal buckling mode: broader LDOS peaks indicate that this pseudo-magnetic field cannot localise electrons efficiently. Similarly, the herringbone mode does not exhibit a gap opening nor flat bands, even for strains as large as 20$\%$. This is a consequence of this buckling mode being the lowest-energy configuration: the average strain in the herringbone mode is generally half of that displayed in the hexagonal mode (and its pseudo-magnetic field is even an order of magnitude smaller).

Buckled graphene superlattices do not localise charge carriers fully. Spatial regions in between those undergoing large pseudo-magnetic fields turn into percolating paths, along which charge carriers propagate. This is seen in the longitudinal conductance, which shows peaks coinciding with the LDOS [146] and confirms the non-localised nature of flat band states.

Electron–electron interactions were introduced by a Hubbard term in a TB model for buckled graphene [147]:

$$\begin{align} H = -t\sum_s\sum_{\langle i,j \rangle}c_{is}^\dagger c_{js}+U\sum_i c_{i\uparrow}^\dagger c_{i \uparrow} c_{i \downarrow}^\dagger c_{i \downarrow}, \end{align} \tag{ 39 }$$

while strain was introduced by modifying the hopping energies connecting the three nearest neighbours $t_n = t+ \delta t_n$ ($n = 1,2,3$) with

$$\begin{align} \delta t_n = - \frac{\sqrt{3}e v_F L_M}{4 \pi}\sin\left(\boldsymbol{b}_n\cdot\boldsymbol{r}\right), \end{align} \tag{ 40 }$$

and the model was solved with the non-collinear mean field formalism for U = t [147].

They found correlated states at values close to $L_M/l_m\approx 6$. While the system remained gapless in the non-interacting case, an emergent magnetic state with a non-homogeneous magnetisation and a correlation gap of $\Delta\approx0.01t$ (consistent with experiment [142]) was obtained when interactions were turned on. The correlated state is not expected to arise for pristine graphene at the same value of U. The strain-induced band flattening is essential for such magnetism to emerge, because it creates a bandwidth $\Delta E$ smaller than the Hubbard interaction U ($\Delta E/U \sim 0.2$). The periodic magnetisation patterns give rise to an emergent honeycomb superlattice with antiferromagnetic ordering. The regions of highest LDOS can be understood as Wannier orbitals, and the sign of the net magnetisation at a given region follows the sign of the pseudo-magnetic field. The magnetic ordering decays as the system is doped away from half filling, a result consistent with experiment [142].

Because of the strain-induced crests and troughs in the superlattice, an electric field impinging along the z-direction can induce non-homogeneous local energy shifts, which are modeled via height-dependent onsite energies $\mu(\boldsymbol{r}) = \mu_0\sum_{i = 1}^3\cos(\boldsymbol{b}_i\cdot\boldsymbol{r})$ in the TB Hamiltonian. This term creates sublattice imbalance in the emergent honeycomb superlattice, and competes with the magnetic ordering, thus holding potential for electrically tuning the magnetic ground state, and for exploring the interplay between strain-induced gauge-fields and electron-electron interactions.

Manesco and Lado [145] derived an effective model for the emergent honeycomb lattice that arises in buckled graphene (figures 13(a) and (b)). The minima and maxima of $B_\mathrm{eff}(\boldsymbol{r})$ contain the highest LDOS, and were considered as Wannier orbitals. The effective TB model for the emergent honeycomb lattice includes sublattice imbalance and a valley-dependent, second-nearest neighbour hopping. Remarkably, such model turns out to be fully analogous to the Kane–Mele model for a spin quantum Hall insulator [148, 149]. Here, the valley isospin plays the role of spin, and the superlattice's mini-valleys play the role of the valleys. Due to the height variation in the buckled graphene, out-of-plane displacement fields increase (or decrease) the energy of states at the $B_\mathrm{eff}(\boldsymbol{r})$ maxima (minima). In the effective model, this is equivalent to increasing the sublattice imbalance, which in turn acts as a knob to tune the system's topology.

Onsite and nearest neighbour Hubbard terms were part of the effective model (with interaction strengths U and V, respectively), and a mean field approach was used. Two different ground states (a charge density wave and an antiferromagnetic state) were found at different values of U and V [145]. Each ground state creates different sublattice imbalances and, as in the Kane–Mele model, this leads to topological transitions.

Indeed, the valley Chern number was calculated for different values of the interaction-induced sublattice imbalance, and both the charge density wave and antiferromagnetic phases in the interacting phase diagram can be further divided into topologically trivial and nontrivial phases. Analogous to the Kane–Mele model, the non-trivial phases are quantum valley Hall insulators that exhibit counter-propagating edge states with opposite valley polarisation. Edge states are further spin-polarised in the topologically non trivial antiferromagnetic phase.

A closely related phenomena is the curvature-induced spin-Hall effect in a graphene Möbius strip. The solution of the Dirac equation shows that despite the absence of a Hall current, a spin-Hall current is a natural consequence for such a topology [53]. Other works have studied topological phases arising from strain and its interplay with interactions in graphene [100, 116, 136, 150–152]. Buckling was also shown to allow control of topological edge states via electric fields in germanene [153].

The results within this subsection place buckled graphene at the interface of strain, flat bands, correlations, and valley topology.

Massive and massless charge carriers were observed on caesium-doped graphene; a system modelled as a strained quantum well [154]. Extended flat bands were reported for this system in a later work [155].

Adatoms within substrates, localised mechanical probes, and electron-electron interactions induce short wavelength lattice modulations on graphene [156]. This subsection is focused on the paradigmatic case of Kekulé patterning observed on graphene deposited over a Cu(111) single crystal (and dubbed Kek-Y) [157] and on Li-intercalated graphene (Kek-O) [158]. As depicted in figure 14 by bonds of two different strengths, Kekulé patterns have few hopping strengths changed within a $\sqrt{3}\times \sqrt{3}$ periodic supercell (shown by dashed lines within the figure). Two Kek-Y structures–related by a reflection with respect to the vertical line–are shown in figure 14 as well [159].

Zoom In Zoom Out Reset image size

A study on strongly-interacting electrons in hexagonal lattices [160] will be highlighted now. Figure 15 presents graphene (a semimetal (SEM): this is, an electronic system with interactions turned off), an antiferromagnetic insulator (AFI), the Kek-O superlattice (renamed as 'dimerized insulator' (DIM) in [160]), and a lattice with a hexagonal distortion (dubbed HEX). Those will turn out to be competing phases on the honeycomb lattice.

Zoom In Zoom Out Reset image size

One sets graphene onto a Hubbard model Hamiltonian with an additional elastic energy term [161]:

$$\begin{align} H\left(\left\{t_{xy}\right\}\right) & = -\sum_{\sigma = \uparrow,\downarrow} \sum_{\langle x,y \rangle} t_{xy}\left(c^{\dagger}_{x,\sigma}c_{y,\sigma}+h.c\right) \notag \\ & \quad + U\sum_{x} \left( n_{x,\uparrow}-\frac{1}{2}\right) \left(n_{x,\downarrow}-\frac{1}{2}\right)+\sum_{\langle x,y \rangle} F\left(t_{xy}\right). \end{align} \tag{ 41 }$$

The $c_{x,\sigma}$ on equation (41) are fermion annihilation operators at lattice site x with spin σ. The occupation number of such site is $n_{x} = c^{\dagger}_{x,\sigma}c_{x,\sigma}$. The on-site repulsion U can have any sign as long as it is the same for all lattice sites. The first nearest-neighbour hopping integrals $t_{x,y}$ between sites x and y are real, as there is no magnetic field, and are assumed positive although not necessarily independent of the $x,y$ pair (x and y here do not denote Cartesian components but different sites within the Honeycomb lattice). The TB parameters t_xy depend on the physical distance r_xy between the lattice sites x and y, and are usually assumed to decay exponentially [1].

The main addition in equation (41) with respect to equation (39) is the elastic distortion energy $F(t_{xy})$, which makes the updated Hamiltonian depend on the distortions of the lattice. The configurations $\{t_{xy}\}$ that minimise the ground state energy, including all periodic, nonperiodic and chaotic structures depend on a Peierls-type of instability of the bond lengths similar to the one employed in the archetypical polyacetylene chain model [162, 163]–also known as the SSH model.

Frank and Lieb demonstrated that only periodic, reflection-symmetric distortions are allowed on the $\sqrt{3}\times \sqrt{3}$ supercell [161] in the thermodynamic limit, and figure 15 presents all the possible structures satisfying such criterion.

A phase diagram based on (quantum) Diffusion Montecarlo (DMC) simulations is provided in figure 16 [160]. The figure presents a comparison between the ground state energies of the competing phases seen in figure 15 for isotropically strained graphene, as obtained from DMC and DFT, and it also includes the magnitude of the enthalpies when the lattice is under stress. DMC produces the lowest enthalpy when compared to unperturbed graphene. For uniform graphene strain in the 7%–15% range (corresponding to a ${\sim} 27$–31 N m⁻¹ stress), the lattice correlation-driven dimerisation freezes Pauling's resonating valence bond into a valence-bond solid realised by an insulating Kekulé-like dimerised phase (DIM) [160, 164, 165] that creates a topological band gap opening in isotropically strained graphene. Kekulé patterning has been observed by STM in strained graphene [166], and it seems to play a role in the superconducting phases of TBG [165, 167, 168].

Zoom In Zoom Out Reset image size

Strong magnetic fields can also create Kekulé distortions in graphene: Using scanning tunnelling spectroscopy (STS), a continuous quantum phase transition from a valley-polarised state onto an intervalley coherent state with a Kekulé distortion was observed under magnetic fields [169]. Importantly, the valley texture extracted from STS measurements of the Kekulé phase revealed valley skyrmion excitations localised near charged defects [169].

Gamayun et al obtained a low-energy four band model from the complete six-band TB model for both Kek-O and Kek-Y patterns without strong electron-electron interactions [159]. Some features of the Y Kekulé pattern will be discussed next.

Pristine graphene and a structure hosting a Kek-Y pattern–indicated by red and black bonds–are shown in figures 17(a) and (b), respectively. As indicated in figure 14, the Kek-Y unit cell contains six atoms and lattice vectors $\tilde{\mathbf{a}}_1$ and $\tilde{\mathbf{a}}_2$ in figure 17(b). The corresponding reciprocal lattice vectors are scaled as $\tilde{\mathbf{K}}_{\pm} = (1/\sqrt{3})\mathbf{K}_{\pm}$ [170].

Zoom In Zoom Out Reset image size

The system is modelled by a TB Hamiltonian with one π-orbital per carbon atom [1]. Going beyond first-nearest neighbour hopping [159], a second-nearest neighbour TB Hamiltonian was employed [172]:

$$\begin{align} H & = -\sum_{\mathbf{r},l} t^{0}_{\mathbf{r},l} \hat{a}_{\mathbf{r}}^{\dagger}\hat{b}_{\mathbf{r+\delta_l}}+\sum_{\mathbf{r}, m\neq n} t_2 \hat{a}_{\mathbf{r}}^\dagger \hat{a}_{\mathbf{r+\delta_m-\delta_n}}\nonumber\\ & \quad + \sum_{\mathbf{r},m\neq n} t^{\left(2\right)}_{\mathbf{r}} \hat{b}_{\mathbf{r}+\delta_m}^\dagger \hat{b}_{\mathbf{r+\delta_n}}+h.c., \end{align} \tag{ 42 }$$

where $\hat{a}_{\mathbf{r}}$ and $\hat{b}_{\mathbf{r}}$ are annihilation operators at site r for a carbon atom in the A or the B sublattice, respectively. Operators $\hat{a}^{\dagger}_{\mathbf{r}}$ and $\hat{b}^{\dagger}_{\mathbf{r}}$ are the corresponding creation counterparts. A bond-density wave modulates the TB parameters and originates the Kekulé distortion: the first-neighbour hopping is given by [159]:

$$\begin{align} t_{r,l}/t_0 = 1+2\text{Re}\left[\Delta e^{i\left(p\mathbf{K}_{+}+q\mathbf{K}_{-}\right)\cdot \mathbf{\delta}_l+i\mathbf{G}\cdot\mathbf{r}}\right], \end{align} \tag{ 43 }$$

where $t_0 = 2.8$ eV is the hopping parameter for first-nearest neighbours in pristine graphene and $\mathbf{G} = \mathbf{K}_+- \mathbf{K}_-$ is a reciprocal lattice vector coupling both valleys (see figure 17(b)). The amount of the bond modulation was $\Delta \approx 0.1$. In addition, p and q are two integers such that the index:

$$\begin{align} \nu = \left(1+q-p\right)\mathrm{mod} \hspace{0.1cm} 3 \end{align} \tag{ 44 }$$

distinguishes the Kekulé O pattern (ν = 0) from the two Y patterns ($\nu = \pm 1$) alluded to in figure 14.

The subsequent, second-nearest neighbour parameters were periodically modulated as:

$$\begin{align} t_{\mathbf{r}}^{\left(0\right)}/t_0 = t_{\mathbf{r}}^{\left(2\right)}/t_2 = 1+2 \Delta \cos\left(\mathbf{G} \cdot \mathbf{r}\right), \end{align} \tag{ 45 }$$

where $t_2 = 0.1t_0$ is the second-nearest neighbour TB parameter of pristine graphene.

A column vector is now defined:

$$\begin{equation*} c_{\mathbf{k}} = \begin{pmatrix} a_{\mathbf{k}}\\ b_{\mathbf{k}}\\ a_{\mathbf{k+G}}\\ b_{\mathbf{k+G}}\\ a_{\mathbf{k-G}}\\ b_{\mathbf{k-G}} \end{pmatrix} \end{equation*}$$

that collects annihilation operators at k and $\mathbf{k+G}$, so that the Hamiltonian can be written as a $6 \times 6$ matrix:

$$\begin{align} H = c_{\mathbf{k}}^{\dagger} \begin{pmatrix} H_{\Gamma} & T \\ T^{\dagger} & H_{G} \end{pmatrix} c_{\mathbf{k}}, \end{align} \tag{ 46a }$$

made from the original graphene Hamiltonian at the $\Gamma-$point:

$$\begin{align} H_{\Gamma} = \begin{pmatrix} f_{\mathbf{k}} & -\epsilon_{\mathbf{k}} \\ -\epsilon_{\mathbf{k}} & f_{\mathbf{k}}\\ \end{pmatrix}, \end{align} \tag{ 46b }$$

the Hamiltonians at the G and $\mathbf{-G}$ points:

$$\begin{align} H_{G} = \begin{pmatrix} f_{\mathbf{k+G}} & -\epsilon_{\mathbf{k+G}} & 0 & -\Delta \epsilon_{\mathbf{k-G}}\\ -\epsilon^*_{\mathbf{k+G}} & f_{\mathbf{k+G}} & -\Delta \epsilon^*_{\mathbf{k+G}} & \Delta f^-_{\mathbf{k-G}}\\ 0 & -\Delta \epsilon_{\mathbf{k+G}} & f_{\mathbf{k-G}} & -\epsilon_{\mathbf{k-G}}\\ -\Delta \epsilon^*_{\mathbf{k-G}} & \Delta f^+_{\mathbf{k+G}} & -\epsilon^*_{\mathbf{k-G}} & f_{\mathbf{k-G}}\\ \end{pmatrix}, \end{align} \tag{ 46c }$$

and the interaction between them:

$$\begin{align} T = \begin{pmatrix} 0 & -\Delta \epsilon_{\mathbf{k+G}} & 0 & -\Delta \epsilon _{\mathbf{k-G}}\\ - \Delta \epsilon^*_{\mathbf{k}} & \Delta f^-_{\mathbf{k+G}} & - \Delta \epsilon^*_{\mathbf{k}} & \Delta f^+_{\mathbf{k-G}} \end{pmatrix}, \end{align} \tag{ 46d }$$

where

$$\begin{align} \epsilon_{\mathbf{k}} & = t_0 \sum_{j = 1}^3 e^{i \mathbf{k} \cdot \boldsymbol{\delta}_j}, \end{align} \tag{ 47a }$$

$$\begin{align} f_{\mathbf{k}} & = t_2 \sum_{m\neq n} e^{i \mathbf{k} \cdot \left(\boldsymbol{\delta}_m-\boldsymbol{\delta}_n\right)},~\text{and} \end{align} \tag{ 47b }$$

$$\begin{align} f\,^{\pm}_{\mathbf{k}} & = t_2 \sum_{m\neq n} e^{i \mathbf{k} \cdot \left(\boldsymbol{\delta}_m-\boldsymbol{\delta}_n\right)}e^{\mp i\mathbf{G} \cdot \boldsymbol{\delta}_n}. \end{align} \tag{ 47c }$$

The spectrum can be found by diagonalisation of the $6 \times 6$ operator or by a projection technique that permits obtaining a low-energy $4\times 4$ matrix effective Hamiltonian. That effective Hamiltonian becomes for the Kek-Y case [172]:

$$\begin{align} \mathcal{H}& = v_F \left(1-\Delta^2\right) \left(\mathbf{p} \cdot \boldsymbol{\sigma}\right) \otimes \boldsymbol{\tau}_0 \nonumber\\ & \quad +v_F \Delta \left(1-\Delta\right) \boldsymbol{\sigma}_0 \otimes \left(\mathbf{p} \cdot \boldsymbol{\tau}\right)\\ &\quad + \frac{\mu}{2}\left(\boldsymbol{\sigma}_0 \otimes \boldsymbol{\tau}_0 + \boldsymbol{\sigma}_x \otimes \boldsymbol{\tau}_x +\boldsymbol{\sigma}_y \otimes \boldsymbol{\tau}_y - \boldsymbol{\sigma}_z \otimes \boldsymbol{\tau}_z\right),\nonumber \end{align} \tag{ 48 }$$

where $\boldsymbol{\sigma} = (\sigma_x,\sigma_y,\sigma_z)$ and $\boldsymbol{\tau} = (\tau_x,\tau_y,\tau_z)$ are two vectors with components defined from Pauli matrices, σ₀ and τ₀ are two $2 \times 2$ identity matrices, and µ is defined as:

$$\begin{align} \mu = \frac{9 \Delta^2 t_0^2 t_2}{t_0^2-9t_2^2}. \end{align} \tag{ 49 }$$

The resulting four low-energy bands turn out to be:

$$\begin{align} E_{D}^{\pm} & = \pm v_D |\mathbf{p}|,~\text{and}\nonumber\\ E_{M}^{\pm} & = \mu \pm \sqrt{v_M^2 \mathbf{p}^2 +\mu^2}, \end{align} \tag{ 50 }$$

were two new velocities have been defined: $v_D = v_F(1-\Delta)$, and $v_M = v_F(1-\Delta)(1+2\Delta) \approx v_F(1+\Delta)$.

As indicated in figure 17(b), the Kek-Y pattern folds the two valleys $\mathbf{K}_{+}$ and $\mathbf{K}_{-}$ into the Γ point when $t_2 = 0$. One valley is a fast Dirac cone with Fermi velocity v_M , and the other a slow Dirac cone with Fermi velocity v_D . Specific signatures on the optical and electrical conductivities, as well as in plasmons, are associated with those two charge carrier flavours [170, 173, 174]: the valley-dependent splitting of the Fermi velocities due to the Kekulé distortion leads to a similar splitting in the dielectric spectrum [173]. Moreover, an absorption phenomenon was found where a resonance peak related to intervalley transport emerges at a 'beat frequency' determined by the difference between characteristic frequencies of each valley [173].

Kekulé-patterning was observed in a work that combined STM and DFT calculations for both graphene monolayers and bilayers [175]. It was shown that strain induced a splitting of the Dirac cones and a band gap (it is not clear if the mechanism is the one predicted in [171] though). Flat bands measured with ARPES were shown to coexist with the bond order in Li-intercalated graphene [176].

The inclusion of second-nearest neighbour hopping changes this picture slightly, as one of the cones gets a small mass [172]. In a similar way, a gap opens in the energy dispersion for the Kek-O pattern [159]. Such band gap was experimentally confirmed by strain-induced Kekulé patterning [166].

Uniform strain has been proposed as a tool to engineer the valley degree of freedom [171]: while valleys are too far away to mix under realistic strain in pristine graphene [1], the folding induced by patterning reduces the separation between them, so that strain can modulate their interactions more easily.

As seen in figure 17(c), uniform strain moves the $\mathbf{K}_{\pm}$ points to $\mathbf{K}^{^{\prime}}_{\pm}$ in graphene [1]. Moreover, strain changes the symmetry of the Bravais lattice, and high-symmetry points of the reciprocal lattice must be labelled differently. As an example, the $\mathbf{K}_{\pm}$ points of the P6/mmm space group are replaced by the F₀ and Δ points in the Cmmm space group when under uniaxial strain, and the Dirac points $\mathbf{K}_{D}^{^{^{\prime}}\pm}$ of the strained lattice do not necessarily coincide neither with $\mathbf{K}_{\pm}^{^{\prime}}$, nor with F₀ or Δ. Figure 17(c) sketches those general observations to avoid confusion.

Uniform strain in Kekulé patterned graphene leads to the situation depicted in figure 17(d): the original valleys are folded into new cones located at $\tilde{\mathbf{K}}^{+}_D$ and $\tilde{\mathbf{K}}^{-}_D$. These new Dirac points are very close to the Γ point, and their separation can be tuned by experimentally accessible strain [171] to enhance valley interactions.

In analogy to electromagnetic fields, strong time-dependent gauge fields can be generated by time-varying strain in graphene. Time-dependent deformations can be induced by local probes such as a transmission electron microscope (TEM), or by SAWs created at a proximal piezoelectric layer [4].

The possibility of controlling electron transport by time-dependent gauge fields was suggested in [1, 177], and phonons were seen to generate effective pseudoelectromagnetic waves [60, 177]. Since then, acoustic waves have been shown to induce a giant synthetic Hall voltage without external magnetic fields [4], and the experimental setup is depicted in figure 18(a): graphene lays on top of LiNbO₃, which itself lays on an insulating SiO₂ substrate, which sits on a p-doped silicon wafer. The carrier concentration is tuned by a voltage V_BG applied to the $p-$doped substrate.

Zoom In Zoom Out Reset image size

The device can interchangeably be used in a magnetotransport Hall configuration under a real magnetic field B, or in an acoustic transport configuration. In the latter configuration, acoustic waves are induced in the graphene by an interdigitated resonator. Figure 18(b) displays a typical Hall measurement with Hall resistivity R_xy and magnetoresistance R_xx . The plateaus in R_xy indicate the presence of a quantum Hall effect. Figure 18(c) presents the acoustic induced Hall effect, seen in the current $I_\mathrm{SAW}$.

SAWs can also provide alternative ways to study topological properties [178–180] and the topological non-adiabaticity at the Dirac point [181]. The issue at hand is that the main assumption of topological phases is the adiabaticity of the time-dependent change of external parameters of a quantum system, which could require a bigger bandgap than the energy provided by external probes. However, there is no bandgap in Dirac materials, and transitions are induced at all frequencies. This produces a region of non-adiabaticity around the cones that requires additional theory to be properly understood [181–183].

Extensions to the strain engineering of other 2D Dirac semimetals–which might be more sensitive than graphene and could thus function as strain sensors–are underway [184].

Strain can also appear in the non-linear conductivity of gapped graphene [185].

When both time reversal and inversion symmetry are present, these higher order phenomena vanish because the Berry curvature Ω is zero [186].

Thus, breaking of inversion symmetry becomes necessary [187] if no magnetic field is present (i.e. when graphene's Hamiltonian is invariant under time reversal). If space symmetries are further reduced by strain, the Berry curvature becomes asymmetrical in reciprocal space, and the integral $\int \mathrm{d}\boldsymbol{k} f_0\partial_{k_x}\boldsymbol{\Omega}$ (the so-called Berry dipole, with f₀ a Fermi–Dirac distribution) turns nonzero. Sodemman and Fu [188] have shown that (under time-dependent bias) the Berry dipole contributes to the second order conductivity in the direction normal to the bias, that is, in the same way as the Hall effect, the resulting second order transversal conductivity tensor being proportional to the Berry dipole.

The stretching of the crystal structure (as a consequence of strain) implies the reduction in the number of crystal symmetries of the lattice. In Bernal-stacked bilayer graphene under a normal electric field, or in graphene on top of hBN, sublattice equivalence is broken and so is inversion symmetry, reducing the symmetry from $C_{6v}$ to $C_{3v}$. Uniaxial strain reduces the symmetry further, from $C_{3v}$ to $C_{2v}$ [189]. When combined with the wrapping of the bands, this leads to an effective Hamiltonian with a nonzero Berry dipole [190].

2D semiconductors are an important platform for the emission of light because of the reduced dielectric screening in two dimensions, the intrinsic absence of total internal reflection within their atomically-thin thickness [217], and because local tensile strain is the only mechanism to move through space (i.e. to funnel) charge neutral electron–hole pairs (i.e. Wannier excitons, or excitons for short) created in monolayers [218, 219]. Furthermore, the deterministic and position-dependent creation of few or single photons is one of the main goals of this field of research (see, e.g. [210]). A recent review on this subject is available [220], and the following paragraphs are mainly an update. A discussion of semiconducting TMDCs monolayers under strain follows, while TMDC monolayers encapsulated within a van deer Waals heterostructure are the subject of section 4.6.

Electron and hole effective masses are comparable in excitons, which makes those oppositely-charged quasiparticles wander for long distances until they break apart as they collide with phonons. The collision processes (and crystal imperfections), combined with exciton-polariton excitations give rise to a broad optical emission spectrum aided by multiple successive exciton creation and reabsorption events. Such broad emission spectrum is customarily known as photoluminiscence [221] and most experimental studies to be described now will be based on this optical technique. Importantly, the photoluminiscence spectrum depends on the quality (impurity density) of the semiconductor or insulating 2D material.

Mukherjee et al placed a WSe₂ monolayer over a regular array of nanopillars separated by 4 µm and having heights of 100 nm. They observed a strain gradient in the vicinity of nanopillars, which periodically localises bright excitons in the WSe₂ monolayer at 4.2 K [29]. As confirmed by auto-correlation measurements with a Hanbury–Brown–Twiss interferometer, localised emission lines turned out to be sources of quantum light. The need for spectral control of the emitted light calls for the encapsulation of the WSe₂ monolayer within a van der Waals stack, and those results will be discussed in section 4.6.

Strain and defect engineering were simultaneously pursued in an effort to increase operating temperatures, photon yield, and purity of single-photon emitters on WSe₂ monolayers by Parto et al [222]. Those emitters show biexciton cascade emission, single-photon purity above 95%, and working temperatures up to 150 K [222].

As depicted in figure 23(a), the spin–orbit coupling in TMDC monolayers splits the conduction band with an optical dipole-allowed bright exciton, but it also leaves the possibility of creating a (spin-forbidden) dark exciton transition at the K and Kʹ points [223–227]. In this context, the bright excitons correspond to a transition among states of the same spin, while the dark excitons correspond to transitions between states with opposite spins. In general terms, a 'dark' exciton is a light excitation requiring additional perturbations–the spin-flips described thus far, and/or momentum transfer [228]—to occur.

Zoom In Zoom Out Reset image size

In the WSe₂ monolayers studied by Hasz et al, the dark exciton has a lower energy (${\sim} 40$ meV) than that of the bright exciton [38]. Dark excitons have a longer radiative lifetime than dipole-allowed bright excitons due to their spin- (and/or momentum) forbidden transitions [38, 229] and hence may function as coherent two-level quantum systems, or as Bose–Einstein excitonic condensates.

Excitonic states can be resonantly enhanced by the use of an AFM metallic tip sampling with a metal surface underneath (figure 23(b)) through the Purcell effect [38], which is the resonant increase of the spontaneous emission rate usually created within a cavity [230]. This technique is especially insightful to scan the photoluminiscence created at defects, grain boundaries, edges, nanobubbles (see figures 23(c) and (d)), and under uniform strain [231]. Additional strain induced during growth on silica has been shown to red-shift Raman features and to create a slowing down of exciton dynamics [33].

Under both bright and dark excitons, and as seen on figures 23(e) and (f), greater tensile strain increases peak intensity and decreases peak energy, while compressive strain increases the intensity and increases peak energy (figures 23(g) and (h)) [38].

Strain can help tune the electronic bandgap of TMDC monolayers [232] and inhomogeneous local strain results in bent band structures with a band gap gradient–the magnitude of the band gap is location-dependent–and bright excitons are driven toward a region of maximum strain [233]; those excitons have been seen to move over distances up to a micrometer.

Working with a WS₂ monolayer placed onto a regular array of nanopillars, momentum-forbidden dark excitions move away from the region of maximum strain in an effect dubbed anti-funnelling. Anti-funnelling is a consequence of a valley-dependent exciton energy within the funnels [30]. In this experiment, a maximum strain of 0.6% takes place at the centre between two pillars, rather than directly atop of pillars. The WS₂ monolayer lies flat there, providing a homogeneous local dielectric environment.

Because of funneling, dark $K\Lambda$ excitons reach spatial regions in which spectral separation of dark and bright excitons turns small, locally enhancing phonon-driven intervalley scattering. This way, originally dark $K\Lambda$ excitons scatter into bright KK states, which turn visible in the photoluminiscence spectra [30]. An additional study showed that excitons created on a MoS₂ monolayer do not feature anti-funnelling, thus emphasising the crucial role of the electronic structure of a WS₂ monolayer (having almost degenerate conduction band edges at the K and Λ k-points) to achieve this effect.

Localised strain was created in encapsulated WSe₂ monolayer nanobubbles, which were demonstrated to operate as single photon sources [37, 234]. Lastly, the conversion of excitons under non-uniform strain to more complex particles (multiexcitons) such as trions has been demonstrated in WS₂ monolayers [235].

The electric polarisation P created by breaking centrosymmetry in graphene (section 3.2.3), or existing within a nanobubble in a hBN monolayer (section 3.4), requires an extrinsic source of energy (strain) to take place.

There are at least three intrinsic ferroic orders: ferroelasticity, ferromagnetism, and ferroelectricity. Ferroelasticity is the ability of a material's unit cell to change its orientation upon the application of strain, and ferroelastic 2D materials will be discussed in section 3.6.2. In analogy with ferromagnets, ferroelectrics lack a crystalline center of inversion and develop an intrinsic electric dipole P that switches orientation by the application of an external electric field [236]. Historically, ferroelectrics have been insulating ternary materials (i.e. they contain at least three chemical elements), and the interplay between electronic properties and ferroelectric ones has been seldom addressed for that reason. (2D ferroelectric monolayers under strain are covered from sections 3.6.3–3.6.8.) The external fields used to control (ferroic) orders are strain (u), electric fields (E), and/or magnetic fields (B) that can switch the crystal structure and/or the magnetic configuration within the available degenerate (or nearly degenerate) quantum states. The conjugated inner fields are the stress (σ), the intrinsic electric dipole (P) and/or the magnetisation (M). Combinations of at least two ferroic orders in a material lead to a multiferroic; incipient 2D multiferroics will be covered in sections 3.8 and 4.9.

Just as it was the case thus far, the Update has been structured in a way that coverage of monolayers will be followed by a discussion of bilayers or few-layer ferroelectrics in sections 4.1, 4.2 and 4.5. Layered antiferroelectrics will be covered in section 4.3.5.

Graphene turns out to be an exception to a trend of structural degeneracies taking place in multiple 2D materials [240]. Even a hBN monolayer, which has a honeycomb structure with two sublattices, is degenerate: its structural energy remains unchanged regardless of whether (i) the A-sublattice contains boron atoms and the B-sublattice contains nitrogen atoms or (ii) the A-sublattice contains nitrogen atoms and the B-sublattice contains boron atoms.

Nevertheless, a hBN monolayer is not ferroic, because ferroic behaviour requires the switching of degenerate structures by strain or temperature. But exchanging all the boron and nitrogen atoms in a hBN monolayer necessitates removing (breaking) and restructuring back all chemical bonds: Trying to switch boron and nitrogen positions will rather burn that 2D crystal down. As it will be next discussed, ferroic 2D materials (having switchable degenerate ground states) do exist.

Ferroic materials are prone to undergo structural phase transitions different from melting at finite temperature, and the structural transformations in figure 24 highlight three ferroic materials: (a) ferroelastic silicene [191, 194], (b) the quantum paraelectric SnO monolayer [237, 241], and (c) a ferroelectric and ferroelastic SnSe monolayer [239, 242–245]. Subplots (i) in figure 24 display energy paths joining two degenerate ground state structures through the smallest possible barrier J, while subplots (ii) in figure 24 show the two degenerate structures and a unit cell with the optimal configuration at the height of the barrier J (which has a higher structural symmetry) explicitly.

Zoom In Zoom Out Reset image size

The following similarities are found: (1) All ferroic materials have energy landscapes with structural degeneracies. In the case of silicene, the degeneracy is created upon a reflection with respect to the z-axis, which swaps the signed vertical separation between atoms in the A and B sublattices. The SnO and SnSe monolayers display a degeneracy upon exchange of their two orthogonal lattice vectors. (2) It is possible to trace paths of steepest descent joining the two degenerate minima, creating Landau-like double-well potentials. At the height of the barrier J, one reaches an atomistic structure with an enhanced symmetry, which is a graphene-like (planar) structure labelled $\Delta_z = 0$ for silicene, and a square unit cell labelled C for the SnO and SnSe monolayers.

Now, if the barriers J are sufficiently large (say, larger than 100 K per primitive unit cell) but also sufficiently small (say, smaller than 1000 K per primitive unit cell), then 2D structural transformations into average structures with an enhanced symmetry become possible, as exemplified in the three subplots labelled (iii) in figure 24, which were obtained by ab initio molecular dynamics [26, 194, 237–239]. The structural transformation of the SnSe monolayer was expressed in terms of the parameter $\Delta \alpha$, a shear strain to be defined in section 3.6.4.

Traversing the landscape from one minima onto the height of the barrier requires the unit cell to vary: the 2D transformations depicted in figure 24 are driven by strain.

Until about a decade ago, most known ferroelectrics were insulating ternary compounds that lost their ferroelectric properties down a few-nanometer-thickness limit. A field of 2D ferroelectrics took off shortly before [1] and, as it will be discussed in what follows, there are ferroelectrics that are only two-atoms thick nowadays (and one or two monolayers thick). Some of those 2D ferroelectrics are binary compounds, but even elemental 2D materials have been proven to be ferroelectrics recently. Ferroelectric materials tend to undergo structural phase transitions at finite temperature, and 2D transitions aided by shear strain will feature in section 3.6.4.

2D ferroelectrics can be semiconductors, semimetals, and even metals. They can display non-trivial topological band structures that couple with P. Those materials command a strong interest for those reasons [246, 247].

Honeycomb lattices are prevalent on many 2D materials such as graphene, hBN monolayers, silicene, germanene, and TMDC monolayers with 1H and 1T atomistic configurations.

Nevertheless, other important and experimentally available 2D materials such as 1T' TMDC monolayers [249, 250], phosphorene [251–253], SnO multiferroic monolayers [237, 241], group-IV monochalcogenide monolayers [238, 242, 243, 246, 248, 254–256], bismuth monolayers [257], and some antiferromagnetic monolayers [258] possess rectangular unit cells. 2D materials with rectangular Bravais lattices are degenerate in energy, as the total energy of the unit cell remains the same upon exchange of x- and y-coordinates.

As shown in figure 25(a), rectangular unit cells can be circumscribed within a rhombus: the short and long diagonals (with magnitudes $2a_2$ and $2a_1$, respectively) form a 90^∘ angle, and the AB and CD sides of the rhombus raise by an angle $\Delta \alpha$ (known as the rhombic distortion angle) away from the horizontal (line AE). Shear strain can turn the rhombus onto a square, whose diagonals are now identical (and set to a_C in figures 24(b)-(ii) and (c)-(ii)) and $\Delta \alpha = 0$.

Zoom In Zoom Out Reset image size

Some materials with rectangular unit cells undergo two-dimensional structural phase transformations at finite temperature, whereby the 2D unit cell suddenly changes from a rectangle onto a square: figure 25(b) depicts $\Delta \alpha$ experimentally determined as a function of temperature [248] for a ferroelectric and ferroelastic SnTe monolayer on epitaxial graphene. Its coalescence to zero is the tell sign of a structural transformation of the two-dimensional SnTe monolayer from a rectangular unit cell onto a square one.

This coupling of the order parameter P_x to strain has been neglected in some theoretical works, where it was instead coupled to an optical phonon on a structure with fixed lattice vectors; see [246] for more details.

As indicated since [242], it is more sensible to couple the magnitude of the intrinsic electric dipole P to the in-plane lattice parameters on 2D ferroelectrics with an in-plane polarisation: $\mathbf{P} = \mathbf{P}(a_1,a_2)$ [238, 243, 246]. Such a direct coupling is the only one consistent with the experimental observation that the rhombic distortion angle (displayed in figure 25(b)) becomes zero. In other words [238, 248]:

$$\begin{align} \frac{a_1}{a_2} = \frac{1+\sin\left(\Delta \alpha\right)}{\cos\left(\Delta \alpha\right)}, \end{align} \tag{ 54 }$$

so $a_1 = a_2=a_C$ when $\Delta \alpha = 0$, and $\mathbf{P}(a_C,a_C) = (0,0,0)$ (figure 25).

The sudden changes in lattice parameters at a critical temperature T_c depicted in figure 24(c) are different from the isotropic thermal expansion seen on graphene, hBN monolayers, or TMDC monolayers. This observation calls for a definition of strain at finite temperature, which is provided now.

Thus far, the theoretical discussion has considered 2D materials at zero temperature. In that context, strain has been implicitly thought of as anything that changes (locally or globally) the atomistic structure of a given 2D with respect to its magnitude at zero temperature (a₀):

$$\begin{align} \varepsilon_{i} = \frac{a_i-a_{0,i}}{a_{0,i}}. \end{align} \tag{ 55 }$$

Nevertheless, the sudden compression of the lattice parameter a₁ and the elongation of a₂ at T_c driving the collapse of $\Delta \alpha$ to zero in figures 24(c) and 25(b) reflect the optimal, lowest-energy atomistic configuration at finite temperature of 2D ferroelectrics with in-plane polarisation, and strain should be measured against the structure in thermal equilibrium. In other words, strain should be temperature-dependent, and equation (55) should be modified as follows:

$$\begin{align} \varepsilon_{i}\left(T\right) = \frac{a_i-\langle a_0\left(T\right)\rangle_i}{\langle a_0\left(T\right)\rangle_i}, \end{align} \tag{ 56 }$$

where $\langle a_0(T)\rangle_i$ is the mean lattice parameter at finite temperature. This modification of strain applies to any 2D material undergoing 2D phase transformations [259].

To go beyond the zero-T paradigm to elasticity, one first obtains an analytical expression for the energy landscapes (subplots (i) in figure 24), and the discussion that follows focuses on the energy landscape depicted in figure 24(c)-(i). That landscape is symmetric with respect to the dashed diagonal line, and a new coordinate system (X, Y) is defined to reflect such symmetry.

The analytical form for the energy landscape $U(X,Y)$ is given in [26]. With it, one can build a function $f(X,Y)$ (for example, the lattice parameter $a_{1}(X,Y)$) and calculate its mean value within the energy landscape $\langle f(U_\mathrm{max})\rangle$ as an average over classically accessible states [26]:

$$\begin{align} \langle\, f\left(U_\mathrm{max}\right)\rangle = \frac{\oint e^{-U\left(X,Y\right)/U_{max}} f\left(X,Y\right) \mathrm{d}X\mathrm{d}Y}{\oint e^{-U\left(X,Y\right)/U_\mathrm{max}} \mathrm{d}X\mathrm{d}Y}, \end{align} \tag{ 57 }$$

with $\mathrm{d}X\mathrm{d}Y$ an area element within the confines of an isoenergy contour $U_\mathrm{max}$ around structure A in figure 24(c)-(i).

$U(X,Y)$ here is a classical potential energy, and one samples accessible crystalline configurations within isoenergy confines. When $U_\mathrm{max}\unicode{x2A7E} J$ nevertheless, the average structure encompasses minima A and B in figure 24(c)-(ii), yielding $\langle a_1\rangle = \langle a_2\rangle$, and it thus is a square (see figure 26(a)). The sampling over independent crystalline unit cells is a limitation of the model, as 2D structural transformations in 2D are driven by disorder [239, 242, 243, 246].

Zoom In Zoom Out Reset image size

Thermal averages of elastic constants C_ij are then determined by [260]

$$\begin{align} \langle C_{ij}\left(U_\mathrm{max}\right)\rangle & = k_\mathrm{B}\left\{\bigg\langle\frac{\partial^2u}{\partial\varepsilon_i\partial\varepsilon_j}\bigg\rangle\right.\\ & \quad \left. -\frac{1}{U_\mathrm{max}}\bigg[\bigg\langle \mathcal{A}\frac{\partial u}{\partial\varepsilon_i}\frac{\partial u}{\partial\varepsilon_j}\bigg\rangle- \langle \mathcal{A}\rangle \bigg\langle \frac{\partial u}{\partial\varepsilon_i}\bigg\rangle \bigg\langle \frac{\partial u}{\partial\varepsilon_j}\bigg\rangle \bigg]\right\},\nonumber \end{align} \tag{ 58 }$$

(with $k_\mathrm{B}$ Boltzmann's constant, $\mathcal{A} = a_1a_2$, and $u = U/\mathcal{A}$) and presented in figure 26(b).

At energies above $U_{max}$ the 2D material transitions from a rectangular onto a square unit cell and $\langle C_{11}\rangle$ and $\langle C_{22}\rangle$ must be identical, but that is not obtained when employing the regular expression for strain based on a zero-temperature value (equation (55)). The need for $\langle C_{11}\rangle = \langle C_{22}\rangle$ for $U_\mathrm{max}\gt J$ dictates the redefinition of equilibrium lattice parameters evolving with temperature (equation (56)), and hence of a thermally-evolving strain.

$\langle C_{12} \rangle$ is the softest elastic modulus and $\langle C_{11}\rangle$ hardens at the transition ($U_\mathrm{max} = J$), while $\langle C_{22}\rangle$ softens at $U_\mathrm{max} = J$. According to figure 26(b), SnSe monolayers are softer than graphene, for which $C_{11}^{(0)} = C_{22}^{(0)} = 336$ N m⁻¹, and $C_{12}^{(0)} = 75$ N m⁻¹.

By direct comparison J and T_C from numerical calculations [239], a correspondence $T\propto 1.42 U_\mathrm{max}$ is established, such that $T_C = 212$ K, and finite−T elastic behaviour can be extracted from figure 26(b).

The discussion now turns into optical ways to modify lattice parameters in a non-thermal manner, but with illumination.

Illumination by light can excite electrons onto the conduction band, changing the electronic density away from its ground state configuration. The modification of the electronic density–in turn–screens the intrinsic electric dipoles of 2D ferroelectrics, producing uniform and anisotropic strain as a result. The process, called photostriction, was theoretically predicted on group-IV monochalcogenide (SnS, SnSe, and GeS) monolayers [7].

The atomistic structure of a SnS monolayer is shown in figure 27(a), and its electronic band structure in figure 27(b). The material has an indirect bandgap, with the top of the valence band at the valley dubbed nX (which stands for near X-point) and the bottom of the conduction band at the valley labelled nY (for near Y-point).

Zoom In Zoom Out Reset image size

Indirect optical transitions require momentum transfer and have low probability of occurrence. For this reason, one considers direct optical transitions instead. Three such transitions are indicated by arrows on figure 27(b). The local valence band valleys at the nX, Γ, and nY k-points have a predominantly s-orbital (spherically-symmetric) hybridisation, and the largest change in intrinsic polarisation occurs upon transitions onto p_x and p_y orbitals. The electronic wavefunction densities at the conduction band edges at the nX k-point is shown in figure 27(c) [7].

Working on a $n_k\times n_k$ k-point grid and considering spin–orbit coupling, one gets a density of charge carriers $n_c = 1/(n_k^2A_0)$, where A₀ is the area of the unit cell before the photoexcitation. Setting $n_k = 41$, a photoexcitation around the nX k-point is simulated by sequentially changing of occupation of the valence band edge by the removal of the occupation at points labelled 1, 2, and 3 within the inset in figure 27(d), promoting the removed electronic density into the conduction band at the same k-point(s), and relaxing the structure. This way, the intrinsic polarisation is shown to decrease as a result of the photoexcitation, which screens P (see its decreasing value in figure 27(d)). Furthermore, while the unit cell remains rectangular with lattice parameters $a_1\gt a_2$, one observes in figure 27(e) that a₂ increases in magnitude while a₁ decreases. This way, one achieves an optomechanical coupling to produce anisotropic strain by illuminating 2D ferroelectrics. The realisation of photostriction in 2D materials is important because the photostriction response time depends on the sample's thickness, ranging from picoseconds in thin films up to a few seconds in the bulk [7].

The effect described in figures 27(d) and (e) persists upon a reflection of the structure depicted in figure 27(a) with respect to the yz plane; an operation that swaps the direction of the intrinsic electric dipole from the $+\hat{x}$ onto the $-\hat{x}$ direction. Then, photostriction is not limited to monolayers, but it persists in bulk layered (orthorhombic) IV–VI compounds with an antipolar dipole configuration [246]. Briefly breaking the logical flow of the update to discuss experiments carried out on a bulk sample now, a technique called megaelectronVolt ultrafast electron diffraction was employed to confirm an expansion of a₂ and a contraction of a₁ on bulk (layered) GeS [8]: measured strains of 0.1% are consistent with those predicted in [7]. Furthermore, and as seen in figure 28, they also report an ultrafast photoexcitation of the order of 10 picoseconds.

Zoom In Zoom Out Reset image size

Similar to graphene, it may be possible to strain 2D ferroelectrics by clamping them onto a substrate and a subsequent bending (figure 29(a)) and one wonders how strain affects the critical temperature T_c at which the intrinsic electric dipole turns zero. To determine the effect of strain on T_c , multiple structural order parameters such as angles and interatomic distances of a ferroelectric SnSe monolayer were tracked as a function of temperature (figure 29(b)) through ab initio molecular dynamics calculations.

Zoom In Zoom Out Reset image size

The results are presented in figure 29(c) for strain applied along the x-direction, and in figure 29(d) when strain is applied vertically. Datapoints are results from molecular dynamics, and solid lines (when present) correspond to a Potts model [238, 261]. The extreme right within the yellow boxes to the left of figures 29(d) is the value of T_c without applied strain.

As seen in figure 29(c), a modest strain of 2% along the direction defined by the largest lattice vector a₁ increases the value of $\Delta\alpha$–from about 2^∘ in the non-clamped SnSe monolayer (figure 24(c)-(iii))–up to about 3.3^∘; see subplot (ii) in figure 29(c).

Still considering subplot (ii) in figure 29(c), one sees that clamping precludes a situation in which $\Delta\alpha = 0$. In this case, and as illustrated in figure 29(b), a transition from a rectangular onto another rectangular unit cell takes place though. In those rectangular cells, the transition takes place through the coalescence of angles α₁ and α₃, and of distances d₂ and d₃. T_c increases regardless of the direction in which uniaxial tensile strain is applied. (Symbols $\langle \cdot \rangle$ stand for thermal averages.)

There are few studies addressing the electronic properties of group-IV monochalcogenide monolayers at domain walls [25, 245, 248, 262], and the amount of strain at those appears to be rather small.

The attention now turns to two atomic layers of bismuth, in which a new frontier in 2D elemental ferroelectrics is arising [263]. (On a bilayer configuration, this material couples a non-trivial electronic topology with ferroelectric behaviour [264].) Similar to SnTe, bulk bismuth has a rhombic structure in the bulk [265, 266]. In the recent past, a structural transformation on ultrathin SnTe was determined whereby it turns into a (layered) orthorhombic phase instead [267, 268].

Consisting of a single chemical element, if a similar rhombic to orthorhombic transition is to take place in bismuth (whose electronic configuration is 6s²6p³), one may first imagine that its 2D structural configuration may resemble that of phosphorene (electronic configuration 3s²3p³), a non-ferroelectric 2D material possessing a center of inversion [251–253].

Nevertheless, and as stated in [263], the larger principal number on bismuth–and its disposition within the Periodic table of the Elements in between insulators and metals–makes its electronic configuration different from that of phosphorus, to the point of developing a net in-plane intrinsic dipole $\mathbf{P}$: while the 3s and 3p orbitals lie at comparable energies in phosphorene, the 6s orbitals lie at a significantly lower energy than the 6p orbitals, which weakens their sp : unlike phosphorene, in which bonds are shared s and $p-$orbitals, atomic bonding on Bi monolayers is mainly due to $p-$orbitals. Furthermore, the originally degenerate p_z orbitals split by occupying the lowermost and highermost atoms within the bilayer, in an electronic configuration resembling that of group-IV monochalcogenide monolayers whereby the chalcogen is more negatively charged: this intrinsic electronic redistribution (polarisation) creates two sublattices and drives the intrinsic electric dipole with an in-plane polarisation in bismuth monolayers.

The atomistic configuration of the bismuth monolayer in figure 30(a) thus resembles that seen in group-IV monochalcogenide monolayers [246] (figures 24(c)-(ii), 27(a), and 29(b)). Unlike phosphorene–similar to the structure B in figure 30(b)—it features different heights for all four atoms in the unit cell, and a twice-degenerate energy landscape reminiscent of that of figure 24(c)-(i) that is depicted in figure 30(c).

Zoom In Zoom Out Reset image size

Even though there are four bismuth atoms in the unit cell, the two ones closest to vacuum receive extra charge from p_z orbitals, giving the unit cell two nonequivalent atoms and a Pmn2₁ space group, just like that taking place in group-IV monochalcogenide monolayers [246].

Structure A' (B) here was labelled B in figure 24(c)-(i). Structure B here is calculated on the same rectangular unit cell used to compute structures A and A', while structure C in figure 24(c)-(i) was calculated on the minimal-energy square unit cell and are thus not equivalent. The charge distribution underpinning the electronically-driven intrinsic dipole is displayed in figure 30(d).

Similar to the experimental realisation of SnTe monolayers [248], bismuth monolayers were also grown in epitaxial graphene. The authors verified its in-plane intrinsic electric polarisation by noticing that the higher-most atom appears bright (dark) in valence (conduction) band sweeps [257], thus confirming the non-equivalence of bismuth atoms in the unit cell.

Furthermore, polarisation switching by the field created by the tip was observed at a domain wall, which verifies the ferroelectric character of 2D bismuth experimentally.

Unlike the case of the SnTe monolayer [262, 267], and similar to SnSe monolayers on epitaxial graphene [245], the bismuth monolayer displays 180^∘ domain walls [257].

The word moiré has an arabic origin and it means 'to wet.' The term is commonly used within the textile industry to describe the appearance of interwoven silk and it should not be capitalised.

The first observation of moirés on turbostatic graphite dates back to 1993 [282]. (This moiré has been revisited quite recently [283].) According to [284], the first observation of unusual electronic properties on rotated bilayer graphene moirés is due to Eva Andrei [285], who determined that the energy separation between van Hove singularities reduces with the relative angle of rotation among monolayers, and posited that their energy separation would be zero at an angle smaller than 1.1^∘ (a magic angle), at which point the linear momentum of charge carriers would turn zero as well. Calculations of a rotated graphene bilayer performed within a continuum model demonstrated the possibility of renormalisation of the Fermi velocity of graphene [286]; velocity renormalisation was confirmed by tight-binding calculations [287, 288]. Bistritzer and MacDonald suggested the possibility of strongly correlated electron–electron quantum phases at the flat bands taking place at magic angles [289]. The subsequent experimental confirmation of unconventional superconductivity and correlated insulating phases by Jarillo-Herrero's group [290] has led to intense and ongoing research in this area. Nowadays it is known that velocity renormalisation also takes place in other moiré 2D materials, such as a twisted hBN bilayer [291]. The experimental challenges in creating precise rotations have been addressed in a recent review [281].

There are multiple ways to create moirés [292] with 2D materials out of the following basic 'ingredients:'

• (i)  
  by a relative twist between two monolayers in a homobilayer
• (ii)  
  by isotropically straining one monolayer in a homobilayer
• (iii)  
  by relative shear strain in a homobilayer
• (iv)  
  by uniaxial simple shear in a homobilayer
• (v)  
  by an unrelaxed hetero-bilayer
• (vi)  
  by stacking more than two monolayers with angular stacking faults

The first three of those alternatives are shown in figure 36.

Zoom In Zoom Out Reset image size

The rotated homo-bilayer moiré is the most commonly described in the theoretical literature, and the process to create a commensurable moiré bilayer by a relative rotation is described now. When commensuration between the two monolayers exists, the lattice vectors of a rotated moiré supercell can be obtained following a process delineated in [282, 292–295] and multiple other publications. Let

$$\begin{align} S = a \begin{pmatrix} -\frac{\sqrt{3}}{2} & \frac{\sqrt{3}}{2}\\ \frac{3}{2} & \frac{3}{2} \end{pmatrix} \end{align} \tag{ 61 }$$

be a matrix built from the lattice vectors written down in equation (1). (Lattice vectors were written in column form in equation (61).)

Let θ be the mismatch angle (this is, the relative rotation angle between the two monolayers). Customarily, the lower monolayer is rotated by $-\theta/2$, while the upper monolayer is rotated by $\theta/2$. The respective rotation matrices are:

$$\begin{align*} R\left(\pm\theta/2\right) = \begin{pmatrix} \cos\left( \frac{\theta}{2} \right) & \mp \sin\left( \frac{\theta}{2} \right)\\[6pt] \pm \sin\left( \frac{\theta}{2} \right) & \cos\left( \frac{\theta}{2} \right) \end{pmatrix}, \end{align*}$$

and the moiré superlattice vectors (figure 37(a)) are obtained as:

$$\begin{align} S_M\left(\theta\right) = \left( S_1\left(\theta\right)^{-1} - S_2\left(\theta\right)^{-1} \right)^{-1}, \end{align} \tag{ 62 }$$

where

$$\begin{align} S_1\left(\theta\right) = R\left(-\theta/2\right)S,~\text{and }S_2\left(\theta\right) = R\left(\theta/2\right)S \end{align} \tag{ 63 }$$

Zoom In Zoom Out Reset image size

are the rotated lattice vectors for the lower and upper monolayer, respectively. Carrying out the math, the moiré lattice vectors (equation (62)) are given by:

$$\begin{align} S_M\left(\theta\right) = \frac{\sqrt{3}a}{2\sin\left(\theta/2\right)} \begin{pmatrix} -\frac{\sqrt{3}}{2} & \frac{\sqrt{3}}{2}\\[6pt] -\frac{1}{2} &-\frac{1}{2} \end{pmatrix}, \end{align} \tag{ 64 }$$

or

$$\begin{align} \mathbf{a}_{M,1}\left(\theta\right) = & a_M\left(\theta\right)\left( -\frac{\sqrt{3}}{2}, -\frac{1}{2} \right),~\text{and }\nonumber\\ \mathbf{a}_{M,2}\left(\theta\right) = & a_M\left(\theta\right)\left( \frac{\sqrt{3}}{2}, -\frac{1}{2} \right), \end{align} \tag{ 65 }$$

with the length of the moiré superlattice vector given by $a_M(\theta) = \frac{\sqrt{3}a}{2\sin(\theta/2)}\simeq \frac{\sqrt{3}a}{\theta}$ for $\theta\lt15^{\circ}$ (when expressed in radians). The area of the moiré is given by $\mathrm{Area}_M(\theta) = [\frac{\sqrt{3}}{2}a_M(\theta)]^2\simeq \frac{3a^2}{2}\theta^{-2}$, and the corresponding reciprocal space is shown in figure 37(b).

Strain and corrugations play a major role in the TBG flat band and energy gaps at the magic angle [298]. As seen in figures 36(a)–(c), different sections of the bilayer form relative conformations that explore every local unit cell configuration seen in figure 33. In principle, the energy cost varies, and certain configurations can be minimised by a process of local atomistic (strain) reconstruction [17, 283, 299–304]; Cazeaux et al have deployed a process to optimise moiré bilayers that minimises a functional of the local strain [292]. It relies on an energy landscape set on a unit cell similar to the one given by equation (60) (but slightly simpler). Configurations (a)–(c) in figure 36 lead to triangular homo-bilayer moirés and to experimentally differentiable patterns.

The energetic preference for AB and BA regions over AA stacking regions (see figure 37(a)) leads to lattice deformations and local strain [17, 19, 292, 298, 305] due to the steric repulsion between carbon p_z orbitals. Therefore, carbon atoms move apart from each other producing a corrugation that tracks the different kinds of stacking regions. As seen in figure 38, DFT calculations indicate that out-of-plane displacements are accompanied along with significant in-plane relaxation [19]. A vortex-like pattern can be observed in figure 38. The direction of circulation on the vortex reverts in between monolayers. Overall, the atomic relaxation results in the shrinking of the AA stacking regions in favor of the $AB/BA$ stacking domains [19].

Zoom In Zoom Out Reset image size

At very small twist angles ($\theta\sim 1^\circ$), the structural relaxation in TBG leads to the formation of triangular domains with Bernal stacking. Under an interlayer bias, these domains are gapped and the electronic transport occurs in one-dimensional networks formed by domain walls. Those states have been visualised in experiments recently [17, 299]. Transport in these one-dimensional networks has been studied using a model where the AA regions in TBG are represented as scattering centers connected by one-dimensional helical channels [306, 307]. Pseudo Landau levels and an effective honeycomb lattice for this model were also discussed [308].

The structural relaxation of AA stacking regions diminishes the tunneling in such regions. Theoretically such tunnelling can be set to zero resulting in a chiral model that contains the basic features of the magic angle physics [296, 297, 309, 310]. In fact, it can be proved that the squared chiral Hamiltonian is basically a Quantum Hall effect hamiltonian [310]. From the chiral model, one can obtain the magic angles for multilayer graphene [311].

The consequences of strain in the interacting phase diagram of TBG were investigated recently [312]: even for small strain ε (within the range of values observed in experiments), strain-induced topological phase transitions occur, supporting the possibility of strain playing an important role in sample-dependent observations. Along those lines, the role played by strain in the fabrication of moiré materials has been recently studied from a reproducibility perspective [281]. The presence of non-uniform strain in TBG might lead to variability in transport results within different regions of a single TBG sample. Models to determine twist angle and strain have been developed [313–315].

Small strains at nonzero integer fillings have been predicted to stabilise an 'incommensurate Kekulé spiral' order [316] (section 3.2.7).

In a study using an atomistic description of TBG [317], the role of a hBN substrate in the structural relaxation and band structure of TBG was shown to be important, introducing large effective masses and large pseudo-magnetic fields. Also, the long-range Coulomb electron–electron interaction affects the distribution of Berry curvature of the bands near charge neutrality of a TBG closely aligned with hBN. Due to the suppressed dispersion of the narrow bands, the band structure is strongly renormalised [54]. Therefore, the topological linear and nonlinear Hall conductivities in TBG/hBN are substantially modified by electron–electron interactions [54].

Figures 39(a)-(i–iv) contain possible relative displacements of the upper monolayer (in orange) with respect to the lower (blue) one. Those displacements, when folded back onto the unit cell, define an order parameter [315].

Zoom In Zoom Out Reset image size

'First order' domain walls obtained by dark field transmission electron microscopy are displayed in figure 39(b)-(i), while 'second order' images (which focus on the domain walls) can be seen in figures 39(b)-(ii–iv). The colors displayed on figures 39(b)-(ii–iii) are consistent with those displayed in figure 36. Figure 39(c) contains domain structures and domain walls. Note how the AA regions have collapsed onto smaller areas; the AA configurations correspond to the points in which two boundary lines in the unit cell in figure 39(d), which are avoided due to their high energy cost (see (red) maxima of the energy landscape at AA onfigurations in figure 33).

The topology of the unit cell can be equivalently observed in the unit cell (figures 39(a)(v), 39(d), and 40(a)) or in the torus displayed as figure 40(b). The crucial point is that the AA regions, being so energetically costly, collapse onto sections with a negligible area, and are explicitly removed ('cut out') in the topological description. Another equivalent description of the topology is obtained by joining the three (green, blue, and red) lines into points AB and BA, as it is done in figure 40(c) [315].

Zoom In Zoom Out Reset image size

To understand the topology of the relative displacements, figure 40(d) depicts a hexagonal circuit, whose center is an (avoided) AA configuration. The direction of a given displacement can be further encoded by the addition of 'inverse' operations, which reverse direction. This way, one works with those elements:

$$\begin{equation*} R,\,R^{-1},\,G,\,G^{-1},\, B,~\text{and }B^{-1}. \end{equation*}$$

An additional simplification calls for the following redefinition of group elements:

$$\begin{align} a\leftrightarrow RG^{-1},~b\leftrightarrow BR^{-1},~\text{and }ba\leftrightarrow BG^{-1}; \end{align} \tag{ 66 }$$

the definition of operations in terms of group generators a and b permits defining a vortex operation:

$$\begin{align} \left[a,b\right] = RG^{-1}BR^{-1}GB^{-1}, \end{align} \tag{ 67 }$$

which, starting at a configuration BA, goes clockwise (from left to right; figure 40(d)). An antivortex operation goes anticlockwise on figure 40 [315]:

$$\begin{align} \left[b,a\right] = BG^{-1}RB^{-1}GR^{-1}. \end{align} \tag{ 68 }$$

Figure 41 combines analyses of topology and geometry of a WSe₂/MoSe₂ moiré hetero-bilayer.

Zoom In Zoom Out Reset image size

A powerful method to describe the structure of twisted 2D materials relies on the cut and projection method used to generate quasicrystals [294, 318–320]. There, the projected structure is interpreted as the set of points of best fit between the two rotated structures [294].

A paradigmatic model for the electronic properties of TBG is described by a $2\times 2$ matrix operator. The starting point is the Bistritzer–MacDonald Hamiltonian [289], obtained by rotating two graphene low-energy effective Dirac Hamiltonians [1], and coupling the two graphene monolayers using the first harmonics of the interlayer interaction potential only [289, 321].

Considering as a basis set the wave vectors $\Phi(\mathbf{r}) = \begin{pmatrix} \psi_1(\mathbf{r}) , \psi_2(\mathbf{r}), \chi_1(\mathbf{r}), \chi_2(\mathbf{r}) \end{pmatrix}^T$, where the index 1 or 2 labels a graphene monolayer and $\psi_j(r)$ and $\chi_j(r)$ are the two Wannier orbitals on each inequivalent site (sublattice) of the graphene's unit cell [1], the chiral Hamiltonian is given by [309, 311]:

$$\begin{align} \mathcal{H} = \begin{pmatrix} 0 & D^{\ast}\left(-\mathbf{r}\right)\\ D\left(\mathbf{r}\right) & 0 \end{pmatrix}. \end{align} \tag{ 69 }$$

Zero-mode operators are defined as:

$$\begin{align} D\left(\mathbf{r}\right) = \begin{pmatrix} -i\bar{\partial} & \alpha U\left(\mathbf{r}\right)\\ \alpha U\left(-\mathbf{r}\right) & -i\bar{\partial} \end{pmatrix}, \end{align} \tag{ 70 }$$

and

$$\begin{align} D^{*}\left(-\mathbf{r}\right) = \begin{pmatrix} -i\partial & \alpha U^{*}\left(-\mathbf{r}\right)\\ \alpha U^{*}\left(\mathbf{r}\right) & -i\partial \end{pmatrix}, \end{align} \tag{ 71 }$$

where $\bar{\partial} = \partial_x+i\partial_y$ is the antiholonomic differential operator and $\partial = \partial_x-i\partial_y$ the holonomic one. The dimensionless coupling (potential) is given by:

$$\begin{align} U\left(\mathbf{r}\right) = e^{-i\mathbf{q}_1\cdot \mathbf{r}}+e^{i\phi}e^{-i\mathbf{q}_2\cdot \mathbf{r}}+e^{-i\phi}e^{-i\mathbf{q}_3\cdot \boldsymbol{r}} \end{align} \tag{ 72 }$$

where $\phi = 2\pi/3$ and the moiré reciprocal lattice vectors are given by $\mathbf{q}_{1} = k_{\theta}(0,-1)$, $\mathbf{q}_{2} = k_{\theta}(\frac{\sqrt{3}}{2},\frac{1}{2})$, $\mathbf{q}_{3} = k_{\theta}(-\frac{\sqrt{3}}{2},\frac{1}{2})$ (see figure 37(b)).

The moiré modulation vector is

$$\begin{align} k_{\theta} = 2k_{D}\sin{\frac{\theta}{2}}, \end{align} \tag{ 73 }$$

where θ is the twist angle between layers and $k_{D} = \frac{4\pi}{3a}$ is the Dirac wave vector norm. The angle and interlayer coupling effects are both captured by the parameter α, defined as:

$$\begin{align} \alpha = \frac{w_1}{v_F k_\theta}, \end{align} \tag{ 74 }$$

with $w_1 = 110$ meV the interlayer coupling for AB stacking and $v_F = \frac{19.81\,\textrm{eV}}{2k_D}$ is the Fermi velocity.

The spectrum can be found by proposing a Bloch wave function with momentum k in the moiré first Brillouin zone:

$$\begin{align} \begin{pmatrix} \psi_{\mathbf{k},1}\left(\mathbf{r}\right)\\ \psi_{\mathbf{k},2}\left(\mathbf{r}\right) \end{pmatrix} = \sum_{mn} \begin{pmatrix} a_{mn}\\ b_{mn}e^{i\mathbf{q}_1\cdot\mathbf{r}} \end{pmatrix}e^{i\left(\mathbf{K}_{mn}+\mathbf{k}\right)\cdot\mathbf{r}}, \end{align} \tag{ 75 }$$

and

$$\begin{align} \begin{pmatrix} \chi_{\mathbf{k},1}\left(\mathbf{r}\right)\\ \chi_{\mathbf{k},2}\left(\mathbf{r}\right) \end{pmatrix} = \sum_{mn}\begin{pmatrix} c_{mn}\\ d_{mn}e^{i\mathbf{q}_1\cdot\mathbf{r}} \end{pmatrix}e^{i\left(\mathbf{K}_{mn}+\mathbf{k}\right)\cdot\mathbf{r}}, \end{align} \tag{ 76 }$$

and solving for the coefficients $a_{m,n}$, $b_{m,n}$, $c_{m,n}$, and $d_{m,n}$ in the Hamiltonian (69).

Figure 42 presents the resulting squared energies E² as a function of α obtained at three k-points: Γ-, K, and $\mathbf{k} = \boldsymbol{q}_1+0.14\mathbf{q}_2+0.23\mathbf{q}_3$, the last one being a k-point of low symmetry. Working with E² instead of E helps emphasise the magnitude of the energy–given that the spectrum is electron–hole symmetric. The band edge is located at the $\Gamma-$point and the dips in figure 42 correspond to the magic angles $\alpha = \alpha(\theta)$, because the band shrinks to zero width for those values of the twist angle θ (see equations (73) and (74)).

Zoom In Zoom Out Reset image size

Chiral symmetry induces an electron–hole symmetry [322–325] that leads into an effective $2 \times 2$ Hamiltonian. Squaring $\mathcal{H}$, it decouples into two identical $2 \times 2$ matrices given by [296]:

$$\begin{align} H^{2} = \begin{pmatrix} -\nabla^{2}+\alpha^{2} |U\left(\mathbf{-r}\right)|^{2} & \alpha F^{\dagger}\left(\mathbf{r}\right)\\ \alpha F\left(\mathbf{r}\right) & -\nabla^{2}+\alpha^{2} |U\left(\mathbf{r}\right)|^{2} \end{pmatrix}. \end{align} \tag{ 77 }$$

The squared norm of $U(\mathbf{r})$ can be understood as an effective confinement potential in the moiré supercell:

$$\begin{align} |U\left(\mathbf{r}\right)|^{2} & = 3+2 \cos\left(\mathbf{b}_1\cdot \mathbf{r}-\phi\right)\nonumber\\ & \quad +2\cos\left(\mathbf{b}_2\cdot \mathbf{r}+\phi \right)+2\cos\left(\mathbf{b}_3\cdot \mathbf{r}+2\phi\right). \end{align} \tag{ 78 }$$

As seen in figure 43(a), $|U(\mathbf{r})|^2$ resembles a Kagome lattice.

The off-diagonal terms of the $2\times 2$ matrix are:

$$\begin{align} F\left(\mathbf{r}\right) = -i\sum_{\mu = 1}^{3}e^{i\mathbf{q}_{\mu}\cdot\mathbf{r}}\left(2\hat{\mathbf{q}}_{\mu}^{\perp}\cdot \mathbf{\nabla}-1\right), \end{align} \tag{ 79 }$$

and,

$$\begin{align} F^{\dagger}\left(\mathbf{r}\right) = -i\sum_{\mu = 1}^{3}e^{-i\mathbf{q}_{\mu}\cdot\mathbf{r}}\left(2\hat{\mathbf{q}}_{\mu}^{\perp}\cdot \mathbf{\nabla}+1\right), \end{align} \tag{ 80 }$$

where $\mathbf{\nabla}^{\dagger} = -\mathbf{\nabla}$ with $\mathbf{\nabla} = (\partial_x,\partial_y)$ and $\mu = 1,2,3$. This is an essential point as eigenvalues must be real (notice that $-F^{\dagger}(\mathbf{-r}) = F(\mathbf{r})$). $\hat{\mathbf{q}}_{\mu}^{\perp}$ is a set of unitary vectors perpendicular to the set $\mathbf{q}_{\mu}$ such that $\hat{\mathbf{q}}_{\mu} \cdot \hat{\mathbf{q}}_{\mu}^{\perp} = 0$.

The contributions from $F^{\dagger}(\mathbf{r})$ and $F(\mathbf{r})$ can be condensed into a new operator

$$\begin{align} \hat{F}\left(\mathbf{r}\right) = \begin{pmatrix} 0 & F^{\dagger}\left(\mathbf{r}\right)\\ F\left(\mathbf{r}\right) & 0 \end{pmatrix} \end{align} \tag{ 81 }$$

proportional to an interlayer current between different bipartite lattices [297] and, eventually, with the angular momentum in the perpendicular direction to the two graphene monolayers: for small twist angles, the squared TBG chiral model Hamiltonian turns out to be akin to a Quantum Hall effect Hamiltonian [310, 326].

The renormalisation simplifies the Hamiltonian considerably and highlights the three main ingredients of charge carrier dynamics on twisted few-layers: (i) a kinetic energy contribution via the $\nabla^{2}$ term, (ii) a confinement, kind of quantum dot, potential energy $|U(\mathbf{r})|^{2}$, and (iii) the interaction operator $\hat{F}(\mathbf{r})$, which is an interlayer current between bipartite lattices that eventually results in a coupling of angular momentum and an effective magnetic field [326].

The expectation values from those three energy terms are presented as a function of the (twist-angle dependent) parameter α at the Γ point in figure 43(b). The kinetic and confinement potential contributions are always positive while the angular contribution is always negative. The first magic angle is different from others, but reachable through perturbation theory [297]. As seen in figure 43(c), the flat band can be understood as a very special condition in which the sum of kinetic and confinement energy is equal to the angular effective confinement contribution achieved through destructive interference paths [310]. Another interesting feature can be seen in figure 43(c): the magic angles follow a quantised rule $\alpha(\theta_{m+1})-\alpha(\theta_m) = 3/2$, where m is the order of the twist angle $\theta$ [310]. The wave functions at magic angles are nearly coherent Landau levels and thus have many important properties and applications [310].

The renormalisation produced by the squaring procedure maps the flat band into a ground state separated by a gap from the rest of the states. The ground state has an antibonding nature in a triangular lattice and then has frustration, associated with a massive degeneration [98, 296, 322–324]. Also, this renormalisation is equivalent to a supersymmetric transformation from a fermion to a bosonic Hamiltonian, a subject of intense research in solid state physics [327, 328].

An intriguing and yet unexplored feature is the equivalence of the flat band with floppy modes on elastic Hamiltonians [296] for mechanical systems with constraints, such as rods or bars between hinges [329–331], or glasses [332, 333].

The chirality of TBG seems to induce spin textures around the K point of the Brillouin zone with alternating vortex-like textures [334]. The helicity of each vortex is inverted by changing the chirality. The spin texture changes as the twist angle evolves.

A unified picture of the quantum phase diagram of TBG in terms of incommensurate Kekulé spiral (IKS) order was proposed recently [316, 335]: this is a time-reversal-symmetric and spatially nonuniform order, which shifts the spatial coordinates and rotates the U(1) angle characterising the spontaneous intervalley coherency simultaneously [316]. Also, many different flavors of fermions can be obtained by making heterostructures of Kekulé over another Kekulé-patterned graphene [336].

Intercalating graphene layers with alkali atoms leads to effects similar to strain [337]. Although those results seem to be relevant for the alkali-intercalated systems where Kekulé patterings have been reported [158, 176], the relation between intercalation and strain has not been discussed in the literature to the best of our knowledge.

Measurements of the second order transversal conductivity have been used to detect topological transitions in strained moiré systems [338], as the latter are closely related to the changes in sign of the Berry dipole. In strained TBG, it has been observed that the nonlinear effects can be even more important than the linear ones [339]. Similar nonlinear effects have been seen in rotated bilayers of WSe₂ [340], where a change in the sign of the Berry dipole across the gap has been identified. The effect is enhanced due to the large spin–orbit gap typical of TMDCs [341]. Apart from charge transport, thermal and thermoelectrical transport have also been investigated in strained TBG [342].

At the 2D limit, the electric field effect can penetrate a bilayer, and ferroelectric metals exist in WTe₂ bilayers [343]. Here, one continues the discussion started in section 4.1 concerning ferroelectric insulating bilayers though. More specifically, figures 34(c) and (d) indicate that the (insulating) hBN bilayer switches its direction of polarisation P as it changes from the AB to the BA configuration. This means that, when a a moiré structure is set with a small angle δ away from the configurations drawn in figures 33(e), (g) and (h), a domain structure like the one shown in figure 42(a) will set up. Additionally, if the types of strain shown in figures 42(b) and (c) are also taking place, those will also give rise to an antiferroelectric disposition that is experimentally depicted by piezo-force microscopy in figure 44. The effect has been experimentally generalised to moiré TMDC bilayers as well [283, 300–302]. That the intrinsic polarisation P increases in discrete steps with additional layers has been demonstrated as well [280].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Flat bands and strong electronic correlations are not exclusive of magic-angle moiré graphene bilayers. As seen in figure 45, charge localisation extends to multiple moiré homo-bilayers such as the hBN one [291, 344–346]. Strong electron–electron correlations on flat-band non-graphene bilayers are being vigorously studied now [347].

Zoom In Zoom Out Reset image size