# On the Surface Fields of a Small Circular Loop Antenna Placed on Plane Stratified Earth

Authors: Mauro Parise
Journal: International Journal of Antennas and Propagation
 (2015)

---

## Abstract

An analytical method is presented which makes it possible to derive exact explicit expressions for the time-harmonic surface fields excited by a small circular loop antenna placed on the top surface of plane layered earth.
The developed procedure leads to casting the complete integral representations for the EM field components into forms
suitable for application of Cauchy’s integral formula. As a result, the surface fields are expressed as sums of Hankel
functions. Numerical simulations are performed to show the validity and accuracy of the proposed solution.


---

## Body

## 1. Introduction

The need to evaluate the electromagnetic fields radiated by horizontal loop antennas located in close proximity to layered media arises in a variety of applications of scientific interest, especially in the areas of electromagnetics [1–9], close-to-the-surface radio communication [10, 11], and geophysical prospecting [11–13]. An excellent illustration of such applications is electromagnetic induction sounding technique [12, 14, 15], which allows acquiring information about the subsurface structure of a terrestrial area. The technique consists of two steps. At first, the fields generated by a current-carrying circular loop antenna lying on the surface of the medium to be explored are measured at a discrete set of frequencies. Next, the presence of shallow buried objects, or the geometrical and EM properties of the structure, is deduced by interpreting the recorded experimental data through matching with the theoretical response curves associated with standard homogeneous or layered earth models [12, 14, 15]. It is easily understood how, for correct data interpretation, high accuracy is required when computing such theoretical response curves.This paper will focus on the problem of an electrically small circular loop antenna (vertical magnetic dipole) lying on the surface ofN-layer plane stratified earth. To date, the fields of a loop source in the presence of a layered ground are not known in closed form, and evaluation of their complete integral representations has been revealed to be challenging from both analytical and numerical points of view. When available, analytical solutions either pertain to simple earth models, like the single-layer configuration [16], or follow from introducing simplifying assumptions. Examples of approximate solutions are the explicit series-form expressions that result from applying standard asymptotic techniques, such as the method of steepest descent (saddle-point method) [11, 12, 17–19]. The obtained formulas are easy to use but are valid under the restrictive assumption that the source-receiver distance is electrically large, that is, in the far-field region of the antenna. Similar considerations can be made with regard to the low-frequency or quasi-static expressions for the fields [7, 12, 13, 20], which are accurate, subject to the condition that the displacement currents are negligible in both the layered material medium and the air space above it.On the other hand, numerical computation of the field integrals through standard numerical approaches, like Gaussian quadrature, is often made difficult and impractical by the oscillating behavior of the integrands and the slow convergence of their envelope, especially when both source and field points lie at the air-medium interface [1, 12, 21]. In an attempt to overcome these drawbacks, some papers document extensive usage of the digital linear filter technique in place of the Gaussian quadrature approach [22, 23]. Unfortunately, it has been recently shown [9, 24, 25] that any digital filter is designed and optimized for a specific known integral transform (the so-called test function) and may lead to uncertain results when applied to the calculation of different integrals.The aim of this work is to derive exact explicit expressions for the fields excited by the loop source at the air-medium interface. The theoretical development is organized as follows. First, the Sommerfeld Integrals describing the fields are split into two terms, that is, the residue contributions from the poles of the integrands and the branch-cut integrals due to the propagation constants in free space and in the bottom layer of the stratified medium. Next, such propagation constants are extracted from the integrands of the branch-cut integrals and replaced with the sums of partial fractions generated by Newton’s iterative method [26, 27]. This makes it possible to apply Cauchy’s integral formula and obtain explicit series representations for the branch-cut integrals. As a result, the Sommerfeld Integrals are converted into infinite sum of Hankel functions. The advantage of the presented approach resides in the fact that the sequence of the sums of partial fractions that are used in place of each propagation constant is convergent, as Newton’s method exhibits quadratic rate of convergence if applied to the computation of a square root [26]. This implies that the derived series representations for the fields are exact.


## 2. Theory

Consider a vertical magnetic dipole of momentmejωt lying on the surface of a flat N-layer lossy medium, as shown in Figure 1. Due to symmetry about the axis of the dipole, a cylindrical coordinate system (r,φ,z) is suitably introduced. The dielectric permittivity and electric conductivity of the nth layer are denoted by ϵn and σn, respectively, while the medium is assumed to have the magnetic permeability of free space μ0. The time-harmonic EM field components generated by the dipole in the air space (the zeroth layer, z<0) may be derived from an electric vector potential which has only z component F. If the time factor ejωt is assumed, it reads(1)Eφ=∂F∂ρ,Hρ=1jωμ0∂2F∂ρ∂z,Hz=-1jωμ0∂2∂ρ2+1ρ∂∂ρF,where F=Fd+Fr, with(2)Fd=jωμ0m8π∫∞ejπ∞eu0zλu0H01λρdλ,(3)Fr=jωμ0m8π∫∞ejπ∞RTEeu0zλu0H01λρdλ,respectively, being the direct and reflected potentials. In (2) and (3), H0(1)(·) is the zeroth-order Hankel function of the first kind, while(4)RTE=Y0-Y^1Y0+Y^1is the plane wave reflection coefficient at z=z0=0, with Y0 and Y^1 being the intrinsic admittance of free space and the surface admittance at z0. The intrinsic admittance of the nth layer is given by(5)Yn=unjωμn,with(6)un=λ2-kn2,Re⁡un>0,(7)kn2=ω2μnϵn-jωμnσn,while, for N layers, the surface admittance at z=zn-1 is obtained through the recurrence relation(8)Y^n=YnY^n+1+Yntanh⁡unzn-zn-1Yn+Y^n+1tanh⁡unzn-zn-1,n=N-1,…,1, with(9)Y^N=YN.The reflected field may be split into two parts, that is, the ideal reflected field (the field of a negative image) and a contribution due to the imperfect conductivity of the earth. Use of the identity(10)RTE=-1+2Y0Y0+Y^1in (3) allows concluding that the direct and image fields exactly cancel out at all points in the air space. Hence, the total potential reads(11)F=jωμ0m4π∫∞ejπ∞Y0eu0zY0+Y^1λu0H01λρdλ=jωμ0m4π∫∞ejπ∞eu0zu0+jωμ0Y^1H01λρλdλ.The fraction under the integral sign has a number of pole singularities, plus the branch point singularities at λ=±k0 and λ=±kN. It is interesting to note that the points λ=±k1,λ=±k2,…,λ=±kN-1 are not branch points of the integrand, as the fraction term is an even function of u1,u2,…,uN-1. We now deform the integration path to the upper infinite semicircle C centered at the origin of the complex λ-plane, leaving out all the singularities located in the first and second quadrants, as indicated in Figure 2. The result of the deformation is a new integration path constituted by C, the two lines wrapped around the branch cuts, that is, Γ0 and ΓN, and a number of infinitesimal circles around the poles of the integrand. It is not difficult to prove that the infinite semicircle does not contribute to the integral. In fact, for large values of λ, the Hankel function may be replaced with its asymptotic form [28](12)H01λρ≅2jπλρejλρ,and the integrand in (11) becomes suitable for application of Jordan’s Lemma [29], which states the following upper bound:(13)∫Cλeu0zu0+jωμ0Y^1ejλρdλ≤πρmaxC⁡λeu0zu0+jωμ0Y^1.Since, in the air space, factor eu0z decays exponentially with increasing λ, the right-hand side of (13) must vanish, and so does the integral along C. As a consequence, F may be expressed as(14)F=jωμ0m4π∫Γeu0zu0+jωμ0Y^1H01λρλdλ,where Γ is the contour consisting of Γ0, ΓN, and the infinitesimal circles. It should be observed that the branch-cut integrals along Γ0 and ΓN describe the above-surface ground wave and the lateral-wave field, respectively, while the integrals around the poles are trapped-surface-wave terms [17].Figure 1
Sketch of a vertical magnetic dipole on a multilayer ground.Figure 2
Deformed integration contour for the Sommerfeld Integrals.The scope of this work is to exactly evaluate the field components at the air-ground interface. To this goal, we substitute (14) into (1) and take the limit as z→0-. It yields (15)Eφ=jωμ0m4π∂S0∂ρ,(16)Hρ=m4π∂S1∂ρ,(17)Hz=-m4π∂2∂ρ2+1ρ∂∂ρS0,with(18)Sh=∫ΓfhλH01λρλdλ,fhλ=u0hu0+jωμ0Y^1.Quantity Sh may be rewritten as(19)Sh=Sh0+ShN+Shp,where(20)Shn=∫ΓnfhλH01λρλdλ,n=0,N,is the integral along Γn, while Sh(p) is the contribution from the residues of the integrand at the poles of fh. The explicit expression for Sh(p) follows from applying Cauchy’s theorem [29]. It yields(21)Shp=2πj∑iRifhλH01λρλ=2πj∑iH01λiρλid/dλfhλ-1λ=λi, where Ri· denotes the residue at the generic pole λi(i=1,2,…). On the other hand, to evaluate Sh(0) and Sh(N), it is convenient to extract the square root terms u0 and uN from fh. This task may be accomplished through letting fhλ=ϕh(u0,uN) and then applying the branch-cut extraction procedure:(22)ϕh=ϕhee+ϕheouN+ϕhoeu0+ϕhoou0uN,where(23)ϕhee=14ϕh+++ϕh-++ϕh+-+ϕh--,(24)ϕheo=14uNϕh+++ϕh-+-ϕh+--ϕh--,(25)ϕhoe=14u0ϕh++-ϕh-++ϕh+--ϕh--,(26)ϕhoo=14u0uNϕh++-ϕh-+-ϕh+-+ϕh--,with(27)ϕh±±=ϕh±u0,±uN.It is straightforward to verify that functions ϕh(ee), ϕh(eo), ϕh(oe), and ϕh(oo) are even in both u0 and uN and, hence, do not exhibit branch cuts. Next, we use Newton’s method to convert u0 and uN into rational functions, so as to replace the branch cuts running from λ=-k0 and λ=-kN with a set of pole singularities. This will make it possible to exactly evaluate the integrals through contour integration. According to Newton’s method, the principal square root of the complex number un2 (i.e., the principal root of the equation ξ2-un2=0) may be obtained by taking the limit of the sequence(28)u~nl=12u~nl-1+un2u~nl-1,u~n0=un2,as l→∞. Recursively applying (28) provides the explicit form of u~n(l); namely,(29)u~nl=2-lPlun2∏i=1l-1Piun2,with(30)Piq=∑k=02i-12i2kqk.Thus, substitution of (29) for un’s in (22) makes ϕh become a meromorphic function of λ, whose poles coincide with those of u~0(l) and u~N(l). The poles of u~n(l)(n=0,N) with positive imaginary part, Ml=2l-1-1 in number, are given by(31)λmnl=-pml+kn2,with pm(l)(m=1,2,…,Ml) being the generic root of the equation(32)∏i=1l-1Piq=0.Notice that 2i-1 zeros of Pi(q) are negative real numbers, and this implies that poles (31) are positioned on the suppressed branch lines, which are described by the equation(33)λ2-kn2=αn,n=0,N,with αn being a negative real parameter. As shown in Figure 3, substitution (29) has the effect of turning the branch-cut integrals Sh(0) and Sh(N) into contour integrals around poles (31). It reads(34)Shn=∑m=1MlΦhmn,n=0,N,where(35)Φhm0=∮Γm0ϕhoeu~0+ϕhoou~0u~NH01λρλdλ,ΦhmN=∮ΓmNϕheou~N+ϕhoou~0u~NH01λρλdλ,with Γmn(n=0,N) being the infinitesimal circle around λmn, and where the superscript l has been omitted for notational simplicity. It should be observed that λmn is not a pole of ϕh(ee) and, hence, it results as(36)∮ΓmnϕheeH01λρλdλ=0.Substituting (29) for u~0 and u~N in (35) provides the expressions(37)Φhm0=ϕhoe+ϕhoouNH01λρλλ=λm0∮Γm0u~0dλ,ΦhmN=ϕheo+ϕhoou0H01λρλλ=λmN∮ΓmNu~Ndλ,where we have taken account of the fact that, far from the relevant poles, quantities u~0 and u~N may be confused with u0 and uN, respectively. Since from (6), (29), and (31) it follows that(38)limλ→λmn⁡u~n=2-lPlpmRmlimλ→λmn⁡1un2-pm=2-lPlpmRmlimλ→λmn⁡1λ2-λmn2,with(39)Rm=∏m′=1m′≠mMlpm-pm′-1,applying Cauchy’s integral formula [29] provides(40)∮Γmnu~ndλ=2-lPlpmRm∮Γmn1λ2-λmn2dλ=πjPlpmRm2lλmn.Thus, (37) are turned into(41)Φhm0=πjPlpmRmvhλm0H01λm0ρ2l,(42)ΦhmN=πjPlpmRmvhλmNH01λmNρ2l,with(43)vhλ=ϕhoe+ϕhoouN,λ=λm0,ϕheo+ϕhoou0,λ=λmN,0,elsewhere.Combining (41) and (42) with (34) finally yields the expression(44)Sh0+ShN=πj2l∑m=1MlPlpmRmvhλm0H01λm0ρ+vhλmNH01λmNρ,where one can clearly distinguish the above-surface ground wave, associated with λm0 terms, and the lateral wave, constituted by λmN terms, the latter traveling in the Nth layer of the medium along its upper boundary.Figure 3
Final integration contour for the Sommerfeld Integrals.Substitution of (44) into (19), in conjunction with (21), provides explicit expressions, in series form, for integrals S0 and S1 appearing in (15)–(17).


## 3. Validation of the Formulation

To test the validity of the developed theory, formulas (15)–(17) and (19), (21), and (44) are applied to the computation of the magnitude of the magnetic field components generated on the top surface of a homogeneous medium at distance ρ from a unit-moment vertical magnetic dipole. In the considered case N=1, and functions fh(λ), which reduce to(45)fhλ=u0hu0+u1,do not exhibit poles. This implies that terms (21) do not contribute to the integrals. Thus, applying (24)–(26) provides(46)ϕ0eo=ϕ1oo=-ϕ0oe=1k02-k12,ϕ1eo=ϕ1oe=ϕ0oo=0,and integrals (19) assume the form(47)S0=πj2lk02-k12∑m=1MlPlpmRmH01λmnρn=0n=1,S1=πj2lk02-k12∑m=1MlPlpmRm-1npm+-1nk02-k12H01λmnρn=1n=0.The numerical results that arise from using (16) and (17), in combination with (47), are depicted in Figures 4–7 and denoted by lines. Such results are compared with those provided by the previously published exact expressions for Hρ and Hz; namely [16],(48)Hρ=mπρα2+β22K1αρI1βρ-αβK2αρI2βρ,Hz=Q1-Q02πk02-k12,where In(·) and Kn(·) are the nth-order modified Bessel functions of the first and second kind, respectively, and(49)α=12jk1+k0,β=12jk1-k0,with(50)Qn=-jkn3ρ3-4kn2ρ2+9jknρ+9e-jknρρ5.Figure 4
Amplitude-frequency spectrum ofHρ, computed at ρ=20m from the source point.Figure 5
Amplitude-frequency spectrum ofHz, computed at ρ=20m from the source point.Figure 6
Profiles ofHρ against the radial distance from the source point.Figure 7
Profiles ofHz against the radial distance from the source point.Figures4 and 5 present amplitude-frequency spectra of Hρ- and Hz-fields generated at ρ=20m, obtained assuming σ1=1mS/m and ϵ1=10ϵ0. Each figure illustrates the behavior of the sequence of spectra generated by the proposed solution as the number of iterations of Newton’s method (l) is increased. As is seen, the sequences converge, and the number of iterations required to achieve perfect agreement with the exact data, denoted by points, is 10 for Hρ-component and 5 for Hz-component. This means that the sums S1 and S0 in (47) must be made up of M10=29-1=511 terms and M5=24-1=15 terms, respectively.It should be observed that the dominant contributions to the sums in (47) come from the residues at the poles located in the portions of the branch cuts which are closer to the real axis, since the magnitude of H0(1)(λρ) is significantly larger at them. In particular, the ground-wave poles (λm0’s) that lie on the negative real λ-axis contribute most to S0 and S1. For fixed l, the number of such poles diminishes as frequency is decreased, and this happens because the branch point -k0 gets closer to the origin of the complex plane. The consequence of this phenomenon is that, at low frequencies, the number of dominant residues may not be sufficient to achieve accurate results. This is the reason why, as illustrated in Figure 4, accurate evaluation of Hρ requires 15 terms if frequency is higher than 10 MHz and at least 511 terms when frequency is lower than 100 kHz.Notice that the onset of oscillations in the frequency spectra is expected, since, when entering the far zone, the total field results from the superposition of two waves, which propagate with wavenumbersk0 and k1 and decay with distance ρ from the source point. Physically, this means that oscillatory interference patterns arise in both space and frequency domains. In the space domain, the period of the oscillations is 2π/(k1-k0), with k1 approaching a real quantity since, at high frequencies, the displacement currents predominate over the conduction currents. In the frequency domain, the patterns oscillate with period (51)Λ=1ρμ0ϵ0ϵr1-1,which, for ρ=20m and ϵr1=10, is equal to(52)Λ=3·1082010-1≅6.94MHz.It may be numerically verified that 6.94 is just the peak-to-peak spacing of the far-zone asymptotic trend of each solid-line curve in Figures 4 and 5.Finally, Figures6 and 7 show profiles of the amplitudes of Hρ and Hz against the source-receiver distance ρ. The calculations have been performed assuming that the loop antenna is positioned on the same medium as in the previous examples and that it operates at 10 MHz. The plotted curves demonstrate that even when taken as functions of ρ, the sequences of partial sums in (47) are convergent and give rise to accurate results after only 8 (Hρ-field) and 9 (Hz-field) iterations of Newton’s method.


## 4. Conclusion

This paper has presented exact series representations for the surface fields excited by a small circular loop antenna placed on the top surface of plane stratified earth. The derived expressions result from a rigorous analytical procedure, which allows casting the complete integral representations for the EM field components into forms suitable for application of Cauchy’s integral formula. Numerical results are presented to show the validity and accuracy of the proposed formulation.



---

*Source: 187806-2015-10-28.xml*