![]() |
deal.II version GIT relicensing-2652-g41e7496fd4 2025-02-17 19:10:00+00:00
|
This tutorial depends on step-6.
Table of contents | |
---|---|
This program was contributed by Manaswinee Bezbaruah (Texas A&M University), and Matthias Maier (Texas A&M University).
A surface plasmon-polariton (SPP) is a slowly decaying electromagnetic wave that is confined near a metal-air (or similar) interface. SPP structures on novel "2D" materials such as graphene, a monoatomic layer of carbon atoms arranged in a hexagonal lattice, typically have wavelengths much shorter than the wavelength of the free-space radiation. This scale separation makes SPPs on 2D materials a promising ingredient in the design of novel optical devices.
In the following, we discuss a method for observing SPPs numerically by solving a suitable electromagnetic model based on time-harmonic Maxwell's equations which incorporate jump conditions on lower-dimensional material interfaces: The conducting sheet is modeled as an idealized hypersurface with an effective electric conductivity and the weak discontinuity for the tangential surface appears naturally in the variational form.
This tutorial presents a direct solver for the time-harmonic Maxwell equations for scattering configurations with lower-dimensional interfaces. In particular, we discuss using complex values, simple first-order absorbing boundary conditions, and a more sophisticated perfectly matched layer (PML) boundary condition for electromagnetic waves.
We start the discussion with a short derivation of the governing equations and some literature references.
In two ( \(d=2\)) or three ( \(d=3\)) spatial dimensions, the time evolution of an electromagnetic wave \((\mathbf{E},\mathbf{H})\) that consists of an electric field component \(\mathbf{E}(t,\mathbf{x})\;:\;\mathbb{R}\times\mathbb{R}^d\to\mathbb{R}^d\) and a magnetic field component \(\mathbf{H}(t,\mathbf{x})\;:\;\mathbb{R}\times\mathbb{R}^d\to\mathbb{R}^d\) is described by Maxwell's equations [195], [159] :
\begin{align*} \frac{\partial}{\partial t} \mathbf{H} + \nabla \times \mathbf{E} &= -\mathbf{M}_a, \\ \nabla \cdot \mathbf{H} &= \rho_m, \\ \frac{\partial}{\partial t} (\varepsilon\mathbf{E}) - \nabla\times(\mu^{-1}\mathbf{H}) &= - \mathbf{J}_a, \\ \nabla\cdot(\varepsilon\mathbf{E}) &= \rho. \end{align*}
Here, \(\nabla\times\) is the curl operator, \(\nabla\cdot\) is the divergence operator, \(\varepsilon\) is the electric permittivity, \(\mu\) is the magnetic permeability, \(\rho\) is the electric charge density, and \(\rho_m\) is a corresponding (hypothetical) magnetic monopole density. \(\mathbf{J}_a\) and \(\mathbf{M}_a\) are the electric and magnetic flux densities which are related to their respective charge densities by the conservation equations [195]
\[ \frac{\partial}{\partial t} \rho + \nabla\cdot\mathbf{J}_a = 0 \text{ and } \frac{\partial}{\partial t} \rho_m + \nabla\cdot\mathbf{M}_a = 0. \]
We now make the important assumption that the material parameters \(\varepsilon\) and \(\mu\) are time-independent and that the fields \(\mathbf{E}\) and \(\mathbf{H}\), the fluxes \(\mathbf{M}_a\) and \(\mathbf{J}_a\), as well as the densities \(\rho\) and \(\rho_m\) are all time-harmonic, i.e., their time evolution is completely described by
\[ \mathbf{F}(\mathbf{x},t) = \text{Re}\{e^{-i\omega t}\tilde{\mathbf{F}}(\mathbf{x})\}, \]
in which \(\omega\) is the temporal angular frequency and \(\tilde{\mathbf{F}}(\mathbf{x})\) is a corresponding complex-valued vector field (or density). Inserting this ansatz into Maxwell's equations, substituting the charge conservation equations and some minor algebra then yields the so-called time-harmonic Maxwell's equations:
\begin{align*} -i\omega \tilde{\mathbf{H}} + \nabla \times \tilde{\mathbf{E}} &= -\tilde{\mathbf{M}}_a, \\ \nabla \cdot \tilde{\mathbf{H}} &= \frac{1}{i\omega}\nabla \cdot \tilde{\mathbf{M}}_a, \\ i\omega\varepsilon\tilde{\mathbf{E}} + \nabla\times(\mu^{-1}\tilde{\mathbf{H}}) &= \tilde{\mathbf{J}}_a, \\ \nabla\cdot(\varepsilon\tilde{\mathbf{E}}) &= \frac{1}{i\omega}\nabla\cdot\tilde{\mathbf{J}}_a. \end{align*}
For the sake of better readability we will now drop the tilde and simply write \(\mathbf{E}(\mathbf{x})\), \(\mathbf{H}(\mathbf{x})\), etc., when referring to the time-harmonic fields.
Graphene is a two-dimensional carbon allotrope with a single atom layer that is arranged in a honeycomb lattice [170]. Due to its atomic thickness it is an example of a so-called 2D material: Compared to the other spatial dimensions (where graphene samples can reach up to several centimeters) the atomistic thickness of graphene typically ranges around 2.5 ångstrom ( \(2.5\times10^{-10}\text{m}\)). We will thus model graphene as a lower-dimensional interface \(\Sigma\) embedded into the computational domain \(\Omega\subset\mathbb{R}^d\). More precisely, \(\Sigma\) is a two-dimensional sheet in three spatial dimensions, or a one-dimensional line in two spatial dimensions. The special electronic structure of graphene gives rise to a current density on the lower-dimensional interface that is modeled with an effective surface conductivity \(\sigma^\Sigma\) obeying Ohm's Law:
\[ \mathbf{J}^\Sigma=\sigma^\Sigma\,\mathbf{E}_T \]
in which \(\mathbf{J}^\Sigma\) is the surface current density, \(\mathbf{E}_T\) denotes the tangential part of the electric field \(\mathbf{E}\), and \(\sigma^\Sigma\) is an appropriately chosen surface conductivity that will be discussed in more detail below. The surface current density gives rise to a jump condition on \(\Sigma\) in the tangential component of the magnetic field. This is best seen by visualizing Ampère's law:
and then taking the limit of the upper and lower part of the line integral approaching the sheet. In contrast, the tangential part of the electric field is continuous. By fixing a unit normal \(\mathbf{\nu}\) on the hypersurface \(\Sigma\) both jump conditions are
\begin{align*} \mathbf{\nu} \times \left[(\mu^{-1}\mathbf{H})^+ - (\mu^{-1}\mathbf{H})^-\right]|_{\Sigma} &= \sigma^{\Sigma}\left[(\mathbf{\nu}\times \mathbf{E}\times \mathbf{\nu})\right]|_{\Sigma}, \\ \mathbf{\nu} \times \left[\mathbf{E}^+ - \mathbf{E}^-\right]|_{\Sigma} &= 0. \end{align*}
The notation \(\mathbf{F}^\pm\) indicates the limit values of the field when approaching the interface from above or below the interface: \(\mathbf{F}^\pm(\mathbf{x})=\lim_{\delta\to0,\delta>0}\mathbf{F}(\mathbf{x}\pm\delta\mathbf{\nu})\).
We will be using a rescaled version of the Maxwell's equations described above. The rescaling has the following key differences:
\[ \mathbf{J}_a = J_0 \mathbf{e}_i\delta(x-a) \]
Accordingly, our electric permittivity and magnetic permeability are rescaled by \(\varepsilon_0\) and \(\mu_0\) as
\[ \mu_r = \frac{1}{\mu_0}\mu \text{ and } \varepsilon_r = \frac{1}{\varepsilon_0}\varepsilon. \]
We use the free-space wave number \(k_0 = \omega\sqrt{\varepsilon_0\mu_0}\) and the dipole strength, \(J_0\) to arrive at the following rescaling of the vector fields and coordinates:
\begin{align*} \hat{x} = k_0x, &\qquad \hat{\nabla} = \frac{1}{k_0}\nabla,\\ \hat{\mathbf{H}} = \frac{k_0}{J_0}\mu^{-1}\mathbf{H},&\qquad \hat{\mathbf{E}} = \frac{k_0^2}{\omega\mu_0 J_0}\mathbf{E},\\ \hat{\mathbf{J}}_a = \frac{1}{J_0}\mathbf{J}_a,&\qquad \hat{\mathbf{M}}_a = \frac{k_0}{\omega\mu_0 J_0}\mathbf{M}_a. \end{align*}
Finally, the interface conductivity is rescaled as
\[ \sigma^{\Sigma}_r = \sqrt{\frac{\mu_0}{\varepsilon_0}}\sigma^{\Sigma}. \]
Accordingly, our rescaled equations are
\begin{align*} -i\mu_r \hat{\mathbf{H}} + \hat{\nabla} \times \hat{\mathbf{E}} &= -\hat{\mathbf{M}}_a, \\ \hat{\nabla} \cdot (\mu_r\hat{\mathbf{H}}) &= \frac{1}{i\omega}\hat{\nabla} \cdot \hat{\mathbf{M}}_a, \\ i\varepsilon_r\hat{\mathbf{E}} + \nabla\times(\mu^{-1}\mathbf{H}) &= \mathbf{J}_a, \\ \nabla\cdot(\varepsilon\mathbf{E}) &= \frac{1}{i\omega}\hat{\nabla} \cdot\hat{\mathbf{J}}_a. \end{align*}
We will omit the hat in further discussion for ease of notation.
Let \(\Omega \subset \mathbb{R}^n\), \((n = 2,3)\) be a simply connected and bounded domain with Lipschitz-continuous and piecewise smooth boundary, \(\partial\Omega\). Let \(\Sigma\) be an oriented, Lipschitz-continuous, piecewise smooth hypersurface. Fix a normal field \(\nu\) on \(\Sigma\) and let \(n\) denote the outer normal vector on \(\partial\Omega\).
In order to arrive at the variational form, we will substitute for \(\mathbf{H}\) in the first equation and obtain
\[ \nabla \times (\mu_r^{-1}\nabla\times\mathbf{E}) - \varepsilon_r \mathbf{E} = i\mathbf{J}_a - \nabla\times (\mu_r^{-1}\mathbf{M}_a). \]
Now, consider a smooth test function \(\varphi\) with complex conjugate \(\bar{\varphi}\). Multiply both sides of the above equation by \(\bar{\varphi}\) and integrate by parts in \(\Omega\backslash\Sigma\).
\[ \int_\Omega (\mu_r^{-1}\nabla\times\mathbf{E})\cdot (\nabla\times\bar{\varphi})\;\text{d}x - \int_\Omega \varepsilon_r\mathbf{E} \cdot \bar{\varphi}\;\text{d}x - \int_\Sigma [\nu \times (\mu_r^{-1}\nabla\times\mathbf{E} + \mu^{-1}\mathbf{M}_a)]_{\Sigma}\cdot \bar{\varphi}_T\;\text{d}o_x\\ \qquad - \int_{\partial\Omega} (\nu \times (\mu_r^{-1}\nabla\times\mathbf{E} + \mu^{-1}\mathbf{M}_a)) \cdot \bar{\varphi}_T\;\text{d}o_x = i\int_\Omega \mathbf{J}_a \cdot \bar{\varphi}\;\text{d}x - \int_\Omega \mu_r^{-1}\mathbf{M}_a \cdot (\nabla \times \bar{\varphi})\;\text{d}x. \]
We use the subscript \(T\) to denote the tangential part of the given vector and \([\cdot]_{\Sigma}\) to denote a jump over \(\Sigma\), i.e.,
\[ \mathbf{F}_T = (\mathbf{\nu}\times \mathbf{F})\times\mathbf{\nu} \text{ and } [\mathbf{F}]_{\Sigma}(\mathbf{x}) = \lim\limits_{s\searrow 0}(\mathbf{F}(\mathbf{x}+s\mathbf{\nu})-\mathbf{F}(\mathbf{x}-s\mathbf{\nu})) \]
for \(\mathbf{x}\in \Sigma\).
For the computational domain \(\Omega\), we introduce the absorbing boundary condition at \(\partial\Omega\), which is obtained by using a first-order approximation of the Silver-Müller radiation condition, truncated at \(\partial\Omega\) [159].
\[ \nu\times\mathbf{H}+\sqrt{\mu_r^{-1}\varepsilon_r}\mathbf{E}=0\qquad x\in\partial\Omega \]
We assume that \(\mu_r^{-1}\) and \(\varepsilon\) have well-defined square roots. In our numerical computation, we combine the above absorbing boundary condition with a PML.
The jump condition can be expressed as a weak discontinuity as follows:
\[ [\nu \times (\mu_r^{-1}\nabla\times\mathbf{E} + \mu^{-1}\mathbf{M}_a)]_{\Sigma} = i\sigma_r^{\Sigma}\mathbf{E}_T,\qquad \text{on }\Sigma\\ \nu \times (\mu_r^{-1}\nabla\times\mathbf{E} + \mu^{-1}\mathbf{M}_a) = -i\sqrt{\mu_r^{-1}\varepsilon_r}\mathbf{E}_T,\qquad \text{on }\partial\Omega. \]
Combining, our weak form is as follows:
\[ \int_\Omega (\mu_r^{-1}\nabla\times\mathbf{E})\cdot (\nabla\times\bar{\varphi})\;\text{d}x - \int_\Omega \varepsilon_r\mathbf{E} \cdot \bar{\varphi}\;\text{d}x - i\int_\Sigma (\sigma_r^{\Sigma}\mathbf{E}_T) \cdot \bar{\varphi}_T\;\text{d}o_x\\ \qquad - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}\mathbf{E}_T) \cdot (\nabla\times\bar{\varphi}_T)\;\text{d}o_x.= i\int_\Omega \mathbf{J}_a \cdot \bar{\varphi}\;\text{d}x - \int_\Omega \mu_r^{-1}\mathbf{M}_a \cdot (\nabla \times \bar{\varphi})\;\text{d}x. \]
Assume that \(\sigma_r^{\Sigma} \in L^{\infty}(\Sigma)^{2\times 2}\) is matrix-valued and symmetric, and has a semidefinite real and complex part. Let \(\varepsilon_r\) be a smooth scalar function with \(-\text{Im}(\varepsilon_r) = 0\), or \(\text{Im}(\varepsilon_r)\ge c > 0\) in \(\Omega\). \(\mu_r^{-1}\) is a smooth scalar such that \(\sqrt{\mu_r^{-1}\varepsilon_r}\) is real valued and strictly positive in \(\partial\Omega\).
\(\mathbf{H}(\text{curl};\Omega)\) is space of vector-valued, measurable and square integrable functions whose weak curl admits a representation by a square integrable function. Define a Hilbert space
\[ X(\Omega) = \{\varphi \in \mathbf{H}(\text{curl};\Omega)\;\;:\;\; \varphi_T|_{\Sigma} \in L^2(\Sigma)^2,\;\varphi_T|_{\partial\Omega} \in L^2(\partial\Omega)^2\} \]
equipped with the norm
\[ \|\varphi\|^2_X = \|\varphi\|^2_{L^2(\Omega)} + \|\nabla\times\varphi\|^2_{L^2(\Omega)} + \|\varphi_T\|^2_{L^2(\Sigma)} + \|\varphi_T\|^2_{L^2(\partial\Omega)}. \]
Define
\[ A(\mathbf{E},\varphi) \dealcoloneq \int_\Omega (\mu_r^{-1}\nabla\times\mathbf{E})\cdot (\nabla\times\bar{\varphi})\;\text{d}x - \int_\Omega \varepsilon_r\mathbf{E} \cdot \bar{\varphi}\;\text{d}x - i\int_\Sigma (\sigma_r^{\Sigma}\mathbf{E}_T) \cdot \bar{\varphi}_T\;\text{d}o_x - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}\mathbf{E}_T) \cdot (\nabla\times\bar{\varphi}_T)\;\text{d}o_x.\\ F(\varphi) \dealcoloneq i\int_\Omega \mathbf{J}_a \cdot \bar{\varphi}\;\text{d}x - \int_\Omega \mu_r^{-1}\mathbf{M}_a \cdot (\nabla \times \bar{\varphi})\;\text{d}x. \]
Then, our rescaled weak formulation is:
\[ \text{Find a unique } \mathbf{E} \in X(\Omega) \text{ such that, for all } \varphi \in X(\Omega), \]
\[ A(\mathbf{E},\varphi) = F(\varphi). \]
Moreover, the above equations are supplemented by the Silver-Müller radiation condition, if the ambient (unbounded) medium is isotropic. This amounts to the requirement that \(\mathbf{E}\) and \(\mathbf{H}\) both approach a spherical wave uniformly in the radial direction for points at infinity and away from the conducting sheet, i.e.,
\[ \lim\limits_{|x|\to\infty} \{\mathbf{H}\times x - c^{-1}|x|\mathbf{E}\} = 0 \text{ and } \lim\limits_{|x|\to\infty} \{\mathbf{E}\times x - c^{-1}|x|\mathbf{H}\} = 0 \text{ for } x \not\in \Sigma. \]
In our case, we eliminate reflection from infinity by implementing a PML, which is described at length below, and avoid the explicit use of the last condition.
The variational form is discretized on a non-uniform quadrilateral mesh with higher-order, curl-conforming Nédélec elements implemented by the FE_NedelecSZ class. This way the interface with a weak discontinuity is optimal, and we get optimal convergence rates.
Consider the finite element subspace \(X_h(\Omega) \subset X(\Omega)\). Define the matrices
\[ A_{ij} = \int_\Omega (\mu_r^{-1}\nabla \times \varphi_j) \cdot (\nabla\times\bar{\varphi}_i)\;\text{d}x - \int_\Omega \varepsilon_r\varphi_j \cdot \bar{\varphi}_i\;\text{d}x - i\int_\Sigma (\sigma_r^{\Sigma}\varphi_{j_T}) \cdot \bar{\varphi}_{i_T}\;\text{d}o_x - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}\varphi_{j_T}) \cdot (\nabla\times \bar{\varphi}_{i_T})\;\text{d}o_x, \]
\[ F_i = i\int_\Omega \mathbf{J}_a \cdot \bar{\varphi_i}\;\text{d}x - \int_\Omega \mu_r^{-1}\mathbf{M}_a \cdot (\nabla \times \bar{\varphi_i}) \;\text{d}x. \]
Then under the assumption of a sufficiently refined initial mesh the discretized variational problem is:
\[ \text{Find a unique } \mathbf{E}_h = \sum_j U_j\mathbf{\varphi}_j \in X_h(\Omega) \text{ such that} \]
\[ \sum_jA_{ij}U_j = F_i\qquad\text{for all }i. \]
The SPP amplitude is negatively affected by the absorbing boundary condition and this causes the solution image to be distorted. In order to reduce the resonance and distortion in our solutions, we are implementing a Perfectly Matched Layer (PML) in the scattering configuration.
The concept of a PML was pioneered by Bérenger [28] and it is is an indispensable tool for truncating unbounded domains for wave equations and often used in the numerical approximation of scattering problems. It is essentially a thin layer with modified material parameters placed near the boundary such that all outgoing electromagnetic waves decay exponentially with no “artificial” reflection due to truncation of the domain.
Our PML is a concentric circle with modified material coefficients ( \(\varepsilon_r, \mu_r, \sigma\)). It is located in a small region near the boundary \(\partial\Omega\) and the transformation of the material coordinates is chosen to be a function of the radial distance \(\rho\) from the origin \(e_r\). The normal field \(\nu\) of \(\Sigma\) is orthogonal to the radial direction \(e_r\), which makes \(\mathbf{J}_a \equiv 0\) and \(\mathbf{M}_a \equiv 0\) within the PML.
Introduce a change of coordinates
\[ x \to \bar{x} = \begin{cases} x + ie_r\int\limits_\rho^r s(\tau)\text{d}\tau,\;\;\;\;\;\;\; r\ge\rho\\ x\;\;\;\;\;\;\;\;\;\text{otherwise} \end{cases} \]
in which \(r = e_r \cdot x\) and \(s(\tau)\) is an appropriately chosen, nonnegative scaling function.
We introduce the following \(2\times2\) matrices
\begin{align*} A &= T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{d}^2}, \frac{1}{d\bar{d}}\right)T_{e_xe_r} \\ B &= T_{e_xe_r}^{-1} \text{diag}\left(d,\bar{d}\right)T_{e_xe_r} \\ C &= T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{d}},\frac{1}{d}\right) T_{e_xe_r} \end{align*}
in which
\begin{align*} d &= 1 + is(r) \\ \bar{d} &= 1 + i/r \int\limits_{\rho}^{r}s(\tau)\text{d}\tau \end{align*}
and \(T_{e_xe_r}\) is the rotation matrix which rotates \(e_r\) onto \(e_x\). Thus, after applying the rescaling, we get the following modified parameters
\begin{align*} \bar{\mu}_r^{-1} &= \frac{\mu_r^{-1}}{d}, \\ \bar{\varepsilon}_r &= A^{-1} \varepsilon_r B^{-1}, \text{ and } \\ \bar{\sigma}^{\Sigma}_r &= C^{-1} \sigma^{\Sigma}_r B^{-1}. \end{align*}
These PML transformations are implemented in our PerfectlyMatchedLayer
class. After the PML is implemented, the electromagnetic wave essentially decays exponentially within the PML region near the boundary, therefore reducing reflection from the boundary of our finite domain. The decay function also depends on the strength of our PML, which can be adjusted to see more or less visible decaying in the PML region.
The set of include files is quite standard. The most notable include is the fe/fe_nedelec_sz.h file which allows us to use the FE_NedelecSZ elements. This is an implementation of the \(H^{curl}\) conforming Nédélec Elements that resolves the sign conflict issues that arise from parametrization (for details we refer to the documentation of the FE_NedelecSZ element).
We begin our actual implementation by declaring all classes with their data structures and methods upfront.
The Parameters class inherits ParameterAcceptor, and instantiates all the coefficients in our variational equations. These coefficients are passed through ParameterAcceptor and are editable through a .prm file. More explanation on the use and inheritance from the ParameterAcceptor can be found in step-60.
epsilon is the Electric Permittivity coefficient and it is a rank 2 tensor. Depending on the material, we assign the i^th diagonal element of the tensor to the material epsilon value (one of the private epsilon_1_ or epsilon_2_ variables).
mu_inv is the inverse of the Magnetic Permiability coefficient and it is a complex number.
sigma is the Surface Conductivity coefficient between material left and material right and it is a rank 2 tensor. It is only changed if we are at the interface between two materials. If we are at an interface, we assign the i^th diagonal element of the tensor to the private sigma_ value.
J_a is the strength and orientation of the dipole. As mentioned in the rescaling,
\[ \mathbf{J}_a = J_0 \mathbf{e}_i\delta(x-a) \]
It is a rank 1 tensor that depends on the private dipole_position, dipole_radius, dipole_strength, dipole_orientation variables.
The PerfectlyMatchedLayer class inherits ParameterAcceptor as well. It implements the transformation matrices used to modify the permittivity and permeability tensors supplied from the Parameters class. The actual transformation of the material tensors will be done in the assembly loop. The radii and the strength of the PML is specified, and the coefficients will be modified using transformation matrices within the PML region. The radii and strength of the PML are editable through a .prm file. The rotation function \(T_{exer}\) is the same as introduced in the perfectly matched layer section of the introduction. Similarly, the matrices A, B and C are defined as follows
\[ A = T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{d}^2},\frac{1}{d\bar{d}}\right)T_{e_xe_r},\qquad B = T_{e_xe_r}^{-1} \text{diag}\left(d,\bar{d}\right)T_{e_xe_r},\qquad C = T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{d}},\frac{1}{d}\right) T_{e_xe_r}.\qquad \]
At this point we are ready to declare all the major building blocks of the finite element program which consists of the usual setup and assembly routines. Most of the structure has already been introduced in previous tutorial programs. The Maxwell class also holds private instances of the Parameters and PerfectlyMatchedLayers classes introduced above. The default values of these parameters are set to show us a standing wave with absorbing boundary conditions and a PML.
The Constructor simply consists of default initialization a number of discretization parameters (such as the domain size, mesh refinement, and the order of finite elements and quadrature) and declaring a corresponding entry via ParameterAcceptor::add_parameter(). All of these can be modified by editing the .prm file. Absorbing boundary conditions can be controlled with the absorbing_boundary boolean. If absorbing boundary conditions are disabled we simply enforce homogeneous Dirichlet conditions on the tangential component of the electric field. In the context of time-harmonic Maxwell's equations these are also known as perfectly conducting boundary conditions.
The Maxwell::make_grid() routine creates the mesh for the computational domain which in our case is a scaled square domain. Additionally, a material interface is introduced by setting the material id of the upper half ( \(y>0\)) to 1 and of the lower half ( \(y<0\)) of the computational domain to 2. We are using a block decomposition into real and imaginary matrices for the solution matrices. More details on this are available under the Results section.
The Maxwell::setup_system() routine follows the usual routine of enumerating all the degrees of freedom and setting up the matrix and vector objects to hold the system data. Enumerating is done by using DoFHandler::distribute_dofs().
This is a helper function that takes the tangential component of a tensor.
Assemble the stiffness matrix and the right-hand side:
\begin{align*} A_{ij} = \int_\Omega (\mu_r^{-1}\nabla \times \varphi_j) \cdot (\nabla\times\bar{\varphi}_i)\text{d}x - \int_\Omega \varepsilon_r\varphi_j \cdot \bar{\varphi}_i\text{d}x - i\int_\Sigma (\sigma_r^{\Sigma}(\varphi_j)_T) \cdot (\bar{\varphi}_i)_T\text{do}x - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\varphi_j)_T) \cdot (\nabla\times(\bar{\varphi}_i)_T)\text{d}x, \end{align*}
\begin{align} F_i = i\int_\Omega J_a \cdot \bar{\varphi_i}\text{d}x - \int_\Omega \mu_r^{-1} \cdot (\nabla \times \bar{\varphi_i}) \text{d}x. \end{align}
In addition, we will be modifying the coefficients if the position of the cell is within the PML region.
Next, let us assemble on the interior of the domain on the left hand side. So we are computing
\begin{align*} \int_\Omega (\mu_r^{-1}\nabla \times \varphi_i) \cdot (\nabla\times\bar{\varphi}_j)\text{d}x - \int_\Omega \varepsilon_r\varphi_i \cdot \bar{\varphi}_j\text{d}x \end{align*}
and
\begin{align} i\int_\Omega J_a \cdot \bar{\varphi_i}\text{d}x - \int_\Omega \mu_r^{-1} \cdot (\nabla \times \bar{\varphi_i}) \text{d}x. \end{align}
In doing so, we need test functions \(\varphi_i\) and \(\varphi_j\), and the curl of these test variables. We must be careful with the signs of the imaginary parts of these complex test variables. Moreover, we have a conditional that changes the parameters if the cell is in the PML region.
Now we assemble the face and the boundary. The following loops will assemble
\begin{align*} - i\int_\Sigma (\sigma_r^{\Sigma}(\varphi_i)_T) \cdot (\bar{\varphi}_j)_T\text{do}x \end{align*}
and
\begin{align} - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\varphi_i)_T) \cdot (\nabla\times(\bar{\varphi}_j)_T)\text{d}x, \end{align}
respectively. The test variables and the PML are implemented similarly as the domain.
We are on an interior face:
We use a direct solver from the SparseDirectUMFPACK to solve the system
The output is written into a vtk file with 4 components
The following main function calls the class step-81(), initializes the ParameterAcceptor, and calls the run() function.
The solution is written to a .vtk file with four components. These are the real and imaginary parts of the \(E_x\) and \(E_y\) solution waves. With the current setup, the output should read
The following images are the outputs for the imaginary \(E_x\) without the interface and with the dipole centered at \((0,0)\). In order to remove the interface, the surface conductivity is set to 0. First, we turn off the absorbing boundary conditions and the PML. Second, we want to see the effect of the PML when absorbing boundary conditions apply. So we set absorbing boundary conditions to true and leave the PML strength to 0. Lastly, we increase the strength of the PML to 4. Change the following in the .prm file :
Following are the output images:
Solution with no interface, Dirichlet boundary conditions and PML strength 0. |
Solution with no interface, absorbing boundary conditions and PML strength 0. |
Solution with no interface, absorbing boundary conditions and PML strength 4. |
We observe that with absorbing boundary conditions and in absence of the PML, there is a lot of distortion and resonance (the real parts will not be generated without a PML). This is, as we stipulated, due to reflection from infinity. As we see, a much more coherent image is generated with an appropriate PML.
Now, let's generate a standing wave by adding an interface at the center. In order to observe this effect, we offset the center of the dipole to \((0, 0.8)\) and set the surface conductivity back to \((0.001, 0.2)\):
Once again, we will visualize the output with absorbing boundary conditions and PML strength 0 and with absorbing boundary conditions and PML strength
Solution with an interface, Dirichlet boundary conditions and PML strength 0. |
Solution with an interface, absorbing boundary conditions and PML strength 0. |
Solution with an interface, absorbing boundary conditions and PML strength 4. |
Solution with an interface, Dirichlet boundary conditions and PML strength 0. |
Solution with an interface, absorbing boundary conditions and PML strength 0. |
Solution with an interface, absorbing boundary conditions and PML strength 4. |
The SPP is confined near the interface that we created, however without absorbing boundary conditions, we don't observe a dissipation effect. On adding the absorbing boundary conditions, we observe distortion and resonance and we still don't notice any dissipation. As expected, the PML removes the distortion and resonance. The standing wave is also dissipating and getting absorbed within the PML, and as we increase the PML strength, the standing wave will dissipate more within the PML ring.
Here are some animations to demonstrate the effect of the PML
Solution with an interface, Dirichlet boundary conditions and PML strength 0. |
Solution with an interface, absorbing boundary conditions and PML strength 0. |
Solution with an interface, absorbing boundary conditions and PML strength 4. |
Solution with an interface, Dirichlet boundary conditions and PML strength 0. |
Solution with an interface, absorbing boundary conditions and PML strength 0. |
Solution with an interface, absorbing boundary conditions and PML strength 4. |
As is evident from the results, we are splitting our solution matrices into the real and the imaginary components. We started off using the \(H^{curl}\) conforming Nédélec Elements, and we made two copies of the Finite Elements in order to represent the real and the imaginary components of our input (FE_NedelecSZ was used instead of FE_Nedelec to avoid the sign conflicts issues present in traditional Nédélec elements). In the assembly, we create two vectors of dimension \(dim\) that assist us in extracting the real and the imaginary components of our finite elements.
As we see in our assembly, our finite element is rotated and scaled as follows:
This \(\phi_i\) variable doesn't need to be scaled in this way, we may choose any arbitrary scaling constants \(a\) and \(b\). If we choose this scaling, the \(\phi_j\) must also be modified with the same scaling, as follows:
Moreover, the cell_rhs need not be the real part of the rhs_value. Say if we modify to take the imaginary part of the computed rhs_value, we must also modify the cell_matrix accordingly to take the imaginary part of temp. However, making these changes to both sides of the equation will not affect our solution, and we will still be able to generate the surface plasmon polariton.
We will create a video demonstrating the wave in motion, which is essentially an implementation of \(e^{-i\omega t}(Re(E) + i*Im(E))\) as we increment time. This is done by slightly changing the output function to generate a series of .vtk files, which will represent out solution wave as we increment time. Introduce an input variable \(t\) in the output_results() class as output_results(unsigned int t). Then change the class itself to the following:
Finally, in the run() function, replace output_results() with
This would generate 100 solution .vtk files, which can be opened in a group on Paraview and then can be saved as an animation. We used FFMPEG to generate gifs.
The example step could be extended in a number of different directions.