IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 56, NO. 1, JANUARY 2008 113 Mixed Finite-Element Time-Domain Method for Transient Maxwell Equations in Doubly Dispersive Media Burkay Donderici and Fernando L. Teixeira, Senior Member, IEEE Abstract—We describe a mixed finite-element time-domain algorithm to solve transient Maxwell equations in inhomogeneous and doubly dispersive linear media where both the permittivity and permeability are functions of frequency. The mixed finite-ele- ment time-domain algorithm is based on the simultaneous use of both electric and magnetic field as state variables with a mix of edge (Whitney 1-form) and face (Whitney 2-form) elements for discretization of the coupled first-order Maxwell curl equations. The constitutive relations are decoupled from the curl equations and cast in terms of (auxiliary) ordinary differential equations in- volving time derivatives. Permittivity and permeability dispersion models considered here are quite general and recover Lorentz, Debye, and Drude models as special cases. The present finite-el- ement time-domain algorithm also incorporates the perfectly matched layer absorbing boundary conditions in a natural way. Index Terms—Dispersive media, finite element time domain. I. INTRODUCTION F INITE-ELEMENT time-domain methods have been suc- cessfully employed over the years for the simulation of transient Maxwell equations in complex geometries [1]–[4]. In contrast to the finite-difference time-domain method [5], the finite-element time-domain method requires the solution of a sparse linear system at every time step; however, when used in conjunction with simplicial grids, the latter is free from staircasing error. Traditionally, the finite-element time-domain method is based on the solution of the second-order vector wave equation for the electric (or magnetic) field after the elim- ination of the magnetic (or electric) field [3]. This facilitates the expansion of the unknown field by a single type of basis function. For the electric field, the basis functions of choice are typically edge elements. The choice of edge elements ensures conformity to a discrete version of the de Rham complex (exact sequence property) [6], [7]. This implies that discrete solutions that are not divergence free necessarily correspond to static fields (gradient-like eigenmodes at zero frequency), which do not pollute the frequency spectrum. Edge elements can be viewed as the natural interpolants of 1-forms and, hence, Manuscript received April 29, 2007; revised August 14, 2007. This paper was supported in part by the Air Force Office of Scientific Research under Grant FA 9550-04-1-0359, by the National Science Foundation under Grant ECCS- 0347502, and by the Ohio Supercomputer Center under Grant PAS-0061 and Grant PAS-0110. The authors are with the ElectroScience Laboratory, Department of Electrical and Computer Engineering The Ohio State University, Columbus, OH 43212 USA (e-mail: donderici.1@osu.edu; teixeira@ece.osu.edu). Digital Object Identifier 10.1109/TMTT.2007.912217 provide the correct discrete representation for the electric field intensity, which in its most general mathematical incarnation is represented as a 1-form in a differential manifold [6]–[18]. In the frequency-domain second-order vector wave equation, the divergence-free nature of the (source free) solutions is enforced in a frequency-dependent fashion (more precisely, quadratically with frequency). Hence, the divergence-free condition becomes progressively weak in the zero-frequency limit, leading to ill-conditioned matrices, which negatively impacts accuracy and convergence [19]. This problem also plagues adaptive mesh refinement because the field approaches the static limit locally for elements in highly refined portions of the grid. In the time domain, a different, but related, type of problem arises: the source-free second-order wave equation admits secular solutions of the form . These are spurious (nonphysical) modes in the null space of the curl operator and in the null space of the second-order time derivative. These spurious modes require elimination using, for example, a tree–cotree decomposition (gauging) [20], [21], grad–div regularization terms, or a posteriori filtering approaches [22]. Alternatively, the finite-element time-domain method can be based directly upon the first-order coupled Maxwell curl equations. In this case, both electric and magnetic fields are employed as unknowns, and mixed finite elements are used. A mixed finite-element time-domain method em- ploying the electric field intensity and magnetic field flux as simultaneous state variables have been considered, e.g., [23]–[27] where edge elements (Whitney 1-forms) are used for and face elements (Whitney 2-forms) are used for . This choice satisfies a discrete version of the de Rham diagram as well. Other desirable characteristics of this approach are: 1) it is free of secular solutions with linear growth; 2) it pro- duces energy-conserving (symplectic) algorithms [28] under an appropriate choice for the time integration scheme; 3) it provides a natural path for hybridization with the finite-differ- ence time-domain method since the latter can be formulated in terms of edge and face elements (now hexahedral) as well; and 4) it is more easily extended to complex media (i.e., with frequency dispersion and/or anisotropy). Although this mixed – finite-element time-domain method utilizes two fields as unknowns, its computer time and memory costs are comparable to those of the second-order vector wave equation formulation. This is because the former requires discretization of first-order time derivatives, while the latter requires second-order deriva- tives. As a result, only one past time-step electric/magnetic field 0018-9480/$25.00 © 2007 IEEE