*2.3. The Transformation from Fictitious Wave Field to Real Di*ff*usion Field*

It should be noted that the fictitious wave field is only used for the convenience of calculation, and it does not exist in reality. Therefore, we need to transform it back to the real diffusion field. The transformations are as follows [23]:

$$J(\omega) = \int\_0^T f'(t') e^{-\sqrt{\omega \omega\_0}t'} e^{i \sqrt{\omega \omega\_0}t'} dt',\tag{25}$$

$$E(\omega) = \sqrt{\frac{-i\omega}{2\omega\_0}} \int\_0^T E'(t')e^{-\sqrt{\omega\omega\_0}t'}e^{i\sqrt{\omega\omega\_0}t'}dt',\tag{26}$$

$$H(\omega) = \int\_0^T H'(t') e^{-\sqrt{\omega \omega\_0} t'} e^{i \sqrt{\omega \omega\_0} t'} dt',\tag{27}$$

$$\mathbf{K}(\omega) = \sqrt{\frac{-i\omega}{2\omega\_0}} \int\_0^T \mathbf{K}'(t') e^{-\sqrt{\alpha\omega\_0}t'} e^{i\sqrt{\alpha\omega\_0}t'} dt',\tag{28}$$

where *T* is the total calculation time in the fictitious wave field, and *t* is the sample time.

#### *2.4. Loading Complex Frequency Shifted Perfectly Matched Layer (CFS-PML) in Fictitious Wave Field*

 Roden (2000) proposed CFS-PML (complex frequency shifted perfectly matched layer) boundary conditions, which significantly saved on memory [27]. The CFS-PML boundary conditions are introduced into the fictitious wave field to improve the absorption effect of the boundary on the low-frequency wave [28].

In x-components, Maxwell's curl equations are as follows:

$$
\nabla\_s \times H'(\omega') + i\omega' \varepsilon' E'(\omega') = 0,\tag{29}
$$

$$\nabla\_{\mathbf{s}} = \frac{1}{\mathfrak{s}\_{\mathbf{x}}} \frac{\partial}{\partial \mathbf{x}} \mathbf{i} + \frac{1}{\mathfrak{s}\_{\mathbf{y}}} \frac{\partial}{\partial \mathbf{y}} \mathbf{j} + \frac{1}{\mathfrak{s}\_{\mathbf{z}}} \frac{\partial}{\partial \mathbf{z}} \mathbf{k}. \tag{30}$$

It should be noted that *i* is the imaginary unit in Maxwell equations, and i, j, and k are the three direction vectors of the coordinate system.  s *<sup>x</sup>*, <sup>s</sup> *<sup>y</sup>*, and  <sup>s</sup> *<sup>z</sup>* are the coordinate stretching factors,  <sup>s</sup> (*xyz*) <sup>=</sup> *<sup>k</sup>*(*x*,*y*,*z*) <sup>+</sup> <sup>σ</sup>(*x*,*y*,*z*) <sup>α</sup>(*x*,*y*,*z*)+*i*ωε <sup>α</sup>(*x*,*y*,*z*) <sup>&</sup>gt; 0, <sup>σ</sup>(*x*,*y*,*z*) <sup>&</sup>gt; 0, and *<sup>k</sup>*(*x*,*y*,*z*) <sup>≥</sup> 1. By setting <sup>α</sup>(*x*,*y*,*z*) <sup>=</sup> <sup>σ</sup>(*x*,*y*,*z*) <sup>=</sup> 0, *k*(*x*,*y*,*z*) = 1, we can get the original stretched coordinate metrics, which is the PML boundary condition. Taking E *x* as an example, and using Equations (29) and (30), we can get

$$
\delta\omega^{\prime}\varepsilon^{\prime}E\_{\infty}^{\prime}(\omega^{\prime}) = -\left(\frac{1}{\mathfrak{E}\_{\mathcal{Y}}}\frac{\partial H^{\prime}\_{z}(\omega^{\prime})}{\partial \mathcal{Y}} - \frac{1}{\mathfrak{E}\_{z}}\frac{\partial H^{\prime}\_{y}(\omega^{\prime})}{\partial z}\right). \tag{31}
$$

Applying Laplace transform to (31) to the time domain, the iterative solution can be obtained.

$$E^{\prime \underline{n} + 1}\_{\underline{\chi}\_{i+1/2, \vert k}} = E^{\prime \underline{n}}\_{\underline{\chi}\_{i+1/2, \vert k}} + \Lambda t^{\prime} \frac{\Delta \underline{n}}{d} \left[ \partial\_{\underline{y}}^{-1} H^{\prime \underline{n} + 1/2}\_{\underline{\chi}\_{i+1/2, \vert k} + 1/2, k} - \partial\_{\underline{x}}^{-1} H^{\prime \underline{n} + 1/2}\_{\underline{\chi}\_{i+1/2, \vert k} + 1/2} \right] + \Lambda t^{\prime} \frac{\Delta \underline{n}}{d} (\Psi^{\underline{n} + 1/2}\_{\underline{\chi}\_{i\vert k} + 1/2, \vert k} - \Psi^{\vert \underline{n} + 1/2}\_{\underline{\chi}\_{i\vert k} + 1/2, \vert k}), \tag{32}$$

$$\mathbf{V}\_{a\_{2\overline{\eta}\_{i+1/2}},k}^{n+1/2} = e^{-(\frac{a\_{\overline{\eta}}}{h\_{\overline{\eta}}}+a\_{\overline{\eta}})\frac{\mathbf{h}\_{\overline{\eta}}^{\rho}}{\mathbf{V}}} \mathbf{V}\_{a\_{2\overline{\eta}\_{i+1/2}},k}^{n+1/2} + \frac{a\_{\overline{\eta}}}{(a\_{\overline{\eta}}k\_{\overline{\eta}} + k\_{\overline{\eta}}^{2}a\_{\overline{\eta}})} (e^{-(\frac{a\_{\overline{\eta}}}{h\_{\overline{\eta}}}+a\_{\overline{\eta}})\frac{\mathbf{h}\_{\overline{\eta}}^{\rho}}{\mathbf{V}}} - 1)(H^{n+1/2}\_{\overline{z}+1/2,\overline{\eta}+1/2,k} - H^{n+1/2}\_{\overline{z}+1/2,\overline{\eta}-1/2,k})/\Delta\_{\overline{y}},\tag{33}$$

$$\Psi\_{\varepsilon\_{\overline{x}\_{i+1/2}},\mathbb{Q}}^{u+1/2} = e^{-(\frac{\varrho\_{\varepsilon}}{\varrho\_{\varepsilon}}+a\_{\varepsilon})\frac{\Delta'}{\overline{r}}} \Psi\_{\varepsilon\_{\overline{x}\_{i+1/2}},\mathbb{Q}}^{u-1/2} + \frac{\varrho\_{\varepsilon}}{(\varrho\_{\varepsilon}k\_{\varepsilon}+k\_{\varepsilon}^{2}a\_{\varepsilon})} (e^{-(\frac{\varrho\_{\varepsilon}}{\varrho\_{\varepsilon}}+a\_{\varepsilon})\frac{\Delta'}{\overline{r}}} - 1)(H^{u+1/2}\_{\frac{\varrho\_{\varepsilon}}{\varrho\_{\varepsilon}}+1/2,\frac{\varrho\_{\varepsilon}}{\overline{r}}+1/2,\frac{\varrho\_{\varepsilon}}{\overline{r}}+1/2,\frac{\varrho\_{\varepsilon}}{\overline{r}}} - H^{u+1/2}\_{\frac{\varrho\_{\varepsilon}}{\varrho\_{\varepsilon}}+1/2,\frac{\varrho\_{\varepsilon}}{\overline{r}}}) / \wedge \Delta\_{\varepsilon}.\tag{34}$$

On the boundary parameters, σ*<sup>i</sup>* and *ki*, which depend on the relative positions of nodes and target areas, are not fixed, whereas <sup>Ψ</sup>*n*+1/2 *exyi*<sup>+</sup>1/2,*j*,*<sup>k</sup>* and <sup>Ψ</sup>*n*+1/2 *exzi*<sup>+</sup>1/2,*j*,*<sup>k</sup>* are two discrete variables.

$$\sigma\_i = \sigma\_{\text{max}} \frac{|k - k\_0|^m}{d^m}, (i = \mathbf{x}, y, z), \tag{35}$$

$$k\_i = 1 + (k\_{\text{max}} - 1) \frac{|k - k\_0|^m}{d^m}, (i = x, y, z), \tag{36}$$

where *k*<sup>0</sup> indicates the intersection of the boundary and the target area boundary, *d* is the boundary thickness, and *m* is the polynomial parameter.
