Electromagnetic surface wave propagation in a metallic wire and the Lambert $W$ function

TL;DR

This paper revisits Sommerfeld's classical solution for electromagnetic surface wave propagation in metallic wires, analyzing the convergence of his iterative method and comparing it with numerical and Lambert W function solutions.

Contribution

It provides an elementary convergence analysis of Sommerfeld's iterative solution and compares it with numerical and Lambert W function approaches.

Findings

Convergence of Sommerfeld's iterative method is established.

Lambert W function offers an alternative solution approach.

Comparison shows the effectiveness of different solution methods.

Abstract

We revisit the solution due to Sommerfeld of a problem in classical electrodynamics, namely, that of the propagation of an electromagnetic axially symmetric surface wave (a low-attenuation single TM$_{01}$ mode) in a cylindrical metallic wire, and his iterative method to solve the transcendental equation that appears in the determination of the propagation wave number from the boundary conditions. We present an elementary analysis of the convergence of Sommerfeld's iterative solution of the approximate problem and compare it with both the numerical solution of the exact transcendental equation and the solution of the approximate problem by means of the Lambert $W$ function.

Figures

Equations

$${\nabla }^2\mathbf{E}-\sigma \mu \frac{\partial \mathbf{E}}{\partial t}-\epsilon \mu \frac{{\partial }^2\mathbf{E}}{\partial t^2}=0,$$

$$\mathbf{E}(\mathbf{r},t)={\mathbf{E}}^0(r,\theta )\varphi (z,t)={\mathbf{E}}^0(r,\theta )\exp [i(hz-\omega t)]$$

$$\nabla \times \mathbf{H}=\sigma \mathbf{E}+\frac{\partial (\epsilon \mathbf{E})}{\partial t}(\text{Ampeˋre’s circuital law})$$

$$\nabla \times \mathbf{E}=-\frac{\partial (\mu \mathbf{H})}{\partial t}(\text{Faraday’s law of induction})$$

$$E_r^0(\hat{\rho })$$

$$\sqrt{\frac{\mu }{\hat{\epsilon }}}{H}_{\theta }^0(\hat{\rho })$$

$${\hat{k}}^2$$

$$\hat{\rho }$$

$$\hat{\epsilon }$$

$$\hat{k}(\omega )=\sqrt{\hat{\epsilon }\mu }\omega$$

$$\frac{d^2{E}_z^0}{d{\hat{\rho }}^2}+\frac{1}{\hat{\rho }}\frac{d{E}_z^0}{d\hat{\rho }}+E_z^0=0.$$

$$E_z^0(\hat{\rho })$$

$$E_z^0(\rho )$$

$$E_z^0(\rho )=B_1{H}_0^{(1)}(\rho )+B_2{H}_0^{(2)}(\rho ),a<r<\infty .$$

$$H_{\nu}^{(1,2)}(\rho)\sim\sqrt{\frac{2}{\pi\rho}}\exp\big{[}\pm i(\rho-\tfrac{1}{2}\nu\pi-\tfrac{1}{4}\pi)\big{]},$$

$$S_{r}(r,t)=-\operatorname{Re}\big{\{}E_{z}^{0}(r)\phi(z,t)\big{\}}\operatorname{Re}\big{\{}H_{\theta}^{0}(r)\phi(z,t)\big{\}}.$$

$$U_{r}(r,t)=4\frac{r}{\rho}\sqrt{\frac{\varepsilon_{0}}{\mu_{0}}}\frac{k}{\sqrt{k^{2}-h^{2}}}\Big{[}B_{1}^{2}\cos^{2}(\rho-\tfrac{1}{4}\pi+\text{phase})-B_{2}^{2}\cos^{2}(\rho-\tfrac{1}{4}\pi-\text{phase})\Big{]}.$$

$$\begin{split}\sqrt{\frac{\mu_{0}}{\varepsilon_{0}}}\,\frac{\sqrt{k^{2}-h^{2}}}{k}\,&\frac{H_{0}^{(1)}\big{(}a\sqrt{k^{2}-h^{2}}\big{)}}{H_{0}^{(1)\,\prime}\big{(}a\sqrt{k^{2}-h^{2}}\big{)}}=\\ &=\sqrt{\frac{\mu}{\hat{\varepsilon}}}\,\frac{\sqrt{\hat{k}^{2}-h^{2}}}{\hat{k}}\,\frac{J_{0}\big{(}a\sqrt{\hat{k}^{2}-h^{2}}\big{)}}{J_{0}^{\,\prime}\big{(}a\sqrt{\hat{k}^{2}-h^{2}}\big{)}},\end{split}$$

$$\hat{k}\approx \frac{1+i}{\sqrt{2}}\sqrt{\sigma \mu \omega }=(1+i)\kappa (\omega ).$$

$$h=20.960+7.516\times 10^{-4}i.$$

$$\frac{\omega }{h}=\frac{\omega }{k}\cdot \frac{k}{h}=c[1-(7.438+3.586i)\times 10^{-5}],$$

$$J_{\nu }(\hat{\rho })\sim \sqrt{\frac{2}{\pi \hat{\rho }}}\cos (\hat{\rho }-\frac{1}{2}\nu \pi -\frac{1}{4}\pi ),$$

$$\sqrt{\frac{\mu}{\hat{\varepsilon}}}\,\frac{\sqrt{\hat{k}^{2}-h^{2}}}{\hat{k}}\,\frac{J_{0}\big{(}a\sqrt{\hat{k}^{2}-h^{2}}\big{)}}{J_{0}^{\,\prime}\big{(}a\sqrt{\hat{k}^{2}-h^{2}}\big{)}}\simeq i\sqrt{\frac{\mu}{\hat{\varepsilon}}}.$$

$$H_0^{(1)}(\rho )$$

$$H_{\nu }^{(1)}(\rho )$$

$$\begin{split}\sqrt{\frac{\mu_{0}}{\varepsilon_{0}}}\,&\frac{\sqrt{k^{2}-h^{2}}}{k}\,\frac{H_{0}\big{(}a\sqrt{k^{2}-h^{2}}\big{)}}{H_{0}^{\prime}\big{(}a\sqrt{k^{2}-h^{2}}\big{)}}\simeq\\ &\simeq\sqrt{\frac{\mu_{0}}{\varepsilon_{0}}}\,\frac{a(k^{2}-h^{2})}{k}\,\ln\Big{(}\frac{e^{\gamma}a\sqrt{k^{2}-h^{2}}}{2i}\Big{)}.\end{split}$$

$$u=\Big{(}\frac{e^{\gamma}a\sqrt{k^{2}-h^{2}}}{2i}\Big{)}^{2}\quad\mathrm{and}\quad v=\frac{e^{2\gamma}}{2i}\sqrt{\frac{\varepsilon_{0}\mu}{\hat{\varepsilon}\mu_{0}}}\,ka,$$

$$u\ln u=v.$$

$$x\ln x=a$$

$$\ln x=\ln (-a)-\ln (-\ln x),$$

$$\ln x_2=\ln (-a)-\ln (-\ln x_1),\text{ or }{x}_2=\frac{a}{\ln x_1}.$$

$$x_1=-a,x_{n+1}=\frac{a}{\ln x_n},n\ge 1.$$

$$\Vert T(y)-T(x)\Vert \le c\Vert y-x\Vert$$

Full text

Keywords: Applied classical electromagnetism; Sommerfeld-Goubau wave; axial cylindrical surface wave; approximate iterative solution; Banach contraction principle

I Introduction

The propagation of electromagnetic (EM) waves across metallic wires was a possibility of considerable interest in the second half of the XIX century, due both to the advent of the modern theory of electromagnetism and the deployment of the early electric telegraph networks.Huurdeman2003 ; Nahim2002 First attempts to investigate EM waves on wires were made, among others, by Hertz during 1886–1887 using spark gaps as sources of EM radiation. Some of the main questions at that time were whether EM waves propagate with a finite velocity, as predicted by Maxwell’s equations, and, if so, whether they propagate at the same speed on air and other materials. Underlying these questions was the larger issue of the similarity of EM waves to light. Hertz could solve neither the theoretical nor the experimental problems involving wires, despite his many other outstanding contributions to the field.Hertz1893 ; Sengupta2003

In 1899, Sommerfeld tackled the problem of the propagation of a circular symmetric EM wave on a thin conducting wire using the full set of Maxwell’s equations, demonstrating that such a cylindrical geometry could support what has become known as a surface wave.Sommerfeld1899 He showed, by means of an example (see Sec. II), that the attenuation of the surface EM waves at high frequencies was too large to make them useful for telecommunications. In the ensuing years, Sommerfeld’s work was extended to include noncircular symmetric modes and to treat coated wires and dielectric cylinders,Harms1907 ; Hondros1909 ; Debye1910 ; Schriever1920 untill by the mid-1930s the general interest shifted to the study of EM cavities (metal boxes) and waveguides (hollow metal tubes) as more practical and efficient devices to store and transmit EM energy across space.Southworth1936 ; Southworth1937 ; Collin1999 ; Carson1936 ; Barrow1937 ; Packard1984

In $1950$, a German physicist by the name of Georg J. E. Goubau, brought to USA from Germany after World War II as part of Operation Paperclip, resumed the investigation of the “Sommerfeld waves” with regard to their practical application to transmission lines.Goubau1950 Goubau was able to show, both theoretically and experimentally, that EM surface waves can be generated on wires by means of coaxial flared horn-shaped launchers, making the otherwise theoretical waves practical. His main result, though, was the demonstration that a sheathed or subwavelength corrugated wire would perform better than a naked wire in the transmission of low frequency ($<1$ THz) EM energy along the line. This is not because of substantially less ohmic or dielectric losses, but because the modified surface of the wire has a much lower effective conductivity than the metal alone, reducing the extension of the radial electric field outside the wire (field confinement), which for low frequencies and large wire diameters can be large. For example, a wire of radius $a=10$ mm carrying a surface wave of $1$ GHz propagates $75\%$ of its power into an area of radius $\sim 1.5$ m around the wire,Goubau1950 not only diverting energy but also potentially coupling with the environment and compromising the entire engineering solution by requiring too much clearance around the transmission line. Surface EM waves on wires have since been known in radio engineering as Sommerfeld-Goubau waves. BarlowCullen1953a ; BarlowCullen1953b ; Roberts1953 ; Barlow1958 Such waves have repeatedly been considered as a cheap alternative to long-distance, beyond-the-horizon transportation of EM energy—for example, for last mile simplex or half-duplex data communication over existing overhead electric power lines—in conditions under which microwave and coaxial links become inefficient or costly.FriedFern2001 ; WangMitt2004 ; JeonZhangGris2005 ; Akalin2006

In this paper we revisit Sommerfeld’s solution of the problem of an EM axially symmetric surface wave (to wit, a pure transverse magnetic, TM01 mode) propagating in a cylindrical metallic wire. Of particular interest to us is the solution of the boundary condition for the EM fields that lead to a transcendental equation that Sommerfeld solved approximately by means of a practical and effective device. We explore the iterative procedure devised by Sommerfeld to solve the transcendental equation and briefly investigate its convergence. As we shall see, the solution of the approximate physical system can immediately be given in terms of the (complex-valued) Lambert $W$ function, a most interesting function that gained a lot of attention in the last few decades. Both solutions (Sommerfeld’s and the one in terms of the Lambert $W$ function) are compared with the exact numerical solution of the boundary condition in a concrete case. Some remarks on the physical significance of the solutions are made.

II Circular symmetric EM wave

in a metallic wire

Comprehensive treatments of EM waves particularly suited to our discussion can be found in Refs. Stratton1941, , Chap. IX and Sommerfeld1952, , $\S{20}$–$\S{25}$, while for definitions and properties of special functions we refer the reader to Refs. SommerfeldVI, ; Arfken7e, ; NIST2018, . Except for one mention to electron-volts in Sec. III, we use SI units throughout.

Maxwell’s equations in a conducting medium without free charges lead to the following wave equation for the electric field $\mathbf{E}=\mathbf{E}(\mathbf{r},t)$,

[TABLE]

with an identical equation for the magnetic field intensity $\mathbf{H}=\mathbf{H}(\mathbf{r},t)$, where $\varepsilon$ and $\mu$ are, respectively, the time-independent electric permittivity and magnetic permeability of the material and $\sigma$ its electrical conductivity. Outside of the conductor $\sigma=0$, and in vacuum (and to a very good approximation also dry air) $\varepsilon=$ $\varepsilon_{0}=$ $1/\mu_{0}c^{2}$ F m*-1* $\simeq 8.854\times 10^{-12}$ F m*-1* and $\mu=$ $\mu_{0}=$ $4\pi\times 10^{-7}$ H m*-1*, with $c=1/\sqrt{\varepsilon_{0}\mu_{0}}=299\,792\,458$ m s*-1* the exact speed of light.

Let us consider a long and thin metallic wire of cylindrical symmetry and constant circular cross section of radius $a$ embedded in a dielectric medium, and assume that there is only a harmonic oscillating electric field along the direction of propagation $z$, i. e., a pure transverse magnetic (TM) mode. This kind of EM wave can be excited on one end of the wire provided that its source is infinitely far removed, or by some device specifically designed for the task, as the coaxial horn-shaped launchers devised by Goubau in the 1950s (essentially mode transducers transforming the TEM wave of a coaxial cable to the TM01 wave of the surface waveguide).Goubau1950 ; RTVNews1955 Following the usual practice, we introduce the complex wave number $h$ (to be determined from the boundary conditions) to express both the propagation and the attenuation of the EM wave across the wire in the $z$-direction through the time-dependent phase factor $\phi(z,t)=\exp[i(hz-{\omega}t)]$ with a real angular frequency $\omega$, such that

[TABLE]

and similarly for $\mathbf{H}(\mathbf{r},t)$. The imaginary part of $h$ must be positive if the wave is to attenuate; we shall also take the real part of $h$ positive, corresponding to the phase of the wave propagating in the positive $z$-direction. Because of the cylindrical symmetry of the problem, the time-independent field vector $\mathbf{E}^{0}(r,\theta)$ depends only on the $r$-coordinate, $\mathbf{E}^{0}(r,\theta)=\mathbf{E}^{0}(r)$, and similarly for $\mathbf{H}^{0}(\mathbf{r})$. Within these settings, Maxwell’s equations

[TABLE]

imply that an electric field $E_{z}$ varying along the $z$-direction induces (by Ampère’s law) a circulating magnetic field $H_{\theta}$ that, in turn, induces (by Faraday’s law) an electric field with components in the $r$ and $z$ directions—and the three fields $E_{z}$, $H_{\theta}$, and $E_{r}$ suffice to describe the physics of the problem. The other components $E_{\theta}$, $H_{z}$, and $H_{r}$ vanish together. The geometry of these fields is depicted schematically in Fig. 1.

The physics of the problem indicates that everything must depend on $E_{z}$, the sole forcing field. Indeed, perusal of Eqs. (3) shows that $E_{r}^{0}$ and $H_{\theta}^{0}$ can be obtained from $E_{z}^{0}$ as

[TABLE]

where we have introduced the complex quantities

[TABLE]

The hats indicate the quantities inside the wire ($0<r<a$), where $\sigma\neq 0$; outside of the wire ($r>a$) they become real and lose their hat. For definiteness, we always take the imaginary part of $\hat{\rho}$ to be positive, which unambiguously defines the sign of the square root in Eq. (5b). The introduction of the complex dielectric constant $\hat{\varepsilon}=\hat{\varepsilon}(\omega)$ allows us to unify the treatment of the problem for dielectric media, conducting or not. In particular, it allows us to write the complex dispersion relation

[TABLE]

everywhere. If the medium is nonconducting, $\hat{\varepsilon}=\varepsilon$ and the usual dispersion relation $k(\omega)=\sqrt{\varepsilon\mu}\,\omega=(n/c)\omega$ recovers, with $n\geq 1$ the refractive index of the medium. Definition (5c) was employed in the simplification of Eq. (4b).

We must find $E_{z}^{0}(\hat{\rho})$ to complete the description of the system. It follows from Eq. (1) that $E_{z}^{0}$ obeys

[TABLE]

Equation (7) is a particular case ($\nu=0$) of Bessel’s differential equation, for which the last term on the left-hand side in the general case would read $(1-\nu^{2}/\hat{\rho}^{2})E_{z}^{0}$, with $\nu$ a complex number. The physical solutions of Eq. (7) are given by

[TABLE]

where $A$ and $B$ are arbitrary real constants and $J_{0}(\hat{\rho})$ and $H_{0}^{(1)}({\rho})$ are, respectively, Bessel and Hankel functions of the first kind. The other field components follow from Eq. (4).

Solution (8a) contains only the regular Bessel function, because the other linearly independent solution—the Neumann function—is singular at the origin. Recently, it has been noticed that solution (8a) (including the phase factor) can actually be split into two counter-propagating conical waves each containing the Neumann function in its description.Correa2017 The net result, however, coincides with Eq. (8a). Solution (8b), however, could in principle read as a linear combination of the Hankel functions of the first and second kinds,

[TABLE]

We dismiss this possibility on the basis of two arguments. First, note that the leading asymptotic form ($\lvert\mkern 1.0mu{\rho}\mkern 1.0mu\rvert\gg 1$, $\lvert\mkern 1.0mu{\nu}\mkern 1.0mu\rvert\ll\lvert\mkern 1.0mu{\rho}\mkern 1.0mu\rvert$) of the Hankel functions is

[TABLE]

where the $+$ ($-$) sign in the exponent belongs to $H_{\nu}^{(1)}$ ($H_{\nu}^{(2)}$). If $\rho=\sqrt{k^{2}-h^{2}}$ is complex, then, as per our sign convention (see discussion following Eq. (5)), we must discard $H_{0}^{(2)}(\rho)$ because it diverges as $\operatorname{Im}{\rho}\to\infty$. If, otherwise, $\rho$ is real, an energy condition rules $H_{0}^{(2)}(\rho)$ out from the solution as follows. For an infinite wire sustaining a longitudinal EM wave that draws its energy from a source at $z=-\infty$, the energy flux across a cylindrical surface coaxial with the wire must vanish. This energy flux per unit length is given by $U_{r}(r,t)=2\pi rS_{r}(r,t)$, where $S_{r}$ is the radial component of the Poynting vector $\mathbf{S}=\mathbf{E}\times\mathbf{H}$, namely,

[TABLE]

Taking Eq. (10) into account and doing the algebra we obtain

[TABLE]

It is always healthy to check once in a while (particularly in lengthy calculations or derivations with many intermediate steps) if the physical dimensions of the expressions obtained are correct. For example, at this point it is a good idea to verify that $U_{r}(r,t)$ indeed has dimensions of energy flux per unit length, [$U_{r}$] $=$ MLT*-3* $=$ $($ML2T${}^{-2}/$T$)/$L (in SI units, J s*-1* m*-1*). Since $r/\rho$ remains finite (in fact, constant) at arbitrarily large $r$ and the phase $hz-\omega t$ can assume any value we conclude that, for real $\rho$, the condition $U_{r}(r\gg a,t)=0$ leads to $B_{1}=0=B_{2}$. Equations (8) are then justified.

Finally, fields $E_{z}^{0}$ and $H_{\theta}^{0}$ must be continuous at the surface ($r=a$) of the wire. This boundary condition leads to the transcendental equation

[TABLE]

where the primes indicate differentiation with respect to the given argument. Equation (13) determines the wave number $h$, from which most of the electromagnetic properties of the system can be derived. In what follows we solve Eq. (13) both numerically and by an approximate method due to Sommerfeld and compare the results.

III Exact numerical solutions of the boundary condition

Before embarking on the solution of Eq. (13), let us look at the magnitude of the physical quantities involved in a typical case.

Consider a microwave propagating in a hard-drawn copper wire of radius $a=1$ mm (approximately wire no. $12$ in the American Wire Gauge system, for which $2a=2.053$ mm). The following parameters are taken from Refs. AshcroftMermin1976, ; Wikipedia2018, . The dc electrical conductivity of copper at $293$ K ($\sim 20$ ∘C) is $5.96\times 10^{7}$ $\Omega^{-1}$ m-1, while for $f=1$ GHz we have $ka\simeq 0.021$, where $k=\omega/c=2\pi f/c=2\pi/\lambda$ is the free-space EM wave number and $\lambda=c/f\simeq 30$ cm its wavelength. Note that for most metals the electrical conductivity is essentially static up to about the infrared region $f\sim 10$ THz or, more meaningfully, photon energies about $E=hf\sim 40$ meV (where Planck’s constant $h\simeq 4.136\times 10^{-15}$ eV s);planck microwaves ($f\sim 1$–$100$ GHz, $\lambda\sim 300$–$3$ mm, $E\sim 0.004$–$0.4$ meV) are far below this energy scale. As we have already mentioned, for copper $\mu\simeq\mu_{0}$ (the more precise current value is $\mu_{\mathrm{Cu}}=0.999\,994\,\mu_{0}$), but the electric permittivity is a complex-valued function of the frequency, $\hat{\varepsilon}=\hat{\varepsilon}(\omega)$. For most metals, however, the conduction current $\mathbf{J}=\sigma\mathbf{E}$ is much larger than the displacement current $\mathbf{J}_{D}=\varepsilon\partial_{t}\mathbf{E}=-i\omega\varepsilon\mathbf{E}$, since $\varepsilon/\sigma\ll 1$. In fact, for most metals $\varepsilon/\sigma\sim 10^{-18}$, and as long as $f\lesssim 10^{16}$ Hz the displacement current can be ignored. This corresponds to the low-frequency limit of an ac Drude model for copper. We thus obtain from Eq. (5a) that $\hat{k}^{2}\approx i\sigma\mu\omega$, or

[TABLE]

The real quantity $\delta(\omega)=1/\kappa(\omega)$ measures the depth in the material at which the current density decreases to a fraction $e^{-1}\simeq 37\%$ of its value at the surface—the so-called skin depth. Note that $\delta(\omega)\sim 1/\sqrt{\omega}$ can become a really small quantity at high frequencies, sometimes just a few nanometers, with the current mostly confined to the periphery of the conductor. If we plug in the figures of our example we obtain $\kappa\simeq 4.851\times 10^{5}$ m*-1*, that is, $\delta\simeq 2~{}\mu$m.

To find the wave number $h$, we solve Eq. (13) numerically with the help of the software package MathematicaWolfram2018 using the figures given above for the physical quantities, as well as approximation (14) for $\hat{k}$. Graphing the two sides of Eq. (13) is unilluminating, so we just quote the solution,

[TABLE]

Two conclusions of physical significance can be immediately drawn from this number. First, that the wave is damped, since the phase factor $\exp(ihz)$ contains an attenuation factor given by $\exp(-7.516\times 10^{-4}z)$, and the amplitude of the wave decays by a factor $e^{-1}$ every $\sim 1330$ m. This attenuation can be considered small, and such a naked copper wire no. 12 at $1$ GHz indeed qualifies as a viable candidate for last mile data carrier! The second conclusion is that the wave travels along the wire at phase speed

[TABLE]

i. e., at about $0.999\,926\,c$, just a tad below the speed of light. The “surface wave” is thus able to travel without much damping for a reasonable distance—but not telegraphic distances, much to the chagrin of the telecommunications engineers at the turn of the XX century.

IV Sommerfeld’s approximate solution

Sommerfeld approached the solution of Eq. (13) by making well informed assumptions about the physical nature of the system at hand and then approximating the functions appearing there accordingly. He was then led to a much simpler transcendental equation that he solved by means of a practical and effective device. In what follows we review the rationale leading to the approximate equation for the boundary condition and Sommerfeld’s iterative solution to the problem.

IV.1 The physical approximations

For a metallic wire, $\sigma$ is large and $\hat{k}$ and $\hat{\rho}$ are large complex numbers. Given that $J_{0}^{\,\prime}(\hat{\rho})=-J_{1}(\hat{\rho})$ and that for large $\hat{\rho}$ in the complex upper half-plane (remember that we are using $\operatorname{Im}\hat{\rho}\geq 0$ by convention) it holds that

[TABLE]

the right-hand side of Eq. (13) can be simplified by taking $\sqrt{\hat{k}^{2}-h^{2}}/\hat{k}\approx 1$ and $J_{0}(\hat{\rho})/J_{0}^{\,\prime}(\hat{\rho})\approx-e^{-i\frac{\pi}{2}}=i$ (we always take the principal branch), such that

[TABLE]

In metals widely used for electric wiring like copper, silver or aluminium, $\lvert\mkern 1.0mu{\mu}\mkern 1.0mu\rvert\simeq\mu_{0}$ to a very good approximation but $\lvert\mkern 1.0mu{\hat{\varepsilon}}\mkern 1.0mu\rvert\gg\varepsilon_{0}$, and we see from Eq. (18) that the left-hand side of Eq. (13) must also be a small number, i. e., $\lvert\mkern 1.0mu{\!\sqrt{k^{2}-h^{2}}}\mkern 1.0mu\rvert$ must be small. The expansion of the Hankel functions about the origin gives

[TABLE]

where $\gamma=0.577\cdots$ is the Euler-Mascheroni constant and $\Gamma(\nu)$ is the usual gamma function. Since $H_{0}^{\,\prime}(\rho)=-H_{1}(\rho)$, the left-hand side of Eq. (13) becomes

[TABLE]

Now we put everything together. Defining the variables

[TABLE]

Eq. (13) for the wave number $h$ of a relatively low frequency pure TM mode propagating in a cylindrical metallic wire becomes, within the aforementioned approximations,

[TABLE]

This equation corresponds to equations ($22$)–($22^{\prime}$) in Sommerfeld’s original paper,Sommerfeld1899 equations ($9$)–($10$) in Goubau’s treatment of the problem,Goubau1950 equations (IX.$23$–$24$) in the textbook by Stratton,Stratton1941 and equation ($22.15$) in in the textbook by Sommerfeld.Sommerfeld1952 Note that while on the left-hand side of Eq. (22) we have the geometric factors of the problem ($a$, $k$, and $h$), the right-hand side concentrates the physical factors ($\hat{\varepsilon}/\varepsilon_{0}$ and $\mu/\mu_{0}$), providing a clear illustration of the role played by boundary conditions in electrodynamics.

IV.2 The iterative solution

Sommerfeld solved Eq. (22) by “a method of successive approximations [that] proved to be very suitable.”Sommerfeld1898 Several references repeat his solution, sometimes second-hand from Stratton,Stratton1941 while Goubau offers an interesting approximate graphical solution of his own.Goubau1950 In what follows we present Sommerfeld’s method according to more modern standards.

Let $f(x)=x\ln{x}$ be a real function of a real argument $x\geq 0$. Routine inspection reveals that $f(x)$ has an absolute minimum at $x_{\text{min}}=e^{-1}\simeq 0.367\,879$, at which $f(x_{\text{min}})=-e^{-1}$, and since $\lim_{x\to 0}f(x)=0$ and $f(1)=0$, it is decreasing for $0\leq x<e^{-1}$ and increasing for $x>e^{-1}$. We thus conclude that the equation

[TABLE]

does not have real solutions for $a<-e^{-1}$, has two solutions for $-e^{-1}\leq a<0$, one in the interval $(0,e^{-1})$ and the other in the interval $(e^{-1},1)$ (the two coinciding when $a=-e^{-1}$), and has a single solution for $a\geq 0$. These facts are summarized in Fig. 2.

IV.2.1 The smaller solution

Based on the previous analysis, Sommerfeld devised two iterative procedures to find the solutions of Eq. (23) depending on the location of the solution. Let us first consider the case $-e^{-1}<a<0$. The smaller solution of Eq. (23), say, $x_{l}^{*}$, is located to the left of $x_{\text{min}}=e^{-1}$. Sommerfeld then observes that if we write Eq. (23) as $x(-\ln{x})=-a$ and take logarithms we obtain

[TABLE]

and since $\ln{x}$ varies slowly in comparison with $x$ (an error of $\delta$ in the value of $x$ gives an error of roughly $\delta/x$ in the value of $\ln{x}$), we can ignore the second term on the right-hand side of Eq. (24) and take $x_{1}\simeq-a$ as a first approximation to $x_{l}^{*}$. Once we have this first approximation to $x_{l}^{*}$, we can substitute it on the right-hand side of Eq. (24) to obtain a new approximation, say, $x_{2}$, hopefully better than the former, given by

[TABLE]

Note that since $a<0$ and $\ln(-a)<-1$ we have $0<x_{2}<x_{1}$. Plugging $x_{2}$ in Eq. (24) leads to $x_{3}=a/\ln{x_{2}}$ with $0<x_{3}<x_{2}$. The iterative procedure is clear:

[TABLE]

Stratton states that “the convergence [of Eqs. (26)] can be demonstrated for complex as well as real values of $a$” (Ref. Stratton1941, , p. 529) but neither offers a proof nor refers to one. He was certainly quoting from Ref. Sommerfeld1899, . The mathematical work of SommerfeldSommerfeld1898 on the iterative solution of Eq. (22) or (23) is rarely mentioned. Sommerfeld proved that $\lim x_{n}=x_{l}^{*}$ from Eqs. (26) based on the ratio test and the fact that $0<x_{n+1}<x_{n}$ to conclude that $(x_{n})$ forms a Cauchy sequence—for every given real $\varepsilon>0$ one can find an integer $N_{\varepsilon}$ such that $\lvert\mkern 1.0mu{x_{n}-x_{m}}\mkern 1.0mu\rvert<\varepsilon$ for all $m,\,n\geq N_{\varepsilon}$. We shall adopt a more general approach.

From elementary mathematical analysis,Donsig2010 ; Shlomo2013 we know that the “discrete dynamical system” (26) has a fixed point $x^{*}=T(x^{*})$ if the map $T(x)=a/\ln{x}$ is a contraction, i. e., if it satisfies the so-called Lipschitz condition

[TABLE]

with $c<1$ for any two points $x,\,y$ in the domain of $T$. The constant $c$ is called the Lipschitz constant of the map. The nice thing about contractions is that a powerful theorem guarantees, under certain conditions, that they have a unique fixed point: the Banach contraction principle. Its statement reads:

A contraction $T\colon X\to X$ on a complete metric space $X$ has a unique fixed point $x^{*}$.

Remember that a complete metric space is a metric space in which every Cauchy sequence is convergent. An easy corollary of the Banach contraction principle is that $\lim_{n\to\infty}T^{(n)}(x)=x^{*}$ for all $x\in X$, where $T^{(n)}(x)$ denotes the $n$-fold composition $(T\circ\cdots\circ T)(x)$.Hille1972

We now revisit Eqs. (26) in the light of the above knowledge. First, note that the interval $X=[0,e^{-1}]$ is a complete metric space (in the usual Euclidean metric $\|y-x\|=\sqrt{(y-x)^{2}}=\lvert\mkern 1.0mu{y-x}\mkern 1.0mu\rvert$) because it is a closed subset of $\mathbb{R}$, which is complete. To check whether $T(x)=a/\ln{x}$ is a contraction on $X$ we take

[TABLE]

The largest possible value of the right-hand side of Eq. (28) occurs when one of its terms becomes arbitrarily small (i. e., when $x$ or $y\to 0$) and the other becomes $1/\ln(e^{-1})=-1$, such that $\lvert\mkern 1.0mu{T(y)-T(x)}\mkern 1.0mu\rvert\leq\lvert\mkern 1.0mu{a}\mkern 1.0mu\rvert$ in $X$. Putting all together we get

[TABLE]

It thus suffices to have $\lvert\mkern 1.0mu{a}\mkern 1.0mu\rvert<e^{-1}$ to make $T$ a contraction. Since $-e^{-1}<a<0$, we conclude that $T(x)$ is indeed a contraction on $X=[0,e^{-1}]$ and that its unique fixed point $x^{*}=T(x^{*})$ is, by construction, the solution of Eq. (23) on $X$.

IV.2.2 The larger solution

If we try to find the second solution of Eq. (23) in the interval $e^{-1}<x<1$ with $-e^{-1}<a<0$ using the recursion relation in Eqs. (26), it may or may not converge and, if it converges, it will do so to the smaller solution. From Eq. (23), however, we see that there are two options to draw an iterative procedure of it: either $x_{n+1}\ln{x_{n}}=a$, corresponding to Eqs. (26), or $x_{n}\ln{x_{n+1}}=a$. Sommerfeld realized that the second option provides a scheme to find the second solution of Eq. (23) in the interval $e^{-1}<x<1$ when $-e^{-1}<a<0$. The recursion relation reads

[TABLE]

We employ again the contraction principle to prove that this recursion relation indeed converges to the right solution.

Let $X=[e^{-1},1]$ and $T\colon X\to X$ be the map $T(x)=\exp(a/x)$. We want to check whether

[TABLE]

with $c<1$. The function $c(a,x,y)=\lvert\mkern 1.0mu{e^{a/y}-e^{a/x}}\mkern 1.0mu\rvert/\lvert\mkern 1.0mu{y-x}\mkern 1.0mu\rvert$ is decreasing in all of its arguments, so its maximum should be at the point $Q=(e^{-1},e^{-1},-e^{-1})$. Unfortunately, $c(a,x,y)$ is undefined for $y=x$. However, since $c(a,x,y)=c(a,y,x)$, we argue that approaching $Q$ along the $x$-axis or the $y$-axis should give the same result, making $\lim_{P\to Q}c(P)$ well defined (if it exists at all). We then fix $x$ and let $y\to x$ to obtain

[TABLE]

where we resolved the indeterminate form $0/0$ by L’Hôpital’s rule. Setting $x=e^{-1}$ and $a=-e^{-1}$ we find that $c(a,x,y)\leq 1$, with strict inequality for $a>-e^{-1}$. We thus conclude that $T(x)$ satisfies the contraction principle on $X$ for $-e^{-1}<a<0$ such that its unique fixed point on $X$ indeed solves Eq. (23).

It is useful to examine an example. Let $a=-0.2$. Application of Eqs. (26) provides the sequence $x_{1}=0.2\to x_{2}=0.124267$ $\to\cdots\to x_{20}=0.0786584\to x_{21}=0.0786584$, converging to seven decimal places in $20$ iterations, with $x_{20}\ln{x_{20}}\simeq-0.200\,000\,061$. The solution is $x_{l}^{*}=0.078\,658\,360\cdots$. Now let us apply Eqs. (30) to find the righmost solution of Eq. (23). We obtain $x_{1}=0.8187308\to x_{2}=0.7832679\to\cdots\to x_{11}=0.7716910\to x_{12}=0.7716910$, and the procedure converges to seven decimal places in just $11$ iterations, with $x_{11}\ln{x_{11}}\simeq-0.199\,999\,981$. The precise value for the rightmost solution is $x_{r}^{*}=0.771\,690\,974\cdots$. The sequence converges faster because the Lipschitz constant for Eqs. (30) is smaller than the one for Eqs. (26), bringing separate points close together faster.

IV.2.3 Cases $a>0$

When $a>0$, Sommerfeld splits the resolution of Eq. (23) again into two cases. For $0<a<e$, the proper convergent iterative scheme reads the same as Eqs. (30), namely,

[TABLE]

while for $a\geq e$ the scheme that converges coincides with Eqs. (26) except that now we start with $x_{1}=a$, to wit,

[TABLE]

We invite the reader to repeat the analysis of the above maps along the lines of what we did for the other cases to prove their convergence, as well as to examine some examples. We will not make use of these procedures in our study.

In his mathematical memoir,Sommerfeld1898 Sommerfeld attempts a generalization of his analyses to an abstract equation of the form

[TABLE]

under mild hypotheses on the continuity of $f$, $\phi$, and their inverse functions. Barring possible factorizations or symmetrizations, a basic iterative scheme to solve Eq. (35) either takes $f(x_{n+1})=a/\phi(x_{n})$ or $\phi(x_{n+1})=a/f(x_{n})$. Analysis of which scheme provides a contraction then suggests the procedure to follow.

IV.3 Solution of the approximate boundary condition

We now apply what we have learned to the solution of Eq. (22). We have a problem, though: equation Eq. (22) involves complex quantities. The Banach contraction principle, however, only requires that the underlying space be complete in some norm, and this is the case of the complex plane—in fact, the complex plane $\mathbb{C}$ is isometric with the real plane $\mathbb{R}^{2}$, with the usual complex norm $\sqrt{z^{*}z}$ in the former corresponding naturally to the Euclidean norm $\sqrt{x^{2}+y^{2}}$ in the latter. Of course, before employing Eqs. (26) or Eqs. (30) one must analyse the moduli and the arguments of $x$ and $a$ to make a better informed trial.

Using the approximations $\mu\simeq\mu_{0}$ and $\varepsilon\omega\ll\sigma$ valid for good metal conductors at low-frequencies (see the discussion preceding Eq. (14)) we obtain for $v$ (Eq. (21))

[TABLE]

and inserting the values of Sec. III for the 1 mm radius copper wire at $1$ GHz we get

[TABLE]

a number with a small negative real part. We shall thus try to apply procedure Eqs. (26) to find the solution to Eq. (22). We obtain, to three decimal places, the following sequence of numbers, where $z_{-e}$ reads $z\times 10^{-e}$: $u_{1}=(7.181+7.181\,i)_{-7}$ $\to u_{2}=(4.892+5.482\,i)_{-8}$ $\to\cdots$ $\to u_{6}=(4.095+4.529\,i)_{-8}$ $\to u_{7}=(4.095+4.529\,i)_{-8}$, and the procedure converged after just $6$ iterations to the solution

[TABLE]

This solution is slightly at variance with that found by Sommerfeld because he uses $v=-7.2\times 10^{-7}(1+i)$ and also because we keep more digits in the intermediate calculations.

From the value for $u$ above we obtain from Eq. (21) that

[TABLE]

Compare this number with the one in Eq. (15). The real parts are virtually identical (the difference starts from the fourth decimal place, not shown), but the imaginary parts differ almost by a factor of two. The root of this difference is likely the somewhat uncontrolled approximation $\sqrt{\hat{k}^{2}-h^{2}}/\hat{k}\approx 1$ that we made in passing from Eq. (13) to Eq. (18); the approximation is valid in absolute value, but not componentwise. Following the analysis presented after Eq. (15), we conclude that the wave is damped by a factor $e^{-1}$ every $\sim 734$ m and that

[TABLE]

i. e., travels along the wire at a phase speed about $0.999\,941\,c$. While the physical conclusions in this case remain the same, namely, that the wave propagates reasonably undamped almost at the speed of light, the attenuation factor appears larger than it really is as a result of the approximations made.

In his textbook,Sommerfeld1952 Sommerfeld compares the above case with that of a very thin platinum wire of radius $a=2\times 10^{-6}$ m forced at $f=300$ MHz. The conductivity of platinum at $293$ K is $\sigma=9.52\times 10^{6}$ $\Omega^{-1}$ m*-1*, approximately $6.3$ times smaller than that of copper.Wikipedia2018 Note that now $\hat{\rho}(a)=\sqrt{\hat{k}^{2}-h^{2}}\,a$ is a small quantity and the approximation Eq. (18) no longer holds—one has to expand the Bessel functions around $\hat{\rho}=0$ instead. The reader is invited to perform the analysis of the EM wave propagation characteristics in this case to conclude that the wire behaves almost like a purely resistive element, something that could be guessed already by a simple estimate of the skin depth $\delta(\omega)\sim\sqrt{2/\sigma\mu_{0}\omega}$.

V The solution in terms of the

Lambert $W$ function

The Lambert $W$ function has been around since at least the XVIII century, when it was introduced in disguise by Lambert in 1758 in the series solution of the trinomial equation $x=q+x^{m}$, a solution which Euler revisited and expanded in 1776. As applied mathematics and mathematical physics became more diverse in the XIX and XX centuries, the Lambert $W$ function started to pop up in several areas, from statistical mechanics and astrophysics to population dynamics, materials science, combinatorics, and the analysis of algorithms, among others.LambertW1993 ; LambertW1996 ; Valluri2000 ; Caillol2003 ; Cranmer2004 ; WarbWang2004 ; Packel2004 ; Hayes2005 ; Scott2006 ; Veberic2010 The fact that numerical routines to evaluate the Lambert $W$ function started to appear in the early 1970s also attests its ubiquitous occurrence in many areas of science and technology.Wright1959 ; Fritsch1973 ; Barry1995 ; Barry743 ; Boyd1998 ; Barry2000 By the mid-1990s, notation—in particular the use of the letter $W$ for the function after its first implementation in the Maple computer algebra software—and basic properties were settled, with Ref. LambertW1996, having become the standard reference on the Lambert $W$ function.

The Lambert $W$ function is defined to be the inverse of the map $z\mapsto ze^{z}$, with $z$ a complex number. Put equivalently, the values of $W(z)$ are given by the solutions of the equation

[TABLE]

Since the map $z\mapsto ze^{z}$ is not injective, the solutions of Eq. (41) are multivalued. The map cannot be inverted easily—in fact, this is the whole point of defining the “special function” $W(z)$. A simple strategy to plot $W(x)$ for real $x$ is to plot $y=xe^{x}$ and interchange the roles of axis $x$ and $y$: look at the plot in a mirror or from behind and rotate it $\pi/2$ clockwise. A plot of $W(z)$ for real $z$ is displayed in Fig. 3. From this picture it is clear that the point $z=-e^{-1}$, at which $W(z)=-1$, is a branching point. The two real branches of $W(z)$ are the principal branch $W_{0}(z)\geq-1$, defined for all $z\geq-e^{-1}$, and $W_{-1}(z)<-1$, defined for $-e^{-1}\leq z<0$ with a singularity $\lim_{z\to 0^{-}}W_{-1}(z)=-\infty$. If $z<-e^{-1}$ or $z$ is complex, then $W(z)$ becomes complex. The asymptotic behavior of $W_{0}(z)$ for small $z$ is given by

[TABLE]

while for large $z$ it is given by

[TABLE]

Back to our mathematical physics problem, if we put $z=\ln{u}$, Eq. (22) becomes $ze^{z}=v$, which immediately furnishes the solution $u=e^{W(v)}$, or, equivalently,

[TABLE]

Unbeknownst to Sommerfeld, the solution of Eq. (22) is “simply” given in terms of the Lambert $W$ function. This is not really helpful unless one has a procedure to compute $W(v)$—exactly what Sommerfeld devised with his iterative method. By comparing Figs. 2 and 3 and noting that the two are connected by the relationship $z=\ln{x}$ one can easily conclude that the leftmost solution of Eq. (23) corresponds to $\exp[W_{-1}(a)]$, while the righmost solution corresponds to $\exp[W_{0}(a)]$. In the software package Mathematica,Wolfram2018 one can calculate $W_{0}(z)$ with the command ProductLog[z], while the calculation of $W_{-1}(z)$ is achieved with the command ProductLog[1,z]. When we insert $v=-7.181\times 10^{-7}(1+i)$ (see Eq. (37)) into Eq. (44) we obtain, using Mathematica, the solution

[TABLE]

which is exactly the same value as before, Eq. (38).

VI Summary and conclusions

We have shown that a certain approximate equation arising from the boundary condition for EM waves between a cylindrical metallic wire and a dielectric can be solved by means of the Lambert $W$ function. Although one can solve the exact condition Eq. (13) numerically without much fuss, the physical rationale leading to the approximate condition Eq. (22) is very illuminating and, as we have seen, provides a reasonably accurate description of the system. In fact, it is quite satisfying to verify that the answer obtained by Sommerfeld through his apparently offhand iterative method and the solution in terms of the Lambert $W$ function agree. While Sommerfeld’s solution and the numerical solution of the exact boundary condition Eq. (13) do not quite match, they are qualitatively compatible and predict the same physics. We can trace the somewhat uncontrolled approximation $\sqrt{\hat{k}^{2}-h^{2}}/\hat{k}\approx 1$ as the probable reason for the mismatch between the exact and the approximate solutions.

Equation (41) is the simplest possible exponential polynomial equation, and this is arguably one of the reasons for its ubiquity. Slightly generalized equations, however, sometimes appear which can also be solved in terms of $W(z)$.Gosper1998 For example, in the study of the capture of a diffusing prey by diffusing predators in one dimension, the authors in Ref. Sidney1999, obtain the equation $ye^{y^{2}}=M$, where $y=x_{\mathrm{last}}/\sqrt{4Dt}$ is the rescaled position of the last predator (the record position of the predators’ one-dimensional random walks) and solve it iteratively as $y=\sqrt{\ln{M}}(1-\frac{1}{4}\ln\ln{M}/\ln{M}+\cdots)$. Clearly, squaring the original equation and then multiplying it by $2$ immediately leads to the solution $y^{2}=\frac{1}{2}W_{0}(2M^{2})$, from which good asymptotics, explicit numbers, and graphs can be obtained with a couple of commands in more or less any modern computer algebra system or scientific library.

The Lambert $W$ function can be computed to arbitrary precision by iterative root-finding of $we^{w}=z$, with the usual trade off between complexity of implementation and number of iterations to converge at the required precision. Newton’s simple first-order method, for example, is appropriate, but converges too slowly for modern standards. The Maple software package implements Halley’s third-order variant of Newton’s method that attains high precision in affordable time.LambertW1993 ; Veberic2018 The investigation of the recursion relations (26) and (30) on the complex plane in the language of dynamical systems is not entirely without interest and may even prove useful in the implementation of numerical routines to calculate $W(z)$. An accessible and comprehensive account of discrete dynamical systems based on the ideas of inverse functions (“root-finding”) and fixed point theorems is given in Ref. Shlomo2013, .

In the introduction to Volume III of his Lectures on Theoretical Physics, Sommerfeld states that “Heinrich Hertz’s greatest paper on the ‘Fundamental Equations of Electrodynamics for Bodies at Rest’ has served as model for my lectures on electrodynamics since my student days.”Sommerfeld1952 This is a remarkable note and advice, namely, that even the masters acknowledge that there is much to be learned from the masters. Maybe it is because they were closer to sparking gaps, vacuum tubes, flaring chemicals, and gears than the average contemporary student will ever be. There is a little bit—usually much more—for everyone. We encourage the reader to browse through Refs. Stratton1941, ; Sommerfeld1952, ; SommerfeldVI, looking for historical perspective, foundations, applications, and classical problems to solve.

Acknowledgements.

The author thanks the LPTMS for kind hospitality during a sabbatical leave in France and FAPESP (Brazil) for partial support through grant no. 2017/22166-9.