# Simulation of scattering from layered spheres with known surface electric
field distributions using Mie theory and modified angular spectrum method:
Applications to corneal sensing

Joel Lamberg1  Faezeh Zarrinkhat1,2  Aleksi Tamminen1  Mariangela Baggio1
Juha Ala-Laurinaho1  Juan Rius2 Jordi Romeu2  Elsayed E. M. Khaled3,4 and
Zachary Taylor1  1 Department of Electronics and Nanoengineering, MilliLab,
Aalto University, Espoo, Finland
2 CommSensLab, Technical University of Catalonia/UPC, Barcelona, Spain
3 Department of Electrical Engineering, Assiut University, Assiut, Egypt
4 High Institute of Engineering and Technology, Sohage, Egypt  *
<EMAIL_ADDRESS>
###### Abstract

Mie theory is a powerful method to evaluate the scattered fields from the
multilayered sphere, where the incident field is expanded to the vector
spherical harmonic (VSH) presentation. Then scattered fields are obtained by
the T-matrix method. However, obtaining the VSH coefficients for arbitrarily
shaped incident fields is difficult and time-consuming. This paper proposes a
novel 3D angular spectrum method (3D ASM) for evaluating the VSH coefficients
for the incident field, which is defined from the required electric field
distribution positioned on the spherical surface. This allows the VSH
expansion and evaluation of the scattered fields from a multilayered sphere
illuminated with an arbitrary incident wavefront in the Mie Scattering range.
This has been computationally challenging with previous methods.

First, the advantage of the beam created with the proposed method compared to
the nominal Gaussian beam illumination is addressed with the spherical
bandstop filter simulation. Then the incident field computed by the proposed
method is compared to the physical-optics simulations showing precise
agreement. As an example of the proposed methodology, the cornea is modeled as
a multilayered spherical structure, and the scattered fields are computed from
the cornea illuminated by the incident field with a spherical top-hat and
tapered top-hat wavefronts. Also, the coupling coefficients of the incident
and scattered fields from the cornea model are computed in the 200 - 400 GHz
frequency range. The results are compared with coupling coefficients obtained
with Gaussian beam illumination and referenced to the reflectivity obtained
from plane wave illumination on an analog planar structure. The interaction
between the layered spherical structures and the spherically matched incident
wavefront shows increased agreement with the planar stratified medium theory
compared to the Gaussian beam illumination. The method enables more accurate
estimations of corneal water content.

††journal: oe††articletype: Research Article

##  1 Introduction

Remote sensing of spherical objects with electromagnetic beams has numerous
applications in non-destructive evaluation (NDE) [1], biology [2], and
medicine [3]. One promising medical application is corneal water content
sensing for detection and management of Fuchs’s dystrophy and assessment of
corneal grafts. Fuchs’ dystrophy causes corneal hyperhydration through
degradation of the cornea’s endothelium and its hydration maintenance
functions. Corneal graft rejection is typically preceded by an inflammatory
response characterized by increased corneal tissue water content (CTWC) [4].
In both cases, excess corneal water content is usually detected via visual
inspection of diffuse light scattering arising from swollen collagen fibers.
In both cases, detection by visually apparent corneal edema is too late in the
disease/rejection process to salvage cornea tissue.

Current diagnostic methods restrict in vivo CTWC quantification to pachymetry,
which infers tissue hydration from thickness quantified via ultrasonic or
optical backscatter measurements. However, the mapping from the pachymetry-
measured central corneal thickness (CCT) to corneal tissue hydration does not
account for physiologic variation and returns significant abnormal hydration
for corneal thicknesses off the nominal 580 $\mu$m [5].

Recent research suggests that terahertz (THz) systems are suitable for mapping
human cornea water content by coupling to the cornea’s lossy longitudinal
modes [6, 7]. When referenced to THz band wavelengths, the cornea presents as
an optically smooth, spherical shell structure with radially varying
dielectric properties arising from the water content gradient in the axial
dimension [4]. Furthermore, natural physiologic variation in radius of
curvature, when referenced to THz band wavelengths, is negligible thus
curvature can be assumed and fixed when designing optical systems.

In many reports describing illumination of the cornea, the incident field is
modeled as a Gaussian beam in which the radius of curvature (RoC) of the
wavefront matches the spherical surface at the optical axis [8, 9]. As the
focused spot size of the incident field is less than the radius of the cornea,
the electromagnetic properties of the cornea can be modeled with stratified
media theory (SMT) [10] to capture the dielectric gradient, and effective
media theories (EMT) [11] to compute the effective permittivity at a given
depth. Application of SMT and EMT implies that a plane wave on a planar film
atop an aqueous half-space is a sufficient analogue for spherical cornea under
the right illumination. This matching criterion is called the plane-wave
condition. Optimization routines can then be applied to estimate corneal water
content and thickness from backscattered THz illumination.

However, as the Gaussian beam is observed further from the optical axis, a
mismatch between the wavefront and spherical surface increases between the
field wavefront and the spherical surface [8], see Figure 1. The phase
deviation between the incident wavefront and the desired wavefront matching
with the surface of the sphere reduces the theoretical CTWC measurement
accuracy when the plane-wave condition and corresponding EMT and SMT are
applied. For improved agreement between THz frequency corneal backscatter and
that predicted by plane wave based analysis, there is a need to synthesize and
simulate a field whose wavefront matches the spherical corneal shape at any
point on the surface. In other words, a method that can create the required
incident field from an arbitrary, complex electric field distribution
initiated at the spherical surface. Furthermore, the method should also be
able to model electromagnetic wave interactions between the incident beam and
the cornea interior modeled as a multilayered spherical structure that
generates complex internal and scattered electric fields from the initial
surface distribution. This leads to an inverse imaging problem; first, the
optimal wavefront distribution should be created at the cornea’s surface, and
then the incident beam should be synthesized with the ability to satisfy this
surface criterion.

Internal and scattered fields arising from a homogenous sphere with arbitrary,
given electric field distributions on the spherical surface can be computed
with conventional methods including full-wave simulations, geometrical optics
(GO) [12], or physical optics (PO) [13, 14]. Conventional methods are also
sufficient for calculating propagating fields from source distributions
defined on a sphere. These techniques are well studied and accurate given the
model fidelity and a suitable wavelength range. However, they cannot assess
the internal and scattered electric fields from multilayered spherical objects
without considerable computational effort. This problem occurs in full-wave
simulation as the volume of the multilayered sphere may require discretization
on the order of $10^{9}$ cells, which would exceed available memory on most
modern computers. The same discretization and computer memory constraint
occurs with the boundary element method, as each layer boundary should be
discretized, to, again, on the order of $10^{9}$ boundary elements. Likewise,
geometrical optics is limited to spheres much larger than the illumination
wavelength, and ray splitting to accommodate a sufficient number of layers may
also exceed available memory. Furthermore, obtaining the steady-state results
with physical optics would need tracking of local reflected and transmitted
fields from each boundary leading to an immense number of data points and
equally infeasible for same computer memory limitations.

The ratio of the corneal radius of curvature ($\sim$ 7.8 mm) and THz band
wavelengths ($\sim$ 3 mm - 0.3 mm) places the problem in the Mie scattering
region thus supporting analysis as a pure scattering problem. Mie theory and
the generalized Lorentz-Mie theory are powerful methods to evaluate the
internal and scattered fields from the multilayered spheres [15]. First, the
known incident field is presented in vector spherical harmonics (VSH)
expansion which is defined by VSH coefficients [16]. However, obtaining the
coefficients for arbitrarily shaped incident fields is difficult [15]. These
VSH coefficients for the incident field can be obtained by the angular
spectrum (AS) method from the planar surface distributions [17]. The AS method
presents the incident field in the angular spectrum domain as the sum of
differently oriented plane waves. When this angular summation of the plane
waves is presented in a direction cosine coordinate system, it can be used
directly to compute the incident VSH coefficients. Then the incident field is
expanded to the VSH presentation. Additionally, the internal and scattered
fields from a multilayered dielectric sphere can be mapped with the T-matrix
method in the Mie scattering region [18]. The limitation of this AS method is
that it is nominally compatible with planar interfaces. In [19], the planar AS
method was expanded to approximate diffractive fields from 2D curved surfaces.
This was obtained by dividing the curved 2D surface into the step-wise
subregions windowed by the Gaussian function. This approach was later expanded
to the 3D surfaces in [20], where the authors used the same planar step-wise
subsections with Gaussian distribution. Both methods approximate well the
diffraction, reflections, and transmission from the curved single-layer
boundaries, when the RoC is much larger than the wavelength [21]. Thus these
methods are not suitable for cornea analysis in the THz region, where the RoC
is comparable to the wavelength.

Figure 1: The incident beam propagating towards the sphere, and the mismatch
between the beam wavefront and the spherical surface is proportional to the
distance from the optical axis.

In this study, we start with the same approach of dividing the spherical
surface into planar subregions. These subregions are presented with equal
amplitude, without Gaussian distribution or windowing. Then these subregions
are shrunk into infinitesimally small areas, which approach the Dirac’s delta
function at the limit. This approach gives exact results of the diffractive
field from the spherical surface on the limit of infinitesimally small
partition. Thus, this study proposes an accurate 3D AS method for the beam
synthesis from the spherical surfaces, which can be expanded to arbitrary
surfaces taking into account surface structures comparable to the wavelength.
Also, with the proposed method the radius of curvature of the surface can be
as small as couple of wavelengths. Moreover, this result is equivalent to the
angular spectrum presentation of the Stratton-Chu method with magnetic dipole
sources [22]. The obtained theory is then utilized to create an incident
electric field with the ideal spherical wavefront [23], and the scattered
field is calculated from the cornea model with the layered spherical structure
of 50 layers.

##  2 Theory

This work presents the synthesis of the incident electromagnetic field from
the arbitrary electric field distribution positioned on the spherical surface
with the modified AS method. The method’s main idea is that Riemann’s surface
integral combines differential surface elements, which can be approximated as
locally planar elements. Then the radiated field from each element is created
with the planar AS method. The total electromagnetic field is the sum of the
fields created from the differential sources by the superposition principle
[17]. In this approach, we make a physical assumption that when the areas of
the elements are small enough, the created field from the differential
elements approaches the field from the original surface. This approximation
works well for a smooth surface with radius of curvature that is larger than
the wavelength.

The derivation starts on an arbitrary point on the unit-sphere, where the
spherical components $r,\theta,\phi$ define the surface. These orthogonal
spherical components
$(\mathbf{e}_{n},\,\mathbf{e}_{\theta},\mathbf{e}_{\phi})$ can be used as base
vectors for the local Cartesian coordinate system, whose origin is positioned
at the same observation point on the spherical surface. Figure 2 presents the
situation, where the spherical components on the surface element are used as
base vectors for the local coordinate system. The main steps of the upcoming
derivation follow the following structure;

  * •

The relation between the local Cartesian coordinates of each surface element
and the global spherical coordinates is derived from the unit-sphere base
vectors $(\mathbf{e}_{n},\,\mathbf{e}_{\theta},\mathbf{e}_{\phi})$.

  * •

Each differential element is considered as a field source and the radiated
field is presented in the angular domain by the Fourier transform of the
electric field distribution.

  * •

The total field is obtained by summing together the fields from the inverse
Fourier transform of angular domain presentation.

  * •

Then the total field with polarization vector parallel to the
$\mathbf{e}_{\theta}$ is expanded into the VSH presentation.

  * •

Each local coordinate system is presented in the global coordinate system with
the origin shift for satisfying the VSH symmetry.

  * •

Finally, the equations of the incident field synthesized from the spherical
surface are presented.

After the derivation, the obtained result is proved to be mathematically exact
in Appendix A and the result can be generalized to the spheres with arbitrary
radii by taking account of the wavelength. Furthermore, the scattered field
equations are presented in Appendix B.

###  2.1 Incident field synthesis

Let us start by defining the base vectors and the relation between coordinate
systems. Position vector in standard form is

$\begin{array}[]{l}\mathbf{r}=x\mathbf{e}_{x}+y\mathbf{e}_{y}+z\mathbf{e}_{z},\end{array}$
(1)

where $\mathbf{e}_{x},\mathbf{e}_{y}$ and $\mathbf{e}_{z}$ are orthonormal
Cartesian base-vectors. Position at unit-sphere is

$\begin{array}[]{l}\mathbf{n}=\mathbf{n}(\theta,\phi)=n_{x}\mathbf{e}_{x}+n_{y}\mathbf{e}_{y}+n_{z}\mathbf{e}_{z},\end{array}$
(2)

where $|\mathbf{n}|=1$, and in spherical coordinates
$n_{x}=\sin{\theta}\cos{\phi},n_{y}=\sin{\theta}\sin{\phi}$ and
$n_{z}=\cos{\theta}$.

Figure 2: Global spherical unit-vectors as local Cartesian base-vectors.

The differential surface element on the unit-sphere in spherical coordinates
is $d\Omega=\sin{\theta}d\theta d\phi$. Let the electric field on the segment
of the sphere $\Omega$ be a continuous function
$E_{0}=E_{0}(\theta,\phi),\\{0\leq\theta\leq\theta_{0},0\leq\phi\leq 2\pi$}.
Let the position $\mathbf{r}$ be fixed and let’s observe the electric field at
$\mathbf{n}(\theta,\phi)$. The differential field element is considered to be
on the plane defined by the unit vectors $\mathbf{e}_{\theta}$ and
$-\mathbf{e}_{\phi}$. Together with their normal $\mathbf{e}_{n}=-\mathbf{n}$,
these determine the right handed orthonormal base
$(\mathbf{e}_{\theta},\mathbf{-e}_{\phi},\mathbf{e}_{n})$. Fixing origin at
$\mathbf{n}$ this base defines a local coordinate system, where coordinates
are denoted by upper-bar like $(\bar{x},\bar{y},\bar{z})$. A well known
relation between spherical and Cartesian coordinates gives

$\begin{array}[]{l}\begin{aligned}
\mathbf{e}_{\theta}&=\frac{d\mathbf{n}}{d\theta}=\cos{\theta}\cos{\phi}\mathbf{e}_{x}+\cos{\theta}\sin{\phi}\mathbf{e}_{y}-\sin{\theta}\mathbf{e}_{z}\\\
-\mathbf{e}_{\phi}&=\frac{-1}{\sin{\theta}}\frac{d\mathbf{n}}{d\phi}=\sin{\phi}\mathbf{e}_{x}-\cos{\phi}\mathbf{e}_{y}.\end{aligned}\end{array}$
(3)

The transformation matrix between the coordinates at the base
$(\mathbf{e}_{\theta},\mathbf{-e}_{\phi,}\mathbf{e}_{n})$ and the coordinates
at the standard base $(\mathbf{e}_{x},\mathbf{e}_{y},\mathbf{e}_{z})$, is
defined as

$\begin{array}[]{l}\Theta(\theta,\phi)=\left[\mathbf{e}_{\theta}\,\,\mathbf{-e}_{\phi}\,\,\mathbf{e}_{n}\right],\end{array}$
(4)

where the columns are the vectors $\mathbf{e}_{\theta},\mathbf{-e}_{\phi}$ and
$\mathbf{e}_{n}$. It is an orthonormal matrix, thus $\Theta^{T}\Theta=I$. The
position $\mathbf{\bar{r}}$ in global coordinate system is

$\begin{array}[]{l}\mathbf{\bar{r}}=\mathbf{r}-\mathbf{n}=(x-\sin{\theta}\cos{\phi})\mathbf{e}_{x}+(y-\sin{\theta}\sin{\phi})\mathbf{e}_{y}+(z-\cos{\theta})\mathbf{e}_{z}.\end{array}$
(5)

Let’s change (5) expressed with respect to the base
$(\mathbf{e}_{\theta},\mathbf{-e}_{\phi},\mathbf{e}_{n}):\,\mathbf{r}-\mathbf{n}=\bar{x}\mathbf{e}_{\theta}+\bar{y}(\mathbf{-e}_{\phi})+\bar{z}\mathbf{e}_{n}$.
By (4) we have position $\mathbf{\bar{r}}$ in local coordinate system as

$\begin{array}[]{l}\left[\bar{x}\,\,\bar{y}\,\,\bar{z}\right]^{T}=\Theta^{T}(\mathbf{r}-\mathbf{n})=\Theta^{T}[x-\sin{\theta}\cos{\phi}\,\,\,\,\,\,y-\sin{\theta}\sin{\phi}\,\,\,\,\,\,z-\cos{\theta}]^{T}.\end{array}$
(6)

Let us divide the segment of the sphere $\Omega$ containing the electric field
distribution to finite partition of small separate sets
$\cup_{c}\Omega_{c}=\Omega$. Select the point
$\mathbf{\bar{n}}=\mathbf{\bar{n}}(\bar{\theta},\bar{\phi})$ on differential
element $\Omega_{c}$, the vectors associated with it are $\mathbf{e}_{\theta}$
and $\mathbf{e}_{\phi}$. Area $\Omega_{c}$ is projected radially to the plane
spanned by the vectors $\mathbf{e}_{\theta}$ and $\mathbf{e}_{\phi}$. The
resulting area at the plane is marked as $\Omega_{t}$. Likewise, function
$E_{0}$ is projected radially into that plane, limited to set $\Omega_{t}$ and
zero elsewhere. This geometry is presented in Figure 3. The resulting
restriction is marked as $E_{0}^{t}$. In the local coordinate system it holds
$\bar{z}=0$ on the piece $\Omega_{t}$, because this is located on the
$\bar{x}\bar{y}-$plane. $E_{0}^{t}$ is valid when $\bar{x}\approx 0$ and
$\bar{y}\approx 0$. Let the Fourier transform of the function $E_{0}^{t}$ on
the $\bar{x}\bar{y}-$plane be

$\begin{array}[]{l}A(k_{\bar{x}},k_{\bar{y}};\bar{0})=\mathcal{F}\left\\{E_{0}^{t}\right\\}(k_{\bar{x}},k_{\bar{y}})=\iint_{\mathbb{R}^{2}}E_{0}^{t}e^{-i(k_{\bar{x}}\bar{x}+k_{\bar{y}}\bar{y})}d\bar{x}d\bar{y},\end{array}$
(7)

where $k=|\mathbf{k}|=\sqrt{k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}=2\pi/\lambda$. The
inverse Fourier transform
$\mathcal{F}\left\\{A(k_{\bar{x}},k_{\bar{y}};\bar{0})\right\\}$ gives back
$E^{t}_{0}(\bar{x},\bar{y})=E(\bar{x},\bar{y},\bar{0})$. Analogously the
angular spectrum $A(k_{\bar{x}},k_{\bar{y}};\bar{z})$ gives back the desired
field $E(\bar{x},\bar{y},\bar{z})$ [18], thus

$\begin{array}[]{l}E(\bar{x},\bar{y},\bar{z})=\frac{1}{4\pi^{2}}\iint_{\mathbb{R}^{2}}A(k_{\bar{x}},k_{\bar{y}},\bar{z})e^{i(k_{\bar{x}}\bar{x}+k_{\bar{y}}\bar{y})}dk_{\bar{x}}dk_{\bar{y}}.\end{array}$
(8)

The relation between $A(k_{\bar{x}},k_{\bar{y}};\bar{0})$ and
$A(k_{\bar{x}},k_{\bar{y}};\bar{z})$ can be obtained by operating the
Helmholtz equation $\nabla^{2}E+k^{2}E=0$ to the function (8) [17]. The
obtained relation is

$\begin{array}[]{l}A(k_{\bar{x}},k_{\bar{y}},\bar{z})=A(k_{\bar{x}},k_{\bar{y}};\bar{0})e^{ik_{\bar{z}}|\bar{z}|}=\mathcal{F}\big{\\{}E_{0}^{t}\big{\\}}(k_{\bar{x}},k_{\bar{y}})e^{i|\bar{z}|\sqrt{k^{2}-k_{\bar{x}}^{2}-k_{\bar{y}}^{2}}},\end{array}$
(9)

where the solution $A(k_{\bar{x}},k_{\bar{y}};\bar{z})=C_{\pm}e^{\pm\
ik_{\bar{z}}\bar{z}}$ allows both propagation directions $\pm z$ as
$e^{ik_{\bar{z}}|\bar{z}|}$. Above the angular spectrum theory is locally
applied on the local coordinate system to function $E_{0}^{t}$ valid on piece
$\Omega_{t}$. Let’s observe the field it creates at point
$\mathbf{\bar{r}}=\bar{x}\mathbf{e}_{\theta}+\bar{y}(\mathbf{-e}_{\phi})+\bar{z}\mathbf{e}_{n}$,
where the coordinates are as in (6). As shown above, it holds at a fixed point
$\mathbf{r}$

$\begin{array}[]{l}E_{t}(\mathbf{r})=\frac{1}{4\pi^{2}}\iint_{\mathbb{R}^{2}}\mathcal{F}\big{\\{}E_{0}^{t}\big{\\}}(k_{\bar{x}},k_{\bar{y}})e^{i(k_{\bar{x}}\bar{x}+k_{\bar{y}}\bar{y})}e^{i|\bar{z}|\sqrt{k^{2}-k_{\bar{x}}^{2}-k_{\bar{y}}^{2}}}dk_{\bar{x}}dk_{\bar{y}}.\end{array}$
(10)

The total field from the spherical surface is obtained by summing the fields
from the differential sources by superposition principle as

$\begin{array}[]{l}\begin{aligned}
\mathbf{E_{\theta}}(\mathbf{r})&\approx\sum_{t}E_{t}(\mathbf{r})\mathbf{e}_{\theta}=\sum_{t}\mathbf{e}_{\theta}\frac{1}{4\pi^{2}}\iint_{\mathbb{R}^{2}}\mathcal{F}\big{\\{}E_{0}^{t}\big{\\}}(k_{\bar{x}},k_{\bar{y}})e^{i(k_{\bar{x}}\bar{x}+k_{\bar{y}}\bar{y})}\\\
&\times
e^{i|\bar{z}|\sqrt{k^{2}-k_{\bar{x}}^{2}-k_{\bar{y}}^{2}}}dk_{\bar{x}}dk_{\bar{y}},\end{aligned}\end{array}$
(11)

where local coordinates $(\bar{x},\bar{y},\bar{z})$ as in (6) and the total
field is polarized along the global $\mathbf{e}_{\theta}$ unit-vector. Because
the operator $\Delta=\nabla^{2}$ is linear, the field (11) also satisfies the
Helmholtz equation as unit vector $\mathbf{e}_{\theta}$ is locally constant.
In order to be precise, if the limit with respect to the areas $|\Omega_{t}|$
of the sets $\Omega_{t}$ exists, we define mathematically precisely

$\begin{array}[]{l}\mathbf{E}_{\theta}(\mathbf{r})=\mathbf{E}_{\theta}(x,y,z)=lim_{|\Omega_{t}|\to
0}\sum_{t}E_{t}(\mathbf{r})\,\mathbf{e}_{\theta}.\end{array}$ (12)

As pointed out at the beginning of the Theory section, this limit also works
physically when considering smooth surfaces, and is even more precise when the
wavelength decreases compared to the radius of the sphere. Finally, a further
examination is made to obtain a closed-form for the clause (12); we postpone
an exact proof at the Appendix A.

Heuristically, when $|\Omega_{t}|$ is getting smaller, the function
$E_{0}^{t}/E_{0}^{t}(\bar{0},\bar{0})|\Omega_{t}|$, normed in volume, is
approaching Dirac delta $\delta_{(0,0)}$ in the sense of the distribution
theory;
$E_{0}^{t}/E_{0}^{t}(\bar{0},\bar{0})|\Omega_{t}|\rightarrow\delta_{(0,0)}$
thus $E_{0}^{t}\rightarrow
E_{0}^{t}(\bar{0},\bar{0})|\Omega_{t}|\delta_{(0,0)}$. As known, the Fourier
transform of Dirac delta is a constant and equals to 1. Thus, by small
$|\Omega_{t}|$

$\begin{array}[]{l}\begin{aligned}
\mathcal{F}\big{\\{}E_{0}^{t}\big{\\}}(k_{\bar{x}},k_{\bar{y}})&=\iint_{\mathbb{R}^{2}}E_{0}^{t}e^{-i(k_{\bar{x}}\bar{x}+k_{\bar{y}}\bar{y})}d\bar{x}d\bar{y}\\\
&\approx\iint_{\mathbb{R}^{2}}E_{0}^{t}(\bar{0},\bar{0})|\Omega_{t}|\delta_{(0,0)}e^{-i(k_{\bar{x}}\bar{x}+k_{\bar{y}}\bar{y})}d\bar{x}d\bar{y}\\\
&=E_{0}(\theta,\phi)|\Omega|\mathcal{F}\big{\\{}\delta_{(0,0)}\big{\\}}=E_{0}(\theta,\phi)|\Omega_{t}|,\end{aligned}\end{array}$
(13)

since $E_{0}(\theta,\phi)=E_{0}^{t}(\bar{0},\bar{0})$ when fixed $\mathbf{n}$.
Then (11) and (13) yield first the form of a Riemann sum, and then an integral
as a limit

$\begin{array}[]{l}\begin{aligned}
\mathbf{E}_{\theta}(\mathbf{r})&=\mathbf{E}_{\theta}(x,y,z)=lim_{|\Omega_{t}|\to
0}\sum_{t}\mathbf{e}_{\theta}\frac{1}{4\pi^{2}}E_{0}(\theta,\phi)|\Omega_{t}|\\\
&\times\iint_{\mathbb{R}^{2}}e^{i(k_{\bar{x}}\bar{x}+k_{\bar{y}}\bar{y})}e^{i|\bar{z}|\sqrt{k^{2}-k_{\bar{x}}^{2}-k_{\bar{y}}^{2}}}dk_{\bar{x}}dk_{\bar{y}}\\\
&=\frac{1}{4\pi^{2}}\iint_{S}E_{0}(\theta,\phi)\mathbf{e}_{\theta}(\theta,\phi)\iint_{\mathbb{R}^{2}}e^{i(k_{\bar{x}}\bar{x}+k_{\bar{y}}\bar{y})}e^{i|\bar{z}|\sqrt{k^{2}-k_{\bar{x}}^{2}-k_{\bar{y}}^{2}}}dk_{\bar{x}}dk_{\bar{y}}d\Omega.\end{aligned}\end{array}$
(14)

The final form can be written as

$\begin{array}[]{l}\mathbf{E}_{\theta}(\mathbf{r})=\frac{1}{4\pi^{2}}\iint_{\Omega}E_{0}(\theta,\phi)\mathbf{E}_{t}(\mathbf{r};\theta,\phi)\sin{\theta}d\theta
d\phi,\end{array}$ (15)

where

$\begin{array}[]{l}\mathbf{E}_{t}(\mathbf{r};\theta,\phi)=\mathbf{e}_{\theta}(\theta,\phi)\iint_{\mathbb{R}^{2}}e^{i(k_{\bar{x}}\bar{x}+k_{\bar{y}}\bar{y})}e^{i|\bar{z}|\sqrt{k^{2}-k_{\bar{x}}^{2}-k_{\bar{y}}^{2}}}dk_{\bar{x}}dk_{\bar{y}}\\\
\end{array}$ (16)

and $\bar{x},\bar{y}$ and $\bar{z}$ are as in (6). It is worthwhile to mention
that equation (16) is similar to the
$\mathbf{E}_{t}=-\nabla\times(\mathbf{M}G)$ after replacing the scalar Green’s
function $G$ by its angular spectrum given by Weyl identity, where
$\mathbf{M}=-2\mathbf{n}\times\mathbf{e}_{\theta}$ is a small magnetic dipole
[24]. This leads equation (15) to be a particular case of the Stratton-Chu
equation for electric field in which only the magnetic sources
$\mathbf{M}=-\mathbf{n}\times\mathbf{E}_{surf}$ are presented and the electric
ones $\mathbf{J}=\mathbf{n}\times\mathbf{H}_{surf}$ have been removed [25].

Figure 3: Spherical segment $\Omega_{c}$ projected into a plane $\Omega_{t}$.

###  2.2 Polarization

In the derived equation (15), the electric field component is along
$\mathbf{e}_{\theta}$ at the spherical surface, creating a globally continuous
polarization. The polarization of each differential source is selected on
$\bar{x}\bar{z}-$plane at source’s local coordinate system. Let us consider
function (16), with the electric field along $\mathbf{e}_{\theta}$

$\begin{array}[]{l}\mathbf{E}_{t}(\bar{x},\bar{y},\bar{z})=\mathbf{e}_{\theta}\iint_{\mathbb{R}^{2}}e^{i\mathbf{\bar{k}}\cdot\mathbf{\bar{r}}_{||}}dk_{\bar{x}}dk_{\bar{y}},\end{array}$
(17)

where $\mathbf{\bar{r}}_{||}=(\bar{x},\bar{y},|\bar{z}|)$. Equation (17) will
give us the local $\bar{x}$-component of the electric field. The local
$\bar{z}-$component for the polarization is obtained from
$\mathbf{E}_{t}\cdot\mathbf{\bar{k}}=E_{\bar{x}}k_{\bar{x}}+E_{\bar{y}}k_{\bar{y}}+E_{\bar{z}}k_{\bar{z}}=0$
for each plane wave [26]. Because of polarization on the
$\bar{x}\bar{z}-$plane as mentioned above, we have $E_{y}=0$ and

$\begin{array}[]{l}E_{\bar{z}}=-\Big{(}\frac{k_{\bar{x}}}{k_{\bar{z}}}\Big{)}E_{\bar{x}}.\end{array}$
(18)

The final form of (17) is obtained as

$\begin{array}[]{l}\mathbf{E}_{t}(\bar{x},\bar{y},\bar{z})=\iint_{\mathbb{R}^{2}}e^{i\mathbf{\bar{k}}\cdot\mathbf{\bar{r}}_{||}}\Big{[}\mathbf{e}_{\bar{x}}-\Big{(}\frac{k_{\bar{x}}}{k_{\bar{z}}}\Big{)}\mathbf{e}_{\bar{z}}\Big{]}dk_{\bar{x}}dk_{\bar{y}},\end{array}$
(19)

where $\mathbf{e}_{\bar{x}}=\mathbf{e}_{\theta},\
\mathbf{e}_{\bar{y}}=-\mathbf{e}_{\phi}$ and
$\mathbf{e}_{\bar{z}}=\mathbf{e}_{n}=-\mathbf{n}$ as defined in (2) and (3).
The total polarization at the spherical surface is thus in line with the
orthodrome lines between the poles of the sphere, and creates a desired
polarization when the middle point of the segment $\Omega$ is located at the
$\theta=\pi/2$. The function (19) is integrated over the real disk
$k_{\bar{x}}^{2}+k_{\bar{y}}^{2}\leq k^{2}$ for obtaining only propagating
waves and those waves outside of the disk are evanescent. Using direction
cosine coordinate system, the above equation can be written as [18]

$\begin{array}[]{l}\mathbf{E}_{t}(\bar{x},\bar{y},\bar{z})=k^{2}\int_{0}^{\pi}\sin{\xi}\Big{\\{}\int_{0}^{\pi}\mathbf{e}_{0}(\xi,\zeta)e^{ik(\mathbf{\bar{s}}\cdot\mathbf{\bar{r}}_{||})}\sin{\zeta}d\zeta\Big{\\}}d\xi,\end{array}$
(20)

where $\xi$ and $\zeta$ are direction cosines related to spherical coordinates
as $\xi=\cos^{-1}{(\sin{\theta}\cos{\phi})}$ and
$\zeta=\cos^{-1}{(\sin{\theta}\sin{\phi})}$. The polarization $\mathbf{e}_{0}$
in terms of direction cosines is given as

$\begin{array}[]{l}\mathbf{e}_{0}(\xi,\zeta)=\mathbf{e}_{\bar{x}}-\frac{\cos{\xi}}{(1-s^{2})^{1/2}}\mathbf{e}_{\bar{z}}.\end{array}$
(21)

Additionally $s^{2}=\cos^{2}{\xi}+\cos^{2}{\zeta}$, and the vector
$\mathbf{\bar{s}}$ is

$\begin{array}[]{l}\mathbf{\bar{s}}=\cos{\xi}\mathbf{e}_{\bar{x}}+\cos{\zeta}\mathbf{e}_{\bar{y}}+(1-s^{2})^{1/2}\mathbf{e}_{\bar{z}}.\end{array}$
(22)

###  2.3 VSH expansion and origin shift

The total incident electric field from the spherical surface distribution is
obtained from the formula (15). When (16) is expanded to the vector spherical
harmonics, the origin of the local coordinate system needs to be shifted to
the global origin due to the spherical symmetry of the VSH presentation. This
means that the dielectric spheres in each local coordinate system will be
centered on the local origin, i.e., at the surface. Thus each local spherical
object needs to be shifted to the global origin while the position of the
surface source is kept fixed. This is done by first shifting the local origins
to the global origins and then shifting the positions of the differential
source elements back to the original places. After this procedure, each local
source element is presented in the global coordinate system with the ability
to add the dielectric sphere to the right location. To do so, equation (20) is
evaluated by numerical integration over a disk of radius $k$ using the
trapezoidal rule with uniform-width $l$ [18]

$\begin{array}[]{l}\mathbf{E}_{t_{ij}}\approx
p^{2}\sum_{i}\sin{\xi_{i}}\sum_{j}\sin{\zeta_{j}}\mathbf{e}_{0_{ij}}(\xi_{t_{i}},\zeta_{t_{j}})e^{ik(\mathbf{\bar{s_{ij}}}\cdot\mathbf{\bar{r}}_{||})},\end{array}$
(23)

where the indices $i$ and $j$ refers to plane waves with propagation constants
$k_{x}^{i}$ and $k_{y}^{j}$ respectively to $\bar{x}$ and $\bar{y}$
directions. $p=\xi_{t+1}-\xi_{t}=\zeta_{t+1}-\zeta_{t}$ is the step-size of
the numerical integration. Integration limit of radius $k$ is selected,
because the electric field outside that area is rapidly decreasing. The
trapezoidal rule method is highly accurate on periodic functions [27], such as
equation (20). The transformation matrix from ($\bar{x}\bar{y}\bar{z}$) to
($xyz$) is $\Theta(\theta,\phi)$ as in (4). We want to express the dot product
$\mathbf{\bar{s}}_{ij}\cdot\mathbf{\bar{r}}_{||}$ in the global coordinates.
Let’s write $\mathbf{\bar{r}}=(\bar{x},\bar{y},\bar{z})$. We have a simple
connection $\mathbf{r}-\mathbf{n}=\Theta\,\mathbf{\bar{r}}$ as in (6); that
will be used. So we define $\mathbf{\bar{s}}_{ij}^{||}=\mathbf{\bar{s}}_{ij}$,
when $\bar{z}>0$, but, when $\bar{z}<0$, the sign of $\bar{z}-$coordinates in
$\mathbf{\bar{s}}_{ij}$ is exchanged. Consequently, we can write
$\mathbf{\bar{s}}_{ij}\cdot\mathbf{\bar{r}}_{||}=\mathbf{\bar{s}}_{ij}^{||}\cdot\mathbf{\bar{r}}$.
Let us define

$\begin{array}[]{l}\mathbf{s}_{ij}=\Theta\,\mathbf{\bar{s}}_{ij}^{||}.\end{array}$
(24)

As a result we get

$\begin{array}[]{l}\mathbf{\bar{s}}_{ij}\cdot\mathbf{\bar{r}}_{||}=\mathbf{\bar{s}}_{ij}^{||}\cdot\mathbf{\bar{r}}=\mathbf{\bar{s}}_{ij}^{||}\cdot\Theta^{T}(\mathbf{r}-\mathbf{n})=\Theta\,\mathbf{\bar{s}}_{ij}^{||}\cdot(\mathbf{r}-\mathbf{n})\\\
=\mathbf{s}_{ij}\cdot(\mathbf{r}-\mathbf{n})=\mathbf{s}_{ij}\cdot\mathbf{r}-\mathbf{s}_{ij}\cdot\mathbf{n}.\end{array}$
(25)

Based on (25), the term
$e^{ik(\mathbf{\bar{s_{ij}}}\cdot\mathbf{\bar{r}}_{||})}$, which is in the
local coordinate system, can be presented in the global coordinate system as

$\begin{array}[]{l}e^{ik(\mathbf{\bar{s_{ij}}}\cdot\mathbf{\bar{r}}_{||})}=e^{ik(\mathbf{s_{ij}}\cdot\mathbf{r})}e^{-ik(\mathbf{s_{ij}}\cdot\mathbf{n})},\end{array}$
(26)

where $e^{-ik(\mathbf{s_{ij}}\cdot\mathbf{n})}$ presents the shift of beam
waist of the differential source from global origin to the unit-sphere
surface. Now each plane wave in the angular spectrum created from the
differential source is presented in a global coordinate system. The total
incident electric field is obtained as the sum of the plane waves in the
angular spectrum, and at an arbitrary point $\mathbf{r}$ in space expressed by
VSHs can be written as [18]

$\begin{array}[]{l}\begin{aligned}
\mathbf{E}_{t}(k\mathbf{r})&=p^{2}\sum_{m}\sum_{n}D_{mn}\Big{[}a^{t}_{\text{emn}}\mathbf{M}^{1}_{\text{emn}}(k\mathbf{\bar{r}})+a^{t}_{\text{omn}}\mathbf{M}^{1}_{\text{omn}}(k\mathbf{\bar{r}})\\\
&+b^{t}_{\text{emn}}\mathbf{N}^{1}_{\text{emn}}(k\mathbf{\bar{r}})+b^{t}_{\text{omn}}\mathbf{N}^{1}_{\text{omn}}(k\mathbf{\bar{r}})\Big{]},\end{aligned}\end{array}$
(27)

where $D_{mn}=\frac{\epsilon_{m}(2n+1)(n-m)!}{4n(n+1)(n+m)!}$ is a
normalization factor. $\epsilon_{m}$ is equal to 1 for $m=0$ and is equal to 2
for $m>0$.
$\mathbf{M}^{1}_{\text{emn}},\mathbf{M}^{1}_{\text{omn}},\mathbf{N}^{1}_{\text{emn}}$
and $\mathbf{N}^{1}_{\text{omn}}$ are the VSH of the first kind and
$a^{t}_{\text{emn}},a^{t}_{\text{omn}},b^{t}_{\text{emn}}$ and
$b^{t}_{\text{omn}}$ are incident beam expansion coefficients given as

$\begin{array}[]{l}a^{t}_{\text{emn}}=\sum_{i}\sin{\xi_{i}}\sum_{j}\sin{\zeta_{j}}a^{t}_{{\text{emn}}_{ij}}e^{-ik(\mathbf{s_{ij}}\cdot\mathbf{n})},\\\
a^{t}_{\text{omn}}=\sum_{i}\sin{\xi_{i}}\sum_{j}\sin{\zeta_{j}}a^{t}_{{\text{omn}}_{ij}}e^{-ik(\mathbf{s_{ij}}\cdot\mathbf{n})},\\\
b^{t}_{\text{emn}}=\sum_{i}\sin{\xi_{i}}\sum_{j}\sin{\zeta_{j}}b^{t}_{{\text{emn}}_{ij}}e^{-ik(\mathbf{s_{ij}}\cdot\mathbf{n})},\\\
b^{t}_{\text{omn}}=\sum_{i}\sin{\xi_{i}}\sum_{j}\sin{\zeta_{j}}b^{t}_{{\text{omn}}_{ij}}e^{-ik(\mathbf{s_{ij}}\cdot\mathbf{n})}.\end{array}$
(28)

The VSH coefficients for each $ij-$ plane wave
$a^{t}_{{\text{emn}}_{ij}},a^{t}_{{\text{omn}}_{ij}},b^{t}_{{\text{emn}}_{ij}}$
and $b^{t}_{{\text{omn}}_{ij}}$ are defined in Appendix B. The equation for
the incident field (15), adapted by polarization (19), can also be written in
vector spherical expansion (27) as

$\begin{array}[]{l}\mathbf{E}(\mathbf{r})_{inc}=\frac{1}{4\pi^{2}}\iint_{\Omega}E_{0}(\theta,\phi)\mathbf{E}_{t}(k\mathbf{r};\theta,\phi)\sin{\theta}d\theta
d\phi.\end{array}$ (29)

A similar scattered field presentation mapped with the T-matrix method is
presented as

$\begin{array}[]{l}\mathbf{E}(\mathbf{r})_{sca}=\frac{1}{4\pi^{2}}\iint_{\Omega}E_{0}(\theta,\phi)\mathbf{E}_{scat}(k\mathbf{r};\theta,\phi)\sin{\theta}d\theta
d\phi,\end{array}$ (30)

where $\mathbf{E}_{scat}$ is introduced in Appendix B.

##  3 Results

Now we present three different simulation scenarios, and they are: 1) Incident
field simulations and comparison to similar physical optics simulations, to
validate the accuracy of the field synthesis. 2) The spherical bandstop filter
simulations to verify the benefits of the proposed method. 3) The cornea
simulations to confirm the increased accuracy of cornea´s thickness and water
content analysis.

###  3.1 Incident field simulations

First, the incident electric field was synthesized from the top-hat field
distribution positioned on the spherical surface with Fibonacci sampling using
a derived equation (29). This field has a matching wavefront at the spherical
surface, with the main polarization along the $\mathbf{e}_{\theta}$ direction.
The incident field was evaluated at the plane transverse to the propagation
x-axis at 40 mm from the origin of the sphere. The position of the evaluation
plane along the propagation axis can be changed without restraints. In our
case, the evaluation plane location is selected to match the approximated
location of the lens in a quasioptical system. At this distance, the
converging beam can be focused with reasonably small focusing elements. Also,
this location is identical to the lens location in a future quasioptical
system, designed to create the desired field distribution matching the
spherical surface. All the following simulations were realized with a constant
electric field within subtended angle of 30∘ on the spherical surface and the
transverse evaluation plane at the yz-plane, see Figure 4. The only exception
is shown in Figure 8, where the evaluation plane is along the xy-plane.

Figure 4: The simulation arrangement of beam propagation between the spherical
surface and evaluating plane, where evaluating plane is marked as blue dots.
This plane has 100 $\times$ 100 mm2 area, located at x = 40 mm from the origin
of the sphere. The source grid with 681 Fibonacci distributed differential
sources is marked as red dots at the spherical surface and at the enlargement
of the surface grid. The subtended angle of the spherical cap enclosing the
differential source grid is 30∘.

Then the Cartesian components of the incident field ($E_{3D\,ASM}$) calculated
by the proposed method were compared to the identical physical-optics
simulation ($E_{PO}$) at 175 GHz frequency. The comparizon provides precise
matching down to less than a $-41$ dB average difference in amplitude and less
than 1∘ average difference in phase, see Figure 5. This difference was
calculated for each component as $||E_{3D\,ASM}-E_{PO}||_{ave}$ for the
amplitude and $||\angle E_{3D\,ASM}-\angle E_{PO}||_{ave}$ for the phase. The
magnitude and phase components of the electric field computed by the physical-
optics are visually identical to the $E_{3D\,ASM}$ results, thus they are not
plotted.

Figure 5: Subfigures $\mathbf{a)}$-$\mathbf{c)}$ are the magnitudes of the
electric field components computed by the proposed 3D ASM method, and
$\mathbf{d)}$-$\mathbf{f)}$ are the magnitude differences between the 3D ASM
and PO method at 175 GHz. Subfigures $\mathbf{g)}$-$\mathbf{i)}$ are the
phases of the electric field components computed by the proposed 3D ASM
method, and $\mathbf{j)}$-$\mathbf{l)}$ are the phase differences between the
3D ASM and PO method at 175 GHz.

###  3.2 Banstop filter backscatter simulations

After the simulation verification of the incident field, let’s investigate the
analogy of the spherical layered structure compared to the planar layered
structure, see Figure 6. To compare these structures properly, we must take
into account two different scenarios. When the spherical layered structure has
a lossless core, it creates a cavity with multiple internal reflections inside
the core. Then the analog planar structure is presented in subfigure 6 b). If
the spherical layered structure has a lossy core to prevent the internal
reflection in the core, the analog planar structure is presented in subfigure
6 c).

Figure 6: Subfigure $\mathbf{a)}$ is stratified spherical shell structure,
subfigure $\mathbf{b)}$ presents the analog planar structure when the
spherical structure have a lossless core, and subfigure $\mathbf{c)}$ presents
the analog planar structure when the spherical structure have a lossy core.

Then, the benefits of illuminating spherical layered structures with the beam
with a spherically matching wavefront compared to the Gaussian beam are
presented. This is done by reproducing a planar 27-layer bandstop filter to an
analog spherical structure. The original filter is 3627 nm thick with
alternating layers of $TiO_{2}$ and $SiO_{2}$, with refractive indexes of 2.5
and 1.45 respectively [28]. A spherically analog structure was created with a
4000 nm radius sphere with a lossless and lossy ($\epsilon=1-0.5i$) core to
prevent the internal reflections. Then spherical bandpass filter was
illuminated in the 450 - 700 nm wavelength range by a normal and tapered top-
hat beam created by the proposed 3D ASM method and by Gaussian beam with
different beam size and position scenarios presented in [29, 30]. These
results are compared to the plane wave illumination of the similar planar
filter structures in Figure 7.

Figure 7: Coupling coefficients from the spherical layered filter model, where
$\mathbf{a)}$ presents the lossy core with $\epsilon=1-0.5i$ permittivity and
$\mathbf{b)}$ presents the lossless core scenario.

In this simulation example, the spherical layers extend the entire radius of
the sphere from the surface to the core, allowing the beam to disperse as it
propagates to the deepest layers. This phenomenon produces a small difference
between the planar stratified medium theory model with PW illumination and the
model of the spherical structure with converging beam illumination.
Nevertheless, the simulation shows an increased and relatively well-matched
coupling coefficient ($K$) between the incident and scattered fields with
normal and tapered top-hat beams compared to the superconfocal Gaussian beam
illumination. The formula for computing $K$ is presented in Appendix B.

###  3.3 Cornea backscatter simulations

The cornea is the transparent front part of the eye covering the iris, pupil,
and anterior chamber. The water percentage of the cornea increases gradually
from the anterior surface to the posterior surface. Thus, the permittivity of
the cornea differs radially and can be modeled by the spherical multilayered
structure. The cornea was modeled as a 7.8 mm radius sphere with a 580 $\mu$m
thick multilayered spherical shell structure with a pure-water core. The
permittivity gradient of the cornea was modeled with 50 layers. Each layer’s
permittivity was obtained with the effective medium theory via the Bruggeman
model [31], and the permittivity of the pure-water core was calculated by the
double-Debye model [32]. The corneal shell consisted of 40 % water on the
anterior surface and 70 % water on the posterior surface. The cornea is
modeled as relatively thin layered structure compared to the radius of the
sphere. Due to this, the propagating beam does not disperse as much as in the
bandpass filter and the difference between the scattering from planar and
spherical structure is decreased. First, the total field
$\mathbf{E}_{tot}=\mathbf{E}_{inc}+\mathbf{E}_{scat}$, including the scattered
field from the 50-layer cornea model was simulated in the plane along the
propagation axis at 175 GHz frequency. This is done for illustrating the beam
behavior and the equal phase wavefront on the spherical surface, where the
internal field $\mathbf{E}_{int}$ is not plotted, see Figure 8. The scattered
field was calculated by the formula (49) from Appendix B. It’s worthwhile to
mention that the 3D ASM method has no singularities at the source points due
to its eigenfunction nature, and the incident and scattered fields can be
evaluated at the source points as well. The source points are included in
Figure 8 at the spherical surface.

Figure 8: $\mathbf{a)}$ presents the radially changing permittivity of the
cornea shell, $\mathbf{b)}$ is amplitude (dB) and $\mathbf{c)}$ is phase (rad)
of the total electric field outside the cornea at xy-plane at 175 GHz. The
cornea is represented with a red circle. The top part in $\mathbf{b)}$ shows a
zoom-in plot in the vicinity of the cornea marked with the rectangle in the
lower part of the red circle.

Each Cartesian component of the scattered field from the same cornea model was
computed at the transverse evaluation plane arrangement, introduced in the
Figure 4. The components are normalized by the main polarization
$E_{z}$-component of the incident field from the Figure 5 (c).

Figure 9: The scattering simulation from the multilayered spherical structure
with 3D ASM method at 175 GHz. Subfigures $\mathbf{a)}$-$\mathbf{c)}$ presents
the magnitudes, and $\mathbf{d)}$-$\mathbf{f)}$ present the phases of the
Cartesian components at the transverse plane at 40 mm from the origin of the
cornea. The cornea was modeled as 50-layer structure.

.

The change in amplitude distribution of the scattered field compared to the
incident field is in an excellent match down to less than $-38$ dB difference
to the change calculated from a similar multilayered stratified medium
structure, i.e., 50 layered planar surfaces illuminated with the plane wave.
This comparison was made as
${E_{3D\,ASM}}_{scat.}/{E_{3D\,ASM}}_{inc.}-{E_{PW}}_{scat.}/{E_{PW}}_{inc.}$.
The plane-wave illumination to the planar surface gives the reference
scattering, which is used to evaluate the thickness and CTWC of the cornea.

The coupling coefficient $K$ between the incident and scattered fields plays a
critical role in the cornea’s CTWC and thickness analysis. In this analysis,
the obtained K in the observed frequency range is fitted to match the cornea’s
CTWC- and thickness-dependent reflectivity curve from the planar reference
structure. In this paper, the K from the cornea model was calculated with two
different illumination scenarios in the 200 GHz - 400 GHz band; by the field
synthesized from the spherical top-hat field distribution and by the Gaussian
beam, whose frequency-dependent waist was located at the origin of the cornea
[8]. Obtained values were compared to the reference values, i.e, plane-wave
illumination of the planar layered structure with the same layer thickness and
permittivity values, see Figure 10. The K of the equal phase top-hat
distribution indicates a high match with amplitude less than 0.8% and tapered
top-hat less than 0.7% average difference compared to the reference plane wave
scenario. Typical Gaussian-beam illumination produces a 4-5% amplitude
difference. This improved matching in K cut down the fitting errors in the
cornea’s CTWC and thickness analysis. It is also worthwhile to mention that
the total reflected energy of the beam with the spherical tapered and
untapered top-hat distribution matches with less than 0.07% average difference
compared to the plane wave scenario, thus we can conclude that the small
difference in the K magnitude of the top-hat beam and reference values are due
to the differences in the diffraction patterns between the incident and
scattered fields.

Figure 10: The magnitude of the coupling coefficient from 200 to 400 GHz
frequency range in four scenarios. The blue line presents the reference value
from plane wave illumination into the planar structure, the red line presents
the top-hat illumination, the yellow line presents tapered top-hat
illumination, and the purple line presents the typical Gaussian beam
illumination for spherical structure.

##  4 Conclusions

This paper proposes a novel method for expanding the nominally planar AS
method to synthesize incident electromagnetic beams from electric field
distributions on spherical surfaces. The derivation of the method was
presented, and the resulting equations were proved mathematically in Appendix
A. The proposed method can be expanded to synthesize beams from the arbitrary
smooth surfaces taking the radii of curvatures and frequency range into
account. The minimum radius of curvature of the surface compared to the
wavelength is $RoC>2\lambda$. The benefit of this method is the ability to
present the synthesized fields in the angular spectrum domain, which can be
used to compute the VSH coefficients of the incident fields. Thus the incident
field can be presented with VSH expansion and scattered fields from
multilayered spherical structures, of an arbitrarily large number of layers,
can be obtained utilizing the Mie theory. The computational effort of
presented 50-layer spherical scattering simulations requires much less memory
and fewer arithmetical operations than, for example, similar full-wave
computation.

First, the advantages of the beams with spherically matching wavefronts
compared to the Gaussian beam illumination for spherical objects were
presented. This was done by comparing the scattering behavior of the top-hat,
tapered top-hat and Gaussian beam illumination from the spherical 27-layer
bandpass filter. These results were compared to the reference scattering from
interaction with plane wave and planar 27-layer bandpass filter. The beams
with spherically matching wavefront had significantly increased coupling
coefficient compared to the traditional Gaussian beam illumination. Next, the
Cartesian components of the incident field were evaluated at the transverse
plane at 40 mm from the origin of the cornea. These results were compared to
the physical optics simulations and demonstrated good agreement. Then the
gradient permittivity of the cornea was modeled as a multilayered spherical
structure of 50 layers, and the scattered beam was evaluated from the ideal
incident field with a spherical top-hat and tapered top-hat wavefront by the
T-matrix method. The coupling coefficients between the incident and scattered
fields of the spherical top-hat distributions were evaluated in the 200 GHz to
400 GHz region and compared to Gaussian-beam illumination and the reference
plane-wave illumination on an analog planar surface. The coupling coefficient
amplitude of the spherical tapered top-hat illumination was closest to the
reference amplitude with less than 0.7% average difference. The good agreement
between the coupling coefficient and the reference validates the proposed
method. The results indicate that the incident beam with a matching spherical
wavefront interacts with the multilayered sphere similar to the plane wave
interacting with the planar structures. Thus, the developed method enables the
synthesis of beams whose interaction with the cornea can be more accurately
predicted/analyzed with standard EMT and SMT and improves the accuracy of CTWC
and thickness measurement in THz band.

The presented method will enable the computation of electric field propagation
between different surfaces. The desired incident field can be synthesized from
the cornea’s surface distribution as presented in this paper and computed at
the optical elements. This electric field distribution at the optical element
can then be used as a source distribution and propagated back to the cornea.
Through this methodology, the minimum optical element apertures can be
evaluated by comparing the surface field distribution from the propagated
field to the original surface distribution.

##  5 Appendix A

######  Theorem 5.1.

Let the function $E_{0}$ be continuous in a closed segment of sphere $\Omega$
and let $\mathbf{r}\in\mathbb{R}^{3}$ be such that
$\mathbf{r}\cdot\mathbf{n}\neq 1$ for all $\mathbf{n}\in\Omega$. The field $E$
created by function $E_{0}$ is presented by (15) and (16) at point
$\mathbf{r}$, where $\bar{x},\bar{y}$ and $\bar{z}$ are as in (6).

###### Proof.

Let us define $\epsilon>0$. Let us show first, that when the partition
$\Omega_{t}$ is small enough, in other words, the areas $|\Omega_{t}|$ are
small enough regardless of t, by replacing the function $E_{0}^{t}$ supported
by each piece $\Omega_{t}$ with a constant
$E_{0}^{t}(\bar{0},\bar{0})=E_{0}(\theta,\phi)$, the error made in formula
(10) is smaller than $\epsilon/4$. Function $E_{0}$ is uniformly continuous on
the compact segment of sphere $\Omega$, and based on this, with sufficiently
small $|\Omega_{t}|$ it holds

$\begin{array}[]{l}|E^{t}_{0}(\bar{x},\bar{y})-E_{0}(\theta,\phi)|\leq\frac{\epsilon}{4I|\Omega|}\end{array}$
(31)

for all $\bar{x},\bar{y}\in\Omega_{t}$ and $t$, where the constant $I$ will be
defined later on. In the Fourier transform (7) on the plane holds

$\begin{array}[]{l}|\mathcal{F}\left\\{E_{0}^{t}\right\\}-\mathcal{F}\left\\{E_{0}(\theta,\phi)\right\\}|\leq\iint_{\Omega_{t}}|E^{t}_{0}(\bar{x},\bar{y})-E_{0}(\theta,\phi)||e^{-i(k_{\bar{x}}\bar{x}+k_{\bar{y}}\bar{y})}|d\bar{x}d\bar{y}\\\
\leq\iint_{\Omega_{t}}\frac{\epsilon}{4I\Omega}|e^{-i(k_{\bar{x}}\bar{x}+k_{\bar{y}}\bar{y})}|d\bar{x}d\bar{y}\leq\iint_{\Omega_{t}}\frac{\epsilon}{4I|\Omega|}d\bar{x}d\bar{y}=\frac{\epsilon|\Omega_{t}|}{4I|\Omega|},\end{array}$
(32)

for all $(k_{\bar{x}},k_{\bar{y}})\in\mathbb{R}^{2}$. Function
$e^{i|\bar{z}|\sqrt{k^{2}-k_{\bar{x}}^{2}-k_{\bar{y}}^{2}}}$ is integrable on
$(k_{\bar{x}},k_{\bar{y}})$ plane, because outside of the circle
$k_{\bar{x}}^{2}+k_{\bar{y}}^{2}=k^{2}$, the exponential becomes real and
negative. Let $I$ be at first a continuous elementary function

$\begin{array}[]{l}I(|\bar{z}|)=\iint_{\mathbb{R}^{2}}|e^{i|\bar{z}|\sqrt{k^{2}-k_{\bar{x}}^{2}-k_{\bar{y}}^{2}}}|dk_{\bar{x}}dk_{\bar{y}},\end{array}$
(33)

where $\bar{z}=\bar{z}(\theta,\phi)$; see (6). A continuous function
$|\bar{z}(\theta,\phi)|$ gets its minimum value $\bar{z}_{min}>0$ on the
compact set
$\left\\{(\theta,\phi)|\theta\in[0,\theta_{0}],\phi\in[0,2\pi]\right\\}$.
Especially $|\bar{z}|\geq\bar{z}_{min}$ for all $t$ regardless of
partitioning. Finally, the constant $I$ is defined

$\begin{array}[]{l}I=I(z_{min})<\infty.\end{array}$ (34)

Due the monotone of integral, $I(|\bar{z}(\theta,\phi)|)\leq I$ on all $t$
regardless of partitions. The upper limit for the error due
$E_{0}(\theta,\phi)$ on the formula (12) is obtained as

$\begin{array}[]{l}\frac{1}{4\pi^{2}}\iint_{\mathbb{R}^{2}}|\mathcal{F}\left\\{E_{0}^{t}\right\\}-\mathcal{F}\left\\{E_{0}(\theta,\phi)\right\\}|(k_{\bar{x}},k_{\bar{y}})|e^{i|\bar{z}|\sqrt{k^{2}-k_{\bar{x}}^{2}-k_{\bar{y}}^{2}}}|dk_{\bar{x}}dk_{\bar{y}}\\\
\leq\frac{1}{4\pi^{2}}\frac{\epsilon|\Omega_{t}|}{4I|\Omega|}\iint_{\mathbb{R}^{2}}|e^{i|\bar{z}|\sqrt{k^{2}-k_{\bar{x}}^{2}-k_{\bar{y}}^{2}}}|dk_{\bar{x}}dk_{\bar{y}}\\\
\leq\frac{1}{4\pi^{2}}\frac{\epsilon|\Omega_{t}|}{4I|\Omega|}I\leq\frac{\epsilon|\Omega_{t}|}{4|\Omega|},\end{array}$
(35)

where (31) and (34) have been used. Because of $||\mathbf{e}_{\theta}||=1$ and
the triangle inequality the error in formula (13) is

$\begin{array}[]{l}err_{1}\leq\sum_{t}\frac{\epsilon|\Omega_{t}|}{4|\Omega|}=\frac{\epsilon}{4}.\end{array}$
(36)

In the above approximation, a constant $E_{0}(\theta,\phi)$ is used on piece
$\Omega$; as a function it is zero on the plane outside of the piece, let us
denote this function simply as $E_{0}(\theta,\phi)$. Next, we consider the
Fourier transform of function $E_{0}(\theta,\phi)$. Because function
$e^{i|\bar{z}|\sqrt{k^{2}-k_{\bar{x}}^{2}-k_{\bar{y}}^{2}}}$ is integrable,
there can be found an origin centered closed disk $\bar{B}$ such that

$\begin{array}[]{l}\left[max_{\Omega}|E_{0}|\right]\iint_{\mathbb{R}^{2}/\bar{B}}|e^{i|\bar{z}|\sqrt{k^{2}-k_{\bar{x}}^{2}-k_{\bar{y}}^{2}}}|dk_{\bar{x}}dk_{\bar{y}}<\frac{\epsilon}{4|\Omega|}\end{array}$
(37)

for all $t$ because $|\bar{z}|\geq\bar{z}_{min}$. For all
$(k_{\bar{x}},k_{\bar{y}})\in\mathbb{R}^{2}$ it holds

$\begin{array}[]{l}|\mathcal{F}\left\\{E_{0}(\theta,\phi)\right\\}(k_{\bar{x}},k_{\bar{y}})|\\\
\leq\left[max_{\Omega}|E_{0}|\right]\iint_{\Omega_{t}}1dk_{\bar{x}}dk_{\bar{y}}=max|E_{0}||\Omega_{t}|.\end{array}$
(38)

Let us now consider the error ($err_{2}$) made in integral (12) when
$\mathcal{F}\left\\{E_{0}(\theta,\phi)\right\\}$ is replaced by the constant
$E_{0}(\theta_{t},\phi_{t})|\Omega_{t}|$. Under estimates (37) and (38) it
follows

$\begin{array}[]{l}\iint_{\mathbb{R}^{2}/\bar{B}}|\mathcal{F}\left\\{E_{0}(\theta,\phi)\right\\}-E_{0}(\theta,\phi)|\Omega_{t}|||e^{i|\bar{z}|\sqrt{k^{2}-k_{\bar{x}}^{2}-k_{\bar{y}}^{2}}}|dk_{\bar{x}}dk_{\bar{y}}\\\
\leq\iint_{\mathbb{R}^{2}/\bar{B}}2\left[max_{\Omega}|E_{0}|\right]|\Omega_{t}||e^{i|\bar{z}|\sqrt{k^{2}-k_{\bar{x}}^{2}-k_{\bar{y}}^{2}}}|dk_{\bar{x}}dk_{\bar{y}}\\\
\leq
2|\Omega_{t}|\frac{\epsilon}{4|\Omega|}=\frac{\epsilon|\Omega_{t}|}{2|\Omega|}.\end{array}$
(39)

For all $(k_{\bar{x}},k_{\bar{y}})\in\bar{B}$ and
($\bar{x},\bar{y})\in\Omega_{t}$, we can estimate by continuity of the
function

$\begin{array}[]{l}\left[max_{\Omega}|E_{0}|\right]|e^{-i(k_{\bar{x}}\bar{x}+k_{\bar{y}}\bar{y})}-1|\leq\frac{\epsilon}{4|\bar{B}|\Omega|}\end{array}$
(40)

always when $|\Omega_{t}|$ is small enough, because then $\bar{x}\approx
0,\bar{y}\approx 0$ and hence $k_{\bar{x}}\bar{x}+k_{\bar{y}}\bar{y}\approx 0\
-$ that is the assignment of $\bar{B}$. Then in the Fourier transform at point
$(k_{\bar{x}},k_{\bar{y}})\in\bar{B}$ it holds

$\begin{array}[]{l}|\mathcal{F}\left\\{E_{0}(\theta,\phi)\right\\}(k_{\bar{x}},k_{\bar{y}})-E_{0}(\theta,\phi)|\Omega_{t}||\\\
=|\iint_{\Omega_{t}}E_{0}(\theta,\phi)e^{-i(k_{\bar{x}}\bar{x}+k_{\bar{y}}\bar{y})}d\bar{x}d\bar{y}-E_{0}(\theta,\phi)1d\bar{x}d\bar{y}|\\\
\leq|E_{0}(\theta,\phi)|\iint_{\Omega_{t}}|e^{-i(k_{\bar{x}}\bar{x}+k_{\bar{y}}\bar{y})}-1|d\bar{x}d\bar{y}\\\
\leq\iint_{\Omega_{t}}\frac{\epsilon}{4|\bar{B}||\Omega|}d\bar{x}d\bar{y}=\frac{\epsilon|\Omega_{t}|}{4|\bar{B}||\Omega|}.\end{array}$
(41)

Replacing by $\mathcal{F}\left\\{E_{0}(\theta,\phi)\right\\}$ with the
constant $E_{0}(\theta,\phi)|\Omega_{t}|$ in formula (12) we make an error

$\begin{array}[]{l}\frac{1}{4\pi^{2}}\Big{(}\iint_{\bar{B}}+\iint_{\mathbb{R}^{2}/\bar{B}}\Big{)}|\mathcal{F}\left\\{E_{0}(\theta,\phi)\right\\}-E_{0}(\theta,\phi)|\Omega_{t}||e^{i|\bar{z}|\sqrt{k^{2}-k_{\bar{x}}^{2}-k_{\bar{y}}^{2}}}|dk_{\bar{x}}dk_{\bar{y}}\\\
\leq\frac{\epsilon|\Omega_{t}|}{4|\bar{B}||\Omega_{t}|}\iint_{\bar{B}}1dk_{\bar{x}}dk_{\bar{y}}+\frac{\epsilon|\Omega_{t}|}{2|\Omega|}=\frac{\epsilon|\Omega_{t}|}{4|\Omega|}+\frac{\epsilon|\Omega_{t}|}{2|\Omega|}=\frac{3\epsilon|\Omega_{t}|}{4|\Omega|},\end{array}$
(42)

where the estimates (39) and (41) have been used. The error for sum (13) is
obtained as

$\begin{array}[]{l}err_{2}\leq\sum_{t}\frac{3\epsilon|\Omega_{t}|}{4|\Omega|}=\frac{3\epsilon}{4}.\end{array}$
(43)

Thus, when the partition $\Omega_{t}$ is dense enough and
$\mathbf{r}\in\mathbb{R}^{3}$ be such that $\mathbf{r}\cdot\mathbf{n}\neq 1$
for all $\mathbf{n}\in\Omega$, we get by (36) and (43)

$\begin{array}[]{l}\Big{|}\Big{|}\sum_{t}E_{t}(\mathbf{r})\mathbf{e}_{\theta}-\frac{1}{4\pi^{2}}\sum_{t}E_{0}(\theta,\phi)\mathbf{e}_{\theta}(\theta,\phi)|\Omega_{t}|\iint_{\mathbb{R}}e^{i(k_{\bar{x}}\bar{x}+k_{\bar{y}}\bar{y})}e^{i|\bar{z}|\sqrt{k^{2}-k_{\bar{x}}^{2}-k_{\bar{y}}^{2}}}dk_{\bar{x}}dk_{\bar{y}}\Big{|}\Big{|}\\\
\leq
err_{1}+err_{2}\leq\frac{\epsilon}{4}+\frac{3\epsilon}{4}=\epsilon.\end{array}$
(44)

Finally, by continuity of the integrand, the integral in (15) exists, and,
because of the estimate (44), for the points $\mathbf{r}\in\mathbb{R}^{3}$
like in Theorem there also exists

$\begin{array}[]{l}\mathbf{E}_{\theta}(\mathbf{r})=\mathbf{E}_{\theta}(x,y,z)=lim_{|\Omega_{t}|\to
0}\sum_{t}E_{t}(\mathbf{r})\mathbf{e}_{\theta}\\\
=\frac{1}{4\pi^{2}}\iint_{\Omega}E_{0}(\theta,\phi)\mathbf{E}_{t}(\mathbf{r};\theta,\phi)\sin{\theta}d\theta
d\phi.\end{array}$ (45)

∎

The proof above is also valid for the sum of the second vector component of
(19). Only a few modifications are needed.

##  6 Appendix B

The VSH beam shape coefficients for each $ij-$ plane wave of the incident
field is defined as [33]

$\begin{array}[]{l}a^{t}_{{\text{emn}}_{ij}}=4i^{n}\mathbf{e}_{o_{ij}}\cdot\big{[}-\mathbf{e}_{\theta}\sin{(m\phi)}\frac{m}{\sin{\theta}}P^{m}_{n}(\cos{\theta})-\mathbf{e}_{\phi}\cos{(m\phi)}\frac{d}{d\theta}P^{m}_{n}(\cos{\theta})\big{]},\\\
a^{t}_{{\text{omn}}_{ij}}=4i^{n}\mathbf{e}_{o_{ij}}\cdot\big{[}\mathbf{e}_{\theta}\cos{(m\phi)}\frac{m}{\sin{\theta}}P^{m}_{n}(\cos{\theta})-\mathbf{e}_{\phi}\sin{(m\phi)}\frac{d}{d\theta}P^{m}_{n}(\cos{\theta})\big{]},\\\
b^{t}_{{\text{emn}}_{ij}}=-4i^{n+1}\mathbf{e}_{o_{ij}}\cdot\big{[}\mathbf{e}_{\theta}\cos{(m\phi)}\frac{m}{\sin{\theta}}P^{m}_{n}(\cos{\theta})-\mathbf{e}_{\phi}\sin{(m\phi)}\frac{d}{d\theta}P^{m}_{n}(\cos{\theta})\big{]},\\\
b^{t}_{{\text{omn}}_{ij}}=-4i^{n+1}\mathbf{e}_{o_{ij}}\cdot\big{[}\mathbf{e}_{\theta}\sin{(m\phi)}\frac{m}{\sin{\theta}}P^{m}_{n}(\cos{\theta})+\mathbf{e}_{\phi}\cos{(m\phi)}\frac{d}{d\theta}P^{m}_{n}(\cos{\theta})\big{]},\end{array}$
(46)

where $P^{m}_{n}$ is the associated Legendre function of the first kind of
degree $n$ and order $m$. The scattered field from differential element in VSH
presentation is obtained as

$\begin{array}[]{l}\begin{aligned}
\mathbf{E}_{scat}(k\mathbf{r})&=p^{2}\sum_{m}\sum_{n}D_{mn}\Big{[}f_{\text{emn}}\mathbf{M}^{(3)}_{\text{emn}}(k\mathbf{\bar{r}})+f_{\text{omn}}\mathbf{M}^{(3)}_{\text{omn}}(k\mathbf{\bar{r}})\\\
&+g_{\text{emn}}\mathbf{N}^{(3)}_{\text{emn}}(k\mathbf{\bar{r}})+g_{\text{omn}}\mathbf{N}^{(3)}_{\text{omn}}(k\mathbf{\bar{r}})\Big{]},\end{aligned}\end{array}$
(47)

where the superscripts (3) presents the vectors spherical harmonics (outgoing
wave) with spherical Hankel function of the first kind
$h_{\text{n}}^{(1)}(kr)$. The $f_{\text{emn}},f_{\text{omn}},g_{\text{emn}}$
and $g_{\text{omn}}$ are vector spherical harmonic coefficient for the
scattered field calculated from the T-matrix method as

$\begin{array}[]{l}\begin{bmatrix}f_{\text{emn}}\\\ f_{omn}\\\
g_{\text{emn}}\\\ g_{omn}\end{bmatrix}=\begin{bmatrix}T_{11}&0&0&0\\\
0&T_{22}&0&0\\\ 0&0&T_{33}&0\\\
0&0&0&T_{44}\end{bmatrix}\begin{bmatrix}a_{\text{emn}}\\\ a_{omn}\\\
b_{\text{emn}}\\\ b_{omn}\end{bmatrix},\end{array}$ (48)

where the T-matrix elements for the multilayered sphere is obtained by the
algorithm defined in [34]. The coupling coefficient between the incident and
scattered field can be obtained as

$K=\frac{\int\int\big{[}\mathbf{E}_{inc}\cdot\mathbf{E}_{sca}\big{]}dxdy}{\int\int\big{[}\mathbf{E}_{inc}\cdot{\mathbf{E}_{inc}}^{*}\big{]}dxdy},$
(49)

where $\mathbf{E}_{inc}$ and $\mathbf{E}_{sca}$ are computed by Eq. (29) and
Eq. (30), and ∗ denotes the complex conjugate.

Funding This research was supported by Academy of Finland $\\#$327640:
Assessment of the Graft Rejection Using Millimeter Waves (AGRUM) and the
Agencia Estatal Investigación PID2019-107885GB-C31/AEI/10.13039,
PRE2018-084326.

Disclosures The authors declare no conflicts of interest.

Data availability Data underlying the results presented in this paper are not
publicly available at this time but may be obtained from the authors upon
reasonable request.

## References

  * [1] Y. Deng and X. Liu, “Electromagnetic imaging methods for nondestructive evaluation applications,” Sensors 11, 11774–11808 (2011).
  * [2] X. Liu and D. Zhao, “The statistical properties of anisotropic electromagnetic beams passing through the biological tissues,” Optics Communications 285, 4152 – 4156 (2012).
  * [3] P. Sowa, J. Rutkowska-Talipska, U. Sulkowska, K. Rutkowski, and R. Rutkowski, “Electromagnetic radiation in modern medicine: Physical and biophysical properties,” Polish Annals of Medicine 19, 139 – 142 (2012).
  * [4] Z. D. Taylor, J. Garritano, S. Sung, N. Bajwa, D. B. Bennett, B. Nowroozi, P. Tewari, J. Sayre, J. P. Hubschman, S. Deng, E. R. Brown, and W. S. Grundfest, “Thz and mm-wave sensing of corneal tissue water content: Electromagnetic modeling and analysis,” IEEE Transactions on Terahertz Science and Technology 5, 170–183 (2015).
  * [5] D. B. Bennett, Z. D. Taylor, P. Tewari, R. S. Singh, M. O. Culjat, W. S. G. M.D., D. J. Sassoon, R. D. Johnson, J.-P. Hubschman, and E. Brown, “Terahertz sensing in corneal tissues,” Journal of Biomedical Optics 16, 1 – 9 (2011).
  * [6] N. W. Bergmann, D. Abbott, A. Hariz, and V. K. Varadan, eds., _Electronics and Structures for MEMS II_ , vol. 4591 (SPIE, 2001).
  * [7] P. Y. Han, G. C. Cho, and X.-C. Zhang, “Time-domain transillumination of biological tissues with terahertz pulses,” Opt. Lett. 25, 242–244 (2000).
  * [8] A. Tamminen, S. V. Pälli, J. Ala-Laurinaho, M. Salkola, A. V. Räisänen, and Z. D. Taylor, “Quasioptical system for corneal sensing at 220–330 ghz: Design, evaluation, and ex vivo cornea parameter extraction,” IEEE Transactions on Terahertz Science and Technology 11, 135–149 (2021).
  * [9] F. Zarrinkhat, J. Lamberg, A. Tamminen, M. Baggio, J. Ala-Laurinaho, E. E. M. Khaled, J. M. Rius, J. R. Robert, and Z. Taylor, “Fourier analysis of submillimeter-wave scattering from the human cornea,” in _2021 15th European Conference on Antennas and Propagation (EuCAP),_ (2021), pp. 1–5.
  * [10] P. Yeh, A. Yariv, and C.-S. Hong, “Electromagnetic propagation in periodic stratified media. i. general theory$\ast$,” J. Opt. Soc. Am. 67, 423–438 (1977).
  * [11] G. A. Niklasson, C. G. Granqvist, and O. Hunderi, “Effective medium models for the optical properties of inhomogeneous materials,” Appl. Opt. 20, 26–30 (1981).
  * [12] M. Gharghi, H. Bai, G. Stevens, and S. Sivoththaman, “Modeling and simulation of spherical solar cells,” in _Conference Record of the Thirty-first IEEE Photovoltaic Specialists Conference, 2005._ , (2005), pp. 1177–1180.
  * [13] N. Keerativoranan, P. Hanpinitsak, K. Saito, and J. I. Tadaka, “Circular mesh based physical optics for scattered field prediction,” IEEE Antennas and Wireless Propagation Letters pp. 1–1 (2021).
  * [14] J. R. Sanford, “Scattering by spherically stratified microwave lens antennas,” IEEE Transactions on Antennas and Propagation 42, 690–698 (1994).
  * [15] G. Gouesbet, J. Lock, and G. Gréhan, “Generalized lorenz–mie theories and description of electromagnetic arbitrary shaped beams: Localized approximations and localized beam models, a review,” Journal of Quantitative Spectroscopy and Radiative Transfer 112, 1–27 (2011).
  * [16] C. Bohren and D. R. Huffman, _Absorption and Scattering of Light by Small Particles_ (Wiley Science Paperback Series, 1998).
  * [17] J. W. Goodman, _Introduction to Fourier optics_ (Englewood, CO: Roberts, 2005).
  * [18] E. E. M. Khaled, S. C. Hill, and P. W. Barber, “Scattered and internal intensity of a sphere illuminated with a gaussian beam,” IEEE Transactions on Antennas and Propagation 41, 295–303 (1993).
  * [19] C.-Y. Hwang, S. Oh, I.-K. Jeong, and H. Kim, “Stepwise angular spectrum method for curved surface diffraction,” Opt. Express 22, 12659–12667 (2014).
  * [20] L. Ebers, M. Hammer, and J. Förstner, “Light diffraction in slab waveguide lenses simulated with the stepwise angular spectrum method,” Opt. Express 28, 36361–36379 (2020).
  * [21] N. G. Worku, R. Hambach, and H. Gross, “Decomposition of a field with smooth wavefront into a set of gaussian beams with non-zero curvatures,” J. Opt. Soc. Am. A 35, 1091–1102 (2018).
  * [22] J.-S. Lee, T.-L. Song, J.-K. Du, and J.-G. Yook, “Near-field to far-field transformation based on stratton-chu fomula for emc measurements,” in _2013 IEEE Antennas and Propagation Society International Symposium (APSURSI),_ (2013), pp. 606–607.
  * [23] J. Lamberg, F. Zarrinkhat, A. Tamminen, E. E. M. Khaled, and Z. Taylor, “Investigation of radial top-hat electric field distributions for corneal reflectometry using modified fourier optics method,” in _2021 46th International Conference on Infrared, Millimeter and Terahertz Waves (IRMMW-THz),_ (2021), pp. 1–2.
  * [24] W. C. Chew, _Planarly Layered Media_ (1995), pp. 45–160.
  * [25] W. Hsu and R. Barakat, “Stratton–chu vectorial diffraction of electromagnetic fields by apertures with application to small-fresnel-number systems,” J. Opt. Soc. Am. A 11, 623–629 (1994).
  * [26] E. E. M. Khaled, S. C. Hill, and P. W. Barber, “Light scattering by a coated sphere illuminated with a gaussian beam,” Appl. Opt. 33, 3308–3314 (1994).
  * [27] L. N. Trefethen and J. A. C. Weideman, “The exponentially convergent trapezoidal rule,” SIAM Review 56, 385–458 (2014).
  * [28] M. A. Butt, S. A. Fomchenkov, and S. N. Khonina, “Multilayer dielectric stack notch filter for 450-700 nm wavelength spectrum,” (2017).
  * [29] F. Zarrinkhat, J. Lamberg, A. Tamminen, M. Baggio, I. Nefedova, J. Ala-Laurinaho, E. E. M. Khaled, J. Rius, J. Romeu, and Z. Taylor, “Vector spherical harmonic analysis and experimental validation of spherical shells illuminated with broadband, millimeter wave gaussian beams: applications to corneal sensing,” Biomed. Opt. Express 13, 3699–3722 (2022).
  * [30] F. Zarrinkhat, M. Baggio, J. Lamberg, A. Tamminen, I. Nefedova, J. Ala-Laurinaho, E. E. M. Khaled, J. M. Rius, J. Romeu, and Z. Taylor, “Calibration alignment sensitivity in corneal terahertz imaging,” Sensors 22 (2022).
  * [31] A. Tamminen, M. Baggio, I. Nefedova, Q. Sun, S. Presnyakov, J. Ala-laurinaho, E. Brown, V. Wallace, E. Macpherson, T. Maloney, N. Kravchenko, M. Salkola, S. Deng, and Z. Taylor, “Extraction of thickness and water content gradients in hydrogel-based, water-backed corneal phantoms via submillimeter wave reflectometry,” IEEE Transactions on Terahertz Science and Technology pp. 1–1 (2021).
  * [32] J. T. Kindt and C. A. Schmuttenmaer, “Far-infrared dielectric properties of polar liquids probed by femtosecond terahertz pulse spectroscopy,” The Journal of Physical Chemistry 100, 10373–10379 (1996).
  * [33] P. W. Barber and S. C. Hill, _Light Scattering By Particles: Computational Methods_ (Advanced Series in Applied Physics, 1990).
  * [34] U. P. O. Peña, “Scattering of electromagnetic radiation by a multilayered sphere,” IEEE Antennas and Propagation Magazine 57, 69–116 (2015).