Dispersive Materials
Nearly all real-world materials exhibit a phenomenon known as dispersion. That is, the speed of light in the medium depends on the EM wave’s frequency. In optics, it manifests as a frequency-dependent refractive index, causing the separation of sunlight when it passes through a prism. In RF/microwave engineering, it appears as a frequency-dependent permittivity and permeability, often causing unexpected changes in a filter’s resonant frequency, or a transmission line’s characteristic impedance. In high-speed digital data lines, dispersion may cause the smearing of signal edges and the broadening of signal pulses. In metamaterial research, one can even deliberately introduce dispersion to control electromagnetic wave propagation in unusual ways.
As a result, while the basic material model with constant permittivity and permeability is sufficient if dispersion is negligible, more demanding simulations call for dispersion models for accurately calculating a material’s wideband response. In fact, strictly speaking, neglecting this frequency dependence violates fundamental physical principles. The assumption of constant material properties implies instantaneous polarization or magnetization in a material, thus it violates the causality requirement as formalized in the Kramers-Kronig relations.
In openEMS, three dispersive material models known as Debye, Drude and Lorentz materials are used to simulate dispersion. These models were originally proposed by solid-state physicists in the early 20th century to explain the microscopic origins of electrical and optical properties of matter.
Due to their compatibility with Maxwell’s equations, FDTD borrows these models for simulating dispersive materials. Because these are analytical models, they do not directly take inputs in the form of a raw frequency-dependent curve or lookup table of material properties. Instead, a material’s permittivity and permeability must be parameterized using terms such as plasma frequency, plasma relaxation time, and Lorentz pole frequency.
Although these solid-state physics parameters may appear unrelated to practical RF/microwave engineering, in FDTD, they mainly act as calculation tools rather than providing physical interpretations. The required parameters are obtained by numerically fitting the model to measured material response curves for the purpose of simulations. Both first-order and higher-order models are supported, using one or more sets of parameters. There are many “degrees of freedom” in the fitting process. If a first-order fit fails to produce the desired curves, more parameters can be added until one obtains a satisfactory fit.
Note
Physics vs. Fitting. FDTD only borrows these solid-state physics models as calculation tools. The fitted parameters are usually purely numerical (i.e. they fit measured data), or phenomenological at best (i.e. they match experimental outcomes without explaining the underlying mechanism). As such, they’re generally not physically meaningful. As von Neumann famously joked, “With four parameters I can fit an elephant, and with five I can make him wiggle his trunk.”
Python. The functions described below are not ported to Python yet.
Complex Permittivity and Permeability
In circuit design, the permittivity of a dielectric material or the permeability of a magnetic material is often given as a real-valued constant. However, this is an incomplete definition that is only suitable for DC or narrowband applications. For wideband EM simulation with dispersion and loss, one needs the full definition.
The complex permittivity and permeability are functions of frequency:
In dielectric materials, the complex permittivity function describes the polarization and associated energy loss. In magnetic materials, the complex permeability describes the magnetization and associated energy loss.
The effects of both dielectric and conductive (ohmic) losses are included in the imaginary parts in the complex representation. Thus, the electric and magnetic loss tangents can be defined as:
Debye Model
Usage
AddDebyeMaterial()
creates a material governed by the Debye model of dielectric
relaxation:
CSX = AddDebyeMaterial(CSX, name, varargin)
The model is defined by the following parameters:
eps_r
: Relative permittivity when frequency approaches infinity (\(\epsilon_{r,\infty}\)).EpsilonDelta_n
: n-th delta dielectric permittivity, a fitted phenomenological term that represents relaxation strength (\(\Delta \epsilon_r\)).EpsilonRelaxTime_n
: n-th relaxation time, inverse of the damping factor, a dissipative term (\(\tau_\mathrm{relax}\)).Kappa
: Electric conductivity, an ohmic loss term (\(\kappa\)). Not to be confused with electric permittivity, sometimes also denoted by the same letter.
Warning
Kappa
(\(\kappa\)) always stands for electric conductivity in openEMS.
It’s not to be confused with electric permittivity \(\epsilon\), which is
sometimes also denoted as \(\kappa\) in the literature (e.g. high-κ
dielectric). This convention is never used in openEMS, its use in simulation
code is strongly discouraged.
For the first-order model, use one set of parameters (i.e. EpsilonDelta_1
,
EpsilonRelaxTime_1
). For higher-order modeling, add more parameters with
increasing indices in their suffixes n
(e.g. EpsilonDelta_2
, EpsilonRelaxTime_2
,
…)
Example:
CSX = AddDebyeMaterial(CSX, 'debye_example');
CSX = SetMaterialProperty(CSX, 'debye_example', 'Epsilon', 5, 'EpsilonDelta_1', 0.1, 'EpsilonRelaxTime_1',1e-9);
% create geometry
CSX = AddBox(CSX, 'debye_example', 10, start, stop);
Calculation
CalcDebyeMaterial()
is a helper function to calculate the numerical values of
permittivity at specified frequency points in an array. It’s useful for fitting
material parameters for simulation preparation as part of a numerical optimization
routine.
Definition:
eps_debye = CalcDebyeMaterial(f, eps_r, kappa, eps_Delta, t_relax)
Parameters:
f
: (vector) Frequency points of interest.eps_r
: Relative permittivity when frequency approaches infinity (\(\epsilon_{r,\infty}\)).kappa
: Electric conductivity, an ohmic loss term (\(\kappa\)).eps_Delta
: (vector) Delta relative permittivity, a fitted phenomenological term that represents oscillator strength (\(\Delta \epsilon_r\))t_relax
: (vector) Relaxation time, the inverse of damping factor, a dissipative term (\(\tau_\mathrm{relax}\)).
Drude Model
Usage
AddLorentzMaterial()
creates a material governed by either the Drude
model or the Lorentz model. The Drude model is obtained as a special case
of Lorentz when relevant parameters are omitted:
CSX = AddLorentzMaterial(CSX, name, varargin)
For dielectric materials, this model is defined by the following parameters:
Epsilon
: Relative permittivity when frequency approaches infinity (\(\epsilon_{r,\infty}\)).EpsilonPlasmaFrequency_n
: n-th electric plasma frequency in hertz (\(f_{p\epsilon}\)).EpsilonRelaxTime_n
: n-th electric plasma relaxation time (\(\tau_\epsilon\)), a dissipative term.Kappa
: Electric conductivity (\(\kappa\)), an ohmic loss term, optional but recommended to improve curve fitting. Not to be confused with electric permittivity, sometimes also denoted by the same letter.
Warning
Kappa
(\(\kappa\)) always stands for electric conductivity in openEMS.
It’s not to be confused with electric permittivity \(\epsilon\), which is
sometimes also denoted as \(\kappa\) in the literature (e.g. high-κ
dielectric). This convention is never used in openEMS, its use in simulation
code is strongly discouraged.
For magnetic materials, this model is defined by the following parameters:
Mue
: Relative permeability when frequency approaches infinity (\(\mu_{r,\infty}\)).MuePlasmaFrequency_n
: n-th magnetic plasma frequency in hertz (\(f_{p\mu}\)).MueRelaxTime_n
: n-th magnetic plasma relaxation time (\(\tau_\mu\)), a dissipative term.Sigma
: Magnetic conductivity (\(\sigma\)), an ohmic-like loss term. This is a non-physical property due to a hypothetical magnetic conduction current, optional but useful to improve fitting by introducing artificial losses.
Both the dielectric and magnetic parameters can be used simultaneously. If the material is not magnetic, magnetic parameters can be omitted.
The terms Epsilon
, Mue
, Kappa
and Sigma
are the effects of the
basic constant-property material model, and are not implemented in
the Drude/Lorentz model per se. But they should be considered in the curve
fitting process.
For the first-order model, use one set of parameters (i.e. EpsilonPlasmaFrequency
,
EpsilonRelaxTime
) without suffix _n
. For higher-order modeling, add more parameters
with increasing indices in their suffixes n (e.g. EpsilonPlasmaFrequency_1
,
EpsilonRelaxTime_1
, EpsilonPlasmaFrequency_2
, EpsilonRelaxTime_2
, …)
Example:
CSX = AddLorentzMaterial(CSX, 'drude_example');
CSX = SetMaterialProperty(CSX, 'drude_example', 'Epsilon', 5, 'EpsilonPlasmaFrequency', 5e9, 'EpsilonRelaxTime', 1e-9);
% optional, for magnetic materials
CSX = SetMaterialProperty(CSX, 'drude_example', 'Mue', 5, 'MuePlasmaFrequency', 5e9, 'MueRelaxTime', 1e-9);
% create geometry
CSX = AddBox(CSX, 'drude_example', 10, start, stop);
Calculation
CalcDrudeMaterial()
is a helper function to calculate the numerical values of
permittivity at specified frequency points in an array. It’s useful for fitting
material parameters for simulation preparation as part of a numerical optimization
routine.
Definition:
eps_drude = CalcDrudeMaterial(f, eps_r, kappa, plasmaFreq, t_relax)
Parameters:
f
: (vector) Frequency points of interest.eps_r
: Relative permittivity when frequency approaches infinity (\(\epsilon_{r,\infty}\)).kappa
: Electric conductivity, an ohmic loss term (\(\kappa\)).plasmaFreq
: (vector) Plasma frequencies in hertz (\(f_p\))t_relax
: (vector) Relaxation time, the inverse of damping factor, a dissipative term (\(\tau_\mathrm{relax}\)).
Since the Drude model can fit both dielectric and magnetic materials with the same formula, this function can be used to calculate permeability as well by symmetry - reinterpreting all inputs as permeability parameters.
Example:
% silver (AG) at optical frequencies (Drude model)
f = linspace(300e12, 1100e12, 201);
eps_model = CalcDrudeMaterial(f, 3.942, 7.97e3, 7e15/2/pi, 0, 1/2.3e13);
figure
plot(f,real(eps_model))
hold on;
grid on;
plot(f,imag(eps_model),'r--')
Lorentz Model
Usage
AddLorentzMaterial()
creates a material governed by the Lorentz
model:
CSX = AddLorentzMaterial(CSX, name, varargin)
The Lorentz model is defined by all Drude parameters as previously described, and two additional Lorentz parameters.
f_eps_Lor_Pole_n
: Electric Lorentz pole frequency in hertz (\(f_{\mathrm{Lor}\epsilon}\)).f_mue_Lor_Pole_n
: Magnetic Lorentz pole frequency in hertz (\(f_{\mathrm{Lor}\mu}\)).
If the material is non-magnetic, f_mue_Lor_Pole_n
is optional.
For the first-order model, use one set of parameters (i.e. f_eps_Lor_Pole
)
without suffix _n
. For higher-order modeling, add more parameters with increasing
indices in their suffixes n (e.g. f_mue_Lor_Pole_1
, f_mue_Lor_Pole_2
, …)
Calculation
CalcLorentzMaterial()
is a helper function to calculate the numerical values of
permittivity at specified frequency points in an array. It’s useful for fitting
material parameters for simulation preparation as part of a numerical optimization
routine.
Definition:
eps_lorentz = CalcLorentzMaterial(f, eps_r, kappa, plasmaFreq, LorPoleFreq, t_relax)
Parameters:
f
: (vector) Frequency points of interest.eps_r
: Relative permittivity when frequency approaches infinity (\(\epsilon_{r,\infty}\)).kappa
: Electric conductivity, an ohmic loss term (\(\kappa\)).plasmaFreq
: (vector) Plasma frequencies in hertz (\(f_p\))LorPoleFreq
: (vector) Lorentz pole frequencies in hertz (\(f_{\mathrm{Lor}\epsilon}\)). If zeroed, Lorentz model reduces to Drude model.t_relax
: (vector) Relaxation time, the inverse of damping factor, a dissipative term (\(\tau_\mathrm{relax}\)).
Since the Lorentz model can fit both dielectric and magnetic materials with the same formula, it can be used to calculate permeability as well by reinterpreting all inputs as permeability parameters.
Example:
% silver (AG) at optical frequencies (Drude+Lorentz model)
f = linspace(300e12, 1100e12, 201);
eps_model = CalcLorentzMaterial(f, 1.138, 4.04e3, [13e15 9.61e15]/2/pi, [0 7.5e15]/2/pi,[1/2.59e13 1/3e14]);
figure
plot(f,real(eps_model))
hold on;
grid on;
plot(f,imag(eps_model),'r--')
Math Notes
Debye Model
The Debye model of dielectric relaxation (not to be confused with Debye model of specific heat) is defined by the following equation ( 1, page 222, equation 37; 3, page 354, equation 9.3; 4, page 294, equation 10.27):
where:
\(n\) is the index of the higher-order model term being evaluated.
\(N\) is the number of model terms in total.
\(\omega = 2\pi f\) is the angular frequency.
\(\epsilon(\omega)\) is the material’s permittivity.
\(\epsilon_0\) is the vacuum permittivity.
\(\epsilon_{r,\infty}\) is the material’s relative permittivity at \(f \to \infty\).
\(\kappa\) is the material’s electric conductivity.
Note that electric conductivity \(\kappa\) is not to be confused with material permittivity, which is sometimes also denoted by \(\kappa\) in the literature.
\(\Delta \epsilon_{r}(n)\) is the n-th oscillator strength.
\(t_\mathrm{relax}(n)\) is the n-th relaxation time (inverse of the damping factor).
This model is only implemented for dielectric materials, not magnetic materials.
Drude Model
The Drude model of dielectric and magnetic materials is defined by the following equations ( 1, page 220, equation 27, 28):
where:
\(n\) is the index of the higher-order model term being evaluated.
\(N\) is the number of model terms in total.
\(\omega = 2\pi f\) is the angular frequency.
\(\omega_\mathrm{p\epsilon}(n)\) is the material’s n-th electric plasma frequency.
\(\tau_\epsilon(n)\) is the n-th electric relaxation time (inverse of the damping factor).
\(\epsilon(\omega)\) is the material’s permittivity.
\(\epsilon_0\) is the vacuum permittivity.
\(\epsilon_{r,\infty}\) is the material’s relative permittivity at \(f \to \infty\).
\(\kappa\) is the material’s electric conductivity.
Note that electric conductivity \(\kappa\) is not to be confused with material permittivity, which is sometimes also denoted by \(\kappa\) in the literature.
Analogously, the Drude model of magnetic materials is defined by a substitution of variables.
\(\omega_\mathrm{p\mu}(n)\) is the material’s n-th magnetic plasma frequency.
\(\tau_\mu(n)\) is the n-th magnetic relaxation time (inverse of the damping factor).
\(\mu(\omega)\) is the material’s permeability.
\(\mu_{r,\infty}\) is the material’s relative permeability at \(f \to \infty\).
\(\sigma\) is the material’s magnetic conductivity. It’s a non-physical property due to a hypothetical magnetic current. In most formulations, \(\sigma = 0\). In openEMS, this term is optional but useful to improve fitting by introducing artificial losses.
Lorentz Model
The Lorentz model of dielectric and magnetic materials is defined by the following equations:
where \(f_\mathrm{\mathrm{Lor}\epsilon}(n)\) and \(f_\mathrm{\mathrm{Lor}\mu}(n)\) are the material’s n-th electric and magnetic Lorentz pole frequencies. Zeroing these terms reduces the Lorentz model to the Drude model.
Other definitions are identical to the Drude model.
Relation to Other Formulations
The Drude/Lorentz formulation used by openEMS is not identical to the common formulations found in textbooks and the literature. In fact, many different variations are in use due to their flexibility. They can be formulated as first-order or higher-order models, with or without the conductivity term, with or without the asymptotic term, with different definitions of the asymptotic terms, and with different symbols and sign conventions.
Thus, it’s necessary to clarify the relationship between openEMS’s formulation with other textbooks to understand how these differences arise.
This section is not exhaustive due to the numerous possible combinations of the mentioned factors. However, after reading this section, readers will hopefully understand how to rewrite one formulation to another by simple algebra.
First-Order Lorentz Model
When there’s only one term, these models reduces to their first-order form. For example, the first-order Lorentz model of dielectric materials is defined by:
First-Order Drude Model
When the Lorentz pole frequency term \(\omega^2_\mathrm{Lor\mu} = 0\), the Lorentz model reduces to the Drude model. For example, the first-order Drude model of dielectric materials is defined by:
First-Order Drude Model w/o Asymptotic and Conductivity Terms
The openEMS’s Drude model contains a conductivity term \(\kappa\). For non-conductive material, the first-order Drude model reduces to:
This formulation is consistent with 3 page 355, equation 9.9; 4 page 292, equation 10.17.
At high frequencies, a material’s relative permittivity decreases, but never to exactly 1. This residue effect is captured by the asymptotic term \(\epsilon_{r,\infty}\). Without the asymptotic term, the Drude model reduces to:
Oscillator Strength Term Definitions
There’s a common parameter \(\Delta \epsilon_{r}\) that appears in all three models. It’s understood as the different between the electrostatic relative permittivity and the infinite-frequency permittivity ( 3, page 254, equation 9.2).
In solid-state physics, it has the physical meaning that ( 5, page 32, equation 2.19):
where \(N\) is the numeber of atoms per unit volume, \(e\) is the electron charge, \(\epsilon_0\) is the vacuum permittivity, \(m_0\) is the electron mass, \(\omega_0\) is the natural frequency of the atoms.
However, in material modeling and simulation rather than solid-state physics, \(\Delta \epsilon_{r}\) is used as a fitting parameter and does not represent physical processes at the atomic level.
Relaxation Time Definitions
Taflove’s formulation of the first-order Lorentz model is ( 3, page 355, equation 9.6):
where \(\epsilon_p\) is the frequency of the pole pair (the undamped resonant frequency of the medium), \(\delta_p\) is the damping coefficient. In comparison, openEMS absorbs the constant term 2 into the time constant \(\tau_p = \frac{1}{2\delta_p}\).
After algebraic manipulations, we obtain the following equation.
This is the standard textbook Lorentz model, rewritten with variable definitions similar to those used in openEMS. However, in addition to the time constant definition, textbooks also define the asymptotic term differently. See the next section.
Asymptotic Term Definitions
The formulation of Drude and Lorentz models are slightly different in comparison to the forms in standard textbooks, due to how the asymptotic term is defined. In most textbooks, the first-order Lorentz model (without the conductivity term) is defined as ( 3, page 355, equation 9.6; 4, page 293, equation 10.25):
But openEMS follows Rennings’s formulation and defines it as ( 1, page 235, equation 58, 59):
The textbook formulation contains two free parameters \(\Delta \epsilon_r\) and \(\epsilon_{r,\infty}\), but openEMS only uses \(\epsilon_{r,\infty}\). The \(\Delta \epsilon_r\) term is implicitly absorbed into \(\omega^{2}_{p\epsilon}\). To convert the textbook formulation to openEMS’s formulation, the plasma frequency must be redefined:
Sign Conventions
Researchers in physics often adopt a different sign convention, in which they use \(e^{-j\omega t}\) for time-harmonic quantities rather than \(e^{j \omega t}\) in engineering.
Example: Taflove’s Drude Model Formulation
The second edition of Taflove’s book defines the Drude model in physics convention, as ( 2, page 230, equation 9.9a):
In the third edition, it’s redefined in engineering convention, for consistency ( 3, page 354, equation 9.1):
These two equations are the same equation, where:
\(j' = -j\) due to the sign convention difference.
\(\gamma = 1/\tau_p\) is the damping factor, which is the inverse of pole relaxation time. It’s sometimes also denoted as the collision frequency of electrons \(\nu_e\) in the context of physics.
From the relationship between electric susceptibility and permittivity, without the asymptotic term:
This is our first-order Drude model without asymptotic and conductivity terms, which matching the description above. Introducing both additional terms according to the aforementioned procedure would reproduce the exact openEMS equations.
Example: Fox and Taflove’s Lorentz Model Formulation
The textbook Optical Properties of Solids by Mark Fox defines the Lorentz model as ( 5, page 36, equation 2.24):
where \(\omega_0\) is the plasmonic resonant frequency, \(f_j\) is a phenomenological oscillator strength parameter (used together with \(\Delta \epsilon_r\)).
This equation uses the opposite sign convention in comparison to most engineering books. By substituting \(j' = -j\), \(\gamma = 1/\tau_p\), \(f_j = \omega_p^2\) and multiplying by \(\frac{-1}{-1}\), we obtain:
This is our first-order Lorentz model without asymptotic and conductivity terms, which matching the description above. Introducing both additional terms according to the aforementioned procedure would reproduce the exact openEMS equations.
Bibliography
- 1(1,2,3)
Andreas Rennings, et, al., Equivalent Circuit (EC) FDTD Method for Dispersive Materials: Derivation, Stability Criteria and Application Examples, Time Domain Methods in Electrodynamics.
- 2
Allen Taflove, Susan. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 2nd ed. Artech House, 1995.
- 3(1,2,3,4,5,6)
Allen Taflove, Susan. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. Artech House, 2005.
- 4(1,2,3)
John B. Schneider. Understanding the FDTD Method.
- 5(1,2)
Mark Fox, Optical Properties of Solids. Oxford, UK: Oxford University Press, 2001.