I Introduction
Photoconductive devices (PCDs) are promising candidates for terahertz (THz) source generation and signal detection because they are compact and frequencystable, and capable of operating at room temperature with low optical input power levels [1, 2, 3, 4, 5]. However, the low opticaltoTHz conversion efficiency of the conventional PCDs has hindered their widespread use in applications of THz technologies. In recent years, significant progress has been made in alleviating this bottleneck. Introduction of metallic/dielectric nanostructures inside or on top of the PCDs’ active region has been shown to increase the conversion efficiency by several orders of magnitude [3, 4, 5, 6]. This significant increase is attributed to several mechanisms. The optical electromagnetic (EM) field is enhanced due to plasmonic [4] or Mie [5] resonances. Nanostructured electrodes reduce the effective distance that the carriers have to travel [7, 6]. Furthermore, recent studies show that the nanostructuretailored bias/static electric field and the carrier screening effect also play important roles [1, 8, 9]. The interplay between these mechanisms makes the relevant device design very complicated. In this context, simulation tools have become indispensable in understanding these physical mechanisms and their coupling and in accelerating the design process.
Due to their geometricallyintricate structure and the complicated EM wave/field interactions they support, simulation of nanostructured PCDs can not be carried out using the methods that have been developed for conventional PCDs and rely on semianalytical approximations/computations [10, 11, 12, 13, 14]. To this end, the finite element method (FEM) has been extensively used in recent years [15, 16, 17, 18]
. The optical EM field is calculated with FEM in frequency domain and is used for predicting the carrier generation rate distribution in space. The time dependency of the generation rate is approximated with an analytical expression. Because the frequencydomain solutions are used, the (nonlinear) coupling between carrier dynamics and EM fields is not fully accounted for and consequently the carrier screening effect cannot be modeled by this approach
[15]. Recently, a multiphysics framework making use of discontinuous Galerkin (DG) methods has been proposed [19, 20]. This framework solves a coupled system of Poisson and stationary driftdiffusion (DD) equations describing the nonequilibrium steady state of the PCD and a coupled system of timedependent Maxwell and DD equations describing the transient stage that involves the optoelectronic response of the PCD. The nonlinear coupling between the electrostatic fields and the carriers and that between the EM fields and the carriers are taken into account by the Gummel method and through the use of an explicit time integration scheme, respectively. Even though this DGbased framework provides higher flexibility and higherorder accuracy in both space discretization and time integration and is more robust in modeling nonlinear coupling mechanisms compared to the FEMbased schemes, it is still computationally demanding especially for practical threedimensional (3D) devices because of the multiple space and time characteristic scales involved in modeling and discretizing a nanostructured PCD (e.g., the Debye length , the optical wavelength , and the device size ) [19, 20].One way to reduce this high computational cost is to make use of the nanostructure’s periodicity, i.e., model and discretize the multiphysical interactions and their coupling on a unit cell to approximate their behavior on the whole device. This approach calls for proper boundary conditions to be enforced on the surfaces of the unit cell. The periodic boundary conditions (PBCs) required by the optical EM field simulation have been wellstudied [21, 7, 6, 15, 22]. However, since a PCD is in a nonequilibrium steadystate under a bias voltage [20, 23], the boundary conditions required by the simulation of the carrier dynamics to be enforced on the unitcell surfaces are not trivial. In [15], it is assumed there is no potentialdrop along the direction perpendicular to the bias electric field and the optical EM field excitation. A PBC is used along this direction, which reduces the computational domain of the carrier dynamics simulation to a strip containing a chain of unit cells. Even with this approach, the reported computational requirement is high [15]. In addition, the nonlinear coupling is not fully considered in [15] since a frequencydomain FEM is used to compute the EM field distribution.
In this work, to reduce the high computational cost of simulating nanostructured PCDs, a unitcell scheme is proposed within the multiphysics DG framework developed in [19, 20]. For Poisson equation, a “potentialdrop” boundary condition (PDBC) is enforced on the opposing surfaces of the unitcell (along the direction of the bias electric field). For carriers and EM fields, PBCs are enforced on the unitcell surfaces. All boundary conditions are “weakly” enforced using the numerical flux of the DG scheme. Numerical results demonstrate that the proposed DGbased unitcell scheme has almost the same accuracy as the DG framework that takes into account the whole device in predicting the THz photocurrent output. It also retains the main advantages of the DG framework [19], while it significantly reduces the computational cost, making it feasible to simulate practical 3D devices on desktop computers.
The rest of paper is organized as follows. In Section IIA, the unitcell model and the relevant boundary conditions are introduced. In Sections IIB and IIC, the coupled systems of Poisson and stationary DD equations and timedependent Maxwell and DD equations are provided with the boundary conditions for the unit cell, respectively. These sections also describe the DG schemes used for discretizing and solving these coupled systems of equations. Numerical results that validate the accuracy of the proposed scheme and demonstrate its computational benefits are provided in Section III. Finally, Section IV concludes the paper and provides some remarks on the future research directions.
Ii The Unitcell Model and Numerical Scheme
Iia The unitcell model
Fig. 1 illustrates an example of a 3D nanostructured PCD. The device consists of a photoconductive layer (LTGaAs), a substrate layer (SIGaAs), two electrodes that are deposited on the photoconductive layer, and a metallic nanograting that is placed between them. The grating is designed to support plasmonic modes that enhance the EM fields induced on the structure upon excitation by an optical EM wave [5].
The operation of PCDs can be analyzed into two stages. Initially, a bias voltage is applied to the electrodes. The resulting (static) electric field changes the carrier distribution. The redistributed carriers in turn affects the electric potential distribution. The device reaches a nonequilibrium steadystate mathematically described by a coupled system of Poisson and stationary DD equations [23]. When an optical EM excitation (i.e., optical pump) is incident on the device, a transient stage starts. The photoconductive material absorbs the EM energy induced on the device due to the optical excitation and generates carriers. The carriers are driven by both the bias electric field and the optical EM fields. The carrier dynamics and EM wave/field interactions are mathematically described by a coupled system of the timedependent Maxwell and DD equations [19, 24].
To accurately model these coupled physical processes that occur on the whole device only within a single unit cell [as illustrated in Fig. 1 (b)], appropriate boundary conditions must be enforced on the unitcell surfaces.
For Poisson equation, one can not simply use PBCs since the potential drops from the anode to the cathode. In the electrostatic problem, the nanostructure generates only local variations in the potential [25] and on average the potential drops linearly between the two electrodes (as in the homogeneous case without the nanostructure) [26]. Furthermore, since all unit cells are the same, the potential drop and the local potential variation within each unit cell should be the same. This analysis suggests that the bias static electric field, which is equal to the gradient of the potential, is periodic. Therefore, the steadystate carrier densities and the field dependent mobility should also be periodic. Furthermore, since the nanostructures and the carrier distributions are periodic, the optical EM fields induced on the structure is the same in all unit cells. The same argument applies to carrier dynamics since the mobility, the static electric field, and the optical EM fields are the same in all unit cells. Therefore, PBCs can be used for stationary DD equations, timedependent Maxwell equations, and timedependent DD equations.
The boundary conditions discussed above are mathematically described as follows. For Poisson equation, the boundary conditions are
(1)  
(2) 
where is the electric potential, and are the widths of the unit cell in  and directions, respectively, and is the potentialdrop between the two faces of the unit cell perpendicular to the direction. In the rest of the text, (1) is termed as PDBC. Note that the PBC (2) is used along the direction because the potential does not drop along this direction (hence is periodic) [15]. The potential drop function in (1) is selected as described next. is positiondependent, e.g., the potential drop becomes smaller away from the electrodes in the direction. Since the height of the device is much smaller than its width [1, 2, 3] and the electrodes extend to the whole width of the device, the potential drop is approximately the same for all and at a given value of (i.e., on a surface perpendicular to direction). Therefore, can be simplified to a single value . Additionally, as discussed before, on average the potential drops linearly between the two electrodes [26, 27]. Consequently,
can be estimated as:
, where the is the bias voltage applied to the electrodes and is the distance between them.For stationary DD equations and timedependent DD and Maxwell equations, PBCs are used:
(3)  
(4) 
where represents the steadystate electron/hole density, the transient electron/hole density, or the transient electric/magnetic field. For all equations, the boundary conditions on the top and bottom surfaces of the unit cell (surfaces perpendicular to direction) are the same as those would be used in the fulldevice simulation [20, 19].
Two comments about the unitcell model introduced in this section are in order: (i) The approximation of using a single value for potential drop can be improved by estimating from the solution of the same device but without the nanostructure (which generates only local variations in the potential). Modeling a simpler device without the nanostructure is easier since a much coarser mesh can be used. (ii) The nanostructure does not have to be metallic for the unitcell model to hold; it is also applicable when nanostructure is made of dielectric material [5].
IiB The unitcell PoissonDD solver
The coupled system of Poisson and stationary DD equations is solved with the Gummel iteration method, which, at every iteration, linearized the Poisson equation. The DD equations are directly treated as linear [23]
. In each iteration, one needs to solve three partial differential equations (PDEs), namely the linearized Poisson equation and two DD equations
[20]. Those equations, together with the boundary conditions described above, are discretized using DG methods.The Poisson equation in the unit cell is expressed as the following boundary value problem (BVP)
(5)  
(6)  
(7)  
(8)  
(9) 
where and are the unknowns to be solved for, and are known coefficients, is the permittivity, is the solution domain, represents the surfaces perpendicular to the direction, ,
is the outward pointing normal vector on
. The homogeneous Neumann boundary condition is used in (9) [19].Discretizing into a set of nonoverlapping elements, DG approximates the unknowns with basis functions (the nodal basis function [28] is used in this work) in each element and applies Galerkin testing to the resulting equations [29, 28, 30]. This process yields a matrix system as
(10) 
Here, and are unknown vectors storing the basis expansion coefficients, and are block diagonal mass matrices, and are block sparse matrices representing the gradient and divergence operators, respectively, and and have contributions from the tested force term and boundary conditions. Detailed expressions of these vectors and matrices can be found in [20].
The continuity of solutions across element interfaces is enforced through a uniquely defined numerical flux. For Poisson equation, the alternate flux is used; it is defined as [29]
(11)  
(12) 
on each element surface in the interior of , where “average” and “jump” operators are defined as and , respectively, with being a scalar or a vector variable. The same definitions of these operators are used throughout this paper. Superscripts “” and “+” refer to variables defined in the interior and exterior of the surface, respectively. determines the upwind direction of and [29, 28, 30]. On , the numerical fluxes are chosen as and [29]. Here, the variables are defined on the surfaces and the explicit dependency on is dropped for the simplicity of the notation.
To enforce PBC, same meshes are created on the opposing surfaces of the unit cell and each element face on a given surface is “connected” to its counterpart on the opposing surface. This means that when calculating the numerical flux on the boundary, (11)(12) are used but, for each element face, the exterior variable is taken from its “connected” face on the opposing surface
(13)  
(14) 
where .
For PDBC, the element faces on the opposing surfaces are “connected” just like it is done for the PBC. The exterior variables in the numerical flux expressions are set as
(15)  
(16)  
(17)  
(18) 
The electron DD equations in the unit cell are expressed as the following BVP [20]
(19)  
(20)  
(21)  
(22)  
(23) 
where and are the unknowns to be solved for, and , and are known coefficients within the Gummel method [20]. The homogeneous Robin boundary condition is used in (23) [19]. The BVP for holes only differs by the sign in front of the drift term.
DG discretization of (19)(23) yields a matrix system as
(24) 
where and are unknown vectors storing the basis expansion coefficients, is same as , and are same as before, is a diagonal matrix, and corresponds to the drift term [20].
IiC The unitcell MaxwellDD solver
The coupled system of timedependent Maxwell and DD equations is integrated in time using an explicit scheme. The nonlinear coupling between the Maxwell and DD equations is accounted for by alternately feeding the updated solutions into each other during the time integration [19]. Each set of equations is discretized using a timedomain DG (DGTD) method [29, 28, 30] that uses its own timestep size [19].
The timedependent electron DD equations in the unit cell are expressed as an initialboundary value problem (IBVP) [19]
(27)  
(28)  
(29)  
(30)  
(31) 
where and are the unknowns to be solved for, and are known coefficients, is the netrecombination, is the transient recombination rate, and is the generation rate [14, 19]. Here, to simplify the notation, the subscript “e” is dropped. The IBVP for holes only differs by the sign in front of the drift term.
Apart from the time dependency and the time derivative term on the left side of (IIC), the system has the same form as (19)(23). Using the same spatial discretization as that used for (19)(23), one can obtain the semidiscrete form [19]
(32)  
(33) 
where and unknown vectors storing timedependent basis expansion coefficients, and are vectors constructed from the netrecombination and boundary conditions (other than PBC), , /, /, and / are the elemental mass, convection, gradient, and divergence matrices, respectively [19].
The numerical fluxes and boundary conditions are the same as those used in the stationary DD equations. The semidiscretized system (32)(33) is integrated in time using an explicit thirdorder totalvariationdiminishing RungeKutta scheme [31].
Likewise, Maxwell equations in the unit cell are expressed as the following IBVP
(34)  
(35)  
(36)  
(37)  
(38) 
where and are the transient electric and magnetic fields, is the permeability, is the transient current density, and . Discretizing (34)(38) with the DGTD method [28, 32, 33, 34, 35, 36] yields
(39)  
(40)  
(41)  
(42)  
(43)  
(44) 
where and are unknown vectors storing timedependent basis expansion coefficients, is the current density vector, and are the constant permittivity and permeability in element , , and , where , , and are the mass, stiffness, and surface mass matrices, respectively [32, 19].
In (39)(44), and are the components of and , respectively, . The upwind flux [28] is used here:
(45)  
(46) 
where and are the wave impedance and wave admittance, respectively. The PBCs (36) and (37) are enforced similarly as (13)(14) and (25)(26), with . The perfect electric conductor (PEC) boundary condition is used on , , and perfectly matched layers (PMLs) [37, 35, 36, 38] are used for absorbing outgoing EM waves. The semidiscrete system (39)(44) is integrated in time using a fivestage fourthorder RungeKutta method [28].
Iii Numerical Results
To validate the proposed scheme, a 3D nanostructured PCD is simulated using the proposed method (with only a unit cell) and the results are compared to those obtained for the full device using the DG framework in [20, 19]. The device is illustrated in Fig. 1. The thickness of the LTGaAs and the SIGaAs layers is m. The nanograting has a periodicity of m in x and ydirections and its height is m. The metal block is a truncated square pyramid with dimensions of m and m on its top and bottom, respectively. The material parameters are same as those used in [19] and are provided in Table I. The DD equations are solved only within the LTGaAs layer while Poisson equation and the timedependent Maxwell equations are solved everywhere in the unit cell (and the full device to obtain the reference results). V, and two continuous lasers with frequencies and excite the devices from top. The laser beam width is in the fulldevice simulation.
C  cm 

cm  
LTGaAs  , , , 
SIGaAs  , 
Mobility  
Recombination 
For the unitcell simulation, the domain is shown in Fig. 1 (b) and the boundary conditions are those explained in Section II. For the fulldevice simulation, the height and width of the device (in  and directions, respectively) are m and m. The nanograting is made of a grid of unit cells. For the DD equations, the Dirichlet boundary condition , is used on electrode/semiconductor interfaces and the homogeneous Robin boundary condition is used on semiconductor/insulator interfaces [23]. For Poisson equation, the Dirichlet boundary condition is enforced on electrode surfaces and the homogeneous Neumann boundary condition is used to truncate the computation domain [23]. For Maxwell equations, the computation domain is truncated with PMLs [37] backed with PEC.
The simulation domains are discretized with tetrahedrons (Fig. 1). The minimum and the maximum edge lengths in the mesh are and , respectively. The numbers of elements are and in the fulldevice and unitcell simulations, respectively. The tolerance of the Gummel iteration is and the solution typically converges after iterations. The linear systems are solved using the GMRES iterative solver and an ILU preconditioner is reused throughout the Gummel iterations [20]. The physical duration of the transient stage is and the timestep sizes for the Maxwell and DD equations are and , respectively [19]. Table II provides the computational times (measured in “corehour”) required by the unitcell and fulldevice simulations to complete the steady state and the transient stage. The unitcell scheme’s steady state and transient stage simulations are and times faster.
steady state  transient stage  

core  hour  corehour  core  hour  corehour  
full device  32  201  6432  4096  35  143360 
unit cell  16  0.22  3.52  32  3  96 

Tested on Shaheen II (https://www.hpc.kaust.edu.sa/), a Cray XC40 system consists of 6174 nodes based on Intel Xeon (R) E52698 v3.
It should be noted here that the unitcell simulation can be carried out on a workstation, while the fulldevice simulation requires of memory and calls for a parallelized solver and a distributedmemory system.
Fig. 2 shows the solutions obtained by the fulldevice simulation. Fig. 2 (a) shows the electric potential distribution on the interface between the nanostructures and the photoconductive layer. One can clearly see that the potential drops equally across each unit cell and the local variations in all unit cells are approximately the same. Figs. 2 (b) and (c) show the electric potential and electric field on several lines along the direction, respectively. The dash lines mark the positions of the unitcell surfaces. Fig. 2 (b) shows that, although the potential distributions at different positions are different, the potential drops are approximately the same. The linear drop estimation agrees with the solution very well on all unitcell surfaces. Fig. 2 (c) shows the bias electric field is periodic as we expected.
Figs. 3 (a) and (c) show the steadystate electric potential and electron density calculated from the unitcell scheme, respectively. For comparison, Figs. 3 (b) and (d) show those calculated using the fulldevice, where the solutions are set transparent (except those on the center unit cell) for better visualization and comparison. Very good agreement between two sets of results is observed. From Figs. 3 (a) and (b), one can see that although only the potential difference between boundaries is used in the unitcell scheme, the potential variation inside the unit cell is same as that obtained from the fulldevice simulation. Since the bias electric fields are the same in all unit cells (the electric potential in different unit cells only differs by a constant), the mobility and the carrier densities are periodic. Fig. 3 (d) shows the electron density distribution in the fulldevice. And, Fig. 3 (c) shows that the solution obtained from the unitcell scheme is same as that obtained from the fulldevice simulation.
Fig. 4 compares the transient current densities obtained from the unitcell and fulldevice simulations. The results agree well. To demonstrate the effect of the nanostructure on the device output, the current density obtained on the same device but without the nanostructure is also shown in the figure. It is clear that the photocurrent density increases significantly with the introduction of the plasmonic nanostructure yielding an enhancement factor of .
It should be noted here that the unitcell scheme assumes an infinitely large pumping aperture (the optical pumping is same for all unit cells). In the fulldevice simulation, the pumping beam and the device have finite widths, which results in two effects: (1) The pumping power near the center cells is higher than that near the boundary cells and (2) the pumping wave is scattered by the electrodes and the / boundaries of the device. These effects lead to a small difference between the transient current densities obtained by the two methods.
Iv Conclusion
The large scale of 3D nanostructured PCDs and various multiscale and nonlinearly coupled physical phenomena involved in their operation render their direct simulation computationally very costly. The unitcell scheme developed in this work dramatically reduces this computational cost, while the complex nonlinear EM/carrier interactions are still accurately accounted for. This scheme solves coupled systems of Poisson and stationary DD and timedependent Maxwell and DD equations in a unit cell which represents one period of the nanostructure. The coupled equations systems are discretized using DG schemes. A PDBC is enforced on the unitcell surfaces for Poisson equation, while PBCs are used for stationary and timedependent DD equations and timedependent Maxwell equations. These BCs are weakly enforced using the numerical flux of the DG methods.
Numerical results demonstrate that the proposed unitcell DG scheme maintains the accuracy of a DG scheme that operates on the full device while it is more than times faster. The unitcell scheme developed in this work can be used in design of PCDs not only with nanostructures but also antireflection layers and substrates. It can also be extended to account for organic devices operating with possibly more than two types of carriers.
Acknowledgment
This research is supported by the King Abdullah University of Science and Technology (KAUST) Office of Sponsored Research (OSR) under Award No 2016CRG52953. The authors would like to thank the KAUST Supercomputing Laboratory (KSL) for providing the required computational resources.
References
 [1] S. Lepeshov, A. Gorodetsky, A. Krasnok, E. Rafailov, and P. Belov, “Enhancement of terahertz photoconductive antenna operation by optical nanoantennas,” Laser Photonics Rev., vol. 11, no. 1, p. 1600199, 2017.
 [2] N. M. Burford and M. O. ElShenawee, “Review of terahertz photoconductive antenna technology,” Opt. Eng., vol. 56, no. 1, p. 010901, 2017.
 [3] J.H. Kang, D.S. Kim, and M. Seo, “Terahertz wave interaction with metallic nanostructures,” Nanophotonics, vol. 7, no. 5, pp. 763–793, 2018.
 [4] N. T. Yardimci and M. Jarrahi, “Nanostructureenhanced photoconductive terahertz emission and detection,” Small, vol. 14, no. 44, p. 1802437, 2018.
 [5] A. E. Yachmenev, D. V. Lavrukhin, I. A. Glinskiy, N. V. Zenchenko, Y. G. Goncharov, I. E. Spektor, R. A. Khabibullin, T. Otsuji, and D. S. Ponomarev, “Metallic and dielectric metasurfaces in photoconductive terahertz devices: a review,” Opt. Eng., vol. 59, no. 6, p. 061608, 2019.
 [6] S.H. Yang, M. R. Hashemi, C. W. Berry, and M. Jarrahi, “7.5% opticaltoterahertz conversion efficiency offered by photoconductive emitters with threedimensional plasmonic contact electrodes,” IEEE Trans. THz Sci. Technol., vol. 4, no. 5, pp. 575–581, 2014.
 [7] C. W. Berry, N. Wang, M. R. Hashemi, M. Unlu, and M. Jarrahi, “Significant performance enhancement in photoconductive terahertz optoelectronics by incorporating plasmonic contact electrodes,” Nat. Commun., vol. 4, p. 1622, 2013.
 [8] K. Moon, I. Lee, J.H. Shin, E. S. Lee, K. N, S. Han, and K. H. Park, “Bias field tailored plasmonic nanoelectrode for highpower terahertz photonic devices,” Sci. Rep., vol. 5, no. 13817, Nov 2012.
 [9] K. Moon, E. S. Lee, J. Choi, D. Lee, I.M. Lee, S.P. Han, H.S. Kim, and K. H. Park, “A comparative study of the plasmon effect in nanoelectrode THz emitters: Pulse vs. continuouswave radiation,” Appl. Phys. Lett., vol. 109, no. 7, p. 071105, 2016.
 [10] P. Kirawanich, S. J. Yakura, and N. E. Islam, “Study of highpower wideband terahertzpulse generation using integrated highspeed photoconductive semiconductor switches,” IEEE Trans. Plasma Sci., vol. 37, no. 1, pp. 219–228, 2008.
 [11] M. Khabiri, M. Neshat, and S. SafaviNaeini, “Hybrid computational simulation and study of continuous wave terahertz photomixers,” IEEE Trans. THz Sci. Technol., vol. 2, no. 6, pp. 605–616, 2012.
 [12] N. Khiabani, Y. Huang, Y.C. Shen, and S. Boyes, “Theoretical modeling of a photoconductive antenna in a terahertz pulsed system,” IEEE Trans. Antennas Propag., vol. 61, no. 4, pp. 1538–1546, 2013.
 [13] J. C. Young, D. Boyd, S. D. Gedney, T. Suzuki, and J. Liu, “A DGFETD port formulation for photoconductive antenna analysis,” IEEE Antennas Wireless Propag. Lett., vol. 14, pp. 386–389, 2014.
 [14] E. Moreno, M. F. Pantoja, S. G. Garcia, A. R. Bretones, and R. G. Martin, “Timedomain numerical modeling of THz photoconductive antennas,” IEEE Trans. THz Sci. Technol., vol. 4, no. 4, pp. 490–500, 2014.
 [15] N. Burford and M. ElShenawee, “Computational modeling of plasmonic thinfilm terahertz photoconductive antennas,” J. Opt. Soc. Am. B, vol. 33, no. 4, pp. 748–759, 2016.
 [16] M. J. MohammadZamani, M. Neshat, and M. K. MoravvejFarshi, “Nanoslit cavity plasmonic modes and builtin fields enhance the cw THz radiation in an unbiased antennaless photomixers array,” Opt. Lett., vol. 41, no. 2, pp. 420–423, 2016.
 [17] M. Bashirpour, S. Ghorbani, M. Kolahdouz, M. Neshat, M. MasnadiShirazi, and H. Aghababa, “Significant performance improvement of a terahertz photoconductive antenna using a hybrid structure,” RSC Advances, vol. 7, no. 83, pp. 53 010–53 017, 2017.
 [18] N. M. Burford, M. J. Evans, and M. O. ElShenawee, “Plasmonic nanodisk thinfilm terahertz photoconductive antenna,” IEEE Trans. THz Sci. Technol., vol. 8, no. 2, pp. 237–247, 2017.
 [19] L. Chen and H. Bagci, “Multiphysics modeling of plasmonic photoconductive devices using discontinuous Galerkin methods,” arXiv preprint arXiv:1912.03639, 2019.
 [20] L. Chen and H. Bagci, “Steadystate simulation of semiconductor devices using discontinuous Galerkin methods,” IEEE Access, vol. 8, pp. 16 203–16 215, 2020.
 [21] B. Heshmat, H. Pahlevaninezhad, Y. Pang, M. MasnadiShirazi, R. Burton Lewis, T. Tiedje, R. Gordon, and T. E. Darcie, “Nanoplasmonic terahertz photoconductive switch on GaAs,” Nano letters, vol. 12, no. 12, pp. 6255–6259, 2012.
 [22] K. Sirenko, Y. Sirenko, and H. Bagci, “Exact absorbing boundary conditions for periodic threedimensional structures: Derivation and implementation in discontinuous Galerkin timedomain method,” IEEE J. Multiscale and Multiphys. Comput. Techn., vol. 3, pp. 108–120, 2018.
 [23] D. Vasileska, S. M. Goodnick, and G. Klimeck, Computational Electronics: semiclassical and quantum device modeling and simulation. CRC press, 2010.
 [24] L. Chen and H. Bagci, “A discontinuous Galerkin framework for multiphysics simulation of photoconductive devices,” in Proc. Int. Appl. Comput. Electromagn. Symp. IEEE, 2019, pp. 1–2.
 [25] L. Chen, M. Dong, and H. Bagci, “Modeling floating potential conductors using discontinuous Galerkin method,” IEEE Access, vol. 8, pp. 7531–7538, 2020.
 [26] A. Menshov and V. I. Okhmatovski, “Surface–volume–surface electric field integral equation for magnetoquasistatic analysis of complex 3D interconnects,” IEEE Trans. Microwave Theory Tech., vol. 62, no. 11, pp. 2563–2573, 2014.
 [27] L. Chen, K. Sirenko, and H. Bagci, “Multiphysics analysis of plasmonic photomixers under periodic boundary conditions using discontinuous Galerkin time domain method,” in Prog. in Electromagn. Res. Symp. August 14, 2018 in Toyama, Japan, 2018.
 [28] J. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. NY, USA: Springer, 2008.
 [29] B. Cockburn and C.W. Shu, “The local discontinuous Galerkin method for timedependent convectiondiffusion systems,” SIAM J. Numer. Anal., vol. 35, no. 6, pp. 2440–2463, 1998.
 [30] C.W. Shu, “Discontinuous Galerkin methods for timedependent convection dominated problems: basics, recent developments and comparison with other methods,” in Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations, Cham: Springer International Publishing, 2016, pp. 371–399.
 [31] C.W. Shu and S. Osher, “Efficient implementation of essentially nonoscillatory shockcapturing schemes,” J. Comput. Phys., vol. 77, no. 2, pp. 439–471, 1988.
 [32] J. Hesthaven and T. Warburton, “Nodal highorder methods on unstructured grids: I. timedomain solution of Maxwell’s equations,” J. Comput. Phys., vol. 181, no. 1, pp. 186 – 221, 2002.
 [33] K. Sirenko, M. Liu, and H. Bagci, “Incorporation of exact boundary conditions into a discontinuous Galerkin finite element method for accurately solving 2d timedependent Maxwell equations,” IEEE Trans. Antennas Propag., vol. 61, no. 1, pp. 472–477, 2012.
 [34] P. Li, L. J. Jiang, and H. Bagci, “A resistive boundary condition enhanced DGTD scheme for the transient analysis of graphene,” IEEE Trans. Antennas Propag., vol. 63, no. 7, pp. 3065–3076, July 2015.
 [35] T. Lu, P. Zhang, and W. Cai, “Discontinuous Galerkin methods for dispersive and lossy Maxwell’s equations and PML boundary conditions,” J. Comput. Phys., vol. 200, no. 2, pp. 549–580, 2004.
 [36] S. D. Gedney, C. Luo, J. A. Roden, R. D. Crawford, B. Guernsey, J. A. Miller, T. Kramer, and E. W. Lucas, “The discontinuous Galerkin finiteelement timedomain method solution of Maxwell’s equations,” Appl. Comput. Electromagn. Soc. J., vol. 24, no. 2, p. 129, 2009.
 [37] J.P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys., vol. 114, no. 2, pp. 185–200, 1994.
 [38] L. Chen, M. B. Ozakin, and H. Bagci, “An efficient implementation of perfectly matched layers within a highorder discontinuous Galerkin time domain method,” December 1720, 2019 in Xiamen, China, 2019.
Comments
There are no comments yet.