The Evaluation of an Asymptotic Solution to the Sommerfeld Radiation Problem Using an Efficient Method for the Calculation of Sommerfeld Integrals in the Spectral Domain

Abstract

A recently developed high-frequency asymptotic solution for the famous “Sommerfeld radiation problem” is revisited. The solution is based on an analysis performed in the spectral domain, through which a compact asymptotic formula describes the behavior of the EM field, which emanates from a vertical Hertzian radiating dipole, located above flat, lossy ground. The paper is divided into two parts. We first demonstrate an efficient technique for the accurate numerical calculation of the well-known Sommerfeld integrals. The results are compared against alternative calculation approaches and validated with the corresponding Norton figures for the surface wave. In the second part, we introduce the asymptotic solution and investigate its performance; we compare the solution with the accurate numerical evaluation for the received EM field and with a more basic asymptotic solution to the given problem, obtained via the application of the Stationary Phase Method. Simulations for various frequencies, distances, altitudes, and ground characteristics are illustrated and inferences for the applicability of the solution are made. Finally, special cases leading to analytical field expressions close as well as far from the interface are examined.

1. Introduction

The Sommerfeld radiation problem is a classical problem in the area of electromagnetic (EM) wave propagation. The original Sommerfeld solution is provided in the spatial domain as an integral expression, utilizing the so-called “Hertz potentials”, but it does not end up in analytic formulas [1]. With the use of the Fourier–Bessel representations in cylindrical coordinates, A. N. Sommerfeld proposed the solution for the Hertz vectors in terms of improper integrals, now known as Sommerfeld integrals. Sommerfeld further derived an asymptotic surface wave expression, assuming high media contrast and large horizontal range. The theory of surface waves was then further developed both by Sommerfeld and J. Zenneck [2,3,4] and mathematically strengthened by H. Ott [5].

Subsequently, K. A. Norton focused on the engineering application of the problem; provided approximate solutions, represented by long algebraic expressions; and described concepts such as the propagating surface wave and its associated “attenuation coefficient”, albeit his definition for the surface wave is not equivalent to Sommerfeld’s definition [6,7]. Various prominent researchers have also dealt with the problem from various points of view. A notable debate was between H. Barlow and J. Wait regarding the existence of surface waves in practical occasions [8]. However, later experimental results verified that a slowly decaying EM field near the interface that decreases approximately as of $1/\sqrt{R}$ can be excited and detected under usual circumstances, providing thus a tangible indication as for the existence of Zenneck waves [9,10]. Related important work is also found in [11,12,13,14,15,16,17,18]. A thorough review of the subject is given in [19,20].

1.1. Previous Contribution by Our Research Group

The problem may also be tackled in the spectral domain. Particularly, in [21,22], we derived the fundamental integral representations for the received EM field by means of a generalized solution to the respective Maxwell equations boundary value problem. This approach has the advantage that no Hertz potentials and their subsequent differentiation are required for the evaluation of the fields.

In [23], the Stationary Phase Method (SPM) [24,25,26] was applied to the EM field’s integral expressions and the well-known analytic formulas for the space wave, defined as the complex summation of the line of sight (LOS) field and a portion of the field emanating from the dipole’s image point (also called reflected field) were obtained as the high frequency asymptotic solution to the complete problem. In [27,28,29], we focused on the numerical evaluation of the field’s integral formulas. It was revealed that accurately evaluating the Sommerfeld integrals is not trivial. The result is sensitive on the position of the singular points in relation to the integration path, an issue that has been also a major problem and matter of debate in various related research works [12].

Then, in [30], the problem was re-examined, for the case where $\mathsf{\sigma }\gg \mathsf{\omega }{\mathsf{\epsilon }}_0$, which is valid for many practical cases of interest in terrestrial communications. As shown there, a special contour integral, called “Etalon integral”, was used to deform the original contour of integration, through the application of the Saddle Point or Steepest Descent Method (SDP). This “Etalon integral” can be expressed in terms of Fresnel integrals and has interesting properties, which can reduce the problem related to the vicinity of the saddle point to the pole point [31,32,33,34,35,36]. The result is a compact formula that better expresses the variation of the field in the high-frequency regime.

1.2. Scope of This Research

As mentioned in [28,29], the accurate evaluation of Sommerfeld integral expressions is not a straightforward task, and this is due to both the presence of singularities along the integration path as well as to the particular complex and rapidly oscillating nature of the integrands. For that reason, various specialized commercial software have been used for obtaining adequate results, for example, the AWAS tool used by Sarkar et al. in [12]. In this paper, we show that, using an appropriate variable transformation, it is possible to convert the generalized integrals of [23] into fast converging formulas, which are rather suitable for numerical calculation, using standard Numerical Integration (NI) techniques. Particularly, the integral expression, describing the received EM field, is broken down into two terms; one relatively easily computed definite integral of finite integation range and another integral of semi-infinite range. However, the latter integrand proves to be a rapidly decaying exponential function, resulting in very fast convergence times. Comparisons against the numerical results, published in [28,29], demonstrate the advantage of the method. Additionally, a validation against Norton’s figures for the well-known surface wave [6,7] is exhibited.

Then, we return to the solution of [30] and elaborate on the derived expressions. Using small and large argument approximations, associated with the Fresnel integrals [37], pure analytic expressions are extracted that better describe the behavior of the EM field, close as well as far away from the ground interface, and provide useful information regarding the propagation mechanism. Moreover, extensive simulations are demonstrated for the purpose of validating these asymptics and a comparison with the solution of [23] (SPM-based solution) is performed. Provided the required conditions are satisfied, our simulations validate the analytic expressions. Related inferences, regarding the field behavior, are made.

1.3. Problem Geometry

The geometry of the problem is shown in Figure 1. A vertical small (Hertzian) dipole (HD), characterized by dipole moment $\underline{\dot{p}}=p·{\widehat{e}}_x,p=\mathrm{const}$, is directed towards the positive x-axis, at altitude X_O above infinite, flat, lossy ground. The dipole radiates time-harmonic EM waves at angular frequency $\mathsf{\omega }=2\pi f$ ($e^{-i\mathsf{\omega }t}$ time dependence is assumed). The relative complex permittivity of the ground is as follows: ${\mathsf{\epsilon }}_r^{{}^{\prime }}={\mathsf{\epsilon }}^{{}^{\prime }}/upvarepsilo{n}_0={\mathsf{\epsilon }}_r+i\mathsf{\sigma }/\mathsf{\omega }{\mathsf{\epsilon }}_0$, where $\mathsf{\sigma }$ is the ground conductivity, f is the carrier frequency, and ${\mathsf{\epsilon }}_0=8.854\times {10}^{-12}$ F/m is the permitivity in vacuum or air. The goal is to evaluate the received field at an arbitrary observation point above the interface, i.e., at point (x,y,z), of Figure 1.

1.4. Structure of the Article

In what follows, Section 2 recaps the fundamentals expressions for the EM field in the spectral domain and the issues associated with their numerical calculation and demonstrates how a simple variable transformation leads to fast converging integral formulas, suitable for evaluation in the computer. Through various simulations, we illustrate the advantages and validate the accuracy of the redefined expressions. Then, in Section 3, we give an overview of the examined asymptotic solution, and through an extended set of simulations and comparisons, we demonstrate its efficiency. Additionally, a discussion regarding the applicability of the resulting closed-form formulas is given. Finally, in Section 4, we summarize on the major findings and propose potential extensions. The analysis is given for the electric field. Expressions for the magnetic field are derived similarly or by suitable use of the duality principle.

2. Efficient Formulation for the EM Field Integral Expressions in Spectral Domain

2.1. Spectral Domain Integral Expressions

The electric field at the observation point of Figure 1 is given by the following integral expression [23],

$$\begin{array}{cc}\hfill \underline{E}& ={\underline{E}}^{LOS}+{\underline{E}}^R=-{\displaystyle \frac{ip}{8\mathsf{\pi }{\mathsf{\epsilon }}_0{\mathsf{\epsilon }}_1}}\left[{\int }_{-\infty }^{\infty }\underline{{\mathrm{f}}_1}\left(k_{\mathsf{\rho }}\right)\mathrm{d}{k}_{\mathsf{\rho }}+{\int }_{-\infty }^{\infty }\underline{{\mathrm{f}}_2}\left(k_{\mathsf{\rho }}\right)\mathrm{d}{k}_{\mathsf{\rho }}\right],\hfill \end{array}$$

(1)

where ${\underline{E}}^{LOS}$ denotes the direct or LOS field, ${\underline{E}}^R$ is for the field scattered by the flat and lossy ground, and the vector functions $\underline{{\mathrm{f}}_1}\left(k_{\mathsf{\rho }}\right)$ and $\underline{{\mathrm{f}}_2}\left(k_{\mathsf{\rho }}\right)$ are given by

$$\begin{array}{cc}\hfill \underline{{\mathrm{f}}_1}\left(k_{\mathsf{\rho }}\right)=& \left({\mathsf{\kappa }}_1\mathrm{sgn}\left(x-x_0\right){\widehat{e}}_{\mathsf{\rho }}-\left|k_{\mathsf{\rho }}\right|{\widehat{e}}_x\right)·{\displaystyle \frac{k_{\mathsf{\rho }}\left|k_{\mathsf{\rho }}\right|}{{\mathsf{\kappa }}_1}}{\mathrm{H}}_0^{\left(1\right)}\left(k_{\mathsf{\rho }}\mathsf{\rho }\right)e^{i{\mathsf{\kappa }}_1|x-x_0|},\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \underline{{\mathrm{f}}_2}\left(k_{\mathsf{\rho }}\right)=& \left({\mathsf{\kappa }}_1{\widehat{e}}_{\mathsf{\rho }}-\left|k_{\mathsf{\rho }}\right|{\widehat{e}}_x\right)k_{\mathsf{\rho }}\left|k_{\mathsf{\rho }}\right|·{\displaystyle \frac{{\mathsf{\epsilon }}_2{\mathsf{\kappa }}_1-{\mathsf{\epsilon }}_1{\mathsf{\kappa }}_2}{{\mathsf{\kappa }}_1\left({\mathsf{\epsilon }}_2{\mathsf{\kappa }}_1+{\mathsf{\epsilon }}_1{\mathsf{\kappa }}_2\right)}}{\mathrm{H}}_0^{\left(1\right)}\left(k_{\mathsf{\rho }}\mathsf{\rho }\right)e^{i{\mathsf{\kappa }}_1(x+x_0)},\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\mathsf{\kappa }}_1=& \sqrt{k_{01}^2-k_{\mathsf{\rho }}^2},{\mathsf{\kappa }}_2=\sqrt{k_{02}^2-k_{\mathsf{\rho }}^2},\hfill \end{array}$$

(4)

with ${\mathrm{H}}_0^{\left(1\right)}$ being the Hankel function of zero order and first kind and where $k_{01}$ are $k_{02}$ the wavenumbers of propagation in the air and lossy medium (ground), respectively. Additionally, $k_{\mathsf{\rho }}$ represents the propagation wavenumber along the horizontal $\mathsf{\rho }$-direction and corresponds to Sommerfeld’s continuous eigenvalues spectrum in the Bessel integral representation of his solution [1,3,4].

Expressions (1)–(4) expose the following difficulties, when coming to the evaluation of the respective integrals through common Numerical Integration (NI) techniques:

–

The range of integration extends from $-\infty$ to $+\infty$, resulting in potential computational errors for large evaluation arguments.

–

The Hankel function, ${\mathrm{H}}_0^{\left(1\right)}$, exhibits a singularity at $k_{\mathsf{\rho }}=0$, and although it is proven that this is a logarithmic singularity [38] and does not break the integral’s convergence (since one can easily show: ${\displaystyle \underset{k_{\mathsf{\rho }}\to 0}{lim}{k}_{\mathsf{\rho }}·{\mathrm{H}}_0^{\left(1\right)}\left(k_{\mathsf{\rho }}\mathsf{\rho }\right)=0}$ [28]), it can affect the accuracy of the NI results.

–

Points $k_{\mathsf{\rho }}=\pm k_{01}$ are also isolated singularities, and despite them being still integrable singularities [38] (it is a square root integrable singularity that applies to Rule 1 of [38]), a sufficiently small range around them must be excluded when evaluating (1) in the computer. As argued in [29], doing so may severely affect the accuracy of the results.

The above accuracy issues are important only regarding the scattered field, ${\underline{E}}^R$, for which no analytic formula exists. For the LOS field, a closed-form expression does exist, which with respect to Figure 1 takes the following full form (near and far field components) in the cylidrical coordinate system [39], with $\mathsf{\zeta }=\sqrt{{\mathsf{\mu }}_1/{\mathsf{\epsilon }}_1}\simeq 120\mathsf{\pi }\mathsf{\Omega }$ being the wave impedance of free space or air; it is used for verification purposes, as shown in Section 2.3 below.

$$\begin{array}{cc}\hfill {\underline{E}}^{LOS}=-{\displaystyle \frac{i\mathsf{\omega }p}{4\mathsf{\pi }}}·e^{i{k}_{01}r}& \left\{\left({\displaystyle \frac{-i\mathsf{\omega }{\mathsf{\mu }}_1}{2r_1}}+{\displaystyle \frac{3\mathsf{\zeta }}{2r_1^2}}-{\displaystyle \frac{3}{2i\mathsf{\omega }{\mathsf{\epsilon }}_1{r}_1^3}}\right)\right.sin2{\mathsf{\theta }}_1·{\widehat{e}}_{\mathsf{\rho }}+\hfill \\ & +\left.\left[{\displaystyle \frac{i\mathsf{\omega }{\mathsf{\mu }}_1}{r_1}}{sin}^2{\mathsf{\theta }}_1+\left({\displaystyle \frac{\mathsf{\zeta }}{r_1^2}}-{\displaystyle \frac{2}{i\mathsf{\omega }{\mathsf{\epsilon }}_1{r}_1^3}}\right)\left(cos2{\mathsf{\theta }}_1+{cos}^2{\mathsf{\theta }}_1\right)\right]{\widehat{e}}_x\right\},\hfill \end{array}$$

(5)

2.2. Reformulated Integral Expressions for the EM Field

We now focus on the scattered field, i.e., the second integral expression of (1), which may be written as

$$\begin{array}{c}\hfill {\underline{E}}^R=-{\displaystyle \frac{ip}{8\mathsf{\pi }{\mathsf{\epsilon }}_0{\mathsf{\epsilon }}_1}}(I_1+\stackrel{I_{23}}{\overbrace{I_2+I_3}}),\end{array}$$

(6)

$$\begin{array}{c}\hfill I_1={\int }_{-k_{01}}^{+k_{01}}\underline{{\mathrm{f}}_2}\left(k_{\mathsf{\rho }}\right)\mathrm{d}{k}_{\mathsf{\rho }},\end{array}$$

(7a)

$$\begin{array}{c}\hfill I_2={\int }_{-\infty }^{-k_{01}}\underline{{\mathrm{f}}_2}\left(k_{\mathsf{\rho }}\right)\mathrm{d}{k}_{\mathsf{\rho }},\end{array}$$

(7b)

$$\begin{array}{c}\hfill I_3={\int }_{+k_{01}}^{+\infty }\underline{{\mathrm{f}}_2}\left(k_{\mathsf{\rho }}\right)\mathrm{d}{k}_{\mathsf{\rho }}.\end{array}$$

(7c)

Starting with (7a), we perform a simple variable transform, $k_{\mathsf{\rho }}=k_{01}sin\xi$, which apparently maps the $\left[-k_{01},+k_{01}\right]$ range to $\left[-\mathsf{\pi }/2,+\mathsf{\pi }/2\right]$. A related variable change is given in [40,41] or in [42] for the horizontal dipole case. However, these are of a different nature and lead to different inferences, as obvious from the expressions that follow, herein. With this transform, (4) is translated to

$${\mathsf{\kappa }}_1=k_{01}cos\xi ,{\mathsf{\kappa }}_2=\sqrt{k_{02}^2-k_{01}^2{sin}^2\xi }.$$

(8)

If we also take into consideration the definition for $\underline{{\mathrm{f}}_2}$, as given by (3), the expression for $I_1$ becomes

$$\begin{array}{c}{I}_1=k_{01}^3{\int }_{-\frac{\mathsf{\pi }}{2}}^{+\frac{\mathsf{\pi }}{2}}\left(cos\xi {\widehat{e}}_{\mathsf{\rho }}-|sin\xi |{\widehat{e}}_x\right)·sin\xi |sin\xi |·R_{\parallel }\left(\xi \right)·{\mathrm{H}}_0^{\left(1\right)}\left(\mathsf{\rho }{k}_{01}sin\xi \right)·e^{i{k}_{01}(x+x_0)cos\xi }\mathrm{d}\xi ,\hfill \\ \hfill \end{array}$$

(9)

$$\begin{array}{c}{R}_{\Vert }\left(\xi \right)=\frac{{\mathsf{\epsilon }}_2{k}_{01}cos\xi -{\mathsf{\epsilon }}_1\sqrt{k_{02}^2-k_{01}^2{sin}^2\xi }}{{\mathsf{\epsilon }}_2{k}_{01}cos\xi +{\mathsf{\epsilon }}_1\sqrt{k_{02}^2-k_{01}^2{sin}^2\xi }},\hfill \end{array}$$

(10)

where $R_{\Vert }$ the reflection coefficient. Equation (10) may be equivalently written as

$$\begin{array}{cc}\hfill I_1& =k_{01}^3\left\{{\int }_0^{\frac{\mathsf{\pi }}{2}}\left[\left(cos\xi {\widehat{e}}_{\mathsf{\rho }}-sin\xi {\widehat{e}}_x\right)·{sin}^2\xi ·R_{\parallel }\left(\xi \right)·{\mathrm{H}}_0^{\left(1\right)}\left(\mathsf{\rho }{k}_{01}sin\xi \right)·e^{i{k}_{01}(x+x_0)cos\xi }\right]\mathrm{d}\xi -\right.\hfill \\ \hfill -& \left.{\int }_{-\frac{\mathsf{\pi }}{2}}^0\left[\left(cos\xi {\widehat{e}}_{\mathsf{\rho }}+sin\xi {\widehat{e}}_x\right)·{sin}^2\xi ·R_{\parallel }\left(\xi \right)·{\mathrm{H}}_0^{\left(1\right)}\left(\mathsf{\rho }{k}_{01}sin\xi \right)·e^{i{k}_{01}(x+x_0)cos\xi }\right]\mathrm{d}\xi \right\}.\hfill \end{array}$$

(11)

We may further elaborate on (11), if we make use of the following properties for the Hankel function [43,44],

$$\begin{array}{cc}\hfill {\mathrm{H}}_0^{\left(1\right)}\left(z\right)+{\mathrm{H}}_0^{\left(2\right)}\left(z\right)& =\phantom{-}2{\mathrm{J}}_0\left(z\right),\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\mathrm{H}}_0^{\left(1\right)}\left(z{e}^{i\mathsf{\pi }}\right)& =-{\mathrm{H}}_0^{\left(2\right)}\left(z\right),\hfill \end{array}$$

(13)

(with the latter implying an analytic continuation of ${\mathrm{H}}_0^{\left(1\right)}$in the upper half plane) and observe from (10) that the reflection coefficient $R_{\Vert }\left(\xi \right)$ is an even function, with respect to $\xi$. Overall, we get

$$\begin{array}{c}\hfill I_1=2k_{01}^3{\int }_0^{\frac{\mathsf{\pi }}{2}}\left[\left(cos\xi {\widehat{e}}_{\mathsf{\rho }}-sin\xi {\widehat{e}}_x\right)·{sin}^2\xi ·R_{\parallel }\left(\xi \right)·{\mathrm{J}}_0\left(\mathsf{\rho }{k}_{01}sin\xi \right)·e^{i{k}_{01}(x+x_0)cos\xi }\right]\mathrm{d}\xi ,\end{array}$$

(14)

where ${\mathrm{J}}_0$ denotes the zero-order Bessel function.

For integrals $I_2$ and $I_3$, we follow a similar approach. Particularly, in (7b), we apply the variable transform $k_{\mathsf{\rho }}=k_{01}cosh\xi$, while in (7c), we set $k_{\mathsf{\rho }}=-k_{01}cosh\xi$. In both cases, the original ranges of integration, $[-\infty ,-k_{01}]$ and $[k_{01},+\infty ]$, are mapped to $[0,+\infty ]$. Moreover, (4) becomes

$${\mathsf{\kappa }}_1=i{k}_{01}sinh\xi ,{\mathsf{\kappa }}_2=\sqrt{k_{02}^2-k_{01}^2{cosh}^2\xi }.$$

(15)

Performing the necessary calculations and using (12) and (13), we combine the results for $I_2$ and $I_3$ as

$$\begin{array}{c}\hfill I_{23}=\frac{2k_{01}^3}{i}{\int }_0^{\infty }\left(isinh\xi {\widehat{e}}_{\mathsf{\rho }}-cosh\xi {\widehat{e}}_x\right)·{cosh}^2\xi ·R_{\parallel }^{\prime }\left(\xi \right)·{\mathrm{J}}_0\left(\mathsf{\rho }{k}_{01}cosh\xi \right)·e^{-k_{01}(x+x_0)sinh\xi }\mathrm{d}\xi ,\end{array}$$

(16)

where $I_{23}=I_2+I_3$ and the reflection coefficient $R_{\Vert }^{\prime }$ is now given by

$$R_{\Vert }^{\prime }\left(\xi \right)=\frac{i{\mathsf{\epsilon }}_2{k}_{01}sinh\xi -{\mathsf{\epsilon }}_1\sqrt{k_{02}^2-k_{01}^2{cosh}^2\xi }}{i{\mathsf{\epsilon }}_2{k}_{01}sinh\xi +{\mathsf{\epsilon }}_1\sqrt{k_{02}^2-k_{01}^2{cosh}^2\xi }}.$$

(17)

Substituting (14) and (16) to (6), we reach the integral formula for the scattered field, ${\underline{E}}^R$. Working similarly for the first integral of (1), we obtain the equivalent expression for the LOS field. Overall, we get

$$\begin{array}{cc}\hfill & {\underline{E}}^{LOS}=\frac{-ip{k}_{01}^3}{4\mathsf{\pi }{\mathsf{\epsilon }}_0{\mathsf{\epsilon }}_1}·\left\{{\int }_0^{\frac{\mathsf{\pi }}{2}}\left[\left(\mathrm{sgn}(x-x_0)·cos\xi {\widehat{e}}_{\mathsf{\rho }}-sin\xi {\widehat{e}}_x\right)·{sin}^2\xi ·{\mathrm{J}}_0\left(\mathsf{\rho }{k}_{01}sin\xi \right)·e^{i{k}_{01}|x-x_0|cos\xi }\right]\mathrm{d}\xi -\right.\hfill \\ \hfill & -i{\int }_0^{{}^{\infty }}\left[\left(i\mathrm{sgn}(x-x_0)·sinh\xi {\widehat{e}}_{\mathsf{\rho }}-cosh\xi {\widehat{e}}_x\right)·{cosh}^2\xi ·{\mathrm{J}}_0\left(\mathsf{\rho }{k}_{01}cosh\xi \right)·e^{-k_{01}|x-x_0|sinh\xi }\right]\mathrm{d}\xi \left\}\phantom{\frac{-ip{k}_{01}^3}{4\mathsf{\pi }{\mathsf{\epsilon }}_0{\mathsf{\epsilon }}_1}}\right.\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill & {\underline{E}}^R=\frac{-ip{k}_{01}^3}{4\mathsf{\pi }{\mathsf{\epsilon }}_0{\mathsf{\epsilon }}_1}·\left\{{\int }_0^{\frac{\mathsf{\pi }}{2}}\left(cos\xi {\widehat{e}}_{\mathsf{\rho }}-sin\xi {\widehat{e}}_x\right)·{sin}^2\xi ·R_{\parallel }\left(\xi \right)·{\mathrm{J}}_0\left(\mathsf{\rho }{k}_{01}sin\xi \right)·e^{i{k}_{01}(x+x_0)cos\xi }\mathrm{d}\xi -\phantom{\frac{-ip{k}_{01}^3}{4\mathsf{\pi }{\mathsf{\epsilon }}_0{\mathsf{\epsilon }}_1}}\right.\hfill \\ \hfill & -i{\int }_0^{{}^{\infty }}\left(isinh\xi {\widehat{e}}_{\mathsf{\rho }}-cosh\xi {\widehat{e}}_x\right)·{cosh}^2\xi ·R_{\Vert }^{\prime }\left(\xi \right)·{\mathrm{J}}_0\left(\mathsf{\rho }{k}_{01}cosh\xi \right)·e^{-k_{01}(x+x_0)sinh\xi }\mathrm{d}\xi \left\}\phantom{\frac{k_1^3}{2}}\right.,\hfill \end{array}$$

(19)

An inspection of (18) and (19) provides useful insights. As required by the problem’s geometry, the LOS field is cylindrically symmetrical (no $\mathsf{\varphi }$–component), and it is expressed as a complex summation of contributions, originating from the dipole’s location, hence its dependence on the horizontal distance, $\mathsf{\rho }$ and the relative height difference, $x-x_0$. Additionally, the x–component of the field is symmetrical, while the $\mathsf{\rho }$–component is antisymetrical, with respect to the dipole’s vertical position $x_0$, in accordance to the conventional solution of the dipole’s problem, as evident from (5) (i.e. $E_{\mathrm{x}}\left(x_0-h\right)=E_x\left(x_0+h\right),E_{\mathsf{\rho }}\left(x_0-h\right)=-E_{\mathsf{\rho }}\left(x_0+h\right),\forall \mathsf{\rho },a$). The expression for ${\underline{E}}^R$ is similar and can be considered the integral generalization of Fresnel’s theory due to the existence of $R_{\Vert }\left(\xi \right)$ and $R_{\Vert }^{\prime }\left(\xi \right)$ in (19) that act as reflection coefficients. Additionally, the field depends on the cumulative distance, $x+x_0$, as if the source is located at the image point ${\mathrm{A}}^{\prime }$ of Figure 1.

Equations (18) and (19) remedy the accuracy issues, mentioned in Section 2.1, above:

–

They utilize the Bessel function, ${\mathrm{J}}_0$, instead of ${\mathrm{H}}_{\phantom{(}0}^{\left(1\right)}$ in (1), which is a bounded function with no singularities.

–

The singularities at points $k_{\mathsf{\rho }}=\pm k_{01}$ have also been removed. Thus, no need to exclude any range around them is required, when using any kind of numerical integration technique, in order to calculate (18) and (19).

–

The result is expressed as the sum of two integrals: one definite integral, in the bounded range $[0,\mathsf{\pi }/2]$ and an improper integral, in which the integration range extends from 0 to ∞. However, due to the presence of $e^{-k_{01}(x+x_0)sinh\xi }$, the second integrand is a fast decaying function, practically making the integral a bound-limits one that quickly converges and is easy to evaluate on a computer.

The above findings are visible in the simulations that follow. The simulation parameters (i.e., T-R heights, ground parameters, operating frequency, etc.) are indicated within the figures.

2.3. Simulation Results and Comparisons

As deduced in [29], SPM results are expected to be accurate in the far field, i.e., at least at distances over 10–15 wavelengths or above 100–150 m, for the 30 MHz case, shown in Figure 2. Therefore, using the SPM data as reference, it is obvious that only the numerical evaluation of (19) achieves the required accuracy and this is noticeably evident for distances larger than the characteristic distance of the so-called Pseudo–Brewster angle, defined as the angle of incidence, ${\mathsf{\theta }}_{\mathrm{B}}$, where the reflected field is minimized [39] (this coincidence is not a general conclusion for every tested scenario; though it is a frequent case, which is why it is mentioned here.) On the contrary, the numerical computation of (1)–(3) fails to describe the electric field behavior, which may be attributed to the reasons analyzed in Section 2.1, above.

In Figure 3, we demonstrate various field types and components for a radiating dipole at 300 KHz, which is regarded as the frontier between the low frequency (LF) and medium frequency (MF) bands [39]. For the LOS field, we used (5), while the space wave was evaluated as in [12], i.e., by using the concept of the Fresnel reflection coefficient for the reflected field. The scattered field (the terms scattered field and reflected field are not equivalent [12]) was numerically computed via (19).

Due to the small antenna heights and the long distances involved (10–20 km), the space wave is expected to diminish [39]. Therefore, the link is established primarily by means of the surface wave, defined as the remaining field, after subtracting the geometrical optics field (or space wave) from the total field [45]. This is actually verified in the top plot of Figure 3, with the total field curve being very close to the surface wave one. As a confirmation of the validity of the results, our surface wave calculations are compared with Norton formulas [6]. The respective curves are almost identical.

The bottom half of Figure 3 displays the behavior of the integrand associated to the second integral of (19), i.e., the generalized integral over the $[0,\infty )$ range. Actually, we deal only with the x–component of this integrand, which is the major field component for this problem [12], denoted as function 𝚐${}_{\mathrm{ex}}\left(\xi \right)$ in Figure 3 (the $\mathsf{\rho }$–component behaves similarly). The integrand is confined in a small window of the integration variable, $\xi$, outside of which and especially for large values of $\xi$, it is practically zeroed. This is due to the fact that $e^{-k_{01}(x+x_0)sinh\xi }$ decreases much faster than the increase rate of ${cosh}^3\xi$ (or alternatively ${lim}_{\xi \to \infty }{e}^{-k_{01}(x+x_0)sinh\xi }·{cosh}^3\xi =0$). The bottom line is that the second integral of (19) essentially becomes a bound-limits integral, which can be easily and accurately evaluated, using common numerical integration techniques. This is an advantage of our formulation.

It is interesting to note the fluctuations of 𝚐${}_{\mathrm{ex}}\left(\xi \right)$. These are an outcome of the oscillating nature of ${\mathrm{J}}_0$. Its effect on 𝚐${}_{\mathrm{ex}}\left(\xi \right)$ is apparent by observing the bold line of Figure 3 (𝚐′${}_{\mathrm{ex}}\left(\xi \right)$ in the figure), which demonstrates how the integrand would behave, if it had not been for ${\mathrm{J}}_0$. Again, the confinement of the integrand within a “narrow-band” of variable $\xi$ is apparent. It also seems that 𝚐′${}_{\mathrm{ex}}\left(\xi \right)$ acts like a slightly shifted envelope function of 𝚐${}_{\mathrm{ex}}\left(\xi \right)$. However, notice that this is a normalized illustration of 𝚐′${}_{\mathrm{ex}}\left(\xi \right)$ to the respective magnitude of 𝚐${}_{\mathrm{ex}}\left(\xi \right)$; the order of magnitude between them is totally different.

The simulations in Figure 3 are now repeated for the VHF/UHF band, in Figure 4. The source and observation points are placed even closer to the ground in an attempt to detect meaningful surface wave values in this higher frequency scenario. Nevertheless, in this case, the space wave almost completely dictates the field behavior. The pursued surface wave becomes quickly negligible, and this is actually in accordance with Norton’s theory, where the large values for the so-called numerical distance (ND) results in very small values for the attenuation coefficient, hence the small surface wave figures [39].

Likewise, notice in the bottom graph of Figure 4 how quickly 𝚐${}_{\mathrm{ex}}\left(\xi \right)$ vanishes (its real part is shown), making the convergence of (19) very fast. Moreover, due to the alternating positive and negative values, it is expected that the effect of 𝚐${}_{\mathrm{ex}}\left(\xi \right)$ is insignificant. The same holds true for the $\mathsf{\rho }$-component of (19), justifying the small observed values, regarding the surface wave field. Put it differently, for the case shown in Figure 4, the major contribution in (19) comes from a narrow area around the stationary point, which in this problem lies within the $\left[-\mathsf{\pi }/2,+\mathsf{\pi }/2\right]$ range [29]. This contribution yields the reflected field in an asymptotic sense, as shown in [23] with the application of the SPM method. In the rest of the integration range, the integrand is related to the surface wave [45] and exposes a behavior similar to Figure 4, thus having minimum impact to the final result. Therefore, the results are also a validation of the SPM method, which as mentioned in Section 1, it emerges as the high-frequency asymptotic solution for the complete problem.

As a last validation, in Figure 5, we demonstrate various field components, for the exact scenario illustrated in Figure 4 of [23]. In [23], only the Norton’s surface wave was evaluated, whereas here, we also compare with the NI results. Moreover, we perform a comparison between the analytic expression for the LOS field and its equivalent integral form (“LOS field NI” in Figure 5), as both are given in Section 2 by (5) and (18), respectively. Again, our numerical evaluation for the surface wave is more or less identical to Norton’s values. Needless to say, we also achieve a perfect match between (5) and (18), essentially meaning that our redefined integral formulation for the EM field, described in Section 2.2, is effective and accurate.

We close this section with a few comments regarding the method’s efficiency. The convergence time depends on four aspects: (a) the utilized HW and SW platform; (b) the selected NI algorithm; (c) the required error tolerance; and (d) the problem parameters, especially the frequency of operation, for a given Transmitter–Receiver (T–R) positioning or on the electric distance, $k_{01}r$, when their combined effect is accommodated. The first three factors seem quite reasonable. Regarding the fourth one, Figure 6 provides a good reasoning. It displays the behavior of the first integrand of (19) (in this case, the real part of the $\mathsf{\rho }$-component), nominated 𝚑${}_{\mathrm{e}\mathsf{\rho }}\left(\xi \right)$ in the figure. It is obvious that higher frequencies contribute to additional oscillations and this is not surprising if one observes that 𝚑${}_{\mathrm{e}\mathsf{\rho }}\left(\xi \right)$ includes a Bessel and a phase function, which are increasingly fluctuating for larger arguments or when the frequency $f=\frac{k_{01}·c}{2\pi }$ increases. Hence, more steps or intervals are required, for the NI algorithm to achieve a target error.

Table 1. Convergence time (ms).

f	Adaptive Simpson’s			Trapezoidal
(MHz)	${\mathbf{10}}^{-\mathbf{3}}$	${\mathbf{10}}^{-\mathbf{6}}$	${\mathbf{10}}^{-\mathbf{9}}$	${\mathbf{10}}^{-\mathbf{3}}$	${\mathbf{10}}^{-\mathbf{6}}$	${\mathbf{10}}^{-\mathbf{9}}$
1	3.88	6.02	18.73	3.43	43.14	969.11
3	4.48	6.98	20.84	4.62	50.94	1272.52
10	5.11	7.92	21.69	7.43	69.27	1921.38
30	6.82	9.02	25.16	14.38	114.05	2923.07
80	9.80	14.93	32.15	31.25	240.26	6748.43
100	11.00	16.20	39.55	39.46	354.70	9560.96
300	21.29	35.77	61.24	57.60	520.40	15,437.25
1000	58.44	103.55	156.89	126.47	973.65	32,240.70

Convergence times for the NI of (19) at various frequencies (f) and different relative error tolerance levels. The horizontal distance was set to ρ = 1 km. The rest of the problem parameters, X, X₀, σ, ε_r, I, and 2 h, were set as in Figure 5. Simulations performed on a 64 bit, Quad Core [email protected] GHz, 16.0 GB RAM platform, using MatLab.

demonstrates the performance of our method at various frequencies, utilizing two common NI techniques to evaluate (19): the Adaptive Simpson’s and the Trapezoidal method [46]. We are able to evaluate the fields at almost an arbitrary accuracy level and at very reasonable convergence times and this is achieved without any need to apply specialized techniques, so as to accurately calculate the Sommerfeld integral tails [40,47,48,49,50,51,52] (note that ${10}^{-12}$ or even lower error is also easily achievable at the expense of computational time). Table 1 also exposes the effectiveness of adaptive quadrature for evaluating such ill-behaved, rapidly fluctuating functions, such as 𝚐${}_{\mathrm{ex}}\left(\xi \right)$ and 𝚑${}_{\mathrm{e}\mathsf{\rho }}\left(\xi \right)$ of Figure 3, Figure 4, and Figure 6 [53,54]. Depending on the allowed error, it seems that, at high frequencies, the selection of an adaptive quadrature technique might be necessary for getting timely results. Similar results are obtained if we apply one of the quadrature integration formulae of Gaussian type, found in [55]. Particularly, to evaluate the two integrals terms of (19), expressions 25.4.30 and 25.4.45 of [55] apply, with total convergence times being equivalent to the adaptive Simpson’s values, shown in Table 1.

3. Evaluating a Novel Asymptotic Solution to the Sommerfeld Problem

Now that we have a solid tool for the calculation of Sommerfeld integrals, we use it to examine a novel asymptotic solution to the Sommerfeld radiation problem. The solution is based on a modified steepest descent method that considers the location of the pole in relation to the saddle point. As such, it may be related with alternative approaches, such as [56]. However, it has distinct characteristics, which, for the cases examined, leads to useful asymptotics and reveals interesting facts for the field’s behavior.

3.1. Outline of the Asymptotic Method

For the usual case, where $\mathsf{\sigma }\gg \mathsf{\omega }{\mathsf{\epsilon }}_0$, the field scattered by the plannar interface can be expressed as [30]

$$\begin{array}{c}\hfill {\underline{E}}^R=-{\widehat{e}}_{{\mathsf{\theta }}_2}{\displaystyle \frac{p{k}_{01}^3}{2{\mathsf{\epsilon }}_0{\mathsf{\epsilon }}_1}}\sqrt{{\displaystyle \frac{-2i}{\mathsf{\pi }{k}_{01}\mathsf{\rho }}}}·e^{i{k}_{01}{r}_{2}cos{\mathsf{\zeta }}_p}·{sin}^{{}^{\frac{3}{2}}}{\mathsf{\theta }}_{2}sin\frac{{\mathsf{\zeta }}_p}{2}{R}_{\Vert }\left({\mathsf{\theta }}_2\right)X\left(k_{01}{r}_2,-{\mathsf{\zeta }}_p\right),\end{array}$$

(20)

where, with respect to Figure 1, ${\widehat{e}}_{{\mathsf{\theta }}_2}={\widehat{e}}_{\mathsf{\rho }}cos{\mathsf{\theta }}_2-{\widehat{e}}_{x}sin{\mathsf{\theta }}_2$ refers to the unit vector along the ${\mathsf{\theta }}_2$–direction of a spherical coordinate system, in which the origin is the dipole’s image (${\mathrm{A}}^{\prime }$) and $R_{\Vert }\left({\mathsf{\theta }}_2\right)$ is given by (10), for $\xi ={\mathsf{\theta }}_2$. Moreover, in (20), ${\mathsf{\zeta }}_p={\xi }_p-{\mathsf{\theta }}_2$, where ${\xi }_p$ is the pole of $R_{\Vert }\left({\mathsf{\theta }}_2\right)$, which in turn may be approximated by

$$\begin{array}{cc}\hfill {\xi }_{\mathrm{p}}& \simeq \frac{\mathsf{\pi }}{2}+\sqrt{\frac{\mathsf{\omega }{\mathsf{\epsilon }}_0{\mathsf{\epsilon }}_1}{2\mathsf{\sigma }}}·\left\{1+\frac{\mathsf{\omega }{\mathsf{\epsilon }}_0\left({\mathsf{\epsilon }}_1+{\mathsf{\epsilon }}_2\right)}{2\mathsf{\sigma }}-i\left[1-\frac{\mathsf{\omega }{\mathsf{\epsilon }}_0\left({\mathsf{\epsilon }}_1+{\mathsf{\epsilon }}_2\right)}{2\mathsf{\sigma }}\right]\right\}.\hfill \end{array}$$

(21)

The most interesting part in (20) is special function X, the so-called “Etalon integral” (EI) [31,32,33,34,35,36]. For parameters $k,\mathsf{\alpha }$, it is defined as the contour integral

$$\begin{array}{ccc}\hfill X\left(k,\mathsf{\alpha }\right)& =\hfill & \frac{1}{4\mathsf{\pi }i}{\int }_S\frac{e^{ik\left(cos\mathsf{\zeta }-cos\mathsf{\alpha }\right)}}{sin\frac{\mathsf{\zeta }+\mathsf{\alpha }}{2}}d\mathsf{\zeta }=\frac{e^{-i\frac{\mathsf{\pi }}{4}}}{\sqrt{2\mathsf{\pi }}}{\int }_{\infty sin\frac{\mathsf{\alpha }}{2}}^{2\sqrt{k}sin\frac{\mathsf{\alpha }}{2}}{e}^{\frac{i{t}^2}{2}}dt\\ & =\hfill & -\frac{1}{2}\mathrm{sgn}\left(\mathrm{Re}\left\{\mathsf{\alpha }\right\}\right)+\frac{1}{2}\mathrm{erf}\left(\sqrt{-2ik}sin\frac{\mathsf{\alpha }}{2}\right),\end{array}$$

(22)

along path S of Figure 7. Regarding the notation in (22): $\infty sin\frac{\mathsf{\alpha }}{2}=\left\{\begin{array}{c}+\infty ,sin\frac{\mathsf{\alpha }}{2}>0\hfill \\ -\infty ,sin\frac{\mathsf{\alpha }}{2}<0\hfill \end{array}\right.$. The “Etalon integral” has useful properties and can be expressed in terms of Fresnel integrals that enable its easy evaluation using the error function. Keep in mind that, to reach (20), we used SDP and then, through the residue theory [57], we deformed the original Sommerfeld contour of integration, $S_z$, into S, so as to apply the above expression for the Etalon integral.

It is possible to elaborate on (22) if we apply the large and small argument approximations for $\mathrm{erf}\left(z\right)$ [37]. This gives

$$\begin{array}{cccc}\hfill X\left(k,\mathsf{\alpha }\right)& \simeq \hfill & \hfill -\sqrt{\frac{i}{2\mathsf{\pi }}}\frac{e^{\left[ik\left(1-cos\mathsf{\alpha }\right)\right]}}{2\sqrt{k}sin\frac{\mathsf{\alpha }}{2}},& \sqrt{2k}|sin\frac{\mathsf{\alpha }}{2}|\gg 1,\end{array}$$

(23)

$$\begin{array}{cccc}\hfill X\left(k,\mathsf{\alpha }\right)& \simeq \hfill & \hfill -\frac{1}{2}\mathrm{sgn}\left(\mathsf{\alpha }\right)+\sqrt{\frac{k}{2\mathsf{\pi }i}}\mathsf{\alpha },& \frac{k{\mathsf{\alpha }}^2}{2}<1,\end{array}$$

(24)

which when used in (20), i.e., for $k=k_{01}{r}_2$ and $\mathsf{\alpha }=-{\mathsf{\zeta }}_p$, yields the following analytic expressions:

$$\begin{array}{ccc}\underline{}{E}^R\hfill & \simeq \hfill & -{\widehat{e}}_{{\mathsf{\theta }}_2}{R}_{\Vert }\left({\mathsf{\theta }}_2\right)\frac{p{k}_{01}^2}{4\mathsf{\pi }{\mathsf{\epsilon }}_0{\mathsf{\epsilon }}_1{r}_2}sin{\mathsf{\theta }}_2·e^{i{k}_{01}{r}_2},\sqrt{2k_{01}{r}_2}·sin\frac{\mathsf{\phi }}{2}\gg 1,\end{array}$$

(25)

$$\begin{array}{c}\hfill \\ {\underline{E}}^R\hfill & \simeq \hfill & {\widehat{e}}_x\mathsf{\delta }\frac{p{k}_{01}^3}{4{\mathsf{\epsilon }}_0{\mathsf{\epsilon }}_1}·\frac{1}{\sqrt{\mathsf{\pi }{k}_{01}\mathsf{\rho }}}{e}^{-\mathsf{\delta }{k}_{01}\left(x+x_0\right)}·e^{i\left(k_{01}\mathsf{\rho }+\mathsf{\pi }/2\right)}·\left[\phantom{\sqrt{\frac{k_{01}\mathsf{\rho }}{\pi }}}\right.\\ & & \left.1+2i\sqrt{\frac{k_{01}\mathsf{\rho }}{\pi }}·\mathsf{\delta }\left(1+k_{01}\mathsf{\rho }{\mathsf{\delta }}^2\right)\right],\begin{array}{c}\hfill \mathsf{\phi }\to 0\\ \hfill k_{01}\mathsf{\rho }{\mathsf{\delta }}^2<1\end{array},\end{array}$$

(26)

with $\mathsf{\delta }=\sqrt{\frac{\mathsf{\omega }{\mathsf{\epsilon }}_0{\mathsf{\epsilon }}_1}{2\mathsf{\sigma }}}\ll 1$ and $\mathsf{\phi }$ being the grazing angle of Figure 1. One may identify that the second condition of (26) is a type of numerical distance [6]; more on it is presented in Section 3.2.

Expression (25) indicates the reflected field emanating from the dipole’s image point (${\mathrm{A}}^{\prime }$ in Figure 1). It should be accurate for a long electric distance, $k_{01}·r_2$, i.e., at the far field region, provided that, at the same time, the grazing angle $\mathsf{\phi }=\mathsf{\pi }/2-{\mathsf{\theta }}_2$ is not very small. In [23], we also reached (25), using the SPM method. However, this method required only the fulfillment of a large electric distance [29]. The effect of the grazing angle was overlooked, and hence, the propagation mechanism for the case of a low height transmission link could not be highlighted. Pay attention to the fact that, if (25) was accurate, even for sliding angles of incidence, just because of a high-frequency transmitting source, the field to be received would essentially be imperceptible, since, in this case, the reflection coefficient, $R_{\Vert }$, approaches $-1$ and ${\underline{E}}^R$ would simply cancel ${\underline{E}}^{LOS}$.

Regarding (26), we are given an expression that describes the behavior of the scattered field for sliding angles of incidence. Due to the existence of the exponentially decaying function, $e^{-\mathsf{\delta }{k}_{01}\left(x+x_0\right)}$, it is confined near the interface as if it is a kind of a surface wave. The field is also spatially limited by the the required conditions, $\mathsf{\phi }\simeq 0,k_{01}\mathsf{\rho }{\mathsf{\delta }}^2<1$, with the implications described in Section 3.2 below. However, pay attention to the fact that this is not a true surface wave, at least when one of the accepted definitions for a type of surface wave is considered [45] (the pole ${\xi }_p$ is located outside the closed contour of Figure 7 and hence its residue is not considered). It simply resembles a surface wave, and this is a consequence of the boundary conditions and the pole’s proximity to the saddle point and the contour of integration (see Figure 7). We choose to call this a pseudo-surface wave.

In the simulations that follow, we compare the introduced closed-form asymptotic solution, (20) against the SPM-based solution of [23]. The reference for our comparisons are the numerical integration results that we obtain using the methodology of Section 2, that is, the numerical integration of (18), (19) and the respective formulas for the magnetic field. Additionally, we examine and comment on the accuracy of analytic expressions (25) and (26).

3.2. Simulation Results

We exhibit two sets of simulations, in Figure 8 and Figure 9, below. Figure 8 demonstrates the effect of the frequency on the received electric field, for a number of scenarios, regarding the T–R horizontal distance, denoted with “d” in the plots. With the exception of Figure 8f, the basic simulation parameters are shown in

Table 2. Simulations parameters.

Symbol	Description	Value
$f_{min}$	minimum frequency	1 MHz
$f_{max}$	maximum frequency	1 GHz
$X_0$	height of transmitting dipole	60 m
$X$	height of receiver’s position	15 m
$I$	dipole’s ipole current source	1 A
$2h$	dipole’s length	0.1 m ${}^{\mathrm{a}}$
$\mathsf{\sigma }$	ground conductivity	4.8 S/m ${}^{\mathrm{b}}$
${\mathsf{\epsilon }}_r$	ground relative permitivity	80 ${}^{\mathrm{b}}$
$\mu$	ground permeability	$4\pi \times {10}^{-7}$ H/m
	numerical integration method	Adapt. Simpsons
	relative error tolerance	${10}^{-6}$

^a much smaller than the wavelength λ = c/f, ^b pertains to the case of sea water.

. The ground parameters, ${\mathsf{\epsilon }}_r$, $\mathsf{\mu }$, and $\mathsf{\sigma }$, are indicative for the case of sea water and do fulfill the requirement, $\mathsf{\sigma }\gg \mathsf{\omega }{\mathsf{\epsilon }}_0$. The altitudes $X_0$ and X are set at 60 m and 15 m, respectively; however by increasing the horizontal distance, d (up to 30 km in Figure 8e), we simulate sliding angles of incidence as well. Finally, for the evaluation of the error function in (22), which according to (20) encompasses a complex argument ($-{\mathsf{\zeta }}_p$), the algorithms of [58,59] were appropriately modified and utilized, since they very accurately evaluate such special functions in the complex plane.

The case of Figure 8a is indicative of non-near-ground-level reception. The T–R relative position is such that the angle of incidence is $\mathsf{\phi }\simeq {15}^{\circ }$. It is evident that there is an almost perfect match between the results obtained numerically, labeled as “NI” in the plots, and what is predicted by (20), depicted via the “Etalon” lines in Figure 8 and Figure 9. It is equally interesting that SPM also yields similar results, particularly for frequencies around 20 MHz and above. Keep in mind that the SPM solution is essentially the expression given by (25), which in turn is derived as a special case of (20). However, the requirement for (25), $\sqrt{k_{01}{r}_2}·sin\frac{\mathsf{\phi }}{2}\gg 1$, is not strictly fulfilled in our case; it goes from 0.31 at 1 MHz to 3.1 at 100 MHz. At 20 MHz, it is about 1.4. Thus, it seems that (25) is an accurate analytic expression for a non-sliding angle of incidence reception in which the validity could be practically extended beyond the strict restrictions imposed for its derivation.

In Figure 8b, the T–R distance is increased to 3 km and as a result the angle of incidence is radically reduced to $\mathsf{\phi }\simeq 1.{43}^{\circ }$. Here, we do observe a discrepancy between the two asymptotic solutions and of both of them with the reference NI results. Of course, this discrepancy appears to be relatively small, and if examined in a broader frequency range, as in Figure 8c, it may be practically regarded as negligible. Nevertheless, it is important to note the tendency of (20) to better follow (19), something that is even more apparent in diagrams (d) and (e) of Figure 8. In these cases, the T–R distance is further increased to 10 km and 30 km, with the incidence angles now being as sliding as $\mathsf{\phi }\simeq 0.{43}^{\circ }$ and $\mathsf{\phi }\simeq 0.{14}^{\circ }$, respectively. Overall, compared with the solution of [23] (SPM-based solution), the Etalon-based solution is a better estimate for the received field. It is also apparent that both methods smoothly converge to (19) in the high frequency regime but that (20) converges faster.

On the contrary, Figure 8f verifies a somehow expected behavior. In lower frequencies, both methods fail to describe the propagation mechanism, being unable to capture the effect of the surface wave, which in this scenario should be rather significant (consider also the lower T–R heights, selected particularly in this case, such that the presence of the surface wave is further exaggerated [39]). Indeed, (20) behaves only marginally better than the respective asymptotic formula of [23]. As already stated, the latter essentially yields only the space wave component. In conclusion, since both solutions are based on the application of high-frequency assymptotic methods (Saddle Point vs SPM), they may therefore yield accurate results only in the high-frequency regime [60].

To confirm and further solidify the above arguments, in Figure 9, we exhibit the field behavior from the perspective of a varying T–R distance. Starting from lower frequencies, in Figure 9a, we once again observe that both asymptotic methods fail to produce accurate results (for reasons of completeness, magnetic field values are given in this case). According to our simulations, this situation holds true almost up to 1 MHz. Moving towards the HF frequency zone, the advantages of the newly introduced asymptotic solution show up in Figure 9b. The difference between (20) and the SPM solution of [23] is more evident for large distances, where the effect of the scattered field is more significant; hence, the improvement that the “Etalon” function, X, yields in (20) becomes visible. Proceeding to the VHF zone of Figure 9c, we realize that both methods begin to converge and ultimately coincide with the complete solution at even larger frequencies, as shown in Figure 9d. At those frequencies, the surface wave is almost negligible. Therefore, there is nothing extra left for special function X to expose and (20) simply yields the reflected field, exactly as the solution of [23] does.

The last two diagrams of Figure 9 are devoted to the investigation of (26), an interesting expression, which, as mentioned in Section 3.1, attributes surface wave characteristics to the near-ground-level scattered field. For this, extended simulations were run, and the outcomes are summarized in Figure 9e,f. Overall, we observe that (26) does converge to (20), from which it was derived when $\mathsf{\phi }\to 0$. In Figure 9e, the convergence occurs at about 6 km, which, for the selected T–R heights, is equivalent to a grazing angle $\mathsf{\phi }$ of less than $1^{\circ }$. As mentioned before, in such frequencies (≤3 MHz), (20) and hence (26) are not very accurate approximations of the complete solution. They are simply better estimates compared to the SPM and hence the visible gap between them and the NI results in Figure 9e.

The same scenario is repeated at 30 MHz and illustrated in the top plot of Figure 9f. This time, convergence is achieved before 1400 m, where $\mathsf{\phi }$ is just above $3^{\circ }$. In general, for the same T–R heights, we observe a tendency of slightly increasing grazing angles to achieve convergence, as the frequency rises. However, bear in mind that the numerical distance, $k_{01}\mathsf{\rho }{\mathsf{\delta }}^2$, quickly approaches the value of 1 for its dependence on the square of the frequency. This eventually shrinks the range for which (26) is accurate. A similar behavior is observed in the bottom diagram of Figure 9f at 80 MHz. Due to the lower antenna heights, we have a match at about 650 m, where $\mathsf{\phi }\simeq 2^{\circ }$. Nevertheless, because of the higher frequency, the validity of (26) is now limited to an even shorter distance range. We thus reach to the conclusion that this pseudo-surface wave has local significance and that its existence is highly sensitive on the value of the numerical distance, $k_{01}\mathsf{\rho }{\mathsf{\delta }}^2$. When the numerical distance exceeds 1, it essentially disappears.

4. Conclusions and Future Research

We demonstrated an efficient method for the numerical calculation of Sommerfeld integrals. The method proves fast and accurate, and when applied to the evaluation of the EM field, radiated by a vertical dipole above flat, lossy ground, it fits very well with existing solutions and Norton’s results.

One limitation of the method is that its efficiency is theoritically limited to non-diminishing antennae heights. This is an outcome of the exponentially decaying factor $e^{-k_{01}(x+x_0)sinh\xi }$ in the integrand of (19), which is “neutralized” when $x+x_0=0$. However, the simulation results showed that, for most practical frequency cases, it takes only a few centimeters of antenna altitudes for the numerical calculation to quickly converge. In fact, the effect of the frequency seems rather more important than the heights of the antennas. In any case, the strict limits of the method’s applicability (probably for rather unexpected practical cases, as an antenna is supposed to be mounted at some distance above the ground) are left for future study.

Using the above numerical method as a reference, we then focused on the evaluation of a high-frequency asymptotic solution. The solution relies on the saddle point method and uses the properties of the so-called “Etalon integral”, as a means to increase the accuracy of the results. We verified that, for the usual case, where $\mathsf{\sigma }\gg \mathsf{\omega }{\mathsf{\epsilon }}_0$, the method does succeed to provide better estimates, as compared with a more basic asymptotic approach, based on the application of the stationary phase method. Moreover, further asymptotic properties for the “Etalon integral” allowed us to reach analytic formulas for the scattered field. Of particular interest is (26) that exposes surface wave characteristics to the field near the interface, provided that the numerical distance $k_{01}\mathsf{\rho }{\mathsf{\delta }}^2<1$.

From the analysis of Section 3.1, one might identify that the so-called lateral waves were not taken into cosideration. However, we still managed to obtain a quite good agreement with NI results. The ultimate goal is to provide asymptotics, applicable for every possible scenario, not only for the usual $\mathsf{\sigma }\gg \mathsf{\omega }{\mathsf{\epsilon }}_0$ case but also for circumstances where the pole is located inside the integration path of Figure 7, where it is possible to analytically extract the surface wave component by means of the residue theory. For that purpose, in the future, we will insist on the investigation of special function $X\left(k,\alpha \right)$ and its properties as well as other special functions that could be used to describe the field behavior. Additionally, the obtained asymptotic solutions can be refined with any accuracy and presented in expanded form, according to known procedures [60] and a comparison with relevant recent research, e.g., [61,62,63,64], is to be made.