Stable pseudoanalytical computation of electromagnetic fields from arbitrarilyoriented dipoles in cylindrically stratified media
Abstract
Computation of electromagnetic fields due to point sources (Hertzian dipoles) in cylindrically stratified media is a classical problem for which analytical expressions of the associated tensor Green’s function have been long known. However, under finiteprecision arithmetic, direct numerical computations based on the application of such analytical (canonical) expressions invariably lead to underflow and overflow problems related to the poor scaling of the eigenfunctions (cylindrical Bessel and Hankel functions) for extreme arguments and/or highorder, as well as convergence problems related to the numerical integration over the spectral wavenumber and to the truncation of the infinite series over the azimuth mode number. These problems are exacerbated when a disparate range of values is to be considered for the layers’ thicknesses and material properties (resistivities, permittivities, and permeabilities), the transverse and longitudinal distances between source and observation points, as well as the source frequency. To overcome these challenges in a systematic fashion, we introduce herein different sets of rangeconditioned, modified cylindrical functions (in lieu of standard cylindrical eigenfunctions), each associated with nonoverlapped subdomains of (numerical) evaluation to allow for stable computations under any range of physical parameters. In addition, adaptivelychosen integration contours are employed in the complex spectral wavenumber plane to ensure convergent numerical integration in all cases. We illustrate the application of the algorithm to problems of geophysical interest involving layer resistivities ranging from 1,000 to , frequencies of operation ranging from 10 MHz down to the low magnetotelluric range of 0.01 Hz, and for various combinations of layer thicknesses.
keywords:
cylindrically stratified media, tensor Green’s function, cylindrical coordinates, electromagnetic radiation[cor1]Corresponding author
1 Introduction
Computation of electromagnetic fields due to arbitrarilyoriented elementary (Hertzian) dipoles in cylindrically stratified media is of interest in a wide range of scenarios, including geophysical exploration, fiber optics, and radar crosssection analysis. Assuming the axis to be the symmetry axis, the analytical formulation of this problem is predicated on the knowledge of the cylindrical eigenfunctions (Bessel and Hankel functions) in the domain transverse to and their modal amplitudes. The derivation of reflection and transmission coefficients at each cylindrical boundary is then ascertained through the use of the proper boundary conditions. Since the eigenfunctions comprise a continuum spectrum in an unbounded domain, a Fouriertype integral along the spectral wavenumber is subsequently necessary to determine the fields (tensor Green’s function) Chew (1983); Lovell and Chew (1987, 1990); Tokgöz and Düral (2000); Hue and Teixeira (2006); Wang et al. (2008); Liu et al. (2012). Unfortunately, numerical computations based on direct use of such (canonical) analytical expressions often lead to underflow and/or overflow problems under finiteprecision arithmetic. These problems are related to the poor scaling of cylindrical eigenfunctions for extreme arguments and/or highorder, as well as convergence problems related to the numerical evaluation of the spectral integral on and truncation of the infinite series over the azimuth mode number . Underflow and overflow problems become especially acute when a disparate range of values needs to be considered for the physical parameters, viz., the layers’ constitutive properties (resistivities, permittivities, and permeabilities) and thicknesses, the transverse and longitudinal distance between source and observation points, as well as the source frequency Yousif and Boutros (1992). In order to stabilize the numerical computation in the case of planewave scattering by highly absorbing layers, Swathi and Tong Swathi and Tong (1988) developed an algorithm based upon scaled cylindrical functions. A stabilization procedure to deal with a very large number of cylindrical layers and disparate radii was proposed in Jiang et al. (2006). Similar issues appear when computing Mie scattering from multilayered spheres Kaiser and Schweiger (1993),Peña and Pal (2009), where continued fractions Lentz (1976) or logarithmic derivatives Toon and Ackerman (1981), for example, can be used to circumvent the recurrence instability of Bessel functions of very large order. When the overall computation cost is not an issue, a more extreme strategy to circumvent this problem is to use arbitraryprecision arithmetic Suzuki and Sandy Lee (2012).
In this work, we extend those efforts by considering the stable computation of electromagnetic fields due to points sources in cylindrically stratified media. A salient feature of our work is that we do not limit ourselves to a particular regime of interest (that is very small or very large radii, or highly absorbing layers) but instead develop a systematic algorithm to enable stable computations in any scenario. Note that, as opposed to Mie scattering or planewave scattering from cylinders considered in the above, the required spectral integration over produces, per se, a large variation on the integrand function arguments. The numerical convergence of such spectral integral, which depends among other factors on the separation between source and observation points, need to be considered in tandem with the numerical stabilization procedure. Our methodology is based on the use of various sets of rangeconditioned, modified cylindrical functions (in lieu of standard cylindrical eigenfunctions), each evaluated in nonoverlapped subdomains to yield stable computations under doubleprecision floatingpoint format for any range of physical parameters. This is combined with different numericallyrobust integration contours that are adaptively chosen in the complex plane to yield fast convergence. We illustrate the algorithm in problems of geophysical interest involving layer resistivities ranging from 1,000 to about , frequencies ranging from 10 MHz to as low as 0.01 Hz, and for various layer thicknesses.
2 RangeConditioned Formulation
Figure 1 shows the geometry of a cylindrically stratified medium. Hereinafter we shall refer to the layer where the point source is present as layer and the region where the fields are computed as layer . As indicated in Fig. 1, each successive layer radius is denoted by .
We use nonprimed coordinates for the observation point and primed coordinates for the source location, denotes the transverse wavenumber in layer and is the longitudinal wavenumber, where .
As usual, represents the angular frequency, the complex permeability, and the complex permittivity, with , denoting the realvalued permittivity and permeability, resp., and , denoting the electric and magnetic conductivities, resp.
(1) 
where is the dipole moment, and
(2) 
is an operator acting on the primed variables at the left, with being a unit vector corresponding to the dipole orientation. The transverse field components can be easily derived from the the knowledge of the components.
Four distinct expressions for the generic factor in the integrand of (1) exist depending on the relative position of the source and observation points, as follows (Chew, 1995, p. 176):
Case 1: and are in the same layer and  
(3a)  
Case 2: and are in the same layer and  
(3b)  
Case 3: and are in different layers and  
(3c)  
Case 4: and are in different layers and  
(3d) 
where and are cylindrical Bessel and Hankel functions of the first kind and order , is the identity matrix, , and with representing the local reflection matrix between two adjacent cylindrical layers, representing the generalized reflection matrix between two adjacent cylindrical layers, and representing generalized transmission matrix, respectively. The reflection and transmission coefficients between cylindrical layers are 22 matrices because both TE and TM waves are in general needed to match the boundary conditions at the interfaces. Expressions for , , are found further down below and are also provided in (Chew, 1995, ch. 3). A note should be made that insofar as the expressions for Cases 3 and 4 are concerned, the leftmost factors found in (Chew, 1995, p. 176) are incorrect and should instead be replaced by , as shown above.
As alluded before, even though (1) provides an exact analytical expression, in many instances the wildly disparate behavior in magnitude of the many different factors that comprise the integrand makes it difficult to obtain accurate results numerically. Furthermore, care should be exercised in choosing the integration path in the complex plane so that a convergent numerical integration is obtained, in a robust fashion. These two challenges are tackled in the remainder of this paper.
2.1 Rangeconditioned cylindrical functions
The evaluation of products of and (and their derivatives) is needed for the computation of the integrals in (3a) – (3d) (and also in similar integrals for the computation of the transverse field components). When , has a very large value whereas has a very small value, and this disparity becomes even more extreme for larger order. On the other hand, when (a condition that arises, for example, at low frequencies in layers with small resistivity), has a very small value while has a very large value. In cylindrically stratified media, these disparities in magnitude can coexist, to a varying extent, in many different layers. In addition, since the magnitude of the arguments for and also depends on , large variations occur in the course of numerical integration over as well. As a result, significant roundoff errors can accumulate unless the inherent characteristics of such functions are taken into consideration beforehand.
Definitions
When , and can be expressed through the following small argument approximations for (Chew, 1995, p. 15).
(4a)  
(4b)  
(4c)  
(4d) 
where . In (4a) – (4d), , , , and are the rangeconditioned cylindrical functions for small arguments. Note that the definition is such that the multiplicative factors and associated with a given cylindrical function and its derivative are the same and the multiplicative factors for and are reciprocal to ones for and . These multiplicative factors play a determinant role in producing the extreme values for the cylindrical eigenfunctions and, as we will see below, should be analytically manipulated (reduced) in an appropriate fashion before any numerical evaluations to stabilize the computation. As seen below, these definitions facilitate subsequent computations.
On the other hand, when , large argument approximations can be used for the cylindrical eigenfunctions (Zhang and Jin, 1996, p. 131, 236). Likewise, rangeconditioned cylindrical functions for large arguments can be defined as
(5) 
(6) 
where , , is the modified Bessel function of the second kind, and and are phase and quadrature polynomial functions with tabulated expressions. In the above, we used the fact that, for , and in (5) can be decomposed into two terms and one of them can be ignored when is large enough such that
(7a)  
(7b) 
Derivatives of rangeconditioned cylindrical functions with large arguments can be obtained through the recursive formulas below
(8)  
(9) 
Note again that the associated multiplicative factors in (5) and (6), as well as in (8) and (9) are chosen to be reciprocal to each other.
When argument is neither too small nor too large, rangeconditioned cylindrical functions are introduced akin to those of small and large arguments, i.e.,
(10a)  
(10b)  
(10c)  
(10d) 
where is determined in B, and its first and second subscript correspond to the radial wavenumber and layer radius, respectively.
In summary, the argument for rangeconditioned cylindrical functions can be classified into three types according to its magnitude: small, moderate, and large. The multiplicative factors that define the respective rangeconditioned functions for each type of argument are summarized in Table 1.
small arguments  moderate arguments  large arguments 






As shown in A, the local reflection and transmission matrices and are written in terms of and matrices. The latter are defined in terms of cylindrical functions so that rangeconditioned versions thereof can be constructed as well. After some algebra, it can be shown that the multiplicative factors associated with the and matrices are the same as for the rangeconditioned cylindrical functions. That is, for small arguments we have
(11a)  
(11b) 
for large arguments we have
(12a)  
(12b) 
and for moderate arguments we have
(13a)  
(13b) 
For a numerical computation, the actual threshold values among the three argument types above need to be determined. This is done in B assuming standard doubleprecision arithmetics.
2.2 Reflection and transmission coefficients for two cylindrical layers
When a medium consists of two cylindrical layers, the expressions for the reflection , and transmission , coefficients at boundary are given in (Chew, 1995, ch. 3) and, for convenience, also presented in A. When one or two of the involved layers are associated with small or large arguments, those coefficients need to be rewritten using rangeconditioned cylindrical functions. Therefore, there are total of nine cases to be considered according to the argument types as listed in Tables 4 – 4. The first group (Cases 1, 2, 3) only incorporates small and moderate arguments; the second group (Cases 4, 5, 6) only incorporates large and moderate arguments; and the third group (Cases 7, 8, 9) consists of remaining combinations.
Case 1  small  small 
Case 2  small  moderate 
Case 3  moderate  small 
Case 4  large  large 
Case 5  large  moderate 
Case 6  moderate  large 
Case 7  small  large 
Case 8  large  small 
Case 9  moderate  moderate 
The canonical expressions for the (local) reflection and transmission coefficients in A can be rewritten using rangeconditioned cylindrical functions for all the above nine cases. The redefined (conditioned) coefficients for all three groups are summarized in Table 5. Note that for the second (Cases 4, 5, 6) and third groups (Cases 7, 8, 9), the redefinition of the coefficients is basically similar to the first group but the associated multiplicative factors are different.
Case 1  

Case 2  
Case 3  
Case 4 

Case 5  
Case 6  
Case 7 

Case 8  
Case 9 
Explicit expressions for , , , and are also provided in A. It should be noted that there is a simple relationship between the original coefficients and the rangeconditioned coefficients in all cases. For any two arbitrarilyindexed cylindrical layers, the relationship can be succinctly expressed as