Total Field/Scattered Field
Section author: Sergei Bastrakov
This section describes the Total Field/Scattered Field (TF/SF) technique for adding incident field to a simulation. TF/SF is used for simulating a field produced by external sources and coming towards the simulation volume from any boundaries and directions. PIConGPU implements this technique as incident field. There are also other ways to generate a laser.
First, we present the general idea of TF/SF and describe its formulation for the standard Yee field solver. Then we show how it is generalized for any FDTD-type solver, this general version is implemented in PIConGPU. Finally, some aspects of using the PIConGPU incident field are covered.
Introduction
When applying TF/SF, we operate on a domain where field values are updated by an explicit FDTD-type Maxwell’s solver. In case a simulation uses a field absorber, the absorber operates externally to this area and independently of TF/SF (no special coupling needed).
This section uses the same Yee grid and notation as the FDTD description. The description and naming is generally based on [Potter2017] [Taflove2005], however there are some differences in notation (also those two sources use different Yee grid layouts).
The field produced in the domain by external sources (such as incoming laser) is called incident field and denoted \(\vec E^{inc}(x, y, z, t)\) and \(\vec B^{inc}(x, y, z, t)\). Total field (TF) is the combined field from the internally occuring field dynamic and the incident field, we denote it \(\vec E^{tot}(x, y, z, t)\) and \(\vec B^{tot}(x, y, z, t)\). Thus, total field is a “normal” result of a simulation with the laser. Finally, the scattered field (SF) is the difference of these two fields:
Scattered field does not have the laser effect, and transparently propagates the waves outgoing from the TF area towards domain boundaries. However, note that Maxwell’s equations hold in both the TF and SF regions. Special handling is needed only near the boundary betweeen the regions.
The field values are represented as a Yee grid of \(\vec E\) and \(\vec B\) values in the domain, TF/SF does not require additionally stored data. However, it changes the interpretation of grid values: the domain is virtually subdivided into the internal area where the stored grid values are TF and external area with SF values.
The subvidision is done with an axis-aligned Huygens surface \(S\). The Huygens surface is chosen so that no Yee grid nodes lay on it. So it could be located at an arbitrary position that does not collide with cell- and half-cell boundaries.
In PIConGPU, the position of the Huygens surface is defined relative to the internal boundary of the field absorber. The surface is located inwards relative to each boundary by a user-defined position in full cells and an additional 0.75 cells. The equations presented in this section hold for this 0.75 shift and the Yee grid layout used in PIConGPU. In principle, a similar derivation can be done in any case, but the resulting expression would not generally match index-by-index.
We use TF/SF formulation to generate fields at the Huygens surface that propagate inwards the simulation area and act as application of \(\vec E^{inc}(x, y, z, t)\), \(\vec B^{inc}(x, y, z, t)\). The fields would mostly propagate in this direction (and not equally to both sides), and as a theoretical limit only inwards. However, in practice due to numerical dispersion a small portion, such as \(10^{-3}\) of the amplitude, would also propagate externally.
The TF/SF technique can be interpreted in two ways. First, as manipulations on discrete FDTD equations given the virtual separation of the domain described above. This “math” view is used in the following description for the Yee and general-case FDTD solvers.
Alternatively, TF/SF can be treated as applying the equivalence theorem on the Huygens surface. Then it is equivalent to impressing the electric and magnetic current sheets, \(\vec J\) and \(\vec K\) respectively, on \(S\):
where \(\vec n\) is a unit normal vector on \(S\) pointing towards the internal volume. Both points of view in principle result in the same scheme.
TF/SF for Standard Yee Solver
Single Boundary
First we describe application of a field source only at the \(x_{min}\) boundary. For this case, suppose the source is given along a plane \(x = x_{min} + \mathrm{position} \Delta x + 0.75 \Delta x\) and acts along the whole domain in \(y\) and \(z\). The source affects transversal field components \(E_y\), \(E_z\), \(B_y\), \(B_z\). Components \(E_x\), \(B_x\) are not affected by TF/SF and we do not consider them in the following description. We also omit the \(\vec J\) term which is applied as usual in the field solver and is also not affected by TF/SF.
The Yee solver update equations for the affected field components and terms are as follows:
When using TF/SF technique, first a usual Yee field solver update is applied to the whole grid, regardless of TF and SF regions. Then a separate stage that we call incident field solver is run to modify the calculated values where necessary. The combined effect of the Yee- and incident field solvers is that Maxwell’s equations hold on the whole grid and the correct incident field is generated. We now proceed to describe how are these values identified and what is the modification necessary.
As mentioned above, values like \(E_y\rvert_{i, j+1/2, k}^{n+1}\) are stored for the whole Yee grid. Whether they represent the total or the scattered field, depends on the position of the node relative to the Huygens surface. To avoid confusion, we use the \(E_y\rvert_{i, j+1/2, k}^{n+1}\) notation for stored grid values, and \(E_y^{other}\left( i \Delta x, (j+1/2) \Delta y, k \Delta z, (n+1) \Delta t \right)\) to denote fields at the same time and space position, but not stored long-term.
Since the Maxwell’s equations hold in both the TF and SF regions, all Yee solver updates involving only grid values from the same region produced correct values that do not need any further modification. A correction is only needed for grid values that were calculated using a mix of TF and SF values. Since the standard Yee solver uses a 2-point central derivative operator, those are a single layer of \(\vec E\) and \(\vec B\) values located near \(S\).
Taking into account the 0.75 shift inwards used by PIConGPU, denote the \(x\) position of \(S\) as \(x_S = (i_S + 0.75) \Delta x\). Then the grid values to be modified by the incident field solver are \(E_y\rvert_{i_S+1, j+1/2, k}^{n+1}\), \(E_z\rvert_{i_S+1, j, k+1/2}^{n+1}\), \(B_y\rvert_{i_S+1/2, j, k+1/2}^{n+3/2}\), and \(B_z\rvert_{i_S+1/2, j+1/2, k}^{n+3/2}\). (All grid values to the right side of those were calculated using only TF values and all grid values on the left side were calculated using only SF values.)
Consider the update of \(E_y\rvert_{i_S+1, j+1/2, k}^{n+1}\) performed by a standard Yee solver for each \(j, k\). All terms but \(B_z\rvert_{i_S+1/2, j+1/2, k}^{n+1/2}\) in this update are in the TF region. Thus, this value has to be modified by the incident field solver in order to preseve the Maxwell’s equations.
To derive the modification necessary, consider a hypothetical Maxwell’s-preserving update at this point where all participating values were TF:
Since \(B_z\rvert_{i_S+1/2, j+1/2, k}^{n+1/2}\) is an SF and by definition of TF and SF,
Substituting it into the update equation and regrouping the terms yields:
Thus, in the incident field stage we have to apply the following update to the grid value calculated by a normal Yee solver :
Grid value \(E_z\rvert_{i_S+1, j, k+1/2}^{n+1}\) is also located in the TF region and with a similar derivation the update for it is
Values \(B_y\rvert_{i_S+1/2, j, k+1/2}^{n+3/2}\), and \(B_z\rvert_{i_S+1/2, j+1/2, k}^{n+3/2}\) are in the SF region. For them the Yee solver update includes one term from the TF region, \(E_z\rvert_{i_S, j, k+1/2}^{n+1}\) and \(E_y\rvert_{i_S, j+1/2, k}^{n+1}\) respectively. Making a similar replacement of an SF value as a difference between a TF value and the incident field value and regrouping, the following update must be applied:
The derivation for the \(x_{max}\) boundary can be done in a similar fashion. Denote the position of \(S\) as \(x_S = (i_{S, max} + 0.25) \Delta x\). Note that our 0.75 cells inwards shift of \(S\) is symmetrical in terms of distance. It implies that the Yee grid incides along \(x\) are not fully symmetric between the two sides of each bondary. The update scheme for \(x_{max}\) is:
Multiple Boundaries
In the general case, \(S\) is comprised of several axis-aligned boundary hyperplanes, 6 planes in 3D, and 4 lines in 2D.
The scheme described above is symmetric for all axes. In case incident field is coming from multiple boundaries, the updates are in principle the same. They can generally be treated as sequential application of the single-boundary case.
Applying TF/SF for each boundary affects the derivatives in the normal direction relative to the boundary. For the standard Yee solver, a single layer of \(\vec E\) and \(\vec B\) values along the boundary is affected. Along other directions, we update all grid values that are internal relative to the Huygens surface. In case a “corner” grid node is near several boundaries, it is updated in all the respective applications of TF/SF.
General Case FDTD
The same principle as for the Yee solver can be applied to any FDTD-type field solver. Same as above, we consider the case of \(x_{min}\) boundary and \(E_y\) field component. The other boundaries and components are treated symmetrically.
We now apply a general-case spatial-only finite-difference operator to calculate derivatives along \(x\). Such operators on the Yee grid naturally have an antisymmetry of coefficients in \(x\) relative to the evaluation point. The antisymmetry is not critical for the following description, but is present in the FDTD solvers implemented and allow simplifying the formulations, and so we assume it. For \(dB_z/dx\) such an operator has the following general form:
Note that there is also typically a symmetry of coefficients along \(y\) and \(z\): \(\alpha_{ii, jj, kk} = \alpha_{ii, -jj, kk}\), \(\alpha_{ii, jj, kk} = \alpha_{ii, jj, -kk}\) but it is not significant for TF/SF. The derivative operator used by the standard Yee solver has \(m_x = m_y = m_z = 0, \alpha_{0, 0, 0} = 1\).
Same as before, denote the \(x\) position of \(S\) as \(x_S = (i_S + 0.75) \Delta x\). In order to stay within the grid, we require that \(i_S \geq m_x\). The incident field solver has to update the grid values of \(E_y\) for which calculating \(dB_z/dx\) involves a mix of TF and SF values. These values can be present in both the TF and SF regions around \(S\):
Take a node in the TF region \(E_y\rvert_{i_0, j+1/2, k}^{n+1} \in E_{TF}\) (\(i_0 = i_S+1+ii_0\) for some \(ii_0 \in [0, m_x]\)). During the FDTD update of this node, the \(dB_z/dx\) operator is calculated:
We split the outer sum over \(ii\) into two parts:
The first sum over \(ii \in [0, i_0-i_S-2]\) only uses \(B_z\) grid values in the TF region (the minimal index in \(x\) used is \(B_z\rvert_{i_S+3/2, j+jj+1/2, k+kk}^{n+1/2}\) for \(ii = i_0-i_S-2\)). Note that if \(i_0-i_S-2 < 0\), this sum has no terms and is equal to 0; the same applies for the following sums. Since the \(E_y\) value in question is also a TF, these terms do not require any action by incident field solver. The remaining sum over \(ii \in [i_0-i_S-1, m_x]\) contains differences of a TF value and an SF value. Each of the latter ones requires a term in the incident field solver update of \(E_y\rvert_{i_0, j+1/2, k}^{n+1}\).
Performing the same kind of substitution and regrouping demonstrated above for the standard Yee solver yields
Thus, we apply the following update for each grid value \(E_y\rvert_{i_0, j+1/2, k}^{n+1} \in E_{TF}\):
For values in SF the treatment is similar. For a node \(E_y\rvert_{i_0, j+1/2, k}^{n+1} \in E_{SF}\) (\(i_0 = i_S+1-ii_0\) for some \(ii_0 \in [1, m_x]\)) we apply \(dB_z/dx\) operator and split the outer sum the same way:
The first sum only has values in the SF region, and the second sum contains differences of TF and SF values. Note that now \(E_y\rvert_{i_0, j+1/2, k}^{n+1}\) is in the SF region and so we express the whole update as for SF:
Thus, we apply the following update for each grid value \(E_y\rvert_{i_0, j+1/2, k}^{n+1} \in E_{SF}\):
Other field components, axes and directions are treated in a similar way.
Example: 4th Order FDTD
For example, consider the 4th order FDTD and \(x_{min}\) boundary. Its derivative operator has \(m_x = 1\), \(m_y = m_z = 0\), \(\alpha_{0, 0, 0} = 27/24\), \(\alpha_{1, 0, 0} = -1/24\). Three layers of \(E_y\) are updated, the first in the SF region and the latter two in the TF region:
Updates of \(E_z\) are done in a similar fashion:
Three layers of \(B_y\) are updated, the first two in the SF region and the last one in the TF region:
Finally, updates of \(B_z\) are as follows:
Calculating Incident B from E
Consider a case when both \(\vec E^{inc}(x, y, z, t)\) and \(\vec B^{inc}(x, y, z, t)\) are theoretically present, but only \(\vec E^{inc}(x, y, z, t)\) is known in explicit form.
When slowly varying elvelope approximation (SVEA) is applicable, one may employ it to calculate the other field as \(\vec B^{inc}(x, y, z, t) = \vec k \times \vec E^{inc}(x, y, z, t) / c\).
Otherwise one may use TF/SF with only the modified known field set as incident and the other one set to 0.
Generally, the interpretation of the result is assisted by the equivalence theorem, and in particular Love and Schelkunoff equivalence principles [Harrington2001] [Balanis2012].
Having \(\vec E^{inc}(x, y, z, t) = \vec 0\) means only electric current \(\vec J\) would be impressed on \(S\).
Taking into account no incident fields in the SF region, the region is effectively a perfect magnetic conductor.
Likewise, having \(\vec B^{inc}(x, y, z, t) = \vec 0\) corresponds to only magnetic current and effectively a perfect electric conductor in the SF region.
Practically, one may try using the known incident field with twice the amplitude and keeping the other incident field zero, as demonstrated in [Rengarajan2000].
Note that using this approach in PIConGPU results in generating pulses going both inwards and outwards of the Huygens surface
(similar to laser profiles with initPlaneY > 0).
Therefore, in this case it is recommended to have no density outside the surface and use a strong field absorber to negate the effect of the artificial outwards-going pulse.
References
M. Potter, J.-P. Berenger A Review of the Total Field/Scattered Field Technique for the FDTD Method FERMAT, Volume 19, Article 1 (2017)
A. Taflove, S.C. Hagness Computational electrodynamics: the finite-difference time-domain method Artech house (2005)
R.F. Harrington Time-Harmonic Electromagnetic Fields McGraw-Hill (2001)
C.A. Balanis Advanced Engineering Electromagnetics John Wiley & Sons (2012)
S.R. Rengarajan, Y. Rahmat-Samii The field equivalence principle: illustration of the establishment of the non-intuitive null fields IEEE Antennas and Propagation Magazine, Volume 42, No. 4 (2000)