Advances in Numerical Analysis

Volume 2010, Article ID 923832, 28 pages

http://dx.doi.org/10.1155/2010/923832

## Finite Element Approximation of Maxwell's Equations with Debye Memory

BICOM, Institute of Computational Mathematics, Brunel University, Uxbridge, Middlesex UB8 3PH, UK

Received 18 May 2010; Accepted 31 August 2010

Academic Editor: Dongwoo Sheen

Copyright © 2010 Simon Shaw. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Maxwell's equations in a bounded Debye medium are formulated in terms of the standard partial differential equations of electromagnetism with a Volterra-type history dependence of the polarization on the electric field intensity. This leads to Maxwell's equations with memory. We make a correspondence between this type of constitutive law and the hereditary integral constitutive laws from linear viscoelasticity, and we are then able to apply known results from viscoelasticity theory to this Maxwell system. In particular, we can show long-time stability by shunning Gronwall's lemma and estimating the history kernels more carefully by appeal to the underlying physical fading memory. We also give a fully discrete scheme for the electric field wave equation and derive stability bounds which are exactly analogous to those for the continuous problem, thus providing a foundation for long-time numerical integration. We finish by also providing error bounds for which the constant grows, at worst, linearly in time (excluding the time dependence in the norms of the exact solution). Although the first (mixed) finite element error analysis for the Debye problem was given by Li (2007), this seems to be the first time sharp constants have been given for this problem.

#### 1. Introduction

The potential for noninvasive electromagnetic detection of biological anomalies in the human body (and of defects in structural artefacts) has been recently noted by Banks et al. in [1]. Such a diagnostic technology has the promise of being quick, painless, cheap, and readily deployable, [2]. Biological tissue is dispersive, and electromagnetic constitutive relationships represent this frequency dependence through “complex moduli” (see, e.g., the data for grey and white brain tissue in [3, 4]). These are manifested though hysteretic fading memory Volterra operators or as an equivalent set of evolution equations, [1, 5]. The basic idea behind such a diagnostic technology is to compare the result of the transmission of electromagnetic waves through a patient's tissue against a datum outcome from a healthy subject. The presence of unwelcome tumours or lesions would be signalled by a change in the lossy response of the biological dielectric.

A possibly more sinister side of the interaction of electromagnetic waves with biological tissue is the concern that wireless communication devices (e.g., cell phones and mobile internet) may have a long-term detrimental effect on health. As far back as 1977 one can find (at http://oai.dtic.mil/oai/oai?verb=getRecord&metadataPrefix=html&identifier=ADA051218 (accessed 26 August 2009)) the following abstract of a workshop dedicated to the concerns that increased use of microwave radiation might raise [6].

The wide application of industrial, commercial and military devices and systems which radiate frequencies in the radiofrequency and microwave portion of the electromagnetic spectrum plus numerous only partially understood indications of microwave effects upon living organisms have raised important questions of the physical basis of the interactions of electromagnetic fields with biological systems. These questions must be answered if the development of regulatory standards and of methods and techniques for controlling radiofrequency and microwave exposure is to be achieved. The same questions must be answered in connection with present and proposed therapeutic applications of these waves. The rapid increase in the use of these frequencies makes these questions matters of imperative concern, particularly in view of the possibilities of cumulative or delayed effects of exposure. The purpose of the Workshop on the Physical Basis of Electromagnetic Interactions with Biological Systems was to bring together the leading investigators in the field to present the results of recent research, to determine the present status of the field and the priority of significant problem areas, and to critically evaluate conflicting theoretical interpretations and experimental techniques.

We can also find more recent contextualising discussions of the potential effects, both good and bad, of the dramatic increase in the use of nonionizing electromagnetic radiation in our world: [7, 8].

The purpose of this paper is to formulate Maxwell's equations with a Debye (i.e., “lossy”) polarization model and give some stability estimates for two continuous models as well as a discrete model, and also some error bounds for the discrete model. We rely on a strong analogy with models of viscoelasticity to provide some high quality estimates in terms of long-time behaviour. To our knowledge this is the first time such an analogy has been used in this context. To elaborate further it is required that we first introduce the mathematical model.

Let be a bounded polyhedral domain in and a finite time interval. We consider Maxwell's equations (e.g., Monk, [9]) in in the form,
where here and throughout the overdot denotes partial time differentiation. In these is the electric field intensity, is the magnetic field intensity, is the electric displacement, and is the magnetic induction. The first equation is *Faraday's law*, the second is *Ampere's modified circuital law*, the last two are *Gauss's law* and *Gauss's magnetic law*, is the total charge density, and is the current density. Conservation of charge implies that and, combined with the above, means that if the Gauss laws hold for one time, , then they hold for all . We therefore build these conditions into the intial data, and take the opportunity to specify the boundary data.
where is the unit outward normal on the boundary.

We assume that where , which is standard, but are concerned here with the Debye model of a dielectric wherein is given as a linear functional of . This will introduce “history integrals” and, therefore, nonlocality in time (but the constitutive law will remain local in space). Our contribution lies in drawing an analogy between the Debye polarization model and the fading memory constitutive laws governing creep and relaxation in linear viscoelastic media. This will allow us to apply known estimates from viscoelasticity theory and derive long-time stability estimates for the solutions to Maxwell's equations (Gronwall's inequality is not used). We also demonstrate a fully discrete time stepping scheme, based spatially on finite elements and temporally on central finite differences, for which these stability estimates carry over exactly. This provides a foundation for long-time numerical integration of the system, and in this spirit we also present an outline error estimate to show that again Gronwall's inequality may be avoided. It is also worth noting that sharp stability estimates for the dual problem are often needed when deriving residual-based *a posteriori* error bounds for time dependent problems; see, for example, [10]. In such cases it is essential to avoid the use of Gronwall's inequality whenever possible.

The original motivations for this work are the papers by Banks et al. [1] and Young and Adams [5], and we will return to these below. We note also that Fabrizio and Morro in [11] have a wealth of information on dielectrics and conductors with memory, for example, general observations on memory effects in Chapter 4, dielectrics and conductors with memory in Chapter 5.2, as well as more general formulations where all three constitutive laws (including Ohm's law) have fading memory in Chapter 9.3.3. They also study the temporally-asymptotic behaviour of dielectrics with memory in Chapter 7.10 using the Fourier transform and in this sense the general flavour of the results given below will come as no surprise. However, we use “energy methods” and explicit representations of the history kernels which are far more revealing in terms of the damping mechanisms at work. These “energy methods” are also well suited to the application to numerical schemes in which the constants can be estimated more carefully (see, e.g., [12–15]).

Besides the temporal decay properties with which we are concerned below, we also would like to briefly remark on the spatial decay exhibited by Debye media. For example, in [16] Roberts and Petropoulos examine the amplitude asymptotics of short and long pulses with application to muscle mimicking materials. They are able to derive tight upper and lower bounds on the exponential spatial decay rates of the temporal norm of the electric field, and they also discuss the nontrivial structure of the “wave hierarchy” induced in the dielectric. This structure was (seemingly) first uncovered in [17] where the Debye “modes” were combined to form a high-order PDE for the electric field.

To provide some more context for both the physical and mathematical aspects of this work, we note that Lakes in [18, Chapter 2] discusses constitutive equations for piezoelectric materials with memory. Here both strain and electric displacement are related to stress, electric field, and temperature through convolution integrals. On the other hand, in terms of numerical approximation schemes, it seems that the first time domain finite element methods (TDFEM, as in [19]) for Maxwell's equations in lossy media were proposed in [20]. They considered cold plasma, Debye and Lorentz media, with single or multiple poles, where the hysteretic effects were captured through a convolution, or “history”, integral. A discontinuous Galerkin TDFEM is given in [21] along with extensive numerics that illustrate also how the lossy medium can be coupled to a perfectly matched layer through the “auxiliary differential [relaxation] equations”, as opposed to a history integral. The first error estimates for Maxwell's equations with loss were developed by Li in [22] for a semidiscrete Galerkin approximation to the vector wave equation for the electric field in cold plasma. This work was extended to an implicit three-stage first-order time stepping scheme in [23], and then to the error analysis of a fully discrete mixed finite element method for cold plasma as well as one-pole Debye and two-pole Lorentz media in [24]. Furthermore, *a posteriori* error bounds for an interior penalty discontinuous Galerkin method are reported in [25].

Our presentation here is intended to complement Li's efforts in that we provide similar (Galerkin) formulations and error bounds but use the strong analogy to viscoelastic fading memory in order to avoid the use of Gronwall's inequality.

This paper is laid out as follows. In Section 2 we briefly review the aspects of linear viscoelasticity that are relevant and then in Section 3 discuss the Debye model and illustrate the similarities to viscoelasticity. This will allow us to make some general assumptions on the Debye kernels which will form the basis of our estimates. In Section 4 we formulate Maxwell's equations with memory as a problem for the electric displacement, , in Section 4.1 and also as a problem for the electric field, , in Section 4.2. In both cases, long-time stability estimates (Theorems 4.3 and 4.5) are derived without recourse to Gronwall's inequality. We give a fully discrete scheme for the second of these problems in Section 5 and show in Theorem 5.4 that the stability estimate continues to hold. A formal *a priori* error bound is stated in Theorem 5.6 which reveals that the constants are not time dependent. It follows from this that the only source of temporal error growth must come from the accumulation effects arising directly from the norms of the exact solution. By making standard assumptions on the approximation properties of the scheme this dependence is made explicit in Corollary 5.9, where we see that the convergence rate is optimal. Finally, Section 6 contains a few concluding remarks.

Everywhere below we will assume that *Ohm's law* holds so that , where is known, but will not exclude the case where the conductivity, , is nonzero since the medium may not be a perfect dielectric and so may exhibit some limited conductivity. At the mathematical level this term is included for the sake of generality and to demonstrate that it can easily be dealt with by using the estimates given below. Also, for the sake of simplicity we assume that all physical parameters are constant in both space and time.

Specific notation will be explained at the point where it is introduced but otherwise our notation is standard. In particular, , , and are the usual Lebesgue, Banach and Hilbert spaces and we set , , and so on. In general terms, if is a Banach space then its norm will be denoted by . The only exception to this will be that will always mean the norm as induced by , the inner product. We will not distinguish norms and inner products between the cases where they act on scalar- or vector-valued functions and time dependence is dealt with in the usual way. Indeed, given a map , then we use standard notation such as to signify the norm of .

The natural space for Maxwell's equations as we will work with them is where with graph norm, . This is a Hilbert space when equipped with the obvious norm-inducing inner product (see, e.g., [9]). We also define, (also a Hilbert space) where is a.e. the unit outward normal on . For Friedrich's-type properties of and associated spaces we refer to Monk in [9, Corollary 4.8].

Finally, in terms of preliminary notation, we recall Young's inequality for products, and Young's inequality for convolutions, and such that .

#### 2. Elements of Viscoelasticity Theory

Before we recall the Debye model, it will be instructive to outline some basic facts from the constitutive theories of linear viscoelasticity theory as described in, for example, [26–28]. For a *synchronous* viscoelastic substance, the stress tensor, , is given as a linear functional of the strain tensor, , by the constitutive law in the following *stress relaxation form*:
This can be inverted to the following *creep form*:
In these is a *stress relaxation function* and is a *creep function*. Also, the subscript denotes partial differentiation so that , for example, where the prime denotes partial differentiation with respect to the *displayed* time argument. Below we will use each of these notations with the aim of clarifying our reasoning in different places and contexts. For example, the -integral of an -derivative in the expressions above is useful later in Lemma 4.1, whereas the prime notation is useful in Assumptions 2.1 below where it neatly characterises the physical assumptions.

Causality and thermodynamics imply rather general conditions on and and, therefore, a good deal of freedom in ascribing specific forms. However, in many practical examples they typically have the Prony form,
subject to normalization, , and where and for each . When we have a *viscoelastic solid* and when we have a *viscoelastic fluid* (see, e.g., [14, 26] but note that the above form for will not apply in this case because as ). The tensor is the Hooke tensor from linear elasticity theory and is subject to some standard symmetry and positivity requirements, but these are not important here.

If we denote the Laplace transformation as and denote by then it is readily shown that the relaxation and creep functions are related by and it follows that, for a solid, . In particular, it is easy to derive the pair, These simple observations illustrate the following more general properties (again, see [26]) which will be our assumptions for all of what follows.

*Assumptions (creep-relaxation pairs). *A normalized creep-relaxation pair subject to satisfy: (i) with , , and ,(ii) with , , and , along with .

We are now in a position to review the Debye model and draw the analogy to viscoelasticity.

#### 3. The Debye Model

In the paper by Banks et al. [1], the electric displacement is written in terms of the electric field and polarization as where is the permittivity of free space and . Here is the instantaneous polarization and is the “relaxation polarization”. With , the constitutive model used in [1] is then
Note that here and throughout the dependence on the spatial variable is suppressed. Integrating this ODE gives
which, by henceforth assuming that (since the instantaneous effects are already taken care of in ), gives us
The connection to viscoelasticity can now be made evident by defining the *creep function *
where , , , , and , and then for , we have,
Immediately we can observe the analogy, or correspondence, between the electric displacement and field, , , and the viscoelastic strain and stress, , . We can also conclude that this relationship can be inverted to give,
for the *relaxation function *. It is clear that and we also have . On physical grounds, we expect that and . We then deduce that is monotone decreasing but bounded above zero by and that is monotone increasing but bounded above by . It therefore follows that this pair satisfy Assumptions 2.1.

Now, under the static conditions where is constant and we have and so, with the Euclidean norm on , we have for all . We can see that this constitutive law represents a monotone creeping process for the electric displacement which is bounded by and .

We note that the Debye model presented in [11, Chapter 4, Section 2] agrees with that just presented if, in [11], we set and .

Just as in viscoelasticity theory the relaxation and creep functions can be represented as a Prony series of decaying exponentials rather than just the “single mode” model just given. This is the model presented by Young and Adams in [5]. They write The connection between this model and the Banks et al. model is clear and we have, Assuming, as before, that each we can therefore write where and defined so that . Drawing again on the viscoelasticity analogy we conclude that again this constitutive law can be inverted with the aid of a relaxation function, , to give Typically, for the creep function in use here we would expect that but this, of course, needs proof. We assume only that the pair satisfy Assumptions 2.1 and give specific examples later in Section 6.

#### 4. Maxwell Systems with Memory

In this section, we introduce the Debye model in each of the forms (3.10) and (3.12) into the field equations, (1.1), the intial and boundary data, (1.2) and (1.3), and to derive partial differential Volterra equations. Firstly, in Section 4.1 we use (3.12) to obtain a PDE for the electric displacement, , and then in Section 4.2 we use (3.10) to derive a PDE for the electric field, . Stability estimates that exploit the “viscoelastic fading memory” in place of Gronwall's lemma are given in each case. We also recall that, for simplicity, all physical parameters are assumed to be constant in space-time.

##### 4.1. A Formulation for the Electric Displacement

From (1.1) and Ohm's law, we obtain
which, recalling (1.5), can be written in weak form as
However, from (3.12) with we obtain after a partial integration
Note that in addition to using an overdot for partial time differentiation, we are using a prime to denote differentiation with respect *to the displayed variable* and subscripts to denote partial differentiation in the usual way.

By introducing these terms into the weak form, we can pose the problem as: find a twice differentiable map, such that, with known and given intial data.

Before we give the stability estimate, we need two lemmas which help us deal with the fading memory in a sharper way than Gronwall's inequality would admit. The first is a coercivity result as used in [29], and the second is a remarkable lemma given originally in [12] by Rivera and Menzala. (The author is very grateful to Professor Igor Bock (Slovak University of Technology, Bratislava) for pointing out Rivera and Menzala's result.) Outline proofs are given for the sake of completeness.

Lemma 4.1 (coercivity of relaxation). *If Assumptions 2.1, hold then
**
for all .*

*Proof. *Since in general terms we need only to observe that from Young's inequality, (1.7). Assumptions 2.1 yield and this completes the proof.

Lemma 4.2 (Rivera and Menzala, [12]). *Let , then
**
for all and for all .*

*Proof. *First we note that,
and substitute for the last term on the right by using
After a rearrangement of terms, this gives
The proof is then completed by integrating over , observing that
and recalling that .

Theorem 4.3 (stability for the equation). * If Assumptions 2.1 hold and if , and are real numbers, then for every we have the following. Firstly, if ,
**
Secondly, if and ,
**
Thirdly, if ,
*

*Proof. *We choose in (4.4) and integrate over to get
We now use Lemmas 4.1 and 4.2 and note that the signs of the and integrands mean that only nonnegative terms contribute to the left-hand side. This completes the proof for . When and , Young's inequality, (1.6), gives,
and we can complete the proof of the second claim. If and **,** we deduce that
and the third claim now follows.

##### 4.2. A Formulation for the Electric Field

We now derive an alternative formulation in which the electric field becomes the “unknown”. If the displacement formulation given above could be called the “relaxation problem”, then this could be called the “creep problem”. We again derive but now use (3.10) to get from which, Introducing this results in the partial differential Volterra equation, or, more usefully, the weak problem: find a twice differentiable map such that, with known and given intial data.

We require a slight adjustment to Lemma 4.1 before giving the stability estimate.

Lemma 4.4 (coercivity of creep). *If Assumptions 2.1 hold, then
**
for all .*

*Proof. *Since and this is similar to Lemma 4.1.

Theorem 4.5 (stability for the equation). * If Assumptions 2.1 hold and if , and are real numbers then for every we have first,
**
and second,
*

*Proof. *We choose in (4.21), integrate over and use Lemma 4.4 to get,
Now by applying the Cauchy-Schwarz and Young inequalities in the following form:
we complete the proof of the first claim. For the second, we easily see that
and with the Hölder and Young inequalities,
we arrive at the second claim.

*Remark 4.6 (growth for ). *If then the fact that as renders the first stability estimate in the above of limited interest.

We now move on to give a numerical scheme for the electric field equation and see that these stability estimates have an exact analogue.

#### 5. A Numerical Scheme for the Electric Field

In this section, we suggest a fully discrete numerical scheme for the electric field problem, (4.21), provide a stability estimate that is essentially identical to that for the continuous problem, and also give an error bound showing that we still do not require Gronwall-type estimates. To prepare we first assume that the domain is “triangulated” with a quasiuniform mesh of hexahedra or tetrahedra, with meshwidth , and that a suitable finite element space, , is specified with respect to this partition. (We would probably use Nédélec elements, [30, 31], but we do not yet need to be specific.) The time interval is then discretized into intervals of equal width so that , where and . We continue to assume that , , and are real numbers but we recognise that some of the working and results given below could be made simpler by fixing on either of the more definite cases or .

Our numerical scheme employs a Galerkin finite element method for spatial discretization and a standard (for the wave equation) two-field centered finite difference approximation in time. For this we define the “velocity” field . Before we begin we just mention that numerical schemes for Maxwell's equations (see, e.g., [32, 33]) are often based on the first two of (1.1) and assume that and are given and (since and *μ* are assumed constant here) satisfy the divergence conditions in (1.2). We will need in the scheme below and also and we note that the latter quantity is available if is because, from (4.3), and from (3.10) and (1.1) we see that and . Based on this, from now on we just assume that we have physically admissible (in terms of (1.2)) intial data such that
Our notation is with
Notice that and represents the trapezium rule.

The numerical scheme is introduced next in Section 5.1 where we establish stability and well-posedness. The error bound is then given in Section 5.2.

##### 5.1. The Numerical Scheme

Our numerical approximation to the electric field wave equation, (4.21), is the problem: for each in turn, find a pair such that, These are supplemented with discrete intial data, and , about which we will say more later (see later below (5.28)), and the discrete convolution, denoted by “”, is defined by, with empty sums set to zero and for any subset . Of course this is merely , the average of at and , for a temporally piecewise constant . Also, for (recall Assumptions 2.1), simple calculations give the alternative form which will be useful below. We next establish a discrete version of (1.7).

Lemma 5.1 (Young's inequality for discrete convolutions). * For sets and of nonnegative real numbers we have for every ,
**
for and the obvious interpretation if .*

*Proof. *Obviously, and so by two Hölder inequalities,
We now interchange the order of summation with
apply Hölder's inequality again and take th roots.

Our first aim is to derive a stability estimate for the discrete problem and thereby establish its well-posedness. To this end, we need a discrete version of Lemma 4.4 but, before this, we need a bound on the discrete convolution.

Lemma 5.2. *If Assumptions 2.1 hold, then
**
for every and every .*

*Proof. *We have from the definition in (5.5) and (5.6) that
and so applying Cauchy-Schwarz and then Lemma 5.1 gives
Now observe that .

We can now state the discrete version of Lemma 4.4.

Lemma 5.3 (coercivity of discrete creep). * If Assumptions 2.1 hold, then
**
for every and every .*

*Proof. *Noting that
the result follows from Lemma 5.2.

Defining the norm we can now state the main result of this subsection.

Theorem 5.4 (discrete well-posedness). * If Assumptions 2.1 hold and we are given approximations to then there exists a unique solution to the discrete problem (5.3) and (5.4). Moreover,
**
for each .*

*Proof. *Since uniqueness implies existence for this finite dimensional problem, we need only prove the estimate and then take . We select in (5.3) and in (5.4), merge them using the common term , and then sum the result over . Noting the usual telescoping sum that arises from products like , we then use Lemma 5.3 to obtain
It is clear from (5.4) that and so if we easily conclude that . Doubling the right-hand side to account for the norm a Young's inequality,
this completes the proof.

The stability estimate given above is an exact analogue of the second estimate given for the continuous problem in Theorem 4.5. We now examine the possibility of deriving an error estimate with a similar control over all of the “constants”.

##### 5.2. An Error Bound

To arrive at an error bound for the discrete approximation, we introduce the projection as the linear continuous map, such that, for every , is the unique solution to It is clear that the best approximation property, is satisfied. Introducing also the projection, , through for all , we then (as usual) split the errors according to so that and along with The first result is a representation formula for the spatially discrete components of the error and for this we define,

Lemma 5.5. *For each , we have
*

*Proof. *We average the continuous problem, (4.21), at the times and and subtract the result from the discrete problem (5.3). Simultaneously adding and subtracting like terms in and introduces the error components, , , , and , and we collect the lower case quantities to the left of the equals sign and all of the other terms to the right. We also notice that the term is zero while can be replaced by . A simple manipulation provides
and so choosing generates the terms on the right and with we can use (5.4), with , to get . Lastly we merge these and sum over .

The first main effort in deriving the error bound is contained in the next lemma. For this we note first the “summation by parts” formula, and, secondly, that on taking we can use (5.4) to get, Thirdly, if the discrete intial data are defined by and , then it follows that .

We can now give a *formal a priori* error bound—formal because it deals only with the time-stepping dynamics of the error rather than with convergence orders. Of course this bound is meaningless until we know that the unresolved terms can be controlled but, in this way, by comparing with the stability estimates in Theorems 4.5 and 5.4 we can see that the error is governed (bounded) by essentially the same temporal dynamics as both the exact and discrete solutions. The “concrete form” of the error bound will follow later in Corollary 5.9.

Theorem 5.6 (“formal” a priori error bound). * If Assumptions 2.1 hold and the discrete intial data are determined by and , we have
**
for every , where the and are given in (5.23) and (5.24).*

*Proof. *Lemmas 5.5 and 5.3 give
and using (5.27) and Young's inequality then gives,
for all . From (5.28), we find
for all . In fact, choosing in the last two estimates gives
and so by applying the triangle inequality in the form , we complete the proof.

To turn this estimate into something more concrete, and to see what the convergence rate is, as well as what the potential growth in time of the error is, we need to estimate all of the terms on the right-hand side and for this we need some preparation. The following arises in a standard way and will be used without proof. We will also need the following estimates for the continuous and discrete convolutions.

Lemma 5.7 (error in convolutions). * If Assumptions 2.1 hold then we have
**
for every and every .*

*Proof. *First, from the definition in (5.6), we get after using (5.36) that
and then summing over , we use Lemma 5.1 and arrive at
which is simplified to give the first result. The proof of the second is similar and begins with,
From (5.34) we define for , and then after a little work, we get
Again we sum over and apply Lemma 5.1 to arrive at,
The same simplification completes the proof.

Corollary 5.9 contains the “concrete” form of the error bound but before we get to it we need to estimate the terms on the right of Theorem 5.6. From now on, we are not so concerned with the precise form of constants that will appear and so we adopt the common practice of writing “” to mean that “” where is a generic constant. In our working will depend only on and (but not on or ).

Lemma 5.8 (auxiliary bounds for and ). * If Assumptions 2.1 hold then for every ,
*

*Proof. *First, from (5.24), (5.35), and (5.36)
and with a Young's inequality, it follows from this that
where the hidden constant is merely a positive real number. Now, with the observation