## 4 Parabolic Equation Modeling
A widely used solution algorithm is the Split-Step Fourier Method (SSFM) [43], an efficient numerical technique for cases where the refractive index \(n(x,z)\) varies slowly. The core idea of the method is to split the original partial differential equation into two parts. One part describes diffraction phenomena, while the other accounts for variations in the medium, specifically the spatial changes in the refractive index. We assume that \(\psi(x,z)\) is sampled on a discrete grid, \(x\) represents the propagation direction, and \(z\) represents the transverse direction. Firstly, for the refractive index step, we solve  
\[\
abla\times\mathbf{E}=-j\omega\mu\mathbf{H}, \tag{7}\]3  
The exact solution over a small step \(\Delta x\) is \(\psi^{\prime}=\psi e^{-j\frac{k_{0}}{2}\left(n^{2}-1\right)\Delta x}\), which is an element-wise multiplication, thus, computationally efficient. Secondly, the equation governing diffraction is as follows:  
\[\
abla\times\mathbf{E}=-j\omega\mu\mathbf{H}, \tag{7}\]4  
Using the Fourier transform, the second derivative with respect to \(z\) takes the form of  
\[\
abla\times\mathbf{E}=-j\omega\mu\mathbf{H}, \tag{7}\]5  
thus, in the frequency domain, the equation becomes:  
\[\
abla\times\mathbf{E}=-j\omega\mu\mathbf{H}, \tag{7}\]6  
This has the exact solution:  
\[\
abla\times\mathbf{E}=-j\omega\mu\mathbf{H}, \tag{7}\]7  
By applying the inverse Fourier transform, the original wave function is obtained effectively: \(\psi(x,z)=\mathcal{F}^{-1}\left\{\hat{\psi}^{\prime}\right\}\). Thanks to the Fast Fourier Transform (FFT) [44], this approach enables a fast and efficient solution. A widely used approach for solving the Parabolic Wave Equation is the Crank-Nicolson Finite-Difference Method (CN-FDM) [45]. This method is regarded as more computationally expensive than SSFM but is numerically stable and effectively handles the effects of rapid spatial variations in the refractive index. The core idea of the method is to approximate the derivatives in the PWE using finite differences. For better understanding, Figure 12 provides an overview of the computational domain's geometry, the introduced variables, and the types of boundary conditions considered in the application of CN-FDM. The first derivative with respect to \(x\) in (20) is approximated using the central difference scheme as follows:  
\[\
abla\times\mathbf{E}=-j\omega\mu\mathbf{H}, \tag{7}\]8where \(\Delta x\) is the step side in the \(x\)-direction, the indexes \(m+1\) and \(m\) refer to the current and next step in \(x\), and \(\partial^{2}\psi/\partial z^{2}\) is the second derivative in \(z\), for which the central difference scheme is:  
\[\
abla\times\mathbf{E}=-j\omega\mu\mathbf{H}, \tag{7}\]9  
where \(\psi_{j}^{m}\) is the wave function at position \(z_{j}\) at step \(m\), and \(\Delta z\) is the transverse grid spacing. By applying this to \(\psi^{m}\) and \(\psi^{m+1}\) as well, we obtain:  
\[\
abla\cdot\mathbf{H}=0,\quad\
abla\cdot\mathbf{E}=0. \tag{8}\]0  
Substituting (28) into (26) and rearranging the terms, a tridiagonal system of equations can be obtained:  
\[\
abla\cdot\mathbf{H}=0,\quad\
abla\cdot\mathbf{E}=0. \tag{8}\]1  
where \(I\) is the identity matrix, and \(D_{z}^{2}\) is the finite-difference operator for the second derivative in \(z\):  
\[D_{z}^{2}=\frac{1}{\Delta z^{2}} \begin{bmatrix}-2&1&0&0&\dots&0\\ 1&-2&1&0&\dots&0\\ 0&1&-2&1&\dots&0\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&\dots&1&-2&1\\ 0&0&\dots&0&1&-2\end{bmatrix} . \tag{30}\]  
Therefore, the tridiagonal system, which should be solved at each step of \(x\), can be written as \(A\psi^{m+1}=B\psi^{m}\), where \(A\) is the left-hand side matrix containing \(D_{z}^{2}\) and \(I\), \(B\) is the right-hand side matrix containing also \(D_{z}^{2}\) but with opposite signs, \(\psi^{m}\) is known, and \(\psi^{m+1}\) is unknown [46; 47]. Since \(A\) is tridiagonal, the system can be efficiently solved by, e.g., LU [48] decomposition or with sparse matrix equation solvers [49; 50].