Question about the Formulation of Drude/Lorentz Model by Taflove and by openEMS

[biergaizi]

While trying to improve documentation of openEMS, I'm trying to understand how openEMS's Drude/Lorentz model is formulated, in comparison to the textbook formulations. It seems that these formulations are not exactly equivalent, even after changing the symbols and notations to the same form.

I'm not sure if I understand it correctly, here's my question.

Drude Model

Naive Model

According to Allen Taflove's Computational Electrodynamics The Finite-Difference Time-Domain Method, Third Edition, page 355, equation 9.7 [1], the first-order Drude model defines a material's electric susceptibility as:

$$ \chi(\omega) = - \frac{\omega_p^2}{\omega^2 \ - j \omega \gamma_p} $$

To be consistent with openEMS's notation, substitute the variable $ v_e = 1 / \gamma_p$.

$$ \chi(\omega) = - \frac{\omega_p^2}{\omega^2 - j \omega \frac{1}{\tau_p}} $$

The naive first-order relationship between susceptibility and permittivity, without the asymptotic and conductivity term, is defined as:

$$ \begin{align} \epsilon(\omega) &= \epsilon_0 \left[ 1 + \chi(\omega) \right] \\ &= \epsilon_0 \left[ 1 - \frac{\omega_p^2}{\omega^2 - j\omega \frac{1}{\tau_p}} \right] \end{align} $$

This equation is consistent with both Taflove's formulation in the book. This is also consistent with Rennings' uncorrected Drude model, as described in the paper Equivalent Circuit (EC) FDTD Method for Dispersive Materials: Derivation, Stability Criteria and Application Examples [2]. The paper gives the uncorrected model in Equation 58.

Asymptotic Models

At high frequencies, the electric susceptibility goes to zero. But real-world materials have a residue term, so we have to add an asymptotic term. But it gets confusing here, as there seems to be two slightly different ways of doing so.

Schneider's Formulation

The first formulation is to replace the constant $1$ term with the asymptotic term $\epsilon_r(\infty)$. In Understanding the FDTD Method by John B. Schneider, page 291, equation 10.8.

$$ \epsilon(\omega) = \epsilon_0 \left[ \epsilon_r(\infty) + \chi(\omega) \right] $$

So the new first-order Drude expression would be:

$$ \begin{align} \epsilon(\omega) &= \epsilon_0 \left[ \epsilon_r(\infty) + \chi(\omega) \right] \\ \epsilon(\omega) &= \epsilon_0 \left[ \epsilon_r(\infty) - \frac{\omega_p^2}{\omega^2 - j\omega \frac{1}{\tau_p}} \right] \\ &= \epsilon_0 \epsilon_r(\infty) - \epsilon_0 \frac{\omega_p^2}{\omega^2 - j\omega \frac{1}{\tau_p}} \end{align} $$

This is the standard expression used by multiple textbooks.

Rennings's Formulation

However, openEMS used Rennings's formulation, and the asymptotic term is applied in a different way. It's applied as a multiplication factor to the overall permittivity $\epsilon(\omega)$, not just the susceptibility $\chi(\omega)$:

$$ \begin{align} \epsilon(\omega) &= \epsilon_0 \epsilon_r(\infty) \left[ 1 + \chi(\omega) \right] \\ &= \epsilon_0 \epsilon_r(\infty) - \epsilon_0 \epsilon_r(\infty) \frac{\omega_p^2}{\omega^2 - j\omega \frac{1}{\tau_p}} \end{align} $$

So Schneider and Renning formulas have a difference by a factor of $\epsilon_r(\infty)$ in the second term.

Lorentz Model

More differences arise if the Drude model is generalized to the Lorentz model, which is described here.

Taflove's Formulation

Allen Taflove's book states that the Lorentz model in page 354, equation 9.4:

$$ \chi(\omega) = \frac{\left[ \epsilon_r(\mathrm{static}) - \epsilon_r(\infty) \right] \omega^2}{\omega_p^2 + 2j \omega \delta_p - \omega^2} $$

Multiply the numerator and denominator by a factor of -1 and substitude the variable $2\delta_p$ with $1 / \tau_p$:

$$ \begin{align} \chi(\omega) &= \frac{-\left[\epsilon_r(\mathrm{static}) - \epsilon_r(\infty) \right] \omega^2}{\omega^2 -\omega_p^2 - j\omega \frac{1}{\tau_p}} \\ &= -\left[ \epsilon_r(\mathrm{static}) - \epsilon_r(\infty) \right] \frac{\omega^2}{\omega^2 -\omega_p^2 - j\omega \frac{1}{\tau_p}} \end{align} $$

Then we substitute it into the permittivity formula:

$$ \begin{align} \epsilon(\omega) &= \epsilon_0 \left[ \epsilon_r(\infty) + \chi(\omega) \right] \\ &= \epsilon_0 \epsilon_r(\infty) -\epsilon_0 \left[ \epsilon_r(\mathrm{static}) - \epsilon_r(\infty) \right] \frac{\omega^2}{\omega^2 -\omega_0^2 - j\omega \frac{1}{\tau_p}} \\ &= \epsilon_0 \epsilon_r(\infty) -\epsilon_0 \Delta \epsilon_r \frac{\omega^2}{\omega^2 -\omega_0^2 - j\omega \frac{1}{\tau_p}} \end{align} $$

This result is found in most textbooks.

Rennings's Formulation

In Rennings's paper, openEMS generalizes the corrected Drude model to the Lorentz model simply by adding an extra Lorentz pole term.

$$ \epsilon(\omega) = \epsilon_0 \epsilon_r(\infty) - \epsilon_0 \epsilon_r(\infty) \frac{\omega_p^2}{\omega^2 - \boxed{\omega^2_\mathrm{Lor}} - j \omega \frac{1}{\tau_p}} $$

Note that Rennings' formulation differs from the standard Taflove formulation. In Taflove's formulation, there are two asymptotic terms $\epsilon_r(\infty)$ and $\Delta \epsilon_r = \epsilon_r(\mathrm{static}) - \epsilon_r(\infty)$. In Rennings's formulation, there's only way asymptotic term $\epsilon_r(\infty)$.

So I think openEMS's Drude/Lorentz model is not exactly equivalent to the textbook Drude/Lorentz model. There's a different in the susceptibility term, so the coefficients in these two formulations will be different. The only way to make them equivalent to is absorb Taflove's $\Delta \epsilon_r$ term into openEMS's plasma frequency term $\omega_p^2$.

Thus, given the same input data, the fitted plasma frequencies $\omega_p$ in the openEMS model will be different from the conventional definition of plasma frequency in standard textbooks by a factor of $\sqrt{\Delta \epsilon_r} = \sqrt{\epsilon_r(\mathrm{static}) - \epsilon_r(\infty)}$. This is the correct meaning of relation to other formulations as documented in openEMS wiki.

Is it correct? If it's correct, I will include this derivation in the new documentation.

Or did I make any mistake?

Thanks for your time.

Bibliography

[1] Allen Taflove, Computational Electrodynamics The Finite-Difference Time-Domain Method, Third Edition.

[2] Andreas Rennings, et, al., Equivalent Circuit (EC) FDTD Method for Dispersive Materials: Derivation, Stability Criteria and Application Examples. https://www.researchgate.net/publication/227133697_Equivalent_Circuit_EC_FDTD_Method_for_Dispersive_Materials_Derivation_Stability_Criteria_and_Application_Examples

[3] John B. Schneider. Understanding the FDTD Method. https://eecs.wsu.edu/~schneidj/ufdtd/index.php

[thliebig]

I'm aware that there are many different formulations for Drude/Lorentz models. In openEMS I have followed closely what Andreas Rennings had in his PhD. He had all the equations, EC-FDTD equations and time-step stability analysis.
If I find the time and find the printed version I will have another look and maybe ask him if I can post some equations and such.

I think there is one matlab example in the openEMS tree showing a simple Lorentz simulation. Not sure if that helps any.

[biergaizi]

I'm aware that there are many different formulations for Drude/Lorentz models. In openEMS I have followed closely what Andreas Rennings had in his PhD.

To my best knowledge, my analysis above is correct, because it's consistent with the paper [2] published by Rennings and with the openEMS Octave source code, and with the wiki's documentation page. The only difference between Rennings's Lorentz model formulation and most textbook formulations is the difference of a factor $\sqrt{\Delta \epsilon_r} = \sqrt{\epsilon_r(\mathrm{static}) - \epsilon_r(\infty)}$. After the correct substitutions, all of these equations give numerically identical results. Unless there's evidence to the contrary, I'm going to include this note in the new Web documentation and submit a Pull Request as part of a large documentation update, with references to papers.

He had all the equations, EC-FDTD equations and time-step stability analysis.

The citation [2] I found also contains a published version of EC-FDTD equations, Drude/Lorentz models, and time-step stability analysis as part of a Springer book Time Domain Methods in Electrodynamics. This is how I finally understood why the asymptotic factor was chosen as described... It's likely not as detailed as the original thesis, but better than none.

[thliebig]

Yeah sound good.