The Sommerfeld half-space problem revisited: from radio frequencies and Zenneck waves to visible light and Fano modes

Abstract

The classical Sommerfeld half-space problem is revisited, with generalizations to multilayer and plasmonic media and focus on the surface field computation. A new ab initio solution is presented for an arbitrarily oriented Hertzian dipole radiating in the presence of a material half-space with arbitrary horizontal stratification. The solution method combines the vector potential approach and the spectral domain transmission line analog of the medium, which results in the most compact formulation and facilitates the inclusion of any number of layers in the analysis. Following Sommerfeld, the solution is first expressed in terms of the Fourier–Bessel transforms, also known as Sommerfeld integrals. Analytical properties of the integrands in the complex plane are then investigated, including the location of the Sommerfeld pole, which gives rise to the Zenneck wave (ZW) or surface plasmon polariton (SPP), and alternative field representations are developed by a deformation of the integration path and analytic continuation of the integrand functions, using hyperbolic and vertical branch cuts. Closed-form expressions for the asymptotic surface fields are also derived and the rôle of the ZW and SPP is elucidated. Numerical examples are included to illustrate the theory, from radio frequencies to visible light.

Keywords:

• Hertzian dipole
• Green function
• layered half-space
• Sommerfeld integral
• Zenneck wave
• surface plasmon polariton

View correction statement:

Correction to: The Sommerfeld half-space problem revisited: from radio frequencies and Zenneck waves to visible light and Fano modes

1. Introduction

The problem of a vertical Hertzian dipole (time-harmonic current element of infinitesimal length) radiating in the presence of a lossy conducting half-space was first rigorously solved by Sommerfeld [1] in 1909 and other dipole configurations were treated in his later works.[2–4]¹ Sommerfeld originally placed the dipole in the interface between the two half-space media, but this caused some unnecessary mathematical difficulties,[3,6] hence most authors – including Sommerfeld in his later monograph [4] – put it above the half-space, as indicated in Figure 1.

In his celebrated 1909 paper, Sommerfeld first expressed the solution in terms of improper integrals with Bessel function kernels – Fourier–Bessel transforms (FBTs) now known as Sommerfeld integrals, which cannot be evaluated in closed form. Sommerfeld then deformed the real-axis integration path into the complex plane, capturing two hyperbolic branch cuts, whose contribution he identified as space waves (Raumwellen), and a pole, whose residue he identified as a radial–cylindrical form of the plane surface wave of Zenneck [7]. Sommerfeld further derived an asymptotic surface field expression, assuming high media contrast and large horizontal range. This formula consisted of a whole-space term multiplied by a factor – now known as the Sommerfeld attenuation function, which accounts for the influence of the ground through a complex parameter – the Sommerfeld numerical distance. Unfortunately, as it later turned out, this groundwave (Bodenwelle) expression was erroneous.² Subsequently, as documented by Maclean and Wu [10], several authors – among them Weyl [11] and Sommerfeld in a paper under the same title as the original one [2] – solved the problem by different methods and obtained the correct asymptotic surface field formula. A reprise of Sommerfeld’s 1926 solution was later given by Kockel [12]. Inexplicably, the discovery of the sign error in Sommerfeld’s original paper unleashed a prolonged controversy in the literature as to its root cause – even though the correct solution was available since at least 1926.[13,14] Various explanations were put forward – all eventually debunked, and doubts were even cast on the very existence of the Sommerfeld pole and the associated ZW.³ From the modern perspective, it is rather surprising that such a relatively simple boundary value problem, which is capable of being solved completely with no room for debate, should have generated so much controversy – with some echoes of the past confusion reverberating to this day.[20,21] The groundwave formulation that eventually became the de facto standard is due to Norton [22], who based it on an earlier work of van der Pol and Niessen [23,24]. This Norton surface wave,[13] or Sommerfeld–Norton groundwave,[25] [10, p.97], [26, p.310], which is applicable at any point in space not too near the antenna, was later derived by Wait [27] using the modified saddle-point method, and by King [28], and more recently by Green [29], using the surface impedance concept and compensation theorem. Perhaps the most recent groundwave formulation is due to King [30], who – eschewing the use of Hertz potentials and complex integral transforms – based it on the real-axis field representation.

Figure 1. Hertzian dipole with the current moment above a half-space medium.

Sommerfeld’s work was motivated by the desire to explain the radio wave propagation along the surface of the earth, before it was realized that reflections from the ionosphere were in fact the decisive factor. However, in recent years there has been a resurgence of interest in the Sommerfeld half-space problem and the Sommerfeld groundwave in particular, in the context of THz applications, near-field optics, plasmonics, and nanophotonics.[31–34] The ZW becomes a SPP in the case of noble metals at optical wavelengths.[35–37] An important extension of the Sommerfeld problem arises when the lower half-space comprises two or more layers. This case was first treated by Grosskopf using the Weyl formulation,[38] and later by Wait [39–42], and more recently by King and Sandler [43].⁴ At radio frequencies, the two-layer case can model the effect of a thin ice layer on top of the ground or body of water,[43] while in the THz regime and plasmonics, it applies in case of a dielectric-coated metal [45] or silver contaminated by a nanolayer of oxide tarnish.[46] The multilayer geometry is relevant, for example, to the ground penetrating radar (GPR) modeling of water spillage from underground pipes.[47]

The solution procedure employed by Sommerfeld may be classified as the method of scattering superposition,[48] since it uses a Fourier–Bessel integral representation of the primary excitation – the Sommerfeld identity, as the point of departure. The reflected and transmitted fields are then expressed in compatible axial transmission representation [49] forms and the total field integrands are matched at the interface, using the boundary conditions to determine the unknown amplitudes. Subsequently, the real-axis integrals may be analytically continued into the complex plane, taking into account any singularities captured in the path deformation process – in this case the branch points associated with the upper and lower half-spaces and the Sommerfeld pole. If the path is confined to the proper Riemann sheet, as in Sommerfeld’s original work, the radial transmission representation [49] of the solution is obtained, comprising integrals around the hyperbolic branch cuts and the pole residue. However, there is a considerable leeway in the choice of branch cuts, and vertical cuts are often employed in surface wave studies, since the path around the branch point associated with the upper half-space then becomes the steepest descent path. The Sommerfeld pole contribution may, or may not, explicitly appear in this nonspectral representation,⁵ depending on the media parameters, but it always exists in the radial transmission representation (with the hyperbolic branch cuts), if there is some loss present in the lower half-space medium.

In the present paper, we revisit the Sommerfeld half-space problem and present a new, ab initio solution – dispensing with the method of scattering superposition – for the more general case of a lower half-space with arbitrary horizontal stratification. Our method combines the vector potential approach and the spectral domain transmission line (TL) network analog of the medium, which results in the most compact formulation and facilitates the inclusion of any number of layers in the analysis. The surface fields are first expressed as FBTs (or Sommerfeld integrals) of the appropriate combinations of the TL Green functions, resulting in the axial transmission representation. The analytical properties of these integrands in the complex plane are then investigated – including the location of the Sommerfeld pole, which gives rise to ZW or SPP, and alternative field representations – termed radial transmission and nonspectral – are developed by suitable deformations of the real-axis integration path around the singularities in the complex plane. Using the nonspectral representation as the point of departure, closed-form expressions of the asymptotic surface fields are derived and the rôle of the Zenneck and surface plasmon waves is discussed. Although from the mathematical standpoint the solution of this problem is identical whether at radio wavelengths or in the optical regime, there are some important nuances in the way these waves contribute to the surface field, and these nuances are elucidated. Further, certain misleading statements recently promulgated in the literature regarding the existence of the Sommerfeld pole and the associated ZW are pointed out. The theoretical developments are supported by numerous numerical examples – from the lossy ground at radio frequencies to plasmonic media in the visible range, including homogeneous and coated half-spaces excited by vertical and horizontal dipoles.

2. Formulation

2.1. Derivation of the Green dyadics – the general case

We begin with the Maxwell’s equations [50]⁶(1)

for a medium with -dependent permeability and permittivity . Using the fact that is divergenceless, we introduce the vector potential and the scalar potential as(2)

to obtain(3)

where is the wavenumber of the medium. Since the number of unknowns in (3(3) ) exceeds by one the number of equations, an additional normalizing equation is needed. Upon inspecting the Cartesian components of (3(3) ), we note that the gage condition [51](4)

results in great simplifications. Using this gauge in (2(2) ), we find that the electric and magnetic fields may be computed from alone as(5)

Next, we use (4(4) ) to eliminate from (3(3) ) and obtain(6) (7) (8)

where the subscript denotes vectors and operators transverse to . It is apparent from these equations that unless the medium is homogeneous, a horizontal - or -directed current source requires not only the corresponding horizontal component of , but also the vertical component, . Further, these equations suggest that we introduce the Green functions , , and as⁷(9) (10) (11)

where is the Dirac delta function. We thus infer that for an arbitrary current distribution, occupying a volume , the vector potential may be found as(12)

where(13)

with . We have thus obtained the simplest possible form of the vector potential Green dyadic , expressed in terms of three scalar functions, for the case of -dependent media.[51–53] Note that the off-diagonal terms in (13(13) ) are absent only for a homogeneous space. For later convenience, we express (13(13) ) as⁸(14)

Turning next to the task of solving (9(9) )–(11(11) ), we first note that the equation for involves in the forcing function, which is cumbersome. To remedy this, we introduce an auxiliary Green function by means of the relation(15)

which upon using (11(11) ) and (9(9) ) leads to(16)

We will thus find indirectly from (15(15) ) by first solving (9(9) ) for and (16(16) ) for . We will also initially assume that the source point lies on the -axis of the cylindrical coordinate system , which results in a much simpler -independent configuration, and subsequently generalize to the off-axis dipole case by the substitutions and , where(17)

Hence, under this assumption, we may express (9(9) ) as(18)

where now(19)

Further, in view of the azimuthal symmetry, we may avail ourselves of the FBT pair⁹ [55, p.21](20)

to simplify the solution. Hence, the application of the direct FBT has the effect(21)

and reduces (18(18) ) to the Sturm–Liouville (SL) problem [56, p.274](22)

subject to the requirement that [57, p.176](23)

Similar treatment of (10(10) ) yields another SL problem(24)

subject to the requirement that(25)

Finally, applying this procedure to (16(16) ), we obtain(26)

subject to (25(25) ), and the direct FBT of (15(15) ) yields(27)

While (22(22) )–(26(26) ) may be solved by a conventional Green function procedure,[57] herein we utilize the TL network method.[49,56,58] The relevant network problems are illustrated in Figure 2(a) and (b), which depict TLs excited by ideal unit-strength current and voltage generators, respectively.

Figure 2. Network problem for the determination of the TL Green functions.

These TLs are characterized by -dependent propagation “constants” and characteristic impedances , where the superscript can assume the values (TM mode) or (TE mode), given as [56, p.190](28)

For the current source excitation, the voltage and current obey the TL equations(29)

where denotes the characteristic admittance, and we find that satisfies the second-order differential equation(30)

For the voltage source excitation, the voltage and current obey the TL equations(31)

and satisfies the second-order differential equation(32)

The voltages and currents in Figure 2 may be considered TL Green functions, since they are excited by unit-strength impulsive sources. The superscript on or indicates the TL type ( for TM line, for TE line), and the subscript or indicates the source type (shunt current generator or series voltage generator, respectively). Note that the TM and TE lines have the same propagation constants, but differ in their characteristic impedances – see (28(28) ). The advantage of the network method is that the solution of any SL problem, such as (22(22) ) or (24(24) ), may be expressed in terms of a suitable TL Green function, by inspection of either (30(30) ) or (32(32) ). Since the behavior of TL voltages and currents is well understood, we thus also gain valuable insight as to the nature of the layered medium Green functions.

Applying the network method to (22(22) ), we note the analogy with (30(30) ) in the TE case, where , and a comparison of these equations immediately yields(33)

Note that the continuity conditions (23(23) ) are automatically satisfied, since the voltage and current are continuous outside the source. Similarly, comparison of (24(24) ) and (32(32) ) in the TM case, where , yields(34)

Turning attention to (26(26) ), we note that it differs in form from the SL problems (22(22) ) and (24(24) ), as it comprises a more complicated driving function, involving a derivative of the Dirac delta. However, in this case, we use the analogy with the equation satisfied by , which follows from (29(29) ):(35)

In the TM case, this reduces to (26(26) ), hence . Finally, from (27(27) ) and (33(33) ) we obtain(36)

We have thus expressed the three spectral domain Green functions required in (13(13) ) in terms of the TL Green functions. We note that a vertical Hertzian dipole contributes a voltage source on a TM TL analog, whereas a horizontal dipole contributes current sources on both the TM and TE lines.

What remains is the simple step of the FBT inversion of the the spectral domain Green functions derived above. For later convenience, we first introduce the shortcut notation(37)

for the Sommerfeld integrals. Hence, on applying the inverse FBT to (33(33) ) and (34(34) ), we arrive at(38) (39)

where we have also invoked (17(17) ). Finally, the FBT inversion of (36(36) ) followed by differentiation with respect to yields¹⁰(40)

whereas differentiation with respect to yields , with in place of .

We have thus formulated the Green dyadic in terms of the TL Green functions of the network analog of the medium along the -axis. Since all TL Green functions can be obtained from and , we have in effect expressed all components in terms of two scalar wave functions. Further, only three Sommerfeld integrals need to be computed to determine in the most general case. Similar representations were previously obtained using different methods by Michalski [53] and Michalski and Mosig [59] – the latter including the uniaxial medium case. The approach presented here combines the direct vector-potential method of Dmitriev [51] and the TL analog of the medium,[56] thus resulting in a very compact, yet general, formulation.

Note that no assumptions have so far been made regarding the variation of and or the source location, hence the formulation is applicable to arbitrary current distributions in media with arbitrary stratification perpendicular to the -axis. In particular, the vector potential due to the Hertzian dipole of Figure 1 may be found as(41)

and the corresponding electric and magnetic fields then readily follow from (5(5) ).

2.2. Dipole above a layered half-space with piecewise-constant parameters

Henceforth, we assume that the dipole resides in air above a layered half-space characterized by piecewise-constant and , which in layer assume the values and , respectively. The total number of layers (including the upper half-space) will be denoted by . The corresponding TL network analog is illustrated in Figure 3, where the th layer is represented by a TL section with propagation constant and characteristic impedance .

Figure 3. Hertzian dipole above a layered half-space and TL network analog.

It will be recalled that in general two TL networks arise, with (TM mode) and (TE mode). The relevant TL Green functions are readily found by a simple network analysis, as described in Appendix 1. Using these results in (38(38) )–(40(40) ), we may express the vector-potential Green functions in the upper half-space as(42) (43) (44)

where and is the reflection coefficient looking down from air into the lower half-space. In the above, we have expressed the direct (whole-space) terms in a closed form, using the Sommerfeld identity [1](45)

The electric and magnetic fields for an arbitrarily oriented dipole above the half-space may now be computed via (41(41) ) and (5(5) ). Henceforth, we assume a unit-strength dipole placed at on the -axis and limit attention to surface fields () in air when h. It will suffice to only consider a vertical (-oriented) and a horizontal (-oriented) dipole. For the vertical dipole, we readily find the nonzero cylindrical components of the surface fields as(46) (47) (48)

For the horizontal dipole, we similarly obtain(49) (50) (51) (52) (53) (54)

Note that some of the expressions above are approximate, meaning that we have retained only the dominant terms, neglecting terms that decay more rapidly with . This is permissible, as we focus on the surface fields far from the dipole.

2.3. Integrand singularities and the Sommerfeld pole

The integrands in (46(46) )–(54(54) ) comprise the propagation constants , which are double valued with branch points at , as discussed in Appendix 2. In the case of multilayer media, the integrands are even functions of associated with all finite-thickness layers, hence and are the only branch points. To show this, we note that the change only affects the reflection and transmission coefficients in sections and . If section is finite, we find that under this sign change turns into its inverse, but and all reflection coefficients to the left of section remain unchanged. The integrands in (46(46) )–(54(54) ) may thus be viewed as defined on a four-sheeted Riemann surface, depicted schematically in Figure 4, where the vertical line segments represent the branch cuts and offer a passage between the sheets.[17, p.56], [60, p.613], [61]. The radiation conditions are satisfied everywhere on sheet , which will be referred to as the top or proper sheet.

Figure 4. Schematic diagram of the four-sheeted Riemann surface with fundamental branch cuts associated with and .

In addition to the branch points, the spectral integrands have poles, which may be found from the transverse resonance condition [62, Section 3.6](55)

Hence, the poles are roots of the denominator of the reflection coefficient (A10(A10) ) looking from air into the layered half-space. The dispersion function is also defined over the four-sheeted Riemann surface and can have roots on any of the sheets.[63] However, only those located on the top sheet are proper and are termed surface wave poles; the ones on the other sheets are improper and are termed leaky wave poles. In the general case of a multilayer half-space, the poles can only be determined by a numerical search for the roots of the dispersion function (55(55) ) in the complex -plane, as no closed-form formula is available. The most effective root location methods are based on Cauchy’s theorem – which, however, requires that the dispersion function be analytic.[64–66] Consequently, prior to the application of these methods, the branch points should be removed by a suitable mapping , which unfolds the four-sheeted Riemann surface, [67–69] and the search is performed in the -plane. The required analytical derivative of may be computed by automatic differentiation, [70] which is readily implemented in modern Fortran by operator overloading.[71] Although we employ these techniques in some of the numerical examples presented here, we omit further details for brevity.

If the lower half-space is homogeneous, as in the original Sommerfeld problem, we may set , and the resonance condition reduces to(56)

where in the TM case () and in the TE case (). Upon squaring both sides of (56(56) ), we readily find the solutions , where(57)

However, since the squaring obliterates any distinction of the square root branches, it is still incumbent upon us to prove if, and under what conditions, these poles exist on the top Riemann sheet. Unfortunately, no simple and general criteria can be given for media with arbitrary permittivity and permeability.[61] We thus limit attention here to the important case of nonmagnetic media (), with no restrictions on , other than passivity. For such media, there are no poles in the TE case (, ), while in the TM case (, ) (57(57) ) reduces to(58)

To test if this pole satisfies the resonance condition (56(56) ), we first note that(59) (60)

assuming real , for simplicity. Using (59(59) )–(60(60) ) in (58(58) ), we obtain(61)

which shows that is located in the 4th quadrant if the positive square root is taken. We next evaluate and at , selecting the square root branches with negative imaginary parts, as required on sheet of the Riemann surface of Figure 4. We thus obtain(62) (63)

If we now substitute (62(62) ) and (63(63) ) into (56(56) ) with , we find that the left- and right-hand expressions match, hence the pole (58(58) ) satisfies the resonance condition on the top sheet .¹¹ Further, this pole also exists on sheet , where the signs of and are reversed, whereas on sheets and the root (58(58) ) satisfies the condition (56(56) ) with the minus sign replaced by plus, i.e. on these sheets it is a root of the numerator of the reflection coefficient (A10(A10) ), and is called the Brewster zero [60, p.613].¹²

It is clear from (58(58) ) that in the high-contrast case, where , the Sommerfeld pole may be arbitrarily close to the branch point . To investigate the pole location and the properties of the associated waves in more detail, we let and consider two important special cases below:

(1)

This condition holds for a lossy ground (or seawater) at radio frequencies (with ), as well as for metals in the infrared (including THz) range (with ). Using the binomial expansion approximation in (58(58) ), we readily obtain(64) which shows that the pole is located to the left or right of , depending on whether or , respectively. Assuming that the associated pole wave behaves as , we now find the propagation length and depth of penetration in air as(65) Note that indicates how tightly is this wave bound to the surface.

(2)

This condition holds for dry ground at high frequencies (with ) and for noble metals in the visible range (with ). In this case, using the binomial approximation in (58(58) ), we obtain(66) which indicates that the pole is located to the left of if and to the right of if .13 The propagation distance and depth of penetration in air now become(67)

Table 1. Taxonomy of surface modes over a homogeneous half-space.

(high-loss)		(low-loss)
ZW	SPP	Brewster mode	Fano mode

It will be shown in due course that the pole location relative to the vertical line passing through the branch point affects the form of the surface field representation. We will thus distinguish between the Sommerfeld case, where the pole is to the left of and the associated pole wave is classified as either ZW or Brewster mode, and the plasmonic case, where the pole is to the right of and the associated pole wave is called either SPP or Fano mode.[82] This taxonomy of surface modes is given in Table 1.

2.4. Sommerfeld integration paths – original and extended

The integration path in (46(46) )–(54(54) ), often referred to as the Sommerfeld integration path (SIP), is along the positive real axis on the top Riemann sheet, as illustrated in Figure 5(a). The poles indicated in Figure 5 are only meant to be representative; their actual location strongly depends on the media parameters. If the upper half-space is lossless (air), the branch point lies on the real axis, and the SIP must be indented above it. Although effective methods for the computation of the Sommerfeld integrals along the SIP are available,[83–87] the rapid oscillations of the Bessel functions on the real axis make them inefficient for large . In this case, more efficient integral representations can be obtained by a suitable deformation of the integration path and analytic continuation of integrand functions into the complex -plane. However, in order to make the integration path deformable, we must first extend it to the whole real axis. The desired form [1] is achieved by first expressing the Bessel function in (37(37) ) in terms of the Hankel functions using [54](68)

Figure 5. Integration paths on the top Riemann sheet. A small loss in the upper half-space is assumed for better drawing clarity.

then introducing in the integrand containing the change of variable , and making use of the circuital relation(69)

whence(70)

provided that is an even (odd) function of for even (odd) order . Since the Hankel function introduces an extraneous logarithmic branch point at the origin with a branch cut along the negative real axis, the extended integration path (EIP) runs just under this branch cut, as illustrated in Figure 5(b). Both SIP and EIP are associated with the axial propagation representation of the layered medium Green functions.[49]

Figure 6. Illustration of the path deformation on the top Riemann sheet.

An important attribute of the axial representation formulation is that the location of any poles of the integrand need not be determined. Also, it is the only formulation applicable in the so-called zero-offset case,[88] where , which may be important in some geophysical and nondestructive testing applications. In practice, the SIP has to be deformed around the branch-point and pole singularities near the real axis, which however does not pose much difficulty.[89,90] The only limitation of this representation stems from the oscillatory behavior of the Bessel function on the real-axis path, which makes it impractical for if the horizontal range exceeds a few hundred wavelengths.¹⁴

2.5. Radial transmission representation

After locating the integrand singularities, it is possible to deform the EIP to obtain the radial transmission representation [49,94,95] (or transverse spectral representation [96]) of the Green functions, as we describe presently. Hence, following standard practice,[97,98] we close the integration path in (70(70) ) down by a semi-circular arc at infinity on the top Riemann sheet, avoiding the branch point and pole singularities, as indicated in Figure 6(a). By Cauchy’s theorem, the integral over this closed path is zero. Also, in view of the large-argument form of the Hankel function, viz.(71)

the semi-circle at infinity contributes nothing to the integral, as long as . Consequently, the real-axis path integral in (70(70) ) may be replaced by integrals around the hyperbolic branch cuts and the residue contributions arising from the presence of the poles, as illustrated in Figure 6(b). Note that on the fundamental branch cut associated with , and , and similarly for the branch cut associated with . It is thus convenient to use the positive real variables or as the variables of integration on the paths around the respective branch cuts. Using these substitutions, we may express the integral in (70(70) ) as(72)

where and denote the integrals around the two hyperbolic branch cuts, is the th pole of and is the corresponding residue. If , where and , the residue may be found as(73)

The first branch cut integral is given as(74)

where the square root branch should be selected for , and where we have introduced the notation(75)

for the folded spectral integrand function. Here, the plus (minus) superscript on indicates that every occurrence of is directly replaced by the variable of integration with a plus (minus) sign. For more analytical details, see Appendix 4.

If the upper half-space is lossless, as we assume here, the branch cut associated with follows the real and imaginary axes, as indicated in Figure 6(b), and we may express as [99, p.134], [83](76)

where is given in (75(75) ), where now the plus or minus superscript on indicates that this function should be evaluated on Riemann sheet or , respectively. In the first integral above, we made a trigonometric substitution to remove the weak singularity at the branch point ,[100] and in the second integral, we used the relationship [54, Equation 10.27.8](77)

to express the Hankel function of imaginary argument in terms of the Macdonald function, which is exponentially decreasing.

The second branch cut integral in (72(72) ) is also given by (74(74) ), except that is replaced by and the roles of and are interchanged. In the case of moist earth or seawater at radio frequencies, as well as for metals in the optical range, this integral is negligible in comparison with the first branch cut integral, hence we omit further details for brevity.

It will be instructive to derive the radial transmission representation of the Sommerfeld integral (45(45) ), which represents the vector potential Green function for an infinite homogeneous space. In this case, the integrand function has branch points at and no other singularities. Consequently, , and we readily obtain(78)

where the last integral follows directly from (74(74) ) and (77(77) ) – with replaced by , to make these expressions applicable to the off-axis case. This Schelkunoff integral [101] is the appropriate representation of the primary excitation when the method of scattering superposition is employed in the presence of cylindrical layers concentric with the -axis, but is not applicable in the presence of planar layers perpendicular to the -axis.¹⁵

The Sommerfeld pole, which is located on the top Riemann sheet, is always captured in the radial transmission representation and its residue contributes a pole wave in (72(72) ). Note that, in view of (71(71) ), pole waves comprise a slow, inverse square root decay factor, and may thus potentially dominate the total surface field in the near and intermediate horizontal range, until the exponential decay eventually takes over. As we have shown in Section 2.3, the Sommerfeld pole is located to the left or to the right of the branch point , depending on whether or , respectively, where . In the former case, which arises for a lossy ground (or seawater) at radio frequencies, the corresponding pole wave is called a (radial–cylindrical) ZW. In the latter case, which arises for metals at THz and optical frequencies, this pole wave is called a SPP. Further, for noble metals in the visible range, the loss tangent may be small, in which case the SPP is also referred to as a Fano mode.[82] The ZW and SPP are thus closely related – they satisfy the same resonance condition [102] – and one may evolve into the other, as the media parameters change with frequency.¹⁶ One difference between ZW and SPP is that the former is a fast wave (meaning that its phase velocity exceeds the speed of light¹⁷), while the latter is slow. Although ZW (or SPP) is always a part of the total surface field in the radial transmission representation, nothing can be said a priori regarding its relative importance. In fact, since the pole in the Sommerfeld case is located just below the hyperbolic branch cut, there is a potential for a substantial cancellation of ZW by the branch cut integral contribution. Note that such a compensation is less likely in the plasmonic case, where the pole is located to the right of the branch point, but it may still occur in the high-contrast case.

In view of the fact that ZW is coupled to a free mode of the layered medium, it is not unreasonable to attempt exciting it by a suitable source distribution.[105,106] However, a pure radial–cylindrical ZW can only be excited by a vertical line current of infinite extent, which is not physically realizable.[107], [108, p.10] With a finite dipole, ZW cannot be isolated from the total surface field – and the latter is what can actually be experimentally observed.[109]

Finally, we note that although (76(76) ) is expressed in a computationally convenient form, the closeness of the Sommerfeld pole to the integration path, as well as the rapid oscillations of the first integrand for large , still adversely affect the efficiency of the radial transmission representation. As a result, this representation is not practical for horizontal ranges exceeding a few hundred wavelengths – similar to the axial transmission representation.¹⁸ These limitations are removed in the nonspectral representation, which will be demonstrated next.

2.6. Nonspectral representation

The nonspectral representation of the integral (70(70) ) has the form(79)

where and denote integrals around the vertical branch cuts, as illustrated in Figure 7, and the summation comprises all poles captured in the deformation, possibly including some on the improper Riemann sheets – the so called leaky wave poles [96]. The shaded areas in Figure 7 are the swaths of the improper Riemann sheets swept when the extended real-axis path is deformed to wrap around these branch cuts. Note that if the Sommerfeld pole is located to the left of – as it would be in the lossy earth case, it is not captured and it does not explicitly contribute to the nonspectral representation of the integral. However, due to its close topological proximity to the downward leg of the integration path around , this pole can still strongly influence the result. In the plasmonic case, this pole is located to the right of , and is thus included in the residue summation in (79(79) ), but it may also strongly influence the branch cut integral. Comparing (79(79) ) and (72(72) ) in the Sommerfeld case, we note that , where is the ZW. Hence, in the nonspectral representation, the branch cut integral of the radial transmission representation associated with the branch point and the contribution from the nearby Sommerfeld pole (i.e. the ZW) are treated collectively.

Figure 7. Integration path wrapped around vertical branch cuts and poles captured. The shaded areas are swaths of improper Riemann sheets swept in the path deformation.

The validity of the contour deformation of Figure 7 may be questioned on grounds that part of the infinite semi-circular arc in Figure 6 now falls within the shaded area, which is on the improper Riemann sheet, potentially causing integrand terms like to grow without bound. However, as the radius of the arc approaches infinity within the finite-width vertical strip on the left side of the branch cut, tends toward a positive real number, whence it follows that must be bounded in this region.[110] Consequently, on this part of the semi-circle the integrand of (79(79) ) will have the form of a bounded function of multiplied by the Hankel function, which decays exponentially in the lower half-plane, thus ensuring the vanishing of the integral over the arc at infinity.

We will presently address the computation of the branch cut integrals in (79(79) ). Hence, upon changing the variable of integration via , we express as [111](80)

with(81)

where now represents the jump in as changes sign across the vertical branch cut associated with . In (81(81) ), we have introduced the normalized Hankel function(82)

which approaches unity for . This substitution makes explicit the asymptotic -dependence of the integrand in (80(80) ). The second branch cut integral, , can similarly be expressed, upon the replacement . As already mentioned in connection with the radial transmission representation, the contribution from this integral is typically negligible.

In order to correctly evaluate in (81(81) ), it is required to keep track of the Riemann sheets traversed as the integration in (80(80) ) proceeds along the vertical branch cut. The correct choices of the square root branches of and are readily determined with reference to Figure A1. Hence, in the Sommerfeld case, where in Figure 7, we enforce and in computing (on the left edge of the cut) and reverse the sign of in computing (on the opposite, right edge of the cut). In the plasmonic case, where , the same rules should be followed in the selection of the branch, but should be enforced in and .

Figure 8. -plane in the vicinity of the origin, which corresponds to the branch point in the -plane.

The mapping of the -plane into the -plane introduced in (80(80) )–(81(81) ) is illustrated in Figure 8, where the contours and are images of the corresponding paths in Figure 7 and is the image of the Sommerfeld pole . In the plasmonic case, this pole can cross the real axis and is then captured when the original path is deformed into the real-axis path . As already mentioned, the Sommerfeld pole may be located arbitrarily close to – either on the left or on the right of the vertical branch cut – with the attendant spike in the integrand of (80(80) ), which makes the integration problematic. In such a case, the most reliable approach is to extract the contribution of the offending pole from the integrand and to add it back integrated in closed form [14,46,112,114]. The analysis is similar to the modified saddle-point nethod of Felsen and Marcuvitz [56, Section 4]. Hence, we express the integral (80(80) ) as¹⁹(83)

with(84)

where is the residue of at the Sommerfeld pole . With the pole term subtracted, the modified integral in (83(83) ) is now amenable to numerical quadrature, while the compensating integral has a closed-form representation, as shown in Appendix 3. Using these results in (80(80) ), we finally obtain(85)

where is the numerical distance and is the attenuation function, given in (C5(C5) ) and (C8(C8) ), respectively. It is important to note that is discontinuous when the sign of changes as the Sommerfeld pole crosses the real-axis integration path in Figure 8 – and the Sommerfeld case turns into the plasmonic case. However, the resulting discontinuity in the branch cut integral (85(85) ) is exactly compensated by the residue term that must now be included in the nonspectral representation (79(79) ).

After the pole subtraction, the integrand of is well behaved and, assuming h, is free of the oscillatory behavior plaguing the axial and radial propagation representations for large . Furthermore, due to the exponential decay factor, the rate of convergence actually improves with increasing radial distance . The latter property makes the nonspectral representation ideal for surface field computations – the only drawback being that leaky wave poles may be captured in the multilayer case.

2.7. Asymptotic evaluation of the branch cut integral

With the pole term extracted, this integrand of in (83(83) ) is smooth and rapidly decaying, due to the fact that the integration path is the steepest descent path (SDP) for h.[112] Further, in view of the exponential factor, most of the contribution to the integral comes from the vicinity of . Hence, we may approximate the integrand function by a Maclaurin expansion [56, p.400], keeping just the first two terms. Since the leading term vanishes, we obtain(86)

where and is given as²⁰(87)

Upon using (86(86) ) in (85(85) ), we obtain the asymptotic form(88)

applicable for large . In the above, we have introduced the modified attenuation function(89)

which, in view of (C10(C10) ), is as .

We may similarly derive the asymptotic form of the second branch cut integral, , where now plays the role of . However, since the Sommerfeld pole is typically located near , the pole extraction is not required in this case. The contribution from this branch cut integral is negligible compared to , except in the case of a low-loss lower half-space.

2.8. Discussion of the surface field

The surface field may be approximated by the asymptotic formula (88(88) ), provided that the second branch cut contribution is negligible and no poles are captured in the nonspectral representation (79(79) ) – as in the Sommerfeld case. If a pole is intercepted – as in the plasmonic case – its contribution must be included in (79(79) ). For convenience, we will refer to the asymptotic branch cut integral (88(88) ) as the lateral wave and the contribution from any proper (top-sheet) pole captured in the nonspectral representation (79(79) ) will be called a trapped wave. We note that the horizontal range dependence of the lateral wave (88(88) ) is influenced by and by , and can be quite complex, due to the interplay of the two terms in the square brackets.²¹ It is clear, however, that sufficiently far from the dipole, the range dependence is , and this behavior is governed by , which is dictated by the limit value of the spectral integrand at the branch point .²² As regards the near- and intermediate-zone behavior, we can be more definitive in the high-contrast case, when the lateral wave assumes the simpler form (D12(D12) ). Since the attenuation function approaches unity for , the field dependence near the dipole is is , and it evolves into dependence with increasing , in view of the asymptotic form (C10(C10) ). We may thus consider the point as a demarcation point where the transition (knee) from the to the asymptote begins, which upon using (D9(D9) ) yields(90)

The asymptotic form (88(88) ) also reveals that the Sommerfeld pole plays a role in the surface field even when not captured in the nonspectral representation (79(79) ). In fact, the exponential term in (C8(C8) ), along with the other factors that appear in (88(88) ), is exactly the ZW.[27,115] This ZW, however, can only be significant in a limited horizontal range, as it is cancelled for large , in view of the asymptotic behavior (C10(C10) ) of .²³ In the plasmonic case, where and the pole is captured and contributes an explicit residue term in (79(79) ), its influence may persist over a larger horizontal range, in view of the slowly decaying amplitude factor – provided that the exponential decay is small.

Figure 9. Integrand function of a typical Sommerfeld integral in the -plane for a lossy ground case. The Sommerfeld pole, despite its location below the branch cut, strongly affects the integrand behavior along the SIP.

3. Numerical results

We include here sample numerical results to elucidate the properties of the alternative field representations and the behavior of the surface fields for several representative cases. However, in view of the recent misguided attempts in the literature to cast doubt on the existence of the Sommerfeld pole on the top sheet,[18,19] it is perhaps appropriate to begin with an example illustrating its effect on the Sommerfeld integrals. Hence, consider a vertical dipole placed at m above a lossy ground with and conductivity mS/m at MHz (m), which results in and the Sommerfeld pole at . For this case, in Figure 9 we plot the logarithmic magnitude of the integrand of at m, where is given by (D5(D5) ), over a square region in the normalized -plane, which includes the branch point. The plot of the integrand along the original SIP, which follows the real axis along the upper edge of the branch cut, is also included. Note that the Sommerfeld pole, which is located very near the branch point on Riemann sheet and is “hidden” just below the branch cut, causes a spike in along the SIP. This spike in the vicinity of the pole, as well as the infinite derivative of at the branch point, make the original SIP problematic – although some authors still use it after annihilating the branch point singularity by either subtraction the singular term or local variable transformation.[117,118] As already mentioned in connection with the axial propagation representation, another simple remedy is a detour into the 1st quadrant of the complex -plane, which skirts the singularities and then rejoins the real-axis tail. In view of the exponential growth of the Bessel function, this detour may not veer off too far into the complex plane, but its shape is otherwise not critical. Here, we have opted for the piecewise linear rooftop contour,[90] with the height inversely proportional to , and on this path the integrals are computed by the adaptive Patterson rule.[119] On the remaining real axis tail, we apply Gaussian quadrature between the asymptotic zeros of the Bessel function and extrapolate the resulting slowly convergent partial sums by the Levin-type Mosig–Michalski transformation.[84], [120, p.102] Note that there is no need to determine the exact location of the poles of the integrand. We use this approach here as a reliable benchmark method for the radial transmission and nonspectral representations, provided that the horizontal range does not exceed a few hundred wavelengths.

As the first example, consider a vertical dipole located at m over dry ground with and mS/m at MHz (m). In this case, and the Sommerfeld pole is at . Note that this is a low-contrast, low-loss Sommerfeld case, hence the pole wave is classified as Brewster mode. Plots of vs. the horizontal distance along the surface, obtained by the radial transmission and nonspectral representations, are shown in Figure 10(a) and (b), respectively.

Figure 10. Results for a vertical dipole over dry ground at MHz.

Here and henceforth, we also include the benchmark result (indicated by triangle symbols) computed using the axial transmission representation with the modified SIP, as explained above. The -branch cut contribution is not negligible in this case and its effect is included in the results. Note that the Brewster mode is strongly excited in the radial propagation representation, but is mostly cancelled by the -branch cut contribution. In the nonspectral representation, the Sommerfeld pole is not captured and the asymptotic lateral wave is practically indistinguishable from the total field when . Both representations being rigorous, the total fields plotted in Figure 10(a) and (b) are, of course, identical. However, the nonspectral representation is much more efficient, especially when the closed-form lateral wave (88(88) ) is employed.

Continuing with the dry ground example, in Figure 11 we present nonspectral representation results showing the field behavior at a fixed distance m from a vertical dipole placed at the surface (), as the frequency is varied over three decades. Similar case was recently considered by Parise [121]. The plots of in Figure 11 illustrate the effect of neglecting the 2nd branch cut integral contribution – confirming that this is not allowed in this low-contrast case. Note that the undulations in the total field are clearly a result of the interference of the two branch cut waves. Plots of and for the same case are shown in Figure 12, where the 2nd branch cut contribution is included. Again, the closed-form lateral wave captures the surface field behavior with remarkable accuracy if the horizontal distance exceeds one wavelength.

In the next example, the vertical dipole is at m above seawater with and S/m, also at MHz. In this case, and the Sommerfeld pole is at . Note that this is a high-contrast, high-loss Sommerfeld case, hence the pole wave is classified as ZW and the -branch cut contribution may be safely neglected. The corresponding radial transmission and nonspectral representation results are presented in Figure 13(a) and (b), respectively.

Figure 11. Plot of vs. frequency at a fixed distance from the source for the dry ground case. Nonspectral representation is employed.

Figure 12. Nonspectral representation results for the dry ground case.

Figure 13. Results for a vertical dipole over seawater at MHz.

We note that the ZW is weakly excited in the radial transmission representation and makes little impact on the total field.²⁴ In the nonspectral representation, the Sommerfeld pole is not captured and the total field is well represented by the lateral wave. The knee point, where the transition from the to the asymptote begins, corresponds to , or km. To conserve space, we do not show the plots of and , which are very similar to those in Figure 13, except that the former is approximately two orders of magnitude lower in strength. For the same example, in Figure 14 we give nonspectral representation results for a horizontal (-oriented) dipole, still at m. Comparing the field plots in Figures 13(b) and 14(a), we note that the latter is much smaller in magnitude, which may be attributed to the fact that a horizontal dipole has a negative image in a conducting half-space. It is also interesting that the horizontal dipole fields exhibit a decay rate in the near zone (up to wavelengths), which in the case of and then turns into a behavior in the intermediate range. The azimuthal component, , is a TE wave, for which there is no Sommerfeld pole, and this component decays as in the entire horizontal range. Mahmoud and Antar [123] state that the TE component is dominant in the far zone, but this is not borne out by the results in Figure 14(b), where both the triangle and square symbols denote benchmark results.

Figure 14. Nonspectral representation results for a horizontal dipole over seawater at MHz.

Staying with the seawater case, we next examine the effect of an ice layer () of varying subwavelength thickness – with the dipole remaining at m above the ice surface. Similar case was recently considered by Mahmoud and Antar [123]. In Figure 15(a) and (b), we plot the total surface fields – obtained using the nonspectral representation, for the vertical and horizontal dipoles, respectively. The -branch cut integral is computed using the asymptotic lateral wave expression (88(88) ), with the contribution from the branch cut associated with seawater safely neglected. Hence, the total field comprises the lateral wave and – if the ice layer is present, the trapped mode. The excellent agreement of the asymptotic and exact benchmark results should be noted. For easy reference, we have also included in the plots the bare seawater case, which was already discussed above. With the ice layer added, the Sommerfeld pole migrates to for cm and to for cm. Hence, as soon as a thin ice layer appears, the Sommerfeld case turns into a plasmonic case and the trapped wave is included in the nonspectral representation (79(79) ). For both the vertical and horizontal dipoles, the ice cover enhances the surface field in the km range – more significantly so in the latter case. This enhancement is due to the trapped wave, which exhibits the decay rate in the range where the lateral wave decays as . For exceeding km, the exponential decay factor in the trapped wave dominates, eventually forcing it to fall below the lateral wave in magnitude – which occurs at a point where the latter follows the asymptote.²⁵ The characteristic wiggles in the total field observed near this point are due to the interference of the lateral and trapped waves when their magnitudes become comparable.

Figure 15. Nonspectral representation results for a dipole over ice-covered seawater at MHz.

Figure 16. Nonspectral representation results for the GPR example and a table of the poles included.

The next example pertains to GPR detection of leaking pipes in the shallow sub-surface at MHz (m).[47] The vertical dipole is in air at the top surface () of a three-layer half-space comprising, from top to bottom, a dry sand bar () of thickness m, an air gap of thickness cm, and water logged soil (, S/m). The nonspectral representation results for this case are plotted in Figure 16, where we also list the poles that were included. The first two poles in the table are proper surface wave poles and the remaining three are improper leaky wave poles. The former are associated with two trapped modes in the sand bar, whose interference creates the undulations observed in Figure 16. The first pole listed, which is close to the branch point, was treated as the Sommerfeld pole in the nonspectral representation. The lateral wave is seen to be negligible, except near the source, and the trapped waves dominate in the horizontal range considered. To bring the nonspectral representation result in agreement with the benchmark, the contributions from the three leaky wave poles, also captured by the deformed integration path, had to be included. These waves are heavily damped, but contribute significantly in the vicinity of the dipole.

Figure 17. Nonspectral representation results for a vertical dipole over aluminum at THz.

As the next example, consider a vertical dipole at m above aluminum with and at THz (m). In this case, and the pole is found at . Note that this is a high-contrast, high-loss plasmonic case, hence the Sommerfeld pole is captured in the nonspectral representation and the pole wave is classified as SPP. We also examine the effect of a thin, m () of polyethylene coating layer with on top of Al.[45] The pole is now found at and it contributes a trapped wave. The nonspectral representation results for the bare and coated Al configurations are presented in Figures 17(a) and (b), respectively. In the former case, the SPP is weakly excited and although its contribution slightly exceeds in magnitude the lateral wave in the intermediate horizontal range, it combines with it destructively, so there is no noticeable enhancement in the total surface field magnitude. The propagation length and depth of penetration in air of the SPP are found from (66(66) ) as m and cm, respectively. Hence, the propagation length is very large and the wave is loosely bound to the surface. However, since the contribution of this mode to the total surface field is insignificant, we contend that it would be futile to try to detect its presence experimentally, as recently endeavored by Jeon and Grischkowsky [31]. The knee point of the lateral wave is found at , or m, whereas the Al sheet used in the surface wave measurements mentioned above was less than one meter long. In the coated Al case, the trapped wave with the dependence is strongly excited and it dominates the total field over a wide horizontal range, until the exponential decay factor makes it fall below the lateral wave, which is in transition to the asymptote at this point. This wave is much more tightly bound to the surface than the SPP in the bare Al case.²⁶ As a result, a larger fraction of the mode energy is concentrated in the lossy conductor, which accelerates the onset of the exponential decay, as compared to the bare metal case.

Figure 18. Nonspectral representation results for a vertical dipole over silver at nm.

Figure 19. Nonspectral representation results for the near-infrared Cr case.

In the next example, the vertical dipole is placed at nm above silver half-space at nm in the visible band. By fitting the experimental dielectric function data for Ag [124] with a low-order partial fraction (PF) model,[125] we find and the Sommerfeld pole at . Note that this is a high-contrast, low-loss plasmonic case, hence the pole is captured in the nonspectral representation and the pole wave is classified as Fano mode. The contribution of the branch cut is negligible, because has a large negative imaginary part. We also consider the effect of a thin, nm silver sulfide tarnish layer, which builds when silver is exposed to air. For AgS, by fitting the tabulated data [126] with the PF model,[125] we find and a numerical search reveals that the pole has now migrated to and it contributes a trapped mode. The nonspectral representation results for the bare and tarnished Ag configurations are presented in Figure 18(a) and (b), respectively. Note that, due to the low losses of Ag in the visible range, the SPP is strongly excited - both in the bare metal and coated cases.²⁷ The propagation length and depth of penetration in air of the Fano mode are found from (67(67) ) as m and m, respectively. Hence, the propagation length is large and the mode is confined to the surface, which is exploited in plasmonic sensor applications.[127] In the coated Ag case, the trapped mode is more heavily damped, due to the additional losses in the tarnish layer, with the concomitant decrease in the propagation range.

Figure 20. PF model fit [125] of the tabulated dielectric function of Cr,[128] including the range where .

As the final example, we consider a vertical dipole located at nm above a chromium half-space, and compute the surface field at a fixed distance m from the dipole as the frequency (photon energy) is varied in the near-infrared band. The nonspectral representation results for this case are shown in Figure 19. Since in the selected frequency range of Cr passes from negative values to positive and back,[35] as indicated in Figure 20, the Sommerfeld pole migrates across the vertical branch cut and there is a transition from the plasmonic case (the pole is captured and contributes a SPP) to the Sommerfeld case (the pole is not captured, but the ZW is implicit in the groundwave) and back. Hence, the total field plot in Figure 19 includes the lateral wave and the SPP in the plasmonic range, and just the lateral wave in the Sommerfeld range. This plot is continuous, which confirms that the jumps in the attenuation function (C8(C8) ) as the Sommerfeld pole crosses the vertical branch cut (or the real axis in the -plane – see Figure 8) are exactly compensated by the SPP, which is either included or not in nonspectral representation (79(79) ). This example thus underscores the close relationship between the ZW and SPP.

4. Summary and conclusions

We have revisited the classical half-space problem of Sommerfeld, with extensions to a layered half-space and plasmonic media and a focus on the surface fields away from the dipole. A new ab initio solution of this problem is presented and alternative field representations – termed axial transmission, radial transmission, and nonspectral – are derived and their relative merits discussed. Using the nonspectral representation as the point of departure, a closed-form asymptotic expression for the surface field is derived. Based on this expression, the surface field behavior is discussed and the rôle of the ZW and SPP is elucidated. The theoretical development is supported by representative numerical examples, which cover a wide spectrum – from radio waves to visible light. In all cases considered, the closed-form asymptotic formulas are remarkably accurate for . The main points of this study may be recapitulated as follows:

• The total surface field of a Hertzian dipole may be decomposed into two branch cut contributions associated with the upper and lower media plus possible pole waves. The branch cut contribution associated with the upper medium (air) is called a lateral wave and it typically dominates the total surface field; the contribution from the second branch cut (associated with the lossy half-space) is negligible if the media contrast is sufficiently high.

• In the case of a homogeneous, lossy half-space, the Sommerfeld pole does exist on the top (proper) Riemann sheet and is always captured in the radial transmission representation (using hyperbolic branch cuts).

• Depending on the real part of the dielectric constant, the pole may be located to the left (Sommerfeld case) or to the right (plasmonic case) of the branch point associated with the upper medium (air).

• In the Sommerfeld case, the pole wave is slow (meaning that its phase velocity exceeds the speed of light) and is called a ZW, and in the plasmonic case, this wave is fast and is called a SPP. In the low-loss plasmonic case, the SPP is also referred to as a Fano mode. In the multilayer half-space case, the proper pole waves are called trapped modes. All pole waves are characterized by a small, , decay rate in the small to moderate horizontal ranges.

• In the nonspectral representation (using vertical branch cuts), the pole is captured in the plasmonic case, but not in the Sommerfeld case.

• The appearance or not of the ZW in a particular surface field representation does not prove or disprove the physical reality of this wave.

• In the Sommerfeld case, the ZW always contributes to the total surface field, either explicitly (in the radial propagation representation) or implicitly (in the nonspectral representation). However, this wave may only be significant in the near to intermediate horizontal range, and is never a dominant contributor to the total surface field. Consequently, it is not feasible to experimentally detect the presence of the ZW in the surface field of a simple antenna.

• The range dependence of the lateral wave is determined by a complex parameter called numerical distance, which depends not only on , but also on the dielectric constant of the lossy half-space. The lateral wave decays as near the dipole and the decay rate increases with increasing , eventually reaching the asymptote. The transition (knee) to the asymptote begins at the value given by .

• In the case of noble metals in the THz (far-infrared) range, the SPP is weakly excited and is weakly bound to the surface, and its contribution to the total surface field is negligible. When the metal is coated with a thin, low-loss dielectric layer, the SPP becomes a trapped mode, which is strongly excited and much more tightly bound to the surface.

• In the case of noble metals in the visible range, where the losses are very small, the SPP (Fano mode) is strongly excited and it may dominate the surface field behavior over a thousand-wavelength range. When a Fano or a trapped mode is present, the total surface field initially exhibits the decay, until the exponential decay factor eventually brings it below the asymptote of the lateral wave.

• The ZW and SPP are closely related – they satisfy the same resonance condition – and one may evolve into the other with the change of frequency.

Additional information

Funding

This work was supported in part by École Polytechnique Fédérale de Lausanne (EPFL), Switzerland.

Notes

No potential conflict of interest was reported by the authors.

Dedicated to the memory of Professor James R. Wait.

1 The horizontal dipole case was first solved by von Hoerschelmann [5] in a PhD thesis supervised by Sommerfeld.

2 The subtle error (just a wrong sign in one of three terms [8]) – never acknowledged by the author – was first publicized by Norton in 1935.[9]

3 One early theory, promulgated in several papers by Kahan and Eckart and soon debunked by Bouwkamp [15] and Ott [16], was that Sommerfeld erroneously included the ZW in the field representation when he deformed the integration path around the hyperbolic branch cuts. This controversy was thoroughly covered by Baños [17, Section 4.10], who ended with a hopeful plea that there would be no more papers attempting to disprove the existence of the ZW term in the Sommerfeld field representation. His wish was nor granted, however, as we again see misguided attempts to deny the existence of the Sommerfeld pole on the top Riemann sheet.[18,19] Another, more long-lasting hypothesis was that Sommerfeld committed a sign error in choosing the wrong branch of the square root of the complex numerical distance parameter that appears in his formula. However, in 2004 Collin [14] found no sign inconsistency in Sommerfeld’s 1909 derivation and concluded that “there are inherent limitations in Sommerfeld’s solution, but they are not caused by a sign error” – where by “limitations” he presumably meant the exceedingly small radii of convergence of the series expansions employed by Sommerfeld, but he did not point out any specific mathematical illegality that lead to the sign error in the asymptotic groundwave formula. Hence, it would appear that the exact cause of the sign error in Sommerfeld’s 1909 paper remains unknown to this day. This is a moot point, however, which now squarely falls within the realm of the history of science, since the correct asymptotic solution has long been available, as indicated above.

4 However, this last formulation was later found to be incomplete, as it does not properly account for the trapped wave.[14,44]

5 So named, because the integration path is not confined to the proper Riemann sheet.

6 The time convention is implied.

7 Source coordinates are distinguished by primes.

8 Unit vectors are distinguished by carets.

9 Note that is the zero-order Bessel function.[54].

10 Here, we use (17(17) ) and the relation , where the prime denotes differentiation with respect to the argument of the Bessel function.

11 Our proof of the existence of the Sommerfeld pole on the top sheet is similar to that originally given by Sommerfeld [1,3] and later by other authors [20,72,73], [74, p.355], except that we have extended it to plasmonic media with . We should also mention here a note by Gogoladze [75], purporting to prove the nonexistence of the Sommerfeld pole on the top sheet. Also, Schelkunoff [76, p.430] gives a proof that there is no pole on the top sheet – however, he seems to assume a lossless half-space. This proof has been uncritically embraced by Sarkar et al. [18,77], who also deny the existence of the Sommerfeld pole on the top sheet. The lossless case is also tacitly assumed by Chew [78, p.99], when he states that only in the plasmonic case is the pole located on the top sheet.

12 This name derives from the fact that the (original) ZW [7] may be set up by assuming a plane wave impinging in air on the lossy half-space at the (complex) Brewster angle,[79] where the latter is determined by the value of the transverse propagation constant , which makes the numerator of the reflection coefficient (A10(A10) ) vanish. Hence, only the incident and transmitted waves exist at this angle, and both must decay in magnitude away from the interface, in order to constitute a surface wave. Since the incident wave propagates toward the interface, with the -dependent factor , and it increases in magnitude as it does so, we must have . The transmitted wave in the lower half-space propagates away from the interface and vanishes for , hence . Consequently, is a root of the numerator of the reflection coefficient (A10(A10) ) on Riemann sheet . Note that in the Sommerfeld half-space problem, the ZW arises as a resonant solution, i.e. it is a mode that exists in the absence of any incident field, and thus the reflection coefficient looking into the half-space must be infinite.[80] To set up such a mode, we postulate waves propagating away from the interface and vanishing for . Hence, the propagation constant is a root of the denominator of the reflection coefficient on sheet , where both and .

13 Incidentally, one should resist the temptation to ascribe any physical significance the the particular threshold numbers , as they are based on approximate expansions (64(64) ) and (66(66) ).[81].

14 Incidentally, the overriding difficulty is not due to the oscillatory integral tail, which can be effectively computed by partition-extrapolation methods,[84,86,91,93] but due to the rapid oscillations of the integrand on the initial, finite-length path segment.

15 The Schelkunoff integral was recently misused in this way by Sarkar et al. [18,19] to derive a spurious radial transmission representation – called the “Schelkunoff formulation” and claimed to be rigorous – which is missing both the ZW and the second branch cut contribution. The fallacy of this formulation was already pointed out by Janaswamy [20], but the criticism was rejected by the authors [21].

16 This is contrary to the recent claim of Dyab et al. [103] that there is no relationship “whatsoever” between SPP and ZW – since, according to these authors, the latter “does not exist.”

17 The group velocity of ZW also exceeds the speed of light,[104] but not its energy transport velocity.[81].

18 Hence, the radial transmission representation is not quite the panacea for the oscillatory integrand behavior it is sometimes purported to be in the literature.[18].

19 Note that the subtracted pole term has been purposely selected to have a zero limit at the origin, where also vanishes.

20 This limit is discussed in more detail in Appendix 4.

21 It should be borne in mind here that and are also -dependent, unless .

22 For TE fields, where the Sommerfeld pole is absent (assuming nonmagnetic media), in (88(88) ) and the lateral wave exhibits only the inverse square -dependence.

23 As a result of the already mentioned elusive error in Sommerfeld’s 1909 paper,[1] there was no asymptotic cancellation of the ZW in his original groundwave formula. Hence, the deleterious effect of this error was that the groundwave seemed to be an important factor in the radio wave propagation along the surface of the earth.[116] It is interesting that Sommerfeld inadvertently arrived at a formula applicable in the plasmonic case. Again, Sommerfeld later gave the correct expression in his 1926 paper.[2].

24 We should mention here the recent attempts of Petrillo et al. [122] to design an efficient ZW launcher, in hopes of extending the range of the maritime surface wave radar. Using a two-dimensional radial transmission representation for that purpose, they arrived at some dipole array configurations, which they describe as “realistic, but still unrealizable.”

25 The knee point, which was at km in the bare seawater case, occurs much sooner when ice is present – it is at km for a cm ice layer, and as close as km when the ice thickness is doubled to cm.

26 Gong et al. [45] call this phenomenon the “THz surface wave collapse.” Note that they refer to the SPP as the “Zenneck THz surface wave.”

27 This is in contrast with the bare Al in the THz range, where the SPP excitation was weak.

29 Henceforth, we omit the superscript for notational simplicity.

30 Note that these branch cuts are in the 1st and 3rd quadrants if the time convention is used.

31 Here, a lossless medium should be treated as a limit of the corresponding lossy medium as losses approach zero.

32 Note that on the hyperbolic extensions of the branch cuts, indicated by the dotted lines.

33 Sommerfeld called this the similarity law (Ähnlichkeitsgesetz) of wireless telegraphy.[1].

34 We henceforth omit the superscript for notational simplicity.

35 Henceforth, square roots with positive real parts are implied.

Appendix 1

Elements of TL analysis

We refer to the layered medium and its TL analog illustrated in Figure 3, where the th TL section – assumed to extend between and along the -axis, has the propagation constant and characteristic impedance , which are given by (28(28) ) with the layer subscript appended, i.e.(A1)

Here, is the wavenumber of layer and its intrinsic impedance is . Since our focus is on the surface fields in air, we only require the TL Green functions corresponding to the upper half-space comprising the source, and these Green functions are readily found from a simple network analysis, which we present next. The general case, where the source and field points may be in any layer, can similarly be treated.[56,129–132]

The voltage and current in the th section that is source free satisfy the TL equations²⁹(A2)

so that may be found from(A3)

and follows from the first of the TL Equations (A2(A2) ). The general solution of (A2(A2) ) may thus be expressed as(A4)

where the upper (lower) expression on the left side corresponds to the upper (lower) sign on the right and are yet-to-be-determined amplitudes of the up- and down-propagating voltage waves, referred to the lower and upper terminals of section , respectively. For observation points down from the source location, it is convenient to eliminate in favor of and express (A4(A4) ) as [132](A5)

where(A6)

is the down-looking voltage reflection coefficient at point and denotes its value at the lower terminal location . The impedance looking down the line at any point in section may now be found as(A7)

and its value at the lower terminal will be denoted by (no argument). Upon enforcing the continuity of the impedance at the interface between sections and , we readily obtain the backward recursion(A8)

with (since there is no reflected wave in the bottom TL section) and(A9)

where is the layer thickness. This recursion allows us to compute, in particular, the reflection coefficient(A10)

looking down from the upper medium into the layered half-space. If the latter is homogeneous, we may use (A10(A10) ) directly, since in this case .

In Section 1, which comprises the source point , we may write (using )(A11) (A12)

where the coefficients depend on the source configuration. For the unit-strength current source excitation (), upon enforcing the jump conditions at , we obtain(A13)

and upon substituting (A13(A13) ) into (A11(A11) )-(A12(A12) ), we find the corresponding voltage and current for as(A14) (A15)

where the signs correspond to and the subscript now refers to the source type. For the unit-strength voltage source excitation (), we similarly find(A16)

and the corresponding current and voltage for are given as(A17) (A18)

The duality between (A14(A14) )–(A15(A15) ) and (A17(A17) )–(A18(A18) ) should be noted: we obtain the latter from the former upon replacing the impedances by the corresponding admittances (which has also the effect of changing the sign of the reflection coefficient).

Figure A1. Riemann sheets associated with in a lossy case.

Appendix 2

Fundamental branch cuts and the Riemann surface

The propagation constants are defined by the square root function (A1(A1) ), which is double valued with branch points at .[17, p.57], [133, Appendix B], [134, Section 1–6] For a unique specification of this square root in the complex -plane, we impose the condition , which defines a two-sheeted Riemann surface [78, Section 2.2.3] with hyperbolic (or fundamental) branch cuts in the 2nd and 4th quadrants.³⁰ We refer to the sheet on which the condition is satisfied as the proper sheet, and the sheet on which this condition does not hold as the improper sheet. This nomenclature derives from the fact that on the proper sheet the function is consistent with the radiation condition, i.e. it vanishes as .³¹ The two Riemann sheets with fundamental branch cuts are illustrated for a lossy medium in Figure A1, where the only way to pass from one sheet to another is by crossing one of the branch cuts, on which .³² If we exclude active [135] and double negative media,[61] the branch point is always located in the 4th quadrant. In the lossless limit, as , approaches the positive real axis and the hyperbolic branch cuts follow the real and imaginary axes. It should be noted here that the choice of branch cuts is wholly arbitrary, as long as they do not cross the real axis path of the Sommerfeld integrals in (46(46) )–(54(54) ). The convenience of the fundamental branch cuts is that the radiation condition is satisfied, or not, on an entire sheet of the Riemann surface. However, there are instances where vertical cuts are preferable – even though the spectral integrands are no longer bounded for all values of on the new top sheet.[4, p.251], [136], [17, p.55], [111]

Appendix 3

Evaluation of the integral

We first transform the integral in (83(83) ) as(C1)

where we have evaluated the first integral in closed form. The second integral, , has been much studied in the literature in the context of the modified saddle-point method.[137, p.240], [56, p.402], [60, Appendix of Chapter 15], [14,138]. Here, using a slightly different and simpler approach, we will express in terms of the Faddeeva (plasma dispersion) function(C2)

for which computer subroutines are readily available.[139–141] This function has the integral forms [54, Equation 7.7.2](C3)

and we note that can be related to the rightmost integral in (C3(C3) ), if we introduce a new variable of integration and set . However, since may be negative (in the Sommerfeld case – see Figure 8), whereas (C3(C3) ) is restricted to , it is necessary to first extend its validity to the lower half-plane by analytic continuation. This is accomplished by adding to the integral times the pole residue when (or half of the pole residue when ). As a result, we obtain(C4)

where the parameter is the Sommerfeld numerical distance [2], [3, p.933], given as(C5)

and is a step function, defined as(C6)

Finally, upon using (C4(C4) ) in (C1(C1) ), we obtain(C7)

where we have introduced(C8)

which, in view of the properties of ,[54, Equation 7.4.3] may also be expressed as(C9)

We note that in the Sommerfeld case reduces to the Sommerfeld attenuation function.[142] It is significant that the behavior of this attenuation function is not determined by or alone, but by the numerical distance parameter , which also includes the effect of the material constants.³³ It can be shown that(C10)

irrespective of the sign of .[27,123,143,144] Note that as .

Appendix 4

Supplementary analytical details

In the main body of this paper, we have presented a general methodology for computing the surface fields (46(46) )–(54(54) ). Here, we fill in some pertinent analytical details. For the sake of brevity, we focus on one representative Sommerfeld integral, associated with the vector potential of a vertical dipole – the case originally treated by Sommerfeld [1]. The treatment of the Sommerfeld integrals in (46(46) )–(54(54) ) is very similar, with some simplifications in the TE () case, where there is no Sommerfeld pole for a homogeneous lower half-space.

Hence, focusing on (43(43) ), we consider the spectral integrand function(D1)

where we have set and . As already mentioned in Section 2.3, in the general multilayer case the poles of this function can only be determined by a numerical search in the complex plane for the roots of . The corresponding residues may then be computed using (73(73) ). However, since analytical evaluation of becomes cumbersome beyond two layers, we resort to the automatic differentiation already implemented in the root search procedure.

The folded , needed in the alternative representations, is readily found as³⁴(D2)

where the branch of is implied and the branch of is selected differently, depending on whether the radial transmission or nonspectral representation is used. In the former case, the branch is always selected, while in the latter case, the or branch should be enforced, depending on whether (Sommerfeld case) or (plasmonic case), respectively – see Figure 7. In deriving (D2(D2) ), we have used the fact that the substitution implies . The branch point limit form, needed to compute in (87(87) ), readily follows as(D3)

where we have used the fact that near and introduced(D4)

In the above, since at , it is convenient to compute as .

If the lower half-space is homogeneous, we substitute , and (D1(D1) ) becomes(D5)

In this case, the Sommerfeld pole is given by the closed-form expression (58(58) ) and the corresponding residue follows from (73(73) ) as³⁵(D6)

where use has been made of (62(62) )–(63(63) ). The homogeneous half-space counterparts of (D2(D2) ) and (D3(D3) ) are readily found as(D7) (D8)

Finally, certain simplifications obtain for high-contrast media (). In this case, we may use the approximations(D9) (D10)

where we have replaced by in the argument of . We also find that(D11)

which may be used to simplify (88(88) ) as(D12)

where we have also replaced by unity, assuming large . If also , we may further replace the attenuation function by its leading asymptotic term as(D13)

which makes explicit the asymptotic range dependence of the lateral wave (D12(D12) ). For and , (D12(D12) ) reduces to the Sommerfeld groundwave formula. [2]