Efficient discontinuous Galerkin scheme for analyzing nanostructured photoconductive devices

: 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.


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][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 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][10][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][13][14][15][16][17][18].To this end, the finite element method (FEM) has been extensively used in recent years [19][20][21][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 ∼ 10 nm, the optical wavelength ∼ 100 nm, and the device size ∼ 10 µ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.

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 φ(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 φ 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 φ drop (y, z) in ( 1) is selected as described next.φ 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][2][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, φ drop (y, z) can be simplified to a single value φ drop .Additionally, as discussed before, on average the potential drops linearly between the two electrodes [31,32].Consequently, φ drop (y, z) can be estimated as: φ drop (y, z) = w x (V bias /w sd ), where the V bias is the bias voltage applied to the electrodes and w 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 φ 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][2][3][4][5]33].

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) where φ(r) and E(r) are the unknowns to be solved for, g(r) and f (r) are known coefficients, ε(r) is the permittivity, Ω is the solution domain, ∂Ω ν represents the surfaces perpendicular to the ν-direction, ν ∈ {x, y, z}, n(r) is the outward pointing normal vector on ∂Ω ν .The homogeneous Neumann boundary condition f N (r) = 0 is used in (9) [23].
Discretizing Ω 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][35][36].This process yields a matrix system as Here, Φ and Ē are unknown vectors storing the basis expansion coefficients, Mg and ME are block diagonal mass matrices, Ḡ and D are block sparse matrices representing the gradient and divergence operators, respectively, and Bϕ and BE 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] is used on each element surface in the interior of Ω, where "average" and "jump" operators are defined as {⊙} = 0.5(⊙ − + ⊙ + ) 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.β = n determines the upwind direction of φ and (εE) [34][35][36].On ∂Ω N , the numerical fluxes are chosen as φ * = φ − and (εE) * = f N [35].Here, the variables are defined on the surfaces and the explicit dependency on r 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 where U ∈ {φ, (ε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 φ * = φ(w x /2, y, z) + 0.5φ drop (y, z) = φ(−w x /2, y, z) − 0.5φ 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] ∇ where n(r) and q(r) are the unknowns to be solved for, and d(r) v(r), and R(r) are known coefficients within the Gummel method [24].The homogeneous Robin boundary condition f R (r) = 0 is used in ( 23) [23].The BVP for holes only differs by the sign in front of the drift term.DG discretization of ( 19)-( 23) yields a matrix system as where N and Q are unknown vectors storing the basis expansion coefficients, Mq is same as ME , Ḡ and D are same as before, d is a diagonal matrix, and C corresponds to the drift term [24].The local Lax-Friedrichs flux [24] (vn is used for the drift term and the alternate flux ( 11)-( 12) is used for the diffusion term (φ and (εE) are replaced with n and q, respectively).The PBC (22) is enforced just like it is done in ( 13)-( 14), with U ∈ {n, (dq)}, 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].

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][35][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] n(x, −w y /2, z, t) = n(x, w y /2, z, t), r ∈ ∂Ω y , ( 30) where n(r, t) and q(r, t) are the unknowns to be solved for, d(r) and v(r, t) are known coefficients, and G(E, H) is the generation rate [18,23].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 (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] Mkk where Nk (t) and Qk (t) unknown vectors storing time-dependent basis expansion coefficients, Bn k (t) and Bq k are vectors constructed from the net-recombination and boundary conditions (other than PBC), Mkk , Ckk / Ckk ′ , Ḡkk / Ḡkk ′ , and Dkk / Dkk ′ are the elemental mass, convection, gradient, and divergence matrices, respectively [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].

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 µm.
The nanograting has a periodicity of 0.18 µm in x and y directions and its height is 0.12 µm.The metal block is a truncated square pyramid with dimensions of 0.06 µm and 0.1 µ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 with ε ∞ = 5.785, ω o = 4.783 × 10 15 rad/s, γ = 4.557 × 10 14 rad/s, ω p = 1.061 × 10 16 rad/s.SI-GaAs is modeled as a dielectric material with ε r = 13.26.The nanogratings and contacts are modeled as gold using the Drude model with parameters ε ∞ = 1.0, ω p = 1.372 × 10 16 rad/s, γ = 8.052 × 10 13 rad/s [23].All materials are considered nonmagnetic µ r = 1.0.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).The bias voltage is fixed at V bias = 10 V.The PCD is operated in the continuous-wave mode [1,4] and excited from top by two continuous-wave lasers operating at 374.5 THz and 375.5 THz, with a frequency difference of 1 THz.The laser beam width is 1.8 µm in the full-device simulation.
Table 1.Semiconductor material parameters used in the PCD example Recombination τ e = 0.3ps, τ h = 0.4ps 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 µm and 0.54 µm.The nanograting is made of a 15 × 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 is used on electrode/semiconductor interfaces and the homogeneous Robin boundary condition n•J e(h) = 0 is used on semiconductor/insulator interfaces [27].For Poisson equation, the Dirichlet boundary condition φ = V el + V T ln (n s e /n i ) is enforced on electrode surfaces and the homogeneous Neumann boundary condition n • ∇φ = 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 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 nm and 200 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 ps and the time-step sizes for the Maxwell and DD equations are 4 × 10 −7 ps and 2 × 10 −6 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 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 ∼ 1 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 tradeoff 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 it N 2 ) and O(N it N), where N it is the number of GMRES iterations and N is the number of unknowns.Note that N 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.

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.

Fig. 1 .
Fig. 1.(a) Geometry of a nanostructured PCD and the corresponding mesh used in the full-device DG scheme.(b) The unit cell used by the proposed scheme for the PCD geometry in (a).

Fig. 3 .
Fig. 3. (a) Electric potential and (c) electron density computed using the unit-cell scheme.(b) Electric potential and (d) electron density computed using the full-device simulation.The displayed solutions are extracted from the interface at z = 0.5 µm.

Fig. 4 .
Fig. 4. The x component of the transient current density obtained from the unit-cell and the full-device simulations.
Ēν k (t) and Hν k (t) are unknown vectors storing time-dependent basis expansion coefficients, Jν k (t) is the current density vector, ε k and µ k are the constant permittivity and permeability in element k, Dν Lk , where Mk , Sν k , and Lk are the mass, stiffness, and surface mass matrices, respectively k = M−1 k Sν k , and Fk = M−1 k