# Efficient perfect matched layer for hybrid MRTD-FDTD computation with lossy dielectric material boundaries

Grigory Ibrahim Massy Sanchez, Nestor Pena, Michel Ney

► To cite this version:

Grigory Ibrahim Massy Sanchez, Nestor Pena, Michel Ney. Efficient perfect matched layer for hybrid MRTD-FDTD computation with lossy dielectric material boundaries. European Physical Journal: Applied Physics, 2012, 57 (01), 10.1051/epjap/2011110113. [hal-00853291](https://hal.science/hal-00853291v1)

HAL Id: hal-00853291

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

Submitted on 18 Feb 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
# Efficient perfect matched layer for hybrid MRTD-FDTD computation with lossy dielectric material boundaries

I. Massy¹,ᵃ, N. Peña¹, and M.M. Ney²

¹ Group of Electronic and Telecommunications Systems, Universidad de Los Andes, Bogotá, Colombia

² Lab-STICC, CNRS, Telecom Institute, Telecom Bretagne, Brest, France

**Abstract.** Existing implementation of perfectly matched layers (PML) constraining the computational volume in an electromagnetic modeling of open structures has been shown to be efficient. However, cases with absorbing materials inserted into PML require substantial amount of memory and calculation time spent on the computation of the field inside the PML layer. In this paper we are addressing these problems by introducing a new computational method. As the procedure requires some substantial amount of memory and operations to compute fields in the PML layer, the paper addresses these issues. The PML split-field equations are revisited yielding a reduction from four-term storage to only three. In addition, the new algorithm is used with a new hybrid MRTD-FDTD (HMRTD) method applied to various lossy structures, conjointly with a higher-order field evaluation. Results demonstrate the efficiency of the new PML algorithm as compared to current schemes.

## 1 Introduction

Special media derived from stretched coordinates [1] are used to limit a computational domain in which volumic approaches such as Finite Difference Time-Domain (FDTD), Transmission-Line Matrix (TLM) or Finite-Element (FEM) techniques are implemented. These media provide good absorption and low level of reflection of electromagnetic waves at arbitrary frequency and incidence on the limits of computational domain. The PML-split version proposed by Berenger [1] is one of the most used methods. Unsplit versions [2-4] were used in FDTD and Multi-Resolution in Time-Domain (MRTD) schemes [5]. The various PML versions differ by their formalism. Their performances have been shown to be similar for MRTD and FDTD implementation. However, Berenger's formulation is slightly simpler than others. In some situations, lossy media extend to the PML region to simulate infinite dimension or matched load. As a result, the PLM algorithm must be modified, which leads to some additional complexity. As Berenger's formalism is simpler, it was selected here to account for lossy heterogeneous media. To the best of our knowledge, this has not been implemented for MRTD simulations. Concerning FDTD computation with lossy PML boundaries, some pioneer work was published by Teixeira et al. [6]. The modified algorithm implies the storage of additional components to perform convolutions to compute spatial derivatives in Maxwell curl equations. This imposes an additional term

to every field component, to store the time integral of the fields. Recently, it was shown [7] that only three electric field components need to be stored instead of four, as done in [6]. In addition, the implementation of the PML update equations has a first-order precision in current techniques as compared to the proposed technique which provides a second-order accuracy.

The new PML algorithm is used with a hybrid FDTD-MRTD (HMRTD) field computation technique [8], to include optimal variable mesh implementation. First, an example in a rectangular waveguide with high level of losses is presented. Then, the case of a Non-Radiating Dielectric (NRD) waveguide is presented and compared with analytical results when relatively high level of losses occurs.

## 2 PML formalism for lossy media

The basic formalism, introduced in [1], is an interpretation of a complex stretched coordinate system, in which a propagating mode becomes amplitude decreasing due to the imaginary parts of the positions. It is possible to implement an arbitrary boundary at the end of the layer neglecting the reflection, because the complex coordinates will absorb all of the energy of the outgoing wave, before it reaches the actual computational domain.

For a lossy medium, Maxwell's curl equations become:

$$(j\omega\varepsilon + \sigma)\mathbf{E} = \nabla \times \mathbf{H}, \quad (1)$$

$$-j\omega\mu\mathbf{H} = \nabla \times \mathbf{E}. \quad (2)$$

<footnote>ᵃ e-mail: gMassy@yahoo.com</footnote>
When inserting stretched coordinates, a split variation of the original Maxwell equations takes the following form:

$$ (j\omega\varepsilon + \sigma)E_{\zeta\xi} = \frac{1}{1 - \frac{j\alpha_{\xi}}{\omega}} \frac{\partial H_{\eta}}{\partial \xi}, \quad (3) $$

$$ (j\omega\varepsilon + \sigma)E_{\zeta\eta} = -\frac{1}{1 - \frac{j\alpha_{\eta}}{\omega}} \frac{\partial H_{\xi}}{\partial \eta}, \quad (4) $$

where $(\xi, \eta, \zeta) \in \{(x, y, z), (y, z, x), (z, x, y)\}$.

According to Teixeira et al. [6], the above equations become in time domain:

$$ \varepsilon \frac{\partial E_{\zeta\xi}}{\partial t} + (\sigma + \varepsilon\alpha_{\xi})E_{\zeta\xi} + \sigma\alpha_{\xi} \int_{0}^{t} E_{\zeta\xi} dt = \frac{\partial H_{\eta}}{\partial \xi}, \quad (5) $$

$$ \varepsilon \frac{\partial E_{\zeta\eta}}{\partial t} + (\sigma + \varepsilon\alpha_{\eta})E_{\zeta\eta} + \sigma\alpha_{\eta} \int_{0}^{t} E_{\zeta\eta} dt = -\frac{\partial H_{\xi}}{\partial \eta}. \quad (6) $$

Computational implementation of equations (5) and (6) needs to store four terms to update the $E_{\zeta}$ field component: two for $E_{\zeta\xi}$ and $E_{\zeta\eta}$ and two others for their integrals specified in equations (5) and (6). As a result, each field sample in the PML medium requires the storage of four terms. However, as suggested in [7], the fourth component can be replaced by only three, two of them following the lossless Berenger's scheme and the third assuming the term corresponding to the conduction current, as shown in equation (7):

$$ j\omega\varepsilon E_{\zeta} = -\sigma E_{\zeta} + \frac{1}{1 - \frac{j\alpha_{\xi}}{\omega}} \frac{\partial H_{\eta}}{\partial \xi} - \frac{1}{1 - \frac{j\alpha_{\eta}}{\omega}} \frac{\partial H_{\xi}}{\partial \eta}. \quad (7) $$

This can be made because in Maxwell's frequency equation only three scaling factors appear in the PML formulation as shown in equation (7). Thus, separating each term one has:

$$ \varepsilon \frac{\partial E_{\zeta\xi}}{\partial t} + \varepsilon\alpha_{\xi} E_{\zeta\xi} = \frac{\partial H_{\eta}}{\partial \xi}, \quad (8.1) $$

$$ \varepsilon \frac{\partial E_{\zeta\eta}}{\partial t} + \varepsilon\alpha_{\eta} E_{\zeta\eta} = -\frac{\partial H_{\xi}}{\partial \eta}, \quad (8.2) $$

$$ \varepsilon \frac{\partial E_{\zeta\zeta}}{\partial t} + \sigma E_{\zeta} = 0, \quad (9) $$

$$ E_{\zeta} = E_{\zeta\zeta} + E_{\zeta\xi} + E_{\zeta\eta}. \quad (10) $$

The scheme to update the magnetic field remains the same as Berenger's algorithm [1]. The reason is that only dielectric losses are concerned here. Thus, there are only three subcomponents to consider in the new scheme and total field is computed via equation (10).

### 3 Implementation of PML including absorbing materials

The above scheme is used to describe the evolution in time of both the electric and magnetic fields. The evolution equations are discretized forms of Maxwell's modified equations (7)–(9) for electric field and original Berenger

<figure><img src="image_2.png" /><figcaption>Fig. 1. (Color online) Field component spatial sampling for both FDTD and MRTD cells.</figcaption></figure>

for magnetic field, depending on field sample location in space and time. Yee scheme [9] approximates the partial derivatives by its corresponding finite differences, a method popularly known as FDTD. Krumpholz and Katehi [5] makes a similar analysis in the MRTD scheme. A useful form, proposed by Taflove and Hagness [10], is to mesh the computational volume into parallelepipedic cells following the scheme of Figure 1.

Differential and/or integral operators are represented by convolution operators in a discrete form. The quality of the discretized operator has some impact on the convergence error. In [9], due to the placement of field components, one can guarantee second-order error. Velocity error (model dispersion) occurs if cell size has maximum dimension greater than one tenth of the minimum wavelength [5].

Denoting boundary cells as integer positions and center cells as half-integer positions, one has for the PML scheme:

$$ \begin{aligned} & \varepsilon \frac{E_{\zeta\xi}|_{n_{\xi}, n_{\eta}, n_{\zeta}+1/2}^{k+1} - E_{\zeta\xi}|_{n_{\xi}, n_{\eta}, n_{\zeta}+1/2}^{k}}{\Delta t} \\ & + \varepsilon\alpha_{\xi}(n_{\xi}) \frac{E_{\zeta\xi}|_{n_{\xi}, n_{\eta}, n_{\zeta}+1/2}^{k+1} + E_{\zeta\xi}|_{n_{\xi}, n_{\eta}, n_{\zeta}+1/2}^{k}}{2} \quad (11) \\ & = \frac{H_{\eta}|_{n_{\xi}+1/2, n_{\eta}, n_{\zeta}+1/2}^{k+1/2} - H_{\eta}|_{n_{\xi}-1/2, n_{\eta}, n_{\zeta}+1/2}^{k+1/2}}{\Delta\xi}. \end{aligned} $$

In equation (11), the k superscript denotes the time evolution and the triple n denotes the volume discrete coordinates. Here, the triple index does not follow the traditional
x, y, z scheme but the permutation (ξ, η, ζ) introduced in the previous section.

Expression (11) is a discrete form of equation (7) and guarantees second-order precision for a uniform mesh. After some manipulation, one expresses the electric field sub-component at the time step k + 1, as a function of other field components at preceding time steps. This is the typical finite difference scheme used for classical FDTD.

An extension for the MRTD yields from equation (7):

$$ 
\begin{aligned}
& \varepsilon \frac{E_{\zeta\xi}|_{n_{\xi},n_{\eta},n_{\zeta}+1/2}^{k+1} - E_{\zeta\xi}|_{n_{\xi},n_{\eta},n_{\zeta}+1/2}^{k}}{\Delta t} \\
& + \varepsilon\alpha_{\xi}(n_{\xi}) \frac{E_{\zeta\xi}|_{n_{\xi},n_{\eta},n_{\zeta}+1/2}^{k+1} + E_{\zeta\xi}|_{n_{\xi},n_{\eta},n_{\zeta}+1/2}^{k}}{2} \quad (12) \\
& = \frac{\sum_{m_{\xi}=-\infty}^{\infty} a_{m_{\xi}-n_{\xi}+1} H_{\eta}|_{m_{\xi}+1/2,n_{\eta},n_{\zeta}+1/2}^{k+1/2}}{\Delta\xi},
\end{aligned}
$$

where the a coefficients are computed as detailed in [5]. Above, summations must be truncated by using coefficient threshold criteria.

One would also like to use variable (structured) meshing. As mentioned in [8], convolution coefficients must be corrected in order to ensure the convergence. An alternative approach is to replace the convolution coefficients by:

$$ b_{n_{\xi},m_{\xi}}^{\xi} = \frac{a_{m_{\xi}-n_{\xi}+1}}{\sum_{\nu=-\infty}^{\infty} a_{\nu-n_{\xi}+1} \xi_{\nu-n_{\xi}+1/2}}. \quad (13) $$

Some problems appear with equation (13), because the partial derivative estimation is of first-order accuracy [8]. As a result, the error increases. However, coefficient correction can be implemented to reach convergence. For example, one can set:

$$ b_{n_{\xi},n_{\xi}-2}^{\xi,\text{CORR}} = b_{n_{\xi},n_{\xi}-2}^{\xi} - \frac{\lambda_{n_{\xi}}}{\xi_{n_{\xi}-1/2} - \xi_{n_{\xi}-3/2}}, \quad (14a) $$

$$ b_{n_{\xi},n_{\xi}-1}^{\xi,\text{CORR}} = b_{n_{\xi},n_{\xi}-1}^{\xi} + \frac{\lambda_{n_{\xi}}}{\xi_{n_{\xi}-1/2} - \xi_{n_{\xi}-3/2}}, \quad (14b) $$

$$ b_{n_{\xi},n_{\xi}}^{\xi,\text{CORR}} = b_{n_{\xi},n_{\xi}}^{\xi} - \frac{\lambda_{n_{\xi}}}{\xi_{n_{\xi}+3/2} - \xi_{n_{\xi}+1/2}}, \quad (14c) $$

$$ b_{n_{\xi},n_{\xi}+1}^{\xi,\text{CORR}} = b_{n_{\xi},n_{\xi}+1}^{\xi} + \frac{\lambda_{n_{\xi}}}{\xi_{n_{\xi}+3/2} - \xi_{n_{\xi}+1/2}}, \quad (14d) $$

$$ b_{n_{\xi},m_{\xi}}^{\xi,\text{CORR}} = b_{n_{\xi},m_{\xi}}^{\xi} |m_{\xi} - n_{\xi} + 1/2| > 3/2, \quad (14e) $$

$$ \lambda_{n_{\xi}} = \frac{\sum_{\nu=-\infty}^{\infty} b_{n_{\xi},\nu}^{\xi}(\xi_{\nu+1/2} - \xi_{n_{\xi}})^2}{(\xi_{n_{\xi}+3/2} + \xi_{n_{\xi}+1/2}) - (\xi_{n_{\xi}-1/2} + \xi_{n_{\xi}-3/2})}. \quad (15) $$

It is worth to note that one can use FDTD a coefficients and replace equation (13) by those of MRTD. The advantage is that one can implement a mixed algorithm by combining Battle-Lemarie [5] and FDTD [9] schemes, as proposed in [8]. This method, called HMRTD, operates on parallelepipedic cells (irregular meshing) with the following strategy: for directions along cell sides which are electrically small, the FDTD scheme applies, while for directions along larger cell sides, MRTD scheme is used instead. The hybrid cell combines the advantage brought by

the fast FDTD algorithm and the low dispersive characteristic of MRTD. Indeed, the latter allows one to use cells with sides as large as 0.4λmin, which is beyond the limit for which FDTD model is at cutoff.

The system in equation (14) is a correction of the convolution coefficients of equation (13). It ensures a second-order accuracy. Following the Taylor polynomial decomposition, the first-order error must be proportional to the second derivative of H in equation (12). The scheme computes approximately this derivative and then removes the first-order error. With the corrections in equation (14), equation (12) becomes:

$$ 
\begin{aligned}
& \varepsilon \frac{E_{\zeta\xi}|_{n_{\xi},n_{\eta},n_{\zeta}+1/2}^{k+1} - E_{\zeta\xi}|_{n_{\xi},n_{\eta},n_{\zeta}+1/2}^{k}}{\Delta t} \\
& + \varepsilon\alpha_{\xi}(n_{\xi}) \frac{E_{\zeta\xi}|_{n_{\xi},n_{\eta},n_{\zeta}+1/2}^{k+1} + E_{\zeta\xi}|_{n_{\xi},n_{\eta},n_{\zeta}+1/2}^{k}}{2} \quad (16) \\
& = \sum_{m_{\xi}=-\infty}^{\infty} b_{n_{\xi},m_{\xi}}^{\xi,\text{CORR}} H_{\eta}|_{m_{\xi}+1/2,n_{\eta},n_{\zeta}+1/2}^{k+1/2}.
\end{aligned}
$$

The updating procedure for the other subcomponents is straightforward. It is only necessary to interchange half-integer indexes by the integer ones, when magnetic field components are updated.

## 4 Comparison of the new model with other methods

Two examples will be discussed. The first one is a finite conductivity waveguide where the performance of the PML is compared to that of Teixeira's work [6]. To the best of our knowledge, it is the most recent paper showing the results for lossy PML. The second example is a NRD waveguide [11] with high level of losses and the results are compared with analytical solution. The NRD technology introduced by Yoneyama [12] is an application for frequency operation around 50 GHz. Modes propagating into the dielectric slab sandwiched between two parallel metallic plates do not radiate laterally and are confined within the slab, as shown in Figure 2.

In the first example, we implemented a finite conductivity waveguide of εr = 12 and σ = 0.2 S/m. The number of cells in the whole volume was 3920 (20 × 2 × 98) for both cases and the number of iterations was 10000. As shown in Figure 3, the comparison between Teixeira's algorithm and the proposed one indicates a better absorption from 20 to 35 dB over the frequency band, with computer time reduced by one half.

In the last experiment, the level of losses was increased with respect to the one used in [7]. In addition, the Longitudinal Section Magnetic (LSM) LSM10 was used instead of the Longitudinal Section Electric (LSE) LSE00. The reason is that the latter produces more attenuation. Therefore, the signal reaching the PML interface would have a low level that would prevent an accurate evaluation of the absorption. Since the problem has an analytical solution [12], results may be trustily compared. Figure 4
<figure><img src="image_3.png" /><figcaption>Fig. 2. (Color online) NRD waveguide [9] under test.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Reflection coefficient for both PMLs. Teixeria [6] with CPU time of 611.68 s and HMRTD [8] with 308.69 s.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. (Color online) Mesh of the NRD structure. FDTD update on small cell and MRTD scheme on large ones.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Reflection coefficient for Figure 2.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Transmission coefficient for Figure 2. Comparison with analytical model.</figcaption></figure>

shows a top view of the NRD guide meshing for which irregular cells were used to demonstrate the efficiency of the HMRTD. Results in Figure 5 show that reflection coefficient is below -70 dB over the whole frequency range. In Figure 6, one can observe that differences with the analytical model are within 0.02 dB only. Note, that dispersion behavior (frequency dependence of the conductivity) was not implemented in any of the models.

## 5 Conclusion

An extension of the PML formalism to include insertion of lossy dielectric materials was presented for HMRTD irregular meshing. The new PML algorithm was shown to yield better performance in terms of computer expenditure as compared to the existing methods developed for
FDTD computation. The algorithm was then tested with lossy NRD structure and very good performances were obtained when compared to an analytic model, even with relatively high conductivity.

Current work concerns the implementation of dispersion behavior of the lossy dielectric.

We thank the PhD program Colciencias, National Doctoral Program, the University of Los Andes, the group of Electronic and Telecommunications Systems of the University of Los Andes, the Research Center of Engineering Faculty CIFI of the University of Los Andes and the Lab-STICC of the TELECOM Bretagne, TELECOM Institute, France.

## References

1. J.P. Berenger, J. Comput. Phys. **114**, 185 (1994)

2. M. Ozkar, G. Lazzi, A. Morttazawi, IEEE Microwave Wireless Compon. Lett. **13**, 220 (2003)

3. Q. Cao, Y. Chen, R. Mittra, in *IEEE Ant. Propag. Soc. Int. Symp. 3* (Orlando, FL, USA, 1999), p. 2010

4. E.M. Tentzeris, R. Robertson, M. Krumpholz, L.P.B. Katehi, in *IEEE Ant. Propag. Soc. Int. Symp. (Baltimore, MD, USA, 1996)*, p. 634

5. M. Krumpholz, L. Katehi, IEEE Trans. Microwave Theory Tech. **44**, 555 (1996)

6. F. Teixeira, C. Moss, W. Chew, J. Kong, IEEE Trans. Microwave Theory Tech. **50**, 30 (2002)

7. I. Massy, N. Peña, M. Ney, in *NUMELEC* (Liège, Belgium, 2008)

8. I. Massy, N. Peña, M. Ney, Int. J. Numer. Model. **23**, 470 (2010)

9. K.S. Yee, IEEE Trans. Ant. Propag. **14**, 302 (1966)

10. A. Tafllove, S.C. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain Method*, 2nd edn. (Artech House, Norwood, MA, 2000)

11. J.A. Monsoriu, B. Gimeno, IEEE Trans. Microwave Theory Tech. **52**, 1693 (2004)

12. T. Yoneyama, IEEE Trans. Microwave Theory Tech. **29**, 1188 (1981)