# Full-Wave Analysis of the Shielding Effectiveness of Thin Graphene Sheets with the 3D Unidirectionally Collocated HIE-FDTD Method

Authors: Arne Van Londersele; Daniël De Zutter; Dries Vande Ginste
Journal: International Journal of Antennas and Propagation
 (2017)

---

## 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.


---

## Body

## 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 Section2, 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 permittivity(1)εrgω=1+σg/jωε01+jω/2Γ,where the volumetric DC conductivity of graphene is (see, e.g., [2])(2)σg=q2kBTπħ2Γdln⁡2cosh⁡μc2kBT,with q being the electron charge, ħ the reduced Planck constant, kB the Boltzmann constant, T the temperature, Γ the scattering rate, μc the chemical potential, and d 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 J. Thus, we are looking for the solution of(3)μ0∂H∂t=-∇×Eε0∂E∂t+J=∇×H12Γ∂J∂t+J=σgE,with μ0 and ε0 being the vacuum permeability and permittivity, respectively. They constitute the vacuum wave impedance and phase velocity(4)Z0=μ0ε0c0=1μ0ε0.The substitution(5)E=Z0E^τ=c0tσ^g=Z0σgpermits symmetrizing Maxwell’s equations, yielding(6)∂H∂τ=-∇×E^∂E^∂τ+J=∇×Hc02Γ∂J∂τ+J=σ^gE^.In 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(7)κyzx∂H¯∂τ+σ^yzxH¯=-∇×E^κzxy∂H∂τ+σ^zxyH=κxyz∂H¯∂τ+σ^xyzH¯∂E̿∂τ+J=∇×Hc02Γ∂J∂τ+J=σ^gE¯¯κyzx∂E¯∂τ+σ^yzxE¯=∂E̿∂τκzxy∂E^∂τ+σ^zxyE^=κxyz∂E¯∂τ+σ^xyzE¯, where a short notation was adopted for the tensors; for example,(8)κyzx=κyκzκxσ^yzx=Z0σyσzσx, with, for the ith PML layer,(9)κu=1+κumax-1inpmlmσu=σumaxinpmlmu∈x,y,z.This is the so-called polynomially graded PML with number of layers npml and power m. In this paper, we always choose the standard values npml=10 and m=4. Also, from [3], we borrow the optimal value(10)σumax=m+1150πΔuu∈x,y,z,where Δu is the spatial step in the u-dimension (i.e., Δx, Δy, or Δz). 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 κu aims to absorb the evanescent waves, whereas the imaginary-stretch parameter σu 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 u-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 Figure1. 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 x-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, Ex and Hx, 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:(11)γy+h¯xi,j+1/2,k+1/2n+1=γy-h¯xi,j+1/2,k+1/2n-1Δye^zi,j+1,k+1/2n+1/2-e^zi,j,k+1/2n+1/2+1Δze^yi,j+1/2,k+1n+1/2-e^yi,j+1/2,kn+1/2γz+hxi,j+1/2,k+1/2n+1=γz-hxi,j+1/2,k+1/2n+γx+h¯xi,j+1/2,k+1/2n+1-γx-h¯xi,j+1/2,k+1/2nfor the magnetic field and(12)1Δτ+σ^g2α+e̿xi,j,kn+1=1Δτ-σ^g2α+e̿xi,j,kn-121+α-α+jxi,j,kn+1Δyhzi,j+1/2,kn+1/2-hzi,j-1/2,kn+1/2-1Δzhyi,j,k+1/2n+1/2-hyi,j,k-1/2n+1/2α+jxi,j,kn+1=α-jxi,j,kn+σ^ge̿xi,j,kn+1+e̿xi,j,knγy+e¯xi,j,kn+1=γy-e¯xi,j,kn+1Δτe̿xi,j,kn+1-e̿xi,j,knγz+e^xi,j,kn+1=γz-e^xi,j,kn+γx+e¯xi,j,kn+1-γx-e¯xi,j,knfor the electric field. Here, we used the notation(13)α±=c0ΓΔτ±1γu±=κuΔτ±Z0σu2u∈x,y,z.Figure 1
The UCHIE cell for implicitization and collocation in thex-dimension. Triangular arrows denote electric field components and square arrows magnetic field components. Blue and red markers are staggered in time.Figure 2
Visualization of the difference method applied to the Maxwell equationε  ∂tEy+Jy=∂zHx-∂xHz. The black dots represent the field samples in the x,t-space, the red arrows are differences, and the green arrows are interpolations. All terms are evaluated with second-order accuracy in the midpoint of the rectangle formed by the four sample points.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 nx cells in the x-dimension, there are nx-1 electric field and nx+1 magnetic field samples in (14)-(15). The bracket notation is used to denote diagonal matrices whose rank is nx-1 or nx+1 depending on whether they act on electric or magnetic field components, respectively.(14)α+-σ^g0000A12A1Δτ00D2000A2γx+00D10-IΔτ0γz+0000-γz+0γy+0000-γy+0γx+jye̿yh¯ze¯yhze^yj+1/2,kn+1/2=α-σ^g0000-A12A1Δτ00-D2000A2γx-00-D10-IΔτ0γz-0000-γz-0γy-0000-γy-0γx-jye̿yh¯ze¯yhze^yj+1/2,kn-1/2+0A1Δzhxj+1/2,k+1/2n-hxj+1/2,k-1/2nA2Δye^xj+1,kn-e^xj,kn000(15)α+-σ^g0000A12A1Δτ00-D2000A2γz+00-D10-IΔτ0γx+0000-γy+0γx+0000-γz+0γy+jze¯¯zh¯ye¯zhye^zj,k+1/2n+1/2=α-σ^g0000-A12A1Δτ00D2000A2γz-00D10-IΔτ0γx-0000-γy-0γx-0000-γz-0γy-jze¯¯zh¯ye¯zhye^zj,k+1/2n-1/2-0A1Δyhxj+1/2,k+1/2n-hxj-1/2,k+1/2nA2Δze^xj,k+1n-e^xj,kn000.Besides, we introduced the x-interpolators(16)A1=111⋱111nx×nx-1A2=11⋱11nx×nx+1,and the x-differentiators(17)D1=Δx1⋱Δxnx-11-11⋱-11-1nx×nx-1D2=Δx1⋱Δxnx-1-11⋱-11nx×nx+1.As 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 M from now on, is partitioned as follows:(18)M11M12M21M22=α+-σ^g0000A12A1Δτ00D2000A2γx+00D10-IIΔτ0γz+0000-γz+0γy+0000-γy+0γx+.Its inverse is given by [4, prop. 2.8.7](19)M-1=S-1-S-1M12M22-1-M22-1M21S-1M22-1+M22-1M21S-1M12M22-1,with(20)M22-1=γz+-1000γy+-10γx+-1γy+γz+-10γx+-1,and Schur complement(21)S=M11-M12M22-1M21,which is, in turn, partitioned as follows:(22)S11S12S21S22=α+-σ^g0A1I2A1IΔτD2γy+-1γz+0D1Δτγx+-1γy+γz+-1A2γx+.Analogously, the inverse of the Schur complement is given by(23)S-1=S11-1+S11-1S12Ss-1S21S11-1-S11-1S12Ss-1-Ss-1S21S11-1Ss-1,such that, in the end, solving the large sparse system (14) only requires the inversion of the rank-2nx matrix(24)Ss=S22-S21S11-1S12=A11Δτ+12α+-1σ^gD2γy+-1γz+D1Δτγx+-1γy+γz+-1A2γx+.The reordering matrix(25)Rij=1i=1,j=nx1i=2k,j=k1i=2k+1,j=nx+k1i=2nx,j=2nx0elsek=1,…,nx-1together with the rank-2nx primary circulant matrix C, that is, the circulant matrix whose first row has second element one and all other elements zero, allows constructing a banded matrix(26)B=CRSsRT,which has four elements per row arranged in a staircase pattern. For example, for nx=5, its sparsity pattern is(27)××××××××××××××××××××××××××××××××××××.The LU factorization with row and column pivoting of B, that is,(28)PBQ=LU,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 v1 and v2 conforming to the partitioning applied in (18).(2)
Computep=v1-M12M22-1v2 and split this vector into two parts p1 and p2 conforming to the partitioning applied in (22).(3)
Use forward and backward triangular sweeps to computee̿y,h¯z=RTQU-1L-1PCRp2-S21S11-1p1.(4)
Computejy=S11-1p1-S12e̿y,h¯z.(5)
Computee¯y,hz,e^y=M22-1v2-M21jy,e̿y,h¯z.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 2nx instead of two tridiagonal matrices of rank nx. 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 hasd=1nm, T=300K, μc=0.05eV, and Γ=0.5THz. Hence, according to (2), the sheet has a DC conductivity σg=6.7S/μm. 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 Ezgt is recorded in one point at the other side of the sheet. The source is a differentiated Gaussian pulse with bandwidth b=0.55THz, temporal width tw=2/πb, and delay td=6tw. The spatial invariance in the y- and z-dimension is fulfilled by imposing periodic boundary conditions (PBCs). Both the transmitted and the back-scattered plane waves are absorbed by UPMLs with κxmax=1. The overall grid is uniformly discretized with steps Δx=Δy=Δz≈27μm, whereas the graphene sheet is locally resolved by two small cells of size d along the x-axis as depicted in Figure 4. The time step equals the Courant limit in the yz-plane, more specifically(29)Δτ=11/Δy2+1/Δz2≈19μm,which is about 19,000 times larger than the time step that would have been used by the classical FDTD method. The recorder is placed 136μm 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 Ez0t. The time-domain data are then Fourier transformed and their ratio determines the shielding effectiveness(30)SEnumω=20log10⁡FEz0tFEzgt.The analytical solution for the SE is known to be (31)SEanaω=20log10⁡1-r2e-jkg2d1-r2ejkgd,with r being the Fresnel reflection coefficient of a single vacuum-graphene interface(32)r=1-εrg1+εrg, and wave number(33)kg=ωc0εrg.Recall that the effective relative permittivity of the graphene sheet εrg 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 0.5THz, the skin depth, hereby meaning the thickness that the graphene sheet should have to reduce the amplitude of a plane wave by 63%, is 360nm and this value increases towards DC.Figure 3
The TFSF plane(blue) excites a plane wave with nonzero components Hy and Ez which normally impinges upon a thin graphene sheet(hexagons) of thickness d=1nm. The electric field Ez is recorded a few cells behind the graphene sheet(red). Infiniteness in the y- and z-dimension is mimicked by PBCs along four of the exterior faces. The remaining two faces are covered by PMLs to absorb outgoing waves in the x-dimension.Figure 4
Discretization along thex-axis. The vacuum part of the grid is uniformly discretized with step Δx≈27μm. The graphene sheet is approximated by two small UCHIE cells of size d=1nm. The graphene conductivity σg=6.7S/μm is solely assigned to the single discretization point enclosed by these two small cells and does not need to be averaged as is the case for the conventional Yee grid.Figure 5
Shielding effectiveness of the1nm graphene sheet illuminated by a normally incident plane wave.


### 5.2. SE for anE- and H-Dipole

A similar simulation is repeated but this time a dipole of lengthΔz is excited, at a distance ds=409μm from the left side of the graphene sheet. All six exterior faces of the grid are covered by UPMLs, this time with κxmax=κymax=κzmax=10 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 E-dipole, the SE is again determined by recording Ez and taking the ratio defined in (30), whereas, for the H-dipole, the dual definition of the SE is adopted where Ez is replaced by Hz. The recorded fields are plotted in Figure 6 together with the plane-wave data from Section 5.1. The resulting SEs for both the E- and the H-dipole are shown in Figure 7. Below 0.12THz, which corresponds to k0ds<1, the graphene sheet experiences the reactive near field of the dipole. For this frequency range, the wave impedances of the E- and H-dipole are approximately ZE=Z0/k0ds>Z0 and ZH=Z0k0ds<Z0, respectively. The graphene sheet behaves like a classical good conductor at low frequencies. As such, it has a very low wave impedance. Consequently, the E-dipole experiences a higher contrast resulting in more reflections, whereas the H-dipole does the opposite. This is exactly what we observe in Figure 7: at low frequencies, the E-dipole is better shielded by the graphene sheet due to strong reflections, whereas the H-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.Figure 6
Fields recorded in the observation point behind the graphene sheet. Only a small part of the total simulated time (0.6428ns) is shown.Figure 7
Shielding effectiveness of the1nm graphene sheet illuminated by a small E- and H-dipole.


## 5.1. SE for a Plane Wave

The graphene sheet under investigation hasd=1nm, T=300K, μc=0.05eV, and Γ=0.5THz. Hence, according to (2), the sheet has a DC conductivity σg=6.7S/μm. 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 Ezgt is recorded in one point at the other side of the sheet. The source is a differentiated Gaussian pulse with bandwidth b=0.55THz, temporal width tw=2/πb, and delay td=6tw. The spatial invariance in the y- and z-dimension is fulfilled by imposing periodic boundary conditions (PBCs). Both the transmitted and the back-scattered plane waves are absorbed by UPMLs with κxmax=1. The overall grid is uniformly discretized with steps Δx=Δy=Δz≈27μm, whereas the graphene sheet is locally resolved by two small cells of size d along the x-axis as depicted in Figure 4. The time step equals the Courant limit in the yz-plane, more specifically(29)Δτ=11/Δy2+1/Δz2≈19μm,which is about 19,000 times larger than the time step that would have been used by the classical FDTD method. The recorder is placed 136μm 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 Ez0t. The time-domain data are then Fourier transformed and their ratio determines the shielding effectiveness(30)SEnumω=20log10⁡FEz0tFEzgt.The analytical solution for the SE is known to be (31)SEanaω=20log10⁡1-r2e-jkg2d1-r2ejkgd,with r being the Fresnel reflection coefficient of a single vacuum-graphene interface(32)r=1-εrg1+εrg, and wave number(33)kg=ωc0εrg.Recall that the effective relative permittivity of the graphene sheet εrg 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 0.5THz, the skin depth, hereby meaning the thickness that the graphene sheet should have to reduce the amplitude of a plane wave by 63%, is 360nm and this value increases towards DC.Figure 3
The TFSF plane(blue) excites a plane wave with nonzero components Hy and Ez which normally impinges upon a thin graphene sheet(hexagons) of thickness d=1nm. The electric field Ez is recorded a few cells behind the graphene sheet(red). Infiniteness in the y- and z-dimension is mimicked by PBCs along four of the exterior faces. The remaining two faces are covered by PMLs to absorb outgoing waves in the x-dimension.Figure 4
Discretization along thex-axis. The vacuum part of the grid is uniformly discretized with step Δx≈27μm. The graphene sheet is approximated by two small UCHIE cells of size d=1nm. The graphene conductivity σg=6.7S/μm is solely assigned to the single discretization point enclosed by these two small cells and does not need to be averaged as is the case for the conventional Yee grid.Figure 5
Shielding effectiveness of the1nm graphene sheet illuminated by a normally incident plane wave.


## 5.2. SE for anE- and H-Dipole

A similar simulation is repeated but this time a dipole of lengthΔz is excited, at a distance ds=409μm from the left side of the graphene sheet. All six exterior faces of the grid are covered by UPMLs, this time with κxmax=κymax=κzmax=10 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 E-dipole, the SE is again determined by recording Ez and taking the ratio defined in (30), whereas, for the H-dipole, the dual definition of the SE is adopted where Ez is replaced by Hz. The recorded fields are plotted in Figure 6 together with the plane-wave data from Section 5.1. The resulting SEs for both the E- and the H-dipole are shown in Figure 7. Below 0.12THz, which corresponds to k0ds<1, the graphene sheet experiences the reactive near field of the dipole. For this frequency range, the wave impedances of the E- and H-dipole are approximately ZE=Z0/k0ds>Z0 and ZH=Z0k0ds<Z0, respectively. The graphene sheet behaves like a classical good conductor at low frequencies. As such, it has a very low wave impedance. Consequently, the E-dipole experiences a higher contrast resulting in more reflections, whereas the H-dipole does the opposite. This is exactly what we observe in Figure 7: at low frequencies, the E-dipole is better shielded by the graphene sheet due to strong reflections, whereas the H-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.Figure 6
Fields recorded in the observation point behind the graphene sheet. Only a small part of the total simulated time (0.6428ns) is shown.Figure 7
Shielding effectiveness of the1nm graphene sheet illuminated by a small E- and H-dipole.


## 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 intrabandand 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].



---

*Source: 5860854-2017-07-03.xml*