## 1.

## Introduction

The fitting of optical constants of metals has given rise to an extensive literature. The optical properties of bulk materials^{1}^{,}^{2} are commonly used for simulations, but recent experimental data for nanostructured materials will be used in the near future.^{3}^{,}^{4} Indeed, the surface plasmon resonance (SPR) setups are able to determine thicknesses of multilayers^{4} and also optical properties.^{5} Therefore, a versatile and efficient method of fitting of experimental data of the optical properties is of interest.

The historical Drude and Lorentz models (DL) were used^{6}^{,}^{7} for gold in a wide domain of wavelengths. Expansion of Lorentzian terms (L4) was proposed by Hao and Nordlander^{8} and tested by comparison between Mie theory and finite time difference domain (FDTD). More recently, the description of gold permittivities by means of critical points model^{9}^{,}^{10} was proved to be efficient for the modeling of spectroscopy with FDTD methods,^{11} finite element method,^{12} and discrete dipole approximation.^{13} Fitting of the optical constants is also useful if eigenvalues of complex structures are computed^{14}^{,}^{15} or if resonances are searched in dispersion curves.^{16}^{,}^{17} For the design and the optimization of nanostructures,^{18} the accuracy of numerical results depends on the quality of the fitting of the relative permittivities.^{19}^{,}^{20}

The models for fitting the relative permittivity are functions of the angular frequency of illumination $\omega $ in a specific range of wavelengths in vacuum ${\lambda}_{0}$. Each model involves parameters considered as degrees of freedom for the fitting of handbooks experimental data from literature. The parameters are searched such that the distance between the reference data^{1}^{,}^{2} and the model is minimal. The use of these models for FDTD requires filling an additional constraint that is included in the proposed method. The novelty of this paper is in the method that consists of using a particle swarm optimization method (PSO) and then in improving the solution by using a Nelder-Mead (NM) algorithm, using the solution of the PSO algorithm as starting point.

The combinations of DL, and Drude and critical points models (DCP) are shortly described in Sec. 2 as well as the criterion that guarantees the convergence of FDTD algorithms. The principle of the proposed PSO method is given in Sec. 3. The results of the fitting of the relative permittivity of some metals are given and discussed in Sec. 4 before concluding.

## 2.

## Two Models for Fitting: the Combination of Drude and Lorentz Models, and the Critical Points Model

## 2.1.

### Combination of Drude and Lorentz Models

The combination of Drude and Lorentz models (DL) describes both the intraband (Drude model) and interband (Lorentz model) electronic transitions.^{21} It enables the imaginary part to be a decreasing function of $\omega $ on the contrary of the single Drude model. The DL model is efficient in the wavelengths range $[\mathrm{500,1000}]\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$.^{7} Nevertheless, the best choice of models of fitting as well as their physical parameters depends on the investigated domain of wavelengths. This study is devoted to the fitting in the visible domain of wavelengths, $[\mathrm{400,800}]\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$.

The limitation to a single Lorentz term limits the memory requirement for FDTD, as it increases linearly with the number of terms used for the dispersion law.^{7} The function of fit ${\u03f5}_{\mathrm{DL}}$ of the relative permittivity of metal is written as the sum of the Drude and the Lorentz models.^{6}^{,}^{10}^{,}^{22}

## (1)

$${\u03f5}_{\mathrm{DL}}(\omega )={\u03f5}_{\infty}-\underset{\text{Drude model}}{\frac{{\omega}_{D}^{2}}{\underbrace{\omega (\omega +i{\gamma}_{D})}}}-\underset{\text{Lorentz model}}{\underbrace{\frac{\mathrm{\Delta}\u03f5{\omega}_{L}^{2}}{{\omega}^{2}-{\omega}_{L}^{2}+i{\gamma}_{L}\omega}}}.$$The pure imaginary number is denoted as $i$ such as ${i}^{2}=-1$. Frictional forces proportional to the velocity of electrons with viscous damping coefficients ${\gamma}_{D}$ and ${\gamma}_{L}$ are introduced and expressed in the same units as $\omega $. The plasma frequency ${\omega}_{L}$ is associated with intraband transitions, $\mathrm{\Delta}\u03f5$ is the oscillator strength, and ${\u03f5}_{\infty}$ is the relative permittivity for high frequencies. For high frequencies, electrons cannot follow excitation, and the illumination does not see free electrons anymore; therefore ${\u03f5}_{\infty}=1$.^{23} Nevertheless, it is commonly admitted that this asymptotic value could also be adjusted for fitting of the relative permittivity in a limited wavelength bandwidth. It is therefore an additional degree of freedom of the fitting. The coefficient ${\omega}_{D}$ is called plasma (or Langmuir) angular frequency (it is not the frequency of volume plasmon^{21}), with ${N}_{e}$ the density of electrons that contribute to the optical properties [Eq. (3), it is not the valence band electron concentration^{21}], $q$ the charge of electron, ${m}_{0}$ the electron mass, and ${\u03f5}_{0}$ the relative permittivity of vacuum.

The volume density of the electrons that contribute to the relative permittivity of gold ${N}_{e}$ (${\mathrm{m}}^{-3}$) is introduced. It can be evaluated from an atomistic model and is actually a function of the density ${\rho}_{m}$ and atomic weight $M$.

with ${N}_{a}=6.022\times {10}^{23}$ the Avogadro constant and ${n}_{e}$ the supposed number of free electrons in a single atom. This value can be compared to that obtained by the fitting in the following. The atomistic volume density of electron ${N}_{e}$ is calculated by using outer electron shell of atoms. Equation (1) could also be written as a product of ${\omega}_{D}^{2}$ by both Drude and Lorentz terms,^{21}but the formulation used in Ref. 7 is preferred to facilitate the comparison.

The physical parameters of the model are $\{{\u03f5}_{\infty},{\gamma}_{D},{N}_{e},\mathrm{\Delta}\u03f5,{\omega}_{L},{\gamma}_{L}\}$. Therefore the dimension of the problem is $\mathrm{dim}(D)=6$ for the fitting.

## 2.2.

### Combination of Drude and Critical Points Models

The critical points model (CP) was introduced in 1998 by Leng et al.^{24} for silicon and used for gold.^{9} This model describes the interband transitions in violet/near-uv region. It satisfies the Kramers-Kronig consistency and is able to reproduce the dispersion of gold with higher accuracy than the previous models. It includes a phase factor ($\varphi $) and corresponds to first-order poles in the complex plane.^{9}^{,}^{22}^{,}^{24} The interband transition angular frequency is $\mathrm{\Omega}$, and the transition broadening is governed by the damping angular frequency $\mathrm{\Gamma}$. It is recommended to take two critical points terms.^{9} This model helps to suppress the nonphysical absorption of the previous fitting in the region of transparency or near-transparency.^{24} In this study, the CP replaces the Lorentz additional term in Eq. (1), and therefore the model is a combination of Drude and CP (DCP).

## (4)

$${\u03f5}_{\mathrm{DCP}}(\omega )={\u03f5}_{\infty}-\frac{{N}_{e}{q}^{2}}{{\u03f5}_{0}{m}_{0}}\frac{1}{\omega (\omega +i{\gamma}_{D})}+\sum _{i=1}^{i=2}\mathrm{\Delta}{\u03f5}_{i}{\mathrm{\Omega}}_{i}[\frac{\mathrm{exp}(i{\varphi}_{i})}{{\mathrm{\Omega}}_{i}-\omega -i{\mathrm{\Gamma}}_{i}}+\frac{\mathrm{exp}(-i{\varphi}_{i})}{{\mathrm{\Omega}}_{i}+\omega +i{\mathrm{\Gamma}}_{i}}].$$This approach was used by Vial et al.^{11} for fitting the relative permittivities of gold and silver from Johnson and Christy.^{2} The CP involves 11 physical parameters: $\{{\u03f5}_{\infty},\text{\hspace{0.17em}}{N}_{e},\text{\hspace{0.17em}}{\gamma}_{D},\text{\hspace{0.17em}}{\mathrm{\Gamma}}_{i},\text{\hspace{0.17em}}{\mathrm{\Omega}}_{i},\text{\hspace{0.17em}}\mathrm{\Delta}{\u03f5}_{i},\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}{\varphi}_{i}\}$ for $i=1$, 2. The dimension of the problem of fitting is $\mathrm{dim}(D)=11$.

## 2.3.

### Criterion of Convergence for FDTD

According to Vial et al.,^{22}^{,}^{25} a criterion must be verified to check the convergence of the FDTD algorithm.

## (6)

$${\chi}_{0}=-{\left(\frac{{\omega}_{D}}{{\gamma}_{D}}\right)}^{2}[1-\mathrm{exp}(-{\gamma}_{D}\mathrm{\Delta}t)]+\frac{{\omega}_{L}^{2}}{{\gamma}_{D}}\mathrm{\Delta}t+\sum \mathfrak{R}(-i\frac{\eta}{\alpha -i\beta}\{1-\mathrm{exp}[(-\alpha +i\beta )\mathrm{\Delta}t]\}),$$In the case of combination of Drude and Lorentz model, ${\chi}_{0}$ is a function of $\alpha ={\gamma}_{L}/2$, $\beta =\sqrt{{\omega}_{L}^{2}-{\alpha}^{2}}$, and $\eta =\mathrm{\Delta}\u03f5{\omega}_{L}^{2}/\beta $. For the CP, ${\chi}_{0}$ is a function of $\alpha ={\mathrm{\Gamma}}_{i}$, $\beta ={\mathrm{\Omega}}_{i}$, and $\eta =2\mathrm{\Delta}{\u03f5}_{i}{\mathrm{\Omega}}_{i}\text{\hspace{0.17em}}\mathrm{exp}(-i{\mathrm{\Phi}}_{i})$.^{25}^{,}^{26} In Eq. (6), $\mathrm{\Delta}t$ is evaluated from the size of the grid $\mathrm{\Delta}x$ used for FDTD: $\mathrm{\Delta}t=\mathrm{\Delta}x/(2c)$, with $c$ the speed of light in vacuum. In the following, we consider $\mathrm{\Delta}x=1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$.^{25}

The criterion in Eq. (5) is a constraint for the method of fitting. Therefore, the random generation of particles $x$ [Eq. (11)] is stressed by this criterion, and the solution given by the NM algorithm is selected among those satisfying this criterion.

## 3.

## Combined Constrained PSO with Nelder-Mead Algorithm

The problem of data fitting belongs to the same class as the resolution of the inverse problem and as the optimization of systems. Classical methods were used for the fitting of relative permittivities: the NM method and the simulated annealing.^{7} The NM method^{27} searches solutions in an unbounded domain. This fact could be a drawback as it could lead to nonphysical solutions or divergence of the algorithm. The NM method was used for problems with small degree of freedom [$\mathrm{dim}(D)$].^{28} The simulated annealing subroutine belongs to the same class as the Monte Carlo method (heuristic methods) and is fully described by Kirkpatrick et al.^{29} and Tarantola.^{30} This method can be bounded or not. Both methods require a starting point for the algorithm initialization and the solution depends on it. The proposed PSO method enables to determine the starting point of the NM method within a space of parameters that ensures the criterion of convergence for FDTD [Eq. (5)].

Among heuristic bounded methods without starting point, the PSO method is widely used for optical applications.^{31}32.^{–}^{33} The PSO method was proposed in 1995 by Kennedy and Eberhart.^{34} It was used for the optimization of gold nanoshells using discrete dipole approximation (DDA) and Mie models for photothermal therapy,^{35} and more general problems of plasmonics including the optimization of nanostructures and the resolution of the inverse problem.^{5}^{,}^{36}^{,}^{37}

The PSO mimics the behavior of a swarm of particles moving to a potential well with analogy to the bees swarming in search of pollen. The nearby of particle with respect to the bottom of the potential well is evaluated through a fitness function $F$. In the present case, the fitness function $F$ is the standard deviation of the data computed from model ${\u03f5}_{M}$ to the ${N}_{\mathrm{ref}}$ values of the relative permittivity of reference ${\u03f5}_{r}$.

## (7)

$$F(x)=\sqrt{\frac{1}{{N}_{\mathrm{ref}}}\sum _{1}^{{N}_{\mathrm{ref}}}{|{\u03f5}_{M}(x)-{\u03f5}_{r}|}^{2}}.$$The fitness function $F(x)$ depends on the physical parameters $x$ that are necessary to compute the model ${\u03f5}_{M}$ for all $\omega $. The models for fitting ${\u03f5}_{M}(x)$ are described in Sec. 2. In the following, ${\u03f5}_{M}(x)={\u03f5}_{\mathrm{DL}}$ [Eq. (1)] or ${\u03f5}_{M}(x)={\u03f5}_{\mathrm{DCP}}$ [Eq. (4)]. According to this definition, for all physical parameters set $x$, the absolute error of fitting can be defined for the real part and the imaginary part of ${\u03f5}_{r}(\omega )$.

## (8)

$${\sigma}_{R}[{\u03f5}_{M}(\omega )]=\sqrt{\frac{1}{{N}_{\mathrm{ref}}}\mathfrak{R}{[{\u03f5}_{M}(\omega )-{\u03f5}_{r}(\omega )]}^{2}},$$## (9)

$${\sigma}_{I}[{\u03f5}_{M}(\omega )]=\sqrt{\frac{1}{{N}_{\mathrm{ref}}}\mathfrak{I}{[{\u03f5}_{M}(\omega )-{\u03f5}_{r}(\omega )]}^{2}}.$$If particles reach the position of the minimum of potential, then the fitness function is minimum; therefore a good position for a particle corresponds to a small value of the fitness function and to a good parameters set of the model. All inputs of the model (physical parameters) form a vector $x(t)$, whose size is the degree of freedom of the model. The simulated cinematic problem for particles depends on a virtual time $t$, which is the iteration step of the algorithm. Each particle moves in the space of physical parameters following Eq. (10).

The velocity $V(t+1)$ is computed from the best position $p(t)$ of each particle over previous generations up to step $t$ and from the swarm global best $g(t)$.^{38}

The inertia weight ${\omega}_{\mathrm{PSO}}$ linearly decreases from 0.9 to 0.4.^{39} The random gene ${U}_{i}(i=1,2)$] always checking the condition given in Eq. (5). The acceleration coefficients are ${c}_{1}={c}_{2}=2$. If the speed is too high, the particle leaves the domain of search. In this case, it is regenerated randomly in this domain. In this study, the number of particles in the swarm is $N=30$: at each step $t$, $N$ particles are moving in the space of search. In this basic version of the algorithm, the stop criterion is a maximum number of iterations of PSO (1000). We checked that the results are hardly dependent on all these exogenous parameters (${\omega}_{\mathrm{PSO}}$, ${c}_{1}$, ${c}_{2}$, and $N$). In Ref. 37, this property and the efficiency of this method was verified by using evolutionary method based on evolution strategies.^{40}41.^{–}^{42}

Then the solution obtained with PSO is used as starting point of an NM algorithm and the success of the method is revealed by the convergence rate to the same best solution. Only solutions that satisfy the condition in Eq. (5) are kept. A hundred realizations of the same algorithm help to check the stability of the solution. If the algorithm gives some identical solutions along realizations, the best solution is assumed to be determined. In the investigated cases, each realization takes a few minutes on a personal computer. This property is directly linked to the simplicity of the models that are detailed in Sec. 2.

## 4.

## Results and Discussion

We first validate the PSO method by comparison of results with those given in the literature. Then we apply the method to the fitting of gold, silver, chromium, aluminum, and titanium in the domain of visible wavelengths. The efficiency of the method is tested by comparisons of our results with published data,^{7}^{,}^{43} where the domain of wavelengths was $[\mathrm{500,1000}]\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$. Of course the influence of the domain of wavelengths for the fitting is crucial to the results.

In this study, the domain of wavelengths is $[\mathrm{400,800}]\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ and the domain of search for the physical parameters is $D$: $\{{\gamma}_{D},{\gamma}_{L}\}b[{10}^{13},{10}^{18}]\text{\hspace{0.17em}}\text{\hspace{0.17em}}(\mathrm{rad}.{\mathrm{s}}^{-1})$, ${\u03f5}_{\infty}\in [-\mathrm{1000,1000}]$, $\mathrm{\Delta}\u03f5\in [\mathrm{0,1000}]$, ${N}_{e}b[{10}^{27},{10}^{30}]\text{\hspace{0.17em}}\text{\hspace{0.17em}}({\mathrm{m}}^{-3})$ [Eq. (3)]. The fitting of bulk data for gold, silver, chromium, aluminum, and titanium are, respectively, investigated.

## 4.1.

### Gold

Gold has mechanical properties and an ease of implementation, which puts it at the forefront of materials used in plasmonics. Moreover, it cannot be oxidized by the common surrounding medium. Table 1 gives the best physical parameters obtained from the proposed method for the fitting of Johnson and Christy data.^{2} Figure 1 shows the results of fitting and the absolute errors ${\sigma}_{R}$ and ${\sigma}_{I}$ [Eqs. (8) and (9)].

## Table 1

Fitting of Johnson and Christy data (ϵr) for gold by the combination of Drude and Lorentz (DL) models [Eq. (1)] and critical points (DCP) model. F is the fitness function, C is the finite difference time domain (FDTD) criterion, σR is the absolute error between real parts of fitting and data, and σI is the absolute error between imaginary parts of fitting and data.

Au2 | ||
---|---|---|

DL | DCP | |

F [Eq. (7)] | 0.55 | 0.15992 |

C [Eq. (5)] | 0.99995 | 0.92761 |

σR | 1.3292 | 0.57097 |

σI | 1.2208 | 1.0747 |

ϵ∞ | 6.15991 | −9.06407 |

γD (rad/s) | 1.66938×1015 | 5.86665×1013 |

ωD (rad/s) | 1.34759×1016 | 1.30423×1016 |

Δϵ | 2.07122 | — |

ωL (rad/s) | 4.66171×1015 | — |

γL (rad/s) | 7.20958×1013 | — |

Ω1 (rad/s) | — | 3.2539×1016 |

Γ1 (rad/s) | — | 5.5350×1015 |

Δϵ1 | — | −10.8876 |

Φ1 (rad) | — | −2.46009 |

Ω2 (rad/s) | — | 3.91172×1015 |

Γ2 (rad/s) | — | 6.95449×1014 |

Δϵ2 | — | 0.718455 |

ϕ2 (rad) | — | −1.13717 |

The absolute error for the fittings using Eqs. (8) and (9) are compared to those in the literature.

• With DL, ${\sigma}_{R}({\u03f5}_{V})=4.5087$ and ${\sigma}_{I}({\u03f5}_{V})=3.6146$.

^{7}The proposed method gives ${\sigma}_{R}=\phantom{\rule{0ex}{0ex}}1.3292$ and ${\sigma}_{I}=1.2208$ (Table 1).• With DCP, ${\sigma}_{R}({\u03f5}_{V})=0.9417$, ${\sigma}_{I}({\u03f5}_{V})=1.0606$,

^{22}${\sigma}_{R}({\u03f5}_{V}^{\prime})=0.7882$, and ${\sigma}_{I}({\u03f5}_{V}^{\prime})=1.2268$.^{43}The proposed method gives ${\sigma}_{R}=0.57097$ and ${\sigma}_{I}=1.0747$ (Table 1).

These results show that the proposed method is competitive with the simulated annealing,^{7}^{,}^{22}^{,}^{43} leaving to smaller errors in the visible range. For DL model, the proposed method produces absolute errors which are less than one third of those found in the literature.

The fitness function for DL is three times that obtained for DCP. The DCP model is therefore more accurate to describe the optical properties of gold in the visible range.

The Palik data for gold come from two series of measurements, leading to a change of slope of the curve. These data are therefore suitable to test the method, by comparison with the fitting of Johnson and Christy data. Table 2 gives the best physical parameters. Figure 2 shows the results of fitting and the absolute errors ${\sigma}_{R}$ and ${\sigma}_{I}$ [Eqs. (8) and (9)]. The results of fitting for DL in Table 1 (Johnson and Christy) are close to those in Table 2 (Palik), on the contrary of DCP fitting. Figure 2 shows that DCP is highly sensitive to the change of slope of the real part of ${\u03f5}_{r}$. The result of fitting for the imaginary part is therefore affected by the change of slope (dashed curve).

## Table 2

Fitting of Palik data (ϵr) for gold by the combination of DL models [Eq. (1)] and DCP model. F is the fitness function, C is the FDTD criterion, σR is the absolute error between real parts of fitting and data, and σI is the absolute error between between imaginary parts of fitting and data.

Au1 | ||
---|---|---|

DL | DCP | |

F [Eq. (7)] | 1.1641 | 0.69043 |

C [Eq. (5)] | 0.99999 | 0.98792 |

σR | 7.7965 | 0.57097 |

σI | 2.7878 | 1.0747 |

ϵ∞ | 6.90939 | −6.23210 |

γD (rad/s) | 1.75628×1015 | 3.17891×1017 |

ωD (rad/s) | 1.38147×1016 | 1.90105×1017 |

Δϵ | 2.31858 | — |

ωL (rad/s) | 4.68266×1015 | — |

γL (rad/s) | 3.60439×1014 | — |

Ω1 (rad/s) | — | 3.37278×1015 |

Γ1 (rad/s) | — | 3.69375×1014 |

Δϵ1 | — | 0.419451 |

Φ1 (rad) | — | −4.44248 |

Ω2 (rad/s) | — | 2.04724×1015 |

Γ2 (rad/s) | — | 1.23992×1015 |

Δϵ2 | — | 28.2908 |

ϕ2 (rad) | — | −3.52016 |

The accuracy of both models given by the fitness function $F$ is less influenced by the chosen model (DL or DCP) than for the fitting of Johnson and Christy data. According to the previous results, the method can be applied to other metals with confidence.

## 4.2.

### Silver

Silver is known to be more efficient than gold for plasmonic nanostructures, although silver suffers from oxidization and changes of optical properties when nanostructured.^{44} Figure 3 shows both fittings and the error of fitting $\sigma (\omega )$. The real and imaginary parts of ${\u03f5}_{V}$^{22} are shown for comparison with those obtained from the proposed method (combination of PSO and NM algorithms). The two errors of fitting are about the same order of magnitude along $\omega $ but exhibit local differences between the reference ^{22} and the PSO results.

The optical properties of silver are known to have smooth behavior in the visible domain. Therefore, the fitness function, the absolute errors of fitting, and the angular frequencies ${\omega}_{D}$ are close together for both models. The amplitudes $\mathrm{\Delta}{\u03f5}_{i}$ of critical points contributions are larger than for gold.

The absolute errors for the fitting with DCP, using Eqs. (8) and (9), are compared to those in the literature for CP model: ${\sigma}_{R}({\u03f5}_{V})=0.66287$, ${\sigma}_{I}({\u03f5}_{V})=0.35788$,^{22} ${\sigma}_{R}({\u03f5}_{V}^{\prime})=0.81118$, and ${\sigma}_{I}({\u03f5}_{V}^{\prime})=0.70911$.^{43} They are roughly two times the values in Table 3: ${\sigma}_{R}=0.36891$ and ${\sigma}_{I}=0.34168$. The parameters for the DL model in Table 3 also give a better fit than that in Ref. 43. The comparison with results of fitting from Refs. 22 and 43 shows a great improvement of the quality of the solution, which is mainly due to a different interval of wavelengths for fitting.

## Table 3

Fitting of Palik data (ϵr) for silver by the combination of DL models [Eq. (1)] and DCP model. F is the fitness function, C is the FDTD criterion, σR is the absolute error between real parts of fitting and data, and σI is the absolute error between imaginary parts of fitting and data.

Ag1 | ||
---|---|---|

DL | DCP | |

F [Eq. (7)] | 0.071605 | 0.068131 |

C [Eq. (5)] | 0.98913 | 0.28263 |

σR | 0.37765 | 0.36891 |

σI | 0.33703 | 0.34168 |

ϵ∞ | 0.114773 | −19.4464 |

γD (rad/s) | 7.05499×1015 | 9.52533×1013 |

ωD (rad/s) | 1.32589×1016 | 1.32472×1016 |

Δϵ | 3.62762 | — |

ωL (rad/s) | 1.58116×1016 | — |

γL (rad/s) | 1.04632×1014 | — |

Ω1 (rad/s) | — | 3.47550×1016 |

Γ1 (rad/s) | — | 5.18324×1016 |

Δϵ1 | — | 1727.14 |

Φ1 (rad) | — | 1.04286 |

Ω2 (rad/s) | — | 4.76055×1016 |

Γ2 (rad/s) | — | 3.23542×1016 |

Δϵ2 | — | −742.061 |

ϕ2 (rad) | — | 1.75438 |

## 4.3.

### Aluminum

Aluminum is a good candidate to plasmon resonance as the real part of its relative permittivity is $<-23$ in the visible domain, even if this material can be easily oxidized.^{37} Table 4 gives the results of fitting for aluminum. These values are used to plot Fig. 4. Figure 4(b) illustrates the absolute error on the fitting of the real part and of the imaginary part. The comparison with results of fitting with DCP model (${\u03f5}_{V}$) from Refs. 22 and 43 are shown for comparison with those obtained from the PSO method. The errors of fitting between the reference ^{22}^{,}^{43} and the proposed fitting are about the same order of magnitude along $\omega $ except for wavelengths $>630\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$. The fitting with DL is much less efficient than DCP as DL cannot handle the change of slope of the real part of ${\u03f5}_{r}$ [see Fig. 4(a)].

## Table 4

Fitting of Palik data (ϵr) for aluminum by the combination of DL models [Eq. (1)] and DCP model. F is the fitness function, C is the FDTD criterion, σR is the absolute error between real parts of fitting and data, and σI is the absolute error between imaginary parts of fitting and data.

Al1 | ||
---|---|---|

DL | DCP | |

F [Eq. (7)] | 2.8752 | 0.29101 |

C [Eq. (5)] | 0.050662 | 0.27149 |

σR | 8.4152 | 1.5629 |

σI | 7.9598 | 1.1941 |

ϵ∞ | 0.0864080 | 5.79805×10−3 |

γD (rad/s) | 4.75280×1015 | 6.78974×1015 |

ωD (rad/s) | 3.91780×1016 | 2.22564×1017 |

Δϵ | −3814.53 | — |

ωL (rad/s) | 5.67222×1014 | — |

γL (rad/s) | 2.76612×1015 | — |

Ω1 (rad/s) | — | 8.10015×1013 |

Γ1 (rad/s) | — | 1.92552×1016 |

Δϵ1 | — | −1.35397×103 |

Φ1 (rad) | — | −0.683343 |

Ω2 (rad/s) | — | 2.25905×1016 |

Γ2 (rad/s) | — | 2.27300×1015 |

Δϵ2 | — | 3.59426 |

ϕ2 (rad) | — | −1.09763 |

Both results of fitting show that the Drude model is not suitable to describe the behavior of the optical properties of aluminum. Consequently, the higher degree of freedom of DCP makes it more efficient for the fitting. From Refs. 22 and 43, ${\sigma}_{R}({\u03f5}_{V})=2.664$ and ${\sigma}_{I}({\u03f5}_{V})=2.178$, where the interval of wavelengths was $[\mathrm{400,1000}]\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$. In Table 4, the corresponding values are ${\sigma}_{R}=1.5629$ and ${\sigma}_{I}({\u03f5}_{V})=1.1941$. The increase of absolute accuracy is around 2.

## 4.4.

### Chromium

Chromium is often used in nanotechnologies as adhesion layer for gold on dielectric substrates. Table 5 gives the results of fitting of Palik data for chromium. The corresponding dispersion curve is plotted in Fig. 5. Figure 5(b) illustrates the absolute error on the fitting of the real part and of the imaginary part. The comparison with results of fitting with DCP model (${\u03f5}_{V}$) from Refs. 22 and 43 are shown for comparison with those obtained from the PSO method. The errors of fitting are about the same order of magnitude along $\omega $ but exhibit local differences between the reference ^{22} and the PSO results.

## Table 5

Fitting of Palik data (ϵr) for chromium by the combination of DL models [Eq. (1)] and DCP model. F is the fitness function, C is the FDTD criterion, σR is the absolute error between real parts of fitting and data, and σI is the absolute error between imaginary parts of fitting and data.

Cr1 | ||
---|---|---|

DL | DCP | |

F [Eq. (7)] | 0.94723 | 0.077607 |

C [Eq. (5)] | 0.99993 | 0.88915 |

σR | 4.9121 | 0.29314 |

σI | 5.625 | 0.48649 |

ϵ∞ | 2.77667 | 0.386846 |

γD (rad/s) | 2.57644×1015 | 1.80690×1015 |

ωD (rad/s) | 1.59078×1016 | 1.39350×1016 |

Δϵ | 13.2908 | — |

ωL (rad/s) | 3.32972×1015 | — |

γL (rad/s) | 2.99660×1015 | — |

Ω1 (rad/s) | — | 3.79297×1015 |

Γ1 (rad/s) | — | 8.01774×1014 |

Δϵ1 | — | 2.12352 |

Φ1 (rad) | — | 0.883949 |

Ω2 (rad/s) | — | 1.75789×1015 |

Γ2 (rad/s) | — | 7.80710×1014 |

Δϵ2 | — | 11.8586 |

ϕ2 (rad) | — | −1.69593 |

The optical properties of a nanometric layer of chromium are known to influence the plasmon resonance of SPR.^{45} Therefore in the visible domain of wavelengths, the absolute errors for the fitting with DCP are significant. From Refs. 22 and 43, ${\sigma}_{R}({\u03f5}_{V})=0.93846$ and ${\sigma}_{I}({\u03f5}_{V})=0.58162$, where the interval of wavelengths was $[\mathrm{400,1000}]\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$. In Table 5, the corresponding values are ${\sigma}_{R}=0.29314$ and ${\sigma}_{I}({\u03f5}_{V})=0.48649$. The comparison with results of fitting from Refs. 22 and 43 shows a great improvement of the quality of the solution as shown in Fig. 5(b). Consequently, extending the domain to far-infrared domain is a drawback for the accurate fitting in the visible domain. The fitting of the optical properties of chromium must be adapted to the domain of wavelengths of interest.

## 4.5.

### Titanium

Titanium is used as an alternative medium for adhesion layer of gold on dielectric samples.^{45} Table 6 gives the results of fitting for titanium. These values are used to plot Fig. 6. Figure 6(b) shows the absolute error on the fitting of the real part and of the imaginary part. The fitting with DL is much less efficient than DCP as DL cannot handle the change of slope of the real part of ${\u03f5}_{r}$ [see Fig. 6(a)]. On the contrary, the imaginary part of the relative permittivity is well described by both models.

## Table 6

Fitting of Johnson and Christy data (ϵr) for titanium by the combination of DL models [Eq. (1)] and DCP model. F is the fitness function, C is the FDTD criterion, σR is the absolute error between real parts of fitting and data, and σI is the absolute error between imaginary parts of fitting and data.

Ti1 | ||
---|---|---|

DL | DCP | |

F [Eq. (7)] | 0.62496 | 0.26951 |

C [Eq. (5)] | 0.99959 | 0.95922 |

σR | 4.418 | 1.2955 |

σI | 1.4696 | 1.483 |

ϵ∞ | 2.17069 | 1.35312 |

γD (rad/s) | 1.08330×1016 | 3.02914×1015 |

ωD (rad/s) | 8.50868×1015 | 2.00924×1016 |

Δϵ | 74.4496 | — |

ωL (rad/s) | 2.93140×1015 | — |

γL (rad/s) | 1.46986×1014 | — |

Ω1 (rad/s) | — | 2.52097×1015 |

Γ1 (rad/s) | — | 2.09204×1015 |

Δϵ1 | — | 9.01823 |

Φ1 (rad) | — | −2.06436 |

Ω2 (rad/s) | — | 1.88081×1015 |

Γ2 (rad/s) | — | 9.05493×1013 |

Δϵ2 | — | 3.90173 |

ϕ2 (rad) | — | 2.76388 |

Both results of fitting show that the Drude model is not suitable to describe the behavior of the optical properties of titanium (Table 6).

## 4.6.

### Concluding Remarks

Both models of dispersion can be used according to the value of the fitness function $F$ (Tables 2 to 6), but Figs. 1 to 6 show local discrepancies of data fitting. In the general case, the error of fitting is spread over the real part or the imaginary part of the relative permittivity. The errors are globally lower compared to that found in the literature, especially for the DCP method. The DL model involves two terms and is less efficient than the CP model using three terms. Therefore, the DL model could be chosen if sparing memory is required for FDTD calculations. However, the DCP model is more accurate for all investigated metals. Aluminum’s optical properties are not well described by the Drude model. Depending on metals, both models are not equivalent, and the inspection of results allows to draw conclusions on the use of models that are summarized in Table 7. The performance of the method is better if the number of “+” increases. This performance is based on the value of the fitness function.

## Table 7

Performance of fitting of ϵr by the combination of DL models and DCP model [Eq. (4)] with the proposed method (particle swarm optimization+Nelder Mead), in the domain of visible wavelengths [400,800] nm. (P) indicates data from Palik1 and (JC) from Johnson and Christy.2

Au (JC) | Au (P) | Ag (P) | Al (P) | Cr (P) | Ti (JC) | |
---|---|---|---|---|---|---|

DCP | +++ | + | +++ | ++ | +++ | ++ |

DL | ++ | − | +++ | − | + | + |

The results of fitting of various metals in the visible bandwidth provide the following remarks:

• The DL model may be considered as competitive with the DCP model in terms of accuracy. The choice between these models for FDTD is based on a compromise between memory requirement that increases linearly with the number of terms and the required accuracy of numerical results.

• The DCP model is able to reproduce the change of slope of the real part of the relative permittivity on the contrary of the DL model.

• The FDTD parameter $\mathrm{\Delta}x=1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ is a strong constraint that could be released easily for applications where the accuracy on the near-field computation is not necessary. The proposed method could generate better fittings in this case.

## 5.

## Conclusion

The particle swarm method under constraint followed with an NM algorithm have been successfully used for the fitting of the relative permittivity of gold, silver, chromium, aluminum, and titanium in the visible wavelength domain. The PSO helps to find a good starting point for the NM method, and the generation of physical parameters that are compatible with their use for FDTD is ensured. The proposed method seems to be efficient for the hard problem of fitting of the relative permittivities of metals as shown by comparison with results in the literature.

Two models are involved for the fitting of the relative permittivities: combinations of the Drude model with Lorentz (DL) and with critical points models (DCP). The efficiency of both models for fitting is compared for each metal. The fitting coefficients of optical properties of metals that are currently used in optical design are given. The use of these results could be particularly useful in plasmonics and design of nanostructured biosensors, but beware of the crude use of these values: the domain of validity of these properties can change dramatically when the size of the metal particles is less than a few tens of nanometers.^{44}^{,}^{46}

The results of fitting can be used directly for any spectroscopic simulation and especially in FDTD codes.^{8}^{,}^{47} These results are complementary to those found in the literature where more than two DL terms are used for fitting. This method could also be applied to absorbing dielectrics, especially for nanowires and nanotube studies, with other models of fitting or other domains of wavelengths.^{48}49.^{–}^{50} The advantage of such a heuristic method lies in its applicability to various problems of fitting, the optimization of complex systems in engineering,^{13}^{,}^{35} the propagation of uncertainties, or the tolerance study of models.^{18}^{,}^{28} The use of more complex models ^{15}^{,}^{16}^{,}^{51}^{,}^{52} can extend the domain of its application.

## Acknowledgments

The authors thank the Agence Nationale de la Recherche for financial support (ANR-2011-NANO-008 NANOMORPH).

## References

## Biography

**Dominique Barchiesi** is a full professor in theoretical physics, applied mathematics, and statistics at the University of Technology of Troyes. He conducts active research in numerical modeling, optimization, and advanced methods for multiphysics, with application to engineering in nanotechnologies, biotechnologies, and plasmonics. He is the author of over 200 articles, conferences, and book chapters since 1992, in the fields of optics, electromagnetism, coupled problems, signal processing, operational research, and cryptography.

**Thomas Grosges** is a full professor in numerical and theoretical physics, and applied mathematics at the University of Technology of Troyes. He conducts active research on modeling the interaction of matter and light and on the optimization of numerical solvers in computer sciences and their applications (multiphysics, cryptography). He is the author of more than 50 publications in physics, applied physics, engineering, and computer science academic journals and books.