Introduction

Moiré patterns formed in twisted bilayer two-dimensional van der Waals materials have led to the observation and exploration of emerging electronic properties1,2,3,4,5, including moiré excitons6, fractional Chern insulators7, and Mott insulators8,9. Along with the surge of research into moiré bilayers of two-dimensional materials, moiré patterns in photonics have also drawn attention from the optical community10,11,12,13,14, with demonstrations of twist-induced flat bands15, energy localization16, and optical soliton17. Moreover, the twist of photonic structures breaks the mirror symmetry and enables chiral optical properties. Here we show that the structural chirality of twisted bilayer photonic crystal (TBPC), through the bound states in the continuum (BICs) in both layers of photonic crystals (PhCs), enables the chirality of light: the optical vortex.

An optical vortex is a light beam with a spiral phase front and an undefined phase at the beam center, resulting in a “doughnut” intensity profile. Optical vortices have been shown to carry orbital angular momentum (OAM), which is a new degree of freedom beyond the spins of light. Since the discovery of light’s OAM18, optical vortices have greatly expanded the horizon of numerous research fields19,20,21,22,23,24, including super-resolution imaging25,26, optical communication27,28,29, micromanipulation30,31, and quantum information processing32,33. As a result, optical vortex generation has attracted extensive attention.

In this work, we discover that the twisted stacking of two-layer PhCs enables the generation of robust optical vortex radiation that is insensitive to either the incident angle or the illumination position of the incoming light. By taking advantage of the special interlayer coupling mechanism between the two twisted PhC slabs, a twist-enabled channel for energy transfer between free space and the bound states in the continuum (BICs)34,35,36,37 in the monolayer PhC is created. Upon arbitrary Gaussian beam illumination onto the TBPC, at- \(\Gamma\) BIC mode in PhC 1 is excited by its coupled guided resonances in PhC 2 via interlayer channeling. The BIC mode in turn excites the radiative guided resonances in PhC 2 for optical vortex emission by inducing Pancharatnam-Berry (PB) phases, as we will show in more detail below. We theoretically and numerically demonstrate the “doughnut” intensity distribution and vortex phase distribution in the generated beam from the TBPC system with a high Q factor, which makes TBPC a promising platform for vortex micro-laser applications. Moreover, optical vortex radiation with both even and odd order and at any wavelength can be realized by adjusting the TBPC design.

Results

Methodology

Figure 1a shows a schematic diagram of the system considered in this work. The device is constructed by stacking two planar PhC slabs with a small twist. Each slab is a honeycomb array of silicon nanodisks. Below we show that upon excitation by an arbitrary pulsed Gaussian beam, the at-Γ BIC modes (herein referred to as BIC for simplicity) in both PhC slabs will be excited, which further radiate energy to the +z and −z directions. Such counterintuitive excitation of BIC modes by Gaussian beam and radiation from BIC modes are realized by a unique interlayer channel created by the moiré structure. Specifically speaking, in a monolayer PhC slab with twofold rotational symmetry, radiative resonant modes can couple to the far-field with frequencies and in-plane wave vectors matching the plane waves in free space. Those radiative resonant modes are usually called guided resonances38. Meanwhile, any singlet at Γ point is an exception since the rotational symmetry of any plane wave with in-plane wave vector \(\left({k}_{x},{k}_{y}\right)=0\) mismatches the rotational symmetry of that at Γ singlet34. Thus, such a singlet is isolated from the free space, and we call this resonant mode which is decoupled from the free space a symmetry-protected BIC. Even with another identical PhC slab introduced to form a bilayer PhC system, the at-Γ mode will remain a BIC mode, due to the protected two-fold symmetry.

Fig. 1: The twisted bilayer photonic crystal (TBPC) system.
figure 1

a Schematic of TBPC with optical vortex emitted from the AA stacking regions to both +z and −z directions. The height and diameter of the disks are 220 nm and 450 nm, respectively. The lattice constant of monolayer honeycomb PhC is 1.0 μm. b Schematic view of the twist-enabled energy transfer between free space light beams and at-Γ BIC modes through guided resonances.

In bilayer PhC systems, when the two PhC layers are perfectly aligned with zero twist angle, the periodicity remains the same as that of monolayers. Thus, the at-Γ BIC modes in both layers are still isolated from guided resonances due to the phase mismatch mechanism. However, the introduction of a small twist between the two layers generates a moiré structure1, which breaks the short-range periodicity of a monolayer PhC. Since the rigorous phase match mechanism rely on the short-range periodicity, it will be broken by introducing moiré structure to the PhCs. The BIC modes in one slab can thus couple to the guided resonances in the other slab which is coupled to the free space15, as shown in the following paragraphs. This interlayer coupling mechanism enabled by the twists in TBPC systems allows energy transfer between BIC modes and the free space via guided resonances, so it is named as a “moiré channel” in this work (Fig. 1b). Moreover, the states of polarization (SOPs) of far-field radiation from these guided resonances are momentum-space variant, and they will form a polarization vortex around the BIC in the momentum space34,35,36. Hence, when BIC modes radiate energy via these guided resonances enabled by the moiré channel, vortices emerge in the real-space radiation.

Moiré channel in TBPC

Here we prove the existence of the moiré channel, i.e., with non-zero twist angles, the coupling between the BIC in one layer and the guided resonances in the other layer is finite. Coupled-mode theory is used as it has been proven to be accurate in analyzing the photonic behaviors of a TBPC system15. Our theory begins with a well-defined transverse electric (TE) mode in the silicon disk15. When two disks are placed closely, the crosstalk between the eigenmodes in either disk is described by:

$$\left\{\begin{array}{c}\frac{{{{{{\rm{d}}}}}}{a}_{1}\left(t\right)}{{{{{{\rm{d}}}}}}t}=(i{\omega }_{1}-{\kappa }_{1}){a}_{1}+i{g}_{12}{a}_{2}\\ \frac{{{{{{\rm{d}}}}}}{a}_{2}\left(t\right)}{{{{{{\rm{d}}}}}}t}=(i{\omega }_{2}-{\kappa }_{2}){a}_{2}+i{g}_{21}{a}_{1}\end{array}\right.$$
(1)

where a, ω, and κ are the mode amplitude, the angular eigenfrequency, and the decay rate, respectively. Note that ω1 = ω2 = ω0, κ1 = κ2 = κ0 for identical disks. g12 and g21 are the coupling coefficient between the two modes. Without loss of generality, we set \({g}_{12}={g}_{21}=g\)15.

For a monolayer PhC consisting of such a disk array, the traveling mode with the wave vector k (mode k in short) can be denoted as a linear combination of each disk modes, and the amplitude \({a}_{k}\left(t\right)\) is: \({a}_{{{{{{\boldsymbol{k}}}}}}}=\frac{1}{\sqrt{N}}\mathop{\sum}\limits_{j}{{\exp }}(i{{{{{\bf{k}}}}}}\cdot {{{{{{\bf{r}}}}}}}_{{{{{{\rm{j}}}}}}}){a}_{j}\), where N is the total number of the disks, and rj is the real-space position of the center of disk j. We set the coupling intensity between BIC in PhC 2 and mode k1 in PhC 1 is ζ(k1). Specifically speaking, mode k1 is a BIC mode when k1=0. We can derive that for non-twist bilayer PhCs, inside the first Brillouin zone of the monolayer, ζ(k1) is always zero except k1 = 0 (Γ point), which means BIC modes in PhC 2 only couple to BIC modes in PhC 1, as we can see in Fig. 2a.

Fig. 2: Interlayer coupling features in bilayer PhCs for (a) perfectly aligned PhCs and (b) twisted bilayer photonic crystals (TBPCs).
figure 2

The color of the upper plane denotes the absolute value of coupling intensity between the BIC in PhC 2 and the guided resonances with wave vector \(({k}_{x1},{k}_{y1})\) in PhC 1. Note that there is no coupling between the BIC in PhC 2 and guided resonances in PhC 1 in a.

However, in a TBPC system, ζ(k1) can be denoted as:

$$\zeta \left({{{{{{\bf{k}}}}}}}_{1}\right)=\frac{i}{N}\mathop{\sum}\limits_{j}{{\exp }}\left(i{{{{{{\bf{k}}}}}}}_{1}\cdot {{{{{{\bf{r}}}}}}}_{1{{{{{\rm{j}}}}}}}\right)\cdot g\left({{{{{{\bf{r}}}}}}}_{1{{{{{\rm{j}}}}}}}\right)$$
(2)

where N is the total number of the disks in one layer, \({{{{{{\bf{r}}}}}}}_{1{{{{{\rm{j}}}}}}}\) is the real-space position of the center of disk j in PhC 1, and \(g\left({{{{{{\bf{r}}}}}}}_{1{{{{{\rm{j}}}}}}}\right)\) is the coupling intensity between disk j in PhC 1 and its nearest neighbor in PhC 2. Since the moiré structure differs the coupling intensity \(g\left({{{{{{\bf{r}}}}}}}_{1{{{{{\rm{j}}}}}}}\right)\) for different j, the coupling strength ζ(k1) is non-zero at non-Γ points. Hence, BIC modes in PhC 1 can couple to the guided resonances in PhC 2, and vice versa, as shown in Fig. 2b. Moreover, we prove that the coupling coefficient ζ(k1) is nearly identical for same |k1| points (Supplementary Fig. 2), which means that the BIC in one layer can excite the guided resonance modes around BIC in the other layer with roughly the same phase and intensity, as illustrated in Fig. 2b.

Optical vortex generation in TBPC

In the following paragraph, we will analytically prove the optical vortex generation, which is radiated from the guided resonances excited by the BIC mode. Since honeycomb lattice is protected by the \({C}_{2}^{Z}T\) symmetry, the SOPs in the far-field will be linearly polarized34, oriented in the direction \(\theta \left({{{{{{\bf{k}}}}}}}_{1}\right)\) to the x axis. To get the explicit expressions of the radiated field, we start from the temporal coupled-mode theory (TCMT) main equations:

$$\left\{\begin{array}{c}\frac{{{{{{\rm{d}}}}}}A({{{{{{\bf{k}}}}}}}_{1})}{{{{{{\rm{d}}}}}}t}=[i\omega ({{{{{{\bf{k}}}}}}}_{1})-\kappa ({{{{{{\bf{k}}}}}}}_{1})]A({{{{{{\bf{k}}}}}}}_{1})+\zeta ({{{{{{\bf{k}}}}}}}_{1}){A}_{{{{{{\rm{BIC}}}}}}}\cdot \exp (i{\omega }_{{{{{{\rm{BIC}}}}}}}t)\\ |{E}_{{{{{{\rm{out}}}}}}}\rangle={{{{{\bf{D}}}}}}A\end{array}\right.$$
(3)

where \(A({{{{{{\bf{k}}}}}}}_{1})\), \(\omega ({{{{{{\bf{k}}}}}}}_{1})\), and \(\kappa ({{{{{{\bf{k}}}}}}}_{1})\) are the amplitude, the eigenfrequency, and the decay rate (due to radiative loss) of mode k1 in PhC 1, respectively. ABIC and ωBIC are the amplitude and the resonant frequency of the BIC mode in PhC 2, respectively. \(\zeta \left({{{{{{\bf{k}}}}}}}_{1}\right)\) is the coupling intensity between mode k1 in PhC 1 and the BIC mode in PhC 2 calculated in Eq. (2). |Eout〉 is the outgoing wave amplitude, and D is the coupling coefficient between the resonance and the outgoing plane wave. The orientation of the SOP should be24:

$$\left(\begin{array}{c}{d}_{x}\\ {d}_{y}\end{array}\right)=\sqrt{{d}_{x}^{2}+{d}_{y}^{2}}\cdot \left(\begin{array}{c}{{\cos }}\left(\theta \left({{{{{{\bf{k}}}}}}}_{1}\right)\right)\\ {{\sin }}\left(\theta \left({{{{{{\bf{k}}}}}}}_{1}\right)\right)\end{array}\right)$$
(4)

where \(\theta \left({{{{{{\bf{k}}}}}}}_{1}\right)\) is the angle of the polarization vector. By dividing this Jones vector in the basis of left- and right-handed circular polarized (LCP and RCP) unit vectors, the radiated field |Eout〉 of mode k1 can be formulated as:

$$|{E}_{{{{{{\rm{out}}}}}}}\rangle=K\zeta ({{{{{{\bf{k}}}}}}}_{1}){e}^{-i\theta ({{{{{{\bf{k}}}}}}}_{1})}|L\rangle+K\zeta ({{{{{{\bf{k}}}}}}}_{1}){e}^{-i\theta ({{{{{{\bf{k}}}}}}}_{1})}|R\rangle$$
(5)

where K is an integrated parameter, which is nearly a constant if we choose a circular loop around the Γ point24. |L〉, |R〉 denote the LCP and RCP unit vectors, respectively. From Fig. 2b, we notice that \(\zeta \left({{{{{{\bf{k}}}}}}}_{1}\right)\) is also roughly a constant in the circular loop around the Γ point. Thus, \(K\zeta \left({{{{{{\bf{k}}}}}}}_{1}\right)\) shall not be able to introduce another phase factor. From Eq. (5) we notice that the emitted light will be composed of two non-trivial parts: LCP and RCP. Both parts of the light will gain a PB phase depending on the orientation of the SOP of the guided resonance: \({\theta }_{{{{{{\rm{LCP}}}}}}}=-\theta \left({{{{{{\bf{k}}}}}}}_{1}\right)\) and \({\theta }_{{{{{{\rm{RCP}}}}}}}=\theta \left({{{{{{\bf{k}}}}}}}_{1}\right)\), where θLCP and θRCP denote the phase of the LCP and RCP portion, respectively, and \(\theta \left({{{{{{\bf{k}}}}}}}_{1}\right)\) is the angle of the polarization vector. Since BIC is the crossing of nodal lines for x direction and y direction portions of the far-field radiation, the SOPs will present a vortex around the BIC34,35,36. Thus, the angle of polarization vector \(\theta \left({{{{{{\bf{k}}}}}}}_{1}\right)\), as well as θLCP and θRCP, will form a vortex structure around the \({{{{{{\bf{k}}}}}}}_{1}{{{{{\boldsymbol{=}}}}}}{{{{{\bf{0}}}}}}\) point. Hence, the emitted beam will gain the desired spiral phase front, whose topological order is the same as that of the BIC mode and can be described as l = ±q. Here q is the topological order of the BIC, minus and plus signs correspond to LCP and RCP portions, respectively.

To verify the above theory, we perform a full-wave numerical simulation of a representative TBPC with a 300 nm separation and a 1.5° twist angle. We excite the system with a pulsed Gaussian beam and measure the field profile on a plane centered at the AA point and is parallel to the TBPC surface after the pulse is over. The AA point denotes the center of AA stacking region where the top and bottom honeycomb PhC layers are well aligned. Simulation details can be found in the Method section. Two peaks are identified in the mode intensity spectrum (Fig. 3a), which are marked as “i” and “ii”. The mode profile in the central cross-section in PhC 1 is shown in Fig. 3c, which matches the BIC mode profile simulated using a monolayer PhC (Fig. 3b). Then we examine the radiated beams corresponding to the two peaks (“i” and “ii”) in the far field. As illustrated by Fig. 3d, e, both the “doughnut” shape in the amplitude profile and the clear vortex in the phase distribution diagram serve as clear evidence of OAM existence. The topological order l equals to −2 and +1 for the RCP portions of peak “i” and peak “ii”, respectively, which is consistent with the theoretical results of the SOPs around the BIC: the lowest topological order of the at-Γ BIC in a honeycomb lattice could be either −2 or 124. Hence, all the full-wave numerical simulation results not only justify our theoretical model, but also directly proves the emission of the optical vortex from the AA region and further demonstrate the important role of BIC mode in the OAM generation. Compared with previous micro-optical vortex generators24, the optical vortex generated by TBPC is not limited to even numbers.

Fig. 3: Optical vortex generation in a representative twisted bilayer photonic crystal (TBPC) system.
figure 3

a The mode intensity spectrum of the TBPC system. The two peaks in the spectrum are marked as “i” and “ii”, with the frequency of 193.4 THz and 211.2 THz, and the Q-factor of 3.1 × 103 and 3.4 × 103, respectively. b Electric-field profile of the BIC modes hosted in monolayer PhC. c Electric-field profiles in PhC 1 of the TBPC system at peak “i” and “ii”. d, e The electric-field distribution of the optical vortex emission 12 μm away from the TBPC surface at peak “i” (d) and “ii” (e). The upper panels depict the normalized amplitude of the LCP and RCP portions of the emission, while the lower panels are the phase distribution of the corresponding upper panels in dashed white circles.

Furthermore, we change the interlayer separation and the twist angle between the two layers. Since the interlayer coupling between two PhC slabs does not change the SOPs around BIC in the monolayer, the topological orders of the emitted vortices will remain the same. We perform a numerical calculation based on the simulated system for Fig. 3 but change its interlayer separation and twist angle following Fig. 4a. For tidiness, only the LCP portion profiles at peak “ii” are included, but all the following conclusions also work for the RCP portion and the peak “i”. As shown in Fig. 4b, c, the topological orders of the vortices are 1 in all cases, proving that the optical vortex emission is robust against disturbance. Therefore, the TBPC system is insensitive to imperfections in experiments, such as fabrication error and thermal expansion effects, making it a stable platform for realizing optical vortices.

Fig. 4: Optical vortex generation as functions of various interlayer separation and twist angles.
figure 4

a The parameter space of interlayer separation and twist angle. The interlayer separation ranges from 200 nm to 350 nm, and the twist angle ranges from 1.0° to 2.5°. The red dots and green circles denote the parameters we choose for the results in b, c, and the dashed lines are guides to the eye. b Normalized amplitude and phase distribution at the red dots in a. c Normalized amplitude and phase distribution at the green circles in a.

Identification of the temporal sequence of the energy transfer in the aforementioned modes provides straightforward evidence of the optical vortex generation mechanism. We demonstrate the time evolution procedure of optical vortex generation where the TBPC system is excited by a pulsed Gaussian beam with a random incident position. At time=0 ps, the Gaussian beam just reaches the TBPC system (Fig. 5a). The Gaussian pulse excites all possible modes, including some low-Q modes in the TBPC system, leading to a chaotic amplitude profile in the monitored plane (Fig. 5c). However, since the low-Q guided resonances dissipate much faster than the BIC mode, the radiation by the BIC mode will gradually appear against the chaotic background as all unwanted low-Q modes vanish after 5 ps. It will form a beautiful optical vortex emission with a zero-intensity point as well as the phase singularity. We notice that the center of the emission is at the AA point rather than the center of the incident beam. It demonstrates that the vortex emission is insensitive to either the incident angle or the illumination position of the incoming light. Consequently, no real-space alignment of the position and incident angle is needed in real applications. This is because the optical vortex emission is excited by the self-contained BIC modes rather than the incident light. More direct demonstration of the robustness of optical vortex generation in TBPC can be found in Supplementary Fig. 6. Thus, in contrast to previous work24, our system has much fewer requirements for the incident light source, making TBPC a promising robust platform for stable vortex generation. Moreover, it is possible to incorporate the TBPC system with micro/nanoelectromechanical systems (MEMS/NEMS) for adjustment of both interlayer separation and twist angle. This could lead to the creation of tunable vortex beams with adjustable quality factors, vortex center positions, divergence angles, etc. Furthermore, with the TBPC system, a tunable OAM laser with adjustable order numbers could be possible.

Fig. 5: Time evolution of the real-space profile.
figure 5

af Depict normalized amplitude and phase distribution of the LCP portion of the emission at different times after excitation. The blue dot in a denotes the position of the AA point, and the white dashed circle denotes the position of the incident Gaussian beam.

In summary, we theoretically show the existence of optical vortex emission in a bilayer photonic crystal system with non-zero twist angles. An interlayer channeling model is formulated based on coupled-mode theory to demonstrate that the optical vortex emission originates from the twist-enabled coupling between the bound state in the continuum (BIC) mode and the guided resonances. Moreover, the optical vortex generation in TBPC is robust against disturbance of geometric parameters, making TBPC a promising platform for well-defined vortex generation, and providing strategies to design stable vortex lasers. Besides the application in vortex generation, the twist-enabled coupling mechanism in TBPC provides a tunable interlayer channel to connect BIC modes to the free space. This will not only benefit the development of BIC study, but also broaden the field of moiré photonics. Due to the intrinsic similarities between photonic systems and condensed matter systems, this work may also guide the exploration of OAM in moiré van der Waals structures.

Methods

Theoretical analysis

Please see the Supplementary Information for the derivations.

Simulations

The electric-field distribution in Fig. 3b, c are calculated using COMSOL Multiphysics from COMSOL Inc. Bloch boundary conditions are applied in the x and y boundaries, while perfect matched layers (PMLs) are applied in the z direction. Without loss of generality, the refractive index of the silicon disks is set to be 3.47. The electric-field amplitude and phase distribution of optical vortices are calculated using commercial finite-difference time-___domain software (Lumerical FDTD solutions from Lumerical Inc.). We place the TBPC system perpendicular to the vertical axis (z axis) and place the AA point at the center of the simulation area. In order to match the experiments, we use a pulsed Gaussian beam as the excitation and apply PMLs in all directions. The electric-field distribution of emission is probed by a frequency-___domain field profile which is 12 μm away from the TBPC.