## I Introduction

In recent years, electromagnetic (EM) systems operating at terahertz (THz) frequencies have found applications in various fields ranging from wireless communication to non-destructive testing [10, 46]. One of the fundamental challenges in developing these systems is implementing efficient THz sources and detectors. Plasmonic photoconductive antennas (PCAs) have become one of the most promising candidates to address this challenge because of their dramatically increased optical-to-THz conversion efficiency [20, 4, 16, 44, 41]. While the experimental research on fabrication of plasmonic PCAs has progressed remarkably in the recent years, there is still room for significant improvement in formulating and implementing numerical methods for rigorous and accurate simulation of this type of devices.

Challenges in simulation of PCAs stem from the fact that their operation relies on strongly coupled and simultaneously occurring multiple physical processes. A PCA incorporates a photoconductive semiconductor device that is biased under an external applied voltage. Upon excitation by an optical EM wave, this semiconductor device generates carriers and produces THz currents [20, 4]. Therefore, the operation of the PCA can be analyzed in two stages: (1) a stationary state corresponding to the non-equilibrium state of the semiconductor device under the bias voltage, and (2) a transient stage describing the dynamic processes after the optical EM wave impinges on the device. The stationary state is a result of the interaction between the static electric field and the carriers and is mathematically described by a coupled system of the Poisson equation and the stationary drift-diffusion (DD) equations. At the transient stage, several coupled processes happen simultaneously: EM field/wave interactions and carrier generation, recombination, drift and diffusion. This stage is mathematically described by a coupled system of the time-dependent Maxwell and DD equations. Because of the strong nonlinear coupling between these two sets of equations, accurate modeling of the transient stage calls for a multiphysics solver operating in time domain.

Characteristic space and time scales of the various physical interactions listed above differ by several orders of magnitude. Space scales change from (Debye length) to m (geometrical dimensions of the device). Furthermore, plasmonic nanostructures and localized plasmon modes generated on them make the space scale even more complicated. Time scales are determined by the carrier dynamics (including the advection and diffusion motions) and the EM wave interactions. These two scales typically differ by orders of magnitude. This multiscale nature of the physical interactions involved in the PCA operation requires highly flexible spatial and temporal discretization techniques.

Several numerical methods have been developed for simulating PCAs. For conventional PCAs (PCAs without plasmonic metallic nanostructures), the finite difference time domain (FDTD) method has been extensively used [19, 27, 28, 25]. These schemes solve the time-dependent Maxwell and DD equations to account for carrier dynamics and THz EM wave interactions on the antenna, respectively, while the carrier generation due to the optical EM wave excitation is taken into account using a closed-form analytical expression. It has been shown in [27, 28, 26] that the results obtained using this approach match experimental results well for conventional PCAs. For a plasmonic PCA, the nanostructures introduced on top of or inside the semiconductor region produce highly localized and strong EM fields resulting in a significantly increased carrier generation and an enhanced optical-to-THz conversion efficiency. However, nanostructures also result in strong scattering of the optical EM wave, making the analytical expression used for carrier generation rate inaccurate. In addition, the geometrically complicated metallic nanostructures and the resulting spatially fast varying plasmon modes cannot be discretized and accounted for accurately using an orthogonal FDTD grid (unless an extremely fine grid, which would prohibitively increase the computational requirements, is used).

Methods other than FDTD are also developed for simulating conventional PCAs [29, 17, 18, 45]. In [29, 17]

, the frequency-domain optical EM field is used to compute the (time-modulated) carrier generation and the transient stage is executed using a coupled Poisson-DD solver. In

[18], an equivalent circuit model is used to predict the optical-to-THz conversion efficiency. In [45], the semiconductor device is treated as a circuit attached to a THz antenna, on which the EM interactions are accounted for using a discontinuous Galerkin time domain (DGTD) scheme. Even though, these approaches have been shown to be efficient for the simulation of conventional PCAs, they cannot be used for simulating plasmonic PCAs since the approximations involved in their formulation and implementation do not allow for accurate modeling of the optical EM wave/field interactions on the nanostructures.To simulate plasmonic PCAs, numerical frameworks that rely on the finite element method (FEM) have been developed [3, 1]

. FEM, thanks to its capability to operate on unstructured meshes, can accurately account for the EM field interactions on the metallic nanostructures. These frameworks compute the space distribution of the carrier generation using the EM field distribution obtained by a frequency-domain FEM. The time dependency of the carrier generation is assumed to be a pulse with a closed-form analytical expression. Then the carrier generation (which is now assumed to be known in space and time) is fed into the time domain simulation of carrier dynamics. Even though the EM field interactions on metallic nanostructures are accurately accounted for with this approach, the nonlinear two-way coupling between the EM interactions and carrier dynamics is not fully considered. Furthermore, the space-charge screening effect, i.e., the fact that separation of electrons and holes induces a polarization vector that cancels the electric field, which results from the time-dependent behavior of the carrier distributions, is not captured by this approach.

In this work, we propose a multiphysics framework that makes use of discontinuous Galerkin (DG) methods to simulate plasmonic PCAs [15, 34, 36]. This framework does not suffer from the shortcomings of the previously-developed methods as briefly described above. Following the two stages of the PCA operation, the semiconductor device is numerically modeled also in two stages: The stationary state is described by a coupled system of the (nonlinear) Poisson equation and the bipolar DD equations. This coupled system is solved iteratively using the Gummel method and the linear systems at each iteration are discretized using (stationary) DG schemes [5]. The transient stage is described by a nonlinearly coupled system of the time-dependent Maxwell and DD equations. This coupled system is discretized using a DGTD scheme, where the nonlinear coupling between the two sets of equations is fully taken into account during the time integration. In this numerical framework, we prefer to use DG-based discretization because (i) it provides higher flexibility in meshing (just like FEM), (ii) it allows for explicit time integration (just like FDTD), and (iii) it permits easy implementation of high-order (spatial) basis functions and high-order time integration schemes. Explicit time integration, which is enabled by the use of DG-based discretization, equips the resulting multiphysics framework with several advantages: (i) The nonlinear coupling between the Maxwell and DD equations can be easily and fully implemented. (ii) Each set of equations can use their own time-step size. Note that using a global time-step size would unnecessarily increase the computational requirements due to the large difference between the characteristic time scales of the Maxwell and DD equations. (iii) Since the explicit time integration does not require any matrix inversion, it can be parallelized easily and executed efficiently on distributed-memory computer clusters. In addition, DG schemes allow for non-conformal meshes, adaptive hp-refinement, and local time stepping. These properties make DG methods very suitable for multiphysics and multiscale simulations [42, 14].

The rest of the paper is organized as follows. Section II-A describes the physical processes involved in the operation of the plasmonic PCAs and presents the systems of equations that are used to mathematically model them. Section II-B details the numerical schemes that are used for discretizing these systems and solving the resulting matrix equations. Some comments about the formulation and implementation are presented in Section II-C. Section III demonstrates the accuracy, the efficiency, and the applicability of the proposed numerical framework via numerical results. Finally, Section IV summarizes the work described in the paper and presents several future research directions.

## Ii Formulation

### Ii-a Physical and Mathematical Description

The semiconductor most commonly used in PCA designs is the low-temperature grown gallium arsenide (LT-GaAs). LT-GaAs has an electron trapping time of less than and a direct optical band gap of 1.42eV. It can “absorb” optical EM energy and generate carriers with a lifetime shorter than . Two conductive electrodes are deposited on the semiconductor substrate. A bias voltage, which is applied to these electrodes, drives the carriers towards them, resulting in a THz current. Consequently, the electrodes act as a current source attached to an antenna. For a conventional PCA, the optical-to-THz conversion efficiency is limited by the amount of absorbed optical EM energy, which is usually small due to the high refractive index of the semiconductor layer.

Plasmonic PCAs utilize metallic nanostructures to increase the optical-to-THz conversion efficiency. These nanostructures are designed to introduce plasmon modes with resonances around the operating frequency of the optical EM wave excitation. Highly localized and strong EM fields associated with these modes significantly increase the carrier generation inside the semiconductor layer. In addition to introducing the plasmon modes, the nanostructures, which are either placed between the electrodes [16] or etched on to them [2], change the static electric field distribution. For the latter case, they act as a part of the electrode they are attached to and effectively reduce the distance between the carriers and the electrodes. Both of these effects significantly change the current [24, 44]. Therefore, the nanostructures are designed not only to generate strong localized EM fields, but also to maximize the overall current response. We also should mention here that the impedance mismatch between the semiconductor device (acting as the current source) and the THz antenna is not as critical as the current generation rate in terms of PCA’s performance [2, 3].

Mathematical modeling of a PCA has to account for carrier dynamics, EM interactions, and their nonlinear coupling. While Maxwell equations are used to describe EM interactions, various models ranging from semi-classical to full-quantum approaches have been developed to describe carrier dynamics. For PCAs, the semi-classical DD model can be used accurately since the device size (m) is far larger than the mean free path of the carriers (m) [8, 39, 28]. In this work, we prefer to use the bipolar DD equations since they can account for the difference in the physical parameters associated with electrons and holes and allow for modeling the space-charge screening effect. Thus, carrier dynamics and EM interactions that happen on a PCA are mathematically described by the following coupled system of equations

(1) | |||

(2) | |||

(3) | |||

(4) |

Here, is the location vector, subscript represents the carrier type and hereinafter the upper and lower signs should be selected for electron () and hole (), respectively, and are the electric and magnetic field intensities, and are the dielectric permittivity and permeability, and are the electron and hole densities, and are the current densities due to electron and hole movement, and are the recombination and generation rates, and are the field-dependent electron and hole mobilities, and are the electron and hole diffusion coefficients, respectively.

The bias voltage is applied to the electrodes continuously throughout the operation of the device. Before the optical EM wave excitation is turned on, the device is under a non-equilibrium stationary state [5]. When the EM wave impinges on the device, carriers are generated. However, the density of these carriers is several orders of magnitude smaller than the doping concentration and the stationary state of the device is assumed to be only weakly perturbed by the optical EM wave excitation. Therefore, field intensities and carrier and current densities are separated into stationary and transient components as , , , , respectively. Here, the superscript “s” and “t” stands for stationary and transient components, respectively.

Similarly, the recombination rate in (3) is decomposed into stationary and transient components as [27]. In this work, we use the trap assisted recombination described by the Shockley-Read-Hall model [8, 39, 5]. The generation rate is defined as [8]

(5) |

Here, is the wavelength of the optical EM wave excitation, is the photon absorption coefficient, is the quantum yield, is the Planck constant, and is the light speed. We should emphasize here that depends only on the transient EM fields and that are generated on the PCA due to the optical EM wave excitation and include the effect of all transient EM interactions.

We should also mention here that one has to account for the field-dependency and velocity-saturation in the mobility model used in (4) for a more accurate prediction of the generated THz current (mobility affects the carrier drift velocity, which in return has a big impact on the current). Here, we use the parallel-field dependent mobility model reported in [28]. Having said that, since the static electric field component is at least two orders of magnitude stronger than its transient counterpart, , it is assumed that mobility is only a function of , i.e., , .

Following the discussion above, the coupled of system of time-dependent Maxwell and DD equations in (1)-(4) is effectively decomposed into two coupled sets of equations. The stationary components of the variables (the ones with superscript “s”) satisfy a coupled system of Poisson and stationary DD equations. This is the first equation system and it reads

(6) | |||

(7) | |||

(8) |

Eliminating the stationary components in (1)-(4) yields a reduced coupled system of Maxwell and DD equations. This is the second equation system and it reads

(9) | |||

(10) | |||

(11) | |||

(12) |

Equations (9)-(II-A) are a strongly-coupled nonlinear system. In (11), leads to an exponential increase in carrier densities. In (II-A), carriers are driven by both and . The carrier motion produces the free current densities and , which contributes to the EM wave/field interactions described by (9)-(10). In return, these EM wave/field interactions, change . We should also note here that is also a nonlinear function, which becomes significant when carrier densities are high, balancing the carrier generation.

### Ii-B The Multiphysics DG Solver

A complete simulator for numerical characterization of PCAs consists of a stationary scheme that solves the Poisson-DD system (6)-(8) and a time domain scheme that solves the Maxwell-DD system (9)-(II-A). The stationary solutions are used as inputs to the time domain scheme. The stationary scheme uses the Gummel iteration method to take the nonlinearity into account. The linearized set of equations that needs to be solved at every iteration of this method is discretized using a stationary DG scheme. For details, we refer the reader to [5]. In the rest of this paper, we focus on the time domain scheme proposed to solve the Maxwell-DD system (9)-(II-A).

#### Ii-B1 Nonlinear Coupling

To integrate the time-dependent Maxwell-DD system in time, we use an explicit scheme. To account for the difference in the characteristic time scales of the EM interactions and the carrier dynamics in an efficient manner, the Maxwell and DD equations are updated using two different time-step sizes. More specifically, since the carrier response is much slower than the time variation of the EM waves/fields, the DD equations are updated using a larger time-step size. The nonlinear coupling between the two equation sets is accounted for by alternately feeding the updated solutions into each other.

This time marching scheme is shown in Fig. 1, where the time-step size for the DD equations is assumed to be twice the step size for the Maxwell equations for illustration purposes. Let us suppose that two equation systems are updated separately using two different time integration schemes (to be discussed in the next section) and denote the time steps and step sizes of the Maxwell and DD systems with {, } and {, }, respectively. We should note here that, in what follows, subscripts and mean that the variables they are attached to are computed/sampled/updated at and , respectively.

Let us assume that the time steps of the two systems are synchronized at time . First, and are used to compute the generation rate using (5). is used in (11) to update the carrier densities (from step to ). Then, are used in (II-A) to compute the current densities and . and are then used in (10) to update the Maxwell equations at steps and to produce {, } and {, }, respectively. The time steps of the two systems match again at time and the process described here is repeated.

#### Ii-B2 Discretization

The time-dependent DD and Maxwell equations [(11)-(II-A) and (9)-(10), respectively] are discretized using DG schemes in time-domain. We start the description of discretization with the DD equations (11)-(II-A). Since the electron and hole DD equations only differ by the sign in front of the drift term, we only discuss the electron DD equation. Also we note that we drop the subscript “e” (meaning electron) and the superscript “t” (meaning transient) from the variables to simplify the notation. Since is provided as an input to the time-domain simulation, variables that are functions of it are assumed to change with only. Furthermore, as the three drift terms in (II-A) are treated in the same way, for brevity, we combine them under one term that is denoted by . Under those notation simplifications, the electron DD equations (11)-(II-A) are expressed as the following initial-boundary value problem (IBVP)

(13) | ||||

(14) | ||||

(15) | ||||

(16) |

Here, is an auxiliary variable introduced to reduce the order of the spatial derivative in the diffusion term, , and represent the surfaces where Dirichlet and Robin boundary conditions are enforced and and are the coefficients associated with these boundary conditions, respectively, and denotes the outward normal vector . For the problems considered in this work, represents the electrode/semiconductor interfaces and ; and represents the semiconductor/insulator interfaces and indicating no carrier spills out those interfaces [32].

To facilitate the numerical solution of the IBVP described by (II-B2)-(16), is discretized into non-overlapping tetrahedrons. The volumetric support of each of these discretization elements is represented by , . Let and denote the surface of and the outward unit vector normal to , respectively. Testing (II-B2) and (14) with Lagrange polynomials, , [15], on element and applying the divergence theorem yield the weak form

(17) | ||||

(18) |

Here,

is the number of interpolation nodes,

is the order of the Lagrange polynomials andis used for identifying the components of the vectors in the Cartesian coordinate system. We note here

and denote the local solutions on element and the global solutions on are the direct sum of the local solutions.In (17) and (18), , , and are numerical fluxes “connecting” element to its neighboring elements. The variables involved in the definition of the numerical fluxes reside on the interfaces between elements and their explicit dependencies on , , and are dropped to simplify the notation. For the diffusion term, the local DG (LDG) alternate flux [9] is used for and . For the drift term, the local Lax-Friedrichs flux [15] is used for . On the surfaces of the element, they are defined as

Here, “average” and “jump” operators are respectively defined as and , where is a scalar or a vector variable. Superscripts “-” and “+” refer to variables defined in element and in its neighboring element, respectively. The vector determines the upwind direction of and . In LDG, it is essential to choose opposite directions for and (note the sign difference in the definitions above), while the precise direction of each variable is not important [9, 15, 34]. In this work, we choose on each element surface. The local Lax-Friedrichs, with [15], mimics the path of the information propagation. On , , and , and on , and . We note that and are not assigned independently on .

and are expanded using Lagrange polynomials [15]

(19) | |||

(20) |

where , , denote the locations of the interpolation nodes, and , , , are the unknown coefficients to be solved for. Substituting (19) and (20) into (17) and (18), we obtain the semi-discretized form

(21) | ||||

(22) |

where the unknown vectors are defined as and , with , .

In (21)-(22), and are mass matrices. is a block diagonal matrix with blocks . is defined as

is a diagonal matrix with entries , where and , , . We note that is assumed isotropic and constant in each element. Matrices and , and and correspond to the gradient and the divergence operators, respectively. For the LDG flux, and . is a matrix and it has contributions from the volume and surface integral terms on the right hand side of (18): , , , where

and

Here, denotes the interface connecting element and and selects the interpolation nodes on the interface

Matrix corresponds to the surface integral term in (18), which involves the unknowns from neighboring elements of element : , where

Similarly, matrix has contributions from the fourth term (the volume integral) and the fifth term (the surface integral) on the right hand side of (17): , where

and

Matrix corresponds to the last surface integral term in (17), which involves the unknowns from neighboring elements of the element :

Matrices and are contributed from the force term and boundary conditions (where element does not exist) and are expressed as

To integrate (21)-(22) in time, an explicit third-order total-variation-diminishing (TVD) Runge-Kutta method [35] is used. The high-order accuracy of this scheme matches that of the spatial discretization. With initial value , time samples of the unknown vector are obtained step by step in time.

The time-dependent Maxwell equations (9) and (10) are discretized using the nodal DG method [15]. First the simulation domain is divided into non-overlapping tetrahedrons with volumetric support , . Just like before, and denote the surface of and the outward unit vector normal to , respectively. Testing (9) and (10) with the Lagrange polynomials , , on element and applying the divergence theorem twice yield the strong form [15]

(23) | ||||

(24) |

Here, and are the local solutions on element , and and are the upwind numerical fluxes “connecting” element to its neighboring elements. They are expressed as [15]

where the average and the jump operators are the same as those defined before, and are the wave impedance and wave admittance, respectively. Same as before, the variables are defined on element surfaces and their explicit dependencies on , , and are dropped to simplify the notation. On boundaries where there are no neighboring elements, the numerical fluxes are assigned using the relevant boundary conditions. That is, on perfect electric conductor (PEC) surfaces, , , and, for absorbing boundary conditions (ABCs), , [15]. We note here that, since the performance of ABCs degrades rapidly with the angle of incidence, they are used to terminate perfectly matched layers (PMLs) wrapping around the simulation domain. Details of the implementation of PMLs can be found in [12].

and are expanded using Lagrange polynomials [15]

(25) | |||

(26) |

where , , denote the locations of the interpolation nodes, denotes the components of vectors in the Cartesian coordinate system, and and are the unknown coefficients to be solved for. Substituting (25) and (26) into (23) and (24) yields [15]

(27) | |||

(28) | |||

(29) | |||

(30) | |||

(31) | |||

(32) |

where and are unknown vectors, is the current density vector, , and and are the components of and , respectively.

Comments

There are no comments yet.