In analogy to Lifshitz transitions in electronic systems, topological transitions have recently attracted widespread attention in photonic metamaterials, metasurfaces, and two-dimensional materials, enabling exotic regimes for light-matter interactions. Here, we discuss and study enhanced photonic local density of states in twisted hyperbolic bilayers, enabled by topological transitions emerging at specific twist angles. Our results enhance the understanding of nanoscale light-matter interactions in stacked optical materials as they are rotated with respect to each other in the context of twistronics and suggest emerging applications of these concepts for photonics, including for the manipulation of radiative heat transfer and the control and harvesting of light at the nanoscale.

Nano-optics studies the interaction of light with nanoscale systems.^{1} In the context of light emission, the nanoscale environment significantly affects the behavior of quantum emitters, and their spontaneous emission rate (SER) can be largely modified by engineered nanophotonic systems, as quantitatively described by the photonic local density of states (LDOS). There has been a significant surge of interest in developing new classes of nanophotonic systems that offer enhanced LDOS. Of particular interest is the emerging class of hyperbolic materials, in which at least two components of the permittivity tensor along specific directions have different signs of their real part.^{2–7} Light dispersion in these materials follows peculiar hyperbolic bands, hence their name, with an open topology, and therefore it can sustain very long wavevectors, corresponding to short wavelengths that boost light matter interactions. Such exotic features in turn facilitate many exciting applications such as hyperlenses with deeply subwavelength resolution,^{8–10} broadband enhancement of photonic LDOS,^{11,12} and enhanced thermal management,^{13,14} just to name a few. To realize hyperbolic materials, a typical approach consists in multiple vertical stacks of metal and dielectrics, yielding out-of-plane hyperbolicity.^{4,9} Metallic or graphene subwavelength gratings can support uniaxial metasurfaces with in-plane hyperbolicity,^{6,15,16} which may be more preferable because they are easier to fabricate, have a thinner profile, and their optical interactions are more easily accessible and less prone to losses since they arise at the interface with dielectrics and not in the bulk. Recently, polar van der Waals nanomaterials with strong contrast of covalent binding and van der Waals bindings have been shown to naturally exhibit hyperbolicity, both out-of-plane and in-plane,^{17} relaxing the need for complex fabrication schemes, in which the minimum features play a major role in fundamentally limiting the overall hyperbolic responses.^{18,19} These naturally hyperbolic materials^{20,21} open totally new opportunities to engineer nanoscale light-matter interactions.

Reconfigurable hyperbolic responses through dispersion engineering have become particularly interesting in this context, as they provide added flexibility in controlling light-matter interactions. For example, varying the filling ratio of alternating dielectric and metal thin layers can turn on and off the hyperbolic response of the resulting metamaterials;^{22} nanostructuring out-of-plane hyperbolic materials can induce in-plane hyperbolicity;^{23,24} stacking hyperbolic metasurface (HMTS) bilayers with a twist angle between them can drastically tune the dispersion contours, transitioning from hyperbolic to elliptical as a function of the rotation angle.^{25–30} In these systems, the emerging hyperbolic-to-elliptic optical transitions, inherently topological in nature, are analogous to Lifshitz transitions^{31} in electronics. In his seminal paper published in 1960, Lifshitz found that the electronic properties of materials strongly depend on the (closed or open) topology of the associated Fermi surface, experiencing an anomaly and a discontinuity of electron density of states^{31} when such topology changes. The photonics analogue corresponds to a modification of the topology of the isofrequency contours of propagating electromagnetic modes, which play the role of Fermi surfaces. At such transition, the photonic LDOS also exhibits an “anomaly,” i.e., a large change in the frequency spectrum at the topological transition of hyperbolic metamaterials^{22} and HMTS.^{16}

In this paper, we study the evolution of the photonic LDOS at topological transitions in twisted HMTS bilayers, occurring at a specific rotation angle known as the photonic magic angle,^{26,27} in analogy with the magic angle in twisted graphene bilayers. In electronic systems, flat bands and dissipation-less transport of electrons (superconductivity) were shown to be supported at a magic rotation angle between two graphene monolayers stacked on top of each other,^{26,27} opening the field of twistronics. The analogue in photonic systems corresponds to the flat dispersion lines and associated diffraction-free propagation of light at the photonic magic angle. Both phenomena share the same origin of the modification of interlayer coupling and hybridization of wavefunctions controlled by the rotation angle. These twisted bilayer systems have received tremendous attention in recent years for on-demand dispersion engineering of light at the nanoscale. Here, we aim at quantitatively exploring the role that these largely reconfigurable extreme optical responses have on the LDOS and demonstrate an anomaly of the photonic LDOS of a localized emitter near twisted metasurface bilayers at the photonic magic-angle topological transitions.^{26,27} In the following, we first derive the dyadic Green's functions of the system. We then numerically study the manipulation of the spontaneous emission rate (SER) at the photonic magic-angle topological transition and discuss promising applications. Our results enable new degrees of freedom to manipulate nanoscale light-matter interactions for various applications enabled by the concept of twistronics for light.^{32}

The geometry under analysis consists of a twisted pair of hyperbolic surfaces that behave as metallic (inductive) for electric fields along one characteristic direction and as dielectric (capacitive) along the orthogonal one. Their implementation can be achieved in various geometries. First, this response can emerge naturally in two-dimensional (2D) black phosphorous^{33} and polar van der Waals (vdW) materials like *α*-MoO_{3} (Refs. 27 and 34–36) and MoTe_{2}.^{37} When the thickness is deeply subwavelength, vdW materials can be treated as surface conductivity sheets, realizing natural hyperbolic surfaces.^{27,34} Another avenue for these in-plane hyperbolic responses, offering more flexibility in terms of the frequency range of interest, is to structure thin materials, such as silver gratings,^{15} graphene nanoribbon arrays,^{16,25,26} and hBN gratings,^{23,24} in which the extreme anisotropy enabled by the pattern geometry showcases hyperbolic responses. In this paper, in order to model twisted hyperbolic bilayers, we study deeply subwavelength graphene nanoribbon arrays and stack them with a twist angle with respect to each other. For simplicity, we use the effective homogeneous surface conductivity model for these nano-arrays neglecting their granularity, which is accurate if the periodicity is smaller than the working wavelength, and the spacing between the two surfaces is larger than their transverse period. Although we use graphene HMTSs, our results can be easily extended to other platforms and do not lose generality.

We adopt an *e ^{-iωt}* time convention, and the global coordinate system is defined such that the

*x*and

*y*directions are perpendicular and parallel to the nanoribbons of the top HMTS. The bottom HMTS has a relative rotation angle ($\theta 2=\Delta \theta $) around

*z*. The HMTSs are located at $z=z1=0$ and $z=z2=\u2212d$, respectively. For simplicity, we assume that the background is air, but our findings can be easily extended to asymmetric backgrounds. As previously reported, the bilayer dispersion can be found by solving the source-free Maxwell equations.

^{26,27}The Fresnel coefficients for different incidence angles will be a useful tool to model the overall response in the presence of an electric dipole source. We assume incident waves with wavevector $k\u2192inc=kx,ky,\u2212kz$ and correspondingly reflected waves with $k\u2192r=kx,ky,kz$. For simplicity, an auxiliary coordinate system ($\rho \u0302,h\u0302,z\u0302$) is defined as

where $k\rho =kx2+ky2$ and $kz=\xb1k02\u2212k\rho 2$. Therefore, the electric field for *s*-waves (transverse electric field, TE) can be written as $s\xb1=h\u0302$, while *p*-waves (transverse magnetic field, TM) have $p\xb1=1k0\xb1kz\rho \u0302\u2212k\rho z\u0302$. Here, the positive sign indicates that the sign of $kz$ is positive (backward, compared to the incident wave), and the negative sign corresponds to the sign of *k _{z}* being negative (forward, compared to the incident wave). The magnetic field can be found as $H\u2192=\u2212i\omega \mu 0\u2207\xd7E\u2192$.

Due to the large anisotropy, the waves in all the dielectric layers have both TE and TM components, independent of the incident polarization (denoted as $E\u2192inc$), as shown in Fig. 1(b), due to cross-polarization conversion between TE and TM modes. To obtain the Fresnel coefficients, we follow the procedure in Refs. 25–27 and derive the equation

Here, *a*_{1} (*a*_{2}) and *b*_{1} (*b*_{2}) denote the coefficients of TE and TM modes for reflected (transmitted) waves, respectively. The weight of TE and TM incidence is denoted as $\delta s$ and $\delta p$, respectively. When $\delta s=0$ and $\delta p=1$, it corresponds to the TM incidence, while the case of $\delta s=1$ and $\delta p=0$ corresponds to the TE incidence. The surface conductivity tensor is then projected onto the ($\rho \u0302,h\u0302,z\u0302$) coordinates:

where $\theta v$ is the rotation angle and $v=1,2$ denotes the top and bottom HMTS, respectively, $R\u033f\theta v$ is the rotation matrix, $\sigma xx=\xi \sigma G$ and $\sigma yy=\xi \sigma G\sigma C\sigma C+1\u2212\xi \sigma G$, where the factor $\sigma C=\u22122i\omega \epsilon 0P\pi ln1sin1\u2212\xi \pi /2$ takes into account nonlocality stemming from the near-field coupling between neighboring graphene nanoribbons, and it can be obtained analytically by mapping the boundary conditions in the electrostatic approximation,^{38}^{,}$\xi =W/P$ is the geometric ratio of the ribbon width (*W*) divided by the periodicity (*P*), and $\sigma G$ is the conductivity of graphene.^{39} We also define the open angle of each HMTS $\alpha =arctan\u2212\sigma xx\sigma yy$, which is the angle between the asymptotic dispersion line of hyperbolic modes supported by each surface individually and the *x* axis, as plotted in [Figs. 1(c) and 1(d)].^{26} For different geometries, the open angle changes [Fig. 1(c)] with frequency with a different disersion, suggesting large turnability in molding the hyperbolic propagation of graphene plasmons in the bilayer system as we describe in the following.

The Fresnel coefficients of the bilayer can be found using Cramer's rule:

^{26,27}In this scenario, we can further simplify the dispersion using a 2 × 2 block matrix, which provides a clear picture of mode hybridization and the role of evanescent coupling between the two HMTSs in controlling the overall dispersion.

^{26}

Equations (2)–(4) go well beyond the eigen-modes of the bilayer and provide the full response of the surface for arbitrary incident waves. Therefore, based on this result, we can calculate the dyadic Green's function of the system as well as the associated photonic LDOS, which measure the overall strength of light-matter interactions for a localized emitter. The reflection dyadic Green's function at the location of the source ($r\u2192$) can be written as^{1,40,41}

Here, the various dyadic tensors can be written as $Ms+s\u2212=s+\u27e9\u2297s\u2212\u27e9$; $Mp+s\u2212=p+\u27e9\u2297s\u2212\u27e9$; $Ms+p\u2212=s+\u27e9\u2297p\u2212\u27e9;$ and $Mp+p\u2212=p+\u27e9\u2297p\u2212\u27e9$, where $\u2297$ denotes the tensor product. This dyadic Green's function denotes the secondary electromagnetic emission from the bilayer at the source location, which can be used to evaluate the LDOS at $r\u2192$. As illustrated in Fig. 1, we can calculate the radiation of a two-level quantum emitter, modeled as an electric dipole, located at $z=D$, and oriented along $n\mu =nx,ny,nz$ with $nx2+ny2+nz2=1$. The spontaneous emission rate (SER) can be obtained as^{1}

while the photonic LDOS is directly proportional to $Im(G\u033fRr\u2192,r\u2192,\omega )$ and thus it is proportional to the SER, if SER is much larger than unity.

Based on this general model and theoretical analysis, we can now explore the response of small emitters placed near twisted HMTS bilayers as they undergo topological transitions at the photonic magic angle. In general, two types of photonic magic angles can exist, corresponding to two different scenarios of topological transitions in momentum space:^{26} hyperbolic to elliptical transitions (hereafter referred to as type I magic angle) and hyperbolic to weak-coupling transitions (hereafter referred to as type II magic angle). Taking into account the open angle of the hyperbolas and the rotation angle of the twisted bilayer, it is possible to predict based on basic geometrical considerations when a magic angle arises. In particular, we can expect a type I magic angle at *π*-2*α*, under the condition that *π*/4<*α* < *π*/2; while a type II magic angle arises at the rotation of 2*α*, when *α* < *π*/4. Figures 2(a)–2(c) showcase how these transitions emerge at these specific angles. The black and blue dashed lines show the dispersion lines of each metasurface when isolated. At their crossing points in momentum space, strong coupling between the modes can arise for sufficiently close surfaces, manifested in an anti-crossing in the band diagram of the coupled system (red lines). Interestingly, the number of anti-crossing points (N_{ACP}) can be simply determined by geometrical arguments, knowing the open angle of the hyperbolas and the rotation angle, as summarized below:^{26}

In an intuitive picture, the anti-crossing will break the dispersion lines of each metasurface, and “broken” dispersion lines will rejoin to form the coupled dispersion lines, supporting closed elliptical dispersion when N_{ACP} = 4 [Fig. 2(a)], open hyperbolic dispersion when N_{ACP} = 2 [Fig. 2(b)], and co-existing, weakly interacting hyperbolic dispersion in a weakly coupled system when N_{ACP} = 0 [Fig. 2(c)]. In this sense, it is easy to recognize that the magic angles are the rotation angles for which the topological invariant N_{ACP} changes, and the two types of photonic magic angles correspond to $\Delta \theta =\pi \u22122\alpha $ (if $\alpha >\pi /4$) with hyperbolic to elliptical topological transitions (type I) and $\Delta \theta =2\alpha $ (if $\alpha <\pi /4$) with bi-hyperbolic to hyperbolic topological transitions (type II). Below, we study the “anomalies” in photonic LDOS emerging at these two types of topological transitions.

We start by studying the effect of enhanced light-matter interactions in a geometry with $\xi =0.8$, *P* = 30 nm, and *d* = 10 nm. Figure 2(d) shows the logscale contour plot of SER with respect to frequency (which controls the open angle *α*) and rotation angle (Δ*θ*), calculated via Eq. (6). From Figs. 2(d) and 2(e), we see a resonant feature near the type I and type II photonic magic angles, where the topological transitions happen. To better illustrate such features, we have shown in Fig. 2(f) the spectrum of SER for tBL HMTS with $\Delta \theta =25\xb0$ and $\Delta \theta =60\xb0$, in comparison with monolayer graphene and graphene HMTSs. Anomalous resonances near 6 THz and 14.5 THz can be observed for $\Delta \theta =25\xb0$ and $\Delta \theta =60\xb0$, respectively, which correspond to type I photonic magic-angle topological transitions, while the spectrum for monolayer graphene or monolayer HMTS are smooth in this spectrum. For $\Delta \theta =60\xb0$, we can observe an additional anomaly around 26 THz [see inset in Fig. 2(c)], due to the type II photonic magic angle. Therefore, these results suggest an anomaly closely related to the photonic magic-angle topological transitions, where the light-matter interaction strength also changes. We also note that at around 32 THz, where the range of hyperbolicity of each individual metasurface ends, or rather topological transitions of individual HMTSs arise [Fig. 1(c)], SER has also a dramatic reduction, which can be seen from Figs. 2(e) and 2(f).

Next, we investigate the flexibility in controlling light-matter interactions near the anomalies of photonic magic angles afforded by the twisted bilayers. As shown in the literature,^{42} the dispersion of graphene plasmons in HMTSs strongly depends on the geometrical parameters (including periodicity and duty ratio), the doping or applied voltage (Fermi levels), and the dielectric background. Here, for simplicity, we focus on the case $P=30\u2009nm$ and $\xi =0.9$. Comparing Figs. 2(d) and 3(a), we see that for larger duty ratio, the hyperbolic frequency of monolayer HMTSs becomes smaller, so does the resonant frequency of SER. If we change the chemical potential of graphene, the topological transition as well as the associated resonant SER also shifts [Figs. 3(a) and 3(b)]. This can be explained by the change of surface conductivity and dispersion of individual graphene HMTSs. Figure 3(c) shows the surface conductivity component perpendicular to the nanoribbons, which clearly shows the shift of the hyperbolic range as the function of the chemical potential and duty ratio. Moreover, Fig. 4 summarizes the change of SER with respect to the distance (*d*) between two HMTS. As the distance becomes larger [Figs. 4(a) and 4(c)], the two HMTSs couple much more weakly to each other and the dispersion distortion due to the anti-crossings becomes less significant. Therefore, the “anomalies” become less obvious. But for smaller distances, shown in Figs. 4(c) and 4(d), the responses are largely enhanced and well visible. Nevertheless, a smaller distance of the source may require finer granularity in the metasurface fabrication to ensure that our homogenized model still works and poses a trade-off between enhanced resonant SERs and fabrication features.

In summary, we have revealed extreme control of photonic LDOS and SER of small optical emitters placed near twisted HMTS bilayers operated near the photonic magic-angle topological transition. This demonstrated response is analogous to the anomalies observed in electronic systems when the Fermi surface changes from hyperbolic to elliptical at a Lifshitz transition. Our results indicate that exotic light-matter interactions can emerge near the photonic magic-angle topological transitions, of great interest for a plethora of applications, including light harvesting, Förster energy transfer, nonlinear optics, and many others. For example, we envision the possibility of manipulating the dipole-dipole interactions in such systems and enhance near-field radiative heat transfer at the nanoscale due to exotic directional energy transport and anomaly-enhanced light-matter interactions, which we will explore in the future. Besides our specific study of twisted HMTS bilayers, future work can also focus on light-matter interactions at substrate-induced topological transitions^{43,44} and in twisted trilayers. More broadly, our findings highlight the interesting opportunities enabled by twistronics for photons, a field that holds the promise to provide exciting potentials for nanophotonics and optoelectronics.

The authors acknowledge funding in part from the Office of Naval Research (Grant No. N00014-19-1-2011), the Air Force Office of Scientific Research, the National Research Foundation Singapore (No. CRP22-2019-0006), and the Advanced Research and Technology Innovation Centre (ARTIC, No. R-261-518-004-720).

## DATA AVAILABILITY

The data that support the findings of this study are available within this article.

## References

_{3}and polar dielectric interfaces