Anisotropic homogenized composite mediums arising from truncated spheres,
spheroids, and ellipsoids

Tom G. Mackay***E–mail<EMAIL_ADDRESS>
School of Mathematics and Maxwell Institute for Mathematical Sciences
University of Edinburgh, Edinburgh EH9 3FD, UK
and
NanoMM — Nanoengineered Metamaterials Group
Department of Engineering Science and Mechanics
The Pennsylvania State University, University Park, PA 16802–6812, USA
Akhlesh Lakhtakia
NanoMM — Nanoengineered Metamaterials Group
Department of Engineering Science and Mechanics
The Pennsylvania State University, University Park, PA 16802–6812, USA

Abstract

Closed-form expressions were recently derived for depolarization dyadics for
truncated spheres and truncated spheroids, and the formalism was extended to
truncated ellipsoids. These results were exploited to develop an
implementation of the Maxwell Garnett homogenization formalism for the
relative permittivity parameters of homogenized composite mediums (HCMs)
arising from an isotropic host medium impregnated with isotropic inclusions
that are truncated spheres, spheroids, and ellipsoids. In so doing, the
anisotropy of the HCM was related to the geometry of the inclusions: in
general, the more the shape of the inclusions deviated from spherical, the
greater was the degree of anisotropy exhibited by the HCM.

Keywords: Maxwell Garnett homogenization formalism; depolarization dyadics;
anisotropy; truncated spheres; truncated spheroids; truncated ellipsoids

##  1 Introduction

A host medium impregnated with particulate inclusions may be regarded as a
homogenized composite medium (HCM) provided that the inclusions are much
smaller than the wavelengths involved [1, 2, 3, 4]. Current interest in HCMs
is intense due to the burgeoning development of nanocomposite materials for
optical applications [5]. In particular, HCMs can exhibit properties not
exhibited by the component mediums from which they arise [6]. For example,
anisotropic HCMs can rise from isotropic component mediums if the geometries
of the inclusions exhibit certain symmetries [7].

Homogenization formalisms are used to estimate the constitutive parameters of
HCMs [1]. The Maxwell Garnett formalism [8] is perhaps the most widely used of
all homogenization formalisms. Most homogenization formalisms (including the
Maxwell Garnett formalism) use _depolarization dyadics_ to represent the
scattering responses of the inclusions [4], with the depolarization dyadic
being the integrated singularity of the corresponding dyadic Green function
[9, 10]. Generally, homogenization formalisms are used for inclusions with
simple shapes, e.g., spherical, spheroidal, and ellipsoidal. This is because
closed-form expressions for depolarization dyadics are only available for a
limited range of inclusion shapes [11, 12, 13, 14, 15, 16, 17], with more
complex inclusion shapes requiring numerical evaluation [18, 19].

Recently, the range of inclusion shapes that can be accommodated in
homogenization formalisms was extended considerably by the derivation of
closed-form expressions for depolarization dyadics for truncated spheres and
truncated spheroids. Furthermore, the methodology was extended to truncated
ellipsoids [20]. In this paper, these recently derived expressions for
depolarization dyadics are used in an implementation of the Maxwell Garnett
formalism to estimate the relative permittivity parameters of HCMs arising
from inclusions shaped as truncated spheres and spheroids. The implementation
is extended to the case of inclusions with truncated ellipsoidal shapes, with
the corresponding depolarization dyadics being evaluated numerically. As
regards notational matters: vectors are boldface; dyadics [4, 10] are double
underlined; $\varepsilon_{\scriptscriptstyle 0}$ is the free-space
permittivity; $\underline{\underline{I}}=\hat{{\bf u}}_{\rm x}\hat{{\bf
u}}_{\rm x}+\hat{{\bf u}}_{\rm y}\hat{{\bf u}}_{\rm y}+\hat{{\bf u}}_{\rm
z}\hat{{\bf u}}_{\rm z}$ is the identity dyadic with $\hat{{\bf u}}_{\rm x}$,
$\hat{{\bf u}}_{\rm y}$, and $\hat{{\bf u}}_{\rm z}$ being the unit vectors
aligned with the coordinate axes of the Cartesian coordinate system $(x,y,z)$.

##  2 Homogenization – theoretical matters

Consider a composite medium comprising a random distribution of identical,
similarly oriented, electromagnetically small inclusions surrounded by a host
medium. The inclusion medium has relative permittivity scalar
${\varepsilon_{\rm in}}$, and the host medium ${\varepsilon_{\rm h}}$. The
volume fraction of the composite material occupied by the inclusions is
denoted by $f_{\rm in}$.

The Maxwell Garnett formalism delivers the following estimate of the relative
permittivity dyadic of homogenized composite material:

$\underline{\underline{\varepsilon}}^{MG}={\varepsilon_{\rm
h}}\,\underline{\underline{I}}+{\frac{f_{\rm
in}}{\varepsilon_{\scriptscriptstyle 0}}}\,\underline{\underline{a}}\mbox{
\tiny{${}^{\bullet}$}
}\left(\underline{\underline{I}}-\frac{1}{3\varepsilon_{\scriptscriptstyle
0}\varepsilon_{\rm h}}f_{\rm in}\,\underline{\underline{a}}\right)^{-1}.$ (1)

Herein the polarizability density dyadic

$\underline{\underline{a}}=\varepsilon_{\scriptscriptstyle
0}\left(\varepsilon_{\rm in}-\varepsilon_{\rm
h}\right)\left[\underline{\underline{I}}+\frac{\varepsilon_{\rm
in}-\varepsilon_{\rm h}}{\varepsilon_{\rm h}}\underline{\underline{L}}({\bf
r}_{\rm o})\right]^{-1};$ (2)

and the depolarization dyadic

$\underline{\underline{L}}({\bf r}_{\rm o})=\frac{1}{4\pi}\iint_{S_{\rm
in}}\hat{{\bf u}}_{\rm n}({\bf r}^{\prime})\frac{{\bf r}^{\prime}-{\bf r}_{\rm
o}}{|{\bf r}^{\prime}-{\bf r}_{\rm o}|^{3}}\,d^{2}{\bf r}^{\prime}$ (3)

is an integral over the surface $S_{\rm in}$ of an inclusion, with $\hat{{\bf
u}}_{\rm n}$ being the unit outward normal vector to $S_{\rm in}$ and ${\bf
r}_{\rm o}$ being the position vector of a point inside the inclusion. Details
leading to Eq. (1) are provided in the Supplementary Material.

##  3 Numerical results and discussion

We used Eq. (1) to calculate the results presented in this section. Since
$\underline{\underline{L}}({\bf r}_{\rm o})$ is diagonalizable, both
$\underline{\underline{a}}$ and $\underline{\underline{\varepsilon}}^{MG}$
also are. In addition, all three dyadics have the same set of eigenvectors.
This fact can be used to obtain simple expressions for the three eigenvalues
of $\underline{\underline{\varepsilon}}^{MG}$, if desired.

The inclusion shapes were taken to be that of truncated spheres, spheroids,
and ellipsoids. For the cases of truncated spheres and spheroids, closed-form
expressions for $\underline{\underline{L}}({\bf r}_{\rm o})$ were recently
derived [20]; these cases are represented schematically in Fig. 1. For each
type of inclusion shown in this figure, $\underline{\underline{L}}({\bf
r}_{\rm o})$ was derived with ${\bf r}_{\rm o}$ specifying the center of the
largest sphere inscribed inside the inclusion.

Without loss of generality, the rotational symmetry axis of the truncated
spheres, spheroids, and ellipsoids was fixed as the $z$ coordinate axis; i.e.,
all truncation planes were parallel to the $xy$ plane. The three eigenvectors
of $\underline{\underline{L}}\equiv\underline{\underline{L}}({\bf r}_{\rm o})$
are $\hat{{\bf u}}_{\rm x}$, $\hat{{\bf u}}_{\rm y}$, and $\hat{{\bf u}}_{\rm
z}$ for all inclusions due to the chosen shapes. Accordingly,

$\left.\begin{array}[]{l}\underline{\underline{L}}={L_{\rm x}}\,\hat{{\bf
u}}_{\rm x}\hat{{\bf u}}_{\rm x}+{L_{\rm y}}\,\hat{{\bf u}}_{\rm y}\hat{{\bf
u}}_{\rm y}+{L_{\rm z}}\,\hat{{\bf u}}_{\rm z}\hat{{\bf u}}_{\rm
z}\vspace{6pt}\\\
\underline{\underline{\varepsilon}}^{MG}=\varepsilon^{MG}_{\rm x}\,\hat{{\bf
u}}_{\rm x}\hat{{\bf u}}_{\rm x}+\varepsilon^{MG}_{\rm y}\,\hat{{\bf u}}_{\rm
y}\hat{{\bf u}}_{\rm y}+\varepsilon^{MG}_{\rm z}\,\hat{{\bf u}}_{\rm
z}\hat{{\bf u}}_{\rm z}\end{array}\right\\};$ (4)

with ${L_{\rm x}}={L_{\rm y}}\equiv L_{\rm t}$ and $\varepsilon^{MG}_{\rm
x}=\varepsilon^{MG}_{\rm y}\equiv\varepsilon^{MG}_{\rm t}$ for truncated
spheres and truncated spheroids.

For the purposes of numerical illustration, the representative values
${\varepsilon_{\rm h}}=3$, ${\varepsilon_{\rm in}}=2+0.5i$, and $f_{\rm
in}=0.3$ were selected, with higher values of $f_{\rm in}$ being incompatible
with the Maxwell Garnett formalism [4, 1]. All inclusions are taken to be
aligned in parallel, but orientational statistics can be easily incorporated
as needed [21].

###  3.1 Inclusions based on spherical geometry

The inclusion shape is specified relative to the unit sphere, centered at the
coordinate origin.

####  3.1.1 Truncated sphere

The inclusion shape is that of the upper part of the sphere which is bounded
below by the plane $z=1-2\kappa$, where $0<\kappa<1$, as schematically
illustrated in Fig. 1(a). Thus, $\kappa$ represents the radius of the largest
sphere that can be inscribed inside the truncated sphere. In this case,
${L_{\rm z}}=1-2L_{\rm t}$ and

$L_{\rm
t}=\frac{\kappa}{6\left(1-\kappa\right)^{3}\tau}\,\left\\{6-\tau\left(3-3\kappa+\kappa^{2}\right)-\kappa\left[11-3\kappa\left(3-\kappa\right)\right]\right\\}\,,$
(5)

where

$\tau=\sqrt{\left(4-3\kappa\right)\kappa}\,.$ (6)

Plots of the real and imaginary parts of the relative permittivity parameters
$\varepsilon^{MG}_{\rm t}$ and $\varepsilon^{MG}_{\rm z}$ of the HCM versus
$\kappa$ are provided in Fig. 2. As $\kappa$ increases, the values of
$\mbox{Re}\left\\{\varepsilon^{MG}_{\rm t}\right\\}$ and
$\mbox{Im}\left\\{\varepsilon^{MG}_{\rm z}\right\\}$ decrease monotonically
whereas the values of $\mbox{Re}\left\\{\varepsilon^{MG}_{\rm z}\right\\}$ and
$\mbox{Im}\left\\{\varepsilon^{MG}_{\rm t}\right\\}$ increase monotonically.
For all $\kappa<1$, $\mbox{Re}\left\\{\varepsilon^{MG}_{\rm
t}\right\\}>\mbox{Re}\left\\{\varepsilon^{MG}_{\rm z}\right\\}$ and
$\mbox{Im}\left\\{\varepsilon^{MG}_{\rm
z}\right\\}>\mbox{Im}\left\\{\varepsilon^{MG}_{\rm t}\right\\}$. The
differences between $\mbox{Re}\left\\{\varepsilon^{MG}_{\rm t}\right\\}$ and
$\mbox{Re}\left\\{\varepsilon^{MG}_{\rm z}\right\\}$ and between
$\mbox{Im}\left\\{\varepsilon^{MG}_{\rm t}\right\\}$ and
$\mbox{Im}\left\\{\varepsilon^{MG}_{\rm z}\right\\}$ decrease monotonically as
$\kappa$ increases, and these differences vanish in the limit $\kappa\to 1$.
That is, the degree of anisotropy that the HCM exhibits decreases
monotonically as $\kappa$ increases and the HCM becomes isotropic in the limit
$\kappa\to 1$.

####  3.1.2 Double-truncated sphere

The inclusion shape is that of the middle part of the sphere which is bounded
below by the plane $z=-\kappa$ and bounded above by the plane $z=\kappa$,
where $0<\kappa<1$, as schematically illustrated in Fig. 1(b). Thus, $\kappa$
represents the radius of the largest sphere that can be inscribed inside the
double-truncated sphere. In this case, ${L_{\rm z}}=1-2L_{\rm t}$ with

$L_{\rm t}=\frac{\left(3-\kappa^{2}\right)\kappa}{6}\,.$ (7)

In Fig. 3, the real and imaginary parts of the relative permittivity
parameters $\varepsilon^{MG}_{\rm t}$ and $\varepsilon^{MG}_{\rm z}$ of the
HCM are plotted against $\kappa$. As is the case for the truncated sphere
presented in §3.1.1, for the double-truncated sphere the HCM is isotropic in
the limit $\kappa\to 1$ and exhibits an increasing degree of anisotropy as
$\kappa$ decreases. The plots of the real and imaginary parts of
$\varepsilon^{MG}_{\rm t}$ and $\varepsilon^{MG}_{\rm z}$ for the double-
truncated sphere in Fig. 3 are similar to the corresponding plots in Fig. 2
for the truncated sphere. Differences are most conspicuous in the regime of
small values of $\kappa$ wherein the real and imaginary parts of
$\varepsilon^{MG}_{\rm t}$ and $\varepsilon^{MG}_{\rm z}$ change less rapidly
as $\kappa$ increases for the double-truncated sphere.

###  3.2 Inclusions based on spheroidal geometry

The inclusion shape is specified relative to the spheroid

$\frac{x^{2}+y^{2}}{\alpha^{2}}+z^{2}\leq 1,$ (8)

with equatorial radius $\alpha>0$.

####  3.2.1 Hemispheroid

The inclusion shape is that of the upper half of the spheroid; this shape lies
between the planes $z=0$ and $z=1$, as schematically illustrated in Fig. 1(c).
The radius of the largest sphere that can be inscribed inside the hemispheroid
is $\kappa$; we have $\kappa=\alpha\sqrt{1-\alpha^{2}}\in(0,1/2)$ for
$\alpha<1/\sqrt{2}$ while $\kappa=1/2$ for $\alpha>1/\sqrt{2}$. For
$\alpha<1/\sqrt{2}$, we have

$\displaystyle L_{\rm t}$ $\displaystyle=$
$\displaystyle\frac{1}{4\nu\gamma^{4}}\left\\{1+\nu+\alpha^{2}\left[\alpha\left(\alpha+\sqrt{2-\alpha^{2}}\right)-\nu-3\right]\right.$
(9)
$\displaystyle\left.-\left(\alpha\gamma\right)^{2}\sqrt{2-\alpha^{2}}\log\displaystyle{\frac{\left(1-\alpha^{2}+\gamma\right)\left(1-\alpha\right)}{\alpha\left(\nu-\gamma\right)}}\right\\},$

where

$\left.\begin{array}[]{l}\nu=\sqrt{2-3\alpha^{2}+\alpha^{4}}\vspace{6pt}\\\
\gamma=\sqrt{1-\alpha^{2}}\end{array}\right\\};$ (10)

for $\alpha>1/\sqrt{2}$, we have

$\displaystyle L_{\rm t}$ $\displaystyle=$
$\displaystyle\displaystyle{\frac{3\left(1+\sqrt{1+4\alpha^{2}}\right)+\alpha^{2}\left(10-8\sqrt{1+4\alpha^{2}}\right)-8\alpha^{4}}{4\gamma^{2}\left(3+8\alpha^{2}-16\alpha^{4}\right)}}$
(11)
$\displaystyle\displaystyle{+\frac{\alpha^{2}\log\displaystyle{\frac{1-2\alpha^{2}+\gamma}{\sqrt{1+3\alpha^{2}-4\alpha^{4}}-1}}}{4\gamma^{3}}}\,.$

For all $\alpha$, ${L_{\rm z}}=1-2L_{\rm t}$.

Plots of the real and imaginary parts of the relative permittivity parameters
$\varepsilon^{MG}_{\rm t}$ and $\varepsilon^{MG}_{\rm z}$ of the HCM versus
$\alpha$ are displayed in Fig. 4. As $\alpha$ increases, the values of
$\mbox{Re}\left\\{\varepsilon^{MG}_{\rm z}\right\\}$ and
$\mbox{Im}\left\\{\varepsilon^{MG}_{\rm t}\right\\}$ decrease monotonically,
whereas the values of $\mbox{Re}\left\\{\varepsilon^{MG}_{\rm t}\right\\}$ and
$\mbox{Im}\left\\{\varepsilon^{MG}_{\rm z}\right\\}$ increase monotonically.
The value $\alpha=0.569$ is noteworthy because for this value
$\varepsilon^{MG}_{\rm t}=\varepsilon^{MG}_{\rm z}$ and the HCM is
consequently isotropic. As $\alpha$ deviates farther from 0.569, the degree of
anisotropy exhibited by the HCM increases. Results for hemispherical
inclusions are obtained for $\alpha=1$.

####  3.2.2 Double-truncated spheroid

The inclusion shape is that of the middle part of the spheroid which is
bounded below by the plane $z=-\kappa$ and bounded above by the plane
$z=\kappa$, with $0<\kappa<1$, as schematically illustrated in Fig. 1(d). The
radius of the largest sphere that can be inscribed inside the double-truncated
spheroid is $\alpha$ for $\kappa>\alpha$ and the radius is $\kappa$ for
$\kappa<\alpha$. In this case

$L_{\rm
t}=\frac{-\nu+\alpha^{2}\tan^{-1}\nu}{2\left(\alpha^{2}-1\right)^{3/2}}$ (12)

with

$\nu=\kappa\sqrt{\frac{\alpha^{2}-1}{\alpha^{2}-\kappa^{2}\left(\alpha^{2}-1\right)}};$
(13)

and again ${L_{\rm z}}=1-2L_{\rm t}$.

In Fig. 5, the real and imaginary parts of the relative permittivity
parameters $\varepsilon^{MG}_{\rm t}$ and $\varepsilon^{MG}_{\rm z}$ of the
HCM are plotted against $\alpha$ and $\kappa$. As $\alpha$ increases,
$\mbox{Re}\left\\{\varepsilon^{MG}_{\rm t}\right\\}$ and
$\mbox{Im}\left\\{\varepsilon^{MG}_{\rm z}\right\\}$ increase monotonically,
and $\mbox{Re}\left\\{\varepsilon^{MG}_{\rm z}\right\\}$ and
$\mbox{Im}\left\\{\varepsilon^{MG}_{\rm t}\right\\}$ decrease monotonically,
for all values of $\kappa$. As $\kappa$ increases,
$\mbox{Re}\left\\{\varepsilon^{MG}_{\rm t}\right\\}$ and
$\mbox{Im}\left\\{\varepsilon^{MG}_{\rm z}\right\\}$ decrease monotonically,
and $\mbox{Re}\left\\{\varepsilon^{MG}_{\rm z}\right\\}$ and
$\mbox{Im}\left\\{\varepsilon^{MG}_{\rm t}\right\\}$ increase monotonically,
for all values of $\alpha$. In the limit $\kappa\to 0$, the values of
$\varepsilon^{MG}_{\rm t}$ and $\varepsilon^{MG}_{\rm z}$ are largely
independent of $\alpha$, except at the smallest values of $\alpha$.

###  3.3 Inclusions based on ellipsoidal geometry

The inclusion shape is specified relative to the ellipsoid

$\frac{x^{2}}{\alpha^{2}}+\frac{y^{2}}{\beta^{2}}+z^{2}\leq 1,$ (14)

with semi-axis lengths $\alpha>0$ and $\beta>0$.

####  3.3.1 Hemi-ellipsoid

The inclusion shape is that of the upper half of the ellipsoid; this shape
lies between the planes $z=0$ and $z=1$. The largest sphere inscribed inside
the hemi-ellipsoid has radius

$\kappa=\left\\{\begin{array}[]{lcr}1/2,&\mbox{for}&\alpha>1/\sqrt{2},\quad\beta>1/\sqrt{2}\vspace{8pt}\\\
\alpha\sqrt{1-\alpha^{2}}&\mbox{for}&\alpha\leq
1/\sqrt{2},\quad\alpha<\beta\vspace{8pt}\\\
\beta\sqrt{1-\beta^{2}}&\mbox{for}&\beta\leq
1/\sqrt{2},\quad\beta<\alpha\end{array}\right..$ (15)

The dyadic $\underline{\underline{L}}$ has to be determined by numerical
integration.

In Fig. 6, plots of the real and imaginary parts of the relative permittivity
parameters $\varepsilon^{MG}_{\rm x}$, $\varepsilon^{MG}_{\rm y}$, and
$\varepsilon^{MG}_{\rm z}$ of the HCM versus $\alpha$ and $\beta$ are
displayed. The real and imaginary parts of each of $\varepsilon^{MG}_{\rm x}$,
$\varepsilon^{MG}_{\rm y}$, and $\varepsilon^{MG}_{\rm z}$ vary smoothly as
$\alpha$ and $\beta$ increase. In particular,
$\mbox{Re}\left\\{\varepsilon^{MG}_{\rm x}\right\\}$ increases markedly and
$\mbox{Im}\left\\{\varepsilon^{MG}_{\rm x}\right\\}$ decreases markedly as
$\alpha$ increases but both are relatively insensitive to changes in $\beta$;
$\mbox{Re}\left\\{\varepsilon^{MG}_{\rm y}\right\\}$ increases markedly and
$\mbox{Im}\left\\{\varepsilon^{MG}_{\rm y}\right\\}$ decreases markedly as
$\beta$ increases but both are relatively insensitive to changes in $\alpha$;
and $\mbox{Re}\left\\{\varepsilon^{MG}_{\rm z}\right\\}$ decreases markedly as
both $\alpha$ and $\beta$ increase whereas
$\mbox{Im}\left\\{\varepsilon^{MG}_{\rm z}\right\\}$ increases markedly as
both $\alpha$ and $\beta$ increase.

####  3.3.2 Double-truncated ellipsoid

The inclusion shape is that of the middle part of the ellipsoid which is
bounded below by the plane $z=-\kappa$ and bounded above by the plane
$z=\kappa$, with $0<\kappa<1$. Again, $\underline{\underline{L}}$ has to be
determined by numerical integration.

For $\kappa=0.3$, plots of the real and imaginary parts of the relative
permittivity parameters $\varepsilon^{MG}_{\rm x}$, $\varepsilon^{MG}_{\rm
y}$, and $\varepsilon^{MG}_{\rm z}$ of the HCM versus $\alpha$ and $\beta$ are
provided in Fig. 7. The plots in Fig. 7 for the double-truncated ellipsoid are
qualitatively similar to those in Fig. 6 for the hemi-ellipsoid. The effects
of varying $\kappa$ are most appreciable at low values of $\alpha$ for both
$\mbox{Re}\left\\{\varepsilon^{MG}_{\rm x}\right\\}$ and
$\mbox{Im}\left\\{\varepsilon^{MG}_{\rm x}\right\\}$, at low values of $\beta$
for both $\mbox{Re}\left\\{\varepsilon^{MG}_{\rm y}\right\\}$ and
$\mbox{Im}\left\\{\varepsilon^{MG}_{\rm y}\right\\}$, and at low values of
both $\alpha$ and $\beta$ for both $\mbox{Re}\left\\{\varepsilon^{MG}_{\rm
z}\right\\}$ and $\mbox{Im}\left\\{\varepsilon^{MG}_{\rm z}\right\\}$. From
further numerical studies (not presented graphically here), the values of
$\varepsilon^{MG}_{\rm x}$, $\varepsilon^{MG}_{\rm y}$, and
$\varepsilon^{MG}_{\rm z}$ are generally much more sensitive to variations in
$\alpha$ and $\beta$ than they are to variations in $\kappa\in(0,1)$.

##  4 Closing remarks

In this paper, the Maxwell Garnett formalism was implemented to estimate the
relative permittivity parameters of HCMs arising from an isotropic host medium
impregnated with isotropic inclusions of relatively complex shapes. The
inclusion shapes considered were truncated spheres, spheroids, and ellipsoids.
In so doing, the anisotropy of the HCM was related to the geometry of the
inclusions: in general, the more the shape of the inclusions deviates from
spherical, the greater is the degree of anisotropy exhibited by the HCM. The
availability of closed-form expressions for depolarization dyadics for
inclusions shaped as truncated spheres and spheroids [20] increases
considerably the reach of the Maxwell Garnett formalism, allowing more complex
HCMs to be theoretically analyzed. Additionally, this formalism can be used
for more realistic modeling of poultry dust [22], house dust [23, 24], cements
[25], volcanic ash [26, 27], mineral dust [28, 29], and interstellar dust
[30], which are presently modeled very simply as agglomerations of complete
ellipsoids [31, 32]. Although we considered all inclusions in a composite
medium to be identical and aligned in parallel, both of those restrictions can
be easily surmounted as shown elsewhere [21, 33].

The Maxwell Garnett formalism is usually implemented for cases involving
spherical, spheroidal, and ellipsoidal inclusions. In such cases, the
derivations of the corresponding polarizability density dyadics rest upon the
assumption that electromagnetic field within each inclusion is spatially
uniform. This assumption is physically reasonable provided that the inclusions
are electrically small. In the preceding sections, this assumption was
extended to truncated spheres, truncated spheroids, and truncated ellipsoids.
For such truncated inclusions, this assumption is also physically reasonable
provided that the inclusions are sufficiently small and that their shapes do
not deviate too much from non-truncated spheres, spheroids, and ellipsoids.

Finally, all homogenization formalisms – including the Maxwell Garnett
formalism – involve various assumptions and approximations. The acid test of
any homogenization formalism can only come from comparison with suitable
experimental data.

## Acknowledgments

TGM was supported by EPSRC (grant number EP/V046322/1). AL was supported by
the Charles Godfrey Binder Endowment at Penn State.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Data availability

All data analyzed in this paper are included in this paper.

## References

  * [1] A. Lakhtakia (Ed.), _Selected Papers on Linear Optical Composite Materials_ (SPIE, 1996).
  * [2] P. S. Neelakanta, Handbook of Electromagnetic Materials: Monolithic and Composite Versions and Their Applications (CRC Press, 1995).
  * [3] T. C. Choy, _Effective Medium Theory: Principles and Applications_ , _2nd Ed._ (Oxford, 2016).
  * [4] T. G. Mackay and A. Lakhtakia, _Modern Analytical Electromagnetic Homogenization with Mathematica_ , _2nd Ed._ (IOP Publishing, 2020).
  * [5] M. P. Mengüç and M. Francoeur (Eds.), _Light, Plasmonics and Particles_ (Elsevier, 2023).
  * [6] A. Lakhtakia and T. G. Mackay, “Meet the metamaterials,” Opt. Photon. News 6(1), 32–39 (2007).
  * [7] T. G. Mackay and W. S. Weiglhofer, “Homogenization of biaxial composite materials: bianisotropic properties,” J. Opt. A: Pure Appl. Opt. 3, 45–52 (2001).
  * [8] J. C. Maxwell Garnett, “Colours in metal glasses and thin films,” Philos. Trans. R. Soc. Lond. A 203, 385–420 (1904).
  * [9] C.-T. Tai, _Dyadic Green Functions in Electromagnetic Theory, 2nd Ed_. Piscataway, NJ, USA: IEEE Press, 1994.
  * [10] M. Faryad and A. Lakhtakia, _Infinite-Space Dyadic Green Functions in Electromagnetism_. San Rafael, CA, USA: Morgan & Claypool, 2018.
  * [11] B. Michel, “A Fourier space approach to the pointwise singularity of an anisotropic dielectric medium,” Int. J. Appl. Electromagn. Mech. 8, 219–227 (1997).
  * [12] A. Moroz, “Depolarization field of spheroidal particles,” J. Opt. Soc. Am. B 26, 517–527 (2009).
  * [13] J. A. Osborn, “Demagnetizing factors of the general ellipsoid,” Phys. Rev. 67, 351–357 (1945).
  * [14] E. C. Stoner, “The demagnetizing factors for ellipsoids,” Phil. Mag. 36, 803–821 (1945).
  * [15] B. Michel and W. S. Weiglhofer, “Pointwise singularity of dyadic Green function in a general bianisotropic medium,” Arch. Elektr. Übertrag. 51, 219–223 (1997). Corrections: 52, 310 (1998).
  * [16] W. S. Weiglhofer, “Electromagnetic depolarization dyadics and elliptic integrals,” J. Phys. A: Math. Gen. 31, 7191–7196 (1998).
  * [17] W. S. Weiglhofer and T. G. Mackay, “Needles and pillboxes in anisotropic mediums,” IEEE Trans. Antennas Propagat. 50, 85–86 (2002).
  * [18] S.-W. Lee, J. Boersma, C.-L. Law, and G. A. Deschamps, “Singularity in Green’s function and its numerical evaluation,” IEEE Trans. Antennas. Propagat. 28, 311–317 (1980).
  * [19] A. Lakhtakia and N. S. Lakhtakia, “A procedure for evaluating depolarization dyadics of polyhedra,” Optik 109, 140–142 (1998).
  * [20] T. G. Mackay and A. Lakhtakia, “Depolarization dyadics for truncated spheres, spheroids, and ellipsoids,” https://doi.org/10.48550/arXiv.2312.08674
  * [21] A. Lakhtakia, “Frequency-dependent continuum electromagnetic properties of a gas of scattering centers,” Adv. Chem. Phys. 85(2), 311–359 (1993).
  * [22] M. Ahaduzzaman, L. Milan, C. L. Morton, P. F. Gerber, and S. W. Walkden-Brown, “Characterization of poultry house dust using chemometrics and scanning electron microscopy imaging,” Poultry Sci. 100, 101188 (2021).
  * [23] C. Lanzerstorfer, “Variations in the composition of house dust by particle size,” J. Environ. Sci. Health A 52, 770–777 (2017).
  * [24] X. Liu, “Understanding semi-volatile organic compounds in indoor dust,” Indoor Built Environment 31, 291–298 (2022).
  * [25] P. Stutzman, “Scanning electron microscopy imaging of hydraulic cement microstructure,” Cement Concrete Compos. 26, 957–966 (2004).
  * [26] A. Vogel, S. Diplas, A. J. Durant, A. S. Azar, M. F. Sunding, W. I. Rose, A. Sytchkova, C. Bonadonna, K. Krüger, and A. Stohl, “Reference data set of volcanic ash physicochemical and optical properties,” J. Geophys. Res.: Atmospheres 122, 9485–9514 (2017).
  * [27] J. I. Farquharson, H. Tuffen, F. B. Wadsworth, J. M. Castro, H. Unwin, and C. I. Schipper, “In-conduit capture of sub-micron volcanic ash particles via turbophoresis and sintering,” Nature Commun. 13, 4713 (2022).
  * [28] P. Knippertz and J.-B. W. Stuut (Eds), Mineral Dust: A Key Player in the Earth System (Springer, 2014).
  * [29] S. Merikallio, T. Nousiainen, M. Kahnert, and A.-M. Harri, “Light scattering by the Martian dust analog, palagonite, modeled with ellipsoids,” Opt. Express 21, 17972–17985 (2013).
  * [30] D. Rogantini, E. Costantini, S. T. Zeegers, C. P. de Vries, W. Bras, F. de Groot, H. Mutschke, and L. B. F. M. Waters, “Investigating the interstellar dust through the Fe K-edge,” Astronom. Astrophys. 609, A22 (2018).
  * [31] R. Stognienko, T. Henning, and V. Ossenkopf, “Optical properties of fluffy particles,” ASP Conf. Series 104, 427–431 (1996).
  * [32] T. Nousiainen, “Impact of particle shape on refractive-index dependence of scattering in resonance domain,” J. Quant. Spectrosc. Radiat. Transf. 108, 464–473 (2007).
  * [33] T. G. Mackay and A. Lakhtakia, “Percolation thresholds in the homogenization of spheroidal particles oriented in two directions,” Opt. Commun. 259, 727–737 (2006).

(a) Truncated unit sphere with largest inscribed sphere of radius
$\kappa\in(0,1)$.

(b) Double-truncated unit sphere with largest inscribed sphere of radius
$\kappa\in(0,1)$.

(c) Hemispheroid with largest inscribed sphere of radius $\kappa\in(0,1/2)$
for $\alpha<1/\sqrt{2}$ (top) and of radius $\kappa=1/2$ for
$\alpha>1/\sqrt{2}$ (bottom).

(d) Double-truncated spheroid with largest inscribed sphere of radius $\alpha$
for $\alpha<\kappa$ (top) and of radius $\kappa$ for $\alpha>\kappa$ (bottom).

Figure 1:  Schematic representations of truncated spheres and spheroids.

Figure 2:  Real and imaginary parts of $\varepsilon^{MG}_{\rm t}$ and
$\varepsilon^{MG}_{\rm z}$ plotted against $\kappa\in(0,1)$ for truncated
spherical inclusions.

Figure 3:  Real and imaginary parts of $\varepsilon^{MG}_{\rm t}$ and
$\varepsilon^{MG}_{\rm z}$ plotted against $\kappa\in(0,1)$ for double-
truncated spherical inclusions.

Figure 4:  Real and imaginary parts of $\varepsilon^{MG}_{\rm t}$ and
$\varepsilon^{MG}_{\rm z}$ plotted against $\alpha\in(0,5)$ for hemispheroidal
inclusions.

Figure 5:  Real and imaginary parts of $\varepsilon^{MG}_{\rm t}$ and
$\varepsilon^{MG}_{\rm z}$ plotted against $\alpha\in(0,4]$ and
$\kappa\in(0,1)$ for double-truncated spheroidal inclusions.

Figure 6:  Real and imaginary parts of $\varepsilon^{MG}_{\rm x}$,
$\varepsilon^{MG}_{\rm y}$, and $\varepsilon^{MG}_{\rm z}$ plotted against
$\alpha\in(0,3]$ and $\beta\in(0,3]$ for hemi-ellipsoidal inclusions.

Figure 7:  Real and imaginary parts of $\varepsilon^{MG}_{\rm x}$,
$\varepsilon^{MG}_{\rm y}$, and $\varepsilon^{MG}_{\rm z}$ plotted against
$\alpha\in(0,3]$ and $\beta\in(0,3]$ for double-truncated ellipsoidal
inclusions with $\kappa=0.3$.

Supplementary Material:
Anisotropic homogenized composite mediums arising from truncated spheres,
spheroids, and ellipsoids

Tom G. Mackay†††E–mail<EMAIL_ADDRESS>
School of Mathematics and Maxwell Institute for Mathematical Sciences
University of Edinburgh, Edinburgh EH9 3FD, UK
and
NanoMM — Nanoengineered Metamaterials Group
Department of Engineering Science and Mechanics
The Pennsylvania State University, University Park, PA 16802–6812, USA
Akhlesh Lakhtakia
NanoMM — Nanoengineered Metamaterials Group
Department of Engineering Science and Mechanics
The Pennsylvania State University, University Park, PA 16802–6812, USA

##  1 Introduction

In this document, the derivations of the polarizability density dyadic and the
Maxwell Garnett formalism are outlined, no novelty being claimed for their
derivations. The electrodynamic provenance of these theoretical constructions
is emphasized.

As regards notational matters: vectors are boldface; dyadics [1, 2] are double
underlined; $i=\sqrt{-1}$; the angular frequency is $\omega$ and an
$\exp(-i\omega t)$ dependence on time $t$ is implicit;
$\varepsilon_{\scriptscriptstyle 0}$ is the free-space permittivity,
$\mu_{\scriptscriptstyle 0}$ is the free-space permeability, and
$k_{\scriptscriptstyle 0}=\omega\sqrt{\varepsilon_{\scriptscriptstyle
0}\mu_{\scriptscriptstyle 0}}$ is the free-space wavenumber;
$\underline{\underline{I}}=\hat{{\bf u}}_{\rm x}\hat{{\bf u}}_{\rm
x}+\hat{{\bf u}}_{\rm y}\hat{{\bf u}}_{\rm y}+\hat{{\bf u}}_{\rm z}\hat{{\bf
u}}_{\rm z}$ is the identity dyadic; ${\bf r}=x\hat{{\bf u}}_{\rm
x}+y\hat{{\bf u}}_{\rm y}+z\hat{{\bf u}}_{\rm z}$ is the position vector; and
$\hat{{\bf u}}_{\rm x}$, $\hat{{\bf u}}_{\rm y}$, and $\hat{{\bf u}}_{\rm z}$
are the unit vectors aligned with the coordinate axes of the Cartesian
coordinate system $(x,y,z)$.

##  2 Polarizability density dyadic

A crucial ingredient in homogenization formalisms for particulate composite
mediums is the polarizability density dyadic of a solitary inclusion immersed
in a host medium. Sometime in 1848, Mossotti [3] conceptualized the existence
of the polarizability density scalar of an electrically small, isotropic-
dielectric sphere in free space. In this section, we present a clear
derivation of the polarizability density dyadic of an electrically small,
isotropic-dielectric inclusion of quite general shape embedded in an isotropic
dielectric host medium. No novelty is implied for this derivation.

Let the closed surface $S_{\rm in}$ separate the electrically small region
$V_{\rm in}$ from the region $V_{\rm h}$ that extends to infinity in all
directions. The shape of $V_{\rm in}$ must be such that the unit outward
normal $\hat{{\bf u}}_{\rm n}({\bf r})$ can be unambiguously identified at
every point ${\bf r}\in S_{\rm in}$; if necessary, wedges and vertexes on
$S_{\rm in}$ can be rounded off slightly to overcome this restriction.

Whereas $V_{\rm in}$ is occupied by a homogeneous dielectric medium of
relative permittivity $\varepsilon_{\rm in}(\omega)$, $V_{\rm h}$ is occupied
by a homogeneous dielectric medium of relative permittivity $\varepsilon_{\rm
h}(\omega)=n_{\rm h}^{2}(\omega)$.

Suppose that the region $V_{\rm in}$ is irradiated by a monochromatic
electromagnetic field phasor ${\bf E}_{\rm s}({\bf r},\omega)$ whose source is
confined to the bounded region $V_{\rm s}\subset V_{\rm h}$ far from $V_{\rm
in}$. Then, the actual electric field phasor at any location ${\bf r}\notin
V_{\rm s}$ is given exactly by [4, 5, 6]

${\bf E}({\bf r},\omega)={\bf E}_{\rm s}({\bf r},\omega)+k_{\scriptscriptstyle
0}^{2}\left[\varepsilon_{\rm in}(\omega)-\varepsilon_{\rm
h}(\omega)\right]\iiint_{V_{\rm in}}\left[\underline{\underline{G}}_{\rm
h}({\bf r},{\bf r}^{\prime};\omega)\mbox{ \tiny{${}^{\bullet}$} }{\bf E}({\bf
r}^{\prime},\omega)\right]\,d^{3}{\bf r}^{\prime}\,,\quad{\bf r}\notin V_{\rm
s}\,,$ (1)

where

$\underline{\underline{G}}_{\rm h}({\bf r},{\bf
r}^{\prime};\omega)=\left[\underline{\underline{I}}+\frac{\nabla\nabla}{k_{\scriptscriptstyle
0}^{2}\varepsilon_{\rm
h}(\omega)}\right]\frac{\exp\left[ik_{\scriptscriptstyle 0}n_{\rm
h}(\omega)|{\bf r}-{\bf r}^{\prime}|\right]}{4\pi|{\bf r}-{\bf
r}^{\prime}|}\,.$ (2)

Setting ${\bf r}={\bf r}_{\rm o}$ in Eq. (1), where ${\bf r}_{\rm o}\in V_{\rm
in}$ but ${\bf r}_{\rm o}\notin S_{\rm in}$, we get

${\bf E}({\bf r}_{\rm o},\omega)={\bf E}_{\rm s}({\bf r}_{\rm
o},\omega)+k_{\scriptscriptstyle 0}^{2}\left[\varepsilon_{\rm
in}(\omega)-\varepsilon_{\rm h}(\omega)\right]\iiint_{V_{\rm
in}}\left[\underline{\underline{G}}_{\rm h}({\bf r}_{\rm o},{\bf
r}^{\prime};\omega)\mbox{ \tiny{${}^{\bullet}$} }{\bf E}({\bf
r}^{\prime},\omega)\right]\,d^{3}{\bf r}^{\prime}\,,\quad{\bf r}_{\rm o}\in
V_{\rm in}\,.$ (3)

The integral on the right side of Eq. (3) is improper due to the appearance of
the term $|{\bf r}_{\rm o}-{\bf r}^{\prime}|^{-3}$ in the integrand [7]. In
order to overcome this difficulty in a “unified and consistent” manner [8],
the term $\underline{\underline{G}}_{\rm P}({\bf r}_{\rm o},{\bf
r}^{\prime};\omega)\mbox{ \tiny{${}^{\bullet}$} }{\bf E}({\bf r}_{\rm
o},\omega)$ is added to and subtracted from that integrand to obtain [9]

$\displaystyle{\bf E}({\bf r}_{\rm o},\omega)={\bf E}_{\rm s}({\bf r}_{\rm
o},\omega)+k_{\scriptscriptstyle 0}^{2}\left[\varepsilon_{\rm
in}(\omega)-\varepsilon_{\rm h}(\omega)\right]\iiint_{V_{\rm
in}}\left[\underline{\underline{G}}_{\rm h}({\bf r}_{\rm o},{\bf
r}^{\prime};\omega)\mbox{ \tiny{${}^{\bullet}$} }{\bf E}({\bf{\bf
r}^{\prime}},\omega)-\underline{\underline{G}}_{\rm P}({\bf r}_{\rm o},{\bf
r}^{\prime};\omega)\mbox{ \tiny{${}^{\bullet}$} }{\bf E}({\bf r}_{\rm
o},\omega)\right]\,d^{3}{\bf r}^{\prime}$
$\displaystyle\qquad+k_{\scriptscriptstyle 0}^{2}\left[\varepsilon_{\rm
in}(\omega)-\varepsilon_{\rm h}(\omega)\right]\left[\iiint_{V_{\rm
in}}\underline{\underline{G}}_{\rm P}({\bf r}_{\rm o},{\bf
r}^{\prime};\omega)\,d^{3}{\bf r}^{\prime}\right]\mbox{ \tiny{${}^{\bullet}$}
}{\bf E}({\bf r}_{\rm o},\omega)\,,\quad{\bf r}_{\rm o}\in V_{\rm in}\,,$ (4)

where

$\underline{\underline{G}}_{\rm P}({\bf r},{\bf
r}^{\prime};\omega)=\frac{1}{k_{\scriptscriptstyle 0}^{2}\varepsilon_{\rm
h}(\omega)}\nabla\nabla\frac{1}{4\pi|{\bf r}-{\bf r}^{\prime}|}\,.$ (5)

In order to treat the second integral on the right side of Eq. (4), first the
identity $\nabla\nabla\left(|{\bf r}-{\bf
r}^{\prime}|^{-1}\right)=-\nabla\nabla^{\prime}\left(|{\bf r}-{\bf
r}^{\prime}|^{-1}\right)$ is invoked and then the Gauss theorem is used to
obtain [8, 4]

${\bf E}({\bf r}_{\rm o},\omega)={\bf E}_{\rm s}({\bf r}_{\rm
o},\omega)+k_{\scriptscriptstyle 0}^{2}\left[\varepsilon_{\rm
in}(\omega)-\varepsilon_{\rm h}(\omega)\right]\left[{\bf M}({\bf r}_{\rm
o},\omega)-\frac{1}{k_{\scriptscriptstyle 0}^{2}\varepsilon_{\rm
h}(\omega)}\underline{\underline{L}}({\bf r}_{\rm o})\mbox{
\tiny{${}^{\bullet}$} }{\bf E}({\bf r}_{\rm o},\omega)\right]\,,\quad{\bf
r}_{\rm o}\in V_{\rm in}\,,$ (6)

where

${\bf M}({\bf r}_{\rm o},\omega)=\iiint_{V_{\rm
in}}\left[\underline{\underline{G}}_{\rm h}({\bf r}_{\rm o},{\bf
r}^{\prime};\omega)\mbox{ \tiny{${}^{\bullet}$} }{\bf E}({\bf{\bf
r}^{\prime}},\omega)-\underline{\underline{G}}_{\rm P}({\bf r}_{\rm o},{\bf
r}^{\prime};\omega)\mbox{ \tiny{${}^{\bullet}$} }{\bf E}({\bf r}_{\rm
o},\omega)\right]\,d^{3}{\bf r}^{\prime}\,,$ (7)

and

$\underline{\underline{L}}({\bf r}_{\rm o})=\frac{1}{4\pi}\iint_{S_{\rm
in}}\hat{{\bf u}}_{\rm n}({\bf r}^{\prime})\frac{{\bf r}^{\prime}-{\bf r}_{\rm
o}}{|{\bf r}^{\prime}-{\bf r}_{\rm o}|^{3}}\,d^{2}{\bf r}^{\prime}\,.$ (8)

Note that $\underline{\underline{L}}({\bf r}_{\rm o})$ is a symmetric dyadic
whose trace equals unity [10].

In the limit ${\bf r}^{\prime}\to{\bf r}_{\rm o}$, the singularity of
$\underline{\underline{G}}_{\rm h}({\bf r}_{\rm o},{\bf r}^{\prime};\omega)$
is offset by the singularity of $\underline{\underline{G}}_{\rm P}({\bf
r}_{\rm o},{\bf r}^{\prime};\omega)$ in the integrand of the volume integral
${\bf M}({\bf r}_{\rm o},\omega)$. Furthermore, this integral is not improper
if it is assumed that the electric field satisfies the Hölder continuity
condition in $V_{\rm in}$ [9, 8]: real constants $\bar{\gamma}>0$ and
$\bar{\nu}>0$ exist such that $|\hat{{\bf u}}_{j}\mbox{ \tiny{${}^{\bullet}$}
}\left[{\bf E}({\bf r},\omega)-{\bf E}({\bf r}_{\rm
o},\omega)\right]|\leq\bar{\gamma}|{\bf r}-{\bf r}_{\rm o}|^{\bar{\nu}}$ for
all ${\bf r}\in V_{\rm in}$ and $j\in\left\\{x,y,z\right\\}$.

Provided that $V_{\rm in}$ is sufficiently small in relation to the
wavelengths in both mediums, we can set ${\bf E}({\bf
r}^{\prime},\omega)\simeq{\bf E}({\bf r}_{\rm o},\omega)$ for all ${\bf
r}^{\prime}\in V_{\rm in}$, as is commonly done [8, 9, 4] following Maxwell
[11] and has been computationally verified [8]. Then, ${\bf M}({\bf r}_{\rm
o},\omega)$ can be factored as $\underline{\underline{M}}({\bf r}_{\rm
o},\omega)\mbox{ \tiny{${}^{\bullet}$} }{\bf E}({\bf r}_{\rm o},\omega)$ and
Eq. (6) simplifies to

${\bf E}({\bf r}_{\rm o},\omega)\simeq{\bf E}_{\rm s}({\bf r}_{\rm
o},\omega)+k_{\scriptscriptstyle 0}^{2}\left[\varepsilon_{\rm
in}(\omega)-\varepsilon_{\rm
h}(\omega)\right]\left[\underline{\underline{M}}({\bf r}_{\rm
o},\omega)-\frac{1}{k_{\scriptscriptstyle 0}^{2}\varepsilon_{\rm
h}(\omega)}\underline{\underline{L}}({\bf r}_{\rm o})\right]\mbox{
\tiny{${}^{\bullet}$} }{\bf E}({\bf r}_{\rm o},\omega)\,,\quad{\bf r}_{\rm
o}\in V_{\rm in}\,.$ (9)

The concept of $\underline{\underline{M}}({\bf r}_{\rm o},\omega)$ must be
credited to Lorenz [12, 13].

Even more commonly, with good justification for really small $V_{\rm in}$ [4,
14], $\underline{\underline{M}}({\bf r}_{\rm o},\omega)$ is ignored in favor
of $\underline{\underline{L}}({\bf r}_{\rm o})/k_{\scriptscriptstyle
0}^{2}\varepsilon_{\rm h}$ [4, 10, 15] to get

${\bf E}({\bf r}_{\rm
o},\omega)\approx\left[\underline{\underline{I}}+\frac{\varepsilon_{\rm
in}(\omega)-\varepsilon_{\rm h}(\omega)}{\varepsilon_{\rm
h}(\omega)}\underline{\underline{L}}({\bf r}_{\rm o})\right]^{-1}\mbox{
\tiny{${}^{\bullet}$} }{\bf E}_{s}({\bf r}_{\rm o},\omega)\,,\quad{\bf r}_{\rm
o}\in V_{\rm in}\,.$ (10)

The electrically small region $V_{\rm in}$ can be said to have acquired an
electric dipole moment [4]

${\bf p}(\omega)={v_{\rm in}}\,\varepsilon_{\scriptscriptstyle
0}\left[\varepsilon_{\rm in}(\omega)-\varepsilon_{\rm h}(\omega)\right]{\bf
E}({\bf r}_{\rm o},\omega)\,$ (11)

located at ${\bf r}_{\rm o}$ as the source of the field scattered into $V_{\rm
h}$, with $v_{\rm in}$ denoting the volume of $V_{\rm in}$. Substituting Eq.
(10) in Eq. (11), we get

${\bf p}(\omega)={v_{\rm in}}\,\underline{\underline{a}}(\omega)\mbox{
\tiny{${}^{\bullet}$} }{\bf E}_{\rm s}({\bf r}_{\rm o},\omega)\,,$ (12)

where the polarizability density dyadic

$\underline{\underline{a}}(\omega)\approx\varepsilon_{\scriptscriptstyle
0}\left[\varepsilon_{\rm in}(\omega)-\varepsilon_{\rm
h}(\omega)\right]\left[\underline{\underline{I}}+\frac{\varepsilon_{\rm
in}(\omega)-\varepsilon_{\rm h}(\omega)}{\varepsilon_{\rm
h}(\omega)}\underline{\underline{L}}({\bf r}_{\rm o})\right]^{-1}\,.$ (13)

The right side of Eq. (13) is definitely an approximation of
$\underline{\underline{a}}(\omega)$, even for spheres, for two reasons. The
first reason is the removal of $\underline{\underline{M}}({\bf r}_{\rm
o},\omega)$ even though $V_{\rm in}$ cannot be infinitesimally small. The
second reason is that the Hölder continuity condition can best be applied
inside the largest ellipsoid [11] centered at ${\bf r}_{\rm o}$ and inscribed
in $V_{\rm in}$. The electric field phasor in the portion of $V_{\rm in}$
outside that inscribed ellipsoid having only been approximately accounted for,
${\bf r}_{\rm o}$ must be carefully chosen [16]. Regardless, Eq. (13) is
implicitly used for cubes [17, 18], rectangular parallelopipeds [19], and
circular cylinders of finite height [10]. Furthermore, the analog of Eq. (13)
exists for quasi-two-dimensional problems [10].

##  3 Maxwell Garnett homogenization formalism

Let us now briefly go through the Maxwell Garnett formalism in this section,
no novelty being implied for the derivation. Suppose that all space $V_{\rm
all}$ be occupied by a composite medium comprising a random distribution of
identical, similarly oriented, electromagnetically small inclusions surrounded
by a host medium. The inclusion medium has relative permittivity scalar
${\varepsilon_{\rm in}}(\omega)$, and the host medium ${\varepsilon_{\rm
h}}(\omega)$. The volume fraction of the composite material occupied by the
inclusions is denoted by $f_{\rm in}$. With $N_{\rm in}$ denoting the number
density of inclusions, $f_{\rm in}/N_{\rm in}$ is the volume of each
inclusion.

The electric field in $V_{\rm all}$ satisfies

$\nabla\times\left[\nabla\times{\bf E}({\bf
r},\omega)\right]-k_{\scriptscriptstyle 0}^{2}\varepsilon_{\rm h}(\omega){\bf
E}({\bf r},\omega)=i\omega\mu_{\scriptscriptstyle 0}{\bf J}({\bf
r},\omega)\,,\quad{\bf r}\in V_{\rm all}\,,$ (14)

where the electric current density phasor ${\bf J}({\bf r},\omega)={\bf 0}$ if
${\bf r}$ lies in a region occupied by the host medium but ${\bf J}({\bf
r},\omega)\neq{\bf 0}$ if ${\bf r}$ lies inside an inclusion [4, 20]. Clearly,
Eq. (14) represents a grainy visualization of the electromagnetic field in the
composite medium. If this medium were to be homogenized, then the HCM will
have

${\bf D}({\bf r},\omega)=\varepsilon_{\scriptscriptstyle 0}\varepsilon_{\rm
h}(\omega){\bf E}({\bf r},\omega)+{\bf P}_{\rm xs}({\bf
r},\omega)=\varepsilon_{\scriptscriptstyle
0}\underline{\underline{\varepsilon}}^{\rm eff}(\omega)\mbox{
\tiny{${}^{\bullet}$} }{\bf E}({\bf r},\omega)\,,$ (15)

as its constitutive relation, with ${\bf D}({\bf r},\omega)$ as the electric
displacement phasor, ${\bf P}_{\rm xs}({\bf r},\omega)$ as the excess
polarization representing the homogenized distribution of the inclusion medium
inside the HCM, and $\underline{\underline{\varepsilon}}^{\rm eff}(\omega)$ as
the effective relative permittivity dyadic of the HCM.

In the Maxwell Garnett homogenization formalism, Eq. (14) is replaced by

$\nabla\times\left[\nabla\times{\bf E}({\bf
r},\omega)\right]-k_{\scriptscriptstyle 0}^{2}\varepsilon_{\rm h}(\omega){\bf
E}({\bf r},\omega)=i\omega\mu_{\scriptscriptstyle 0}{\bf J}_{\rm avg}({\bf
r},\omega)\,,\quad{\bf r}\in V_{\rm all}\,,$ (16)

where ${\bf J}_{\rm avg}({\bf r},\omega)=-i\omega{\bf P}_{\rm xs}({\bf
r},\omega)$ is the spatial average of ${\bf J}({\bf r},\omega)$ in some small
neighborhood of ${\bf r}$ called the Lorentzian cavity [21]. Accordingly, the
electric field phasor in the HCM is given by

${\bf E}({\bf r},\omega)={\bf E}_{\rm cf}({\bf
r},\omega)+\omega^{2}\mu_{\scriptscriptstyle 0}\iiint_{V_{\rm
all}}\underline{\underline{G}}_{\rm h}({\bf r},{\bf r}^{\prime};\omega)\mbox{
\tiny{${}^{\bullet}$} }{\bf P}_{\rm xs}({\bf r}^{\prime},\omega)\,d^{3}{\bf
r}^{\prime}\,,\quad{\bf r}\in V_{\rm all}\,,$ (17)

where ${\bf E}_{\rm cf}({\bf r},\omega)$ is the complementary function.

Let the electrically small spherical region $V_{\ell}$ centered at ${\bf r}$
be identified as the Lorentzian cavity. Then, we can define

${\bf E}_{\rm loc}({\bf r},\omega)={\bf E}_{\rm cf}({\bf
r},\omega)+\omega^{2}\mu_{\scriptscriptstyle 0}\iiint_{V_{\rm
all}-V_{\ell}}\underline{\underline{G}}_{\rm h}({\bf r},{\bf
r}^{\prime};\omega)\mbox{ \tiny{${}^{\bullet}$} }{\bf P}_{xs}({\bf
r}^{\prime},\omega)\,d^{3}{\bf r}^{\prime}$ (18)

as the electric field phasor present at ${\bf r}$ if the excess polarization
were to be null-valued in $V_{\ell}$. Equations (17) and (18) together yield

${\bf E}({\bf r},\omega)={\bf E}_{\rm loc}({\bf
r},\omega)+\omega^{2}\mu_{\scriptscriptstyle
0}\iiint_{V_{\ell}}\underline{\underline{G}}_{\rm h}({\bf r},{\bf
r}^{\prime};\omega)\mbox{ \tiny{${}^{\bullet}$} }{\bf P}_{\rm xs}({\bf
r}^{\prime},\omega)\,d^{3}{\bf r}^{\prime}\,.$ (19)

The conventional assumption is that the excess polarization is uniform inside
a Lorentzian cavity, so that [21]

${\bf E}({\bf r},\omega)\simeq{\bf E}_{\rm loc}({\bf
r},\omega)-\frac{1}{3\varepsilon_{\scriptscriptstyle 0}\varepsilon_{\rm
h}(\omega)}{\bf P}_{\rm xs}({\bf r},\omega)\,$ (20)

emerges from Eq. (19).

In the Maxwell Garnett formalism, we set ${\bf P}_{\rm xs}=N_{\rm in}{\bf p}$
and use Eq. (12) with ${\bf E}_{\rm loc}({\bf r},\omega)$ fulfilling the role
of ${\bf E}_{\rm s}({\bf r},\omega)$; accordingly,

${\bf P}_{\rm xs}({\bf r},\omega)=f_{\rm
in}\,\underline{\underline{a}}(\omega)\mbox{ \tiny{${}^{\bullet}$} }{\bf
E}_{\rm loc}({\bf r},\omega)\,.$ (21)

We then get

${\bf E}_{\rm loc}({\bf
r},\omega)\simeq\left[\underline{\underline{I}}-\frac{1}{3\varepsilon_{\scriptscriptstyle
0}\varepsilon_{\rm h}(\omega)}f_{\rm
in}\,\underline{\underline{a}}(\omega)\right]^{-1}\mbox{ \tiny{${}^{\bullet}$}
}{\bf E}({\bf r},\omega)$ (22)

from Eq. (20), followed by

${\bf P}_{\rm xs}({\bf r},\omega)\simeq f_{\rm
in}\,\underline{\underline{a}}(\omega)\mbox{ \tiny{${}^{\bullet}$}
}\left[\underline{\underline{I}}-\frac{1}{3\varepsilon_{\scriptscriptstyle
0}\varepsilon_{\rm h}(\omega)}f_{\rm
in}\,\underline{\underline{a}}(\omega)\right]^{-1}\mbox{ \tiny{${}^{\bullet}$}
}{\bf E}({\bf r},\omega)$ (23)

from Eq. (21). The Maxwell Garnett estimate

$\underline{\underline{\varepsilon}}^{MG}(\omega)={\varepsilon_{\rm
h}}(\omega)\,\underline{\underline{I}}+{\frac{f_{\rm
in}}{\varepsilon_{\scriptscriptstyle
0}}}\,\underline{\underline{a}}(\omega)\mbox{ \tiny{${}^{\bullet}$}
}\left[\underline{\underline{I}}-\frac{1}{3\varepsilon_{\scriptscriptstyle
0}\varepsilon_{\rm h}(\omega)}f_{\rm
in}\,\underline{\underline{a}}(\omega)\right]^{-1}$ (24)

of the relative permittivity dyadic
$\underline{\underline{\varepsilon}}^{eff}(\omega)$ then emerges from Eqs.
(15) and (23).

Note that the Maxwell Garnett formalism does not really depend on any explicit
formula for $\underline{\underline{a}}(\omega)$, although more than a century
ago, Maxwell Garnett himself [22] and others [23] had used the
Mossotti–Clausius expression [13] for the polarizability scalar of an
electrically small, isotropic, dielectric sphere. Parenthetically, even for
that simple inclusion, other expressions have also been used [13, 12, 24, 25].
Finally, Eq. (24) for $\underline{\underline{\varepsilon}}^{MG}$ readily
generalizes to cases involving inclusions of more than one type as well as to
cases involving anisotropic inclusion materials [1, 26].

## References

  * [1] T. G. Mackay and A. Lakhtakia, _Modern Analytical Electromagnetic Homogenization with Mathematica_ , _2nd Ed._ (IOP Publishing, 2020).
  * [2] M. Faryad and A. Lakhtakia, _Infinite-Space Dyadic Green Functions in Electromagnetism_. San Rafael, CA, USA: Morgan & Claypool, 2018.
  * [3] O. F. Mossotti, “Discussione analitica sull’influenza che l’azione di un mezzo dielettrico ha sulla distribuzione dell’elettricità alla superficie di più corpi electrici disseminati in esso,” Mem. Math. Fisica Modena 24, 49–74 (1850).
  * [4] A. Lakhtakia, “On two numerical techniques for light scattering by dielectric agglomerated structures,” J. Res. Nat. Inst. Stand. Technol. 98, 699–716 (1993).
  * [5] W. C. Chew, Waves and Fields in Inhomogeneous Media (IEEE, 1995); Sec. 8.9.2.
  * [6] M. Sancer, K. Sertel, J. Volakis, and P. Alstine, “On volume integral equations,” IEEE Trans. Antennas Propagat. 54, 1488–1495 (2006).
  * [7] O. D. Kellogg, Foundations of Potential Theory (Springer, 1967); Sec. VI.2.
  * [8] J. J. H. Wang, “A unified and consistent view on the singularities of the electric dyadic Green’s function in the source region,” IEEE Trans. Antennas Propagat. 30, 463–468 (1982).
  * [9] J. G. Fikioris, “Electromagnetic field inside a current-carrying region,” J. Math. Phys. 6, 1617–1620 (1965).
  * [10] A. D. Yaghjian, “Electric dyadic Green’s function in the source region,” Proc. IEEE 68, 248–263 (1980).
  * [11] J. C. Maxwell, A Treatise on Electricity and Magnetism. _3rd Ed._ (Clarendon Press,1904); Arts. 437 and 438.
  * [12] L. V. Lorenz, “Experimentale og theoretiske Undersøgelser over Legemernes Brydningsforhold, II,” K. Dan. Vidensk. Selsk. Forh. 10, 485–518 (1875).
  * [13] M. T. Prinkey, A. Lakhtakia, and B. Shanker, “On the extended Maxwell-Garnett and the extended Bruggeman approaches for dielectric-in-dielectric composites,” Optik 96, 25–30 (1994).
  * [14] R. C. Robertson, “Computation of the electric field within a region of uniform current density,” Int. J. Electron. 69, 707–716 (1990).
  * [15] J. Van Bladel, “Some remarks on Green’s dyadic for infinite space,” IRE Trans. Antennas Propagat. 9, 563–566 (1961).
  * [16] T. G. Mackay and A. Lakhtakia, “Depolarization dyadics for truncated spheres, spheroids, and ellipsoids,” https://doi.org/10.48550/arXiv.2312.08674
  * [17] D. E. Livesay and K.-M. Chen, “Electromagnetic field induced inside arbitrarily shaped biological bodies,” IEEE Trans. Microwave Theory Tech. 22, 1273–1280 (1974).
  * [18] J. Avelin, A. N. Arslan, J. Brännback, M. Flykt, C. Icheln, J. Juntunen, K. Kärkkäinen, T. Niemi, O. Nieminen, T. Tares, C. Toma, T. Uusitupa, and A. Sihvola, “Electric fields in the source region: the depolarization dyadic for a cubic cavity,” Elec. Eng. 81, 199–202 (1998).
  * [19] M. Cauterman, P. Degauque, B. Demoulin, and R. Gabillard, “Remote sensing of buried resistive inhomogeneities by electromagnetic transmission measurements between the ground surface and a borehole—Theory and experiment,” Radio Sci. 13, 371–378 (1978).
  * [20] R. F. Harrington, Time-Harmonic Electromagnetic Fields (IEEE, 2001); Sec. 3-11.
  * [21] D. E. Aspnes, “Local-field effects and effective-medium theory: a microscopic perspective,” Am. J. Phys. 50, 704–709 (1982).
  * [22] J. C. Maxwell Garnett, “Colours in metal glasses and thin films,” Philos. Trans. R. Soc. Lond. A 203, 385–420 (1904).
  * [23] A. Sihvola, Electromagnetic Mixing Formulas and Applications (IEE, 1999).
  * [24] C. E. Dungey and C. F. Bohren, “Light scattering by nonspherical particles: a refinement to the coupled-dipole method,” J. Opt. Soc. Am. A 8, 81–87 (1991).
  * [25] P. Mallet, C. A. Guérin, and A. Sentenac, “Maxwell-Garnett mixing rule in the presence of multiple scattering: Derivation and accuracy,” Phys. Rev. B 72, 014205 (2005).
  * [26] A. Lakhtakia, “Frequency-dependent continuum electromagnetic properties of a gas of scattering centers,” Adv. Chem. Phys. 85(2), 311–359 (1993).


# Building astroBERT, a language model for Astronomy & Astrophysics

Felix Grezes1, Sergi Blanco-Cuaresma, Alberto Accomazzi, Michael J. Kurtz,
Golnaz Shapurian, Edwin Henneken, Carolyn S. Grant, Donna M. Thompson, Roman
Chyla, Stephen McDonald, Timothy W. Hostetler, Matthew R. Templeton, Kelly E.
Lockhart, Nemanja Martinovic, Shinyi Chen, Chris Tanner, and Pavlos Protopapas

###### Abstract

The existing search tools for exploring the NASA Astrophysics Data System
(ADS) can be quite rich and empowering (e.g., similar and trending operators),
but researchers are not yet allowed to fully leverage semantic search. For
example, a query for "results from the Planck mission" should be able to
distinguish between all the various meanings of Planck (person, mission,
constant, institutions and more) without further clarification from the user.
At ADS, we are applying modern machine learning and natural language
processing techniques to our dataset of recent astronomy publications to train
astroBERT, a deeply contextual language model based on research at Google.
Using astroBERT, we aim to enrich the ADS dataset and improve its
discoverability, and in particular we are developing our own named entity
recognition tool. We present here our preliminary results and lessons learned.

1Harvard-Smithsonian Center for Astrophysics, Cambridge, MA, USA;
<EMAIL_ADDRESS>
##  1 Introduction

The NASA ADS (https://ui.adsabs.harvard.edu) has been publicly servicing the
astrophysical research community since Kurtz et al. (1993). It maintains a
searchable bibliographic collection of more than 15 million records. The
database can be explored via complex queries using filters (e.g. authors, year
range, text in title and/or body, refereed or not) combined with Boolean
operators to explore the citations and references, and more advanced operators
such as paper similarity, co-reads, and reviews.

As powerful as the existing tools are, they do not make use of recent advances
in natural language processing and semantic search. Ideally, if a user wanted
to filter by papers that use data from the ESA Planck space observatory, a
simple query of "Planck mission" should suffice, without requiring manual
input of specific keywords such as year:2009-2022 or NOT author:"Planck, Max".
Such queries highlight the inherent ambiguity of natural language and how
remarkably flexible humans are at understanding semantics from contextual
clues. Currently, ADS handles many of the most common ambiguities with
carefully crafted rules (e.g. using SIMBAD objects).

To automate the process of identifying, disambiguating, and tagging named
entities within the ADS database, we are building astroBERT, a language model
for astronomical content. This astroBERT is based on Google’s BERT deep neural
network transformer architecture (Devlin et al. 2018), inspired by the
AllenAI’s SciBERT (Beltagy et al. 2019) efforts, built using tools from
Huggingface (Wolf et al. 2019), and is trained is trained using a large
collection of recent astronomy ADS papers.

Our preliminary results show that astroBERT outperforms BERT on the named
entity recognition task on ADS data. Based on this encouraging result, we can
now work on improving and integrating astroBERT into ADS services, and
providing access to it for the astrophysics community.

##  2 Related Work: BERT and SciBERT

Devlin et al. (2018) introduced the Bidirectional Encoder Representations from
Transformers (BERT) language model, based on the Transformer technology
(Vaswani et al. 2017). It was trained with text from Wikipedia and the Brown
Corpus, and then fine-tuned to achieve state-of-the-art performance on 11 NLP
tasks, including sentiment analysis, semantic role labeling, sentence
classification and word disambiguation.

Building off the success of BERT, Beltagy et al. (2019) created SciBERT,
Focusing on scientific papers from the Semantic Scholar corpus, SciBERT is
part of a growing list of models that adapt BERT to specific domains and
tasks. Its state-of-the-art performance on a wide range of scientific domain
NLP tasks motivates our work on astroBERT. The SciBERT paper shows that pre-
training BERT with domain-specific language data improves its performance when
compared to the original BERT.

##  3 Technical Details

The work presented here was implemented using the open-source python-based
Huggingface (Wolf et al. 2019) library, which provides easy access to state-
of-the-art NLP models, and tools. Runtimes are reported as run on a machine
with two NVIDIA V100 GPUs, two 12-core Intel Xeon Gold 6246 CPUs, and 768GB of
RAM. Unless specified, all parameters were left to Huggingface v4.10.2
defaults .

###  3.1 Data Preparation

In order to train our language model, our first task was to build a large
dataset of clean English language text from astronomy papers. From the ADS
database, we selected papers for which the source was provided by the
publishers in XML format. The final dataset is made up of 395,499 documents
(16GB on disk). Per best practices, this was evenly split into training,
evaluation and testing datasets. Despite our efforts to select the best
quality data sources, we required to execute a clean-up to remove unprintable
unicode characters and rare unicode characters (those appearing 50 times or
less). We also replaced accented characters by the non-accented counterparts
when possible. In future re-iterations we will consider a more exhaustive
cleaning to filter out meaningless strings representing data tables,
mathematical formulas, and appendices of raw data.

###  3.2 Tokenizer

An important part of the text processing pipeline is the tokenizer, which
converts text data into inputs compatible with the model To stick as close as
possible to the original BERT design choices. we trained a WordPiece tokenizer
maximum vocabulary size of 30,000 tokens. Due to the extensive use of
accronyms in the astrophysics literature, and its tendency to correspond to
English words (e.g., AURA, BLAST, CLUSTER, INTEGRAL, WISE), we decided to keep
the case information in our model. the tokenizer did not automatically convert
every character to the same case, however accents were removed to simplify
downstream tasks. Training the tokenizer over the entire data took
approximately 45 minutes. The astroBERT vocabulary overlaps with that of BERT
and SciBERT by 24.5% and 35.3% respectively (over total number of tokens),
while BERT and SciBERT overlapped by 27.0%. The total size of the data
corresponds to 3,819,322,591 tokens, or 2,977,635,680 words and 121,207,934
sentences when parsed by NLTK (Bird et al. 2009) .

###  3.3 The MLM, NSP and NER Tasks

Following the methodology of the BERT paper, astroBERT is pre-trained on two
tasks that do not require hand labeled data: first the masked language model
(MLM) task in which the model predicts tokens that were randomly masked by
using the contextual text; and second MLM jointly with the next sentence
prediction (NSP) task in which the model predicts if a sentence is subsequent
to another in the original document. The original random masking and shuffling
probabilities were kept unchanged (imitating the original BERT setup).

Following astroBERT pre-training, we transferred the models to our first
downstream task: named entity recognition (NER), in which the model tries to
identify organizations. Using a list of 908 astronomy organization acronyms
and their full names, we scanned the acknowledgement sections of 44K astronomy
manuscripts and built a dataset of 1856 sentences with 6279 annotated
organizations. The models were evaluated using 10-fold cross-validation.

##  4 Training astroBERT and Preliminary Results

We trained astroBERT(MLM) for 15 epochs on the MLM task, and astroBERT(NSP
+MLM) for 3 epochs on the joint NSP and MLM tasks. astroBERT(MLM) was
initialized with BERT weights, and astroBERT(NSP+MLM) with the trained
astroBERT (MLM). Training astroBERT takes approximately 8 hours per epoch for
the MLM task, and 22 hours for the joint NSP+MLM task. Training was done using
mixed floating point precision, which has been shown to reduce memory
consumption and lead to faster learning (Micikevicius et al. 2017). The NER
models were fine-tuned for 3 epochs, each epoch taking approximately 90
seconds.

In table 1 we report perplexity scores our astroBERT, BERT, and SciBERT models
on our cleaned ADS dataset, as well across various standard English texts
dataset provided by Huggingface. Because these models use different
vocabularies, scores should be compared within model and across datasets, and
not across models. These scores show that BERT, which is trained on Wikipedia
English, can be re-trained into astroBERT to performs best on ADS text.
SciBERT, which was trained on data from Semantic Scholar (not tested here)
performs similarly on all 3 scientific text datasets. As expected, adding the
NSP task to astroBERT(NSP+MLM) training causes it to perform worse in terms of
perplexity on the MLM task.

However as shown in table 2, in which we report the average f1-scores and
standard deviations across cross-validation runs for BERT and astroBERT after
fine-tuning, astroBERT(NSP+MLM) slightly outperforms astroBERT(MLM) on the NER
task, and both outperform the original BERT model.

Table 1.: Comparisons of perplexity scores across BERT models | BERT | SciBERT | astroBERT | astroBERT
---|---|---|---|---
|  |  | (MLM) | (NSP+MLM)
ADS | 38.4 | 4.05 | 4.16 | 5.71
wikitext-2-raw-v1 | 29.7 | 8.64 | 18.5 | 25.9
ptb_text_only | 82.8 | 9.63 | 20.0 | 27.4
scientific_papers (arxiv) | 41.3 | 4.23 | 4.82 | 6.87
scientific_papers (pubmed) | 34.0 | 3.96 | 8.41 | 12.6
Table 2.: Comparisons of NER scores across BERT models on ADS data | BERT | astroBERT | astroBERT
---|---|---|---
|  | (MLM) | (NSP+MLM)
f1-score | 0.859 | 0.893 | 0.902
standard deviation | 0.014 | 0.028 | 0.014

##  5 Future Work

With astroBERT showing encouraging results on the preliminary ADS data, we are
currently building a more complete NER dataset covering more named entities
and papers. Our goal is to extend our services to automatically identify
entities such as facilities, and allow ADS users to unambiguously find papers
that mention them.

Beyond improving ADS services, astroBERT will be made publicly available to
the astrophysics community through the Huggingface Model Hub , which provides
free access to many state-of-the-art language models.

## References

  * Beltagy et al. (2019) Beltagy, I., Lo, K., & Cohan, A. 2019, arXiv e-prints. 1903.10676
  * Bird et al. (2009) Bird, S., Klein, E., & Loper, E. 2009, Natural language processing with Python: analyzing text with the natural language toolkit (" O’Reilly Media, Inc.")
  * Devlin et al. (2018) Devlin, J., Chang, M.-W., Lee, K., & Toutanova, K. 2018, arXiv e-prints. 1810.04805
  * Kurtz et al. (1993) Kurtz, M. J., Karakashian, T., Grant, C. S., Eichhorn, G., Murray, S. S., Watson, J. M., Ossorio, P. G., & Stoner, J. L. 1993, in Astronomical Data Analysis Software and Systems II, edited by R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, vol. 52 of Astronomical Society of the Pacific Conference Series, 132
  * Micikevicius et al. (2017) Micikevicius, P., Narang, S., Alben, J., Diamos, G., Elsen, E., Garcia, D., Ginsburg, B., Houston, M., Kuchaiev, O., Venkatesh, G., & Wu, H. 2017, arXiv e-prints. 1710.03740
  * Vaswani et al. (2017) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, Ł., & Polosukhin, I. 2017, in Advances in neural information processing systems, 5998
  * Wolf et al. (2019) Wolf, T., Debut, L., Sanh, V., Chaumond, J., Delangue, C., Moi, A., Cistac, P., Rault, T., Louf, R., Funtowicz, M., Davison, J., Shleifer, S., von Platen, P., Ma, C., Jernite, Y., Plu, J., Xu, C., Le Scao, T., Gugger, S., Drame, M., Lhoest, Q., & Rush, A. M. 2019, arXiv e-prints. 1910.03771