minor improvements to documentation

doc/docs/C++_Developer_Information.md

@@ -2,7 +2,7 @@
 # C++ Developer Information
 ---
 
-An overview of Meep's inner workings is provided in [Computer Physics Communications, Vol. 181, pp. 687-702, 2010](http://ab-initio.mit.edu/~oskooi/papers/Oskooi10.pdf). This page is a supplement which provides a description of the source code.
+An overview of Meep's inner workings is provided in [Computer Physics Communications, Vol. 181, pp. 687-702 (2010)](http://ab-initio.mit.edu/~oskooi/papers/Oskooi10.pdf). This page is a supplement which provides a description of the source code.
 
 For additional details, see [Chunks and Symmetry](Chunks_and_Symmetry.md)
 

doc/docs/Chunks_and_Symmetry.md

@@ -2,7 +2,7 @@
 # Chunks and Symmetry
 ---
 
-As described in [Computer Physics Communications, Vol. 181, pp. 687-702, 2010](http://ab-initio.mit.edu/~oskooi/papers/Oskooi10.pdf),
+As described in [Computer Physics Communications, Vol. 181, pp. 687-702 (2010)](http://ab-initio.mit.edu/~oskooi/papers/Oskooi10.pdf),
 Meep subdivides geometries into **chunks**. Each chunk is a contiguous region of space — a line, rectangle, or parallelepiped for 1d/2d/3d Cartesian geometries, or an annular section in a cylindrical geometry—whose sizes are automatically determined by `libmeep.` In [parallel calculations](Parallel_Meep.md), each chunk is assigned, in its entirety, to precisely one process — that is, no chunk exists partly on one processor and partly on another.
 
 Many internal operations in Meep consist of looping over points in the
@@ -105,7 +105,7 @@ each own 4 chunks, while processes 1, 3, 4, and 6 each own 2 chunks.
 
 Meep's approach to handling symmetries is discussed
 from the user's perspective in [Exploiting Symmetry](Exploiting_Symmetry.md)
-and from a high-level algorithmic perspective in [Computer Physics Communications, Vol. 181, pp. 687-702, 2010](http://ab-initio.mit.edu/~oskooi/papers/Oskooi10.pdf).
+and from a high-level algorithmic perspective in [Computer Physics Communications, Vol. 181, pp. 687-702 (2010)](http://ab-initio.mit.edu/~oskooi/papers/Oskooi10.pdf).
 
 The following is a brief synopsis of the implementation of this feature.
 

doc/docs/Materials.md

@@ -90,7 +90,7 @@ Conductivity and Complex ε
 
 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.
 
-One approach to this problem would be allowing you to specify a constant, frequency-independent, imaginary part of ε, but this has the disadvantage of requiring the simulation to employ complex fields which double the memory and time requirements, and also tends to be numerically unstable. Instead, the approach in Meep is for you to set the conductivity $\sigma_D$ (or $\sigma_B$ for an imaginary part of μ), chosen so that $\mathrm{Im}\, \varepsilon = \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 μ.
+One approach to this problem would be allowing you to specify a constant, frequency-independent, imaginary part of $\varepsilon$, but this has the disadvantage of requiring the simulation to employ complex fields which double the memory and time requirements, and also tends to be numerically unstable. Instead, the approach in Meep is for you to set the conductivity $\sigma_D$ (or $\sigma_B$ for an imaginary part of $\mu$), chosen so that $\mathrm{Im}\, \varepsilon = \varepsilon_\infty \sigma_D / \omega$ is the correct value at your frequency $\omega$ of interest. Note that, in Meep, you specify $f = \omega/2\pi$ instead of $\omega$ for the frequency, however, so you need to include the factor of $2\pi$ when computing the corresponding imaginary part of $\varepsilon$. Conductivities can be implemented with purely real fields, so they are not nearly as expensive as implementing a frequency-independent complex $\varepsilon$ or $\mu$.
 
 For example, suppose you want to simulate a medium with $\varepsilon = 3.4 + 0.101i$ at a frequency 0.42 (in your Meep units), and you only care about the material in a narrow bandwidth around this frequency (i.e. you don't need to simulate the full experimental frequency-dependent permittivity). Then, in Meep, you could use `meep.Medium(epsilon=3.4, D_conductivity=2*math.pi*0.42*0.101/3.4)` in Python or `(make medium (epsilon 3.4) (D-conductivity (* 2 pi 0.42 0.101 (/ 3.4))))` in Scheme; i.e. $\varepsilon_\infty = \mathrm{Re}[\varepsilon] = 3.4$ and $\sigma_D = \omega \, \mathrm{Im}[\varepsilon / \varepsilon_\infty] = (2\pi \, 0.42) \, 0.101 / 3.4$.
 
@@ -99,7 +99,7 @@ For example, suppose you want to simulate a medium with $\varepsilon = 3.4 + 0.1
 Nonlinearity
 ------------
 
-In general, ε can be changed anisotropically by the **E** field itself, with:
+In general, $\varepsilon$ can be changed anisotropically by the $\mathbf{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$$
 
@@ -109,18 +109,18 @@ Meep supports instantaneous, isotropic Pockels and Kerr nonlinearities, correspo
 
 $$\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$\times$3 matrix with the components of **E** along the diagonal.
+Here, "diag($\mathbf{E}$)" indicates the diagonal 3$\times$3 matrix with the components of $\mathbf{E}$ along the diagonal.
 
-Normally, for nonlinear systems you will want to use real fields **E**. This is 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.
+Normally, for nonlinear systems you will want to use real fields $\mathbf{E}$. This is 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 interpretation 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**|<sup>2</sup> is the squared magnitude of the *real* part of **E** for when computing the real part of **D**, and conversely for the imaginary part.
+For complex fields in nonlinear systems, the physical interpretation 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 $\mathbf{E}$ and $\mathbf{D}$) separately: e.g., $|\mathbf{E}^2|$ is the squared magnitude of the *real* part of $\mathbf{E}$ for when computing the real part of $\mathbf{D}$, and conversely for the imaginary part.
 
 For a discussion of how to relate $\chi^{(3)}$ in Meep to experimental Kerr coefficients, see [Units and Nonlinearity](Units_and_Nonlinearity.md).
 
 Magnetic Permeability μ
 -----------------------
 
-All of the above features that are supported for the electric permittivity ε are also supported for the magnetic permeability μ. That is, Meep supports μ with dispersion from magnetic conductivity and Lorentzian resonances, as well as magnetic $\chi^{(2)}$ and $\chi^{(3)}$ nonlinearities. The description of these is exactly the same as above, so we won't repeat it here &mdash; just take the above descriptions and replace ε, **E**, **D**, and σ<sub>D</sub> by μ, **H**, **B**, and σ<sub>B</sub>, respectively.
+All of the above features that are supported for the electric permittivity $\varepsilon$ are also supported for the magnetic permeability $\mu$. That is, Meep supports $\mu$ with dispersion from magnetic conductivity and Lorentzian resonances, as well as magnetic $\chi^{(2)}$ and $\chi^{(3)}$ nonlinearities. The description of these is exactly the same as above, so we won't repeat it here &mdash; just take the above descriptions and replace $\varepsilon$, $\mathbf{E}$, $\mathbf{D}$, and $\sigma_D$ by $\mu$, $\mathbf{H}$, $\mathbf{B}$, and $\sigma_B$, respectively.
 
 Saturable Gain and Absorption
 -----------------------------
@@ -139,7 +139,7 @@ where $\Gamma_{ij}$ is the non-radiative decay or pumping rate from level $i$ to
 
 In Meep, one can specify an arbitrary number of atomic levels with any number of lasing transitions between them, enabling one to realize common two- and four-level saturable media, as well as entire manifolds of levels and transitions found in realistic models of saturable media. When assigning the necessary transition frequencies, $\omega_n$, and widths, $\gamma_n$, of the atomic transition in Meep, these are specified in units of 2π$c$/$a$. However, the pumping and decay rates, $\Gamma_{ij}$, are instead specified in units of $c/a$. Finally, as part of initializing a saturable medium, the total atomic density, $N_0$, must be specified, and Meep will ensure that $\sum_i N_i = N_0$.
 
-In principle, there are a few different ways to alter the effective strength of the gain or absorption experienced by the electric field. For example, to increase the gain of a two-level medium, one could increase the density of the gain medium, $N_0$ (assuming $\Gamma_{12} > \Gamma_{21}$), further increase $\Gamma_{12}$, or increase the relevant components of the coupling strength tensor, $\bar{\sigma}_n$. In practice though, there are two considerations to keep in mind. First, for atomic and molecular gain media, the medium density, coupling strength, and non-radiative decay rates are all fixed by the intrinsic properties of the gain medium -- physically what can be changed is the pumping rate. However, this can occasionally be an inconvenient parameter to tune because it also effects the relaxation rate of the gain medium (see discussion about $\Gamma_\parallel$ below), which can change the stability of the system, see [Applied Physics Letters, Vol. 117, pp. 051102, 2020](https://aip.scitation.org/doi/abs/10.1063/5.0019353?journalCode=apl). Second, Meep will not automatically check to ensure that the non-radiative decay times, $1/\Gamma_{ij}$, are larger than the timestep, $\Delta t$. If the non-radiative decay times are too short, this will lead to unphysical behavior and can result in diverging fields. As such, we would generally recommend tuning the effective gain or loss by changing $N_0$.
+In principle, there are a few different ways to alter the effective strength of the gain or absorption experienced by the electric field. For example, to increase the gain of a two-level medium, one could increase the density of the gain medium, $N_0$ (assuming $\Gamma_{12} > \Gamma_{21}$), further increase $\Gamma_{12}$, or increase the relevant components of the coupling strength tensor, $\bar{\sigma}_n$. In practice though, there are two considerations to keep in mind. First, for atomic and molecular gain media, the medium density, coupling strength, and non-radiative decay rates are all fixed by the intrinsic properties of the gain medium -- physically what can be changed is the pumping rate. However, this can occasionally be an inconvenient parameter to tune because it also effects the relaxation rate of the gain medium (see discussion about $\Gamma_\parallel$ below), which can change the stability of the system, see [Applied Physics Letters, Vol. 117, Num. 051102 (2020)](https://aip.scitation.org/doi/abs/10.1063/5.0019353?journalCode=apl). Second, Meep will not automatically check to ensure that the non-radiative decay times, $1/\Gamma_{ij}$, are larger than the timestep, $\Delta t$. If the non-radiative decay times are too short, this will lead to unphysical behavior and can result in diverging fields. As such, we would generally recommend tuning the effective gain or loss by changing $N_0$.
 
 ### Frequency Domain Susceptibility in Steady-State Regime
 
@@ -220,7 +220,7 @@ The second gyrotropy model is a linearized [Landau-Lifshitz-Gilbert equation](ht
 
 $$\frac{d\mathbf{P}_n}{dt} = \mathbf{b}_n \times \left( - \sigma_n \mathbf{E} + \omega_n \mathbf{P}_n + \alpha_n \frac{d\mathbf{P}_n}{dt} \right) - \gamma_n \mathbf{P}_n$$
 
-Note: although the above equation is written in terms of electric susceptibilities, this model is typically used for magnetic susceptibilities. Meep places no restriction on the field type that either gyrotropy model can be applied to. As usual, electric and magnetic susceptibilities can be swapped by substituting ε with μ, **E** with **H**, etc.
+Note: although the above equation is written in terms of electric susceptibilities, this model is typically used for magnetic susceptibilities. Meep places no restriction on the field type that either gyrotropy model can be applied to. As usual, electric and magnetic susceptibilities can be swapped by substituting ε with μ, $\mathbf{E}$ with $\mathbf{H}$, etc.
 
 The Landau-Lifshitz-Gilbert equation describes the precessional motion of a saturated point magnetic dipole in a magnetic field. In the above equation, the variable $\mathbf{P}_n$ represents the linearized deviation of the polarization from its static equilibrium value (assumed to be much larger and aligned parallel to $\mathbf{b}_n$). Note that this equation of motion is completely different from the [Drude-Lorentz equation](Materials.md#material-dispersion), though the constants σ$_n$, ω$_n$, and γ$_n$ play analogous roles (σ$_n$ couples the polarization to the driving field, ω$_n$ is the angular frequency of precession, and γ$_n$ is a damping factor).
 

doc/docs/Perfectly_Matched_Layer.md

@@ -2,7 +2,7 @@
 # Perfectly Matched Layer
 ---
 
-The [perfectly matched layer](https://en.wikipedia.org/wiki/Perfectly_matched_layer) (**PML**) approach to implementing absorbing boundary conditions in FDTD simulation was originally proposed in [J. Computational Physics, Vol. 114, pp. 185-200 (1994)](http://dx.doi.org/10.1006/jcph.1994.1159). The approach involves surrounding the computational cell with a medium that in theory absorbs without any reflection electromagnetic waves at *all* frequencies and angles of incidence. The original PML formulation involved "splitting" Maxwell's equations into two sets of equations in the absorbing layers, appropriately defined. These split-field equations produce wave attenuation but are unphysical. It was later shown in [IEEE Transactions on Antennas and Propagation, Vol. 43, pp. 1460-3, 1995](https://ieeexplore.ieee.org/abstract/document/477075) that a similar reflectionless absorbing medium can be constructed as a lossy anisotropic dielectric and magnetic material with "matched" impedance and electrical and magnetic conductivities. This is known as the uniaxial PML (UPML).
+The [perfectly matched layer](https://en.wikipedia.org/wiki/Perfectly_matched_layer) (**PML**) approach to implementing absorbing boundary conditions in FDTD simulation was originally proposed in [J. Computational Physics, Vol. 114, pp. 185-200 (1994)](http://dx.doi.org/10.1006/jcph.1994.1159). The approach involves surrounding the computational cell with a medium that in theory absorbs without any reflection electromagnetic waves at *all* frequencies and angles of incidence. The original PML formulation involved "splitting" Maxwell's equations into two sets of equations in the absorbing layers, appropriately defined. These split-field equations produce wave attenuation but are unphysical. It was later shown in [IEEE Transactions on Antennas and Propagation, Vol. 43, pp. 1460-3 (1995)](https://ieeexplore.ieee.org/abstract/document/477075) that a similar reflectionless absorbing medium can be constructed as a lossy anisotropic dielectric and magnetic material with "matched" impedance and electrical and magnetic conductivities. This is known as the uniaxial PML (UPML).
 
 The finite-difference implementation of PML requires the conductivities to be turned on *gradually* over a distance of a few grid points to avoid numerical reflections from the discontinuity. It is also important when using PMLs to make the computational cell sufficiently large so as not to overlap the PML with [evanescent fields](https://en.wikipedia.org/wiki/Evanescent_field) from resonant-cavity or waveguide modes (otherwise, the PML could induce artificial losses in such modes). As a rule of thumb, a PML thickness comparable to half the largest wavelength in the simulation usually works well. The PML thickness should be repeatedly doubled until the simulation results are [sufficiently converged](FAQ.md#checking-convergence).