Introduction

Surface waves (SW) on circular conducting wires have been of theoretical interest since their discovery by Sommerfeld1. However, due to the large lateral extent of the fields at low frequencies, practical applicability seemed limited at first. Goubau discovered that by coating wires in dielectric or corrugating the wire’s surface, the fields’ lateral confinement could be drastically enhanced2. These coated wires, named Goubau lines in later years, showed low loss and weak dispersion. Hence, they were discussed as an alternative to traditional, two-conductor transmission lines. Recently, interest in surface wave technology has re-emerged at the GHz-THz frequency range where it presents a promising alternative to current waveguide technology3,4. Furthermore, by introducing sub-wavelength corrugations, spoof surface plasmon polaritons emerge which have tunable properties and can exhibit sub-wavelength lateral confinement5,6,7. These new technologies have been discussed as solutions to problems such as signal integrity in integrated circuits and backhaul solutions for the network standard 5G8,9.

In many cases it is favourable to print a conductor design on a substrate using established printed circuit board fabrication processes such as etching. This led to the invention of a planar Goubau line (PGL) consisting of a thin rectangular conducting strip on a substrate10,11,12,13. It represents an alternative to standard substrate transmission lines namely the microstrip line and as such multiple essential electronic components have been devised for the PGL including broadband loads, power dividers and frequency selective filters14,15,16. Substrate integrated leaky wave antennas have been designed based on the PGL as well17,18,19. Additionally, application of PGLs in terahertz spectroscopy has been established20,21.

Despite these advances, only numerical and experimental studies have been published on the PGL to date and no analytic theory or model exists22,23. A difficulty in the exact treatment is presented by the presence of sharp corners which introduce lightning rod effects12. Therefore, a simplified model of the Goubau line in which the rectangular conductor is exchanged for a conductor with circular cross section will be considered in this paper. For ease of notation the simplified system will be referred to as a PGL as well. A related system - the single conductor above a semi-infinite conducting earth - has been extensively studied (see24 and references therein). We will draw parallels to this system where appropriate.

The paper is structured as follows: First, we discuss the wave created by an infinitesimally small current filament above a substrate. Then, the finite thickness of wire is incorporated to derive a characteristic equation for the system and applicability criteria for this approach are discussed. Using the derived equation, the dependence of the propagation constant on input parameters is established and some field patterns are reported. We validate our results against numerical data obtained through the finite element method and experiment. Finally, we draw parallels to the Goubau line with rectangular conductor.

Derivation of the electromagnetic fields

Figure 1
figure 1

Cross sectional illustration of the simplified planar Goubau line. A wire of radius a is at a height h over a dielectric slab with permittivity \(\epsilon _2\) and thickness b. The surrounding medium has permittivity \(\epsilon _1\). A Cartesian coordinate system (x, y, z) is defined with its origin placed on the substrate in line with the centre of the conductor.

Full size image

The system we are investigating consists of a perfectly conducting wire of radius a with its centre located at a height h above a substrate of thickness b. Figure 1 shows the cross section of the system and defines the coordinate system (x, y, z) located on the surface of the substrate with its origin in line with the centre of the wire. We assume the wire and substrate are uniform in the z-direction which will be the direction of wave propagation. Additionally, the substrate extends infinitely far in the x-plane. It has a dielectric constant of \(\epsilon _2\) and is immersed in a medium with dielectric constant \(\epsilon _1\) with \(\epsilon _1<\epsilon _2\). For most practical applications the surrounding medium is air whose dielectric constant can be approximated as the dielectric permittivity of vacuum \(\epsilon _0\). All materials are assumed to be non-magnetic and have permeability equal to the magnetic constant \(\mu _0\).

We start our derivation by postulating a time harmonic current density \(\vec {J}\) that replaces the wire. Its form, which is motivated by the current found in regular surface waveguides such as the Goubau line, is given by:

$$\begin{aligned} \vec {J}=I\delta (x)\delta (y-h)e^{i\omega t - i\beta z}{\hat{z}} \end{aligned}$$
(1)

with angular frequency \(\omega\), time t, propagation constant \(\beta\), current amplitude I, unit vector \({\hat{z}}\) and Dirac delta function \(\delta\). This approach is only strictly valid in the case where the wire is infinitesimally small but it can also give reasonable results for thin wires. The physical conditions under which wires may be considered thin will be given later when discussing the characteristic equation. The presented current density can be interpreted as the exact current distribution averaged across the cross section of the wire. For now, we note that as we are not resolving the exact current distribution inside the wire, the near field close to the wire will deviate from an exact solution. However, at distances much greater than the wire radius, the distribution of current inside the wire should have an insignificant effect on the electromagnetic fields.

Next, we assume that the total field can be separated into a transverse magnetic (TM) and transverse electric (TE) component which are characterised by having no longitudinal magnetic or electric field, respectively. Thus, we may express the total electric and magnetic field, \(\vec {E}\) and \(\vec {H}\), as

$$\begin{aligned} \vec {E}=\vec {E}^{TM}+\vec {E}^{TE} \quad &\vec {H}=\vec {H}^{TM}+\vec {H}^{TE} \end{aligned}$$
(2)

We expect the fields to be in phase with the current density and should also contain a factor \(e^{i\omega t - i \beta z}\) which will be implicitly assumed but omitted for clarity. Ampere’s and Faraday’s law then take the form

$$\begin{aligned} \vec {\nabla } \times \vec {H}^{TM}&=-i\omega \epsilon _j\vec {E}^{TM}+I\delta (x)\delta (y-h){\hat{z}}&\vec {\nabla } \times \vec {E}^{TM}&=i \omega \mu _0 \vec {H}^{TM} \end{aligned}$$
(3)
$$\begin{aligned} \vec {\nabla } \times \vec {H}^{TE}&=-i\omega \epsilon _j \vec {E}^{TE}&\vec {\nabla } \times \vec {E}^{TE}&=i \omega \mu _0 \vec {H}^{TE} \end{aligned}$$
(4)

where \(\epsilon _j\) takes either the value \(\epsilon _1\) or \(\epsilon _2\). From the defining property of TE and TM modes, it may be shown that all field components may be calculated from the longitudinal component of the magnetic and electric field, \(H_z\) and \(E_z\), respectively25. Taking the curl of the second equation in (3) and using Faraday’s and Gauss’ law in combination with the continuity equation as well as standard vector calculus identities, we arrive at

$$\begin{aligned} \Big (\frac{\partial ^2}{\partial x^2}+\frac{\partial ^2}{\partial y^2}+\gamma _j^2\Big )E_z^{TM}=-i\omega \mu _0I\frac{\gamma _j^2}{k_j^2}\delta (x) \delta (y-h) \end{aligned}$$
(5)

where we have introduced \(k_j^2=\omega ^2 \epsilon _j \mu _0\) and the new variable \(\gamma _j^2=k_j^2-\beta ^2\) with \(j\in \{1,2\}\). In a similar fashion, we may manipulate Eq. (4) to arrive at

$$\begin{aligned} \Big (\frac{\partial ^2}{\partial x^2}+\frac{\partial ^2}{\partial y^2}+\gamma _j^2\Big )H_z^{TE}=0 \end{aligned}$$
(6)

As the boundaries along the dielectric substrate extend infinitely along the x-direction, it is convenient to introduce a Fourier transform and its inverse as

$$\begin{aligned} E^{TM}_z=\frac{1}{2\pi }\int _{-\infty }^{\infty }{\tilde{E}}^{TM}_z e^{-i \xi x}d\xi \quad &{\tilde{E}}^{TM}_z=\int _{-\infty }^{\infty }E_z^{TM} e^{i \xi x} dx \end{aligned}$$
(7)

where \(\xi\) represents a spatial frequency. Similarly, \({\tilde{H}}_z^{TE}\) is the Fourier transform of \(H_z^{TE}\). Generally, we signify functions in Fourier space by a tilde. The transformed Eqs. (5) and (6) become:

$$\begin{aligned}&\Big (\frac{\partial ^2}{\partial y^2}+u_j^2\Big ){\tilde{E}}_z^{TM}=-i\omega \mu _0I\frac{\gamma _j^2}{k_j^2} \delta (y-h) \end{aligned}$$
(8)
$$\begin{aligned}&\Big (\frac{\partial ^2}{\partial y^2}+u_j^2\Big ){\tilde{H}}_z^{TE}=0 \end{aligned}$$
(9)

with \(u_j^2=\gamma _j^2-\xi ^2\). Solutions to these equations are readily available26. We impose the condition that fields should decay towards infinity and find

$$\begin{aligned} {\tilde{E}}^{TM}_z={\left\{ \begin{array}{ll} -\omega \mu _0 I \frac{\gamma _1^2}{k_1^2} \frac{e^{i u_1 |y-h|}}{2u_1}+C_1 e^{iu_1y} &\qquad y\ge 0\\ C_2 e^{i u_2 y}+ C_3 e^{-i u_2 y} &\quad-b\le y \le 0 \\ C_4 e^{-i u_1 y} &\quad y\le -b \end{array}\right. } \end{aligned}$$
(10)
$$\begin{aligned} {\tilde{H}}^{TE}_z={\left\{ \begin{array}{ll} \frac{C_5}{\eta _0} e^{iu_1y} &\quad y\ge 0\\ \frac{C_6}{\eta _0}e^{i u_2 y}+ \frac{C_7}{\eta _0} e^{-i u_2 y} &\quad -b\le y \le 0 \\ \frac{C_8}{\eta _0} e^{-i u_1 y} &\quad y\le -b \end{array}\right. }. \end{aligned}$$
(11)

where \(u_1\) is defined to have positive imaginary part such that the fields’ energy remains finite. \(C_1\) to \(C_8\) are constants yet to be determined. We introduced the free space impedance \(\eta _0=\sqrt{\mu _0/\epsilon _0}\) so all constants have the same dimensions.

At the interface between the dielectric substrate and air, the tangential components of \(\vec {E}\) and \(\vec {H}\) must be continuous. Introducing \(A=-\omega \mu _0 I \frac{\gamma _1^2}{k_1^2}\frac{e^{i u_1 h}}{2u_1}\), these boundary conditions may be expressed in the following matrix equation

$$\begin{aligned} \underline{\underline{Q}}\begin{pmatrix} C_1 \\ C_2\\ C_3\\ C_4\\ C_5\\ C_6\\ C_7\\ C_8 \end{pmatrix}= \begin{pmatrix} A\\ 0\\ 0\\ 0\\ -\frac{k_1 \sqrt{\varepsilon _{r,1}}u_1}{\gamma _1^2}A\\ 0\\ \frac{\beta \xi }{\gamma _1^2}A\\ 0 \end{pmatrix} \end{aligned}$$
(12)

with relative dielectric permittivities \(\varepsilon _{r,j}=\epsilon _j/\epsilon _0\) and matrix \(\underline{\underline{Q}}\) which is given in the supplementary material. \(\underline{\underline{Q}}\) may be inverted to find expressions for the constants \(C_1\) to \(C_8\). We note here that the TE and TM modes are coupled through the continuity of \(H_x\) and \(E_x\) at the boundaries between substrate and air. Pure TM solutions with \(C_5, C_6, C_7\) and \(C_8\) equal to zero cannot fulfil the matrix equation. Hence, the resultant electromagnetic field will be hybrid in nature contrary to the Goubau mode for the dielectric coated cylinder.

Characteristic equation

So far we have discussed the exact solution for an arbitrary, infinitesimal current filament carrying a known current wave at a height h above a substrate. However, in most practical cases, we want to find the propagation constant of the wave carried by an extended wire of finite size. This is still a formidable task even with the possibility of formally expressing the electromagnetic fields given any current distribution by convoluting the calculated fundamental solution with a given source term26.

As a means of characterising the propagating mode, we introduce the effective refractive index \(n_{eff}\) which is related to the free space wave number \(k_0\) via \(\beta =n_{eff}k_0\). In order to formulate an approximate characteristic equation, we assume the wire is a perfect electrical conductor (PEC). This is valid for many metals in the GHz to THz frequency range if the radius is much larger than the skin depth. As a PEC, the tangential electric field should be zero at its surface. In particular, the z-component of the electric field must be zero. Imposing this condition at any point on the wire’s surface, gives an equation for the approximate propagation constant if the wire is sufficiently thin27. Hence, the characteristic equation may be expressed as

$$\begin{aligned} E_z(x=a, y=h)=0 \end{aligned}$$
(13)

In a study on the validity of this approach for a wire in air above a semi-infinite earth, Pogorzelski and Chang showed that reasonable results are obtained if the contribution of azimuthally varying currents in the wire can be neglected28. Furthermore, it was shown in the same work that when the wire is close to the earth signified by \(k_2 h \ll 1\) and \(|\gamma _1| h\ll 1\), the contribution due to the first order azimuthal terms in the effective refractive index of the wave scales as

$$\begin{aligned} \frac{\Delta n_{eff}}{n_{eff}^{(0)}}=\frac{g(n_{eff})}{n_{eff}^{(0)}\partial E_z/\partial n_{eff}|_{x^2+(y-h)^2=a^2}}\Big |_{n_{eff}=n_{eff}^{(0)}} \end{aligned}$$
(14)

where \(\Delta n_{eff}\) is the correction due to higher order terms, \(n_{eff}^{(0)}\) is the zero order effective index and \(g(n_{eff})\) is a function of the effective refractive index and the geometry given by

$$\begin{aligned} g(n_{eff})=n_{eff}^2\Big (\frac{a}{2h}\Big )^2\Big (\frac{\epsilon _2-\epsilon _0}{\epsilon _2+\epsilon _0}\Big )\Big [\Big (\frac{a}{2h}\Big )^2+\frac{1}{1-2i\frac{ \epsilon _2-n_{eff}^2\epsilon _0}{(1-n_{eff}^2)(\epsilon _2+\epsilon _0)}}\Big ]^{-1} \end{aligned}$$
(15)

Hence, the correction to the effective refractive index \(\Delta n_{eff}\) due to azimuthal currents can be neglected if the absolute value of the right hand side in Eq. (14) is small. In our case the higher order contributions should be even smaller because the dielectric is only of finite thickness. Therefore, Eq. (14) provides an estimate on the obtainable accuracy when using the thin wire approximation to calculate the propagation constant.

In order to express and solve the characteristic equation, we focus on the amplitude \(C_1\) which is required for calculating \(E_z\) in real space, in contrast to Fourier space, above the substrate via the inverse Fourier transform (7). Solving the matrix Eq. (12), we find that it may be written in the form

$$\begin{aligned} C_1=A(-1+F(\xi )) \end{aligned}$$
(16)

where F is a function of \(\xi\) whose complete form is given in the supplementary material. It has some noteworthy properties. First it only depends on \(\xi ^2\) reflecting the mirror symmetry of the system with respect to the plane \(x=0\). Furthermore, for \(|\xi |\gg |\gamma _2|\) and \(|\xi b| \gg 1\) it asymptotically behaves as

$$\begin{aligned} F(\xi )\sim \frac{2k_1^2}{\gamma _1^2}\frac{\xi ^2 + u_1u_2}{k_1^2u_2 + k_2^2u_1}u_1 \sim \frac{k_1^2}{\gamma _1^2}\Big (1-\frac{2\beta ^2}{k_1^2+k_2^2}\Big ) \end{aligned}$$
(17)

which can be shown as all exponential terms \(e^{iu_2b}\) in \(F(\xi )\) will be very small. This expression is identical to the one obtained by Wait for the case of a wire above a semi-infinite earth27. In fact, in the corresponding limit of \(b\rightarrow \infty\) the conditions on \(\xi\) can be relaxed. This is due to all contributions involving \(e^{iu_2b}\) becoming infinitely fast oscillating or zero so that they can be neglected in the inverse Fourier transform in Eq. (7).

Let us reiterate the Fourier transform of \(E_z\), which is of the following form

$$\begin{aligned} \tilde{E_z}=-\omega \mu _0I\frac{\gamma _1^2}{k_1^2}\Big (\frac{e^{iu_1|y-h|}}{2u_1}-\frac{e^{iu_1(y+h)}}{2u_1}+\frac{e^{i u_1 (y+h)}}{2u_1}F(\xi )\Big ) \end{aligned}$$
(18)

To calculate the inverse Fourier transform, we use the identity

$$\begin{aligned} \int _{-\infty }^{\infty }\frac{e^{iu_1|(y\pm h)|}}{2u_1}e^{-i\xi x}d\xi =-iK_0\big (-i\gamma _1 \sqrt{x^2+(y\pm h)^2}\big ) \end{aligned}$$
(19)

where \(K_0\) is the zeroth order modified Bessel function of the second kind29. Furthermore, for a real number V for which \(V b \gg 1\) and \(V\gg |\gamma _2|\) we may replace \(u_1\) with \(i\xi\) and \(F(\xi )\) with its asymptotic form. Hence, the inverse transform gives

$$\begin{aligned} \int _{V}^{\infty }F(\xi )\frac{e^{i u_1 (y+h)-i \xi x}}{2u_1} d\xi \sim \frac{-ik_1^2}{2\gamma _1^2}\Big (1-\frac{2\beta ^2}{k_1^2+k_2^2}\Big )\Gamma \big (0,(ix+y+h)V\big ) \end{aligned}$$
(20)

where \(\Gamma\) is the incomplete Gamma function30. A similar expression is obtained for the integration from \(-\infty\) to \(-V\) as \(F(\xi )\) only depends on \(\xi ^2\). In the region from -V to V, no analytic expression for the integral was found but it may be calculated numerically. Note that, in general, the Cauchy principal value of the integral needs to be taken because the integrand may contain poles.

The field may then be expressed as

$$\begin{aligned} E_z= & {} -i\omega \mu _0 I \frac{\gamma _1^2}{2\pi k_1^2}\Big [K_0\big (-i\gamma _1\sqrt{x^2+(y+h)^2}\big )-K_0\big (-i\gamma _1\sqrt{x^2+(y-h)^2}\big )\nonumber \\&-\frac{k_1^2}{\gamma _1^2}\big (1-\frac{2\beta ^2}{k_1^2+k_2^2}\big )\mathfrak {R}\{\Gamma (0,(ix+y+h)V)\}-2i\int _{0}^{V}\cos (\xi x)F(\xi )\frac{e^{iu_1(y+h)}}{2u_1}d\xi \Big ] \end{aligned}$$
(21)

where \(\mathfrak {R}\) indicates that the real part of the expression in brackets should be taken.

Finally, let us consider the range of values which the free parameters of the problem, namely angular frequency \(\omega\), wire radius a, height h, substrate thickness b, current amplitude I and subtrate dielectric constant \(\varepsilon _2\) can take. The frequency \(\omega\) for which the model can be applied is capped by the possibility to carry out the numerical integration necessary for determining the longitudinal electric field in Eq. (21). We found that calculations above 30 GHz proved increasingly difficult. Additionally, higher order modes and azimuthally varying currents will play an increasing role with frequency deteriorating the accuracy of the model. The conductor radius is limited in size by the thin-wire approximation (14). As a rule of thumb \(a\ll \lambda\) with \(\lambda\) being the free space wavelength. Clearly, this also limits the height which cannot be smaller than the conductor radius. Similarly, for the substrate permittivity Eq. (14) provides a limit to the applicability. Generally, \(\varepsilon _2\) and \(\varepsilon _1\) should be of the same order of magnitude. The larger their difference, the stronger the azimuthally dependent currents are which make the theory inaccurate. The substrate thickness b will influence this as well and therefore should be small compared to the wavelength \(b \ll \lambda\). Especially, for \(b \sim \lambda\) or \(b\gg \lambda\) the theories in24 may be more applicable. Lastly, the theory is independent of the current amplitude I as long as non-linear effects can be neglected.

Solutions to the characteristic equation

Figure 2
figure 2

Effective refractive index \(n_{eff}\) of the propagating mode as a function of various parameters: (a) height, (b) frequency, (c) substrate permittivity and (d) substrate thickness. Results were obtained with the presented model and the finite element solver COMSOL Multiphysics. While sweeping one parameter, all other parameters were kept at their nominal values \(a=0.1\) mm, \(b=1.6\) mm, \(\varepsilon _{r,2}=3\), \(h=0.1\) mm and \(f=10\) GHz.

Full size image

In general, the characteristic Eq. (13) must be solved numerically due to the integral containing \(F(\xi )\). However, if the wire is located very far above the substrate such that \(|\gamma _1 h| \gg 1\) , the integral may be neglected due to the strong exponential damping. In fact, the characteristic equation is then dominated by a single term

$$\begin{aligned} K_0(-i \gamma _1 a)=0 \end{aligned}$$
(22)

This is the characteristic equation for a surface wave on a perfectly conducting cylinder surrounded by air which has been shown not to support any bound solutions31.

For all other cases we have solved Eq. (13) using Wolfram Mathematica. Due to the system being lossless, any bound mode will have a real \(n_{eff}\) with \(\sqrt{\varepsilon _{r,1}}<n_{eff}<\sqrt{\varepsilon _{r,2}}\) . Note that a large value of \(n_{eff}\) generally indicates a stronger confinement of the wave to the wire and substrate. We examine the effects of the substrate thickness and dielectric constant, signal frequency and the wire’s height above the substrate on the propagation constant by varying their values but keeping all other parameters constant. Nominal parameter values are \(a=0.1\) mm, \(b=1.6\) mm, \(\varepsilon _{r,1}=1\), \(\varepsilon _{r,2}=3\), \(h=0.1\) mm and frequency \(f=10\) GHz (cf. Fig 1). Our results are validated against finite element numerical solutions obtained with COMSOL Multiphysics32. Details on the simulation are given in the supplementary material. Figure 2 shows the results of the parameter sweeps. It can be seen from the figure that the effective refractive index and in turn the propagation constant crucially depends on all input parameters.

For instance, the height of the wire above the substrate influences how much electromagnetic energy can travel inside the dielectric substrate. In general, if the conductor is further away from the substrate, less energy travels in the dielectric. This means that the propagating mode has an effective refractive index closer to that of the surrounding medium. Consequently, as the conductor approaches the substrate, the effective refractive index increases as shown in Fig. 2a. Good agreement between our model and results obtained with COMSOL can be seen. However, at small heights COMSOL produces an effective index which is slightly higher than predicted by our method.

In fact, the results obtained with COMSOL in Fig. 2 seem to systematically lie above the results of our model. On the one hand, this deviation may be explained by our model neglecting the exact current distribution in the wire leading to errors such as those predicted in Eq. (14) which were on the order of 1-3% throughout the sweep. On the other hand, as the height becomes very small, COMSOL is forced to use elongated mesh elements between the wire and the substrate which is generally not recommended.

The effect of frequency on the PGL mode is shown in Fig. 2b. At high frequencies the wave localises close to the conductor similar to the classical Sommerfeld or Goubau line. This leads to an increase in \(n_{eff}\) which ensures a fast transverse decay in the surrounding air. Additionally, one can think of increasing the frequency as localising more of the wave energy inside the substrate which slows down the wave. Thus, the PGL is generally dispersive. At very low frequencies, the wire cannot be approximated as a PEC anymore as the skin depth becomes similar to the wire radius. Hence, only values of the refractive index above 1 GHz are reported. Results obtained with our model and COMSOL are within a few percent and both show the same general trend although COMSOL predicts a slightly higher mode index.

Increasing the dielectric constant of the substrate slows the wave down. This effect looks to be nearly linear in the magnitude of the effective dielectric constant as seen in Fig. 2c. Agreement between the model and COMSOL is good at small permittivities but the results deviate increasingly with the dielectric constant of the substrate.

Finally, varying the substrate thickness influences the amount of energy that travels inside the substrate. Increasing the thickness slows the wave down leading to a higher effective refractive index. Fig. 2d also includes the effective refractive index of the \(\hbox {TE}_0\) substrate mode as Gacemi et al. reported that with increasing substrate thickness the PGL mode mixes and ultimately merges into this mode23. Indeed the COMSOL result aligns nicely with the \(\hbox {TE}_0\) mode and past a thickness of 7.5 mm only the substrate mode was detected in COMSOL. On the other hand, our model produces results beyond this thickness. However, a noticeable change in the behaviour \(n_{eff}\) with the substrate thickness is observed after intersecting with the \(n_{eff}\) of the substrate mode at around 6 mm. In fact, we were not able to obtain physical field patterns for values of the refractive index after this point. Hence, we believe that while results for larger substrate thicknesses can be calculated they do not hold any physical relevance.

Due to the dependence on geometrical parameters, the planar Goubau line can be designed to have a high or low effective refractive index signifying a strongly or weakly confined mode respectively. Clearly, a mode which is more localised near the substrate will experience increased dielectric loss. Hence, a trade-off between loss and field extent will need to be made. This can to some degree be mitigated by using low-loss substrates for instance quartz or plastics.

Field pattern

Figure 3
figure 3

Plot of \(E_z\) obtained through (a) COMSOL and (b) calculated with the presented analytic model . Excellent agreement between the two field profiles is observed. (c) Contour plot showing \(\log (|E_z|)\) calculated with the analytic model. The field is elongated along the substrate and exponentially decays away from the wire.

Full size image

Once the characteristic equation has been solved, Eq. (7) may be computed at every point in space to calculate the distribution of the z-component of the electric field. Plots of the resulting field are shown in Fig 3 using the parameters given in the previous section. The top two plots show the longitudinal electric field calculated in COMSOL (Fig. 3a) and with our model (Fig. 3b). Both plots are in excellent agreement with each other. It is interesting to observe that the field changes sign between opposite sides of the dielectric. This behaviour was also found in the simulations presented by Horestani et al. for the PGL but was not discussed16. We emphasise that the sign change is unique to the PGL and not found for the classic Goubau line where the mode is cylindrically symmetric. This shows that despite some similarities between the PGL and the classic Goubau line such as an exponential decay at large distances, the presence of a single sided substrate substantially alters the Goubau mode. It breaks the cylindrical symmetry of the system which in turn means no pure TM mode can propagate. As a result only a hybrid mode exists on the PGL.

The logarithmic contour plot in Fig. 3c shows that \(E_z\) decays exponentially away from the wire at distances much greater than the wire radius. This can also be shown directly from our expression for \(E_z\) in Eq. (21). There, we can neglect the integral for \(|\gamma _1 y|\gg 1\) due to the exponential damping. As the incomplete Gamma function is small for large argument, the field is dominated by the modified Bessel functions which have an exponential decay for large, real argument. Thus, the field drops off exponentially with a decay constant \(|\gamma _1|\). This is consistent with previously reported field profiles13,16,17.

Experimental validation

Figure 4
figure 4

(a) Measured and calculated effective refractive index \(n_{eff}\) of the PGL mode over a FR4 substrate as a function of frequency. Error bars correspond to the uncertainty in the theory due to first order azimuthal currents as given in Eq. (14). (b) \(n_{eff}\) of the PGL over a FR4 substrate as a function of height at a frequency of 3 GHz. Here, error bars are included in the measurement data and represent the error in measuring the height in experiment.

Full size image

To validate our results experimentally, we measured the effective refractive index of the PGL mode as a function of frequency. This is achieved with a simple setup. Using scaled versions of the planar launchers discussed by Akalin et al., we excite a Sommerfeld surface wave on a single copper wire1,10. S-Parameters are obtained with a Vector Network Analyser (VNA). Then, we introduce a dielectric slab of finite length l into the path such that the dielectric represents the substrate discussed in our model. In transmission, the substrate will cause a phase delay \(\Delta \varphi\) due to the increased refractive index of the now propagating PGL mode relative to the Sommerfeld mode. The delay can be measured with the VNA. Under the assumption that the Sommerfeld wave travels approximately at the speed of light, the phase delay is given by

$$\begin{aligned} \Delta \varphi =k_0 l \big (n_{eff}-1\big ) \end{aligned}$$
(23)

This is now easily solved to give the effective refractive index of the PGL mode. Note that as we only measure the difference in phase that is introduced by the substrate, this measurement is independent of SW launching as long as a SW is propagating. Further details of the measurement such as wire radius and substrate size are found in the methods section at the end of the paper.

Figure 4a shows the measured refractive index of the PGL with the substrate touching the wire, i.e. \(h=a\), as a function of frequency together with theoretically predicted values. To ensure that the wire touches the substrate, it was taped onto it. Error bars were added to the theoretical values according to Eq. (14) estimating the higher order azimuthal current effects. Good agreement between theory and measurement is observed. Above 20 GHz the effect of azimuthal currents increases drastically reducing the accuracy of the presented theory. Error estimates for the measurement were omitted in the figure for clarity as they were significantly smaller than the theoretical ones.

A further experiment is shown in Fig. 4b. It shows the effective refractive index of the PGL mode against the height of the wire above the substrate. The frequency of this measurement was fixed at 3 GHz. The error bars in the measurement represent the uncertainty in setting the height of the wire above the substrate which was estimated as \(\pm 1\) mm. The accuracy was mainly limited by a slight bend in the substrate. Note that this error does not appear at \(h=a\) as the wire was taped to the substrate eliminating the uncertainty in height. Uncertainty in the theory as given by Eq. (14) was omitted as it becomes negligible at large heights \(h\gg a\). Given the accuracy of the measurement, the agreement between theory and experiment is considered good. However, measurement values are systematically below their theoretical predictions. We interpret this as a systematic error in our measurement of height due to the difficulties described previously.

Considerations for realistic planar Goubau lines

The presented simplified system has many similarities with the standard planar Goubau line such as the reported scaling behaviours with geometrical parameters and frequency. One difference is that in a real system, conductor and dielectric will be lossy leading to a complex effective refractive index. However, the main differences are that the conductor of a standard planar Goubau line has a rectangular cross section and the substrate is of finite width. In many cases the width of the substrate can be neglected as the electromagnetic fields of the bound Goubau mode decay exponentially. Thus, the edges do not have a strong influence on the field pattern. However, the shape of the conductor has been shown to strongly influence the effective refractive index12.

It is beyond the scope of this work to derive a complete theory to incorporate these effects. However, we will try to give some qualitative arguments to describe the observed trends within our framework. If we consider a conductor of finite thickness but variable width, then most of the electric fields and currents will be localised near the edges due to the lightning rod effect. Hence, a natural model is two parallel wires on the substrate which are located at the edges of the rectangular conductor carrying coupled surface waves. This type of coupling has recently been studied for Sommerfeld wires and it was shown to reduce the effective refractive index33. This behaviour is consistent with the observations in Ref.12.

Conclusion

This paper presents a theoretical investigation of a planar Goubau line consisting of a cylindrical wire above an infinitely wide substrate. To this end, the electromagnetic field of an infinitesimal current filament above a substrate has been derived in Fourier space showing that the resultant field will be a hybrid mode. By incorporating the finite width of a realistic wire, an approximate characteristic equation was derived. Estimates to the applicability of this equation were given before exploring the influence of geometrical parameters and frequency on the wave characteristics. It was found that, depending on the operating frequency and setup, both weakly and strongly confined surface waves can propagate allowing one to tune the geometry depending on the desired application. Furthermore, the field profile for a particular set of parameters was calculated. All results agreed well with numerical simulations and previously reported experimental studies. Finally, the derived model was experimentally tested and excellent agreement between theory and measurement was found. The limitations of the presented model with respect to a realistic planar Goubau line with rectangular conductor were discussed in the last section and some qualitative arguments were made to incorporate the effect of the conductor geometry.

To our knowledge this work is the first attempt at analytically modelling the planar Goubau line. Although only a version with circular conductor was discussed, the behaviour with frequency and other parameters was found to agree with the previous reports for a standard planar Goubau line with rectangular conductor. Hence, it may be considered as an important step towards a complete model of the PGL. In conclusion, we expect these insights to help better understand and utilise planar Goubau lines in printed circuit board designs for challenging applications such as terahertz spectroscopy or high frequency circuitry.

Methods

In our experimental setup a 15 cm long, 1.6 mm thick, dielectric slab made from FR4 epoxy in conjunction with a 0.5 mm annealed copper wire were used. The relative permittivity of the substrate is given as 4.55 by the manufacturer. For the measurement where wire and substrate are touching, the wire was taped onto the substrate to ensure contact. Contrary, in the measurement with different heights, we used ROHACELL as support for the FR4 epoxy substrate. ROHACELL has a relative permittivity close to one. Therefore, we assume it does not influence the electromagnetic fields. S-parameters were measured by a 8722D VNA from Agilent Technologies.

Calculations based on the presented theoretical model were carried out in Wolfram Mathematica. The principal value of the integral in (21) is taken explicitly with a value of \(\epsilon =10^{-5}\) after finding the roots of the denominator using a numerical root finding algorithm. Equation (13) is solved by calculating \(E_z\) for different \(n_{eff}\) until it changes sign. Then, we interpolate the data to find the value of \(n_{eff}\) which solves (13). For details on the COMSOL simulation, please refer to the supplementary material.