Frequency-Selective Surface Based on Negative-Group-Delay Bismuth–Mica Medium

Abstract

Negative group delay may be observed in dispersive media with anomalous dispersion in a certain frequency range. The fact that an outgoing wave packet precedes an incoming one does not violate the causality principle but is only a consequence of a waveform reshaping. This effect is observed in media such as photonic crystals, hyperbolic and epsilon-near-zero metamaterials, undersized waveguides, subwavelength apertures, side-by-side prisms, and resonant circuits at various frequencies. The current work is devoted to the design of a simple negative-group-delay medium with tunable properties in the THz frequency range. This medium consists of a bismuth-based frequency-selective surface on a dielectric substrate and may be tuned both statically and dynamically. While a geometry variation defines a main form of an effective permittivity dispersion and group delay/group velocity spectra, an external voltage allows one to adjust them with high precision. For the configuration proposed in this work, all frequency regions with noticeable change in group delay/group velocity lie within atmospheric transparency windows, which are to be used in 6G communications. This medium may be applied to THz photonics for a tunable phase-shift compensation, dispersion management in systems of THz signal modulation, and for encoding in next-generation wireless communication systems.

1. Introduction

Negative group delay (NGD) is observed for well-shaped waves (pulses) traveling through a medium with anomalous dispersion in a certain frequency range, i.e., group velocity may become superluminal [1,2,3]. In this case, an outgoing wave packet precedes an incoming one, which does not violate the causality principle but is only a consequence of a waveform reshaping associated with the fact that high-frequency components of a wave move faster than low-frequency components [4]. These frequency bands with anomalous dispersion are observed in any dispersive media and are accompanied by strong absorption as a consequence of the principle described by Kramers–Kronig relations [5]. However, attenuation may be prevented during wave propagation in an active medium [6,7]. Group velocity describes only smoothly varying electromagnetic wave packet features such as position of the peak. It should be noted that the superluminal group velocity of a wave packet transmitting through a lossy medium is not equal to superluminal energy/information propagation [8] and may become infinite or negative within anomalous dispersion regions and near such frequency bands.

Group delay (${\tau }_g$) is related to group velocity ($v_g$) as ${\tau }_g=d/v_g$, where d is a medium length (thickness). It may be extracted as

$${\tau }_g=\frac{d}{c}\left[Re\left(\widehat{n}\left(\omega \right)\right)+\omega \frac{\partial Re\left(\widehat{n}\left(\omega \right)\right)}{\partial \omega }\right],$$

(1)

where c is the speed of light in vacuum, $\widehat{n}\left(\omega \right)$ is a refractive index dispersion of a medium under study, $\omega =2\pi f$ is an angular frequency, and f is an ordinary radiation frequency [9]. The concept of group delay is applicable to both spatially extended and spatially negligible media, such as optical fibers and thin filters, respectively. Superluminal group velocity has been observed in anomalous dispersion regions of positive-index media, as well as in photonic crystals, hyperbolic and epsilon-near-zero (ENZ) metamaterials, undersized waveguides, subwavelength apertures, side-by-side prisms, and resonant circuits [10,11,12,13,14,15,16,17].

There are many works devoted to the study of this phenomenon from kHz to GHz and at optical frequencies [18,19,20,21]. Due to the development of THz next-generation communications, it is of great interest to study superluminal propagation in the THz region located in the frequency range of 0.1–10 THz, i.e., between the microwave and mid-infrared regions. Terahertz radiation has a non-ionizing nature and selectively interacts with many dielectric media being absorbed by conductive materials [22]. It has attracted the attention of many researchers and industries due to a recent development of fast-acting and cost-effective sources and receivers of such a radiation and due to its applicability in imaging (medicine, pharmaceutics, security systems, and contactless diagnostics), next-generation communications (6G systems), and spectral analysis [23,24,25,26,27,28]. For research and development in the THz domain, both radiophysical and optical approaches are used. At the moment, important work is underway to create devices and materials to control THz radiation parameters, e.g., for 6G communication systems [29,30]. An application of NGD media in the THz domain allows for tunable phase-shift compensation, as well as dispersion management in systems of THz signal modulation and encoding. To date, several media have been proposed to reach NGD at THz frequencies: periodic bandgap structures (photonic crystals), waveguides close to cutoff frequency, and metamaterials [31,32,33,34]. In the field of metamaterials and metasurfaces, a theoretical work [35] proposed a graphene-based structure that provides a tunable group delay (up to 37 ps) in the THz frequency range; however, these results are have not been experimentally validated, especially considering difficulties with the fabrication of fully graphene structures. Another work [36] proposed a relatively complex frequency-selective surface based on aligned nematic liquid crystal cells, for which an effective permittivity dispersion was analyzed.

In this work, we demonstrate statically/dynamically tunable amplitude transmission and phase-shift spectra, effective permittivity dispersion, and group delay and group velocity spectra both numerically and experimentally on the basis of a simple medium consisting of only two layers: a dielectric substrate and a conductive thin film. We show that all the frequency regions with noticeable change in group delay/group velocity lie within the atmospheric transparency windows, which are to be used in 6G communications. We investigate cases of both continuous and structured media with optical properties tuned both statically and dynamically in the frequency range of 0.2–1.0 THz. A solid mica substrate is used as a dielectric layer, and a thin-film bismuth is used as a conductive layer since its permittivity depends on the thickness, substrate material, temperature, etc. Despite the fact that this medium provides negative group delay in the case of continuous bismuth film (low transmittance), additional features may be introduced by making apertures in a bismuth film with a period and size comparable to the radiation wavelength. Such a structuring allows one to reduce losses and to reach a dual-band or a multiband NGD or to make a medium work in ENZ mode. Finally, we investigate a possibility to adjust an effective permittivity dispersion and group delay/group velocity spectra using an external voltage source. As a result, the structure combines a simple manufacturing process with the possibility of dynamic tuning using printed electrodes and operation in frequency bands that are well-matched with atmospheric transparency windows within the THz frequency range.

2. Methods

2.1. Sample Preparation

The thickness of the bismuth film was chosen toachieve high conductivity and to provide transmission at the same time [37,38]. A thin, almost transparent mica substrate was chosen for THz waves due to its stable geometry and its ability to grow a high-quality bismuth film with a large average crystallit size. Bi film with a thickness of 120 nm on 21 μm mica substrate was synthesized by thermal evaporation (spraying) in a vacuum of ${10}^{-5}$ Torr. A certain portion of a metered substance converted into a gaseous state was deposited on a substrate located above and heated up to 120 °C. Subsequent annealing was performed for 30 min at 250 °C. The thickness of the bismuth film was controlled using a Linnik MII-4 interferometer with an error of 5 nm. The synthesis conditions mentioned above are optimal and therefore allow one to obtain a high-quality thin-film bismuth [39] with a conductivity comparable to that of bulk Bi. It should be mentioned that a classical size effect arises in bismuth films with a thickness of less than about 200 nm, leading to a smooth decrease in their conductivity with decreasing thickness (a quantum size effect occurs under temperatures on the order of liquid nitrogen) [40]. A photo of the continuous-film sample is presented in Figure 1 (inset).

2.2. Measurements, Numerical Simulation, and Material Parameter Extraction

The conventional method of THz time-domain spectroscopy was applied to investigate the optical properties of the sample in the 0.2–1.0 THz frequency range using experimental equipment from Tydex, LLC [41] (see Figure 1). Using a setup (TERA K8, Menlo Systems, Germany) operating in a transmission mode (dynamic range of 65 dB), waveforms of THz pulses that passed through the air, the mica substrate, and the bismuth–mica structure were measured. The instrumental measurement error is estimated to be no more than 5%, including errors in substrate and thin film thickness measurements. To ensure high-frequency resolution of 3.62 GHz, the length of waveforms was set to 276 ps. No filtering was applied in order to save all spectral features of the medium under study. Then, waveforms were converted into complex amplitude spectra, providing information about both amplitude and the phase spectra of the THz wave. Rationing the spectrum of complex amplitudes to a reference spectrum (THz signal passed through the air) yields the transmission spectrum of the structure and the corresponding phase delay. These data are used to calculate dispersions of a complex refractive index ($\widehat{n}\left(f\right)$) and a complex permittivity ($\widehat{\epsilon }\left(f\right)={\widehat{n}}^2\left(f\right)$), whether it is a dielectric substrate (homogeneous medium) or a substrate-film effective medium (heterogeneous) [42,43]:

$$n^{\prime }\left(f\right)=\frac{c\left[{\varphi }_{medium}\left(f\right)-{\varphi }_{air}\left(f\right)\right]}{2\pi fd}+n_{air},$$

(2a)

$$T\left(f\right)=\frac{4n^{\prime }\left(f\right)n_{air}}{{[n^{\prime }\left(f\right)+n_{air}]}^2},$$

(2b)

$$n^{″}\left(f\right)=\frac{cln\left[T\left(f\right)/\left|{\widehat{E}}_{medium}\left(f\right)/{\widehat{E}}_{air}\left(f\right)\right|\right]}{2\pi fd},$$

(2c)

$$\widehat{n}\left(f\right)=n^{\prime }\left(f\right)-i{n}^{″}\left(f\right),$$

(2d)

where ${\widehat{E}}_{medium}\left(f\right)$ and ${\widehat{E}}_{air}\left(f\right)$ are complex amplitudes ($\widehat{E}=\left|\widehat{E}\right|e^{i\varphi }$) of a wave passed through a medium and the air, respectively correspondingly; ${\varphi }_{medium}\left(f\right)$ and ${\varphi }_{air}\left(f\right)$ are corresponding phases; $n_{air}=1$ is the air refractive index; $T\left(f\right)$ is a power transmission coefficient; and i represents the imaginary unit. However, this equation may result in an error for an absorbing medium, since it does not take into account the mutual influence of real and imaginary parts of $\widehat{n}\left(f\right)$ on each other. To obtain an exact value of $\widehat{n}\left(f\right)$, the iterative method that takes into account a complex power transmission coefficient is applied [43,44]. Briefly, an initial value of a complex refractive index (extracted by Equation (2)) is used to calculate a complex transfer function

$${\widehat{H}}_{model}\left(f\right)=\frac{4\widehat{n}\left(f\right)n_{air}}{{\left[\widehat{n}\left(f\right)+n_{air}\right]}^2}·exp\left(-i\left[\widehat{n}\left(f\right)-n_{air}\right]\frac{2\pi fd}{c}\right)·\widehat{FP}\left(f\right),$$

(3)

where the first term is a complex power transmission coefficient, and the last term is responsible for multiple reflections inside the sample (Fabry–Pérot effect):

$$\widehat{FP}\left(f\right)={\left(1-{\left[\frac{\widehat{n}\left(f\right)-n_{air}}{\widehat{n}\left(f\right)+n_{air}}\right]}^{2}exp\left[-2i\widehat{n}\left(f\right)\frac{2\pi fd}{c}\right]\right)}^{-1}.$$

(4)

This transfer function is compared with that measured in an experiment or obtained as a result of simulation (${\widehat{H}}_{measured}\left(f\right)={\widehat{E}}_{medium}\left(f\right)/{\widehat{E}}_{air}\left(f\right)$), and a corresponding correction is applied to values of real and imaginary partsof a complex refractive index:

$$E{R}_m\left(f\right)=\left|{\widehat{H}}_{measured}\left(f\right)\right|-\left|{\widehat{H}}_{model}\left(f\right)\right|,$$

(5a)

$$E{R}_p\left(f\right)=\angle {\widehat{H}}_{measured}\left(f\right)-\angle {\widehat{H}}_{model}\left(f\right),$$

(5b)

$$Re\left[{\widehat{n}}_{new}\left(f\right)\right]=Re\left[{\widehat{n}}_{old}\left(f\right)\right]-s·E{R}_p\left(f\right),$$

(6a)

$$Im\left[{\widehat{n}}_{new}\left(f\right)\right]=Im\left[{\widehat{n}}_{old}\left(f\right)\right]-s·E{R}_m\left(f\right),$$

(6b)

where $E{R}_m\left(f\right)$ and $E{R}_p\left(f\right)$ are magnitude and phase errors, respectively; $\left|{\widehat{H}}_{measured}\left(f\right)\right|$ and $\left|{\widehat{H}}_{model}\left(f\right)\right|$ are measured and modeled magnitude dispersions, respectively; $\angle {\widehat{H}}_{measured}\left(f\right)$ and $\angle {\widehat{H}}_{model}\left(f\right)$ are corresponding phase angle dispersions; ${\widehat{n}}_{new}\left(f\right)$ ${\widehat{n}}_{old}\left(f\right)$ are complex refractive index values on current and previous iterations, respectively; and s is a correction factor typically set to 0.01 to ensure algorithm convergence. The procedure is repeated until the difference between measured and calculated (or simulated) transfer functions is less than a specified error. This function can be used to extract a complex permittivity dispersion of both a substrate and a substrate-film effective medium. Only after this stage is the roughness of the resulting dispersion smoothed using high-order polynomials.

Since a thin bismuth film is placed on a substrate and its permittivity dispersion cannot be measured directly, Tinkham’s method for thin films is applied to extract a Bi film in-plane complex conductivity dispersion [45]:

$${\widehat{\sigma }}_{Bi}\left(f\right)=[({\widehat{n}}_{mica}\left(f\right)+1){\widehat{E}}_{mica}\left(f\right)/{\widehat{E}}_{Bi+mica}\left(f\right){\widehat{n}}_{mica}\left(f\right)-1]/Z_0,$$

(7)

where ${\widehat{E}}_{mica}\left(f\right)$ and ${\widehat{E}}_{Bi+mica}\left(f\right)$ are complex amplitudes of a THz wave transmitted through the mica substrate and the bismuth–mica structure, respectively, and $Z_0$ is the free space impedance. This conductivity dispersion is then converted into an in-plane permittivity dispersion [46,47]:

$${\widehat{\epsilon }}_{Bi}\left(f\right)=1+i{\widehat{\sigma }}_{Bi}\left(f\right)/\left(2\pi f{d}_{Bi}{\epsilon }_0\right),$$

(8)

where $d_{Bi}$ is the thickness of the bismuth film, and ${\epsilon }_0$ is the absolute dielectric constant of the vacuum.

Complex permittivity dispersions of a mica and a bismuth film are used in numerical simulation of a heterogeneous layered structure as described further. The result of this simulation, namely a scattering matrix containing scattering parameters, also provides information on the amplitude and phase components of a THz pulse transmitted through the structure in a frequency range of interest. While amplitude information is obtained directly, a corresponding free-space phase shift is subtracted from a phase component to obtain the phase delay of the medium under study. This simulation is carried out using the frequency domain solver of CST Microwave Studio [48]. The periodic boundary conditions are applied in two directions in the xy plane (sample plane), and linearly polarized THz radiation propagates in the z direction. Results of numerical simulation are converted into a complex effective permittivity dispersion by the same method described above.

3. Results and Discussion

Optical properties of the mica substrate (with a thickness of 21 μm) and the Bi film (with an area of 15 × 15 mm and a thickness of 120 nm) on this substrate were measured using THz TDS in the 0.2–1.0 THz frequency range. Amplitude and phase spectra (Figure 2) were used to extract complex permittivity dispersions of 120 nm Bi film and mica (Figure 3). Figure 2 shows that that Bi-on-mica medium transmits a very small part of radiation in comparison with mica substrate in the entire frequency range under study. The phase delay is positive for mica (normal dispersion) and negative for Bi-on-mica medium (broadband anomalous dispersion) in the frequency range of 0.2–1.0 THz.

Complex permittivity dispersions of 120 nm Bi film and mica were used to calculate an effective permittivity dispersion (${\widehat{\epsilon }}_{eff}\left(f\right)$) of the heterogeneous bismuth–mica structure by numerical simulation in CST Microwave Studio and using the transfer matrix method, obtaining the same results as those reported in [47]. In the case of normal (to the sample plane) incidence of THz radiation, the dispersion relation for this layered structure yields a similar solution for both TE and TM modes:

$$k_0=\frac{2\pi f}{c},$$

(9a)

$${\widehat{k}}_{z,Bi}=k_0\sqrt{{\widehat{\epsilon }}_{Bi}},$$

(9b)

$${\widehat{k}}_{z,mica}=k_0\sqrt{{\widehat{\epsilon }}_{mica}},$$

(9c)

$${\widehat{\eta }}_{Bi}=2\pi f{\epsilon }_0{\widehat{\epsilon }}_{Bi}/{\widehat{k}}_{z,Bi},$$

(9d)

$${\widehat{\eta }}_{mica}=2\pi f{\epsilon }_0{\widehat{\epsilon }}_{mica}/{\widehat{k}}_{z,mica},$$

(9e)

$$\begin{array}{c}\hfill {\widehat{k}}_z=\frac{cos\left({\widehat{k}}_{z,Bi}{d}_{Bi}\right)cos\left({\widehat{k}}_{z,mica}{d}_{mica}\right)}{d_{Bi}+d_{mica}}-\frac{1}{2(d_{Bi}+d_{mica})}\left(\frac{{\widehat{\eta }}_{Bi}}{{\widehat{\eta }}_{mica}}+\frac{{\widehat{\eta }}_{mica}}{{\widehat{\eta }}_{Bi}}\right)·\\ \hfill ·sin\left({\widehat{k}}_{z,Bi}{d}_{Bi}\right)sin\left({\widehat{k}}_{z,mica}{d}_{mica}\right),\end{array}$$

(9f)

$${\widehat{\epsilon }}_{eff}={\widehat{k}}_z^2/k_0^2,$$

(9g)

where $d_{mica}$ is the thickness of the mica substrate, ${\widehat{\epsilon }}_{mica}$ is its permittivity calculated according to Equations (2)–(6), $k_0$ is a free-space wave number, and ${\widehat{k}}_z$ is a Bloch wave number with real and imaginary parts responsible for propagating and evanescent waves, respectively.

The effective permittivity dispersion (${\widehat{\epsilon }}_{eff}\left(f\right)$) was also calculated on the basis of complex amplitude spectra measured for the effective bismuth–mica medium using THz TDS. The results are depicted in Figure 4. There is almost exact coincidence between simulated and experimental dispersions.

In order to impart a frequency-dependent response, periodical cross-like apertures were made in the bismuth film in accordance with one of the standard configurations of frequency-selective surfaces with a polarization-independent electromagnetic response [49]. The scheme of this FSS is depicted in Figure 5a. At fixed film and substrate thicknesses, its electromagnetic response is defined by three geometrical parameters, namely unit cell side length (G) and the length (L) and width (K) of a cross-like slot.

Numerical simulation of the bismuth-based FSS was performed at $G=280$ μm, $L=250$ $\mathsf{\mu }$m, and $K=80$ $\mathsf{\mu }$m to obtain at least two resonances within the frequency range under study. Effective amplitude transmission (Figure 6a) and phase delay (Figure 6b) spectra were extracted to calculate the FSS effective permittivity (Figure 7a,b). The experimental samples (see Figure 5b) were structured using the laser engraving technique with a spot diameter of 30 µm, which was taken into account in both numerical and experimental studies. The experimental curves are presented in these figures, along with simulated curves.

There are two frequency bands with negative phase shift and two bands with positive phase shift associated with corresponding values of an effective permittivity and positive/negative group delay. At $G=280$ $\mathsf{\mu }$m, $L=250$ $\mathsf{\mu }$m, and $K=80$ $\mathsf{\mu }$m, there are three frequency points (at 0.32, 0.68, and 0.83 THz) corresponding to zero effective permittivity.

Excepting the case of continuous Bi film on a mica substrate presented in Figure 4, the experimental and modeling results slightly differ from each other. Laser radiation possibly affects bismuth film permittivity on borders of cross-like slots. The effect may be associated with solid-state regrowth during laser annealing, as discussed in detail in [50]. In addition, the permittivity tensor of anisotropic Bi film [51] may be changed due to the implemented geometrical discontinuity, but in numerical simulation, only the in-plane Bi film permittivity component was used, since it was measured in the experiment. The resulting difference may relate to the component of the dielectric permittivity tensor, which is perpendicular to the sample plane, which was not taken into account in the simulation.

A comparison of THz waveforms measured by THz TDS for the cases of solid Bi film and Bi-based FSS with $G=280$ $\mathsf{\mu }$m, $L=250$ $\mathsf{\mu }$m, and $K=80$ $\mathsf{\mu }$m is presented in Figure 8. The group delay of a wave in mica is positive. In the case of continuous Bi film on a mica substrate, the group delay is negative in relation to both substrate and air signals, which corresponds to superluminal group velocity (the real part of the effective permittivity is negative in the whole frequency range under study, in accordance with Figure 4). In the case of Bi-based FSS on mica, the amplitude transmission is higher, but the group delay is negative in relation to the substrate signal only. In this case, the response of the medium is more complex and should be analyzed on the basis of the effective permittivity dispersion presented in Figure 7a.

Static tuning is achieved by geometric variation. The influence of G, L, and K on amplitude and phase spectra, effective permittivity dispersion, and group delay and group velocity spectra was studied further. In all the investigated cases, a value of the parameter under study was shifted towards both larger and smaller values by 20 $\mathsf{\mu }$m, while the other two parameters remained fixed.

With rising square unit cell size (G), the Q factor should increase, and amplitude transmission should decrease in the entire frequency range, enabling a higher difference between extreme values of the real part of the effective permittivity dispersion (stretching in both directions along the vertical axis), but also leading to an increase in its imaginary part in the entire frequency range under study. The numerical simulation and experimental results are presented in Figure 9 and Figure 10. According to the experimental results, with G rising from 260 $\mathsf{\mu }$m to 300 $\mathsf{\mu }$m, the ENZ point at 0.27 THz moves to 0.34 THz, the second point at 0.70 THz shifts to 0.67 THz, and the third point at 0.83 THz shifts to 0.86 THz.

When the group delay (${\tau }_g$) is higher than that in free space, a wave passes through a medium with subluminal group velocity ($v_g$), and its amplitude fades weakly. When ${\tau }_g$ is less than that in free space, a wave propagates with superluminal $v_g$ exceeding c. If ${\tau }_g<0$, the group velocity ($v_g$) becomes negative, i.e., negative with respect to the phase velocity. In both latter cases, the transmission coefficient falls significantly.

The spectra of group delay and group velocity are presented in Figure 11, along with horizontal lines denoting free-space group delay ($(d_{Bi}+d_{mica})/c$) and the speed of light (in the latter case, velocity is normalized as $v_g/c$). With increasing G, the frequency regions with negative group delay expand, while those with positive group delay narrow. Variation of the structural unit cell size (G) affects group delay in the entire frequency range under study. However, the main effect observed in this case with rising G is the rapidly increasing superluminal $v_g$ near 0.20–0.26 THz (the main jump is between $G=280$ $\mathsf{\mu }$m and $G=300$ $\mathsf{\mu }$m) and its shift near 0.8–0.84 THz.

These frequencies lie within atmospheric transparency windows proposed for use in 6G wireless communication systems [52], which are about 0.2–0.3 THz, 0.4 THz, 0.6 THz, 0.65 THz, 0.8–0.9 THz, etc. [53]. For example, current efforts in the field of 6G communications are focused but not limited on frequencies such as 0.14 THz, 0.22 THz, and 0.34 THz [52].

Variation of the cross-like aperture side length (L) should affect the frequency position of the first resonant peak: with rising L, it should shift towards the low-frequency region (both amplitude and phase spectra). Both real and imaginary parts of the effective permittivity dispersion should also stretch into the low-frequency region, which will lead to a strong change in their values near some fixed frequencies. In addition, the transmission amplitude at the second resonant peak should slightly increase, since rising slot length produces an increase in the cross-sectional area at a constant value of K. The results are presented in Figure 12 and Figure 13.

By varying value of the aperture length (L), it is possible to obtain more broadband ENZ response near 0.8 THz. With rising L, the low-frequency region with NGD tends to narrow, while the region with positive group delay expands, stretching noticeably into the low-frequency range. According to experimental results, with L rising from 230 $\mathsf{\mu }$m to 270 $\mathsf{\mu }$m, the ENZ point at 0.38 THz moves to 0.26 THz, the second point at 0.70 THz shifts to 0.68 THz, and the third point at 0.87 THz does not shift.

The frequency spectra of group delay and group velocity with a variable value of L are presented in Figure 14. With rising L, the group delay spectrum expands mainly to lower frequencies and changes its value depending on the specific frequency. For example, the group velocity is negative when $L=$ 230 $\mathsf{\mu }$m, superluminal for $L=$ 250 $\mathsf{\mu }$m, and subluminal for $L=$ 270 $\mathsf{\mu }$m at 0.24 THz. At the same time, with L rising from $L=$ 230 $\mathsf{\mu }$m to $L=$ 270 $\mathsf{\mu }$m, the operation mode changes from subluminal to near-luminal at 0.87 THz, and a noticeable shift is observed near 0.6 THz. All these frequencies also lie within transparency windows of the atmosphere.

Next, the influence of the cross-like aperture side width (K) was investigated. With rising K, the second resonant peak should shift towards lower frequencies, and at the same time, amplitude transmission should increase at both resonances, since a cut area of the cross-like slot does increase. In other words, with rising K, FSS tends to work as a high-pass filter instead a band-pass filter. The results are presented in Figure 15 and Figure 16. Obviously, increasing the aperture width (K) reduces the Q factor, eliciting a relatively broadband ENZ response near 0.7–0.8 THz by switching from positive to negative group delay in the high-frequency region above 0.63 THz. According to the experimental results, with K rising from 60 $\mathsf{\mu }$m to 100 $\mathsf{\mu }$m, the ENZ point at 0.32 THz moves to 0.36 THz, and the second third points at 0.67 THz and 0.97 THz disappear, since the dispersion in this frequency range becomes positive.

The frequency spectra of group delay and group velocity with variable value of K are presented in Figure 17. With rising K, the group delay changes significantly from negative to positive near 0.8–0.9 THz and falls smoothly near 0.25–0.5 THz. For example, with increasing L, the group velocity switches from negative to superluminal and then to subluminal mode near 0.8–0.9 THz. At the same time, the noticeable shift occurs near such frequencies as 0.2–0.3 THz, and 0.6 THz. These frequencies also lie within transparency windows of the atmosphere.

The dynamical tuning was also studied. The influence of voltage applied between two opposite electrodes on the shift of the effective permittivity dispersion was investigated in the range of 0–10 V (see Figure 18). There is a possibility to adjust effective permittivity dispersion with high precision. When a voltage of 10 V is applied, the maximum range of such a change along vertical axis reaches 1. For example, it allows one to move ENZ point at 0.705 THz to 0.724 THz or to change the effective permittivity real part at 0.304 THz from 0.23 to −0.26. The similar shift of group delay is presented in Figure 19: these dependencies confirm the possibility of a group delay precise adjustment at any desired frequency in the range under study. For example, in the case of $G=280$ $\mathsf{\mu }$m, $L=250$ $\mathsf{\mu }$m, and $K=80$ $\mathsf{\mu }$m, the most noticeable change in this parameter is observed near 0.4 THz, 0.6 THz, 0.8 THz, and 0.95 THz. Most of these frequencies (excluding 0.95 THz) lie within atmospheric transparency windows.

4. Conclusions

In conclusion, statically/dynamically tunable amplitude transmission and phase shift spectra, effective permittivity dispersion, and group delay and group velocity spectra were demonstrated both numerically and experimentally on the basis of FSS consisting of a thin Bi film on a mica substrate. FSS with a polarization-independent configuration is produced by laser engraving, which provides a high-quality pattern. Static tuning is achieved by geometric variation; FSS unit cell size affects its Q factor and transmission coefficient, expanding or narrowing the frequency regions with NGD; cross-like slot length variation shifts the first resonant peak and stretches the corresponding effective permittivity into a low-frequency region; and cross-like slot width influences the second peak position (affecting the operation of the structure as a band-pass or a high-pass filter) and allows one to obtain a wide-band epsilon-near-zero dispersion. Dynamic tuning provides precise adjustment when voltage applied between two opposite electrodes influences the shift of an effective permittivity dispersion and a group delay spectrum. All the frequency regions with noticeable change in group delay/group velocity lie within atmospheric transparency windows, which are to be used in 6G communications. The proposed structure combines a simple manufacturing procedure with the possibility of dynamic tuning using printed electrodes and operation in frequency bands that are well-matched with atmospheric transparency windows and may be applied for tunable phase shift compensation, temporal modulation, dispersion management in systems of THz signal modulation, and for encoding in communications.