Graphene Hybrid Metasurfaces for Mid-Infrared Molecular Sensors

We integrated graphene with asymmetric metal metasurfaces and optimised the geometry dependent photoresponse towards optoelectronic molecular sensor devices. Through careful tuning and characterisation, combining finite-difference time-domain simulations, electron-beam lithography-based nanofabrication, and micro-Fourier transform infrared spectroscopy, we achieved precise control over the mid-infrared peak response wavelengths, transmittance, and reflectance. Our methods enabled simple, reproducible and targeted mid-infrared molecular sensing over a wide range of geometrical parameters. With ultimate minimization potential down to atomic thicknesses and a diverse range of complimentary nanomaterial combinations, we anticipate a high impact potential of these technologies for environmental monitoring, threat detection, and point of care diagnostics.


Introduction
Distinctive mid-infrared (MIR) molecular vibrations, ranging between~2 and 12 µm, act as characteristic 'molecular fingerprints' for label-free identification of a wide range of chemicals and biomolecules. Usefully, atmospheric transparency windows, at 3-5 µm and 8-13 µm, enable a range of applications such as CO 2 gas sensing and alcohol detection. Within the same MIR spectral range, black body radiation can also be utilised for photodetection and thermal imaging technologies. Taken together, MIR technologies have a substantial role in environmental monitoring, medical diagnosis, and security. However, in comparison to other wavelength regions, there exists a relative lack of MIR sources, detectors, and methodologies.
The current state-of-the-art technologies for MIR sensors are based on semiconductor bulk or quantum structures, such as gallium arsenide (GaAs), indium antimonide (InSb), and mercury-cadmium-telluride (MCT) [1]. Whilst these detectors offer very high sensitivity performance, their wider implementation and dissemination is significantly limited by a combination of high costs, limited spectral range, scarce materials, temperature sensitivity, and need for cooling. As a result, commercially available MIR sensing technologies are typically bulky and expensive, requiring specially controlled operating conditions. Whilst MIR technologies continue to develop, alternative materials and approaches are in high demand to meet the multiple challenges of device sensitivity, spectral range, size, cost, and power consumption [1][2][3].
A new family of low dimensional nanomaterials (graphene, transition metal dichalcogenides, topological insulators) offer unique optoelectronic functionalities and new technological solutions beyond those attainable with conventional semiconductors [4][5][6][7][8][9]. Graphenebased devices, in particular, attracted extensive research and attention due to their unique optoelectronic properties, broadband absorption, and electronic tuneability. For example, graphene shows good potential as atomically thin transparent conductive electrodes, combining high optical transparency (over 97%) for a wide range of wavelengths with

Integration of Graphene with Metal Arrays to Form Hybrid Metasurfaces
Two approaches were implemented for integrating the metal metastructures with monolayer CVD graphene (Graphenea, San Sebastián, Spain): (1) graphene transfer onto pre-fabricated and characterised metal arrays and (2) direct metal deposition on the graphene surface. Transferring graphene allowed direct comparison of optically characterized surfaces and reduced risks related to graphene-metal adhesion. Nominally single layer graphene was acquired from and transferred on to the fabricated and optically characterised metal arrays on SiO 2 substrates by an adhesive transfer method [24]. The presence of graphene layers was confirmed by a combination of contrast enhanced optical microscopy [44], electrical conductivity, and confocal micro-Raman analysis (S&I Spectroscopy & Imaging, Warstein, Germany).
Graphene optoelectronic device structures were patterned by electron beam lithography, graphene plasma ashing, metallisation, and lift-off as described in [45]. Material and device conductivity were measured at room temperature in an electronic probe station and 2450 SourceMeter (Keithley, Cleveland, OH, USA). In the absence of graphene, macroscopic surface conductivity was too low to measure for the Au metasurfaces on SiO 2 , with resistivities of a few Ω on the Au contacts, indicating correctly isolated Au nanostructures.
With the addition of graphene, two-dimensional sheet resistivities of the devices were estimated as~0.7 ± 0.3 kΩ/square. By considering commercial material standardisation and comparing to previous characterisation data, the graphene device mobilities were estimated as~1000 cm 2 ·V −1 ·s −1 , with two dimensional charge carrier densities of

Fourier-Transform Infrared Characterisation of the Metasurfaces
Micro reflection and transmission Fourier-transform infrared spectroscopy (FTIR) measurements were performed by a VERTEX 80v FT-IR spectrometer attached to a HYPER-ION 2000 FTIR microscope (Bruker, Ettlingen, Germany), which allowed to measure both transmission and reflection spectra of the samples. The FTIR spectra were recorded in the range of 5000-500 cm −1 (2-20 µm), with a resolution of 2 cm −1 , measured over areas from 40 × 40 µm 2 to 1 × 1 cm 2 . Polarization measurements used an additional model P03 IR Polarizer (Bruker, Ettlingen, Germany).

Finite-Difference Time-Domain Studies
The applied FDTD code was developed in-house, as described in several previous studies [46][47][48]. The purpose of the development and application of our own FDTD codes was to tackle the considerable difference between the mid-infrared wavelengths (2~10 µm) of interest and the thickness (d = 0.335 nm) of the graphene sheet. Dielectric coefficients of gold were numerically described using Lorentz equation with two poles by fitting refractive index and extinction coefficient data obtained from [49]. Periodic boundary conditions were applied in the x and y directions, while the ends of the calculation domain in the z dimension were simulated by perfectly matched layers (PMLs). Time step duration in the FDTD calculations was 1.83 × 10 −13 s, while a mesh size of 0.012 nm was used, so that detailed features of the electromagnetic wave within the single graphene sheet were properly resolved. The total number of simulated time steps was 16 × 10 5 ; it was sufficient that all fields propagated away from the calculation domain of the FDTD study. The FDTD code required up to 200 GB RAM in a mini-supercomputer Intel(R) Xeon(R) 144 cores, 500 GB RAM, 20 TB HDD.

Geometric Tuneability of Metal Metasurfaces on SiO 2
Metasurfaces, comprised of asymmetric metal nanoantenna arrays, were initially designed based on FDTD studies, which indicated strong interactions with electromagnetic waves, even for sparse metal arrays, with significantly enhanced reflectance (85%), a substantial diffraction (10%), and a much-reduced transmittance (5%) for an array of only 15% surface metal coverage [48]. Importantly, the propagating electromagnetic fields were estimated to be transiently concentrated around the surface nanolayer (e.g., graphene) in a time duration on the order of tens of nanoseconds, suggesting a novel efficient near-field optical coupling.
Since the direct interaction of unpatterned graphene in the midinfrared can be relatively weak, in this approach, we began with fabrication and analysis of metal nanoantenna metasurface on SiO 2 -Si substrates, followed by description and analysis of the integrated hybrid graphene metasurfaces. Figure 1 shows the geometry of one of the studied metastructure arrays measured by scanning electron microscopy. For all structures presented in this study, the designed metal thickness (L z = 50 nm ≈ λ MIR /80), width (L y = 200 nm), and nanogap lengths (g x = P x − L x = 200 nm) remained fixed, whilst the metal length (L x ) and lateral pitch (P y ) were varied. Metal coverages were estimated as M % = (L x × L y )/(P x × P y ).
The MIR photoresponse of the fabricated metal metasurfaces was investigated by micro-FTIR spectroscopy. Figure 2a shows FTIR spectra comparing metasurfaces with widely varying metal coverages, defined by the lateral pitch P y . Two main transmittance minima were observed, defined here as λ 1 and λ 2 . The second transmittance minima, around 9-10 µm (λ 2 ), was also observed for unpatterned reference regions of Si-SiO 2 (SiO 2 layer thickness 192 nm), which was attributed to asymmetric Si-O stretching modes for the sample substrate [50]. However, the interaction strength appeared locally enhanced with increasing coverage density of the metasurface arrays. In contrast, the transmittance minimum, λ 1 , displayed around 4 µm, was absent for the unpatterned Si-SiO 2 regions. This minimum depended strongly on the IR polarization and metal geometry, and can be fully attributed to the patterned metasurfaces. Peak interaction strengths (transmittance minima, reflectance maxima) were observed to increase with the metal coverage in the range of 3% (P y = 10 µm) − 20% (P y = 1.2 µm). For higher density arrays of nanoantenna, the intensity increases appeared to saturate, with a slight broadening of the peak. Nanomaterials 2023, 13, x FOR PEER REVIEW 5 of 15 minimum, λ1, displayed around 4 µm, was absent for the unpatterned Si-SiO2 regions. This minimum depended strongly on the IR polarization and metal geometry, and can be fully attributed to the patterned metasurfaces. Peak interaction strengths (transmittance minima, reflectance maxima) were observed to increase with the metal coverage in the range of 3% (Py = 10 µm) − 20% (Py = 1.2 µm). For higher density arrays of nanoantenna, the intensity increases appeared to saturate, with a slight broadening of the peak.  Here, the metasurface photoresponses λ1 were targeted at the 4.25 µm absorption peak for CO2, with wavelength dependent peak interaction intensities increasing with the metal coverage. (b) Measured dependence of the peak-response wavelengths (λ1, λ2) on the designed metal length (Lx), for different array pitch (Py). λ1 can be attributed to the metal metasurface geometry and increases with metal length Lx whilst decreasing with sub-wavelength metal coverage densities (Pitch, %Metal). The dashed line follows the expected trend for isolated nanoantenna, derived in Section 3.3.
The peak response wavelengths of the main MIR features were investigated as a function of the metasurface geometry. Figure 2b displays FTIR analysis of the peak mid-IR , and nanogap distance (200 nm) were fixed, whilst the lateral pitch P y and horizontal metal lengths L x were varied. From left to right, the density of metal antennas and coverage % are increased, whilst the horizontal metal lengths L x were also slightly increased to compensate for P y -induced blueshifts and targeting CO 2 (4.25 µm) with their peak photoresponse wavelengths. minimum, λ1, displayed around 4 µm, was absent for the unpatterned Si-SiO2 regions. This minimum depended strongly on the IR polarization and metal geometry, and can be fully attributed to the patterned metasurfaces. Peak interaction strengths (transmittance minima, reflectance maxima) were observed to increase with the metal coverage in the range of 3% (Py = 10 µm) − 20% (Py = 1.2 µm). For higher density arrays of nanoantenna, the intensity increases appeared to saturate, with a slight broadening of the peak.  Here, the metasurface photoresponses λ1 were targeted at the 4.25 µm absorption peak for CO2, with wavelength dependent peak interaction intensities increasing with the metal coverage. (b) Measured dependence of the peak-response wavelengths (λ1, λ2) on the designed metal length (Lx), for different array pitch (Py). λ1 can be attributed to the metal metasurface geometry and increases with metal length Lx whilst decreasing with sub-wavelength metal coverage densities (Pitch, %Metal). The dashed line follows the expected trend for isolated nanoantenna, derived in Section 3.3.
The peak response wavelengths of the main MIR features were investigated as a function of the metasurface geometry. Figure 2b displays FTIR analysis of the peak mid-IR Figure 2. (a) Linear polarized micro-FTIR spectra of the fabricated metasurfaces shows two clear mid-infrared photoresponses peaks, λ 1 , λ 2 . Here, the metasurface photoresponses λ 1 were targeted at the 4.25 µm absorption peak for CO 2 , with wavelength dependent peak interaction intensities increasing with the metal coverage. (b) Measured dependence of the peak-response wavelengths (λ 1 , λ 2 ) on the designed metal length (L x ), for different array pitch (P y ). λ 1 can be attributed to the metal metasurface geometry and increases with metal length L x whilst decreasing with sub-wavelength metal coverage densities (Pitch, %Metal). The dashed line follows the expected trend for isolated nanoantenna, derived in Section 3.3.
The peak response wavelengths of the main MIR features were investigated as a function of the metasurface geometry. Figure 2b displays FTIR analysis of the peak mid-IR photoresponse for 38 unique fabricated and characterised metal metasurface arrays on SiO 2 substrates. Both λ 1 and λ 2 displayed geometry dependant tuneability. It was found that the primary geometric variable was the nanoantenna metal length L x , whilst an important perturbation was induced by the lateral pitch P y . The shorter infrared peak wavelength response λ 1 exhibited the strongest dependence on the array geometries. As with λ 1 , transmittance minimum λ 2 was also observed to shift with varying metal length and lateral pitch (L x , P y ). However, the shift in λ 2 was significantly less sensitive to changes in the geometry, only ranging between~9.2 and 10.0 µm for the fabricated structures. The smaller observed blueshift in λ 2 , at high metal coverage, was understood to be the effect of the non-symmetric density-of-states of the Si-O vibration modes in the three-dimensional SiO 2 layer (thickness 192 nm) [51,52]. The enhanced transmittance at λ 2 was the result of the transiently concentrated electromagnetic field around the surface nanolayer which strongly activates the Si-O vibration modes in the SiO 2 layer. Midinfrared photoresponses were qualitatively similar for both peaks from room temperature down to 10 Kelvin, suggesting temperature resilient performance, however with an apparent blueshift of approximately 150 nm at low temperatures.
The characterised photoresponse for both features were also observed to be less critically sensitive to other geometric factors, varied within~50%, such as metal thickness, metal width, longitudinal pitch, nanogaps, and SiO 2 dielectric thickness. Combining this information, enabled reliable reproduction of metasurfaces with optimised and well-defined photoresponse for integration with graphene, and towards devices for targeted molecular sensing applications.

Graphene Metasurface Device Integration and Photoresponse
Two technological approaches were investigated to integrate graphene with the geometrically optimised metal arrays to form hybrid metasurface devices. In the first approach, monolayer CVD graphene was transferred onto the pre-characterised SiO 2 -metal metasurfaces, enabling direct comparison of the photoresponse [24]. Alternatively, the nanofabrication of the metasurfaces was replicated on substrates where graphene already covered the full surface (Figure 3a), allowing simplified and reproducible device processing [45]. With the addition of graphene layers, FTIR spectra were qualitatively similar to uncovered metal arrays on bare SiO 2 , with the persistent presence of λ 1 and λ 2 (Figure 3b). However, with the addition of graphene monolayers, an additional blueshift was observed at λ 1 for each of the measured array structures (Table 1). Significant signal enhancements were observed, of Si-O (at 9.5 µm) and PMMA (at 5.8 µm) peaks of +16-19% and +11%, respectively. Considering that the total metal coverages were only 20%, and that the most active focus area around and between the poles of antenna elements was close to~1% coverage, this represents an order of magnitude enhancement of the local molecular signals. Table 1. Graphene-induced blueshifts in the FTIR peak photoresponses of metasurface arrays following transfer of monolayer CVD graphene (G).

Time-Resolved FDTD Study of the Infrared Pulse Transmission of Graphene
To understand our experimental results, we theoretically investigated infrared pulse transmission through hybrid graphene metasurfaces by the FDTD method. The wavelength range of interest was within the midinfrared range, 2-10 µm, so the intraband conductivity of graphene sheet is adopted [27,53]: where is the Fermi level, in the order of approximately 0.2 eV and τ = 100 fs, and the real and imaginary parts of the relative dielectric coefficient are for the graphene plasmonic resonance to the IR radiation as functions of the wavelength of the IR radiation, where d = 0.335 nm is the thickness of the graphene sheet. Note that ε′ is negative, and the absolute value increases quickly following the increase of the IR wavelength. ε′′ increases with the IR wavelength, and is positive, indicating that the transmission of the IR waves through the graphene sheet is very lossy. Note that in Equations (1) and (2), = 0.2 eV and τ = 100 fs are adopted from [27,53]. In [53], the graphene mobility was 2700 cm 2 V −1 s −1 , whilst the experimental graphene properties in this study were estimated as 1000 cm −2 /V s (at ~1 × 10 12 cm −2 ). By scaling τ proportional to mobility, we obtain a new τ = 100/2.7 fs. Such a modification does not affect Equations (1) and (2), since 1/τ = 1.0/100.0 10 = 1.0 10 s −1 , 1/(100 10 /2.7) = 2.7 10 s −1 , while in the wavelength range of interest, ω = 2π c/λ = 6.28 × 3 10 /4.0 10 = 4.7 10 s −1 is much larger than the τ −1 factor.
The first theoretical study subject was to understand the almost identical FTIR spectra of graphene and removed graphene in Figure 3b. In other words, a single graphene sheet does not affect the optical properties of SiO2-Si. This sounds reasonable, since the wavelengths of the IR waves were much larger than the sub-nano-feature sizes of the

Time-Resolved FDTD Study of the Infrared Pulse Transmission of Graphene
To understand our experimental results, we theoretically investigated infrared pulse transmission through hybrid graphene metasurfaces by the FDTD method. The wavelength range of interest was within the midinfrared range, 2-10 µm, so the intraband conductivity of graphene sheet is adopted [27,53]: where E f is the Fermi level, in the order of approximately 0.2 eV and τ = 100 fs, and the real and imaginary parts of the relative dielectric coefficient are for the graphene plasmonic resonance to the IR radiation as functions of the wavelength of the IR radiation, where d = 0.335 nm is the thickness of the graphene sheet. Note that ε is negative, and the absolute value increases quickly following the increase of the IR wavelength. ε increases with the IR wavelength, and is positive, indicating that the transmission of the IR waves through the graphene sheet is very lossy. Note that in Equations (1) and (2), E f = 0.2 eV and τ = 100 fs are adopted from [27,53]. In [53], the graphene mobility was 2700 cm 2 V −1 s −1 , whilst the experimental graphene properties in this study were estimated as 1000 cm −2 /V s (at~1 × 10 12 cm −2 ). By scaling τ proportional to mobility, we obtain a new τ = 100/2.7 fs. Such a modification does not affect Equations (1) and (2), since 1/τ = 1.0/100.0 × 10 −15 = 1.0 × 10 13 s −1 , 1/(100 × 10 −15 /2.7) = 2.7 × 10 13 s −1 , while in the wavelength range of interest, ω = 2π c/λ = 6.28 × 3 × 10 8 /4.0 × 10 −6 = 4.7 × 10 14 s −1 is much larger than the τ −1 factor. The first theoretical study subject was to understand the almost identical FTIR spectra of graphene and removed graphene in Figure 3b. In other words, a single graphene sheet does not affect the optical properties of SiO 2 -Si. This sounds reasonable, since the wavelengths of the IR waves were much larger than the sub-nano-feature sizes of the graphene sheet (0.335 nm); so, macroscopically, the IR waves should transmit through without perceivable perturbation. However, previous studies reported significant crests and troughs in the IR transmission spectra through single graphene sheet, e.g., [54].
We started with a single graphene sheet in vacuum. The transmittance of the IR plane wave at normal incidence through this single graphene sheet extended in the xy plane positioned at z = 0 is easily calculated by applying Fresnel's equations where E r , E i , and E t denote the electric field of the reflected, incident, and transmitted wave, respectively, between medium 1 denoted by complex refractive index ñ 1 and medium 2 having ñ 2 . Let ñ = n + iκ be the complex refractive index of the thin graphene film (denoted as medium 2), and the measurement is performed in air so that the refractive indices of the spaces above the upper interface of the graphene sheet (medium 1) and below the lower interface of the graphene sheet (medium 3) are 1.0. By Equation (3), the reflection and refraction coefficients at the upper and lower interfaces are for a single reflection and transmission. Here, r 12 is the reflection of the plane wave from medium 1 back to medium 1 reflected by the upper interface of the graphene sheet, t 12 is the refraction from medium 1 into medium 2 at the upper interface. r 23 and t 23 are likewise defined but at the lower interface. Note that r 12 = −r 23 . The series of the transmitted waves are E t E i = e iδ t 12 t 23 1 + β + β 2 + · · · = e iδ t 12 t 23 lim n→∞ n ∑ i=0 β i where δ = ωñd/c 0 and β = e iδ r 23 r 21 . It is easy to see that the result of the above infinite summation is E t E i = e iδ t 12 t 23 1 + e 2iδ r 12 r 23 (6) Since the optical power of the transmitted light is while the incident optical power is S i = 2c 0 ε 0 E i 2 from which we obtain the transmittance T through the thin graphene film The numerically calculated transmittance, by Equation (5) with limited summation over n, then Equation (6) for n = ∞ (the black line marked with "∞") are presented in Figure 4a. For n = ∞, the transmittance is very close to 1.0 due to the extremely thin layer thickness of the graphene sheet (d = 0.335 nm) and appears featureless as a function of the IR wavelength. However, when n is limited, strong oscillation in the transmittance spectrum is observed.  Next, we performed FDTD numerical calculations using in-house FDTD codes [46][47][48]. Numerical results were carefully examined in both space and time domains.
FDTD-calculated transmission spectrum of the graphene sheet is presented in Figure  4b as a function of the number of FDTD simulation steps, where the graphene sheet was placed on the xy plane, and the IR pulse impinged on the graphene structure along the z axis (see Figure 4c). As compared with the black ∞ line in Figure 4a, which is also presented in Figure 4b for direct comparison, the FDTD transmission spectrum oscillated strongly along the simulation time.
The critical aspect about results of Equation (6) and FDTD calculation is that Equation (6) is derived for a time interval of infinite length, i.e., n → ∞ in Equation (5). When we calculated the transmittance as a function of n, as a means to emulate measurements of finite time intervals, Equations (5) and (6) produced crests and troughs in the transmission spectrum (see Figure 4a), similar to the oscillations in Figure 4b. For the single graphene sheet, there existed a very substantial difference between n = 100 and n = ∞. This can be expected since the intensity of the IR wave in the graphene sheet will reduce gradually due to both the absorption (loss, ε′′ ≠ 0, but the real loss was negligibly small due to thin thickness) but mainly the transmission.
This also explains the time-dependence of the FDTD-calculated transmission spectrum presented in Figure 4b. The Ex field of the transmitting electromagnetic field is shown in Figure 4d, indicating that the major electromagnetic field passed through the graphene sheet already at approximately the time step of 10 5 (the time duration of each time step is δt = 1.83 × 10 −13 sec in FDTD); we still observed significant EM field passing through the transmission detector, which caused the strong oscillation in the FDTD-calculated transmission spectrum in Figure 4b. This effect was very significant for the graphene sheet since the dielectric coefficient of the graphene sheet was very large, so that the effective light speed there was much slowed.

FDTD Analysis of the Hybrid Graphene Metasurfaces
Here, we studied periodic gold nanoantenna metasurfaces embedded in SiO2 then covered with a single graphene sheet schematically shown in Figure 5. The two- Next, we performed FDTD numerical calculations using in-house FDTD codes [46][47][48]. Numerical results were carefully examined in both space and time domains.
FDTD-calculated transmission spectrum of the graphene sheet is presented in Figure 4b as a function of the number of FDTD simulation steps, where the graphene sheet was placed on the xy plane, and the IR pulse impinged on the graphene structure along the z axis (see Figure 4c). As compared with the black ∞ line in Figure 4a, which is also presented in Figure 4b for direct comparison, the FDTD transmission spectrum oscillated strongly along the simulation time.
The critical aspect about results of Equation (6) and FDTD calculation is that Equation (6) is derived for a time interval of infinite length, i.e., n → ∞ in Equation (5). When we calculated the transmittance as a function of n, as a means to emulate measurements of finite time intervals, Equations (5) and (6) produced crests and troughs in the transmission spectrum (see Figure 4a), similar to the oscillations in Figure 4b. For the single graphene sheet, there existed a very substantial difference between n = 100 and n = ∞. This can be expected since the intensity of the IR wave in the graphene sheet will reduce gradually due to both the absorption (loss, ε = 0, but the real loss was negligibly small due to thin thickness) but mainly the transmission.
This also explains the time-dependence of the FDTD-calculated transmission spectrum presented in Figure 4b. The E x field of the transmitting electromagnetic field is shown in Figure 4d, indicating that the major electromagnetic field passed through the graphene sheet already at approximately the time step of 10 5 (the time duration of each time step is δt = 1.83 × 10 −13 s in FDTD); we still observed significant EM field passing through the transmission detector, which caused the strong oscillation in the FDTD-calculated transmission spectrum in Figure 4b. This effect was very significant for the graphene sheet since the dielectric coefficient of the graphene sheet was very large, so that the effective light speed there was much slowed.

FDTD Analysis of the Hybrid Graphene Metasurfaces
Here, we studied periodic gold nanoantenna metasurfaces embedded in SiO 2 then covered with a single graphene sheet schematically shown in Figure 5. The two-dimensional metasurface array of gold nanorod antenna was placed on the surface of an insulating SiO 2 /Si substrate. The size of the metal patch was denoted as L x × L y × L z , where L x = 1.4 µm is the metal length in the x direction, L y = 0.2 µm is the metal width in the y direction, and L z = 0.05 µm is the metal thickness in the z direction. The periods of the array are denoted as P x = 1.6 µm and P y = 1.2 µm in the x and y direction, respectively.   (Table 1). Figure 5 shows the transmittance spectra of the three sub-micron structures at three different simulation times. Because of the large dielectric coefficients of both the graphene sheet and the gold nanoantenna as well as the long wavelengths of interest (1∼10 µm), the transmittance spectra shown in Figure 5a strongly oscillated because the transmission detector in the simulation was placed within the residual electromagnetic fields captured around the graphene sheet and gold nanoantenna, diffracted from the z propagation direction to propagate along the x and y directions, which persist a long time since the submicron structures are periodic along the xy plane. Only after ~12 × 10 5 simulation steps, all transmitted EM field passed through the detector and the transmittance spectra converge (see Figure 5c). By comparison, here, we emphasised the technical importance of avoiding numerical artifacts for the system (Figure 5a,b). Similar to Figure 4, the transmittance through the single graphene sheet was almost perfect, save a small reduction due to the reflection by the SiO2-air interface. The addition of the graphene sheet to the gold nanoantenna array blueshifted the transmission minimum from 4.10 µm to 3.63 µm. This blueshift of approximately −470 nm is in close agreement with the experimental observations ( Table 1). Note that due to the simplified wavelength-independent dielectric coefficients of the theoretical models, additional vibrational mode features of SiO2 and Si substrate materials appearing in the experimental spectra from (2-10 µm) are not displayed in the simulated spectra.

Discussion
Similar blueshifts, of around −500 nm (~10%), were observed both from (i) direct comparison by graphene transfer onto pre-characterised metal metasurfaces and (ii) indirect comparison by patterning identical metal metasurfaces on graphene-covered-or bare-SiO2 substrates. A weaker blueshift, but much enhanced reflectance was also observed for the second photoresponse peak λ2 due to the Si-O vibration mode in the SiO2 layer (Table 1).
For the application relevant PMMA encapsulated (200 nm) metal metasurfaces and hybrid-graphene devices (Figure 3), FTIR revealed only a limited additional absorption in Figure 5. FDTD transmittance spectra of parallel (E x , H y ) polarized EM fields. The gold metasurface geometry was L x × L y × L z = 1.4 × 0.2 × 0.05 µm 3 , with pitch P x × P y = 1.6 × 1.2 µm 2 . The number of simulation steps was increased from (a) 4 × 10 5 , (b) 8 × 10 5 , up to (c) 16 × 10 5 to allow all transmitted EM field through the detector. The addition of graphene is associated with a blueshift in the photoresponse, in close agreement with experimental FTIR spectra (Table 1). Figure 5 shows the transmittance spectra of the three sub-micron structures at three different simulation times. Because of the large dielectric coefficients of both the graphene sheet and the gold nanoantenna as well as the long wavelengths of interest (1∼10 µm), the transmittance spectra shown in Figure 5a strongly oscillated because the transmission detector in the simulation was placed within the residual electromagnetic fields captured around the graphene sheet and gold nanoantenna, diffracted from the z propagation direction to propagate along the x and y directions, which persist a long time since the sub-micron structures are periodic along the xy plane. Only after~12 × 10 5 simulation steps, all transmitted EM field passed through the detector and the transmittance spectra converge (see Figure 5c). By comparison, here, we emphasised the technical importance of avoiding numerical artifacts for the system (Figure 5a,b). Similar to Figure 4, the transmittance through the single graphene sheet was almost perfect, save a small reduction due to the reflection by the SiO 2 -air interface. The addition of the graphene sheet to the gold nanoantenna array blueshifted the transmission minimum from 4.10 µm to 3.63 µm. This blueshift of approximately −470 nm is in close agreement with the experimental observations ( Table 1). Note that due to the simplified wavelengthindependent dielectric coefficients of the theoretical models, additional vibrational mode features of SiO 2 and Si substrate materials appearing in the experimental spectra from (2-10 µm) are not displayed in the simulated spectra.

Discussion
Similar blueshifts, of around −500 nm (~10%), were observed both from (i) direct comparison by graphene transfer onto pre-characterised metal metasurfaces and (ii) indirect comparison by patterning identical metal metasurfaces on graphene-covered-or bare-SiO 2 substrates. A weaker blueshift, but much enhanced reflectance was also observed for the second photoresponse peak λ 2 due to the Si-O vibration mode in the SiO 2 layer (Table 1).
For the application relevant PMMA encapsulated (200 nm) metal metasurfaces and hybrid-graphene devices (Figure 3), FTIR revealed only a limited additional absorption in the MWIR region of interest, with some sharp characteristic features around 5.8 µm, and at longer wavelengths. However, a substantial redshift of around 400 nm was observed of λ 1 for these devices.
These effects can have a significant practical impact on the design and functionality of such devices, where the optical response may require precise tuning for suitable operational efficiency. For example, for the application of molecular sensing where the absorption wavelengths can be relatively narrow. The demonstrated influence of the graphene layer, with indications of the range of sensitivity and application, also suggests the possibility of a further dynamic tunability of the system by modifying the graphene properties by electrostatic or electrochemical gating [55]. For a molecular sensing mechanism by peak wavelength shifts, the sharpness of the metasurface-induced MIR peaks could be further sharpened by replacing the metal nanoantenna arrays with high-q dielectrics and by introducing chirality to the geometrical design [36].

Empirical Metasurface Photoresponse Calculator for Reliable Precision Design
The main factors determining the peak photoresponse wavelengths and interaction strength were experimentally observed to be the metal length Lx, the lateral pitch Py, and the material properties at the metasurface. From [56], we can expect the resonance condition for a single one dimensional metal antenna to occur at where m N is an integer denoting the order of the resonance condition, n eff is the effective mode refractive index, and δ L represents the deviation of the effective metal length from the geometrical length L x . Figure 6 shows the dependence of (λ N /L x ) vs. (1/P y ), revealing the overall trend for the metasurface arrays to be approximated by with extracted fitting parameters of A (SiO2-Au) ≈ 4.12, b (SiO2-Au) ≈ 1.07 for metasurfaces on SiO 2 substrates. With the addition of graphene, we find A (SiO2-Au-G) ≈ 3.72, b (SiO2-Au-G) ≈ 1.32. From this, we can estimate n eff(SiO2-Au) ≈ 2.06, and n eff(SiO2-Au-G) ≈ 1.86. This pitch related shift effect can be quite substantial for subwavelength periodicities, e.g., for P y = 0.6 µm, The extraction of accurate fitting relations can further be used to improve the design precision of the metasurface geometries. This simple photoresponse design method was used to define the geometric parameters within this study to target specific wavelengths, for example in Figures 2 and 3. Although we focused here on 4.25 µm, CO 2 , Table 2 indicates the expected design parameters to target a range of other MIR wavelengths and molecules.
The physics of the blueshift due to the insertion of the graphene sheet into the metal metasurface can be interpreted from two perspectives: (1) The electromagnetic wave is strongly perturbed, transiently, by the graphene sheet. This, in turn, affects the free electrons in the metal nanoantenna, resulting in a stronger plasmonic effect and a blueshift; (2) by viewing the graphene and the metal nanoantenna independently, the nanoantenna encloses a space in the form of a resonant cavity. The insertion of the graphene sheet reduces this cavity volume so that the resonant frequency is increased, resulting in a blueshift. Since this affect is expected to be influenced by the graphene charge carrier density, it is anticipated that the spectral photoresponse wavelength of the hybrid graphene metasurfaces can be further tuned or modulated by electrostatic or electrochemical gating [38][39][40].
with extracted fitting parameters of A(SiO2-Au) ≈ 4.12, b(SiO2-Au) ≈ 1.07 for metasurfaces on SiO2 substrates. With the addition of graphene, we find A(SiO2-Au-G) ≈ 3.72, b(SiO2-Au-G) ≈ 1.32. From this, we can estimate neff(SiO2-Au) ≈ 2.06, and neff(SiO2-Au-G) ≈ 1.86. This pitch related shift effect can be quite substantial for subwavelength periodicities, e.g., for Py = 0.6 µm, ∆λ1 ≈ −40%. Figure 6. Dependence of the peak-response wavelengths on the metal length Lx and lateral pitch Py. Strong blueshifts are evident for Py < λ 1 , and where the gold metasurfaces are integrated with graphene. Analysis of these trends enables simple and effective targeted design of the metasurface geometries. Table 2. Estimated metal metasurface geometries for midinfrared molecular targeting, for a range of metal lengths L x and lateral pitch P y , with and without the presence of graphene.

Conclusions
In summary, we demonstrated simple, precise, and reproducible geometrical photoresponse tuning of hybrid graphene metasurface towards targeted molecular sensing. Systematic midinfrared photoresponse characterisation was enabled by a combination of electron beam lithography-based nanofabrication, micro-FTIR spectroscopy, and FDTD studies. Peak-response wavelengths were found to depend most critically on two geometrical parameters; the longitudinal metal nanoantenna length and the lateral pitch between the antenna. Substantial blueshifts were observed and characterised upon the integration of graphene with the metal metasurfaces, and for high density metasurface structures, observed up to~40%. The careful interpretation of more than 100 unique structures enabled the development of a simple and precise set of design tools to ensure geometrically fine-tuned photoresponses for reproducible mid-infrared molecular targeting. The combi-nation of hybrid metasurfaces and their detailed characterisation is important for the next generation of smart portable sensors and lab-on-chip technologies.

Data Availability Statement:
The data presented in this study are available on request from the corresponding authors.