Recent progress on mathematical analysis and numerical simulations for Maxwell's equations in perfectly matched layers and complex media: a review

1.   Introduction

The well-known Maxwell's equations are a system of coupled partial differential equations (PDEs), which form the foundation of classical electromagnetism, classical optics, and electric circuits. The equations provide a mathematical model for countless applications involving electric, optical, and radio technologies, such as lenses, radar, power generation, electric motors, wireless communication, etc. Due to the importance of Maxwell's equations, many robust and efficient numerical methods have been developed and implemented in solving Maxwell's equations over the past several decades (cf., monographs ^{[1,2,3,4,5,6,7,8]} and references therein).

In recent years, many new electromagnetic materials (e.g., metamaterials ^[9,10,11], graphene ^[12,13]) have been successfully constructed. To efficiently simulate wave propagation in these new media, many new numerical methods have been developed in recent years. In this paper, we give a review on some recent progress on the mathematical analysis and numerical simulations for Maxwell's equations in these media.

The rest of the paper is organized as follows. In Section 2, we present a review on metamaterials. In Section 3, we give the review on graphene. In Section 4, our review goes to the perfectly matched layer (PML), which can be treated as a special medium. Finally, we conclude our review in Section 5.

2.   Metamaterials

Electromagnetic metamaterials are artificially structured materials with some exotic properties that are not observed in natural materials. Here we are interested in a specific type of metamaterials with both permittivity and permeability being negative. Study of such a negative index metamaterial (NIM) was initiated by Veselago ^[14] back in 1968. In this seminar paper, Veselago found that simultaneously negative permittivity $\epsilon$ and permeability $\mu$ of a material lead to a negative refractive index in the material, and waves in such metamaterials travel backward toward the source wave location. Since NIMs do not exist in nature, Veselago's work didn't get much attention until the successful construction of such a material in 2000 ^[15] and successful demonstration of the negative refraction index in 2001 ^[16]. Another catalyst for great interest in studying metamaterials is caused by Pendry's perfect lens paper ^[17]. According to [11, p.317], these four seminar papers together made the birth of the subject of metamaterials.

Since 2000, there has been a tremendous growing interest in studying the design of metamaterials and the investigation of their potential applications ^[10,18] ranging from electronics and telecommunications to sensing, radar and defense, nanolithography with light, subwavelength imaging, data storage, and electromagnetic cloaking device to achieve invisibility. Detailed overviews of the state of the art in the progress on metamaterials can be consulted from some recent books published in this area (e.g., ^{[10,11,19,20]}).

The appearance of metamaterials provides many challenging problems that require novel incorporations of complex material models (e.g., ^[21,22,23]) into the Maxwell's equations, and it offers many new problems in terms of mathematical and numerical issues (cf. ^{[24,25,26,27,28,29,30,31,32,33,34,35]}). ^[36,37] provided some reviews on the mathematical models for metamaterials and their numerical analysis. Due to the special interest in invisibility cloaking, below we will focus on this subject.

The idea of designing invisibility cloaks with metamaterials was originated in 2006 by Leonhardt ^[38] and Pendry et al. ^[39] independently. In addition to a huge amount of publications in physics and engineering (cf. ^[20,40]), there are many excellent works published in the mathematics community, too. In ^[41], Ammari et al. used the scattering coefficients vanishing approach to consider near-cloaking for the full Maxwell equations in the frequency domain. Near cloaking for Maxwell's equations is also considered in ^[42]. When considering near-cloaking for the Helmholtz equation ^[43], Ammari et al. proved that cloaking is increasingly difficult as the cloaked object becomes bigger or the operating frequency becomes higher. Cloaking by anomalous localized resonance was proposed and studied intensively, too (e.g., ^[44,45]). However, cloaking by anomalous localized resonance has some serious limitations, and it can take place if and only if the dipole type source lies inside critical radii determined by the radii of the core and the shell. Abstract mathematical analysis of cloaking phenomena has been investigated in many papers (e.g., ^{[41,46,47,48,49]}), and numerical analysis and computer simulations of cloaking phenomena have been carried out by the Finite-Difference Time Domain (FDTD) method (e.g., ^[19]), finite element method (FEM) (e.g., ^[50,51,52]), and the spectral element method (e.g., ^[53,54,55]).

Broadband cloaking ^[48,56,57] encourages us to pursue the development and analysis of the finite element time-domain (FETD) method for simulating invisibility cloaks, which involves solving time-dependent Maxwell's equations supplemented with some auxillary equations resulting from the constitutive relations for the complex media ^{[21,25,26,27,58]}. In 2013, Wu and Li ^[59] showed that total reflection and total transmission (cloaking) can be achieved through a zero index metamaterial (ZIM) waveguide embedded with rectangular dielectric defects. In 2014, Li ^[60] analyzed a time-domain spherical cloaking model. Later, Li and his collaborators studied some time-domain cylindrical and elliptical cloaks ^[61] and several arbitrary star-shaped 2D electromagnetic cloaking models ^[62]. Time-domain finite element schemes were developed and cloaking phenomena were simulated by edge elements.

Here we just illustrate an interesting "carpet cloak" originally proposed by Li and Pendry in 2008 ^[63] by using the quasi-conformal mapping technique. The idea is to transform a bulging reflecting surface into a flat one, rendering anything under the bulging surface invisible from outside observers. Experimental realizations of carpet cloaking were successfully demonstrated from microwave regime to terahertz and optical frequencies. The governing equations for modeling the wave propagation in the carpet cloak were derived in ^[64] and given as follows (cf. [65, (2.1)–(2.3)]):

where we denote $H$ for the magnetic field, and ${{\boldsymbol{D}}}$ and ${{\boldsymbol{E}}}$ for the 2D electric displacement and 2D electric field, respectively. Moreover, $\partial_{t^k}u$ denotes for the $k$-th derivative $\partial^k u/\partial t^k$ of a function $u$, and we adopt the following 2D curl operators:

where $E_x$ and $E_y$ are the $x$ and $y$ components of ${{\boldsymbol{E}}}$, respectively. Here, ${\epsilon}_0$ and $\mu_0$ are, respectively, the permittivity and permeability in vacuum, $\mu$ is the relative permeability, $\lambda_2$ and $\omega_{p}$ are some positive constants, $M_A^{-1}$ is the inverse of a two-by-two positive definite matrix $M_A$, and $M_C$ is a two-by-two positive semi-definite matrix. Detailed expressions can be seen in ^[65].

Theorem 1. The following energy identity holds true for the solution $({{\boldsymbol{D}}}, H, {{\boldsymbol{E}}})$ of (2.1)–(2.3):

where we denote

Moreover, we have the following stability:

where the positive constant $C_*$ depends on the model physical parameters.

In ^[65], we developed two finite element methods to solve the carpet cloak model (2.1)–(2.3) with edge elements. Stability and convergence estimates were established also in ^[65]. Figure 1 shows some snapshots of the numerical magnetic fields $H$ obtained with an incident Gaussian wave

imposed on a 45-degree angle line segment. This figure shows that the incident wave is totally reflected from the flat ground. This produces the invisibility cloaking phenomenon, and our result is similar to the simulation obtained by software COMSOL for the acoustic carpet cloak [66, Figure 7.7].

3.   Graphene

Graphene is an atom-thick planar sheet of $sp^2$-bonded carbon atoms packed in a honeycomb lattice. Graphene was first isolated and characterized by Novoselov, Geim and co-workers ^[67] at the University of Manchester in 2004. The Nobel Prize in Physics 2010 was awarded jointly to Andre Geim and Konstantin Novoselov "for groundbreaking experiments regarding the two-dimensional material graphene." Since then, the study of graphene has attracted great interest of researchers (e.g., ^{[12,13,68,69]}) due to its many interesting properties. For example, graphene exhibits high electronic mobility as a result of its unique gapless conical band structure, which promotes the quantum Hall effect at low temperatures and generates high magnetic fields for both electrons and holes. Graphene has low absorption of light over a wide wavelength range, an extremely high thermal conductivity, and saturable absorption (an interesting phenomenon where the light absorption decreases with an increasing light intensity). Graphene also exhibits some unique mechanical properties, such as low mass density, high surface area-to-volume ratio, and Young's modulus. The excellent properties of graphene bring a broad range of promising applications in various fields. Researchers have been successfully utilizing graphene to generate picosecond laser pulses due to its wide absorption range, fast decay, and high stability properties. Graphene can be used in sensors to simultaneously sense gas, mass, charge, tension, diseases, and explosives; graphene can be used in many high-performance products, including low-cost display screens of mobile devices, lithium-ion batteries with fast recharge capacity, ultracapacitors with improved performance relative to that of batteries, hydrogen storage for fuel cell-powered cars, and low-cost fuel cells and water desalination, etc. Graphene and its derivatives have potential usages in bio-imaging, ultraviolet lens, frequency multiplier, Hall effect sensors, conductive ink, optoelectronics, spintronics, charge conductor, sound transducers, radio wave absorption, catalyst, waterproof coating, condenser coating, and coolant additive ^[70].

Though physical research on graphene and graphene-based devices has progressed a lot, the study of graphene in the mathematics community still lags behind. In recent years, some papers have published on the mathematical analysis of graphene and modeling (e.g., ^{[71,72,73,74,75,76,77]} and references therein). However, except our recent papers ^[58,78,79], to our best knowledge, we are unaware of other papers on time-domain finite element analysis and applications for graphene models.

In 2020, the author and his collaborators initiated a study on graphene based devices. Starting from the Kubo formula for the graphene conductivity, we derived the following graphene model ^[78]:

where ${{\boldsymbol{E}}} = (E_x, E_y)'$ and $H = H_z$ are the electric field and magnetic field, repsectively, ${\epsilon}_0$ and $\mu_0$ are the permittivity and permeability in vacuum, $\omega_{pe}, \gamma, b_2, b_1, a_2, a_1, a_0$ are physical constants, and ${{\boldsymbol{J}}}_d$ and ${{\boldsymbol{J}}}_p$ are the induced current densities through the intraband and interband conductivities, respectively. Note that this model is formed as a system of 7 coupled differential equations. Developing an accurate and efficient numerical method for solving it is quite challenging. In ^[78], we further proposed and analyzed an FETD method for solving the coupled system (3.1)–(3.4). The surface plasmon polaritons (SPP) propagating along the graphene sheet are achieved by our proposed FETD scheme, which treats the graphene as a thin plate with finite thickness.

The technique by treating graphene with some thickness can become impracticable due to the high requirements in terms of memory storage and computer time, since a very fine mesh is usually needed to adequately discretize the thin graphene layer in order to achieve a good numerical accuracy. In 2022, we proposed a new finite element scheme by treating the graphene as a zero-thickness sheet. In this case, the govening equations are Maxwell's equations coupled with an ordinary differential equation on the graphene interface ^[79]:

where we denote $K_s$ for an imposed magnetic source function, ${{\boldsymbol{J}}}: = {{\boldsymbol{J}}}_d$ (as denoted in ^[78]) for the induced intraband surface current in graphene, the positive constant $\tau_0$for the relaxation time, and the positive constant $\sigma_0$ for the graphene surface conductivity. Moreover, we denote $\Gamma$ for the graphene sheet buried in the physical domain ${\Omega}$. In a 2D domain, the graphene sheet appears as a line (cf. Figure 2 shown below). Comparing to the modeling equations with a finite thickness (3.1)–(3.3), the new model (3.5)–(3.7) is much simpler. We now need to impose the following graphene interface boundary conditions:

i.e., across the graphene interface, the tangential electric field is continuous, and the jump of the magnetic field is equal to the tangential surface current. Here we denote $H_1$ and $H_2$ for the magnetic field above and below the interface, respectively, $\hat{n}: = (n_x, n_y)'$ for the unit normal vector pointing upward, and $\hat{n}_1$ and $\hat{n}_2$ are the unit outward normal vectors from top and bottom subdomains of the interface, respectively.

Compared to treating graphene with some thickness, this way uses less meshes and avoids the stringent time step size caused by the fine mesh needed in the graphene region. Similar SPP phenomenon is achieved by our new method ^[79]. In Figure 2, we present the SPP simulation along a bifurcated graphene sheet to demonstrate the flexibility of our FETD scheme in handling the complex geometry. Some snapshots of the obtained numerical magnetic fields $H$ are presented in Figure 2, which shows that our new method captures the SPPs very well.

For graphene simulation, how to choose the right interface conditions is challenging. Engineers have proposed many different interface conditions ^[80] such as impedance transmission boundary condition, impedance matrix boundary condition, surface current boundary condition, and surface impedance boundary condition. How to incorporate these interface conditions to time-domain grahene simulations and theoretically justify them from mathematical point of view would be very interesting. Graphene is known to be spatially dispersive (also known as nonlocal), i.e., the graphene conductivity tensor varies with wavevectors. For example, when considering the graphene ribbon in the nanometer scale, the graphene's surface conductivity becomes spatially dispersive ^[81]. Hence, accurate modeling of the spatial dispersion effects is vital for the simulation of both passive and active 2D nano-device. However, the nonlocal conductivity makes the modeling equations more complicated and challenging. In 2018, A Discontinuous Galerkin time domain (DGTD) framework was developed to study the effects of the spatial dispersion due to the nonlocality of the graphene's conductivity in [81, Eqs (10)–(12)]). Here, a third-order approximate nonlocal model leads to four extra time-dependent PDEs (cf., [81, Eqs (10)–(12)]). When the thermal energy is much less than the chemical potential of graphene, the surface current ${{\boldsymbol{J}}}(t)$ of graphene becomes a highly nonlinear function of the electric field. In ^[82], a novel FDTD scheme has been proposed to model the nonlinear electrodynamic properties of graphene at THz frequencies.

4.   Perfectly matched layers

In practice, many interesting wave propagation and radiation/scattering problems happen in unbounded domains, and due to limited computer memory we usually have to truncate the unbounded domain to a bounded one by either using the so-called absorbing boundary conditions ^[83,84,85], or PMLs.

The PML concept was originally introduced by Bérenger in 1994 ^[86] for efficiently solving the unbounded electromagnetic problems with the finite-difference time-domain method. The PML idea is to surround the computational domain with a specially designed artificial (nonphysical) absorbing layer of finite thickness to absorb the outgoing waves (leaving the physical domain) without any spurious reflections.

Inspired by Bérenger's work, many different PML models have been proposed and analyzed for solving Maxwell's equations in the past two decades. In ^[87], a Maxwellian material based PML model was developed for the time-dependent Maxwell's equations. Later, we carried out the stability analysis for this PML model in its original form ^[88] and in integro-differential form ^[89]. Finite element methods are proposed and its wave absorbing performances are demonstrated in ^[88,89]. In ^[90], Cohen and Monk developed another PML model for the time-dependent Maxwell's equations by using the stretched coordinates approach. The PML performance is tested by a finite element method with mass-lumped edge elements. Later in ^[91], we estabilshed a stability result for the Cohen-Monk model, and a new finite element scheme was proposed, analyzed, and implemented. In ^[92], a new type of absorbing layer for Maxwell's equations and the linearized Euler equations in the time domain were proposed and analyzed. The PML model for the Maxwell's equations is augmented from the free space Maxwell's equations with six additional unknowns and three damping functions. One major interest from the mathematics community is in the convergence analysis of the PMLs. For time-harmonic electromagnetic scattering problems, the exponential convergence has been established in terms of the thickness of the PML layer for the circular or spherical PML method (e.g., ^{[93,94,95,96,97,98]}) and in ^{[57,15,17,19,20]} for the uniaxial (or Cartesian) PML method (e.g., ^[99,100,101]). On the other hand, there are fewer works on the convergence analysis of PML models in the time domain. The exponential decays and convergence of the PML solutions are established for the one-dimensional time-dependent Maxwell system, acoustic equations and hyperbolic systems in ^[102], which is the first result on exponential decays and convergence of PMLs for time-dependent systems. For the transient electromagnetic scattering problems ^[103,104], similar exponential convergences have been established for the spherical PML method ^[105] and the uniaxial PML method ^[106]. A good survey on PML models for Maxwell's equations can be seen in the reivew article by Teixeira and Chew ^[107].

Due to its excellent performance demonstrated for electromagnetic wave simulations (e.g., ^{[108,109,110,111]}), the PML technique originally proposed for the Maxwell's equations has been quickly extended to solve other wave problems. For example, PML models have been studied for the linearized Euler equations (e.g., ^[112,113]), for elastic wave problems (cf., a recent review paper by Pled and Desceliers ^[114] and references therein), for general hyperbolic systems ^[115], for hyperbolic-parabolic systems ^[116], for the wave equation ^[117], for the 2D Helmholtz equation ^[118], and for the elastodynamic problems ^{[119,120,121]}.

4.1.   Berenger's PML

Consider the two-dimensional Bérenger's $TE_z$ PML model (Transverse Electric with the electric field lying in the $(x, y)$ plane), which involves only three unknowns, $E_x, E_y, H_z$, and is given as [86, (3.a)–(3.d)]: i

where $(E_x, E_y)$ and $H_z = H_{zx}+H_{zy}$ (splitting $H_z$ into two components) are the electric and magnetic fields, respectively; ${\epsilon}_0 = 8.854\cdot 10^{-12}F/m$ and $\mu_0 = 4\pi\cdot 10^{-7}H/m$ are the vacuum permittivity and permeability, respectively; and $\sigma_x, \sigma_x, \sigma_x^*, \sigma_y^*$ are the electric and magnetic conductivities. To avoid the reflection across the vacuum-medium interface, we require the following impedance matching conditions to be satisfied:

The Bérenger's PML model, such as (4.1a)–(4.1d), works very well in practical simulations. However, the mathematical study by Abarbanel and Gottlieb ^[122] pointed out that Bérenger's PML model is weakly well-posed but not strongly well-posed, i.e., it may suffer the long time instability. Later, a careful study by Becache and Joly ^[123] proved that for the Cauchy problem made of Eqs (4.1a)–(4.1d) with $\sigma_y(y) = \sigma_y^*(y) = 0$ (i.e., in the PML layer parallel to the $y$ axis) and ${\epsilon}_0 = \mu_0 = 1$, the following stability [123, Theorem 1.2]:

holds true for any initial function $U^{0}\in (H^1(\mathcal{R}^2))^4$, where $U = (E_x, E_y, H_{zx}, H_{zy})'$ denotes the solution of (4.1a)–(4.1d) with the initial condition

Here, the constant $C_{\sigma} = C\min\left(1, 1/\sigma\right)$ and $\sigma = \sigma_x = \sigma_{x}^{*} > 0$. The stability (4.3) shows that the weak well-posedness appears in the loss of regularity between the initial data and the solution at time $t$. They proved further [123, Theorem 1.3] that the loss of regularity only affects the split fields, while the physical fields ${{\boldsymbol{E}}} = (E_x, E_y)'$ and $H_z = H_{zx}+H_{zy}$ still satisfy the usual estimate, as for the standard Maxwell's equations, i.e.,

where constant $C > 0$ is independent of $\sigma$.

In ^[123], Becache and Joly also used the energy techniques to prove some nice stability results for the PML model based on Zhao-Cangellaris's formulation ^[124] for general positive variable conductivities. It is a pity that their proof of Theorem 2.7 in ^[123] has some typos. Inspired by Becache and Joly's work ^[123], we rederived the TEz PML model (by correcting their typos and putting the physical permittivity and permeability back), which is given as follows [125, Eqs (5)–(9)]: For any $({{\boldsymbol{x}}}, t) \in \Omega\times (0, T]$,

where we denote ${{\boldsymbol{E}}} = (E_x, E_y)'$ and $H$ for the electric field and magnetic field, and $\widetilde{{{\boldsymbol{E}}}} = (\widetilde{E}_x, \widetilde{E}_y)'$, $\widetilde{H}$ and $H^*$ for the auxiliary variables. Moreover, we denote

where $\sigma_x(x), \sigma_y(y)\geq 0$ are the damping functions in the $x, y$ directions, respectively. Here, ${\Omega}$ is assumed to be an open bounded Lipschitz polygon in $\mathcal{R}^2$ with boundary ${\partial}{\Omega}$ and outward unit normal vector ${{\boldsymbol{n}}}$. Furthermore, we assume that the model problem (4.5a)–(4.5e) satisfies the perfect electric conductive (PEC) boundary condition:

and the initial conditions:

where ${{\boldsymbol{E}}}_0, \widetilde{{{\boldsymbol{E}}}}_0, H_0, H^{*}_0$ and $\widetilde{H}_0$ are some given functions.

Following the technique developed by Becache and Joly in ^[123], we established the following stability for the model (4.5a)–(4.5e).

Theorem 2. [125, Theorem 1] Denote the energy at any time $t$ for the solution $({{\boldsymbol{E}}}, \widetilde{{{\boldsymbol{E}}}}, H, \widetilde{H})$ of (4.5a)–(4.5e) subject to the PEC boundary condition as:

Under the PEC boundary condition, we have

Once the continuous stability is establiished, we can develop various numerical methods to solve this PML model (4.5a)–(4.5e). In ^[125], both an FDTD method (with 4th order in space and 2nd order in time convergence) and an FEM using edge elements are developed and analyzed. A novel explicit unconditionally stable finite element scheme using edge elements was proposed and analyzed recently in ^[126].

4.2.   PML models for metamaterials

It is known that the classical Berenger's PML is unstable when it is used to interact with metamaterials directly due to the presence of backward waves. This phenomenon was pointed out in several papers from the physicist community (cf., ^{[127,128,129,130,131]}). Furthermore, they proposed some stable PML models to overcome the unstability. In 2015, Bécache et al. ^[132] systematically derived a new PML model with 16 unknowns to overcome this unstability issue, and they carried out some simulations by the FDTD method to demonstrate the stability of the PML model. Unfortunately they did not provide any details about the specific numerical scheme in ^[132] and their continuous work ^[133].

Inspired by ^[132], recently we proposed, analyzed and implemented both an FDTD method ^[134] and a finite element method ^[135] for a special case of the Drude PML model established by Bécache et al. (cf. [132, Eq (48)] and [134, (2.1)]):

where ${{\boldsymbol{E}}} = (E_x, E_y)$ and $H = H^x+H^y$ are the electric field and magnetic field (in split form) respectively, ${{\boldsymbol{J}}} = (J_x, J_y)$ and ${{\boldsymbol{K}}} = (K^x, K^y)$ are the auxiliary variables, $\sigma_x(x)\geq 0$ and $\sigma_y(y)\geq 0$ are the damping functions in the $x$ and $y$ directions, $\omega_e$ and $\omega_m$ are the electric and the magnetic plasma frequencies given in the following Drude model:

where $\omega$ denotes the general wave frequency. In Figure 3, we presented our simulation result for a transmission problem originally proposed in ^[132] and simulated by the FDTD method without much details. This transmission problem is given between vacuum and a Drude metamaterial surrounded by $Berenger's$ PML and the metamaterial PML, respectively. Snapshots of the magnetic field $H$ obtained by our finite element method ^[135] with edge elements are presented in Figure 3, which is similar to what Bécache et al. obtained by the FDTD method [132, Figure 12].

In ^[134], we established the following energy identity and stabiltiy result for the solution to the PML model (4.9)–(4.16) subject to a PEC boundary condition.

Theorem 3. [134, Theorem 2.1] Define the energy

For any nonnegative functions $\sigma_x(x)$ and $\sigma_y(y)$, we have

When $\sigma_x = \sigma_y = \sigma \geq 0$ (i.e., a positive constant), the energy is decreasing:

We like to remark that the above stability results are established only when the damping functions $\sigma_x = \sigma_y$ are a positive constant. It is still an open issue whether or not the stability holds true when the damping functions vary with spatial variables, which happens in practical applications. Though some nice analytical results are established for PML models in metamaterials or general dispersive media ^{[132,133,136]}, few works have been done for those models in terms of numerical method developments and applications.

5.   Conclusions

Here we reviewed some recent progress on mathematical analysis and numerical simulations for Maxwell's equations in metamaterials, graphene, and perfectly matched layers. Some open issues are mentioned in the review. Actually many PML models used in practice are quite challenging to establish their stabilities. More mathematical works are needed for analyzing those differential equations developed to simulate wave interactions with other complex media such as metasurfaces (e.g., ^{[137,138,139]}) and time-varying media (e.g., ^{[140,141,142,143,144]}).

Acknowledgments

This work was partially supported by NSF grant DMS-2011943. I would like to thank Prof. Jia Zhao's kind invitation for me to be a plenary speaker for the 2nd Northern States Section Meeting held at Utah State University, Logan, UT, April 15-16, 2023. I very much appreciate his kindly hosting me and hard work on editing this special issue. Finally, I am very grateful to the anonymous reviewers for their helpful comments on improving this paper.

Conflict of interest

The author declares no conflicts of interest in this paper.