Introduction

Terahertz (THz) photonics is an emerging field with significant potential for the 6th generation wireless systems (6G), enabling broadband features and good directivity.^[1-3] The large capacity and high-speed 6G communication need THz devices for high-performance beam steering,^{[4, 5]} modulation,^{[6, 7]} and multiplexing. More recently, a series of studies on THz multiplexing, including polarization division multiplexing (PDM), wavelength division multiplexing (WDM), and mode division multiplexing (MDM), have received considerable interest.^[8-11] Among these, MDM that explores the orthogonality of modes^[12-14] can establish multiple distinguishable data paths using the orthogonal orbital angular momentum (OAM) modes.^[15] The OAM determines the spiral phase front of the vortex beam proportional to exp(ilφ), here φ is the azimuthal angle, l represents the topological charge, and the value of l can take any integer. OAM-MDM has been extensively discussed in optical and quantum communications.^[16-19] In recent years, it has gradually garnered attention as an efficient candidate to increase the spectral efficiency of communication systems in the THz regime.^[20]

The realization of high-quality MDM relies on the convenient and efficient (de)multiplexer for the orthogonal modes.^[21-25] Spatial light modulators and spiral phase plates are two straightforward methods for OAM mode conversion.^{[26, 27]} However, multiple converters and additional beam combiners are required to realize coaxial transmission of orthogonal modes for (de)multiplexing. Metasurfaces composed of sub-wavelength meta-atoms exhibit extraordinary wavefront manipulation applications such as abnormal refraction, beam modulation, perfect absorption, and mode conversion,^[28-32] providing a customizable design method for the high-performance OAM mode generators and multiplexers.^[33-38] Notably, by taking advantage of shared-aperture arrangement^[39] and polarization-dependent meta-atoms,^[40] a single metasurface can implement the 2D multiplexing of polarization and mode to increase the spectral efficiency immensely.^[41-44] Meanwhile, the deflection and focusing of the converted beam can also be integrated into this metasurface, which can significantly benefit the compactness and efficiency of the MDM system.^[45] Recent pioneering reports in the THz regime have shown encouraging advances in metasurfaces for multichannel vortex beam generation^[46-48] and OAM mode multiplexing.^{[49, 50]} Moreover, combining OAM-based MDM and PDM can enhance spectral efficiency and system integration, enabling THz communications to support massive device connections and achieve ultra-high-rate data transmission. Significantly, the working bandwidth of the multi-functional metasurfaces used in multi-dimensional multiplexing and the crosstalk between adjacent channels have not been investigated, which is vital to compact THz wireless communication systems with extremely high data capacity.

Here, we propose all-dielectric metasurfaces for multi-dimensional multiplexing and demultiplexing in the THz regime. The multi-functional phase-modulated metasurfaces function as versatile (de)multiplexers in the MDM and PDM systems. As shown in Figure 1, the Gaussian beams with different incident angles and polarizations are first converted to coaxially transmitted vortex beams with different OAM modes, and each incident beam serves as a data-carrying channel for information transmission. At the receiving terminal, the demultiplexer recovers the coaxial OAM modes to the corresponding Gaussian beams and focuses the demultiplexed beams on different positions. In the experiment, we demonstrated the demultiplexing process of 4-channel and 6-channel coaxial beams using multi-functional metasurfaces. The low crosstalk and high isolation of the demultiplexed channels show good performance of the fabricated metasurface samples at 0.75 THz. Notably, both two types of metasurfaces with linearly polarized (LP) responses and circularly polarized (CP) responses exhibit broadband responses and have the potential to be compatible with integrated WDM. The multi-functional metasurfaces demonstrated here can benefit the development of compact and broadband platforms for multi-dimensional multiplexing in the THz communication system.

Results

2.1 Preparation of Multi-Functional Metasurfaces

Figure 1 illustrates 4-channel 2D multiplexing and demultiplexing based on a pair of multi-functional metasurfaces. Two horizontally polarized (HP) and two vertically polarized (VP) Gaussian beams with various incident angles are converted into high-order OAM beams by the designed metasurface. Here, the state of the four orthogonal channels can be described by $$|H\rangle |l = - 2\rangle $$, $$|H\rangle |l = 2\rangle $$, $$|V\rangle |l = - 2\rangle $$ and $$|V\rangle |l = 2\rangle $$. The coaxial OAM modes with orthogonal polarizations can serve as independent channels to carry digital signals for high-capacity and high-speed data transmission, the intensity distributions of these four channels are shown in Figure 1a. After coaxial transmission over a specific distance, the vortex beams with different topological charges can be recovered to the Gaussian beams ($$|H\rangle |l = 0\rangle $$ and $$|V\rangle |l = 0\rangle $$) and deflected to different directions using the demultiplexing metasurface. It is worth mentioning that the demultiplexed metasurface can focus the output beam and thus optimize the beam divergence angle to obtain the high-quality Gaussian beam at the receiving terminal. As shown in Figure 1b, the demultiplexed beam corresponding to Channel 4 (Ch4) can be focused on the aim terminal. As the interleaved method is adopted in our metasurface design, this demultiplexed beam is the superposition of correct-demultiplexed $$|V\rangle |l = 0\rangle $$ and incorrect-demultiplexed $$|V\rangle |l = - 4\rangle $$. The detailed analysis is shown in Section S.1 (Supporting Information). However, we can efficiently receive the energy of the fundamental Gaussian mode by adjusting the diameter of the receiving device because the spot diameter of the OAM beam with the topological charge of l = −4 is much larger than that of the fundamental-mode Gaussian beam and the intensity distribution of the OAM beam is a hollow ring, as shown in the dotted ring in the figure. Similarly, the electric fields of the other three channels are also a superposition of beams with l = 0 and |l| = 4. We can identify the energy of the correct-demultiplexed fundamental mode without additional mode filters. The designed two metasurfaces can replace many bulky devices, such as q-plates, beam deflectors, and Fourier transform lenses, thus providing an efficient, versatile, and compact platform for multi-dimensional (de)multiplexing.

The same design method can be used to design multiplexers and demultiplexers according to the reversibility of the light path. Two different phase distributions corresponding to the HP and VP components are required to implement polarization and mode multiplexing/demultiplexing simultaneously. As shown in Figure 2a,b, for the HP component, the transmission function can be expressed as:

$$$$\begin{equation}{\Phi }_{{\mathrm{HP}}} = A{\mathrm{exp}}\left[ {{\mathrm{i(}}l\theta + m\frac{{2{{\pi}}x}}{T} + k{\mathrm{(}}{\mathrm f} - \sqrt {{{\mathrm f} }^2 + {x}^2 + {y}^2} {\mathrm{))}}} \right]\end{equation}$$$$(1)

where A is the transmission amplitude of the metasurface and is assumed to be 1 when only phase modulation is considered, l is the topological charge, θ = arctan(y/x) is the azimuthal coordinate, T is the period of phase gradient, m is the diffraction order and m = ±1 for the 4-channel MDM. Here, we adopt the checkboard-like interleaving method to spatially sample the phase distribution (abbreviated to PD) for OAM mode multiplexing with the topological charge of l = +2 (PD1+PD3 in Figure 2a) and l = −2 (PD2+PD4 in Figure 2). The last term in the phase is for focusing, corresponding to the PD5 in the figure. Here, k = 2π/λ is the wavevector, λ represents the wavelength, and f is the focal length. Similarly, for the VP component, the transmission function is:

$$$$\begin{equation}{\Phi }_{{\mathrm{VP}}} = A{\mathrm{exp}}\left[ {{\mathrm{i(}}l\theta + m\frac{{2{{\pi}}y}}{T} + k{\mathrm{(}}f - \sqrt {{{\mathrm f} }^2 + {x}^2 + {y}^2} {\mathrm{))}}} \right]\end{equation}$$$$(2)

To obtain the different transmission responses for HP and VP components, the phase modulation of the metasurface at each point should satisfy the Φ_HP and Φ_VP in the two formulas above. As shown in Figure 2e, we used polarization-dependent silicon meta-atoms to compose the metasurface. The high-resistive silicon wafer exhibits a non-dispersive refractive index (∼3.45) and low absorption loss in the THz band. Sixty-four meta-atoms with a period of 160 µm (0.4 λ at 0.75 THz) were selected to achieve independent phase modulation for HP and VP with a phase interval of 45°. Details of the chosen pillars can be seen in Table S1 (Supporting Information). The metasurface contains three layers: a silicon pillar with a thickness of h = 260 µm, a 688-µm-thick silicon substrate, and a periodic anti-reflection layer with h_a = 52 µm thickness. As shown in the right image of Figure 2e, the anti-reflection layer^[51] consists of the periodic quadrate pillars with side length l_x = l_y = 48 µm and the period p = 63 µm. The phase accumulation of the selected meta-atoms under HP and VP incidence simulated by CST Microwave Studio is shown in Figure 2g. Figure 2f shows the simulated transmission amplitude of the selected meta-atoms. The average t_x(t_y) of the chosen meta-atoms is 0.815, and the transmissivity of the ordinary silicon wafer is ≈0.7, which indicates that the anti-reflection layer can effectively improve transmission efficiency. We fabricated the proposed metasurface by double-sided registration lithography and deep reactive ion etching. The details can be found in the “Method” section. For the fabricated MS1 sample, the optical image of the upper layer is shown in the lower right corner of Figure 2c.

2.2 Two-Dimension Demultiplexer with LP Responses

The performance of the multi-functional metasurface was first investigated by the Rayleigh–Sommerfeld (RS) diffraction integration^{[52, 53]} in theory. The RS diffraction formula we used in the calculation is described by:

$$$$\begin{equation}E{\mathrm{(}}{P}_1{\mathrm{) = }}\sum {\sum {\frac{d}{{{\mathrm{i}}\lambda }}} } {E}_0{\mathrm{(}}{P}_0{\mathrm{)}}\frac{{{e}^{ - {\mathrm{i}}k{r}_{01}}}}{{r_{01}^2}}dxdy\end{equation}$$$$(3)

where E(P₁) and E₀(P₀) represent the electric field at point P₁ on the objective plane and point P₀ on the metasurface, respectively. d is the distance between the metasurface and the objective plane; r₀₁ is the distance between the points P₀ and P₁. For OAM mode multiplexing, the complex amplitude E₀ can be expressed as E₀ = E_in·exp(iΦ_m). Here, E_in = A₀·exp(iΦ_g) is defined as the initial incident Gaussian beam, A₀ = exp[−(x² + y²)/w²] is the amplitude, and Φ_g is the phase gradient due to oblique incidence, Φ_m represents the phase distribution for metasurface serving as OAM mode multiplexer. For OAM mode demultiplexing, the E₀ in Equation 3 can be defined by: E₀ = E_OAM·exp(iΦ_dem), E_OAM is the complex amplitude of coaxial OAM beams, and Φ_dem represents the phase distribution for metasurface serving as mode demultiplexer. The period of phase gradient in Φ_dem is set to T = 2.02 mm for separating the different channels with a refractive angle of 11.42°. The calculated intensity distributions of demultiplexed beams at the focal plane of f = 20 mm are shown in Figure 3a. The focus points corresponding to the correct-demultiplexed fundamental mode can be observed at four different positions, indicating that the four coaxial channels have been demultiplexed and separated successfully.

The characteristics of the fabricated metasurfaces were experimentally investigated by the THz near-field time-domain spectroscopy (NF-TDS) system. Details of the experimental setup are shown in Figure S2 (Supporting Information). Since it was difficult to tilt the transmit antenna in our NF-TDS system, we fabricated the mode converter to replace the multiplexer used for oblique incidence, which can convert the normal-incident THz beam into OAM beams with the topological charge of l = ±2, the details of the design method are shown in Section S.4 (Supporting Information). The intensity and phase distributions of converted superimposed vortex beams with HP and VP can be seen in Figure 2d. We use the fabricated multi-functional metasurfaces as the 2D demultiplexer. Figure 3 shows the measured and calculated results of demultiplexed channels at four receive positions, verifying the 2D modulation of the metasurface. A 4-bit binary code is used to simulate incident beams that carry different information. Out of the demultiplexer, the spatially co-propagated OAM beams carried different information that could be recovered to Gaussian beams and be focused to four positions, as shown in Figure 3a. The result scanned by the NF-TDS system in Figure 3a is the x-y plane perpendicular to the z-axis. But in practice, the detector should be placed in the plane perpendicular to the wavevector of the demultiplexed beam. Here, we process the measured data by RS integration to obtain the intensity and phase distributions in the plane perpendicular to the wavevector of the demultiplexed channels. The results for each channel are shown in Figure 3b, and the data processing details are in Section S.5 (Supporting Information). When the incident beams carry 1100, the clear focus spots and uniform phase distributions can be observed in two symmetric positions on the x-axis. If the incident beams carry 0011, the two focus spots can be seen on the y-axis.

Based on the measurement results in Figure 3b, we can quantitatively analyze the crosstalk and isolation of four channels, where the aperture of the receiver is set to 1.36 mm (white dashed circles in the figure), the same as the diameter (calculated by 1/e²) of the focused spot. The intensity matrix of the four demultiplexed channels is shown in Figure 3c, where the values are calculated by N_dB = 10·lg (I_r/I_i), I_r is the intensity of the correct-demultiplexed Gaussian beam within the receiving aperture, and I_i represents the intensity of the corresponding incident beam. The diagonal of the matrix is the value of the correct-demultiplexed beam, and the others are crosstalk. The maximum crosstalk among all channels is −16.88 dB. The isolation of each channel is shown in Figure 3d, which is calculated by N_is = N_r – N_c, N_r is the intensity of the correct-demultiplexed beam in the diagonal of the intensity matrix. N_c = 10·lg (I_c/I_i), I_c represents the sum of the intensities of all crosstalk aroused by incorrect-demultiplexed beam within the receive aperture, and the measured isolation of the demultiplexed channels at 0.75 THz can be greater than 15.08 dB. We also analyzed the performance of the designed devices over the broadband range of 0.65–0.85 THz based on the measured results. It is worth mentioning that we used broadband mode convertors in the experiment, and the normalized intensity of all incident OAM modes in the measured frequency range is greater than 83%. The details of incident vortex beams can be found in Figure S5 (Supporting Information). The isolation of the four channels in the whole interested range is shown in Figure 3e, while the crosstalk is analyzed in Section S.7 (Supporting Information). Due to the dispersion, the wavelength deviation will decrease in the received intensity. However, the proposed metasurface still exhibits wideband performance. The 3 dB bandwidth of all correct-demultiplexed channels is 0.13 THz, and the isolation of each channel can be greater than 10 dB over a range of 0.15 THz. Considering that isolation and insertion loss are two critical factors in multiplexing systems, we defined the adequate working bandwidth (AWB) of the devices as the frequency range in which the received energy attenuates less than 3 dB with respect to the center frequency and the isolation of each channel is greater than 10 dB. Thus, the AWB of the designed MS1 is 0.13 THz. The results indicate that the proposed 2D demultiplexer is promising for THz communication compatible with WDM. As two cascaded metasurfaces are adopted in our MDM strategy, the misalignment issues may influence the performance of the devices. We provide a theoretical analysis and discussion on this issue in Section S.8 (Supporting Information). The calculated results show that as long as the displacement is less than ±2 mm (31% of the side length of the metasurface), the decrease of received intensity and isolation can be less than 3 dB and 8 dB respectively, indicating that the multiplexing system can still perform well under minor displacement. Meanwhile, the calculated results shown that the devices can be robust to the minor rotation if the diffraction orders are designed to be symmetric.

2.3 Two-Dimension Demultiplexer with CP Responses

The HP and VP are two essential and commonly used orthogonal polarization states for the PDM technique used in optical communications. Besides HP and VP, the left-circular polarized (LCP) and right-circular polarized (RCP) components also can be used as two orthogonal polarization channels. Using the circularly polarized components can decrease the alignment difficulty of the polarization coordinates of the transmitter and the receiver, and the polarization degree of a circularly polarized beam can maintain better performance in the free space transmission compared with that of a linearly polarized beam.^[54] Therefore, we designed MS2 with CP responses. Similar to MS1, the designed MS2 needs to provide two independent phase profiles for the LCP and RCP components, and the transmission functions of Φ_LCP and Φ_RCP are the same as Φ_HP and Φ_VP we mentioned above, respectively. Here, we obtain the target phase profiles by combining geometric phase and propagation phase modulation^[55] in the design method. The phase modulation of the metasurface at each point should satisfy Φ_LCP = exp[i(φ_l)] = exp[i(φ_p + φ_g)] and Φ_RCP = exp[i(φ_r)] = exp[i(φ_p − φ_g)], here φ_p and φ_g represent propagation phase and geometric phase, respectively. To design a metasurface with arbitrary phase distribution φ_l and φ_r, the propagation phase at each point should satisfy φ_p(x, y) = (φ_l + φ_r)/2, and the orientation angle ?? is defined by ??(x, y) = φ_g/2 = (φ_l − φ_r)/4. As shown in Figure 4c, we used the anisotropic silicon pillars to compose the metasurface. The image of the fabricated MS2 sample is shown in the bottom right corner of Figure 4a. We selected eight pillars with a propagation phase interval of 45°, four of them are displayed in Figure 4d. The other four structures can be obtained by rotating these four pillars with an angle of 90°. For the chosen meta-atoms, the phase difference between the geometric phase φ_x and φ_y is very close to π, thereby implementing efficient conversion of the two orthogonally circular polarizations, and the average conversion efficiency η_lr(η_rl) is 0.8244. Figure 4e shows the geometric phase of eight meta-atoms, which illustrates that the phase delay φ_l can cover the range of 0–2π when rotating the pillars from 0 to π. The details of the selected meta-atoms are shown in Table S2 (Supporting Information).

Figure 4b illustrates the measured intensity and phase distributions of the coaxial beams with superimposed OAM modes of l = ±2 for LCP and RCP, respectively, which was obtained by another mode converter. As illustrated in Figure 4a, MS2 is used as a 2D demultiplexer with CP responses, and the intensity distributions of the demultiplexed channels are shown in Figure 5a. The right of Figure 5a shows the intensity and phase distribution at the focused positions of the four channels, demonstrating that the OAM modes with different polarizations are successfully demultiplexed into spatially separated fundamental Gaussian modes by the MS2. Figure 5b shows the measured intensity distribution in demultiplexed focal plane illuminated by different OAM beams, respectively, and the intensity distributions at four focal positions in the plane perpendicular to the wavevector of correct-demultiplexed beams are shown on the right of Figure 5b. Figure 5c represents the intensity matrix of the four demultiplexed channels. It can be seen that the crosstalk between the orthogonal polarization channels is greater than that between different OAM mode channels with the same polarization. The former of crosstalk can be attributed to fabrication errors because the small deviation of element size will cause the reduction of polarization conversion efficiency, and the incomplete conversion of energy can increase the crosstalk. However, the fabricated MS2 still shows good performance. The maximum crosstalk among all channels is −14.48 dB, and the measured isolation of the channels is above 13.06 dB (as shown in Figure 5d). Figure 5e is obtained by superposing the intensity distribution in the focal plane of demultiplexed beams under four incident vortex beams carrying the codes 1000, 0100, 0010, and 0001, respectively. We also evaluate the performance of the designed MS2 over the broadband range of 0.65–0.85 THz, and details of the broadband incident vortex beams are shown in Figure S6 (Supporting Information). As shown in Figure 5f, the 3 dB bandwidth of the demultiplexed channels is 0.13 THz, and the measured isolation of each channel can be greater than 10 dB over a range of 0.12 THz, which means that the AWB of MS2 can exceed 0.1 THz. These indicate that the designed metasurface can perform 2D demultiplexing in a broadband regime.

Discussion

The design method we used here can be extended to realize multiplexing and demultiplexing for more OAM channels or other orthogonal mode channels, such as Herimite-Gaussian modes and vector modes. As a proof of concept, we designed MS3 for 6-channel multi-dimensional demultiplexing, and the states of the six coaxial channels can be described by $$|{\mathrm{LCP}}\rangle {\mathrm{|}}l = \pm {\mathrm{2,}}l = - 4\rangle $$ and $$|{\mathrm{RCP}}\rangle l = \pm 2,l = - 4\rangle $$. Figure 6a shows the intensity distribution of the demultiplexed channels illuminated by different vortex beams. The superposition results of all the LCP and RCP components in Figure 6c indicate that MS3 could also show good performance under the incidence with multiple modes. As interleaving is a kind of shared-aperture design method, the increase of demultiplexed channels inevitably leads to decreased energy received by each channel and the introduction of extra crosstalk.^[56] Designing the anti-reflection layer is an effective way to improve device efficiency, as demonstrated here before. In addition, using the optimization algorithm to optimize the phase distribution of the meta-devices can increase the power received in each channel.^[57] In the experiment, we adjusted the aperture of the receiving terminal to decrease crosstalk, and the diameter of the receiving aperture was set to be the same as the full width at half maximum of correct-demultiplexed Gaussian mode in the top row of Figure 6d. The bottom row of Figure 6d shows the superposed intensity distribution of all the crosstalk aroused by the five incorrect-demultiplexed beams. A significant intensity difference between the Gaussian mode and crosstalk within the receiving aperture can be observed. The measured isolation of Ch1-Ch6 is shown in Figure 6b. Except for Ch4, whose isolation is 9.77 dB, the isolation of other channels is greater than 10 dB. According to our definition of device bandwidth, the MS3 can only be used as a narrow-band device. The results in Figure 6d show that the isolation of Ch2 and Ch4 is greater than that of other channels as the minimum interval of the topological charges between these two channels and adjacent channels is |4| comparing |2| of other channels, which indicates that the crosstalk can be reduced by increasing the intervals of the topological charges between adjacent channels. In a multiplexing system, isolation and insertion loss are two critical factors. In theory, it is difficult to accurately predict the deterioration of these two factors as the number of channels increases due to the complexity of the 2D multiplexing system. Our measured results demonstrated that, under the given experimental conditions, a 6-channel metasurface-based demultiplexing setup can achieve an isolation of nearly 10 dB and an insertion loss of ≈–15 dB. This indicates that the few-channel (6 or less) MDM+PDM strategy holds high promise to be integrated into broadband THz modulation systems. Moreover, we analyze the sensitivity of the MDM system to the angle error of the transmitter and receiver in Section S.10 (Supporting Information). The results indicate that the error of the angle should be less than ±1°. Reconfigurable and intelligent metasurfaces^[58-61] can be utilized to design multiplexers that response rapidly to wide-range incident angles and realize self-adaptive beamforming, which hold promise to increase the flexibility of THz demultiplexer.

Conclusion

In conclusion, we proposed a multi-functional metasurface to realize the 2D MDM-PDM in the THz domain. 4-channel LP demultiplexing and 6-channel CP demultiplexing were experimentally demonstrated by the NF-TDS system. For the 4-channel LP demultiplexing, a low crosstalk of −16.88 dB between different channels is achieved. Moreover, the 3 dB bandwidth of the metasurfaces is 0.13 THz, and the measured isolation of each channel can be greater than 10 dB over a range of 0.1 THz, which indicates that they have great potential to be compatible with WDM. This work provides an efficient and compact platform for integrated THz MDM and PDM, and the results will help to pave the way for multi-dimensional multiplexing in high-capacity THz communication systems.

Experimental Section