The beam splitter is a fundamental component in optics and electromagnetic science, which can change the propagation path of an incident beam and distribute it into multiple output channels. Initially, beam splitters were realized using the birefringence effect of prisms, and later, some more compact implementations based on planar structures appeared [1,2]. So far, beam splitters have played important roles in various applications such as interferometers [3], quantum optics [4], 3D laser printing [5], optical signal processing [6], and optical wireless communications [7]. Beam splitters are also expected to contribute to terahertz (THz) technology. Given the abundant spectral resources, the THz band is expected to be fully explored in future 6G mobile communication systems to provide higher data transmission rates [8], which can accommodate a large number of users within a single communication link. To this end, it is critical to split THz signals and focus them onto multiple receivers by developing multi-port THz beam splitters. At the same time, the number and positions of output channels are preferred to be reconfigurable to respond to different working statuses flexibly.

In traditional THz systems, the common devices that manipulate THz waves are dielectric lenses and metal reflectors, which are bulky and thus hinder the compact integration of the system. Metasurface-based devices (meta-devices) have attracted wide attention recently due to their ultra-thin and compact properties [9–11]. Metasurfaces consist of sub-wavelength artificial structures (meta-atoms) that can manipulate the wavefront of electromagnetic waves with sub-wavelength precision. Therefore, using different materials to design meta-atoms of different geometric sizes offers additional degrees of freedom for highly flexible manipulation of electromagnetic waves. In this way, the functions of existing optical components can be extended to develop new devices, such as deflectors [12], metalenses [13], holographic displayers [14,15], three-dimensional (3D) imagers [16], orbital angular momentum generators [17], meta-sensors [18], vector beam generators [19], and surface wave generation devices [20,21]. The metasurfaces can also facilitate advances in 6G THz communication application fields because of their abilities in mode, polarization, and frequency multiplexing [22,23].

One limitation of common metasurfaces is that once fabricated, their functions are determined and difficult to adjust flexibly. The researchers put forward plenty of schemes to realize dynamic reconfigurability for metasurfaces. In general, this function can be achieved by introducing active components such as liquid crystals [24], phase change materials [25], graphene [26], micro-electro-mechanical systems (MEMS) [27], and complementary metal oxide semiconductor (CMOS) transistors [28], but these components are costly and complex to develop in the THz band. Fortunately, there is an alternative approach to achieving reconfigurable functions without introducing these expensive and complex active components, namely, stacking multiple metasurfaces in space and adjusting their relative spatial positions. This reconfiguration method attracts vast attention due to its simple and low-loss features [29,30].

Most existing stacked metasurfaces are designed using traditional methods guided by physical intuition [31–33]. However, using this design method, it is challenging to multiplex various functions into them. Recently, inverse design methods have provided more flexibility in designing metasurfaces [34,35]. These methods introduce computational optimization to design metasurfaces according to the functions they are supposed to achieve. They can encode multiple functions into a single meta-device, thus allowing the realization of many interesting optical phenomena, such as diffractive neural networks [36–39], optical encryption/decryption [40], and dynamically switchable diffraction patterns [41]. Therefore, using inverse-designed stacked metasurfaces is a promising way to achieve a reconfigurable THz beam splitter.

In this paper, we propose a reconfigurable terahertz beam splitter consisting of two inverse-designed metasurfaces cascaded together, as schematically shown in Fig. 1. Adjusting the relative rotation angle between the two metasurfaces allows the beam splitting direction and focusing area of THz waves to be dynamically controlled and switched. The focal spot arrays of THz waves can be projected to any position on the focal plane (i.e., any multi-port projection), and the number and positions of these focal spot arrays can be flexibly switched by rotating the meta-device. The phase distribution is realized by tuning the height of meta-atoms, and this meta-device was fabricated using the 3D printing technique. We introduce the gradient descent method to inversely design the phase distribution, ensuring that THz waves can be effectively split and focused to form focal spot arrays. When the rotation angles are set to be 0°, 90°, 180°, and 270°, the meta-device can produce four different focal spot arrays. In this way, the four working modes can be flexibly switched by rotation. Our measurement results agree well with simulations, validating a good performance of this meta-device. Our scheme provides an alternative for designing the signal splitters for future 6G THz communications.

We note that various far-field diffraction patterns represent different beam-splitting directions. To obtain these far-field diffraction patterns, we use an inverse design method based on gradient-descent optimization (shown in Fig. 2). We first consider the supercells with phase distributions of only 2×15×15 pixels, which are the optimization variables. The relationship between the phase distributions ϕ and the far-field diffraction intensities I at different rotation angles θ can be described as Iθ(u,v)=|FT(exp(jϕ1(x′,y′))×rotationθ(exp(jϕ2(x′,y′))))|2,(1)where FT(·) denotes the Fourier transform, rotationθ(·) represents the matrix rotation operation at rotation angles θ, and (x′,y′) denotes the Cartesian coordinate of the supercells. This meta-device can produce distinct far-field diffraction patterns for different θ, which are supposed to form our desired far-field intensities. Such a relationship between the phase distributions and their far-field diffraction patterns is realized using the optimization process. The objective of the optimization process is to minimize the difference between diffraction intensities I and target intensities I^ as much as possible. The loss function is defined by calculating the mean square error (MSE) between I and I^ to evaluate their difference: L=1N∑i=1N(Ii(u,v)−I^i(u,v))2,(2)where N=4×15×15=900 is the total pixel number of the four intensity pictures and (u,v) denotes the Cartesian coordinate of the spatial frequency plane. Then the gradient of the loss ∇ϕL is automatically calculated using the machine learning library of PyTorch (see more details in Appendix A). Forward and backward propagations are the two basic ideas in PyTorch. Here, we define the forward propagation as the far-field calculation described by Eq. (1), while the backward propagation is the gradient of the loss function ∇ϕL. Once the forward propagation is defined, the backward propagation can be automatically calculated to reduce the loss, and thus, the variables can be optimized. In this way, the entire optimization process can be expressed as finding the variables that meet ϕ^=argminϕL(ϕ). After some iterations until the loss converges, we can obtain the final updated phase distributions.

Compared with traditional phase retrieval methods like the Gerchberg-Saxton (GS) algorithm, this scheme is based on modern machine learning libraries (like PyTorch). Specifically, as we have defined the forward propagation in Eq. (1), the far-field intensities can be linked to the phase distributions of the meta-device, which are the variables that need to be optimized. Thus, the forward propagation can be easily realized. PyTorch is a uniform and easy-to-use framework to perform optimization tasks and is adopted here. It can automatically calculate the backward propagation (i.e., the gradient of the loss). Since the loss function has been well-defined [see Eq. (2)], PyTorch can find the maximum gradient descent direction in the variable space, and then the phase variables are updated along this descent direction. In this process, the value of loss decreases and finally converges, and we can get the final optimized phase variables that make the loss minimum, i.e., achieving the desired far-field intensities. Moreover, various optimizers [such as stochastic gradient descent (SGD), root mean square propagation (RMSProp), adaptive moment estimation (Adam), and limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS)] in the machine-learning-based framework have been proposed by former researchers to realize advanced automatic differentiation methods, thus making optimization tasks more accurate and faster.

The optimized phase profiles are chosen to form the supercells of the two metasurfaces, each composed of N0×N0 (N0=15) meta-atoms with a pixel size of P=2mm, as shown in Fig. 3(a). The corresponding far-field intensities at four different rotation angles (0°, 90°, 180°, and 270°) are shown in Fig. 3(b). Theoretically, there are 15×15=225 locations for beam-splitting channels in each rotation angle. According to the Nyquist sampling theorem, the maximum spatial frequency (fx,y) is half of the spatial sampling frequency (1/P), i.e., |fx,y|≤12P. Next, we need to consider the interference effects caused by periodically arranged supercells. Mathematically, suppose we periodically pad the supercells m times. In that case, the corresponding intensities after the Fourier transform will be interpolated by n zeros, as shown in Figs. 3(c) and 3(d) (with m=2, Nrep=m+1=3, and N1=Nrep×N0=45). Thus, after the periodic arrangement, the splitting direction (spatial frequency) of each beam remains unchanged, but the resolution improves, indicating more accurate focal spot arrays. Before the periodic arrangement, the resolution is Δf0=1N0P [see Figs. 3(a) and 3(b)] in the spatial frequency plane; by contrast, after that, the resolution increases to Δf1=1N1P [see Figs. 3(c) and 3(d)], with a resolution improvement of Δf0Δf1=N1N0=Nrep=3. Such a resolution improvement can also be physically explained by the increased numerical aperture (NA). After periodic arrangement, the NA improves three times, and the resolution also increases accordingly.

In this way, we obtain the beam-splitting phase profiles of the two metasurfaces. Then, we introduce the beam-focusing phase profile (ϕf) to make the beams with different splitting directions focus on the focal plane [Fig. 3(e)]: ϕf=−π(x2+y2)λz,(3)where (x,y) is the position in the Cartesian coordinate system, λ=1mm is the working wavelength, and z=150mm is the focal length between the meta-device and the focal plane. This phase profile works as a Fourier lens, relating its object plane (meta-device) to the image plane (focal plane) in the spatial Fourier domain. As shown in Fig. 3(e), we set the pixel size to p=0.5mm (half wavelength) to provide a higher spatial resolution. The corresponding pixel number is N2=4×N1=180, ensuring that the dimensions of the two metasurfaces remain unchanged (L=90mm). These phase profiles are then employed to map the two metasurfaces [Figs. 3(f) and 3(g)], which are fabricated using the 3D printing technique (see detailed sample fabrication process in Appendix B). The 3D printing technique has the advantages of high fabrication speed and low manufacturing cost. We print the meta-atoms with different heights h using photo-curable high-temperature resin (n˜=n+jκ=1.6311+j0.0245) to match the phase profiles of the meta-device [Fig. 3(h)]. The meta-atom has the shape of a square pillar with varying height, whose cross section is 0.5mm×0.5mm. The period size p is also 0.5 mm (half wavelength), which makes it capable of manipulating THz waves at a subwavelength resolution. The meta-atom’s phase modulation is realized by changing the equivalent THz wave propagation length [i.e., (n−1)h] along the z-direction, and this relationship between the corresponding phase and height can be described as ϕh=2πλ(n−1)h. Thus, the phase profiles are matched by tuning h. We also note that no polarization-related factors are involved. So, this meta-device has the advantages of polarization insensitivity, which can modulate the THz waves with any specific polarization. In this way, the 3D-printed meta-device is encoded with the pre-designed phase profiles.

A THz planar scanning measurement setup is assembled to characterize the performance of the meta-device, as shown in Fig. 4(a) (see detailed experimental setup in Appendix C). A vector network analyzer (VNA) generates radio frequency (RF) signals, which are then upconverted to the working frequency of 0.3 THz by a frequency extender. The y-polarized THz waves are radiated via a horn antenna (Tx) and then collimated by a dielectric lens. The collimated THz waves go through the sample of the meta-device, generating the required field distributions (in y-polarization) on the detection plane, which are scanned by a probe (Rx). The Rx is mounted on a three-axis translation stage to facilitate scanning. Finally, the THz signals are downconverted back to RF by another frequency extender, then analyzed and recorded by the VNA.

We next demonstrate the performance of the meta-device working at four various rotation angles (θ=0°, 90°, 180°, and 270°). Figure 4(b) presents the field distributions on the focal plane (focal length z=150mm), including both numerically simulated and experimentally measured results (see detailed numerical simulation in Appendix D). We note that the distance between neighboring focal spots (Δx,y) is 5 mm, calculated by Δx,y=z(N0P)−1λ−1. It is because the focusing plane is also called the Fourier plane, with the relationship of Δx,yz=Δkx,yk0=(N0P)−1λ−1. This distance is verified by simulated and measured results [shown in Fig. 4(b)]. The four rotation angles (four working modes) correspond to different scenarios. As for working mode 1 (θ=0°), the input THz signal is split and focused to form nine uniformly distributed output THz signals in the center. By contrast, in working mode 2 (θ=90°), two groups of focal spot arrays are positioned at the lower-left and upper-right corners, with varying distances from the center, each composed of four closely spaced focal spots. As for working mode 3 (θ=180°), the eight focal spots are located around the center, with the shape of the focal spot array differing from that of working mode 1. Interestingly, when switched to working mode 4 (θ=270°), the meta-device generates a single intensive bright spot at the center of the focal plane, which well suits the scenario of a single receiver. Therefore, the focal spots’ number and locations can be flexibly switched by rotating this meta-device. The measured results agree well with the simulations, validating the good performance of the proposed meta-device.

In the optimization and simulation, we neglect the distance between the two metasurfaces, i.e., Δz=0mm. However, in practice, this distance will inevitably be introduced. We evaluate the performance of the meta-device affected by Δz, as shown in Fig. 5. As Δz increases, the performance deteriorates. Because of the diffraction effect, the wavefront cannot remain unchanged after propagating through this distance, and the diffraction effect becomes more obvious as the distance grows. Thus, the wavefront’s deformation degree also increases. Besides, we notice that the deformation degree depends on the phase gradient of the original wavefront before diffraction. A smaller phase gradient may cause a lower deformation degree. At the same time, to minimize undesired effects (such as the near-field coupling), Δz should not be too small. Therefore, Δz=2−3mm (2λ−3λ) may be more suitable.

The dimensions of the meta-device will also influence its performance (Fig. 6). We use the repetitive number Nrep to represent the dimension (L=N0×Nrep×p) of the meta-device. As discussed in Fig. 3, a larger Nrep corresponds to a higher resolution of the focal spot arrays (smaller spot size). It is important to note that the resolution is subject to the Rayleigh diffraction limit. The meta-device should not be too large, as it would take up too much space. In practical use, the incident THz beam may be a Gaussian beam, and the approximate plane wave incidence cannot be effectively ensured across the entire range of the meta-device. It is the main reason for the discrepancies between the simulated and measured results in Fig. 4. The limited aperture of the dielectric lens prevents the approximated plane wave from fully covering the meta-device, thereby compromising its performance.

We also notice there are several sidelobes around the main spots. Such a phenomenon is caused by the limited aperture and uniform amplitude (i.e., all-one amplitude) modulation of the meta-device. The meta-devices with finite apertures always produce sidelobes around the main spots. Besides, according to the array antenna theory, the units with uniform amplitude can also lead to sidelobes. The approaches to suppress them are by increasing the spatial sizes and introducing amplitude-modulated meta-atoms, and the amplitude distribution may also be optimized with well-defined objective functions. As shown in Fig. 6, the meta-device with a larger spatial size (Nrep=4) can generate weaker sidelobes than that with a smaller spatial size (Nrep=3). However, in the experiment, we cannot generate such a wide THz beam to cover the entire meta-device with Nrep=4. Thus, we still use the smaller one (Nrep=3) to conduct the experiment. In the present design, we do not consider realizing complex-amplitude modulation. Other well-designed metasurfaces may help realize this [42], thus achieving sidelobe suppression.

Compared with other existing reconfigurable terahertz metasurfaces, this mechanically reconfigurable meta-device has higher efficiency. Because it is made of passive dielectric materials, the losses are lower than metal materials. For example, no ohmic loss is introduced. To illustrate this, we compare the focusing efficiency of previously reported reconfigurable terahertz metasurfaces, as shown in Table 2. Here, we consider the efficiency of working mode 4 (single focusing spot) to be the focusing efficiency of this meta-device (see Table 1).Table 2.

Efficiency Comparison of Our Proposed Meta-Device with Previously Reported Works

The simulated field intensities on the focal plane at three different frequencies (0.27 THz, 0.3 THz, and 0.33 THz) are calculated, as shown in Fig. 8. The location of focal spots changes as frequency varies. To ensure the meta-device’s normal operation, each focal spot’s spatial deviation at different frequencies should be less than 2.5 mm (half the distance between neighboring focal spots at 0.3 THz). To obey this, we estimate the normal working range is from 0.27 THz to 0.33 THz. Thus, the bandwidth is 0.33THz−0.27THz=0.06THz, and the relative bandwidth is (0.06THz/0.3THz)×100%=20%. Although the meta-device is designed at a single working frequency of 0.3 THz, it still shows a relatively broadband behavior (0.27–0.33 THz). This feature makes it promising for future 6G mobile communication system applications.

To conclude, we propose a scheme to generate reconfigurable THz focal spot arrays with a rotatable doublet meta-device. By adjusting the relative rotation angles between the two metasurfaces to 0°, 90°, 180°, and 270°, the number and position of focal spot arrays on the focal plane can be flexibly switched among four distinct working modes. The gradient descent method is introduced to inversely design the phase distributions, ensuring the THz waves are precisely split and focused to form the desired focal spot arrays. To physically achieve this, we fabricate the meta-device using the 3D printing technique, and the phase distribution is realized by tuning the heights of meta-atoms. The measurement results and simulations are in good agreement, validating the excellent performance of this doublet meta-device. This approach paves a way to design inexpensive components for future 6G THz communication systems.