## 2 Materials and Methods
### Numerical Modelling  
Equations (19)-(21) are discretised using the classical FDTD method of second order in space and time, obtaining the algebraic Equations (22)-(24). In Figure 1, the calculation molecule of \(E_{z}\), \(H_{x}\), and \(H_{y}\) in space and time is depicted. \[Hx_{(i,j-1/2)}^{n+1/2}=Hx_{(i,j-1/2)}^{n-1/2}\left(\frac{\mu_{(i-1/2)}^{n-1/2}} {\mu_{(i,j-1/2)}^{n+1/2}}\right)-\frac{\Delta t}{\mu_{(i,j-1/2)}^{n+1/2}}\left( \frac{Ez_{(i,j)}^{n}-Ez_{(i,j-1)}^{n}}{\Delta y}\right)\, \tag{22}\]  
\[Hy_{(i-1/2,j)}^{n+1/2}=Hy_{(i-1/2,j)}^{n-1/2}\left(\frac{\mu_{(i-1/2,j)}^{n-1/2 }}{\mu_{(i-1/2,j)}^{n+1/2}}\right)+\left(\frac{\Delta t}{\mu_{(i-1/2,j)}^{n+1/ 2}}\right)\left(\frac{Ez_{(i,j)}^{n}-Ez_{(i-1,j)}^{n}}{\Delta x}\right) \tag{23}\]  
Figure 1: Mesh configuration for the propagation of electrical transverse modes in the \(x,y,t\) space. and  
\[Ez_{(i,j)}^{n+1}=Ez_{(i,j)}^{n}\left(\frac{4\epsilon_{(i,j)}^{n}- \Delta t\Big{(}\sigma_{(i,j)}^{n+1}+\sigma_{(i,j)}^{n}\Big{)}}{4\epsilon_{(i,j) }^{n+1}+\Delta t\Big{(}\sigma_{(i,j)}^{n+1}+\sigma_{(i,j)}^{n}\Big{)}}\right)+ \frac{\Delta t}{\epsilon_{(i,j)}^{n+1}+\frac{\Delta t}{4}\Big{(}\sigma_{(i,j)}^ {n+1}+\sigma_{(i,j)}^{n}\Big{)}} \tag{24}\] \[\left[\frac{Hy_{(i+1/2,j)}^{n-1/2}-Hy_{(i-1/2,j)}^{n-1/2}}{\Delta x }-\frac{Hx_{(i,j+1/2)}^{n-1/2}-Hx_{(i,j-1/2)}^{n-1/2}}{\Delta y}-\frac{Jz_{(i, j)}^{n+1}+Jz_{(i,j)}^{n}}{2}\right]\,,\]  
where \(i,j,n\) denote the discrete position corresponding to the coordinates \(x,y,\) and \(t\), respectively, while \(\Delta x\), \(\Delta y\), and \(\Delta t\) are the FDTD grid size in the \(x,y,\) and \(t\) variables, respectively. The time step must satisfy the Courant-Friedrichs-Lewy (CFL) criterion \(\Delta t\leq\frac{\Delta x}{c_{0}}\sqrt{\frac{\varepsilon_{r}H_{r}}{2}}\) to ensure numerical stability [52; 53], where \(c_{0}\) is the speed of light in a vacuum, and \(\varepsilon_{r}\) and \(\mu_{r}\) are the propagating medium's relative permittivity and permeability, respectively. Several numerical experiments in the context of electromagnetic propagation have shown that a relation of \(\Delta x=\Delta y=\
icefrac{{\lambda}}{{20}}\), where \(\lambda\) is the propagated wavelength, allows the quality of the numerical solution to be maintained [53]. That is, when setting the computational grid, there should be at least 20 spatial grid cells of the FDTD method for an accurate resolution of one wavelength of the signal. The absorbing boundary condition implemented is the Convolutional PML (CPML) [51] in its simplified version by [PERSON] and [PERSON] [54]; this version of the Convolutional PML saves memory and is relatively straightforward to code. To implement the CPML boundary conditions, it is necessary to introduce the Luebbers convolutional memory variables [55] for each spatial derivative. The convolutional term \(\Psi\) is introduced in Equations (22)-(24) as  
\[\Psi_{i}^{n}=\sum_{m=0}^{n-1}\frac{\partial}{\partial x}\bigg{|}_{i\Delta x}^ {\left|n-(m+\frac{1}{2})\right|\Delta t}Z(m)\;, \tag{25}\]  
where  
\[Z(m)=\left(\frac{d}{\kappa d+\alpha x^{2}}\right)\left\{e^{-\left(\frac{d}{x }+\alpha\right)(m+1)\Delta t}-e^{-\left(\frac{d}{x}+\alpha\right)m\Delta t} \right\}\,, \tag{26}\]  
and where \(\alpha(x),\kappa(x),\) and \(d(x)\) are the CPML parameters.