The electromagnetic wave scatterings by a spherically layered media (SLM) is a canonical theoretical problem in electromagnetics. Its solution has wide applications in science and engineering, such as material and optical physics [1,2], chemistry [3,4], atmospheric physics [5], antenna designs [6], etc. For the case of plane wave incidence, Mie theory (or Lorentz-Mie theory) can be used to solve the scattering problem from a single-layered sphere (or multi-layered sphere) [7–9]. The key point in the Lorenz-Mie theory is to solve for the scattering parameters, which are the coefficients in front of the spherical Bessel and Hankel functions defined in each dielectric layer. The scattering parameters are just the amplitudes of the standing and outgoing spherical waves and are solved directly without using the concepts of reflection and transmission. The disadvantage is that the expression for the scattering parameter is complicated and is only applicable to the case of plane wave incidence.

The SLM theory is similar to the Lorentz-Mie theory but uses the concepts of reflection and transmission to solve electromagnetic propagation and scattering problems in SLM [10]. The SLM theory is relatively simple to compute and can handle arbitrary incident waves. For example, it is used to calculate the dyadic Green’s function for SLM, which is an important tool for solving electromagnetic scattering problems in SLM backgrounds [11–13].

Both Lorentz-Mie theory and SLM theory suffer from numerical instability problems because of the involvement of Bessel and Hankel functions [14]. First, the Hankel function becomes very large and the Bessel function goes to zero when the function’s order goes to infinity. Second, when the argument goes to zero, the Hankel (Bessel) function grows (decreases) at the rate of power functions. Third, when the argument has a large imaginary part, which corresponds to high-loss media, the Bessel function is exponentially large. All of these features may lead to numerical overflows or underflows in numerical computations.

Recently, we have proposed the renormalized spherically layered media (RSLM) theory, where the reflection and transmission coefficients are renormalized so that their magnitudes are kept at ordinary level [15]. Based on that theory, a method that cancels divergent factors in the asymptotic formulas of Bessel functions has been proposed to solve the first and second numerical instability issues mentioned above [16], while the third issue remains unsolved.

This paper addresses the third issue by introducing scaled Bessel functions into RSLM theory. Scaled Bessel functions are Bessel functions with the asymptotic factors under large arguments extracted out so that scaled Bessel functions are not divergent in high-loss case. The computation of scaled Bessel functions has already been supported by open source mathematical libraries such as SLATEC¹ or software such as Matlab and GNU Octave² [17]. The asymptotic factors extracted from the Bessel functions can be canceled out in the RSLM theory so that the whole theory can be computed stably even for large lossy cases. In the following, we will first review the RSLM theory for the convenience of discussion. Then, stable formulas employing scaled Bessel functions will be derived. After that, numerical tests will be shown to validate the derived stable formulas. The time factor $ {\mathrm{e}}^{-\mathrm{i}\omega t} $ is adopted throughout this paper.

Unlike Cartesian coordinates, the unit vectors $ \hat{r} $, $ \hat{\theta } $, and $ \hat{\phi } $ of the spherical coordinates are variable vectors. As a result, the associated electric field and magnetic field components do not satisfy the scalar homogeneous Helmholtz equation. Instead, Debye potentials $ {\pi }_{e} $ and $ {\pi }_{m} $ are commonly used to express TE and TM spherical wave functions in free space respectively [10].

Debye potentials satisfy the Helmholtz equation

$ \left({\nabla }^{2}+{k}^{2}\right)\pi =0 $ (1)

where $ \pi $ denotes either $ {\pi }_{e} $ or $ {\pi }_{m} $, and $ k $ is the wave number. The solution to (1) is of the form

$ \left[a{j}_{n}\right(kr)+b{h}_{n}^{\left(1\right)}(kr\left)\right]{P}_{n}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)[c\mathrm{c}\mathrm{o}\mathrm{s}m\phi +d\mathrm{s}\mathrm{i}\mathrm{n}m\phi ] $ (2)

where $ {j}_{n}\left(kr\right) $ is the spherical Bessel function, $ {h}_{n}^{\left(1\right)}\left(kr\right) $ is the first kind spherical Hankel function, $ {P}_{n}^{m}(\mathrm{c}\mathrm{o}\mathrm{s}\theta ) $ is the associated Legendre function, and $ a\mathrm{,}\ b\mathrm{,\ }c, $ and $ d $ are the linear superposition coefficients. It can be shown that for the $ n $-th spherical harmonic potentials $ {\pi }_{mn} $ and $ {\pi }_{en} $, the corresponding radial field components are

$ {H}_{rn}=\frac{1}{{\mathrm{i}}\omega \mu }\frac{n(n+1)}{r}{\pi }_{mn} $ (3a)

$ {E}_{rn}=-\frac{1}{{\mathrm{i}}\omega \textit{ϵ}}\frac{n(n+1)}{r}{\pi }_{en} $ (3b)

and the transverse to $ r $ components are

$ {\boldsymbol{H}}_{s}=-\boldsymbol{r}\times {\nabla }_{s}{\pi }_{e}+\frac{1}{{\mathrm{i}}\omega \mu }\frac{1}{r}\frac{\partial }{\partial r}{r}^{2}{\nabla }_{s}{\pi }_{m} $ (4a)

$ {\boldsymbol{E}}_{s}=-\boldsymbol{r}\times {\nabla }_{s}{\pi }_{m}-\frac{1}{{\mathrm{i}}\omega \textit{ϵ} }\frac{1}{r}\frac{\partial }{\partial r}{r}^{2}{\nabla }_{s}{\pi }_{e} $ (4b)

where the subscript $ s $ denotes the transverse components, i.e. $ \theta $ and $ \phi $ components. (3) and (4) imply that once the $ n $-th harmonic radial fields are known, the $ n $-th harmonic transverse fields are also known. Solving the scalar wave equation in (1) for the Debye potentials is the key to the problem.

Under spherical coordinates, the incident wave is usually expressed as an out-going wave at the field point of radius larger than the radius of the source point. The wave propagations in a multi-layered sphere can be derived by first considering the propagation properties for the single-layered or single-interface case.

Given an out-going wave $ {A}_{1}{h}_{n}^{\left(1\right)}\left({k}_{1}r\right){P}_{n}^{m}\left(\mathrm{cos}\theta \right){{\rm{e}}}^{{\mathrm{i}}m\phi } $ impinging on a single interface of two different media as shown in Fig. 1, the electric or magnetic total potentials in region 1 ($ r < a $) and region 2 ($ r > a $) can be uniformly written as

$ {\pi }_{1}= {A}_{1}\left[{h}_{n}^{\left(1\right)}\left({k}_{1}r)+{R}_{12}{j}_{n}({k}_{1}r\right)\right]{P}_{n}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right){{\rm{e}}}^{{\mathrm{i}}m \phi }= {A}_{1}\left[{h}_{n}^{\left(1\right)} ({k}_{1}r)+{R}_{12}^\prime\frac{{h}_{n}^{\left(1\right)}\left({k}_{1}a\right)}{{j}_{n}\left({k}_{1}a\right)}{j}_{n}({k}_{1}r )\right] {P}_{n}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right){{\rm{e}}}^{{\mathrm{i}}m\phi } $ (5a)

$ {\pi }_{2}={A}_{1}{T}_{12}{h}_{n}^{\left(1\right)}\left({k}_{2}r\right){P}_{n}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right){{\rm{e}}}^{{\mathrm{i}}m\phi }={A}_{1}{T}_{12}^\prime\frac{{h}_{n}^{\left(1\right)}\left({k}_{1}a\right)}{{h}_{n}^{\left(1\right)}\left({k}_{2}a\right)}{h}_{n}^{\left(1\right)}\left({k}_{2}r\right){P}_{n}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right){{\rm{e}}}^{{\mathrm{i}}m\phi } $ (5b)

where $ {R}_{12} $ ($ {R}_{12}^\prime $) and $ {T}_{12} $ ($ {T}_{12}^{\,\prime} $) are the single-interface (renormalized) reflection and transmission coefficients respectively, and $ {A}_{1} $ is the amplitude of the incident out-going wave. Using the field continuity condition across the interface, it can be shown that

$ R^{{\mathrm{TM}}}_{12}=\frac{\sqrt{\textit{ϵ}_2\mu _1}\hat{H}^{(1)}_{n}(k_2a)\hat{H}{^{(1)}_n}^\prime(k_1a)-\sqrt{\textit{ϵ} _1\mu _2}\hat{H}^{(1)}_{n}(k_1a)\hat{H}{^{(1)}_n}^\prime(k_2a)}{\sqrt{\textit{ϵ} _1\mu _2}\hat{J}_n(k_1a)\hat{H}{^{(1)}_n}^\prime(k_2a)-\sqrt{\textit{ϵ} _2\mu _1}\hat{H}^{(1)}_{n}(k_2a)\hat{J}_n^\prime(k_1a)} $ (6a)

$ T^{{\mathrm{TM}}}_{12}=\frac{i\textit{ϵ} _2\sqrt{\dfrac{\mu _2}{\textit{ϵ} _1} }}{\sqrt{\textit{ϵ} _1\mu _2}\hat{J}_{n}(k_1a)\hat{H}{^{(1)}_n}^\prime(k_2a)-\sqrt{\textit{ϵ} _2\mu _1}\hat{H}^{(1)}_{n}(k_2a)\hat{J}^\prime_n(k_1a)} $ (6b)

for the TM case, where $ {\hat{J}}_{n}\left(x\right)=x{j}_{n}\left(x\right) $ and $ {\hat{H}}_{n}^{\left(1\right)}\left(x\right)=x{h}_{n}^{\left(1\right)}\left(x\right) $ are known as the Riccati-Bessel functions [14], and the superscript prime means derivative with respect to the argument. $ {R}_{12}^{\mathrm{T}\mathrm{E}} $ and $ {T}_{12}^{\mathrm{T}\mathrm{E}} $ can be obtained by the duality principle.

In addition, we have

$ {R}_{12}^\prime=\frac{{j}_{n}\left({k}_{1}a\right)}{{h}_{n}^{\left(1\right)}\left({k}_{1}a\right)}{R}_{12} $ (7a)

$ {T}_{12}^{\,\prime}=\frac{{h}_{n}^{\left(1\right)}\left({k}_{2}a\right)}{{h}_{n}^{\left(1\right)}\left({k}_{1}a\right)}{T}_{12} $ (7b)

where the superscript prime indicates that the quantity is renormalized. Note that from (5), $ {R}_{12}^\prime $ and $ {T}_{12}^{\,\prime} $ are the amplitude ratios of the reflected and transmitted waves to the incident wave at the interface, so that they have ordinary magnitudes. However, it can be seen from (7) that since the renormalized reflection and transmission coefficients have ordinary magnitudes, the conventional definition of reflection and transmission coefficients could lead to too large or too small values if the Bessel functions in (7) have high orders, small arguments, or large loss.

With the single-interface reflection and transmission coefficients, we can derive the generalized reflection and transmission coefficients for a spherically multi-layered medium. It can be shown that the generalized and renormalized reflection and transmission coefficients are [15,16]

$ {\tilde{R}}_{i{\mathrm{,}}i+1}^\prime={R}_{i{\mathrm{,}}i+1}^\prime+\frac{{T}_{i+1{\mathrm{,}}i}^\prime{\tilde{R}}_{i+1{\mathrm{,}}i+2}^\prime{T}_{i{\mathrm{,}}i+1}^{\,\prime}{\alpha }_{i+1}}{1-{R}_{i+1{\mathrm{,}}i}^\prime{\tilde{R}}_{i+1{\mathrm{,}}i+2}^\prime{\alpha }_{i+1}} $ (8a)

$ {S}_{i{\mathrm{,}}i+1}^{\,\prime}=\frac{{T}_{i{\mathrm{,}}i+1}^{\,\prime}}{1-{R}_{i+1{\mathrm{,}}i}^\prime{\tilde{R}}_{i+1{\mathrm{,}}i+2}^\prime{\alpha }_{i+1}} $ (8b)

where the parameter $ {\alpha }_{i} $ is defined as

$ {\alpha }_{i}\equiv \frac{{j}_{n}\left({k}_{i}{a}_{i-1}\right)}{{j}_{n}\left({k}_{i}{a}_{i}\right)}\frac{{h}_{n}^{\left(1\right)}\left({k}_{i}{a}_{i}\right)}{{h}_{n}^{\left(1\right)}\left({k}_{i}{a}_{i-1}\right)} $ (9)

which is only relevant to the $ i $-th layer’s parameter. The superscript prime in (8) means that the quantity is renormalized.

Thus, the renormalized total transmission coefficient, that is the amplitude ratio of the out-going wave leaving the interface of radius $ {a}_{m-1} $ and the out-going wave entering the interface of radius $ {a}_{1} $, is defined and calculated as

$ {\tilde{T}}_{1m}^{\,\prime}\equiv \frac{{A}_{m}{h}_{n}^{\left(1\right)}\left({k}_{m}{a}_{m-1}\right)}{{A}_{1}{h}_{n}^{\left(1\right)}\left({k}_{1}{a}_{1}\right)}=\left[\prod _{i=1}^{m-2}\frac{{h}_{n}^{\left(1\right)}\left({k}_{i+1}{a}_{i+1}\right)}{{h}_{n}^{\left(1\right)}\left({k}_{i+1}{a}_{i}\right)}{S}_{i{\mathrm{,}}i+1}^{\,\prime}\right]{S}_{m-1{\mathrm{,}}m}^{\,\prime} $ (10a)

when $ 3\le m\le N $, and

$ {\tilde{T}}_{12}^{\,\prime}={S}_{12}^{\,\prime} $ (10b)

when $ m=2 $.

Under spherical coordinates, the incident wave is usually expressed as a standing wave at the field point of radius smaller than the radius of the source point. Similarly, the wave propagations in a multi-layered sphere can be derived by first considering the propagation properties of single-interface case.

Given a standing wave $ {A}_{2}{j}_{n}\left({k}_{2}r\right){P}_{n}^{m}\left(\mathrm{cos}\theta \right){{\rm{e}}}^{{\mathrm{i}}m\phi } $ impinging on a single interface as illustrated in Fig. 2, the total Debye potentials in each region can be expressed as

$ {\pi }_{1}={A}_{2}{T}_{21}{j}_{n}\left({k}_{1}r\right){P}_{n}^{m}(\mathrm{c}\mathrm{o}\mathrm{s}\theta ){{\rm{e}}}^{{\mathrm{i}}m\phi }={A}_{2}{T}_{21}^{\,\prime}\frac{{j}_{n}\left({k}_{2}a\right)}{{j}_{n}\left({k}_{1}a\right)}{j}_{n}\left({k}_{1}r\right){P}_{n}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right){{\rm{e}}}^{{\mathrm{i}}m\phi } $ (11a)

$ \begin{split} \pi_2 =& A_2 \left[j_n(k_{2}r)+R_{21} h_n^{(1)}(k_{2}r)\right]P_n^m(\cos\theta){\rm{e}}^{{\mathrm{i}}m\phi} =\\ & A_2\left[j_n(k_{2}r) + R_{21}^\prime\frac{j_{n}(k_{2}a)}{h_{n}^{(1)}(k_{2}a)}h_{n}^{(1)}(k_{2}r)\right]P_n^m(\cos\theta){\rm{e}}^{{\mathrm{i}}m\phi} \end{split} $ (11b)

For TE case,

$ R^{\rm TE}_{21}=\frac{\sqrt{\textit{ϵ} _1\mu _2}\hat{J}_n(k_2a)\hat{J}^\prime_n(k_1a)-\sqrt{\textit{ϵ} _2\mu _1}\hat{J}_n(k_1a)\hat{J}^\prime_n(k_2a)}{\sqrt{\textit{ϵ} _2\mu _1}\hat{J}_n(k_1a)\hat{H}{^{(1)}_n}^\prime(k_2a)-\sqrt{\textit{ϵ} _1\mu _2}\hat{H}^{(1)}_n(k_2a)\hat{J}^\prime_n(k_1a)} $ (12a)

$ T^{\rm TE}_{21}=\frac{i\mu _2\sqrt{\dfrac{\textit{ϵ} _1}{\mu_2} }}{\sqrt{\textit{ϵ} _2\mu _1}\hat{J}_{n}(k_1a)\hat{H}{^{(1)}_n}^\prime(k_2a)-\sqrt{\textit{ϵ} _1\mu _2}\hat{H}^{(1)}_{n}(k_2a)\hat{J}^{\prime}_{n}(k_1a)} $ (12b)

By the duality principle, $ {R}_{21}^{\mathrm{T}\mathrm{M}} $ and $ {T}_{21}^{\mathrm{T}\mathrm{M}} $ can be obtained. Also, we have

$ {R}_{21}^\prime=\frac{{h}_{n}^{\left(1\right)}\left({k}_{2}a\right)}{{j}_{n}\left({k}_{2}a\right)}{R}_{21} $ (13a)

$ {T}_{21}^\prime=\frac{{j}_{n}\left({k}_{1}a\right)}{{j}_{n}\left({k}_{2}a\right)}{T}_{21} $ (13b)

Following the same idea, we can obtain the renormalized and generalized coefficients’ expressions as

$ {\tilde{R}}_{i{\mathrm{,}}i-1}^\prime={R}_{i{\mathrm{,}}i-1}^\prime+\frac{{T}_{i-1{\mathrm{,}}i}^\prime{\tilde{R}}_{i-1{\mathrm{,}}i-2}^\prime{T}_{i{\mathrm{,}}i-1}^\prime{\alpha }_{i-1}}{1-{R}_{i-1{\mathrm{,}}i}^\prime{\tilde{R}}_{i-1{\mathrm{,}}i-2}^\prime{\alpha }_{i-1}} $ (14a)

$ {S}_{i{\mathrm{,}}i-1}^\prime=\frac{{T}_{i{\mathrm{,}}i-1}^\prime}{1-{R}_{i-1{\mathrm{,}}i}^\prime{\tilde{R}}_{i-1{\mathrm{,}}i-2}^\prime{\alpha }_{i-1}} $ (14b)

Thus, the renormalized total transmission coefficient, that is the amplitude ratio of the standing wave leaving the interface of radius $ {a}_{m} $ and the standing wave entering the interface of radius $ {a}_{N-1} $, is defined and calculated as

$ {\tilde{T}}_{Nm}^{\,\prime}\equiv \frac{{A}_{m}{j}_{n}\left({k}_{m}{a}_{m}\right)}{{A}_{N}{j}_{n}\left({k}_{N}{a}_{N-1}\right)}=\left[\prod _{i=N}^{m+2}\frac{{j}_{n}\left({k}_{i-1}{a}_{i-2}\right)}{{j}_{n}\left({k}_{i-1}{a}_{i-1}\right)}{S}_{i{\mathrm{,}}i-1}^{\,\prime}\right]{S}_{m+1{\mathrm{,}}m}^{\,\prime} $ (15a)

when $ 1\le m\le N-2 $, and

$ {\tilde{T}}_{N{\mathrm{,}}N-1}^{\,\prime}={S}_{N{\mathrm{,}}N-1}^{\,\prime} $ (15b)

when $ m=N-1 $.

Although the renormalized reflection and transmission coefficients in the RSLM theory reviewed above have ordinary values, the computation of RSLM is still unstable due to the Bessel functions involved in the formulas. When the media are high lossy or the radius of a lossy layer is quite large, then the arguments of Bessel functions in the theory will have a large imaginary part. It is well known that when $ {x''}=\mathfrak{I}\left(x\right)\to +{\infty } $,

$ {j}_{n}\left(x\right)\approx \frac{1}{x}\mathrm{c}\mathrm{o}\mathrm{s}\left(x-\frac{n+1}{2}\pi \right)\sim \frac{{\rm e}^{{x''}}}{x} $ (16a)

$ {h}_{n}^{\left(1\right)}\left(x\right)\approx \frac{1}{x}{\rm e}^{{\mathrm{i}}\left(x-\tfrac{n\pi }{2}-\tfrac{\pi }{2}\right)}\sim \frac{{\rm e}^{-{x''}}}{x} $ (16b)

This asymptotic behavior will lead to the divergence of $ {j}_{n}\left(x\right) $ or division by zero problem if $ {h}_{n}^{\left(1\right)}\left(x\right) $ is in the denominator. In order to avoid the tough asymptotic properties of Bessel functions, scaled Bessel functions are defined by extracting the asymptotic factors out of the Bessel functions as

$ {\mathring{h}}_{n}^{\left(1\right)}\left(x\right)\equiv {h}_{n}^{\left(1\right)}\left(x\right){\rm e}^{-{\mathrm{i}}x} $ (17a)

$ {\mathring{j}}_{n}\left(x\right)\equiv {j}_{n}\left(x\right){\rm e}^{-\left|{x''}\right|} $ (17b)

Hence, if the Bessel functions in RSLM theory are written as the product of scaled Bessel functions and asymptotic factors, then the asymptotic factors can be canceled out whereas the scaled Bessel functions have good numerical behavior under large lossy case.

Single-interface reflection and transmission coefficients are the building blocks of SLM or RSLM theory. Hence, we start with deriving the high lossy stable formulas for single-interface coefficients, and then consider the multi-interface case.

Applying the scaled Bessel functions to the renormalized reflection and transmission coefficients in (6) and (7), we get

$ {R'}_{12}^{\mathrm{T}\mathrm{M}}=\frac{{\mathring{j}}_{n}\left({k}_{1}a\right)}{{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{1}a\right)}\frac{\sqrt{{\textit{ϵ} }_{2}{\mu }_{1}}{k}_{2}a{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{2}a\right){B}_{n}\left({k}_{1}a\right)-\sqrt{{\textit{ϵ} }_{1}{\mu }_{2}}{k}_{1}a{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{1}a\right){B}_{n}\left({k}_{2}a\right)}{\sqrt{{\textit{ϵ} }_{1}{\mu }_{2}}{k}_{1}a{\mathring{j}}_{n}\left({k}_{1}a\right){B}_{n}\left({k}_{2}a\right)-\sqrt{{\textit{ϵ} }_{2}{\mu }_{1}}{k}_{2}a{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{2}a\right){A}_{n}\left({k}_{1}a\right)} $ (18a)

$ {T\,'}_{12}^{\mathrm{T}\mathrm{M}}=\frac{{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{2}a\right)}{{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{1}a\right)}\frac{i{\textit{ϵ} }_{2}\sqrt{\dfrac{{\mu }_{2}}{{\textit{ϵ} }_{1}}}}{\sqrt{{\textit{ϵ} }_{1}{\mu }_{2}}{k}_{1}a{\mathring{j}}_{n}\left({k}_{1}a\right){B}_{n}\left({k}_{2}a\right)-\sqrt{{\textit{ϵ} }_{2}{\mu }_{1}}{k}_{2}a{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{2}a\right){A}_{n}\left({k}_{1}a\right)}{\rm e}^{-({\mathrm{i}}{k}_{1}+|{k}_{1}''\left|\right)a} $ (18b)

where

$ {A}_{n}\left(x\right)={\mathring{j}}_{n}\left(x\right)-x{\mathring{j}}_{n+1}\left(x\right)+n{\mathring{j}}_{n}\left(x\right) $ (19a)

$ {B}_{n}\left(x\right)={\mathring{h}}_{n}^{\left(1\right)}\left(x\right)-x{\mathring{h}}_{n+1}^{\left(1\right)}\left(x\right)+n{\mathring{h}}_{n}^{\left(1\right)}\left(x\right) $ (19b)

The above formulas are stable in high lossy case since only scaled Bessel functions are employed and $ |\mathrm{e}\mathrm{x}\mathrm{p}[-(i{k}_{1}+|{k}_{1}''\left|\right)a\left]\right| $ is equal to unity since $ {k}_{1}'' > 0 $ in lossy case.

The generalized transmission and reflection coefficients in (8) are functions of the single-interface transmission and reflection coefficients as well as the parameter $ {\alpha }_{i} $. The stable formulas of $ {\alpha }_{i} $ can be easily shown as

$ {\alpha }_{i}=\frac{{\mathring{j}}_{n}\left({k}_{i}{a}_{i-1}\right)}{{\mathring{j}}_{n}\left({k}_{i}{a}_{i}\right)}\frac{{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{i}{a}_{i}\right)}{{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{i}{a}_{i-1}\right)}{\rm e}^{({\mathrm{i}}{k}_{i}-|{k}_{i}''\left|\right)({a}_{i}-{a}_{i-1})} $ (20)

Note that the modulus of the exponential factor is smaller than 1 in lossy case so that the formula is numerically stable.

Applying the scaled Bessel functions to the total transmission coefficient in (10a) yields

$ {\tilde{T}}_{1m}^{\,\prime}=\left[\prod _{i=1}^{m-2}\frac{{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{i+1}{a}_{i+1}\right)}{{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{i+1}{a}_{i}\right)}{\rm e}^{{\mathrm{i}}{k}_{i+1}\left({a}_{i+1}-{a}_{i}\right)}{S}_{i{\mathrm{,}}i+1}^{\,\prime}\right]{S}_{m-1{\mathrm{,}}m}^{\,\prime} $ (21)

for $ 3\le m\le N $, which is stable under high lossy case obviously.

Similar derivations can be done for the standing wave incidence case. The stable expressions for the transmission and reflection coefficients understanding wave incidence case are

$ {R'}_{21}^{\mathrm{T}\mathrm{M}}=\frac{{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{2}a\right)}{{\mathring{j}}_{n}\left({k}_{2}a\right)}\frac{\sqrt{{\varepsilon }_{2}{\mu }_{1}}{k}_{2}a{\mathring{j}}_{n}\left({k}_{2}a\right){A}_{n}\left({k}_{1}a\right)-\sqrt{{\varepsilon }_{1}{\mu }_{2}}{k}_{1}a{\mathring{j}}_{n}\left({k}_{1}a\right){A}_{n}\left({k}_{2}a\right)}{\sqrt{{\varepsilon}_{1}{\mu }_{2}}{k}_{1}a{\mathring{j}}_{n}\left({k}_{1}a\right){B}_{n}\left({k}_{2}a\right)-\sqrt{{\varepsilon }_{2}{\mu }_{1}}{k}_{2}a{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{2}a\right){A}_{n}\left({k}_{1}a\right)} $ (22a)

$ {T\,'}_{21}^{\mathrm{T}\mathrm{M}}=\frac{{\mathring{j}}_{n}\left({k}_{1}a\right)}{{\mathring{j}}_{n}\left({k}_{2}a\right)}\frac{i{\varepsilon }_{1}\sqrt{\dfrac{{\mu }_{1}}{{\varepsilon }_{2}}}}{\sqrt{{\varepsilon }_{1}{\mu }_{2}}{k}_{1}a{\mathring{j}}_{n}\left({k}_{1}a\right){B}_{n}\left({k}_{2}a\right)-\sqrt{{\varepsilon }_{2}{\mu }_{1}}{k}_{2}a{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{2}a\right){A}_{n}\left({k}_{1}a\right)}{{\mathrm{e}}}^{-({\mathrm{i}}{k}_{2}+|{k}_{2}''\left|\right)a} $ (22b)

The stable formula of the total transmission coefficient in (15a) is

$ {\tilde{T}\,'}_{ Nm}=\left[\prod _{i=N}^{m+2}\frac{{\mathring{j}}_{n}\left({k}_{i-1}{a}_{i-2}\right)}{{\mathring{j}}_{n}\left({k}_{i-1}{a}_{i-1}\right)}{{\mathrm{e}}}^{\left|{k}_{i-1}''\right|\left({a}_{i-2}-{a}_{i-1}\right)}{S}_{i{\mathrm{,}}i-1}^{\,\prime}\right]{S}_{m+1{\mathrm{,}}m}^{\,\prime} $ (23)

for $ 1\le m\le N-2 $.

An important application of SLM theory is to compute the dyadic Green’s function, which is the response of a point dipole source in SLM background. For a point dipole source in the $ j $-th layer, the Debye potentials can be written in an unified form as

$ \pi ={D}_{j}^\prime\sum _{n=1}^{+ {\infty }}{F}_{n}\left(r{\mathrm{,}}\;{r}^\prime\right)\frac{{A}_{n}\left(\theta {\mathrm{,}}\;\phi ;{\theta }^\prime{\mathrm{,}}{\phi }^\prime\right)}{n\left(n+1\right)} $ (24)

where $ {D}_{j}^\prime=-\omega {\mu }_{j}{k}_{j}{\hat{\alpha }}^\prime\cdot {\nabla }^\prime\times {\hat{r}}^\prime{r}^\prime $ corresponds to $ {\pi }_{m} $, and $ {D}_{j}^\prime=i{k}_{j}{\hat{\alpha }}^\prime\cdot {\nabla }^\prime\times {\nabla }^\prime\times {\hat{r}}^\prime{r}^\prime $ corresponds to $ {\pi }_{e} $. The sign ‘$ \mathrm{\text{'}} $’denotes the coordinates related to the source. $ {\hat{\alpha }}^\prime $ denotes the polarization direction of the point dipole. The spherical harmonic function

$ \begin{split}{A}_{n}\left(\theta {\mathrm{,}}\;\phi ;{\theta }^\prime{\mathrm{,}}\;{\phi }^\prime\right) = & \displaystyle\sum _{m=-n}^{n}{Y}_{nm}\left(\theta {\mathrm{,}}\;\phi \right){Y}_{nm}^{*}\left({\theta }^\prime{\mathrm{,}}\;{\phi }^\prime\right)=\sum _{m=-n}^{n}\frac{\left(n-m\right)!}{\left(n+m\right)!}\frac{2n+1}{4\pi }{P}_{n}^{m}\left(\mathrm{cos}\theta \right){{\rm{e}}}^{{\mathrm{i}}m\phi }{P}_{n}^{m}\left(\mathrm{cos}{\theta }^\prime\right){{\rm{e}}}^{-{\mathrm{i}}m{\phi }^\prime} \\ = & \sum _{m=0}^{n}{\varepsilon }_{m}\frac{\left(n-m\right)!}{\left(n+m\right)!}\frac{2n+1}{4\pi }{P}_{n}^{m}\left(\mathrm{cos}\theta \right){P}_{n}^{m}\left(\mathrm{cos}{\theta }^\prime\right)\mathrm{cos}\left(\phi -{\phi }^\prime\right) \end{split} $ (25)

where $ {\varepsilon }_{m}=1 $ when $ m=0 $ and $ {\varepsilon }_{m}=2 $ when $ m > 0 $.

The function $ {F}_{n}(r{\mathrm{,}}\;{r}^\prime) $ describes wave propagations along the radial direction [10], which is also numerically unstable in the high loss case. By applying the scaled Bessel functions to it, stable formulations can be obtained. When both source and observation are in the $ j $-th layer, we have

$ \begin{split} {F}_{n}(r{\mathrm{,}}\;{r}^\prime)= & \left[{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{j}{r}_{ > }\right)+{\tilde{R}}_{j{\mathrm{,}}\;j+1}^\prime\frac{{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{j}{a}_{j}\right)}{{\mathring{j}}_{n}\left({k}_{j}{a}_{j}\right)}{\mathring{j}}_{n}\left({k}_{j}{r}_{ > }\right){e}^{\left(i{k}_{j}-|{k}_{j}''|\right)\left({a}_{j}-{r}_{ > }\right)}\right] \times \\ & \left[{\mathring{j}}_{n}\left({k}_{j}{r}_{ < }\right)+{\tilde{R}}_{j{\mathrm{,}}\;j-1}^\prime\frac{{\mathring{j}}_{n}\left({k}_{j}{a}_{j-1}\right)}{{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{j}{a}_{j-1}\right)}{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{j}{r}_{ < }\right){e}^{\left(i{k}_{j}-|{k}_{j}''|\right)\left({r}_{ < }-{a}_{j-1}\right)}\right]{\tilde{M}}_{j}^\prime{e}^{\left(i{k}_{j}{r}_{ > }+|{k}_{j}''|{r}_{ < }\right)} \end{split} $ (26a)

where $ {\tilde{M}}_{j}^\prime=(1-{\tilde{R}}_{j{\mathrm{,}}\;j-1}^\prime{\tilde{R}}_{j{\mathrm{,}}\;j+1}^\prime{\alpha }_{j}{)}^{-1} $ and $ {r}_{ > ( < )} $ denotes the lager (smaller) one between $ r $ and $ {r}^\prime $. When source is in the $ j $-th layer and observation is in the $ i $-th layer, for the case of $ i > j $ we have

$ \begin{split} {F}_{n}(r{\mathrm{,}}\;{r}^\prime)= & \left[{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{i}r\right)+{\tilde{R}}_{i{\mathrm{,}}\;i+1}^\prime\frac{{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{i}{a}_{i}\right)}{{\mathring{j}}_{n}\left({k}_{i}{a}_{i}\right)}{\mathring{j}}_{n}\left({k}_{i}r\right){\rm e}^{\left({\mathrm{i}}{k}_{i}-|{k}_{i}''|\right)({a}_{i}-r)}\right]\times \\ & \left[{\mathring{j}}_{n}\left({k}_{j}{r}^\prime \right)+{\tilde{R}}_{j{\mathrm{,}}\;j-1}^\prime\frac{{\mathring{j}}_{n}\left({k}_{j}{a}_{j-1}\right)}{{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{j}{a}_{j-1}\right)}{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{j}{r}^\prime\right){\rm e}^{\left({\mathrm{i}}{k}_{j}-|{k}_{j}''|\right)\left({r}^\prime-{a}_{j-1}\right)}\right]{\tilde{M}}_{j}^\prime{\tilde{T}}_{ji}^\prime\frac{{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{j}{a}_{j}\right)}{{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{i}{a}_{i-1}\right)}{\rm e}^{\left[{\mathrm{i}}{k}_{i}\left(r-{a}_{i-1}\right)+{\mathrm{i}}{k}_{j}{a}_{j}+|{k}_{j}''|{r}^\prime\right]} \end{split} $ (26b)

and for the case of $ i < j $ we have

$ \begin{split}{F}_{n}(r{\mathrm{,}}\;{r}^\prime)= & \left[{\mathring{j}}_{n}\left({k}_{i}r\right)+{\tilde{R}}_{i{\mathrm{,}}\;i-1}^\prime\frac{{\mathring{j}}_{n}\left({k}_{i}{a}_{i-1}\right)}{{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{i}{a}_{i-1}\right)}{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{i}r\right){\rm e}^{\left({\mathrm{i}}{k}_{i}-|{k}_{i}''|\right)\left(r-{a}_{i-1}\right)}\right]\times \\ &\qquad \left[{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{j}{r}^\prime \right)+{\tilde{R}}_{j{\mathrm{,}}\;j+1}^\prime\frac{{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{j}{a}_{j}\right)}{{\mathring{j}}_{n}\left({k}_{j}{a}_{j}\right)}{\mathring{j}}_{n}\left({k}_{j}{r}^\prime\right){\rm e}^{\left({\mathrm{i}}{k}_{j}-|{k}_{j}'' |\right)\left({a}_{j}-{r}^\prime\right)}\right]{\tilde{M}}_{j}^\prime{\tilde{T}}_{ji}^\prime\frac{{\mathring{j}}_{n}\left({k}_{j}{a}_{j-1}\right)}{{\mathring{j}}_{n}\left({k}_{i}{a}_{i}\right)}{\rm e}^{\left[ |{k}_{i}''|\left(r-{a}_{i}\right)+{\mathrm{i}}{k}_{j}{r}^\prime+|{k}_{j}''|{a}_{j-1}\right]} \end{split} $ (26c)

All terms in (26), including the exponential functions, are not divergent with high lossy media.

The plane wave scatterings of a layered sphere is also a very common problem. It will be shown in this section that, under the high lossy case, its computation is not stable if no special treatment is adopted, whereas stable computation can be done with the use of scaled Bessel functions.

Assuming an incident plane wave

$ {\boldsymbol{E}}^{i}=\hat{x}{E}_{x}=\hat{x}{\rm e}^{-{\mathrm{i}}{k}_{0}z} $ (27)

illuminating a layered sphere, as shown in Fig. 3. The $ r $ component of the incident plane wave is

$ {E}_{r}^{i} = {\rm e}^{-{\mathrm{i}}{k}_{0}r\mathrm{c}\mathrm{o}\mathrm{s}\theta }\mathrm{s}\mathrm{i}\mathrm{n}\theta \mathrm{c}\mathrm{o}\mathrm{s}\phi =\frac{1}{{\mathrm{i}}{k}_{0}r}\frac{\partial }{\partial \theta }{\rm e}^{-{\mathrm{i}}{k}_{0}r\mathrm{c}\mathrm{o}\mathrm{s}\theta }\mathrm{c}\mathrm{o}\mathrm{s}\phi $ (28)

where $ \theta $ and $\phi $ obey the standard definitions in spherical coordinates. According to the plane-wave-to-spherical-wave transformation

$ {{\mathrm{e}}}^{-{\mathrm{i}}{k}_{0}r\mathrm{c}\mathrm{o}\mathrm{s}\theta }=\sum _{n=0}^{+ {\infty }}(-i{)}^{n}\left(2n+1\right){j}_{n}\left({k}_{0}r\right){P}_{n}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right) $ (29)

we have

$ {E}_{r}^{i}=\frac{1}{{\mathrm{i}}{k}_{0}r}\sum _{n=1}^{+ {\infty }}(-i{)}^{n}\left(2n+1\right){j}_{n}\left({k}_{0}r\right){P}_{n}^{1}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)\mathrm{c}\mathrm{o}\mathrm{s}\phi $ (30)

Therefore, the electric Debye potential of the incident plane wave is

$ {\pi }_{en}=-\frac{{\mathrm{i}}\omega {\textit{ϵ} }_{0}r}{n(n+1)}{E}_{rn}=-\frac{\omega {\textit{ϵ} }_{0}}{{k}_{0}}(-i{)}^{n}\frac{2n+1}{n(n+1)}{j}_{n}\left({k}_{0}r\right){P}_{n}^{1}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)\mathrm{c}\mathrm{o}\mathrm{s}\phi $ (31a)

Similarly, it can be shown that the magnetic Debye potential of the incident plane wave is

$ {\pi }_{mn}=-(-i{)}^{n}\frac{2n+1}{n(n+1)}{j}_{n}\left({k}_{0}r\right){P}_{n}^{1}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)\mathrm{s}\mathrm{i}\mathrm{n}\phi $ (31b)

When such a plane wave is incident onto a layered sphere, the Debye potentials in the $ i $-th layer are

$ {\pi }_{en}=-\frac{\omega {\textit{ϵ} }_{0}}{{k}_{0}}(-{\mathrm{i}}{)}^{n}\frac{2n+1}{n(n+1)}{F}_{n}^{e}\left(r\right){P}_{n}^{1}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)\mathrm{c}\mathrm{o}\mathrm{s}\phi $ (32a)

$ {\pi }_{mn}=-(-{\mathrm{i}}{)}^{n}\frac{2n+1}{n(n+1)}{F}_{n}^{m}\left(r\right){P}_{n}^{1}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)\mathrm{s}\mathrm{i}\mathrm{n}\phi $ (32b)

where

$ {F}_{n}^{e\left(m\right)}\left(r\right)=\left[{j}_{n}\left({k}_{i}r\right)+{\tilde{R}}_{i{\mathrm{,}}i-1}^{\mathrm{T}\mathrm{M}\left(\mathrm{T}\mathrm{E}\right)}{h}_{n}^{\left(1\right)}\left({k}_{i}r\right)\right]{\tilde{T}}_{Ni}^{e\left(m\right)} $ (33)

In order to compute the problem stably with large lossy media, the function $ {F}_{n}\left(r\right) $ should be first written in the renormalized form as

$ {F}_{n}\left(r\right)=\left[{j}_{n}\left({k}_{i}r\right)+{\tilde{R}}_{i{\mathrm{,}}\; i-1}^\prime\frac{{j}_{n}\left({k}_{i}{a}_{i-1}\right)}{{h}_{n}^{\left(1\right)}\left({k}_{i}{a}_{i-1}\right)}{h}_{n}^{\left(1\right)}\left({k}_{i}r\right)\right]{\tilde{T}}_{Ni}^\prime\frac{{j}_{n}\left({k}_{N}{a}_{N-1}\right)}{{j}_{n}\left({k}_{i}{a}_{i}\right)} $ (34)

where the renormalized reflection and transmission coefficients are used. Then, to avoid numerical instability issues under the high lossy case, we apply the scaled Bessel function to (34) to get

$ {F}_{n}\left(r\right)=\left[{\mathring{j}}_{n}\left({k}_{i}r\right)+{\tilde{R}}_{i{\mathrm{,}}\; i-1}^\prime\frac{{\mathring{j}}_{n}\left({k}_{i}{a}_{i-1}\right)}{{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{i}{a}_{i-1}\right)}{\mathring{h}}_{n}^{\left(1\right)}\left({k}_{i}r\right){\rm e}^{\left({\mathrm{i}}{k}_{i}-|{k}_{i}''|\right)\left(r-{a}_{i-1}\right)}\right]{\tilde{T}}_{Ni}^\prime\frac{{\mathring{j}}_{n}\left({k}_{N}{a}_{N-1}\right)}{{\mathring{j}}_{n}\left({k}_{i}{a}_{i}\right)}{\rm e}^{\left[|{k}_{N}''|{a}_{N-1}+|{k}_{i}''|\left(r-{a}_{i}\right)\right]} $ (35)

In order to verify the numerical stability of the proposed method, numerical tests have been carried out. We first consider the case of $ \phi $-polarized point dipole source excitation. Fig. 4 illustrates the model in the test, with the model parameters listed in Table 1. The observation point is located at $ (r=10.9{\mathrm{,}}\;\theta =\pi /2{\mathrm{,}}\; \phi =0.4\pi ) $, and the $ \phi $-polarized dipole source (marked by the blue arrow) is located at $ ({r}^\prime=11.5{\mathrm{,}}\;{\theta }^\prime=\pi /2{\mathrm{,}}\;{\phi }^\prime=\pi /2) $. The frequency is chosen so that $ {k}_{0}=1 $. Fig. 5(a) shows the numerical results of the magnetic field $ {H}_{r} $ as a function of medium loss $ {x}'' $ of layer 2. As seen, the conventional formulation breaks down at around $ {x}'' =63 $ when the medium loss of layer 2 increases, while the proposed formulation works stably.

As mentioned above, even though the media loss is not large, a numerical breakdown may still happen if the layers’ radii are very large, so that the Bessel functions’ arguments have a large imaginary part. Fig. 5(b) compares the stability of conventional and proposed formulations when the radii of the layers are very large but with the ordinary level of media loss. The observation point is located at $ (r={a}_{2}-10{\mathrm{,}}\; \theta =\pi /2{\mathrm{,}}\; \phi =0.4\pi ) $, and the source point is located at $ ({r}^\prime={a}_{2}+10{\mathrm{,}}\; {\theta }^\prime=\pi /2{\mathrm{,}}\; {\phi }^\prime=\pi /2) $. We set the frequency so that $ {k}_{0}=1 $. The model parameters for this test are listed in Table 2. As seen in Fig. 5(b), when the radius of layer 2 increases, the conventional method suffers from numerical breakdown at around a₂ = 350.8 m, whereas the proposed method is quite stable.

Next, we performed the numerical test with a $ r $-polarized dipole as shown in Fig. 6. The observation point is located at $ (r=16{\mathrm{,}} \; \theta =\pi /2{\mathrm{,}} \phi =0.4\pi ) $, and the $ r $-polarized dipole source (marked by the blue arrow) is located at $ ({r}^\prime=10{\mathrm{,}}\;{\theta }^\prime=\pi /2{\mathrm{,}}\;{\phi }^\prime=\pi /2) $. The frequency is $ f=1 $ GHz. The other parameters can be found in Table 1. Fig. 7(a) shows the numerical results of the electric field $ {E}_{r} $ as a function of the medium loss $ {x}'' $ of layer 2. It can be seen that as the loss of the medium increases, the traditional formulation becomes unstable and breaks down at $ {x}'' =3.04 $, while the formulation proposed in this paper can still work stably.

Fig. 7(b) compares the stability of the conventional and proposed formulations when the radii of the layers are very large but with ordinary loss. The observation point is located at $ (r=1100{\mathrm{,}}\;\theta =\pi /2{\mathrm{,}}\;\phi =0.4\pi ) $, and the $ r $-polarized dipole source (marked by the blue arrow) is located at $ ({r}^\prime=300{\mathrm{,}}\;{\theta }^\prime=\pi /2{\mathrm{,}}\;{\phi }^\prime=\pi /2) $. The frequency is chosen so that $ {k}_{0}=1 $. The other parameters are given in Table 2. As shown in the figure, as the radius $ {a}_{2} $ increases, the proposed method works stably while the conventional method breaks down around $ {a}_{2}=350.8 $ m.

The performance of the conventional method and the proposed method with plane wave incidence were also compared. The model parameters for this test are given in Table 3. Fig. 8 shows the numerically computed $ \left|{E}_{r}\right| $ in the case of a plane wave incidence along the $ -z $ direction. The frequency is 1 GHz. The observation point is located at $ (r=4{\mathrm{,}}\;\theta =\pi /2{\mathrm{,}}\;\phi =0) $. In Fig. 8(a), $ {x}''\in (8{\mathrm{,}}\; 9) $ and $ {a}_{3}=3 $; in Fig. 8(b), $ {x}''=8 $ and $ {a}_{3}\in (4.15{\mathrm{,}}\; 4.16) $. We can see that as the layer’s medium loss or radius increases, the conventional method breaks down at around $ {x}''=8.3 $ or $ {a}_{3}=4.154 $ m, while the proposed method still works stably.

We have also compared the computational efficiencies exhibited by the proposed formulation and the traditional formulation in the numerical tests above. As seen in Table 4, the computational efficiencies of the proposed method and the traditional method are on the same order. It is because our code based on the traditional formulation has been well optimized that it is a little bit more time efficient. In principle, these two formulations should have the same computational efficiency since they have the same complexity. The absolute computational time depends on the stopping precision set in the code.

Moreover, the accuracy of the proposed method is the same as that of the traditional method since no approximation has been made in the proposed stable formulations employing scaled Bessel functions.

In the computation of SLM theory, when the media’s loss or radius is large, the involved Bessel function’s argument will have a large imaginary part, which could lead to numerical overflow or underflow and computation failure. However, scaled Bessel functions do not diverge in these cases. In this paper, scaled Bessel functions are introduced to the RSLM theory, and stable formulations of the theory are derived. Numerical tests have shown that the proposed approach works stably with high media loss or large layer radii, whereas the conventional formulation breaks down. The proposed approach can be combined with those stabilization methods for small argument cases so as to achieve fully stable computation of SLM theory under various situations. In addition, the proposed approach can be extended to cases of eccentric layered-sphere or layered-cylinder. The stable and analytical formulas derived are also beneficial for training artificial intelligence techniques in computational electromagnetics [18,19].

The authors declare no conflicts of interest.