Surface plasmons (SPs) refer to the collective electron oscillations at the conductor-dielectric interface (Low and Avouris, 2014). Their outstanding capabilities of manipulating and localizing electromagnetic fields at the subwavelength scale has triggered various applications in nanophotonics in the visible spectral region (Shalaev, 2007; Luk’yanchuk et al., 2010), including photovoltaic device (Atwater and Polman, 2010), biosensing (Xu et al., 1999; Kabashin et al., 2009), integrated photonic devices, (Gramotnev and Bozhevolnyi, 2010; Schuller et al., 2010; Novotny & van Hulst, 2011), thermal radiation control, (Basov et al., 2016) etc., However, when it comes to IR and especially THz bands, SPs in noble metals experience a severe reduction in optical confinement, which hampers its practical applications (Maier, 2007). In contrast, graphene has attracted considerable attention recently as an excellent SPs materials in THz spectral regions owing to the capability of supporting highly confined THz plasmonic modes. As a two-dimensional material, graphene can in addition be electrically tuned by injecting charge carriers so that their plasmonic responses can be controlled dynamically, (Huard et al., 2007), leading to a variety of potential applications in active optoelectronic devices (Luo et al., 2013; Low and Avouris, 2014; Otsuji et al., 2014; Cui et al., 2021).

Despite its importance, the atomic thickness and low free-carrier concentration of graphene still lead to its weak interactions with electromagnetic fields in comparison with its metal counterparts. Inspired by the design concept of metamaterials and metasurfaces, enormous efforts have been devoted to pattern and shape graphene sheets into arrays for enhancing their interaction cross-sections with THz waves (Liu et al., 2012; Thongrattanasiri et al., 2012; Wang, 2012; Yan et al., 2012; Fan et al., 2013; Shen et al., 2014; Fan et al., 2015; Guo et al., 2018). Accordingly, several analytical models including sheet retrieval method (Fan et al., 2015), plasmon wave functions (PWFs), (García de Abajo, 2014; Silveiro et al., 2015; Yu et al., 2017), and coupled-mode theory (CMT) (Guo et al., 2018) have been proposed to optimize THz wave absorption of graphene. In particular, graphene cut-wire array (GCWA) has been recognized as the simplest THz plasmonic metasurface and its absorption is reported to reach the maximum value of 50% (Fan et al., 2015). However, as too many parameters including geometry, carrier concentration, wavelength, and substrate effect are involved in these models, the optimization of a specific GCWA is still case-dependent, and a numerical data bank with lengthy numerical calculations is still needed. For the efficient design of GCWA with optimized THz absorption, an intuitive understanding of their THz absorption mechanism with all parameters incorporated is strongly desired.

In this study, we propose a universal theoretical approach to design GCWAs with maximized THz wave absorption. Specifically, we combine temporal CMT with quasistatic PWFs to derive a formula for calculating THz absorption of the GCWAs. This approach, which we denote as quasistatic coupled-mode theory (QCMT), enables an efficient optimization of the GCWAs absorption by matching its absorption (Γ_abs) and radiative (Γ_rad) decay rates. When reaching the optimal condition, the maximum absorption is only determined by the substrate given by (1−r ₀)/2, with r ₀ the substrate reflection coefficient without the GCWA structure. In particular, we identify that this optimal condition can be fulfilled simply by tuning the duty ratio of the array, regardless of operating frequency or geometrical parameters, such as the rectangle length, width, or array period. The theoretical results are further verified numerically using finite element method (FEM) simulations. It is noted that our results can be applied throughout the entire IR and THz spectral regions, which therefore can guide the design of graphene nano-micro structure arrays for applications in high-performance optoelectronic devices.

The GCWA is schematically shown in Figure 1 (left), which is placed onto a dielectric substrate with permittivity of ε ₂. Considering that our study is focused on the THz spectral region, we assume that the array period (P _x and P _y as shown in Figure 1) is much smaller than the wavelength of interest. Accordingly, only specular transmission and reflection are allowed. The excitation THz waves, which are denoted as | s + 〉 = [ s 1 + s 2 + ] T , illuminate perpendicularly towards the GCWA surface. Similarly, the outgoing THz waves are defined by | s − 〉 = [ s 1 − s 2 − ] T . Symbols s ₁₊ (s ₁₋) and s ₂₊ (s _2-) are complex incident (outgoing) wave amplitudes at the air and the dielectric side. The time-variant mode amplitude, a, of the GCWA can be expressed according to CMT as (Haus, 1984; Fan et al., 2003): { d a d t = j ( ω 0 + j Γ tot / 2 ) a + ( 〈 d | ∗ ) | s + 〉 | s − 〉 = C | s + 〉 + a | d 〉 (1) where Γ_tot represents the total decay rate, which is the sum of radiative decay rate Γ_rad and absorption rate Γ_abs. | d 〉 = [ d 1 d 2 ] T is the coupling vector with 〈 d | d 〉 = Γ rad . Parameter C is the background scattering matrix without the GCWA, which is stated as, C = [ r 0 j t 0 j t 0 r 0 ] (2) where r 0 = ε 2 − ε 1 ε 2 + ε 1 and t 0 = 2 ε 1 ε 2 ε 2 + ε 1 are the port reflection and transmission coefficients, respectively (Haus, 1984) ε ₁ and ε ₂ are permittivities of the substrate and surrounding environment (set as air in our current study), respectively. Note that from the time-reversal symmetry and conservation of energy, we have CC ⁺ = I and C | d 〉 ∗ = − | d 〉 . Thereafter, d ₁ and d ₂ can be calculated as, { d 1 = − Γ rad ( 1 − r 0 ) 2 d 2 = j Γ rad ( 1 + r 0 ) 2 (3)

By solving the differential. Eqs. 1–3, one can obtain, a = ( 〈 d | ∗ ) | s + 〉 j ( ω − ω 0 − j Γ tot / 2 ) (4)

When the THz wave is illuminated from the air side, i.e., | s + 〉 = [ s 1 + 0 ] T , the reflection and transmission coefficients for the GCWA can be written as, { r = s − 1 s + 1 = r 0 + Γ rad ( 1 − r 0 ) / 2 j ( ω − ω 0 ) + Γ tot / 2 t = s − 2 s + 1 = j t 0 − j Γ rad t 0 / 2 j ( ω − ω 0 ) + Γ tot / 2 (5)

Consequently, the total absorption of the GCWA is given by, A = 1 − | r | 2 − | t | 2 = | a | 2 Γ abs | s + 1 | 2 = Γ abs Γ rad ( 1 − r 0 ) / 2 ( ω − ω 0 ) 2 + ( Γ tot / 2 ) 2 (6)

The absorption rate of the graphene array is determined by electron scattering rates, impurity and defect densities, which can be approximated as a constant.(Jablan et al., 2009). Therefore, the maximum absorption of the GCWA can be achieved at the condition of ∂ A / ∂ Γ rad = 0 , yielding the well-known critical coupling condition Γ_rad = Γ_abs. As a result, the theoretical maximum absorption for graphene array becomes (at oscillation frequency ω = ω ₀) A max = ( 1 − r 0 ) / 2 , which only depends on the reflection of the substrate. For GCWA in vacuum, the reflection coefficient r ₀ is 0 and the maximum absorption rate is A _max,vac = 0.5. This result is consistent with literature where the graphene array is suspended in the air (Fan et al., 2015). When graphene array is placed onto silicon oxide substrate ( n S i O 2 = 1.955 and r 0 , S i O 2 = 0.323 ), the maximum absorption is reduced to A m a x , S i O 2 = 0.338 . Specifically, if the graphene array is supported onto a perfect metal, the QCMT reduces to single port condition (r ₀ = −1), and the perfect absorption A _max = 1 can be achieved (Guo et al., 2018).

In the QCMT model, ω ₀, Γ_rad, and Γ_abs are phenomenological parameters that are associated with a specific structure. For the GCWA, we can deduce their explicit expressions based on quasistatic PWFs. (Garcıía de Abajo, 2013; Silveiro et al., 2015; Yu et al., 2017). Specifically, in the long-wave limit (q << k _F, k _F is Fermi wave vector) and high doping condition (E _F >> ℏ ω , E _F is the Fermi energy), the surface conductivity of graphene can be approximated by Drude model, (Hanson, 2008; Low and Avouris, 2014), σ ( ω ) = e 2 π ℏ 2 j E F ω + j τ − 1 (7) where e is the charge of an electron, ℏ is the reduced Planck constant, and τ is the electron relaxation time. When the size of the graphene cut-wire is much smaller than the incidence wavelength, it is convenient to describe its response in the electrostatic limit (Fang et al., 2013; Garcıía de Abajo, 2013; García de Abajo, 2014; Silveiro et al., 2015; Yu et al., 2017) According to the quasistatic PWFs, the transmission of the GCWA can be expressed as, (García de Abajo, 2007; Abd El-Fattah et al., 2019) t = t 0 ( 1 + j β α − 1 − G ) (8) where G = g / S 3 / 2 + j β (g is the dimensionless parameter and for square lattice, g = 4.52), (Thongrattanasiri et al., 2012; García de Abajo, 2014) β = 2 π ω / c n ¯ S is the radiation coupling term, ω is the incidence frequency, c is the speed of light in vacuum, S is the lattice area with S = P _x P _y , n ¯ is the average refractive index of the surrounding medium, α = ζ n 2 ε ¯ D 3 / [ 1 / η ( ω ) − 1 / η n ] is the polarizability of individual graphene cut-wire, ζ n is the mode dipole moment, η n is the mode eigenvalue, η ( ω ) = j σ ( ω ) / ( ω D ε ¯ ) , and n is the order of the plasmon mode. Parameter ε ¯ = ( ε 1 + ε 2 ) / 2 is the average dielectric constant of the surrounding medium and D = L W is the characteristic size of the graphene cut-wire (García de Abajo, 2014; Yu et al., 2017; Fang et al., 2013). For the lowest-order dipolar plasmon mode with electron oscillations along x-axis, n = 1, and only η 1 and ζ 1 are considered in Eq. 8 (Yu et al., 2017) Eq. 8 can be further written as, t = t 0 [ 1 + j Γ rad ω 0 2 / ω − ω − j ( Γ rad + Γ abs ) ] (9) where ω 0 = e 2 E F 4 π 2 ε 0 ℏ 2 ( − η 1 ) ε ¯ D − c n ¯ g 2 π S ε ¯ Γ rad (10) Γ rad = ζ 1 2 E F e 2 2 π ε 0 c n ¯ ℏ 2 D 2 S (11) Γ abs = τ − 1 (12)

For coupled-mode with narrow linewidth (high quality factor Q) and when ω ≈ ω ₀, Eq. 9 is reduced to the same form as Eq. 5 except a universal phase factor. Therefore, Eqs. 10–12 are explicit expressions of the phenomenological parameters used in QCMT model, which can be applied to calculate the THz wave absorption of the GCWAs.

Several observations can be made from Eqs. 10–12. First, according to Eq. 10, the first term inside the square root represents the oscillation frequency of the individual graphene cut-wire, while the second term is accounted for interactions between different cut-wires in the array. Therefore, when the period of the lattice decreases, plasmon frequency undergoes a small redshift due to such interactions. Second, the absorption rate of GCWA keeps constant as Γ_abs = τ ⁻¹, which is independent of the specific geometry and frequency. This is a direct consequence of electrostatic approximation. Third, the radiative decay rate Γ_rad is proportional to E _F, the duty ratio of the array ξ = D 2 / S and the square of the mode dipole moment ζ 1 2 . Importantly, ζ 1 is solely determined by the shape of the cut-wire (Silveiro et al., 2015; Yu et al., 2017; Abd El-Fattah et al., 2019). As will be shown in the following discussion, ζ 1 can be approximated as a constant for the longitudinal plasmonic dipole mode. Finally, by applying Equations (11) and (12) to the critical coupling condition (Γ_rad = Γ_abs), one can obtain, ξ = τ − 1 2 π ε 0 c n ¯ ℏ 2 ζ 1 2 E F e 2 (13)

Eq. 13 is the central result of this work. That is, with substrate refractive index n ¯ , Fermi energy E _F, and electron relaxation time τ at hand, the optimum absorption condition only depends on the duty ratio of the array, regardless of other geometrical parameters including period of the array, length and width, as well as oscillation frequency of the GCWA. Additionally, such a conclusion can be scaled throughout the entire infrared and terahertz band once the electrostatic approximation is valid.

We then carried out FEM simulations (COMSOL MULTIPHYSIC) ¹ to validate the theoretical predictions. To that end, we first verify the critical coupling condition by modeling the optical response of GCWA with different duty ratios (ξ) suspended in the air (n _air = 1) and supported onto SiO₂ substrate ( n S i O 2 = 1.955 ). E _F and the graphene damping rate are set to 0.3 eV and τ ⁻¹ = 0.249 THz as phenomenological parameters (Jablan et al., 2009). The L and W of a graphene cut-wire are 10 and 5 μm, respectively. The periods of the GCWA are P x = L / ξ , P y = W / ξ . The other geometrical parameters used in the simulations are listed in Table 1.

The simulated reflection |r|², transmission |t|² and absorption A spectra of the two GCWAs are plotted in Figure 2A–C and Figures 2E–G. It can be seen clearly that |r|², |t|², and A are strongly dependent on the ξ of the GCWAs. In particular, there exists optimized ξ for the THz absorption of the GCWAs. This can be seen more clearly by plotting the evolution of A at the oscillation frequency of the GCWA as a function of the corresponding ξ (Figures 2D,H, orange squares and lines). Specifically, maximum A can be obtained at ξ = 0.26 and 0.38 for GCWA suspended in the air and supported onto SiO₂ substrate, respectively. The two ξ values for optimum absorption condition are exactly the same as those calculated using Eq. 13, where n ¯ is one and 1.478 = ( 1 + n SiO 2 ) / 2 are employed for air and SiO₂ substrate, respectively. Furthermore, in such a circumstance, Γ_rad is equal to Γ_abs.

To further compare the simulation and QCMT model results quantitatively, Eq. 5 obtained from the QCMT model were then employed to fit the simulated spectra. As shown by the dashed lines corresponding to each ξ , the QCMT model corroborates with the FEM simulation very well. Figures 2D,H depict the fitted absorption rate (Γ_abs, purple) and radiative decay rate (Γ_rad, green), where the graphene damping rate τ ⁻¹ is plotted as a dashed line for comparison. As expected from the QCMT model, Γ_rad is proportional to ξ, while Γ_abs remains as a constant that is close to the intrinsic damping rate of graphene (Eq. 12). At the maximum absorption, the Γ_rad coincides with Γ_abs exactly (Figures 2D,H, black dashed lines). These results indicate that as described by QCMT model, when critical coupling condition is fulfilled the GCWA will exhibit maximum absorption in response to THz wave illumination.

We proceed to validate the analytical expressions of ω ₀ and Γ_rad (Equation (10) and (11)). These two parameters can be calculated based on the mode eigenvalue η ₁ and dipole moment ζ ₁. For graphene cut-wires, they are only determined by the aspect ratio defined as κ = L/W. We first calculated η ₁ and ζ ₁ associated with different κ through mode analysis on an isolated graphene cut-wire (more details on the calculations of η ₁ and ζ ₁ can be found in Supporting Information, Note 1). As shown in Figure 3A and 3D, η ₁ decreases against κ, while ζ ₁ first increases as κ increases and then keeps approximately constant when κ is larger than 1. Because the plasmonic mode is associated with electron oscillations along the x-axis, it becomes a longitudinal mode when κ > 1. Therefore, the dependence of ζ ₁ on κ indicates that the effective dipole moment of the longitudinal dipole mode of the cut-wire stays unchanged. With the knowledge of η ₁ and ζ ₁, we then calculated Γ_rad and ω ₀ using Eqations (10) and (11). Different GCWAs with L and W ranging from 5 to 15 μm and 2–8 μm, respectively, are considered. The lattice period is fixed as P _x = P _y = 24.5 μm. As shown in Figure 3B, Γ_rad increases as either L or W becomes larger. This is reasonable because the area of the cut-wire D ² is equal to L×W. A larger D ² means that more electrons will be driven by the external field, and consequently more energy will be radiated back to the vacuum due to the stronger electron oscillations. The dependence of working frequency ω ₀ on the geometrical parameters is shown in Figure 3E. Clearly, when fixing L (W) while increasing W (L), the ω ₀ increases (decreases) monotonically.

We further simulated the reflection and transmission spectra of the GCWAs with different L and W using the FEM. QCMT model (Equation (5) and (6)) was then employed to fit the simulated spectra and extract the corresponding Γ_rad and ω ₀ (Detailed simulations and fitting results can be found in Supporting Information, Note two and Supplementary Figure S1). As shown in Figures 3C,F, the Γ_rad and ω ₀ obtained from the QCMT model agree well with those deduced from Eqations (10) and (11), which unambiguously confirms the accuracy of the QCMT model we developed.

It should be emphasized that for graphene cut-wires of mico-meter scale characteristic sizes, our QCMT model and the duty ratio-oriented design can in principle cover a vast geometric space and spectral region ranging from the mid-infrared to THz frequencies, as long as the electrostatic limit is still valid. Specifically, given a dielectric substrate of refractive index n _die and a targeting ω ₀, the geometrical parameters of GCWA with maximum absorption can be obtained by using Eqs 10–13 and the numerical results shown in Figures 3A,D. To simplify the discussion, we set ζ ₁ = 0.925 according to Figure 3A, and consequently Eq. 13 can be simplified as, ξ = α n ¯ ℏ τ − 1 E F (14) where α = 80.13 is a dimensionless coefficient. To validate Eq. 14, we consider GCWA with E _F = 0.5 eV, ℏ τ − 1 = 1 meV (τ ⁻¹ = 0.249 THz), and substrate refractive index n = 1. The optimal ξ of GCWA is calculated to be 0.165 . We have numerically simulated three types of GCWAs with different combinations of L, W, P _x , and P _y : a) fixed lattice period (P _x = 4.93 μm and P _y = 2.46 μm) and characteristic size (D = 1.414 μm), but varying graphene cut-wire aspect ratio κ; b) fixed graphene cut-wire (L = 2 μm and W = 1 μm) and lattice area (S = 12.13 μm²), but varying lattice period ratio P _x :P _y ; c) varying all L, W, P _x and P _y parameters. All the GCWAs share fixed duty ratio ( ξ = 0.165 ) with detailed geometrical parameters are listed in Table 2–4. Their corresponding absorption spectra are depicted in Figure 4 . As shown in Figure 4A, the oscillation frequency of GCWAs decreases with the increasing κ. This can be explained by the decrease in η ₁ of individual graphene cut-wire (the first term in Eq. 10). In contrast, for graphene cut-wire with fixed size, the oscillation frequency of GCWA experiences a small redshift when increasing the lattice period ratio. This can be attributed to the mutual interactions between different cut-wires (the increase of g in the second term in Eq. 10 (Thongrattanasiri et al., 2012; García de Abajo, 2014)). Figure 4C shows the absorption spectra for cut-wires arrays with the different combinations of P _x :P _y and κ. Their oscillation frequencies are deliberately chosen to cover 3–28 THz. Nevertheless, all the spectra exhibit maximum absorptions in spite of the geometry diversity, confirming the applicability of the QCMT.

The analyses in Figures 2–4 have validated the accuracy of our developed QCMT model. This model can be employed as an effective tool to design GCWAs with optimized electromagnetic wave absorption. To demonstrate the capability of the QCMT model, we employed it to design three types of GCWAs with optimum electromagnetic wave absorptions. We consider three types of GCWAs: I) E _F = 0.4 eV, ω ₀ = 30 THz, n _die = 1; II) E _F = 0.7 eV, ω ₀ = 6 THz, n _die = 1.955; III) E _F = 1 eV, ω ₀ = 2 THz, n _die = 3.42. The GCWAs I to III are supported in the air, onto SiO₂, and silicon substrates, respectively. In addition, we set the aspect ratios κ of all the three GCWAs to 2. Based on Equations (10) and (11), it is able to deduce the feature size of an individual graphene cut-wire to be L×W = 0.088 μm × 0.044 μm (Type I), 1.63 μm × 0.81 μm (Type II), and 7.92 μm × 3.96 μm (Type III), respectively. According to Eq. 14 the duty ratio of each GCWA can be determined as ξ = 0.2, 0.17, and 0.18, respectively. Afterwards the periods P x = P y = L W / ξ of the GCWAs can be obtained.

To validate the above designs, we numerically simulated the absorption spectra of the three types of GCWAs with different duty ratios. The GCWAs are respectively suspended in the air or supported onto SiO₂ and silicon substrates. The lengths and widths of the three types of GCWAs are employed according to the QCMT designs. As shown in Figure 5A, for the GCWA suspended in the air (Type I), a working frequency at 31.3 THz (mid-infrared regime) can be found with a duty ratio of 0.2. Additionally, such a duty ratio leads to a maximum peak absorption of 0.5, which is consistent with the QCMT calculation and previous study (Fan et al., 2015). For the GCWAs supported onto SiO₂ (Type II) and silicon (Type III) substrates, the maximum absorption occurs for duty ratios of 0.17 and 0.18, respectively (Figures 5B,C). The corresponding working frequencies are 6.23 THz and 2.02 THz, respectively. These results are all corroborated with the QCMT predictions, except for small deviations on the oscillation frequencies (the discrepancies are 4.3% for Type I, 3.8% for Type II and 1% for Type III). We attribute such discrepancies to the deviation of GCWA from quasistatic limitation as well as the high-Q approximation when reducing Equation (9) to Equation (5). Nevertheless, our results unambiguously demonstrate the effectiveness of the QCMT model in designing of GCWAs with optimum electromagnetic absorption in mid-infrared to THz spectral regime.

In summary, we have developed an analytical model, i.e., the QCMT model, for designing GCWA with optimum electromagnetic wave absorption in the mid-infrared to THz spectral regime. In the QCMT model, all of the phenomenological parameters employed to calculate the electromagnetic wave reflection and transmission of a specific GCWA can be calculated analytically, which on one hand provides an intuitive physical insight on these parameters, and on the other hand can help optimize the electromagnetic wave absorption of the GCWA. As a result, by utilizing the QCMT model it is able to facilely obtain the geometrical parameters associated with maximum electromagnetic absorption of the GCWA, which otherwise requires arduous parameter securitizing by numerical simulations. We show that a maximum THz wave absorption can be achieved by choosing proper duty ratio of the GCWA, regardless of its working frequencies and other geometrical parameters. The model developed in our study can in principle generalized to other types of polaritonic 2D crystals, which therefore provide a guideline for the design of 2D micro-/nanostructure arrays operating in the long-wave regime with optimized electromagnetic wave absorption.