Spectral curve fitting of dielectric constants

Optical constants are crucial in determining how a substance responds to different frequencies of incident light. As an intrinsic property of the material,¹ optical constants contain information on the composition of material,^2–4 its anisotropy,^5–7 and its optical activity.⁸ The optical properties of a material can be inferred from experimental spectra from X-ray frequencies through to the Far-Infrared/Terahertz regions, in either transmittance or absorbance or reflectance mode.^2,9–13 From the measured spectra, it is possible to extract both the real and imaginary components of the complex optical constants by employing the Kramer-Kronig (K-K) relation, a mathematical expression that connects the real and imaginary parts.¹⁴ However, due to the limited range of available spectral data, the optical constants extracted by direct K-K analysis contain truncation errors.¹⁵ 

Curve fitting is an alternate method to recover optical constants from spectral data.^16,17 By employing physically meaningful line shape functions that are K-K consistent, it is possible to retrieve optical constants to high accuracy. In addition, it is possible to conveniently express these optical constants just using a few parameters, instead of supplementing entire spectral data tables.^14,17–22

The Lorentzian line shape is an obvious choice as it is a simple curve that has theoretical origin from optical transitions,^23–25 but it is not always the most appropriate. Inhomogeneous effects which are intrinsically Gaussian in nature convolve the Lorentzian into a Voigt line shape. It has been realised that the Lorentzian line shape does not satisfactorily describe spectral features of amorphous solids^4,5 or interband electronic transitions of noble metals,²⁶ whereas the Voigt line shape can represent their dielectric constants. A more versatile but less physically intuitive way to model optical constants is to use K-K consistent functions which are capable of fitting various line shapes.^{7,17,22,25,27–30} However, there is no comparison of these methods in a systematic and comprehensive manner. The purpose of this manuscript is to scrutinize four commonly employed line shape functions using representative molecular crystal systems, and provide insight on their applicability to optical constant fitting and K-K consistency.

Refractive indices m_v (m_v = n_v + ik_v), along with the dielectric constant $εv(εv=ε˙v+iε¨v)$ and, the polarizability $αv(αv=α˙v+iα¨v)$⁠, are the wavelength dependent optical constants of a material. These optical constants are interrelated,¹ i.e.

Polarizability is related to the dielectric constant and thereby the refractive index by the Lorentz-Lorentz relation,^1,31

In certain instances, these optical constants are directly linked to measured spectra. The imaginary part of refractive indices, k_v, is proportional to the absorbance of thin films¹ while the imaginary part of polarizability, $α¨v$⁠, is proportional to the absorbance of small particles.^32,33 The area under the curve bounded by $α¨v$ and wavenumber v (⁠$Area=∫a¨vdv$⁠) is directly related to the transition dipole moment of the material.¹⁴ Therefore, macroscopic behaviour of the material is connected via polarizability to its molecular level property.

The real and imaginary parts of optical constants (whether refractive indices, polarizability or dielectric constants) are related through K-K transformation.^1,34 For consistency with most previous methodology papers, this work focusses on dielectric constants. The real part of the dielectric constant may be computed by performing a K-K transformation on the imaginary part,

where P denotes the Cauchy principal value,³⁵ and $ε˙(v∞)$ is the real part of dielectric constant at infinite frequency; which can be approximated as the real part of dielectric constant at visible wavelength for colourless materials.³⁶ 

The complex dielectric constant and polarizability can be expressed as a function of band parameters using the classical damped harmonic oscillator (CDHO) model.¹⁴ In the CDHO model, if the j^th oscillator intensity is S_j (cm^-2), the damping constant is $γj$ (cm^-1), and the band position is v_j (cm^-1) then the dielectric constant is:^14,19

Note that as the CDHO line shape is a Lorentzian function, it can not correctly reproduce Gaussian or asymmetric line shapes.³⁷ 

The Voigt line shape is a convolution of Lorentzian and Gaussian functions that has been shown to accurately describe the spectral profiles of amorphous solids.^5,18,20 The Voigt function cannot be expressed analytically in closed form, however, it can be calculated using various approximations.^38,39 Specifically, the Voigt function, K(x,y), is the real part of Faddeeva function, w(z):³⁸ 

Here, z = x + i y, where x and y are related to peak position, and the Lorentzian and Gaussian peak width.³⁸ The complex optical constants can be described by a slightly modified version of the Faddeeva function¹⁸ that satisfies the causality rule. In this work we evaluate the Faddeeva function using the exponential series approximation as employed by Abrarov.³⁸ 

If the imaginary part of dielectric constant, $ε¨(v)$⁠, is significantly non-zero only in the range -w₁ < v < w₁, then the imaginary and real parts can be expressed as Fourier sine-cosine series:^29,30

As optical constants are rarely non-zero outside the measured range, the optical data can be remapped so that the Equations (8) and (9) are applicable. The idea is to extrapolate the imaginary part linearly until its value is zero, zero-fill both above and below the available range, and then form a new data set that is antisymmetric about the origin. The negative axis values can be set to $ε¨(−v)=−ε¨(v)$ to allow the new data set to be fitted to Equation (8). The p_k values returned can then be inserted into Equation (9) to calculate the real part. The remapping is illustrated in Figure 1 for dielectric constants of water ice⁴⁰ in the mid-IR range.

It has been shown that triangular basis functions can be used to model dielectric constants.^7,17,27 Note that a triangular shape is B-spline of degree one but B-splines of various degrees can also be employed.¹⁷ The idea is to find a set of triangular basis functions $ε¨j(v)$ anchored at $νj$ so that the imaginary part of the dielectric constant can be expressed as a linear combination of these basis functions.^6,27 One such triangular function is:¹⁷ 

The real part of the above triangular function is:¹⁷ 

The bases can be decided variationally^17,27 or manually chosen, with the simplest scheme being to assign a basis at every other point.⁶ Once these basis functions are determined, the dielectric constants can be expressed as:²⁷ 

Here N is the number of basis functions, which is at most half the total number of data points, and the coefficients of first and last basis function (A₁ and A_N), should be set to zero.²⁷ 

The fitting methods described here are generally applicable to other spectral regions, subject to the conditions previously outlined concerning peak widths and separation of the region of interest from other absorptions at high or low frequency, but we have focussed on the IR spectra of molecular ices according to our particular research interests. In this study we have selected example substances with optical constants to represent a range of dielectric materials. In the context of the tests carried out here, it is better to use available optical data so that: i) the fitted data can be compared to the original data allowing the absolute error for each method to be quantified, ii) optical data of various line profiles can be tested to establish the best fitting model for each of the example materials. The data chosen for this study includes the refractive indices for water ice⁴⁰ at 266 K, acetonitrile ice⁴¹ at 95 K, and the nitric acid dihydrate (NAD)⁴² at 160 K. Water ice spectral bands are very broad (>100 cm^-1) while acetonitrile peaks are very narrow in comparison (∼2 cm^-1), and nitric acid dihydrate spectra have a mixture of sharp and diffuse bands. We have converted the reference refractive indices to dielectric constants (Equation 1), which is a more generalised way to represent optical constants. The dielectric constants of these ices in the mid IR region are shown in Figure 2. The $ε˙(∞)$ value taken for these ices are 1.712,⁴⁰ 2.045,⁴¹ and 2.35 respectively.⁴³ 

The fitting of dielectric constants for various ices are performed by minimising the sum of squared residuals (SSR), i.e.:

Here $εi$ and $εtrue$ are fitted and true (experimental) dielectric constants, respectively. When evaluating the K-K consistency of various methods (i.e. Figs. 4–6), the ‘true’ value of $ε˙$ is derived from K-K transformation of the $ε¨$ values between 0 and 10,000 cm^-1 at intervals of 0.1 cm^-1. Note that in this range the ‘truncation’ error is in the order of 10^-6 which is at least two orders of magnitude smaller than other errors, so the truncation error is negligible.

The object function to be minimised can take various forms,⁴ i.e. SSR of the imaginary part, real part or some combination of both. Depending on whether the spectral data is absorbance/transmittance or reflectance mode that nonlinearly depends on the both parts of dielectric constant, a carefully chosen objective function can minimise the errors, for details see references 4, 26. We chose SSR of the imaginary part as the object function, because (1) under certain experimental conditions the imaginary part of some optical constants are directly proportional to measured spectra (see introduction) and (2) the findings outlined in this paper do not change if other forms of object function are chosen. All the calculations are performed using Matlab (R2014b)⁴⁴ on a PC. The object function minimisation is done using a built-in nonlinear least square algorithm. Fitting times are reported to provide an approximate comparison of computational efficiency. Coefficients of determination, adjusted R², are reported to aid assessment of the fit quality as R² accounts for the number of parameters used and allows for better comparison between spectra in which the dielectric constants have quite different magnitude.

The error associated with K-K transformations within a limited data range is thoroughly discussed in Reference 13 by Segal-Rosenheimer and Linker.¹⁵ The key conclusions are: By neglecting the imaginary component of low frequency bands, the K-K transformation overestimates the real part and the associated error is largest at the lowest available frequency; By neglecting the imaginary component of high frequency bands, the K-K transformation underestimates the real part and the associated error is largest at the highest available frequency.

We also find that if the data interval is not small enough to cover the details of a peak (at least 10 points per FWHM), the K-K transformation also introduces errors that are comparable to or larger than truncation errors, especially around the peak maximum. To illustrate this point we compute dielectric constants using asymmetric Lorentzian³⁵ profile with v₀=2500 cm^-1, $ε¨max=1$⁠, and full width at half maximum (FWHM) = 10 cm^-1. Shown in Figure 3 (a) are the simulated dielectric constants, and (b-d) are the errors in the real part upon K-K transformation of the imaginary part between 500 – 4500 cm^-1, at intervals of 0.5 cm^-1 (Fig. 3b), 2 cm^-1 (Fig. 3c), and 5 cm^-1 (Fig. 3d solid line). The error in Fig. 3(b) is due to truncation error, as demonstrated by Segal-Rosenheimer¹⁵ while errors in Fig. 3 (c) and (d) are clearly due to the sampling interval. We suggest the oscillatory errors originate from the method used,³⁵ and its application of odd or even parity of input data indices to perform the K-K transformation. However, similar errors are observed upon Hilbert transformations, which approximate K-K transformations.³⁴ 

When the data interval is too large (smaller than 2.5 points per FWHM), the errors can be improved over the whole spectral range by interpolating between every data points. Upon interpolation, the errors on the wings of the peak decrease most significantly while around peak maximum the error can be almost halved. In Figure 3 (d) we interpolated (cubic spline) the imaginary part between every point (halving the interval to $Δν$ = 2.5 cm^-1), and the error after K-K transformation shown in dotted line is considerably reduced. Interpolation at smaller data intervals does not further decrease the overall error significantly, presumably due to the approximations inherent in the interpolation. Thus the data interval requirement for accurate K-K transformation is more strict than the peak resolving power requirement in spectroscopic measurements.⁴⁵ In FTIR spectroscopy in case the spectral data are measured at lower resolution, the data interval can be decreased by post-processing the original interferogram with increased zero filling or as demonstrated here, by interpolation.

In this section, we test K-K consistency between real and imaginary dielectric constants derived by the various methods. None of the methods besides the CDHO model⁴⁶ strictly satisfy the K-K relation due various approximations, as will be shown later. Shown in Figure 4(A) is the Voigt modelling for the imaginary component of the dielectric constant, and (B) its real part residual. As can be seen from Figure 4 and Table I, the errors increase as the peaks become broader and more Lorentzian. The errors are due to the K-K inconsistency of the real and imaginary parts,¹⁸ and to lesser extent, the errors in the evaluation of the Feddeeva function.³⁸ Therefore, the Voigt line shape may be applicable only to species that exhibit narrow spectral bands.

The dielectric constants associated with a Fourier series fitting procedure is shown in Figure 5 (a), computed for a Lorentzian band (FWHM = 10 cm^-1) using 1000 terms. It is clear from Figure 5 (b) that compared to the benchmark K-K values, the real part $ε˙$ is underestimated across the whole spectral range. The source of the errors are due to i) the finite number of points in the spectrum, meaning the maximum number of terms are limited,²⁹ and ii) the sine functions are all zero at zero wavenumber but are not necessarily zero beyond the high frequency limit of the spectrum. If spectral range is limited, which is common, the error associated with the Fourier series fitting method becomes similar to the truncation error in K-K transformation,¹⁵ i.e. the values are underestimated at high frequency and overestimated at low frequency.

Shown in Figure 6(a) are the dielectric constants obtained for a Triangular function of FWHM 10 cm^-1 centred at 2000 cm^-1, and in 6(b and c) its real part residual compared to the benchmark K-K values. The ‘bumps’ around 2000 cm^-1 (Figure 6(c)) coincide with the triangle basis function’s anchor position. It is clear that the real part is underestimated at low wavenumber and overestimated at high wavenumber side of the spectrum. The calculated errors do tend to decrease as the band becomes narrower. However, if a spectral feature is broad, then more triangular bases are needed for accurate fitting, thus the errors can sum up to larger values.

We now evaluate how well each method reproduces the dielectric constants of crystalline water ice. Shown in Figure 7 are the dielectric constant residuals for (a) CDHO fitting using 16 bands, (b) Voigt fitting using 16 bands, (c) Fourier series fitting using 100 terms, and (d) Triangular basis fitting using 500 triangles. It is clear that the fits to the imaginary part $ε¨$ are all reasonably good except for an edge artefact of the Triangular fitting at lower wavenumber end (not shown here, the imaginary and real part residuals go to -0.8 and 1.3, respectively). The real part $ε˙$ errors are much more diverse, though all are worst at the low wavenumber limit due to the non-zero value of the libration band at the edge of the spectrum. The CDHO and Voigt fits can be slightly improved by increasing the number of bands but, as can be seen from Table II, the improvement after more than 16 bands is not significant and the fitting procedure becomes inefficient. Note that the real part error in Voigt fitting is large and does not decrease much even after using 30 bands. Also note the slow convergence of Voigt fitting.

The Fourier series fit of the imaginary part can be improved by increasing the number of terms (Table II), but the real part residuals do not decrease significantly after inclusion of 200 terms. By increasing the series terms, the oscillatory part of the residuals can be decreased but the general trend in the real part residual plot is not improved, since it is due to the intrinsic error of the Fourier method (see Figure 5).

The error in Triangular basis fitting for using 500 triangular bases is shown in Figure 7 (d). By using more triangle bases, the errors at the edge of spectra (∼700 cm^-1) can be decreased significantly but throughout the other regions of the spectra, the errors stay the same. The errors can also be decreased by linear extrapolation of the spectra to lower wavenumber prior to the fitting, similar to the findings of a previous study.¹⁷ 

As crystalline acetonitrile ice peaks are very sharp, all the methods can be readily employed to fit the optical constants with reasonable accuracy. Shown in Figure 8 are the residuals from fitting the acetonitrile dielectric constants by the (a) CDHO (28 bands), (b) Voigt (28 bands), (c) Fourier series (4000 terms), and (d) triangle fitting (4281 basis) procedures. All except Voigt fitting can be improved upon increasing the number of bands or terms (Table III). Also note that the fitting is fastest (a few seconds) for the CDHO method due to its analytic nature and availability of analytic derivatives.³⁷ The triangular fitting converges similarly fast. The Fourier series fitting is also fast (tens of seconds), and it may be possible to accelerate further by evaluating series coefficients by Fast Fourier Transform.²⁹ The Voigt fitting converges very slowly (a few hundreds of seconds). Because there are many peaks and most peaks are narrow many terms and bases are required for Fourier series and Triangular fitting.

Nitric acid dehydrate (NAD) has a combination of narrow peaks as well as broad bands, meaning some of its spectral features are comparable to those of water, while others are more similar to acetonitrile. The bands above 1500 cm^-1 are very broad (∼100 cm^-1), which translates to large errors associated with the Voigt lineshape (see Figure 4 and Table I). The residuals of the NAD dielectric constants are shown in Figure 9 for (a) CDHO fit (30 bands), (b) Voigt fit (30 bands), (c) Fourier series fit (1000 terms), and (d) Triangle fit (500 terms). The CDHO and Fourier series fitting are of high accuracy and both converge within a few seconds. Note that in the Triangular basis fitting procedure, even though the error in imaginary part is the smallest among all methods (see Table IV), the error in the real part is substantial. This is expected considering the intrinsic error of the Triangular basis method (see Section IV B and Figure 6).

The accuracy of the four curve fitting methods that are commonly employed to extract optical constants from spectral data are evaluated, and in our study applied to model the dielectric constants of water, acetonitrile, and nitric acid dihydrate ices. These three materials have spectra with quite different characteristics. For water ice where spectral bands are very broad, the Fourier series fitting is most accurate and computationally rapid, while the CDHO fitting can be as fast but the accuracy is lower. This is due to the nonlocal nature of CDHO line shape where, for broad bands like water ice, a CDHO line has non-negligible contribution even far from the line centre. For broad bands, therefore, CDHO fitting is best suited for a species where the spectral intensity is significantly nonzero over the whole spectral range, as in the case for nitric acid dihydrate. For spectra where peaks are narrow (of the order of a few cm^-1), like acetonitrile, the CDHO fitting as well as the Triangular fitting are the most accurate and efficient of the methods. The Voigt fitting procedure is particularly slow even though the accuracy is comparable to CDHO fitting. The Fourier series fitting seems universally applicable to narrow as well as broad spectral bands with reasonable accuracy and speed. Note that in the cases of water and nitric acid dihydrate ice, whose spectral features are not necessarily narrow, the Fourier series fitting is superior to the remaining methods when considering accuracy and speed of convergence.

One advantage that is unique to CDHO fitting is the dielectric constants on the whole spectral range can be expressed using just a few parameters with reasonable accuracy.⁴⁷ Another advantage of CDHO fitting is the derivatives of the parameters to be optimised are readily available.³⁷ Thus CDHO fitting is one convenient method when further analysis is needed following the optical constants fitting. For example, when both refractive indices and particle size parameters are to be extracted from large aerosol particle spectra, the CDHO fitting can be employed where only a handful of parameters are need optimisation.^37,48

Based on the quantitative analysis in this study, we conclude that the Fourier series fitting is the most desirable method, especially when spectral features are not very narrow (> 10 cm^-1). In cases where the spectral features are narrow (< 10 cm^-1) the CDHO and the Triangular fitting are the most accurate and fast of the four methods. The Voigt fitting can be accurate for narrow peaks but the very slow convergence speed is discouraging.

This project has been supported by the Australian Research Council via a Discovery Early Career Researcher Award (DE150100301). M. Ruiz thanks La Trobe University for funding via Postgraduate Scholarship. We like to acknowledge the helpful discussions with Dr. Alexey Kuzmenko. We also kindly thank the Australian Synchrotron and staff at the THz/Far-IR beamline.