How to derive the Debye series expansion of the Mie scattered field? The usual Mie scattering theory is just a series expansion of the scattered field and the field inside the scattering sphere, done using a cleverly chosen basis. This lets us calculate the full solution, but provides little insight into the physics of scattering.
OTOH, the Debye series is a decomposition of each term of the Mie series into (an infinite number of) physically-meaningful parts, each corresponding to number of internal reflections of the wave, which lets us examine the scattering process in detail.
The derivation of the Mie solution is relatively easy to find, e.g. in [1]. But where I could find discussion of Debye expansion, I only got either qualitative description, or, like in [2], the algorithm without explanation of how it was derived.
So, how is the Debye series derived? Is there a general procedure that one could apply to any solution of a scattering problem, in order to obtain a similar decomposition?

References

*

*C.F.Bohren, D.R. Huffman, "Absorption and Scattering of Light by Small Particles", Wiley 1983.

*Jianqi Shen and Huarui Wang, "Calculation of Debye series expansion of light scattering", Appl. Opt. 49, 2422-2428 (2010).

 A: $\newcommand{\h}[1]{{h^{(#1)}_\ell}}$
$\newcommand{\j}{j_\ell}$
For simplicity, in this answer I'll consider scattering of a scalar wave (e.g. acoustic), rather than vector-valued EM wave. The case of EM is a bit more complicated because of the two polarizations (TM and TE), but the logic is mostly the same.
Mie solution
This solution can be derived in the same way as described in the ref. 1 in the OP. Below I'll summarize the result, to use it as a reference in the next section detailing the derivation of the Debye expansion.
An incident scalar plane wave can be expanded into a series of standing spherical waves as
\begin{equation}
\operatorname{inc}(r,\theta)=\sum_{\ell=0}^\infty s_\ell \j(kr) P_\ell(\cos\theta),
\tag1
\label{inc}
\end{equation}
where $s_\ell=i^\ell (2\ell+1),$ $\j$ are the spherical Bessel functions of the first kind, and $P_\ell$ are the Legendre polynomials.
Assuming interface conditions of continuity of the wave function and its derivative, we can calculate, separately for each $\ell$, the amplitudes $u_\ell$ of the wave inside the sphere and the $v_\ell$ for the wave outside. For a sphere of radius $a$ and refractive index $n$ the Mie solution for the scattered field (outside the sphere; doesn't include incident wave) is then
\begin{equation}
\operatorname{scat}(r,\theta)=\sum_{\ell=0}^\infty v_\ell s_\ell \h1(kr)P_\ell(\cos\theta),
\tag2
\label{scat}
\end{equation}
where $\h1$ are the spherical Hankel functions of the first kind, and
\begin{equation}
v_\ell=\frac{n\j(ka)\j'(nka)-\j(nka)\j'(ka)}{M_\ell}
\tag3
\label{vl}
\end{equation}
with
\begin{equation}
M_\ell=\j(nka)\h1'(ka)-n\h1(ka)\j'(nka).
\tag4
\end{equation}
The Mie solution for the field inside the sphere is
\begin{equation}
\operatorname{int}(r,\theta)=\sum_{\ell=0}^\infty u_\ell s_\ell \j(nkr)P_\ell(\cos\theta),
\tag5
\label{int}
\end{equation}
with
\begin{equation}
u_\ell=\frac{i}{k^2a^2M_\ell}.
\tag6
\label{ul}
\end{equation}
Debye expansion components
I used Appendix A from ref. I as a guide in the following derivation.
A standing spherical wave (of which $\operatorname{inc}$ in \eqref{inc} consists) is a superposition of a converging and an expanding spherical waves:
\begin{equation}
\j(r)=\frac12\left(\h1(kr)+\h2(kr)\right).
\tag7
\label{jAsH1H2}
\end{equation}
Each of these two running spherical waves are the solutions of wave equations with respectively a sink and a source at the center. The superposition of them can be viewed as a converging spherical wave that crosses the center $r=0$, and then continues propagation, now expanding away from the center.
When a spherical object is present around $r=0$, the converging wave is partially reflected at the surface of the sphere with reflection coefficient $R^{22}_\ell$ determined by the interface conditions. Solving the resulting equations yields
\begin{equation}
R^{22}_\ell=\frac{n\h2(ka)\h2'(nka) - \h2(nka)\h2'(ka)}{D_\ell}
\tag8
\label{R22}
\end{equation}
with the denominator common to all such coefficients for given $\ell$ being
\begin{equation}
D_\ell=\h2(nka)\h1'(ka) - n\h1(ka)\h2'(nka).
\tag9
\end{equation}
Thus, instead of unit amplitude of the outgoing wave that is present in the original incident wave, the reflected term should have amplitude $R^{22}_\ell$, while the original reflected term should be canceled, thus we have $R^{22}_\ell-1$. Inserting the coefficient $\frac12$ from \eqref{jAsH1H2}, we get the $p=0$ terms of the Debye series—the ones representing externally-reflected wave and the shadow of the sphere:
\begin{equation}
\operatorname{debye}_0(r,\theta)=\sum_{\ell=0}^\infty \frac12(R^{22}_\ell-1) s_\ell \h1(kr) P_\ell(\cos\theta).
\tag{10}
\end{equation}
The remaining part of the incident wave is transmitted into the sphere with transmission coefficient $T^{21}_\ell$. (The superscript "21" here designates transition from region 2 (medium) into region 1 (sphere).) The transmission coefficient $T^{21}_\ell$ is determined from the same equations as $R^{22}_\ell$ in \eqref{R22}, with the result being
\begin{equation}
T^{21}_\ell=\frac{2i}{k^2a^2D_\ell}.
\tag{11}
\end{equation}
This transmitted wave crosses the center of the sphere and expands out until it hits the surface of the sphere from inside, at which point it once again partially reflects with reflection coefficient $R^{11}$ and partially transmits out with transmission coefficient $T^{12}$. These coefficients are given by
\begin{equation}
R^{11}_\ell=\frac{n\h1(ka)\h1'(nka)-\h1(nka)\h1'(ka)}{D_\ell},
\tag{12}
\end{equation}
\begin{equation}
T^{12}_\ell=\frac{T^{21}}n.
\tag{13}
\end{equation}
Now every outgoing wave of higher order in the Debye series corresponds to a single transmission from the medium into the sphere, $(p-1)$ internal reflections and a single transmission outwards. Thus corresponding terms in the Debye expansion have the amplitude of
\begin{equation}
T^{21}_\ell(R^{11}_\ell)^{p-1}T^{12}_\ell,
\tag{14}
\end{equation}
which yields corresponding Debye expansion terms
\begin{equation}
\operatorname{debye}_p(r,\theta)=\sum_{\ell=0}^\infty \frac12\left(T^{21}_\ell(R^{11}_\ell)^{p-1}T^{12}_\ell\right) s_\ell \h1(kr) P_\ell(\cos\theta).
\tag{15}
\end{equation}
Connection of Debye series with Mie series
Summing all the Debye series as
\begin{equation}
\sum_{p=0}^\infty \operatorname{debye}_p(r,\theta)
\tag{16}
\end{equation}
by exchanging the order of summation over $\ell$ and over $p$, we'll find that the $\ell$-dependent expansion coefficients
\begin{equation}
\frac12\left(R^{22}_\ell-1+\sum_{p=1}^\infty T^{21}_\ell(R^{11}_\ell)^{p-1}T^{12}_\ell\right)=
\frac12\left(R^{22}_\ell-1+\frac{T^{21}_\ell T^{12}_\ell}{1-R^{11}_\ell}\right)=
v_\ell,
\end{equation}
which means that the Debye series is indeed an exact representation of the Mie solution \eqref{scat}–\eqref{vl}.
Similarly, simplification of the amplitudes of internal waves yields
\begin{equation}
\sum_{p=0}^\infty T^{21}_\ell \left(R^{11}_\ell\right)^p=\frac{T^{21}_\ell}{1-R^{11}}=u_\ell,
\end{equation}
which means that Debye expansion also works for internal field as given by the Mie solution \eqref{int}–\eqref{ul}.
References
I. Edward A. Hovenac and James A. Lock, "Assessing the contributions of surface waves and complex rays to far-field Mie scattering by use of the Debye series", J. Opt. Soc. Am. A 9, 781-795 (1992)
