On the terahertz response of metal-gratings on anisotropic dielectric substrates and its prospective application for anisotropic refractive index characterization

Terahertz (THz) frequencies have been drawing increasing attention among the research community to bridge the gap between microwaves and optics. THz radiation can be applied to biomedical imaging,¹ security detection,² medical diagnosis,^3,4 remote sensing,⁵ and nondestructive evaluation applications.^6,7 Terahertz spectroscopy has also been widely used for accurately measuring optical constants of materials in the terahertz region.^8–10 

Since Ebbesen et al.¹¹ proposed the extraordinary optical transmission (EOT) through subwavelength hole arrays, many researchers have extensively studied the electromagnetic wave transmission properties of metallic gratings with subwavelength slits.^12–16 With the terahertz wave impinging the metallic grating, surface currents are induced, which flow toward the slits, causing charge accumulation at the slit sides.¹⁷ Such grating structures over dielectric substrates were analyzed by Minot et al., above and below the limit in which radiative diffraction occurs, showing in the later frequency range a feature-rich transmission spectrum with extraordinary optical transmission explained by a pure interference mechanism.¹⁸ Ferraro et al. also studied such resonances, where particular emphasis was put on a resonant mode (guided-mode resonance) arising from the coupling between the first-order diffraction mode from the grating with modes guided in the dielectric substrate.¹⁹ 

In this work, we analyze the peculiarities of such resonances in gratings with subwavelength slits on electromagnetically thick dielectric substrates. As a difference from previous investigations, particular focus is put on anisotropic dielectric substrates and the effect of this anisotropy on the excited modes. In particular, we show that the resonance associated with the waveguide (WG) TM_1,1 mode is a narrow linewidth resonance, dictated mainly by the out-of-plane index of the substrate, which becomes nearly independent of substrate thickness. Therefore, from the spectral position of this resonance, it is theoretically possible to extract the out-of-plane component of the substrate refractive index with very small uncertainty. Although this method, because of relying on the spectral position of a narrow linewidth resonance, formally enables the extraction of refractive index only at one frequency, for the case of materials with no dispersion, as experimentally studied in this work, the extracted refractive index will represent the “THz refractive index” of the material. Overall, this approach allows extracting either the permittivity for the case of materials with no significant dispersion or otherwise the permittivity at discrete frequencies for dispersive materials under a low-loss approximation. Furthermore, from the transmission amplitude at the resonance, it is generally possible to extract the imaginary part of the refractive index.

Based on this principle, we demonstrate that metallic gratings with sub-millimeter periodicity fabricated on dielectric substrates can be used to aid the extraction of anisotropic refractive index. We perform measurements of such gratings on a series of insulating non-dispersive dielectric substrates, namely, isotropic undoped silicon, (0001) anisotropic sapphire, and anisotropic (010) beta gallium oxide (β-Ga₂O₃), to illustrate this effect. The parallel metallic gratings can be utilized to determine the refractive index and thickness of the silicon substrate and to directly extract the out-of-plane refractive indices of the sapphire and β-Ga₂O₃ substrates. In the case of these two anisotropic substrates, in-plane refractive index extraction can be performed through direct transmission measurements through bare substrates.

The schematic of the analyzed grating structure is shown in Fig. 1(a). The width of the parallel metallic stripes is w, and the period of the stripes is p. The substrates used in this analysis are a single-side polished undoped 4 in. silicon wafer, a double-side polished 2 in. (0001) sapphire wafer, and a single-side polished 1 in. (010) β-Ga₂O₃ substrate. The semi-insulating Fe-doped CZ-grown (Czochralski) β-Ga₂O₃ substrate was provided by Northrop Grumman SYNOPTICS. The chosen grating sizes of each substrate are listed in Table I. The substrate thickness h varies with different substrates, with nominal thickness listed in Table I. The grating thickness d is set to be 340 nm in the simulation. Full-wave electromagnetic simulations were performed utilizing the commercial finite element software COMSOL Multiphysics. Results were validated by simulations in Lumerical FDTD. The simulations were carried out under 2D (x–y) condition. Periodic boundary conditions were applied in the x direction, and perfectly matched layer conditions were applied in the y direction. For the sapphire and β-Ga₂O₃ substrate, their refractive indices were set as a second-rank tensor,

The incident wave, with electric field polarized along the x direction (TM mode), propagates normal to the surface. The metallic gold gratings are assumed to have a finite conductivity of σ_m= 4.1 × 10⁷ S/m.²⁰ 

Metal stacks with target thickness Cr(10 nm)/Au(340 nm) were deposited on the silicon and sapphire substrates using DC sputtering. The same stack was deposited on the β-Ga₂O₃ substrate using e-beam evaporation. Cr was used to improve the adhesion of the metal stack. The metal gratings were patterned using standard photolithography, followed by lift-off. The measured grating thicknesses were 388, 367 and 358 nm for silicon, sapphire, and β-Ga₂O₃ samples, respectively, verified using profilometry and confocal laser scanning microscopy (CLSM).

The transmission spectra of the samples were measured by frequency-domain spectroscopy employing a linearly polarized terahertz excitation at normal incidence. The setup is a commercial diode-laser-driven photo-mixing spectrometer using InGaAs photo-mixers at 1550 nm (TOPTICA photonics). Frequency was swept in the range from 0.2 to 0.6 THz. The frequency scan step size was set to be 0.04 GHz to ensure adequate frequency scan points. All the samples were measured with a terahertz collimated beam with a smaller spot than the sample. Reference measurements through air in the absence of samples were performed before each sample transmission measurement.

We first study the Si sample to showcase the structure's fundamental physics and general properties. Si is a widely studied isotropic material. Then, we analyze a sapphire sample, which constitutes a widely known material that is anisotropic. Finally, we study β-Ga₂O₃, a not very well-known anisotropic material, intending to characterize its permittivity along the b axis. This direction is a principal axis for permittivity in β-Ga₂O₃, wherein refractive index values reported in the literature show considerable variations.

The measured terahertz transmission spectrum for the Si sample is shown in Fig. 2, displayed by the black curve. Five measured transmission peaks with resonant frequencies 0.2323, 0.3168, 0.4112, 0.4344, and 0.5210 THz are selected for refractive index and thickness of Si substrate extraction. The position of the selected resonant peaks is verified by Lorentzian fitting with, in all cases, <0.1 GHz confidence bound interval width. Simulations were performed for this structure to guide our analysis.

Assuming a refractive index of 3.38 and grating dimensions and substrate thickness as listed in Table I, we performed electromagnetic simulations. The mode profiles depicted in Fig. 3(a) indicate that peak 1 and peak 2 are third and fourth-order Fabry–Pérot (FP) 1D cavity TEM modes, for which simplified analytical equations can be obtained. The overall theoretical analysis discussed in the following paragraphs is based on a generic anisotropic substrate and, therefore, applies to all the materials studied in our work. The following analysis requires the knowledge of the reflection and transmission coefficient at the grating interface for normal incident light and scattering light, which do not have simple analytical equations and are strongly dependent on frequency and geometry.²¹ For simplicity, the following analysis only considers pure wave-guided modes without mode coupling through the grating. However, the analysis can explain the phenomenon observed in experimental results and the full-wave simulations to some extent. For a TEM mode in the substrate, the phase change from reflection at the substrate/grating interface can be estimated to be $wπp$⁠, as shown in Fig. S4 in the supplementary material, and the phase change from reflection at the substrate/air interface is zero. At the TEM mode resonance, the round-trip phase change in the substrate $(2πfTEM,N/c0)nxx2h+wπ/p$ should be a multiple of $2π$⁠, where $c0$ is the speed of light in vacuum. Thus, the N-order TEM mode resonant frequency, $fTEM,N$⁠, is

Black dashed lines in Fig. 4(a) represent the third, fourth, and fifth FP TEM modes calculated from Eq. (1), which agree with the full-wave electromagnetic simulations.

For normal incident light with only the E_x field, the gold grating can induce the E_y field (perpendicular to the bottom surface of the grating) and standing waveguide (WG) TM modes with electric field in the y direction. A simplified analytical formula can be calculated from the boundary conditions of the structure. The periodic boundary condition of the grating requires the $Nx$ mode wave-vector in the x direction $kx,Nx$ to follow

The boundary conditions at substrate/air and substrate/grating interfaces impose the $Ny$ order even/odd mode wave-vector in the y direction $ky,even,Ny/ky,odd,Ny$ to follow

The $(Nx,Ny)$ order WG TM even/odd mode resonant frequencies, $fTM,even,Nx,Ny/fTM,odd,Nx,Ny$ thus follow²² 

For $nyy=3.38$ and $p=220μm$⁠, the frequency of higher-order $(Nx≥2)$ modes will be larger than 0.8 THz; thus, only $fTM,even/odd,1,Ny$ modes can be observed in the THz frequency range measured in our study. For convenience, WG mode order refers to the $Ny$ mode order in the manuscript.

As discussed in the supplementary material, mode field amplitude $(Ex,TM,m,1,Ny,Ey,TM,m,1,Ny)$ and wave vectors are related as

The ratio of sensitivity of $fTM,m,1,Ny$ $(m=even,odd)$ with respect to $nyy$ and $nxx$ can be obtained as

The $Ny$ order even modes can accumulate and thus generate maximum transmission due to the constructive interference after the $Ny×2π$ round-trip phase shift.²³ Red dashed lines in Fig. 4(a) represent the first- and second-order WG TM even modes calculated from Eq. (4), and gray dashed lines represent the zero- and first-order WG TM odd modes, matching well with the full-wave electromagnetic simulations. Anti-crossing effect of the coupled FP modes and WG modes,²⁴ as illustrated, e.g., by the peak splitting in the black circle, is observed in Fig. 4(a). As shown in the mode profiles in Fig. 3(a), the E_y distributions for peak 3 and peak 4 can be approximated by $sin(2πpx)sin(πhy)$ and $sin(2πpx)sin(2πhy)$⁠, respectively; furthermore, there are some E_x components in the mode profiles. This suggests that peak 3 is dominantly a first-order WG even mode $fTM,even,1,1$ and peak 4 is dominantly a second-order WG even mode $fTM,even,1,2$⁠, and they are both weakly coupled with fifth-order FP modes to have asymmetrical line shapes as shown in Fig. 2.

At the odd mode resonance, electric fields and corresponding coupling efficiency at the transmitted interface are much weaker than those at the input interface, and as a result, the WG mode energy cannot be radiated through the output port. The simulated mode profiles for the first two transmission dips (at 0.405 and 0.424 THz in the transmission spectra in Fig. 2) are shown in Fig. 3(b). As shown in the mode profiles in Fig. 3(b) and Fig. S5(c) in the supplementary material, the E_y distributions for peak 3 and peak 4 can be approximated by $sin(2πpx)sin(π2hy)$ and $sin(2πpx)sin(3π2hy)$⁠, respectively, suggesting that dip 1 and dip 2 are zero- and first-order WG odd modes. As the substrate thickness is increased, the resonant frequency of all the even and odd waveguide modes $fTM,even/odd,1,Ny$ will converge to $c0nyy×p=403.1GHz$⁠, which agrees with the transmission dip frequency of 403.4 GHz for a grating on a semi-infinite silicon substrate calculated by the full-wave electromagnetic simulation in Fig. 4(b). $c0nyy×p$ can accurately predict the dip frequency for a semi-infinite substrate with different $nyy,nxx$ values as shown in Fig. S7 in the supplementary material. Because of this, for low order modes and electromagnetically thick substrates, resonance frequencies become little affected by the substrate thickness.

As shown in Figs. S2 and S3 in the supplementary material, FP modes are generated by the WG TM modes diffraction at the grating, and WG TM modes are generated by the FP mode diffraction. Using the mode coupling theory, the transmission peaks correspond to coupled modes mixed with both FP and TM modes.²³ Equation (6) shows that the resonance frequency of a pure mode is more sensitive to $nyy$ than $nxx$ when the mode has more $Ey$ field energy. Since the FP mode has no $Ey$ field, the FP mode depends only on $nxx$ and not on $nyy$⁠. We expect the ratio of sensitivity of the coupled FP and WG mode resonant frequency to $nxx$ and $nyy$ can be indicated by the ratio of total $Ex$ field energy $Ux$ and total $Ey$ field energy $Uy$⁠, where

Among the transmission peaks, peak 3 has the highest $Ey$ field energy ratio $αy=Uy(Ux+Uy)$⁠. It is to be noticed that dip 1 has an even higher energy ratio than peak 3, as shown in Fig. 3(b). From this perspective, among the analyzed resonances, peak 3 and dip 1 are dominated by the $Ey$ field and are most sensitive to the perpendicular (out-of-plane) index $nyy$⁠.

Due to the anti-crossing and grating diffraction effects, full-wave electromagnetic simulations are still needed for accurate refractive index extraction. In these simulations, the thickness of the Si substrate was varied from 475 to 525 μm with a step size of 5 μm. The refractive index was varied from 3.20 to 3.60 with a step size of 0.025. Curves representing the combinations of thickness and refractive index for the Si sample matching the spectral position of the measured resonant peaks were extracted from the simulations for the five selected resonant frequencies and are shown in Figs. 5(a)–5(e). Figure 5(f) shows the combined curves arising from the five selected resonant peaks. As shown in Fig. 5(f), the curves intersect at a point corresponding to an Si thickness of 497 μm and a refractive index of 3.390. The simulated transmission spectrum of the Si sample is also shown in Fig. 2, displayed by the red dashed line. From Fig. 2, the resonant peaks in the measured transmission spectrum match well with the simulations. The lower transmission amplitude for the narrow resonances can be attributed to a large sensitivity of the structure to imperfection in fabrication (e.g., lift-off process) as well as to metal loss; in general, losses can be enhanced at resonant peaks with narrow linewidths, thus associated with large-quality factors.¹⁹ The obtained refractive index 3.390 in the THz frequency range is close to the reported value of 3.418 in Refs. 25 and 26. The slight difference between these extracted values, 0.8%, might be the result of various experimental and computational factors, which, although each of them have a very small role, might all together add and lead to error. These factors may include non-perfect normal incidence of the beam and other misalignments, interpolation errors, and inhomogeneities in the sample (to a much less extent).

It is to note in Fig. 5(c) that the position of the first narrow linewidth peak, peak 3, which occurs at 0.4112 THz, is not affected significantly by the substrate thickness. This peak's position mainly depends on the refractive index of the substrate. Therefore, by accurately determining in measurements its spectral position, it is theoretically possible to extract the terahertz refractive index with minimal uncertainty. This is possible because of the narrow spectral linewidth of this resonance in the analyzed structure (which has a relatively large period to gap ratio). For this narrow linewidth peak, one could define sensitivity as the change in the resonance frequency with a change in the substrate's refractive index (⁠$S=Δf/Δn$⁠, where $Δf$ is the change in the resonance frequency and $Δn$ is the change of the refractive index²⁷). The calculated sensitivity is ∼0.1 Thz per refractive index unit (THz/RIU) for peak 3, which is associated with a figure-of-merit $(FOM)$ ∼ 12 (⁠$FOM=S×Q$⁠, where quality factor Q was estimated to be ∼120 for peak 3). Considering a ∼ 0.1 GHz uncertainty in resonance peak frequency determination, this translates to uncertainty levels in index extraction on the order of Δn ∼ 0.001. In practice, numerical errors in the simulations and interpolation errors when generating the contour plots in Fig. 5 are sources of error/uncertainty in this index extraction (together with experimental errors, as mentioned before).

The measured terahertz transmission spectrum for the sapphire sample is shown in Fig. 6, indicated by the black curve. In Fig. 6, peaks 1, 2, and 3 with resonant frequencies 0.2871, 0.4160, and 0.4495 THz are selected to guide the sapphire in-plane n(//) and out-of-plane $n(⊥)$ refractive indices fitting. Here, the in-plane refractive index n(//) indicates the one parallel to the incident wave polarization (perpendicular to the c axis), and the out-of-plane refractive index $n(⊥)$ indicates the perpendicular (parallel to the c axis). Peaks 1, 3, and 4 with resonant frequencies 0.2871, 0.4495, and 0.5695 THz are selected for sapphire thickness and in-plane refractive index n(//) extraction.

Simulated spectral positions of peaks 1, 2, and 3 when varying the in-plane refractive index n(//) and the out-of-plane refractive index $n(⊥)$ for the sapphire substrate are shown in Figs. 7(a)–7(c). The colored lines correspond to the measured resonant frequencies 0.2871, 0.4160, and 0.4495 THz, respectively. The in-plane refractive index n(//) varies from 2.90 to 3.40 with a step size of 0.05, and the out-of-plane refractive index $n(⊥)$ varies from 3.15 to 3.70 with the same step size of 0.05 in the simulations. The thickness of the sapphire substrate was assumed to be 430 μm. Figure 7(d) shows the combined curves associated with the three selected resonances. It is observed that these lines intersect at a point with an in-plane refractive index of 3.130 and an out-of-plane refractive index of 3.380. Figures 7(e) and 7(f) show the simulated curves for the same three selected resonances assuming a sapphire substrate thickness of 420 [Fig. 7(e)] and 440 μm [Fig. 7(f)], respectively. Once again, the lines intersect at a point. The intersection occurs for an in-plane refractive index of 3.205 and an out-of-plane refractive index of 3.380 in Fig. 7(e), while in Fig. 7(f), the intersection point corresponds to an in-plane refractive index of 3.060 and an out-of-plane refractive index of 3.380. From this analysis, we can conclude that regardless of the assumed thickness, the out-of-plane refractive index is ∼3.380. It is to note that the spectral position of the first narrow linewidth peak, labeled as peak 2 in Figs. 6 and 7, is what sets the value of this out-of-plane index. The above simulation analysis verifies that the position of this peak is not significantly altered when varying either the substrate thickness or the in-plane index (as theoretically discussed in the Sec. A). Therefore, we can conclude that accurate out-of-plane refractive index extraction is possible through the proposed approach.

Next, we employ an out-of-plane refractive index value of 3.380, from the previous analysis, for thickness and in-plane refractive index n(//) extraction. Shown in Figs. 8(a)–8(c) are the simulated spectral positions for peak 1, peak 3, and peak 4, respectively, when varying sapphire substrate thickness and in-plane refractive index n(//). The colored curves in Figs. 8(a)–8(c) correspond to the measured resonance frequencies, namely, 0.2871, 0.4495, and 0.5695 THz. In the simulations, the in-plane refractive index was varied from 3.00 to 3.30 with a step size of 0.05, and the thickness of the sapphire substrate was varied from 410 to 450 μm with a step size of 5 μm. The combined simulated curves associated with the measured spectral positions of the three resonances are shown in Fig. 8(d). As observed in Fig. 8(d), the curves rather than intercepting are almost parallel to each other. This has its physical explanation on these resonances arising either as pure Fabry–Pérot modes in the substrate (Peak 1) or resulting from coupled modes wherein the n_xx dependence follows a FP mode (peaks 3 and 4). Therefore, in either case, index and thickness are coupled. So, the thickness of the sapphire sample was mechanically measured, from where a value of $440±4μm$ was obtained. Consequently, an in-plane refractive index of 3.07 was roughly estimated. The final simulated transmission response of the sapphire sample with thickness 440 μm, in-plane refractive index 3.07, and out-of-plane refractive index 3.380 is shown in Fig. 6, displayed by the dashed red line. From Fig. 6, the measured transmission response of the sapphire sample matches well with the simulated transmission response with the obtained anisotropic refractive indices and thickness. The obtained in-plane refractive index (perpendicular to the c axis) 3.07 and out-of-plane refractive index (parallel to the c axis) 3.380 are in good agreement with the refractive indices 3.067 and 3.407 reported in Ref. 28, and 3.064 and 3.406 reported in Ref. 25.

Single crystalline monoclinic structure beta gallium oxide (β-Ga₂O₃) is being extensively studied due to its wide variety of potential applications, such as power switching,^29,30 radio frequency (RF) components,³¹ and UV optoelectronic devices.^32,33 Therefore, it is important to study its refractive indices at varied frequencies ranging from far infrared to deep ultraviolet.^32–34 β-Ga₂O₃ has strongly anisotropic refractive indices along different crystal directions due to its low symmetry monoclinic structure.³⁴ Its photoluminescence as well as THz emission features inherit this anisotropy.

For the analyzed β-Ga₂O₃ sample, which is Fe-doped and thus behaves as a lossless non-dispersive substrate in the frequency range analyzed in this study,³⁵ metallic gratings were fabricated along the a* axis, so the incident electric field is set to be parallel to the c axis under the TM mode.

The measured terahertz transmission spectrum of the β-Ga₂O₃ sample is shown in Fig. 9. As indicated by the red arrow, only peak 1 with resonant frequency 0.3006 THz is selected for crystallographic b axis refractive index extraction following our previous discussions in silicon and sapphire samples. The effects of β-Ga₂O₃ substrate thickness, a*, c, and b axes refractive indices are studied, and the corresponding results are shown in Fig. 10. As shown in Fig. 10(a), the resonant peaks generally shift to a lower frequency range as the substrate thickness increases. However, the selected peak 1 resonant frequency remains almost unchanged. As the incident electric field is parallel to the c axis (perpendicular to the a* axis), the simulated transmission spectra remain unchanged when varying a* axis refractive index values, as depicted in Fig. 10(b). Figure 10(c) displays the simulated transmission spectra with varied c axis refractive index values n(c). The resonant peaks shift to lower frequencies with an increasing n(c) value. However, like in Fig. 10(a), the selected peak 1 resonant frequency remains almost unchanged. Hence, the b axis refractive index can be obtained by matching the position of the simulated peak 1 position to the measured one. Very strong sensitivity of this resonance position to n(b) is observed in Fig. 10(d), which enables once again accurate out-of-plane index extraction. After running a b axis refractive index n(b) sweep with a fine step size of 0.005, the transmission spectrum with obtained n(b) value, 3.195, was found to provide the best match to the selected peak 1 resonant frequency. It is to note that this value is associated with a quasi-static permittivity of 10.21 in good agreement with extractions in bulk samples from direct TDS transmission (10.4),³⁶ polarized CW transmission (10.09 ± 0.23),³⁴ and smaller than those extracted from global ellipsometry measurements (10.6 ± 0.06).³⁴ An important feature of the method employed here is that the extraction of this out-of-plane index component is based on the position of a very narrow linewidth resonance, which, as discussed through this manuscript, is nearly independent of substrate thickness and in-plane indices. Because of this, we believe that the value extracted here is a very reliable representation of the index along the b axis.

To sum up, this paper discussed the general electromagnetic properties of a grating on top of an anisotropic dielectric substrate. From simulations and theoretical analysis, it is shown that various resonant modes can be excited in this structure. For the case of the TM_1,1 guided mode in electromagnetically thick substrates, which manifests as a narrow linewidth transmission peak, the resonance position is mainly dictated by the grating period and out-of-plane refractive index. Neither substrate thickness nor in-plane indices significantly affect the position of this resonance. Thus, we propose a method for refractive index characterization with the assistance of parallel metallic gratings. This method was explored to analyze the properties of three substrates (Si, sapphire, and β-Ga₂O₃).

See the supplementary material for further discussion on the modes of grating diffracted fields, reflection coefficient of grating, effect of finite and infinite thick substrate, and Gaussian beam illumination.

W.J., M.L, P.G., and B.S.-R. acknowledge the AFOSR (No. FA9550-18-1-0332) for primary funding support and the National Science Foundation (ECCS No. 1936729) for partial funding support. A.B. and S.K. acknowledge funding support from the II-VI Foundation Block Gift Program. B.S.-R., A.B., and S.K. also acknowledge support from the AFOSR under Award No. FA9550-21-0078 (Program Manager: Dr. Ali Sayir). Any opinions, findings, and conclusions or recommendations in this material are those of the author(s) and do not necessarily reflect the views of the United States Air Force. This work was performed in part at the Utah Nanofab sponsored by the College of Engineering, Office of the Vice President for Research, and the Utah Science Technology and Research (USTAR) initiative of the State of Utah. The authors appreciate the support of the staff and facilities that made this work possible.

The authors have no conflicts to disclose.

W.J. and M.L. contributed equally to this work.

The data that support the findings of this study are available from the corresponding author upon reasonable request.