{'fraction_non_alphanumeric': 0.02814005404770297, 'fraction_numerical': 0.01186699565268476, 'mean_word_length': 4.827798699075659, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
FDTD subcell graphene model beyond the thin-film approximation

A subcell technique for calculation of optical properties of graphene with the finite-difference time-domain (FDTD) method is presented. The technique takes into account the surface conductivity of graphene which allows the correct calculation of its dispersive response for arbitrarily polarized incident waves interacting with the graphene.The developed technique is verified for a planar graphene sheet configuration against the exact analytical solution. Based on the same test case scenario, we also show that the subcell technique demonstrates a superior accuracy and numerical efficiency with respect to the widely used thin-film FDTD approach for modeling graphene. We further apply our technique to the simulations of a graphene metamaterial containing periodically spaced graphene strips (graphene strip-grating) and demonstrate good agreement with the available theoretical results.


Introduction
Graphene has emerged as a candidate for a large number of novel applications in the field of nano-optics. Unique optical properties of graphene are determined by a 2D nature of its conductivity, with one of the most notable consequences being the existence of ultra-subwavelength graphene surface plasmon polaritons on its surface in the THz frequency range [1]. Coupled with the fact that the conductivity of graphene is easily controlled by electrical gating, this effect can be utilized to create tunable graphene-based THz metamaterials [2].
Though the electromagnetic problem for the continuous graphene sheet can be solved analytically [3,4], in most cases graphene-based metamaterials with complicated geometry require numerical methods. Nowadays, there are three popular computational electromagnetic methods for graphene simulation: finite-difference time-domain method (FDTD), finite element method (FEM), and method of moments (MoM). Each method has its own advantages and disadvantages depending on the specific problem (see the review [5]). Compared to frequency-domain methods FDTD being time-domain has several advantages, among which less computing memory and simulation time for obtaining spectra in a broad frequency range. There are dif-ferent approaches to model graphene in FDTD. In the most common approach, 2D graphene sheet is represented on a 3D mesh by a thin film with its thickness of the order of a single mesh step. For the wavelengths of interest in the THz range (100 microns) the mesh step should be of the order of several microns. However, it was shown that in order to provide a satisfactory approximation of a 2D sheet conducting properties, the mesh step should be superfine (tens of nanometers or less) [6,7,8], which makes the task computationally challenging. To overcome this problem some modified FDTD methods for graphene have been developed, among which the locally one-dimensional method [9], various subcell methods [10,11], and the surface boundary condition method [12].
In this paper we present a subcell technique for conducting 2D materials which employs a tensor representation of the frequency-dependent dielectric permittivity. This technique was developed for interfaces of dispersive materials [13] and implemented in the framework of EMTL library [14]. Unlike other subcell methods it formally allows for arbitrary orientation of the interface normal, not requiring the interface to be aligned with the FDTD grid. We verify the proposed technique by calculating the optical properties of a free-standing continuous graphene sheet at oblique incidence against the analytical solution. We then perform a head-to-head comparison of the subcell technique with the thin-film approach in identical simulation setups. We show that the proposed technique accurately accounts for the interaction of arbitrarily polarized light with the graphene sheet, in contrast to the thin-film approach. We further compute the transmittance, reflectance and absorbance of a graphene strip-grating metamaterial at oblique incidence and demonstrate an excellent agreement of the results with the available semi-analytic solution [15]. where e is the electron charge, c is the speed of light in vacuum, E F is the Fermi energy of electrons in graphene, ω is the frequency, and Γ is the phenomenological electron relaxation rate (damping) [16]. E F can be controlled by applying the gate voltage to the graphene, which in turn leads to the change of the conductivity σ graph [2].
Let us consider an oblique incidence of a plane wave on a continuous graphene sheet. The graphene sheet sandwiched between two half-spaces 1 and 2 filled with dielectric media ε 1 and ε 2 respectively, as shown in Fig. 1. The respective wavevectors of the wave in the half-spaces 1 and 2 are k 1,2 = ω c √ ε 1,2 . The incident wave gives rise to surface currents at the sheet, which can be expressed as: where J graph is the surface current density at the graphene sheet, and E(0) is the electric field at the graphene sheet. Taking into account the relation (2), we can write the boundary conditions for the electric (E) and magnetic (H) fields in case of TE-and TM-polarized incident waves respectively: These boundary conditions are then used to derive the Fresnel coefficients r g and t g of a graphene sheet sandwiched between the two dielectric half-spaces [3,4]. In case of the TE-polarized incident wave we obtain: where θ 1 is an incident angle, and θ 2 is a transmittance angle.
For TM-polarized wave these coefficients take form: (5) and (6) can be used to obtain analytically the reflectance and transmittance properties for the considered geometry for any given angle of incidence and the conductivity of the graphene.
An alternative so called 'thin-film' approach to the the problem depicted in Fig. 1 consists in substituting the graphene surface with a thin film. In this case, the thin film possesses an effective dielectric function, which mimics the properties of the graphene [6]: where d f ilm is a film thickness. The reflectance and transmittance can then be calculated analytically with Airy formula.
In order to estimate the range of applicability of the thin film approach we now consider a special case of a free-standing continuous graphene sheet possessing a typical value of E F = 0.3 eV. In our case the electron relaxation rate Γ ≪ E F is of the order of several meV [16] and it slightly affects the results, so we can set Γ = 0. Let us compare the exact analytical solution (5) and (6) for this case with the thin-film approach. The results are presented in Fig. 2 for a range of the incident angles. In our calculations we took a typical thickness d f ilm = 5 nm, which is used in FDTD simulations of the graphene [6].
We immediately see that while the thin-film approach gives accurate results in case of TE-polarization, it yields incorrect results in case of TM-polarization. Particularly, an unphysical resonance appears which becomes stronger as the incident angle is increased. The reason for this resonance to occur lies in unphysical normal current components arising within the thin film when the electric field in the incident wave has non-zero normal components (which is the case in TM-polarized wave). In the subsequent sections we describe the subcell technique which allows to avoid this effect.

Subcell graphene model
As shown in the previous section, the thin-film representation of graphene leads to substantial differences in the reflection and transmission of the sheet compared to the actual 2D material. In [13] we used the surface normal direction to construct an effective dielectric tensor for mesh cells in the vicinity of an interface. In this approach, the graphene sheet may be regarded as such an interface, introducing additional dispersion along the surface directions. Following the notation from [13], the effective inverse dielectric permittivity tensor has the following form [17,18,19,20]: where P is the projection matrix P ij = n i n j onto the normal n to the graphene sheet. Here the averaging <> should be performed over the characteristic mesh cell volume δv, and assuming uniform bulk medium it may be obtained analytically:ǫ where ǫ || is the 2D graphene permittivity, and δv graph is the volume fraction of the mesh cell, "occupied" by the graphene sheet. The value of δv graph depends on the discretization thickness of the graphene sheet and the surface fraction of the sheet residing inside the cell δv graph = d graph s cell . Note that the discretization thickness d graph may be chosen smaller than the typical mesh step, thus reducing δv graph . The solution of discrete Maxwell's equations with the tensor (9) is then obtained by the auxiliary differential equation (ADE) method in the way described in [13].
Applying the subcell approach (9) to the graphene sheets aligned with the (rectangular) Yee mesh is straightforward and leads to accurate solutions, as will be shown by the results presented in the next Section. In this case we select δv graph = δv and the electric field components aligned with the graphene sheet are updated with the graphene dispersive ADEs, whereas the normal components are nondispersive. Note that in the thin-film approach the normal components of the electric field inside the sheet (represented by one or more mesh cells) are also governed by the dispersive equations, which appears to lead to inaccuracy in the solution.
Here we must remark that although the tensor method formally allows for arbitrary direction of the 2D interface, its application to the inclined 2D sheets requires special care. Our preliminary tests have shown, that variations in δv graph between adjacent mesh cells, typical for the inclined geometry, lead to large discontinuities in the polarization currents drastically reducing the accuracy. Therefore we leave the study of the inclined subcell model for future publications, concentrating currently on the comparison of the thin film and the subcell approaches in the meshaligned case, but for the oblique wave incidence.

Continuous graphene
In this section the FDTD simulation of oblique incidence of a plane wave on a graphene sheet in vacuum is described. The material parameters of graphene are the same as in the previous sections (E F = 0.3 eV, Γ = 0). The compari-  5) and (6)).
son of the subcell technique with the thin-film approach is carried out. The results are also compared to the analytical solution ( (5) and (6)).
In the simulations, the incoming oblique broadband plane wave impulse spans the near-IR to THz frequency range (500 -3000 cm −1 ). We use the previously developed iterative oblique incidence technique [21] to obtain a broadband response for a given angle of incidence in a single simulation run.
In simulations with the developed subcell technique, the choice of the FDTD mesh step value is governed by the frequency range of interest as well as the geometrical features of considered structures. In our particular case of a planar graphene sheet the mesh step ∆ value is determined solely by the frequency range. We found that the reflectance at a given frequency converges to within 2 · 10 −3 at ∆ subcell = 50 nm with respect to the result for ∆ subcell = 10 nm. The results for TM-polarization are shown in Fig. 3 (circles). The relative error of the reflectance calculations with respect to the analytical result is within 4 · 10 −2 for the considered frequencies (below 3000 cm −1 ), which corresponds to the absolute error less than 10 −5 in the calculation of reflectance. Stability is en-sured by performing prolonged simulations for 10 5 time steps (1.7 ps). Typical simulations were run for 10 4 time steps, which took about 10 s per one iteration of the oblique incidence method [21] at Intel Core i5 6200 processor.
On the other hand, the thin-film approach dictates the choice of a much more refined mesh step along the dimension normal to the graphene surface at least. In the reference simulations with the thin-film approach we used ∆ thin−f ilm = 5 nm, with d f ilm = ∆ thin−f ilm . This leads to a dramatic slow-down of the simulations by the factor of 240 with respect to the subcell technique (about 2400 s per one oblique iteration at Intel Core i5 6200 processor). The results and comparison with the subcell technique as well as the analytical solution are shown in Fig. 3 (crosses). The relative error of the simulations with respect to the analytical result is greater than 0.1 for θ = 10 degrees, rising to the order of 1 and even more for greater angles (Fig. 3).
The magnitude of error in case of the thin-film approach is due to two factors: a) the unphysical normal component of polarization current as discussed in previous section, b) the used mesh step is insufficient for the considered film thickness (single mesh step per film). In order to analyze the second factor, we computed the reflectance spectra of  the graphene sheet assuming a thicker (d f ilm = 20 nm) film in the thin-film approach, but leaving ∆ thin−f ilm = 5 nm. At d f ilm /∆ thin−f ilm = 4 mesh steps per film thickness, results appear to reproduce the analytical result for the thin film obtained with Airy formula (Fig. 4). This result also indicates that a finer mesh is required in order to accurately reproduce the analytical thin-film result. For TE polarization, the accuracy of the both methods is comparable. The subcell method is still much more preferable due to much higher computational performance.

Graphene strip-grating
We further demonstrate the application of the subcell method to the modeling of a oblique incidence on a graphene metamaterial. We consider a graphene stripgrating (Fig. 5) possessing multiple surface plasmon resonances in THz range. These resonances can only be excited with a TM-polarized light (when the electric field has a component perpendicular to the strip axis) [2]. Dependence of the resonance positions on the width of the strips makes this structure suitable to serve as a tunable polarizationselective filter. The grating has the period of 70 µm and the strip width of 20, 40, and 60 µm. The ambient medium is vacuum. The TM-polarized plane wave is incident on the structure at an angle of 30 degrees.
While the width of the strips is of order of the wavelength of the impinging light, the arising resonances are strongly localized near the grating surface. Moreover, they possess very sharp field gradients at the edges of the strips [15]. This places strong requirements on the FDTD mesh step size, which should be as fine as 100 mesh steps per strip width at least [22]. The convergence tests for 60 µm strips at the normal incidence indicate that the transmittance spectra convergence to within 1% relative error in 0-10 THz range at ∆ = 250 nm. This mesh step was used in the oblique incidence calculations.
The obtained spectra (Fig. 6) reveal the transmittance dips and absorbance peaks associated with the surface plasmon resonances. The results are consistent with the exact semi-analytical solution for this structure [15] based on the Fourier expansion method. However, this semianalytical method assumes the numerical solution of the integral equations and requires using regularization techniques [15,23]. Although such semi-analitical numerical methods allow modeling of more complicated structures than the strip-grating, such as corrugated graphene [24] or even biperiodic graphene metasurfaces [25], for each new configuration they demand to rewrite all the equations and to rebuild the numerical solution. In contrast, the subcell method allows using FDTD as a universal package for the analysis of arbitrary geometrical configurations and propagation directions, including analysis of propagation along graphene sheets.

Conclusions
This paper presents a subcell technique for modeling graphene within the framework of FDTD method. The proposed subcell technique accurately reproduces the surface conductivity of graphene in frame of the FDTD calculations, which is proved by simulating an oblique incidence of a plane wave on a free-standing continuous graphene sheet. It is shown that the subcell technique demonstrates a superior accuracy with respect to the thin-film approach. In particular, analysis of TM-polarization case shows that thinfilm approach becomes inaccurate due to the unphysical polarization current component normal to the graphene surface. It is also shown, that the subcell method does not require a fine mesh leading to a drastic performance improvement of FDTD simulations of graphene containing structures as compared to the thin-film approach. The subcell technique is further verified by simulating a graphene stripgrating metamaterial at oblique incidence. It is shown that the results are consistent with the available semi-analytical solution. We conclude that the subcell technique presents a flexible and efficient approach to modeling of graphene metamaterials and can be applied to optimize their optical properties for various applications.