Time-domain discrete-dipole approximation for simulation of temporal response of plasmonic nanoparticles

Kwang-Hyon Kim; Maxim A. Yurkin

doi:10.1364/OE.23.015555

1. Introduction

Over the past decade, nonlinear plasmonic nanostructures have attracted great interest in the community of nanophotonics [1,2]. Such systems include metal nanostructures with ultrafast nonlinearity [3,4] and spaser-based nanolasers [5], loss-compensated metamaterials [6], and surface plasmon polariton amplifiers [7,8] in which active material embedded in the nanostructure exhibits saturation behavior. One needs to account for the temporal evolution of the system in response to external field for the study of these problems. Frequency-domain methods, such as the discrete dipole approximation [9,10], frequency-domain finite element method (FEM) [11], boundary integral method [12] and many others can account for the temporal response only after extensive calculation of the responses of the systems over the whole frequency range under consideration (see e. g [13,14].). Such circumstance becomes more limiting when one needs to account for the nonlinearity of nanostructure taking part in the interactions with light whose intensity depends on time. For solving these problems, time-domain methods are evidently superior to the frequency-domain methods. As a popular time-domain method, the finite-difference time-domain (FDTD) method [15] is widely applied for numerical simulations of temporal responses of plasmonic nanostructures.

The discrete-dipole approximation (DDA) has the advantages of small memory requirement for calculation and its simplicity compared to the other methods (see e. g [13].) such as FDTD, FEM and others, leading to wide applications for the investigation of electromagnetic responses of nanostructures. Nevertheless, DDA cannot directly be applied for studying time-evolution problems, since it is a frequency-domain method. Recently reported DDA in time-domain [13,14,16] performs numerical calculation in frequency-domain, and the temporal response of the system is calculated by the inverse Fourier or Laplace transform. By using such indirect method, dynamically accounting for the temporal evolution of nonlinear or intensity-dependent response is impossible.

In this paper we develop a direct time-domain DDA (TDDDA) which enables evaluating the temporal evolution of plasmonic nanostructures. The main equation for the TDDDA is obtained by taking the inverse Fourier transform of the Taylor expansion of the frequency-domain DDA in terms of frequency deviation from the central frequency. The method can be considered an intermediate between the FDTD and the frequency-domain DDA. In particuar, it solves for the oscillations amplitudes of dipole moments but does it through the time-domain differential equation. Due to the latter it requires an analytic frequency dependence of electric permittivity. We obtain the differential equation for dipole moments from its frequency-domain expression for complete description of temporal evolution of local field. We compare the simulation results for silver nanoparticles with that of the frequency-domain DDA and of the FDTD to demonstrate the correctness of the proposed TDDDA. Moreover, we discuss the numerical efficiency of the latter and its possible extension to simulation of nonlinear processes.

2. Time-domain DDA

The idea for obtaining the equation describing the temporal evolution of local field is as follows. First, we obtain a modified DDA equation in terms of frequency deviation from central frequency. The central frequency can be, for example, of a pulse incident onto the medium containing metal nanoparticles. Taking the inverse Fourier transform, we can get the TDDDA equation containing the time derivatives.

We begin with DDA equation describing distribution of local field (also denoted as exciting field) in inhomogeneus media [9,10,17]:

Ε_{m} = Ε_{m}^{i n} + \sum_{n \neq m} G_{m n} p_{n},

where E_m and

E_{m}^{in}

are the discretized local and incident fields, respectively, p_n = α_nE_n is the dipole moment of n-th dipole, and α_n = (3v/4π)(ε_n − ε_h)/(ε_n + 2ε_h) is the polarizability, where v is the volume of the dipole, ε_n = ε_n(ω) is the dielectric function of material at the position of n-th dipole, and ε_h is the permittivity of host material. We chose the simplest formulation (Claussius–Mossotti) for polarizability [10], since we are aiming at smaller nanoparticles, for which difference between polarizability formulations is negligibly small. G_mn in Eq. (1) is given by [10]

G_{m n} = \frac{\exp (i k R_{m n})}{R_{m n}} [k^{2} (1 - {\hat{r}}_{m n} \otimes {\hat{r}}_{m n}) + \frac{i k R_{m n} - 1}{R_{m n}^{2}} (1 - 3 {\hat{r}}_{m n} \otimes {\hat{r}}_{m n})],

where R_mn = |r_m – r_n|,

{\hat{r}}_{m n} = (r_{m} - r_{n}) / R_{m n}

, and k = ω/c₀, where c₀ is the speed of light in vacuum.

Let us expand Eq. (2) into the Taylor series:

G_{m n} (k) = G_{m n} (k_{0}) + {\frac{\partial G_{m n}}{\partial k} |}_{k = k_{0}} (k - k_{0}) + \frac{1}{2} {\frac{\partial^{2} G_{m n}}{\partial k^{2}} |}_{k = k_{0}} {(k - k_{0})}^{2} + \dots,

where

k - k_{0} = Δ ω ε_{h}^{1 / 2} / c_{0}

, ∆ω = ω – ω₀, and ω₀ is the central angular frequency. We consider the incident field in Eq. (1) as a plane wave given by

E_{m}^{in} = E^{in} (ω) \exp (i k r_{m ⊥}),

where Eⁱⁿ(ω) is the amplitude and r_m_⊥ = (k∙r_m)/k is the projection of r_m onto the direction of wavevector k. The above incident field can be expanded as follows:

E_{m}^{in} = E^{in} (ω) \exp (i k_{0} r_{m ⊥}) [1 + i(k - k_{0}) r_{m ⊥} - \frac{1}{2} {(k - k_{0})}^{2} r_{m ⊥}^{2} + \dots] .

To make the numerical scheme practical we have to truncate the Taylor expansions in Eqs. (3) and (5) down to a few terms. This results in the main limitation of the proposed TDDDA – (k – k₀)R << 1, where R is the characteristic particle size, assuming that the origin of the reference frame in Eq. (5) is placed inside the particle. In other words, we assume that incident wavefronts of different frequencies accumulate relatively small phase difference when passing the particle. In particular, this approximation is reasonable for broadband response of nanoparticles much smaller than the wavelength, which is the main motivation of this paper. The applicability domain of the method in terms of R somewhat increases with the number of accounted Taylor term; still, the case of (k – k₀)R ~1 cannot be efficiently handled. Further on, we consider only the second-order approximation, under which Eq. (1) can be rewritten as

E_{m} = E^{in} (ω) \exp (i k_{0} r_{m ⊥}) (1 + i \frac{ε_{h}^{1 / 2}}{c_{0}} r_{m ⊥} Δ ω - \frac{ε_{h}}{2 c_{0}^{2}} r_{m ⊥}^{2} Δ ω^{2}) + \sum_{n \neq m} (G_{m n}^{(0)} + G_{m n}^{(1)} Δ ω + G_{m n}^{(2)} Δ ω^{2}) p_{n},

where

G_{m n}^{(0)} = {G_{m n} |}_{k = k_{0}},

G_{m n}^{(1)} = \frac{ε_{h}^{1 / 2}}{c_{0}} \frac{\exp (i k_{0} R_{m n})}{R_{m n}} [k_{0} (2 + i k_{0} R_{m n}) (1 - {\hat{r}}_{m n} \otimes {\hat{r}}_{m n}) - k_{0} (1 - 3 {\hat{r}}_{m n} \otimes {\hat{r}}_{m n})],

G_{m n}^{(2)} = \frac{ε_{h}}{2 c_{0}^{2}} \frac{\exp (i k_{0} R_{m n})}{R_{m n}} [(- k_{0}^{2} R_{m n}^{2} + 4 i k_{0} R_{m n} + 2) (1 - {\hat{r}}_{m n} \otimes {\hat{r}}_{m n}) - (i k_{0} R_{m n} + 1) (i k_{0} R_{m n} + 1)] .

To validate the applicability of the second-order approximation in a quantitative manner, we further evaluate their magnitudes. For the Green’s tesor it is sufficient to compare the magnitudes of the following scalar function which constitutes the tensor:

g = \frac{\exp (i k R)}{R} [k^{2} (δ - s) + R^{- 2} (i k R - 1) (δ - 3 s)]

with its first- and second-order approximations:

g_{1} = {g |}_{k = k_{0}} + \frac{ε_{h}^{1 / 2}}{c_{0}} \frac{\exp (i k_{0} R)}{R} [k_{0} (2 + i k_{0} R) (δ - s) - k_{0} (δ - 3 s)] Δ ω,

g_{2} = g_{1} + \frac{ε_{h}}{2 c_{0}^{2}} \frac{\exp (i k_{0} R)}{R} [(- k_{0}^{2} R^{2} + 4 i k_{0} R + 2) (δ - s) - (i k_{0} R + 1) (δ - 3 s)] Δ ω^{2} .

In Eqs. (9)–(11), δ = 1 and 0 for diagonal and nondiagonal elements, respectively. For the incident field we compare the position-dependent part in Eq. (5) f = exp(ikr_⊥) with its approximations

f_{1} = \exp (i k_{0} r_{⊥}) (1 + i \frac{ε_{h}^{1 / 2}}{с_{0}} r_{m ⊥} Δ ω),

f_{2} = \exp (i k_{0} r_{⊥}) (1 + i \frac{ε_{h}^{1 / 2}}{с_{0}} r_{m ⊥} Δ ω - \frac{ε_{h}}{2 c_{0}^{2}} r_{m ⊥}^{2} Δ ω^{2}) .

Figure 1(a) presents the comparison of absolute values of g, g₁, and g₂ as the functions of wavelength λ for R = 50 nm, δ = 1, s = 0.1, and λ₀ = 2πс₀/ω₀ = 400 nm, as an example. Figures 1(b) and 1(c) presents similar results for real and imaginary parts of f, f₁, and f₂ for r_⊥ = 50 nm. Those results confirm that the second approximations perfectly agree with the corresponding original functions (maximum relative errors are 2.5% for g and 0.4% for f) for the considered wavelength ranges and sizes of nanoparticles.

Fig. 1 Validation the second-order approximation: absolute values of g, g₁, and g₂ (a) and real (b) and imaginary (c) parts of f, f₁, and f₂ versus the wavelength for the values of parameters R = r_⊥ = 50 nm, δ = 1, s = 0.1, and λ₀ = 2πс/ω₀ = 400 nm.

Download Full Size | PDF

To obtain time-domain equation for the local field, we take the inverse Fourier transform of both sides in Eq. (5), when the frequency difference is transformed to the time derivatives: ∆ω → i∂/∂t and ∆ω² → −∂²/∂t². As a result, we obtain the DDA equation describing the temporal evolution of amplitude of local field A_m(t):

A_{m} = (1 - \frac{ε_{h}^{1 / 2}}{c_{0}} r_{m ⊥} \frac{\partial}{\partial t} + \frac{ε_{h}}{2 c_{0}^{2}} r_{m ⊥}^{2} \frac{\partial^{2}}{\partial t^{2}}) A^{in} \exp (i k_{0} r_{m ⊥}) + \sum_{n \neq m} (G_{m n}^{(0)} +i G_{m n}^{(1)} \frac{\partial}{\partial t} - G_{m n}^{(2)} \frac{\partial^{2}}{\partial t^{2}}) q_{n},

where Aⁱⁿ, A_m, and q_m are the amplitudes defined by Eⁱⁿ(t) = Aⁱⁿ(t)exp(−iω₀t), E_m(t) = A_m(t)exp(−iω₀t), and p_n(t) = q_n(t)exp(−iω₀t), respectively.

3. Temporal evolution of dipole moments

Using Eq. (14), we can directly obtain the temporal evolution of optical response of metal nanostructures, unlike the indirect method presented in [13,14]. However, the drawback of operating in time domain is the need to relate the time dependences of q_n and A_n, which requires analytic expression for ε_n(ω), analogous to the FDTD [15]. In this paper we consider the extended Drude model for the dielectric function of metal:

ε_{n} (ω) = ε_{\infty} - \frac{ω_{p}^{2}}{ω^{2} + i γ ω},

where ε_∞ is the dielectric function for infinite frequency, ω_p is the plasma frequency, and γ is the electron collision frequency in metal. In principle, more complicated expressions for ε_n(ω) can also be incorporated, but we leave it for a future research. For convenience, we write the dipole moment as p_n(ω) = p_n₀(ω) + p_n₁(ω), where

p_{n 0} (ω) = α_{0} E_{n} (ω) = \frac{3 v}{4 π} \frac{ε_{\infty} - ε_{h}}{ε_{\infty} + 2 ε_{h}} E_{n} (ω),

p_{n 1} (ω) = - \frac{3 v}{4 π} \frac{w E_{n} (ω)}{ω^{2} + i ω γ - ω_{p}^{2} / (ε_{\infty} + 2 ε_{h})},

and

w = 3 ε_{h} ω_{p}^{2} / {(ε_{\infty} + 2 ε_{h})}^{2}

. Equation (17) can be rewritten as

(ω^{2} + i ω γ - \frac{ω_{p}^{2}}{ε_{\infty} + 2 ε_{h}}) p_{n 1} (ω) = - \frac{3 v}{4 π} w E_{n} (ω) .

Taking the inverse Fourier transform of both sides in the Eq. (18), we obtain the differential equation for the amplitude q_n₁(t), defined by p_n₁(t) = q_n₁(t)exp(−iω₀t):

{\ddot{q}}_{n 1} + c {\dot{q}}_{n 1} + b q_{n 1} = a A_{n},

where a = 3wv/(4π),

b = ω_{p}^{2} / (ε_{\infty} + 2 ε_{h}) - i γ ω_{0} - ω_{0}^{2}

, and c = γ – 2iω₀.

Considering p_n(t) = p_n₀(t) + p_n₁(t) = [α₀A_n(t) + q_n₁(t)]exp(−iω₀t) and Eq. (14) we obtain the final equation for the TDDDA:

A_{m} = (1 - \frac{ε_{h}^{1 / 2}}{c_{0}} r_{m ⊥} \frac{\partial}{\partial t} + \frac{ε_{h}}{2 c_{0}^{2}} r_{m ⊥}^{2} \frac{\partial^{2}}{\partial t^{2}}) A^{in} e^{i k_{0} r_{m ⊥}} + \sum_{n \neq m} (G_{m n}^{(0)} +i G_{m n}^{(1)} \frac{\partial}{\partial t} - G_{m n}^{(2)} \frac{\partial^{2}}{\partial t^{2}}) (α_{0} A_{n} + q_{n 1}) .

Equations (19) and (20) are the central result of this paper. Solving them simultaneously, we can obtain both q_n₁(t) and A_n(t) inside the nanoparticle. If required, Eq. (20) can further be used to evaluate evolution of the electric field in the vicinity of the nanoparticle as well – in full analogy with the frequency-domain DDA [10]. However, we do not currently pursue this opportunity.

4. Difference equation for the time-domain DDA

For the numerical implementation of the TDDDA, we apply the implicit difference scheme for Eqs. (19) and (20). From Eq. (19) we obtain

q_{n 1, N} = \frac{a Δ t^{2} A_{n, N} + (2 + c Δ t) q_{n 1, N - 1} - q_{n 1, N - 2}}{1 + c Δ t + b Δ t^{2}},

where q_n_1,_N = q_n₁(N∆t). Its first and the second derivatives are given by

{\dot{q}}_{n 1, N} = \frac{q_{n 1, N} - q_{n 1, N - 1}}{Δ t} = \frac{a Δ t^{2} A_{n, N} - (b Δ t^{2} - 1) q_{n 1, N - 1} - q_{n 1, N - 2}}{Δ t (1 + c Δ t + b Δ t^{2})},

{\ddot{q}}_{n 1, N} = \frac{q_{n 1, N} - 2 q_{n 1, N - 1} + q_{n 1, N - 2}}{Δ t^{2}} = \frac{a Δ t A_{n, N} - (c + 2 b Δ t) q_{n 1, N - 1} + (c + b Δ t) q_{n 1, N - 2}}{Δ t (1 + c Δ t + b Δ t^{2})},

respectively. Substituting Eqs. (21)–(23) into Eq. (20), we obtain

\begin{array}{l} A_{m, N} - (α_{0} + \frac{a Δ t^{2}}{l}) \sum_{n \neq m} (G_{m n}^{(0)} + \frac{i}{Δ t} G_{m n}^{(1)} - \frac{1}{Δ t^{2}} G_{m n}^{(2)}) A_{n, N} = \\ = [(1 - \frac{ε_{h}^{1 / 2} r_{m ⊥}}{c_{0} Δ t} + \frac{ε_{h} r_{m ⊥}^{2}}{2 c_{0}^{2} Δ t^{2}}) A_{m, N}^{in} + (\frac{ε_{h}^{1 / 2} r_{m ⊥}}{c_{0} Δ t} - \frac{ε_{h} r_{m ⊥}^{2}}{c_{0}^{2} Δ t^{2}}) A_{m, N - 1}^{in} + \frac{ε_{h} r_{m ⊥}^{2}}{2 c_{0}^{2} Δ t^{2}} A_{m, N - 2}^{in}] \exp (i k_{0} r_{m ⊥}) \\ + \sum_{n \neq m} (G_{m n}^{(0)} \frac{(2 + c Δ t) q_{n 1, N - 1} - q_{n 1, N - 2}}{l} - i G_{m n}^{(1)} {\frac{α_{0}}{Δ t} A_{n, N - 1} + \frac{1}{l Δ t} [(b Δ t^{2} - 1) q_{n 1, N - 1} + q_{n 1, N - 2}]} \\ - G_{m n}^{(2)} {\frac{α_{0}}{Δ t^{2}} (- 2 A_{n, N - 1} + A_{n, N - 2}) + \frac{1}{l Δ t} [- (c + 2 b Δ t) q_{n 1, N - 1} + (c + b Δ t) q_{n 1, N - 2}]}), \end{array}

where l = 1 + с∆t + b∆t². Using Eqs. (24) and (21) alternatingly we can numerically calculate the evolution of local field.

The enhanced field intensity $E_{m}^{enh}$ , also denoted as total internal field [10], is given by the product in the frequency domain: $E_{m}^{enh} (ω) = f E_{m} (ω)$ , where f = 3ε_h/(ε_m + 2ε_h) is the field enhancement factor. From the extended Drude model (9) for the dielectric function of metal, we have

E_{m}^{enh} (ω) = \frac{f_{0} (ω^{2} + i ω γ)}{ω^{2} + i ω γ - ω_{p}^{2} / (ε_{\infty} + 2 ε_{h})} E_{m} (ω),

where f = 3ε_h/(ε_∞ + 2ε_h). Taking the inverse Fourier transform and through the derivation similar to the case of Eq. (19), we obtain the differential equation for the amplitude of enhanced field

A_{m}^{enh} (t)

given by

A_{m}^{enh} (t) = E_{m}^{enh} (t) \exp (- i ω_{0} t)

:

{\ddot{A}}_{m}^{enh} + c {\dot{A}}_{m}^{enh} + b A_{m}^{enh} = f_{0} ({\ddot{A}}_{m} + c {\dot{A}}_{m} - h A_{m}),

where b and c are the same as in Eq. (19) and

h = ω_{0}^{2} + i γ ω_{0}

. Applying the implicit finite-difference method, we obtain the formula for numerical calculation of

A_{m}^{enh}

:

A_{m, N}^{enh} = \frac{f_{0}}{l} (1 + c Δ t - h Δ t^{2}) A_{m, N} - \frac{2 + c Δ t}{l} (f_{0} A_{m, N - 1} - A_{m, N - 1}^{enh}) + \frac{1}{l} (f_{0} A_{m, N - 2} - A_{m, N - 2}^{enh}) .

5. Numerical results and discussions

We perform the numerical simulation for the temporal response of silver nanoparticles illuminated by a pulse with planar wavefront. The result is compared with preceding approaches: the FDTD, the Mie theory, and the frequency-domain DDA. We have used the FFT for both time- and frequency-domain DDA calculations in order to reduce the calculation time [18]. The solution of the linear system has been performed by using stabilized biconjugate gradient method. All the numerical simulations presented in the paper have been performed for cubic grid with size of ∆x = 2 nm.

Figure 2 shows the temporal evolution of enhanced field in silver nanosphere with a diameter of 70 nm surrounded by air, in comparison with results of the FDTD. The Drude parameters has been taken to be ε_∞ = 5.9809, ω_p = 14.624 fs⁻¹, and γ = 0.3333 fs⁻¹ which are obtained by fitting the dielectric function of silver [19] in the wavelength range from 330 to 500 nm. The amplitude A^enh is calculated by Eqs. (21), (24), and (27) with the time step of 0.2 fs. This time step is much larger than that for the FDTD because the TDDDA deals only with the envelope or amplitude of the field. The amplitude of incident Gaussian pulse is Aⁱⁿ(t) = e_x exp(−[(t − t₀)/τ]²), where e_x is the unit vector along the x-axis, We considered parameter values t₀ = 5 fs and τ = 1.6 fs, leading to the pulse duration of 1.88 fs. The central wavelength of incident pulse is λ₀ = 390 nm. In the FDTD calculation, the total domain is 81 × 81 × 81, PML thickness is 5 grid elements. The calculation time step is ∆t = ∆x/(2c) = 3.3 × 10⁻³ fs, which is 60 times smaller than that for the TDDDA.

Fig. 2 Temporal evolution of enhanced field in Ag nanosphere with a diameter 70 nm surrounded by air. Grid size is 2 nm. The central wavelength and the duration of incident pulse are 390 nm and 1.88 fs, respectively. Figure 2(a) presents the x-component of real part of enhanced field calculated by using TDDDA in comparison with the result of FDTD at the central position of the nanosphere. Figure 2(b) illustrates the amplitudes of enhanced field, averaged over the space occupied by the nanosphere, for TDDDA and FDTD.

Download Full Size | PDF

It is improtant to note that the DDA operates with complex-valued fields and amplitudes (including all those mentioned above), implying that the actual physical field is the real part of the computed one. By contrast, the FDTD operates directly with real-valued fields. In particular the incident field in the FDTD is given as Re(Eⁱⁿ) = Aⁱⁿ(t)cos(ω₀t). To match the two methods we used two options. First, Re(E^enh) computed by the TDDDA is exactly the physical field, so its x-component at the center of nanosphere is compared with that of the FDTD in Fig. 2(a). Second, we computed the envelope amplitude with the FDTD. For that we repeated the FDTD calculation with incident field Im(Eⁱⁿ) = Aⁱⁿ(t)sin(ω₀t), which resulted in the computed field equal to Im(E^enh). Then, we used a trivial relation

‖ A^{enh} ‖ = \sqrt{{‖ Re (E^{enh}) ‖}^{2} + {‖ Im (E^{enh}) ‖}^{2}},

where ||.|| denotes the Euclidian norm of the vector. In Fig. 2(b) we compare ||A^enh|| averaged over the volume of nanosphere between the TDDDA and the FDTD. Overall, the results in Fig. 2 show good agreement between the methods and thus reveal the validity of the proposed TDDDA.

To further validate our method in the frequency domain, we present in Fig. 3 the extinction spectra for the nanosphere shown in Fig. 2, in comparison with the result of Mie theory and the conventional frequency-domain DDA. The latter implementation has been validated by comparison with the ADDA code [20] for several frequencies. The extinction efficiency Q_ext has been determined in the TDDDA through the intermediate calculation of frequency-domain local fields. The plasmon resonance wavelength of about 385 nm is nearly the same for the three methods, however, the maximum value of Q_ext calculated with both variants of the DDA is about 10% smaller than that of the Mie theory. This difference, as well as discrepancies at larger wavelengths, are due to relatively coarse discretization of the nanosphere [21]. We stress, however, that the results of time- and frequency-domain DDA agree perfectly. The maximum relative difference is 2.5%, but it decreases rapidly when approaching the central wavelength. The latter is explained by the fact that the only source of this difference is the truncation of Taylor expansions (3) and (5). The same dependence of differences on k – k₀ is expected to apply to any quantity, including near- and internal fields. While we have not performed point-by-point comparison of the latter in the frequency-domain, it should be accurate on average, since both Q_ext (Fig. 3) and time-domain fields (Fig. 2) can be expressed in terms of the internal fields.

Fig. 3 Extinction spectra for the same nanosphere as in in Fig. 2 calculated from the data obtained by the TDDDA, in comparison with that directly computed with the Mie theory and the frequency-domain DDA.

Download Full Size | PDF

Let us further compare the numerical performance of the TDDDA with that of the FDTD and the DDA, taking the results shown in Figs. 2 and 3 as examples. The TDDDA required 9 iterations of the iterative solver on average for each time step to reach the relative error of less than 10⁻⁵. Since the total simulation domain has the size of M = 35³ and the number of time steps is 100, the total computational complexity (in flops) is the order of 100 × 9 × MlogM ≈4 × 10⁸. By contrast, the frequency-domain DDA required about 160 iterations on average for the same error threshold, and 50 different frequencies were considered, leading to the computational complexity of 50 × 160 × MlogM ≈4 × 10⁹, which is about 10 times larger than that of the TDDDA. The FDTD required calculation time about 5 times longer than that of the TDDDA, mainly due to 60-times smaller time step. The only drawback of the TDDDA in this case is relatively large memory required for the time-domain fields, which are further processed into frequency-domain quantities. However, this drawback is not relevant when only the temporal responses of nanostructures is computed. In the latter case the TDDDA is shown to be superior to the other two methods. The relative acceleration of the TDDDA in comparison with the FDTD may become more significant for simulations over longer time interval, e.g., picoseconds, since the required number of time steps in the TDDDA for accurate description of the envelope is fairly independent of its width. While the frequency-domain DDA is probably more practical for simulation of linear response to long pulses, it can hardly be applied for nonlinear simulations (see below).

As examples of application of the TDDDA to nonspherical particles, we consider silver nanorod and nanodisk in Fig. 4. The nanorod has a diameter of 40 nm and a length of 70 nm and is embedded in silica glass. The incident wave propagates along the direction perpendicular to the long axis of the rod with polarization parallel to this axis. We set λ₀ = 680 nm, the pulse duration 4.71 fs, ∆t = 0.40 fs, and permittivity of silica ε_h = 2.1229. The disk has a diameter of 70 nm and a thickness of 42 nm and is also embedded in silica. Both propagation and polarization directions of incident wave are perpendicular to the symmetry axis of the disk. In this case λ₀ = 550 nm, the pulse duration is 4.12 fs, ∆t = 0.35 fs, and ε_h = 2.1414 corresponding to this value of λ₀. In both cases we have used the extended Drude model for the dielectric function of silver with parameters ε_∞ = 4.3378, ω_p = 13.385 fs⁻¹, and γ = 0.1264 fs⁻¹, obtained by fitting the data presented in [19] in the wavelength range from 300 to 900 nm. In particular, the computed extinction spectra agree well with that of frequency-domain DDA (maximum relative differences are less than 8%, but much smaller near the central wavelength).

Fig. 4 TDDDA simulations of time- and frequency-domain responses of silver nanorod with a diameter of 40 nm and a length of 70 nm (a, b) and nanodisk with a diameter of 70 nm and a thickness of 42 nm (c, d), embedded in silica glass. The time-domain response (a, c) is given in terms of the space-averaged amplitude of enhanced field in comparison with the incident one. The frequency-domain response (c, d) is given in terms of extinction spectra in comparison with the direct DDA simulations.

Download Full Size | PDF

Being a single-pass method the TDDDA proposed in this paper has the potential for the study of temporal behaviors of plasmonic systems with nonlinear responses. One example is ultrafast nonlinear response of metal nanostructures illuminated by laser pulses shorter than around 30 fs, which cannot be solved without considering the dispersion of field enhancement factor [22]. Spasers or plasmonic nanosystems with saturable gain also require the numerical calculation of temporal evolution of the system. Even for the steady state operation, frequency-domain methods can provide accurate result only after time-consuming iterative self-consistent calculations [23,24]. The nonlinearity can be incorporated through the following equation [4,23]

\frac{d Δ ε_{n}}{d t} = - \frac{Δ ε_{n}}{τ_{ep}} + \frac{χ_{n}^{(3)}}{τ_{ee} τ_{ep}} \int_{- \infty}^{t} d t^{'} {‖ A_{n}^{enh} (t^{'}) ‖}^{2} \exp (\frac{t - t^{'}}{τ_{ee}}),

where ∆ε_n is the nonlinear increment of dielectric function (to be added to the linear part),

χ_{n}^{(3)}

is the degenerate third-order nonlinear susceptibility, τ_ee and τ_ep are the electron–electron and electron–photon response times which have the magnitudes of the order of 0.1 and 1 ps, respectively. Equation (29) is a reasonable approximation for moderate laser power, i.e. when one needs to account only for the lowest-order nonlinear susceptibility [4]. Moreover, similar nonlinear modification of dielectric function was used in the framework of the frequency-domain DDA [25]. Unfortunately, such extension incurs significant technical complexity, as the whole formulation need to be modified starting from Eq. (15). Therefore, we leave it for a future research.

Although we have considered only plasmonic nanoparticles, one can apply the TDDDA also for dielectric or semiconductor ones by using the differential equation for dipole moment obtained from the Lorentz model, describing the dispersion of dielectric function of those materials, as in the FDTD. Moreover, the TDDDA can be extented to anisotropic and/or magnetodielectric scatterers with the same additional complexity as that for the frequency-domain and previous implementations of the time-domain DDA [10,14]. In particular, the magnetodielectric case requires one to double the number of equations. The only nontrivial additional component is then the cross-term Green’s tensor (between the electric and magnetic) parts, which need to be expanded in Taylor series, similar to Eq. (4).

6. Conclusion

We have proposed the TDDDA, directly describing the temporal evolution of electric field in plasmonic nanostructures (a single-pass approach). The main equation is obtained by taking the inverse Fourier transform of the Taylor expansion of the frequency-domain DDA in terms of frequency deviation from the central frequency. The differential equation for dipole moments is obtained from its frequency-domain expression for complete description of temporal evolution of local field. Compared to the frequency-domain DDA, the TDDDA has only two limitations. First, it requires an analytic frequency dependence of permittivity, e.g. the Drude model (similar to the FDTD). Second, the incident wavefronts of different frequencies should accumulate relatively small phase difference when passing the particle. In particular, the latter condition is always satisfied for nanoparticles much smaller than the wavelength. The results obtained with the TDDDA agree well with its frequency-domain counterpart and with the FDTD, taking silver nanoparticles as examples. Moreover, the TDDDA is the fastest of the three methods for incident pulses of several-femtoseconds width. Thus, it can be effectively applied for modeling the temporal responses and optimizing the structural and material parameters of plasmonic nanostructures. Also, it can potentially be extended to simulate ultrafast nonlinear processes in nanostructures.

Acknowledgments

This work was partly supported by Russian Science Foundation (Grant No. 14-15-00155).

References and links

1. M. Kauranen and A. V. Zayats, “Nonlinear plasmonics,” Nat. Photonics 6(11), 737–748 (2012). [CrossRef]

2. J. Renger, R. Quidant, N. van Hulst, and L. Novotny, “Surface enhanced nonlinear four-wave mixing,” Phys. Rev. Lett. 104(4), 04683 (2010). [CrossRef]

3. Y. Zhu, X. Hu, Y. Fu, H. Yang, and Q. Gong, “Ultralow-power and ultrafast all-optical tunable plasmon-induced transparency in metamaterials at optical communication,” Sci. Rep. 3, 2338 (2013).

4. K.-H. Kim, U. Griebner, and J. Herrmann, “Theory of passive mode locking of solid-state lasers using metal nanocomposites as slow saturable absorbers,” Opt. Lett. 37(9), 1490–1492 (2012). [CrossRef] [PubMed]

5. M. A. Noginov, G. Zhu, A. M. Belgrave, R. Bakker, V. M. Shalaev, E. E. Narimanov, S. Stout, E. Herz, T. Suteewong, and U. Wiesner, “Demonstration of a spaser-based nanolaser,” Nature 460(7259), 1110–1112 (2009). [CrossRef] [PubMed]

6. S. Wuestner, A. Pusch, K. L. Tsakmakidis, J. M. Hamm, and O. Hess, “Overcoming losses with gain in a negative refractive index metamaterial,” Phys. Rev. Lett. 105(12), 127401 (2010). [CrossRef] [PubMed]

7. D. Y. Fedyanin, A. V. Krasavin, A. V. Arsenin, and A. V. Zayats, “Surface plasmon polariton amplification upon electrical injection in highly integrated plasmonic circuits,” Nano Lett. 12(5), 2459–2463 (2012). [CrossRef] [PubMed]

8. K.-H. Kim, A. Husakou, and J. Herrmann, “Theory of plasmonic femtosecond pulse generation by mode-locking of long-range surface plasmon polariton lasers,” Opt. Express 20(1), 462–473 (2012). [CrossRef] [PubMed]

9. B. T. Draine and P. J. Flatau, “Discrete-dipole approximation for scattering calculations,” J. Opt. Soc. Am. A 11(4), 1491–1499 (1994). [CrossRef]

10. M. A. Yurkin and A. G. Hoekstra, “The discrete dipole approximation: An overview and recent developments,” J. Quant. Spectrosc. Radiat. Transf. 106(1-3), 558–589 (2007). [CrossRef]

11. J. L. Volakis, A. Chatterjee, and L. C. Kempel, Finite Element Method for Electromagnetics (IEEE, New York, 1998).

12. U. Hohenester and J. Krenn, “Surface plasmon resonances of single and coupled metallic nanoparticles: A boundary integral method approach,” Phys. Rev. B 72(19), 195429 (2005). [CrossRef]

13. P. C. Chaumet, K. Belkebir, and A. Rahmani, “Coupled-dipole method in time domain,” Opt. Express 16(25), 20157–20165 (2008). [CrossRef] [PubMed]

14. P. C. Chaumet, K. Belkebir, and A. Rahmani, “Discrete dipole approximation for time-domain computation of optical forces on magnetodielectric scatterers,” Opt. Express 19(3), 2466–2475 (2011). [CrossRef] [PubMed]

15. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech House, Boston, 2000).

16. P. C. Chaumet, T. Zhang, A. Rahmani, B. Gralak, and K. Belkebir, “Discrete dipole approximation in time domain through the Laplace transform,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 88(6), 063303 (2013). [CrossRef] [PubMed]

17. G. H. Goedecke and S. G. O’Brien, “Scattering by irregular inhomogeneous particles via the digitized Green’s function algorithm,” Appl. Opt. 27(12), 2431–2438 (1988). [CrossRef] [PubMed]

18. J. J. Goodman, B. T. Draine, and P. J. Flatau, “Application of fast-Fourier-transform techniques to the discrete-dipole approximation,” Opt. Lett. 16(15), 1198–1200 (1991). [CrossRef] [PubMed]

19. E. D. Palik, Handbook of Optical Properties of Solids (Academic, New York, 1985).

20. M. A. Yurkin and A. G. Hoekstra, “The discrete-dipole-approximation code ADDA: capabilities and known limitations,” J. Quant. Spectrosc. Radiat. Transf. 112(13), 2234–2247 (2011). [CrossRef]

21. M. A. Yurkin, D. de Kanter, and A. G. Hoekstra, “Accuracy of the discrete dipole approximation for simulation of optical properties of gold nanoparticles,” J. Nanophoton. 4(1), 041585 (2010). [CrossRef]

22. K.-H. Kim, A. Husakou, and J. Herrmann, “Saturable absorption in composites doped with metal nanoparticles,” Opt. Express 18(21), 21918–21925 (2010). [CrossRef] [PubMed]

23. K.-H. Kim, U. Griebner, and J. Herrmann, “Theory of passive mode-locking of semiconductor disk lasers in the blue specral range by metal nanocomposites,” Opt. Express 20(15), 16174–16179 (2012). [CrossRef]

24. S.-L. Chua, Y. Chong, A. D. Stone, M. Soljacić, and J. Bravo-Abad, “Low-threshold lasing action in photonic crystal slabs enabled by Fano resonances,” Opt. Express 19(2), 1539–1562 (2011). [CrossRef] [PubMed]

25. M. S. Dhoni and W. Ji, “Extension of discrete-dipole approximation model to compute nonlinear absorption in gold nanostructures,” J. Phys. Chem. C 115(42), 20359–20366 (2011). [CrossRef]

Time-domain discrete-dipole approximation for simulation of temporal response of plasmonic nanoparticles

Abstract

1. Introduction

2. Time-domain DDA

3. Temporal evolution of dipole moments

4. Difference equation for the time-domain DDA

5. Numerical results and discussions

6. Conclusion

Acknowledgments

References and links

Cited By

Figures (4)

Equations (29)

Optics Express