# Materials in Meep

(Difference between revisions)
 Revision as of 01:13, 20 July 2008 (edit)Stevenj (Talk | contribs) (Dielectric materials in Meep moved to Materials in Meep: more generic name)← Previous diff Revision as of 01:47, 20 July 2008 (edit)Stevenj (Talk | contribs) (new features)Next diff → Line 1: Line 1: {{Meep}} {{Meep}} - The material structure in Maxwell's equations is determined by the dielectric function ε('''x'''), but ε is not only a function of position — in general, it also depends on frequency (''material dispersion'') and on the electric field '''E''' itself (''nonlinearity''). Material dispersion, in turn, is generally associated with absorption loss in the material, or possibly ''gain''. All of these effects can be simulated in Meep, with certain restrictions. + The material structure in Maxwell's equations is determined by the relative permittivity ε('''x'''), but ε is not only a function of position — in general, it also depends on frequency (''material dispersion'') and on the electric field '''E''' itself (''nonlinearity''). Material dispersion, in turn, is generally associated with absorption loss in the material, or possibly ''gain''. All of these effects can be simulated in Meep, with certain restrictions. + + Similarly for the relative permeability μ('''x'''), although the current version of Meep does not support arbitrary dispersion in μ, nor does it support nonlinearity in μ. In this section, we describe the form of the equations and material properties that Meep can simulate. The actual interface with which you specify these properties is described in the [[Meep reference]]. In this section, we describe the form of the equations and material properties that Meep can simulate. The actual interface with which you specify these properties is described in the [[Meep reference]]. Line 11: Line 13: :$\mathbf{D} = \varepsilon_\infty \mathbf{E} + \mathbf{P}$ :$\mathbf{D} = \varepsilon_\infty \mathbf{E} + \mathbf{P}$ - where $\varepsilon_\infty$ is the ''instantaneous'' dielectric function (the infinite-frequency response) and '''P''' is the ''polarization'' density in the material. '''P''', in turn, has its ''own'' time-evolution equation, and the exact form of this equation determines the frequency-dependence ε(ω). In particular, Meep supports any material dispersion of the form of a sum of harmonic resonances: + where $\varepsilon_\infty$ is the ''instantaneous'' dielectric function (the infinite-frequency response) and '''P''' is the ''polarization'' density in the material. '''P''', in turn, has its ''own'' time-evolution equation, and the exact form of this equation determines the frequency-dependence ε(ω). In particular, Meep supports any material dispersion of the form of a sum of harmonic resonances, plus terms from a frequency-independent electric conductivity $\sigma_D$: - :$\varepsilon(\omega,\mathbf{x}) = \varepsilon_\infty(\mathbf{x}) + \sum_n \frac{\sigma_n(\mathbf{x}) \cdot \omega_n^2 \Delta\varepsilon_n}{\omega_n^2 - \omega^2 - i\omega\gamma_n} ,$ + :$\varepsilon(\omega,\mathbf{x}) = \varepsilon_\infty(\mathbf{x}) + \frac{i \varepsilon_\infty(\mathbf{x}) \cdot \sigma_D(\mathbf{x})}{\omega} + \sum_n \frac{\sigma_n(\mathbf{x}) \cdot \omega_n^2 \Delta\varepsilon_n}{\omega_n^2 - \omega^2 - i\omega\gamma_n} ,$ where $\omega_n$, $\gamma_n$, and $\Delta\varepsilon_n$ are user-specified constants and $\sigma_n(\mathbf{x})$ is a user-specified function of position (usually 0 or 1). This corresponds to evolving '''P''' via the equations: where $\omega_n$, $\gamma_n$, and $\Delta\varepsilon_n$ are user-specified constants and $\sigma_n(\mathbf{x})$ is a user-specified function of position (usually 0 or 1). This corresponds to evolving '''P''' via the equations: Line 22: Line 24: That is, we must store and evolve a set of auxiliary fields $\mathbf{P}_n$ along with the electric field in order to keep track of the polarization '''P'''. Essentially any ε(ω) could be modeled by including enough of these polarization fields — Meep allows you to specify any number of these, limited only by computer memory and time (which must increase with the number of polarization terms you require). That is, we must store and evolve a set of auxiliary fields $\mathbf{P}_n$ along with the electric field in order to keep track of the polarization '''P'''. Essentially any ε(ω) could be modeled by including enough of these polarization fields — Meep allows you to specify any number of these, limited only by computer memory and time (which must increase with the number of polarization terms you require). - To implement a ''Drude model'' of $\varepsilon(\omega)$, in which $\omega_n$ is zero in the denominator, we just set $\omega_n to be a very small number (e.g. 1e-20) and make [itex]\Delta\varepsilon_n$ large to compensate in the numerator. + Note that the conductivity $\sigma_D$ corresponds to an imaginary part of ε given by $i \varepsilon_\infty \sigma_D / \omega$. When you specify frequency in Meep units, however, you are specifying ''f'' without the 2π, so the imaginary part of ε is $i \varepsilon_\infty \sigma_D / 2\pi f$. == Loss and gain == == Loss and gain == - If γ above is nonzero, then the dielectric function ε(ω) becomes ''complex'', where the imaginary part is associated with absorption loss in the material if it is positive, or gain if it is negative. + If γ above is nonzero, then the dielectric function ε(ω) becomes ''complex'', where the imaginary part is associated with absorption loss in the material if it is positive, or gain if it is negative. Alternatively, a dissipation loss or gain may be added by a positive or negative conductivity, respectively. - If you look at Maxwell's equations, then $d\mathbf{P}/dt$ plays exactly the same role as a current $\mathbf{J}$. Just as $\mathbf{J} \cdot \mathbf{E}$ is the rate of change of mechanical energy (the power expended by the electric field on moving the currents), therefore, the rate at which energy is lost to absorption is given by: + If you look at Maxwell's equations, then $d\mathbf{P}/dt$ plays exactly the same role as a current $\mathbf{J}$. Just as $\mathbf{J} \cdot \mathbf{E}$ is the rate of change of mechanical energy (the power expended by the electric field on moving the currents), therefore, the rate at which energy is lost to absorption for is given by: :absorption rate $\sim \frac{d\mathbf{P}}{dt} \cdot \mathbf{E}$ :absorption rate $\sim \frac{d\mathbf{P}}{dt} \cdot \mathbf{E}$ - Meep can keep track of this energy, which for gain gives the amount of energy expended in amplifying the field. Using this energy, Meep supports the idea of a ''saturable gain'' (e.g. a situation in which there is a depletable population inversion causing the gain). For more information, see [[saturable gain in Meep]]. + Meep can keep track of this energy for the Lorentzian polarizability terms (but not for the conductivity terms), which for gain gives the amount of energy expended in amplifying the field. Using this energy, Meep supports the idea of a ''saturable gain'' (e.g. a situation in which there is a depletable population inversion causing the gain). For more information, see [[saturable gain in Meep]]. + + == Conductivity == + + Often, you only care about the absorption loss in a narrow bandwidth, where you just want to set the imaginary part of ε (or μ) to some known experimental value, in the same way that you often just care about setting a dispersionless real ε that is the correct value in your bandwidth of interest. + + Some FDTD packages approach this problem by allowing you to specify a constant (frequency-independent) imaginary part of ε, but this has the disadvantage of requiring the simulation to employ complex fields (doubling the memory and time requirements). Instead, the approach recommended in Meep is for you to set the conductivity $\sigma_D$ (or $\sigma_B$ for an imaginary part of μ), chosen so that $\Im \varepsilon = i \varepsilon_\infty \sigma_D / \omega$ is the correct value at your frequency ω of interest. (Note that, in Meep, you specify $f = \omega/2\pi$ instead of ω for the frequency, however, so you need to include the factor of 2π when computing the corresponding imaginary part of ε!) (Conductivities can be implemented with purely real fields, so they are not nearly as expensive as implementing a frequency-independent complex ε or μ.) + + However, note that the "conductivity" in Meep is slightly different from the conductivity you might find in a textbook, because (for computational convenience) it appears as $\sigma_D \mathbf{D}$ in our Maxwell equations rather than the more-conventional $\sigma \mathbf{E}$; this just means that our definition is different from the usual electric conductivity by a factor of ε. Also, just as Meep uses the dimensionless relative permittivity for ε, it uses nondimensionalized units of 1/''a'' (where ''a'' is your unit of distance) for the conductivities $\sigma_{D,B}$. If you have the electric conductivity σ in SI units and want to convert to $\sigma_D$ in Meep units, you can simply use the formula: $\sigma_D = (a/c) \sigma / \varepsilon_r \varepsilon_0$ (where ''a'' is your unit of distance in meters, ''c'' is the vacuum speed of light in m/s, $\varepsilon_0$ is the SI vacuum permittivity, and $\varepsilon_r$ is the relative permittivity). == Nonlinearity == == Nonlinearity ==

## Revision as of 01:47, 20 July 2008

The material structure in Maxwell's equations is determined by the relative permittivity ε(x), but ε is not only a function of position — in general, it also depends on frequency (material dispersion) and on the electric field E itself (nonlinearity). Material dispersion, in turn, is generally associated with absorption loss in the material, or possibly gain. All of these effects can be simulated in Meep, with certain restrictions.

Similarly for the relative permeability μ(x), although the current version of Meep does not support arbitrary dispersion in μ, nor does it support nonlinearity in μ.

In this section, we describe the form of the equations and material properties that Meep can simulate. The actual interface with which you specify these properties is described in the Meep reference.

## Material dispersion

Physically, material dispersion arises because the polarization of the material does not respond instantaneously to an applied field E, and this is essentially the way that it is implemented in FDTD. In particular, $\mathbf{D} = \varepsilon\mathbf{E}$ is expanded to:

$\mathbf{D} = \varepsilon_\infty \mathbf{E} + \mathbf{P}$

where $\varepsilon_\infty$ is the instantaneous dielectric function (the infinite-frequency response) and P is the polarization density in the material. P, in turn, has its own time-evolution equation, and the exact form of this equation determines the frequency-dependence ε(ω). In particular, Meep supports any material dispersion of the form of a sum of harmonic resonances, plus terms from a frequency-independent electric conductivity σD:

$\varepsilon(\omega,\mathbf{x}) = \varepsilon_\infty(\mathbf{x}) + \frac{i \varepsilon_\infty(\mathbf{x}) \cdot \sigma_D(\mathbf{x})}{\omega} + \sum_n \frac{\sigma_n(\mathbf{x}) \cdot \omega_n^2 \Delta\varepsilon_n}{\omega_n^2 - \omega^2 - i\omega\gamma_n} ,$

where ωn, γn, and $\Delta\varepsilon_n$ are user-specified constants and $\sigma_n(\mathbf{x})$ is a user-specified function of position (usually 0 or 1). This corresponds to evolving P via the equations:

$\mathbf{P} = \sum_n \mathbf{P_n}$
$\frac{d^2\mathbf{P}_n}{dt^2} + \gamma_n \frac{d\mathbf{P}_n}{dt} + \omega_n^2 \mathbf{P}_n = \sigma_n(\mathbf{x}) \omega_n^2 \Delta\varepsilon_n \mathbf{E}$

That is, we must store and evolve a set of auxiliary fields $\mathbf{P}_n$ along with the electric field in order to keep track of the polarization P. Essentially any ε(ω) could be modeled by including enough of these polarization fields — Meep allows you to specify any number of these, limited only by computer memory and time (which must increase with the number of polarization terms you require).

Note that the conductivity σD corresponds to an imaginary part of ε given by $i \varepsilon_\infty \sigma_D / \omega$. When you specify frequency in Meep units, however, you are specifying f without the 2π, so the imaginary part of ε is $i \varepsilon_\infty \sigma_D / 2\pi f$.

## Loss and gain

If γ above is nonzero, then the dielectric function ε(ω) becomes complex, where the imaginary part is associated with absorption loss in the material if it is positive, or gain if it is negative. Alternatively, a dissipation loss or gain may be added by a positive or negative conductivity, respectively.

If you look at Maxwell's equations, then $d\mathbf{P}/dt$ plays exactly the same role as a current $\mathbf{J}$. Just as $\mathbf{J} \cdot \mathbf{E}$ is the rate of change of mechanical energy (the power expended by the electric field on moving the currents), therefore, the rate at which energy is lost to absorption for is given by:

absorption rate $\sim \frac{d\mathbf{P}}{dt} \cdot \mathbf{E}$

Meep can keep track of this energy for the Lorentzian polarizability terms (but not for the conductivity terms), which for gain gives the amount of energy expended in amplifying the field. Using this energy, Meep supports the idea of a saturable gain (e.g. a situation in which there is a depletable population inversion causing the gain). For more information, see saturable gain in Meep.

## Conductivity

Often, you only care about the absorption loss in a narrow bandwidth, where you just want to set the imaginary part of ε (or μ) to some known experimental value, in the same way that you often just care about setting a dispersionless real ε that is the correct value in your bandwidth of interest.

Some FDTD packages approach this problem by allowing you to specify a constant (frequency-independent) imaginary part of ε, but this has the disadvantage of requiring the simulation to employ complex fields (doubling the memory and time requirements). Instead, the approach recommended in Meep is for you to set the conductivity σD (or σB for an imaginary part of μ), chosen so that $\Im \varepsilon = i \varepsilon_\infty \sigma_D / \omega$ is the correct value at your frequency ω of interest. (Note that, in Meep, you specify f = ω / 2π instead of ω for the frequency, however, so you need to include the factor of 2π when computing the corresponding imaginary part of ε!) (Conductivities can be implemented with purely real fields, so they are not nearly as expensive as implementing a frequency-independent complex ε or μ.)

However, note that the "conductivity" in Meep is slightly different from the conductivity you might find in a textbook, because (for computational convenience) it appears as $\sigma_D \mathbf{D}$ in our Maxwell equations rather than the more-conventional $\sigma \mathbf{E}$; this just means that our definition is different from the usual electric conductivity by a factor of ε. Also, just as Meep uses the dimensionless relative permittivity for ε, it uses nondimensionalized units of 1/a (where a is your unit of distance) for the conductivities σD,B. If you have the electric conductivity σ in SI units and want to convert to σD in Meep units, you can simply use the formula: $\sigma_D = (a/c) \sigma / \varepsilon_r \varepsilon_0$ (where a is your unit of distance in meters, c is the vacuum speed of light in m/s, $\varepsilon_0$ is the SI vacuum permittivity, and $\varepsilon_r$ is the relative permittivity).

## Nonlinearity

In general, ε can be changed anisotropically by the E field itself, with:

$\Delta\varepsilon_{ij} = \sum_{k} \chi_{ijk}^{(2)} E_k + \sum_{k\ell} \chi_{ijk\ell}^{(3)} E_k E_\ell + \cdots$

where the ij is the index of the change in the 3×3 ε tensor and the χ terms are the nonlinear susceptibilities. The χ(2) sum is the Pockels effect and the χ(3) sum is the Kerr effect. (If the above expansion is frequency-independent, then the nonlinearity is instantaneous; more generally, Δε would depend on some average of the fields at previous times.)

Currently, Meep supports instantaneous, isotropic Pockels and Kerr nonlinearities, corresponding to a frequency-independent $\chi_{ijk}^{(2)} = \chi^{(2)} \cdot \delta_{ij} \delta_{jk}$ and $\chi_{ijk\ell}^{(3)} = \chi^{(3)} \cdot \delta_{ij} \delta_{k\ell}$, respectively. Thus,

$\mathbf{D} = \left( \varepsilon_\infty(\mathbf{x}) + \chi^{(2)}(\mathbf{x})\cdot \mathrm{diag}(\mathbf{E}) + \chi^{(3)}(\mathbf{x}) \cdot |\mathbf{E}|^2 \right) \mathbf{E} + \mathbf{P}$

Here, "diag(E)" indicates the diagonal 3×3 matrix with the components of E along the diagonal.

Normally, for nonlinear systems you will want to use real fields E. (This is usually the default. However, Meep uses complex fields if you have Bloch-periodic boundary conditions with a non-zero Bloch wavevector k, or in cylindrical coordinates with $m \neq 0$. In the C++ interface, real fields must be explicitly specified.)

For complex fields in nonlinear systems, the physical interpretration of the above equations is unclear because one cannot simply obtain the physical solution by taking the real part any more. In particular, Meep simply defines the meaning of the nonlinearity for complex fields as follows: the real and imaginary parts of the fields do not interact nonlinearly. That is, the above equation should be taken to hold for the real and imaginary parts (of E and D) separately. (e.g. |E|2 is the squared magnitude of the real part of E for when computing the real part of D, and conversely for the imaginary part.)

Note: The behavior for complex fields was changed for Meep 0.10. Also, in Meep 0.9 there was a bug: when you specified χ(3) in the interface, you were actually specifying $\chi^{(3)}/\varepsilon_\infty^4$. This was fixed in Meep 0.10.

For a discussion of how to relate χ(3) in Meep to experimental Kerr coefficients, see Units and nonlinearity in Meep.