## 4 Parabolic Equation Modeling  
Parabolic equation modeling is a widely used method for the numerical study of electromagnetic wave propagation. The application of this method dates back to the 1940s, and with the development of numerical techniques and digital elevation models of the Earth [38], it has become one of the most important tools for solving electromagnetic wave propagation problems. This method can be considered a deterministic model because it is based on an approximate numerical solution of [PERSON]'s equations [14]. Consider the frequency-domain forms of the equations \((\partial/\partial t=j\omega t)\) and assuming a homogeneous isotropic medium \((\mathbf{D}=\epsilon\mathbf{E}\) and \(\mathbf{B}=\mu\mathbf{H})\), and also assuming that there is no free current source \((\mathbf{J}=0)\) and no charge density \((\rho=0)\) in the propagation medium, the following simplified equations can be written [39]:  
\[\
abla\times\mathbf{H}=j\omega\epsilon\mathbf{E}, \tag{6}\]  
\[\
abla\times\mathbf{E}=-j\omega\mu\mathbf{H}, \tag{7}\]  
\[\
abla\cdot\mathbf{H}=0,\quad\
abla\cdot\mathbf{E}=0. \tag{8}\]  
By taking the curl of both sides of (7) and substituting (6) into the obtained equation, we get:  
\[\
abla\times(\
abla\times\mathbf{E})=-j\omega\mu(\
abla\times\mathbf{H})=-j \omega\mu(j\omega\epsilon\mathbf{E})=\omega^{2}\mu\epsilon\mathbf{E}. \tag{9}\]  
Applying the identity  
\[\
abla\times(\
abla\times\mathbf{E})=\
abla(\
abla\cdot\mathbf{E})-\Delta \mathbf{E} \tag{10}\]  
and considering (8), we get  
\[\Delta\mathbf{E}+\omega^{2}\mu\epsilon\mathbf{E}=0. \tag{11}\]  
By defining the wave number as \(k=\omega\sqrt{\mu\epsilon}\), we get the so-called Helmholtz equation [40] for the electric field strength \(\mathbf{E}\) as  
\[\Delta\mathbf{E}+k^{2}\mathbf{E}=0. \tag{12}\]  
If the electromagnetic wave propagates not in a vacuum but in an inhomogeneous medium, the wave number varies in space and depends on the refractive index:  
\[k(x,y,z)=k_{0}n(x,y,z), \tag{13}\]  
where \(k_{0}\) is the wave number in vacuum and \(n(x,y,z)=\sqrt{\epsilon_{r}\mu_{r}}\) is the refracting index in function of the relative permittivity and permeability. In this case, the Helmholtz equation can be written as  
\[\Delta\mathbf{E}+k_{0}^{2}n^{2}\mathbf{E}=0. \tag{14}\]  
A similar equation can be obtained for the magnetic field strength if we start from (6). For horizontal polarization, \(\mathbf{E}\) has only one non-zero component \(E_{y}(x,z)\), thus the Helmholtz equation can be written as  
\[\frac{\partial^{2}E_{y}}{\partial x^{2}}+\frac{\partial^{2}E_{y}}{\partial z ^{2}}+k_{0}^{2}n^{2}E_{y}=0. \tag{15}\]Considering that most of the wave energy propagates in the forward (positive \(x\)) direction and that the field varies slowly in the transverse direction, it is common to separate the fast-oscillating and slowly varying components of the electric field. Therefore, a reduced field representation is introduced in the following form:  
\[\
abla\times\mathbf{E}=-j\omega\mu\mathbf{H}, \tag{7}\]0  
where \(\psi(x,z)\) is a slowly varying envelope function and \(n_{0}\) is a constant reference refractive index, e.g., \(n_{0}=1\) for air. Two examples of this type of function can be seen in Figure 11, where the difference between the transverse beam widths, the attenuation coefficients, and the wave numbers can be clearly observed, with the same exponential character. These demonstrate how the wave varies more rapidly along the propagation axis than in the transverse direction, highlighting the paraxial nature of the approximation. Substituting (16) back into (15) and determining the derivatives gives:  
\[\
abla\times\mathbf{E}=-j\omega\mu\mathbf{H}, \tag{7}\]1 \[\
abla\times\mathbf{E}=-j\omega\mu\mathbf{H}, \tag{7}\]2  
This equation is referred to as the Parabolic Wave Equation (PWE) [41], which is solved numerically in simulations conducted for wave propagation studies [42]. Among the numerous available algorithms, this paper presents two fundamentally different types.