Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Accurate Time-Domain Modeling of Arbitrarily Shaped Graphene Layers Utilizing Unstructured Triangular Grids

Accurate Time-Domain Modeling of Arbitrarily Shaped Graphene Layers Utilizing Unstructured... axioms Article Accurate Time-Domain Modeling of Arbitrarily Shaped Graphene Layers Utilizing Unstructured Triangular Grids 1, 2 1 Stamatios Amanatiadis * , Theodoros Zygiridis and Nikolaos Kantartzis Department of Electrical & Computer Engineering, Aristotle University of Thessaloniki, GR-54124 Thessaloniki, Greece; kant@auth.gr Department of Electrical & Computer Engineering, University of Western Macedonia, GR-50100 Kozani, Greece; tzygiridis@uowm.gr * Correspondence: samanati@auth.gr Abstract: The accurate modeling of curved graphene layers for time-domain electromagnetic simula- tions is discussed in the present work. Initially, the advanced properties of graphene are presented, focusing on the propagation of strongly confined surface plasmon polariton waves at the far-infrared regime. Then, the implementation of an unstructured triangular grid was examined, based on the De- launay triangulation method. The electric-field components were placed at the edges of the triangles, while two different techniques were proposed for the sampling of the magnetic ones. Specifically, the first one suggests that the magnetic component is placed at the triangle’s circumcenter providing more accurate results, although instability may occur for nonacute triangles. On the other hand, the magnetic field was sampled at the triangle’s centroid, considering the second technique, ensuring the algorithm’s stability, but further approximations were required, leading to a slight accuracy reduction. Moreover, the updating equations in the time-domain were extracted via an appropriate approximation of Maxwell equations in their integral form. Finally, graphene was introduced in the computational domain as an equivalent surface current density, whose location matches the Citation: Amanatiadis, S.; corresponding electric components. The validity of our methodology was successfully performed via Zygiridis, T.; Kantartzis, N. Accurate the comparison of graphene surface wave propagation properties to their theoretical values, whereas Time-Domain Modeling of the global error determination indicates the minimal triangle dimensions. Additionally, an instructive Arbitrarily Shaped Graphene Layers setup comprising a circular graphene scatterer was analyzed thoroughly, to reveal our technique’s Utilizing Unstructured Triangular advantages compared to the conventional staircase discretization. Grids. Axioms 2022, 11, 44. https:// doi.org/10.3390/axioms11020044 Keywords: auxiliary differential equation; curved modeling; equilateral triangle; FDTD; finite inte- Academic Editor: Nhon gration; surface current density; surface wave Nguyen-Thanh Received: 7 December 2021 Accepted: 19 January 2022 1. Introduction Published: 22 January 2022 In previous years, graphene attracted rapid scientific attention due to its extraordinary Publisher’s Note: MDPI stays neutral properties in a broad field of applications. Particularly, this naturally two-dimensional with regard to jurisdictional claims in carbon allotrope presents significant and useful characteristics, in spite of its negligible published maps and institutional affil- thickness [1]. From an electromagnetic point of view, a finite surface conductivity is ob- iations. served, thus enabling several exotic phenomena, such as gyrotropy and the propagation of strongly confined surface plasmon polariton (SPP) waves at the far-infrared frequen- cies [2–4]. The latter paved the way for the design of a wide class of plasmonic devices at Copyright: © 2022 by the authors. the THz regime [5–8]. Licensee MDPI, Basel, Switzerland. The majority of these devices require efficient simulation techniques for their pre- This article is an open access article cise characterization. Therefore, the development of suitable numerical algorithms that distributed under the terms and model graphene accurately has been triggered. A very popular numerical scheme is the conditions of the Creative Commons finite-difference time-domain (FDTD) technique, which facilitates broadband simulations Attribution (CC BY) license (https:// of complex structures [9] owing to its explicit time-domain nature. For this reason, sev- creativecommons.org/licenses/by/ eral algorithms have been reported that modify the conventional FDTD method properly, 4.0/). Axioms 2022, 11, 44. https://doi.org/10.3390/axioms11020044 https://www.mdpi.com/journal/axioms Axioms 2022, 11, 44 2 of 11 while open-source implementations are available [10]. In general, graphene is treated as a surface boundary, and its inherent frequency dispersion can be modeled through an auxiliary differential equation [11,12] or by exploiting the magnetic-field boundary condi- tions [13]. Additionally, more advanced properties have been implemented, such as the case of magnetically biased graphene [14–16], its third-order response [17,18], and its complex conductivity at optical frequencies [19,20]. Despite its apparent flexibility, the main draw- back of the original FDTD algorithm is the discretization of the computational domain with a structured orthogonal mesh. On the other hand, it should be kept in mind that various de- vices feature optimized designs exploiting curved surfaces. As a consequence, an extremely fine mesh is required for the staircase approximation of such configurations [21–23], which degrades the simulation’s performance. Alternatively, frequency-domain numerical solvers can be employed [24–27], sacrificing the advantages of time-domain solutions. In this paper, we propose the accurate numerical modeling of arbitrarily shaped graphene layers in the time domain via the discretization of the computational space with unstructured triangular meshes. The analysis starts with the presentation of graphene main attributes, such as its surface conductivity and the fundamental propagation properties of the supported SPP waves. Furthermore, the basic aspects of the unstructured triangular grids were investigated for time-domain electromagnetic simulations. Specifically, two different sampling techniques for the magnetic components were examined, at the triangle circumcenter and at its centroid, while the electric components were placed at the edges of the triangles. Then, the graphene contribution is introduced as an equivalent surface current density, while a recursive convolution method (RCM) was developed for the consistent modeling of its frequency dispersion. The proposed algorithm was successfully validated through the direct comparison of the numerically extracted SPP wave propagation properties with the analytical values, while the global error calculation indicated certain discretization limitations. Finally, an interesting configuration that includes a circular graphene scatterer was designed, and the comparison with the conventional staircase approximation highlighted the superiority of the unstructured grids, in case of problems with nonorthogonal features. 2. Theoretical Formulation 2.1. Graphene Surface Conductivity and Surface Wave Propagation Properties Throughout this work, graphene is considered to be a naturally infinitesimally thin material, characterized via its surface conductivity s (w, m , G, T). The latter depends on gr c the radial frequency w, the chemical potential m that can be controlled either via chemical doping or an electrostatic bias application, the scattering rate G that constitutes the main loss mechanism of the material, and the temperature T. Graphene’s surface conductivity is evaluated utilizing the Kubo formula [28], leading to the following compact expression in the absence of an external magnetostatic field [29]: e (jw + 2G) 1 ¶ f (#) ¶ f (#) d d s (w, m , G, T) =  # d# gr c 2 2 (jw + 2G) ¶# ¶# ph ¯ f (#) f (#) d d + d# , (1) 2 2 (jw + 2G) + 4(#/h ¯ ) where e is a single electron charge, h ¯ the reduced Planck’s constant and f (#) is the Fermi– Dirac distribution f (#) = , (2) (#m )/k T (e + 1) with k the Boltzmann’s constant. It is evident that (1) consists of two terms that are connected to the major electron transitions in graphene, namely the intraband and the interband. The former dominates at lower frequencies, approximately until the far-infrared Axioms 2022, 11, 44 3 of 11 regime, while the latter becomes considerable beyond the mid-infrared frequencies. Al- though the proposed methodology can be extended to incorporate more complex graphene properties, we focus on the lower regime, where the interband term is negligible. Conse- quently, the graphene’s surface conductivity expression can be significantly simplified to a Debye dispersion relation: e k T m m B m /k T c s (w, m , G, T) = + 2 ln(e + 1) = , (3) gr,intra c k T jw + 2G ph ¯ (jw + 2G) B where A is the frequency-independent term. It is important to mention that isotropic graphene conductivity is considered instead of the anisotropic one since our analysis is performed for a two-dimensional computational domain. The capability of graphene to support strongly confined SPP waves at the far-infrared regime is well established. The propagation properties of such waves depend on graphene’s surface conductivity, and their characterization can be performed via the propagation constant k . The latter has been evaluated analytically for a free-standing graphene layer of infinite dimensions by means of [29] k = k 1 , (4) s h gr 0 where k and h are the free-space wavenumber and wave impedance, respectively. Now, 0 0 the main propagation attributes, such as the SPP wavelength l , the propagation length SPP L and the confinement z, were straightforwardly extracted from the propagation con- SPP stant k via 2p 1 1 nq o l = , L = , z = , (5) SPP SPP <fk g 2=fk g 2 r r 2 < k k with <fg and =fg denoting the real and the imaginary part of the argument, respectively. 2.2. Unstructured Triangular Grids for Electromagnetic Analysis in the Time Domain The FDTD approach constitutes one of the most popular methods for electromagnetic simulations in the time domain [30]. The conventional algorithm solves the differential form of Maxwell equations by applying a finite-difference scheme on a structured rectangular grid. Despite the algorithm’s powerful characteristics, the use of a rectangular grid enforces staircase approximations to curved interfaces, which may lead to inaccurate results. For this reason, unstructured triangular grids have been proposed to effectively model such surfaces [31,32]. In our approach, the integral form of Maxwell equations is utilized for isotropic media I ZZ ¶H E dl = m m  dA, (6) ¶t C A I ZZ ¶E H dl = # # + sE  dA, (7) 0 r ¶t C A where A is the area limited by the closed-loop C with length l, # and m are the electric 0 0 permittivity and magnetic permeability of free-space, respectively, and # , m are their r r relative equivalents. Our work focuses on a transverse electric two-dimensional (2D) case. The computa- tional domain is discretized into triangular elements using an efficient Delaunay triangula- tion algorithm, as depicted in Figure 1. Here, the electric-field components are tangential to the 2D domain and located along the edges of the triangles. It should be mentioned that the sampling of the magnetic components, which are considered normal to the domain, is not a trivial issue. In particular, two different approaches, depicted in Figure 1, are investigated. Axioms 2022, 11, 44 4 of 11 The initial one (Figure 1a) places the magnetic components at the triangle circumcenter, namely the center of the circle that is defined by its vertices creating thus, a Voronoi diagram [31]. The advantage of this approach is that the magnetic-field connection between two adjacent triangles is perpendicular to their common edge. This is important, since the inner product at the right-hand side of (7) suggests that the electric component is vertical to the virtual area of the connecting magnetic fields. Unfortunately, these points may lie outside the element for an obtuse triangle case, potentially influencing the algorithm’s stability. To circumvent this hindrance, a second technique can be used (Figure 1b), where the magnetic components are located at the triangle centroid, namely the averaged coordinates of its three vertices [32]. Although this approach ensures that each magnetic component lies inside the corresponding element, the magnetic-field connection between adjacent triangles is not necessarily vertical to the common electric component. Hence, an angular correction is required for the accurate implementation of (7). Note that the two techniques are identical for the case of equilateral triangles. In this manner, an optimal triangulation can be considered to consist of such elements. k k l l E E i,k H i,k i,l E ,J E E ,J i,j gr[i,j] i,l i,j gr[i,j] graphene graphene j j (a) Circumcenter sampling (b) Centroid sampling Figure 1. Topology of the unstructured triangular grid and the electromagnetic field assignment concerning the two magnetic field sampling techniques. In particular, the arrangement of the electric components is at the edges of the triangles, while the magnetic components lie inside the triangle (a) at its circumcenter or (b) at its centroid. Now, a finite integration was applied to (6) and (7). First of all, the temporal discretiza- tion of fields was defined to coincide with the conventional FDTD algorithm, i.e., E was updated at nDt and H at (n + 1/2)Dt instants, where Dt is the time-step. Then, considering element i, (6) was developed into n+1/2 n1/2 H H n n n i i E l + E l + E l = m m A, 0 r i,j i,k i,l i,j i,k i,l Dt Dt n+1/2 n1/2 n n n H = H (E l + E l + E l ). (8) i,j i,k i,l i,j i,k i,l i i A m m i 0 r The finite integration of (7) depends on the applied magnetic field component sam- pling technique. The circumcenter-sampling case provides the following straightforward update equation for E i,j n1 n1 n n E E E + E i,j i,j i,j i,j n1/2 n1/2 H H = # # + s d 0 r i,j i j Dt 2 2# # sDt 2Dt 1 0 r n n1 n1/2 n1/2 E = E + (H H ) . (9) i,j i,j i j 2# # + sDt 2# # + sDt d 0 r 0 r i,j i,j i,j i,l i,l Axioms 2022, 11, 44 5 of 11 For the centroid-sampling case, an angular correction must be included to model the finite integration accurately. Explicitly, 2# # sDt 2Dt 1 0 r n n1 n1/2 n1/2 E = E + (H H ) . (10) i,j i,j i j 2# # + sDt 2# # + sDt d sin q r r 0 0 i,j The electromagnetic properties (described by # , m , s) of each element are assumed r r to be constant within each element. Consequently, the tangential continuity of the electric component is satisfied a priori. Each set of d and l is shared between two adjacent i,j i,j elements. Additionally, E is shared between triangles i and j, and a unique definition of i,j its orientation must be selected. Herein, the positive direction is counterclockwise with respect to the element with a lower index (i in Figure 1). Thus, a negative contribution is assumed for the element with a larger index (j in Figure 1). Finally, a stability criterion should be established, similar to the conventional FDTD algorithm. In essence, the time increment should satisfy the relation [31] 0 1 B C 1 2A B C u i Dt 6 minB C, (11) @ t l A i,fg i,fg where A is the area of the element i, l the length of the common edge with triangle i,fg fg, d the magnetic-field distance from the corresponding triangle (Figure 1), and c the i,fg free-space speed of light. Particularly, the inner expression of (11) was calculated for all i,fg triangles, where the ratio corresponds to the cell’s dimensions since it is connected i,fg to all three triangle edges. Then, the stable time increment was defined via the mini- mum value of the aforementioned calculations. It is worth mentioning that this criterion converges to the popular Courant–Friedrichs–Lewy condition of the conventional FDTD structured discretization. 2.3. Graphene Modeling within Unstructured Triangular Grids Graphene was embedded in the computational domain in terms of an equivalent surface current J = s E, as demonstrated in Figure 1. Observe that the locations of gr gr the surface current components are identical to the corresponding electric ones. The contribution of this surface current is included into the computational domain via Ampere’s law, namely (7) via I ZZ ZZ ¶E H dl = # # + sE  dA + J  dA. (12) 0 r gr ¶t C A A Consequently, the electric-field updating Equation (9), for the circumcenter sampling case, was properly modified to 2# # sDt 2Dt 1 0 r n n1 n1/2 n1/2 E = E + (H H ) i,j i,j i j 2# # + sDt 2# # + sDt d 0 r 0 r i,j 2Dt 1 n1/2 J . (13) gr[i,j] 2# # + sDt d 0 r i,j One may observe the d parameter in the denominator, which is attributed to the i,j fact that graphene is modeled as a surface current density, instead of a volume density. Therefore, the infinitesimally thin nature of the material was retained. Similarly, (10), for the centroid sampling case, was modified to Axioms 2022, 11, 44 6 of 11 2# # sDt 2Dt 1 0 r n n1 n1/2 n1/2 E = E + (H H ) i,j i,j i j 2# # + sDt 2# # + sDt d sin q 0 r 0 r i,j 2Dt 1 n1/2 J . (14) gr[i,j] 2# # + sDt d sin q 0 r i,j Now, the main concern is the modeling of graphene’s frequency dispersion. Firstly, a transition of J = s E and the surface conductivity (3) into the time domain should be gr gr performed using an inverse Fourier transform (IFT) m IFT 2Gt J (w) = s (w)E (w) = E (w) ! J (t) = A e  E (t). (15) gr i,j i,j m i,j gr[i,j] gr[i,j] c jw + 2G Although a procedure that involves a convolution should be performed for the surface current updating, the efficient RCM concept is utilized, i.e., n n o n+1/2 2G(nm)Dt m J = A e E Dt , å c i,j gr[i,j] m=0 n+1/2 2GDt n1/2 n J = A e J + A DtE . (16) m m c c i,j gr[i,j] gr[i,j] Interestingly, there is no connection between the frequency-dispersion modeling and the computational-domain discretization with unstructured triangular meshes. Therefore, more advanced schemes, such as interband conductivity approximation, third-order re- sponse, or anisotropic conductivity due to a magnetostatic-bias application, for 3D domains, can be straightforwardly combined with the triangular grids. Additionally, other complex dielectrics are easily importable via equivalent frequency-dispersion schemes, such as the auxiliary differential equation of conventional FDTD, since they are, also, independent of the computational-domain discretization. The entire procedure of the proposed algorithm is summarized in the block diagram of Figure 2. Here, every individual step is included, such as the initialization and the electromagnetic component updating via the well-known leap-frog scheme. Initialization Simulation n n+1/2 H locations H updating eq. (7) Mesh E locations E updating generation eq. (12) or (13) J updating gr eq. (15) Element size Figure 2. Block diagram of the proposed scheme. 3. Proposed Algorithm Validation The validity of the featured algorithm is pursued through the well-established setup of surface wave propagation onto a straight free-standing graphene layer of infinite size. The advantage of this example is its direct comparison to the analytical expressions (4) and (5). Moreover, the conventional FDTD method is able to provide accurate results for this orthog- onal setup, in order to calculate reliably the global error within computational domain. To this aim, graphene chemical potential was selected m = 0.2 eV, and its scattering rate was set to G = 0.11 meV, while the analysis was performed at the far-infrared spectrum, namely from 1 to 3 THz. The dimension of the computational domain is 600 m  150 m, with the graphene layer located along the x-axis. As discussed, an unstructured Delaunay trian- gulation was utilized with a maximum edge size D = 3 m, resulting in approximately max 32,000 elements. The time-step was selected 3.7 fs to ensure the algorithm’s stability, while open boundaries were terminated via an 8-cell thick perfectly matched layer (PML). It is Axioms 2022, 11, 44 7 of 11 important to mention that PML is selected instead of an Absorbing Boundary Condition since it is possible to extend graphene inside the PML region. This leads to the effective ab- sorption of the strongly confined SPP waves. Finally, 8-cells are adequate for the considered applications due to the inherent losses of graphene surface waves. The numerically extracted propagation properties of the supported surface wave on graphene are compared to the aforementioned theoretical calculations in Figure 3. It is clear that the matching of the results is evident, especially until 2.5 THz. Beyond this frequency, the SPP wavelength is less than 30 m, and the conventional l/10 discretization limit is exceeded, leading to small discrepancies. It is important to mention that the simulated results represent both magnetic-field sampling techniques, since their differences is negligible. Furthermore, the distribution of the magnetic field at 2 THz is demonstrated in Figure 4, clearly indicating the propagation of the strongly confined surface wave on the graphene surface. Consequently, the proposed integration of the infinitesimally thin material into an unstructured triangular grid was validated successfully. −3 SPP −4 SPP −5 Theoretical Simulation −6 1 1.5 2 2.5 3 Frequency [THz] Figure 3. Numerical validation of surface wave propagation characteristics for a graphene layer with m = 0.2 eV and G = 0.11 meV. −50 −200 −150 −100 −50 0 50 100 150 200 x-axis [mm] Figure 4. Distribution of the magnetic-field component for a graphene layer with m = 0.2 eV and G = 0.11 meV at 2 THz. Finally, the global error in the computational domain was investigated for different element sizes, to check the influence of the element size. For this reason, a numerical simulation through the conventional orthogonal grid of the FDTD method was conducted using a considerably fine mesh (Dx = Dy < l /20). The mean global error was calculated for the SPP magnetic-field component, while it was normalized to the mean magnetic-field value of the computational domain. The calculated results are depicted in Figure 5, and it is apparent that smaller triangle dimensions lead to more precise results. Note that the shape of the curves is, generally, altered with D since different triangular meshes are generated. In addition, the max magnetic-field sampling at the triangle circumcenter (technique 1) presents a slightly lower global error. In any case, the satisfaction of the well-known D < l /10 rule secures max SPP reliable numerical results. Length (mm) y-axis [mm] Axioms 2022, 11, 44 8 of 11 0.5 Technique 1 Technique 2 0.4 0.3 D =10 mm max D =5 mm 0.2 max 0.1 D =3 mm max 1 1.5 2 2.5 3 Frequency [THz] Figure 5. Global error of the magnetic field, normalized to the mean computational domain value, for different triangle size discretization of the two proposed techniques. 4. Incident Field toward a Circular Graphene Scatterer Although the previous example validated the accuracy of the proposed methodology, the computational efficiency is not optimal compared to the conventional FDTD algorithm. The main reason for the degraded performance is the requirement for the storage of ad- ditional parameters, such as triangle dimensions as well as pointers to identify adjacent elements. For this reason, a second example is discussed in this section to exploit the advantages of the triangular computational domain. Particularly, a plane wave propagates toward a periodic structure of circular graphene scatterers of radius R = 20 m and period- icity of w = 100 m, as illustrated in Figure 6. The graphene chemical potential is set to 0.2 eV, and its scattering rate is 0.11 meV at the frequency range 1–3 THz. The computational domain is divided into 23,000 triangular elements of maximum edge size D = 3 m, max while the time-step is selected 3.7 fs. The side boundaries satisfy periodic conditions, and the extracted quantity of interest is the configuration’s transmission coefficient. periodic boundaries graphene Figure 6. Setup of an incident plane wave toward a periodic structure of circular graphene scatterers. The extracted numerical results of the proposed algorithm are compared to the conven- tional staircase approximation of the FDTD algorithm, as well as the commercial software COMSOL Multiphysics commercial package [33] for reference. Figure 7 proves the excel- lent behavior of our methodology. On the other hand, the staircase approximation fails to track the resonance frequencies accurately, even when using a very fine mesh (D = 1 m). This performance is explained via the distribution of the electric field in the scatterer in Figure 8. Specifically, the edge effects of the staircase approximation degrade the plasmonic resonances on the graphene, in contrast to the proposed method, where the electric field is smoothly distributed on the material surface. Relative global error Axioms 2022, 11, 44 9 of 11 −5 −10 COMSOL −15 Proposed method Staircase FDTD (D=5mm) Staircase FDTD (D=1mm) −20 1 1.5 2 2.5 3 Frequency [THz] Figure 7. Transmission coefficient comparison between the proposed method and the conventional staircase FDTD algorithm for a circular graphene scatterer of 20 m radius. [V/m] 0.8 0.6 0.4 0.2 −0.2 −0.4 −0.6 −0.8 −1 (a) Staircase (b) Proposed Figure 8. Comparison of the electric field distribution between the proposed method and the conventional staircase FDTD algorithm for a circular graphene scatterer of 20 m radius at 1.72 THz. 5. Conclusions The accurate time-domain modeling of curved graphene layers was accomplished in the present work using an unstructured triangular mesh. Initially, the fundamental characteristics of the two-dimensional material were investigated, as well as the finite integration scheme for the proposed computational grid. Particularly, two techniques were examined that coincide for equilateral triangle discretization. The strength of the first one relies on its more physical representation to achieve increased accuracy, but acute triangles are required to ensure stability. This problem is resolved via the second technique, although the necessity of further approximations slightly decreases the algorithm’s accuracy. Then, graphene was introduced into the numerical algorithm as an equivalent surface current density, while its frequency dispersion was treated via an RCM technique. The proposed methodology was thoroughly examined in terms of the surface wave propagation properties, validating the proposed implementation. Finally, a circular graphene scatterer setup was designed to exploit the algorithm’s applicability and advantages, compared to the conventional staircase approximation. Transmission coefficient [dB] Axioms 2022, 11, 44 10 of 11 Author Contributions: Conceptualization, S.A.; methodology, S.A.; validation, S.A.; writing— original draft preparation, S.A.; writing—review and editing, T.Z. and N.K.; supervision, T.Z. and N.K. All authors have read and agreed to the published version of the manuscript. Funding: This research was cofinanced by Greece and the European Union (European Social Fund- ESF) through the Operational Programme «Human Resources Development, Education and Lifelong Learning» in the context of the project “Reinforcement of Postdoctoral Researchers—2nd Cycle” (MIS-5033021), implemented by the State Scholarships Foundation (IKY). Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: The data presented in this study are available via a straightforward reproduction of the provided simulation details. Conflicts of Interest: The authors declare no conflict of interest. Abbreviations The following abbreviations are used in this manuscript: FDTD Finite-Difference Time-Domain PML Perfectly Matched Layer RCM Recursive Convolution Method SPP Surface Plasmon Polariton References 1. Novoselov, K.; Geim, A.K.; Morozov, S.; Jiang, D.; Zhang, Y.; Dubonos, S.; Grigorieva, I.; Firsov, A. Electric field effect in atomically thin carbon films. Science 2004, 306, 666–669. [CrossRef] [PubMed] 2. Mikhailov, S.A.; Ziegler, K. New electromagnetic mode in graphene. Phys. Rev. Lett. 2007, 99, 016803. [CrossRef] [PubMed] 3. Grigorenko, A.N.; Polini, M.; Novoselov, K. Graphene plasmonics. Nat. Photonics 2012, 6, 749–758. [CrossRef] 4. Bludov, Y.V.; Ferreira, A.; Peres, N.M.; Vasilevskiy, M.I. A primer on surface plasmon-polaritons in graphene. Int. J. Mod. Phys. B 2013, 27, 1341001. [CrossRef] 5. Ansell, D.; Radko, I.; Han, Z.; Rodriguez, F.; Bozhevolnyi, S.; Grigorenko, A. Hybrid graphene plasmonic waveguide modulators. Nat. Commun. 2015, 6, 8846. [CrossRef] [PubMed] 6. Xiao, S.; Zhu, X.; Li, B.H.; Mortensen, N.A. Graphene-plasmon polaritons: From fundamental properties to potential applications. Front. Phys. 2016, 11, 117801. [CrossRef] 7. Li, Y.; Tantiwanichapan, K.; Swan, A.K.; Paiella, R. Graphene plasmonic devices for terahertz optoelectronics. Nanophotonics 2020, 9, 1901–1920. [CrossRef] 8. Zhang, J.; Hong, Q.; Zou, J.; He, Y.; Yuan, X.; Zhu, Z.; Qin, S. Fano-resonance in hybrid metal-graphene metamaterial and its application as mid-infrared plasmonic sensor. Micromachines 2020, 11, 268. [CrossRef] 9. Choroszucho, A. Analysis of the influence of the complex structure of clay hollow bricks on the values of electric field intensity by using the FDTD method. Arch. Electr. Eng. 2016, 65, 745–759. [CrossRef] 10. Oskooi, A.F.; Roundy, D.; Ibanescu, M.; Bermel, P.; Joannopoulos, J.D.; Johnson, S.G. MEEP: A flexible free-software package for electromagnetic simulations by the FDTD method. Comput. Phys. Commun. 2010, 181, 687–702. [CrossRef] 11. Bouzianas, G.D.; Kantartzis, N.V.; Yioultsis, T.V.; Tsiboukis, T.D. Consistent study of graphene structures through the direct incorporation of surface conductivity. IEEE Trans. Magn. 2014, 50, 161–164. [CrossRef] 12. Ramadan, O. Improved direct integration auxiliary differential equation FDTD scheme for modeling graphene drude dispersion. Optik 2020, 219, 165173. [CrossRef] 13. Nayyeri, V.; Soleimani, M.; Ramahi, O.M. Modeling graphene in the finite-difference time-domain method using a surface boundary condition. IEEE Trans. Antennas Propag. 2013, 61, 4176–4182. [CrossRef] 14. Wang, X.H.; Yin, W.Y.; Chen, Z. Matrix exponential FDTD modeling of magnetized graphene sheet. IEEE Antennas Wirel. Propag. Lett. 2013, 12, 1129–1132. [CrossRef] 15. Amanatiadis, S.A.; Kantartzis, N.V.; Ohtani, T.; Kanai, Y. Precise modeling of magnetically biased graphene through a recursive convolutional FDTD method. IEEE Trans. Magn. 2017, 54, 7201504. [CrossRef] Axioms 2022, 11, 44 11 of 11 16. Ramadan, O. Simplified FDTD Formulations for Magnetostatic Biased Graphene Simulations. IEEE Antennas Wirel. Propag. Lett. 2020, 19, 2290–2294. [CrossRef] 17. Mock, A. Modeling passive mode-locking via saturable absorption in graphene using the finite-difference time-domain method. IEEE J. Quantum Electron. 2017, 53, 1300210. [CrossRef] 18. Amanatiadis, S.A.; Ohtani, T.; Zygiridis, T.T.; Kanai, Y.; Kantartzis, N.V. Modeling the Third-Order Electrodynamic Response of Graphene via an Efficient Finite-Difference Time-Domain Scheme. IEEE Trans. Magn. 2019, 56, 6700704. [CrossRef] 19. Mock, A. Padé approximant spectral fit for FDTD simulation of graphene in the near infrared. Opt. Mater. Express 2012, 2, 771–781. [CrossRef] 20. Amanatiadis, S.; Zygiridis, T.; Ohtani, T.; Kanai, Y.; Kantartzis, N. A Consistent Scheme for the Precise FDTD Modeling of the Graphene Interband Contribution. IEEE Trans. Magn. 2021, 57, 1600104. [CrossRef] 21. Xiao, S.; Wang, T.; Liu, Y.; Xu, C.; Han, X.; Yan, X. Tunable light trapping and absorption enhancement with graphene ring arrays. Phys. Chem. Chem. Phys. 2016, 18, 26661–26669. [CrossRef] [PubMed] 22. Xu, H.; Zhang, Z.; Wang, S.; Liu, Y.; Zhang, J.; Chen, D.; Ouyang, J.; Yang, J. Tunable graphene-based plasmon-induced transparency based on edge mode in the mid-infrared region. Nanomaterials 2019, 9, 448. [CrossRef] [PubMed] 23. Liu, T.; Zhou, C.; Jiang, X.; Cheng, L.; Xu, C.; Xiao, S. Tunable light trapping and absorption enhancement with graphene-based complementary metasurfaces. Opt. Mater. Express 2019, 9, 1469–1478. [CrossRef] 24. Liu, P.; Cai, W.; Wang, L.; Zhang, X.; Xu, J. Tunable terahertz optical antennas based on graphene ring structures. Appl. Phys. Lett. 2012, 100, 153111. [CrossRef] 25. Ye-xin, S.; Jiu-sheng, L.; Le, Z. Graphene-integrated split-ring resonator terahertz modulator. Opt. Quantum Electron. 2017, 49, 350. [CrossRef] 26. Wu, L.; Liu, H.; Li, J.; Wang, S.; Qu, S.; Dong, L. A 130 GHz electro-optic ring modulator with double-layer graphene. Crystals 2017, 7, 65. [CrossRef] 27. Bao, Z.; Tang, Y.; Hu, Z.D.; Zhang, C.; Balmakou, A.; Khakhomov, S.; Semchenko, I.; Wang, J. Inversion method characterization of graphene-based coordination absorbers incorporating periodically patterned metal ring metasurfaces. Nanomaterials 2020, 10, 1102. [CrossRef] 28. Gusynin, V.; Sharapov, S.; Carbotte, J. Magneto-optical conductivity in graphene. J. Phys. Condens. Matter 2006, 19, 026222. [CrossRef] 29. Hanson, G.W. Dyadic Green’s functions and guided surface waves for a surface conductivity model of graphene. J. Appl. Phys. 2008, 103, 064302. [CrossRef] 30. Taflove, A.; Hagness, S.C.; Piket-May, M. Computational electromagnetics: The finite-difference time-domain method. In The Electrical Engineering Handbook; Elsevier Inc.: New York, NY, USA, 2005; Volume 3. 31. Lee, C.; McCartin, B.; Shin, R.; Kong, J.A. A triangular-grid finite-difference time-domain method for electromagnetic scattering problems. J. Electromagn. Waves Appl. 1994, 8, 449–470. [CrossRef] 32. Liu, Y.; Sarris, C.D.; Eleftheriades, G.V. Triangular-mesh-based FDTD analysis of two-dimensional plasmonic structures supporting backward waves at optical frequencies. J. Light. Technol. 2007, 25, 938–945. [CrossRef] 33. COMSOL Multiphysics v. 5.5, COMSOL AB. 2019. Available online: www.comsol.com (accessed on November 2021). http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Axioms Multidisciplinary Digital Publishing Institute

Accurate Time-Domain Modeling of Arbitrarily Shaped Graphene Layers Utilizing Unstructured Triangular Grids

Download PDF
11 pages

Loading next page...
 
/lp/multidisciplinary-digital-publishing-institute/accurate-time-domain-modeling-of-arbitrarily-shaped-graphene-layers-MptC5NYkR1
Publisher
Multidisciplinary Digital Publishing Institute
Copyright
© 1996-2022 MDPI (Basel, Switzerland) unless otherwise stated Disclaimer The statements, opinions and data contained in the journals are solely those of the individual authors and contributors and not of the publisher and the editor(s). MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Terms and Conditions Privacy Policy
ISSN
2075-1680
DOI
10.3390/axioms11020044
Publisher site
See Article on Publisher Site

Abstract

axioms Article Accurate Time-Domain Modeling of Arbitrarily Shaped Graphene Layers Utilizing Unstructured Triangular Grids 1, 2 1 Stamatios Amanatiadis * , Theodoros Zygiridis and Nikolaos Kantartzis Department of Electrical & Computer Engineering, Aristotle University of Thessaloniki, GR-54124 Thessaloniki, Greece; kant@auth.gr Department of Electrical & Computer Engineering, University of Western Macedonia, GR-50100 Kozani, Greece; tzygiridis@uowm.gr * Correspondence: samanati@auth.gr Abstract: The accurate modeling of curved graphene layers for time-domain electromagnetic simula- tions is discussed in the present work. Initially, the advanced properties of graphene are presented, focusing on the propagation of strongly confined surface plasmon polariton waves at the far-infrared regime. Then, the implementation of an unstructured triangular grid was examined, based on the De- launay triangulation method. The electric-field components were placed at the edges of the triangles, while two different techniques were proposed for the sampling of the magnetic ones. Specifically, the first one suggests that the magnetic component is placed at the triangle’s circumcenter providing more accurate results, although instability may occur for nonacute triangles. On the other hand, the magnetic field was sampled at the triangle’s centroid, considering the second technique, ensuring the algorithm’s stability, but further approximations were required, leading to a slight accuracy reduction. Moreover, the updating equations in the time-domain were extracted via an appropriate approximation of Maxwell equations in their integral form. Finally, graphene was introduced in the computational domain as an equivalent surface current density, whose location matches the Citation: Amanatiadis, S.; corresponding electric components. The validity of our methodology was successfully performed via Zygiridis, T.; Kantartzis, N. Accurate the comparison of graphene surface wave propagation properties to their theoretical values, whereas Time-Domain Modeling of the global error determination indicates the minimal triangle dimensions. Additionally, an instructive Arbitrarily Shaped Graphene Layers setup comprising a circular graphene scatterer was analyzed thoroughly, to reveal our technique’s Utilizing Unstructured Triangular advantages compared to the conventional staircase discretization. Grids. Axioms 2022, 11, 44. https:// doi.org/10.3390/axioms11020044 Keywords: auxiliary differential equation; curved modeling; equilateral triangle; FDTD; finite inte- Academic Editor: Nhon gration; surface current density; surface wave Nguyen-Thanh Received: 7 December 2021 Accepted: 19 January 2022 1. Introduction Published: 22 January 2022 In previous years, graphene attracted rapid scientific attention due to its extraordinary Publisher’s Note: MDPI stays neutral properties in a broad field of applications. Particularly, this naturally two-dimensional with regard to jurisdictional claims in carbon allotrope presents significant and useful characteristics, in spite of its negligible published maps and institutional affil- thickness [1]. From an electromagnetic point of view, a finite surface conductivity is ob- iations. served, thus enabling several exotic phenomena, such as gyrotropy and the propagation of strongly confined surface plasmon polariton (SPP) waves at the far-infrared frequen- cies [2–4]. The latter paved the way for the design of a wide class of plasmonic devices at Copyright: © 2022 by the authors. the THz regime [5–8]. Licensee MDPI, Basel, Switzerland. The majority of these devices require efficient simulation techniques for their pre- This article is an open access article cise characterization. Therefore, the development of suitable numerical algorithms that distributed under the terms and model graphene accurately has been triggered. A very popular numerical scheme is the conditions of the Creative Commons finite-difference time-domain (FDTD) technique, which facilitates broadband simulations Attribution (CC BY) license (https:// of complex structures [9] owing to its explicit time-domain nature. For this reason, sev- creativecommons.org/licenses/by/ eral algorithms have been reported that modify the conventional FDTD method properly, 4.0/). Axioms 2022, 11, 44. https://doi.org/10.3390/axioms11020044 https://www.mdpi.com/journal/axioms Axioms 2022, 11, 44 2 of 11 while open-source implementations are available [10]. In general, graphene is treated as a surface boundary, and its inherent frequency dispersion can be modeled through an auxiliary differential equation [11,12] or by exploiting the magnetic-field boundary condi- tions [13]. Additionally, more advanced properties have been implemented, such as the case of magnetically biased graphene [14–16], its third-order response [17,18], and its complex conductivity at optical frequencies [19,20]. Despite its apparent flexibility, the main draw- back of the original FDTD algorithm is the discretization of the computational domain with a structured orthogonal mesh. On the other hand, it should be kept in mind that various de- vices feature optimized designs exploiting curved surfaces. As a consequence, an extremely fine mesh is required for the staircase approximation of such configurations [21–23], which degrades the simulation’s performance. Alternatively, frequency-domain numerical solvers can be employed [24–27], sacrificing the advantages of time-domain solutions. In this paper, we propose the accurate numerical modeling of arbitrarily shaped graphene layers in the time domain via the discretization of the computational space with unstructured triangular meshes. The analysis starts with the presentation of graphene main attributes, such as its surface conductivity and the fundamental propagation properties of the supported SPP waves. Furthermore, the basic aspects of the unstructured triangular grids were investigated for time-domain electromagnetic simulations. Specifically, two different sampling techniques for the magnetic components were examined, at the triangle circumcenter and at its centroid, while the electric components were placed at the edges of the triangles. Then, the graphene contribution is introduced as an equivalent surface current density, while a recursive convolution method (RCM) was developed for the consistent modeling of its frequency dispersion. The proposed algorithm was successfully validated through the direct comparison of the numerically extracted SPP wave propagation properties with the analytical values, while the global error calculation indicated certain discretization limitations. Finally, an interesting configuration that includes a circular graphene scatterer was designed, and the comparison with the conventional staircase approximation highlighted the superiority of the unstructured grids, in case of problems with nonorthogonal features. 2. Theoretical Formulation 2.1. Graphene Surface Conductivity and Surface Wave Propagation Properties Throughout this work, graphene is considered to be a naturally infinitesimally thin material, characterized via its surface conductivity s (w, m , G, T). The latter depends on gr c the radial frequency w, the chemical potential m that can be controlled either via chemical doping or an electrostatic bias application, the scattering rate G that constitutes the main loss mechanism of the material, and the temperature T. Graphene’s surface conductivity is evaluated utilizing the Kubo formula [28], leading to the following compact expression in the absence of an external magnetostatic field [29]: e (jw + 2G) 1 ¶ f (#) ¶ f (#) d d s (w, m , G, T) =  # d# gr c 2 2 (jw + 2G) ¶# ¶# ph ¯ f (#) f (#) d d + d# , (1) 2 2 (jw + 2G) + 4(#/h ¯ ) where e is a single electron charge, h ¯ the reduced Planck’s constant and f (#) is the Fermi– Dirac distribution f (#) = , (2) (#m )/k T (e + 1) with k the Boltzmann’s constant. It is evident that (1) consists of two terms that are connected to the major electron transitions in graphene, namely the intraband and the interband. The former dominates at lower frequencies, approximately until the far-infrared Axioms 2022, 11, 44 3 of 11 regime, while the latter becomes considerable beyond the mid-infrared frequencies. Al- though the proposed methodology can be extended to incorporate more complex graphene properties, we focus on the lower regime, where the interband term is negligible. Conse- quently, the graphene’s surface conductivity expression can be significantly simplified to a Debye dispersion relation: e k T m m B m /k T c s (w, m , G, T) = + 2 ln(e + 1) = , (3) gr,intra c k T jw + 2G ph ¯ (jw + 2G) B where A is the frequency-independent term. It is important to mention that isotropic graphene conductivity is considered instead of the anisotropic one since our analysis is performed for a two-dimensional computational domain. The capability of graphene to support strongly confined SPP waves at the far-infrared regime is well established. The propagation properties of such waves depend on graphene’s surface conductivity, and their characterization can be performed via the propagation constant k . The latter has been evaluated analytically for a free-standing graphene layer of infinite dimensions by means of [29] k = k 1 , (4) s h gr 0 where k and h are the free-space wavenumber and wave impedance, respectively. Now, 0 0 the main propagation attributes, such as the SPP wavelength l , the propagation length SPP L and the confinement z, were straightforwardly extracted from the propagation con- SPP stant k via 2p 1 1 nq o l = , L = , z = , (5) SPP SPP <fk g 2=fk g 2 r r 2 < k k with <fg and =fg denoting the real and the imaginary part of the argument, respectively. 2.2. Unstructured Triangular Grids for Electromagnetic Analysis in the Time Domain The FDTD approach constitutes one of the most popular methods for electromagnetic simulations in the time domain [30]. The conventional algorithm solves the differential form of Maxwell equations by applying a finite-difference scheme on a structured rectangular grid. Despite the algorithm’s powerful characteristics, the use of a rectangular grid enforces staircase approximations to curved interfaces, which may lead to inaccurate results. For this reason, unstructured triangular grids have been proposed to effectively model such surfaces [31,32]. In our approach, the integral form of Maxwell equations is utilized for isotropic media I ZZ ¶H E dl = m m  dA, (6) ¶t C A I ZZ ¶E H dl = # # + sE  dA, (7) 0 r ¶t C A where A is the area limited by the closed-loop C with length l, # and m are the electric 0 0 permittivity and magnetic permeability of free-space, respectively, and # , m are their r r relative equivalents. Our work focuses on a transverse electric two-dimensional (2D) case. The computa- tional domain is discretized into triangular elements using an efficient Delaunay triangula- tion algorithm, as depicted in Figure 1. Here, the electric-field components are tangential to the 2D domain and located along the edges of the triangles. It should be mentioned that the sampling of the magnetic components, which are considered normal to the domain, is not a trivial issue. In particular, two different approaches, depicted in Figure 1, are investigated. Axioms 2022, 11, 44 4 of 11 The initial one (Figure 1a) places the magnetic components at the triangle circumcenter, namely the center of the circle that is defined by its vertices creating thus, a Voronoi diagram [31]. The advantage of this approach is that the magnetic-field connection between two adjacent triangles is perpendicular to their common edge. This is important, since the inner product at the right-hand side of (7) suggests that the electric component is vertical to the virtual area of the connecting magnetic fields. Unfortunately, these points may lie outside the element for an obtuse triangle case, potentially influencing the algorithm’s stability. To circumvent this hindrance, a second technique can be used (Figure 1b), where the magnetic components are located at the triangle centroid, namely the averaged coordinates of its three vertices [32]. Although this approach ensures that each magnetic component lies inside the corresponding element, the magnetic-field connection between adjacent triangles is not necessarily vertical to the common electric component. Hence, an angular correction is required for the accurate implementation of (7). Note that the two techniques are identical for the case of equilateral triangles. In this manner, an optimal triangulation can be considered to consist of such elements. k k l l E E i,k H i,k i,l E ,J E E ,J i,j gr[i,j] i,l i,j gr[i,j] graphene graphene j j (a) Circumcenter sampling (b) Centroid sampling Figure 1. Topology of the unstructured triangular grid and the electromagnetic field assignment concerning the two magnetic field sampling techniques. In particular, the arrangement of the electric components is at the edges of the triangles, while the magnetic components lie inside the triangle (a) at its circumcenter or (b) at its centroid. Now, a finite integration was applied to (6) and (7). First of all, the temporal discretiza- tion of fields was defined to coincide with the conventional FDTD algorithm, i.e., E was updated at nDt and H at (n + 1/2)Dt instants, where Dt is the time-step. Then, considering element i, (6) was developed into n+1/2 n1/2 H H n n n i i E l + E l + E l = m m A, 0 r i,j i,k i,l i,j i,k i,l Dt Dt n+1/2 n1/2 n n n H = H (E l + E l + E l ). (8) i,j i,k i,l i,j i,k i,l i i A m m i 0 r The finite integration of (7) depends on the applied magnetic field component sam- pling technique. The circumcenter-sampling case provides the following straightforward update equation for E i,j n1 n1 n n E E E + E i,j i,j i,j i,j n1/2 n1/2 H H = # # + s d 0 r i,j i j Dt 2 2# # sDt 2Dt 1 0 r n n1 n1/2 n1/2 E = E + (H H ) . (9) i,j i,j i j 2# # + sDt 2# # + sDt d 0 r 0 r i,j i,j i,j i,l i,l Axioms 2022, 11, 44 5 of 11 For the centroid-sampling case, an angular correction must be included to model the finite integration accurately. Explicitly, 2# # sDt 2Dt 1 0 r n n1 n1/2 n1/2 E = E + (H H ) . (10) i,j i,j i j 2# # + sDt 2# # + sDt d sin q r r 0 0 i,j The electromagnetic properties (described by # , m , s) of each element are assumed r r to be constant within each element. Consequently, the tangential continuity of the electric component is satisfied a priori. Each set of d and l is shared between two adjacent i,j i,j elements. Additionally, E is shared between triangles i and j, and a unique definition of i,j its orientation must be selected. Herein, the positive direction is counterclockwise with respect to the element with a lower index (i in Figure 1). Thus, a negative contribution is assumed for the element with a larger index (j in Figure 1). Finally, a stability criterion should be established, similar to the conventional FDTD algorithm. In essence, the time increment should satisfy the relation [31] 0 1 B C 1 2A B C u i Dt 6 minB C, (11) @ t l A i,fg i,fg where A is the area of the element i, l the length of the common edge with triangle i,fg fg, d the magnetic-field distance from the corresponding triangle (Figure 1), and c the i,fg free-space speed of light. Particularly, the inner expression of (11) was calculated for all i,fg triangles, where the ratio corresponds to the cell’s dimensions since it is connected i,fg to all three triangle edges. Then, the stable time increment was defined via the mini- mum value of the aforementioned calculations. It is worth mentioning that this criterion converges to the popular Courant–Friedrichs–Lewy condition of the conventional FDTD structured discretization. 2.3. Graphene Modeling within Unstructured Triangular Grids Graphene was embedded in the computational domain in terms of an equivalent surface current J = s E, as demonstrated in Figure 1. Observe that the locations of gr gr the surface current components are identical to the corresponding electric ones. The contribution of this surface current is included into the computational domain via Ampere’s law, namely (7) via I ZZ ZZ ¶E H dl = # # + sE  dA + J  dA. (12) 0 r gr ¶t C A A Consequently, the electric-field updating Equation (9), for the circumcenter sampling case, was properly modified to 2# # sDt 2Dt 1 0 r n n1 n1/2 n1/2 E = E + (H H ) i,j i,j i j 2# # + sDt 2# # + sDt d 0 r 0 r i,j 2Dt 1 n1/2 J . (13) gr[i,j] 2# # + sDt d 0 r i,j One may observe the d parameter in the denominator, which is attributed to the i,j fact that graphene is modeled as a surface current density, instead of a volume density. Therefore, the infinitesimally thin nature of the material was retained. Similarly, (10), for the centroid sampling case, was modified to Axioms 2022, 11, 44 6 of 11 2# # sDt 2Dt 1 0 r n n1 n1/2 n1/2 E = E + (H H ) i,j i,j i j 2# # + sDt 2# # + sDt d sin q 0 r 0 r i,j 2Dt 1 n1/2 J . (14) gr[i,j] 2# # + sDt d sin q 0 r i,j Now, the main concern is the modeling of graphene’s frequency dispersion. Firstly, a transition of J = s E and the surface conductivity (3) into the time domain should be gr gr performed using an inverse Fourier transform (IFT) m IFT 2Gt J (w) = s (w)E (w) = E (w) ! J (t) = A e  E (t). (15) gr i,j i,j m i,j gr[i,j] gr[i,j] c jw + 2G Although a procedure that involves a convolution should be performed for the surface current updating, the efficient RCM concept is utilized, i.e., n n o n+1/2 2G(nm)Dt m J = A e E Dt , å c i,j gr[i,j] m=0 n+1/2 2GDt n1/2 n J = A e J + A DtE . (16) m m c c i,j gr[i,j] gr[i,j] Interestingly, there is no connection between the frequency-dispersion modeling and the computational-domain discretization with unstructured triangular meshes. Therefore, more advanced schemes, such as interband conductivity approximation, third-order re- sponse, or anisotropic conductivity due to a magnetostatic-bias application, for 3D domains, can be straightforwardly combined with the triangular grids. Additionally, other complex dielectrics are easily importable via equivalent frequency-dispersion schemes, such as the auxiliary differential equation of conventional FDTD, since they are, also, independent of the computational-domain discretization. The entire procedure of the proposed algorithm is summarized in the block diagram of Figure 2. Here, every individual step is included, such as the initialization and the electromagnetic component updating via the well-known leap-frog scheme. Initialization Simulation n n+1/2 H locations H updating eq. (7) Mesh E locations E updating generation eq. (12) or (13) J updating gr eq. (15) Element size Figure 2. Block diagram of the proposed scheme. 3. Proposed Algorithm Validation The validity of the featured algorithm is pursued through the well-established setup of surface wave propagation onto a straight free-standing graphene layer of infinite size. The advantage of this example is its direct comparison to the analytical expressions (4) and (5). Moreover, the conventional FDTD method is able to provide accurate results for this orthog- onal setup, in order to calculate reliably the global error within computational domain. To this aim, graphene chemical potential was selected m = 0.2 eV, and its scattering rate was set to G = 0.11 meV, while the analysis was performed at the far-infrared spectrum, namely from 1 to 3 THz. The dimension of the computational domain is 600 m  150 m, with the graphene layer located along the x-axis. As discussed, an unstructured Delaunay trian- gulation was utilized with a maximum edge size D = 3 m, resulting in approximately max 32,000 elements. The time-step was selected 3.7 fs to ensure the algorithm’s stability, while open boundaries were terminated via an 8-cell thick perfectly matched layer (PML). It is Axioms 2022, 11, 44 7 of 11 important to mention that PML is selected instead of an Absorbing Boundary Condition since it is possible to extend graphene inside the PML region. This leads to the effective ab- sorption of the strongly confined SPP waves. Finally, 8-cells are adequate for the considered applications due to the inherent losses of graphene surface waves. The numerically extracted propagation properties of the supported surface wave on graphene are compared to the aforementioned theoretical calculations in Figure 3. It is clear that the matching of the results is evident, especially until 2.5 THz. Beyond this frequency, the SPP wavelength is less than 30 m, and the conventional l/10 discretization limit is exceeded, leading to small discrepancies. It is important to mention that the simulated results represent both magnetic-field sampling techniques, since their differences is negligible. Furthermore, the distribution of the magnetic field at 2 THz is demonstrated in Figure 4, clearly indicating the propagation of the strongly confined surface wave on the graphene surface. Consequently, the proposed integration of the infinitesimally thin material into an unstructured triangular grid was validated successfully. −3 SPP −4 SPP −5 Theoretical Simulation −6 1 1.5 2 2.5 3 Frequency [THz] Figure 3. Numerical validation of surface wave propagation characteristics for a graphene layer with m = 0.2 eV and G = 0.11 meV. −50 −200 −150 −100 −50 0 50 100 150 200 x-axis [mm] Figure 4. Distribution of the magnetic-field component for a graphene layer with m = 0.2 eV and G = 0.11 meV at 2 THz. Finally, the global error in the computational domain was investigated for different element sizes, to check the influence of the element size. For this reason, a numerical simulation through the conventional orthogonal grid of the FDTD method was conducted using a considerably fine mesh (Dx = Dy < l /20). The mean global error was calculated for the SPP magnetic-field component, while it was normalized to the mean magnetic-field value of the computational domain. The calculated results are depicted in Figure 5, and it is apparent that smaller triangle dimensions lead to more precise results. Note that the shape of the curves is, generally, altered with D since different triangular meshes are generated. In addition, the max magnetic-field sampling at the triangle circumcenter (technique 1) presents a slightly lower global error. In any case, the satisfaction of the well-known D < l /10 rule secures max SPP reliable numerical results. Length (mm) y-axis [mm] Axioms 2022, 11, 44 8 of 11 0.5 Technique 1 Technique 2 0.4 0.3 D =10 mm max D =5 mm 0.2 max 0.1 D =3 mm max 1 1.5 2 2.5 3 Frequency [THz] Figure 5. Global error of the magnetic field, normalized to the mean computational domain value, for different triangle size discretization of the two proposed techniques. 4. Incident Field toward a Circular Graphene Scatterer Although the previous example validated the accuracy of the proposed methodology, the computational efficiency is not optimal compared to the conventional FDTD algorithm. The main reason for the degraded performance is the requirement for the storage of ad- ditional parameters, such as triangle dimensions as well as pointers to identify adjacent elements. For this reason, a second example is discussed in this section to exploit the advantages of the triangular computational domain. Particularly, a plane wave propagates toward a periodic structure of circular graphene scatterers of radius R = 20 m and period- icity of w = 100 m, as illustrated in Figure 6. The graphene chemical potential is set to 0.2 eV, and its scattering rate is 0.11 meV at the frequency range 1–3 THz. The computational domain is divided into 23,000 triangular elements of maximum edge size D = 3 m, max while the time-step is selected 3.7 fs. The side boundaries satisfy periodic conditions, and the extracted quantity of interest is the configuration’s transmission coefficient. periodic boundaries graphene Figure 6. Setup of an incident plane wave toward a periodic structure of circular graphene scatterers. The extracted numerical results of the proposed algorithm are compared to the conven- tional staircase approximation of the FDTD algorithm, as well as the commercial software COMSOL Multiphysics commercial package [33] for reference. Figure 7 proves the excel- lent behavior of our methodology. On the other hand, the staircase approximation fails to track the resonance frequencies accurately, even when using a very fine mesh (D = 1 m). This performance is explained via the distribution of the electric field in the scatterer in Figure 8. Specifically, the edge effects of the staircase approximation degrade the plasmonic resonances on the graphene, in contrast to the proposed method, where the electric field is smoothly distributed on the material surface. Relative global error Axioms 2022, 11, 44 9 of 11 −5 −10 COMSOL −15 Proposed method Staircase FDTD (D=5mm) Staircase FDTD (D=1mm) −20 1 1.5 2 2.5 3 Frequency [THz] Figure 7. Transmission coefficient comparison between the proposed method and the conventional staircase FDTD algorithm for a circular graphene scatterer of 20 m radius. [V/m] 0.8 0.6 0.4 0.2 −0.2 −0.4 −0.6 −0.8 −1 (a) Staircase (b) Proposed Figure 8. Comparison of the electric field distribution between the proposed method and the conventional staircase FDTD algorithm for a circular graphene scatterer of 20 m radius at 1.72 THz. 5. Conclusions The accurate time-domain modeling of curved graphene layers was accomplished in the present work using an unstructured triangular mesh. Initially, the fundamental characteristics of the two-dimensional material were investigated, as well as the finite integration scheme for the proposed computational grid. Particularly, two techniques were examined that coincide for equilateral triangle discretization. The strength of the first one relies on its more physical representation to achieve increased accuracy, but acute triangles are required to ensure stability. This problem is resolved via the second technique, although the necessity of further approximations slightly decreases the algorithm’s accuracy. Then, graphene was introduced into the numerical algorithm as an equivalent surface current density, while its frequency dispersion was treated via an RCM technique. The proposed methodology was thoroughly examined in terms of the surface wave propagation properties, validating the proposed implementation. Finally, a circular graphene scatterer setup was designed to exploit the algorithm’s applicability and advantages, compared to the conventional staircase approximation. Transmission coefficient [dB] Axioms 2022, 11, 44 10 of 11 Author Contributions: Conceptualization, S.A.; methodology, S.A.; validation, S.A.; writing— original draft preparation, S.A.; writing—review and editing, T.Z. and N.K.; supervision, T.Z. and N.K. All authors have read and agreed to the published version of the manuscript. Funding: This research was cofinanced by Greece and the European Union (European Social Fund- ESF) through the Operational Programme «Human Resources Development, Education and Lifelong Learning» in the context of the project “Reinforcement of Postdoctoral Researchers—2nd Cycle” (MIS-5033021), implemented by the State Scholarships Foundation (IKY). Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: The data presented in this study are available via a straightforward reproduction of the provided simulation details. Conflicts of Interest: The authors declare no conflict of interest. Abbreviations The following abbreviations are used in this manuscript: FDTD Finite-Difference Time-Domain PML Perfectly Matched Layer RCM Recursive Convolution Method SPP Surface Plasmon Polariton References 1. Novoselov, K.; Geim, A.K.; Morozov, S.; Jiang, D.; Zhang, Y.; Dubonos, S.; Grigorieva, I.; Firsov, A. Electric field effect in atomically thin carbon films. Science 2004, 306, 666–669. [CrossRef] [PubMed] 2. Mikhailov, S.A.; Ziegler, K. New electromagnetic mode in graphene. Phys. Rev. Lett. 2007, 99, 016803. [CrossRef] [PubMed] 3. Grigorenko, A.N.; Polini, M.; Novoselov, K. Graphene plasmonics. Nat. Photonics 2012, 6, 749–758. [CrossRef] 4. Bludov, Y.V.; Ferreira, A.; Peres, N.M.; Vasilevskiy, M.I. A primer on surface plasmon-polaritons in graphene. Int. J. Mod. Phys. B 2013, 27, 1341001. [CrossRef] 5. Ansell, D.; Radko, I.; Han, Z.; Rodriguez, F.; Bozhevolnyi, S.; Grigorenko, A. Hybrid graphene plasmonic waveguide modulators. Nat. Commun. 2015, 6, 8846. [CrossRef] [PubMed] 6. Xiao, S.; Zhu, X.; Li, B.H.; Mortensen, N.A. Graphene-plasmon polaritons: From fundamental properties to potential applications. Front. Phys. 2016, 11, 117801. [CrossRef] 7. Li, Y.; Tantiwanichapan, K.; Swan, A.K.; Paiella, R. Graphene plasmonic devices for terahertz optoelectronics. Nanophotonics 2020, 9, 1901–1920. [CrossRef] 8. Zhang, J.; Hong, Q.; Zou, J.; He, Y.; Yuan, X.; Zhu, Z.; Qin, S. Fano-resonance in hybrid metal-graphene metamaterial and its application as mid-infrared plasmonic sensor. Micromachines 2020, 11, 268. [CrossRef] 9. Choroszucho, A. Analysis of the influence of the complex structure of clay hollow bricks on the values of electric field intensity by using the FDTD method. Arch. Electr. Eng. 2016, 65, 745–759. [CrossRef] 10. Oskooi, A.F.; Roundy, D.; Ibanescu, M.; Bermel, P.; Joannopoulos, J.D.; Johnson, S.G. MEEP: A flexible free-software package for electromagnetic simulations by the FDTD method. Comput. Phys. Commun. 2010, 181, 687–702. [CrossRef] 11. Bouzianas, G.D.; Kantartzis, N.V.; Yioultsis, T.V.; Tsiboukis, T.D. Consistent study of graphene structures through the direct incorporation of surface conductivity. IEEE Trans. Magn. 2014, 50, 161–164. [CrossRef] 12. Ramadan, O. Improved direct integration auxiliary differential equation FDTD scheme for modeling graphene drude dispersion. Optik 2020, 219, 165173. [CrossRef] 13. Nayyeri, V.; Soleimani, M.; Ramahi, O.M. Modeling graphene in the finite-difference time-domain method using a surface boundary condition. IEEE Trans. Antennas Propag. 2013, 61, 4176–4182. [CrossRef] 14. Wang, X.H.; Yin, W.Y.; Chen, Z. Matrix exponential FDTD modeling of magnetized graphene sheet. IEEE Antennas Wirel. Propag. Lett. 2013, 12, 1129–1132. [CrossRef] 15. Amanatiadis, S.A.; Kantartzis, N.V.; Ohtani, T.; Kanai, Y. Precise modeling of magnetically biased graphene through a recursive convolutional FDTD method. IEEE Trans. Magn. 2017, 54, 7201504. [CrossRef] Axioms 2022, 11, 44 11 of 11 16. Ramadan, O. Simplified FDTD Formulations for Magnetostatic Biased Graphene Simulations. IEEE Antennas Wirel. Propag. Lett. 2020, 19, 2290–2294. [CrossRef] 17. Mock, A. Modeling passive mode-locking via saturable absorption in graphene using the finite-difference time-domain method. IEEE J. Quantum Electron. 2017, 53, 1300210. [CrossRef] 18. Amanatiadis, S.A.; Ohtani, T.; Zygiridis, T.T.; Kanai, Y.; Kantartzis, N.V. Modeling the Third-Order Electrodynamic Response of Graphene via an Efficient Finite-Difference Time-Domain Scheme. IEEE Trans. Magn. 2019, 56, 6700704. [CrossRef] 19. Mock, A. Padé approximant spectral fit for FDTD simulation of graphene in the near infrared. Opt. Mater. Express 2012, 2, 771–781. [CrossRef] 20. Amanatiadis, S.; Zygiridis, T.; Ohtani, T.; Kanai, Y.; Kantartzis, N. A Consistent Scheme for the Precise FDTD Modeling of the Graphene Interband Contribution. IEEE Trans. Magn. 2021, 57, 1600104. [CrossRef] 21. Xiao, S.; Wang, T.; Liu, Y.; Xu, C.; Han, X.; Yan, X. Tunable light trapping and absorption enhancement with graphene ring arrays. Phys. Chem. Chem. Phys. 2016, 18, 26661–26669. [CrossRef] [PubMed] 22. Xu, H.; Zhang, Z.; Wang, S.; Liu, Y.; Zhang, J.; Chen, D.; Ouyang, J.; Yang, J. Tunable graphene-based plasmon-induced transparency based on edge mode in the mid-infrared region. Nanomaterials 2019, 9, 448. [CrossRef] [PubMed] 23. Liu, T.; Zhou, C.; Jiang, X.; Cheng, L.; Xu, C.; Xiao, S. Tunable light trapping and absorption enhancement with graphene-based complementary metasurfaces. Opt. Mater. Express 2019, 9, 1469–1478. [CrossRef] 24. Liu, P.; Cai, W.; Wang, L.; Zhang, X.; Xu, J. Tunable terahertz optical antennas based on graphene ring structures. Appl. Phys. Lett. 2012, 100, 153111. [CrossRef] 25. Ye-xin, S.; Jiu-sheng, L.; Le, Z. Graphene-integrated split-ring resonator terahertz modulator. Opt. Quantum Electron. 2017, 49, 350. [CrossRef] 26. Wu, L.; Liu, H.; Li, J.; Wang, S.; Qu, S.; Dong, L. A 130 GHz electro-optic ring modulator with double-layer graphene. Crystals 2017, 7, 65. [CrossRef] 27. Bao, Z.; Tang, Y.; Hu, Z.D.; Zhang, C.; Balmakou, A.; Khakhomov, S.; Semchenko, I.; Wang, J. Inversion method characterization of graphene-based coordination absorbers incorporating periodically patterned metal ring metasurfaces. Nanomaterials 2020, 10, 1102. [CrossRef] 28. Gusynin, V.; Sharapov, S.; Carbotte, J. Magneto-optical conductivity in graphene. J. Phys. Condens. Matter 2006, 19, 026222. [CrossRef] 29. Hanson, G.W. Dyadic Green’s functions and guided surface waves for a surface conductivity model of graphene. J. Appl. Phys. 2008, 103, 064302. [CrossRef] 30. Taflove, A.; Hagness, S.C.; Piket-May, M. Computational electromagnetics: The finite-difference time-domain method. In The Electrical Engineering Handbook; Elsevier Inc.: New York, NY, USA, 2005; Volume 3. 31. Lee, C.; McCartin, B.; Shin, R.; Kong, J.A. A triangular-grid finite-difference time-domain method for electromagnetic scattering problems. J. Electromagn. Waves Appl. 1994, 8, 449–470. [CrossRef] 32. Liu, Y.; Sarris, C.D.; Eleftheriades, G.V. Triangular-mesh-based FDTD analysis of two-dimensional plasmonic structures supporting backward waves at optical frequencies. J. Light. Technol. 2007, 25, 938–945. [CrossRef] 33. COMSOL Multiphysics v. 5.5, COMSOL AB. 2019. Available online: www.comsol.com (accessed on November 2021).

Journal

Axioms – Multidisciplinary Digital Publishing Institute

Published: Jan 22, 2022

Keywords: auxiliary differential equation; curved modeling; equilateral triangle; FDTD; finite integration; surface current density; surface wave

There are no references for this article.