Simulation of ageing and wear effect on graphene THz passive components using finite element method

In the growing scenario of 2D material-based metamaterials and metasurfaces for Terahertz (THz) applications, assessing the impact of ageing and wear due to environmental stressors on the components’ performance is becoming mandatory to understand the long-term reliability of novel technologies. This paper introduces approaches to assess the ageing and wear effects on THz passive components through numerical simulations. For this purpose, common techniques for introducing 2D materials and thin metal layers in numerical models are discussed. As a case study, this work explores the effects of graphene degradation and reflective metal ageing on the electromagnetic response of a graphene-enhanced reflective grating for THz absorption and modulation by finite element (FE) analysis. The developed FE model is validated against experimental data obtained through THz Time-Domain Spectroscopy. By computing the device’s transmission, reflection, and absorption spectra for degrading graphene and metal conductive properties, this work provides insights into the influence of ageing and wear on THz passive components.


Introduction
The historical lack of materials suitable for the realization of solid-state components for applications in the terahertz (THz) spectral range (0.1-10 THz frequencies, corresponding to 3 mm-30 µm wavelengths), has hindered the advancement of THz technologies, despite the high potential impact that this range of the electromagnetic (EM) spectrum has in various Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
applications.The development of materials [1] and components such as sources [2,3], modulators [4,5], absorbers [6][7][8][9], detectors [10][11][12], and shields [13][14][15] operating at THz frequencies has recently found new propelling strength thanks to the discovery of metamaterials and metasurfaces specifically designed to be responsive in the THz frequency range.By exploiting resonant microstructures, especially in form or arrays of sub-wavelength unit cells, these artificial materials with geometry-based EM properties enable the interaction with the THz waves, envisaging applications to the detection of the THz spectral fingerprint of chemicals and biomolecules for advanced molecular sensing as well as to the high-frequency signal modulation.
THz metamaterials and metasurfaces are possibly developed through microfabrication of dielectric substrates [16,17], with the addition of metals [18] or semimetals that allow attributing tailored EM response, thus creating passive THz components with fixed spectral response.Including materials with tunable electrical properties, such as graphene or other bi-dimensional materials (2DMs), allows providing the capability to change the THz response by varying the applied electric potential, hence creating active THz components [19][20][21].
Indeed, the forthcoming THz components, both passive and active, allow designing new applications in sensing [22], biosensing [23][24][25], communications [26][27][28], security [29,30], imaging [31], and more.It is therefore plausible that, in the next future, THz components will be used for a wide range of fields, including systems oriented to critical sectors facing demanding environmental conditions, such as aerospace and aviation [32,33].For this reason, it is expected that they will need to withstand critical conditions such as abrupt temperature changes, humidity fluctuations, mechanical stresses, and eventually ionising radiation.As a consequence, the issue of ageing and wear caused by the prolonged exposure to external agents, cannot be neglected.These external agents impair the structural integrity of the components, inducing, for instance, the formation of structural defects and cracks and the delamination of layers, ultimately impairing the THz device performance.
This paper aims at providing the basics of numerical modeling of THz passive components for the computer-based assessment of the effects of ageing and wear, reporting a case study of a graphene-enhanced absorber made by a dielectric grating encompassing a reflective metal coating.In section 2, the models for numerical analysis of ageing and wear effect on THz devices are introduced.In section 3, the finite element (FE) model of a graphene-enhanced reflective grating for the computation of the THz absorption is described, and its performance validated through comparison with THz Time Domain Spectroscopy (TTDS) measurements of a fabricated sample.In section 4, the evaluation of the effects of ageing and wear on the THz performance of the graphene-enhanced reflective grating is proposed.Finally, the conclusions are reported in section 5.

Models for numerical analysis of ageing and wear effect on THz devices
Numerical analyses represent a cost-effective and timeefficient tool for the design and optimization of metamaterials for THz applications, in place of expensive and timeconsuming experimental prototyping [34].Numerical analyses can be performed both at the material-level, to design the unit cell tailoring its THz response, and at the devicelevel, to analyse its performance during operation.Due to the complex structure of metamaterials, and to the dependence of their EM properties on the layout and geometry of the unit cell [35], the numerical analyses of metamaterials and metadevices for THz applications must be based on models with parametrized geometrical properties.For this reason, instead of classical analytical models, space-dependent models for FE and finite difference time domain (FDTD) techniques are generally adopted for this purpose.Both FE and FDTD discretize the volume of the component under analysis by dividing into small elements of size much smaller than the wavelength, and solve linearized Maxwell's equations.FE models can be partitioned by means of different element shapes, enabling the discretization of complex structures, thus providing a great flexibility for modelling arbitrary geometries.The frequencydependent expressions of material properties, such as electrical conductivity, make analysis in the frequency domain easier and faster than time domain analysis.These reasons make FE method more suitable than FDTD for the simulation of devices and components for THz applications.
In the case of 2DM-enhanced metamaterials, the 2DM must be also included in the FE model, with consequent complications related to its very small thickness with respect to that of other dielectric and metal layers involved in the metamaterial structure.Because 2DMs are significantly thinner than any other layer, these numerical methods model the 2DM layer as a zero-thickness layer, i.e. as a surface current boundary condition (SCBC) applied at the interface between the two media surrounding the 2DM [36].The surface current J generated at the interface follows where σ is the surface conductivity of the 2DM and E is the electric field impinging on the surface.For what concerns pristine single-layer graphene, it was demonstrated that its DC conductivity can be obtained through the solution of a Boltzmann quantum equation, and it behaves as a conductive film with collision-dominated transport of free carriers in the low frequency range [37].By increasing the frequency up to the sub-THz range, graphene exhibits Drudelike electrical conductivity, dominated by intraband transitions of charge carriers.Overall, graphene's electrical conductivity, σ g , can be exchanged with optical conductivity and is sum of two contributions: one related to intraband transitions (Drudelike behaviour) and one related to interband transitions [38].
In the far-infrared (IR) range, intraband transitions are initially dominant (∼300 GHz < f < ∼ 3 THz), and graphene's conductivity is well represented by the Drude model [39].By increasing the frequency (∼3 THz < f < ∼ 30 THz), the contribution of interband transitions becomes appreciable, causing slight deviations from the traditional Drude behaviour.In this frequency range, it was demonstrated that the application of a static bias in gated graphene can tune the Fermi level to change carrier doping, thereby shifting the onset of interband transitions [4].When transitioning from far to mid-IR (f ∼ 30 THz), plasmonic effects appear [40].
From mid-to-near-IR (∼30 THz < f < ∼ 300 THz), the contribution of interband transitions to electrical conductivity becomes more relevant, developing increasing absorption that can be modulated from 0 to ∼2.3% through Fermi level modification, by means of an applied electric field [41,42].Finally (f > 300 THz), single-layer graphene exhibits constant absorption of ∼2.3% through all visible range [43].On this basis, in the numerical model of passive graphene components operating in the sub-THz range, an ideally continuous single-layer of graphene can be modelled by a SCBC with the conductivity σ g based on the intraband transitions only, and well described by the Kubo formula [38] where k B is the Boltzmann constant, T is the temperature, e is the electron charge, ℏ is the reduced Planck constant, µ is graphene's chemical potential, i is the imaginary unit, ω is the angular frequency, and γ is graphene's scattering rate.The validity of this formulation of is limited to lowtemperature operation; therefore, it cannot be used to simulate devices for temperature-critical applications.However, its flexibility make it suitable to fit simulations (in which graphene is assumed a perfect continuous single layer) to experiments (in which graphene presents superposition of layers in correspondence of grain boundaries and other defects, such as inclusions and vacancies).
Metals can be modelled with the same approach by using their conductivity σ m expressed by the Drude formula in which ω p is the plasma frequency and τ is the relaxation time.
The deterioration of film integrity, continuity, and homogeneity due to device ageing and prolonged exposure to external factors adversely affects the device's performance by compromising its conductivity.For this reason, the effects of ageing and wear can be included in models based on the thin-layer approximation by changing the conductivity of the thin conductive films.The simplest approach involves varying the conductivity across two extremes: the 'best case', in which the conductivity is not impaired at all (i.e., holds its nominal value), and the 'worst case', in which the conductivity is completely compromised (i.e., reaches zero).More refined approaches include introducing parameters that represent deterioration-induced physical phenomena occurring in the degraded material, and affect specific terms in the mathematical expression of conductivity.Both approaches are followed in the proposed case study, to account for the degradation of the graphene layer and the reflective metal layer of the considered THz metamaterial.

FE modelling of a graphene-enhanced reflective grating for THz absorption
A FE model of a multilayer metamaterial for THz absorbers and modulators was developed for the simulation-assisted design of THz passive components based on graphene by The modelled component is a reflective grating, made by a patterned SiO 2 substrate covered by a gold layer, and coupled to single-layer graphene mechanically supported on both sides by polymethylmethacrylate (PMMA).The pattern size is chosen to interact with THz radiation in the range of frequency f = 0.1-1 THz (λ = 3 mm-300 µm).The metal layer reflects the sub-THz EM radiation towards the graphene layer, to improve the absorption of the metamaterial in the sub-THz frequency range, making it suitable for the design of THz absorbers.
The pattern is composed by straight parallel grooves, obtained by femtosecond laser micromachining of the fused silica substrate, replicated with sub-wavelength spatial period p = 167 µm; the grooves' width is w = 100 µm and the depth is h = 80 µm.A SEM image of the machined grating is displayed in figure 1.
The gold layer, deposited by sputtering technique, has a thickness d g = 30 nm.The PMMA top and bottom layers, with thickness t ps = 150 nm and t pb = 500 nm respectively, are transparent to THz radiation.The pattern profile has an elliptical shape due to the micromachining process.Details about the geometry are in figure 2(a), and the 3D sketch of the different layers composing the considered metamaterial is in figure 2(b).
The fabricated samples have shown great mechanical robustness in time, thanks to the good adhesion between the PMMA layer and plane gold [44] as well as between PMMA and graphene [45].The EM response of the fabricated samples was experimentally analyzed using a commercially available THz time-domain (TD) spectrometer (T-Spec, EKSPLA, Vilnius, Lithuania).The experimental setup is schematically represented in figure 3. The THz-TD setup includes a pumping laser that generates the 145 fs pulses with 1064 nm wavelength and output power of approximately 85 mJ.The laser radiation is split by a beam splitter with a 55:45 power ratio into pump and probe beams.The pump beam, after traversing a fast delay line, excites the gap of the photoconductive antenna (low-temperature grown GaBiAs).Meanwhile, the probe beam traverses a slow delay line and is directed to  For what concerns the FE modelling of the device, due to the sample lateral size much greater than the wavelength in the frequency range of interest, the structure could be considered infinite in the lateral directions, and the FE modelling could be approached by a 2D model of the cross-section of one period of the grating's pattern [46].Symmetry boundary conditions were applied to the model left and right edges, to represent the pattern replication.On top of the structure and below the structure, air gaps were included to represent the propagating medium.The meshing operation was usercontrolled.Specifically, mapped meshing by rectangular elements was chosen for the PMMA layers, and free triangular meshing was used for the remaining part of the structural model.The maximum element size was limited to one tenth of the smallest wavelength considered (i.e. 30 µm), resulting in 4934 vertices and 8376 elements (7176 triangles, 1200 rectangles) of 0.8964 average element quality (scalar value between 0 and 1, where higher is better).
As previously mentioned, graphene was described by a SCBC, as in equation (1), in which the conductivity follows equation (2).Gold was included using the same approach of zero-thickness modelling adopted for the introduction of graphene.However, due to the micromachined substrate roughness, the continuity and homogeneity of the gold layer was impaired.For this reason, the gold layer electrical conductivity was described by refining the Drude model of equation ( 3) with an additional term, i.e. the first term of the Drude-Smith generalization of the Drude model [47], which takes into account the backscattering of carriers resulting from the imperfectness of the reflective surface [48].The equation adopted for gold conductivity, σ Au , is in which c = <cosϑ>, i.e. the average cosine of the scattering angle ϑ formed between the propagation direction of the incoming wave and the wave scattered by the gold grain.
The developed FE model is shown in figure 4.
The finite element analysis (FEA) was performed by using the COMSOL Multiphysics ® RF module for EM waves field analysis.The software allows direct computation of the scattering (S-)parameters by defining input and output ports, where the electric field is either provided or observed.In particular, the software numerically solves the linearized equations to determine the electric field across all nodes of the model, resulting from the meshing of areas.The distribution of the electric field depends on the boundary conditions imposed along model edges and at interfaces between different materials, as well as the material properties attributed to the involved media.Hence, the conductivity defined during the application of the boundary condition described by equation ( 2) determines the transmitted and reflected electric field at the output port.
FEA was carried out on the described model by exciting the input port placed on the top air domain with an electric field with fixed input power, representing the incoming THz radiation, and observing the electric field at the output port placed below the bottom air gap.Then, the frequency-dependent scattering parameters in the range f = 0.1−1 THz were obtained.The power transmission, reflection, and absorption spectra were finally computed according to Following, the model parameters {µ, γ, ω p , τ } were tuned and refined based on the comparison with the measured T and R coefficients, in both cases of electric field polarised parallel and perpendicular to the grating's pattern.Figure 4 shows the results of the FE model simulation compared to the experimental data obtained by THz Time-Domain Spectroscopy (TTDS) of the fabricated sample.
As can be seen in figure 5, the device is mainly an absorber of THz radiation across all the investigated frequency range.The device's absorbance does not demonstrate substantial frequency dispersion in the entire THz range, being 50% at 200 GHz, 55% at 600 THz for parallel polarization and around 30% at 200 and 600 GHz for perpendicular one.These frequencies are chosen to demonstrate the device performance at state-of-the art industrial and R&D communication units, respectively.
The resonating behaviour of the structure visible in figure 5 is caused by the multilayer composition of the metamaterial, presenting distinct patterns or fringes in the EM response due to the constructive and destructive interference of multiple reflections happening at the different interfaces between the layers.The spikes visible between 0.9 and 1 THz are numerical artifacts due to the volume discretization in mesh elements of finite dimension, that cause imperfect model performance at smaller wavelengths.
The simulation results show very good agreement with the measurements, proving the validity of the proposed approach to the modelling of the graphene-based metamaterial for THz passive components.Specifically, the agreement is excellent in the case of parallel polarization.In the case of perpendicular polarization, the simulation underestimates the transmission up to ∼10% in correspondence of the peak values, and the variability of reflection in the higher part of the investigated frequency range, thereby underestimating also the variation of the resulting absorption.Nevertheless, the model well catches the peak values of absorption also at high frequency.The unmatched performance in perpendicular polarization can be ascribed to the microfabrication-related variability of the grooves' lateral size.The model, indeed, assumes periodic grooves with equal aperture and spacing.Deviations from this ideal condition are more relevant when the electric field is oriented orthogonal to the grating (i.e. to the grooves non-uniform lateral size).For this reason, the assessment of ageing and wear effect on the grating's response was performed by assuming parallel polarization of the incoming electric field, when the model has proven greater accuracy.

Ageing and wear effect on the reflective grating's frequency response
The developed FE model was used to evaluate the impact of metal degradation, graphene removal, and structural damage on the EM behaviour of the reflective grating.In this way, the possible effects of material ageing and wear caused by harsh environmental conditions were assessed.The operating condition considered in the simulation is that of input electric field polarized parallel to the grating's pattern, as an exemplary case, because simulations in this condition have a better agreement with the performed measurements.
The effect of graphene degradation was evaluated by computing the T, R, and A spectra in the two extreme conditions of pristine graphene and graphene completely removed by external agents.The effect of gold degradation was assessed by acting on the scattering parameter to account for the increased roughness and inhomogeneity caused by the structural integrity impairment.

Effect of graphene degradation
Pristine graphene can have defects such as vacancies, grain boundaries, and edge irregularities.Environmental conditions, especially ionising radiation, mechanical stress, extreme temperature, can make existing defects more pronounced or introduce new ones.For instance, the edge roughness can be worsened, grain boundaries can move, and contamination by foreign elements can occur.These phenomena generate structural defects that act as scattering centres for charge carriers, significantly influencing graphene's electronic properties by negatively affecting its electrical conductivity.
Extreme environmental conditions and external factors, in addition to normal degradation due to use over the life of the component, can ultimately lead to the complete removal of the graphene layer.In the zero-thickness modelling approach, this condition corresponds to the removal of the surface current boundary condition (SCBC), i.e. the achievement of zero conductivity where graphene was applied.
In between the nominal conductivity condition and the zero conductivity condition, the EM response of the grating changes progressively, moving from the 'pristine graphene' case to the 'removed graphene' case.The power transmission, reflection, and absorption of the metamaterial for THz absorbers was simulated by FEA in these two extreme conditions, and is shown in figure 6.The component's behaviour in the intermediate conditions of 'damaged graphene' is confined within the area delimited by the two extreme cases.
As can be seen in figure 6(a), degradation of the graphene layer causes a progressive increase of the transmission of THz incoming waves by the metamaterial, due to the reduction of the electric conductivity at the interface.Therefore, more energy is capable of passing through graphene and ultimately reaching the other side of the metamaterial.Both T and R exhibit a shift of the resonance peaks (see figures 6(a) and (b)).The frequency of the observed peaks changes as the electrical properties of graphene change.As a result of this different balance between transmitted and reflected energy, the absorption significantly decreases across all the considered frequency range spanning roughly between −22% and −29%, as can be seen in figure 6(c).This is coherent with the concept that graphene cannot absorb the radiation reflected by the goldcovered grating, due to the decreasing electrical conductivity.

Effect of reflective metal degradation
After the destruction of the PMMA/graphene/PMMA multilayer structure on top of the reflective grating, the gold layer is exposed to external agents.Environmental factors, thus, reach the metal and impact also on the structural integrity of gold.
When subject to external factors, gold undergoes morphological modifications.Grain boundaries re-define as the material ages and particles migrate from and towards the substrate, affecting the long-range connections between the islands.Dislocations move due to the applied external stress, and determine local deformation leading to non-uniform electrical properties and eventually interrupting the charge flow.Vacancies and impurities can appear, increasing the scattering of electrons and hindering the efficacy of change transport.Delamination of the thin film may also occur.Overall, the electrical conductivity changes, and the reflective properties of the gold film ultimately worsen.
In the proposed model of the gold-covered grating, the conductivity of gold is modelled by using a Drude-Smith model in place of traditional Drude model, as described in equation ( 4), to take into account the scattering of carriers in opposite direction with respect to the propagation of the incoming wave.For this reason, the modification of the electrical conductivity of the gold layer can be represented by the variation of the backscattering coefficient c.
Ideally, metals have infinite conductivity (i.e. they are perfect electrical conductors).In reality, thin conductive films at high frequency have a conductivity that follows the Drude model (equation ( 3)).This case can be seen as a Drude-Smith model (equation ( 4)) when c = 0, and represents the isotropic scattering condition.In fact, c represents the statistical average of the cosine of the scattering angle.This means that, when carriers are equally scattered in all directions, the average is null and c = 0. Conversely, the stronger the scattering of carriers in the direction opposite to the carrier original direction, the closer the value of c to −1.Indeed, the sign change represents a 180 • phase inversion, meaning that the electron touching the gold interface turns in the opposite direction, hindering the film conductivity.Moving from c = 0 to c = − 1, i.e. using progressively more negative values of c, allows modelling a metal layer with worsening conduction properties.
To observe the effect of the morphological structure modification caused by ageing and wear, the FE model of the grating was simulated by varying the c coefficient around the initial value of c = − 0.97 obtained by fitting the experiments.The initial value of c was obtained by fitting the simulation results to the experimental data of T and R. It is close to c = − 1 due to the irregular surface of the micromachined substrate, which causes the thin film of gold to be already very rough and from the beginning.The rearrangement of the gold grains can either increase or reduce the backscattering, causing the c coefficient to vary in both directions.
Figure 7 shows the change of T, R, and A, computed when the c coefficient varies between c = −0.90 and c = −1.By comparing the simulation results with experiments, it can be noticed that the value of the c coefficient that better fits the measurement is c = −0.97.However, the computed R does not perfectly match the measured reflection, that falls outside of the simulated curves.This mismatch is due to the unavoidable discrepancies between the simplified virtual model and the real component, this latter being affected by imperfections and nonuniformities both in the grating's geometry and the graphene quality.However, the simulations can provide an outlook on the impact of specific parameters, delivering a general overview on the effect of gold imperfectness in terms of increase or decrease of the quantities observed.
As can be seen, increased backscattering (c approaching −1) causes an increase of T (figure 7(a)), a reduction of R (figure 7(b)), and, consequently, a reduction of A across all the observed frequency range, spanning from a −13.5% to a −30.5% decrease (figure 7(c)).This behaviour can be explained by the progressive deterioration of the homogeneity of the thin film of gold, which becomes less conductive as the backscattering increases.As a consequence, the metal layer is less capable to reflect the impinging power, leading to the increase of T. Aside, it can be noticed that the fringes related to the multiple reflections at the different interfaces become more pronounced.

Conclusions
This paper presented a numerical analysis-based approach to the evaluation of the effect of ageing and wear on the spectral response of 2DM-enhanced metamaterials and metasurfaces for THz absorbers and modulators.Possible methods for the computer-based simulation of such THz passive components were discussed, motivating the use of FEA to better address the peculiarities of 2DM-based metastructures.A zero-thickness modeling technique was suggested for representing both 2DM and thin metal layers with Drude-like electrical conductivity by means of SCBC.Two approaches for including the 2DM and metal degradation in the FE model were described, and applied to the case study of a graphene-enhanced reflective grating.
For the considered component, the removal of graphene was modeled by simulating the device in the two extreme conditions of pristine graphene (nominal electrical conductivity) and fully-removed graphene (zero conductivity).The simulation results suggest that the removal of graphene translates in a worsening of the power absorption across all the investigated frequency range.This worsened absorption is mainly due to the increase of transmission caused by the progressive reduction of the surface current intensity at the boundary, dictated by the impairment of graphene electrical conductivity.Removing graphene does not affect the reflection, which is therefore ruled by the gold layer covering the dielectric grating.The degradation of the gold layer was modeled by representing gold electrical conductivity with a Drude-Smith model for the electrical conductivity, and considering fluctuation of the related scattering parameter around the nominal value.The findings show that gold degradation also determines a significant reduction of the overall absorption at all frequencies in the considered interval.
The proposed study provides suggestions for introducing the degradation of 2DM and metal layers in numerical models of metamaterials and metasurfaces for THz passive components.Furthermore, the application of the proposed techniques presented in this paper gives insight into the performance impairment of a graphene-enhanced gold-covered dielectric grating conceived for THz modulation and absorption, as related to the ageing and wear effects on the graphene and gold layer.

Figure 1 .
Figure 1.SEM image of the grating obtained by femtosecond machining of the silicon oxide substrate.

Figure 2 .
Figure 2. (a) Cross-section of the modelled metamaterial, with the description of the geometrical parameters.(b) Exploded 3D view of the layers composing the modelled metamaterial.Both pictures are not to scale.

Figure 3 .
Figure 3. Schematic representation of the THz-TD experimental setup for the characterization of the fabricated samples.

Figure 4 .
Figure 4. FE model developed for the considered multilayer metamaterial.The air domain extends above the structural model for more than ten times the maximum wavelength, to ensure plane wave propagation according to the 2D representation of the problem.As can be seen in the enlarged section, gold and graphene are included in form of surface current boundary condition (SCBC) at the interfaces pointed by the arrows.

Figure 5 .
Figure 5.Comparison between the power transmission, reflection, and absorption spectra, measured by TTDS and simulated by the developed FE model.

Figure 6 .
Figure 6.(a) Transmission, (b) reflection, (c) absorption spectra in the sub-THz frequency range computed in correspondence of the two extreme conditions of pristine graphene and removed graphene.All intermediate conditions of damaged graphene lead to an electromagnetic behaviour within the grey area.