Fitting the optical constants of gold, silver, chromium, titanium, and aluminum in the visible bandwidth

Abstract. The fitting of metal optical properties is a topic that has applications in advanced simulations of spectroscopy, plasmonics, and optical engineering. In particular, the finite difference time domain method (FDTD) requires an analytical model of dispersion that verifies specific conditions to produce a full spectrum in a single run. Combination of Drude and Lorentz models, and Drude and critical points models, are known to be efficient, but the number of parameters to be adjusted for fitting data can prevent accurate results from simulated annealing or Nelder-Mead. The complex number relative permittivities of Au, Ag, Al, Cr, and Ti from either Palik or Johnson and Christy experimental data in the visible domain of wavelengths are successfully fitted by using the result of the particle swarm optimization method with FDTD constraint, as a starting point for the Nelder-Mead method. The results are well positioned compared to those that can be found in the literature. The results can be used directly for numerical simulations in the visible domain. The method can be applied to other materials, such as dielectrics, and to other domain of wavelengths.

Fitting the optical constants of gold, silver, chromium, titanium, and aluminum in the visible bandwidth

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 ω in a specific range of wavelengths in vacuum λ 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

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 ω on the contrary of the single Drude model. The DL model is efficient in the wavelengths range ½500;1000 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, ½400;800 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 ϵ DL of the relative permittivity of metal is written as the sum of the Drude and the Lorentz models. 6,10,22 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 γ D and γ L are introduced and expressed in the same units as ω. The plasma frequency ω L is associated with intraband transitions, Δϵ is the oscillator strength, and ϵ ∞ is the relative permittivity for high frequencies. For high frequencies, electrons cannot follow excitation, and the illumination does not see free electrons anymore; therefore ϵ ∞ ¼ 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 ω 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 ϵ 0 the relative permittivity of vacuum.
The volume density of the electrons that contribute to the relative permittivity of gold N e (m −3 ) is introduced. It can be evaluated from an atomistic model and is actually a function of the density ρ m and atomic weight M.
with N a ¼ 6.022 × 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 ω 2 D 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 fϵ ∞ ; γ D ; N e ; Δϵ; ω L ; γ L g. Therefore the dimension of the problem is dimðDÞ ¼ 6 for the fitting.

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 (ϕ) and corresponds to first-order poles in the complex plane. 9,22,24 The interband transition angular frequency is Ω, and the transition broadening is governed by the damping angular frequency Γ. 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).
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: fϵ ∞ ; N e ; γ D ; Γ i ; Ω i ; Δϵ i ; and ϕ i g for i ¼ 1, 2. The dimension of the problem of fitting is dimðDÞ ¼ 11.

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. with where Rð·Þ is the real part of a complex number.
In the case of combination of Drude and Lorentz model, 25,26 In Eq. (6), Δt is evaluated from the size of the grid Δx used for FDTD: Δt ¼ Δx∕ð2cÞ, with c the speed of light in vacuum. In the following, we consider Δx ¼ 1 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.

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 [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 ϵ M to the N ref values of the relative permittivity of reference ϵ r .
The fitness function FðxÞ depends on the physical parameters x that are necessary to compute the model ϵ M for all ω. The models for fitting ϵ M ðxÞ are described in Sec. 2. In the following, (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 ϵ r ðωÞ.
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 ω 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 (ω 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.

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 ½500;1000 nm. Of course the influence of the domain of wavelengths for the fitting is crucial to the results.

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 σ R and σ 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.
The absolute error for the fittings using Eqs. (8) and (9) are compared to those in the literature.
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 σ R and σ 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 ϵ r . The result of fitting for the imaginary part is therefore affected by the change of slope (dashed curve).
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.

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 σðωÞ. The real and imaginary parts of ϵ 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 ω 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 ω D are close together for both models. The amplitudes Δϵ 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: 43 They are roughly two times the values in Table 3: σ R ¼ 0.36891 and   Fig. 1 (a) Real and imaginary part of the experimental data for gold 2 and fitting with the Drudecritical points models (DCP) ( Table 1) Fig. 3 (a) Real and imaginary part of the experimental data for silver 1 and fitting with the DL and the DCP (Table 3) 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.

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 (ϵ 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 ω except for wavelengths >630 nm. The fitting with DL is much less efficient than DCP as DL cannot handle the change of slope of the real part of ϵ r [see Fig. 4(a)].
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, σ R ðϵ V Þ ¼ 2.664 and σ I ðϵ V Þ ¼ 2.178, where the interval of wavelengths was ½400;1000 nm. In Table 4, the corresponding values are σ R ¼ 1.5629 and σ I ðϵ V Þ ¼ 1.1941. The increase of absolute accuracy is around 2. 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.

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 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, σ R ðϵ V Þ ¼ 0.93846 and σ I ðϵ V Þ ¼ 0.58162, where the interval of wavelengths was ½400;1000 nm. In Table 5, the corresponding values are σ R ¼ 0.29314 and σ I ðϵ 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.

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 ϵ r [see Fig. 6(a)]. On the contrary, the imaginary part of the relative permittivity is well described by both models.
Both results of fitting show that the Drude model is not suitable to describe the behavior of the optical properties of titanium (Table 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,  Fig. 4 (a) Real and imaginary part of the experimental data for aluminum 1 and fitting with the DL and the CP (Table 4). (b) The error between the fitting and the reference data. The fitting with parameters from Refs. 22 and 43 (ϵ V ) are plotted. 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.   (Table 5). (b) The error between the fitting and the reference data. The fitting with parameters from Refs. 22 and 43 (ϵ V ) are plotted. 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.   (Table 6). (b) The error between the fitting and the reference data.
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.
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 Δx ¼ 1 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.

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.