A modular implementation of dispersive materials for time-domain simulations with application to gold nanospheres at optical frequencies

D. Baumann; C. Fumeaux; C. Hafner; E. P. Li

doi:10.1364/OE.17.015186

1. Introduction

Present research in computational electromagnetics (CEM) is tackling complex problems, which have become solvable only in recent years thanks to widespread availability of powerful computers with ample memory. Complex problems can be typically defined as those that require prohibitively high computer resources or demand modelling capabilities going beyond standard simulation tools. Those categories might include electrically large problems, multi-scale structures, or devices which involve materials with demanding electrical characteristics. Among the areas of research that challenge currently available EM simulation tools, the following can be mentioned: Biomedical applications for health monitoring and non-invasive diagnosis techniques, the blooming area of metamaterials or the full-wave analysis of optical nano-structures, including carbon nanotubes, photonic bandgap structures and surface plasmon devices [1, 2]. As different as they are, these examples share the common ground of depending on materials with strongly frequency-dependent characteristics. The modeling of such materials is a crucial element which must not be neglected in order to obtain accurate simulation results.

The choice of an appropriate CEM tool is a first important step for successful simulations. Obviously, frequency-domain (FD) methods can handle dispersive materials naturally as they solve Maxwell’s equations at fixed frequencies. However, time-domain (TD) techniques become the methods of choice if transient problems are investigated, or if nonlinear phenomena need to be included in the simulation, as it might be the case, e.g. when nano-structures produce a strong localized field enhancement. The most prominent TD methods are applied commonly in structured, hexahedral (orthogonal) meshes, and build the basis of many commercial CEM tools. While permitting very efficient and fast solvers, structured meshes might hamper the computation of complex problems for several well-known reasons. Submeshing techniques can alleviate the problem of a 3-D explosion of the number of cells in multi-scale problems, however up to date, submeshing remains a rather delicate matter. In contrast, a discretization with an unstructured, polyhedral mesh, inherently enables multi-scaling because it can be strongly inhomogeneous without the drawback of requiring hanging nodes (as is the case in structured meshes) [3–5]. The Finite-Element TD (FETD) method is a relatively common technique applied in tetrahedral meshes [6]. However, as a drawback FETD tends to be inefficient for large-scale problems, since large matrixes have to be inverted at every time step. This is why recently the Discontinuous Galerkin TD (DGTD) method has received increasing attention. The DGTD algorithm is a weak FE formulation based uniquely on local interactions between adjacent cells, which enables an explicit time-stepping scheme [7]. As an alternative to this method which is still in its infancy, the Finite-Volume TD (FVTD) method, which can be considered as low-order DGTD, is by now well investigated and matured. The explicit nature of the FVTD method facilitates the computation of very large problems due to a linear scaling of memory with the number of elements. However, compared to e.g. the Finite-Difference TD (FDTD) method [8–12], the development of dispersive materials for FVTD has not been investigated deeply. In this paper, the three most common implementations of frequency-dependent materials are adapted to FVTD, namely the recursive convolution (RC) [8], piecewise-linear recursive convolution (PLRC) [9], and the auxiliary differential equation (ADE) [10, 11] techniques. Their derivation for FVTD is accomplished in a general and modular way, where only one additional loss term is added to the standard FVTD update equation. This approach ensures that frequency-dependent materials can be easily added to existing finite-volume (FV) codes as well as to other TD methods without changing the logic in the algorithm’s main iteration loop.

The paper is organized as follows: In section 2 the derivation for the RC, PLRC and ADE models of Debye, Lorentz and Drude type frequency-dependent electric and magnetic materials are given together with a description of the necessary parameters. The implementation will be verified by means of backscattering calculations and the necessary total-field scattered-field (TF/SF) technique is described for FVTD in section 3. The analytical solution for backscattering, computed by a Mie series is briefly introduced in section 4 as it includes extended definitions of Mie coefficients for electric and magnetic materials. Several validation examples proving the high accuracy of the simulation of dispersive materials in FVTD can be found in section 5. The modular approach presented here, enables a flexible approximation of complex material characteristics. For example, a proper choice of material models (here a two-pole Lorentz/one-pole Drude combination) can accurately represent the permittivity of gold in the optical frequency range. This is exploited in a TD investigation of two gold nanospheres coupled in close proximity, as presented in section 6. Although other numeric tools can solve this problem as well, the attractive properties of FVTD are well expressed for this type of geometries, as it is not deterred by very small gap widths and the unstructured inhomogeneous mesh can easily tackle also non-canonical shapes of particles.

2. Dispersive Materials for the FVTD Method

Since the introduction of FV methods for electromagnetics at the end of the 1980s [₁₃], FVTD has demonstrated attractive features for the solution of the Maxwell’s equations for complex, real-world problems [14–17]. The FVTD method’s versatility arise from two main characteristics: On the one hand, FV can be implemented in an explicit TD scheme, and on the other hand, it is applied in an unstructured, polyhedral mesh. Generally, FV methods are used for numerical solution of hyperbolic systems in conservation-law form. Since the Maxwell’s equations are of hyperbolic nature, FV is suitable for the application to CEM. However, in order for FVTD to be applied to Maxwell’s equation, Faraday’s and Ampere’s law have to be cast into a conservative form [3], which lead to the semi-discrete form of the FVTD update equation

Λ \partial_{t} U_{i} = - \frac{1}{V_{i}} Σ_{k = 1}^{4} Ψ_{U_{k}} - σ U_{i} - L_{i}^{PML} - L_{i}^{DM} .

In this equation, U _i=(E_x,E_y,E_z,H_x,H_y,H_z)^T denotes the collocated electromagnetic field vector in unit cell i with volume V_i. A tetrahedral mesh is used in the present implementation. The flux $Ψ_{U_{k}}$ through the triangular face k of cell i is constructed from the tangential field components located in the barycenter of k. Details of the standard flux evaluation can be found in [3]. Based on a local plane wave approach, the fluxes can be split into outgoing $Ψ_{U_{k}}^{+}$ and incoming $Ψ_{U_{k}}^{-}$ contributions which commonly is exploited for several tasks, as e.g. for S-parameter extraction [18], boundary conditions or for excitation schemes, as explained in section 3. The matrix σ=diag{σ_e,σ_e,σ_e,0,0,0} is the conductivity matrix. The material matrix Λ contains the permittivity and permeability of the medium. L^PML _i corresponds to the loss vector for an unsplit perfectly matched layer (PML) formulation, which is described in detail for radial and conformal formulations in [19, 20]. The remainder of the paper focuses on the last term in Eq. (1), namely the dielectric-material (DM) loss vector $L_{i}^{DM} = {(L_{i}^{{DM}_{e}}, L_{i}^{{DM}_{m}})}^{T}$ , which consists of two three-element vectors, representing the electric and magnetic losses of the dispersive materials.

Employing a Lax-Wendroff time discretization scheme with a predictor/corrector time stepping [21], the semi-discrete FVTD update equation, Eq. (1), can be written in a compact form as

U_{i}^{n + κ} = U_{i}^{n} - {(Λ^{κ})}^{- 1} κ Δ t_{ℓ} (\frac{1}{V_{i}} Σ_{k = 1}^{4} Ψ_{U_{k}}^{n + κ - 0.5} - σ U_{i}^{n} + L_{κ, i}^{PML, n} + L_{κ, i}^{DM, n})

where n denotes the time step index. The inhomogeneity of a tetrahedral mesh can be exploited in terms of efficiency by employing a local-time stepping (LTS) scheme, where different time steps Δt_ℓare used in different domains of the mesh Ω_ℓ depending on the cell size [22]. Application of a LTS commonly provides a computational speed-up between 2 and ten, depending on the inhomogeneity of the mesh. In Eq. (2), the two possible values of the parameter κ translate into a time advancement of a half step 0.5Δt_ℓ in the predictor step and accordingly a full step 1.0Δt_ℓ in the corrector step. Obviously, a Lax-Wendroff time-stepping algorithm makes the implementation of dispersive materials more challenging compared to e.g. the standard leap-frog time iteration of the FDTD approach.

In the following, the material matrix and dielectric loss term are defined. The material matrix Λ^κ=diag{ε^κ,ε^κ,ε^κ,µ^κ,µ^κ,µ^κ} contains

ε^{κ} = ε_{0} ε_{\infty} + \frac{κ Δ t_{ℓ}}{2} σ_{e} + \frac{1}{2} Σ_{p = 1}^{p} ℜ {{\overset{˘}{β}}_{p}^{κ}}

μ^{κ} = μ_{0} μ_{\infty} + \frac{1}{2} Σ_{p = 1}^{p} ℜ {{\overset{˘}{α}}_{p}^{κ}},

with parameters β^κ_p=ℜ{β̆^κ_p} and α^κ_p=ℜ{ᾰ^κ_p} defined as the real part of complex parameters β̆^κ_p and ᾰ^κ _p. These parameters specifically depend both on the material model (Debye, Lorentz and Drude) as well as on the implementation scheme (RC, PLRC and ADE) and will be explicitly given in paragraphs 2.1 to 2.3. The second term of the permittivity equation in Eq. (3) contains the electric losses σ_e and appears due to a semi-implicit approximation of the electro-magnetic field in Eq. (2) at time step n+0.5. The loss vector L ^DM,n_κ can be generally written as

L_{κ}^{DM, n} = ℜ {l_{1} {\overset{˘}{J}}^{n} + l_{2} {\overset{˘}{J}}^{n - κ} + l_{3} U^{n} + l_{4} U^{n - κ}}

and depends on current and previous field vectors U as well as on a (yet to be defined) loss current J̆ⁿ. The update parameters l^ν=ℜ{l̆_ν} with l̆_ν=(l̆^e_ν,l̆^m_ν) (ν=1…4) will be determined in paragraphs 2.1 to 2.3 for each material model and implementation scheme. The update parameters are different for the electric field components l̆eν and the magnetic field components l̆^m_ν, and may be of complex nature. The complex loss current J̆^ν conveniently can be computed iteratively as

{\overset{˘}{J}}^{n + κ} = {\overset{˘}{j}}_{1} {\overset{˘}{J}}^{n} + {\overset{˘}{j}}_{2} J^{n - κ} + {\overset{˘}{j}}_{3} U^{n + κ} + {\overset{˘}{j}}_{4} U^{n} .

Again, the complex-valued parameters j̆_ν=(j̆^e_ν,j̆_m_ν) (ν=1…4) will be determined in paragraphs 2.1 to 2.3. The emphasis of the FVTD approach presented here is placed on a modular implementation of dispersive materials with the aim of achieving the highest possible flexibility for generating and investigating lossy and dispersive materials. For example, in order to investigate surface plasmons of gold nanospheres, a mixed double-pole Lorentz/single-pole Drude model is able to approximate with reasonable accuracy the measured behavior of gold’s permittivity at optical frequencies. Therefore, in order to describe the properties of dispersive materials in terms of permittivity and permeability, the following general approach is chosen to define the material parameters

{\overset{˘}{ε}}_{r} = ε_{\infty} + {\overset{˘}{χ}}_{e}^{De} + {\overset{˘}{χ}}_{e}^{Lo} + {\overset{˘}{χ}}_{e}^{Dr}

{\overset{˘}{μ}}_{r} = μ_{\infty} + {\overset{˘}{χ}}_{m}^{De} + {\overset{˘}{χ}}_{m}^{Lo} + {\overset{˘}{χ}}_{m}^{Dr},

where ε_∞and µ_∞ are the relative permittivity and permeability at infinite frequency and χ̆_e and χ̆_m are the (complex) electric and magnetic susceptibilities. In the following subsections, the electric susceptibilities χ̆_e for Debye (De), Lorentz (Lo) and Drude (Dr) materials will be specified in detail. On this basis, the update parameters β^κ_p, l̆^e_ν and j̆^e_ν for the RC, PLRC and ADE method can be derived. The magnetic susceptibilities χ̆_m and the parameters α^κ_p, l̆^m_ν and j̆^m_ν can be easily derived accordingly, and therefore will not be specified explicitly here.

In the formalism of Eq. (5) and Eq. (6) different sets of update parameters enable the implementation of the RC, PLRC and ADE methods in a common framework. The RC and PLRC methods [8, 9] are quite similar: Both methods rely on a numerical convolution of the timedependent permittivity with time-dependent electric fields. The difference arises since RC assumes a constant electric field during one time step, whereas PLRC assumes a linear dependency. Obviously, PLRC is able to approximate the time-dependency of the electromagnetic field with a higher accuracy than RC [23]. The ADE method follows a very different approach as it does not compute a convolution, but solves a differential equation for an auxiliary current that describes the dielectric losses in frequency domain. ADE subsequently exploits the differentiation theorem of the Fourier transform in order to achieve a TD formulation [10, 11]. In terms of accuracy, PLRC and ADE are comparable.

The modular implementation of the different methods as proposed in this paper allows on one hand a unified realization in FVTD and on the other hand gives the user the possibility to choose between RC, PLRC and ADE formulations. While mostly PLRC and ADE are of practical use, the RC parameters are also given for the sake of completeness. The proposed implementation is achieved in a way which preserves the standard FVTD update equation as it only requires an additional loss term. The presented approach should be easily adaptable to other time-domain methods. In the following subsections, the three material models are introduced. As the theory behind RC, PLRC and ADE methods can be found in literature [23], only the resulting update parameters for the FVTD method are presented.

2.1. Debye material

A P-pole Debye medium can be characterized by a complex susceptibility ε̆^De as a function of angular frequency w and pole relaxation time γ ^De _p

{\overset{˘}{χ}}_{e}^{De} = Σ_{p = 1}^{p} \frac{Δ ε_{p}^{De}}{1 + j ω γ_{p}^{De}}

with Δ ε ^De _p=ε ^De _s,p-ε ^De _s,p-1, where ε ^De _s,p is the static relative permittivity of the p-th pole and γ De _p is the pole relaxation time. It is important to note that according to the definition employed here, ε ^De _s,0=ε_∞. After derivation of the RC, PLRC and ADE method for a Debye medium in the formalism of Eq. (2), the necessary parameters for the permittivity equation, Eq. (3), as well as for the update equations, Eq. (5) and Eq. (6), can be determined. The resulting parameters are summarized in Table 1.

2.2. Lorentz material

The complex susceptibility ε̌^Lo of a Lorentz media is commonly defined as

{\overset{˘}{χ}}_{e}^{Lo} = Σ_{p = 1}^{p} \frac{Δ ε_{p}^{Lo} \cdot {(ω_{p}^{Lo})}^{2}}{{(ω_{p}^{Lo})}^{2} + 2 j ω γ_{p}^{Lo} - ω^{2}}

where Δε ^Lo _p=ε ^Lo _s,p-ε ^Lo _s,p-1 (with ε ^Lo _s,0=ε_∞), γ ^Lo _p is the damping coefficient and wLo p is the undamped resonant angular frequency of the pole pair. Following parameters are introduced for the sake of simplicity:

ρ_{p}^{Lo} = γ_{p}^{Lo} - j \sqrt{{(ω_{p}^{Lo})}^{2} - {(γ_{p}^{Lo})}^{2}} and δ_{p}^{Lo} = ε_{0} Δ ε_{p}^{Lo} \cdot {(ω_{p}^{Lo})}^{2} ⁄ \sqrt{{(ω_{p}^{Lo})}^{2} - {(γ_{p}^{Lo})}^{2}}

Table 2 lists the necessary update parameters required in Eq. (3), Eq. (5) and Eq. (6).

Table 1. Real update parameters for Debye materials

View Table | View all tables in this article

Table 2. Complex update parameters for Lorentz materials.

View Table | View all tables in this article

Table 3. Real update parameters for Drude materials. The coefficients l₂=l₄=j₂=0 are equal to zero for all methods RC, PLRC and ADE.

View Table | View all tables in this article

2.3. Drude material

At a macroscopic scale, the Drude model is widely used to treat electromagnetic wave interaction with metals at optical wavelengths. Generally in frequency domain, the P-pole susceptibility function can be written as

{\overset{˘}{χ}}_{e}^{Dr} = - Σ_{p = 1}^{p} \frac{{(ω_{p}^{Dr})}^{2}}{ω^{2} - j ω γ_{p}^{Dr}}

where ω ^Dr _p is the Drude pole angular frequency and γ ^Dr _p is the reciprocal of the pole relaxation time. The update parameters for Eq. (3), Eq. (5) and Eq. (6) are given in Table 3.

2.4. Note on time-dependency convention

Commonly, physicists and engineers employ different conventions to describe harmonic time dependency: While physicists broadly use exp(-iωt), engineers commonly employ exp(+ jωt). Although this is well known, it is still a source of confusion. As a consequence, the dispersive models for Debye, Lorentz and Drude materials can either exhibit a positive imaginary part (when using the physicists’ models) or a negative imaginary part of χ̆, as it is the case in Eq. (9)–Eq. (12), which stick to the engineers’ convention. Of course no matter which convention is used, the physics behind the models does not change. The crucial point is that the material model employed must describe an exponentially decaying behavior of the electro-magnetic field after a transient excitation. This constraint dictates the sign of the imaginary part of the permittivity (and permeability respectively). If this is not satisfied, an active medium is created and instabilities in the simulation are very likely to arise.

2.5. Implementation

In terms of computational costs, a tetrahedral FVTD cell in a dispersive medium requires approximately 3.5 times more memory compared to a standard FVTD cell. However, the CPU time is only marginally increased since all necessary parameters can be computed in a preprocessing step.

The experience gained from the examples shown in section 5 suggests that PLRC can handle Debye and Lorentz models of any kind with a high accuracy. However, certain sets of Drude parameters seem to require an extremely small time step in PLRC in order to avoid instabilities. In contrast, ADE does not appear to require a reduction of the time step for Drude models, but it seems to be slightly less accurate for Lorentz materials (in the implementation presented in this paper). It was therefore empirically found that in the example 5.6, a double-pole PLRC Lorentz/single-pole ADE Drude model is the best choice for a successful simulation of the dispersive behavior of gold.

In oder to obtain a reliable validation for the implementation of dispersive materials in FVTD test examples are advantageously chosen where analytical solutions are available. A prominent test setup for CEM is electromagnetic scattering from spheres, where the analytical backscattering coefficients can be obtained by Mie series. Hence, in the following sections, two numerical tools necessary for the validation are introduced. First the total-field scattered-tield (TF/SF) formulation for FVTD is summarized, followed by a generalized definition of the Mie coefficients for electric and magnetic materials.

3. Total-Field Scattered-Field (TF/SF) Formulation

In numerical computations of scattering phenomena, the scattered field has to be separated from the incident field, as their sum (the total field) is the actual product of the simulation. A straightforward approach is the TF/SF technique, which performs the computation of the total field only inside of a bounded Huygens’ source surface in the center of the computational domain, in the total-field (TF) region. The scattered-field is directly available outside of this region in the scattered-field (SF) region. The theory behind TF/SF is widely known for the FDTD method [23], however, implementation of this excitation scheme for FVTD is not well documented. In a tetrahedral FVTD mesh, the source terms are introduced through the incoming fluxes on a TF/SF source surface, which can have an arbitrary shape. More precisely, the incident source terms are added to the incoming fluxes on the inner source boundary (towards the TF domain), and subtracted from the incoming fluxes at the outer source boundary (towards the SF domain)

Ψ^{-} = Ψ_{U_{k}}^{-} \pm Ψ_{i} .

Figure 1(a) displays the relation between outgoing Ψ⁺ and incoming Ψ^- fluxes for the cells adjacent to the Huygens’ source surface. The location of the depticted boundary cells is arbitrary, in the sense that they could as well be adjacent cells. For an incident plane wave, the source flux on the source surface is defined as

Ψ_{i} = \frac{1}{2} {- [{\vec{n}}_{k} \times \frac{1}{ε} {\vec{H}}_{i} - {\vec{n}}_{k} \times ({\vec{n}}_{k} \times c {\vec{E}}_{i})], [{\vec{n}}_{k} \times \frac{1}{μ} {\vec{E}}_{i} + {\vec{n}}_{k} \times ({\vec{n}}_{k} \times {c \vec{H}}_{i})]}^{T}

where c=(εµ)^-1/2 is the velocity of light in the medium. As source terms are always added to the incoming flux, the source surface can even be placed on the computational boundary and e.g. combined with a radiation boundary condition. As an example, Fig. 1(b) depicts four snapshots of the incident and scattered field around a perfectly electric conducting (PEC) sphere (R=7.5mm) located in the center of a spherical TF region (r=27.0mm) and illuminated by a pulsed plane wave in the frequency range of 10 to 50 GHz.

Fig. 1. (a) Source fluxes at the boundary of a TF/SF region (gray). The source flux is added (inside cell) or subtracted (outside) from the original flux term. (b) Snap shots of the near field around a PEC scatterer.

Download Full Size | PDF

4. Mie Scattering

The analytical solution for the scattering parameters of a material spheres can be computed using the Mie scattering algorithm [24]. Mostly, Mie scattering coefficients for electric materials with ε_r≠1 and µ_r=1 can be found in literature, but rarely for simultaneously electric ε_r≠1 and magnetic µ_r≠1 materials [25, 26]. Although the theory of Mie scattering is well-known and investigated since the late 1950s, accurate computation of the Mie coefficients is still not straightforward since the calculation of Riccati-Bessel functions becomes problematic if their argument has a large imaginary part, which is typical for Drude materials. A very clever and efficient approach for an iterative computation of the Mie coefficients can be found in [27]. Here, this approach is borrowed and extended for magnetic materials. The resulting parameters are

a_{n} = \frac{[u r_{n} (m x) ⁄ m + n (1 - u ⁄ m^{2}) ⁄ x] ψ_{n} (x) - ψ_{n - 1} (x)}{[u r_{n} (m x) ⁄ m + n (1 - u ⁄ m^{2}) ⁄ x] ζ_{n} (x) - ζ_{n - 1} (x)}

b_{n} = \frac{ψ_{n} (x) [m r_{n} (m x) ⁄ u + n (1 - 1 ⁄ u) ⁄ x] - ψ_{n - 1} (x)}{ζ_{n} (x) [m r_{n} (m x) ⁄ u + n (1 - 1 ⁄ u) ⁄ x] - ζ_{n - 1} (x)},

where the Riccati-Bessel functions ψ_n(x)=xj_n(x) and ζ_n=xh_n(x) are derived from the Spherical Bessel functions of first j_n(x) and third kind h_n(x). The ratio of the Riccati-Bessel functions are defined as r_n(mx)=ψ _n-1(mx)/ψ_n(mx) and the ratio of the material parameters are $m = \sqrt{ε_{2} μ_{2}} ⁄ \sqrt{ε_{1} μ_{1}}$ . The material of the scattering sphere is denoted with ε2 and µ2, and the ε1 and µ1 represents the background material (commonly vacuum). The normalized backscattering radar cross section σ_n used in the following section is computed using

σ_{n} = 4 ∣ S_{1} (π) ∣ ⁄ {(kR)}^{2}

where R is the radius of the sphere and k=2π/λ the angular wave number. The function S₁(π) is given by

S_{1} (π) = - Σ_{n = 1}^{\infty} (n + \frac{1}{2}) {(- 1)}^{n} (a_{n} - b_{n}) .

The interested reader is referred to [24] for a more detailed explanation on the general theory of Mie scattering.

Fig. 2. (a) Frequency characteristic of permittivity ε̆_r of a single-pole electric Debye medium resembling water. (b) Backscattering: The computed FVTD results (solid line) are compared to the analytical Mie series (line with diamond-shaped markers).

Download Full Size | PDF

5. Validation

The validation examples presented in this paper are inspired by [4, 28], where the simulated backscattering from single spheres, consisting of various dispersive materials, is compared to the analytical Mie series solution. However, here the modular implementation of dispersive materials in FVTD is put under test by considering combination of different materials featuring single and multiple poles, electric and magnetic dispersions, as well as Debye, Lorentz and Drude material models. The discretizations used for all examples range from approximately λ _rmin/15 on the surface of the scattering sphere to λ _0min/7 at the computational boundary, where λ _rmin and λ 0min are the wavelengths at the highest frequency of interest in the material and in vacuum, respectively.

5.1. Backscattering from single-pole electric Debye sphere

The electric dispersive behavior of water can be conveniently modeled as a Debye material. In this example, the backscattering from a water droplet with a radius of R=420 µm is calculated at microwave frequencies. The parameters which describe the dielectric characteristic of water are ε _∞=5.9, ε ^De _s,1=80.2, and γ ^De ₁=9.5ps [29]. Figure 2(a) depicts the value of the complex permittivity in the frequency range from 10 to 50 GHz, as well as the comparison between the analytical Mie solution (line with diamond-shaped markers) and the numerical results obtained with FVTD (solid line) in Fig. 2(b). The agreement between the Mie series and FVTD is excellent over the whole frequency range.

5.2. Backscattering from single-pole electric Lorentz sphere

The material model of Lorentz media is tested through the computation of backscattering from a sphere with radius R=900 µm. The employed Lorentz parameters are ε∞=4.1, ε ^Lo _s,1=5.8, ω ^Lo ₁=40π · 10⁹ s^-1 and γ ^Lo ₁=30π · 10⁸ s^-1. The characteristic of the permittivity is shown in Fig. 3(a) and the comparison between analytical Mie scattering (line with diamond-shaped markers) and FVTD simulation results (solid line) is plotted in Fig. 3(b). Again, the agreement between the Mie series and FVTD is very good in the frequency range from 10 to 50 GHz where the Lorentz pole shows its resonance.

5.3. Backscattering from single-pole electric Drude sphere

The third example provides the verification of the Drude material implementation. Here the backscattering from a large plasma sphere (R=3000 µm) is computed. The parameters of the Drude model are ε∞=1.0, ω ^Dr ₁=80π · 10⁹ s^-1 and γ ^Dr ₁=1.5 · 10¹⁰ s^-1. This example is especially demanding since the radius of the sphere is very large in terms of wavelength for the observed frequency range of 1 to 100 GHz. A discretization corresponding to approximately λ_{r min}/10 at 100 GHz is deployed on the surface of the sphere. Figure 4(a) plots the frequency dependent behavior of the permittivity of the Drude medium. The normalized radar cross section σ_n is shown in Fig. 4(b), where a very good agreement between the analytical Mie scattering series (line with diamond-shaped markers) and the simulated FVTD results (solid line) can be observed. The small discrepancies observed at the highest frequencies are explained by the relatively coarse discretization.

Fig. 3. (a) Frequency characteristic of permittivity $\overset{˘}{ε}$ _r of a single-pole electric Lorentz medium. (b) Backscattering: The computed FVTD results (solid line) are compared the analytical Mie series results (line with diamond-shaped markers).

Download Full Size | PDF

Fig. 4. (a) Frequency characteristic of permittivity of a single-pole electric Drude medium. (b) Backscattering: The computed FVTD results (solid line) are compared to the analytical Mie series (line with diamond-shaped markers).

Download Full Size | PDF

5.4. Backscattering from double-pole electric Debye sphere

After the verification of single-pole media of all three types (Debye, Lorentz and Drude), the general adaptability of the FVTD approach is demonstrated by the example of a double-pole material. In order to emphasize the practical relevance of such a material, a double-pole Debye model is chosen, which is commonly found to describe the frequency-dependent dielectric characteristic of human tissue. This is a relevant feature of a time-domain CEM tool, e.g. for microwave imaging devices for breast-cancer detection. The parameters of the two Debye poles representing human muscle tissue are ε∞=11.05, ε ^De _s,1=83.97, ε ^De _s,2=43.35, τ ^De ₁=8.56·10-¹² s, τ ^De ₂=2.33 ·10^-10 s. These values have been obtained by a least-square fitting with a Cole-Cole model for muscle tissue in the licensed frequency range of ultra-wideband (UWB) signals (3.1 to 10.6 GHz) [30]. Figure 5(a) depicts the frequency-dependent characteristic of the permittivity of human muscle tissue, and in (b) the comparison between FVTD results (solid line) and analytical Mie series (line with diamond-shaped markers) of the normalized RCS from a sphere composed of human muscle tissue (R=1800 µm).

Fig. 5. (a) Frequency characteristic of permittivity of a double-pole electric Debye medium resembling human muscle tissue. (b) Backscattering: The FVTD results (solid line) are compared to the analytical Mie series (line with diamond-shaped markers).

Download Full Size | PDF

Fig. 6. (a) Frequency characteristic of permittivity ε̆_r and permeability µ̆_r of a mixed Lorentz medium. (b) Backscattering: The computed FVTD results (solid line) are compared to the analytical Mie series (line with diamond-shaped markers).

Download Full Size | PDF

5.5. Backscattering from mixed single-pole electric and single-pole magnetic Lorentz sphere

As another example, the backscattering from a mixed single-pole electric and single-pole magnetic Lorentz sphere (R=900µm) is calculated. The frequency characteristic of permittivity and permeability of the Lorentz medium are depicted in Fig. 6(a). The electric Lorentz pole is determined by ε∞=4.1, ε ^Lo _s,1=8.0, ω ^Lo _e,1=40π ·10⁹ s^-1, and γ ^Lo _e,1=30π ·10⁹ s^-1. The magnetic pole is characterized by µ∞=4.1, µ ^Lo _s,1=5.0, ω ^Lo _m,1=60π · 10⁹ s^-1 and γ ^Lo _m,1=40π · 10⁹ s^-1. The normalized radar cross section σ_n is plotted in Fig. 6(b) and a good agreement between the FVTD simulation (solid line) results and the analytical Mie solution (line with diamond-shaped markers) is found. This example nicely demonstrates the flexibility of the FVTD implementation of dispersive materials for simultaneous deployment of electric and magnetic dispersive poles in a single medium.

5.6. Backscattering from a gold nanosphere

As a final validation example, the backscattering from a gold nanosphere (R=40nm) is calculated, illustrating the relevance of the technique for optical nano-structure simulations. The permittivity of gold is approximated by an electric mixed double-pole Lorentz/single-pole Drude model, which parameters have been obtained by a least-square fitting of experimental data [31]. The material model is characterized by ε∞=3.65, two Lorentz poles by ε ^Lo _s,1=7.01, ε ^Lo _s,2=5.54, ω ^Lo ₁=4.79 · 10¹⁵ s^-1, ω ^Lo ₂=6.45 · 10¹⁵ s^-1, γ ^Lo ₁=9.08 · 10¹⁴ s^-1, γ ^Lo ₂=1.39 · 10¹⁵ s^-1; the Drude pole is modeled with ω ^Dr ₁=1.28 · 10¹⁶ s^-1 and γ ^Dr ₁=2.76 · 10¹³ s^-1. The wavelength-dependent characteristic of the fitted permittivity values is plotted in Fig. 7(a) in the visible spectrum of wavelengths, where plasma oscillations occur. The measured permittivity characteristic is indicated with diamonds at discrete frequencies, which suggests that the fitted parameters approximate the experimental data very well. The simulated results (solid line) of the backscattering is shown in Fig. 7(b) where they are compared to the reference solution obtained by an analytical Mie series (line with diamond-shaped markers).

Fig. 7. (a) Frequency characteristic of permittivity of gold. The measured permittivity values are indicated with diamonds. (b) Backscattering: The FVTD results (solid line) are compared to the analytical Mie series (line with diamond-shaped markers).

Download Full Size | PDF

6. Plasmonic Resonance of Two Coupled Gold Nanospheres

As experimentally shown in [1], two gold nanospheres in close proximity can produce a large field enhancement that can be exploited for nonlinear optical mixing. A similar arrangement has been simulated using FDTD in [32], employing an extended Debye model to approximate the permittivity of gold. The setup investigated here with FVTD is depicted in Fig. 8(a) and exhibits a very narrow gap (d=1nm) in between the two spheres, where a strong plasmonic resonance of the EM field can be achieved. The permittivity of gold is approximated with the same parameters as given in section 5.6, as they provide a very good agreement with experimental data. From a numerical point of view, such an arrangement becomes increasingly challenging as the gap between the spheres becomes smaller. However, an unstructured and strongly inhomogeneous mesh can tackle such multi-scale problems, as illustrated in Fig. 8(b), where the surface and volume meshes of the computational model employed here are depicted. The gap region is magnified in order to display the fine mesh in between the spheres. Overall, the model heavily exploits the inhomogeneity of a tetrahedral mesh and hence minimizes the total number of cells to about 270’000. The LTS scheme allows to apply a time steps 64 times higher in the larger free-space cells compared to the time step in the tiny cells in the gap.

A reference solution for this example can be obtained using a quasi-analytical frequency-domain method such as the Multiple-Multipole Method (MMP) [33], which is able to generate very accurate results for such a configuration. The two nanospheres are illuminated by a plane wave with polarization parallel to the sphere-sphere axial direction. The plasmonic resonance occurs in the center of the visible range of the electromagnetic spectrum:Fig. 8(c) depicts the normalized electric field in the center of the gap between the two spheres, where the solid black line represents the reference MMP solution and the dashed red line shows the FVTD results. Although the results do not match perfectly, the difference between theMMPand FVTD solutions is lower than typical experimental measurement uncertainties at nanometer scale.

Fig. 8. (a) Sketch of two gold nanospheres separated by a narrow gap. (b) Surface and volume mesh employed in the simulation (only TF domain is shown). The gap is shown in a magnified view. (c) Normalized electric field in the center of the gap. The FVTD result (dashed line) is compared to the reference MMP solution (solid line).

Download Full Size | PDF

This example proves the applicability of FVTD for accurate simulation of optical structures: The advantage of a general-purpose simulation method such as FVTD for such problems is that there is in principle no limitations to any extension toward more complex geometries. In particular, the present problem could be extended through the use of non-canonical particle shapes or the inclusion of a substrate.

7. Conclusion

This paper has presented a modular implementation of dispersive materials for the FVTD method. The chosen scheme is flexible and general, allowing to choose between RC, PLRC and ADE methods for the modeling of P-pole electric and magnetic Debye, Lorentz and Drude materials, which can be combined in any fashion. This results in a very flexible treatment of dispersive materials, as even complicated frequency-dependent characteristics of permittivities and permeabilities can be approximated very well. As a challenging example, the plasmonic resonance of two gold nanospheres coupled in close proximity has been simulated successfully. Since FVTD is relying on an unstructured mesh, arbitrarily-shaped particles and structures can be simulated without restrictions and without increase in computational costs for a comparable size. This extension has demonstrated the capability of FVTD as a promising CEM tool for complex optical problems exhibiting dispersive characteristics.

References and links

1. M. Danckwerts and L. Novotny, “Optical Frequency Mixing at Coupled Gold Nanoparticles,” Phys. Rev. Lett. 98, 1–4 (2007). [CrossRef]

2. D. Pinto and S. Obayya, “Accurate Perfectly Matched Layer Finite-Volume Time-Domain Method for Photonic Bandgap Devices,” IEEE Photon. Technol. Lett. 20, 339–341 (2008). [CrossRef]

3. P. Bonnet, X. Ferrieres, B. Michielsen, P. Klotz, and J. Roumiguieres, Time Domain Electromagnetics, chap. 9 Finite-Volume Time-Domain Method, pp. 307–367 (Academic Press, 1999).

4. F. Edelvik and B. Strand, “Frequency Dispersive Materials for 3-D Hybrid Solvers in Time Domain,” IEEE Trans. Antennas Propag. 51, 1199–1205 (2003). [CrossRef]

5. P. Sewell, T. Benson, C. Christopoulos, D. Thomas, A. Vokuvic, and J. Wykes, “Transmission-Line Modeling (TLM) Based Upon Unstructured Tetrahedral Meshes,” IEEE Trans. Microwave Theory Tech. 53, 1919–1928 (2005). [CrossRef]

6. J.-F. Lee, R. Lee, and A. Cangellaris, “Time-Domain Finite-Element Methods,” IEEE Trans. Antennas Propag. 45, 430–442 (1997). [CrossRef]

7. J. S. Hesthaven and T. Warburton, Nodel Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, vol. 54 of Texts in Applied Mathematics (Springer New York, 2008).

8. R. Luebbers, F. P. Hunsberger, K. S. Kunz, R. B. Standler, and M. Schneider, “A Frequency-Dependent Finite-Difference Time-Domain Formulation for Dispersive Materials,” IEEE Trans. Electromagn. Compat. 32, 222–227 (1990). [CrossRef]

9. D. F. Kelley and R. J. Luebbers, “Piecewise Linear Recursive Convolution for Dispersive Media using FDTD,” IEEE Trans. Antennas Propag. 44, 792–797 (1996). [CrossRef]

10. M. Okoniewski, M. Mrozowski, and M. A. Stuchly, “Simple Treatment of Multi-Term Dispersion in FDTD,” IEEE Microwave Guided Wave Lett. 7, 121–123 (1997). [CrossRef]

11. Y. Takayama and W. Klaus, “Reinterpretation of the Auxiliary Differential Equation Method for FDTD,” IEEE Microwave Wireless Comp. Lett. 12, 102–104 (2002). [CrossRef]

12. F. L. Teixeira, “Time-Domain Finite-Difference and Finite-Element Methods for Maxwell Equations in Complex Media,” IEEE Trans. Antennas Propag. 56, 2150–2166 (2008). [CrossRef]

13. V. Shankar, W. Hall, and A. Mohammadian, “A time-domain differential solver for electromagnetic scattering problems,” Proc. IEEE 77, 709–721 (1989). [CrossRef]

14. D. Baumann, C. Fumeaux, P. Leuchtmann, and R. Vahldieck, “Finite-Volume Time-Domain (FVTD) Modeling of a Broadband Double-Ridged Horn Antenna,” Int. J. Numer. Model. 17, 285–298 (2004). [CrossRef]

15. D. K. Firsov and J. LoVetri, “FVTD - Integral Equation Hybrid for Maxwell’s Equations,” Int. J. Numer. Model. 21, 29–42 (2007). [CrossRef]

16. C. Fumeaux, D. Baumann, and R. Vahldieck, “Finite-Volume Time-Domain Analysis of a Cavity-Backed Archimedean Spiral Antenna,” IEEE Trans. Antennas Propag. 54, 844–851 (2006). [CrossRef]

17. Y. Shi and C.-H. Liang, “The Finite-Volume Time-Domain Algorithm using Least Square Method in Solving Maxwell’s Equations,” J. Comput. Phys. 226, 1444–1457 (2007). [CrossRef]

18. D. Baumann, C. Fumeaux, and R. Vahldieck, “Field-Based Scattering-Matrix Extraction Scheme for the FVTD Method Exploiting a Flux-Splitting Algorithm,” IEEE Trans. Microwave Theory Tech. 53, 3595–3605 (2005). [CrossRef]

19. C. Fumeaux, K. Sankaran, and R. Vahldieck, “Spherical Perfectly Matched Absorber for Finite-Volume 3-D Domain Truncation,” IEEE Trans. Microwave Theory Tech. 55, 2773–2781 (2007). [CrossRef]

20. D. Baumann, C. Fumeaux, R. Vahldieck, and E. P. Li, “Conformal Perfectly Matched Absorber for Finite-Volume Simulations,” in 19th Intern. Zurich Symposium on EMC (Asia Pacific EMC Week) (Singapore, 2008).

21. R. F. Warming and R. M. Beam, “Upwind second-order difference schemes and applications in aerodynamic flows,” AIAA J. 14, 1241–1249 (1976). [CrossRef]

22. C. Fumeaux, D. Baumann, and R. Vahldieck, “A generalized local time-step scheme for efficient FVTD simulations in strongly inhomogeneous meshes,” IEEE Trans. Microwave Theory Tech. 52, 1067–1076 (2004). [CrossRef]

23. A. Taflove and S. C. Hagness, Computational Electrodynamics, The Finite-Difference Time-Domain Method, 3rd ed. (Artech House, 2005).

24. H. C. Van de Hulst, Light Scattering by Small Particles (John Wiley & Sons, Inc., 1957).

25. R. F. Harrington, Time-Harmonic Electromagnetic Fields (John Wiley & Sons, Inc., 2001). [CrossRef]

26. M. Kerker, W. D.S., and C. Giles, “Electromagnetic Scattering by Magnetic Spheres,” J. Opt. Soc. Am. 73, 765–767 (1983). [CrossRef]

27. H. Du, “Mie-Scattering Calculation,” Appl. Opt. 43, 1951–1956 (2004). [CrossRef] [PubMed]

28. R. Luebbers, D. Steich, and K. Kunz, “FDTD Calculation of Scattering from Frequency-Dependent Materials,” IEEE Trans. Antennas Propag. 41, 1249–1257 (1993). [CrossRef]

29. F. S. Barnes and B. Greenebaum, eds., Biological and Medical Aspects of Electromagnetic Fields, 3rd ed. (CRC Press, 2007).

30. S. Gabriel, R. W. Lau, and C. Gabriel, “The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues,” Phys. Med. Biol. 41, 2271–2293 (1996). [CrossRef] [PubMed]

31. P. B. Johnson and R. W. Christy, “Optical Constants of the Noble Metals,” Phys. Rev. B 6, 4370–4379 (1972). [CrossRef]

32. A. Dhawan, S. J. Norton, M. D. Gerhold, and T. Vo-Dinh, “Comparison of FDTD numerical computations and analytical multipole expansion method for plasmonics-active nanosphere dimers,” Opt. Express 17, 9688–9703 (2009). [CrossRef] [PubMed]

33. C. Hafner, Generalized Multipole Techniques for Electromagnetic and Light Scattering, chap. 3, pp. 21–38 (Elsevier Science B.V., 1999).

A modular implementation of dispersive materials for time-domain simulations with application to gold nanospheres at optical frequencies

Abstract

2. Dispersive Materials for the FVTD Method

2.1. Debye material

2.2. Lorentz material

2.3. Drude material

2.4. Note on time-dependency convention

2.5. Implementation

3. Total-Field Scattered-Field (TF/SF) Formulation

4. Mie Scattering

5. Validation

5.1. Backscattering from single-pole electric Debye sphere

5.2. Backscattering from single-pole electric Lorentz sphere

5.3. Backscattering from single-pole electric Drude sphere

5.4. Backscattering from double-pole electric Debye sphere

5.5. Backscattering from mixed single-pole electric and single-pole magnetic Lorentz sphere

5.6. Backscattering from a gold nanosphere

6. Plasmonic Resonance of Two Coupled Gold Nanospheres

7. Conclusion

References and links

Cited By

Figures (8)

Tables (3)

Equations (19)

Optics Express