Photoconductive devices (PCDs) are promising candidates for terahertz (THz) source generation and signal detection because they are compact and frequency-stable, and capable of operating at room temperature with low optical input power levels [1, 2, 3, 4, 5]. However, the low optical-to-THz 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  or Mie  resonances. Nanostructured electrodes reduce the effective distance that the carriers have to travel [7, 6]. Furthermore, recent studies show that the nanostructure-tailored 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 geometrically-intricate 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 semi-analytical 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 frequency-domain 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. 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 drift-diffusion (DD) equations describing the nonequilibrium steady state of the PCD and a coupled system of time-dependent 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 DG-based framework provides higher flexibility and higher-order accuracy in both space discretization and time integration and is more robust in modeling nonlinear coupling mechanisms compared to the FEM-based schemes, it is still computationally demanding especially for practical three-dimensional (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 multi-physical 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 well-studied [21, 7, 6, 15, 22]. However, since a PCD is in a non-equilibrium steady-state under a bias voltage [20, 23], the boundary conditions required by the simulation of the carrier dynamics to be enforced on the unit-cell surfaces are not trivial. In , it is assumed there is no potential-drop 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 . In addition, the nonlinear coupling is not fully considered in  since a frequency-domain FEM is used to compute the EM field distribution.
In this work, to reduce the high computational cost of simulating nanostructured PCDs, a unit-cell scheme is proposed within the multiphysics DG framework developed in [19, 20]. For Poisson equation, a “potential-drop” boundary condition (PDBC) is enforced on the opposing surfaces of the unit-cell (along the direction of the bias electric field). For carriers and EM fields, PBCs are enforced on the unit-cell surfaces. All boundary conditions are “weakly” enforced using the numerical flux of the DG scheme. Numerical results demonstrate that the proposed DG-based unit-cell 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 , 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 II-A, the unit-cell model and the relevant boundary conditions are introduced. In Sections II-B and II-C, the coupled systems of Poisson and stationary DD equations and time-dependent 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 Unit-cell Model and Numerical Scheme
Ii-a The unit-cell model
Fig. 1 illustrates an example of a 3D nanostructured PCD. The device consists of a photoconductive layer (LT-GaAs), a substrate layer (SI-GaAs), 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 .
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 re-distributed carriers in turn affects the electric potential distribution. The device reaches a non-equilibrium steady-state mathematically described by a coupled system of Poisson and stationary DD equations . 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 time-dependent 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 unit-cell 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  and on average the potential drops linearly between the two electrodes (as in the homogeneous case without the nanostructure) . 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 steady-state 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, time-dependent Maxwell equations, and time-dependent DD equations.
The boundary conditions discussed above are mathematically described as follows. For Poisson equation, the boundary conditions are
where is the electric potential, and are the widths of the unit cell in - and -directions, respectively, and is the potential-drop 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) . The potential drop function in (1) is selected as described next. is position-dependent, 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 time-dependent DD and Maxwell equations, PBCs are used:
where represents the steady-state 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 full-device simulation [20, 19].
Two comments about the unit-cell 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 unit-cell model to hold; it is also applicable when nanostructure is made of dielectric material .
Ii-B The unit-cell Poisson-DD 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 
. In each iteration, one needs to solve three partial differential equations (PDEs), namely the linearized Poisson equation and two DD equations. 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)
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) .
Discretizing into a set of non-overlapping elements, DG approximates the unknowns with basis functions (the nodal basis function  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
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 .
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 
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 . 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
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
The electron DD equations in the unit cell are expressed as the following BVP 
where and are the unknowns to be solved for, and , and are known coefficients within the Gummel method . The homogeneous Robin boundary condition is used in (23) . The BVP for holes only differs by the sign in front of the drift term.
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 .
The local Lax-Friedrichs flux  is used for the drift term
Ii-C The unit-cell Maxwell-DD solver
The coupled system of time-dependent 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 . Each set of equations is discretized using a time-domain DG (DGTD) method [29, 28, 30] that uses its own time-step size .
The time-dependent electron DD equations in the unit cell are expressed as an initial-boundary value problem (IBVP) 
where and are the unknowns to be solved for, and are known coefficients, is the net-recombination, 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 (II-C), the system has the same form as (19)-(23). Using the same spatial discretization as that used for (19)-(23), one can obtain the semi-discrete form 
where and unknown vectors storing time-dependent basis expansion coefficients, and are vectors constructed from the net-recombination and boundary conditions (other than PBC), , /, /, and / are the elemental mass, convection, gradient, and divergence matrices, respectively .
The numerical fluxes and boundary conditions are the same as those used in the stationary DD equations. The semi-discretized system (32)-(33) is integrated in time using an explicit third-order total-variation-diminishing Runge-Kutta scheme .
Likewise, Maxwell equations in the unit cell are expressed as the following IBVP
where and are unknown vectors storing time-dependent 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].
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 semi-discrete system (39)-(44) is integrated in time using a five-stage fourth-order Runge-Kutta method .
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 LT-GaAs and the SI-GaAs layers is m. The nanograting has a periodicity of m in x- and y-directions 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  and are provided in Table I. The DD equations are solved only within the LT-GaAs layer while Poisson equation and the time-dependent 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 full-device simulation.
|LT-GaAs||, , ,|
For the unit-cell simulation, the domain is shown in Fig. 1 (b) and the boundary conditions are those explained in Section II. For the full-device 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 . 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 . For Maxwell equations, the computation domain is truncated with PMLs  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 full-device and unit-cell 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 . The physical duration of the transient stage is and the time-step sizes for the Maxwell and DD equations are and , respectively . Table II provides the computational times (measured in “core-hour”) required by the unit-cell and full-device simulations to complete the steady state and the transient stage. The unit-cell scheme’s steady state and transient stage simulations are and times faster.
|steady state||transient stage|
Tested on Shaheen II (https://www.hpc.kaust.edu.sa/), a Cray XC40 system consists of 6174 nodes based on Intel Xeon (R) E5-2698 v3.
It should be noted here that the unit-cell simulation can be carried out on a workstation, while the full-device simulation requires of memory and calls for a parallelized solver and a distributed-memory system.
Fig. 2 shows the solutions obtained by the full-device 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 unit-cell 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 unit-cell surfaces. Fig. 2 (c) shows the bias electric field is periodic as we expected.
Figs. 3 (a) and (c) show the steady-state electric potential and electron density calculated from the unit-cell scheme, respectively. For comparison, Figs. 3 (b) and (d) show those calculated using the full-device, 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 unit-cell scheme, the potential variation inside the unit cell is same as that obtained from the full-device 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 full-device. And, Fig. 3 (c) shows that the solution obtained from the unit-cell scheme is same as that obtained from the full-device simulation.
Fig. 4 compares the transient current densities obtained from the unit-cell and full-device 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 unit-cell scheme assumes an infinitely large pumping aperture (the optical pumping is same for all unit cells). In the full-device 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.
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 unit-cell 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 time-dependent 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 unit-cell surfaces for Poisson equation, while PBCs are used for stationary and time-dependent DD equations and time-dependent Maxwell equations. These BCs are weakly enforced using the numerical flux of the DG methods.
Numerical results demonstrate that the proposed unit-cell DG scheme maintains the accuracy of a DG scheme that operates on the full device while it is more than times faster. The unit-cell 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.
This research is supported by the King Abdullah University of Science and Technology (KAUST) Office of Sponsored Research (OSR) under Award No 2016-CRG5-2953. The authors would like to thank the KAUST Supercomputing Laboratory (KSL) for providing the required computational resources.
-  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.
-  N. M. Burford and M. O. El-Shenawee, “Review of terahertz photoconductive antenna technology,” Opt. Eng., vol. 56, no. 1, p. 010901, 2017.
-  J.-H. Kang, D.-S. Kim, and M. Seo, “Terahertz wave interaction with metallic nanostructures,” Nanophotonics, vol. 7, no. 5, pp. 763–793, 2018.
-  N. T. Yardimci and M. Jarrahi, “Nanostructure-enhanced photoconductive terahertz emission and detection,” Small, vol. 14, no. 44, p. 1802437, 2018.
-  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.
-  S.-H. Yang, M. R. Hashemi, C. W. Berry, and M. Jarrahi, “7.5% optical-to-terahertz conversion efficiency offered by photoconductive emitters with three-dimensional plasmonic contact electrodes,” IEEE Trans. THz Sci. Technol., vol. 4, no. 5, pp. 575–581, 2014.
-  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.
-  K. Moon, I. Lee, J.-H. Shin, E. S. Lee, K. N, S. Han, and K. H. Park, “Bias field tailored plasmonic nano-electrode for high-power terahertz photonic devices,” Sci. Rep., vol. 5, no. 13817, Nov 2012.
-  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. continuous-wave radiation,” Appl. Phys. Lett., vol. 109, no. 7, p. 071105, 2016.
-  P. Kirawanich, S. J. Yakura, and N. E. Islam, “Study of high-power wideband terahertz-pulse generation using integrated high-speed photoconductive semiconductor switches,” IEEE Trans. Plasma Sci., vol. 37, no. 1, pp. 219–228, 2008.
-  M. Khabiri, M. Neshat, and S. Safavi-Naeini, “Hybrid computational simulation and study of continuous wave terahertz photomixers,” IEEE Trans. THz Sci. Technol., vol. 2, no. 6, pp. 605–616, 2012.
-  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.
-  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.
-  E. Moreno, M. F. Pantoja, S. G. Garcia, A. R. Bretones, and R. G. Martin, “Time-domain numerical modeling of THz photoconductive antennas,” IEEE Trans. THz Sci. Technol., vol. 4, no. 4, pp. 490–500, 2014.
-  N. Burford and M. El-Shenawee, “Computational modeling of plasmonic thin-film terahertz photoconductive antennas,” J. Opt. Soc. Am. B, vol. 33, no. 4, pp. 748–759, 2016.
-  M. J. Mohammad-Zamani, M. Neshat, and M. K. Moravvej-Farshi, “Nanoslit cavity plasmonic modes and built-in fields enhance the cw THz radiation in an unbiased antennaless photomixers array,” Opt. Lett., vol. 41, no. 2, pp. 420–423, 2016.
-  M. Bashirpour, S. Ghorbani, M. Kolahdouz, M. Neshat, M. Masnadi-Shirazi, 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.
-  N. M. Burford, M. J. Evans, and M. O. El-Shenawee, “Plasmonic nanodisk thin-film terahertz photoconductive antenna,” IEEE Trans. THz Sci. Technol., vol. 8, no. 2, pp. 237–247, 2017.
-  L. Chen and H. Bagci, “Multiphysics modeling of plasmonic photoconductive devices using discontinuous Galerkin methods,” arXiv preprint arXiv:1912.03639, 2019.
-  L. Chen and H. Bagci, “Steady-state simulation of semiconductor devices using discontinuous Galerkin methods,” IEEE Access, vol. 8, pp. 16 203–16 215, 2020.
-  B. Heshmat, H. Pahlevaninezhad, Y. Pang, M. Masnadi-Shirazi, 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.
-  K. Sirenko, Y. Sirenko, and H. Bagci, “Exact absorbing boundary conditions for periodic three-dimensional structures: Derivation and implementation in discontinuous Galerkin time-domain method,” IEEE J. Multiscale and Multiphys. Comput. Techn., vol. 3, pp. 108–120, 2018.
-  D. Vasileska, S. M. Goodnick, and G. Klimeck, Computational Electronics: semiclassical and quantum device modeling and simulation. CRC press, 2010.
-  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.
-  L. Chen, M. Dong, and H. Bagci, “Modeling floating potential conductors using discontinuous Galerkin method,” IEEE Access, vol. 8, pp. 7531–7538, 2020.
-  A. Menshov and V. I. Okhmatovski, “Surface–volume–surface electric field integral equation for magneto-quasi-static analysis of complex 3-D interconnects,” IEEE Trans. Microwave Theory Tech., vol. 62, no. 11, pp. 2563–2573, 2014.
-  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 1-4, 2018 in Toyama, Japan, 2018.
-  J. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. NY, USA: Springer, 2008.
-  B. Cockburn and C.-W. Shu, “The local discontinuous Galerkin method for time-dependent convection-diffusion systems,” SIAM J. Numer. Anal., vol. 35, no. 6, pp. 2440–2463, 1998.
-  C.-W. Shu, “Discontinuous Galerkin methods for time-dependent 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.
-  C.-W. Shu and S. Osher, “Efficient implementation of essentially non-oscillatory shock-capturing schemes,” J. Comput. Phys., vol. 77, no. 2, pp. 439–471, 1988.
-  J. Hesthaven and T. Warburton, “Nodal high-order methods on unstructured grids: I. time-domain solution of Maxwell’s equations,” J. Comput. Phys., vol. 181, no. 1, pp. 186 – 221, 2002.
-  K. Sirenko, M. Liu, and H. Bagci, “Incorporation of exact boundary conditions into a discontinuous Galerkin finite element method for accurately solving 2d time-dependent Maxwell equations,” IEEE Trans. Antennas Propag., vol. 61, no. 1, pp. 472–477, 2012.
-  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.
-  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.
-  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 finite-element time-domain method solution of Maxwell’s equations,” Appl. Comput. Electromagn. Soc. J., vol. 24, no. 2, p. 129, 2009.
-  J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys., vol. 114, no. 2, pp. 185–200, 1994.
-  L. Chen, M. B. Ozakin, and H. Bagci, “An efficient implementation of perfectly matched layers within a high-order discontinuous Galerkin time domain method,” December 17-20, 2019 in Xiamen, China, 2019.