Abstract
Graphene-based electrical components are inherently multiscale, which poses a real challenge for finite-difference time-domain (FDTD) solvers due to the stringent time step upper bound. Here, a unidirectionally collocated hybrid implicit-explicit (UCHIE) FDTD method is put forward that exploits the planar structure of graphene to increase the time step by implicitizing the critical dimension. The method replaces the traditional Yee discretization by a partially collocated scheme that allows a more accurate numerical description of the material boundaries. Moreover, the UCHIE-FDTD method preserves second-order accuracy even for nonuniform discretization in the direction of collocation. The auxiliary differential equation (ADE) approach is used to implement the graphene sheet as a dispersive Drude medium. The finite grid is terminated by a uniaxial perfectly matched layer (UPML) to permit open-space simulations. Special care is taken to elaborate on the efficient implementation of the implicit update equations. The UCHIE-FDTD method is validated by computing the shielding effectiveness of a typical graphene sheet.
1. Introduction
Graphene, a carbon monolayer with honeycomb pattern, is a promising material to manufacture high-speed transistors and other switching devices thanks to its thinness, its mechanical strength, and its extraordinarily high electron mobility. Consequently, it has attracted lots of research interests, in particular with regard to the development of computational electromagnetic design tools. The finite-difference time-domain (FDTD) method is amid the most popular electromagnetic full-wave solvers. Classically, it discretizes both the electric and magnetic fields in a staggered fashion on a finite cubic grid and approximates the derivatives occurring in Maxwell’s equations by second-order accurate central differences. The algorithm marches on in time by alternately updating the electric and magnetic fields with a time step that is bounded by the cell size. Hence, the thinness of graphene poses a huge computational burden on the conventional FDTD method as it invokes a dense temporal discretization. Also, the spatial staggering gives rise to material-averaging errors near the graphene boundary.
To resolve both problems at once, the unidirectionally collocated hybrid implicit-explicit (UCHIE) FDTD method was proposed in 2D in [1] for good conductors in a multiscale environment and is extended here to 3D problems including graphene sheets. This 3D UCHIE-FDTD method implicitizes the direction perpendicular to the graphene boundary, hereby eliminating the small cell sizes from the time step limit and, consequently, speeding up the simulation. At the same time, it collocates the fields in this direction, which permits an unambiguous definition of the graphene boundary, strongly enhancing the accuracy. Furthermore, the proposed UCHIE-FDTD method wisely combines interpolations and differences as to preserve second-order accuracy even for nonuniform discretizations. Hence, different cell sizes can be adopted inside and outside the graphene sheet without any form of accuracy trade-off. In contrast to [1], where the UCHIE method was confined to a small part in the interior of the simulation domain, open-space simulations now require the construction of a perfectly matched layer (PML) for the UCHIE method. Moreover, the dispersive properties of graphene necessitate an auxiliary differential equation (ADE) formulation of the UCHIE method. These are two important issues that will be tackled here.
In the remainder of this article, the set of continuous equations that need to be solved are first presented in Section 2, including a uniaxial PML. In Section 3, the corresponding update equations are listed. An efficient direct solver based on the twofold use of the Schur complement is next presented in Section 4. Finally, in Section 5, the proposed UCHIE-FDTD technique is validated by computing the shielding effectiveness (SE) of a practical graphene sheet illuminated by a plane wave, an electric dipole, and a magnetic dipole. Section 6 gives a short conclusion and lists some possible future research topics.
2. ADE Formulation
At microwave and THz frequencies, the conductivity of graphene is mainly attributed to intraband electron transitions. In that case, the graphene sheet is typically characterized by a Drude model with effective relative permittivitywhere the volumetric DC conductivity of graphene is (see, e.g., [2])with being the electron charge, the reduced Planck constant, the Boltzmann constant, the temperature, the scattering rate, the chemical potential, and the sheet thickness. The Drude model, which takes the dispersive nature of graphene into account, is incorporated into Maxwell’s equations by means of an auxiliary differential equation describing the time-domain behavior of the graphene current . Thus, we are looking for the solution ofwith and being the vacuum permeability and permittivity, respectively. They constitute the vacuum wave impedance and phase velocityThe substitutionpermits symmetrizing Maxwell’s equations, yieldingIn order to perform simulations in open space, a UPML is added [3]. Thereto, the electric and magnetic field each require one auxiliary equation to model the dispersive PML medium. Moreover, the electric field requires an additional auxiliary equation to prolong the graphene sheet inside the PML. We highlight these auxiliary unknowns by single and double bars. The final set of equations, embracing all these considerations, is where a short notation was adopted for the tensors; for example, with, for the th PML layer,This is the so-called polynomially graded PML with number of layers and power . In this paper, we always choose the standard values and . Also, from [3], we borrow the optimal valuewhere is the spatial step in the -dimension (i.e., , , or ). The UPML is mathematically perceived as a complex stretching of the Cartesian coordinates (and the fields) ensuring impedance matching across its interface. In other words, it is a reflectionless lossy medium. The real-stretch parameter aims to absorb the evanescent waves, whereas the imaginary-stretch parameter is responsible for the absorption of the traveling waves. They are, respectively, one and zero inside the simulation region of interest, and they are polynomially graded according to (9) inside the PMLs normal to the -dimension. In the next section, (7) is discretized, giving rise to the UCHIE-FDTD update equations.
3. Update Equations
In a nutshell, the UCHIE-FDTD method organizes the electromagnetic field components in unidirectionally collocated cells such as the one depicted in Figure 1. Unlike traditional collocated finite-difference methods, the proposed scheme ensures second-order accuracy, albeit at the cost of computationally more involved update equations. Thereto, it retains the conventional central differences, but it evaluates the Maxwell equations that contain a derivative with respect to the direction of collocation (here the -dimension) right in the middle between two adjacent samples in that direction by means of linear interpolations. In order to end up with a consistent update scheme, similar field interpolations are needed in time. This principle is exemplified in Figure 2. In the end, the two field components along the direction of collocation, that is, and , are updated explicitly, whereas the remaining four field components are updated implicitly. In line with conventional ADE schemes, the electric current components have the exact same discretization as their associated electric field components. Hence, we get the following set of explicit update equations:for the magnetic field andfor the electric field. Here, we used the notation
There are two implicit updates (14) and (15). The PEC boundary conditions at the back of the PMLs require the exterior tangential electric fields to be zero. Consequently, for cells in the -dimension, there are electric field and magnetic field samples in (14)-(15). The bracket notation is used to denote diagonal matrices whose rank is or depending on whether they act on electric or magnetic field components, respectively.Besides, we introduced the -interpolatorsand the -differentiatorsAs suggested by the time-indices of the field samples, the time stepping algorithm alternately performs sets of explicit and implicit updates. As one possible valid implementation, we opted to treat the update equations in the same order as they occur in this section. For more insights on the numerical dispersion and stability, the interested reader is referred to [1].
4. Efficient Solution Method
In this section, an efficient strategy to solve the sparse linear system (14) is pointed out. First, the matrix in the l.h.s., which we will call from now on, is partitioned as follows:Its inverse is given by [4, prop. 2.8.7]withand Schur complementwhich is, in turn, partitioned as follows:Analogously, the inverse of the Schur complement is given bysuch that, in the end, solving the large sparse system (14) only requires the inversion of the rank- matrixThe reordering matrixtogether with the rank- primary circulant matrix , that is, the circulant matrix whose first row has second element one and all other elements zero, allows constructing a banded matrixwhich has four elements per row arranged in a staircase pattern. For example, for , its sparsity pattern isThe LU factorization with row and column pivoting of , that is,is obtained at negligible cost. Once all necessary matrices are precomputed, (14) is updated via forward-backward substitutions and some sparse-matrix multiplications. More specifically, the final algorithm to update (14) inside the time stepping loop goes as follows: (1)Compute the vector in the r.h.s. of (14) and split it into two parts and conforming to the partitioning applied in (18).(2)Compute and split this vector into two parts and conforming to the partitioning applied in (22).(3)Use forward and backward triangular sweeps to compute (4)Compute .(5)Compute .For clarity, we omitted the transposition superscripts above, because it should be quite clear that all occurring vectors are column vectors. Note that some matrix products can be computed prior to time stepping. Compared to the conventional HIE-FDTD method [5], the proposed UCHIE-FDTD method requires the inversion of two pentadiagonal matrices of rank instead of two tridiagonal matrices of rank . Hence, the improved accuracy of the UCHIE-FDTD method comes at a slightly higher computational cost.
5. Numerical Examples
5.1. SE for a Plane Wave
The graphene sheet under investigation has , , , and . Hence, according to (2), the sheet has a DC conductivity . In order to validate the stability and accuracy of the proposed unidirectionally collocated HIE-FDTD method, the shielding effectiveness is determined for normal plane-wave incidence on an infinite graphene sheet. Thereto, a similar configuration as the one in [6, 7] is adopted and depicted in Figure 3. A plane wave is generated by a total-field scattered-field (TFSF) surface placed at one side of the graphene sheet, whereas the transmitted field is recorded in one point at the other side of the sheet. The source is a differentiated Gaussian pulse with bandwidth , temporal width , and delay . The spatial invariance in the - and -dimension is fulfilled by imposing periodic boundary conditions (PBCs). Both the transmitted and the back-scattered plane waves are absorbed by UPMLs with . The overall grid is uniformly discretized with steps , whereas the graphene sheet is locally resolved by two small cells of size along the -axis as depicted in Figure 4. The time step equals the Courant limit in the -plane, more specificallywhich is about 19,000 times larger than the time step that would have been used by the classical FDTD method. The recorder is placed beyond the graphene sheet. The simulation performs 10,000 iterations. An auxiliary simulation uses the same configuration but replaces the graphene sheet by a vacuum layer as to record the reference field . The time-domain data are then Fourier transformed and their ratio determines the shielding effectivenessThe analytical solution for the SE is known to be with being the Fresnel reflection coefficient of a single vacuum-graphene interface and wave numberRecall that the effective relative permittivity of the graphene sheet was defined previously in (1). The resulting SEs are plotted in Figure 5. Despite the enormous refinement ratio of 27,000, the analytical and numerical curves are indistinguishable, confirming the accuracy of the UCHIE-FDTD method. The shielding is mainly ascribed to reflections from the two interfaces rather than to the skin effect. At , the skin depth, hereby meaning the thickness that the graphene sheet should have to reduce the amplitude of a plane wave by 63%, is and this value increases towards DC.
5.2. SE for an - and -Dipole
A similar simulation is repeated but this time a dipole of length is excited, at a distance from the left side of the graphene sheet. All six exterior faces of the grid are covered by UPMLs, this time with in order to absorb the evanescent waves radiated by the dipole source as well as the evanescent waves inside the graphene sheet. Indeed, the graphene sheet extends inside the PMLs in order to exclude reflections from its outer edges, such that we are again modeling an infinite graphene sheet. For the -dipole, the SE is again determined by recording and taking the ratio defined in (30), whereas, for the -dipole, the dual definition of the SE is adopted where is replaced by . The recorded fields are plotted in Figure 6 together with the plane-wave data from Section 5.1. The resulting SEs for both the - and the -dipole are shown in Figure 7. Below , which corresponds to , the graphene sheet experiences the reactive near field of the dipole. For this frequency range, the wave impedances of the - and -dipole are approximately and , respectively. The graphene sheet behaves like a classical good conductor at low frequencies. As such, it has a very low wave impedance. Consequently, the -dipole experiences a higher contrast resulting in more reflections, whereas the -dipole does the opposite. This is exactly what we observe in Figure 7: at low frequencies, the -dipole is better shielded by the graphene sheet due to strong reflections, whereas the -dipole does not even notice the graphene sheet because their wave impedances have the same order of magnitude. At higher frequencies, both dipoles are electrically further removed from the graphene sheet and, consequently, they resemble a plane wave as demonstrated by the convergence of the three curves in Figure 7.
6. Conclusion
A novel HIE-FDTD method is proposed that features collocation in the direction of implicitization. This technique, combined with the ADE formalism to incorporate Drude media, is demonstrated to accurately capture the multiscale and dispersive intricacies of graphene sheets in the microwave and THz regime. At higher frequencies such as the near infrared, the contributions from both the intraband and interband interactions to the overall conductivity of the graphene sheet should be taken into account as was done in [7] by means of Padé fitting. Future work focuses on the extension to gyrotropic graphene sheets in the presence of a biasing magnetic field, which is preferably treated with a fully collocated FDTD method such as the one described in [8]. However, an intensive study is required to reduce the computational cost of the occurring matrix inversion, which will be critical for the success of this approach. Also, it should be possible to combine the presented UCHIE-FDTD method with classical Yee-FDTD by means of a correct generalization to three dimensions of the interface condition derived in [1].
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.