1. Introduction

The computation of electromagnetic (EM) wave scatterings of spherically layered media is a canonical problem encountered in both science and engineering research, such as material and optical physics [1], chemistry [2], atmospheric physics [3], geophysics [4], antenna engineering [5], etc. Both Lorentz-Mie theory and spherically layered media theory deal with EM wave scatterings of a sphere with inhomogeneous media distributions along radial direction [6–10]. The major difference is that Lorentz-Mie theory deals with plane wave incidence, whereas spherically layered media theory can handle arbitrary incident waves. One important application of spherically layered media theory is to compute the dyadic Green’s function of spherically layered media [11–15].

Both Lorentz-Mie theory and spherically layered media theory suffer from the notorious numerical instability issue due to the involved Bessel and Hankel functions with small argument, large order or high loss [16–18]. The logarithmic derivative method has been proposed to successfully solve the numerical instability issues for Lorentz-Mie theory first [19–21] and then spherically layered media theory very recently [18,22],

It is well known that asymptotic formulas exist for Bessel and Hankel functions when the argument goes to zero or infinity or the order is quite large. Hence, if the Bessel and Hankel functions are represented by the asymptotic formulas under these cases, the diverging terms in spherically layered media theory can be canceled out so that the calculation will be stable. In addition, the asymptotic formula of Bessel (and Hankel) function under the case of small argument is the same with that under the case of large order [23,24], hence a unified method can be achieved to solve the numerical issue. It is the purpose of this paper to solve both the small-argument issue and the large-order issue of the spherically layered media theory by adopting these asymptotic formulas. It will be shown that the proposed asymptotic method enjoys the feature of simpler formulas and higher efficiency compared to Ref. [18]. The issue under large lossy case can be solved similarly or by employing the scaled Bessel and Hankel functions provided by the popular softwares like Matlab and GNU Octave [25]. As an application example, the stable formulas are used to derive the stable field formulations for a point dipole source excitation in stratified media, which is related to the dyadic Green’s function. Finally, numerical tests are performed to verify the derived stable formulas. Time factor $e^{-i\omega t}$ is used throughout this paper.

2. Stable formulations of the spherically layered media theory for large orders and small arguments

In Ref. [18], we proposed the renormalized spherically layered media where the reflection and transmission coefficients are renormalized so as to have ordinary magnitudes. This definition can avoid the numerical diverging of the coefficients themselves. However, their calculations are still numerically instable due to the involved spherical Bessel functions with small argument or large order. This issue was solved by using the logarithmic derivative method in Ref. [18]. To solve this problem alternatively, the spherical Bessel functions can be approximated by their asymptotic formulas first, and then the diverging terms in the asymptotic formulas are canceled out to obtain the stable formulas. These diverging terms can always be canceled out because the renormalized coefficients have finite values.

From the ascending series of spherical Bessel functions [26], it is known that for integer $n>0$, the asymptotic expansions of the spherical Bessel functions for fixed $n$ and $\lvert x \rvert \rightarrow 0$ are [23,24]

The asymptotic expansions of the spherical Bessel functions for fixed $x$ and $n\rightarrow +\infty$ are

As shown in Eqs. (1) and (2), the asymptotic behaviors of the spherical Bessel functions under small argument case and under large order case can be uniformly represented by the same leading order term. Hence, the numerical instability problems due to small argument or large order can be treated in a unified manner. When the ratio $n/|x|$ is large, either $n$ is large, $|x|$ is small or both, so that the leading term in the asymptotic expansion dominates.

2.1 Stable formulas of the renormalized reflection and transmission coefficients

As shown in Ref. [18], the renormalized and generalized reflection and transmission coefficients (i.e. $\tilde {R}^\prime$ and $\tilde {T}^\prime$) are the functions of the renormalized single-interface reflection and transmission coefficients (i.e. $R^\prime$ and $T^\prime$). Hence, we only need to derive the stable formulas for the renormalized single-interface coefficients. When the single-interface coefficients can be computed stably, the generalized ones can also be computed stably.

2.1.1 Out-going wave incidence

Applying the above asymptotic formulas, the stable formulas can be derived for different cases. To produce a non-zero electromagnetic field, the order $n$ of the Bessel functions in Debye potentials shall be larger than $0$. When $n/\lvert k_1a \rvert \to +\infty$, we have

Here, we use ’$+\infty$’ to mean that the ratio is sufficiently large so that the asymptotic approximation is valid, rather than that the ratio actually goes to infinity. $R^\prime _{12}$ and $T^\prime _{12}$ are the renormalized single-interface reflection and transmission coefficients with the out-going incident wave from layer 1 to layer 2 [18]. $\epsilon _i (\mu _i),i=1,2$ are the permittivity (permeability) of region 1 and region 2 respectively. The wavenumber $k_i=\omega \sqrt {\epsilon _i\mu _i}$. $a$ is the radius of the interface. $\hat {J}_n(x)=xj_n(x)$ and $\hat {H}^{(1)}_n(x)=xh^{(1)}_n(x)$ are known as the Riccati-Bessel functions [10]. When $n/\lvert k_2a \vert \to +\infty$,

When both $n/\lvert k_1a \rvert \to +\infty$ and $n/\lvert k_2a \rvert \to +\infty$,

The justification for defining the renormalized transmission and reflection coefficients can be easily seen from Eq. (5) because it is always stable for any ordinary medium parameters. According to the duality principle, $R^{\prime TE}_{12}$ and $T^{\prime TE}_{12}$ can be obtained by swapping $\epsilon$ and $\mu$.

2.1.2 Standing wave incidence

Using the same idea, the stable formulas for the case of standing wave incidence can be derived. When $n/\lvert k_1a \rvert \to +\infty$, we have

When $n/\lvert k_2a \rvert \to +\infty$,

When both $n/\lvert k_1a \rvert \to +\infty$ and $n/\lvert k_2a \rvert \to +\infty$,

Applying the duality principle, the formulas of $R_{21}^{\prime TE}$ and $T_{21}^{\prime TE}$ can be obtained by swapping $\epsilon$ and $\mu$.

2.2 Stable formulas of $\alpha$

As shown in Ref. [18], the renormalized and generalized transmission and reflection coefficients contain the parameter $\alpha _i$, which is expressed as

The appearance of $\alpha _i$ is due to the use of the renormalized transmission and reflection coefficients. The arguments of $\alpha _i$ involve $k_ia_i$ and $k_ia_{i-1}$. Since $\lvert k_ia_i \rvert > \lvert k_ia_{i-1} \rvert$, when $n/\lvert k_ia_i \rvert$ satisfies the approximation conditions, $n/\lvert k_ia_{i-1} \rvert$ also does. Hence, when $n/\lvert k_ia_i \rvert \rightarrow +\infty$, we have

When only $n/\lvert k_ia_{i-1} \rvert \rightarrow +\infty$,

To avoid the numerical overflow of the Gamma function with large $n$, by $\Gamma \left (n + \frac {1}{2}\right )=\frac {\left (2n - 1\right )!!}{2^n}\Gamma \left (\frac {1}{2}\right )$,

2.3 Stable formulas of the renormalized total transmission coefficient

The arguments of the renormalized total transmission coefficient involve $k_i$, as well as the outer radius $a_i$ and inner radius $a_{i-1}$ of the $i$th layer. Hence, the fact $a_{i}>a_{i-1}$ is used to distinguish different approximation cases.

2.3.1 Out-going wave incidence

It is shown in Ref. [18] that the renormalized total transmission coefficient from layer 1 to layer $m$ is

When only $n/\lvert k_{i+1}a_i \rvert \rightarrow +\infty$,

2.3.2 Standing wave incidence

It has been shown in Ref. [18] that

When $n/\lvert k_{i-1}a_{i-1}\rvert \rightarrow +\infty$, we have

When only $n/\lvert k_{i-1}a_{i-2}\rvert \rightarrow +\infty$, Eq. (11a) is stable so that it can be calculated directly.

3. Point dipole source in a spherically layered medium

An important application of the spherically layered media theory is to calculate the EM fields of a point dipole source with the background of spherically layered media. When a point dipole source is in the $j$th layer, i.e. $r^\prime \in \mathrm {region}~j$, the Debye potentials can be written in a unified form as

The canonical expressions of $F_n(r,r^\prime )$, as shown in Chapter 3 of Ref. [10], are also numerically instable. The renormalized transmission and reflection coefficients as well as asymptotic formulas can be applied to solve this issue.

3.1 Stable formulas when source and observation are in the same layer

When the source and observation are in the same layer, $F_n(r,r^\prime )$ can be modified as

When $n/\lvert k_{j}a_{j-1} \rvert \rightarrow +\infty$, we have

When $n/\lvert k_{j}r_< \rvert \rightarrow +\infty$,

When $n/\lvert k_{j}r_> \rvert \rightarrow +\infty$,

When $n/\lvert k_{j}a_{j} \rvert \rightarrow +\infty$,

It should be mentioned that the summation of $h^{(1)}_n(k_jr_>)j_n(k_jr_<)$ over $n$ in Eq. (13a) leads to the dipole radiation in free space, which can be computed by using the scalar Green’s function instead of the series summation so that the numerical issue with this term automatically disappears. By doing so, only the scattering fields, i.e. the terms containing $\tilde {R}^\prime$ in Eq. (13), need to be calculated through the series summation, which will improve the computation efficiency. Also, conventional singularity integral techniques can be used in this way when the spherically layered media theory is applied to method of moments.

3.2 Observation is in the $i$th layer and the source is in the $j$th layer with $i>j$

When $r$ is located in the $i$th layer with $i>j$, $F_n(r,r^\prime )$ is modified as

The approximation formulas can be derived similarly with the fact $k_ia_{i-1} < k_i r < k_ia_i$ and $k_ja_{j-1} < k_jr^\prime < k_ja_j$ being noted to distinguish the different cases of approximation.

We define $F_i$ and $F_j$ as

When $n/\lvert k_ia_{i-1} \rvert \rightarrow +\infty$, we have

When $n/\lvert k_ir \rvert \rightarrow +\infty$,

In Eqs. (16a) and (16b), in order to avoid numerical overflow with large $n$, the factor $\frac {\left (k_ia_{i-1}\right )^n}{\left (2n-1\right )!!}$ should be expanded first as is done in the above. When $n/\lvert k_ia_i \rvert \rightarrow +\infty$,

Similarly for $F_j$, when $n/\lvert k_ja_{j-1} \rvert \rightarrow +\infty$, we have

The factor $\left [\cdot \right ]^2$ should be handled as above. When $n/\lvert k_jr^\prime \rvert \rightarrow +\infty$,

When $n/\lvert k_ja_j \rvert \rightarrow +\infty$,

3.3 Observation is in the $i$th layer and the source is in the $j$th layer with $i<j$

For the case that $r$ is in the $i$th layer with $i<j$, $F_n(r,r^\prime )$ is modified as

We define $F_i$ and $F_j$ as

For $F_i$, when $n/\lvert k_ia_{i-1} \rvert \rightarrow +\infty$,

When $n/\lvert k_ir \rvert \rightarrow +\infty$,

When $n/\lvert k_ia_i \rvert \rightarrow +\infty$,

As for $F_j$, when $n/\lvert k_ja_{j-1} \rvert \rightarrow +\infty$, Eq. (19b) is numerically stable. When $n/\lvert k_jr^\prime \rvert \rightarrow +\infty$, we have

When $n/\lvert k_ja_j \rvert \rightarrow +\infty$,

Finally, the first and second derivatives of $F_n(r,r^\prime )$ will be encountered in EM fields calculation [18]. The differentiation can be directly applied on the asymptotic forms of $F_n(r,r^\prime )$. The resultant formulas are much simpler than those shown in [18].

4. Numerical results

Numerical tests have been performed to characterize the derived stable formulas with the results shown below. More numerical results can be found in the supplemental document. The programming was done by using GNU Octave [25].

4.1 Validations with a point dipole source

In this validation, we examine the $H_r$ and $E_r$ components produced by a unit point dipole source, using both asymptotic and conventional formulas. The derivations of the field expressions can be found in the Supplement 1 of this paper or Ref. [18]. This validation is relevant to the dyadic Green’s function that is an important application of the spherically layered media theory. The media distribution is shown in Fig. 1 with the parameters listed in Table 1. The frequency is chosen so that $k_0=1$.

 

Fig. 1. The layered dielectric sphere with a point dipole source for the validation.

Download Full Size | PDF

Table 1. The parameters of Fig. 1

View Table | View all tables in this article

The observation point of $H_r$ is $(r=0.95\,\text {m}, \theta = \pi /2-\pi /20, \phi = \pi /2)$, and the source point is $(r^\prime =1.01\,\text {m}\sim 1.29\,\text {m}, \theta ^\prime =\pi /2, \phi ^\prime =\pi /2)$. The real and imaginary parts of $H_r$ are plotted in Figs. 2(a) and 2(b) respectively. As $r^\prime$ approaches to around 1.1 m, the distance between the source and observation points gets smaller so that larger order Bessel functions must be included in the series summation, which leads to the numerical overflow of the conventional formulas. The asymptotic formulas are more stable than the canonical formulas. In this test, $n/|x|>10$ is the threshold to trigger the use of asymptotic formulas. The error of the asymptotic result relative to the canonical result can be decreased if we increase the threshold.

 

Fig. 2. The $H_r$ and $E_r$ calculated by the asymptotic method and the conventional method, as well as the relative errors. PMCHW result is also shown. (a) $\text {Re}(H_r)$, (b) $\text {Im}(H_r)$, (c) $\text {Re}(E_r)$ and (d) $\text {Im}(E_r)$.

Download Full Size | PDF

In Fig. 2(a), the results obtained with Poggio-Miller-Chang-Harrington-Wu (PMCHW) method are given to validate the results of layered media theory [27]. Since the imaginary part of $H_r$ shown in Fig. 2(b) is much smaller than the real part of $H_r$ shown in Fig. 2(a), the comparison with the PMCHW results makes no sense so that it is not done in Fig. 2(b).

The observation point of $E_r$ is $(r=11.1\,\text {m}, \theta = \pi /2, \phi =\pi /2-\pi /10)$, and the source point is $(r^\prime =7.3\,\text {m}\sim 10.3\,\text {m}, \theta ^\prime =\pi /2, \phi ^\prime =\pi /2)$. It is seen in Fig. 2(c) and 2(d) that as $r^\prime \rightarrow r$, larger order Bessel functions are summed so that the conventional formulas suffer from numerical overflow at $r^\prime =9.3 \text {m}$. However, the asymptotic formulas works normally within the whole parameter range. The relative error of the asymptotic and canonical results is zero in this test. The PMCHW results are in good agreement with the results of the layered media theory.

We have also compared the averaged computational times for the data in Fig. 2 among the asymptotic formulas, the conventional formulas and the method in [18]. In the computation, the iteration of series summation is stopped when the increment is smaller than $1e-6$ of the series sum. It can be seen in Table 2 that the asymptotic formulas takes a little bit longer computation time than the canonical formulas, which may be because the renormalized and asymptotic formulas are a little more complicated than the conventional formulas, and need to evaluate the order-to-argument ratio and do logical judgment in programming in order to trigger the correct asymptotic formulas. However, the proposed method is faster than the method in Ref. [18] in the test. The main reason is that for the method in Ref. [18], the stable iteration of $j_{n-1}(x)/j_n(x)$ should be carried out in the direction of decreasing $n$ so that some preprocessing is required in the algorithm.

Table 2. Comparison of the averaged computational time

View Table | View all tables in this article

4.2 Plane-wave scattering by a layered sphere

Plane-wave scatterings by a layered sphere, as illustrated in Fig. 3(a), is a canonical problem in electromagnetics. When the computation region is large, higher-order Bessel functions have to be included in order to approximate a plane wave, which may lead to numerical instabilities. Here, we calculate the backward scattering of a layered sphere illuminated by an $x$-polarized incident plane wave as illustrated. The frequency is 1 GHz, $\epsilon _r=\mu _r=3$, and the radii of the medium interfaces are 0.1 m, 0.2 m and 0.3 m respectively. The observation points are located on the positive $z$ axis. Figure 3(b) shows the $x$ component of the scattered electric field at different observation points. It can bee seen that as the observation point gets farther, the canonical formulas break down because the order of the Bessel functions gets higher and higher, whereas the asymptotic formulas are still stable. The results obtained with asymptotic formulas agree with the PMCHW results also. The derivation details can be found in the Supplement 1.

 

Fig. 3. (a) The plane wave scatterings by a layered sphere. (b) The amplitude of the scattering $E^s_x$ component.

Download Full Size | PDF

4.3 Schumann resonance frequency estimation

Schumann resonance occurs because the space between the Earth and the conductive ionosphere forms a closed electromagnetic resonance cavity at extremely low frequencies. Observations of Schumann resonances can be used to track global lightning activity, study the lower ionosphere on Earth, predict short-term earthquake, etc. The Schumann resonance frequencies can be approximated as $f_n=c\sqrt {n(n+1)}/2\pi a$, where $c$ is the light speed in vacuum, $a$ is the Earth radius and $n$ is resonance mode number.

By employing the derived stable formulas, we modeled Schumann resonance as shown in Fig. 4(a). The relative permittivity of the ionosphere is $\epsilon _r=1-\left (\frac {\omega _p}{\omega }\right )^2$ with the plasma resonance frequency $\omega _p = e\sqrt \frac {n_e}{m_e \epsilon _0}$, where $e$ is the electron charge, $n_e=10^9\,\text {m}^{-3}$ is the free electron density and $m_e$ is the electron mass. The earth is modeled as a conductor with the conductivity of $10^7$ S/m. In this model, an electric dipole source is put at the height of 10 km from the Earth ground, and the observer is at 35 km from the Earth ground. The calculated strength of $E_r$ component versus frequency by using the asymptotic theory is plotted in Fig. 4(b). As shown, the first two Schumann resonances occur at 10.5 Hz and 18.25 Hz, which are close to 10.59 Hz and 18.35 Hz given by the above formula. The electric field integral equation (EFIE) method was also used for validation. Since the permittivity of the ionosphere and the earth is very large in magnitude, it is difficult to do the meshing. Hence, both the ionosphere and the earth are modeled as perfect electric conductors in the EFIE validation. The good agreement has been achieved as shown.

 

Fig. 4. Schumann resonance frequency estimation. (a) An electric dipole source is put in the model to excite the fields. (b) The amplitude of the $E_r$ component and the maximum order of the Bessel functions in the calculation at each frequency.

Download Full Size | PDF

Figure 4(b) also shows the iteration number, which is equal to the maximum order of the Bessel function to achieve the required calculation precision at each frequency. For a cavity, when the excitation frequency is not equal to the resonance frequency, a large number of modes will be excited; whereas when the excitation frequency is close to the eigen-frequency of one mode, that mode will be excited very much stronger than the other modes. For the typical parameter values of this problem, $n/|x| \gg 1$ so that conventional formulations of the layered medium theory lead to numerical instability and computation failure.

4.4 Luneberg lens design

As a final test, we compared the numerical stabilities of the proposed method and the conventional method by designing a Luneberg lens. Luneberg lens can generate highly directive beam with a dipole source located at the surface of the lens. As shown in Fig. 5(a), the radius of the lens is $R=3$ m, $\mu _r=1$ and $\epsilon _r=2-\left (\frac {r}{R}\right )^2$. A $z$-polarized dipole source is put at $(0,-\text {3 m},0)$. The electric field $E_z$ at $xoy$ plane is computed by using the proposed method, the conventional method and a commercial solver separately, and the results are plotted in Figs. 5(b)$\sim$5(d) for comparison. As shown, the result of the proposed method agrees with that of the commercial solver, whereas the conventional method breaks down (the white region) due to numerical instabilities.

 

Fig. 5. The $E_z$ field at $xoy$ plane of (a) the Luneberg lens calculated by using (b) the proposed method, (c) the conventional method, and (d) a commercial solver.

Download Full Size | PDF

5. Discussions and conclusions

We used the measure $n/|x|$ to trigger the asymptotic formulas for Bessel and Hankel functions because if $n/|x|$ is sufficiently large, then either $n$ is large or $|x|$ is small so that the asymptotic formulas work. However, $n/|x^2|$ could be a better choice since it is a measure of the accuracy of the asymptotic formulas, as seen in the remainder term in Eqs. (1) and (2). In other words, the stable asymptotic formulas are triggered as long as they reach a certain accuracy.

In Table 3, we have compared the accuracies of the asymptotic formulas with the thresholds $n/|x|=10$ and $n/|x^2|=10$ respectively. Note that the case of $n=10$ and $x=1$ applies to both thresholds. It can be seen that for the threshold $n/|x^2|=10$, the relative error remains at the same order, whereas the relative error changes greatly for the threshold $n/|x|=10$. A uniform error is obtained with the threshold $n/|x^2|$.

Table 3. Accuracy comparison of the asymptotic formulas with thresholds $n/|x|=10$ and $n/|x^2|=10$

View Table | View all tables in this article

Under the limiting case that all the Bessel and Hankel functions in the algorithm may numerically underflow or overflow, the stable asymptotic formulas reduce to very simple and numerically stable form with only physical parameters of the problem, such as dielectric constants and geometric parameters, which can be seen from Eqs. (5), (8), (9b), (10b), (11b), (13e) for example.

In conclusion, numerically stable asymptotic formulas of the renormalized spherically layered media theory have been proposed uniformly for the cases of large orders and small arguments. To avoid numerical instabilities under these cases, the Bessel and Hankel functions are first expressed in asymptotic forms, and then the diverging terms are canceled out. The derived asymptotic formulas for the single-interface transmission coefficient $T^\prime$ and reflection coefficient $R^\prime$, the parameter $\alpha$, the renormalized total transmission coefficient $\tilde {T}^\prime$, and the function $F$ are numerically stable. Numerical tests have shown that the derived formulas are good approximations to the conventional formulas, but are more stable when the Bessel (and Hankel) function’s order is large or its argument is small. The numerical issue under large lossy case can be solved similarly or by employing the scaled Bessel and Hankel functions, which will be addressed in another paper.