## 8 Numerical Simulations  
This section is dedicated to the presentation of some numerical simulations, obtained by applying the schemes described in Section 6. The computational domain is the square \(Q=[0,2L]\times[-L,L]\), where \(L>0\) is a parameter. We chose by default \(L=1\). The first variable is related to \(r\) and the second one to \(z\), so that we are actually operating on a cylinder with a vertical axis. In \(Q\), we consider a squared grid, spaced by \(h:=1/N\), where \(N\) is an integer. The total number of nodes is then proportional to \(N^{2}\). The speed of light in a vacuum is normalized to \(1\). The development of the process is studied in the time interval \([0,T]\), for a certain \(T>0\). In order to guarantee stability, the time-step \(\Delta t\) will be chosen to be suitably small. We set up in advance the velocity field \(\mathbf{V}\) in order to have \(|\mathbf{V}|=1/\sqrt{\varepsilon_{r}}<1\) inside the dielectric, where we suppose that the vectors are emanating from the origin. We have \(|\mathbf{V}|=1\) on the remaining part of the domain where the vectors are horizontal (perpendicular to the \(z\)-axis). The deployment of the velocity field is displayed in Figure 5, for the cases when the dielectric is missing (left) and when it is present (right). In a more rigorous framework, \(|\mathbf{V}|\) should not be fixed in advance, but computed as time evolves according to Equations (18) and (19). As we specified in Section 6, this approach is not easily implementable from the numerical viewpoint, especially when the time advancing procedure is the trivial Euler method. Numerical solutions, in this circumstance, are too much affected by viscosity, so that the results are not reliable. Since the development of more sophisticated techniques is outside of the goals of this paper, we do not insist further on these generalizations. For the above reasons, we decided to assume that \(|\mathbf{V}|\) is known a priori. As we are interested in the evolution within the region determined by the two cones of the antenna, suitable constraints are placed in the computational domain. For simplicity, we assume that these guides are perfect conductors, though the situation might be not realistic. Additionally, we should distinguish between the nodes inside the dielectric and those outside. In the transition between the two media the fields are updated by enforcing that \(\mathbf{B}\) remains continuous, whereas \(\mathbf{E}\) has a  
Figure 4: Shapes of \(\partial\Omega\) for different values of the relative permittivity coefficient: from left to right, \(\varepsilon_{r}=2,3,4,5,6,7\). jump (these are consequences of the Fresnel's conditions; see, e.g., (28), p. 40). In practice, we truncate to zero the longitudinal component \(E_{r}\) at the first node just outside the dielectric. As an initial guess, we take a portion of spherical wave. Based on (24), we set \(g(\theta)=(\sin\theta)^{2}\) for \(\pi/4\leq\theta\leq 3\pi/4\). This modulates the displacement along the transverse direction. Concerning the longitudinal direction, we take \(f(R)=(\cos\omega R)^{2}\) with \(\omega=7\) and \(\pi/2\omega\leq R\leq 3\pi/2\omega\). The resulting electric field is shown in Figure 6 (left). The corresponding shapes of \(|\mathbf{E}|\) can be seen in the very first plot of Figures 7 and 8. Note that this displacement has non-zero divergence (\(\rho\
eq 0\)). The corresponding initial magnetic field is orthogonal to the plane of the figure and we have \(|\mathbf{B}|=|\mathbf{E}|\). Successively, the wave evolves according to the model equations. The single initial bump displays a wave-length comparable to the dimensions of the antenna. Of course, one may also consider a train of pulses lined up one after the other. This option provides however results perfectly similar to those pertaining to the single pulse. The boundaries are placed at a 45 degrees angle; thus, the imposition of the constraints comes very naturally. Neumann boundary conditions for the electric field are imposed on the antenna guide.