# Implementation of dispersion models in the split-field-finite-difference-time-domain algorithm for the study of metallic periodic structures at oblique incidence

A. Belkhir, O. Arar, S.S. Benabbes, O. Lamrous, F.I. Baida

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

## To cite this version:

A. Belkhir, O. Arar, S.S. Benabbes, O. Lamrous, F.I. Baida. Implementation of dispersion models in the split-field-finite-difference-time-domain algorithm for the study of metallic periodic structures at oblique incidence. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics [2001-2015], 2010, 81, pp.046705. 10.1103/PhysRevE.81.046705. hal-00479461

HAL Id: hal-00479461

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

Submitted on 9 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_2.png" /></figure>
# Implementation of dispersion models in the split-field-finite-difference-time-domain algorithm for the study of metallic periodic structures at oblique incidence

A. Belkhir,¹ O. Arar,² S. S. Benabbes,² O. Lamrous,¹ and F. I. Baida²,*

¹Laboratoire de Physique et Chimie Quantique, Université Mouloud Mammeri, Tizi-Ouzou, Algeria

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

(Received 8 October 2009; revised manuscript received 22 February 2010; published 29 April 2010)

We extend here the finite-difference-time-domain (FDTD) algorithm working in oblique incidence to dispersive media. The split-field method (SFM) is used and adapted for taking into account the metal dispersion. The additional equations to the FDTD algorithm are given. Instead of the 24 field components usually used in the SFM, 38 and 112 field components are needed to implement the cases of Drude and Drude-Lorentz dispersion models, respectively. Some tests are presented to validate our code as long as the angle of incidence is lower than 76° in addition to an application dealing with enhanced transmission and showing original results.

DOI: 10.1103/PhysRevE.81.046705

PACS number(s): 02.70.-c, 95.75.Pq

## I. INTRODUCTION

The study of periodic structures (photonic crystals) or samples exhibiting collective effects (enhanced transmission) is still a challenge for theoreticians working in the domain of nano-optic even if a lot of methods can hardly treat such a problem. On the contrary of the determination of the band diagram of a three-dimensional (3D) photonic crystal, the calculation of reflection and transmission through a finite structure is still difficult especially if this later is composed of metallic parts that present large dispersion in the considered electromagnetic domain. Moreover, this difficulty believes when the symmetry of the problem is destroyed such as the off-plane propagation inside a photonic crystal when this last is illuminated at oblique incidence. So, it does not exist only identified and rigorous method that is worth its used for this large panel of calculations. In these circumstances, the finite-difference-time-domain (FDTD) method seems a natural and competitive tool considering its large domain of application [1]. Despite its popularity, it is necessary to operate a development of the FDTD method that follows the scientific evolution and allowing an extension of its domain of application. In what concerns us here, this development involves with the oblique incidence problem in the case of light interaction with dispersive metallic periodic structures.

More precisely, and after recently developing the split-field method (SFM) in the FDTD algorithm [2] and use it to make original study dealing with the enhanced transmission [3], it was necessary to integrate dispersion models to describe metal dispersion in the visible range. Let us note that other techniques adapting the FDTD to the oblique incidence case exist (see Ref. [4]). However, our aim here is to extend the implementation of the SFM with integration of the metal dispersion considering the well-known Drude and Drude-Lorentz models. Thus this metal dispersion is implemented with such a technique. Although, the Valuev team [5] recently

proposed and demonstrated an iterative technique that can include dispersion models used in conventional FDTD algorithm. However, this technique needs the consideration of more than one spatial period of the structure and seems to be very time consuming because of the large number of iterations necessary for its convergence. Let us mention that a comparison between the SFM and the iterative method is done in Ref. [4].

In Ref. [4] a very brief description of the method which includes only the Debye dispersion model with using the auxiliary differential equation (ADE) method is presented. However, this paper is devoted to the implementation of the Drude and Drude-Lorentz models in the context of the SFM-FDTD by using the ADE and recursive convolution methods, respectively.

## II. THEORETICAL DEVELOPMENTS

Let us begin by introducing the SFM technique: the FDTD algorithm, in its basic version, is unable to treat the oblique incidence problem in the case of a periodic structure. In fact, for such periodic structure, only one period must be considered in addition to the periodic conditions that enable to process the infinite structure. When analytically expressed, these periodic conditions let appear a term that depends explicitly on the frequency of the wave. Because of the temporal character of the FDTD, this term cannot be directly incorporated into the algorithm. Some solutions were developed in order to solve this problem [1]. One of them, called the SFM, was recently implemented by us [2] in the case of nondispersive materials such as perfect conductors or pure dielectrics. The SFM technique was chosen according to its ability to be easily extended to the case of dispersive media.

As it is well known, in the visible range, real metals are no more perfect conductors and they have frequency-dependent dielectric properties. To be taken into account by the FDTD method, the dispersion must be analytically described. Drude model [6] is one of the most used and efficient way to describe the dispersion of some noble metals (aluminum and silver, for example) in the visible range. For

<footnote>*Corresponding author. FAX: 00 (33) 3 81 66 64 16; fbaida@univ-fcomte.fr</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Scheme of a 2D periodic structure with the geometrical parameters needed for the description of the incident field and for the expansion of the periodic conditions.</figcaption></figure>

other noble metal, such as gold, the Drude model is not efficient because of the contribution of the bound electrons to the dielectric function in the visible range. So, a Drude-Lorentz model is considered.

In order to give a clear insight of our algorithm, we should explicitly give some equations.

First, let us consider a plane wave impinging on the bidimensional [two-dimensional (2D)] periodic structure at oblique incidence (see Fig. 1). The electromagnetic field of this plane wave can be expressed as

$$ \vec{E}(x,y,z,t) = \vec{E}_0(t)e^{j(k_x x + k_y y + k_z z)}, \quad (1) $$

$$ \vec{H}(x,y,z,t) = \vec{H}_0(t)e^{j(k_x x + k_y y + k_z z)}, \quad (2) $$

where $k_x = \frac{\omega}{v_i} \sin \theta \cos \phi$, $k_y = \frac{\omega}{v_i} \sin \theta \sin \phi$, and $k_z = \frac{\omega}{v_i} \cos \theta$ are the three Cartesian components of the incident wave vector and $v_i$ is the light speed in the medium of incidence.

As mentioned above, only one period of the structure is considered for the FDTD calculations. The periodic conditions can then be written as

$$ \vec{E}(x + p_x, y, z, t) = \vec{E}(x, y, z, t)e^{jk_x p_x}, \quad (3a) $$

$$ \vec{E}(x, y + p_y, z, t) = \vec{E}(x, y, z, t)e^{jk_y p_y}, \quad (3b) $$

$$ \vec{H}(x, y, z, t) = \vec{H}(x + p_x, y, z, t)e^{-jk_x p_x}, \quad (3c) $$

$$ \vec{H}(x, y, z, t) = \vec{H}(x, y + p_y, z, t)e^{-jk_y p_y}. \quad (3d) $$

As described in Ref. [2], a set of variables ($\vec{P}, \vec{Q}$), also called the global variables, are necessary to implement these last equations;

$$ \vec{P} = \vec{E}e^{-j(k_x x + k_y y)}, \quad (4a) $$

$$ \vec{Q} = \vec{H}e^{-j(k_x x + k_y y)}. \quad (4b) $$

According to this transformation, additional variables [$\vec{P}_{ia}, \vec{Q}_{ka}$ with $(i, k) = x, y$ or $z$], named intermediate variables, are also introduced and the Maxwell equations lead to two

separable systems of equations. The first system allows the determination of the intermediate variables versus the global ones. It is given by

$$ \varepsilon \frac{\partial P_{xa}}{\partial t} = \frac{\partial Q_z}{\partial y} - \frac{\partial Q_y}{\partial z}, \quad (5) $$

$$ \varepsilon \frac{\partial P_{ya}}{\partial t} = \frac{\partial Q_x}{\partial z} - \frac{\partial Q_y}{\partial x}, \quad (6) $$

$$ \varepsilon \frac{\partial P_{za}}{\partial t} = \frac{\partial Q_y}{\partial x} - \frac{\partial Q_x}{\partial y}, \quad (7) $$

$$ \mu \frac{\partial Q_{xa}}{\partial t} = \frac{\partial P_y}{\partial z} - \frac{\partial P_z}{\partial y}, \quad (8) $$

$$ \mu \frac{\partial Q_{ya}}{\partial t} = \frac{\partial P_z}{\partial x} - \frac{\partial P_x}{\partial z}, \quad (9) $$

$$ \mu \frac{\partial Q_{za}}{\partial t} = \frac{\partial P_x}{\partial y} - \frac{\partial P_y}{\partial x}. \quad (10) $$

In the same way, the second system grants the global variables versus the intermediate ones:

$$ P_z = \left(\frac{\varepsilon}{\varepsilon - \beta}\right) \left[ P_{za} + \frac{\sin(\theta)\cos(\phi)}{\varepsilon v_i} Q_{ya} - \frac{\sin(\theta)\sin(\phi)}{\varepsilon v_i} Q_{xa} \right], \quad (11) $$

$$ Q_z = \left(\frac{\varepsilon}{\varepsilon - \beta}\right) \left[ Q_{za} + \frac{\sin(\theta)\sin(\phi)}{\mu v_i} P_{xa} - \frac{\sin(\theta)\cos(\phi)}{\mu v_i} P_{ya} \right], \quad (12) $$

$$ P_x = P_{xa} + \frac{\sin(\theta)\sin(\phi)}{\varepsilon v_i} Q_z, \quad (13) $$

$$ P_y = P_{ya} - \frac{\sin(\theta)\cos(\phi)}{\varepsilon v_i} Q_z, \quad (14) $$

$$ Q_x = Q_{xa} - \frac{\sin(\theta)\sin(\phi)}{\mu v_i} P_z, \quad (15) $$

$$ Q_y = Q_{ya} + \frac{\sin(\theta)\cos(\phi)}{\mu v_i} P_z, \quad (16) $$

with $\beta = \frac{\sin^2 \theta}{\mu v_i^2}$.

The explicit dependence of these last equations on the permittivity $\varepsilon$ requires a specific treatment in the case of dispersive material. One can verify through Eqs. (11)–(16) that, at normal incidence, the intermediate variables are equal to the global ones and consequently, the second system of equations disappears.

### A. Drude model adaptation

Because of the dispersion [$\varepsilon = \varepsilon(\omega)$], these two systems [Eqs. (5)–(16)] cannot be directly incorporated in the FDTD
algorithm. The resolution needs a new set of variables that can be expressed similarly to the constitutive equation connecting the electric displacement field ($\tilde{D}$) to the electric field ($\tilde{E}$). To illustrate this relation, let us consider the x component of the intermediate variables $P_{xa}$ for which the introduced variable is $L_{xa} = \varepsilon(\omega)P_{xa}$. By injecting this equality into Eq. (5), one obtains

$$ \frac{\partial L_{xa}}{\partial t} = \frac{\partial Q_z}{\partial y} - \frac{\partial Q_y}{\partial z}. \qquad (17) $$

The discretization of this last equation allows the calculation of the variable $L_{xa}$ as follows:

$$ L_{xa}^{n+1}(i,j,k) = L_{xa}^{n}(i,j,k) + \frac{\Delta t}{\Delta y}[Q_{z}^{n}(i,j,k) - Q_{z}^{n}(i,j-1,k)] \\ - \frac{\Delta t}{\Delta z}[Q_{y}^{n}(i,j,k) - Q_{y}^{n}(i,j,k-1)], \qquad (18) $$

where $\Delta x$ and $\Delta y$ are the spatial steps and $\Delta t$ is the temporal one.

By injecting the dielectric function given by the Drude model, we can determine the intermediate variable $P_{xa}$ versus $L_{xa}$ through:

$$ L_{xa} = \varepsilon_0 \left( 1 - \frac{\omega_p^2}{\omega^2 + i\omega\Gamma_D} \right) P_{xa}, \qquad (19) $$

where $\varepsilon_0$ is the permittivity of vacuum. Assuming time dependence in $e^{-i\omega t}$, a simple Fourier transform ($\omega \to t$) of this last equation [7] leads to

$$ \frac{\partial^2 L_{xa}}{\partial t^2} + \Gamma_D \frac{\partial L_{xa}}{\partial t} = \varepsilon_0 \left( \frac{\partial^2 P_{xa}}{\partial t^2} + \Gamma_D \frac{\partial P_{xa}}{\partial t} + \omega_p^2 P_{xa} \right). \qquad (20) $$

The partial derivatives of this equation are then replaced by their expression through the centered finite difference schema. A discretized equation is then obtained:

$$ \xi P_{xa}^{n+1} = -\chi P_{xa}^{n-1} + 4\varepsilon_0 P_{xa}^{n} + L_{xa}^{n+1}[\Gamma_D \Delta t + 2] - 4L_{xa}^{n} \\ + [-\Gamma_D \Delta t + 2]L_{xa}^{n-1}, \qquad (21) $$

where $\xi = \varepsilon_0[\omega_p^2\Delta t^2 + \Gamma_D\Delta t + 2]$ and $\chi = \varepsilon_0[\omega_p^2\Delta t^2 - \Gamma_D\Delta t + 2]$. This last equation connects the $P_{xa}$ component at the moment $(n+1)\Delta t$ to the same component at the moments $n\Delta t$ and $(n-1)\Delta t$ in addition to the new component $L_{xa}$ that is given by Eq. (18) at the two moments $n\Delta t$ and $(n-1)\Delta t$. The same strategy of calculation is then used to determine the two other components $P_{ya}$ and $P_{za}$.

The next step consists of the calculation of the global variables $\tilde{P}$ and $\tilde{Q}$. For this purpose, we have to inject the dielectric function given by the Drude model into Eqs. (11)–(16). After some algebraic operations we obtain the following four equations:

$$ \beta P_z + \frac{\sin \theta \cos \phi}{v_i} Q_{ya} - \frac{\sin \theta \sin \phi}{v_i} Q_{xa} = \varepsilon_0 [1 - \omega_p^2/(\omega^2 + i\omega\Gamma_D)] (P_z - P_{za}), \qquad (22) $$

$$ \beta Q_z = \varepsilon_0 [1 - \omega_p^2/(\omega^2 + i\omega\Gamma_D)] (Q_z - Q_{za} - \frac{\sin \theta \sin \phi}{\mu v_i} P_{xa} + \frac{\sin \theta \cos \phi}{\mu v_i} P_{ya}), \qquad (23) $$

$$ \varepsilon_0(P_x - P_{xa})[1 - \omega_p^2/(\omega^2 + i\omega\Gamma_D)] = \frac{\sin \theta \sin \phi}{v_i} Q_z, \qquad (24) $$

$$ \varepsilon_0(P_y - P_{ya})[1 - \omega_p^2/(\omega^2 + i\omega\Gamma_D)] = - \frac{\sin \theta \cos \phi}{v_i} Q_z. \qquad (25) $$

To illustrate the rest of the calculations, we present the step to be followed for the calculation of the $P_z$ component from Eq. (22). We define

$$ L_z = \varepsilon(\omega)P_z = \varepsilon_0[1 - \omega_p^2/(\omega^2 + i\omega\Gamma_D)]P_z. \qquad (26) $$

According to Eq. (22), we have

$$ L_z = L_{za} + \beta P_z + \frac{\sin \theta \cos \phi}{v_i} Q_{ya} - \frac{\sin \theta \sin \phi}{v_i} Q_{xa}, \qquad (27) $$

where $L_{za} = \varepsilon(\omega)P_{za}$, is one of the components calculated above and used to determine $P_{za}$.

Once $L_z$ is calculated, we determine $P_z$ through Eq. (26). Using an inverse Fourier transform to the time domain followed by a centered finite difference discretization of the partial derivatives with respect to time, we get

$$ \xi P_z^{n+1} = -\chi P_z^{n-1} + 4\varepsilon_0 P_z^{n} + L_z^{n+1}[\Gamma_D \Delta t + 2] - 4L_z^{n} \\ + [-\Gamma_D \Delta t + 2]L_z^{n-1}. \qquad (28) $$

Replacing $L_z^{n+1}$ given in expression (27) into Eq. (28), we finally obtain the equation for updating $P_z$ at time $(n+1)$:

$$ \eta P_z^{n+1} = -\chi P_z^{n-1} + 4\varepsilon_0 P_z^{n} + (\Gamma_D \Delta t + 2) \left[ L_{za}^{n+1} + \frac{\sin \theta \cos \phi}{v_i} Q_{ya}^{n+1} - \frac{\sin \theta \sin \phi}{v_i} Q_{xa}^{n+1} \right] \\ - 4L_z^{n} + (-\Gamma_D \Delta t + 2)L_z^{n-1}, \qquad (29) $$

where $\eta = [\xi - (\Gamma_D \Delta t + 2)\beta]$.

The other components, namely, $Q_z$, $P_x$, and $P_y$ are obtained in a similar way.

## B. Drude-Lorentz model adaptation

The dielectric function given by the Drude-Lorentz model, which corrects the Drude model for some noble metals in the visible range, is written as follows [8]:

$$ \varepsilon(\omega) = \varepsilon_{\infty} - \frac{\omega_p^2}{\omega^2 + i\omega\Gamma_D} - \frac{\Delta\varepsilon\Omega_L^2}{\omega^2 - \Omega_L^2 + i\Gamma_L\omega}. \qquad (30) $$

The components which need a particular treatment in the dispersive media described by this model are the same as in the Drude model case.
The determination of the intermediate variables $P_{xa}$, $P_{ya}$, and $P_{za}$ is done through equations that are analog to the ones used for normal incidence case. For example, the calculation of the $P_{xa}$ component needs replacing the Eq. (5) by the set of the three following equations:

$$ \psi_{Dxa}^{n} = C_{D}^{p} \psi_{Dxa}^{n-1} + C_{D}^{\delta} P_{xa}^{n}, \quad (31a) $$

$$ \psi_{Lax}^{n} = C_{L}^{p} \psi_{Lxa}^{n-1} + C_{L}^{\delta} P_{xa}^{n}, \quad (31b) $$

$$ P_{xa}^{n+1} = C^{\alpha} P_{xa}^{n} + C^{\beta} \left( \frac{\partial Q_{z}^{n+1/2}}{\partial y} - \frac{\partial Q_{y}^{n+1/2}}{\partial z} \right) + C^{\gamma} \mathrm{Re}(\psi_{Dxa}^{n} + \psi_{Lxa}^{n}), \quad (31c) $$

where the constants $C^{\alpha}$, $C^{\beta}$, $C^{\gamma}$, $C_{D}^{p}$, $C_{D}^{\delta}$, $C_{L}^{p}$, and $C_{L}^{\delta}$ are defined in [8]. The $\psi_{Dxa}$ and $\psi_{Lxa}$ are recursive accumulators functions.

The calculation of $P_{ya}$ and $P_{za}$ is carried out exactly in the same manner. The total field components $\vec{P}$ and $\vec{Q}$ requires the introduction of other variables leading to new equations. To illustrate the manner in which these components are calculated, let us consider the $P_z$ component given by Eq. (11). By introducing $M_z=P_z-P_{za}$ into this later, we obtain

$$ \varepsilon M_z = \beta P_z + \frac{\sin \theta \cos \phi}{v_i} Q_{ya} - \frac{\sin \theta \sin \phi}{v_i} Q_{xa}. \quad (32) $$

To use the same technique used with $P_{xa}$, let us transform this Eq. (32) in such a way it looks like Eq. (5). For this purpose, we differentiate the Eq. (32) with respect to time and obtain

$$ \varepsilon \frac{\partial M_z}{\partial t} = \beta \frac{\partial P_z}{\partial t} + \frac{\sin \theta \cos \phi}{v_i} \frac{\partial Q_{ya}}{\partial t} - \frac{\sin \theta \sin \phi}{v_i} \frac{\partial Q_{xa}}{\partial t}. \quad (33) $$

Similarly to the scheme used above, Eq. (33) is replaced by

$$ \psi_{Dz}^{n} = C_{D}^{p} \psi_{Dza}^{n-1} + C_{D}^{\delta} M_{z}^{n}, \quad (34a) $$

$$ \psi_{Lz}^{n} = C_{L}^{p} \psi_{Lz}^{n-1} + C_{L}^{\delta} M_{z}^{n}, \quad (34b) $$

$$ M_{z}^{n+1} = C^{\alpha} M_{z}^{n} + C^{\beta} \left( \beta \frac{\partial P_{z}}{\partial t} + \frac{\sin \theta \cos \phi}{v_{i}} \frac{\partial Q_{ya}}{\partial t} - \frac{\sin \theta \sin \phi}{v_{i}} \frac{\partial Q_{xa}}{\partial t} \right) + C^{\gamma} \mathrm{Re}(\psi_{Dz}^{n} + \psi_{Lz}^{n}). \quad (34c) $$

Time derivatives of the Eq. (34c) are then transformed into centered finite differences which leads to

$$ (1 - C^{\beta}/\Delta t \beta) P_{z}^{n+1} = P_{za}^{n+1} + C^{\alpha} M_{z}^{n} + C^{\beta}/\Delta t \left( -\beta P_{z}^{n} \right. \\ \left. \qquad + \frac{\sin \theta \cos \phi}{v_{i}} (Q_{ya}^{n+1} - Q_{ya}^{n}) \right. \\ \left. \qquad - \frac{\sin \theta \sin \phi}{v_{i}} (Q_{xa}^{n+1} - Q_{xa}^{n}) \right) \\ + C^{\gamma} \mathrm{Re}(\psi_{Dz}^{n} + \psi_{Lz}^{n}). \quad (35) $$

Finally, in addition to Eqs. (34a), (34b), and (35), a fourth equation is necessary to the determination of the $P_z$ component:

$$ M_{z}^{n+1} = C^{\alpha} M_{z}^{n} + C^{\beta}/\Delta t \left( \beta(P_{z}^{n+1} - P_{z}^{n}) + \frac{\sin \theta \cos \phi}{v_{i}}(Q_{ya}^{n+1} - Q_{ya}^{n}) \right. \\ \left. - Q_{ya}^{n} \right) - \frac{\sin \theta \sin \phi}{v_{i}}(Q_{xa}^{n+1} - Q_{xa}^{n}) \\ + C^{\gamma} \mathrm{Re}(\psi_{Dz}^{n} + \psi_{Lz}^{n}), \quad (36) $$

where $P_x$, $P_y$, and $Q_z$ are then obtained by following a similar procedure.

As in Ref. [2], the PML technique is also adapted and used to solve the boundary conditions in the z direction.

Finally, we remind that the temporal step $\Delta t$ must satisfy the stability criterion given by

$$ \Delta t \le \frac{\min(\Delta x, \Delta y, \Delta z)(v_i^2 \mu Y - \sin^2 \theta)}{v_i[|\sin \theta \cos \phi| + |\sin \theta \sin \phi| + \sqrt{3v_i^2 \mu Y - 2 \sin^2 \theta(1 - |\sin \theta \cos \phi|)}]}, \quad (37) $$

where $Y=\min(\varepsilon, \varepsilon_{\infty})$ with $\varepsilon$ as the minimal dielectric permittivity of the nondispersive media in the calculation domain and $\varepsilon_{\infty}$ as the smallest permittivity of the considered dispersive media at infinite frequencies. $\mu$ is the smallest magnetic permeability of these different materials and $v_i$ is the phase velocity of the wave in the incident medium.

According to this stability criterion, one can see that the SFM-FDTD algorithm is limited to angles of incidence less than a virtual value $\Theta$ given by $\sin \Theta = v_i \sqrt{\mu Y}$. Note that this

limit disappears in the case where the incident media is vacuum and $\varepsilon_{\infty}=\varepsilon_0$.

### III. DEMONSTRATION OF THE ACCURACY OF OUR CODE

We choose to present both 2D and full 3D calculations in order to validate the developed algorithm. In addition, and in order to check the accuracy of our code according to inci-
<figure><img src="image_4.png" /></figure>

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

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

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

<figure><img src="image_8.png" /><figcaption>FIG. 2. (Color online) Transmission coefficients through 30-nm-thick silver [(b) and (c)] and gold [(e) and (d)] films suspended in air (a). The circles correspond to the FDTD calculations ($T_{\text{FDTD}}$), the solid curves are analytical data ($T_{\text{ana}}$), and the dashed curves correspond to the error signals. The error signals are calculated as $|T_{\text{FDTD}} - T_{\text{ana}}|/T_{\text{ana}}$. The wavelength is fixed to $\lambda=650$ nm. The silver and gold dispersions are described by the Drude and Drude-Lorentz models, respectively. The spatial step is set to 1 nm.</figcaption></figure>

dence angle, we will start presenting one-dimensional (1D) tests with Drude and Drude-Lorentz models. A quantitative comparison is established between our SFM-FDTD results and those obtained analytically. This is shown in Fig. 2. Let us mention that the wavelength is fixed to $\lambda=650$ nm and the incidence angle $\theta$ varies from 0° to 76°. For incidence angles greater than 76°, the temporal step $\Delta t$ given by the stability criterion of Eq. (37) becomes very small requiring considerable computing time.

As seen in Fig. 2, the relative error, which remains lower than 10%, is maximal for angles close to 76° in TE polarization for both cases (gold and silver). In all cases the calculated relative error for $\theta \le 60°$ remains lower than 1%. This

last value can be decreased by increasing the PML layer number and/or by decreasing the spatial step. Since the most experimental procedures in nano-optics frequently occurs with incidence angle $\theta < 60°$, we can conclude that our algorithm is sufficiently accurate for simulating such experiences.

The 2D test deals with enhanced transmission through a 2D metallic grating made in silver. The considered structure [see Fig. 3(a)] was theoretically and experimentally studied by Blaikie et al. in Ref. [9]. The structure is composed of three layers deposited on a glass substrate (n=1.512). The first one is a 35-nm-thick silver layer followed by 200 nm of a dielectric layer (n=1.454). The third one is the grating
<figure><img src="image_9.png" /><figcaption>FIG. 3. (Color online) Results corresponding to the 2D test: (a) schematic of the studied structure with all the geometrical parameters. (b) and (c) show the transmission spectra calculated by both SFM-FDTD and FMM methods for $\theta=0^\circ$ and $\theta=5^\circ$, respectively.</figcaption></figure>

itself which is composed of slits engraved into a silver film of 35 nm of thickness. The grating period is set to 485 nm and the slit width is of 220 nm. In order to compare our results with the ones published in Ref. [9], we first determine the Drude-Lorentz model that matches faithfully the experi-

mental data of Ref. [10] used in Ref. [9]. Thus, using a genetic algorithm we find the constants of Eq. (30): $\epsilon_\infty = 3.600$, $\omega_p = 1.406 \times 10^{16}$ rad/s, $\Gamma_D = 2.794 \times 10^{13}$ rad/s, $\Delta\epsilon = 0.663$, $\Omega_L = 1.632 \times 10^{16}$ rad/s, and $\Gamma_L = 1.827 \times 10^{16}$ rad/s. The obtained Drude-Lorentz model is then

<figure><img src="image_10.png" /><figcaption>FIG. 4. (Color online) Comparison between the FMM and the SFM-FDTD in the case of a TM polarized incident plane wave. (a) Scheme of the square annular aperture array with all the geometrical parameters. (b) The transmission spectra versus the angle of incidence obtained with the FMM. (c) The same spectra calculated with the SFM-FDTD developed algorithm. (d) Comparison between two specific spectra at $\theta=20^\circ$. In all simulations the azimuthal angle $\phi$ is fixed to $0^\circ$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIG. 5. (Color online) (a) Transmission spectra in the case of an annular aperture arrays made in silver or gold films. For both the tow metals, the inner and outer radii of each aperture are set to 50 and 100 nm, respectively, and the structure is supposed to be self-suspended in air and has 240 nm thickness. [(b)–(d)] The three cylindrical components of the electric field [$E_r(x,y)$ in (b), $E_\phi(x,y)$ in (c), and $E_z(x,y)$ in (d)] inside one aperture for the case of silver film illuminated at $\theta=45^\circ$ in TM polarization. All amplitudes are normalized to the maximum of the radial component $E_r$. As shown in (c) and (d), $E_\phi$ and $E_z$ are very small compared to $E_r$ but they do not completely vanish because the metal is not perfectly conductor. As in Fig. 4, the azimuthal angle $\phi$ is set to zero.</figcaption></figure>

used in both homemade Fourier modal method (FMM) and SFM-FDTD codes to determine the transmission through the structure when it is illuminated by a TM-polarized plane wave (the slits are perpendicular to the electric field of the incident wave). The plane wave impinges on the structure at $\theta=0^\circ$ or $\theta=5^\circ$ from the substrate side.

The obtained results are presented in Fig. 3 where a comparison is done between our own results obtained by SFM-FDTD and FMM codes. A good agreement is found when comparing these results to the one of Figs. 3(b) (dashed red curve) and Fig. 4(a) (full blue curve) of Ref. [9].

A full 3D calculations in the case of a square annular apertures is performed in order to compare the FDTD results with the ones published in Ref. [11] obtained with FMM. The structure consists on a 100-nm-thick silver film deposited on a glass substrate and engraved with square annular grooves according to the parameters of Fig. 4(a). The dispersion of silver is described by a Drude model with $\omega_p = 1.374 \times 10^{16}$ rad/s and $\gamma = 3.21 \times 10^{13}$ rad/s. The incident plane wave illuminates the silver film from the air side. The transmission spectra versus the angle of incidence in the case of a TM polarization are depicted in Figs. 4(b) and 4(c) for both the FMM and the FDTD methods, respectively. These spectra agree qualitatively very well with the ones of Fig. 3a of Ref. [11]. A more qualitative comparison is shown in Fig.

4(d) where the same spectrum (for $\theta=20^\circ$) is calculated by the two methods. According to us, the difference between these two curves is directly linked to the spatial discretization of the sample that is performed in the FDTD method. On the contrary of the FMM, this spatial description leads to a smoothing of the structure edges. On the other hand, the convergence of the FMM greatly depends on the truncation order that was fixed here to $15 ([15 \times 2 + 1]^2 = 961$ Fourier harmonics).

A third application which is an original study is also presented: it consists on the extension of the TEM-like mode excitation [12] of an annular aperture arrays to the case of a real metal. It is well known that this mode does not have cut-off frequency and consequently, it can be of big importance for applications such as enhanced transmission. In this context, this guided mode was evoked in some papers [13–15] but its excitation was often doubtful. In fact, it was recently analytically demonstrated [3] that, to excite this peculiar guided mode, the incident plane wave must fulfill two conditions. Namely, it must be TM polarized and impinges the metallic film under oblique incidence. The excitation of such a mode was also demonstrated in the same paper [3] but only in case of a perfectly conducting metal.

In this paper, we only point out the contribution of this mode in the process of the enhanced transmission through
annular aperture arrays. A complete study including the variation in the geometrical parameters together with the metal nature is to be performed if we attempt to optimize the excitation and the propagation of this specific mode.

So, let us consider an annular aperture arrays [apertures with circular sections instead of a square ones as shown in Fig. 4(a)] engraved into a self-suspended (in vacuum) silver film. First, we set the values of the inner and the outer radius of each aperture to 50 and 100 nm, respectively, and the metal thickness to 240 nm. For this specific structure, the excitation of the TEM-like mode is clearly demonstrated on Fig. 5(a) where additional peak of transmission occur only in the TM polarization and for oblique incidence but for both silver and gold. The light distribution inside one annular aperture is plotted in Figs. 5(b)–5(d) at the TEM peak (λ=900 nm here) in order to confirm the TEM-like characteristics of this mode.

The case of a gold structure is also presented in the same figure and shows, as expected, redshift of all the transmission peaks accompanied with a transmission decreasing. This property is directly linked to the metal dispersion, i.e., to its absorption in the considered spectral domain [13]. For this simulation, the gold dispersion is obtained through a Drude-

Lorentz model [Eq. (30)] with parameters given in Ref. [8] that fit the experimental data of Ref. [16].

## IV. CONCLUSION

In summary, we demonstrate the implementation of the Drude and the Drude-Lorentz dispersion models in a SFM-FDTD algorithm working at oblique incidence. The SFM is demonstrated to be very compatible with this implementation and allows a competitive efficiency compared to recent results published in [5] based on an iterative technique. Our algorithm is very explicit and can be used for several metals over a large spectral domain. The integration of a dispersion model that includes a sum of Lorentzian terms is conceivable. This opens the way to numerically perform quasirigorous calculations involving with metallodielectric photonic crystals that are now widely used in nano-optics.

## ACKNOWLEDGMENT

This work was partly supported by the project MAGNETO-PHOT n° BLAN06-2-135594 funded by the Agence National de la Recherche (ANR).

[1] A. Taflove and S. C. Hagness, *Computational Electrodynamics: The Finite-Difference-Time-Domain Method*, 2nd ed. (Artech House, Norwood, MA, 2005).

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

[3] F. I. Baida, *Appl. Phys. B: Lasers Opt.* **89**, 145(R) (2007).

[4] F. I. Baida and A. Blekhir, *Opt. Lett.* **34**, 2453 (2009).

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

[6] The Drude model is given by $\epsilon(\omega)=1-\frac{\omega_p^2}{\omega(\omega+i\gamma)}$, where $\omega_p$ is the plasma frequency and $\gamma$ is associated to the losses.

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

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

[9] R. J. Blaikie, L. Lin, and R. J. Reeves, *Proceedings of the International Conference on Nanoscience and Nanotechnology (ICONN 2006)*, Brisbane, Australia, (IEEE Press, Piscataway,

New Jersey, 2006), pp. 458–461.

[10] P. B. Johnson and R. W. Christy, *Phys. Rev. B* **6**, 4370 (1972).

[11] D. Van Labeke, D. Gérard, B. Guizal, F. I. Baida, and L. Li, *Opt. Express* **14**, 11945 (2006).

[12] For a perfect conductor, this mode is called TEM mode ($E_\phi$ and $H_z$ are zero). In the case of a real metal, these two components of the guided electromagnetic field can not completely vanish because of the penetration of the field inside the metal.

[13] F. I. Baida and D. Van Labeke, *Opt. Commun.* **209**, 17 (2002).

[14] J. Rybczynski, K. Kempa, A. Herczynski, Y. Wang, M. J. Naughton, Z. F. Ren, Z. P. Huang, D. Cai, and M. Giersig, *Appl. Phys. Lett.* **90**, 021104 (2007).

[15] T. Thio, *Nat. Nanotechnol.* **2**, 136 (2007).

[16] E. Palik, *Handbook of Optical Constants of Solids* (Academic Press, New York, 1985).