Electron–matter interaction, which has been studied for a long time ^[1–14], is of great significance in many fields, such as particle detection ^[15,16], radiation emission ^[17–20], and scanning electron microscopes ^[21–24]. It is known that radiation emits when electrons traverse the interface of the different mediums due to the mismatch of the boundary conditions, corresponding to transition radiation (TR) ^[1,2], which distributes in a wide angular range. While the electrons move faster than the light velocity in the medium, the accompanying fields of the moving electrons are thrown away to radiate in the medium in a cone determined by both the electron velocity and the medium permittivity, known as Cherenkov radiation (CR) ^[4,5]. According to these mechanisms, both the TR and CR require a vast accelerating voltage to get relativistic electrons, and the radiation angle is difficult to be controlled, which limits the further development of the amplification of electron–matter interaction.

The progress of artificial materials paves a promising route to explore novel phenomena of the conventional electron–matter interaction to address these problems ^[25–29]. It is reported that the threshold-free inverse CR is found in hyperbolic metamaterials (HMMs) consisting of periodically arranged Au and dielectric layers ^[26], as well as angle-controllable CR by resonance TR in designed photonic crystals ^[27]. Besides, the inverse Doppler and Smith–Purcell effects ^[28,29] have also been presented by employing kinds of metamaterials.

Among these kinds of artificial materials, HMMs exhibit electromagnetic features unattainable in conventional mediums due to the unique hyperbolic iso-frequency dispersion induced by the coupling of propagating or localized surface plasmons (SPs) ^[30–34]. The coupling enables the propagation of evanescent waves and enhances the field, for which HMMs have attracted great interest in nanophotonics and light emissions, for instance, in hyperlenses ^[35] and enhanced spontaneous and thermal emissions ^[36,37]. Due to the excellent electronic properties, graphene, a two-dimensional honeycomb of carbon atoms, has attracted great interest in modern science and technologies such as nonlinear photonics, absorbers for ultrafast pulsed lasers, broadband polarizers, and optical cavities ^[38–41]. And owing to the SPs’ sustainability induced by the dynamic conductivity and tunability resulting from the bias voltage gating, the graphene sheet is also considered a promising candidate for HMMs in the terahertz region ^[42,43]. For a monolayer graphene sheet, it is reported that swift electrons excite both TR and graphene plasmons (GPs) on the surface ^[44–46], and a conventional TR distributed in a wide angular range based on graphene HMMs excited by swift electrons was also presented ^[47].

Here we present a mechanism of direction controllable inverse transition radiation (DCITR) based on the interaction between low-velocity electrons and graphene HMMs by a rigorous method of transfer matrix and GP excitation. The mechanism reveals that the DCITR originates from a spatial dispersion depending on wavevectors induced by the GP coupling and the superposition of the individual coupled bulk GP modes in graphene HMMs. Thus, distinct from the conventional TR, the DCITR propagates along a special inversed radiation angle, which is normal to the Poynting vector and can be controlled by adjusting the chemical potential of the graphene sheets. Further, due to the excitation of the GPs, only hundreds of volts for the accelerating voltage are required and the field of DCITR is greatly enhanced in the HMMs.

When electrons traverse the interface between two different mediums, the waves accompanied with the moving electrons cannot satisfy the boundary conditions, which results in the excitation of additional fields, including the TR propagating outward from the incident spot, and the evanescent waves attenuating rapidly on the interface ^[1,2]. For an interface loaded with a doped graphene sheet, the additional evanescent waves could also excite the GPs further, causing them to propagate and become strongly confined on the interface ^[24].

For the stacked graphene consisting of periodically arranged alternative graphene sheets and dielectric buffers shown in Fig. 1, there are more distinctive characteristics of the GPs in the stack due to the strong coupling between the adjacent graphene sheets. On the one hand, a number of GP modes are excited, which indicates a much more complex frequency dispersion induced by the coupling. On the other hand, the strong coupling makes the GPs become bulk modes, which introduces a spatial dispersion for the stack depending on the wavevector. Thus, the superposition of the bulk modes and the induced spatial dispersion give rise to a DCITR.

To analyze the DCITR, the stack with N layers of graphene sheets is divided into N−1 regions for the different buffers. The free electrons traverse normally to the surface of the stack along the Z direction. Because of the low velocity of the electrons, we only consider the electrons traversing the first graphene sheet with a constant velocity u0 and charge quality q. The conductivity of graphene sheet is σ ^[44,45], and the thickness and permittivity of the dielectric buffer are h and εd, respectively. Then the fields accompanied by the moving electrons are expressed as ^[48] where the subscript i denotes the incident region or the first buffer region, ⊥ denotes the perpendicular direction (the X and Y directions), and k0 is the wavevector in a vacuum.

By employing the additional fields, including the TR and GP fields, the boundary conditions on the first graphene sheet are expressed by the continuity of the fields as where the superscript a denotes the additional fields. In the stack, the fields in the first buffer region are coupled with those in the adjacent regions, which can be obtained by a transfer matrix as where the matrices Mip and Mib are the propagating and boundary transfer matrices in region i, respectively, and where k//,i is the wavevector along the Z direction in region i, k//,i2=εdk02−k⊥2, and η0 is the wave impedance in a vacuum.

Solving Eq. (3), the coupled fields in region 1 are modified by factors ϑa and ϑb, corresponding to the waves propagating along the positive and negative Z direction, respectively. Then the additional fields in the substrate are given by Eq. (2) through an integral in wavevector domain ^[24,47,48] as where χinc=(ω/u0)2+k⊥2−k02εinc and χ1=(ω/u0)2+k⊥2−k02εd.

It is found that Eq. (5) represents a propagating wave when k⊥ is smaller than k0, corresponding to the conventional TR. For the case of k⊥ being larger than k0, the denominator of Eq. (5) is just the frequency dispersion relation of the coupled GPs in the stack, indicating the excitation of the GPs when it reaches zero value at the plasmon poles. The GPs confined in the stack’s energies are much stronger than those of TR. Accordingly, the excited GPs could be obtained by the pole approximation of plasmons and the Taylor expansion ^[24]. We mark the wavevectors at the poles as kGSPsi, and then the DCITR fields are expressed by the superposition of the N modes as ^[24,48] where r is the distance between the electrons spot and observation points, g(x) is the molecule, and f(x)′ is the derivatives of the denominator of Eq. (5), respectively. Then all the fields in the stack can be obtained by the transfer matrix as shown in Eq. (3). For the excitation of GPs, only TM modes are considered in this paper.

Based on the presented mechanism, to study the features of DCITR, the frequency dispersion and excitation probabilities of GPs in the stack are studied first. Assuming a stack consisting of 10 layers of the graphene sheet, the dispersion curves for 10 individual modes are demonstrated in Fig. 2(a). Due to the multi-boundary in the stack, there could be two fundamental GP modes corresponding to the vacuum/buffer and buffer/buffer interfaces shown as the green dashed lines in Fig. 2(a). It is found that the dispersion curves deviate around the fundamental modes, since all 10 individual modes originate from the modification of the fundamental modes induced by the coupling between the adjacent graphene sheets. In the strong coupling region, because of the rising difference between the two wavevectors of the fundamental modes, these dispersion curves are distributed in an increasingly wide region with the increase of frequency. As shown in Fig. 2(b), the plasmon poles for 15 THz are more widely distributed than those for 12 THz. With the further increase of frequency, the coupling is reduced due to the decrease of the penetration depth of each GP mode on the interfaces. Accordingly, it can be seen in Fig. 2(a) that most of the deviated dispersion lines are again gradually approaching the fundamental modes at high frequencies. As shown in Fig. 2(b), the poles for 17 THz present a narrower distribution than for 15 THz; for 22 THz, the poles would almost be concentrated at the fundamental modes.

The normalized excitation probabilities, defined by |f′|−1 according to Eq. (6), are drawn as the gradually changing color in Fig. 2(a), which shows that the GP modes could be excited efficiently at the poles with low group velocities. Thus, it is found that the amplitudes of the modes reach the maximum values around 15 THz and drop sharply with the increase or decrease of the frequency. The detailed excitation probabilities are shown in Fig. 2(b) at the fixed frequencies 12, 15, and 17 THz. It is demonstrated that most of the GP modes are excited efficiently in the strong coupling region, like those at 12, 15, and 17 THz. For 22 THz, only the modes contiguous to the fundamental modes are excited.

To confirm the details of the presented mechanism of the DCITR, the contour maps of the individual bulk GP modes are also studied at a fixed frequency, which is chosen to be 12 THz for a more obvious DCITR at this frequency component. As shown in Fig. 3(a), the contour map of the fourth mode at 12 THz demonstrates that the fields propagate along the R (X and Y) direction, and they distribute along the Z direction as a standing wave, shaping a TM4-like bulk mode in the stack. From a detailed viewpoint, it can be found that the amplitude reaches the maximum on the interfaces of the graphene sheets, indicating that the GPs, which are supposed to be confined on the interfaces, are coupled into a bulk mode by the field interactions between the adjacent graphene sheets. Accordingly, the GP fields can be considered as a propagating wave with a wavevector larger than k0 along the R direction determined by the plasmon poles and an effective one along the Z direction induced by the coupling, which strongly indicates a hyperbolic-like dispersion of the stack. For a higher-order mode, like the tenth mode shown in Fig. 3(b), the fields are mainly distributed on the top several graphene sheets and attenuate rapidly along the Z direction. This is induced by the weak coupling due to the raising plasmon poles, which results in a shorter penetration depth for the high-order mode and makes the stack present a plasmonic-like dispersion. Therefore, due to the GP coupling, the graphene stack presents a spatial dispersion depending on the wavevector, which is the physics of the DCITR.

Next we move on to the synthesized fields of these individual bulk GP modes at 12 THz. As demonstrated in Fig. 4(a), the overlying modes make the fields appear to have unique characteristics that are distinct from the conventional TR or GPs on a monolayer graphene sheet, as shown in Figs. 4(b) and 4(c).

First, parts of the fields are confined around the first several graphene sheets to form the GP-like fields, which are induced by the superposition of the high-order GP modes at large wavevectors with plasmonic-like spatial dispersion. The others are leaked from the upper surface into the graphene stack in a special upwards radiation angle and reflect in the stack to form interference spots corresponding to the propagating-like fields. This is because the hyperbolic-like spatial dispersions of low-order GP modes with small wavevectors enable the GP modes, which are supposed to be confined on the graphene surface, to propagate in the stack. Meanwhile, the composition of these low-order GP modes also makes the propagating fields appear at a unique radiation angle further, which is determined by the plasmon pole distribution of these modes. This is confirmed further by the contour at 22 THz in the stack. It is found in the inset of Fig. 4(a) that the fields are all confined on the surface of the stack, resulting from the plasmonic-like spatial dispersion at a large wavevector due to weak GP coupling.

Second, further calculation shows that the Poynting vector of the GP-like fields lies in the R direction, but that of the propagating-like fields is in the direction of θ (normal to the direction of propagation), which greatly distinguishes it from those of conventional TR or individual GP modes as shown in Figs. 4(b) and 4(c). Furthermore, the dependence of the field intensity on the number of graphene sheets presented in Fig. 4(d) also indicates that the fields probed on the last graphene sheet decrease in the form of an inverse proportional function of the graphene sheet number, which is consistent with the radiation feature. But the fields for a monolayer graphene sheet attenuate exponentially with distance as shown in Fig. 4(e). Accordingly, the DCITR, resulting from the hyperbolic-like spatial dispersion and overlying of the coupled bulk GP modes, is realized.

Based on the mechanism presented above, it can be concluded that the radiation angle and field intensity of DCITR strongly depend on the pole distribution and mode excitation. The dependences of the radiation angle and field intensity on the frequency probed on the lower surface are shown in Fig. 5(a). The stack is chosen to be 30 layers to exclude the GP-like part from fields for the attenuation features of the GP-like and DCITR fields. For a more concentrated pole distribution, the DCITR would be closer to the GPs, which induces a larger radiation angle, until it propagates along the R direction. As the red line shows in Fig. 5(a), the radiation angle decreases with the rising frequency until 15 THz, and then it increases, which is in accordance with the frequency dispersion given in Fig. 2(a). For the field intensity, the maximum is also obtained around 15 THz as the solid blue line shows in Fig. 5(a), which is induced by the largest excitation probability given in Fig. 2(a).

Further, due to the excitation of GPs and propagation of the DCITR, the field intensity is greatly enhanced by about 10 orders of magnitude compared with that of conventional TR or GPs on a monolayer graphene sheet as shown in Fig. 5(a). Therefore, the field intensity of the DCITR can be optimized by adjusting the velocity of these traversing electrons. As shown in Fig. 5(b), the fields reach the maximum for the relatively slow electrons corresponding to an accelerating voltage of about 150 V, which are much smaller than those of conventional TR or CR and are in accordance with those of GPs excited by perpendicular moving electrons ^[44–46].

The pole distributions and excitation probabilities are strongly determined by the coupling of the graphene sheets, which depends on the penetration depths of the GP modes on the individual graphene sheets, and the features of the DCITR can be controlled by adjusting the chemical potentials as well as the thickness and dielectric permittivity of the dielectric buffer. For instance, the rising chemical potential shifts the GP dispersion into a higher frequency region, and then the poles of the fixed frequency are shifted from the weak coupling region into the strong region. Thus, both the excitation probabilities and the range of the pole distribution increase first and then decrease, which indicates thresholds of chemical potential for both the field intensity and radiation angle of the DCITR. As shown in Fig. 5(c), the fields of DCITR at 12 THz reach the maximum at 0.112 eV. Meanwhile, the increasing GP pole distribution range leads to an increasing radiation angle, and a decreasing range induces a decreasing angle. Accordingly, the radiation angle at 12 THz is deduced from 67° to 38.7° by adjusting the chemical potential from 0.05 to 0.1 eV and is raised from 38.7° to 64.6° for the chemical potential ranging from 0.1 to 0.3 eV as shown in Fig. 5(d).

In summary, the mechanism of DCITR in the graphene stack excited by low-velocity electrons is presented. The detailed analyses show that the DCITR originates from the hyperbolic-like spatial dispersion induced by the GP coupling between the adjacent graphene sheets, and the excited bulk GP mode superposition makes the DCITR propagate along an inverse radiation angle normal to the Poynting vector. Therefore, the features of the DCITR strongly depend on the plasmon poles’ distribution and excitation probabilities. Thus, the radiation angle of DCITR can be controlled by adjusting the chemical potential as well as the design of the stack. Further, owing to the excitation of the GP modes, only hundreds of volts for the accelerating voltage are required for the traversing electrons, which is much lower than that of the relativistic electrons used in conventional TR or CR. The field intensities of DCITR are also dramatically enhanced by the excited GP modes.

On the other hand, the mechanism revealed here by rigorous electromagnetic theory also shows how the evanescent waves in a normal medium are transformed into radiation fields in HMMs. Thus, the mechanism can also be applied to HMMs based on the stacked semiconductors in the infrared region as well as stacked noble metals in the visible and ultraviolet region. Even a DCITR in non-periodical stacked plasmonic materials can be analyzed by this mechanism, for instance, the non-periodical graphene chemical potential or dielectric buffer. Accordingly, the presented mechanism of DCITR is of great significance in particle detection, radiation emission, and so on.