# Split-field FDTD method for oblique incidence study of periodic dispersive metallic structures

F.I. Baida, A. Belkhir

► To cite this version:

F.I. Baida, A. Belkhir. Split-field FDTD method for oblique incidence study of periodic dispersive metallic structures. Optics Letters, 2009, 34 (16), pp.2453-2455. 10.1364/OL.34.002453. hal-00440782

HAL Id: hal-00440782

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

Submitted on 18 Apr 2021

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
# Split-field FDTD method for oblique incidence study of periodic dispersive metallic structures

F. I. Baida¹,* and A. Belkhir²

¹Département d'Optique P.M. Duffieux, Institut FEMTO-ST UMR 6174 CNRS Université de Franche-Comté, 25030 Besançon Cedex, France

²Laboratoire de Physique et Chimie Quantique, Université Mouloud Mameri, Tizi-Ouzou, Algeria

*Corresponding author: fbaida@univ-fcomte.fr

The study of periodic structures illuminated by a normally incident plane wave is a simple task that can be numerically simulated by the finite-difference time-domain (FDTD) method. On the contrary, for off-normal incidence, a widely modified algorithm must be developed in order to bypass the frequency dependence appearing in the periodic boundary conditions. After recently implementing this FDTD algorithm for pure dielectric materials, we here extend it to the study of metallic structures where dispersion can be described by analytical models. The accuracy of our code is demonstrated through comparisons with already-published results in the case of 1D and 3D structures.

The extension of the finite-difference time-domain (FDTD) algorithm to the study of infinitely biperiodic structures illuminated by plane waves at oblique incidence is a very difficult task that was partially done in the 1990s [1]. Several techniques were used to solve the implementation of the periodic boundary conditions that are explicitly frequency dependent. In fact, the FDTD algorithm, which works in the temporal domain, does not allow such a dependence. One of the most popular methods to bypass this problem is the split-field method (SFM), which was introduced by Roden et al. [2]. The implementation of this method needs the consideration of only one unit cell of the periodic structure. It is based on the well-known $\{\vec{E}, \vec{H}\} \to \{\vec{P}, \vec{Q}\}$ transformation, leading to new boundary conditions similar to those obtained in the case of normal incidence. Moreover, the discretization of the obtained equations, to calculate the total fields ($\vec{P}$ and $\vec{Q}$), is done by introducing new variables $\vec{P}_a$ and $\vec{Q}_a$ called “a-fields,” in order to conserve the stability of the Yee algorithm. We already implemented and discussed the temporal criterion of this method in [3], where dielectric or perfect conductors are considered. For these two cases, where there is no dispersion, the number of variables in the FDTD grows from six (the conventional components of the electromagnetic field) to 24 components [3]. Compared with the iterative technique recently published by Valuev et al. [4], the SFM seems to be less computer-time consuming, because it requires the discretization of only one period in addition to the fact that the convergence of the iterative technique is obtained for more than 20 iterations. On the other hand, the SFM is perfectly adapted to the integration of analytical dispersion models necessary to describe metals in the visible and IR domains. Consequently, we demonstrate here the implementation of three dispersion models i.e., Debye, Drude, and Drude-Lorentz ones, and show its accuracy by presenting results in comparison with analytical or published ones.

The mathematical demonstration of the implementation is briefly presented in the case of a generalized model that can handle both Debye dispersion and metal with constant conductivity. This is done using the auxiliary differential equation (ADE) method. Assuming a time dependence $e^{-i\omega t}$, the dielectric function that describes both these two cases can be expressed as

$$\epsilon(\omega) = \epsilon_r + \frac{\Delta\epsilon}{\alpha - i\omega\tau}, \quad (1)$$

where $\alpha$ is set to 0 for metal with constant conductivity or $\alpha=1$ to find the Debye model. According to Eqs. (7)–(9) of [3], the three components of $\vec{P}$ and those of $\vec{P}_a$, in addition to $Q_z$, depend explicitly on $\epsilon(\omega)$ and cannot be directly calculated by FDTD algorithm. Their resolution needs a new set of variables that can be expressed similarly to the constitutive equation connecting the electric displacement vector to the

<figure><img src="image_2.png" /><figcaption>Fig. 1. (Color online) Comparison of the transmission of the monolayer structure depicted in the inset when illuminated at 45° calculated by the three methods mentioned in the legend. The transmission is defined as the total flux of the transmitted Poynting vector divided by the incident one. In the SFM-FDTD simulation, the spatial step is fixed to the period over 80.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Fig. 2. (Color online) Transmission coefficients through 20-nm-thick silver [(a)–(d)] and gold [(e)–(h)] films suspended in air. In all eight subfigures, the circles correspond to the FDTD calculation ($T_{\text{FDTD}}$), the solid curves are analytical data ($T_{\text{Ana}}$), and the dotted curves correspond to the error signal. (a), (b) and (e), (f) correspond to TE polarization, while (c), (d) and (g), (h) correspond to TM one. In (a), (c) and (e), (g) the wavelength is fixed to $\lambda=700$ nm. The error signal is calculated as $|T_{\text{FDTD}}-T_{\text{Ana}}|/T_{\text{Ana}}$. In (b), (d) and (f), (h) the incidence angle is set to $\theta=45^{\circ}$. The parameters of the Drude model are $\omega_p=1.374 \times 10^{16}$ rad/s and $\gamma=3.21 \times 10^{13}$ rad/s. For the Drude–Lorentz model we consider the parameters given in Table I of [6]. In all simulations, the spatial step is set to 0.5 nm.</figcaption></figure>

electric field. For example, to calculate the components of $\vec{P}_a$, we define $\vec{L}_a = \epsilon_0 \epsilon(\omega) \vec{P}_a$. By injecting this new relation into the modified Maxwell equations we connect this new vector $\vec{L}_a$ to the $\vec{Q}$ one by

$$ \frac{\partial \vec{L}_a}{\partial t} = \vec{\nabla} \wedge \vec{Q}. \qquad (2) $$

After discretization of the differential Eq. (2) by means of centered finite differences, we first calculate $\vec{L}_a$, and second we determine the $\vec{P}_a$ components, as in the case of conventional FDTD algorithm [5], using the constitutive equation. After Fourier transform of this last equation, we get the temporal relation

$$ \epsilon_0 \epsilon_r \tau \frac{\partial \vec{P}_a}{\partial t} + \epsilon_0 (\epsilon_r \alpha + \Delta \epsilon) \vec{P}_a = \tau \frac{\partial \vec{L}_a}{\partial t} + \alpha \vec{L}_a. \qquad (3) $$

The same strategy of calculation is then used to determine the other components of $\vec{P}$ and $\vec{Q}_z$.

The accuracy of this code is first demonstrated through the study of the same 3D structure considered in [4]. This later consists of a square lattice photonic crystal monolayer of metallic spheres, where the metal is considered with constant conductivity ($\sigma=2$) and $\epsilon_r=4$. In our case, this metal dispersion is taken into account by fixing $\alpha=0$ and $\Delta\epsilon/\tau=\sigma/\epsilon_0$ in the Debye dielectric function [Eq. (1)]. The sphere radius is set to $0.375a$, where $a$ is the period of the structure. To have an accurate description of the spheres, no subgriding technique is used, but a spatial step of $a/80$ is considered. As presented in Fig. 1, the obtained spectrum presents a very good agreement with the layer Korringa–Kohn–Rostocker method, especially for the low values of frequency. Note here that, for applications in nano-optics where the ratio $a/\lambda$ is generally less than 2, the SFM-FDTD results seem to be more efficient than the ones obtained with the iterative technique (see the dashed zone in Fig. 1). In addition, the iterative technique deprives the FDTD method of its main advantage, i.e., the obtaining of a wideband spectral response with only one simulation.
<figure><img src="image_4.png" /><figcaption>Fig. 3. (Color online) Relative error deduced from analytical data (dashed curves) for a thickness uncertainty of 0.25 nm in comparison with the SFM-FDTD one (solid curves) in the case of $h=20$-nm-thick gold film. (a) Transmission angular spectrum at $\lambda=700$ nm and (b) transmission versus wavelength at $\theta=45^{\circ}$.</figcaption></figure>

The Drude model, which is generally accurate for metals such aluminum or silver in the visible and IR spectral domains, is also implemented with the same technique of ADE. In this case, the total number of variables (field components) is then of about 38 components, instead of the 12 components used to implement the same Drude model in the case of normal incidence.

For other noble metals, such gold, the Drude model is not efficient. A Drude–Lorentz model is then used and implemented with the same parameters given in Table I of [6]. In this case of dispersion, the SFM-FDTD algorithm, together with the recursive convolution method, becomes more complex and the number of variables increases to exceed 110 components instead of the 18 electromagnetic field components needed in the case of normal incidence. Note here that multiple Lorentz terms can be considered in order to faithfully describe the metal dispersion [7].

From the several tests we have done we choose here to present results demonstrating the good agreement with analytical data by considering transmission through $h=20$-nm-thick silver (for the Drude

model) or gold (for the Drude–Lorentz model) film in vacuum. Results of these comparisons are presented in Fig. 2 where both spectral and angular (versus the angle of incidence) responses are studied. These simulations were performed for both TE and TM polarizations, while the case of conical incidence (case of an arbitrary plane of incidence) is also implemented.

To estimate the error induced by the spatial meshing, which is generally linked to the geometrical description of the structure, we have compared our SFM-FDTD relative error to the analytical one induced by an uncertainty of 0.25 nm on the metallic film thickness (the order of convergence of the standard FDTD is supposed to be linear with the spatial step [8]). So, two analytical calculations were done for $h=19.75$ nm and $h=20.25$ nm to obtain the relative analytical errors presented in Fig. 3 in comparison with the FDTD ones. The good agreement between the results presented in Figs. 2 and 3, confirmed by a maximum relative error less than almost $10^{-2}$, demonstrates the efficiency of the developed code and opens up the path for us to the modeling of several recent experiments dealing with metallodielectric periodic structures, especially in the domain of nano-optics.

In summary, the developed SFM-FDTD code is a powerful tool allying at once the FDTD advantages (especially the obtention of wideband spectral response through only one calculation) to the rather faithful description of metal dispersion. Its main drawback is the manipulation of a very large number of variables and, sometimes, the need of a posttreatment of the recorded data to extract the real electromagnetic signal. Practically, this algorithm is restricted to incidence angles less than 75°, which is not a real limitation compared with what is performed in most experiences.

## References

1. A. Taflove and S. C. Hagness, *Computational Electrodynamics, the Finite-Difference Time-Domain Method*, 2nd ed. (Artech House, 2005).

2. J. A. Roden, S. D. Gedney, M. P. Kesler, J. G. Maloney, and P. H. Harms, IEEE Trans. Microwave Theory Tech. **46**, 420 (1998).

3. A. Belkhir and F. I. Baida, Phys. Rev. E **77**, 056701 (2008).

4. I. Valuev, A. Deinega, and S. Belousov, Opt. Lett. **33**, 1491 (2008).

5. R. M. Joseph, S. C. Hagness, and A. Taflove, Opt. Lett. **16**, 1412 (1991).

6. A. Vial, A.-S. Grimault, D. Macías, D. Barchiesi, and M.-L. de la Chapelle, Phys. Rev. B **71**, 085416 (2005).

7. T. O. Körner and W. Fichtner, Opt. Lett. **22**, 1586 (1997).

8. A. Farjadpour, D. Roundy, A. Rodriguez, M. Ibanescu, P. Bermel, J. D. Joannopoulos, S. G. Johnson, and G. W. Burr, Opt. Lett. **31**, 2972 (2006).