## Abstract

Incorporation of plasmonic nanostructures in the design of photoconductive devices (PCDs) has significantly improved their optical-to-terahertz conversion efficiency. However, this improvement comes at the cost of increased complexity for the design and simulation of these devices. Indeed, accurate and efficient modeling of multiphysics processes and intricate device geometries of nanostructured PCDs is challenging due to the high computational cost resulting from multiple characteristic scales in time and space. In this work, a discontinuous Galerkin (DG)-based unit-cell scheme for efficient simulation of PCDs with periodic nanostructures is proposed. The scheme considers two physical stages of the device and models them using two coupled systems: a system of Poisson and drift-diffusion equations describing the nonequilibrium steady state, and a system of Maxwell and drift-diffusion equations describing the transient stage. A “potential-drop” boundary condition is enforced on the opposing boundaries of the unit cell to mimic the effect of the bias voltage. Periodic boundary conditions are used for carrier densities and electromagnetic fields. The unit-cell model described by these coupled equations and boundary conditions is discretized using DG methods. Numerical results demonstrate that the proposed DG-based unit-cell scheme has the same accuracy in predicting the THz photocurrent as the DG framework that takes into account the whole device, while it significantly reduces the computational cost.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

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–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–6]. This significant increase is attributed to several mechanisms: The optical electromagnetic (EM) field is enhanced due to plasmon [4] or Mie [5] resonances; nanostructured electrodes reduce the effective distance that the carriers have to travel [6,7]; and tailoring bias electric field using nanostructures can improve the device efficiency under low carrier densities [8]. Furthermore, for both conventional and nanostructured PCDs, physical effects resulting from the coupling between EM fields and carriers alter the device performance. For example, carrier screening causes saturation in the output THz current at high levels of optical pump power [9–11]. The interplay between all these different physical mechanisms makes the relevant device design very complicated. In this context, rigorous multiphysics simulation tools have become indispensable in understanding these physical mechanisms and their coupling and in enabling and/or accelerating the design process.

Due to their geometrically-intricate structure and the complicated EM wave/field interactions they support, simulation of nanostructured PCDs cannot be carried out using the methods that have been developed for conventional PCDs and rely on semi-analytical approximations/computations [12–18]. To this end, the finite element method (FEM) has been extensively used in recent years [19–22]. 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 effects cannot be modeled by this approach [19].

Recently, a multiphysics framework making use of discontinuous Galerkin (DG) methods has been proposed [23,24]. 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. This is simply because of the multiple space and time characteristic scales involved in the physical processes, e.g., the Debye length $\sim 10~\mathrm {nm}$, the optical wavelength $\sim 100~\mathrm {nm}$, and the device size $\sim 10~\mu \mathrm {m}$ [23].

One way to reduce this high computational cost is to make use of the nanostructure’s periodicity, i.e., model and discretize the multiphysics 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 [6,7,19,25,26]. However, since a PCD is in a non-equilibrium steady-state under a bias voltage [24,27], the boundary conditions required by the simulation of the carrier dynamics to be enforced on the unit-cell surfaces are not trivial. In [19], it is assumed that 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 computation 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 [19]. In addition, the nonlinear coupling is not fully considered in [19] since a frequency-domain FEM is used to compute the EM field distribution.

In this work, to reduce the high computational cost of nanostructured PCD simulations, a unit-cell scheme is proposed within the multiphysics DG framework developed in [23,24]. 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 practically 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 [23] while significantly reducing the computational cost and making it feasible to simulate practical 3D devices on desktop computers.

The rest of this paper is organized as follows. In Section 2.1, the unit-cell model and the relevant boundary conditions are introduced. Sections 2.2 and 2.3 describe the coupled systems of Poisson and stationary DD equations and time-dependent Maxwell and DD equations, respectively. Also, the DG schemes used for discretizing and solving these coupled systems of equations are provided in these sections. In Section 3, numerical results that validate the accuracy of the proposed scheme and demonstrate its computational benefits are provided. Finally, Section 4 concludes the paper and provides some remarks on the future research directions.

## 2. Mathematical model and discretization

#### 2.1 Unit-cell model

Figure 1 illustrates an example of a 3D nanostructured PCD. The device consists of a photoconductive region (LT-GaAs), a substrate layer (SI-GaAs), two electrodes that are deposited on the photoconductive region, 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 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 [27]. 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 [23,28].

To accurately capture these coupled physical processes that occur on the whole device using only 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. A critical observation that makes modeling the biased state using a unit cell possible is, in the electrostatic problem: The nanostructure generates only local variations in the potential [29,30] and on average the potential drops linearly between the two electrodes (as in the homogeneous case without the nanostructure) [31]. 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 are 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 $\varphi (x,y,z)$ is the electric potential, $w_x$ and $w_y$ are the widths of the unit cell in $x$ and $y$ directions, respectively, and $\varphi _{\mathrm {drop}}(y,z)$ is the potential-drop between the two faces of the unit cell perpendicular to the $x$ direction. In the rest of the text, (1) is termed as PDBC. Note that the PBC (2) is used along the $y$ direction because the potential does not drop along this direction (hence is periodic) [19]. The potential drop function ${\varphi _{\mathrm {drop}}}(y,z)$ in (1) is selected as described next. ${\varphi _{\mathrm {drop}}}(y,z)$ is position-dependent, e.g., the potential drop becomes smaller away from the electrodes in the $-z$ direction. Since the height of the device is much smaller than its width [1–3] and the electrodes extend to the whole width of the device, the potential drop is approximately the same for all $y$ and $z$ at a given value of $x$ (i.e., on a surface perpendicular to $x$ direction). Therefore, ${\varphi _{\mathrm {drop}}}(y,z)$ can be simplified to a single value ${\varphi _{\mathrm {drop}}}$. Additionally, as discussed before, on average the potential drops linearly between the two electrodes [31,32]. Consequently, $\varphi _{\mathrm {drop}}(y,z)$ can be estimated as: ${\varphi _{\mathrm {drop}}}(y,z) = {w_{x}}({{V_{\mathrm {bias}}}}\left / w_{\mathrm {sd}} \right . )$, where the ${V_{\mathrm {bias}}}$ is the bias voltage applied to the electrodes and $w_{\mathrm {sd}}$ is the distance between them.For stationary DD equations and time-dependent DD and Maxwell equations, PBCs are used:

where $U(x,y,z)$ 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 $z$ direction) are the same as those would be used in the full-device simulation [23,24].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 ${\varphi _{\mathrm {drop}}}(y,z)$ 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 the nanostructure is made of a dielectric material [1–5,33].

#### 2.2 Unit-cell Poisson-DD solver

The coupled system of Poisson and stationary DD equations is solved using the Gummel iteration method [24,27]. In each iteration, three partial differential equations (PDEs), namely the linearized Poisson equation and two DD equations are solved [24]. These equations, together with the boundary conditions described above, are discretized using DG methods.

The Poisson equation in the unit cell is expressed as a boundary value problem (BVP)

Discretizing $\Omega$ into a set of non-overlapping elements, DG approximates the unknowns with basis functions (the nodal basis function [34] is used in this work) in each element and applies Galerkin testing to the resulting equations [34–36]. This process yields a matrix system as

Here, $\bar \Phi$ and $\bar E$ are unknown vectors storing the basis expansion coefficients, ${{\bar M}^g}$ and ${\bar M}^{\mathbf {E}}$ are block diagonal mass matrices, $\bar G$ and $\bar D$ are block sparse matrices representing the gradient and divergence operators, respectively, and $\bar B^\varphi$ and $\bar B^{\mathbf {E}}$ have contributions from the tested force term and boundary conditions. Detailed expressions of these vectors and matrices can be found in [24].

The continuity of solutions across element interfaces is enforced through a uniquely defined numerical flux. For Poisson equation, the alternate flux [35]

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

where $U \in \{\varphi , (\varepsilon \mathbf {E})\}$.For PDBC, the element faces on the opposing surfaces are “connected” just like it is done for the PBC. An intermediate state is defined as $\varphi ^* = {\varphi (w_x/2,y,z)} + 0.5{\varphi _{\mathrm {drop}}}(y,z) = {\varphi (-w_x/2,y,z)} - 0.5{\varphi _{\mathrm {drop}}}(y,z)$, and the exterior variables in the numerical flux expressions are set as

The electron DD equations in the unit cell are expressed as a BVP [24]

DG discretization of (19)–(23) yields a matrix system as

The local Lax-Friedrichs flux [24]

The PBC (22) is enforced just like it is done in (13)–(14), with $U \in \{n, (d\mathbf {q})\}$, and for (21)

Solutions of the matrix systems (10) and (24) are computed using linear solvers at every iteration of the Gummel method. The steady-state solutions are obtained after the Gummel iteration converges and are used as inputs for the transient solver [23].

#### 2.3 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 their updated solutions into each other during the time integration [23]. Each set of equations is discretized using a time-domain DG (DGTD) method [34–36] that uses its own time-step size [23].

The time-dependent electron DD equations in the unit cell are expressed as an initial-boundary value problem (IBVP) [23]

Apart from the time dependency and the time derivative term on the left side of (27), 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 [23]

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 [37].

Likewise, Maxwell equations in the unit cell are expressed as the following IBVP

In (39)–(44), $F_k^{{\mathbf {E}},\nu }(t)$ and $F_k^{{\mathbf {H}},\nu }(t)$ are the $\nu$ components of $\hat {\mathbf {n}} \times [{\mathbf {H}_k}(t) - {\mathbf {H}}_k^*(t)]$ and $\hat {\mathbf {n}} \times [{\mathbf {E}_k}(t) - {\mathbf {E}}_k^*(t)]$, respectively, $\nu \in \{x,y,z\}$. The upwind flux [34] is used here:

## 3. 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 [23,24] (which has been validated against experimental data). The device is illustrated in Fig. 1. The thickness of the LT-GaAs and the SI-GaAs layers is $0.5~\mu$m. The nanograting has a periodicity of $0.18~\mu$m in $x$ and $y$ directions and its height is $0.12~\mu$m. The metal block is a truncated square pyramid with dimensions of $0.06~\mu$m and $0.1~\mu$m on its top and bottom, respectively. The semiconductor material parameters are same as those used in [23] and are provided in Table 1. The permittivity of LT-GaAs is modeled using the Lorentz dispersion relation

For the unit-cell simulation, the domain is shown in Fig. 1(b) and the boundary conditions are those explained in Section 2. For the full-device simulation, the height and width of the device (in $x$ and $y$ directions, respectively) are $3.1~\mu$m and $0.54~\mu$m. The nanograting is made of a $15 \times 3$ grid of unit cells, which is smaller than practical devices but is large enough for the purpose of validation in this work. For the DD equations, the Dirichlet boundary condition ${n_e} = (C + \sqrt {{C^2} + 4{n_i^2}} )/2$, ${n_h} = {n_i^2}/{n_e}$ is used on electrode/semiconductor interfaces and the homogeneous Robin boundary condition $\hat { \mathbf {n} } \cdot { \mathbf {J}_{e(h)}} = 0$ is used on semiconductor/insulator interfaces [27]. For Poisson equation, the Dirichlet boundary condition $\varphi = {V_{el}} + {V_T} \ln {(n_e^s/n_i)}$ is enforced on electrode surfaces and the homogeneous Neumann boundary condition $\hat {\mathbf {n}} \cdot \nabla \varphi = 0$ is used to truncate the computation domain [27]. Here, ${V_{el}}$ is the electric potential on the electrode, and ${V_{el}}=0$ for the negative electrode and ${V_{el}}=V_{\mathrm {bias}}$ for the positive one. For Maxwell equations, the computation domain is truncated with PMLs backed with PEC [44,45]. The simulation domains are discretized using tetrahedrons (Fig. 1). The minimum and the maximum edge lengths in the mesh are $10~\mathrm {nm}$ and $200~\mathrm {nm}$, respectively. The numbers of elements are ${1\,107\,866}$ and $13\,591$ in the full-device and unit-cell simulations, respectively. The tolerance of the Gummel iteration is $10^{-5}$ and the solution typically converges after $150$ iterations. The linear systems are solved using the GMRES iterative solver and an ILU preconditioner is reused throughout the Gummel iterations [24]. The physical duration of the transient stage is $8~\mathrm {ps}$ and the time-step sizes for the Maxwell and DD equations are $4\times {10^{-7}}~\mathrm {ps}$ and $2 \times {10^{-6}}~\mathrm {ps}$, respectively. These time-step sizes are chosen based on the Courant-Friedrichs-Lewy (CFL) conditions for the Maxwell and the DD equations [23], and are much smaller than the relaxation time ($10^{-3}~\mathrm {ps}$) of the carrier response on the PCD [27]. Table 2 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. Note that the unit-cell simulation can be carried out on a workstation, while the full-device simulation requires $\sim 1~\mathrm {TB}$ of memory and calls for a parallelized solver and a distributed-memory system. The unit-cell scheme’s steady state and transient stage simulations are $1800$ and $1500$ times faster.

The linear systems solved for the steady-state simulation are sparse however there is a trade-off between the ILU preconditioner sparsity and the number of GMRES iterations (denser preconditioner means smaller number of iterations). Therefore, the computational cost of the steady-state simulation is estimated to be between $O(N_{\mathrm {it}}N^2)$ and $O(N_{\mathrm {it}}N)$, where $N_{\mathrm {it}}$ is the number of GMRES iterations and $N$ is the number of unknowns. Note that $N_{\mathrm {it}}$ is larger for the full-device steady-state simulation since the mesh is more non-uniform (worsening the conditioning of the matrix systems). This short computational complexity analysis explains the large difference between the computation times required by the full-device and the unit-cell steady-state simulations.

The computation domain of the full-device transient-stage simulation contains many elements in the air background and in the PMLs. During time marching, only EM fields are updated for these elements while both EM fields and carrier densities are updated for the elements in the photoconductive region. Note that the exponentially decaying boundary layer of carrier densities require fine meshes on the boundaries of the photoconductive region. These fine meshes are required on almost all six faces of the photoconductive region for the full-device computation domain while in the unit-cell computation domain, fine meshes are only required on the top and bottom boundaries. These differences in the meshes used for full-device and unit-cell computation domains make the parallel load-balancing of the full-device transient-stage simulation significantly less efficient than that of the unit-cell transient-stage simulation. Therefore, even though the computational cost of the transient-stage simulation theoretically scales with $O(N)$ (linear in the number of unknowns $N$), the measured computation time comparison shows that the full-device transient-stage simulation is much slower than 80 times (the ratio of the numbers of unknowns in full-device and unit-cell computation domains).

Figure 2 shows the solutions obtained by the full-device simulation. Figure 2(a) shows the electric potential distribution on the interface between the nanostructures and the photoconductive region. 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. Figures 2(b) and (c) show the electric potential and electric field on several lines along the $x$ direction, respectively. The dash lines mark the positions of the unit-cell surfaces. Figure 2(b) shows that, although the potential distributions at different $z$ 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. Figure 2(c) shows the bias electric field is periodic as we expected.

Figures 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. Figure 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.

Figure 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 $5.9$.

Note that the unit-cell scheme assumes an infinitely large optical pump aperture (the optical pump is same for all unit cells). In the full-device simulation, the pump beam and the device have finite widths, which results in two effects: (1) The pump power near the center cells is higher than that near the boundary cells and (2) the optical field is scattered by the electrodes and the $x$/$y$ boundaries of the device. These effects lead to a small difference between the transient current densities obtained by the two methods.

## 4. Conclusion

The large scale of 3D nanostructured PCDs and various multiscale/multiphysics 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 the DG scheme that operates on the full device while it is more than $1500$ times faster.

The unit-cell scheme developed in this work can be used in the 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.

## Funding

King Abdullah University of Science and Technology (2016-CRG5-2953); Okawa Foundation Research Grant.

## Acknowledgments

The authors would like to thank the KAUST Supercomputing Laboratory (KSL) for providing the required computational resources.

## Disclosures

The authors declare no conflicts of interest.

## Data Availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## 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. **11**(1), 1600199 (2017). [CrossRef]

**2. **N. M. Burford and M. O. El-Shenawee, “Review of terahertz photoconductive antenna technology,” Opt. Eng. **56**(1), 010901 (2017). [CrossRef]

**3. **J.-H. Kang, D.-S. Kim, and M. Seo, “Terahertz wave interaction with metallic nanostructures,” Nanophotonics **7**(5), 763–793 (2018). [CrossRef]

**4. **N. T. Yardimci and M. Jarrahi, “Nanostructure-enhanced photoconductive terahertz emission and detection,” Small **14**(44), 1802437 (2018). [CrossRef]

**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. **59**(06), 1 (2019). [CrossRef]

**6. **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. **4**(5), 575–581 (2014). [CrossRef]

**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. **4**(1), 1622 (2013). [CrossRef]

**8. **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. **5**(1), 13817 (2015). [CrossRef]

**9. **J. T. Darrow, X. C. Zhang, D. H. Auston, and J. D. Morse, “Saturation properties of large-aperture photoconducting antennas,” IEEE J. Quantum Electron. **28**(6), 1607–1616 (1992). [CrossRef]

**10. **G. Rodriguez and A. J. Taylor, “Screening of the bias field in terahertz generation from photoconductors,” Opt. Lett. **21**(14), 1046–1048 (1996). [CrossRef]

**11. **M. Khorshidi and G. Dadashzadeh, “Plasmonic photoconductive antennas with rectangular and stepped rods: a theoretical analysis,” J. Opt. Soc. Am. B **33**(12), 2502–2511 (2016). [CrossRef]

**12. **M. Sirbu, S. B. P. Lepaul, and F. Aniel, “Coupling 3-D Maxwell’s and Boltzmann’s equations for analyzing a terahertz photoconductive switch,” IEEE Trans. Microw. Theory Tech. **53**(9), 2991–2998 (2005). [CrossRef]

**13. **M. Neshat, D. Saeedkia, L. Rezaee, and S. Safavi-Naeini, “A global approach for modeling and analysis of edge-coupled traveling-wave terahertz photoconductive sources,” IEEE Trans. Microw. Theory Tech. **58**(7), 1952–1966 (2010). [CrossRef]

**14. **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. **37**(1), 219–228 (2009). [CrossRef]

**15. **M. Khabiri, M. Neshat, and S. Safavi-Naeini, “Hybrid computational simulation and study of continuous wave terahertz photomixers,” IEEE Trans. THz Sci. Technol. **2**(6), 605–616 (2012). [CrossRef]

**16. **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. **61**(4), 1538–1546 (2013). [CrossRef]

**17. **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. **14**, 386–389 (2015). [CrossRef]

**18. **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. **4**(4), 490–500 (2014). [CrossRef]

**19. **N. Burford and M. El-Shenawee, “Computational modeling of plasmonic thin-film terahertz photoconductive antennas,” J. Opt. Soc. Am. B **33**(4), 748–759 (2016). [CrossRef]

**20. **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. **41**(2), 420–423 (2016). [CrossRef]

**21. **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 Adv. **7**(83), 53010–53017 (2017). [CrossRef]

**22. **N. M. Burford, M. J. Evans, and M. O. El-Shenawee, “Plasmonic nanodisk thin-film terahertz photoconductive antenna,” IEEE Trans. THz Sci. Technol. **8**(2), 237–247 (2018). [CrossRef]

**23. **L. Chen and H. Bagci, “Multiphysics simulation of plasmonic photoconductive devices using discontinuous Galerkin methods,” IEEE J. Multiscale Multiphys. Comput. Tech. **5**, 188–200 (2020). [CrossRef]

**24. **L. Chen and H. Bagci, “Steady-state simulation of semiconductor devices using discontinuous Galerkin methods,” IEEE Access **8**, 16203–16215 (2020). [CrossRef]

**25. **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 Lett. **12**(12), 6255–6259 (2012). [CrossRef]

**26. **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 Multiphys. Comput. Tech. **3**, 108–120 (2018). [CrossRef]

**27. **D. Vasileska, S. M. Goodnick, and G. Klimeck, * Computational Electronics: semiclassical and quantum device modeling and simulation* (CRC, 2010).

**28. **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.

**29. **L. Chen, M. Dong, and H. Bagci, “Modeling floating potential conductors using discontinuous Galerkin method,” IEEE Access **8**, 7531–7538 (2020). [CrossRef]

**30. **L. Chen, M. Dong, P. Li, and H. Bagci, “A hybridizable discontinuous Galerkin method for simulation of electrostatic problems with floating potential conductors,” Int. J. Numer. Model.: Electron. Networks, Device Fields p. e2894 (2020).

**31. **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. Microw. Theory Tech. **62**(11), 2563–2573 (2014). [CrossRef]

**32. **L. Chen and H. Bagci, “A unit-cell discontinuous Galerkin scheme for analyzing plasmonic photomixers,” in Proc. IEEE Int. Symp. Antennas Propag., (2019), pp. 1069–1070.

**33. **M. Fang, K. Niu, Z. Huang, W. E. I. Sha, X. Wu, T. Koschny, and C. M. Soukoulis, “Investigation of broadband terahertz generation from metasurface,” Opt. Express **26**(11), 14241–14250 (2018). [CrossRef]

**34. **J. Hesthaven and T. Warburton, * Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications* (Springer, 2008).

**35. **B. Cockburn and C.-W. Shu, “The local discontinuous Galerkin method for time-dependent convection-diffusion systems,” SIAM J. Numer. Anal. **35**(6), 2440–2463 (1998). [CrossRef]

**36. **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, G. R. Barrenechea, F. Brezzi, A. Cangiani, and E. H. Georgoulis, eds. (Springer, 2016), pp. 371–399.

**37. **C.-W. Shu and S. Osher, “Efficient implementation of essentially non-oscillatory shock-capturing schemes,” J. Comput. Phys. **77**(2), 439–471 (1988). [CrossRef]

**38. **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. **61**(1), 472–477 (2013). [CrossRef]

**39. **P. Li, Y. Shi, L. J. Jiang, and H. Bagci, “DGTD analysis of electromagnetic scattering from penetrable conductive objects with IBC,” IEEE Trans. Antennas Propag. **63**(12), 5296–5304 (2015). [CrossRef]

**40. **P. Li, L. J. Jiang, Y. J. Zhang, S. Xu, and H. Bagci, “An efficient mode-based domain decomposition hybrid 2-D/Q-2D finite-element time-domain method for power/ground plate-pair analysis,” IEEE Trans. Microw. Theory Tech. **66**(10), 4357–4366 (2018). [CrossRef]

**41. **T. Lu, P. Zhang, and W. Cai, “Discontinuous Galerkin methods for dispersive and lossy Maxwell’s equations and PML boundary conditions,” J. Comput. Phys. **200**(2), 549–580 (2004). [CrossRef]

**42. **J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. **114**(2), 185–200 (1994). [CrossRef]

**43. **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. **24**, 129 (2009).

**44. **L. Chen, M. B. Ozakin, S. Ahmed, and H. Bagci, “A memory-efficient implementation of perfectly matched layer with smoothly-varying coefficients in discontinuous Galerkin time-domain method,” IEEE Trans. Antennas Propag. p. 1 (2020).

**45. **P. Li, Y. Shi, L. J. Jiang, and H. Bagci, “A hybrid time-domain discontinuous Galerkin-boundary integral method for electromagnetic scattering analysis,” IEEE Trans. Antennas Propag. **62**(5), 2841–2846 (2014). [CrossRef]