Introduction

Forty years ago1 Michael Berry introduced what are now routinely called “Berry phase” and “Berry curvature”2,3,4. Soon after5, Joshua Zak realized that these geometric concepts could be generalized to Bloch-periodic systems, where the parameters (quasi-momenta) are varied in closed loops (bands or Fermi surfaces) by applying electric fields.

Berry phase and Berry curvature have been powerful tools for understanding a plethora of intrinsic (i.e., geometric) contributions to crystals’ properties. These properties of the bulk Bloch bands have been identified to play a pivotal role in a wide range of physical phenomena, including electric polarization6, magnetic oscillations in metals7,8,9, anomalous Hall effects10,11,12, orbital magnetism13, quantum Hall effects of various kinds14,15,16,17,18,19, and quantized charge pumping20. Additionally, they are crucial also in systems where nontrivial real space (as opposed to momentum space) topology emerges, e.g., when spin skyrmion lattices are present21,22. Originally proposed in nuclear physics by Skyrme in 196223, skyrmions have become fundamental in understanding spin structures in condensed matter systems. Initially conceptualized as vortices in the spatial distribution of spin magnetization in crystals24, skyrmions have been experimentally observed in magnetic materials25,26,27 and other physical systems, including focused photonic vector beams28 and suitably-designed space-coiling metastructures for localized spoof plasmons29.

In this Letter, we focus on a different solid-state platform where skyrmion lattices play a crucial role, i.e., twisted transition metal dichalcogenide (TMD) bilayers, such as twisted MoTe2 homobilayers. These materials have recently attracted a great deal of attention because of the experimental discovery of fractional Chern insulating states in zero magnetic field30,31,32,33 and, more recently, superconductivity34,35. Available continuum-model Hamiltonians36,37 describing their single-particle topological moiré bands—which have been recently dubbed skyrmion Chern-band models38,39—harbor orbital (rather than spin) skyrmion lattices.

Intriguingly, because of the presence of this orbital skyrmion lattice, it is possible to approximately map these models into the problem of Landau levels subject to a periodic potential38,40. More precisely, by means of an adiabatic approximation on the layer-pseudospin degree of freedom, the twisted TMD Hamiltonian transforms into one for holes under the effect of a periodic potential and a periodic magnetic field38,40. This effective periodic magnetic field has a non-zero average, and its strength is related to the topological charge of the effective skyrmion lattice36,37,38,40.

In this Letter, we first present a thorough theoretical study of the optical and plasmonic properties of these skyrmion Chern-band models. Although numerical results are presented for the case of twisted MoTe2, the theoretical presentation is carried out in general. We then argue that plasmons in twisted TMD bilayers are a probe of orbital skyrmion textures (see Fig. 1) in the sense that the gap at zero wave number is approximately related to the uniform component of the skyrmion effective magnetic field. Our findings suggest that the effective magnetic field description can be used in order to qualitatively understand the main properties of this class of systems, while an accurate quantum treatment is instead needed for obtaining quantitative results. In addition, we perform a derivation of the mapping onto a single layer-pseudospin sector, which is an alternative to the of refs. 38,40 and based on semiclassical techniques. This gives rise to a formal series expansion, providing us with a systematic way of constructing all higher-order terms in the adiabatic parameter.

Fig. 1: The energy loss function \({\mathcal{L}}({\boldsymbol{q}},\omega )\) of twisted MoTe2 as a function of the in-plane wave vector q and frequency ω.
figure 1

Results in this plot refer to filling factor ν = − 1 and temperature T = 5 K (chemical potential μ ≈ 41 meV). The twist angle is fixed at θ = 3. 1∘. a shows the loss function for energies up to 60 meV. Note that the lowest-energy mode, which is a slow inter-band plasmon, is gapped. The long-wavelength gap measures the q = 0 uniform component of the Skyrme pseudo-magnetic field \({B}_{{\rm{eff}}}^{z}({\boldsymbol{r}})\) defined in Eq. (36). b shows the inter-band mode between the first and second valence bands. The dashed black line represents a thermally activated intra-band plasmon (Eq. (27)). The dotted blue line is an analytical approximation for the gapped inter-band plasmon—see Eq. (30).

Full size image

Results

Plasmons are collective excitations of the electron density in a crystal41, and their behavior can indeed be influenced by various factors, including the underlying band structure and, in principle, its topological properties. To investigate the impact of Berry curvature on plasmons42, it is convenient to set up collisionless hydrodynamic equations, which describe the long-wavelength collective motion of the electron fluid. These can then be coupled to Hamilton equations of motion describing an electron wave packet in a crystal13:

$$\dot{{\boldsymbol{r}}}=\frac{1}{\hslash }{{\boldsymbol{\nabla }}}_{{\boldsymbol{k}}}{\epsilon }_{\lambda }({\boldsymbol{k}})+{{\mathbf{\Omega }}}_{\lambda }({\boldsymbol{k}})\times \dot{{\boldsymbol{k}}},$$
(1)

and

$$\dot{{\boldsymbol{k}}}=-\frac{1}{\hslash }{{\boldsymbol{\nabla }}}_{{\boldsymbol{r}}}V({\boldsymbol{r}},t).$$
(2)

Here, λ is a band index, ϵλ(k) is the band energy, V(r, t) is the instantaneous Hartree potential41, and r, k are the position and quasi-momentum of the Bloch wave packet, respectively. The second term on the right-hand-side of Eq. (1) is the so-called anomalous velocity and affects the motion of the wave packet when the band possesses a non-zero Berry curvature Ωλ(k), defined as:

$${{\mathbf{\Omega }}}_{\lambda }({\boldsymbol{k}})\equiv i\langle {\partial }_{{\boldsymbol{k}}}{u}_{{\boldsymbol{k}},\lambda }| \times | {\partial }_{{\boldsymbol{k}}}{u}_{{\boldsymbol{k}},\lambda }\rangle .$$
(3)

Here, \(| {u}_{{\boldsymbol{k}},\lambda }\left.\right\rangle\) is the periodic part of the Bloch wave function corresponding to a given k and λ.

We now write down the (collisionless) hydrodynamic equations describing plasmons, which are the Euler equation41,

$${\partial }_{t}{\boldsymbol{j}}+{c}_{{\rm{s}}}^{2}{{\boldsymbol{\nabla }}}_{{\boldsymbol{r}}}\delta n-\frac{e{n}_{0}}{{m}^{* }}{{\boldsymbol{\nabla }}}_{{\boldsymbol{r}}}V({\boldsymbol{r}},t)=0,$$
(4)

and the continuity equation,

$${\partial }_{t}\delta n+{{\boldsymbol{\nabla }}}_{{\boldsymbol{r}}}\cdot {{\boldsymbol{j}}}_{{\rm{tot}}}=0.$$
(5)

In Eq. (4), j = p/m* is the current density, cs is the sound velocity41, δn(r, t) ≡ n(r, t) − n0 is the deviation of the carrier density n(r, t) from the uniform density n0 at equilibrium, and m* is the carrier effective mass. Notice that the total, physical current density

$${{\boldsymbol{j}}}_{{\rm{tot}}}\equiv {\boldsymbol{j}}+{{\mathbf{\Omega }}}_{\lambda }({\boldsymbol{k}})\times \dot{{\boldsymbol{k}}}$$
(6)

appears in the continuity equation (5). Finally, V(r, t) in Eq. (4) is related to the density fluctuation δn through the Poisson equation:

$${\nabla }_{{\boldsymbol{r}}}^{2}V({\boldsymbol{r}},t)=-4\pi {e}^{2}\delta n({\boldsymbol{r}},t),$$
(7)

where, for the sake of simplicity, the dielectric environment has been set to the vacuum.

Combining Eqs. (6) and (2), we immediately conclude that the anomalous velocity contribution drops out of the continuity equation, since what matters there is the divergence of the physical current:

$$\begin{array}{lll}{{\boldsymbol{\nabla }}}_{{\boldsymbol{r}}}\cdot \left[{{\mathbf{\Omega }}}_{\lambda }({\boldsymbol{k}})\times \dot{{\boldsymbol{k}}}\right]&=&-\frac{1}{\hslash }{{\boldsymbol{\nabla }}}_{{\boldsymbol{r}}}\cdot \left[{{\mathbf{\Omega }}}_{\lambda }({\boldsymbol{k}})\times {{\boldsymbol{\nabla }}}_{{\boldsymbol{r}}}V({\boldsymbol{r}})\right]\\ &=&-\frac{1}{\hslash }{{\mathbf{\Omega }}}_{\lambda }({\boldsymbol{k}})\cdot \left[{{\boldsymbol{\nabla }}}_{{\boldsymbol{r}}}\times {{\boldsymbol{\nabla }}}_{{\boldsymbol{r}}}V({\boldsymbol{r}})\right]\\ &=&0.\end{array}$$
(8)

This is one of the main results of the pioneering paper by Song and Rudner42, i.e., in the long-wavelength (hydrodynamic) limit, a Bloch momentum-space Berry curvature does not modify the bulk plasmon dispersion. It does instead produce edge plasmons—which were dubbed chiral Berry plasmons by the authors of ref. 42—at the boundary between a topologically non-trivial phase and a trivial one.

In what follows, we study the impact of real-space topology—as induced by an orbital skyrmion lattice—on the optical and plasmonic properties of the hosting material. As mentioned above, we carry out these investigations in the realm of the continuum models that have been introduced36,38,39,40 to describe the topological moiré bands of twisted TMD homobilayers. In particular, we focus on the longitudinal and Hall conductivity and the energy loss function of twisted MoTe2.

Skyrmion Chern-band models

We consider a twisted homobilayer TMD constructed by stacking two TMD monolayers at a relative distance d and rotating one with respect to the other by a twist angle θ. Focusing attention on the K valley and using a 2 × 2 representation in the layer-pseudospin degrees of freedom, the single-particle low-energy valence-band model Hamiltonian of such a system is given by the following expression36:

$$\begin{array}{l}{\hat{{\mathcal{H}}}}_{\uparrow ,v}^{K}({\boldsymbol{r}})=\\ \left(\begin{array}{cc}-\frac{{(\hat{{\boldsymbol{p}}}-\hslash {{\boldsymbol{K}}}_{{\rm{b}}})}^{2}}{2{m}^{* }}+{\Delta }_{{\rm{b}}}({\boldsymbol{r}})&{\tilde{\Delta }}_{{\rm{T}}}({\boldsymbol{r}})\\ {\tilde{\Delta }}_{{\rm{T}}}{({\boldsymbol{r}})}^{\dagger }&-\frac{{(\hat{{\boldsymbol{p}}}-\hslash {{\boldsymbol{K}}}_{{\rm{t}}})}^{2}}{2{m}^{* }}+{\Delta }_{{\rm{t}}}({\boldsymbol{r}})\end{array}\right),\end{array}$$
(9)

where

$${\tilde{\Delta }}_{{\rm{T}}}({\boldsymbol{r}})\equiv w\left[1+{e}^{i{{\boldsymbol{G}}}_{2}(\theta )\cdot {\boldsymbol{r}}}+{e}^{i{{\boldsymbol{G}}}_{3}(\theta )\cdot {\boldsymbol{r}}}\right]$$
(10)

is an inter-layer potential and

$${\Delta }_{\ell }({\boldsymbol{r}})\equiv 2V\sum _{j=1,3,5}\cos \left[{{\boldsymbol{G}}}_{j}(\theta )\cdot {\boldsymbol{r}}\pm \psi \right]$$
(11)

is an intra-layer potential, the + sign applying to ℓ = t (top layer) and the − sign to the ℓ = b (bottom layer) case. These intra-layer potentials, which are periodic over the moiré unit cell, come from different alignments of the metal and chalcogen atoms of the top and bottom TMD monolayers36. The reciprocal lattice vectors Gj(θ) appearing on the right-hand side of Eq. (11) belong to the first shell of the moiré reciprocal lattice and are defined by

$${{\boldsymbol{G}}}_{j}(\theta )=\frac{8\pi }{\sqrt{3}a}\sin \left(\frac{\theta }{2}\right)\left(\cos \left(j\frac{\pi }{3}\right),\sin \left(j\frac{\pi }{3}\right)\right).$$
(12)

Here, a is the lattice parameter of the underlying monolayer TMD of interest. In Eq. (9), the momentum shifts Kℓ correspond to the position of the K valley of the top and bottom layers, after rotation—see Fig. 2a:

$${{\boldsymbol{K}}}_{\ell }=\frac{8\pi }{\sqrt{3}a}\sin \left(\frac{\theta }{2}\right)\left(-\frac{1}{2},\pm \frac{\sqrt{3}}{2}\right).$$
(13)

As in the case of Eq. (11), the + sign applies to ℓ = t while the − sign to the ℓ = b (bottom layer) case.

Fig. 2
figure 2

a The Brillouin zones (BZs) of two monolayer TMDs rotated relative to each other (blue and red big hexagons) are drawn together with the moiré BZ of a twisted homobilayer TMD (small black hexagons), constructed around the K and \({K}^{{\prime} }\) valleys. The high symmetry path Γm-Mm-Km-Γm of the moiré BZ is highlighted by green arrows. b The first five valence bands ϵk,λ of twisted homobilayer MoTe2. The bands are obtained from Eq. (9), for a twist angle θ = 3. 1∘. The dependence on k is displayed along the high symmetry path Γm-Mm-Km-Γm of the moiré BZ. c shows a zoom of the first two valence bands. In b and c, blue vertical arrows connect pairs of nested bands, i.e., mark interband transitions with the highest spectral weight. In b, ℏω02 is between the first and third valence bands while in c, ℏω01 involves the first and second valence bands. The energies associated with these transitions are: ℏω01 ≈ 8 meV and ℏω02 ≈ 33 meV, respectively. d shows the density of states relative to the set of bands displayed in (b). The dashed dot line represents the chemical potential, fixed at μ ≈ 41 meV, corresponding to a filling factor of ν = −1 at T = 5 K.

Full size image

Since TMD monolayers display, as a consequence of strong spin-orbit coupling, spin-valley locking43, the single-particle valence-band Hamiltonian relative to the other valley, \({K}^{{\prime} }\), is obtained from \({\hat{{\mathcal{H}}}}_{\uparrow ,v}^{K}({\boldsymbol{r}})\) in Eq. (9) by applying to it the time-reversal operator36. In practice, this corresponds to changing the sign of the momentum shifts Kℓ in Eq. (13), while simultaneously taking the complex conjugate of the Hamiltonian (9).

The energy bands ϵk,λ can be easily found by diagonalizing the Hamiltonian, provided that, for a fixed value of the twist angle θ, one chooses a physically sound set of parameters \(\left(V,\psi ,w\right)\). In Fig. 2b and c we have reported the first few valence bands of twisted homobilayer MoTe2, as obtained from Eq. (9) at a twist angle θ = 3. 1∘, with the following choice of microscopic parameters44: a = 0.352 nm, m* = 0.6me, me being the bare electron mass in vacuum, and \(\left(V,\psi ,w\right)=(20.8\,{\rm{meV}},107.{7}^{\circ },-23.8\,{\rm{meV}})\). The corresponding density of states has been reported in Fig. 2d.

We now comment on the origin of the name “skyrmion Chern bands” for the bands of the model defined by Eq. (9). The key point is that, upon carrying out a suitable gauge transformation38, the valence-band Hamiltonian (9) reduces to the following physically transparent form:

$$\hat{{\mathcal{H}}}=-\frac{{\hat{{\boldsymbol{p}}}}^{2}}{2{m}^{* }}{\tau }_{0}+{\mathbf{\Delta }}({\boldsymbol{r}})\cdot {\boldsymbol{\tau }}+{\Delta }_{0}({\boldsymbol{r}}){\tau }_{0},$$
(14)

where we have omitted the indices ↑ and K and introduced the following quantities:

$${\mathbf{\Delta }}({\boldsymbol{r}})\equiv \left({\rm{Re}}\,{\Delta }_{{\rm{T}}}({\boldsymbol{r}}),{\rm{Im}}\,{\Delta }_{{\rm{T}}}({\boldsymbol{r}}),\frac{{\Delta }_{{\rm{b}}}({\boldsymbol{r}})-{\Delta }_{{\rm{t}}}({\boldsymbol{r}})}{2}\right),$$
(15)
$${\Delta }_{0}({\boldsymbol{r}})\equiv \frac{{\Delta }_{{\rm{b}}}({\boldsymbol{r}})+{\Delta }_{{\rm{t}}}({\boldsymbol{r}})}{2},$$
(16)

and

$${\Delta }_{{\rm{T}}}({\boldsymbol{r}})\equiv w\left({e}^{i{{\boldsymbol{q}}}_{1}\cdot {\boldsymbol{r}}}+{e}^{i{{\boldsymbol{q}}}_{2}\cdot {\boldsymbol{r}}}+{e}^{i{{\boldsymbol{q}}}_{3}\cdot {\boldsymbol{r}}}\right),$$
(17)

with

$${{\boldsymbol{q}}}_{1}=-\frac{8\pi }{\sqrt{3}a}\sin \left(\frac{\theta }{2}\right)\left(0,\frac{1}{\sqrt{3}}\right),$$
(18)

q2 = G2 + q1, and q3 = G3 + q1. In Eq. (14), the matrices τ and τ0 are Pauli matrices acting on the layer-pseudospin (rather than spin) degrees of freedom.

The Hamiltonian in Eq. (14) is the analog of a spin-skyrmion Hamiltonian, with the skyrmion texture defined by Δ(r)/∣Δ(r)∣36,38. The only difference is that real spin is here replaced by the layer-pseudospin degree of freedom.

Local optical conductivity

We now proceed to calculate the local optical conductivity tensor σαβ(ω) of a twisted homobilayer TMD. The quantity σαβ(ω) is the linear response function relating the electrical current flowing in the Cartesian direction α in response to a uniform electric field applied in the direction β41. By using the Kubo formula41, it can be shown that in a 2D crystal, i.e., in a Bloch translationally-invariant system in which the single-particle eigenstates are of the Bloch type, σαβ(ω) is given by45

$${\sigma }_{\alpha \beta }(\omega )=\sum _{\xi =K,{K}^{{\prime} }}{\sigma }_{\alpha \beta ,\xi }(\omega ),$$
(19)

where

$${\sigma }_{\alpha \beta ,\xi }(\omega )={\sigma }_{\alpha \beta ,\xi }^{{\rm{intra}}}(\omega )+{\sigma }_{\alpha \beta ,\xi }^{{\rm{inter}}}(\omega )$$
(20)

is the sum of intra-band and inter-band contributions. Note that, in Eq. (19), we have explicitly divided the physical optical conductivity σαβ(ω) into its valley-resolved contributions σαβ,ξ(ω), with \(\xi =K,{K}^{{\prime} }\). We stress that, due to the underlying giant spin-orbit coupling at the monolayer TMD level, also in the case of the twisted homobilayer TMD, the spin degree of freedom is locked to the valley one. Therefore, the quantity σαβ,ξ(ω) physically represents a single-valley and single-spin optical conductivity tensor.

The final expressions for the single-spin and single-valley intra- and inter-band optical conductivities are:

$${\sigma }_{\alpha \beta ,\xi }^{{\rm{intra}}}(\omega )=\frac{\hslash }{\pi }\frac{i{{\mathcal{D}}}_{\alpha \beta ,\xi }}{\hslash \omega +i\eta },$$
(21)

where

$$\begin{array}{lll}{{\mathcal{D}}}_{\alpha \beta ,\xi }&=&\frac{\pi {e}^{2}}{{\hslash }^{2}}\sum\limits_{\lambda }\int\frac{{d}^{2}{\boldsymbol{k}}}{{(2\pi )}^{2}}[-{f}_{{\boldsymbol{k}},\lambda }^{{\prime} }]\times \\ &&{\times }_{\xi }\left\langle {\boldsymbol{k}},\lambda \right\vert \frac{\hslash }{{m}^{* }}{\hat{p}}_{\alpha }{\left\vert {\boldsymbol{k}},\lambda \right\rangle }_{\xi }{\langle {\boldsymbol{k}},\lambda | \frac{\hslash }{{m}^{* }}{\hat{p}}_{\beta }| {\boldsymbol{k}},\lambda \rangle }_{\xi },\end{array}$$
(22)

is the Drude weight, and

$$\begin{array}{rcl}{\sigma }_{\alpha \beta ,\xi }^{{\rm{inter}}}(\omega )&=&\frac{i{e}^{2}}{\hslash }\sum\limits_{\lambda \ne {\lambda }^{{\prime} }}\int\frac{{d}^{2}{\boldsymbol{k}}}{{(2\pi )}^{2}}\left[-\frac{{f}_{{\boldsymbol{k}},\lambda }-{f}_{{\boldsymbol{k}},{\lambda }^{{\prime} }}}{{\epsilon }_{{\boldsymbol{k}},\lambda }-{\epsilon }_{{\boldsymbol{k}},{\lambda }^{{\prime} }}}\right]\times \\ &&\times \frac{{\atop{\xi}} \left\langle {\boldsymbol{k}},\lambda \right\vert \frac{\hslash }{{m}^{* }}{\hat{p}}_{\alpha }{\left\vert {\boldsymbol{k}},{\lambda }^{{\prime} }\right\rangle }{\atop{\xi} }{\langle {\boldsymbol{k}},{\lambda }^{{\prime} }| \frac{\hslash }{{m}^{* }}{\hat{p}}_{\beta }| {\boldsymbol{k}},\lambda \rangle }_{\xi }}{{\epsilon }_{{\boldsymbol{k}},\lambda }-{\epsilon }_{{\boldsymbol{k}},{\lambda }^{{\prime} }}+\hslash \omega +i\eta }.\end{array}$$
(23)

Here, ϵk,λ is the eigenvalue associated with the Bloch eigenstate \({\left\vert {\boldsymbol{k}},\lambda \right\rangle }_{\xi }\), fk,λ is the Fermi-Dirac distribution, \({f}_{{\boldsymbol{k}},\lambda }^{{\prime} }\) is its derivative with respect to the energy ϵk,λ, \({\hat{p}}_{\alpha }\) is the α-th Cartesian component of the momentum operator, and η = 0+ is a positive infinitesimal. We point out that the integrals in Eqs. (22) and (23) are performed over the first moiré BZ. Energy bands ϵk,λ and Bloch states \({\left\vert {\boldsymbol{k}},\lambda \right\rangle }_{\xi }\) are found numerically, by diagonalizing the Hamiltonian in Eq. (9).

Figure 3 shows the single-valley and single-spin optical conductivity of MoTe2 evaluated at a twist angle θ = 3. 1∘ and for an integer filling factor ν = −1—see Fig. 2d. The longitudinal part, which is displayed in Fig. 3a, shows two peaks, the first one, at ℏω = ℏω01 ≈ 8 meV, which is associated to an interband transition between the first two valence bands and the second one, at energy ℏω = ℏω02 ≈ 33 meV, which is associated to an interband transition between the first and third valence bands.

Fig. 3: Single-valley and single-spin local optical conductivity σαβ,ξ=K(ω) of MoTe2 as a function of ω in units of G0 = 2e2/h.
figure 3

Results in this plot refer to filling factor ν = −1 and temperature T = 5 K. The twist angle is fixed at θ = 3. 1∘. a shows the longitudinal component σxx,K(ω); b shows the Hall component σxy,K(ω). In both panels, the solid (dashed) line refers to the real (imaginary) part. The colored vertical dashed lines highlight the transition energies ℏω01 ≈ 8 meV (green) and ℏω02 ≈ 33 meV (blue).

Full size image

Figure 3b shows instead the transverse Hall-like component of the optical conductivity tensor, which is nonzero in the single-valley and single-spin sector. In the following, we provide a physical interpretation of the appearance of such a finite Hall-like conductivity, showing that its origin stems from the orbital magnetic moment of the Bloch states (see Section I of ref. 46). Another possibility to explain the finiteness of σxy,ξ(ω) at the level of the single-spin and single-valley model is to invoke a mapping into a Landau-level problem38. In such a Landau-level mapping, σxy,ξ(ω) ≠ 0 is a direct consequence of the effective magnetic field due to the presence of the orbital skyrmion lattice. Finally, the fact that σxy,ξ ≠ 0 ultimately stems from the broken \({C}_{2}{\mathcal{T}}\) symmetry (inversion times time-reversal)3. This is due to the z component of Δ3,47 in Eq. (14), defined in Eq. (15). This contribution allows for a nonzero orbital magnetic moment in each valley47. This is in stark contrast to the case of twisted bilayer graphene, which is \({C}_{2}{\mathcal{T}}\) symmetric (and this is why Dirac points around the moiré K and \({K}^{{\prime} }\) points are protected in this material) and σxy,ξ ≡ 0.

Despite the finiteness of σxy,ξ in our study of twisted homobilayer TMDs, in the case in which time-reversal symmetry is not explicitly broken, the transverse Hall conductivity is odd, i.e., \({\sigma }_{xy,K}=-{\sigma }_{xy,{K}^{{\prime} }}\). The total physical conductivity σαβ defined in Eq. (19) is therefore purely longitudinal. Indeed, the longitudinal single-valley and single-spin optical conductivity is even under the exchange of the valley index, i.e., \({\sigma }_{\alpha \alpha ,K}={\sigma }_{\alpha \alpha ,{K}^{{\prime} }}\).

Although the single-valley and single-spin optical conductivity calculations described in this Section are performed at the single-particle (i.e., non-interacting) level, we expect that they qualitatively capture the physics of valley and spin-polarized states when electron-electron interactions are taken into account44. Indeed, strong interactions tend to lift up the valley and spin degeneracies, leading to a flavor-polarized ground state44,48,49. On the other hand, since spin-valley locking in the twisted homobilayer TMDs is inherited from the underlying monolayer TMD, circularly polarized light could be used to measure the Hall response at finite ω43.

Collective modes: loss function

We now turn to discuss collective modes in twisted homobilayer TMDs. In particular, we focus on the electron energy loss function \({\mathcal{L}}({\boldsymbol{q}},\omega )\), which represents the probability of exciting the electron system by applying a scalar perturbation with wave vector q and energy ℏω that couples to the total electron density. This quantity therefore also describes self-sustained charge oscillations (plasmons), which appear as sharp peaks, and contains information about the spectral density of incoherent electron-hole pairs, which yield a broadly distributed background in the q-ω plane. The energy loss function can be measured in principle via electron energy loss spectroscopy50 and scattering-type near-field optical spectroscopy51,52,53,54,55,56. In crystals, it is defined as:

$${\mathcal{L}}({\boldsymbol{q}},\omega )=-{\rm{Im}}\left\{{\left[\varepsilon {({\boldsymbol{q}},\omega )}^{-1}\right]}_{{\boldsymbol{G}} = {\boldsymbol{0}},{{\boldsymbol{G}}}^{{\prime} } = {\boldsymbol{0}}}\right\},$$
(24)

where \(\varepsilon ({\boldsymbol{q}},\omega )={[\varepsilon ({\boldsymbol{q}},\omega )]}_{{\boldsymbol{G}},{{\boldsymbol{G}}}^{{\prime} }}\) is the dynamical dielectric function, which is a matrix41 with respect to the reciprocal lattice vectors G and \({{\boldsymbol{G}}}^{{\prime} }\).

In the random phase approximation (RPA)41 and in the local limit (i.e., q → 0 and \({\boldsymbol{G}},{{\boldsymbol{G}}}^{{\prime} }={\boldsymbol{0}}\)), this matrix can be expressed as45:

$$\mathop{\lim }\limits_{{\boldsymbol{q}}\to {\boldsymbol{0}}}{\left[\varepsilon ({\boldsymbol{q}},\omega )\right]}_{{\boldsymbol{0}},{\boldsymbol{0}}}\approx 1+i\frac{L({\boldsymbol{q}},\omega )}{\omega }{q}_{\alpha }{q}_{\beta }{\sigma }_{\alpha \beta }(\omega ),$$
(25)

where L(q, ω) is the Coulomb propagator relating the charge density fluctuations δnq(ω) to the self-induced electrical potential, i.e., W(q, ω) = e2L(q, ω)δnq(ω). Notice that, in Eq. (25), the total physical optical conductivity (19) appears.

For a 2D system embedded in a homogeneous and isotropic dielectric environment with a dielectric permittivity \(\bar{\varepsilon }\), the non-retarded Coulomb propagator is (In writing Eq. (26) we have neglected finite-thickness effects. In a semiconducting system of thickness d, the Coulomb propagator \(2\pi /(\bar{\varepsilon }q)\) in Eq. (26) is replaced by the famous Keldysh interaction L(q) = 2πK(q)/q (Keldysh, L. V. Coulomb interaction in thin semiconductor and semimetal films. http://jetpletters.ru/ps/1458/article_22207.pdfJETP Lett. 29, 658 (1979) and Rösner, M., Şaşıoğlu, E., Friedrich, C., Blügel, S., & Wehling, T. O. Wannier function approach to realistic Coulomb interactions in layered materials and heterostructures. https://doi.org/10.1103/PhysRevB.92.085102Phys. Rev. B 92, 085102 (2015)), where \(K(q)\equiv \frac{2[{\epsilon }_{{\rm{i}}}\cosh (qd)+{\epsilon }_{{\rm{b}}}\sinh (qd)]}{{\epsilon }_{{\rm{i}}}({\epsilon }_{{\rm{a}}}+{\epsilon }_{{\rm{b}}})\cosh (qd)+({\epsilon }_{{\rm{i}}}^{2}+{\epsilon }_{{\rm{a}}}{\epsilon }_{{\rm{b}}})\sinh (qd)}.\)

Here, ϵi is the intrinsic dielectric constant of the semiconducting system and ϵa, ϵb are the top and bottom dielectric constants, respectively. This interaction has been used e.g., to calculate the plasmon dispersion in monolayer TMDs (see, for example, Torbatian, Z. & Asgari, R. Plasmonic physics of 2D crystalline materials. https://doi.org/10.3390/app8020238Appl. Sci. 8, 238 (2018)). However, we highlight that such finite-thickness effects do not alter the long-wavelength dispersion relation because, in this limit, the plasmon wavelength is much larger than the thickness d of the TMD. Indeed, taking the limit qd ≪ 1 of K(q), one recovers the Coulomb propagator \(2\pi /(\bar{\varepsilon }q)\) in Eq. (26), with \(\bar{\varepsilon }=({\epsilon }_{{\rm{a}}}+{\epsilon }_{{\rm{b}}})/2\).):

$$L({\boldsymbol{q}},\omega )\equiv \frac{2\pi }{q\bar{\varepsilon }}.$$
(26)

Note that we have neglected the frequency dependence of \(\bar{\varepsilon }\) on purpose, with the aim of highlighting intrinsic features of the plasmonic spectrum of twisted homobilayer TMDs. Also, in our numerical calculations, we have considered an isolated system in vacuum (i.e., a suspended twisted homobilayer TMD), for which \(\bar{\varepsilon }=1\). Extrinsic effects, which may modify the plasmonic spectrum, such as hyperbolic phonon polaritons57,58,59 that appear, e.g., when twisted homobilayer TMDs are encapsulated in hexagonal Boron Nitride (hBN), can be easily taken into account60 but are not of our interest here.

A summary of our main results for the loss function \({\mathcal{L}}({\boldsymbol{q}},\omega )\) of twisted homobilayer MoTe2 is reported in Fig. 1. Results in this figure have been obtained for the exact same choice of parameters as in Fig. 3. In the explored range of values of ℏω, we clearly see two sharp and slow collective modes, the lower energy one falling in the Terahertz spectral range. (Therefore, this low-energy mode will not be modified by the phonon polariton modes of the hBN slabs encapsulating MoTe2, which are at much higher energies, in the mid-infrared range.) The sharpness and slowness of these two collective modes is linked (as in the case of plasmons in twisted bilayer graphene45,61,62) to the flatness of the moiré bands. Their existence can be inferred from the above-mentioned single-particle optical inter-band transitions with the highest spectral weight. Indeed, considering that we have chosen a hole density such that the filling factor is ν = − 1, the lowest-energy single-particle electron-hole excitation is the vertical inter-band transition between the first two flat valence bands, depicted as a blue arrow in Fig. 2c, with associated transition energy ℏω = ℏω01 ≈ 8 meV. Upon including electron-electron interactions, as in our RPA calculations, these single-particle excitations merge into an inter-band collective plasmon mode, whose dispersion relation terminates at a finite value of ℏω given by ℏω = ℏω01 in the long-wavelength q = 0 limit. Similarly, the higher-energy plasmon branch occurring at higher energies, which is clearly visible in Fig. 1, stems from a coalescence of single-particle inter-band transitions involving the first and third valence bands, depicted as a blue arrow in Fig. 2b, with an associated transition energy ℏω = ℏω02 ≈ 33 meV. Notice that, due to inter-band Landau damping (see Section II of ref. 46), this higher-energy plasmon branch appears blurred at long wavelengths. Finally, the gapped behavior of these modes resembles that of magnetoplasmons63, which are collective modes of 2D electron liquids subject to a strong perpendicular magnetic field. We will come back to this point below.

We now comment on the role of finite-temperature effects. Since our calculations are performed at finite T, we need, at least in principle, to check whether the results shown in Fig. 1 can be explained in terms of a thermally activated intra-band plasmon. The latter would have a long-wavelength dispersion given by45

$${\omega }_{{\rm{th}}}(q)=\sqrt{\frac{2{{\mathcal{D}}}_{{\rm{L}}}}{\bar{\varepsilon }}q},$$
(27)

where \({{\mathcal{D}}}_{{\rm{L}}}\equiv {{\mathcal{D}}}_{xx,K}+{{\mathcal{D}}}_{xx,{K}^{{\prime} }}={{\mathcal{D}}}_{yy,K}+{{\mathcal{D}}}_{yy,{K}^{{\prime} }}\) is the longitudinal part of the Drude weight tensor defined in Eq. (22). Such intra-band mode is depicted in Fig. 1(b) with a black dashed line. We conclude that our findings in Fig. 1 cannot be explained in terms of trivial, thermally-activated, intra-band plasmons.

Interpretation of the numerical results

In what follows, we provide a physical interpretation of our main numerical results in Figs. 1 and 3b. Starting from the Kubo formula for the optical conductivity, we first discuss the orbital character of the Hall-like component of the conductivity tensor and then provide a simple analytical formula for the gapped inter-band plasmon dispersion.

Since the chemical potential lies in the gap between the first two flat valence bands and we work at relatively low temperatures (kBT ≪ ℏω01), we neglect the intra-band contribution to the conductivity (21) and concentrate only on the inter-band term in Eq. (23). In the analysis below, we further restrict our study to the contribution to \({\sigma }_{\alpha \beta ,\xi }^{{\rm{inter}}}(\omega )\) stemming from the lowest inter-band transition at ℏω = ℏω01.

Manipulating Eq. (23) in the case of the Hall-like term \({\sigma }_{xy,\xi }^{{\rm{inter}}}(\omega )\) we find (see Section I of46):

$$\begin{array}{l}{\sigma }_{xy,\xi }^{{\rm{inter}}}(\omega )=\\ -ec\sum\limits_{\lambda \ne {\lambda }^{{\prime} }}\int\frac{{d}^{2}{\boldsymbol{k}}}{{(2\pi )}^{2}}\frac{{f}_{{\boldsymbol{k}},\lambda }-{f}_{{\boldsymbol{k}},{\lambda }^{{\prime} }}}{{\epsilon }_{{\boldsymbol{k}},\lambda }-{\epsilon }_{{\boldsymbol{k}},{\lambda }^{{\prime} }}+\hslash \omega +i\eta }{{\boldsymbol{m}}}_{{\boldsymbol{k}},\lambda ,{\lambda }^{{\prime} }}^{\xi }\cdot {\boldsymbol{z}},\end{array}$$
(28)

where

$$\begin{array}{l}{{\boldsymbol{m}}}_{{\boldsymbol{k}},\lambda ,{\lambda }^{{\prime} }}^{\xi }\equiv \frac{ie}{2\hslash c}\times \\ \frac{{\atop{\xi}} \left\langle {u}_{{\boldsymbol{k}},\lambda }\right\vert {\partial }_{{\boldsymbol{k}}}{\hat{{\mathcal{H}}}}^{\xi }({\boldsymbol{k}}){\left\vert {u}_{{\boldsymbol{k}},{\lambda }^{{\prime} }}\right\rangle }_{\xi }{\times}{\atop{\xi} }\left\langle {u}_{{\boldsymbol{k}},{\lambda }^{{\prime} }}\right\vert {\partial }_{{\boldsymbol{k}}}{\hat{{\mathcal{H}}}}^{\xi }({\boldsymbol{k}}){\left\vert {u}_{{\boldsymbol{k}},\lambda }\right\rangle }_{\xi }}{{\epsilon }_{{\boldsymbol{k}},\lambda }-{\epsilon }_{{\boldsymbol{k}},{\lambda }^{{\prime} }}}\end{array}$$
(29)

is the inter-band orbital magnetic moment of the Bloch states4. Here, \({\left\vert {u}_{{\boldsymbol{k}},\lambda }\right\rangle }_{\xi }\) denotes the periodic part of the Bloch state46. Eq. (28) shows that a non-zero orbital magnetic moment of the Bloch states can result in a Hall-like conductivity \({\sigma }_{xy,\xi }^{{\rm{inter}}}(\omega )\). In our case, since we are dealing with a single-valley and single-spin model where the \({C}_{2}{\mathcal{T}}\) symmetry is broken, we obtain a non-trivial orbital magnetic moment texture in momentum space4,47. In Section I of ref. 46 we further manipulate Eq. (28) in order to capture analytically the behavior of \({\sigma }_{xy,\xi }^{{\rm{inter}}}(\omega )\) near ℏω = ℏω01, and show numerical results for \({{\boldsymbol{m}}}_{{\boldsymbol{k}},\lambda ,{\lambda }^{{\prime} }}^{\xi }\) in the first moiré BZ.

As far as the plasmonic modes are concerned, we stress that Eq. (25) is controlled only by the longitudinal part σxx(ω) of the physical optical conductivity since σxy(ω) = 0 due a cancellation of the two valley-resolved contributions σxy,ξ(ω). Imposing the condition \(\mathop{\lim }\limits_{{\boldsymbol{q}}\to {\boldsymbol{0}}}{\left[\varepsilon ({\boldsymbol{q}},\omega )\right]}_{{\boldsymbol{0}},{\boldsymbol{0}}}=0\) for the existence of a plasmon mode, we find the following analytical expression for the long-wavelength dispersion relation (see Section III of Ref. 46):

$${\omega }_{{\rm{pl}}}(q)=\sqrt{{\omega }_{01}^{2}+\frac{4\pi {e}^{2}}{\bar{\varepsilon }\hslash {\omega }_{01}}({\alpha }_{K}+{\alpha }_{{K}^{{\prime} }})q},$$
(30)

where

$${\alpha }_{\xi }\equiv \sum _{\lambda \ne {\lambda }^{{\prime} }\in \{0,1\}}{(-1)}^{{\lambda }^{{\prime} }}\int\frac{{d}^{2}{\boldsymbol{k}}}{{(2\pi )}^{2}}{f}_{{\boldsymbol{k}},\lambda }{| }_{\xi }\left\langle {\boldsymbol{k}},\lambda \right\vert \frac{{\hat{p}}_{x}}{{m}^{* }}{\left\vert {\boldsymbol{k}},{\lambda }^{{\prime} }\right\rangle }_{\xi }{| }^{2}.$$
(31)

Eq. (30) has been obtained by retaining only two valence bands and assuming ϵk,0 − ϵk,1 ≈ ℏω01. Indeed, \(\lambda ,{\lambda }^{{\prime} }=0,1\) in Eq. (31) refer to the first and second valence bands, respectively. As we can see in Fig. 1b, the inter-band plasmon dispersion (30)—represented by a dotted blue line—properly captures the long-wavelength behavior of the numerically calculated dispersion and, in particular, the value of the plasmon gap at q = 0, i.e., ωpl(0) = ω01.

Discussion

Figure 2b and c show that the first and second valence bands are rather flat and separated by a gap, akin to two adjacent Landau levels. This energy structure suggests the use of techniques that have been developed in the context of the theory of the quantum Hall effect41. In these cases, indeed, it is often convenient to project the full Hamiltonian onto a restricted Hilbert space comprising one or two Landau levels64,65.

For twisted MoTe2, the mapping onto a Landau-level model was first proposed in ref. 38. The first step in this approach is to project the Hamiltonian (14) onto a single layer-pseudospin sector21,38,40,66,67. This projection is justified when relative energy splitting between the two pseudospin sectors is large, i.e., when the condition \({p}_{0}^{2}/{m}^{* }\ll\, \Delta\), where p0 is the typical electron momentum and Δ is the typical value of ∣Δ(r)∣. Alternatively, one can write this condition as ℏ2/(m*AM) ≪ Δ, where AM is the area of the moiré unit cell.

In Section IV of ref. 46, we discuss a formalism based on semiclassical techniques68,69,70, with which one can perform such projection to any desired order in the dimensionless parameter ζ ≡ ℏ2/(m*AMΔ). This approach is rooted into the relation between quantum mechanical operators on Hilbert spaces and classical observables on phase space71,72,73. Its main advantage over previously considered methods21,38,40,66,67 is that it provides a systematic way to construct the projection as an expansion in powers of ζ, including all higher-order terms.

In Section IV.C of ref. 46 we show that, to lowest order in the parameter ζ, our effective scalar Hamiltonian for the upper layer-pseudospin sector coincides with the previously obtained result21,38,40,66,67, i.e.,

$$\hat{H}=-\frac{1}{2{m}^{* }}{\left[\hat{{\boldsymbol{p}}}+\frac{e}{c}{\boldsymbol{A}}({\boldsymbol{r}})\right]}^{2}-D({\boldsymbol{r}})+{\Delta }_{0}({\boldsymbol{r}})+| {\mathbf{\Delta }}({\boldsymbol{r}})| ,$$
(32)

where \(\hat{{\boldsymbol{p}}}=-i\hslash {\partial }_{{\boldsymbol{r}}}\) is the usual momentum operator. The vector potential A in Eq. (32) is given by

$${\boldsymbol{A}}({\boldsymbol{r}})=\frac{\hslash c}{2e}\frac{1}{1-{n}_{z}}\left(-{n}_{x}\frac{\partial {n}_{y}}{\partial {\boldsymbol{r}}}+{n}_{y}\frac{\partial {n}_{x}}{\partial {\boldsymbol{r}}}\right)$$
(33)

while

$$D({\boldsymbol{r}})\equiv \frac{{\hslash }^{2}}{8{m}^{*}}\left[{({\partial }_{x}{\boldsymbol{n}})}^{2}+{({\partial }_{y}{\boldsymbol{n}})}^{2}\right].$$
(34)

In Eqs. (33) and (34), n(r) ≡ Δ(r)/|Δ(r)| is the so-called Skyrme texture. We remind the reader that the quantities Δ(r) and Δ0(r) have been introduced above in Eq. (15) and (16), respectively.

Since the Skyrme texture has topological charge36

$$N=\frac{1}{4\pi }{\int}_{{A}_{{\rm{M}}}}{\boldsymbol{n}}({\boldsymbol{r}})\cdot \left[{\partial }_{x}{\boldsymbol{n}}({\boldsymbol{r}})\times {\partial }_{y}{\boldsymbol{n}}({\boldsymbol{r}})\right]{d}^{2}{\boldsymbol{r}}=-1,$$
(35)

the induced real-space magnetic field21,38

$${B}_{z}^{{\rm{eff}}}({\boldsymbol{r}})=\frac{\hslash c}{2e}{\boldsymbol{n}}({\boldsymbol{r}})\cdot \left[{\partial }_{x}{\boldsymbol{n}}({\boldsymbol{r}})\times {\partial }_{y}{\boldsymbol{n}}({\boldsymbol{r}})\right]$$
(36)

carries a total magnetic flux Φ0 = − hc/e through the unit cell of the moiré lattice of area AM.

As discussed in ref. 38, the magnetic field (36) has a non-zero average given by \({B}_{0z}^{{\rm{eff}}}\equiv {\Phi }_{0}/{A}_{{\rm{M}}}\). We can therefore split the vector potential A into a linear part A0, which gives rise to this constant magnetic field, and a remainder with a zero-average over the moiré unit cell. One subsequently introduces \(\hat{{\boldsymbol{\Pi }}}=\hat{{\boldsymbol{p}}}+e{{\boldsymbol{A}}}_{0}/c\) and the raising and lowering operators38,41

$${\hat{a}}^{\dagger }=\frac{\ell }{\sqrt{2}\hslash }\left({\hat{\Pi }}_{y}-i{\hat{\Pi }}_{x}\right),\quad \hat{a}=\frac{\ell }{\sqrt{2}\hslash }\left({\hat{\Pi }}_{y}+i{\hat{\Pi }}_{x}\right),$$
(37)

where ℓ is the effective magnetic length, defined by the condition \(2\pi {\ell }^{2}{B}_{0z}^{{\rm{eff}}}={\Phi }_{0}\). Using the definition of \({B}_{0z}^{{\rm{eff}}}\), we can also write this relation as 2πℓ2 = AM.

In terms of the ladder operators (37), we can write the Hamiltonian (32) as38

$$\hat{H}=-\hslash {\omega }_{{\rm{c}}}\left({\hat{a}}^{\dagger }\hat{a}+\frac{1}{2}\right)+{\hat{H}}_{{\rm{cor}}},$$
(38)

where we have introduced the effective cyclotron frequency

$${\omega }_{{\rm{c}}}=\frac{e{B}_{0z}^{{\rm{eff}}}}{{m}^{* }c}$$
(39)

and \({\hat{H}}_{{\rm{cor}}}\) is a “correction” term, which was explicitly computed in Ref. 38.

In absence of \({\hat{H}}_{{\rm{cor}}}\), the Hamiltonian (38) is simply the Landau-level Hamiltonian. Its spectrum corresponds to perfectly flat bands located at the energies \({E}_{n}=-(n+\frac{1}{2})\hslash {\omega }_{{\rm{c}}}\). In the vicinity of the magic angle, the corrections arising from \({\hat{H}}_{{\rm{cor}}}\) can be neglected for the first valence band (n = 0), as shown explicitly in ref. 38. Figure 2 shows that the second valence band (n = 1) is also rather flat (with a bandwidth smaller than 10 meV) at θ = 3.1∘, allowing us to neglect the corrections for this band as well. We therefore conclude that the original Hamiltonian (32) can be approximately mapped onto a Landau-level Hamiltonian, at least when the aim is to discuss physics involved with the first and second valence bands.

The mapping onto a Landau-level model allows us to analytically compute the contribution to the inter-band (single-spin and single-valley) optical conductivity due to transitions between the first and second valence bands. In Section V of ref. 46, we explicitly perform this calculation, starting from an expression similar to Eq. (23). Introducing the Landau-level eigenstates, we obtain

$${\sigma }_{\alpha \beta }^{{\rm{inter}}}(\omega )=\frac{{e}^{2}}{{m}^{* }{A}_{{\rm{M}}}}\frac{1}{{\omega }_{{\rm{c}}}^{2}-{(\omega +i\eta )}^{2}}\left(\begin{array}{cc}-i\omega &{\omega }_{{\rm{c}}}\\ -{\omega }_{{\rm{c}}}&-i\omega \end{array}\right),$$
(40)

where η = 0+. The Landau-level model therefore qualitatively explains the existence of a non-zero (single-spin and single-valley) Hall conductivity \({\sigma }_{xy}^{{\rm{inter}}}(\omega )\) as seen in Fig. 3b. In this context, \({\sigma }_{xy}^{{\rm{inter}}}(\omega )\) arises from the effective magnetic field \({B}_{0z}^{{\rm{eff}}}\) in the Hamiltonian (32).

We now turn to discuss how gapped plasmon modes emerge in this Landau-level picture. We can calculate the long-wavelength dispersion of inter-band plasmons from the diagonal elements of the conductivity (40). Imposing that the real part of the dielectric function (25) vanishes, we obtain

$${\rm{Re}}{\left[\varepsilon ({\boldsymbol{q}}\to {\boldsymbol{0}},\omega )\right]}_{{\boldsymbol{0}},{\boldsymbol{0}}}={\rm{Re}}\left[1+i\frac{2\pi }{q\bar{\varepsilon }}\frac{{q}^{2}{\sigma }_{{\rm{L}}}(\omega )}{\omega }\right]=0,$$
(41)

where \({\sigma }_{{\rm{L}}}(\omega )=2{\sigma }_{xx}^{{\rm{inter}}}(\omega )=2{\sigma }_{yy}^{{\rm{inter}}}(\omega )\). Here, the factor of 2 is a degeneracy factor that stems from the fact that the conductivity in Eq. (40) is a single-spin and single-valley one. Note that the Hall conductivity does not appear in the plasmon dispersion, because of the usual cancellation of the two valley-resolved contributions in the case that time-reversal symmetry is not broken. Using our result (40) for the conductivity, we find that the long-wavelength dispersion of collective modes in our system is given by the famous magnetoplasmon formula63

$${\omega }_{{\rm{pl}}}(q)=\sqrt{{\omega }_{{\rm{c}}}^{2}+{\omega }_{{\rm{2D}}}^{2}(q)},$$
(42)

where \({\omega }_{{\rm{2D}}}(q)=\sqrt{2\pi {n}_{0}{e}^{2}q/({m}^{* }\bar{\varepsilon })}\) is the usual 2D gapless plasmon dispersion in a parabolic-band electron liquid in zero external magnetic field41. Equation (42) only gives us the long-wavelength dispersion of the plasmon, since it stems from a local approximation to the optical conductivity tensor. Given the Landau-level-type Hamiltonian (38) one can go beyond the local limit, in principle, by computing the full non-interacting (bubble or) polarization function χ0(q, ω)41. This calculation is performed explicitly in Section V of ref. 46, where we show that it yields the same result as in Eq. (42) for the plasmon dispersion in the long-wavelength limit.

The analytical results (40) for the conductivity and (42) for the plasmon dispersion are textbook results41 for a 2D electron liquid subject to a homogeneous perpendicular magnetic field. Their main merit is to qualitatively explain the numerical results presented above. More precisely, Eq. (40) shows that the (single spin and single valley) Hall conductivity \({\sigma }_{xy}^{{\rm{inter}}}(\omega )\) is non-zero, in accordance with Fig. 3b. Equation (42) predicts a gapped plasmon mode, in accordance with Fig. 1.

However, a quantitative comparison reveals that the predictive power of the Landau-level model is quite limited. Let us first take a closer look at the gap in the plasmon dispersion. According to Eq. (42), in the Landau-level-mapping setting discussed in this Section, the plasmon gap at q = 0 coincides with the cyclotron frequency, i.e., ℏωpl(0) = ℏωc. Using the definition of \({B}_{0z}^{{\rm{eff}}}\) introduced above, i.e., \({B}_{0z}^{{\rm{eff}}}={\Phi }_{0}/{A}_{{\rm{M}}}\), we can estimate the cyclotron frequency (39) for small twist angles θ as following:

$$\hslash {\omega }_{{\rm{c}}}\approx \frac{2\pi {\hslash }^{2}}{{m}^{* }}\frac{{\theta }^{2}}{{a}^{2}},$$
(43)

Using the typical values for the parameters m* and a given above (for MoTe244), we find ℏωc ≈ 1.94θ2 meV, with θ expressed in degrees. For θ = 3.1∘, this implies ℏωc ≈ 19 meV. Comparing this value to the plasmon gap obtained analytically above in Eq. (30), i.e., ℏωpl(0) = ℏω01 ≈ 8 meV, we immediately see that they disagree by more than a factor of 2. The gap in the plasmon dispersion predicted by the Landau level model, i.e., Eq. (42), is far too large in comparison with the analytical gap (30) calculated above, which compares well with the full numerical results shown in Fig. 1.

The Landau-level model has another, more qualitative difficulty. The second term under the square-root sign in Eq. (42) is linear in the wave number q as in the case of the analytical calculation (30). The coefficients controlling these linear terms are however dramatically different. In the Landau-level model, this coefficient is determined by the f-sum rule41. on reasonable (30), it is instead related to the inter-band Berry connection46,74. Eq. (30) also takes into account the fact that the bands are not exactly flat.

In this article, we have calculated the optical properties and plasmonic spectrum of twisted homobilayer TMDs, focusing, for purely illustrative properties, on twisted MoTe2. Our main numerical results are shown in Figs. 1 and 3. Results in Fig. 1 are amenable to experimental testing via scattering-type near-field optical spectroscopy51,52,53,54,55,56. We have then gone one step further trying to understand the link between the real-space topological properties of the skyrmion Chern-band models introduced in refs. 38,39 and our numerical predictions. With the caveats discussed in the previous section, we have shown that the bulk plasmon dispersion in the long-wavelength limit roughly probes the uniform, i.e., spatially-constant, part of the skyrmion magnetic field (36). Real-space variations of the latter, which are encoded in the correction term \({\hat{H}}_{{\rm{cor}}}\) in Eq. (38), are expected to appear in the non-local corrections to the plasmon dispersion relation at finite q. In this work, we have included electron-electron interactions in the calculation of collective modes at the level of the RPA41, which well describes weakly correlated electron systems in the long-wavelength limit. In systems with flat bands and strong electron-electron interactions, however, one may want to calculate collective modes with non-perturbative approaches. One of these is certainly brute-force exact diagonalization75. Another one is the Bijl-Feynmann single-mode approximation (SMA), which was first used76 in the calculation of the phonon-roton spectrum of 4He. This approach was successfully generalized to the case of strongly correlated fractional quantum Hall fluids by Girvin, MacDonald, and Platzman (GMP)65,77 leading to the prediction of gapped magneto-rotons in the lowest Landau level. The GMP-SMA has also been applied to the present problem78. We emphasize, however, that a qualitative difference between fractional quantum Hall fluids and interacting Chern insulators in twisted MoTe2 needs to be taken into account. In the standard GMP theory65,77, Landau-level mixing can indeed be safely neglected. In the present problem, however, Landau-level mixing is not negligible. In fact, the electron-electron interaction energy scale Vee greatly exceeds the effective cyclotron energy. Indeed the length scale \(\sqrt{{A}_{{\rm{M}}}}\) as a typical value of the inter-electron distance, we find:

$$\frac{{V}_{{\rm{ee}}}}{\hslash {\omega }_{{\rm{c}}}} \sim \frac{\sqrt{{A}_{{\rm{M}}}}}{{a}_{{\rm{B}}}^{* }}\approx \frac{65}{\bar{\varepsilon }},$$
(44)

where \({a}_{{\rm{B}}}^{* }=\bar{\varepsilon }{\hslash }^{2}/({m}^{* }{e}^{2})\) is the material’s Bohr radius. The last approximation holds for the choice of parameters we made throughout this work, i.e., a = 0.352 nm, m* = 0.6me, and θ = 3. 1∘.

In the future, it will be interesting to study metallic magnetic materials79,80, where itinerant electrons coexist with real-spin topological Skyrme textures81,82. As suggested by the present results, bulk plasmons could also serve in this case as a very useful probe of these interesting real-space topological structures. The strong plasmon-magnon interaction recently predicted in 2D magnetic materials83 points towards this direction.

Methods

The band structure and DOS results presented in Fig. 2b–d have been carried out by employing a plane-wave expansion of the Hamiltonian in Eq. (9). The chosen basis consists of 37 plane waves at each wave vector k in the moiré BZ (corresponding to the first three hexagonal shells spanned by the moiré reciprocal lattice vectors). This results in a total of 148 valence bands. We have then retained 30 of them for the evaluation of the sums over the band indices appearing in Eqs. (21) and (23).

For the optical conductivity (Fig. 3), numerical results have been obtained by performing the integrals over a mesh of 90 × 90 = 8100 equally-spaced points in the moiré BZ. The value of the small parameter η appearing in Eqs. (21) and (23) were taken to be η = 0.25 meV.

The loss function (Fig. 1) at small values of q has been obtained by combining Eqs. (24)–(26) into the following expression:

$${\mathcal{L}}(q,\omega )\approx -{\rm{Im}}\left[\frac{1}{1+\frac{2\pi iq}{\omega \bar{\epsilon }}{\sigma }_{xx}(\omega )}\right].$$
(45)

In the range of values of q between q = 0 nm−1 and q = 0.08 nm−1 we have included 150 points.