Stable pseudoanalytical computation of electromagnetic fields from arbitrarily-oriented dipoles in cylindrically stratified media

TL;DR

This paper presents a stable computational method for electromagnetic fields from arbitrarily-oriented dipoles in cylindrically stratified media, overcoming numerical issues in traditional approaches through range-conditioned functions and adaptive integration.

Contribution

It introduces range-conditioned cylindrical functions and adaptive spectral integration techniques for stable, accurate electromagnetic field computations in complex layered media.

Findings

Successfully computes fields across wide parameter ranges

Prevents overflow and underflow in numerical calculations

Applicable to geophysical problems with diverse layer properties

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 finite-precision 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 high-order, 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 range-conditioned, modified cylindrical functions (in lieu of standard cylindrical eigenfunctions), each associated with non-overlapped subdomains of (numerical) evaluation to allow for stable computations under any range of physical parameters. In addition adaptively-chosen 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 1000 $\Omega \cdot$m to 10$^{-8} \Omega \cdot$m, frequencies of operation ranging from 10 MHz down to the low magnetotelluric range of 0.01 Hz, and for various combinations of layer thicknesses.