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

This is an update of a previous review (Naumis et al 2017 Rep. Prog. Phys. 80 096501). Experimental and theoretical advances for straining graphene and other metallic, insulating, ferroelectric, ferroelastic, ferromagnetic and multiferroic 2D materials were considered. We surveyed (i) methods to induce valley and sublattice polarisation (P) in graphene, (ii) time-dependent strain and its impact on graphene’s electronic properties, (iii) the role of local and global strain on superconductivity and other highly correlated and/or topological phases of graphene, (iv) inducing polarisation P on hexagonal boron nitride monolayers via strain, (v) modifying the optoelectronic properties of transition metal dichalcogenide monolayers through strain, (vi) ferroic 2D materials with intrinsic elastic (σ), electric (P) and magnetic (M) polarisation under strain, as well as incipient 2D multiferroics and (vii) moiré bilayers exhibiting flat electronic bands and exotic quantum phase diagrams, and other bilayer or few-layer systems exhibiting ferroic orders tunable by rotations and shear strain. The update features the experimental realisations of a tunable two-dimensional Quantum Spin Hall effect in germanene, of elemental 2D ferroelectric bismuth, and 2D multiferroic NiI2. The document was structured for a discussion of effects taking place in monolayers first, followed by discussions concerning bilayers and few-layers, and it represents an up-to-date overview of exciting and newest developments on the fast-paced field of 2D materials.


Motivation and Update's structure
One Review of strain and the electronic and optical properties of graphene and other 2D materials appeared six years ago [1].Exciting advances since then include (1) strain on lateral superlattices; (2) non-linear Hall effects, Berry dipoles, and exotic quantum phase diagrams due to strain in graphene; (3) the experimental verification of piezoelectricity in hexagonal boron nitride (hBN) monolayers; (4) the creation of strain by a sudden thermal quenching during 2D material growth and the effects of strain on the optical properties of transition metal dichalcogenide (TMDC) monolayers; (5) the role of strain on 2D phase transitions in 2D ferroelectrics and ferromagnets; (6) twisted (moiré) 2D materials exhibiting flat electronic bands; (7) the role of local and global (and even time-dependent) strain on superconductivity and topological phases of twisted few-layer graphene; (8) magnetoelectric couplings on 2D magnetic bilayers; among others.Previous list of topics contains strained monolayers (topics 1 through 5) and few layers (topics 6 through 8) under strain.Strain is an integral part of the field of 2D materials, and recently synthesised 2D materials such as (Quantum Spin Hall) germanene, (elemental ferroelectric) bismuth monolayers, and (magnetoelectric multiferroic) NiI 2 monolayers-all adding functionalities and additional 2D material platforms-are featured here as well.
The Update is written to reflect that order: section 2 lists experimental techniques to strain 2D materials, section 3 contains a discussion of strained graphene, piezoelectric hBN monolayers, strained TMDC monolayers, ferroelectric and magnetic 2D materials under strain, and section 4 is devoted to strained bilayers and multilayers (including moirés, ferroelectric and antiferroelectric bilayers, as well as engineered multiferroic bilayers).Recently synthesised 2D materials appear within those sections.Conclusions are presented in section 5.
Prior to undertaking the main task of this Update, it may be necessary to reassess the physical meaning of mechanical strain and to disclose the definition followed here.Before the atomistic hypothesis was born, strain was a measure of mechanical stress and of the elasticity of bulk materials at finite temperature, without concern for their microstructure.In this sense, strain is an out-of-equilibrium mechanical modification of a material.Examples of this definition include topics (i) through (vi) in the list above.
But, sometimes strain is also thought of as a modification of a 'pristine' (crystalline) zero-temperature lattice.In other words, though a lattice may be in mechanical and thermal equilibrium (and hence unstrained in the strictest sense of the word), its lattice parameters or atomic positions may have changed with respect to those in the pristine scenario.Examples of that definition of (permanent-or quasi-permanent in case a nano-bubble were to explode, or temperature changed through a 2D phase transition) strain include processes (vii) through (xiii).
Both uses of the word strain are employed by the community, and it is important to make the point that those definitions are not identical.
The effect of strain on materials' properties is studied in what follows.

Strained monolayers
Strain affects the electro-optical properties of 2D materials due to the stretching, bending and torsion of chemical bonds.Stretching has the largest effect because it alters electronic hopping between atoms the most [39].

Updates on graphene strain engineering
The most common theoretical approaches for graphene's strain engineering are tight-binding (TB) methods, the low energy approximation in a continuum or within a discrete lattice, and density-functional theory (DFT) calculations [1].
Aiming for continuity of the theoretical descriptions, the notation and definitions used in [1] will be maintained.For example, graphene's zero-temperature lattice vectors are still given by: ) , and a 2 = a ) , (1) where a = 1.42 Å is the distance between carbon atoms.Standard notation for atomic displacements induced by strain is kept unchanged as well, so that atomic positions under deformation are given by this field: with in-plane displacements given by u = u(r) = (u x (r), u y (r)) and out-of-plane displacements denoted by h = h(r), where r = (x, y) is a 2D position vector.The strain tensor ε (with components ε µν ) is related to the deformation field as follows: (3)

Correcting second-and third-nearest neighbour interactions on empirical tight-binding
Hamiltonians.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]: where r ′ runs over all sites of the deformed honeycomb lattice and the hopping integral t r ′ ,d ′ n (r) varies due to changing interatomic distances.The operators c † r ′ and c r ′ +d ′ n (r) create an electron at r ′ and annihilate another at r ′ + d ′ n (r).Vectors d ′ n (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.

Corrections of empirical TB Hamiltonians due to spinorbit coupling.
In-plane strain leads to quadratic corrections on spin-orbit coupling [42].Out of plane deformations, on the other hand, hybridise π-orbitals with σ−electrons and produce a first-order contribution in the spin-orbit interaction strength [43] with three contributions [42]: The colour map corresponds to the angle θ, which varied from 0 to π. Second-and third-nearest neighbour interactions require a larger relative correction when compared to first-nearest neighbour interactions.Reprinted with permission from [40].
where the labels A 1 , B 2 , and G ′ represent the irreducible representations of the group C ′ ′ 6 , which results from considering a √ 3 × √ 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]).
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 π − σ 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]: where the position-dependent Fermi velocity tensor v(r) is given by: with v F the Fermi velocity (v F /c ≈ 1/300, where c is the speed of light in vacuum), σ = (ησ x , σ y ) is the Pauli matrix vector and η = ±1 labels the valley, and σ 0 is the 2 × 2 identity matrix.In turn, the (deformation) scalar and (pseudomagnetic) vector potentials are given by: 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 Γ s (r) (l = x, y, z) is: Unlike A s , Γ 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 ψ = Mψ to retain hermiticity, and for self-adjoint extensions to preserve currents [72][73][74][75].M is a matrix containing microscopic details and symmetries, and ψ is the electron/hole wavefunction at the boundary.Here, Γ s (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 nonhomogeneous 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]: where g µν is the metric tensor associated with strain, δ ai is the Kronecker delta, √ g the square root of the determinant of the metric tensor, and Tr(g) = ∑ µ g µµ 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 = ((∆R)/R)/ε xx , where R was the resistance and ∆R the change of R due to strain [77].Modifications of optical properties under uniaxial strain are detailed elsewhere [78,79].

Electrostatic effects due to substrates.
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 sub (r)σ 0 [47] and, if the substrate produces a bipartite symmetry breaking, an extra term ∆V(r)σ z = (V A (r) − V B (r))σ z [44], which measures the difference of the electrostatic potential in the A and B sublattices [81].Localised electronic states at folds.From [83].Reprinted with permission from AAAS.

Cut and folded graphene: sublattice and valley polarisation.
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 dI/dV measurements taken at locations marked by circles in figures 2(e) and (f) [85].The splitting of dI/dV peaks (and their weight redistribution) in figure 2(g) are manifestations of localised states on folded graphene [82,83].
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][89][90] were theoretically studied using continuous models with effective pseudomagnetic fields and TB Reprinted with permission from [84].Copyright (2018) American Chemical Society.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 valleypolarised currents [84], opening an avenue for graphene folds as straintronic quantum wires [92,93].
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. u(r) = 0. From equations (3) and ( 9), one gets: 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: 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, A s (r) has only one non-zero component [95]: Equation (17) implies that ∇ • A s (r) = 0, so that A s (r) is in the Coulomb gauge.Thus it can be derived from a pseudoscalar magnetic potential: with i, j = x, y and ε ij a 2D Levi-Civita tensor.Expressing Φ(y) as a Fourier series: with and The associated pseudomagnetic field becomes: Φ 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 A s (r) is derived from Φ(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(r).Using equation ( 6) for E = 0, two wave functions are obtained: 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 ∆ A .The average of the coefficients in equation ( 19) can be expressed as: 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]: from where it follows that [97]: with 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][102][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 : where E is the electron's energy, v F the Fermi velocity, Ψ A (Ψ 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.
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][116][117][118]. Depicted in figure 5, one method involves the combination of straininduced 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, Hexagons and the two sublattices were overlaid.The LDOS is higher in the sublattice corresponding to the red dots at higher currents.Reprinted with permission from [101].Copyright (2017) American Chemical Society.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 ± : Thus, energy levels of the K + valley are pushed up, while those of the K − valley are pushed down.This effect was observed in rippled graphene grown on rhodium foil [112,[119][120][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 B R = 0 T (figure 5(c)) show pseudo-magnetic Landau levels only at points where B s ̸ = 0.The (n ̸ = 0) 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.

Graphene superlattices
The focus on this section are periodic patterns on graphene.

Strain-modulated superlattices.
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 B s (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 ∼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).
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 ∝ √ n, but to a quantisation appearing linear in energy (E n ∝ 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 = √ h/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  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: 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 ∼ 100 T) that can be experimentally created: graphene ripples might provide a platform for valley-dependent electron transport.
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.Estimations of strain with a STM are usually derived from dI/dV 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 dI/dV 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.
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 2 nanospheres [124], or over black phosphorus substrates [128].

Second harmonic generation on sublattice-polarised graphene.
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]: where Ē = (E s,A + E s,B )/2 and ∆ E = (E s,A − E s,B )/2.Second-harmonic generation has also been experimentally employed to sense strain in MoS 2 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.

Flat-band superconductivity in strained graphene.
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 ρ(E F ), that theory indicates a critical temperature depending exponentially on the coupling strength λ, T c ∼ e −1/(|λ|ρ(EF)) [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 ∼ λ [132][133][134] for large ρ(E F ).

Review
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][134][135][136][137][138][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 0 (see equation ( 9)): which can be realised by the following in-plane displacement fields: (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 = |∇ × A s | is shown in figures 11(c) and (d).(The pseudomagnetic field strength (B s ≈ 100 T) and period (d ≈ 14 nm) visualised in experiments [142] would lead to u 0 ≈ 5 for the displacement fields defined here.) 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 ∆ A/B is sublatticedependent.∆ 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 A 1D cos 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.
The SC order parameter ∆ exceeds the flat-band bandwidth for high enough values of the strain strength u 0 and the interaction strength λ, and its dependence is linear in λ [132][133][134].At lower values of u 0 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 ∆ ∼ λ occurs) and it would require increasing the strain strength u 0 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 spatiallydependent antiferromagnetic coupling [138], showing that a chiral, d-wave SC can be stabilised via strain as well.

Flat electronic bands and electron-electron correlations in buckled graphene.
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 2 [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 dI/dV spectra exhibited peaks that were identified as the pseudo-Landau levels induced by a pseudomagnetic field B s (E n ∼ sgn(n) √ |n|B s for n = 0, ±1, ±2, . ..), 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).
The dI/dV 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]: (37) where B is the amplitude of the pseudo-magnetic field, a b the superlattice period, b The choice of coordinates is consistent with C 6 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 ∆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]: at which ∆E = E I for the nth band.Φ 0 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]: while strain was introduced by modifying the hopping energies connecting the three nearest neighbours t n = t + δt n (n = 1, 2, 3) with 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 ≈ 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 ∆ ≈ 0.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 ∆E smaller than the Hubbard interaction U (∆E/U ∼ 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 µ(r) = µ 0 ∑ 3 i =1 cos(b i • 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 eff (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 minivalleys 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 eff (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 spinpolarised 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][151][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.

Adatom-induced graphene superlattices.
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].

Short-wavelength periodic lattice deformations: Kekulé and other
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 Liintercalated graphene (Kek-O) [158].As depicted in figure 14 by bonds of two different strengths, Kekulé patterns have few hopping strengths changed within a √ 3 × √ 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].
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.
One sets graphene onto a Hubbard model Hamiltonian with an additional elastic energy term [161]: The c x,σ on equation ( 41) are fermion annihilation operators at lattice site x with spin σ.The occupation number of such site is n x = c † x,σ c x,σ .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, reflectionsymmetric distortions are allowed on the √ 3 × √ 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 ∼27-31 N m −1 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].
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.√ 3)K ± [170].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]: where âr and br are annihilation operators at site r for a carbon atom in the A or the B sublattice, respectively.Operators â † r and b † r are the corresponding creation counterparts.A bonddensity wave modulates the TB parameters and originates the Kekulé distortion: the first-neighbour hopping is given by [159]: where t 0 = 2.8 eV is the hopping parameter for first-nearest neighbours in pristine graphene and G = K + − K − is a reciprocal lattice vector coupling both valleys (see figure 17(b)).
The amount of the bond modulation was ∆ ≈ 0.1.In addition, p and q are two integers such that the index: distinguishes the Kekulé O pattern (ν = 0) from the two Y patterns (ν = ±1) alluded to in figure 14.
The subsequent, second-nearest neighbour parameters were periodically modulated as: where t 2 = 0.1t 0 is the second-nearest neighbour TB parameter of pristine graphene.A column vector is now defined: that collects annihilation operators at k and k + G, so that the Hamiltonian can be written as a 6 × 6 matrix: made from the original graphene Hamiltonian at the Γ−point: the Hamiltonians at the G and −G points: and the interaction between them: where The spectrum can be found by diagonalisation of the 6 × 6 operator or by a projection technique that permits obtaining a low-energy 4 × 4 matrix effective Hamiltonian.That effective Hamiltonian becomes for the Kek-Y case [172]: where σ = (σ x , σ y , σ z ) and τ = (τ x , τ y , τ z ) are two vectors with components defined from Pauli matrices, σ 0 and τ 0 are two 2 × 2 identity matrices, and µ is defined as: The resulting four low-energy bands turn out to be: and were two new velocities have been defined: As indicated in figure 17(b), the Kek-Y pattern folds the two valleys K + and 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 Liintercalated 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 K ± points to K ′ ± 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 K ± points of the P6/mmm space group are replaced by the F 0 and ∆ points in the Cmmm space group when under uniaxial strain, and the Dirac points K ′ ± D of the strained lattice do not necessarily coincide neither with K ′ ± , nor with F 0 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 K+ D and 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.

Time-dependent strain.
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 timedependent 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 3 , which itself lays on an insulating SiO 2 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.
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.SAWs can also provide alternative ways to study topological properties [178][179][180] and the topological nonadiabaticity 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][182][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].

Berry dipole and nonlinear transport from strain.
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 ´dkf 0 ∂ kx Ω (the so-called Berry dipole, with f 0 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].

Experimental verification of a tunable electronic topology in germanene
The discussion continues with germanene, another monoelemental 2D material with a honeycomb lattice, because of the very recent experimental verification of its tunable topological properties [153].Its availability opens possibilities for topological field-effect transistors whereby a back or top gate turns dissipation-less transport along its edge on and off [148].
There are advantages in growing germanium on Ge 2 Pt(101) though [153].The most salient one is the creation of a buffer layer reminiscent of that created in epitaxial graphene grown on silicon carbide [196], which renders a graphene-like dispersion on the next subsequent carbon layer.Here, a second layer above appears electronically and chemically decoupled, thus featuring the electronic and topological properties predicted in freestanding germanene [191].In particular, a spin-orbitcoupling-induced band gap of 23.9 meV, which is sufficiently large to observe the Quantum Spin Hall effect.Experimentally, the observed gap was much larger at 70 meV; this discrepancy may be due to strain, buckling [43], stacking sequence, and effects due to the proximity to the substrate.

Piezoelectricity by strain in hexagonal boron nitride monolayers
An essential feature of noncentrosymmetric crystals is piezoelectricity, which permits exchanging strain and intrinsic (internal) electric fields.The induced charge density ρ couples to strain through variations in the intrinsic polarisation P through [197]: The ith-component of the intrinsic polarisation P couples to the crystal strain through: where γ ijk is the 3rd rank piezoelectric tensor, whose only non-zero entry is γ ≡ γ yyy = −γ yxx = −γ xyx = −γ xxy for 2D materials with D 3h symmetry.This way, it can be shown that P(r) = γA(r) × ẑ, where A(r) = (u xx − u yy )x − 2u xy ŷ, with x = (1, 0, 0), ŷ = (0, 1, 0), and ẑ = (0, 0, 1).The piezoelectric coefficient is estimated to be 2.91 × 10 −10 C m −1 using the modern theory of polarisation, working with a twoband model of a massive Dirac equation with a 5.97 eV band gap [35].Electrical measurements performed by electrostatic force microscopy on top of bubbles and creases of hBN monolayers (figure 19) show a net electric field, which in turn represents another method to quantify local strain [35].(Electrostatic force microscopy is a non-contact technique that maps the local electrostatic interaction a tip and a sample.)No such contrast was observed on hBN bilayers nor in bulk samples, which recover a centre of inversion despite strain.

Quantifying strain on transition metal dichalcogenide (TMDC) monolayers
TMDCs transition from an indirect onto a direct band gap at monolayer thickness.They display multiple electro-optical properties [198], including strong exciton binding energies and photoluminiscence, that can be tuned by strain.
Quantitatively assessing the degree of strain, nevertheless, is not a trivial task.The difficulty stems from the fact that x-ray diffraction (XRD)-the most powerful technique to measure precise lattice constants-is not readily applicable to monolayers not having a macroscopically uniform in-plane orientation.This is because the incident x-ray will spread over a macroscopic distance (of the order of millimeters) for a grazing-incidence XRD geometry required to probe any in-plane strain (figure 20), thus leading to a small averaged signal for a sample with small size and/or random orientation.
Substrates can impose a uniform strain in the film grown on the surface.In standard epitaxial growth, substrates are chosen so that their surfaces have the same in-plane lattice symmetry as the film to be grown.In this case, strain is caused by the difference in the in-plane (super)lattice spacing of the substrate and the film.Synchrotron-based grazing incidence XRD was used to measure such strain in a WSe 2 monolayer epitaxially grown on graphene, where a lattice compression of −0.19% was found in WSe 2 [32] (figure 20).The 3 × 3 unit cell of the slightly compressed WSe 2 was perfectly commensurate with a 4 × 4 graphene lattice with a experimental mismatch below 0.03%, which could explain why the WSe 2 monolayer was compressed on graphene.
As indicated in section 3.2, non-uniform strain may be realised by patterned structures manufactured on substrates.For the sake of quantitative evaluation, XRD can be applied more accurately under uniform strain conditions.Non-uniform strain, on the other hand, can be characterised with local probes such as Raman spectroscopy.In that case, the expected shift of phonon frequencies as a function of strain should be available, at least theoretically, to convert observed shifts in Raman modes into a magnitude of strain.As seen in figure 21 and similar to graphene (section 3.2.3),non-uniform strain can be measured through second harmonic generation when inversion symmetry is removed by strain [130] as well.
Using methods such as local heating, strain induced on amorphous substrates by a thermal expansion mismatch has been studied [207].CVD-grown MoS 2 monolayer is also seen to be strained after the cool down from growth at high temperatures [208].The (per cent) strain induced in the TMDC monolayer in the latter scenario can be estimated by: where α TMD is the thermal expansion coefficient of a TMDC monolayer, α sub is that of the substrate, ∆T = T growth − T room (where T growth and T room denote the growth temperature and room temperature, respectively).Recently, strain in TMDC monolayers was enhanced by utilising a high growth temperature in PVD (T growth = 1200 • C) and a systematic change in strain levels for Si 3 N 4 and SiO 2 substrates, consistent with the above scenario, was observed by using PL and Raman spectroscopy.The resultant strain caused large shifts in the excitonic transition energy of ∼80 meV for MoS 2 monolayers on SiO 2 (figure 22).Similar to graphene and hBN monolayers, TMDC monolayers can also be strained by polymer encapsulation [15,16], dielectric nanopillars [29,31,209,210], nanobubbles [36][37][38], and nanoindentation [10].Nevertheless, and in contrast with graphene, metal nanogaps [211], metal nanostructures [212], and optical waveguides [213][214][215][216] are also employed to create strain and/or to enhance their optical response.An additional pathway for strain and the tuning of optical properties is the creation of lateral heterostructures [20][21][22][23][24].

Tuning optical properties of TMDC monolayers by strain.
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  [203,205] and theory [206], while filled squares denote those for PVD-grown monolayers on Si 3 N 4 (blue square) and SiO 2 (orange square) in which strain is due to thermal expansion mismatch induced after the high-temperature growth.Grey lines are linear fits to data from strained PVD monolayers.
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 2 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 2 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 2 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 2 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][224][225][226][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 spinflips described thus far, and/or momentum transfer [228]-to occur.
In the WSe 2 monolayers studied by Hasz et al, the dark exciton has a lower energy (∼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 2 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 valleydependent 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 2 monolayer lies flat there, providing a homogeneous local dielectric environment.
Because of funneling, dark KΛ 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Λ excitons scatter into bright KK states, which turn visible in the photoluminiscence spectra [30].An additional study showed that excitons created on a MoS 2 monolayer do not feature anti-funnelling, thus emphasising the crucial role of the electronic structure of a WS 2 monolayer (having almost degenerate conduction band edges at the K and Λ k-points) to achieve this effect.
Localised strain was created in encapsulated WSe 2 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 2 monolayers [235].
3.6.Ferroic behaviour and strain in monolayers with degenerate structural configurations 3.6.1.Context.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.

Structural degeneracies and ferroicity of some 2D
materials.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][243][244][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.
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 ∆ 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][238][239].The structural transformation of the SnSe monolayer was expressed in terms of the parameter ∆α, 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.

2D ferroelectrics within group-IV monochalcogenide
monolayers.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].

Two-dimensional phase transformations in 2D ferroelectric and ferroelastic materials driven by shear strain.
Honeycomb lattices are prevalent on many 2D materials such as graphene, hBN monolayers, silicene, germanene, and TMDC monolayers with 1H and 1T atomistic configurations.
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 ∆α (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 ∆α = 0. Some materials with rectangular unit cells undergo twodimensional structural phase transformations at finite temperature, whereby the 2D unit cell suddenly changes from a rectangle onto a square: figure 25(b) depicts ∆α 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: P = 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]: so a 1 = a 2 = a C when ∆α = 0, and P(a C , a C ) = (0, 0, 0) (figure 25).

Elastic properties of 2D ferroelectrics through 2D structural transformations.
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 0 ): Nevertheless, the sudden compression of the lattice parameter a 1 and the elongation of a 2 at T c driving the collapse of ∆α 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: where ⟨a 0 (T)⟩ 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 ⟨f(U max )⟩ as an average over classically accessible states [26]: with dXdY an area element within the confines of an isoenergy contour U 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 max ⩾ J nevertheless, the average structure encompasses minima A and B in figure 24(c)-(ii), yielding ⟨a 1 ⟩ = ⟨a 2 ⟩, 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].
Thermal averages of elastic constants C ij are then determined by [260] (with k B Boltzmann's constant, A = a 1 a 2 , and u = U/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 ⟨C 11 ⟩ and ⟨C 22 ⟩ 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 ⟨C 11 ⟩ = ⟨C 22 ⟩ for U max > J dictates the redefinition of equilibrium lattice parameters evolving with temperature (equation ( 56)), and hence of a thermally-evolving strain.
⟨C 12 ⟩ is the softest elastic modulus and ⟨C 11 ⟩ hardens at the transition (U max = J), while ⟨C 22 ⟩ softens at U max = J.According to figure 26 By direct comparison J and T C from numerical calculations [239], a correspondence T ∝ 1.42U max is established, such that T C = 212 K, and finite−T elastic behaviour can be extracted from figure 26(b).
3.6.6.Photostriction: anisotropic, non-thermal strain created by illumination.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).
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 × n k k-point grid and considering spinorbit coupling, one gets a density of charge carriers n c = 1/(n 2 k A 0 ), where A 0 is the area of the unit cell before the photoexcitation.Setting n k = 41, a photoexcitation around the Experimental demonstration of photostriction on orthorhombic, bulk GeS.The sample was excited by a 400 nm femtosecond laser with a photon energy above GeS' bandgap and probed by ∼100 fs ultrashort electron pulses at normal incidence.The reversible photoinduced mechanical deformation of the lattice was obtained from the shift of Bragg peak center positions in reciprocal space (Q).The {020} AC ({200} ZZ ) direction corresponds to the x (y) direction in figure 27(a).As seen in figure 27(e), an anisotropic strain in which the x direction compresses and the y direction elongates is experimentally verified.Reprinted with permission from [8].Copyright (2023) American Chemical Society.
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 > a 2 , one observes in figure 27(e) that a 2 increases in magnitude while a 1 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 +x onto the −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 2 and a contraction of a 1 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.

Using uniaxial strain to tune the critical temperature of 2D ferroelectrics.
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.
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 1 increases the value of ∆α-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 ∆α = 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 α 1 and α 3 , and of distances d 2 and d 3 .T c increases regardless of the direction in which uniaxial tensile strain is applied.(Symbols ⟨•⟩ 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.
3.6.8.2D bismuth: an elemental ferroelectric monolayer.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 nontrivial 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 2 6p 3 ), one may first imagine that its 2D structural configuration may resemble that of phosphorene (electronic configuration 3s 2 3p 3 ), a nonferroelectric 2D material possessing a center of inversion [251][252][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 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).
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 1 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.

Tuning the magnetic configuration of 2D magnets by strain
Strain can modify the spin configuration of magnetic 2D materials and make them undergo ferromagnetic (FM) to antiferromagnetic (AFM) phase transitions [269][270][271].As a first example-obtained via Monte Carlo calculationsuniaxial compressive strain was applied to CrGeTe 3 monolayers to switch in between FM and AFM phases.The switching of magnetic phases happens by crossing a region of parameter space with a high degree of magnetic frustration (f ) [269].Uniaxial strain breaks the equivalency Cr−Cr bonds and modifies exchange coupling parameters, resulting in a change in magnetic configuration [269].
A zero-temperature, quantum phase transition from a FM to an AFM takes place at about 2.8% uniaxial compressive strain applied along the zigzag direction.The critical temperature increases for larger compresive strain (figure 31(a)).On the other hand, FM and AFM phases are separated by 2.9% compressive strain at zero temperature when strain is applied along the armchair direction (figure 31(b)) [269].

Magnetoelectric multiferroic monolayers
This section is written to highlight a new frontier on 2D materials that is quickly converging with another area of research dealing with bulk materials having combined magnetic and electric ferroic orders.This section will also serve as preamble for a discussion of engineered multiferroics on magnetic bilayers.
NiI 2 is a layered material with a screw spin helix in an otherwise centrosymmetric lattice [272].As it turns out, the electromagnetic coupling takes place within a generalised Katsura-Nagaosa-Balatsky interaction that is written as follows [273]: where P ij is the intrinsic electric polarisation vector (at point i in the lattice, created by a spin at lattice positions i and i), M is the electromagnetic coupling, S i is the spin at location i, and S j the corresponding spin at location j.Such polar helimagnetic phase persists down to a monolayer, making NiI 2 the first experimentally-verified 2D multiferroic hosting an improper electronic ferroelectricity, where P is driven by an inversion-symmetry-breaking screw spin helix, in an otherwise centrosymmetric lattice.Spin interactions beyond nearest-neighbour Ni sites stabilise the screw spin helix.
The discussion of monolayers ends at this point, to give way to a discussion of novel phenomena taking place in few-layer materials.

Strained multilayers (including moirés)
The emphasis turns into stacks of 2D materials in what follows.We start with a discussion of non-rotated bilayers and few-layers, and end with a discussion of few-layers with relative rotations (moirés).

Energy landscape of slid homo-bilayers without relative rotations
Homo-bilayers are created out of monolayers made from the same 2D material.A hetero-bilayer, on the other hand, is made from two dissimilar 2D materials and it will always produce a moiré.The discussion in the rest of this review is mostly on homo-bilayers and on stacks made from the same 2D material.Hetero-layers will only be discussed when dealing with encapsulated devices here.
To understand the atomistic reconstructions taking place in moirés, it is illustrative to first realise that there is an optimal stacking of two monolayers on a homo-bilayer without relative rotations.The optimal stacking can be calculated using ab initio techniques that yield (sliding) energy landscapes similar to those displayed in figures 24(subplots (i)) and 30(c) (later on, figure 43(a) will feature an energy landscape as well).
Indicated in figure 32(a), graphene's preferred bilayer stacking is labeled AB [275] (or BA) and it is two-fold degenerate.Binary compounds expand the number of possible stacking sequences, and a hBN bilayer can be either set on its lowermost energy configuration AA' (figure 32(b)), or on the AB configuration seen in figure 32(c) (which is related by a 60 • rotation-or a reflection of the upper monolayer in the AA' structure given the three-fold symmetry).
The energy landscapes on figures 33(a), (c), (e), (g) and (h) indicate the optimal relative placement of a homo-bilayer without moiré.The important observation is the existence of two degenerate energy minima per unit cell in figures 33(a), (e) and (h).Those minima correspond to the structures shown in figures 34(a), (b), (c), (d) and (f).The (degenerate) energy landscapes can be fitted to the following function [274]: Recent experimental work showed a transition from a ferroelectric configuration on WSe 2 bilayers (one in which the intrinsic electric dipole moment P is finite) onto a paraelectric one (in which P = 0) at finite temperature [279].The atomistic phenomena leading to a paraelectric behaviour on those homo-bilayers is quite different to that observed in ferroelectrics before, where a two-fold degenerate energy landscape (figure 35(a)) suffices to describe such behaviour.
The first observation to justify a novel description of ferroelectricity on bilayers lacking a center of inversion is the presence of multiple degenerate minima (two per unit cell) that are all accessible by sliding (figure 35(b)).
The second observation is that the paraelectric phase is created by sliding events on the honeycomb lattice, with the outof-plane polarisation P changing direction at every discrete sliding step shown in figures 35(d)-(f) [276].
A third and most remarkable experimental observation is that layered ferroelectrics can sustain an intrinsic electric polarisation that increases in discrete steps with the number of added monolayers [280].
As it turns out, it is extremely challenging to realise bilayers with lattice vectors of each monolayer forming exact multiple of 30 • angles with respect to one another [281].This fact will give rise to a wealth of physical behaviours that will be surveyed in the next section.The word moiré has an arabic origin and it means 'to wet.'The term is commonly used within the textile industry to describe  33(e) for the two degenerate minima within the unit cell).Adapted with permission from [277].Copyright (2017) American Chemical Society.This phenomenology extends to other readily available bilayers, such as TMDC (MX 2 ) bilayers, which have a center of inversion in their 'hexagonal' (H) conformation, but lack a point of inversion in their 'rhombohedral' (R) configuration, and develop an out-of-plane intrinsic electric dipole P as well.Reprinted with permission from [278].Copyright (2018) American Chemical Society.Sliding events are observed at temperatures above 490 K; the magnitude of those displacements is consistent with the vectors drawn as an inset in plot (c) which furnish a honeycomb lattice.The critical temperature T C is assigned to that for which the first sliding event occurs within the full 1 µs simulation time.Reprinted with permission from [276].Copyright (2022) American Chemical Society.
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].

Different ways to create moirés.
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.
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][293][294][295] and multiple other publications.Let 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 Reproduced from [292], with permission from Springer Nature.monolayer is rotated by −θ/2, while the upper monolayer is rotated by θ/2.The respective rotation matrices are: , and the moiré superlattice vectors (figure 37(a)) are obtained as: where 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: or ) , and ) with the length of the moiré superlattice vector given by a for θ < 15 • (when expressed in radians).The area of the moiré is given by Area , 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][300][301][302][303][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].
At very small twist angles (θ ∼ 1 • ), 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 sampledependent 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][314][315].
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].

Real-space topology and strain on moiré bilayers.
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].'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].
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: R, R −1 , G, G −1 , B, and B −1 .
the definition of operations in terms of group generators a and b permits defining a vortex operation: which, starting at a configuration BA, goes clockwise (from left to right; figure 40(d)).An antivortex operation goes anticlockwise on figure 40 [315]: Figure 41 combines analyses of topology and geometry of a WSe 2 /MoSe 2 moiré hetero-bilayer.

Magic angle twisted bilayer graphene.
A powerful method to describe the structure of twisted 2D materials relies on the cut and projection method used to generate quasicrystals [294,[318][319][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 × 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 Φ(r) = ( ψ 1 (r), ψ 2 (r), χ 1 (r), χ 2 (r) ) T , where the index 1 or 2 labels a graphene monolayer and ψ j (r) and χ 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]: Zero-mode operators are defined as: and where ∂ = ∂ x + i∂ y is the antiholonomic differential operator and ∂ = ∂ x − i∂ y the holonomic one.The dimensionless coupling (potential) is given by: where ϕ = 2π/3 and the moiré reciprocal lattice vectors are given by q ) (see figure 37(b)).
The moiré modulation vector is where θ is the twist angle between layers and k D = 4π 3a is the Dirac wave vector norm.The angle and interlayer coupling effects are both captured by the parameter α, defined as: with w 1 = 110 meV the interlayer coupling for AB stacking and v F = 19.81eV 2kD is the Fermi velocity.The spectrum can be found by proposing a Bloch wave function with momentum k in the moiré first Brillouin zone: and 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 2 as a function of α obtained at three k-points: Γ-, K, and k = q 1 + 0.14q 2 + 0.23q 3 , the last one being a k-point of low symmetry.Working with E 2 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 Γ−point and the dips in figure 42 correspond to the magic angles α = α(θ), because the band shrinks to zero width for those values of the twist angle θ (see equations ( 73) and ( 74)).
The contributions from F † (r) and F(r) can be condensed into a new operator 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 ∇ 2 term, (ii) a confinement, kind of quantum dot, potential energy |U(r)| 2 , and (iii) the interaction operator F(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 α(θ m+1 ) − α(θ m ) = 3/2, where m is the order of the twist angle θ [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][323][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][330][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 2 [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].
4.3.5.Antiferroelectric moiré bilayers.At the 2D limit, the electric field effect can penetrate a bilayer, and ferroelectric metals exist in WTe 2 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][301][302].That the intrinsic polarisation P increases in discrete steps with additional layers has been demonstrated as well [280].

Flat bands, charge localisation, and electronic correla-
tions in non-graphene moiré homo-bilayers.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][345][346].Strong electron-electron correlations on flat-band non-graphene bilayers are being vigorously studied now [347].

An argument for charge localisation in graphene and hexagonal boron nitride moiré homo-bilayers
The domain walls separating AB and BA regions develop a limit length [292].At that limit, electronic states belonging to AB and BA domains happen to be orthogonal.This is proven for a graphene bilayer, and it can be proven for a hBN bilayer along similar lines.
To prove the orthogonality of electronic states belonging to either and AB or a BA domains in bilayer graphene, the πelectron Hamiltonian in [275] will be employed.In the absence of external out-of-plane electric field or dopants, it reads as: with f(k) the usual factor introduced in [1], γ 0 = 3.16 eV, γ 1 = 0.381 eV, γ 3 = 0.38 eV, γ 4 = 0.14 eV, and ∆ ′ = 0.022 eV [275,348].
The π-orbitals used to build the Hamiltonian in equation ( 82) are ordered as A 1 , B 1 , A 2 , and B 2 , where the letter A or B indicates the sublattice within a given layer, and the subindex (1 or 2) indicates the layer to which a given orbital belongs.
Starting from an AB bilayer, the Hamiltonian for a BA bilayer requires an inversion along the z-direction (see figures 34(a) and (b)), or a monolayer relabeling leading to: Focusing on the energy dispersion at the Dirac points, the question is whether an electronic state at zero energy at a AB domain at one K-point can be transmitted onto the BA domain at the same K-point.This question is addressed by looking at whether the relevant states are orthogonal: if they were, states at the AB domain cannot be transmitted onto a BA domain.

Correlated physics in triple-layer graphene
More recently, triple layer twisted graphene has been found to be the most strongly interacting correlated material [349,350], and quasi-crystalline electronic phases have been uncovered on those systems [351].Given that they share similar quantum phase diagrams, such discoveries unveiled the possibility of using 2D materials to understand unconventional superconductivity in cuprates and heavy fermion systems [290,312,349,352,353].
The search for superconductivity and other correlated phases in graphene-based materials was extended to nontwisted graphene multilayers.Remarkably, superconductivity was reported in ABC (or rhombohedral) trilayer graphene [354] and Bernal-stacked blayer graphene.
For one of the superconducting phases found in ABCtrilayer graphene, multiple studies have focused on a purely electronic mechanism [355].A proposal by Ghazaryan et al [356] pointed out a general mechanism for superconductivity induced by an annular Fermi surface.Remarkably, an annular Fermi surface has not only been observed in ABC-trilayer graphene, but in other graphene superlattices [154] as well.

Modified optical properties of transition metal dichalcogenide monolayers by strain
Single-photon emitters from 2D materials provide unique advantages over existing quantum light sources and may have potential applications in quantum cryptography and quantum computing.Recently, a van der Waals heterostructure consisting of a WSe 2 monolayer, hBN, and few-layer graphene (FLG) was fabricated to produce electric fieldtunable single photon emission [29] by optical excitation.As seen in figure 46, demonstrations of electrical excitation of photoemitters on strain-tunable van der Waals heterostructures have been provided as well [357].
By applying an electric field, the electron and hole wave functions become spatially separated, resulting in a change in the total energy of the trapped exciton.Electric field-tuned SPE can be observed at positions where strain is induced by a nanopillar.These strain-induced emitters have higher emission energies (blueshifts) than previously reported quantum emitters in WSe 2 monolayers [29,358].

Ferroelectric non-centrosymmetric few-layer elemental stacks
In connection with section 3.6.8,elemental, non-centrosymmetric four-layer graphite has been shown to turn ferroelectric theoretically (figure 47) [359,360] and experimentally [361].Polar and non-polar configurations can be reached by the relative sliding of specific monolayers within the stack.Besides it being yet another example of elemental ferroelectricity, this recent finding of polar behaviour on few-layer graphite is quite relevant from a historic perspective.
Indeed, there was a strong focus on ways to open a gap in the early days of graphene research (this in order to create transistors out of graphene).One way to do so involves placing a graphene bilayer under an external electric field [275,362,363].The ability to create an intrinsic or internal electric field in tetralayer graphene results in the creation of a (small) band gap.

Tuning magnetism on few-layer magnets by tensile uniaxial strain
Strain has been put forward as a way to control the magnetic properties of thin layered materials.To give an example, a strain-induced antiferromagnetic-to-ferromagnetic phase transition was realised on ∼20 nm thick CrSBr flakes; magnetoelastic couplings can give rise to skirmion patterns [364].

Engineered bilayer multiferroics by rotations and shear strain
Creating intrinsic electric dipoles can also be combined with magnetism, providing a novel method for magnetoelectric multiferroics in rotated 2D materials.
By rotating the upper monolayer by 60 • , the centre of inversion of the CrI 3 AB configuration is removed, resulting in an intrinsic electric dipole.A rotated bilayer (labelled s 1 in figure 48) belongs to the crystalline space group Cm.The upper layer of the bilayer s 1 can be further displaced by approximately 2a 1 /3 (labelled s 2,1 ) or 2a 2 /3 (labelled s 2,2 ), resulting degenerate structures (s 2,1 and s 2,2 ) that still belongs to the same crystalline space group, Cm.As a result, it creates an out-of-plane electric dipole moment and an even stronger in-plane intrinsic electric dipole.

Conclusion
It may be fair to say that a natural state of atomistic membranes is one of local strain.Here, experimental and theoretical avenues to strain graphene and other metallic, insulating, ferroelectric, ferroelastic, ferromagnetic and multiferroic 2D materials were reviewed.A broad emphasis was placed on: • valley and sublattice polarisation (P) in graphene • time-dependent strain on graphene • local and/or global strain and superconductivity and other highly correlated and/or topological phases of graphene • piezoelectric P hBN monolayers • strain and the optoelectronic properties of TMDCs • ferroic 2D materials with intrinsic elastic (σ), electric (P) and magnetic (M) polarisation under strain, as well as incipient 2D multiferroics • moiré few-layers exhibiting flat electronic bands and exotic quantum phase diagrams, and other bilayer or few-layer systems exhibiting ferroic orders tunable by rotations and shear strain.
Realisations of a tunable two-dimensional Quantum Spin Hall effect in germanene, of elemental 2D ferroelectric bismuth, and 2D multiferroic NiI 2 were highlighted as well.
The comprehensive discussion of effects of strain facilitated within this Update can be seen as an integral and necessary component for the holistic understanding of 2D materials.The authors hope that the breadth and the number of extremely recent works reviewed help jumpstart research avenues for a broad range of researchers on this burgeoning field.
CONACyT, México for a PhD Scholarship.S B L was funded by the U.S. Department of Energy (Award DE-SC0022120).H N acknowledges support from the MonArk NSF Quantum Foundry supported by the National Science Foundation Q-AMASE-i program under NSF Award No. DMR-1906383 and the Department of Defense under award FA9550-23-1-0500.We thank E Andrade-Amezcua, R Carrillo-Bastos, L Navarro, P Pantaleon, F Guinea, J Roll, J Ferrer, A Garcia-Fuente, A Pacheco-Sanjuan, P Kumar, M Asmar, A Huamán, S Patel, A Singh, and A Ceferino for scientific discussions, help with figures, or for their critical reading of the manuscript.

Figure 1 .
Figure 1.Effect of a uniform (15%) tensile deformations at different angles θ on TB hopping parameters.Three sets of lattice vectors for different values of θ are shown in the main plot, while the tensile deformation is shown as an inset.The radiuses of the circles indicate the strength of hopping parameters from an atom centred at the origin.The colour red (blue) indicates a positive (negative) hopping parameter.The colour map corresponds to the angle θ, which varied from 0 to π. Second-and third-nearest neighbour interactions require a larger relative correction when compared to first-nearest neighbour interactions.Reprinted with permission from [40].Copyright (2018) American Chemical Society.

Figure 2 .
Figure 2. (a) Tip-induced tearing and folding at a step edge.Inset: zoom-in of the folded edge.(b)-(d) Folded graphene on a graphitic surface.Insets are snapshots taken prior to folding.Reprinted with permission from [82].Copyright (2022) by the American Physical Society.(e), (f) Graphene folding with specific chirality.(g) Localised electronic states at folds.From [83].Reprinted with permission from AAAS.

Figure 3 .
Figure 3. (a) Scanning electron microscope image of a graphene fold (dark line within cyan dashed rectangle).(b) Hall bar device created within the cyan rectangle in (a): graphene is encapsulated with hBN, and Hall bar contacts are made.(c) Device image: side contacts to graphene are labelled G1 and G2, while contacts to the graphene fold (nanowire) are N1 and N2.The red dashed line between blue lines highlights the linearly-shaped strain region.Reprinted with permission from[84].Copyright (2018) American Chemical Society.

Figure 4 .
Figure 4. Sublattice symmetry breaking by gaussian deformations in graphene.(a) Gaussian deformation induced by an STM tip (amplitude H = 1 Å, width b = 5 Å).(b) The deformation (black dashed line) follows the STM tip (whose atomistic structure is depicted by red atoms and a blue apex), leading to the apparent STM image (blue curve) which is lifted with respect to the relaxed one (red curve).Yellow bars represent tunnelling currents.(c) Pseudomagnetic field pattern of the gaussian deformation of panel (a).The brightness corresponds to the LDOS calculated from a first nearest neighbour TB Hamiltonian.Insets (white squares) are areas of maximum Bs. (d)-(g) Constant-current STM images of the same graphene area.Hexagons and the two sublattices were overlaid.The LDOS is higher in the sublattice corresponding to the red dots at higher currents.Reprinted with permission from[101].Copyright (2017) American Chemical Society.

Figure 5 .
Figure 5.Valley polarisation and inversion in a graphene ripple.(a) STM image.Dots mark positions where tunnelling spectra was measured.The pseudo-magnetic field measured in the region between the two dashed curves was zero and valley-polarisation was not detected in Landau levels there.(b) Pseudo-magnetic field Bs along the arrow in panel (a).(c) Three tunnelling spectra measured at different positions in panel (a) at B R =0 T. Pseudo-Landau level indexes were marked.(d) Tunnelling spectra curves measured at different positions in panel (a) at B R =11 T. The indexes for the Landau levels are indicated.The subscript + (−) represents Landau levels at the K+ (K−) valley.Reprinted (figure) with permission from [112].Copyright (2020) by the American Physical Society.

Figure 6 .
Figure 6.(a) STM topography showing draped graphene over a step separating two terraces.The draped section forms ripples due to strain.(b) STM topography around a step edge, highlighting locations at which measurements were taken.Reprinted with permission from [122].Copyright (2020) American Chemical Society.

Figure 7 .
Figure 7. (a) A model of Dirac electrons in an alternating magnetic field with period L B is used to understand the dI/dV features of graphene on a draped step (figure 6).(b) LDOS as a function of L B for |B| = 160 T. The black curve depicts Landau levels in graphene with energy peaks at En ∝ √ n.As L B is reduced (blue trace) until it has the same order of magnitude of the magnetic length (red curve), the LDOS turns drastically different from that created by Landau levels, with peaks becoming almost equally spaced.Plots were offset for clarity.(c) and (d) Additional calculations reproduce the spatial variation of the quantised states.(c) Strained and (d) unstrained plus pseudopotential TB models capture features observed in the experimental spectroscopy along the ripple (blue arrow in figure 6(b)).Reprinted with permission from [122].Copyright (2020) American Chemical Society.

Figure 9 .
Figure 9. (a) Scanning electron microscopy image of graphene over gold nanopillars (white spots) and strain-induced ripples (lines in between nanopillars).(b) STM topography of graphene at a region away from the nanopillars.The dI/dV curve shown as an inset displays the characteristic 'V-shape' of unstrained graphene.(c) Fourier transform of the image shown in (b).The superlattice reflects a moiré pattern between graphene and an hBN flake.(d) Measured moiré pattern on strained graphene over the nanopillars (top) and simulated model (bottom).The supercell is shown in black in both panels.(e), (f) Moiré pattern and its Fourier transform on the strained graphene lattice and the hBN substrate at the position marked by the red square in (a).Arrows are the lattice vectors of the distorted graphene.Reprinted with permission from [28].Copyright (2017) American Chemical Society.

Figure 10 .
Figure 10.(a) Strained graphene on a nanopillar array gives rise to second harmonic generation.(b) The shift of Raman 2D peak verifies the strain on graphene.(c) Temperature-dependent emission at 1035/2 nm.Since the wavelength of the pump laser is 1035 nm is, emission at half-wavelength confirms the creation of SHG.Reproduced from [129].CC BY 4.0.

Figure 11 .
Figure 11.(a), (b) Displacement fields in graphene leading to pseudo vector potentials A 1D cos and A 2D cos (the amplitudes were increased for clarity).(c), (d) Corresponding pseudomagnetic fields Bs = ∇ × As. (e), (f) Corresponding (sublattice-dependent) SC order parameter ∆ A/B (A orange, B blue), which peaks at the extreme values of Bs.Adapted from [139].© IOP Publishing Ltd.All rights reserved.

Figure 12 .
Figure 12.Low-energy electronic dispersion in the normal state of unstained (orange) and strained (blue) graphene for the 1D cosine potential A 1D cos (subplots (a) and (b)) and the 2D cosine potential A 2D cos (subplot (c)).The dispersion arising from a 1D strain is shown in both the mixed and reduced zone schemes for better visualisation.Adapted from [139].© IOP Publishing Ltd.All rights reserved.

Figure 13 .
Figure 13.(a) Buckled graphene with pseudomagnetic field Bs indicated in false colour.(b) The local density of states shows the formation of an effective, large-scale honeycomb lattice.Valley-projected band structure of the graphene superlattice across a path in the mini-Brillouin zone (c) without strain, and (d) with buckling-induced strain.Reprinted with permission from [145].© IOP Publishing Ltd.All rights reserved.

Figure 14 .
Figure 14.The Kek-O (subplot (a)) and two possible Kek-Y honeycomb structures (subplots (b) and (c)): as schematically indicated by thicker or thinner lines, bond strengths are different.The √ 3 × √ 3 unit cells were drawn with dashed lines.The two Kek-Y structures are related by a reflection with respect to a vertical line.

Figure 15 .
Figure 15.Competing phases when electron-electron interactions and deformations are included: (a) semimetallic (SEM) graphene, (b) honeycomb antiferromagnetic insulator (AFI), (c) Kek-O or dimerised Kekulé-like insulator (DIM), and (d) distorted hexagonal (HEX) insulator.The labels t A , t B and t C are hopping integrals for the corresponding TB models used in the kinetic part of a Hubbard Hamiltonian that includes the elastic energy due to elastic distortions given in equation (41).Reprinted (figure) with permission from [160], Copyright (2018) by the American Physical Society.

Figure 16 .
Figure 16.(a) Ground state energy E relative to graphene's (SEM) energy E − E SEM versus strain ε.The trends were obtained by DMC and DFT for the insulating Kekulé like dimerised (DIM), antiferromagnetic insulator (AFI), and the distorted hexagonal insulator (HEX) phases.(b) Stress (σ)-strain (ε) relation for strained graphene obtained by fitting DMC energies.Dashed lines mark the transition stress values for the SEM-DIM transition.(c) Enthalpy H of strained graphene relative to the SEM phase H SEM versus tensile stress σ.The blue-shaded region indicates the error bars on the enthalpies for DIM and AFI phases by DMC.Reprinted (figure) with permission from [160], Copyright (2018) by the American Physical Society.

Figure 17 .
Figure 17.Effect of strain on a Kekulé distorted lattice.Lattices in real space, first Brillouin zones, and energy dispersion for: (a) pristine graphene, (b) Kek-Y distorted graphene, (c) strained graphene, and (d) Kek-Y distorted graphene with strain.The Kekulé and strained Kekulé vectors G and G ′ are indicated in the upper side of the hexagonal first Brillouin zones.In (a), the Dirac cones are shown in dark gray.In (b), the gray Dirac cones are folded into a degenerate (blue) Dirac cone to the Γ point.In subplot (c), the original Dirac cones (in gray) are deformed and translated to the (red) points K ± D .Dirac cones (gray) are folded into the Γ point in panel (d).Strain breaks the degeneracy resulting in an overlap of the two shifted Dirac cones, indicated in purple.Reprinted (figure) with permission from [171].Copyright (2019) by the American Physical Society.
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 SAW .

Figure 18 .
Figure 18.Synthetic Hall voltages in graphene by acoustic waves.(a) Top: Raman spectra of graphene (MLG).Inset: diagram of graphene on substrate.Carrier concentration was tuned by a voltage V BG applied to the p-doped substrate.Bottom: two measurement configurations.In the magnetotransport configuration, a constant current I is passed through the Hall bar.In the acoustic transport configuration, a microwave frequency is applied to the inter-digitated top contacts (IDT) and the resulting acoustic current, I SAW , is measured in a short-circuited configuration.(b) A typical magnetotransport measurement for V BG = −20 V for the Hall resistivity Rxy and magnetoresistance Rxx versus magnetic field B. The inset is carrier density versus V BG .(c) In the acoustic measurement, I SAW as a function of frequency, power and V BG for the fundamental resonance of the IDT.Data were obtained at 4.2 K. Reprinted (figure) with permission from [4].Copyright (2022) by the American Physical Society.

Figure 19 .
Figure 19.(a) Topography and (b) electric images of a hexagonal boron nitride monolayer triangular bubble.(c) Profiles along the lines in (a) and (b) and in a corresponding dielectric image.(d) Simulated topography image and corresponding (e) electric-field energy density and (f) polarisation.Reproduced from [35].CC BY 4.0.

Figure 21 .
Figure 21.Determination of strain by second harmonic generation on a MoS 2 monolayer: the arrows indicate the strain field.A SEM image is shown as an inset.Reproduced from [130].CC BY 4.0.

Figure 22 .
Figure 22.Strain dependence of exciton transition energy as well as Raman shifts (out-of-plane A1g mode) for MoS 2 and MoSe 2 monolayers.Solid and dashed lines are linear fits to the data from mechanically-strained monolayers[203,205] and theory[206], while filled squares denote those for PVD-grown monolayers on Si 3 N 4 (blue square) and SiO 2 (orange square) in which strain is due to thermal expansion mismatch induced after the high-temperature growth.Grey lines are linear fits to data from strained PVD monolayers.

Figure 23 .
Figure 23.(a) Band diagram of a WSe 2 monolayer with (spin-allowed) bright (X 0 ) and (spin-forbidden) dark exciton (X D ) states.∆E CB is the spin-split energy difference at the conduction band.(b) Schematics of tip-enhanced photoluminescence with a cantilever gold tip.The bottom and tip metals create a Purcell (optical resonance) effect.(c) Bright exciton and (d) dark exciton intensities.(e) Average cross-section of intensities and strain for blue dotted regions in (c) and (d).(f) Increasing bright (X 0 ) and dark (X D ) emission under tensile strain from data shown in (e).(g) Average cross-section of intensities and strain for red dotted regions in (c) and (d).(h) Increasing bright (X 0 ) and dark (X D ) emission under compressive strain.Adapted with permission from [38].Copyright (2023) American Chemical Society.

Figure 24 .
Figure 24.Highlighting similar features of ferroic 2D materials: (i) energy landscapes featuring energy degeneracies and small energy barriers J < 1000 K per primitive unit cell, (ii) atomistic structures with enhanced symmetry at the top of the barrier, and (iii) thermally driven 2D structural transformations.(a) Ferroelastic silicene.Reprinted (figure) with permission from [194].Copyright (2023) by the American Physical Society.(b) Quantum paraelectric SnO monolayer.Reprinted (figure) with permission from [237].Copyright (2019) by the American Physical Society.(c) Ferroelectric and ferroelastic SnSe monolayer.Reprinted (figure) with permission from [26], Copyright (2022) by the American Physical Society.Reprinted (figure) with permission from [238], Copyright (2018) by the American Physical Society.Reprinted (figure) with permission from[239], Copyright (2020) by the American Physical Society.Note that traversing the landscape from one minima up to the height of the barrier J requires interatomic distances to change, and hence, that the transformations along the landscape can also be realised by strain.

Figure 25 .
Figure 25.(a) Schematics of the relation the lattice parameters of the rectangular unit cell a 1 and a 2 , and ∆α leading to equation (54).Reprinted (figure) with permission from [238].Copyright (2018) by the American Physical Society.(b) Experimental coalescence of the rhombic distortion angle ∆α on a ferroelectric SnTe monolayer grown on epitaxial graphene: both ∆α and the in-plane intrinsic electric dipole Px become zero at 270 K. From [248].Adapted with permission from AAAS.

Figure 26 .
Figure 26.(a) Thermal evolution of lattice parameters of a ferroelectric SnSe monolayer as a function of an effective temperature Umax.(b) Evolution of lattice constants.The evolution is calculated using a normalised Boltzmann distribution over the energy landscape.Reprinted (figure) with permission from [26].Copyright (2022) by the American Physical Society.
(b), SnSe monolayers are softer than graphene, for which C

Figure 27 .
Figure 27.(a) Structural snapshot of a SnS monolayer.The four atoms in the unit cell are indicated.The direction of the intrinsic electric dipole are shown by red arrows in the top view, and the side view features a 0.3 eÅ −3 isosurface.(b) Three direct optical transitions between valence and conduction bands.(c) Monolayer states at the nX conduction local valley minima have mostly px orbital symmetry.(d) Reduction of the magnitude of the intrinsic electric dipole P as a result of photoexcitation to occupy the states in the nX conduction valley indicated within inset.(e) The longer lattice parameter a 1 elongates and the shorter lattice parameter a 2 shrinks as a result of the photoexcitation.Reprinted (figure) with permission from [7].Copyright (2017) by the American Physical Society.

Figure 28 .
Figure 28.Experimental demonstration of photostriction on orthorhombic, bulk GeS.The sample was excited by a 400 nm femtosecond laser with a photon energy above GeS' bandgap and probed by ∼100 fs ultrashort electron pulses at normal incidence.The reversible photoinduced mechanical deformation of the lattice was obtained from the shift of Bragg peak center positions in reciprocal space (Q).The {020} AC ({200} ZZ ) direction corresponds to the x (y) direction in figure27(a).As seen in figure27(e), an anisotropic strain in which the x direction compresses and the y direction elongates is experimentally verified.Reprinted with permission from[8].Copyright (2023) American Chemical Society.

Figure 29 .
Figure 29.(a) Ferroelectric SnSe monolayer clamped onto a strained substrate.(b) Clamping prevents ∆α from becoming zero (panel (i)), which is no longer a good order parameter to tell the structural transformation.Nevertheless, a 2D structural phase transition can still be registered if angles α 1 and α 3 turn identical (panel (ii)), and distances d 2 and d 3 turn equal as well (panel (iii)).(c) Thermal evolution of (i) lattice parameters a 1 and a 2 , (ii) ∆α, (iii) angles α 1 , α 2 , α 3 , and α 4 , (iv) distances d 2 and d 3 , and in-plane intrinsic dipole (Px, Py, 0) under a 2% strain along the x-direction (the 2D material is still able to change a 2 , which is not fixed).As a result of strain, Tc increased up to 500 K.(d) Similar to (c), but for a strain applied along the y-direction: strain raises Tc to 250 K. T C increases regardless of the clamping direction, and it can swap the intrinsic electric dipole from being oriented along the x-direction (plot (c)) onto the y-direction (plot (d)).Solid curves in plots (c) and (d) are the result of a Potts model.Reprinted (figure) with permission from [238].Copyright (2018) by the American Physical Society.

Figure 30 .
Figure 30.As, Sb, and Bi monolayers undergo a structural distortion onto a ferroelectric phase with an in-plane intrinsic electric polarisation P. (a) Generic atomistic structure.(b) Ferroelectric phases (A and A') and a rectangular paraelectric phase (B) with a higher symmetry.(c) Two-fold degenerate elastic energy landscape as a function of the height of the red atoms in (b); blue atoms were left unchanged with respect to their positions in structure B. (d) An electronic reconfiguration of pz orbitals onto outermost atoms drives the creation of a purely-electronic in-plane intrinsic dipole P. Figure taken from [263].John Wiley & Sons.© 2018 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim.

Figure 31 .
Figure 31.(a) Magnetic phase diagram versus percent strain and temperature for a CrGeTe 3 monolayer under uniaxial compressive strain along the zigzag (a ∥ ) direction.(b) Similar to (a), but for compressive strain along the armchair (a ⊥ ) direction.The color coding corresponds to the level of magnetic frustration (f ).Insets depict the AFM phases, with black and white circles representing Cr ions with opposite spin directions.Reprinted (figure) with permission from [269].Copyright (2021) by the American Physical Society.

Figure 33 .
Figure 33.(a) Energy landscape of a graphene bilayer; the AB stacking is the preferred one, and it is degenerate (note the dark blue regions).On the other hand, the AA-stacked configurations (shown in red) are the most energetically unfavorable.(b) Corresponding optimal relative height of the graphene bilayer.(c) Optimal energy and (d) relative height of the hexagonal boron nitride bilayer, originally on the AA' configuration (darkest blue).(e) Energy landscape and optimal relative height of the non-centro-symmetric hexagonal boron nitride bilayer: note the existence of degenerate two minima within the unit cell.Reprinted (figure) with permission from [274].Copyright (2015) by the American Physical Society.(g) Energy landscape of a WSe 2 2H bilayer; it has only one minima in the unit cell.(h) Energy landscape of the non-centro-symmetric rhombic WSe 2 bilayer.Reprinted with permission from [276].Copyright (2022) American Chemical Society.

Figure 34 .
Figure 34.(a) AB bilayer graphene[275].(b) BA bilayer graphene can be obtained by a reflection or by a relative sliding (see degenerate minima in figure33(a).(c) Substitution of carbon atoms in the AB configuration by boron and nitrogen removes the point of inversion, creating an intrinsic electric dipole P that (d) switches direction upon reflection and/or sliding (see figure33(e) for the two degenerate minima within the unit cell).Adapted with permission from[277].Copyright (2017) American Chemical Society.This phenomenology extends to other readily available bilayers, such as TMDC (MX 2 ) bilayers, which have a center of inversion in their 'hexagonal' (H) conformation, but lack a point of inversion in their 'rhombohedral' (R) configuration, and develop an out-of-plane intrinsic electric dipole P as well.Reprinted with permission from[278].Copyright (2018) American Chemical Society.

Figure 35 .
Figure 35.(a) Ferroelectrics are traditionally described by a polynomial energy landscape with two degenerate minima, an energy barrier J, and energy-confining potential energy walls.(b) Rhombohedrally-stacked (R) TMDC bilayers are unusual ferroelectrics having an infinite number of degenerate minima.The paraelectric state is the time average of P over long times; P swaps sign at random times and it averages down to zero.(c)-(f) Demonstrating the temperature-activated relative sliding of the R WSe 2 bilayer: The vector r M−M = (r 1,M−M , r 2,M−M , r 3,M−M ) tracks a pair of W atoms belonging to opposite monolayers on the same unit cell as a function of time.Sliding events are observed at temperatures above 490 K; the magnitude of those displacements is consistent with the vectors drawn as an inset in plot (c) which furnish a honeycomb lattice.The critical temperature T C is assigned to that for which the first sliding event occurs within the full 1 µs simulation time.Reprinted with permission from[276].Copyright (2022) American Chemical Society.

Figure 36 .
Figure 36.Strain from (a) twisting, (b) isotropic strain, or (c) pure shear have signatures that allow identification from one another.Reproduced from[292], with permission from Springer Nature.

Figure 37 .
Figure 37. (a) Wigner-Seitz supercell of a twisted graphene moiré bilayer (atomic positions were not relaxed).Regions with AA, AB and BA stacking configurations, as well as the moiré lattice vectors a 1 and a 2 are indicated.(b) Reciprocal lattice of the twisted bilayer graphene.High-symmetry points are indicated along with moiré reciprocal lattice vectors b 1 , b 2 , b 3 and a set of vectors q 1 , q 2 , q 3 .(a) Reprinted (figure) with permission from [296].Copyright (2021) by the American Physical Society.(a) Reprinted (figure) with permission from [297].Copyright (2022) by the American Physical Society.

Figure 38 .
Figure 38.TBG in plane relaxation for (a) the top and (b) bottom graphene layers for a mismatch angle of θ = 1.08 • .The vectors represent the displacement field.Its magnitude is also highlighted by colour.Perfect AA stacking points are at the corners.The minimally strained regions are AA stacking-like.Vortices of the strain field are seen around AA like regions.One dimensional channels are also seen.Reproduced from [19].CC BY 4.0.

Figure 39 .
Figure 39.(a) Subplots (i)-(iv): definition of possible atomic displacements; the direction of relative monolayer motion is coded by color.Subplot (v) folds the displacements back into a single unit cell, and it defines an order parameter.(b) Transmission electron microscopy and Fourier transforms permit identifying the type of moiré created.For instance, the domains seen in subplots (i)-(v) correspond to a moiré created away from a bilayer near an AB configuration, while those in subplots (iv) and (vi) take place in a moiré near an AA' configuration.(c) Identification of domains and domain walls (see figure 36).(d) Probability of displacements folded into a single unit cell: coloring is consistent with that in subplot (c).Reprinted (figure) with permission from [315].Copyright (2023) by the American Physical Society.

Figure 40 .
Figure 40.(a) The order parameter observed in figure 39(d) can be redrawn (b) as a cut torus that avoids configuration AA, or (c) as three lines joining two points AB and BA.(d) A loop round a configuration AA, which is thus avoided.(e)-(h) Real-space dislocations corresponding to clockwise paths in configuration space, generating vortices ((e) and (f)) and antivortices ((g) and (h)).Domain walls were coloured based on the R, G, or B move in configuration space.The direction of the moves is shown by black arrows.The structures were created from: (e) isotropic, (f) twist, (g) uniaxial, and (h) shear displacements.Reprinted (figure) with permission from [315].Copyright (2023) by the American Physical Society.

Figure 41 .
Figure 41.(a) Dark-field TEM image of a WSe 2 /MoSe 2 hetero-bilayer.The scale bar is 100 nm.The inset contains a diffraction pattern and Burgers vectors.(b) Vortex density.(c) Isotropic strain map.(d) Uniaxial strain map.(e) The shear strain has magnitudes around 1% in highly elongated domains.(f) Twist is the largest contributor to the moiré.Reprinted (figure) with permission from [315].Copyright (2023) by the American Physical Society.

Figure 42 .
Figure 42.Numerically found lowest-band logarithm of the squared energy of the twisted bilayer graphene versus the parameter α (related to θ): The upper red dots corresponds to the Γ point; black dots to the K−point, while purple dots correspond to another generic k point, i.e. chosen no to be a high symmetry point, in this case k = q 1 + 0.14q 2 + 0.23q 3 .At magic α's, indicated by vertical red lines, the bandwidth goes to zero, as indicated by a drop in the band energy values.The K point is always a zero solution for any α.Reprinted (figure) with permission from [297].Copyright (2022) by the American Physical Society.

Figure 43 .
Figure 43.(a) Dimensionless confinement potential |U(r)| 2 for a twisted bilayer graphene model in which the tunnelling at AA regions is set to zero.The degenerate minima of |U(r)| 2 occurs at AA (green dot) and AB (red dots) stacking points.The maximum are located at BA stacking points (yellow dots).The Wigner-Seitz cell of the moiré lattice is indicated (see figure 37(a)).(b) Expectation values of the components of H 2 at the Γ point, as a function of α near the first magic angle (vertical line).Numerical results are indicated with dashed lines and points.−⟨∇ 2 ⟩ is shown in blue, α 2 ⟨|U(r)| 2 ⟩ is green.The off-diagonal contribution ⟨A⟩ = α⟨ Â(r)⟩ is seen in black and orange, corresponding to the first and second terms in equation (79).The solid lines are analytical, perturbative solutions[297].(c) Kinetic plus confinement energy expectation value without the absolute value of the angular momentum contribution at the Γ point.Flat-bands arise at magic angles (vertical lines) whenever the kinetic, confinement and angular momenta contributions add up to zero.The kinetic and confinement potential contributions are positive, while the angular contribution is negative.Reprinted (figure) with permission from[297].Copyright (2022) by the American Physical Society.

Figure 44 .
Figure 44.(a) Small-angle twisted hexagonal boron nitride bilayer.The red circled dot and red circled X represent up and down polarisation, respectively.(b), (c) Vertical piezo-force microscopy phase and amplitude images.From [300].Reprinted with permission from AAAS.

Figure 45 .
Figure 45.(a) Hexagonal boron nitride bilayer with a small moire away from configuration BN/BN-2 in figure 32(c).(b) Hexagonal boron nitride bilayer with a small moire away from configuration BN/BN-1 (lower in energy) in figure 32(b).(c) Flat band on a moiré hexagonal boron nitride bilayer with a rotation angle of 2.64 • away from the BN/BN-2 configuration.Reproduced from [291].CC BY 4.0.

Figure 46 .
Figure 46.(a) Vertical heterostructure (top graphene/WSe 2 monolayer/h-BN/bottom graphene on a polymer substrate).The optical emission site is located below the AFM tip.(b) Band bending of the strained heterostructure under an external bias.E B , E V , E D , and E Defect are the excitonic bright band, valence band, dark excitonic band, and defect state of the WSe 2 monolayer.Solid and dashed red lines in the right panel are spatial variations of the band gap due to local strain.(c) Fabrication process: (i) A PMMA layer was spin-coated on top of the TMDC on the SiO 2 /Si substrate, and peeled off.(ii) The detached PMMA layer was turned upside down.(iii) AFM indentation was performed afterwards.(d) Optical microscope image of the heterostructure.Scale bar: 10 µm.(e) AFM image of the fabricated heterostructure.Seven indents were made to the WSe 2 monolayer.Scale bar: 5 µm.Reproduced with permission from [357].CC BY-NC 4.0.

Figure 48 .
Figure 48.(a) Energy landscape of CrI 3 bilayer with FM (uuuu) and AFM (uudd) magnetic configurations, and optimal vertical separation between monolayers (D).Structure s 1 sits at the point from which lattice vectors a 1 and a 2 join, and it is periodic upon lattice translations.Structure s 2,1 (s 2,2 ) sits 0.6389a 1 (0.6806a 2 ) away from s 1 locations.(b) One-dimensional cuts of the energy landscape along (i) a 1 and (ii) a 2 .(c) Structural stability of s 1 bilayer, and vibrations at saddle point p, displaying imaginary modes.(d) Three-dimensional electric dipole along the straight path joining structures s 1 and s 2,1 .Inset: change in magnetisation at a nearly degenerate crossover AFM and FM configurations at point p (see inset labeled`along a 1 ' in subplot (b)).Reprinted (figure) with permission from [370].Copyright (2023) by the American Physical Society.