# Improving Yee's scheme with asymptotic dispersion correction for time-harmonic Maxwell's equations

Pierre-Henri Cocquet, Martin Gander

► To cite this version:

Pierre-Henri Cocquet, Martin Gander. Improving Yee's scheme with asymptotic dispersion correction for time-harmonic Maxwell's equations. Journal of Computational Physics, 2025, pp.114602. 10.1016/j.jcp.2025.114602. hal-05441004

HAL Id: hal-05441004

https://hal.science/hal-05441004v1

Submitted on 5 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Improving Yee's scheme with asymptotic dispersion correction for time-harmonic Maxwell's equations

Pierre-Henri Cocquet

Laboratoire des Sciences pour l'Ingénieur Appliquées à la Mécanique et au génie Électrique (SIAME), Université de Pau et des Pays de l'Adour, Avenue de l'université, Pau, 64000, France

Martin J. Gander

Section de Mathématiques de Genève, Rue du Conseil-Général 9, CP 64, Geneva, 1211, Suisse

**Abstract**

In this paper, we show how to reduce the dispersion error associated to Yee's finite difference scheme applied to time-harmonic Maxwell's equations in one, two and three spatial dimensions. Our method, called asymptotic dispersion correction, is based on the introduction of a shifted angular frequency depending on a free parameter in the Yee stencil. The optimal parameter, called the asymptotically optimal shift, is next explicitly determined by minimizing the dispersion error for small enough meshsize or, equivalently, for large enough number of grid points per wavelength. Numerical experiments are provided and show that the relative error is reduced when using the optimal shifted angular frequency as soon as the number of grid points per wavelength is large enough.

*Keywords:*

Time-harmonic Maxwell's equations, Yee's finite difference method, Dispersion correction.

**1. Introduction**

Maxwell's equations are used to describe the propagation of time-harmonic electromagnetic waves. For its simplicity, the finite difference method introduced by K. Yee in [54] is widely used for its numerical simulation. Like most
numerical approximations of time-harmonic wave propagation problems, the Yee's scheme is subject to dispersion error. The latter is the fact that plane waves at the discrete level propagate at a speed different from the continuous one. This phase difference produces the so-called *pollution effect* (see [2, 20, 32, 36, 42, 57] for references related to the Helmholtz equation and [7, 6, 26, 25] for the time-harmonic Maxwell's equations) which indicates that the relative error grows with the wavenumber $k$ when a fixed number of grid points per wavelength $G = 2\pi/(kh)$ is used while the wavenumber increases. Because of the pollution effect, considering meshes for which the number of grid points per wavelength is large enough, hence $kh \ll 1$, does not ensure the numerical solution to be a good approximation of the continuous one unless the meshsize complies with additional restrictions (e.g. $k^3h^2 \ll 1$ for P1 FEM [4, 57, 25]) resulting in having very large linear systems to be solved.

The pollution effect is currently well-understood for the Helmholtz equation which is a model problem for scalar time-harmonic wave propagation. For instance, if one consider polynomials with degree $q$ such that $q \ge C \log(k)$ for some $C > 0$ and mesh sizes so that $kh/q$ is small enough, then the pollution effect can be overcome (see [36, 42, 57]). In addition, for 1d Helmholtz equation, it is possible to suppress the dispersion error which also removes the pollution effect (see [4, 49, 31, 14]) by suitably modifying the numerical schemes.

We emphasize that, in dimension greater than 1, the dispersion error (and thus the pollution effect as well) can not be suppressed. Nevertheless, there still exist some methods aiming at reducing the dispersion error which will be referred to as *dispersion correction*. It is worth noting that one of the main advantages of performing dispersion correction for a given numerical scheme is that the resulting scheme has smaller relative error compared to their un-corrected version (see e.g. all the subsequent references).
Dispersion correction for the Helmholtz equation is currently widely studied. In [56], the authors show that the dispersion error for hp-CIP-FEM on tensor meshes can be reduced by suitably choosing a penalty parameter. A similar technique has been used in [5] for a second-order IPDG method on structured meshes with equilateral triangles. Second and sixth-order finite-difference stencils have been derived in [43] where some free parameters are then obtained by minimizing phase and amplitude errors. Here the coefficients are sought to be polynomials of $kh$. A slightly similar technique has been used in [9] where a second order 27-point finite difference stencil with free parameters is considered and where the coefficients are then obtained by
minimizing the dispersion error over some intervals of number of grid points per wavelength. The same approach has been derived in [8, 10] for a second order 9-point stencil, in [19] for a 4th-order 17- and a 27-point stencil, in [50] for a 4th-order 9-point stencil and in [52] for a 4th-order 27-point finite difference method. Reducing the dispersion has also been done for a fixed number of grid points per wavelength by considering finite difference stencils with free-parameters that are next obtained by minimizing the phase error in [51, 34, 44] for a second-order 9-point stencil and in [53] for a 6th-order 9-point stencil. In [24, 31, 37], the dispersion error has also been reduced by considering a 9-point stencil with free parameters that are defined as polynomials of $kh$ whose coefficients are determined by minimizing the average truncation error of plane waves. We emphasize here that the truncation order of the stencil applied to plane waves is directly linked to the dispersion error as indicated by [14, Theorem 3.2, Theorem A.1]. We also would like to mention [35] where the dispersion correction is done by factoring the Helmholtz operator with 2 pseudo-differential operators that are then approximated using different approaches (e.g. Padé approximants) to reduce the dispersion error. Finally, in [21], the authors present a correction of the eigenvalues of the finite difference matrices associated to the 5-point (in 2d) and 7-point (in 3d) stencils. Their results indicate that the pollution effect can be greatly reduced using this approach.

Regarding dispersion correction for time-harmonic Maxwell's equation, to the best of our knowledge, we only found 2 references. In [41], the centered finite-difference formula for the first-order derivative is modified though the introduction of some constant that is next slightly perturbed thanks to some numerical optimization (see e.g. [41, Table II]) with the objective of reducing the dispersion error. In [33], the authors considered a least-squares finite element method for the Maxwell's equations where the additional least-squares term depends on some constant. For the lowest order hexahedral Nédélec elements, this constant is next explicitly computed to reduce the dispersion error for small enough meshsize. One drawback of this approach is that the dispersion analysis can not be readily extended to higher order Nédélec elements since, in that case, the additional least-squares term can no longer be seen as a constant times the mass matrix.

According to the above literature review on dispersion correction (or dispersion minimizing) schemes for time harmonic wave propagation problems, it is worth noting that most approaches consider numerical schemes with free parameters that are next determined though numerical optimization by
minimizing the dispersion error. In addition, to the best of our knowledge, the case of time-harmonic Maxwell's equation does not seem to have been considered yet.

In this paper, we would like to address dispersion correction for the Yee finite difference scheme. Since the Yee's stencil does not have any free parameter in it, we are going to rely on the technique we developed in [17, 14, 16, 15] for finite difference methods applied to the Helmholtz equation. Our method, called asymptotic dispersion correction is based first on the introduction of a shifted wavenumber, which depends on some parameter, in the finite difference stencil. Next, we compute explicitly the shift that minimizes the dispersion error either for small enough meshsize or for large enough number of grid points per wavelength. The resulting scheme is thus very similar to the original one but with reduced dispersion error. The major advantage of this approach is that it can be used for numerical schemes for which we do not have access to parameters to optimize with.

In the present paper, we extend our method to Maxwell problems in the following ways : first we are going to apply the asymptotic dispersion correction to full 1d, 2d and 3d Maxwell's equations, we also define asymptotically optimal shifts for media with complex permittivity and permeability that can also depend on the angular frequency. Until now (see e.g. [14]), the asymptotically optimal shift was defined by minimizing the $L^\infty$ dispersion error and we now defined it also for $L^p$-dispersion error, for $p \in [1, +\infty]$, although only the cases $p = 2, +\infty$ permit for closed form formulas. We also show, with numerical experiments, that our explicit formula can be used for piecewise-constant parameters.

This paper is organized as follows: In the next section, we present the Yee finite differences scheme, we then recall some of its standard properties and compute its dispersion relation. Next, we show how to perform the asymptotic dispersion correction using a shifted angular frequency in the Yee scheme. We then give some numerical results to show the influence of the asymptotic dispersion correction on the relative error. We finally end with some conclusions and outlook.

## 2. Maxwell's systems and their dispersion relations

The 3d time-harmonic Maxwell equations in dielectric media are given by

$$ \begin{cases} i\omega\varepsilon\mathbf{E} - \nabla \times \mathbf{H} = -\mathbf{J}, \\ i\omega\mu\mathbf{H} + \nabla \times \mathbf{E} = \mathbf{0}, \end{cases} \qquad (1) $$
where **E** is the electric field, **H** the magnetic field, ε the permittivity, μ the permeability and **J** is the electric current. Let **e**<sub>j</sub> ∈ ℝ<sup>3</sup> be the vectors of the cartesian basis that are defined component-wise as

$$
\mathbf{e}_{j,i} = \delta_{j,i},
$$

where $\delta$ is the Kronecker symbol. We can write

$$
\nabla \times \mathbf{E} = \sum_{j=1}^{3} [\mathbf{e}_j \times] \partial_j \mathbf{E},
$$

where [e_j x] are 3 x 3 matrices related to the cross product x of some vector
X ∈ C^3 as [e_j x] X = e_j x X.

Since this paper deals with dispersion correction for the Maxwell’s equa-
tion which is a system of first order partial differential equations, we recall
below how the dispersion relation and the scalar wavenumber are defined.
Consider the first-order differential operator with constant coefficients

$$
\mathcal{L}(\omega; \nabla) := i\omega A_0 + \sum_{j=1}^{d} A_j \partial_j,
$$

where $A_0, A_j$ are square matrices. For any $\mathbf{k} \in \mathbb{C}^d$, its symbol is defined as

$$
\mathcal{L}(\omega, \mathbf{i}\mathbf{k}) = i\omega A_0 + \sum_{j=1}^{d} \mathbf{i}k_j A_j,
$$

and the dispersion relation is then

$$
D_c := \{ \mathbf{k} \in \mathbb{C}^d \mid \det (\mathcal{L}(\omega, \mathbf{i}\mathbf{k})) = 0 \}.
$$

In addition, for any wave-vector **k** ∈ *D*<sub>*c*</sub> there is some **θ** ∈ *S*<sup>*d*−1</sup> so that
**k** = ||**k**|| **θ** and the scalar wavenumber *k*<sub>s</sub>(*ω*, **θ**) is then given as

$$
k_s (\omega, \boldsymbol{\theta}) = \| \boldsymbol{k} \| .
$$

Below, we give the Maxwell’s equations for one, two and three dimensions
and also compute their dispersion relation.
**Remark 2.1.** When working on a bounded domain $\Omega \subset \mathbb{R}^3$ with outward unit normal $\mathbf{n}$, the Maxwell's equation (1) is supplemented with boundary conditions like e.g.

Perfect Electric Conductor : $\mathbf{E} \times \mathbf{n} = 0$ on $\partial\Omega$,

Impedance : $\mathbf{n} \times \mathbf{E} - \mathbf{n} \times \Lambda(x) (\mathbf{n} \times \mathbf{H}) = 0$ on $\partial\Omega$,

for some symmetric positive definite matrix $\Lambda(x)$. For all $\omega > 0$ and any $\mathbf{J} \in L^2(\Omega)$, the existence and uniqueness of a solution, $(\mathbf{E}, \mathbf{H}) \in L^2(\Omega)^6$ with $(\nabla \times \mathbf{E}, \nabla \times \mathbf{H}) \in L^2(\Omega)^6$, to the Maxwell's equation with one of the above boundary conditions has been proved for example in [38] (see also [18]). These results hold for coercive dielectric parameters $\varepsilon, \mu \in L^\infty(\Omega, \mathbb{C}^{3 \times 3})$ and for dispersive media, hence having frequency dependent parameters, as well. These results also hold for 1d and 2d Maxwell's equations.

## 2.1. 1d Maxwell system

The one-dimensional Maxwell's equations can be obtained (see e.g. [28]) from (1) by considering electric and magnetic fields given by

$$ \mathbf{E} = (0, E_2, 0)^T, \quad \mathbf{H} = (H_1, 0, 0)^T, $$

where $E_2$ and $H_1$ only depend on the $z$ variable. Using these simplifications in (1) yields

$$ \left\{ \begin{array}{l} i\omega\varepsilon E_2 - \partial_z H_1 = -J_2, \\ i\omega\mu H_1 - \partial_z E_2 = 0, \end{array} \right. \qquad (2) $$

which will be referred to as the 1d Maxwell's equation.

Its continuous symbol is

$$ \mathcal{L}^{1d}(k_3) = i\omega \begin{pmatrix} \varepsilon & 0 \\ 0 & \mu \end{pmatrix} + ik_3 \begin{pmatrix} 0 & -1 \\ -1 & 0 \end{pmatrix}, $$

so that the dispersion relation is

$$ D_c^{1d} = \{k_3 \in \mathbb{C} \mid -\omega^2\varepsilon\mu + k_3^2 = 0\}. $$

The scalar wavenumber (being also the wave vector for one-dimensional problems) is thus

$$ k_s^{1d} = \omega\sqrt{\varepsilon\mu}. $$

We emphasize that there exist other ways to derive the one-dimensional Maxwell equation. Note nevertheless that their dispersion properties are exactly the same as those of (2).
## 2.2. 2d Maxwell's equations: Transverse Electric (TE) and Transverse Magnetic (TM) modes

The transverse electric Maxwell equations can be obtained from (1) by considering

$$ \mathbf{E} = (0, 0, E_3)^T, \quad \mathbf{H} = (H_1, H_2, 0)^T, $$

where each component depends on $(x, y)$. The resulting system of equations is thus

$$ \left\{ \begin{array}{l} i\omega\varepsilon E_3 + \partial_x H_2 - \partial_y H_1 = J_3, \\ i\omega\mu H_1 - \partial_y E_3 = 0, \\ i\omega\mu H_2 + \partial_x E_3 = 0. \end{array} \right. \qquad (3) $$

Its continuous symbol is

$$ \mathcal{L}^{2d}(k_1, k_2) = i\omega \begin{pmatrix} \varepsilon & 0 & 0 \\ 0 & \mu & 0 \\ 0 & 0 & \mu \end{pmatrix} + ik_1 \begin{pmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \end{pmatrix} + ik_2 \begin{pmatrix} 0 & -1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}. $$

The dispersion relation is

$$ D_c^{2d} = \{ (k_1, k_2) \in \mathbb{C}^2 \mid i\omega\mu(-\varepsilon\mu\omega^2 + \|\mathbf{k}\|^2) = 0 \}, $$

and the scalar wavenumber is thus

$$ k_s^{2d} = \|\mathbf{k}\| = \omega\sqrt{\varepsilon\mu}. $$

**Remark 2.2.** Some transverse-magnetic (TM) Maxwell's equations can also be defined from (1) by considering

$$ \mathbf{E} = (E_1, E_2, 0)^T, \quad \mathbf{H} = (0, 0, H_3)^T, $$

where each components depend on $(x, y)$. This gives

$$ \left\{ \begin{array}{l} i\omega\varepsilon E_1 - \partial_y H_3 = -J_1, \\ i\omega\varepsilon E_2 + \partial_y H_3 = -J_2, \\ i\omega\mu H_3 + \partial_x E_2 - \partial_y E_1 = 0. \end{array} \right. \qquad (4) $$

Some computations can then show that the dispersion relation of (4) is the same as the one of (3).
### 2.3. 3d Maxwell's equations

The continuous symbol of the 3d-Maxwell's equations (1) is

$$ \mathcal{L}(\omega, \mathbf{k}) = i\omega \begin{pmatrix} \varepsilon\mathbb{I} & 0 \\ 0 & \mu\mathbb{I} \end{pmatrix} + \sum_{j=1}^{3} i k_j \begin{pmatrix} 0 & [\mathbf{e}_j \times] \\ -[\mathbf{e}_j \times] & 0 \end{pmatrix}. $$

Since

$$ \det(\mathcal{L}(\omega, \mathbf{k})) = -\omega^2 \varepsilon \mu (\varepsilon \mu \omega^2 - \|\mathbf{k}\|^2)^2, $$

the dispersion relation is

$$ D_c^{3d} = \{(k_1, k_2) \in \mathbb{C}^2 \mid \varepsilon \mu \omega^2 - \|\mathbf{k}\|^2 = 0\}, $$

and the scalar wavenumber is thus

$$ k_s = \|\mathbf{k}\| = \omega \sqrt{\varepsilon \mu}. \qquad (5) $$

We emphasize that $k_s$ actually does not depend on $\theta$.

**Remark 2.3.** For $\varepsilon, \mu \in \mathbb{C}$, the wave-vector $\mathbf{k} \in D_c$ is also a complex vector. Nevertheless, the scalar wavenumber can still be defined by (5) where we use the principal square root of a complex number.

## 3. The Yee scheme and its dispersion properties

The Yee scheme has been introduced in [54] and has then been extensively studied (see for example the review [45], [11, Chapter 37] and [39] for the convergence when applied to time-domain Maxwell's equations).

We consider here a uniform grid with meshsize $h$. Writing the first three equations of (1) as $\nabla \times \mathbf{H} - i\omega\varepsilon\mathbf{E} - \mathbf{J} = 0$, the Yee finite difference approximation at $(x, y, z)$ reads as (see e.g. [46, Eq. (2.8)] for a similar stencil
involving only integer index)

$$
\begin{align*}
& \frac{1}{h} \left( H_3 \left( x + \frac{h}{2}, y + \frac{h}{2}, z \right) - H_3 \left( x + \frac{h}{2}, y - \frac{h}{2}, z \right) \right) \\
& + \frac{1}{h} \left( H_2 \left( x + \frac{h}{2}, y, z - \frac{h}{2} \right) - H_2 \left( x + \frac{h}{2}, y, z + \frac{h}{2} \right) \right) \\
& - i\omega\varepsilon E_1 \left( x + \frac{h}{2}, y, z \right) - J_1 \left( x + \frac{h}{2}, y, z \right) = 0, \\
& \frac{1}{h} \left( H_1 \left( x, y + \frac{h}{2}, z + \frac{h}{2} \right) - H_1 \left( x, y + \frac{h}{2}, z - \frac{h}{2} \right) \right) \\
& + \frac{1}{h} \left( H_3 \left( x - \frac{h}{2}, y + \frac{h}{2}, z \right) - H_3 \left( x + \frac{h}{2}, y + \frac{h}{2}, z \right) \right) \\
& - i\omega\varepsilon E_2 \left( x, y + \frac{h}{2}, z \right) - J_2 \left( x, y + \frac{h}{2}, z \right) = 0, \\
& \frac{1}{h} \left( H_2 \left( x + \frac{h}{2}, y, z + \frac{h}{2} \right) - H_2 \left( x - \frac{h}{2}, y, z + \frac{h}{2} \right) \right) \\
& + \frac{1}{h} \left( H_1 \left( x, y - \frac{h}{2}, z + \frac{h}{2} \right) - H_1 \left( x, y + \frac{h}{2}, z + \frac{h}{2} \right) \right) \\
& - i\omega\varepsilon E_3 \left( x, y, z + \frac{h}{2} \right) - J_3 \left( x, y, z + \frac{h}{2} \right) = 0.
\end{align*}
$$

The second part of (1) is $\nabla \times \mathbf{E} + i\omega\mu\mathbf{H} = 0$ and the Yee finite difference
scheme at $(x, y, z)$ is given as

$$ 
\begin{aligned}
& \frac{1}{h} \left( E_2 \left( x, y + \frac{h}{2}, z \right) - E_2 \left( x, y + \frac{h}{2}, z + h \right) \right) \\
& + \frac{1}{h} \left( E_3 \left( x, y + h, z + \frac{h}{2} \right) - E_3 \left( x, y, z + \frac{h}{2} \right) \right) \\
& + i\omega\mu H_1 \left( x, y + \frac{h}{2}, z + \frac{h}{2} \right) = 0, \\
& \frac{1}{h} \left( E_3 \left( x, y, z + \frac{h}{2} \right) - E_3 \left( x + h, y, z + \frac{h}{2} \right) \right) \\
& + \frac{1}{h} \left( E_1 \left( x + \frac{h}{2}, y, z + h \right) - E_1 \left( x + \frac{h}{2}, y, z \right) \right) \\
& + i\omega\mu H_2 \left( x + \frac{h}{2}, y, z + \frac{h}{2} \right) = 0, \\
& \frac{1}{h} \left( E_1 \left( x + \frac{h}{2}, y, z \right) - E_1 \left( x + \frac{h}{2}, y + h, z \right) \right) \\
& + \frac{1}{h} \left( E_2 \left( x + h, y + \frac{h}{2}, z \right) - E_2 \left( x, y + \frac{h}{2}, z \right) \right) \\
& + i\omega\mu H_3 \left( x + \frac{h}{2}, y + \frac{h}{2}, z \right) = 0.
\end{aligned}
$$

In the following subsections, we review first some classical properties of the Yee scheme for the 3d-Maxwell's equations. It is worth noting that these properties can also be obtained for the 1d and 2d Maxwell's equations since the latter are derived from the 3d system by setting specific components of both the electric and magnetic fields to zero. After this, we will compute the discrete wavenumbers associated to the Yee scheme for the 1d, 2d and 3d Maxwell system.

### 3.1. Classical properties of the Yee scheme

We review here some properties of the Yee scheme. First, we introduce the forward and backward differential operators in the x-direction,

$$ \delta_h^{+, (1)} \varphi = \frac{1}{h} (\varphi_{i+1,j,k} - \varphi_{i,j,k}), \quad \delta_h^{-, (1)} \varphi = \frac{1}{h} (\varphi_{i,j,k} - \varphi_{i-1,j,k}). $$

The discrete derivatives can be defined similarly for the y and z directions and are denoted by $\delta_h^{\pm, (2)}$ and $\delta_h^{\pm, (3)}$. As a result, the forward and backward
discrete gradients at a grid point $(x_i, y_j, z_k)$ are

$$ \nabla_h^+ \varphi = \frac{1}{h} \left( \begin{array}{l} \varphi_{i+1,j,k} - \varphi_{i,j,k} \\ \varphi_{i,j+1,k} - \varphi_{i,j,k} \\ \varphi_{i,j,k+1} - \varphi_{i,j,k} \end{array} \right) = \sum_{j=1}^{3} \delta_h^{+, (j)} \varphi e_j, $$

$$ \nabla_h^- \varphi = \frac{1}{h} \left( \begin{array}{l} \varphi_{i,j,k} - \varphi_{i-1,j,k} \\ \varphi_{i,j,k} - \varphi_{i,j-1,k} \\ \varphi_{i,j,k} - \varphi_{i,j,k-1} \end{array} \right) = \sum_{j=1}^{3} \delta_h^{-, (j)} \varphi e_j. $$

Note that we do not highlight the dependence of these discrete operators on the grid point, but these formulas are valid for any $\varphi = (\varphi_{i,j,k})_{i,j,k \in \mathbb{Z}}$. We can then define the discrete divergence and curl operators of a vector field $\mathbf{\Phi} = (\mathbf{\Phi}_1, \mathbf{\Phi}_2, \mathbf{\Phi}_3)$ as

$$ \mathrm{div}_{h}^{\pm} \mathbf{\Phi} = \nabla_{h}^{\pm} \cdot \mathbf{\Phi} := \sum_{j=1}^{3} \delta_{h}^{\pm, (j)} \mathbf{\Phi}_{j}, $$

$$ \mathrm{curl}_{h}^{\pm} \mathbf{\Phi} = \nabla_{h}^{\pm} \times \mathbf{\Phi} := \sum_{j=1}^{3} [e_j \times] \delta_{h}^{\pm, (j)} \mathbf{\Phi}, $$

where everything is computed at a grid point $(x_i, y_j, z_k)$. Introducing the *shifting* operators

$$ (S_h E)_{i,j,k} = E_{1;i+\frac{1}{2},j,k} e_1 + E_{2;i,j+\frac{1}{2},k} e_2 + E_{3;i,j,k+\frac{1}{2}} e_3, $$

$$ (\tilde{S}_h H)_{i,j,k} = H_{1;i+\frac{1}{2},j,k+\frac{1}{2}} e_1 + H_{2;i+\frac{1}{2},j,k+\frac{1}{2}} e_2 + H_{3;i+\frac{1}{2},j+\frac{1}{2},k} e_3, $$

the Yee scheme can then be written at each grid point $(x_i, y_j, z_k)$, in compact form, as

$$ i\omega\varepsilon S_h E - \nabla_h^- \times \tilde{S}_h H = -S_h J, \quad (6) $$

$$ i\omega\mu \tilde{S}_h H + \nabla_h^+ \times S_h E = 0. $$

The discrete differential operators also have some properties verified by the continuous operators.
**Theorem 3.1.** At each grid point, we have<fnref n="1" />

$$ \mathrm{div}_{h}^{\pm}(\nabla_{h}^{\pm} \times \Phi) = 0, $$

$$ \nabla_{h}^{\pm} \times (\nabla_{h}^{\pm} \varphi) = 0, $$

$$ \nabla_{h}^{+} \times (\nabla_{h}^{-} \times \Phi) = \nabla_{h}^{-} (\mathrm{div}_{h}^{+} \Phi) - \Delta_{h} \Phi, $$

$$ \nabla_{h}^{-} \times (\nabla_{h}^{+} \times \Phi) = \nabla_{h}^{+} (\mathrm{div}_{h}^{-} \Phi) - \Delta_{h} \Phi, $$

where $\Delta_h \Phi = \sum_j \Delta_h \Phi_j e_j$ is the discrete vector Laplace operator defined as $\Delta_h = \mathrm{div}_h^+ \nabla_h^- = \mathrm{div}_h^- \nabla_h^+$.

The properties given in Theorem 3.1 can be checked by direct computations. In addition, we can see that

$$ \Delta_h \varphi = \frac{\varphi_{i+1,j,k} + \varphi_{i-1,j,k} + \varphi_{i,j,k+1} + \varphi_{i,j,k-1} + \varphi_{i,j+1,k} + \varphi_{i,j-1,k} - 6\varphi_{i,j,k}}{h^2}, $$

which is thus the standard 7-point stencil for the 3d-Laplace operator. Using Theorem 3.1, we can show that any solution to (6) also satisfies, at each grid point, the discrete vector Helmholtz equations

$$ \begin{aligned} -\Delta_h S_h E - \omega^2 \varepsilon \mu S_h E &= -i\omega\mu S_h J + \frac{1}{i\omega\varepsilon} \nabla_h^+ (\mathrm{div}_h^- S_h J), \\ -\Delta_h \tilde{S}_h H - \omega^2 \varepsilon \mu \tilde{S}_h H &= \nabla_h^+ \times (S_h J), \end{aligned} \quad (7) $$

$$ \mathrm{div}_h^+ (\tilde{S}_h H) = 0, $$

$$ i\omega\varepsilon\mathrm{div}_h^- (S_h E) = -\mathrm{div}_h^- (\tilde{S}_h J), $$

which is the discrete counterpart of the vector wave equations satisfied by the electric and magnetic fields at the continuous level. We end this section by noting that, in what follows, we will not use the second order formulation of the continuous and the discrete Maxwell's equations. Nevertheless, it is worth noting that (7) shows the link between the Yee scheme and the Helmholtz equation discretized using the 7-point stencil for the 3d-Laplace operator. We can therefore expect the Yee scheme has dispersion properties that are similar to those of this stencil.

Regarding the truncation order of the Yee scheme, although it is defined using only first order finite difference operators, the resulting scheme is nevertheless second order accurate as one can prove using Taylor expansion.

**Theorem 3.2.** The Yee's scheme has second order truncation error.

<footnote n="1">¹These properties are also known as discrete De Rham complex.</footnote>
### 3.2. Dispersion analysis for the Yee scheme applied to the 1d Maxwell equation

Starting from the 3d-Maxwell's equations (1), its 1d counterpart (2) can be derived by setting

$$ \mathbf{E} = (0, E_2, 0)^T, \quad \mathbf{H} = (H_1, 0, 0)^T, $$

which also works for the Yee scheme for the 1d-Maxwell's equations. As a result, at any grid points $x_j = jh$, the 1d Yee scheme is

$$ \left\{ \begin{array}{l} i\omega\varepsilon E_j - \frac{1}{h}(H_{j+1/2} - H_{j-1/2}) = -J(x_j), \\ i\omega\mu H_{j+1/2} - \frac{1}{h}(E_{j+1} - E_j) = 0. \end{array} \right. \qquad (8) $$

We next show how to compute the discrete symbol associated with the 1d Yee stencil given by

$$ \left( \mathcal{M}(h) \begin{pmatrix} E \\ H \end{pmatrix} \right)_i = \left( \begin{array}{l} i\omega\varepsilon E_j - \frac{1}{h}(H_{j+1/2} - H_{j-1/2}) \\ i\omega\mu H_{j+1/2} - \frac{1}{h}(E_{j+1} - E_j) \end{array} \right). $$

Inserting a plane wave $Xe^{-ixk}$ into the stencil yields

$$ \left( \mathcal{M}(h) \begin{pmatrix} X_1 \\ X_2 \end{pmatrix} e^{-ixk} \right)_i = \left( \begin{array}{l} i\omega\varepsilon X_1 e^{-ix_j k} - X_2 \frac{2i}{h} (e^{-ix_{j+1/2} k} - e^{-ix_{j-1/2} k}) \\ i\omega\mu X_2 e^{-ix_{j+1/2} k} - X_1 \frac{2i}{h} (e^{-ix_{j+1} k} - e^{-ix_{j} k}) \end{array} \right). $$

Multiplying by the diagonal matrix $D_i = \text{diag}(e^{ix_j k}, e^{ix_{j+1/2} k})$, we can remove the dependence on the grid point where the approximation is computed. We obtain

$$ D_i \left( \mathcal{M}(h) \begin{pmatrix} X_1 \\ X_2 \end{pmatrix} e^{-ixk} \right)_i = \left( \begin{array}{l} i\omega\varepsilon X_1 - X_2 \frac{2i}{h} \sin\left(\frac{kh}{2}\right) \\ i\omega\mu X_2 - X_1 \frac{2i}{h} \sin\left(\frac{kh}{2}\right) \end{array} \right) = \sigma_{\mathcal{M}}(h; k) \mathbf{X}, $$

where

$$ \sigma_{\mathcal{M}}(h, k) = \begin{pmatrix} i\omega\varepsilon & -2i\frac{\sin\left(\frac{kh}{2}\right)}{h} \\ -2i\frac{\sin\left(\frac{kh}{2}\right)}{h} & i\omega\mu \end{pmatrix} $$

is the discrete symbol. Computing the determinant of the discrete symbol shows that the discrete dispersion relation is

$$ \left( -\omega^2\varepsilon\mu + \frac{4}{h^2}\sin\left(\frac{kh}{2}\right)^2 \right) = 0. \qquad (9) $$
The discrete wavenumber is therefore

$$k_d(\omega; h) = \frac{2}{h} \arcsin \left( \frac{\omega \sqrt{\varepsilon \mu h}}{2} \right), \qquad (10)$$

where $\arcsin$ here is the principal branch of the inverse sine function that is defined for $\omega\sqrt{\varepsilon\mu h}/2 \in \mathbb{C} \setminus ((-\infty, -1) \cup (1, +\infty))$. Note that the discrete wavenumber can therefore be expressed with (10) as soon as $|\omega\sqrt{\varepsilon\mu h}/2| \le 1$, which amounts to consider $h$ small enough, or equivalently a number of grid points per wavelength such that $G \ge \pi$. A Taylor expansion for $\omega\sqrt{\varepsilon\mu h} = k_s h$ small gives

$$k_d(\omega; h) = k_s + \frac{k_s}{24} (k_s h)^2 + O((k_s h)^4), \qquad (11)$$

and we then get that the discrete wavenumber is a second order approximation of the continuous one.

### 3.3. Dispersion analysis of the Yee scheme for 2d Maxwell equation

The two-dimensional TE-Maxwell's equation (3) are derived from (1) by considering

$$\mathbf{E} = (0, 0, E_3)^T, \quad \mathbf{H} = (H_1, H_2, 0)^T,$$

and the 2d-Yee's scheme is then

$$\left\{ \begin{array}{l} i\omega\varepsilon E_{3,i,j} + \frac{1}{h}(H_{2,i+1/2,j} - H_{2,i-1/2,j}) - \frac{1}{h}(H_{1,i,j+1/2} - H_{1,i,j-1/2}) = -J_{i,j}, \\ i\omega\mu H_{1,i,j+1/2} - \frac{1}{h}(E_{3,i,j+1} - E_{3,i,j}) = 0, \\ i\omega\mu H_{2,i+1/2,j} + \frac{1}{h}(E_{3,i+1,j} - E_{3,i,j}) = 0. \end{array} \right. \qquad (12)$$

The discrete symbol is now

$$\sigma_M(h, \mathbf{k}) = \begin{pmatrix} i\omega\varepsilon & -\frac{2i}{h}\sin\left(\frac{k_2h}{2}\right) & \frac{2i}{h}\sin\left(\frac{k_1h}{2}\right) \\ -\frac{2i}{h}\sin\left(\frac{k_2h}{2}\right) & i\omega\mu & 0 \\ \frac{2i}{h}\sin\left(\frac{k_1h}{2}\right) & 0 & i\omega\mu, \end{pmatrix},$$

and the discrete dispersion relation is thus

$$\mathcal{D}_d^{2d} := \left\{ \mathbf{k} \in \mathbb{C}^2 \mid -\omega^2 \varepsilon \mu + \frac{4}{h^2} \left( \sin\left(\frac{k_1 h}{2}\right)^2 + \sin\left(\frac{k_2 h}{2}\right)^2 \right) = 0 \right\}. \quad (13)$$

For any $\theta = (\cos(\theta), \sin(\theta)) \in S^1$, we can write $\mathbf{k} \in \mathcal{D}_d^{2d}$ as $\mathbf{k} = k_d(\omega; \theta, h)\theta$ to get the discrete wavenumber $k_d$. Note that it is difficult to find an explicit expression for $k_d$. Nevertheless, we can still compute an approximate
expression as $k_s h$ goes to 0 assuming it admits the expansion

$$k_d(\omega; \theta, h) = k_0 + (k_s h)k_1 + (k_s h)^2 k_2 + \cdots$$

Inserting this ansats into the dispersion relation and identifying terms with same degree in $k_s h$ yields that

$$k_d(\omega; \theta, h) = k_s + \frac{k_s}{24} (\cos(s)^4 + \sin(s)^4) (k_s h)^2 + O((k_s h)^4), \quad (14)$$

and, as in the 1d-setting, the discrete wavenumber is again a second order approximation of the continuous one.

### 3.4. Dispersion analysis of the Yee scheme for 3d Maxwell equation

The discrete symbol associated to the 3d-Yee scheme is given as (see Appendix B for computations using Maple)

$$\sigma_M(h; \mathbf{k}) = \mathrm{diag} \left( \begin{array}{c} i\omega\epsilon e^{i\mathbf{h}k_1/2} \\ i\omega\epsilon e^{i\mathbf{h}k_2/2} \\ i\omega\epsilon e^{i\mathbf{h}k_3/2} \\ i\omega\mu e^{i\omega\mathbf{h}k_2/2}e^{i\omega\mathbf{h}k_3/2} \\ i\omega\mu e^{i\omega\mathbf{h}k_1/2}e^{i\omega\mathbf{h}k_2/2} \\ i\omega\mu e^{i\omega\mathbf{h}k_1/2}e^{i\omega\mathbf{h}k_2/2} \end{array} \right) + \begin{pmatrix} 0 & N_h \\ M_h & 0 \end{pmatrix},$$

where

$$N_h = \begin{pmatrix} 0 & 2i \frac{\sin(hk_3/2)}{h} e^{ihk_1/2} & -2i \frac{\sin(hk_2/2)}{h} e^{ihk_1/2} \\ -2i \frac{\sin(hk_3/2)}{h} e^{ihk_2/2} & 0 & 2i \frac{\sin(hk_1/2)}{h} e^{ihk_2/2} \\ 2i \frac{\sin(hk_2/2)}{h} e^{ihk_3/2} & -2i \frac{\sin(hk_1/2)}{h} e^{ihk_3/2} & 0 \end{pmatrix},$$

and

$$M_h = \begin{pmatrix} 0 & -2i \frac{\sin(\omega h k_3/2)}{h} e^{i\omega h k_2/2} e^{i\omega h k_3/2} & 2i \frac{\sin(\omega h k_2/2)}{h} e^{i\omega h k_2/2} e^{i\omega h k_3/2} \\ 2i \frac{\sin(\omega h k_3/2)}{h} e^{i\omega h k_1/2} e^{i\omega h k_3/2} & 0 & -2i \frac{\sin(\omega h k_1/2)}{h} e^{i\omega h k_1/2} e^{i\omega h k_3/2} \\ -2i \frac{\sin(\omega h k_2/2)}{h} e^{i\omega h k_1/2} e^{i\omega h k_2/2} & 2i \frac{\sin(\omega h k_1/2)}{h} e^{i\omega h k_1/2} e^{i\omega h k_2/2} & 0 \end{pmatrix}.$$

The determinant of $\sigma_M(h; \mathbf{k})$ can be computed using Maple (see Appendix B) and is

$$\det(\sigma_M(h; \mathbf{k})) = -\frac{\omega^2 \mu \epsilon e^{3i\mathbf{h}(k_1+k_2+k_3)/2}}{h^4} \times \left( h^2 w^2 \epsilon \mu - 4 \sin\left(\frac{k_3 h}{2}\right)^2 - 4 \sin\left(\frac{k_2 h}{2}\right)^2 - 4 \sin\left(\frac{k_1 h}{2}\right)^2 \right)^2,$$
where $\mathbf{k} = (k_1, k_2, k_3)$. The discrete dispersion relation is then defined by

$$ \mathcal{D}_d^{3d} = \left\{ \mathbf{k} \in \mathbb{C}^3 \mid h^2 w^2 \varepsilon \mu - 4 \sin \left( \frac{k_3 h}{2} \right)^2 - 4 \sin \left( \frac{k_2 h}{2} \right)^2 - 4 \sin \left( \frac{k_1 h}{2} \right)^2 = 0 \right\}. $$

**Remark 3.3.** We emphasize that the discrete dispersion relation $\mathcal{D}_d^{3d}$ is actually the discrete dispersion relation associated to the discretization, using the 7-point stencil for the Laplace operator, of the (scalar) Helmholtz equation $-\Delta - k^2$ with $k = \omega\sqrt{\varepsilon\mu}$. This can be explained by noting that the solutions to the discretized Maxwell's equations also satisfy the vector Helmholtz equation (7) where each component is discretized using the 7-point stencil.

It is also worth noting that the same holds for the dispersion relation associated to the 1d and 2d Yee scheme. Indeed, the dispersion relation for the 1d Yee scheme (9) is the same as the dispersion relation of the 3-point stencil and, for the 2d Yee scheme, the discrete dispersion relation (13) is the same as the one of the 5-point stencil.

The discrete wavenumber is defined for each $\boldsymbol{\theta}$ by $\mathbf{k} = k_d(\omega; \boldsymbol{\theta}, h)\boldsymbol{\theta}$ where any $\boldsymbol{\theta} \in S^2$ is expressed as

$$ \boldsymbol{\theta} = \begin{pmatrix} \cos(\phi) \sin(\theta) \\ \sin(\phi) \sin(\theta) \\ \cos(\theta) \end{pmatrix}, \quad \phi \in [0, 2\pi], \theta \in [0, \pi]. $$

Since the equation satisfied by $k_d$ seems to hardly have a closed form solution, we look for $k_d$ such that

$$ k_d = k_0 + (k_s h)^2 k_2 + \dots \quad (15) $$

We emphasize that $k_1 = 0$ since $\det(\sigma_M(-h; \mathbf{k})) = \det(\sigma_M(h; \mathbf{k}))$. Inserting (15) into the discrete dispersion relation and identifying each term, we obtain (see Appendix C for computations using Maple)

$$ k_0(\omega) = \omega\sqrt{\varepsilon\mu}, \quad k_2(\omega, \boldsymbol{\theta}) = \frac{\omega\sqrt{\varepsilon\mu}}{24} F(\phi, \theta), $$

where $F: (\phi, \theta) \in [0, 2\pi] \times [0, \pi] \mapsto F(\phi, \theta) \in \mathbb{R}$ is defined as

$$ F(\phi, \theta) = \cos(\phi)^4 \sin(\theta)^4 + \sin(\phi)^4 \sin(\theta)^4 + \cos(\theta)^4. \quad (16) $$
The discrete wavenumber then admits, for $h\omega\sqrt{\varepsilon\mu} = k_s h$ small, the expansion

$$k_d(\omega; \boldsymbol{\theta}, h) = \omega\sqrt{\varepsilon\mu} + (k_s h)^2 \frac{k_s}{24} F(\phi, \theta) + \dots \quad (17)$$

where the dots represent higher order terms. In addition, the dispersion error satisfies

$$\max_{\phi, \theta} \left| \frac{k_d(\omega; \boldsymbol{\theta}, h) - k_s}{k_s} \right| = O(k_s^2 h^2), \quad (18)$$

and thus the discrete wavenumber is a second order approximation of the continuous one. We emphasize that (18) can be recast as

$$\max_{\phi, \theta} \left| \frac{k_d(\omega; \boldsymbol{\theta}, h) - k_s}{k_s} \right| = O((k_s h)^2) = O(G^{-2}),$$

where $G = 2\pi/(k_s h)$ is the number of grid points per wavelength. Finally, we emphasize that an estimate similar to (18) (see (11), (14)) also holds for the 1d and 2d Yee scheme.

## 4. Asymptotic dispersion correction with the shifted angular frequency

We introduce now the dispersion correction for the Yee scheme. The proposed approach is based on our previous work (see [15, 16, 17, 14]) related to dispersion correction for finite difference schemes for Helmholtz equation. The starting point is to introduce, in the scheme, a *shifted angular frequency* $\hat{\omega}$ which depends on a free parameter that is next determined by minimizing the dispersion error for small enough meshsize $h$ (or identically for small enough $k_s h$ or even for a large enough grid points per wavelength given by $G = 2\pi/|k_s h|$). Since our method reduces the dispersion error for small enough meshsize only, we called it *asymptotic dispersion correction*.

Since the discrete wavenumber is a second order approximation of the continuous one (see (11),(14),(17)), the shifted angular frequency will be defined as

$$\hat{\omega} = \omega + h^2\omega_2.$$

The shifted angular frequency being a second order perturbation of $\omega$, it is worth noting that the Yee scheme with shifted angular frequency will still be second order accurate owing to Theorem 3.2. With the previously computed
expansions of the discrete wavenumbers (see (11),(14),(17)), they can all been written as

$$k_d(\omega; \boldsymbol{\theta}, h) = k_s(\omega) + (k_s(\omega) h)^2 k_s(\omega) F(\boldsymbol{\theta}) + \dots, \quad (19)$$

for some specific functions $F$. Using this expansion, we can obtain a second order expansion of the discrete wavenumber with shifted angular frequency

$$\begin{aligned} \hat{k}_d(\omega; \boldsymbol{\theta}, h) &:= k_d(\hat{\omega}; \boldsymbol{\theta}, h) \\ &= k_s(\omega) + (k_s(\omega) h)^2 \left( k_s(\omega) F(\boldsymbol{\theta}) + \omega_2 \frac{\partial_{\omega} k_s(\omega)}{k_s(\omega)^2} \right) + \dots. \end{aligned}$$

The asymptotically optimal shift $\omega_2$ can then be defined as

$$\omega_2^{*,p} := \arg \min_{\omega_2 \in \mathbb{C}} \left\| k_s(\omega) F(\boldsymbol{\theta}) + \omega_2 \frac{\partial_{\omega} k_s(\omega)}{k_s(\omega)^2} \right\|_{L^p(\Omega)}, \quad (20)$$

where $\partial_{\omega} k_s(\omega) = \sqrt{\varepsilon\mu}$ for material with dielectric parameters that do not depend on $\omega$, $\Omega = [0, 2\pi]$ in 2d and $\Omega = [0, 2\pi] \times [0, \pi]$ in 3d. We emphasize that $\omega_2^{*,p}$ aims at reducing the dispersion error for small enough meshsize (or equivalently for a large enough number of grid points per wavelength). It is worth noting that we could have also defined the shifted angular frequency as $\hat{\omega} = \omega + (k_s h)^2 \tilde{\omega}_2$ or as $\hat{\omega} = \omega + G^{-2} \tilde{\omega}_2$. Since $\tilde{\omega}_2$ and $\omega_2$ are defined through $\hat{\omega}$, they are related to each other by $(k_s h)^2 \tilde{\omega}_2 = G^{-2} \omega_2^*$, and in these cases, we can get the asymptotically optimal shift directly from $\omega_2^{*,p}$.

This section is dedicated to compute exact formulas for the asymptotically optimized shifts. We begin with some general formula for the solution to (20). Next, we show how the case of materials with frequency-dependent parameters can be reduced to the case of materials with $\varepsilon, \mu$ that do not depend on $\omega$. We then compute the exact formula of the asymptotically optimal shift for the 1d-Yee scheme for which we can actually suppress the dispersion error. Finally, we compute the shift for the 2d and 3d Yee scheme.

### 4.1. Main results for the asymptotic dispersion correction

As seen above, the asymptotically optimal shift is the solution to the optimization problem (20) in any dimensional setting. The latter can be recast as

$$\omega^{*,p} := \arg \min_{\omega \in \mathbb{C}} \| z_1 F(\cdot) + \omega z_2 \|_{L^p(\Omega)} \quad (21)$$
for two given complex numbers $z_1, z_2$ and a domain $\Omega \subset \mathbb{R}^d$. This optimization problem can actually be solved in closed form for $p = \{2, +\infty\}$ as one can see in the next result.

**Proposition 4.1.** Assume that $F: x \in \Omega \mapsto F(x) \in \mathbb{R}$ is continuous on $\bar{\Omega}$. Then, for any $p \in [1, +\infty]$, the solution to the optimization problem (21) is unique. In addition, for $p \in [1, +\infty)$, it satisfies

$$\int_{\Omega} |z_1 F(x) + \omega^{*,p} z_2|^{p-2} \overline{(z_2 (z_1 F(x) + \omega^{*,p} z_2))} dx = 0.$$

For $p = 2$, this translates to a linear equation whose solution is

$$\omega^{*,2} = - \frac{(F, 1)_{L^2(\Omega)}}{|\Omega|} \times \frac{z_1}{z_2},$$

and the minimal value is

$$\min_{\omega} \|z_1 F(\cdot) + \omega z_2\|_{L^2(\Omega)}^2 = |z_1|^2 \left( \|F\|_{L^2(\Omega)}^2 - \left( \frac{1}{|\Omega|} (F, 1)_{L^2(\Omega)} \right)^2 \right).$$

For $p = \infty$, the optimal solution to (21) is

$$\omega^{*,\infty} = - \left( \frac{F_{\max} + F_{\min}}{2} \right) \frac{z_1}{z_2},$$

and

$$\min_{\omega} \|z_1 F(\cdot) + \omega z_2\|_{L^{\infty}} = |z_1|^2 \frac{|F_{\max} - F_{\min}|}{2},$$

where $F_{\max} = \max_{\Omega} F$ and $F_{\min} = \min_{\Omega} F$.

*Proof*. First, note that for any $p \in [1, +\infty]$ the function $G_p : \omega \in \mathbb{C} \mapsto \|z_1 F(\cdot) + \omega z_2\|_{L^p(\Omega)} \in \mathbb{R}$ is convex and thus has a unique minimum.

For $p \in [1, +\infty)$, we have $\|z_1 F(\cdot) + \omega z_2\|_{L^p(\Omega)} \neq 0$ for any $\omega \in \mathbb{C}$ and thus the function $G_p$ is smooth. As a result, $\omega^{*,p}$ is a critical point of $G_p$. Setting $\omega = \omega_R + i\omega_I$, we can compute the gradient of $G_p$ with respect to the real and imaginary parts of $\omega$. For simplicity, we rather compute the derivatives of $G_p$ which are

$$\partial_{\omega_R} G_p^p(\omega_R, \omega_I) = p \int_{\Omega} |z_1 F(x) + \omega z_2|^{p-2} \Re(\overline{z_2} (z_1 F(x) + \omega z_2)) dx, (22)$$

$$\partial_{\omega_R} G_p^p(\omega_R, \omega_I) = -p \int_{\Omega} |z_1 F(x) + \omega z_2|^{p-2} \Im(\overline{z_2} (z_1 F(x) + \omega z_2)) dx.$$
The critical point of $G_p$ verifies

$$ 
\begin{aligned}
0 &= \partial_{\omega_R} G_p(\omega_R^{*,p}, \omega_I^{*,p}) + i \partial_{\omega_I} G_p(\omega_R^{*,p}, \omega_I^{*,p}) \\
  &= \frac{1}{p G_p^p(\omega_R^{*,p}, \omega_I^{*,p})} \left( \partial_{\omega_R} G_p^p(\omega_R^{*,p}, \omega_I^{*,p}) + i \partial_{\omega_I} G_p^p(\omega_R^{*,p}, \omega_I^{*,p}) \right),
\end{aligned}
$$

from which we obtain the non-linear equation satisfied by $\omega^{*,p}$.

For $p = 2$, the non-linear system (22) reduces to a linear system from which we obtain that $\omega^{*,2}$ satisfies

$$ 
\begin{aligned}
0 &= \int_{\Omega} \Re(\bar{z}_2 (z_1 F(x) + \omega^{*,2} z_2)) - i \Im(\bar{z}_2 (z_1 F(x) + \omega^{*,2} z_2)) dx \\
  &= \int_{\Omega} \overline{(\bar{z}_2 (z_1 F(x) + \omega^{*,2} z_2))} dx,
\end{aligned}
$$

and thus

$$ \omega^{*,2} = - \frac{z_1 (F, 1)_{L^2(\Omega)}}{z_2 |\Omega|}. $$

For $p = \infty$, we start by noting that

$$ \min_{\omega} \max_{x \in \Omega} |F(x)z_1 + \omega z_2|^2 = |z_1|^2 \min_{\omega} \left( \Im(\omega z)^2 + \max_{x \in \Omega} \{F(x) + \Re(\omega z)\}^2 \right), $$

where $z = z_2/z_1$. Introducing the real variables $X = \Re(\omega z)$ and $Y = \Im(\omega z)$, we now have to solve the optimization problem

$$ \min_{X,Y} \left( Y^2 + \max_{x \in \Omega} \{F(x) + X\}^2 \right) = \min_{Y} Y^2 + \min_{X} \left( \max_{x \in \Omega} \{F(x) + X\}^2 \right), $$

from which we readily infer that $Y^* = \Im(\omega^{*,\infty} z) = 0$. To get the optimal $X$, note that $F : x \in \bar{\Omega} \mapsto F(x) \in \mathbb{R}$ is continuous over a compact set and thus reaches its minimal value $F(x_{\min}) = F_{\min} := \min_{x \in \Omega} F(x)$, and its maximal value $F(x_{\max}) = F_{\max} := \max_{x \in \Omega} F(x)$. As a result, for any $X \in \mathbb{R}$, the function $x \to (F(x) + X)^2$ reaches its maximal value either at $x = x_{\min}$ or at $x = x_{\max}$, and thus

$$ \max_{x \in \Omega} \{F(x) + X\}^2 = \max \{ (F_{\max} + X)^2, (F_{\min} + X)^2 \}, $$

which reaches its minimum value at $X^*$ such that $F_{\max} + X^* = -(F_{\min} + X^*)$ hence $X^* = -(F_{\max} + F_{\min})/2$. Using then that $X^* = \Re(\omega^{*,\infty}) =$
$\Re(\omega^{*,\infty})\Re(z) - \Im(\omega^{*,\infty})\Im(z)$ and $Y^* = \Re(\omega^{*,\infty})\Im(z) + \Im(\omega^{*,\infty})\Re(z)$ and doing some computations, we obtain the optimal solution to (21) when using the infinity-norm. $\square$

In what follows, we are going to use Proposition 4.1 to compute the asymptotically optimal shift for the 2d and 3d Yee schemes.

## 4.2. Dispersion correction for dispersive media: Standard vs Frozen approaches

We show here how to perform the asymptotic dispersion correction for media having frequency-dependent parameters. In this case, the scalar wavenumber also depends on $\omega$ and we can introduce the shifted angular frequency in the following two ways:

1. Standard: We replace $\omega$ by $\hat{\omega} = \omega + h^2\omega_2$ everywhere. Using the standard method, the 1d Yee scheme (8) with shifted angular frequency then reads as

$$ \left\{ \begin{array}{l} i\hat{\omega}\varepsilon(\hat{\omega}) E_j - \frac{1}{h}(H_{j+1/2} - H_{j-1/2}) = -J(x_j), \\ i\hat{\omega}\mu(\hat{\omega}) H_{j+1/2} - \frac{1}{h}(E_{j+1} - E_j) = 0, \end{array} \right. $$

and the extension to higher dimensions is straightforward.

2. Frozen: We replace $\omega$ by $\hat{\omega} = \omega + h^2\omega_2$ only in the term $i\omega$ but not in the dielectric parameters $\varepsilon(\omega)$ and $\mu(\omega)$. Using the frozen method, the 1d Yee scheme (8) is now

$$ \left\{ \begin{array}{l} i\hat{\omega}\varepsilon(\omega) E_j - \frac{1}{h}(H_{j+1/2} - H_{j-1/2}) = -J(x_j), \\ i\hat{\omega}\mu(\omega) H_{j+1/2} - \frac{1}{h}(E_{j+1} - E_j) = 0, \end{array} \right. $$

whose extension to higher dimensions is also straightforward.

In the standard case, the asymptotically optimal shift satisfies (20), with $k_s(\omega) = \omega\sqrt{\varepsilon(\omega)\mu(\omega)}$, whose solution can be computed from Proposition 4.1. For the frozen case, we set

$$ k_s^f(\hat{\omega}, \omega) = \hat{\omega}\sqrt{\varepsilon(\omega)\mu(\omega)}, $$

which satisfies $k_s^f(\omega, \omega) = k_s(\omega)$. Note that the expansion (19) with which the discrete wavenumber complies can be written as

$$ k_d(\omega; \boldsymbol{\theta}, h) = k_s^f(\omega, \omega) + (k_s^f(\omega, \omega)h)^2 k_s^f(\omega, \omega)F(\boldsymbol{\theta}) + \dots $$
The frozen dispersion correction now amounts to insert the shifted frequency only in the first argument of $k_s^f$. The discrete wavenumber associated to the frozen dispersion correction then has the next expansion

$$ 
\begin{aligned}
\hat{k}_d^f(\omega; \boldsymbol{\theta}, h) &= k_s^f(\hat{\omega}, \omega) + (k_s^f(\hat{\omega}, \omega)h)^2 k_s^f(\hat{\omega}, \omega)F(\boldsymbol{\theta}) + \dots \\
&= k_s^f(\omega + h^2\omega_2, \omega) + (k_s^f(\omega + h^2\omega_2, \omega)h)^2 k_s^f(\omega + h^2\omega_2, \omega)F(\boldsymbol{\theta}) + \dots \\
&= k_s(\omega) + (k_s(\omega)h)^2 \left( k_s(\omega)F(\boldsymbol{\theta}) + \omega_2 \frac{\partial_{\hat{\omega}}k_s(\hat{\omega}, \omega)}{k_s(\omega)^2} \right) + \dots,
\end{aligned}
$$

where $\partial_{\hat{\omega}}k_s(\hat{\omega}, \omega) = \sqrt{\varepsilon(\omega)\mu(\omega)}$. The asymptotically optimal shift in the frozen case can then be defined as

$$ \omega_2^{*,f,p} := \arg \min_{\omega_2 \in \mathbb{C}} \left\| k_s(\omega) F(\cdot) + \omega_2 \frac{\sqrt{\varepsilon(\omega)\mu(\omega)}}{k_s(\omega)^2} \right\|_{L^p(\Omega)} . \quad (23) $$

Note that both standard (hence given by (20)) and frozen (hence given by (23)) asymptotically optimal shifts can be computed from Proposition 4.1 both with

$$ z_1 = k_s(\omega), \quad F = F(\boldsymbol{\theta}). $$

The value of $z_2$ is however different, since

$$ z_2^{\text{std}} = \frac{\partial_{\omega} k_s(\omega)}{k_s(\omega)^2}, $$

for the standard case and

$$ z_2^f = \frac{\sqrt{\varepsilon(\omega)\mu(\omega)}}{k_s(\omega)^2}, $$

for the frozen case and the shifts are then going to be different regarding the approach chosen. Nevertheless, note that the values of $\min_{\omega} \|z_1F(\cdot) + \omega z_2\|_{L^p(\Omega)} = \|z_1F(\cdot) + \omega^{*,p}z_2\|_{L^p(\Omega)}$ for $p = 2, +\infty$ do not depend on $z_2$ (see Proposition 4.1). As a result, the reduction of the dispersion error brought by the asymptotically optimal shift is going to be the same for both approaches. In what follows we will only rely on the frozen approach to give explicit formulas for the asymptotically optimal shift because this does not require to compute the derivative of the dielectric parameters with respect to the frequency.
### 4.3. No dispersion error for the 1d-Yee scheme

The dispersion relation associated to the 1d-Yee scheme is given by

$$ \omega^2 \varepsilon \mu - \frac{4}{h^2} \sin \left( \frac{kh}{2} \right)^2 = 0, $$

from which we get that the discrete wavenumber is

$$ k_d(\omega; h) = \frac{2}{h} \arcsin \left( \frac{\omega \sqrt{\varepsilon \mu h}}{2} \right), $$

where $\arcsin$ is the principal branch of the inverse sine function defined for $\omega\sqrt{\varepsilon\mu h}/2 \in \mathbb{C} \setminus ((-\infty, -1) \cup (1, +\infty))$. In this one-dimensional setting, we can find a shifted angular frequency $\hat{\omega}$ for which there is no dispersion error. This amounts to find $\hat{\omega}$ such that

$$ k_d(\hat{\omega}; h) = k_s(\omega). $$

Using that $\sin(\arcsin(z)) = z$ for any $z \in \mathbb{C}$ with $|z| \le 1$, we obtain that

$$ \hat{\omega} = \frac{2}{h\sqrt{\varepsilon\mu}} \sin \left( \frac{\omega\sqrt{\varepsilon\mu h}}{2} \right), \qquad (24) $$

where the derivation of this formula is valid as soon as $G \ge \pi$. A Taylor expansion gives

$$ \hat{\omega} = \omega - h^2 \frac{\omega^3 \varepsilon \mu}{24} + \dots = \omega - (k_s h)^2 \frac{\omega}{24} + \dots, $$

from which we see that the shifted angular frequency is a second order perturbation (in $k_s h$) of $\omega$. This ensures that the Yee scheme with shifted angular frequency, replacing $\omega$ by $\hat{\omega}$, is still second order accurate.

### 4.4. Asymptotic dispersion correction for the 2d Yee scheme

We compute next the asymptotically optimal shift $\omega_2^{*,p}$ for the 2d Yee scheme whose expansion of its discrete wavenumber can be found in (14). We then have to solve the optimization problem (20) with

$$ F(\boldsymbol{\theta}) = \frac{1}{24} (\cos(s)^4 + \sin(s)^4), $$

which fits in the framework of Proposition 4.1 with

$$ z_1 = k_s(\omega) = \omega\sqrt{\varepsilon\mu} \text{ and } z_2 = \frac{1}{\omega^2\sqrt{\varepsilon\mu}}. $$

We then have the following result.
**Proposition 4.2.** The asymptotically optimal shifts for the 2d-Yee's scheme are

$$ \omega_2^{*,\infty} = \omega_2^{*,2} = -\frac{\omega^3 \varepsilon \mu}{32}. $$

*Proof*. We apply Proposition 4.1 with the above $z_1$ and $z_2$. For $\omega_2^{*,2}$, we only need to compute the mean value of $F$ over $[0, 2\pi]$,

$$ \frac{1}{2\pi} \int_{0}^{2\pi} \frac{1}{24} (\cos(s)^4 + \sin(s)^4) \, ds = \frac{1}{32}. $$

This gives

$$ \omega_2^{*,2} = -\frac{z_1}{z_2} \frac{1}{2\pi} \int_{0}^{2\pi} F(\theta) \, ds = -\frac{\omega^3 \varepsilon \mu}{32}. $$

To get $\omega_2^{*,\infty}$, we have to find the extremal values of $F$ over $[0, 2\pi]$. Setting $t = \sin(s)^2$ and using that $\cos(s)^2 = 1 - \sin(s)^2$, we get

$$ F(\theta) = Y(t) = \frac{1}{24}((1-t)^2 + t^2), $$

and we now have to find the extremal values of $Y$ over $[0, 1]$. Direct computations give that $F_{\min} = 1/48$ and $F_{\max} = 1/24$ and thus

$$ \omega_2^{*,\infty} = -\left(\frac{F_{\min} + F_{\max}}{2}\right) \frac{z_1}{z_2} = -\frac{\omega^3 \varepsilon \mu}{32}. $$

<figure><img src="image_2.png" /></figure>

In this case, both shifts are equal and are the same as the one obtained for the 5-point stencil in [14]. In [14, Section 4.1], we compared the relative dispersion errors for the schemes with and without dispersion correction and showed how they are reduced with the shift (see [14, Figure 4.1]). Since these results are going to hold without any changes for the 2d Yee scheme, we do not reproduce them here.

### 4.5. Asymptotic dispersion correction for the 3d Yee scheme

We compute now the asymptotically optimal shift $\omega_2^{*,p}$ for the 3d Yee scheme whose discrete wavenumber has been expanded for $k_s h$ small in (17). We then have to solve the optimization problem (20) with

$$ F(\theta) = \frac{1}{24} (\cos(\phi)^4 \sin(\theta)^4 + \sin(\phi)^4 \sin(\theta)^4 + \cos(\theta)^4), $$
which fits the framework of Proposition 4.1 with

$$z_1 = k_s(\omega) = \omega\sqrt{\varepsilon\mu} \text{ and } z_2 = \frac{1}{\omega^2\sqrt{\varepsilon\mu}}.$$

**Proposition 4.3.** We have $\omega_2^{*,\infty} = -\frac{\omega^3\varepsilon\mu}{36}$, $\omega_2^{*,2} = -\omega^3\varepsilon\mu\frac{7}{256}$.

*Proof*. From Proposition 4.1, we have

$$\omega_2^{*,2} = -\frac{z_1 \int_0^{2\pi} \int_0^{\pi} F(\phi, \theta) \, d\phi \, d\theta}{z_2 2\pi^2}.$$

Computing the above integral (which gives $21\pi^2/(16 \times 24)$), we obtain $\omega_2^{*,2} = -\omega^3\varepsilon\mu\frac{7}{256}$.

For the asymptotically optimal shifts with the $L^\infty$-norm, we need the extremal values of $F$ over $\Omega = [0, 2\pi] \times [0, \pi]$. Fortunately, we already computed the extremal values of this function in [14, Theorem 4.3] when doing dispersion correction for the 7-point stencil applied to the 3d-Helmholtz equation. They are given by $F_{\max} = \frac{1}{24}$, $F_{\min} = \frac{1}{72}$ and we finally obtain $\omega_2^{*,\infty} = -\frac{\omega^3\varepsilon\mu}{36}$.

<figure><img src="image_3.png" /></figure>

We emphasize that

$$\frac{|\omega_2^{*,\infty} - \omega_2^{*,2}|}{|\omega^3\varepsilon\mu|} = \left| \frac{1}{36} - \frac{7}{256} \right| \approx 0.0004340277778, \quad (25)$$

and thus using one or the other should yields similar results. We can also compute the (asymptotic) reduction factors which indicate how much the dispersion error is reduced when using the shifted angular frequency. The asymptotic reduction factors are the ratio of dispersion errors with and without shifted angular frequency,

$$R_f^p = \frac{\|k_d(\omega; \cdot, h) - k_s\|_{L^p}}{\|k_d(\hat{\omega}^{*,p}; \cdot, h) - k_s\|_{L^p}},$$

where $\hat{\omega}^{*,p} = \omega + h^2\omega_2^{*,p}$ is the shifted angular frequency with asymptotically optimal shift. Using the expansion of the discrete wavenumber as $k_s h \to 0$ together with Proposition 4.1, we get

$$\begin{align}
R_f^\infty &= \frac{\max\{F_{\max}, F_{\min}\}}{\frac{|F_{\max} - F_{\min}|}{2}} + \dots = 3 + \dots, \nonumber \\
R_f^2 &= \frac{\|F\|_{L^2}}{\|F - \frac{7 \times 24}{256}\|_{L^2}} + \dots = \sqrt{\frac{643}{55}} + \dots. \tag{26}
\end{align}$$
<figure><img src="image_4.png" /><figcaption>Figure 1: $k_d(\omega; \theta, h)/k_s(\omega)$ for $\varepsilon = \mu = 1$, $\omega = 100$ and $h = 0.01$. Left: No dispersion correction. Right: with dispersion correction.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 2: Log-log plot of $\max_{\theta}|k_d(\omega; \theta, h) - k_s(\omega)|/k_s(\omega)$ for $\varepsilon = \mu = 1$, $h = 0.01$ and various $\omega$.</figcaption></figure>

Since $\sqrt{643}/55 \approx 3.42$, the asymptotic dispersion correction in the $L^2$-norm reduces the dispersion error a bit more than the asymptotic dispersion in the $L^\infty$-norm.

We show in Figure 1 the ratio of the discrete wavenumber over the continuous one with and without dispersion correction using $\hat{\omega}^{*,2} = \omega + h^2\omega_2^{*,2}$ for the shifted frequency. In this setting, we have roughly $G = \frac{2\pi}{\omega\sqrt{\varepsilon\mu h}} \approx 6.28$ points per wavelength. We see that, without dispersion correction, the discrete wavenumber is always above the continuous one while $k_d$ oscillates around $k_s$ when using dispersion correction. We show in Figure 2 the dispersion error, with and without dispersion correction, for various values of $G$. These numerical results show that the dispersion error is reduced when using the shift and that it behaves like $O(G^{-2})$ as expected from (18).
<figure><img src="image_6.png" /><figcaption>Figure 3: Top row: Discrete wavenumber $k_d(\omega; \theta, h)$. Bottom row: discrete wavenumber $k_d(\hat{\omega}^{\text{asy},2}; \theta, h)$ using dispersion correction. Left column: real part. Right column: imaginary part.</figcaption></figure>

We now illustrate the effect of the $L^2$-dispersion correction (hence using $\hat{\omega} = \omega + h^2\hat{\omega}^{*,2}$) for the complex dielectric parameters

$$ \varepsilon = 50 - 12i, \quad \mu = 1, \qquad (27) $$

which, from [1, Table 1], aim at modeling the dielectric properties of the human body. Note that the continuous wavenumber is given here as $k_s = \omega\sqrt{\varepsilon\mu} \approx \omega(7.1211 - 0.8426i)$. We show in Figure 3 the real and imaginary parts of the discrete wavenumber for $\omega = 10$ and $h = 0.01$ which yields $G = \frac{2\pi}{|k_s|h} \approx 8.76$, grid points per wavelength.

Note that both real and imaginary parts of the discrete wavenumber are above the continuous one without dispersion correction but both oscillate around $k_s$ when using the shifted frequency. To further test the effect of the asymptotic dispersion correction, we now use $\omega = 5, 10, 20, 40$, so that the number of grid points per wavelength is approximately $G = 17.5245, 8.7622, 4.3811, 2.1906$.
<figure><img src="image_7.png" /><figcaption>Figure 4: Log-log plot of $\max_{\theta} |k_d(\omega; \theta, h) - k_s(\omega)| / |k_s(\omega)|$ for $\varepsilon = 50-12i$, $\mu = 1$, $h = 0.01$ and various $\omega$.</figcaption></figure>

We show in Figure 4 the dispersion error for varying frequencies which results in varying the number of grid point per wavelength. These numerical results shows that the dispersion error behaves like $O(G^{-2})$ for large enough $G$ and that the dispersion correction indeed reduces the dispersion error. Note that the dispersion correction fails if too few grid points per wavelength are used. This is expected since the dispersion correction has been done for small enough $h$ hence for large enough $G$.

We end this section by computing numerically the optimal shift, denoted by $\omega_2^{\text{opt}}$, which minimizes the $L^2$-relative dispersion error

$$ \omega_2^{\text{opt}} = \arg\min_{\omega_2} \left\| \frac{k_d(\omega + h^2\omega_2; \theta, h) - k_s(\omega)}{k_s(\omega)} \right\|_{L^2}. $$

These numerical experiments check how far the (explicit) asymptotically optimal shift is from the optimal one. The results are given in Table 1, using the dielectric parameters from Eq. (27) and various frequencies $\omega$. In Table 1, we compute the relative dispersion error

$$ \text{Err}_d(\omega_2) = \left\| \frac{k_d(\omega + h^2\omega_2; \theta, h) - k_s(\omega)}{k_s(\omega)} \right\|_{L^2}, $$

when no dispersion correction is performed, hence using $\omega_2 = 0$, when using the optimal shift and when using the asymptotically optimal shift. These results show that the asymptotically optimal shift is actually close to the optimal one even if a small number of grid points per wavelength is considered.
<table><caption>Table 1: Comparisons between the asymptotically optimal shift and the optimal one for the frequencies $\omega = 5, 10, 20, 40$ with the dielectric parameters from (27) and meshsize $h = 0.01$.</caption><tr><td>G</td><td>17.5245</td><td>8.7622</td><td>4.3811</td><td>2.1906</td></tr><tr><td>Err<sub>d</sub>(0)</td><td>0.0038</td><td>0.0155</td><td>0.0730</td><td>0.4142</td></tr><tr><td>Err<sub>d</sub>(ω<sub>2</sub><sup>opt</sup>)</td><td>0.0011</td><td>0.0045</td><td>0.0204</td><td>0.1295</td></tr><tr><td>Err<sub>d</sub>(ω<sub>2</sub><sup>asy</sup>)</td><td>0.0011</td><td>0.0045</td><td>0.0205</td><td>0.1328</td></tr><tr><td>|ω<sub>2</sub><sup>opt</sup> - ω<sub>2</sub><sup>asy</sup>| / |ω<sub>2</sub><sup>opt</sup>|</td><td>0.0082</td><td>0.0135</td><td>0.0252</td><td>0.0587</td></tr></table>

### 4.6. Which $L^p$-norm should be chosen to minimize the asymptotic dispersion error?

As shown in Proposition 4.2 and 4.3, the $L^2$ and $L^\infty$ shifts are either equal (in 2d) or very close to each other (in 3d). The only difference appears in the value of the asymptotic reduction factors $R_f^p$ but the two norms being different, such a discrepancy is to be expected. In this subsection, we further study the influence of $p$ on the asymptotically optimal shift.

Assuming that $\varepsilon, \mu$ are real-valued, we get that $z_1, z_2 \in \mathbb{R}$ and we can then transform the minimization problem defining $\omega_2^{*,p}$, namely

$$ \omega_2^{*,p} = \arg \min_{\omega_2} \|z_1 F(\cdot) + z_2 \omega_2\|_{L^p(\Omega)}^p $$

into the minimization problem

$$ \tilde{\omega}_2^{*,p} = \arg \min_{\omega_2} \|F(\cdot) + \tilde{\omega}_2\|_{L^p(\Omega)}^p. \qquad (28) $$

The two are then related by $\tilde{\omega}_2^{*,p} = \omega_2^{*,p} z_2 / z_1$. This change of variables allows us to eliminate the physical parameters and to work only with $\omega = 1$ and $\varepsilon = \mu = 1$.

In (28), the function $F$ is defined in 2d by $F(s) = (\cos(s)^4 + \sin(s)^4)/24$ with $\Omega = [0, 2\pi]$ and in 3d by $F(\phi, \theta) = (\cos(\phi)^4 \sin(\theta)^4 + \sin(\phi)^4 \sin(\theta)^4 + \cos(\theta)^4)/24$ with $\Omega = [0, 2\pi] \times [0, \pi]$. Proposition 4.1 then ensures that $\tilde{\omega}_2^{*,p}$ is the unique solution to the non-linear equation

$$ C(\tilde{\omega}_2^{*,p}) = \int_{\Omega} (F(x) + \tilde{\omega}_2^{*,p}) |F(x) + \tilde{\omega}_2^{*,p}|^{p-2} dx = 0. \qquad (29) $$

We compute numerically the solution to the minimization problem (28) using the Nelder-Mead algorithm. We approximate the $L^p$-norms with Riemann
<figure><img src="image_8.png" /><figcaption>Figure 5: Top left: $\omega_{\text{approx}}^p$, the asymptotically optimal shifts computed with the Nelder-Mead algorithm; Top right: $C(\omega_{\text{approx}}^p)$; Bottom left: $L^p$ reduction factors defined by $R_f^p = \|F\|_{L^p(\Omega)} / \|F(\cdot) + \tilde{\omega}_{\text{approx}}^p\|_{L^p(\Omega)}$; Bottom right: $|\tilde{\omega}_{\text{approx}}^p - \tilde{\omega}_2^{*,\infty}|$.</figcaption></figure>

sums using a regular mesh of $\Omega$ composed of rectangles with $10^5$ points in 2d and $10^3 \times 10^3$ points in 3d. We then verify the numerical solution $\tilde{\omega}_{\text{approx}}^p$ is close to $\tilde{\omega}_2^{*,p}$ by computing $C(\tilde{\omega}_{\text{approx}}^p)$.

The results are shown in Figure 5, where we used $p = 2^j$ for $j = 0, 1, \dots, 8$. One can see from the values of $C(\tilde{\omega}_{\text{approx}}^p)$ that the optimality conditions are well satisfied by the $\tilde{\omega}_{\text{approx}}^p$. Also, in the 2d and 3d settings, both reduction factors reach their maximal value for $p = 1$. In addition, we see that the asymptotically optimal shifts in the 2d setting do not depend on the $L^p$ norm chosen to minimize the asymptotic dispersion error, whereas they slightly vary in 3d. We can therefore conclude that, in 2d, the choice of $L^p$-norm used for the dispersion correction does not matter.

For the 3d case, for every considered $p$, the values of $\tilde{\omega}_{\text{approx}}^p$ are actually close (up to 3 digits except for $p = 1$ where only 2 digits match) to $\tilde{\omega}_2^{*,\infty}$. To
<table><caption>Table 2: Reduction factor using either the numerically optimized shift or $\tilde{\omega}_2^{*,\infty} = -1/36$.</caption><thead><tr><th>p</th><th>1</th><th>2</th><th>4</th><th>8</th><th>16</th><th>32</th><th>64</th><th>128</th></tr></thead><tbody><tr><td>R<sub>f</sub><sup>p</sup> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">||F||<sub>L<sup>p</sup>(&Omega;)</sub></div><div>||F(&middot;)+&omega;<sub>approx</sub><sup>p</sup>||<sub>L<sup>p</sup>(&Omega;)</sub></div></div></td><td>3.7742</td><td>3.4183</td><td>3.2129</td><td>3.1259</td><td>3.0827</td><td>3.0535</td><td>3.0333</td><td>3.0202</td></tr><tr><td>R&#771;<sub>f</sub><sup>p</sup> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">||F||<sub>L<sup>p</sup>(&Omega;)</sub></div><div>||F(&middot;)-1/36||<sub>L<sup>p</sup>(&Omega;)</sub></div></div></td><td>3.7097</td><td>3.4141</td><td>3.2079</td><td>3.1073</td><td>3.0606</td><td>3.0343</td><td>3.0189</td><td>3.0101</td></tr></tbody></table>

justify choosing $p = \infty$ in 3d, we compare in Table 2 the reduction factors when using $\tilde{\omega}_{approx}^p$ with the reduction factors when using $\tilde{\omega}_2^{*,\infty} = -1/36$. This shows that using $\omega_2^{*,\infty}$ instead of the asymptotically optimal one for $p \neq \infty$ still reduces the asymptotic dispersion error with a reduction factor that is very close to the optimal one. Hence, the choice of $p$ does not matter and using $p = +\infty$ is sufficient for an efficient dispersion correction.

## 5. Numerical experiments

This section is devoted to numerical results to show the effect of the dispersion correction (using the shifted angular frequency in the Yee scheme) on the relative error for the one and the three dimensional setting<fnref n="2" />. For each case, we will consider manufactured solutions that oscillate proportionally to $1/k_s$. In addition, we work with homogeneous media first, a material with piecewise constant parameters next, and finally with a media having slowly varying parameters. Finally, for the 3d results, we only consider the $L^\infty$-asymptotically optimized shift since it is very close to the $L^2$-one (see (25)).

### 5.1. Numerical results for the 1d Maxwell's equation

This section gives some numerical experiments to show the effect of the shifted angular frequency on the performance of the 1d Yee scheme. We consider first a homogeneous material, next a media with piecewise constant parameters, and then a media with slowly varying parameters.

#### 5.1.1. Homogeneous media

We consider here the 1d-Maxwell's equation on $\Omega = [0, 1]$ with homogeneous Dirichlet boundary conditions $E|_{\partial\Omega} = 0$. A manufactured solution to

<footnote n="2">²Since the dispersion properties of the 2d Yee scheme are the same as those of the 5-point stencil applied to the Helmholtz equation, we do not provide numerical experiments for this case that has already been investigated in [14].</footnote>
(2) can be found as

$$E^{ex}(x) = x(1-x) \cos(\omega\sqrt{\varepsilon\mu}x), \quad H^{ex}(x) = \cos(\omega\sqrt{\varepsilon\mu}x) + 3,$$

with source terms<fnref n="3" />

$$J(x) = \partial_z H^{ex} - i\omega\varepsilon E^{ex}, \quad M(x) = \partial_z E^{ex} - i\omega\mu H^{ex}.$$

In these numerical experiments, we set $\varepsilon = 1$, $\mu = 1$, together with the angular frequencies $\omega = 25, 50, 100, 200$. We consider a uniform mesh with $N+1$ interior points so that the meshsize is $h = 1/N$. Finally, to highlight the pollution effect as $\omega$ increases, we work with the numbers of interior points

$$N = 2^j \times 10, \quad j = 0, 1, \dots, 7.$$

As a result, the number of grid points per wavelength, namely $G = 2\pi/(\omega\sqrt{\varepsilon\mu}h)$, decreases as $\omega$ increases although the mesh sizes are fixed for each $\omega$ considered.

The numerical results are in Figure 6. From these numerical experiments, we clearly see the pollution effect together with its suppression when using the shifted angular frequency. In this 1d-setting, the shifted angular frequency completely suppresses the dispersion error and we can thus conclude that the pollution effect is driven by the dispersion error.

### 5.1.2. Media with piecewise constant parameters

We next show numerical results for media with piecewise constant parameters. To do this, we consider $\Omega = [-1, +1]$ and the parameters

$$\varepsilon(x) = \varepsilon_1 \mathbf{1}_{x<0}(x) + \varepsilon_2 \mathbf{1}_{x>0}(x), \quad \mu(x) = \mu_1 \mathbf{1}_{x<0}(x) + \mu_2 \mathbf{1}_{x>0}(x),$$

together with the manufactured solution to (A.1)

$$E^{ex}(x) = (1 - x^2) \cos(Kx), \quad H^{ex}(x) = \cos(Kx) + 3,$$

with $K = \max\{|\omega\sqrt{|\varepsilon_i\mu_i|}\} = \max\{k_i\}$. The source terms are again given by

$$J(x) = \partial_z H^{ex} - i\omega\varepsilon E^{ex}, \quad M(x) = \partial_z E^{ex} - i\omega\mu H^{ex},$$

<footnote n="3">³We consider here Maxwell's equation with magnetic current $M$ that are involved in Maxwell's equations as $i\omega\mu H + \nabla \times E = -M$.</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 6: Relative error with and without shifted angular frequency for the 1d Maxwell's equation. Top left: $\omega = 25$; Top right: $\omega = 50$; Bottom left: $\omega = 100$; Bottom right: $\omega = 200$.</figcaption></figure>
and the permittivity and permeability will be

$$ \varepsilon_1 = \mu_1 = 1, \varepsilon_2 = \mu_2 = \alpha, $$

for varying values of the contrast $\alpha$. Again, we consider a uniform mesh with $N+1$ interior points so the meshsize is $h = 2/N$. Since we work on $[-1, +1]$, we now chose

$$ N = 2 \times 2^j \times 10, j = 0, 1, \dots, 7, $$

so that the meshsizes are the same as those from the previous numerical experiments.

The Yee scheme (8) can be applied for piecewise-constant parameters as well with the only exception at the grid point $x = 0$ for which the parameters are not defined. Since the permeability is computed at intermediate grid points $jh + h/2$, there is no problem with the discontinuity but the permittivity needs to be computed at $x = 0$ for which we simply use its mean value $\varepsilon(0) = (\varepsilon_1 + \varepsilon_2)/2$. The shifted angular frequency is then applied directly on each sub-interval and, for the grid point corresponding to $x = 0$ we also use the mean value of the left and right shifted angular frequencies.

The numerical results can be found in Figure 6 where the number of grid point per wavelength is computed as $G = 2\pi/(Kh)$. We see that the dispersion correction once again greatly reduces the pollution effect. Indeed, when using dispersion correction, we recover the expected second order accuracy for a number of grid points per wavelength far smaller than the one required for the un-corrected scheme to converge. In addition, for large enough $G$, the relative error when using the shifted angular frequency is always smaller than the relative error of the standard Yee's scheme.

We end this subsection with additional numerical experiments using now 4 sub-intervals and higher contrasts. We consider the parameters

$$ \varepsilon(x) = \varepsilon_1 \mathbf{1}_{x<-0.5}(x) + \varepsilon_2 \mathbf{1}_{-0.5<x<0}(x) + \varepsilon_3 \mathbf{1}_{0<x<0.5}(x) + \varepsilon_4 \mathbf{1}_{0.5<x}(x), $$

$$ \mu(x) = \mu_1 \mathbf{1}_{x<-0.5}(x) + \mu_2 \mathbf{1}_{-0.5<x<0}(x) + \mu_3 \mathbf{1}_{0<x<0.5}(x) + \mu_4 \mathbf{1}_{0.5<x}(x), $$

with values

$$ \varepsilon_1 = \mu_1 = \varepsilon_4 = \mu_4 = 1, \varepsilon_2 = \mu_2 = \frac{3}{2}\alpha, \varepsilon_3 = \mu_3 = 2\alpha. $$

At the discontinuities located at $x \in \{-0.5, 0, 0.5\}$ we use the mean value of the permittivity. As done above, the shifted angular frequency is again
<figure><img src="image_10.png" /><figcaption>Figure 7: Relative error with and without shifted angular frequency for the 1d Maxwell's equation. Top row: Contrast $\alpha = 2$. Bottom row: Contrast $\alpha = 4$. Left column: $\omega = 25$. Right column: $\omega = 50$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 8: Relative error with and without shifted angular frequency for the 1d Maxwell's equation using 4 sub-intervals and $\alpha = 4$ hence with the largest contrast being 8. Left: $\omega = 25$. Right: $\omega = 50$.</figcaption></figure>

applied on each sub-interval and the mean value is used at the discontinuities. Figure 8 shows the numerical results using the same mesh sizes as before for $\omega = 25, 50$. We also set $\alpha = 4$ and thus have a maximal contrast of 8. Similar conclusions can be drawn from these results, namely that the 1d Yee scheme with shifted angular frequency greatly reduces the relative error as soon as a large enough number of grid points is considered.

### 5.1.3. Media with slowly varying smooth parameters

We consider the same numerical setting as in Section 5.1.2, but now the media has slowly varying smooth parameters given by

$$ \varepsilon(x) = \mu(x) = 1 + \delta |\cos(\omega_0 x)|. \qquad (30) $$

The manufactured solution is now defined with $K = \omega \|\sqrt{\varepsilon\mu}\|_{L^{\infty}}$. We emphasize that the shifted angular frequency is defined in (24) and involves the values of $\varepsilon$ and $\mu$. In the case of spatially varying $\varepsilon, \mu$, we simply used formula (24) at each grid point to define the Yee stencil with dispersion correction.

The numerical results are shown in Figure 9, where the number of grid points per wavelength is now computed as $G = 2\pi / (\omega \|\sqrt{\varepsilon\mu}\|_{L^{\infty}} h)$. These numerical results show that the dispersion correction also reduces the relative error for large enough $G$ even for media with spatially varying parameters. More precisely, on the finest grid, the ratios of relative error without dispersion correction are 9.13 for $\delta = 0.1$ and 4.6 for $\delta = 0.5$, indicating that, as $\delta$ increases, the reduction factor decreases.
<figure><img src="image_12.png" /><figcaption>Figure 9: Relative error with and without shifted angular frequency for the 1d Maxwell's equation using the spatially varying parameter (30) with $\omega = 50$ and $\omega_0 = 20\pi$. Left: $\delta = 0.1$. Right: $\delta = 0.5$.</figcaption></figure>

## 5.2. Numerical results for the 3d Maxwell's equations

Next we solve the Maxwell's equations (1) on $\Omega = (0, \pi)^3$ with Perfect Electric Conductor (PEC) boundary conditions on $\partial\Omega$,

$$ \boldsymbol{n} \times \mathbf{E} = 0 \text{ on } \partial\Omega, $$

with $\boldsymbol{n}$ being the outward unit normal.

We use a manufactured solution to Maxwell equations on $\Omega$ with PEC boundary conditions. Considering the electric field

$$ E_1^{ex} = \sin(z) \sin(y) \sin(kx), \quad E_2^{ex} = \sin(z) \sin(x) \sin(ky), \quad E_3^{ex} = \sin(y) \sin(x) \sin(kz), $$

we can compute the associated magnetic field using $H^{ex} = -(i\omega\mu)^{-1}\nabla \times \mathbf{E}^{ex}$ which is thus

$$ \begin{aligned} H_1^{ex} &= \frac{1}{\omega\mu} (-i\sin(x)(\cos(z)\sin(ky) - \cos(y)\sin(kz))), \\ H_2^{ex} &= \frac{1}{\omega\mu} (i\sin(y)(\sin(kx)\cos(z) - \cos(x)\sin(kz))), \\ H_3^{ex} &= \frac{1}{\omega\mu} (i\sin(z)(\cos(x)\sin(ky) - \sin(kx)\cos(y))). \end{aligned} \quad (31) $$

The electric current can then be obtained by computing $\mathbf{J}^{ex} = -i\omega\varepsilon\mathbf{E}^{ex} + \nabla \times \mathbf{H}^{ex}$.
We emphasize that, du to the way it has been computed, the above manufactured solution can also be used for spatially varying permittivity but only for constant permeability. In all numerical experiments, we set

$$k = \omega \sqrt{\|\varepsilon\|_{L^{\infty}(\Omega)} |\mu|},$$

so that our explicit solution oscillates with period $1/k$. In addition, we are going to compute the relative error

$$\mathbf{Err} := \frac{\max \{\|\mathbf{E}_h - \mathbf{E}^{\text{ex}}\|_{l^{\infty}}, \|\mathbf{H}_h - \mathbf{H}^{\text{ex}}\|_{l^{\infty}}\}}{\max \{\|\mathbf{E}^{\text{ex}}\|_{l^{\infty}}, \|\mathbf{H}^{\text{ex}}\|_{l^{\infty}}\}},$$

where the $l^{\infty}$-norm is defined as the maximal value over all grid points.

In what follows, we study the two following settings: we consider first the case of media with homogeneous and real permittivity and permeability. Next, we test the dispersion correction on Maxwell's equations with piece-wise constant permittivity having constant permeability.

In all our numerical experiments, we use a uniform $N \times N \times N$ grid with $N = 8, 16, 32, 64$ point in each direction which results in having 512, 4096, 32768 and 262144 degrees of freedom for each of the six component of the electric and magnetic fields. We work with fixed values for the permeability and the permittivity but with varying angular frequencies ranging from 1 to 10 with $\omega = 1, 2.5, 5, 10$. All our convergence results are presented using the number of grid points per wavelength $G = 2\pi/(kh)$, with $h = \pi/(N-1)$, instead of using the meshsize $h$, since $G$ is the relevant parameter for time harmonic wave propagation problems. We actually have the equivalence between the asymptotic dispersion correction as $h \to 0$ and the one as $G \to +\infty$ (see also [14, Appendix A]). Finally, to see the influence of the asymptotic dispersion correction, we are going to compute the reduction factor defined as $Err(\omega)/Err(\hat{\omega})$, corresponding to the ratio of relative errors without dispersion correction over relative error when using the shifted angular frequency.

### 5.2.1. Experiments with homogeneous real parameters

We consider the homogeneous and real parameters $\varepsilon = 1, \mu = 1$. We show in Figure 10 the relative errors without dispersion correction and using the $L^{\infty}$ dispersion correction (see Proposition 4.3 for the formula). We also give in Table 3 the reduction factors corresponding to the ratio of the relative error without and with dispersion correction.
<figure><img src="image_13.png" /><figcaption>Figure 10: Relative error with and without shifted angular frequency. Top left: $\omega = 1$; Top right: $\omega = 2.5$; Bottom left: $\omega = 5$; Bottom right: $\omega = 10$.</figcaption></figure>

<table><caption>Table 3: Ratio of the errors without to with dispersion correction for varying meshsize and angular frequency for media with constant real parameters using $L^\infty$ dispersion correction.</caption><thead><tr><th>N</th><th>8</th><th>16</th><th>32</th><th>64</th></tr></thead><tbody><tr><td>ω = 1</td><td>1.8093</td><td>1.7972</td><td>1.7935</td><td>1.7922</td></tr><tr><td>ω = 2.5</td><td>0.4146</td><td>1.3256</td><td>1.7444</td><td>1.8696</td></tr><tr><td>ω = 5</td><td>0.5656</td><td>1.6861</td><td>2.0678</td><td>2.1842</td></tr><tr><td>ω = 10</td><td>0.0192</td><td>0.0922</td><td>2.8673</td><td>2.7275</td></tr></tbody></table>
<table><caption>Table 4: Ratio of the errors without dispersion correction for varying meshsize and angular frequency for media with constant real parameters using $L^\infty$ dispersion correction.</caption><thead><tr><th>N</th><th>5</th><th>10</th><th>20</th><th>40</th></tr></thead><tbody><tr><td>ω = 1</td><td>1.8325</td><td>1.8189</td><td>1.8154</td><td>1.8141</td></tr><tr><td>ω = 2.5</td><td>0.3759</td><td>1.2572</td><td>1.6657</td><td>1.7883</td></tr><tr><td>ω = 5</td><td>0.5442</td><td>1.4380</td><td>2.0130</td><td>2.2769</td></tr><tr><td>ω = 10</td><td>0.0126</td><td>0.0930</td><td>3.3351</td><td>5.5244</td></tr></tbody></table>

We see from the numerical results given in Figures 10 that the dispersion correction allows to reduce the relative error as soon as a large enough number of grid point per wavelength is considered. This behavior is expected since the asymptotically optimal shift is defined through the minimization of the dispersion error for small enough meshsize (or equivalently for $G$ large enough). Table 3 gives the reduction factor from which we see that the error is reduced when using the asymptotically optimal shift and $G$ is large enough. Note that the relative error can be divided by 2 using this dispersion correction for no additional cost!

### 5.2.2. Experiments with piecewise constant parameters

We finally consider a material with constant permeability and piecewise constant permittivity defined as

$$ \epsilon = \begin{cases} 1 & \text{if } x \in \Omega \setminus C \\ 0.7 - 0.1i & \text{if } x \in C \end{cases}, \mu = 1, $$

with $C =]1/3, 2/3]^3$. We show the relative errors in Figure 11 and the reduction factors in Table 4.

From these results, we can see that the shift reduces the relative error for large enough number of grid points per wavelength. In addition, the error can be roughly divided by 2 up to 5.

## 6. Conclusions and Outlook

We provided a dispersion correction technique for Yee finite difference scheme applied to the time-harmonic Maxwell's equations. Our method is based on a modification of the angular frequency in the Yee stencil and can thus be easily implemented in numerical codes. We showed, through numerical experiments that this dispersion correction reduces the relative error as long as one has a large enough number of grid points per wavelength.
<figure><img src="image_14.png" /><figcaption>Figure 11: Relative error with and without shifted angular frequency. (Top left) $\omega = 1$; (Top right) $\omega = 2.5$; (Bottom left) $\omega = 5$; (Bottom right) $\omega = 10$.</figcaption></figure>
It would be interesting to get more theoretical insight on the effect of the dispersion correction on the convergence of the scheme. Some ideas might be borrowed from [29] where Fourier analysis is used to derive precise error bounds for the 1d Helmholtz equation. Another interesting future work is related to extending this dispersion correction method to materials modelled by other constitutive relations like e.g. chiral media [12, 55] or bi-anisotropic materials [27].

We would also like to extend our dispersion correction method based on the introduction of a shifted angular frequency to other numerical methods like the fourth-order non-compact stencil from [23], the fourth-order compact stencil from [47, 48] and also to finite element discretizations involving Nédélec edge elements [40, 3, 2].

We finally emphasize that our dispersion correction method is based on a modification of the shifted angular frequency, which thus also modifies the wave speed at the discrete level. As a result, this technique is likely to be difficult to generalize to time-dependent wave propagation problems, and finding a dispersion correction method for such problems will probably require a different approach.

# Appendix A. Eigenvalues analysis for the 1d Maxwell's equation with PEC boundary condition

The (continuous and discrete) Maxwell's equations with Perfect Electric Conductor (PEC) boundary conditions are singular for some values of $\omega$. In this appendix, we show that, when using the shifted angular frequency $\hat{\omega}$ given in (24), both continuous and discrete problems are singular for the same values of $\omega$. Note that the shift has a similar effect for the 1d Helmholtz equation discretized with the centered 3-point stencil (see [14]).

We compute the eigenvalues and eigenfunctions associated to the continuous 1d-Maxwell's operator with homogeneous PEC boundary conditions and the eigenvalues and eigenvectors of the discretized 1d Maxwell equation using the Yee scheme. In addition to the aforementioned result regarding the invertibility of the discrete Maxwell's equation, the explicit formulas for both the eigenvalues and the eigenvectors will be useful for future research. They could indeed be used to get some precise error bounds as those from [29] or help to understand the convergence behavior of some iterative methods like multigrid (see e.g. [13, 22, 30]).
Appendix A.1. Eigenvalues associated to the continuous Maxwell's equations

The continuous problem (2) on $\Omega = (0, 1)$ with Perfect Electric Conductor (PEC) boundary conditions is

$$ \left\{ \begin{array}{l} i\omega\varepsilon E - \partial_z H = -J, \\ i\omega\mu H - \partial_z E = 0, \\ E(0) = 0, \\ E(1) = 0. \end{array} \right. \qquad (A.1) $$

Setting $E = \frac{1}{\sqrt{\varepsilon}}\tilde{E}$, $H = \frac{1}{\sqrt{\mu}}\tilde{H}$, $\tilde{J} = \sqrt{\mu}J$, we get

$$ \left\{ \begin{array}{l} ik\tilde{E} - \partial_z\tilde{H} = -\tilde{J}, \\ ik\tilde{H} - \partial_z\tilde{E} = 0, \\ \tilde{E}(0) = 0, \\ \tilde{E}(1) = 0, \end{array} \right. \qquad (A.2) $$

with the wavenumber $k = \omega\sqrt{\varepsilon\mu}$.

We compute the eigenvalues and eigenfunctions associated to the continuous problem (A.2), looking for $\lambda \in \mathbb{C}$ and $(\varphi(z), \psi(z))$ such that

$$ \left\{ \begin{array}{l} ik\varphi - \partial_z\psi = \lambda\varphi, \\ ik\psi - \partial_z\varphi = \lambda\psi, \\ \varphi(0) = 0, \\ \varphi(1) = 0. \end{array} \right. \qquad (A.3) $$

Using the first equation of (A.3), we note that the homogeneous Dirichlet boundary condition on $\varphi$ gives a homogeneous Neumann boundary condition, hence

$$ \partial_z\psi(0) = \partial_z\psi(1) = 0. $$

Therefore, any solution to (A.3) also satisfies

$$ \left\{ \begin{array}{l} (ik - \lambda)^2 \varphi - \partial_{zz} \varphi = 0, \qquad \varphi(0) = \varphi(1) = 0, \\ (ik - \lambda)^2 \psi - \partial_{zz} \psi = 0, \qquad \partial_z \psi(0) = \partial_z \psi(1) = 0. \end{array} \right. \qquad (A.4) $$

Any non-zero solution to (A.4) can be obtained by setting $\lambda = i\sigma$ and some computations give the eigenvalues $k - \sigma = j\pi$, $j \in \mathbb{Z}$, hence $\lambda_j^\pm = ik \pm ij\pi$, $j \in \mathbb{N}$. In addition, the associated eigenfunctions are necessarily eigenfunctions of the operator $-\partial_{zz}$ with homogeneous Dirichlet or Neumann boundary
conditions given by $\varphi = \sin(j\pi z)$ and $\psi = \cos(j\pi z)$. The eigenfunctions of (A.3) are then

$$ \forall j \in \mathbb{N}: \mathbf{v}_j^{\pm}(z) = \begin{pmatrix} \mp i \sin(j\pi z) \\ \cos(j\pi z) \end{pmatrix}, $$

and this is an orthonormal basis of $L^2((0,1), \mathbb{C}^2)$ as eigenfunctions of the unbounded skew-symmetric operator $\tilde{M} = \begin{pmatrix} ik & -\partial_z \\ -\partial_z & ik \end{pmatrix}$ acting on $H_0^1((0,1), \mathbb{C}) \times L^2((0,1), \mathbb{C})$. Finally, note that (A.2) has a unique solution for any $k$ such that $\lambda_j^{\pm} \neq 0$ hence $k \neq j\pi$ for $j \in \mathbb{Z}$.

**Remark Appendix A.1.** We can compute as well the eigenvalues and eigenfunctions of the operator

$$ \mathcal{M} = \begin{pmatrix} i\omega\epsilon & -\partial_z \\ -\partial_z & i\omega\mu \end{pmatrix}. $$

Indeed, looking for $\lambda \in \mathbb{C}$ and $\mathbf{v}(y) = (\varphi(y), \psi(y))$ such that $\varphi \in H_0^1(0,1)$ and

$$ \mathcal{M}\mathbf{v} = \tilde{\lambda}\mathbf{v}, $$

gives that necessarily

$$(i\omega\epsilon - \tilde{\lambda})(i\omega\mu - \tilde{\lambda})\varphi - \partial_{zz}\varphi = 0, \quad \varphi(0) = \varphi(1) = 0,$$

$$(i\omega\epsilon - \tilde{\lambda})(i\omega\mu - \tilde{\lambda})\psi - \partial_{zz}\psi = 0, \quad \partial_z\psi(0) = \partial_z\psi(1) = 0.$$

As a result, any eigenvalue $\lambda_j$ of $\mathcal{M}$ satisfies the quadratic equation

$$ (i\omega\epsilon - \tilde{\lambda}_j)(i\omega\mu - \tilde{\lambda}_j) = -(j\pi)^2, \quad j \in \mathbb{N}. \qquad (A.5) $$

Solving (A.5) gives

$$ \tilde{\lambda}_j^{\pm} = \frac{i}{2} \left( \omega\epsilon + \omega\mu \pm \sqrt{\omega^2(\epsilon - \mu)^2 + 4j^2\pi^2} \right), \quad j \in \mathbb{N}. $$

In addition, the associated eigenfunctions are $\varphi(y) = \alpha_j^{\pm} \sin(j\pi y)$, $\psi_j(y) = \beta_j^{\pm} \cos(j\pi y)$ for some specific constants $\alpha_j^{\pm}, \beta_j^{\pm} \in \mathbb{C}$ such that

$$ \begin{pmatrix} \alpha_j^{\pm} \\ \beta_j^{\pm} \end{pmatrix} \in \ker (A_j^{\pm}) \quad \text{with} \quad A_j^{\pm} = \begin{pmatrix} (i\omega\epsilon - \tilde{\lambda}_j^{\pm}) & j\pi \\ -j\pi & (i\omega\mu - \tilde{\lambda}_j^{\pm}) \end{pmatrix}. $$

Problem (A.1) then has a unique solution for any angular frequency $\omega$ such that $\tilde{\lambda}_j^{\pm} \neq 0$, that is $\omega \neq \pm \frac{\pi j}{\sqrt{\epsilon\mu}}$, $j \in \mathbb{Z}$.
Appendix A.2. Eigenvalues associated to the discrete problem

We consider a uniform mesh of $[0, 1]$ with $N + 1$ points so that

$$x_j = jh, \quad j = 0, 1, \dots, N, \quad h = \frac{1}{N}.$$

Denoting by $E_j, H_{j+1/2}$ the unknowns associated to the electric and magnetic fields, the 1d Yee scheme is then

$$ \left\{ \begin{array}{l} \hat{i}k E_j - \frac{1}{h}(H_{j+1/2} - H_{j-1/2}) = -J(x_j), \quad j = 1, 2, \dots, N-1, \\ \hat{i}k H_{j+1/2} - \frac{1}{h}(E_{j+1} - E_j) = 0, \quad j = 0, 1, \dots, N-1, \\ E_0 = 0, \\ E_N = 0, \end{array} \right. \qquad (A.6) $$

where $\hat{k}$ is either $k$ for the usual Yee scheme or

$$ \hat{k} = \hat{\omega}\sqrt{\varepsilon\mu} = \omega\sqrt{\varepsilon\mu} \left( \frac{\sin\left(\frac{h\omega\sqrt{\varepsilon\mu}}{2}\right)}{h\omega\sqrt{\varepsilon\mu}} \right) = k \left( \frac{\sin(hk/2)}{hk/2} \right) $$

when using dispersion correction.

We keep only the interior points and introduce the notation

$$ \mathbf{E}_h = (E_1, E_2, \dots, E_{N-1})^T \in \mathbb{C}^{N-1}, \quad \mathbf{H}_h = (H_{1/2}, H_{3/2}, \dots, H_{N-1/2}) \in \mathbb{C}^N, \\ \mathbf{J}_h = (J_1, J_2, \dots, J_{N-1})^T \in \mathbb{C}^{N-1}. $$

The linear system associated to (A.6) is

$$ \left( \begin{array}{cc} \hat{i}kI_{N-1} & L \\ -L^T & \hat{i}kI_N \end{array} \right) \left( \begin{array}{c} \mathbf{E}_h \\ \mathbf{H}_h \end{array} \right) = \left( \begin{array}{c} -\mathbf{J}_h \\ O \end{array} \right), \qquad (A.7) $$

where $L \in \mathbb{C}^{(N-1) \times N}$ is the bi-diagonal rectangular matrix given by

$$ L = \frac{1}{h} \begin{pmatrix} 1 & -1 & 0 & 0 & 0 & 0 \\ 0 & 1 & -1 & 0 & 0 & 0 \\ 0 & 0 & \cdots & \cdots & 0 & 0 \\ 0 & 0 & 0 & \cdots & \cdots & 0 \\ 0 & 0 & 0 & 0 & 1 & -1 \end{pmatrix}. $$
We are going to diagonalize the matrix

$$ \mathcal{M}_h = \begin{pmatrix} i\hat{k}I_{N-1} & L \\ -L^T & i\hat{k}I_N \end{pmatrix}, $$

appearing in (A.7). It is worth noting that $\mathcal{M}_h = i\hat{k}I_{2N-1} + K$, where $K$ satisfies $K^* = -K$ which translates to $(iK)^* = iK$, and we get that the matrix $iK$ is Hermitian and thus it admits a family of eigenvectors that is an orthonormal basis of $\mathbb{C}^{2N-1}$. We compute below both this orthonormal family of eigenvectors and the associated eigenvalues as well.

**Theorem Appendix A.2.** The spectrum of $\mathcal{M}_h$ is given by $\sigma(\mathcal{M}_h) = \{i\hat{k}, \lambda_{h,j}^{\pm}\}$ where $\lambda_{h,j}^{\pm}$ are

$$ \lambda_{h,j}^{\pm} = i\hat{k} \pm i\sqrt{\sigma_{h,j}}, \quad \sigma_{h,j} = \frac{4}{\hbar^2} \sin\left(\frac{j\pi h}{2}\right)^2, \quad j = 1, \dots, N-1. $$

The eigenvector associated to the eigenvalue $i\hat{k}$ is

$$ \mathbf{v}_0 = \sqrt{\hbar} \begin{pmatrix} 0 \\ \mathbf{1}_N \end{pmatrix}, $$

and each other eigenvector can be written as

$$ \mathbf{v}_j^{\pm} = \frac{1}{\sqrt{2}} \begin{pmatrix} \mp i \mathbf{v}_{E,j} \\ \mathbf{v}_{H,j} \end{pmatrix}, \quad j = 1, \dots, N-1, $$

where

$$ \mathbf{v}_{E,j} = \sqrt{2\hbar} [\sin(jl\pi h)]_{l=1}^{N-1}, \quad \mathbf{v}_{H,j} = \sqrt{2\hbar} \left[ \cos\left(\left(\frac{h}{2} + lh\right)\pi j\right) \right]_{l=0}^{N-1} . \quad (A.8) $$

*Proof:* We start by considering the eigenvalue problem $\mathcal{M}_h\mathbf{v} = \lambda\mathbf{v}$ which, setting $\mathbf{v} = (\mathbf{v}_E, \mathbf{v}_H)^T \in \mathbb{C}^{(N-1)+N}$, can be recast as

$$ i\hat{k}\mathbf{v}_E + L\mathbf{v}_H = \lambda\mathbf{v}_E, \quad i\hat{k}\mathbf{v}_H - L^T\mathbf{v}_E = \lambda\mathbf{v}_H. $$

Inserting $(i\hat{k} - \lambda)\mathbf{v}_E = -L\mathbf{v}_H$ into the second equation and $(i\hat{k} - \lambda)\mathbf{v}_H = L^T\mathbf{v}_E$ into the first one, we obtain

$$ LL^T \mathbf{v}_E = -(i\hat{k} - \lambda)^2 \mathbf{v}_E, \quad L^T L \mathbf{v}_H = -(i\hat{k} - \lambda)^2 \mathbf{v}_H. $$
This shows that it is necessary for $\mathbf{v}_E$ to be an eigenvector of $LL^T$ and $\mathbf{v}_H$ to be an eigenvector of $L^T L$. Since

$$LL^T = \frac{1}{h^2} \text{tridiag}(-1, 2, -1) \in \mathbb{R}^{(N-1) \times (N-1)},$$

$$L^T L = \frac{1}{h^2} \text{tridiag}(-1, \mathbf{u}_N, -1) \in \mathbb{R}^{N \times N}, \quad \mathbf{u}_N = (1, 2, \dots, 2, 1),$$

we get that $\mathbf{v}_E$ and $\mathbf{v}_H$ are given as in (A.8). In addition, any eigenvalue $\lambda$ of $\mathcal{M}_h$ must satisfy

$$\forall \gamma \in \sigma(LL^T) \cup \sigma(L^T L) : - (i\hat{k} - \lambda)^2 = \gamma.$$

Using the following trigonometric identities valid for any integer $l$,

$$\frac{\sin((l+1)x) - \sin(lx)}{\cos((l+\frac{1}{2})x)} = 2\sin(\frac{x}{2}),$$

$$\frac{-\cos((l+\frac{1}{2})x) + \cos((l-\frac{1}{2})x)}{\sin(lx)} = 2\sin(\frac{x}{2}),$$

with $x = j\pi h$, we obtain the formulas

$$L \mathbf{v}_{H,j} = \sqrt{\sigma_j} \mathbf{v}_{E,j}, \quad L^T \mathbf{v}_{E,j} = \sqrt{\sigma_j} \mathbf{v}_{H,j},$$

from which we infer that

$$LL^T \mathbf{v}_{E,j} = \sigma_j \mathbf{v}_{E,j}, \quad j = 1, \dots, N-1 \text{ and } L^T L \mathbf{v}_{H,j} = \sigma_j \mathbf{v}_{H,j}, \quad j = 1, \dots, N-1.$$

Using that $L^T L \mathbf{v}_0 = \mathbf{0}_N$, we obtain that

$$\sigma(LL^T) \cup \sigma(L^T L) = \{0\} \cup \{\sigma_j, \quad j = 1, \dots, N-1\},$$

and thus the spectrum of $\mathcal{M}_h$ is

$$\sigma(\mathcal{M}_h) = \{i\hat{k}\} \cup \{\lambda_{h,j}^{\pm} := i\hat{k} \pm i\sqrt{\sigma_{h,j}}, \quad j = 1, \dots, N-1\}.$$

Since the eigenvector $\mathbf{v}_j^{\pm}$ associated to $\lambda_{h,j}^{\pm}$ has to be written as $(\alpha_j^{\pm} \mathbf{v}_E, \beta_j^{\pm} \mathbf{v}_H)^T$, it now remains to get the constants $\alpha_j^{\pm}$ and $\beta_j^{\pm}$. The constant satisfy the following $2 \times 2$ linear system

$$A \begin{pmatrix} \alpha_j^{\pm} \\ \beta_j^{\pm} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} \quad \text{with} \quad A^{\pm} = \begin{pmatrix} (i\hat{k} - \lambda_j^{\pm}) & \sqrt{\sigma_j} \\ -\sqrt{\sigma_j} & (i\hat{k} - \lambda_j^{\pm}) \end{pmatrix} = \sqrt{\sigma_j} \begin{pmatrix} \mp i & 1 \\ -1 & \mp i \end{pmatrix}.$$
Since this linear system is singular, one gets $(\alpha_j^{\pm}, \beta_j^{\pm})^T \in \ker(A^{\pm})$ hence

$$ \begin{pmatrix} \alpha_j^{\pm} \\ \beta_j^{\pm} \end{pmatrix} = \begin{pmatrix} \mp i \\ 1 \end{pmatrix}, $$

which completes the proof. □

We emphasize that the families $\{v_{E,j}\}_{j=1}^{N-1}$ and $\{v_0\} \cup \{v_{H,j}\}_{j=1}^{N-1}$ are orthonormal bases of $\mathbb{R}^{N-1}$ and $\mathbb{R}^N$. Therefore the family of eigenvectors of $M_h$ computed in Theorem Appendix A.2 defined as

$$ v_j^{\pm} = \frac{1}{\sqrt{2}} \begin{pmatrix} \mp i v_{E,j} \\ v_{H,j} \end{pmatrix}, \quad j = 1, \dots, N-1 $$

is an orthonormal basis of $\mathbb{C}^{2N-1}$.

We end this section by showing the effect of the shifted angular frequency on the invertibility of the linear system (A.7). Without dispersion correction, this gives

$$ k \neq \pm \sqrt{\sigma_{h,j}}, \quad j = 1, \dots, N-1, $$

whereas with dispersion correction, one has

$$ \hat{k} \neq \pm \sqrt{\sigma_{h,j}} = \frac{2}{h} \sin \left( \pm \frac{j\pi h}{2} \right), \quad j = 1, \dots, N-1. $$

Using the expression for $\hat{k}$, this yields

$$ \frac{2}{h} \sin \left( \frac{kh}{2} \right) \neq \frac{2}{h} \sin \left( \pm \frac{j\pi h}{2} \right), \quad j = 1, \dots, N-1, $$

which translates to

$$ k \neq \pm j\pi, \quad j = 1, \dots, N-1, $$

and hence the shifted angular frequency ensures that the linear system (A.7) is invertible whenever the continuous problem (A.1) is.

**Remark Appendix A.3.** Theorem Appendix A.2 can also be obtained for the discrete problem associated to (A.1) whose associated matrix is

$$ \tilde{M}_h = \begin{pmatrix} i\hat{\omega}\epsilon I_{N-1} & L \\ -L^T & i\hat{\omega}\mu I_N \end{pmatrix}. $$
In this setting, we have that any $\tilde{\lambda} \in \sigma(\tilde{\mathcal{M}}_h)$ satisfies $-(i\hat{\omega}\varepsilon - \tilde{\lambda})(i\hat{\omega}\mu - \tilde{\lambda}) = \sigma_j$, for all $j = 1, \dots, N-1$, which gives $2N-2$ eigenvalues. The last one can be obtained by using that

$$ \begin{pmatrix} i\hat{\omega}(\varepsilon - \mu)I_{N-1} & L \\ -L^T & 0 \end{pmatrix} \mathbf{v}_0 = \mathbf{0}, $$

from which we see that $\tilde{\lambda} = i\hat{\omega}\mu$ is also an eigenvalue.

Following the proof of Theorem Appendix A.2, any eigenvector $\tilde{\mathbf{v}}_j^\pm$ of $\tilde{\mathcal{M}}_h$ associated to $\tilde{\lambda}_j^\pm$ can be expressed as $(\tilde{\alpha}_j^\pm \mathbf{v}_E, \tilde{\beta}_j^\pm \mathbf{v}_H)^T$ where $(\tilde{\alpha}_j^\pm, \tilde{\beta}_j^\pm)^T \in \ker(\tilde{\mathcal{A}}^\pm)$ with

$$ \tilde{\mathcal{A}}^{\pm} = \begin{pmatrix} (i\hat{w}\varepsilon - \tilde{\lambda}_j^{\pm}) & \sqrt{\sigma_j} \\ -\sqrt{\sigma_j} & (i\hat{w}\mu - \tilde{\lambda}_j^{\pm}) \end{pmatrix}. $$

Without dispersion correction, the matrix $\tilde{\mathcal{M}}_h$ is invertible whenever $\omega$ satisfies

$$ \omega \neq \pm \frac{\sigma_{h,j}}{\sqrt{\varepsilon\mu}}, \quad j = 1, \dots, N-1. $$

When using dispersion correction, the matrix $\tilde{\mathcal{M}}_h$ is invertible as soon as

$$ \hat{\omega} \neq \pm \frac{\sigma_{h,j}}{\sqrt{\varepsilon\mu}}, \quad j = 1, \dots, N-1 $$

from which we see that

$$ \omega \neq \pm \frac{\pi j}{\sqrt{\varepsilon\mu}}, $$

which are exactly the frequencies for which the continuous problem (A.1) has a unique solution.

## Appendix B. Computations with Maple: Continuous and discrete symbols

We give here the Maple code to compute the continuous and discrete symbols. The continuous symbol $\mathcal{L}(\omega, \mathbf{k})$ is computed as
restart;
with(LinearAlgebra);
Nx := Matrix([[0, 0, 0], [0, 0, 1], [0, -1, 0]]);
Ny := Matrix([[0, 0, -1], [0, 0, 0], [1, 0, 0]]);
Nz := Matrix([[0, 1, 0], [-1, 0, 0], [0, 0, 0]]);
Z := Matrix([[0, 0, 0], [0, 0, 0], [0, 0, 0]]);
Id := Matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]]);
Gx := Matrix([[Z, Nx], [-Nx, Z]]);
Gy := Matrix([[Z, Ny], [-Ny, Z]]);
Gz := Matrix([[Z, Nz], [-Nz, Z]]);
Gw := Matrix([[epsilon*Id, Z], [Z, mu*Id]]);
G := w*Gw*I + kx*Gx*I + ky*Gy*I + kz*Gz*I;

The discrete symbol is denoted by $\mathcal{L}(h, \omega, \mathbf{k})$ and can be computed as

$$
\begin{align*}
\mathrm{Ex} &:= \exp(kx\cdot h\cdot I/2); \quad \mathrm{Ey} := \exp(ky\cdot h\cdot I/2); \quad \mathrm{Ez} := \exp(kz\cdot h\cdot I/2); \\
\mathrm{sx} &:= 2\cdot I\cdot \sin(h/2\cdot kx)/h; \quad \mathrm{sy} := 2\cdot I\cdot \sin(h/2\cdot ky)/h; \quad \mathrm{sz} := 2\cdot I\cdot \sin(h/2\cdot kz)/h; \\
\mathrm{Gd} &:= \begin{bmatrix}
\text{w}\cdot\text{epsilon}\cdot\text{Ex}\cdot I, & 0, & 0, & 0, & \text{sz}\cdot\text{Ex}, & -\text{sy}\cdot\text{Ex}, \\
0, & \text{w}\cdot\text{epsilon}\cdot\text{Ey}\cdot I, & 0, & -\text{sz}\cdot\text{Ey}, & 0, & \text{sx}\cdot\text{Ey}, \\
0, & 0, & \text{w}\cdot\text{epsilon}\cdot\text{Ez}\cdot I, & \text{sy}\cdot\text{Ez}, & -\text{sx}\cdot\text{Ez}, & 0, \\
0, & -\text{sz}\cdot\text{Ey}\cdot\text{Ez}, & \text{sy}\cdot\text{Ey}\cdot\text{Ez}, & \text{w}\cdot\text{mu}\cdot\text{Ey}\cdot\text{Ez}\cdot I, & 0, & 0, \\
[\text{sz}\cdot\text{Ex}\cdot\text{Ez}, 0, & -\text{sx}\cdot\text{Ex}\cdot\text{Ez}, 0, & \text{w}\cdot\text{mu}\cdot\text{Ex}\cdot\text{Ez}\cdot I, 0, & 0, & 0, \\
[-\text{sy}\cdot\text{Ex}\cdot\text{Ey}, \text{sx}\cdot\text{Ex}\cdot\text{Ey}, 0, 0, 0, 0, \text{w}\cdot\text{mu}\cdot\text{Ex}\cdot\text{Ey}\cdot I]]
\end{bmatrix};
\end{align*}
$$

The continuous (denoted as $S_c = \det(\mathcal{L}(\omega, \mathbf{k}))$ and discrete ($S_d = \det(\mathcal{L}(h, \omega, \mathbf{k}))$)
dispersion relations are then obtained with

Sc := simplify(Determinant(G));
Sd := simplify(Determinant(Gd));

## Appendix C. Computations with Maple: Discrete wavenumber for 3d Yee scheme

Here are the Maple commands to get the expansion of the discrete wavenumber for small enough meshsize:

kx := kd*cos(phi)*sin theta;
ky := kd*sin(phi)*sin theta;
kz := kd*cos theta;
kd := k0 + k2*h^2;
Expansion_Sd := convert(simplify serie(Sd, h = 0, 10)), polynom, h):
k0 sol := solve(coeff(Expansion_Sd, h, 0), k0); k0 := k0 sol [1];
k2 sol := solve(coeff(Expansion_Sd, h, 4), k2);
k2:=simplify(k2 sol [1]);

In these commands, we chose to keep $k_0 = k_{0sol}[1]$ since this gives $w\sqrt{\varepsilon\mu}$ which is the (positive) continuous wavenumber. To get $k_2$, we have to get it using the fourth order coefficient in the expansion, as $h \to 0$, of $S_d$ = det($\mathcal{L}(h, k_d\theta)$) since the other lower order terms vanish once $k_0$ is assigned.

## Appendix D. Manufactured solution

Here are the Maple commands used to compute the manufactured solution $E^{\text{ex}}$, $H^{\text{ex}}$:

E1 := sin(z)*sin(y)*sin(k*x);
E2 := sin(z)*sin(x)*sin(k*y);
E3 := sin(y)*sin(x)*sin(k*z);
H1 := simplify(1/(w*mu*I)*(diff(E2, z) - diff(E3, y)));
H2 := simplify(1/(w*mu*I)*(diff(E3, x) - diff(E1, z)));
H3 := simplify(1/(w*mu*I)*(diff(E1, y) - diff(E2, x)));

The associated electric current $J^{\text{ex}}$ can then be computed as

JE1 := -simplify(w*e*E1*I + (-diff(H3, y) + diff(H2, z)));
JE2 := -simplify(w*e*E2*I + (-diff(H1, z) + diff(H3, x)));
JE3 := -simplify(w*e*E3*I + (-diff(H2, x) + diff(H1, y)));

## References

[1] Ahmad, F., Amin, M. G., and Dogaru, T. (2014). Partially sparse imaging of stationary indoor scenes. *EURASIP Journal on Advances in Signal Processing*, 2014:1–15.

[2] Ainsworth, M. (2004a). Discrete dispersion relation for hp-version finite element approximation at high wave number. *SIAM Journal on Numerical Analysis*, 42(2):553–575.

[3] Ainsworth, M. (2004b). Dispersive properties of high-order Nédélec/edge element approximation of the time-harmonic Maxwell equations. *Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences*, 362(1816):471–491.
[4] Babuska, I. M. and Sauter, S. A. (1997). Is the pollution effect of the fem avoidable for the Helmholtz equation considering high wavenumbers? *SIAM Journal on Numerical Analysis*, 34(6):2392-2423.

[5] Bendali, A. (2023). Dispersion reduction in Feng and Wu's IPDG method. *Computer Methods in Applied Mechanics and Engineering*, 416:116388.

[6] Bonazzoli, M., Dolean, V., Graham, I., Spence, E., and Tournier, P.-H. (2019). Domain decomposition preconditioning for the high-frequency time-harmonic Maxwell equations with absorption. *Mathematics of Computation*, 88(320):2559-2604.

[7] Bouajaji, M. E., Dolean, V., Gander, M. J., and Lanteri, S. (2012). Optimized schwarz methods for the time-harmonic Maxwell equations with damping. *SIAM Journal on Scientific Computing*, 34(4):A2048-A2071.

[8] Chen, Z., Cheng, D., Feng, W., and Wu, T. (2013). An optimal 9-point finite difference scheme for the Helmholtz equation with PML. *International Journal of Numerical Analysis & Modeling*, 10(2).

[9] Chen, Z., Cheng, D., and Wu, T. (2012). A dispersion minimizing finite difference scheme and preconditioned solver for the 3d Helmholtz equation. *Journal of Computational Physics*, 231(24):8152-8175.

[10] Cheng, D., Tan, X., and Zeng, T. (2017). A dispersion minimizing finite difference scheme for the Helmholtz equation based on point-weighting. *Computers & Mathematics with Applications*, 73(11):2345-2359.

[11] Chew, W. C. (2020). Lectures on electromagnetic field theory. USA, Purdue University.

[12] Cho, K. (2019). Dispersion relation in chiral media: Credibility of Drude-Born-Fedorov equations. In *Electromagnetic Metamaterials*, pages 189-199. Springer.

[13] Cocquet, P.-H. and Gander, M. J. (2017). How large a shift is needed in the shifted Helmholtz preconditioner for its effective inversion by multi-grid? *SIAM Journal on Scientific Computing*, 39(2):A438-A478.

[14] Cocquet, P.-H. and Gander, M. J. (2024). Asymptotic dispersion correction in general finite difference schemes for Helmholtz problems. *SIAM Journal on Scientific Computing*, 46(2):A670-A696.
[15] Cocquet, P.-H., Gander, M. J., and Xiang, X. (2018). A finite difference method with optimized dispersion correction for the Helmholtz equation. In *Domain Decomposition Methods in Science and Engineering XXIV 24*, pages 205–213. Springer.

[16] Cocquet, P.-H., Gander, M. J., and Xiang, X. (2020). Dispersion correction for Helmholtz in 1d with piecewise constant wavenumber. In *Domain Decomposition Methods in Science and Engineering XXV 25*, pages 359–366. Springer.

[17] Cocquet, P.-H., Gander, M. J., and Xiang, X. (2021). Closed form dispersion corrections including a real shifted wavenumber for finite difference discretizations of 2d constant coefficient helmholtz problems. *SIAM Journal on Scientific Computing*, 43(1):A278–A308.

[18] Cocquet, P.-H., Mazet, P.-A., and Mouysset, V. (2012). On the existence and uniqueness of a solution for some frequency-dependent partial differential equations coming from the modeling of metamaterials. *SIAM Journal on Mathematical Analysis*, 44(6):3806–3833.

[19] Dastour, H. and Liao, W. (2019). A fourth-order optimal finite difference scheme for the Helmholtz equation with pml. *Computers & Mathematics with Applications*, 78(6):2147–2165.

[20] Deraemaeker, A., Babuška, I., and Bouillard, P. (1999). Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions. *International journal for numerical methods in engineering*, 46(4):471–499.

[21] Dwarka, V. and Vuik, C. (2021). Pollution and accuracy of solutions of the Helmholtz equation: A novel perspective from the eigenvalues. *Journal of Computational and Applied Mathematics*, 395:113549.

[22] Ernst, O. G. and Gander, M. J. (2013). Multigrid methods for Helmholtz problems: A convergent scheme in 1d using standard components. *Direct and Inverse Problems in Wave Propagation and Applications. De Gruyer*, pages 135–186.

[23] Fathy, A., Wang, C., Wilson, J., and Yang, S. (2008). A fourth order difference scheme for the Maxwell equations on Yee grid. *Journal of Hyperbolic Differential Equations*, 5(03):613–642.
[24] Feng, Q., Han, B., and Michelle, M. (2021). Sixth order compact finite difference method for 2d Helmholtz equations with singular sources and reduced pollution effect. *arXiv preprint arXiv:2112.07154*.

[25] Feng, X., Lu, P., and Xu, X. (2016). A hybridizable discontinuous Galerkin method for the time-harmonic Maxwell equations with high wavenumber. *Computational Methods in Applied Mathematics*, 16(3):429–445.

[26] Feng, X. and Wu, H. (2014). An absolutely stable discontinuous Galerkin method for the indefinite time-harmonic Maxwell equations with large wavenumber. *SIAM Journal on Numerical Analysis*, 52(5):2356–2380.

[27] Fernandes, P. and Raffetto, M. (2009). Well-posedness and finite element approximability of time-harmonic electromagnetic boundary value problems involving bianisotropic materials and metamaterials. *Mathematical Models and Methods in Applied Sciences*, 19(12):2299–2335.

[28] Fezoui, L. F. and Lanteri, S. (2018). *Finite volume scheme for the 1D Maxwell equations with GSTC conditions*. PhD thesis, Inria.

[29] Gander, M. J. and Zhang, H. (2025). Fourier analysis of finite difference schemes for the Helmholtz equation: Sharp estimates and relative errors. *arXiv preprint arXiv:2501.16696*.

[30] Hackbusch, W. (2013). *Multi-grid methods and applications*, volume 4. Springer Science & Business Media.

[31] Han, B., Michelle, M., and Wong, Y. S. (2021). Dirac assisted tree method for 1d heterogeneous Helmholtz equations with arbitrary variable wavenumbers. *Computers & Mathematics with Applications*, 97:416–438.

[32] Ihlenburg, F. and Babuška, I. (1995). Finite element solution of the Helmholtz equation with high wavenumber part i: The h-version of the FEM. *Computers & Mathematics with Applications*, 30(9):9–37.

[33] Jagalur-Mohan, J., Feijoo, G., and Oberai, A. (2013). A Galerkin least squares method for time harmonic Maxwell equations using Nédélec elements. *Journal of Computational Physics*, 235:67–81.
[34] Liu, Z.-l., Song, P., Li, J.-s., Li, J., and Zhang, X.-b. (2015). An optimized implicit finite-difference scheme for the two-dimensional Helmholtz equation. *Geophysical Journal International*, 202(3):1805–1826.

[35] Lytaev, M. S. (2023). Reducing the numerical dispersion of the one-way Helmholtz equation via the differential evolution method. *Journal of Computational Science*, 71:102057.

[36] Melenk, J. M. and Sauter, S. (2011). Wavenumber explicit convergence analysis for Galerkin discretizations of the Helmholtz equation. *SIAM Journal on Numerical Analysis*, 49(3):1210–1243.

[37] Michelle, M. (2022). Numerical study of the Helmholtz equation with large wavenumbers.

[38] Monk, P. (2003). *Finite element methods for Maxwell's equations*. Oxford university press.

[39] Monk, P. and Süli, E. (1994). A convergence analysis of Yee's scheme on nonuniform grids. *SIAM Journal on Numerical Analysis*, 31(2):393–412.

[40] Nédélec, J.-C. (1980). Mixed finite elements in $\mathbb{R}^3$. *Numerische Mathematik*, 35:315–341.

[41] Rewienski, M. and Mrozowski, M. (2000). An iterative algorithm for reducing dispersion error on Yee's mesh in cylindrical coordinates. *IEEE microwave and guided wave letters*, 10(9):353–355.

[42] Spence, E. A. (2023). A simple proof that the hp-FEM does not suffer from the pollution effect for the constant-coefficient full-space Helmholtz equation. *Advances in Computational Mathematics*, 49(2):27.

[43] Stolk, C. C. (2016). A dispersion minimizing scheme for the 3-d Helmholtz equation based on ray theory. *Journal of Computational Physics*, 314:618–646.

[44] Stolk, C. C., Ahmed, M., and Bhowmik, S. K. (2014). A multigrid method for the Helmholtz equation with optimized coarse grid corrections. *SIAM Journal on Scientific Computing*, 36(6):A2819–A2841.
[45] Taflove, A. (1988). Review of the formulation and applications of the finite-difference time-domain method for numerical modeling of electromagnetic wave interactions with arbitrary structures. *Wave Motion*, 10(6):547-582.

[46] Tsuji, P. and Ying, L. (2012). A sweeping preconditioner for Yee's finite difference approximation of time-harmonic Maxwell's equations. *Frontiers of Mathematics in China*, 7:347-363.

[47] Turkel, E. and Yefet, A. (2000). On the construction of a high order difference scheme for complex domains in a cartesian grid. *Applied Numerical Mathematics*, 33(1-4):113-124.

[48] Versano, I., Turkel, E., and Tsynkov, S. (2024). Fourth-order accurate compact scheme for first-order Maxwell's equations. *Journal of Scientific Computing*, 100(2):31.

[49] Wang, K. and Wong, Y. S. (2014). Pollution-free finite difference schemes for non-homogeneous Helmholtz equation. *International Journal of Numerical Analysis & Modeling*, 11(4).

[50] Wu, T. (2017). A dispersion minimizing compact finite difference scheme for the 2d Helmholtz equation. *Journal of Computational and Applied Mathematics*, 311:497-512.

[51] Wu, T. and Chen, Z. (2014). A dispersion minimizing subgridding finite difference scheme for the Helmholtz equation with PML. *Journal of Computational and Applied Mathematics*, 267:82-95.

[52] Wu, T., Sun, Y., and Cheng, D. (2021). A new finite difference scheme for the 3d Helmholtz equation with a preconditioned iterative solver. *Applied Numerical Mathematics*, 161:348-371.

[53] Wu, T. and Xu, R. (2018). An optimal compact sixth-order finite difference scheme for the Helmholtz equation. *Computers & Mathematics with Applications*, 75(7):2520-2537.

[54] Yee, K. (1966). Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. *IEEE Transactions on antennas and propagation*, 14(3):302-307.
[55] Zhao, R., Koschny, T., and Soukoulis, C. M. (2010). Chiral metamaterials: retrieval of the effective parameters with and without substrate. *Optics express*, 18(14):14553–14567.

[56] Zhou, Y. and Wu, H. (2023). Dispersion analysis of CIP-FEM for the Helmholtz equation. *SIAM Journal on Numerical Analysis*, 61(3):1278–1292.

[57] Zhu, L. and Wu, H. (2013). Preasymptotic error analysis of CIP-FEM and FEM for Helmholtz equation with high wavenumber. part ii: hp version. *SIAM Journal on Numerical Analysis*, 51(3):1828–1852.