Optimizing the Drude-Lorentz model for material permittivity: Examples for semiconductors

TL;DR

This paper presents an optimized method for fitting the Drude-Lorentz model to experimental permittivity data of semiconductors, improving the accuracy of optical property modeling across different frequency ranges.

Contribution

It introduces a systematic optimization approach for the Drude-Lorentz model parameters tailored for semiconductors, enhancing the modeling of their permittivity dispersion.

Findings

Optimized parameters for Si, GaAs, and Ge across various frequency ranges.

Up to five pairs of Lorentz poles can be effectively fitted.

The method reduces the deviation between model and measured permittivity.

Abstract

Approximating the frequency dispersion of the permittivity of materials with simple analytical functions is of fundamental importance for understanding and modeling their optical properties. Quite generally, the permittivity can be treated in the complex frequency plane as an analytic function having a countable number of simple poles which determine the dispersion of the permittivity, with the pole weights corresponding to generalized conductivities of the medium at these resonances. The resulting Drude-Lorentz model separates the poles at frequencies with zero real part (Ohm's law and Drude poles) from poles with finite real part (Lorentz poles). To find the parameters of such an analytic function, we minimize the error weighted deviation between the model and measured values of the permittivity. We show examples of such optimizations for various semiconductors (Si, GaAs and Ge), for different frequency ranges and up to five pairs of Lorentz poles accounted for in the model.

Figures

Tables

Table 1. Table 1: Optimized model parameters for different semiconductors, using the fit function with five pairs Lorentz pole and optimization energy ranges corresponding to the data shown in Figs.  2 , 4 , and 6 . The number of data values 2 ​ N 2 𝑁 2N , the number of fit parameters 1 + 4 ​ L 1 4 𝐿 1+4L , and the resulting error S 𝑆 S are also given.

Material	Si	GaAs	Ge
	$0.81568$	$-0.54651$	$0.79842$
	$3.3736$	$1.4377$	$1.168$
${\Omega }_1^{\prime \prime }\left(\mathrm{eV}\right)$	$-0.11402$	$-0.05948$	$-0.33778$
${\sigma }_1^{\prime }\left(\mathrm{eV}\right)$	$1.6934$	$0.01981$	$0.47159$
${\sigma }_1^{\prime \prime }\left(\mathrm{eV}\right)$	$2.084$	$0.01122$	$0.01002$
	$3.6519$	$2.7229$	$2.174$
${\Omega }_2^{\prime \prime }\left(\mathrm{eV}\right)$	$-0.52378$	$-1.2972$	$-0.28077$
${\sigma }_2^{\prime }\left(\mathrm{eV}\right)$	$5.2573$	$7.8336$	$3.2926$
${\sigma }_2^{\prime \prime }\left(\mathrm{eV}\right)$	$8.0106$	$8.3274$	$4.1239$
	$4.2877$	$2.8922$	$3.781$
${\Omega }_3^{\prime \prime }\left(\mathrm{eV}\right)$	$-0.21116$	$-0.23992$	$-1.1461$
${\sigma }_3^{\prime }\left(\mathrm{eV}\right)$	$-1.7164$	$2.706$	$0.86584$
${\sigma }_3^{\prime \prime }\left(\mathrm{eV}\right)$	$5.9939$	$1.616$	$18.898$
	$5.3188$	$4.5222$	$4.3232$
${\Omega }_4^{\prime \prime }\left(\mathrm{eV}\right)$	$-0.18434$	$-0.42072$	$-0.20006$
${\sigma }_4^{\prime }\left(\mathrm{eV}\right)$	$-0.00528$	$2.1137$	$-1.7377$
${\sigma }_4^{\prime \prime }\left(\mathrm{eV}\right)$	$0.32911$	$4.6445$	$2.5278$
	$5.5064$	$4.9278$	$5.6442$
${\Omega }_5^{\prime \prime }\left(\mathrm{eV}\right)$	$-1.7892$	$-0.19972$	$-0.41214$
${\sigma }_5^{\prime }\left(\mathrm{eV}\right)$	$-3.8438$	$-1.243$	$0.10451$
${\sigma }_5^{\prime \prime }\left(\mathrm{eV}\right)$	$6.9298$	$1.4424$	$1.0292$
	$1.0$	$1.3$	$0.5$
$\hslash {\omega }_N\left(\mathrm{eV}\right)$	$7.6$	$6.0$	$6.0$
$2N$	$662$	$190$	$222$
$1+4L$	$21$	$21$	$21$
$S$	$0.01016$	$0.01157$	$0.01327$

Equations

$$\hat{\mathbf{\epsilon }}(\omega )={\hat{\mathbf{\epsilon }}}_{\infty }+\sum _j\frac{i{\hat{\mathbf{\sigma }}}_j}{\omega -{\Omega }_j},$$

$$\epsilon (\omega )={\epsilon }_{\infty }+\sum _{k=1}^L\left(\frac{i{\sigma }_k}{\omega -{\Omega }_k}+\frac{i{\sigma }_k^{\ast }}{\omega +{\Omega }_k^{\ast }}\right),$$

$$E=\sum _{j=1}^N{\left(\frac{{\epsilon }^{\prime }({\omega }_j)-{\epsilon }_j^{\prime }}{\Delta {\epsilon }_j^{\prime }}\right)}^2+{\left(\frac{{\epsilon }^{\prime \prime }({\omega }_j)-{\epsilon }_j^{\prime \prime }}{\Delta {\epsilon }_j^{\prime \prime }}\right)}^2,$$

$$\begin{gathered} \epsilon ({\omega }_j)={\epsilon }_{\infty }+\sum _{k=1}^L[{\sigma }_k^{\prime }\left(\frac{i}{\omega -{\Omega }_k}+\frac{i}{\omega +{\Omega }_k^{\ast }}\right) \\ +{\sigma }_k^{\prime \prime }\left(\frac{-1}{\omega -{\Omega }_k}+\frac{1}{\omega +{\Omega }_k^{\ast }}\right)] \end{gathered}$$

$$S=\sqrt{\frac{\tilde{E}}{2N}}.$$

$${\Omega }_1={\omega }_1{\left(\frac{{\omega }_N}{{\omega }_1}\right)}^Y-i({\omega }_N-{\omega }_1)N^{Y^{\prime }-1}$$

Full text

I Introduction

In this paper we expand on our recent work [Sehmi16, ] on fitting a generalized Drude-Lorentz (DL) model to the experimental permittivity of metals (gold, silver and copper), by applying the same method to common semiconductors. The fitting procedure employs an error minimization using a semi-analytic approach in which we first exactly determine the pole weights for given pole frequencies by solving a set of linear equations. This reduces the number of remaining fit parameters by about a factor of two. These parameters, which are the pole frequencies, are then determined numerically using a gradient decent.

The use of an analytical model of $\hat{\boldsymbol{\varepsilon}}(\omega)$, which contains only simple poles, is motivated by physical arguments, such as the presence of resonances in the material self energy and response functions. Furthermore, this form of the permittivity can be efficiently implemented in numerical methods, such as the finite difference in time domain (FDTD) method  [VialJPD07, ], and in more analytic and rigorous approaches, such as the dispersive resonant-state expansion (RSE) [MuljarovPRB16, ]. While real metals are described well by the Drude model and as a result have a dominating contribution of Ohm’s law and Drude poles, the permittivity of semiconductors is dictated by the interband transitions which are represented by Lorentz components. We will therefore use here only Lorentz terms for $\hat{\boldsymbol{\varepsilon}}(\omega)$, and call the model a Lorentz model.

We use an efficient algorithm of fitting experimental data with the Lorentz model with an arbitrary number of Lorentz poles, which takes into account experimental errors if available. This algorithm combines an exact analytical approach for the linear parameters of the model, with a numerical solver for the optimization of the nonlinear parameters of the model. This method, which is discussed in detail in our previous paper [Sehmi16, ], results in a significant speed-up of the process of finding a complete fit for a set of measured data. We have also developed methods of quickly finding the non-linear parameters with MATLAB’s fminunc function through suited selection of start values. We show examples of fitted data for Si, GaAs and Ge, for different frequency ranges and up to 5 pairs of Lorentz poles. We illustrate the resulting pole positions and their weights in the complex plane to give some physical insight how the model approximates the electronic transitions in real materials.

The paper is organized as follows. Sec. II introduces a generalized DL model of the permittivity and explains why we use in the present work a reduced version of this model. The fit procedure, including the analytical and numerical optimization of parameters of the DL model and the algorithm for determining appropriate starting values in the gradient descent minimization is summarized in Sec. III. Results of the fit are provided in Sec. IV for the semiconductor crystals silicon, gallium arsenide, and germanium.

II Lorentz model

Quite generally, the permittivity $\hat{\boldsymbol{\varepsilon}}(\mathbf{r},\omega)$ can be treated as an analytic function in the complex frequency plane, having a countable number of simple poles. Then, according to Mittag-Leffler theorem, it can be expressed as

[TABLE]

where $\hat{\boldsymbol{\varepsilon}}_{\infty}$ is the high-frequency value of the permittivity and $\Omega_{j}$ are the resonance frequencies – poles of the permittivity, determining its dispersion, with the weight tensors $\hat{\boldsymbol{\sigma}}_{j}$ corresponding to generalized conductivities of the medium at these resonances. The Lorentz reciprocity theorem requires that all tensors in Eq. (1) are symmetric, and the causality principle requires that $\hat{\boldsymbol{\varepsilon}}(\omega)$ has no poles in the upper half $\omega$ plane and that $\hat{\boldsymbol{\varepsilon}}^{\ast}(\omega)=\hat{\boldsymbol{\varepsilon}}(-\omega^{\ast})$ [LandauLifshitzV8Book84, ]. Therefore, for a physically relevant dispersion, each pole of the permittivity with a positive real part of $\Omega_{j}$ has a partner at $\Omega_{-j}=-\Omega^{\ast}_{j}$ with $\hat{\boldsymbol{\sigma}}_{-j}=\hat{\boldsymbol{\sigma}}^{\ast}_{j}$. Poles with zero real part of $\Omega_{j}$ have real $\hat{\boldsymbol{\sigma}}_{j}$. For simplicity, we assume in the following an isotropic response, such that the conductivities and thus the permittivity are described by scalars. We note however that it is straightforward to extend the presented treatment to anisotropic response.

We can separate the poles with zero real part of the frequency, which describe the conductivity of materials in the Drude model. There is a pole at zero frequency which represents Ohm’s law, corresponding to the $\omega^{-1}$ low-frequency limit of the dispersion. Together with a second pole which has a negative imaginary part, it provides the $\omega^{-2}$ high-frequency asymptotics, originating from the non-zero mass of the charge carriers. In real materials, the carrier mass and the damping can show a frequency dependence, which is not included in the Drude model. To describe such effects, the DC conductivity can be split  [AllenPRB77, ] into several Drude contributions.

In this work we will describe semiconductors which have a negligible free carrier density, and thus a negligible Drude pole weight. Their susceptibility in the visible and ultraviolet range is dominated by interband transitions. We thus use a version of $\varepsilon(\omega)$ which contains only pairs of Lorentz poles,

[TABLE]

where $L$ is the number of such pairs. Both the pole positions $\Omega_{k}=\Omega^{\prime}_{k}+i\Omega^{\prime\prime}_{k}$ and generalized conductivities $\sigma_{k}=\sigma_{k}^{\prime}+i\sigma_{k}^{\prime\prime}$ are complex. We denote real and imaginary parts of complex quantities with prime and double prime, respectively, and keep using this notation throughout the paper.

The model of the permittivity $\varepsilon(\omega)$ given by the analytic function Eq. (2) with $\Omega^{\prime\prime}_{k}\leqslant 0$ respects the constrain of causality by construction. The parameters of the model, which are the conductivities and the resonance frequencies, have to be determined from the experimentally measured data.

III Optimization

With the analytic model Eq. (2) of the permittivity, the task of fitting the experimental data reduces to finding the parameters of the model which minimize the error weighted deviation $E$ between the analytic and the measured values of $\varepsilon$, as this maximizes the probability of the model given the data. Assuming Gaussian errors, we use the squared deviation, weighted with the root-mean-square (RMS) errors:

[TABLE]

where $\varepsilon_{j}$ are experimental values, $\Delta\varepsilon_{j}$ are the corresponding errors, and $N$ is the total number of data points taken into account in the fit. Considering that typical experimental data sets consist of tens to hundreds of points, and $\varepsilon(\omega)$ is an analytic function of $\omega$ with a large number of parameters, typically in the order of a few tens, a robust and efficient algorithm is needed. To achieve this goal, we first make use of an exact, analytical minimization with respect to the parameters in which $\varepsilon$ is linear – these are all the conductivities and $\varepsilon_{\infty}$. This is the reason why it is advantageous to fit $\varepsilon$ instead of the complex refractive index $n+i\kappa$, as for the latter none of the parameters is linear. Then for the rest of the parameters, in which $\varepsilon$ is nonlinear – these are the pole frequencies – we use an iterative minimization with a gradient decent and a suited selection of starting points.

To make this linear dependence on parameters more clear, we write the permittivity as

[TABLE]

in which one can see $1+2L$ linear real parameters: $\varepsilon_{\infty}$, $\sigma_{k}^{\prime}$ and $\sigma_{k}^{\prime\prime}$. Minimization of the total error $E$, given by Eq. (3), with respect to the linear parameters can be done analytically by setting to zero all the first derivatives of $E$ with respect to these parameters. This provides a set of $1+2L$ linear equations which can be solved numerically exactly using standard linear algebra packages. We can fix the value of $\varepsilon_{\infty}$ if requested, removing it from the set of linear parameters, by subtracting our chosen value $\varepsilon_{\infty}$ from $\varepsilon(\omega_{j})$.

Using the values of the linear parameters determined by exact minimization of $E$, we now define, via Eq. (3), a new error function $\widetilde{E}$, which has been already minimized with respect to the linear parameters and depends only on the nonlinear parameters, which are the complex frequencies $\Omega_{k}$ of the Lorentz poles. Overall, there are $2L$ remaining real parameters over which $\widetilde{E}$ has to be minimized. To represent the average deviation of the model from the measured data points relative to their experimental RMS error, we introduce the relative error of our fit,

[TABLE]

To minimize $\tilde{E}$ over the $2L$ nonlinear parameters we use known minimization algorithms based on the gradient decent (implemented in MATLAB as function ‘fminunc’). The main challenge is to select suited starting points for the parameters, from which the algorithm finds local minima. The starting points should be selected in a way that the global minimum is amongst the local minima found. Remaining abrupt changes of $S$ with $N$ can, but do not have to, indicate that the global minimum was not yet obtained, and more starting values should be employed.

For $L=1$, we have a pair of Lorentz poles given by a single complex parameter $\Omega_{1}$. For the starting value of $\Omega_{1}$, we use a random logarithmic distribution within the range of the measured data, specifically

[TABLE]

where $Y$ and $Y^{\prime}$ are random numbers with a uniform distribution between 0 and 1. The minimization is repeated with different starting points until at least ten resulting $S$ values are equal within 10%, and the parameters for the lowest $S$ are accepted as global minimum.

The parameter space volume to be covered in such a procedure increases exponentially with $L$, making it computationally prohibitive to use this approach for large $L$. Increasing $L$, we therefore revert to a different strategy. Instead of guessing all $\Omega_{k}$ randomly, we use the optimized values for $\Omega_{1},...,\Omega_{L-1}$ of the model with $L-1$ poles as starting values for the simulation for $L$ poles, and choose the starting value for the additional pole as $\Omega_{L}=[1-{i}/{(L+1)}]\omega_{N}$. This method is fast but can result in not finding the global minimum. However, we can vary the range of the experimental data to be fitted in order to provide more starting points. Here, we choose to keep the lowest frequency $\omega_{1}$ fixed but vary $\omega_{N}$ and consequently $N$. Increasing or decreasing $N$ by one, we use as starting point the optimized values for $N$.

Furthermore, going back, from $L+1$ to $L$, just removing one pair of Lorentz poles provides $L+1$ additional starting values for the simulation with $L$ poles. It is also possible to go back multiple steps, e.g., from $L+2$ to $L$ provides $(L+2)(L+1)/2$ starting values – this however has not been used to produce the $S$ values in this paper.

IV Results

Here we discuss examples of the Lorentz model optimized for measured material dispersions. We show results for crystals of silicon (Si), gallium arsenide (GaAs) and germanium (Ge) using data from [Sopra, ]. In our previous work we considered metals (gold, silver and copper) and showed the full range of the experimental permittivity fitted with one Drude pole and three or four pairs of Lorentz poles giving a satisfactory fit. As semiconductors do not have a significant free carrier density, they do not require a Drude pole, so we have fitted the experimental permittivity with the background term ($\varepsilon_{\infty}$) and up to five pairs of Lorentz poles. We use a variable optimization range, from the lowest measured photon energy $\hbar\omega_{1}$ to a variable upper boundary of the photon energy $\hbar\omega_{N}$ taking all available measured values. We show the resulting $S$ values for Si in Fig. 1(b) for different numbers of poles taken into account. We can see that with an increasing number of poles, the error is decreasing, as expected considering the increasing number of parameters. Also, increasing $\hbar\omega_{N}$ results in larger $S$ values, since a model of a given number of parameters is used to describe an increasing number of data.

As we have no experimental errors for the permittivity we use in this work $\Delta\varepsilon_{j}=\varepsilon_{j}$, which means that the $S$ values we use are the normalized relative error. We show a dashed line in Fig. 1(b) at 2% relative error as a guide to a satisfactory fit.

Concerning the relation of the poles to interband transitions in solids, it is important to emphasize that in microscopic theory the optical response is due to a large number of transitions, often described by a continuum. This continuum, however, can be represented by an infinite or a finite number of poles of the self-energy describing the effects of screening and frequency dispersion. Therefore, the model with a limited number of Lorentz oscillators presents a fully physical though approximate approach, collecting the oscillator strength and transition energies of the continuum into a finite number of poles. The resulting pole positions and weights depend on the energy range to be described and represent microscopic transitions in the material.

For Si, adding the first pair of Lorentz poles ($L=1$) we see the effect of the interband transitions can be described up to the first peak in $\varepsilon^{\prime\prime}$, around 3.4 eV. Adding the second pair of Lorentz poles ($L=2$), the effect of the interband transitions can be described up to about 4 eV, as $\hbar\omega_{N}$ approaches another peak in $\varepsilon^{\prime\prime}$. Three pairs ($L=3$) describe both peaks in $\varepsilon^{\prime\prime}$ up to 4.7 eV and finally four pairs of Lorentz poles adequately describes the full range of measured data up to 7.6 eV, with $S=0.0162$. We do see a significant improvement in $S$ when going to five pairs of poles.

The optimized model for Si with $L=5$ is compared with the experimental data in Fig. 2(a). The refraction and absorption indices are shown as a function of the photon energy, with the measured data as circles and lines representing the fit functions of the Lorentz model. The poles of the model [see Eq.(2)] are shown as circles in Fig. 2(b), centered at their pole positions $\Omega_{j}$ in the complex photon energy plane, with the complex pole weight represented by the circle area proportional to $\sqrt{|\sigma_{j}|}$ and the color giving the phase. We find $S=0.0102$ for this fit, with other parameters given in Table 1. We can see that the Lorentz poles are properly positioned to model the interband transitions of silicon. The phases of all 5 poles are close to $\pi/2$, corresponding to a classical damped Lorentz oscillator, such as a mass on a spring. The resonances $\Omega_{1}^{\prime}\sim 3.4$ eV and $\Omega_{3}^{\prime}\sim 4.3$ eV are around the centers of the two interband transitions well seen in Fig1.(a) within the optimization range, and the half-width of the resonances, $-\Omega_{1}^{\prime\prime}\sim 0.1$ eV and $-\Omega_{3}^{\prime\prime}\sim 0.2$ eV, are approximately covering the half-width of these transitions.

We obtained similar results for GaAs and Ge which can be seen in Fig. 3 and Fig. 5, respectively. As in Fig. 1, we use a variable upper limit $\hbar\omega_{N}$ of the optimization range and show the resulting $S$ values for different numbers of poles. For GaAs we find that the model with one pair of Lorentz oscillators works well up to 1.6 eV, two up to 2.9 eV, three up to 4.6 eV, four up to 5 eV, and five up to a value above the upper limit of 6 eV. For Ge we find that the one pair of Lorentz oscillators is a good approximation up to 1.1 eV, two up to 2.2 eV, three up to 4.2 eV, four up to 5.6 eV, and five beyond the upper limit of 6 eV.

In Fig. 4, the first Lorentz pole conductivity $\sigma_{1}$ is multiplied by a factor of 10 for clarity. We find phases of all poles close to $\pi/2$, corresponding to a classical damped Lorentz oscillator. There are similar results for Ge which can be seen in Fig. 6, most poles having a phase close to $\pi/2$ with the exception of the first pole which has a phase very close to zero. We note that the indirect band gap of both Si and Ge is less suited for the modeling by simple poles, due to the phonon-assisted absorption leading to a weak tail in the absorption spectrum. We can see this clearly in the absorption index for $\hbar\omega<1$ eV, shown in Fig. 6(a).

All parameters and values of $S$ for the fits shown in Fig. 2, Fig. 4 and Fig. 6 are given in Table 1.

V Conclusion

In conclusion we have demonstrated the performance of our optimization algorithm to determine the parameters of a generalized Drude-Lorentz model for the permittivity of semiconductors, which can be modeled with Lorentz poles only. For $L$ pairs of Lorentz poles taken into account, the developed algorithm uses an analytic minimization over the $2L+1$ linear parameters of the model (the generalized conductivities and high frequency value $\varepsilon_{\infty}$), and a gradient decent method for determining the $2L$ nonlinear parameters of the model (the Lorentz pole frequencies), with a suited choice of the starting values, resulting in fast and reliable determination of the best global fit. Examples of the fit using up to 5 pairs of Lorentz poles are provided for Si, GaAs and Ge. The optimization program implementing the described algorithm to model any measured data for the refractive index or permittivity is also provided [Program, ].

Acknowledgements.

This work was supported by the Sêr Cymru National Research Network in Advanced Engineering and Materials.