In recent years, terahertz (THz) devices have become indispensable components of electromagnetic (EM) systems in applications ranging from wireless communication to non-destructive testing. The main obstacle in the way of widespread industrial use of THz technologies is the difficulty of the implementing compact, efficient, frequency-stable and tunable THz source generators capable of operating at room temperatures. Plasmonic photoconductive devices (PCDs) have become one of the most promising candidates for THz source generation because of their dramatically increased optical-to-THz conversion efficiency and polarization insensitivity [1, 2, 3]. While the experimental research on plasmonic PCDs has made rather significant progress, rigorous numerical modeling of this type of devices is still challenging.
Challenges in the simulation of PCDs reside in the multiphysics nature of the involved physical process. PCDs are based on photoconductive semiconductor devices, which generate carriers under optical pumping and produce THz currents under the driving force of an applied bias voltage [4, 5]. Physically, a PCD has two operation stages: (1) a steady-state corresponding to the non-equilibrium state of the semiconductor device under a bias voltage; (2) a transient-stage describing the dynamic process after a pumping laser impinges upon the device. The steady-state involves the interaction between electric potential and carrier dynamics, which is usually described by a coupled system of Poisson equation and drift-diffusion (DD) equations. In the transient-stage, several processes happen simultaneously: EM wave scattering, carrier generation and recombination, carrier drift driven by electric fields, and carrier diffusion. This stage can be described by a coupled system of Maxwell equations and time-dependent DD equations. Due to the strong nonlinear coupling between carriers and EM waves, especially the exponential growth of carrier densities upon the absorption of light, accurate modeling of this process calls for a time domain multiphysics solver.
Moreover, due to the involvement of different physics, characteristic scales in this device also differ by several orders of magnitude. Length scales range from nm (Debye screening length) to m (the device geometry features). Additionally, the plasmonic nanostructures and the localized plasmonic modes make the length scale even more complicated. Time scales are signified by carrier dynamics (including the advection and diffusion motions), the EM wave propagation, and the response time between their coupling. Those scales typically differ by orders of magnitude. The multiscale nature requires highly flexible discretization schemes. Furthermore, the long physical time, regarding to the ratio between the THz wave period and the optical wave period, requires high order accuracy in time integration.
Several numerical approaches have been proposed for PCDs simulations. For conventional PCDs, the finite difference time domain (FDTD) method has been extensively used [6, 7, 8, 9, 10]. In this approach, an analytic carrier generation rate is used to model the optical pumping. The carrier dynamics is modeled with time-dependent DD equations and Maxwell equations are used to model the THz radiation. All equations are solved with finite difference methods. For conventional PCDs, this approach has been shown to agree with experiments very well [7, 8, 9]. However, for plasmonic PCDs, this approach turns out to be inadequate. In plasmonic PCDs, a crucial design is the nanostructure, which can produce enhanced local EM fields and hence enhance the carrier generation. The nanostructure leads strong scattering of the pumping EM wave, making the previously used analytic generation rate invalid. Meanwhile, the plasmonic structures and the plasmonic mode patterns are rather complicated, which require prohibitively dense meshes if using orthogonal FDTD grids. Several other methods are also reported, e.g., in [11, 12]13] an equivalent circuit model is used to estimate the optical-to-THz conversion efficiency, in  the semiconductor device is treated as a one dimensional port model and the EM radiation is solved with the discontinuous Galerkin finite element time domain method. Although those approaches are efficient for the conventional PCDs, the adopted approximations make them inadequate for plasmonic PCDs due to the complicated scattering of the pumping EM wave.
To model plasmonic PCDs, recent works utilize the finite element method (FEM) [15, 16]. Taking advantages of the unstructured meshes in FEM, this approach can easily solve the complex EM field in the plasmonic nanostructures. The carrier generation is calculated from the EM field and used in time domain carrier dynamics simulations. However, because of using frequency domain FEM, the space-distribution of the carrier generation is considered as time-invariant and the time-dependency is approximated by an analytic pulse shape. This means the nonlinear coupling between EM waves and carriers is not fully considered in this approach. In practice, the EM scattering, the carrier generation, and carrier movements happen simultaneously. Due to the scattering of the pumping EM wave, the carrier generation rate is space-time dependent. The generated nonuniform carriers in turn lead to scattering of the optical EM wave, which will again change the carrier generation. Since carrier densities increase exponentially as the pumping wave energy is absorbed, the whole process is strongly nonlinear. Furthermore, the space-charge screening effect, i.e., the separation of electrons and holes induces polarization that cancels the electric field, which depends on time-varying carrier distributions, is not captured in the FEM-based approach. This effect is important for understanding the high-power saturation effect in PCDs and has been extensively studied [17, 18, 4].
In this work, we propose a multiphysics framework for modeling plasmonic PCDs using the discontinuous Galerkin (DG) methods [19, 20, 21]. The device is modeled in two stages as the physical operation stages described above. The steady-state is solved by a nonlinear solver solving a coupled system of Poisson equation and bipolar DD equations. The transient-stage is modeled by solving time-dependent Maxwell equations and DD equations simultaneously, with the nonlinear couplings between carriers and EM waves being fully taken into account. The two solvers are connected by separating static and time-dependent variables. All equations are solved with DG methods. DG methods provide high flexibility in meshing (as FEM does) and allow explicit time integration (like FDTD). By using high-order basis functions and high-order time integration schemes, they can easily achieve high-order accuracy. In addition, DG methods provide advanced features such as non-conformal meshes, adaptive hp-refinement, and local time stepping. When using explicit time integration, they also show high parallel efficiencies, making them efficient on clusters. Those properties make DG methods promising for multiphysics and multiscale simulations [22, 23]. In the context of plasmonic PCDs, plasmonic structures and plasmonic modes can be modeled easily using unstructured meshes and the nonlinear couplings in the optical-to-THz conversion process can be captured naturally with explicit time integration. The flexible discretization also allows using different meshes and time step sizes for different equations, which can alleviate the computational effort a lot when dealing with the multiscale features.
This paper is organized as follows. Section II describes the physical processes involved in plasmonic PCDs. Two nonlinearly coupled systems are introduced to describe the steady-state and the transient-stage. Section III presents the numerical approach dealing with the nonlinear couplings and the DG discretization for each equation. Two examples are presented in Section IV, where a conventional PCD is modeled and compared with references and a plasmonic PCD simulation is presented to demonstrate the capability of the proposed framework. Finally, Section V concludes this paper with some remarks on applications of the multiphysics framework.
Ii-a Physical Model
PCDs are based on photoconductive semiconductors, of which the low temperature grown gallium arsenide (LT-GaAs) is the most popular one. LT-GaAs has electron trapping time less than 1 picosecond (ps) and a direct optical band gap of 1.42eV. They can absorb optical wave energy and generate pulsed photocarriers with pulse width ps, providing THz components. A typical PCD device has two electrodes deposited on the photoconductive semiconductor substrate. The electrodes are biased by an external voltage, which can drive the photocarriers into the electrodes. In the meantime, the electrodes act as parts of a radiating THz antenna. In conventional PCDs, the optical-to-THz conversion efficiency is limited by the amount of the absorbed optical energy, which is usually small due to the high refractive index of LT-GaAs.
Plasmonic PCDs utilize metallic nanostructures to increase the optical-to-THz conversion efficiency. The nanostructures are designed to introduce plasmonic resonances at the operating optical frequency, which can greatly enhance the local EM fields and increase the light absorption. Aside from introducing plasmonic resonances, the nanostructures, which are either placed in-between the gap of the electrodes or etched into the electrodes, also change the electric potential distribution and the distance between photocarriers and electrodes. The latter effects also significantly influence the photocurrent . Thus, the design of the nanostructures is not only to optimize the plasmonic EM modes, but to optimize the overall photocurrent response. Another factor of the device design is the impedance matching between the semiconductor layer and the THz antenna. However, this factor seems not so critical as the photocurrent generation, as is reported in several experimental reports [1, 15], where the output THz radiation power is observed to be proportional to the photocurrent. Thus, in this work, we focus on modeling the photocurrent generation and the THz radiation is not discussed.
Physical description of this device includes carrier dynamics and electrodynamics (including the bias electric field and EM waves corresponding to the pumping light and the THz radiation). While Maxwell equations are used to describe electrodynamics, various models exist for describing carrier transport in semiconductor devices, ranging from semi-classical to quantum approaches. For PCDs, the semi-classical DD model is widely used [25, 26, 8] since the device size (m) is far larger than the mean free path of the carriers (m). Regarding the difference of parameters for electrons and holes, and the space-charge screening effect, the bipolar DD equations are considered. Thus, the coupled system of electrodynamics and carrier dynamics is described by the following coupled equations
represents the location vector, the subscriptrepresents the carrier type and hereafter the upper and lower signs should be selected for electron () and hole (), respectively, and are the total electric and magnetic field intensities, and are the dielectric permittivity and permeability, and are the total electron and hole densities, and are the total electron and hole current densities, and are the recombination rate and the generation rate, respectively. We note that Gauss law is implicitly included in (1)-(3). In (4), and are the field-dependent electron and hole mobilities, and are the electron and hole diffusion coefficients, respectively.
In PCDs, a biased voltage is added on the electrodes throughout operation process. Before turning on the optical pump, the device is under a non-equilibrium steady-state . When the pump pulse enters the device, the EM wave energy is absorbed by the semiconductor and photocarriers are generated. Since the photocarrier density is several orders of magnitude smaller than the doping concentration, the biased steady-state is only weakly perturbed by the pumping. Thus, the total fields, the total carrier densities, and the total current densities can be separated into steady-state DC components and transient components as , , , , respectively. Here, the superscript “s” and “t” stand for steady-state and transient (time-dependent) variables, respectively.
Similarly, the total recombination rate in (3) is divided as , where describes the recombination rate of steady-state carriers and describes the recombination rate of photocarriers. In this work, we consider the two most common processes, namely the trap assisted recombination described by the Shockley-Read-Hall (SRH) model [26, 27]. The generation rate is defined as
where is the pumping wavelength, is the photon absorption coefficient , is the quantum yield, is the Planck constant, is the light speed, and is the energy flux density of the optical wave
which only depends on the time-dependent EM fields of the optical wave. Here, using the energy flux density allows us to model the carrier generation in the complicated space-time varying EM fields (resulted from the scattering from the plasmonic nanostructures and the nonuniformly distributed carriers). It is justified by the fact that the optical pump dominates the energy flux density and the contribution from THz frequency components is ignorable in the semiconductor substrate. If monochromatic wave is considered, the energy flux density returns to the optical intensity, which is widely used in analytical estimations [25, 12, 8, 15, 16].
In general, when using the DD model to simulate nonequilibrium state of semiconductor devices, one should consider field-dependency and velocity-saturation in the mobility model. Because mobility influences the carrier drift velocity, it is important to consider appropriate models for the photocurrent modeling. Here, we use the parallel-field dependent mobility model reported in . Since is two orders of magnitude stronger than , is ignored in the field-dependent mobility model.
Regarding to the fact that the steady-state components satisfy the coupled system of Poisson and DD equations
Here we should note that the same separation of variables has been done in .
Equations (11)-(II-A) is a strongly-nonlinear system. In (12), leads to exponential increasing of carrier densities. In (II-A), carriers are driven by both and . The carrier motion produces free current densities and , which results in EM wave scattering. This is described by the current term in (11). The EM wave scattering, in turn, changes . is also a nonlinear function, which becomes significant when carrier densities are high, balancing the carrier generation.
Ii-B The Multi-physics DG Solver
A complete simulator for PCD comprises a steady-state solver solving the Poisson-DD system (7)-(9) and a transient solver solving the Maxwell-DD system (10)-(II-A). The steady-state solutions are used as input parameters of the transient solver. The steady-state is solved by a nonlinear solver making use of the Gummel iteration method. For details, we refer the reader to . Here, we focus on the transient solver.
Ii-B1 Nonlinear Coupling
To solve the time-dependent Maxwell-DD system, we use an explicit time integration approach. Regarding the difference of time scales in EM wave propagation and carrier dynamics, Maxwell and the DD equations are updated separately as two subsystems. Such that we can use different time step sizes in the two subsystems. In PCDs, carriers response much slower than the variation of EM fields. Hence we can update the DD equations with a larger time step size without sacrificing accuracy. The nonlinear couplings between two subsystems are handled by alternately feeding updated solutions into each other.
The updating procedure is shown in Fig. 1, where the time step size for the DD equations is assumed to be twice the step size for Maxwell equations for illustration. Let us suppose the two systems are updated separately with certain time integration schemes (to be discussed in the next subsection) and denote the time steps of the Maxwell and DD systems as and , respectively.
Assuming that both subsystems are synchronized at step and (with the same physical time), first, and are used to calculate the generation rate using (5). and are used to update the DD equations from step to using (12). Then, the new carriers densities are used to calculate the current densities and (), which are used in (11) to update Maxwell equations at steps and , respectively. The physical time of two subsystems matches again at step and . One can again update DD equations to step with the EM fields and generation rate at time step .
Both the time-dependent DD equations (12)-(II-A) and Maxwell equations (10)-(11) are solved with DG methods in time-domain. We start with the time-dependent DD equations (12)-(II-A). Since the electron and hole DD equations only differ in the sign in front of the drift term and the values of physical parameters, we only discuss the electron DD equation and, to simplify the notation, we drop the subscript “e” (denoting electron) and the superscript “t” (denoting transient) when there is no confusion. Furthermore, as the three drift terms in (II-A) are treated in the same way, for brevity, in the following formulation we only discuss one drift term and denote it as . The electron DD model (12)-(II-A) can be written as the following initial-boundary value problem (IBVP)
where 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, denotes the outward normal vector of the surface, and and are the coefficients associated with these boundary conditions, respectively. For the problems considered in this work, is set on electrode/semiconductor interfaces and is set as (note that and are imposed on the steady-state components based on local charge neutrality [28, 27]). The homogeneous Robin boundary condition is used on semiconductor/insulator interfaces, indicating no carrier spills out those interfaces .
To facilitate the numerical solution of the IBVP described by (II-B2)-(17), is discretized into non-overlapping tetrahedrons. The volumetric support of each of these elements is represented by , . Let denote the element surface of and denote the outward unit vector normal to . Testing (II-B2) and (15) with Lagrange polynomials , , , on element and applying the divergence theorem yield the weak form
is the number of interpolating nodes,is the order of the Lagrange polynomials and
is used for identifying the components of the vectors in the Cartesian coordinate system. We note hereand are used to denote the local solutions on element and the global solutions on are the direct sum of the local solutions.
, , and are numerical fluxes “connecting” element to its neighboring elements. Here, the variables are defined on the interface between elements and the explicit dependencies on and are dropped for simplicity of notation. For the diffusion term, the LDG alternate flux  is used for the primary variable and the auxiliary variable . For the drift term, the local Lax-Friedrichs flux  is used for . In the interior of , they are defined as
where average and “jump” operators are defined as and , respectively, with could be 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 , while the precise direction of each variable is not important [29, 19, 20]. In this work, we choose the same as on each element surface. The local Lax-Friedrichs, with , mimics the information propagation path of the directional propagation. On , the numerical fluxes are chosen as , and . On , and . We note that and are not assigned independently on .
Expanding and with Lagrange polynomials 
where the unknown vectors are defined as and , with , .
In (22), and are mass matrices. is a block diagonal matrix with each block the same as . is defined as
is a diagonal matrix with diagonal entries , where , , , , and . We note that is assumed isotropic and constant in each element.
and are matrices corresponding to gradient operator and and are matrices corresponding to the divergence operator, respectively. Using the LDG flux, one finds and . is a matrix and it has contribution from the volume integral term and the surface integral term on the right hand side of (II-B2): , , , where
Here, denotes the interface connecting element and , is defined to select nodes on the interface
corresponds to the surface integral term [in (II-B2)] involving unknowns from neighboring elements of the current element , , where
Similarly, is contributed by the fourth term (the volume integral) and the fifth term (the surface integral) on the right hand side of (II-B2): , where
corresponds to the last surface integral term [in (II-B2)] involving unknowns from neighboring elements,
and are contributed from the force term and boundary conditions (where element does not exist) and are expressed as
For the time discretization of (22), an explicit third-order total-variation-diminishing (TVD) Runge-Kutta method  is used to incoparate with the high-order accuracy of the spatial discretization. With intial values being zero, is iterated step by step in time. By using the explicit scheme, we implement an MPI-parallelized solver which does not involve solving any linear system in the time iteration. Meanwhile, from (22), we can see is a local variable which only needs to be defined as a temporary array of size and is flushed repeatly when looping over elements. For this reason, this scheme is called local DG method.
and the flux vector is defined as
Discretizing the simulation domain into elements, testing (23) with Lagrange polynomials on each element, and applying the divergence theorem twice yield the strong form
where the subscript “” denotes that the corresponding variable is defined locally on element . To recover a meaningful global solution, the local solutions are required to be identical on the interface between two neighboring elements, which is guaranteed by the uniquely defined numerical flux on the interface. For Maxwell equations, which belongs to hyperbolic problems (with real and distinct characteristic paths), solving the Riemann problem using the Rankine-Hugoniot condition yields the upwind flux, which leads to the expression
where the average and jump operators are the same as those defined before, and are wave impedance and wave admittance, respectively. Here, the variables are defined on element surfaces and their explicit dependency on and is dropped for simplicity of notation. On boundaries where the neighboring unknowns are not defined, the numerical flux is assigned according to the boundary conditions. That is, on perfect electric conductor (PEC) boundaries, , , and, for absorbing boundary conditions (ABC), , . We note that, since the performance of ABC degrades rapidly with the angle of incidence, in addtion to ABC, we also use the perfectly matched layers (PML) to truncate the simulation domain. Details of the PML implementaion are reffered to .
We expand and with Lagrange polynomials
where denotes the location of interpolating nodes, denotes the components of vectors in the Cartesian coordinate system. and are the unknown coefficients to be solved for. Substituting (25) and (26) into (II-B2) yields
where and are elemental unknown vectors, , are the current density vectors, and are the components of and , respectively. The mass matrix , stiff matrices , and lift operator are defined as