# Charge trapping and recombination in dipolar field of charged defect cluster
in silicon

Darius Abramavicius  Institute of Chemical Physics, Physics Department,
Vilnius University, Sauletekio al. 9-III  Juozas Vidmantis Vaitkus  Institute
of Photonics and Nanotechnology, Physics Department, Vilnius University,
Sauletekio al. 9-III

###### Abstract

Extensive irradiation of silicon crystal sensors by high energy particles in
e. g. accelerators yield defect clusters of different types. Trapping of
electrons and holes result in extended internal electric fields that drive
remaining free charges. The question whether these internal electric fields
affect the experimental observables, e.g. recombination process and lifetime
of free charges is the main focus of this paper. Including the drift and
diffusion of electrons and holes we calculate the recombination rate in a
cubic sample with a single dipolar cluster of defects. It is shown that the
large effect on charge lifetime is to be expected when charge diffusion length
during the charge lifetime is comparable to the dimensions of the cluster. If
the diffusion length exceeds the cluster size, the cluster barely affects the
recombination rate.

##  I Introduction

Silicon-based semiconductor sensors and detectors are widely used as high
energy (HE) particle detectors [1]. The fundamental principal design of
detectors relies on generation of mobile charges by HE particles. These mobile
charges are being extracted and detected as electrical signals in e. g.
avalanche diode configurations, such as low gain avalanche detectors
(LGADs)[2]. However, in specific experiments the irradiation energy may soar
as high as 13-14 TeV (in e. g. Large Hadron Colliders - LHCs at CERN); the HE
particles strongly interact with the crystal lattice and their kinetic energy
is sufficient to knock-out atoms from their rest positions. In the simplest
case two types of structural defects are then generated in the bulk: vacancies
and interstitials, which shorten the devices’ lifetime and the detectors
degrade [3, 4]. It has been proposed that atoms, knocked out from their
lattice positions have sufficient energy to further on knock out other atoms,
thus creating an avalanch-like defect generation of vacancies and interstitial
atoms [5]. Additional to that, as interstitials are created in the lattice
they could relax into several distinct lattice positions. The interstitial
defects tend to migrate in the crystal, recombine with vacancies until the
defects relax into quasi-stable configuration after impact with the HE
particle. As a result, the complex pattern of extended defect clusters should
be highly probable [5, 6]. Clustering of defects is also induced by ablation
process [7]. It is highly probable that groups of interstitials and vacancies
organize into nano-size effective clusters where high content of vacancies and
interstitials are distributed non homogeneously [5]. Resulting defect clusters
may demonstrate dipolar-like characteristics.

At least three types of interstitials are stable in Si crystal [8, 9, 6] with
their unique electronic properties. Electronic properties of point defects are
essentially determined by position of electron energy at the defect [10, 11].
Some of the point defects can trap electrons, others can trap holes, while
some are recombination centers [12]. Theoretical determination of isolated
point defects can be calculated by quantum methods at high accuracy at the
atomic level [13]. Clustering of various defects becomes impossible to analyze
at the same level. Nonhomogeneous distribution of point defects in defect
cluster accompanied by trapping of charges creates long-range internal lattice
fields that affects dynamics of the remaining charges on long distances [14].
Strong correlation between the type of defect clustering and the type of
irradiation has been established [15, 7]. Defect clusters are most probable
after irradiation by neutrons, while irradiation by charged particles may lead
to sparse point defects [5]. Wide distribution of types of defect clusters
should be thus expected. It has been estimated that only about a third of
point defects become polarized, i. e. trap one or other type of charges, while
the point defects that effectively trap both charges, become recombination
centers [16]. Hence, the largest amount of defects turn into recombination
centers, which can be assumed to be distributed uniformly in the vicinity of
(and within) the cluster. In the simplified picture, the charged defects thus
create polarized medium for free charges and it thus becomes a significant
factor, which drives the free charges that become captured by neutral
recombination centers. The question, whether such internal fields affect the
charge recombination and at what conditions, is targeted in this paper.

Experimental observations of recombination kinetics after photoexcitation
allow estimation of charge recombination rates and demonstrate that the charge
recombination time becomes directly proportional to the overall doze of
irradiation [17, 14]. It is still debatable whether this variation of charge
lifetime is due to increase of overall defect density in the bulk or due to
spatial charge accumulated dipolar-like defect clusters because presumably
they are responsible for variation of charge lifetime as the function of
irradiation intensity [18]. In this work we focus on this question by computer
simulations of the recombination process by considering diffusion and drift
processes of free charges in the vicinity of the simplest possible defect
cluster with the restriction that the density of point defects is assumed to
be not too high so that the bulk properties, such as mobility and diffusion
coefficients of charges are not affected. This is certainly an approximation,
however, it allows to highlight the effect that the dipolar field of the
defect cluster makes on the charge recombination. Simulations of the
recombination process in a box with a single cluster demonstrate that due to
large diffusion distances of free electrons and holes, a nano-sized defect
clusters with dipolar internal fields weakly affect charge recombination rates
when the density of recombination centers is low. Influence of the internal
fields becomes considerable when the charge lifetime is artificially reduced
so that the diffusion length becomes comparable to the size of the cluster.

##  II Theory of charge density evolution

Consider a single defect cluster that is isolated in the crystal. Electron-
hole recombination process in the defect cluster involves (at least) two
opposite free charges, both mobile, both affected by the internal field of
trapped charges. Their dynamics has to be followed simultaneously. They must
interact via Coulomb field and with the field of trapped charges in the defect
cluster. We restrict the problem to the semiclassical approach in the
overdamped limit (relaxation-time approximation), where the charges are
characterized by their effective masses and relaxation times. The main
parameters characterizing their dynamics are the diffusion coefficients and
mobilities.

A pair of free charges can be described in space by a distribution function,
$\rho(\bm{r}_{e},\bm{r}_{h},t)d^{6}r$ in six-dimensional space, denoting the
probability that electron is at position $\bm{r}_{e}$ in the volume element
$d^{3}r$, while the hole is at $\bm{r}_{h}$ in volume element $d^{3}r$.
Dynamics is governed by diffusion and drift processes described by Fokker-
Planck equation of motion [19]:

$\displaystyle\dot{\rho}(\bm{r}_{e},\bm{r}_{h},t)$ $\displaystyle=$
$\displaystyle
D_{e}\Delta_{e}\rho(\bm{r}_{e},\bm{r}_{h},t)+D_{h}\Delta_{h}\rho(\bm{r}_{e},\bm{r}_{h},t)$
(1)
$\displaystyle-\nabla_{e}\cdot\left(\rho(\bm{r}_{e},\bm{r}_{h},t)\bm{f}_{e}(\bm{r}_{e},\bm{r}_{h},t)\right)$
$\displaystyle-\nabla_{h}\cdot\left(\rho(\bm{r}_{e},\bm{r}_{h},t)\bm{f}_{h}(\bm{r}_{e},\bm{r}_{h},t)\right)$
$\displaystyle-S(\bm{r}_{e},\bm{r}_{h})\rho(\bm{r}_{e},\bm{r}_{h},t).$

Here the dot denotes time derivative, we also assume the homogeneous diffusion
with diffusion coefficients $D_{e}$ for electrons and $D_{h}$ for holes,
$\bm{f}_{e}(\bm{r}_{e},\bm{r}_{h},t)$ is the drift force acting on electrons,
while $\bm{f}_{h}(\bm{r}_{e},\bm{r}_{h},t)$ is the corresponding force acting
on holes. Correspondingly, $\Delta_{i}$ and $\nabla_{i}$ are the 3D Laplace
and _nabla_ differential operators acting on electron and hole coordinates,
respectively. $S(\bm{r}_{e},\bm{r}_{h})$ describes the recombination. All
these parameters are described below.

The combined probability density can be used to define separate electron and
hole distributions. We can define the electron distribution function by
integrating out the hole variables

$\rho_{e}(\bm{r}_{e},t)=\int_{V}d^{3}\bm{r}_{h}\rho(\bm{r}_{e},\bm{r}_{h},t);$
(2)

and similarly for the hole distribution function

$\rho_{h}(\bm{r}_{h},t)=\int_{V}d^{3}\bm{r}_{e}\rho(\bm{r}_{e},\bm{r}_{h},t),$
(3)

while the total probability (or density) of electrons and holes being in the
system is given by

$n(t)=\int_{V}d^{3}\bm{r}_{h}\int_{V}d^{3}\bm{r}_{e}\rho(\bm{r}_{e},\bm{r}_{h},t).$
(4)

The total probability density $n(t)$ is the quantity of interest when
characterizing the charge recombination process. Following time evolution of
$n(t)$ allows to define the recombination rate.

Above definitions yield separate, while coupled, equations for electron and
hole distributions by integrating equation 1. Using notation $(i,-i)=(e,h)$ or
$(h,e)$ we find

$\displaystyle\dot{\rho}_{i}(\bm{r},t)$ $\displaystyle=$ $\displaystyle
D_{i}\Delta\rho_{i}(\bm{r},t)$ (5)
$\displaystyle-\nabla_{r}\cdot\left(\int_{V}d^{3}\bm{r}^{\prime}\rho_{i}(\bm{r},t)\frac{\rho_{-i}(\bm{r}^{\prime},t)}{n(t)}\bm{f}_{i}(\bm{r},\bm{r}^{\prime},t)\right)$
$\displaystyle-\int_{V}d^{3}\bm{r}^{\prime}S(\bm{r},\bm{r}^{\prime})\frac{\rho_{-i}(\bm{r}^{\prime},t)}{n(t)}\rho_{i}(\bm{r},t).$

Here the first term on the right side of equation is due to diffusion, second
is responsible for electron-hole interaction and interaction with other
trapped charges, while the third term defines the recombination process.

To define explicit expression for the drift forces we recall that the drift
force in Fokker-Planck equation corresponds to equations of motion for a
particle with coordinate $x_{i}(t)$ in the form

$\dot{x}_{i}(t)=f_{i}(x_{1}...x_{N},t).$ (6)

The force function $f_{i}(x_{1}...x_{N},t)$ can be defined by considering all
internal Coulomb fields. For electrons and holes in the semiclassical
approximation we have to include inter-particle Coulomb attraction as well as
internal field-induced attraction or repulsion by the trapped charges. All
these fields may be expressed as gradients of internal potential created by
charged species:

$\displaystyle U(\bm{r}_{e},\bm{r}_{h})$ $\displaystyle=$
$\displaystyle-\chi\frac{1}{\left|\bm{r}_{e}-\bm{r}_{h}\right|}$ (7)
$\displaystyle-\chi\sum_{d}Q_{d}\left(\frac{1}{\left|\bm{r}_{e}-\bm{R}_{d}\right|}-\frac{1}{\left|\bm{r}_{h}-\bm{R}_{d}\right|}\right)$

Here the first term corresponds to the electron-hole attraction of mobile
charges with $\chi=e^{2}/(4\pi\epsilon\epsilon_{0})$ being the amplitude of
the Coulomb potential, other two terms are due to static trapped charges:
$Q_{i}$ is the charge of the defect $d$ given in terms of elementary electron
charge $e$, while $\bm{R}_{d}$ \- its position. Considering the times longer
than the scattering time, averaged motion of the charges in the overdamped
limit can be described by equation in the form of eq 6 with $i=e,h$, what
defines the drift for $i-$ th particle [20]

$\dot{\bm{r}}_{i}=-\frac{\tau_{i}}{m_{i}^{\ast}}\nabla_{i}U(\bm{r}_{i},\bm{r}_{-i})\equiv\bm{f}_{i}(\bm{r}_{e},\bm{r}_{h}).$
(8)

Here we find $m_{i}^{\ast}$ \- the effective mass, $\tau_{i}$ is the
relaxation time. Also notice that the parameters, which enter the drift force
in Eq. 8 can be related via Drude expression, relating the relaxation time,
effective mass, and mobility: $\mu_{i}=e\tau_{i}/m_{i}^{\ast}$; $e$ is the
electron charge. Consequently, the ratio $\tau_{i}/m_{i}^{\ast}$ reduces to
$\tau_{i}/m_{i}^{\ast}=\mu_{i}/e$. Then Einstein relation
$D_{i}=\mu_{i}k_{B}T$ allows to reduce the drift amplitude to the thermal
energy. Further on we use the standard notation for the inverse thermal energy
$\beta=(k_{B}T)^{-1}$ and can write $\tau_{i}/m_{i}^{\ast}=\beta D_{i}/e$.

The last term in eq. 5 characterizes the recombination. Assuming that _the
process_ is time independent and local in space, i. e. the electron and hole
must meet in space on a specific recombination center in order the
recombination to occur; and that the probability of finding a recombination
center is distributed uniformly in space, we obtain a simple expression for
the recombination function
$S(\bm{r}_{e},\bm{r}_{h})\equiv\bar{R}\delta(\bm{r}_{e}-\bm{r}_{h})$, with
$\bar{R}$ \- the mean on-site recombination rate.

Combining the recombination rate and drift force into Eq. 5 we obtain the
real-space equations of motion describing drift, diffusion and recombination
of electrons and holes in the cluster:

$\displaystyle\dot{\rho}_{i}(\bm{r},t)=D_{i}\Delta\rho_{i}(\bm{r},t)-\bar{R}\rho_{i}(\bm{r},t)\frac{\rho_{-i}(\bm{r},t)}{n(t)}$
$\displaystyle\qquad-\frac{\beta
D_{i}\chi}{e}\int_{V}d^{3}\bm{r}^{\prime}\frac{\rho_{-i}(\bm{r}^{\prime},t)}{n(t)}\frac{\bm{r}-\bm{r}^{\prime}}{\left|\bm{r}-\bm{r}^{\prime}\right|^{3}}\cdot\nabla\rho_{e}(\bm{r},t)$
$\displaystyle\qquad+\frac{\beta
D_{i}\chi}{e}q_{i}Q_{d}\sum_{d}\frac{\bm{r}-\bm{R}_{d}}{\left|\bm{r}-\bm{R}_{d}\right|^{3}}\cdot\nabla\rho_{e}(\bm{r},t)$
(9)

where $q_{i}$ is the charge of the $i-$th particle in terms of elementary
electron charge. Integration volume $V$ corresponds to the system under
consideration. Notice that resulting equations are linear in electron and hole
density so the absolute density does not play any role.

In the following we study the recombination process in the simplest neutral
system, where the internal field is of the dipolar type that can be
represented by positioning two opposite stationary charges at the given
positions. It is convenient to set orientation of the coordinate system with
respect to the dipole: the symmetry axis of the dipole is oriented along $z$
axis, positive charge of the dipole is positioned at
$\bm{z}_{0}\equiv(0,0,z_{0})$, while negative - at $-\bm{z}_{0}$. The absolute
charge center is then at the origin.

##  III Stationary equation in reciprocal space

The obtained equations of motion depend on spatial as well as on time
derivatives. Numerically the high order mixed derivatives are poorly
represented by a finite difference. Certain derivatives can be eliminated by
transforming the problem into the reciprocal $\bm{k-}$ space by applying
Fourier transform. Using definition

$\rho_{i}(\bm{r},t)=\int\frac{d^{3}\bm{k}}{(2\pi)^{3}}e^{i\bm{k}r}\sigma_{i}(\bm{k},t)$
(10)

we have a simple conversion for the core of the Coulomb potential:

$\frac{1}{|\bm{r}|}\to\frac{4\pi}{|\bm{k}|^{2}}.$ (11)

Additionally, at zero wavevector $\sigma_{i}(0,t)=\sigma_{-i}(0,t)=n(t)$ is
the total charge density. Consequently, evaluating the recombination kinetics
becomes equivalent to following the zero-value wave vector amplitude. The
resulting equation in $\bm{k}-$ space becomes simplified:

$\displaystyle\dot{\sigma}_{i}(\bm{k},t)=-\bm{k}^{2}D_{i}\sigma_{i}(\bm{k},t)$
$\displaystyle\ +\frac{CD_{i}L}{8\pi^{3}n(t)}\int
d^{3}\bm{k}^{\prime}\frac{\bm{k}\cdot\bm{k}^{\prime}}{|\bm{k}^{\prime}|^{2}+\eta^{2}}\sigma_{-i}(\bm{k}^{\prime},t)\sigma_{i}(\bm{k}-\bm{k}^{\prime},t)$
$\displaystyle\ +\frac{iq_{i}QCD_{i}L}{4\pi^{3}}\int
d^{3}\bm{k}^{\prime}\frac{\bm{k}\cdot\bm{k}^{\prime}}{|\bm{k}^{\prime}|^{2}+\eta^{2}}\sin\left(\bm{k}^{\prime}\bm{z}_{0}\right)\sigma_{i}(\bm{k}-\bm{k}^{\prime},t)$
$\displaystyle\ -\frac{\bar{R}}{8\pi^{3}n(t)}\int
d^{3}\bm{k}^{\prime}\sigma_{-i}(\bm{k}^{\prime},t)\sigma_{i}(\bm{k}-\bm{k}^{\prime},t).$
(12)

Here we denoted $C=4\pi\beta\chi/(eL)$ \- the dimensionless temperture-
dependant amplitude of the Coulomb interaction. The first term in the right
hand side of Eq. 12 is the diffusion-related damping, then the integrals are
due to Coulomb drift force and final is the recombination term. The drift-
related integral vanishes for $\bm{k}=0$ or $\bm{k}^{\prime}=0$ since the
corresponding amplitudes represent homogeneous charge densities having
translational invariance. To emphasize this and to avoid divergencies, an
additional $\eta^{2}$ shift with $\eta\ll\Delta k$ is added in the
denominators (the Coulomb potential function in real space becomes equal to
the Yukawa potential).

While equations do not involve spatial derivatives (i. e. diffusion has been
transformed into derivative-less form), the tripple integrations over momentum
$\bm{k}^{\prime}$ cannot be reduced, so full numerical propagation at every
time is computationally very expensive. Additionally, notice that the
diffusion-related terms in the equation are proportional to $\bm{k}^{2}$. So
the diffusion-related term diverges for large wave vectors requiring
decreasing of time step. To avoid explicit time propagation we consider the
stationary dynamics, when the charge distribution relaxed in space for the
given steady decay of nonstationary populations. Then the recombination rate
can be evaluated with the Ansatz that the electron (and hole) density decays
exponentially. We thus request

$\sigma_{i}(\bm{k},t)=\bar{\sigma}_{i}(\bm{k})e^{-\gamma t}$ (13)

With the initial condition that at time zero, the probability to find electron
and hole is equal to $1$, we thus have the initial condition corresponding to
the boundary value $\bar{\sigma}_{i}(\bm{k}=0)\equiv n(t=0)=1$. For
$\bm{k}\neq 0$ from Eq. 12 we obtain the integral equation for the spatial
charge distribution in discretized form

$\displaystyle\bar{\sigma}_{i}(\bm{k})=\frac{C}{(\bm{k}L)^{2}}\sum_{\bm{k}^{\prime}}\frac{\bm{k}\cdot\bm{k}^{\prime}}{|\bm{k}^{\prime}|^{2}}\bar{\sigma}_{-i}(\bm{k}^{\prime})\bar{\sigma}_{i}(\bm{k}-\bm{k}^{\prime})$
(14)
$\displaystyle\qquad+\frac{2iq_{i}QC}{(\bm{k}L)^{2}}\sum_{\bm{k}^{\prime}}\frac{\bm{k}\cdot\bm{k}^{\prime}}{|\bm{k}^{\prime}|^{2}}\sin\left(\bm{k}^{\prime}\bm{z}_{0}\right)\bar{\sigma}_{i}(\bm{k}-\bm{k}^{\prime})$
$\displaystyle\qquad+\frac{\gamma}{\bm{k}^{2}D_{i}}\bar{\sigma}_{i}(\bm{k})-\frac{\gamma_{h}}{\bm{k}^{2}D_{i}}\sum_{\bm{k}^{\prime}}\bar{\sigma}_{-i}(\bm{k}^{\prime})\bar{\sigma}_{i}(\bm{k}-\bm{k}^{\prime}),$

while $\bm{k}=0$ in Eq. 12 yields the decay rate of the charge density

$\displaystyle\gamma$ $\displaystyle=$
$\displaystyle\gamma_{h}\sum_{\bm{k}^{\prime}}\bar{\sigma}_{-i}(\bm{k}^{\prime})\bar{\sigma}_{i}(-\bm{k}^{\prime}).$
(15)

We also introduced $\gamma_{h}=\bar{R}/V$, being the recombination rate for
the homogeneous system (see next section). In the discretized form all terms
are dimensionless. These equations can be solved iteratively starting from e.
g. the homogeneous distribution.

##  IV Results

Homogeneous distribution of charges constitutes the reference point which
represents the case when the internal defect field is not present. In that
case the combined 6D normalized density can be written as

$\rho^{(hom)}(\bm{r}_{e},\bm{r}_{h},t)=\frac{n^{(hom)}(t)}{V^{2}},$ (16)

where $V$ is the volume under consideration. Consequently

$\rho_{e}^{(hom)}(\bm{r},t)=\rho_{h}^{(hom)}(\bm{r},t)=\frac{n^{(hom)}(t)}{V}$
(17)

and in reciprocal space $\sigma_{i}^{(hom)}(\bm{k},t)=n^{(hom)}(t)$ for
$\bm{k}=0$ and zero otherwise. For the homogeneous system the translational
invariance holds (all Coulomb interactions vanish). Equation 12 can then be
integrated what yields the recombination rate, $\gamma_{h}=\bar{R}/V$, of the
homogeneous distribution

$\frac{\dot{n}^{(hom)}(t)}{n^{(hom)}(t)}=-\gamma_{h}.$ (18)

Notice that it rescales the on-site rate $\bar{R}$ by the volume of the
system. This is related to the fact that when the volume of the system is
increased, the probability for the two carriers to meet on the single site,
which is proportional to $V^{-2}$, is reduced. The scaling in Eq. 18 is
$V^{-1}$ because for the constant density of recombination centers, their
number increases linearly with volume $V$. However, in our simulations the
rate $\gamma_{h}$ should be considered as tuning parameter since the density
of recombination centers in defect clusters is not known.

The homogeneous distribution function is taken as the initial condition for
numerical solution of Eq. 14. Starting from
$\bar{\sigma}_{i}^{(hom)}(\bm{k}=0)=1$ (and zero otherwise) we solved Eq. 14
iteratively together with Eq. 15. For calculations we assume the widely
available set of parameters of pure Si semiconductor. We used a length scale
of 1$\mu m$ and time scale of $1\mu s$. The modelling system is taken as a box
of volume $V=L^{3}=1\,(\mu m)^{3}$. The parameters, that enter the drift force
in Eq. 8, specifically, the effective masses and relaxation times, have been
reduced to diffusion constants and mobilities and finally to thermal energy in
Eq. 14. We use the standard set of Si characteristics at room temperature, i.
e. the experimental mobility values, $\mu_{e}=1400\ cm^{2}/(Vs)$,
$\mu_{h}=450\,cm^{2}/(Vs)$, and diffusion constants, $D_{e}=36\,cm^{2}/s$,
$D_{h}=12\,cm^{2}/s$. Additionally the Coulomb field is scaled by Si
dielectric constant $\epsilon=11.68$. The box is discretized into $40\times
40\times 40$ points. The resulting effective cubic lattice constant is
$a=1/40$ $\mu m$. In the reciprocal space the integration step $\Delta
k=2\pi/L=2\pi$ $(\mu m)^{-1}$.

The last parameter that has to be considered is the trapped charge in the
defect $Q$. For the defect cluster we assume that the cluster contains net
charge $Q>e$, i. e. the defect cluster can trap more then a single charge. The
total trapped charge is the property of the defect and is equal to the number
of trapping sited in the cluster. A single point defect can in principle
accommodate only 1-2 e. However, in the extended defect cluster we may
consider accumulation of several point defects in a single region.
Consequently, the charge $Q$ becomes a unique number for each particular
cluster (Eq. 14). In present study we chose this number arbitrarily and set
$Q=4e$.

Figure 1: Convergence of decay rate with iterations. Dipole charges $Q=4e$,
dipole separation $2z_{0}=4a\equiv 100$ nm, homogeneous recombination rate
$\gamma_{h}=1$ ns-1. Figure 2: Distribution of electron density in real space
$\bar{\sigma}_{e}(\bm{r})$ at a specific plain. Dipole charges $Q=4e$, dipole
separation $2z_{0}=4a\equiv 100$ nm, homogeneous recombination rate
$\gamma_{h}=1$ ns-1. Figure 3: Distribution of electron density
$\bar{\sigma}_{e}(\bm{r})$ along $z$ axis at $x,y=0$. Dipole charges $Q=4e$,
dipole separation $2z_{0}=4a\equiv 100$ nm, homogeneous recombination rate
$\gamma_{h}=1$ ns-1.

First, we consider the slow recombination regime when the probability of
recombination centers is low. This is the case when for our chosen size of the
cluster we take $\gamma_{h}<(\Delta k)^{2}D_{i}$ in Eq. 14. This corresponds
to the case when the diffusion-related charge travel time over the system
length $L$ is shorter than the recombination time $\gamma_{h}^{-1}$. For the
electron and hole the diffusion coefficients we find $(\Delta k)^{2}D_{e}=$
$142$ ns-1 for electrons and $(\Delta k)^{2}D_{h}=$ $47$ ns-1 for holes. To
represent this case we assume that the homogeneous recombination rate
parameter $\gamma_{h}=1$ ns-1. Fig. 1 demonstrates iterative convergence of
Eq. 14 towards the solution for this type of system. The convergence is mostly
exponential and can be safely truncated after $\sim$20-th iteration. The
resulting electron distribution is presented in Fig. 2 at the plains $xy$,
$xz$, while along $z$ axis at $(x,y)=0$ is shown in Fig 3. As expected, the
positive trapped charge attracts the electron density, while the negative
trapped charge repels it. The distribution becomes rotationally symmetric
around $z$ axis. It is worth noting that due to relatively large diffusion
constants, the distribution is much wider at distances larger than the dipolar
charge separation. The integral charge around the dipole (while with low
amplitude) makes significant contribution to the recombination process, but
the conditions are essentially quite similar to the homogeneous conditions and
the effect of the dipolar field of the trapped charges is small. The
calculated recombination rates, presented in Fig. 4, demonstrate that while
the effect of the dipole is weak, there is still dependence on the separation
of charges in the dipole: as the separation is increased, the recombination
rate is reduced. This is due to the fact that the electrons and holes become
condensed in separate spatial regions so their overlap becomes smaller.

Similar result is obtained for the homogeneous recombination rate
$\gamma_{h}=$10 ns-1 (100 ps decay time), which also satisfies
$\gamma_{h}<(\Delta k)^{2}D_{i}$. Since the rate converges similarly and the
distributions obtained are similar to the previous case, we do not present
charge distributions. Only the final recombination rates are presented in Fig.
4 by red crosses. We find that the relative dependence on dipole separation is
slightly stronger, but still the effect is weak.

Figure 4: Normalized variation of recombination rates for dipole charges
$Q=4e$ as a function of separation between dipole charges when $\gamma_{h}=1$
ns-1 (green dots) and (100 ps)-1 (red crosses).  Figure 5: Convergence of
charge distribution rate with iterations. Dipole charges $Q=4e$, dipole
separation $2z_{0}=4a\equiv 100$ nm, homogeneous recombination rate
$\gamma_{h}=100$ ns-1. Figure 6: Distribution of electron density
$\bar{\sigma}_{e}(\bm{r},t)$ at a specific plain. Dipole charges $Q=4e$,
dipole separation $2z_{0}=4a\equiv 100$ nm, homogeneous recombination rate
$\gamma_{h}=100$ ns-1.  Figure 7: Distribution of electron (red) and hole
(blue) density $\bar{\sigma}_{e}(\bm{r},t)$ along $z$ axis at $x,y=0$. Dipole
charges $Q=4e$, dipole charge separation $2z_{0}=4a\equiv 100$ nm, homogeneous
recombination rate $\gamma_{h}=100$ ns-1. Figure 8: Calculated recombination
rates for dipole charges $Q=4e$ as a function of separation between dipole
charges when $\gamma_{h}=100$ ns-1.

Completely different result is obtained when the homogeneous decay is fast:
$\gamma_{h}=100$ ns-1 (10 ps decay time). In that case we have (for holes) the
diffusion term $\bm{k}^{2}D_{h}$ that can be larger or smaller than
$\gamma_{h}$ depending on the amplitude of wave vector $\bm{k}$. In Fig 5 we
show the convergence of the solution. First of all notice additional
complexity in convergence. The dependence is non-exponential - the rate
evolves in step-like fashion as $\gamma$ value traverses specific discrete
values of $\bm{k}^{2}D$. Due to numerical approach, the result does not relax
to the complete stable solution and small fluctuating behavior remains
visible. However, the rate can be estimated. For estimation of the rate we
take the lower limit, where the solution is mostly stabilized.

Charge distribution becomes more localized compared to the previous case: in
Fig. 2 the distribution was homogeneous except at the positions of trapped
charges. In the present case in Fig. 6 we find that the long-range charge
attraction (due to positive charge) and charge repulsion (due to negative
charge) create depleted regions. That corresponds to the situation where the
charges recombine faster than they escape the dipolar wells. Very large
variations of charge density are obtained along $z$ axis as shown in Figs. 6
and 7. The competition between diffusion, recombination and drift can be
observed. At $z>0$ the positive defect charge attracts the electron density
from the long distances $z>0$ as the Coulomb potential is long-range. This
shifts electron and hole distributions along $z$ axis in opposite directions
so the recombination rate is reduced. Very close to the defect positive charge
(at $z\sim 2a$), while charge density has spikes, the opposite charges are
repelled out from the nearest neighborhood and the recombination happens
slowly.

Such large variation of electron density has a strong effect on the
recombination rate. In Fig. 8 we show that the recombination rate is reduced
from the homogeneous case by 20% - much larger value compared to Fig. 4. The
effect starts from smallest $z_{0}$ values, while for $2z_{0}>150$ nm
dependence of the rate is almost constant.

##  V Discussion and conclusions

Presented model allows determining the effect of internal defect cluster
fields on the recombination rate (or time). The model allows to assign a decay
rate constant to the highly non-stationary charge distribution decay process.
This becomes possible by using the finite box with periodic boundary
conditions ($k$-space solution). The approach, thus, allows to characterize
the vicinity of the defect cluster. Considering the cluster without periodic
boundary conditions it would be impossible to obtain a single rate constant by
this method because the system would be described by non-exponential decay
process as the decay rate close to the cluster and far from it are essentially
different. So the present model is approximate while allows to demonstrate the
effect that local fields induce on the recombination process.

According to the presented calculations, the most significant parameter
describing charge evolution and recombination in the defect cluster is the
diffusion coefficient. It relates the diffusion length $l$ and diffusion time
$\tau$ via relation $D_{i}=l_{i}^{2}/\tau_{i}$. By taking the diffusion time
to be equal to the lifetime, we obtain the total diffusion length, i. e. the
distance a charge can traverse during its lifetime by diffusion. The total
diffusion length is the parameter that should be correlated with the size of
the defect cluster. In our analysis we used a box with edge length $L=$1 $\mu
m$. If the total diffusion length is much longer than the defect cluster size,
the charge will be able to escape the cluster (and the Coulomb field of the
dipole): the cluster will have negligible effect on the charge recombination.
Diffusion constant for holes in silicon are of order of $10$ cm2/s$\equiv
1000$ ($\mu$m)2/$\mu$s, consequently, the hole covers $\sim 10$ ($\mu$m)
before recombination takes place if the recombination time is $100$ ns. In
this case, if the defect cluster extends only over few hundred nm as it was
assumed in simulations, its effect on the recombination is almost vanishing,
as is observed in Fig. 4. In the other limit, if the lifetime is short - few
ps, the total diffusion length corresponds to few tens of nm and that becomes
comparable to the size of the defect cluster, and the effect of the cluster
dipolar field on the recombination time becomes considerable.

Defining the size of the defect cluster is also not trivial. As Coulomb
potential is long range, size of the cluster does not necessarily corresponds
with the distance between point defects in the cluster. For the dipolar
internal fields, the Coulomb potential energy of the dipole decays
approximately as $R^{2}$, where $R$ is the distance from the dipole. Comparing
the cluster-induced Coulomb potential energy (last two terms in Eq. 7) to
$k_{B}T$ we can define the effective dipolar capture distance, and the
corresponding cluster size
$L_{d}=4z_{0}\sqrt{Q(8\pi\epsilon\epsilon_{0}z_{0}k_{B}T)^{-1}}$, where
$2z_{0}$ is the separation of dipolar charges. In our model for dipole charge
separation of $100$ nm and $Q=4e$ we can find $L_{d}\approx 100$ nm. The
relevant size of the cluster that should be considered for the diffusion
process is therefore even smaller than we have discussed above.

In this discussion we are not restrictive about homogeneous system
recombination rate. The silicon crystals are never pure and various types of
impurities create recombination centers. The homogeneous recombination rate
therefore varies in various crystal preparation conditions. Additionally, in
order to isolate the effect of dipolar fields we neglected the possibility
that HE irradiation which creates trapping centers (responsible for creation
of internal fields), also generates additional recombination centers in the
same area. The local “homogeneous” recombination rate in the vicinity of the
defect clusters therefore could be quite different compared to the bulk values
of unaffected silicon. Therefore, high recombination rates in the vicinity of
clusters and strong influence of the internal cluster fields on the
recombination process is highly probable.

Diffusion process is quite easy to evaluate and estimate, while the drift
force depends on various additional parameters: the amount of trapped charge,
$Q$, effective separation of dipole charges (which could be associated with
the size of the defect cluster, as discussed above). Consider holes as their
diffusion constant is smaller. The quantity
$\frac{4\pi\chi\mu_{i}}{V}\frac{Q}{e}$, which scales the Coulomb fields in Eq.
9 (notice that $d^{3}k=(2\pi)^{3}/V$) and equal to $\frac{4\pi\chi
D_{i}}{Vk_{B}T}\frac{Q}{e}=278\mu s^{-1}$ with $Q=4$ and at room temperature
has the dimension on $s^{-1}$. Parameter $\bm{k}^{2}D_{i}$ has the same
dimension $s^{-1}$ and can be directly compared. Taking the smallest
$\bm{k}\equiv dk=2\pi/L$ we find $\bm{k}^{2}D_{h}=47400\ \mu s^{-1}$.
Diffusion coefficient of electrons is even larger. The Coulomb drift forces
acting on electrons and holes are thus considerably weaker than the diffusion-
related gradient forces (including the large dielectric constant). This
demonstrates that the diffusion process is extremely fast in Si compared to
microscopic internal field-induced gradient forces due to trapped charges by
defect clusters. Additionally it should be taken into account that the Coulomb
field is long-range but quickly decays on small scale. The thermal activation
energy $k_{B}T$ overcomes the Coulomb attraction energy for trapped charge
$Q=4e$ at a distance of $20$ nm, hence, only at such small distances the
mobile charges cannot escape the potential well of trapped charge. In our
modeling the box size is $\sim 1$ $\mu$m, consequently, the large amount of
charge is barely affected by Coulomb attraction forces.

However, this study of a single dipolar internal field applies when defect
concentration is relatively small. If the defect concentration is high, the
internal fields are no more dipolar, bulk properties like mobility, diffusion
constant, effective masses may be extensively affected. As the parameters
become different, the eventual recombination rate can be strongly affected.
However, our model with the bulk diffusion constant being independent on the
defect concentration allows to highlight how the dipolar internal field
redistributes the mobile charges so that charge recombination becomes altered
solely due to internal field effects.

Considering realistic samples we should keep in mind that high energy
irradiation creates clusters of defects with various configurations of trapped
charges. In Figs. 4 and 8 we demonstrated how the rate of recombination
changes with the simplest model of trapped dipole. When the distance between
trapped charges increases, the whole picture would finally reflect the case of
independent positively and negatively charged point defects scattered in the
bulk. Since the governing parameter is the diffusion length of a charge, the
transition to independent charged defects is directly related to the total
diffusion length.

The homogeneous recombination rate, that we consider in our simulations is
barely observable quantity in experiments due to variability of types of
recombination centers in irradiated silicon. As we have demonstrated the large
effect of the dipolar field on recombination rate emerges when the diffusion
length is comparable with the dimension of the defect cluster. This can be
achieved when the size of the defect cluster is large. As we have discussed,
for the recombination rates of $1$ $n$s, the effective diffusion length is $1$
$\mu$m. If the defect clusters are of the same dimensions we will have
conditions that the charges cannot escape the internal defect fields during
their lifetime and the defect cluster is the determining factor governing the
recombination rate. Theoretically this is easily demonstrated by making order-
by-order solution of the dynamic equation as demonstrated in Appendix A. When
the dimensions of the cluster become comparable to the diffusion length (by
varying parameter $\gamma_{h}$) the effect is amplified as difference
$\bm{k}^{\prime 2}D_{i}-\gamma_{h}$ can get to zero in the denominator. Our
results thus demonstrate that the ratio of the charge diffusion length and the
size of clusters is the most important parameter that governs the sensitivity
of the recombination rate to the defect density.

Having this in mind, consider the case of independent homogeneously
distributed monopole trapped charges without clasterization, $i.$e.
$z_{0}=\infty$ to the case where defects are organized into clusters, i. e.
$z_{0}$ finite. Figures 4 and 8 have to be compared starting from large
$z_{0}$ and going to smaller values. We then observe _increasing_
recombination rate with dipolar-type cluster formation. That is how the
present model should be conveyed. Totally pure semiconductor crystals do not
exist and recombination usually takes place via point defects even before
irradiation. In that case, the uniform distribution of point defects would
yield the homogeneous recombination rate of such sample $\gamma_{hom}$. This
parameter is essentially proportional to the density of the recombination
defects, $N_{r}$, while a single defect is characterized by a charge capture
radius, $R_{c}$. Irradiation of semiconductor by HE particles creates clusters
of defects with additional recombination centers $M_{r}$ and trapping centers.
Since no foreign atoms are introduced, additional recombination defects are of
the same type as the original recombination defects, while their density is
larger in the clusters. We therefore should expect increase of $\gamma_{hom}$
in the vicinity of the clusters, i. e. we can consider the “homogeneous”
recombination rate to take very different values in various spatial positions.

Another questionable point is whether bulk crystal parameters should be
applicable to the defect clusters. Charge mobility in Si decays with the
concentration of dopants (down to $\sim$100 cm2/Vs for electrons when defect
concentration reaches $10^{19}$ cm-3)[21]. While overall concentration of
defects may not be very high, local concentration in the vicinity of defect
cluster may be very high. As we have discussed the lower diffusion
coefficients means that variation of the recombination rate at the cluster
could be appreciable even at low homogeneous relaxation rates $\gamma_{hom}$.
Our simulations thus correspond to the limiting case when density of defects
inside the cluster is low. We find that internal electric field of the cluster
becomes playing considerable effect at high concentration of point
recombination centers. If we assume that the charge mobility decreases
tenfold, the concentration of point recombination centers can be ten times
lower to observe the same strong effect. Electron mobilities of $\sim$100
cm2/Vs corresponds to the lower limit and thus should represent the case when
electrons diffuse over the point defects via tunneling. So these conditions
would correspond to the high concentration of trapping centers and that would
reflect the other limiting case. However, in the conditions of high
concentration of defects the structure of internal fields should be also
considered more complicated.

It should be stressed that such clustering of defects may be important for
various types of semiconductor high energy irradiation sensors. Defect
clustering creates local internal fields and conditions that are different
from the bulk. If they determine the bulk properties (e.g. by creating sinks
of charge density) this will be observed by bulk experiments.

Concluding, we have presented theoretical analysis of charge recombination in
Si crystal in the vicinity of trapped charge creating dipolar field starting
from the real space Fokker-Planck equation which describes charge evolution
including diffusion and drift forces as well as charge recombination. Using
the reciprocal space representation we have the governing equation in the
compact form, what allows iterative solution of the equation. Using the
stationary conditions we define the experimentally observed recombination rate
of the homogeneous system. We found that for slow recombination process in the
case of homogeneous conditions, inserting dipolar field minimally affects the
rate of homogeneous recombination. Large variation of the recombination rate
occurs when the homogeneous recombination rate is fast (on the order of tens
of picoseconds). The reason behind this effect is that the diffusion length of
charges during their lifetime becomes comparable to the extensions of defect
clusters.

##  Appendix A Order-by-order solution

In order to demonstrate the sensitivity of the effect on various parameters it
is instructive to perform analytically few steps. Hence, take Eq. 14 and start
from the homogeneous distribution

$\bar{\sigma}_{i}(\bm{k})=\left\\{\begin{array}[]{c}1,\ \bm{k}=0\\\ 0,\
\mathrm{otherwise}\end{array}\right.$ (19)

The first iteration starting from homogeneous distribution we find an
additional term in distribution function at $\bm{k}\neq 0$:

$\displaystyle\bar{\sigma}_{i}(\bm{k})$ $\displaystyle=$
$\displaystyle\frac{(2\pi)^{3}}{V}\frac{iq_{i}Q\chi\tau_{i}}{\pi^{2}m_{i}^{\ast}\left(\bm{k}^{2}D_{i}-\gamma_{h}\right)}\sin\left(\bm{k}\bm{z}_{0}\right)$
(20)

which gives the recombination rate

$\displaystyle\frac{\gamma-\gamma_{h}}{\gamma_{h}}$ $\displaystyle=$
$\displaystyle\frac{8Q^{2}\chi^{2}\tau_{h}\tau_{e}}{\pi
Vm_{h}^{\ast}m_{e}^{\ast}}\int\frac{-d\bm{k}^{\prime}\sin^{2}\left(\bm{k}^{\prime}\bm{z}_{0}\right)}{\left(\bm{k}^{\prime
2}D_{h}-\gamma_{h}\right)\left(\bm{k}^{\prime 2}D_{e}-\gamma_{h}\right)}.$
(21)

We thus find that at first iteration we obtain correction to the recombination
rate. Notice that $\gamma<\gamma_{h}$ if $D_{e}=D_{h}$. However, the situation
may become different if $D_{h}<D_{e}$ as in the case of Si crystal. In that
case the integral may yield negative value if $\bm{k}^{\prime
2}D_{h}<\gamma_{h}<\bm{k}^{\prime 2}D_{e}$. In this case the overall
recombination rate may possibly become larger than $\gamma_{h}$.

## Acknowledgements

This work was performed in the framework of the CERN RD50 collaboration. It
was funded by the Lithuanian Academy of Sciences, grant No CERN-VU-2021-2022.
Computations were performed on resources at the High Performance Computing
Center, “HPC Sauletekis” in Vilnius University Faculty of Physics.

## Author Declarations

The authors have no conflicts to disclose.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Hartmann [2017] F. Hartmann,  _Evolution of Silicon Sensor Technology in Particle Physics_ , 2nd ed. (Springer, 2017).
  * Moffat  _et al._ [2018] N. Moffat, R. Bates, M. Bullough, L. Flores, D. Maneuski, L. Simon, N. Tartoni, F. Doherty, and J. Ashby, Low Gain Avalanche Detectors (LGAD) for particle physics and synchrotron applications, J. Instrum. 13, C03014 (2018).
  * Lindström [2003] G. Lindström, Radiation damage in silicon detectors, Nucl. Instruments Methods Phys. Res. Sect. A Accel. Spectrometers, Detect. Assoc. Equip. 512, 30 (2003).
  * Gallrapp  _et al._ [2017] C. Gallrapp, M. Fernández García, S. Hidalgo, I. Mateu, M. Moll, S. Otero Ugobono, and G. Pellegrini, Study of gain homogeneity and radiation effects of Low Gain Avalanche Pad Detectors, Nucl. Instruments Methods Phys. Res. Sect. A Accel. Spectrometers, Detect. Assoc. Equip. 875, 27 (2017).
  * Radu  _et al._ [2015] R. Radu, I. Pintilie, L. C. Nistor, E. Fretwurst, G. Lindstroem, and L. F. Makarenko, Investigation of point and extended defects in electron irradiated silicon - Dependence on the particle energy, J. Appl. Phys. 117, 164503 (2015).
  * Centoni  _et al._ [2005] S. A. Centoni, B. Sadigh, G. H. Gilmer, T. J. Lenosky, T. Díaz de la Rubia, and C. B. Musgrave, First-principles calculation of intrinsic defect formation volumes in silicon, Phys. Rev. B 72, 195206 (2005) .
  * Aye  _et al._ [2021] M. M. Aye, E. Rivasto, M. Z. Khan, H. Rijckaert, E. Salojärvi, C. Haalisto, E. Mäkilä, H. Palonen, H. Huhtinen, I. Van Driessche, and P. Paturi, Control of the nanosized defect network in superconducting thin films by target grain size, Sci. Rep. 11, 6010 (2021).
  * Newman [1982] R. C. Newman, Defects in silicon, Reports Prog. Phys. 45, 1163 (1982).
  * Leung  _et al._ [1999] W.-K. Leung, R. J. Needs, G. Rajagopal, S. Itoh, and S. Ihara, Calculations of Silicon Self-Interstitial Defects, Phys. Rev. Lett. 83, 2351 (1999).
  * McCluskey and Haller [2018] M. D. McCluskey and E. E. Haller,  _Dopants and Defects in Semiconductors 2nd Edition_ , 2nd ed. (CRC Press, 2018).
  * McCluskey and Janotti [2020] M. D. McCluskey and A. Janotti, Defects in Semiconductors, J. Appl. Phys. 127, 190401 (2020).
  * Cottom  _et al._ [2016] J. Cottom, G. Gruber, P. Hadley, M. Koch, G. Pobegen, T. Aichinger, and A. Shluger, Recombination centers in 4H-SiC investigated by electrically detected magnetic resonance and ab initio modeling, J. Appl. Phys. 119, 181507 (2016).
  * Alkauskas  _et al._ [2016] A. Alkauskas, M. D. McCluskey, and C. G. Van de Walle, Tutorial: Defects in semiconductors—Combining experiment and theory, J. Appl. Phys. 119, 181101 (2016).
  * Ža̧sinas and Vaitkus [2016] E. Ža̧sinas and J. V. Vaitkus, Disordered small defect clusters in silicon, Lith. J. Phys. 55, 330 (2016).
  * Huhtinen [2002] M. Huhtinen, Simulation of non-ionising energy loss and defect formation in silicon, Nucl. Instruments Methods Phys. Res. Sect. A Accel. Spectrometers, Detect. Assoc. Equip. 491, 194 (2002).
  * Vaitkus  _et al._ [2021] J. V. Vaitkus, A. Mekys, and Š. Vaitekonis, Electron mobility dependence on neutron irradiation fluence in heavily irradiated silicon, Lith. J. Phys. 61, 91 (2021).
  * Gaubas  _et al._ [2008] E. Gaubas, A. Kadys, A. Uleckas, and J. Vaitkus, Investigation of Carrier Recombination in Si Heavily Irradiated by Neutrons, Acta Phys. Pol. A 113, 829 (2008).
  * Gaubas  _et al._ [2018] E. Gaubas, T. Ceponis, L. Deveikis, D. Meskauskaite, J. Pavlov, V. Rumbauskas, J. Vaitkus, M. Moll, and F. Ravotti, Anneal induced transformations of defects in hadron irradiated Si wafers and Schottky diodes, Mater. Sci. Semicond. Process. 75, 157 (2018).
  * Breuer and Petruccione [2007] H. P. Breuer and F. Petruccione,  _Theory Open Quantum Syst._ (Oxford University Press, 2007).
  * Ashcroft and Mermin [1976] N. W. Ashcroft and N. D. Mermin,  _Solid state physics_ (Thomson Learning Inc., 1976).
  * Arora  _et al._ [1982] N. Arora, J. Hauser, and D. Roulston, Electron and hole mobilities in silicon as a function of concentration and temperature, IEEE Trans. Electron Devices 29, 292 (1982).


2021
MCM/ICM
Summary Sheet

Asian Giant Hornet Control based on Image Processing and Biological Dispersal

Changjie Lu, Shen Zheng, Hailu Qiu
College of Science and Technology, Wenzhou-Kean University

Summary

The Asian giant hornet (AGH) appeared in Washington State appears to have a
potential danger of bioinvasion. Washington State has collected public photos
and videos of detected insects for verification and further investigation. In
this paper, we analyze AGH using data analysis, statistics, discrete
mathematics, and deep learning techniques to process the data to control AGH
spreading.

First, we visualize the geographical distribution of insects in Washington
State. Then we investigate insect populations to varying months of the year
and different days of a month. Third, we employ wavelet analysis to examine
the periodic spread of AGH. Fourth, we apply ordinary differential equations
to examine AGH numbers at the different natural growth rate and reaction speed
and output the potential propagation coefficient. Next, we leverage cellular
automaton combined with the potential propagation coefficient to simulate the
geographical spread under changing potential propagation. To update the model,
we use delayed differential equations to simulate human intervention. We use
the time difference between detection time and submission time to determine
the unit of time to delay time. After that, we construct a lightweight CNN
called SqueezeNet and assess its classification performance. We then relate
several non-reference image quality metrics, including NIQE, image gradient,
entropy, contrast, and TOPSIS to judge the cause of misclassification.
Furthermore, we build a Random Forest classifier to identify positive and
negative samples based on image qualities only. We also display the feature
importance and conduct an error analysis. Besides, we present sensitivity
analysis to verify the robustness of our models. Finally, we show the
strengths and weaknesses of our model and derives the conclusions.

We find meaningful conclusions from our analysis. (1) Most Asian Giant Hornet
(AGH) is reported at West Washington from May to December. However, only 14
AGH proves to exist. (2) Local insects have a monthly active cycle of about
7.71 days. AGH breed at a rapid rate every Spring. Meanwhile, the geographical
spread of AGH heavily depends on the potential propagation coefficient. (3)
Earlier artificial intervention would efficiently control the AGH population
and protect other insects. (4) Light-weight convolutional neural network is
efficient and robust for classifying AGH. (5) High image Gradient Pictures
tend to contain AGH, while high image Contrast and background interference
will misguide classification.

We finally provide recommendations based on our analysis result. (1) AGH has a
tiny population and is nearly eradicated in Washington State. Therefore is no
need for massive panic. (2) AGH has a predicted spread period of approximately
one and a half months. When their growth rate is high, we should take earlier
actions and eliminate them before expansion. (3) We should prioritize image
gradient over image contrast and ignore background noise for subjective
classification. (4) We recommend a lightweight CNN for accurate
classification.

Key Words: SqueezeNet; Cellular Automaton; Wavelet Analysis; Delay
Differential Equation; Random Forest; Non-reference Image Quality Metrics

###### Contents

  1. 1 Introduction
    1. 1.1 Background
    2. 1.2 Related Works
    3. 1.3 Solutions
  2. 2 Assumptions and Notations
  3. 3 Data Preprocessing
    1. 3.1 Data Acquisition
    2. 3.2 Data Visualization
  4. 4 Expansion Model
    1. 4.1 Differential Equations
      1. 4.1.1 Ordinary Differential Equations
      2. 4.1.2 Delay Differential Equation
    2. 4.2 Cell Automaton
    3. 4.3 Wavelet Analysis
  5. 5 Classification Model
    1. 5.1 Data Preparation
    2. 5.2 Model Construction
    3. 5.3 Model Evaluation
      1. 5.3.1 Evaluation Metrics
      2. 5.3.2 Visual Demonstration
      3. 5.3.3 Quantitative Comparison
    4. 5.4 Image Evaluation
      1. 5.4.1 Image Quality Metrics
      2. 5.4.2 Model Implementation
      3. 5.4.3 Feature Importance & Error Analysis
  6. 6 Sensitivity Analysis
    1. 6.1 Expansion Model
    2. 6.2 Classification Model
  7. 7 Model Reflection
    1. 7.1 Strengths
    2. 7.2 Weaknesses
  8. 8 Memorandum

##  1 Introduction

###  1.1 Background

Asian Giant Hornet (Vespa mandarinia), a hornet species native to eastern
Asia’s temperate and tropical areas, is a predator to pollinating insects, the
world’s most enormous hornets.

On Vancouver Island in British Columbia, Canada, Sep 2019, an AGH colony was
detected and immediately destroyed. In May 2020, an AGH individual was
discovered in Washington state, which means AGH survived a winter in northern
America.

In the 21st century, biological invasions have become one of the most
significant problems to biodiversity, crop production, and ecological balance.
As one kind of social Hymenoptera, AGH perfectly fits in a spreading pattern
by human transportation and hidden skills from standard border inspection.
They hide in transported goods and can quickly establish a new colony with
only one fertilized queen. Furthermore, high fertility, comprehensive diet and
habitat, great defense, and competitive abilities also make them dominant in
the area[1]. A species under the same genus with AGH, which is called V.
velutina, invaded throughout the European continent with a rapid speed, and
their control costs can be very high[2]. A laboratory researched AGH spreading
pattern by building habitat suitability models and dispersal simulations,
concluded that AGH has a high possibility of spreading across western North
America with no human intervention[3].

To investigate the current situation of AGH and find the government of
Washington State offered a solution defending its bioinvasion, hotlines, and a
website for people to report the hornets’ emerges in different areas. The
public reports need to be processed through multiple analyses for further
study and giving solutions.

###  1.2 Related Works

Convolutional neural network (CNN) has been widely applied to image
classification tasks because of their excellent accuracy, capability, and
flexibility. Krizhevsky [4] introduces AlexNet, which uses max-pooling layers
and drop-out techniques. Simonyan [5], with others in VGG (Visual Geometry
Group), designed very deep convolutional networks for large-scale image
recognition, which grows deeper layers with a small kernel size of 3 * 3.
Lightweight image classification models have been designed to fit mobile
devices that have limited memory capacity. MobileNet [6] uses depth-wise
separable convolutions to a trade-off between speed and accuracy efficiently.
SqueezeNet [7] adopts small filters, late downsampling, and fire block
(squeeze and expand operations) for flexible and efficient model compression.

###  1.3 Solutions

In this paper, we propose a novel approach for analyzing, predicting, and
giving recommendations towards Asia Hornet. Our analysis uses different
discrete mathematics, deep learning, and data visualization tools. First, we
import and preprocess the data, including mapping images, videos to text, and
numbers. Second, we use dispersal models. Our dispersal model mainly includes
different equations, cell automation (CA), and wavelet analysis. Third, we
bring a CNN for classifying AGH. We leverage statistical evaluation metrics
and activation maps to verify our model performance. Fourth, we explore the
cause for misclassification using a random forest [8] classifier with several
non-reference image quality metrics.

Figure 1: Paper Workflow

The rest of the paper is organized as below. Section 2 introduces basic
notations and assumptions. Data acquisition and data visualization are
conducted in Section 3. Section 4 and section 5 construct the expansion model
and the classification model, respectively. Section 6 perform sensitivity
analysis and Section 7 analyzes the strengths and weaknesses. Conclusion and
recommendations will be addressed in the Memorandum.

##  2 Assumptions and Notations

| Abbreviation | Full form |
---|---|---|---
| AGH | Asian Giant Hornet |
| CNN | Convolutional Neural Network |
| CA | Cellular Automaton |
| RF | Random Forest |
| ODE | Ordinary Differential Equation |
| DDE | Delay Differential Equation |
| TOPSIS | Technique for Order of Preference by Similarity to Ideal Solution |
| NIQE | Naturalness Image Quality Evaluator |
| ReLU | Rectified Linear Unit |
| Conv2d | 2D Convolution Layer |
| OOB | Out Of Bag |
Table 1: Abbreviation of Relevant Terms | Symbol | Explanation |
---|---|---|---
| A(t) | Amount of AGH at any given time |
| B(t) | Amount of local species at any given time |
| a | Natural growth rate of AGH |
| b | Internal competition coefficient of AGH |
| c | Aggression coefficient of AGH to local species |
| d | Natural growth rate of local species |
| e | Internal competition coefficient of local species |
| f | Resistance coefficient of local species to AGH |
| g | Effects of manual intervention on AGH |
| h | Effects of manual intervention on local species |
Table 2: Symbols and Explanations

  * •

We assume that AGH only competes with local species and receive mininum impact
from altitude, weather, and human activities.

  * •

We assume that the map of Washington State can divides into 128*128 squares.

  * •

We assume that the spread pattern of AGH in Washington State is relatively
fixed.

  * •

We assume that the proportional relationship between AGH and other dwelling
species is similar to the number of reports submitted within a specific range
of latitude and longitude.

##  3 Data Preprocessing

###  3.1 Data Acquisition

We have several data files for analysis. The first is a pdf file named
Vespamandarinia, which we have used in the introduction section. The second is
a pair of a spreadsheet. We use GlobalID to combine two spreadsheets. We first
convert all images into jpg format. We use the strategy for word and pdf file
that contains only images. We then take a screenshot of the video and convert
it to jpg images. We take the screenshot at the first explicit appearance of
insects.

###  3.2 Data Visualization

Figure 2: Geographical Visualization

We perform classification using lab comments and lab status. We count each
insect species’ number under the same lab status in different years by
detection date, from January to December, integrating them into one
spreadsheet. The Geographical Visualization result in Fig.1 shows most
reported cases appear in west Washington. Besides, all positive ID appears in
northwest Washington.

Figure 3: Change of Insect population over Time. X-axis: the month where a
insect is observed. Y-axis: Corrected log value of insect amount in a specific
month. The category positive, negative, unprocessed, and unverified from lab
status data means AGH species insects, non-AGH species insects, unprocessed
insects, and unverified insects respectively.

We have several findings from Fig.3. First, each category begins a rapid
growth in May and peaked in around September. Second, the negative category
has the highest numbers. Third, the positive and unprocessed category has the
least numbers with the similar trend, where as the unverified category has a
similar trend with the negative category.

In Fig.4(b), we have categories including symphyta (sawfly), wasp, Sphex
ichneumoneus (great golden sand digger), Bombus (bumblebee) and Apis (honey
bee), Dolichovespula maculata (bald-head hornet) and Vespa crabro (European
hornet), and other species.

Figure 4: (a): Change of AGH population over time. X-axis: the year and month
when a AGH is observed. Y-axis: Numbers of AGH at a specific time. (b): Change
of Non-AGH population over time. X-axis: the month where a insect is observed.
Y-axis: the numbers of reported insects that is not AGH. We roughly classified
insects into 5 categories and others.

We have several findings from Fig.4. First, there are 14 AGHs detected on
specific dates. Second, no AGHs is detected during 2020-01 to 2020-05. Third,
only one AGD is reported in one detection. Besides, all non-AGH insects have a
rapid increase in June and a signficant decline in September. Last, Sphex
ichneumoneus and wasp have peak values of around 400 during July and August,
whereas others category of insects does not.

Figure 5: Total Insect Numbers Detected by Day in a Month. X-axis: day in a
month, from 1 to 31. Y-axis: lab status comment for indicating insect
categories. Legend: different lab status Figure 6: Insect Periodic Activity
with Wavelet Analysis. Left X-axis: time in days since simulation. Left
Y-axis: cycle time period. Right X-axis: power spectrum frequency. Right
Y-axis: intensity

We introduce wavelet analysis [9] in Fig.6 to investigate if the insects have
monthly activity patterns. Wavelet analysis figure show that inserts have a
7.71-day period, which explains Washington State insects’ monthly activity as
7 days.

##  4 Expansion Model

###  4.1 Differential Equations

####  4.1.1 Ordinary Differential Equations

We use systems of ordinary differential equations (ODE) to simulate how insect
populations evolve. Our different equations seek equilibrium points given
initial values and constraints. In Fig.7, we describe AGH as competing species
to local species.

We set up a differential equation to study the population change of AGH in
Washington. The system of differential equations established is as follows:

$\left\\{\begin{array}[]{l}\frac{dA(t)}{dt}=aA(t)-bA^{2}(t)-cA(t)B(t)\\\ \\\
\frac{dB(t)}{dt}=dB(t)-eB^{2}(t)-fA(t)B(t)\end{array}\right.$ (1)

(a) a=0.5

(b) a=0.55

Figure 7: The solution of our ODE systems. (a): Solution when Natural Growth
Rate of Asian Giant Hornet (AGH) is 0.5. (b): Solution when growth rate is
0.5. For both conditions, we set b = 0.3, c = 0.2, d = 0.4, e = 0.1, and f =
0.8

The number of AGH found in northwest Washington is roughly 1:10 for the number
of Local species. Therefore, we set A = 0.4, B = 4 as a starting point. The
result of our ODE appears in Fig.7. When the natural growth rate of AGH is
0.5, the amount of AGH decreases steadily from 0.5 to 0, whereas the amount of
local species decreases sharply to 2.5 before slowly increasing to 4.0. When
the natural growth rate of AGH is 0.55, the amount of AGH increases to 2,
whereas the number of local species decreases to 0. When time is greater than
25, they arrive at equilibrium.

####  4.1.2 Delay Differential Equation

We use systems of delay differential equations (DDE) to analyze the problem
further. Unlike ODE, DDE has the derivative of the unknown function in terms
of previous time points. The formula of our delayed differential equation is
below. Note the last terms at each line comes from an early reaction. Others
are the same as the ODE.

We calculated the average difference between image submission time and
detection time. The parameters g and h represent the influence of the manual
intervention on the population. (g = 0.05, h = 0.001) Delays in manual
intervention can lead to dramatic changes in population evolution.

$\left\\{\begin{array}[]{l}\frac{dA(t)}{dt}=aA(t)-bA^{2}(t)-cA(t)B(t)-gA(t-10)\\\
\\\ \frac{dB(t)}{dt}=dB(t)-eB^{2}(t)-fA(t)B(t)-hB(t-10)\end{array}\right.$ (2)

Figure 8: Solution Comparison between ODE and DDE. We set the natural growth
rate of AGH as 0.55. We use red for AGH and blur for other species.

Fig.8 shows the result from earlier and on-time actions towards AGH. In the
ODE solution, there is quick increase of AGH and a sharp decrease of other
species. In the DDE solution, there is quick increase of other species and
steadily decline of AGH. Our result demonstrate earlier action help control
AGH numbers.

###  4.2 Cell Automaton

Cellular Automaton (CA) is a discrete computational model in automation
theory. It consists of finite-dimensional cells with states (e.g.,
True/False). A new generation is created according to specific algorithms to
iterative update the neighboring cells.

From the the ODE, We get values of A(t) and B(t) with respect to time.
Therefore, we define the potential propagation coefficients:

$\alpha(t)=\frac{A(t)}{A(t)+B(t)}$ (3)

We apply CA to simulate the spread of AGH in Washington and simplify the map
into a square, dividing into 128*128 cells. In the upperleft 15*15 grid, we
use red to represent AGH existence and black to represent native species.
Since most AGH cases is reported in northwest Washington, we set initial
probability ratio as 4:1. The probability of each grid occupied by AGH is (if
P ¡ 1):

$P(\text{cells
}=1)=\frac{\operatorname{sum}}{2}\cdot\alpha\left(\frac{t}{10}\right)$ (4)

Where sum(x,y) represents AGH’s amount around the grid with coordinate
position (x,y), and t is the time unit of CA. After every 50 units of time,
there will be 10 units where AGH cannot reproduce and has a death probability
of 0.3:

$\text{ P(cells }=0)=0.3$ (5)

(a) t=1

(b) t=11

(c) t=21

(d) t=50

(e) t=64

Figure 9: Cellular Automaton Using the First Potential Propagation Coefficient

(a) t=1

(b) t=12

(c) t=67

(d) t=85

(e) t=153

(f) t=163

(g) t=241

Figure 10: Cellular Automaton Using the Second Potential Propagation

In Fig.9, we visualize AGH’s evolution based on the first differential
equation, from t = 1 to 64, AGH gradually decrease and finally disappear.

In Fig.10, we visualize AGH’s evolution based on the first differential
equation, from t = 1 to 64, AGH gradually decrease and finally disappear. , we
visualize AGH’s evolution based on the second differential equation. From t =
1 to 67, AGH has small and fluctuating numbers. Starting from t = 85, the AGH
population grow explosively and finally occupy Washington.

###  4.3 Wavelet Analysis

Wavelet transform is a decomposition of a continuous signal into different
scale components. Wavelet analysis is better than Fourier analysis in that it
can capture the time information of frequency change. Therefore, wavelet
analysis is suitable for complicated movements such as the spread of AGH. Our
Morlet wavelet is a continuous plane wave modulated by Gauss function:

$\Psi(t)=e^{icl}\left(\mathrm{e}^{-\frac{t^{2}}{2}}-\sqrt{2}\mathrm{e}^{-\frac{c^{2}}{4}}\mathrm{e}^{-l^{2}}\right)$
(6)

Its discrete mathematical expression is:

$W_{f}(a,b)=|a|^{-\frac{1}{2}}\sum_{i=1}^{N}f(i\delta
t)\psi^{\star}\left(\frac{i\delta t-b}{a}\right)$ (7)

where * represents complex conjugate; a is scale factor (related to period and
frequency); b is the translation factor (time position); I is time and
position label of data series; f(t) is the variable time series; W is the
wavelet coefficient; $\delta t$ is the time interval of the variable series.
Since the wavelet is in complex from, the coefficients after the wavelet
transform are also complex.

Figure 11: Red Grid Numbers in Cellular Automaton over Time. X-axis: Time
since starting simulation. Y-axis: Red Grid amount. We select the condition of
the second differential equation for plotting. More Grid counts represents
more AGH occupancy.

Fig.11 shows a general increasing trend for AGH numbers. However, there is a
sharp decrease for after every 50 time units. Therefore, eliminating the Asian
bumblebee in winter is a better option.

Figure 12: Wavelet Analysis of Red Grids Numbers. Left X-axis: Time in days
since simulation. Left Y-axis: Cycle Time Period. Right X-axis: Power
Spectrum. Right Y-axis: Intensity

Fig.12 shows that the amount of AGH in CA has a period of 46.75 time unit,
which meets the stagnation phase of 50 days. Besides, the only difference
between two ODE is the reproduction ability of AGH.

The main factors affecting the reproductive ability of AGH are temperature,
humidity, and manual intervention. AGH has high fertility, and a little
deviation may cause the amount of the AGH has divergent results. Moreover, AGH
is intolerant to frigidity; First, AGH almost disappear In the 67th unit which
is in winter. However, it has fast outbreak afterwards. Therefore, the
Washington Department of Agriculture should pay attention to the number of AGH
at next spring. Second, 3 unidentified samples were classified as AGH at
spring or summer on the west part of Washington. Therefore, we should monitor
the potentially spreading at west Washington.

##  5 Classification Model

###  5.1 Data Preparation

Figure 13: Batch of Cropped Dataset

We use an AGH dataset from a Github repository [10] to make 509 positive
samples. We then select 569 negative examples from the given dataset (See
Fig.13 for a sample batch). The total amount of pictures for the model
training and testing is 1078. We first resize the cropped dataset to the shape
of 224*224. Then we add Random Horizontal Flip to flip with a probability of
0.5 horizontally. Next, we standardize the data using ImageNet’s [11] mean and
standard deviation. After that, we split the data into training, validation,
and testing set with the ratio of 0.70, 0.15, and 0.15, respectively. We have
757 images for training, 160 for validation, and 161 for testing. Ultimately,
we prepare training and testing labels, labeling “1” for positive (AGH) and
“0” for negative (non-AGH).

###  5.2 Model Construction

We choose a light-weight model called SqueezeNet [7]. Since our dataset is
small ( 1000 images), using a light-weight CNN can prevent overfitting and
achieve better accuracy. We use pretrained weight from ImageNet. We then
change the final classification block output from 1000 to 2 for binary
classification. We use Dropout [12] with 0.5 probability to prevent
overfitting. We also use ReLU [13] activation layers to mitigate vanishing
gradient problem.

Figure 14: Model Architecture of SqueezeNet. SqueezeNet consists of features
layers and classifier layers. We use Conv2d, ReLU, MaxPool2d, Fire, Dropout
and AvgPool2d for convolutional layer, ReLU activation layer, maximum pooling
layer, Fire Block, Dropout layer and average pooling layer, respectively.
Figure 15: Geographical Visualization after SqueezeNet classification. We use
red for positive samples and green for negative samples

Fig.15 show most positive samples (AGH) are found at Northwest Washington.
Only 1 case is found in Southwest Washington. We use the Pytorch framework for
training our model. We set the epoch as 3 and batch size as 8. We adopt cross-
entropy loss as the cost function. We utilize stochastic gradient descent with
a learning rate of 1e-3 and momentum of 0.9. Since our training only consists
of 3 epochs, there is no need for learning rate decay.

(a) Loss (b) Accuracy

Figure 16: Loss-Accuracy Plot. (a): Training and Validation Loss. (b):
Training and Validation Accuracy.

Fig.16 shows the loss and accuracy during training and validation. Our
training loss quickly converges and remains stable between 0 to 0.03. Our
validation loss has a spike at the iteration of 30 but has a general trend of
declining from 0.3 to 0.1. Except for the sharp decline in the middle, our
training and validation accuracy quickly increase and stay between 0.95 and
1.00. Since both accuracies are similar at a high level, there is neither
overfitting nor underfitting. We obtain the final testing accuracy of 0.975.

###  5.3 Model Evaluation

####  5.3.1 Evaluation Metrics

We use a group of evaluation metrics to evaluate our models. Specifically, we
include loss, total accuracy, class accuracy, recall, specificity, precision,
f1 score, and confusion metric for evaluating our image classification model.

Confusion matrix displays a table of true positive (TP), False Positive (FP),
False Negative (FN) and True Negative (TN). The formula for other evaluation
metrics is given as below.

$\mathrm{Accuracy}=\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}}$
(8) $\operatorname{Recall}=\frac{TP}{TP+FN}$ (9) $\text{ Specificity
}=\frac{TN}{TN+FP}$ (10) $\text{ Precision }=\frac{TP}{TP+FP}$ (11) $F1\text{
Score }=\frac{2*\text{Precision}*\text{Recall}}{\text{ Precision }+\text{
Recall }}$ (12)

####  5.3.2 Visual Demonstration

(a) Pos1

(b) Pos2

(c) Pos3

(d) Pos4

(e) Act1

(f) Act2

(g) Act3

(h) Act4

Figure 17: Activation Layers Result with samples classified as positive (AGH).
The upper row is the input image. The lower row is the output from the first
convolution layer. We use ”Pos” for positive samples and ”Act“ for activation
layers, respectively.

(a) Neg1

(b) Neg2

(c) Neg3

(d) Neg4

(e) Act1

(f) Act2

(g) Act3

(h) Act4

Figure 18: Activation Layers Result with samples classified as Negative (non-
AGH). The upper row is the input image. The lower row is the output from the
first convolution layer. We use ”Neg” for positive samples and ”Act“ for
activation layers, respectively.

We visualize the first convolutional layer result in Fig.17 and Fig.18. Our
rows display the positive sample, positive sample activation result, negative
sample, and negative sample activation result, respectively. As shown in the
figure, our neural network can filter background interference and separate the
given insects’ pattern. Besides, our CNN displays more evident edges for
positive samples. That means our model can efficiently capture the distinct
feature of the AGH.

####  5.3.3 Quantitative Comparison

In this part, we quantitatively analyze our model classification performance.
We first display the confusion matrix. According to Fig.19, we have 82 true
positives, 75 true negatives, 3 false positives, and 1 false negative. We also
note from Table 3 that all evaluation metrics have a score greater than 0.96.
That means our model has great classification accuracy. Besides, our model has
better predictive accuracy at positive samples than negative samples. That is
reasonable because the negative examples contain numbers of insects that look
similar to AGH (e.g., European Hornet). Table 5 shows the image quality.
Positive samples win all evaluation metrics except image contrast. That means
image contrast is a misleading metric for classifying AGH.

| Accuracy | Class Accuracy | Recall | Specificity | Precision | F1 Score
---|---|---|---|---|---|---
Value | 0.975 | 0.965/0.987 | 0.988 | 0.962 | 0.965 | 0.976
Table 3: Quantitative result of our image classification model. Class Accuracy
has negative (not AGH) on the left and positive (AGH) on the right. Figure 19:
Confusion Matrix of our image classification model. Left: Actual Label.
Bottom: Predicted Label. from upperleft to lowerright cell: True Positive,
False Positive, False Negative and True Negative,

###  5.4 Image Evaluation

####  5.4.1 Image Quality Metrics

We use a set of non-reference image quality metrics to understand why people
would misclassify other insects as Asian hornets. That group of image quality
metrics would also help investigators screen out a more likely picture to
contain AGH. Our image quality metrics include Naturalness Image Quality
Evaluator (NIQE) [14] , image gradient, entropy, and contrast. We also use
TOPSIS to summarize the indicators above. NIQE is a statistical measure
derives from natural, undistorted images. Image Gradient is a vectorized
change in image intensity. Image Entropy is a statistical measure of
randomness for assessing the texture of an image. Image Contrast is the
difference in luminance that makes a picture visually distinguishable. TOPSIS
is a multi-criteria decision analysis method that uses weighted and normalized
geometric distance for assessing image quality. We use TOPSIS to give a
comprehensive evaluation of image quality.
NIQE:

$\displaystyle D\left(\nu_{1},\nu_{2},\right.$
$\displaystyle\left.\Sigma_{1},\Sigma_{2}\right)$ (13)
$\displaystyle=\sqrt{\left(\left(\nu_{1}-\nu_{2}\right)^{T}\left(\frac{\Sigma_{1}+\Sigma_{2}}{2}\right)^{-1}\left(\nu_{1}-\nu_{2}\right)\right)}$

Image Gradient:

$\nabla f=\left[\begin{array}[]{l}g_{x}\\\
g_{y}\end{array}\right]=\left[\begin{array}[]{l}\frac{\partial f}{\partial
x}\\\ \frac{\partial f}{\partial y}\end{array}\right]$ (14)

Entropy:

$S=k_{b}\ln\Omega$ (15)

Contrast:

$\sum_{i=0}^{N-1}\sum_{j=0}^{M-1}\left(I_{ij}-\bar{I}\right)^{2}$ (16)

We acquire the positive and negative ideal points of each indicator in TOPSIS:

$d_{i}^{*}=\sqrt{\sum_{j=1}^{4}\left(z_{ij}-z_{j}^{*}\right)^{2}},d_{i}^{-}=\sqrt{\sum_{j=1}^{4}\left(z_{ij}-z_{j}^{-}\right)^{2}}$
(17)

where Z is the weighted evaluation matrix.

| NIQE | Gradient | Entropy | Contrast | TOPSIS
---|---|---|---|---|---
Negative | 6.98 | 0.057 | 7.385 | 3092 | 0.1912
Positive | 10.10 | 0.075 | 7.525 | 2838 | 0.2655
Table 4: Image Quality Comparison for positive and negative samples. We use
bold to indicate higher values.

####  5.4.2 Model Implementation

We use four non-reference image quality metrics values as features and the
other one as the label. We then split the training and testing data with a
train/test ratio of 2:1 and used random forest classifiers to classify
positive and negative samples.

Random Forest classifier is an ensembled supervised machine learning
algorithm. It is a meta estimator that leverages numbers of decision trees and
takes the average of those predictors to improve forecasting accuracy and
reduce predictability variability.

| n_estimators | max_depth | min_samples_split | min_samples_leaf
---|---|---|---|---
Value | 644 | 110 | 10 | 2
Table 5: Best Hyperparameters of the random forest classifier. n_estimators:
numbers of decision tree in random forest. max_depth: maximum depth of the
tree.min_samples_split: minimum samples for splitting internal nodes.
min_samples_leaf: minimum samples for a leaf node.

We initially set the min sample split as two and the numbers of estimators as
150. We then select four hyperparameters for random grid search cross-
validation. Specifically, We run 3-fold cross-validation for 100 iterations
and find the hyperparameter group that minimizes the loss function. The best
hyperparameters summarize in the Table 5. Our Testing Accuracy from the
optimized hyperparameters is approximately 0.686.

####  5.4.3 Feature Importance & Error Analysis

Figure 20: (a): Feature Importance Score. X-axis: feature index. From left to
right, NIQE, gradient, entropy and contrast. Y-axis: Feature Importance Score,
from 0 to 1. (b): Error Analysis. X-axis: number of estimators. Y-axis: out-
of-bag (OOB) error rate

We visualize the feature importance in the Fig.20(a). The feature importance
score describes how significant a specific feature is for image
classification. We observe from the plot that the image gradient has dominant
importance (around 0.40) for classification, whereas other factors contribute
similarly.

We also examine the source of classification errors in Fig.20(b). In Fig.20
(b), we note different types of maximum features: the number of features for
calculating the best split. When max features are ’sqrt’ (the blue line), the
classifier will consider maximum features as the total features’ square root.
When max features are ’None’ (the yellow line), it equates maximum features to
total features.

For both max features types, we have decreasing OOB error rate with increasing
numbers of features. However, we note that when numbers of estimators grow
greater than 600, the OOB error rate remains stable. Therefore, it is
efficient to stop at around 600.

##  6 Sensitivity Analysis

###  6.1 Expansion Model

(a) 0.475¡a¡0.525

(b) 0.38¡c¡0.42

(c) 0.19¡b¡0.21

(d) 0.76¡e¡0.84

Figure 21: Expansion Model Response to Input Variation

Fig.21 shows how our expansion model responds to input variation. For AGH, Our
model has reasonable stability under input variation. For other insects, The
internal competition coefficient of AGH results in the most significant
fluctuation, whereas the internal competition coefficient of local species
leads to the least fluctuation.

###  6.2 Classification Model

(a) Ground Truth

(b) Averaging Blur

(c) Gaussian Blur

(d) Median Blur

Figure 22: Clean and Blurred Images. All blur has kernel size of 5.

In this subsection, we conduct sensitivity analysis to the classification
model. We first apply disruptive factors, including averaging, Gaussian Blur,
and Median Blur, to the groudtruth. The resulting images are in Fig.22.

(a) Average Blur (b) Gaussian Blur (c) Median Blur

Figure 23: Confusion matrix under different blur conditions

In Fig.23, We display the confusion matrix under different blur conditions.
Under all blur conditions, our model has 0 or only 1 false negatives.
Meanwhile, it has some false positives, but its number is significantly less
than true positives and true negatives. That means our model can recognize
Asia Wasps under disturbance.

Figure 24: Radar Plot under different blur conditions

We then build a radar plot in Fig.24 to assess model performance under
disturbance. All blurring has no significant influence on positive accuracy
and recall. However, Averaging Blur and Gaussian Blur will slightly decrease
direct accuracy, specificity, and f1 score and significantly reduce negative
accuracy and precision. Meanwhile, median blur results in less degradation.

##  7 Model Reflection

###  7.1 Strengths

  * •

Consistent classification and predictive ability.

  * •

Quick and Convenient Model Update from additional data.

  * •

High Simulation Accuracy with seasonal periodic change, Brownian motion and
potential propagation coefficient

  * •

Efficient Analysis of complex images and videos

  * •

Great classification performance and robustness

  * •

Novel binary classification from pure image quality.

###  7.2 Weaknesses

  * •

Idealize inter-species competition and override human intervention factors.

  * •

Possibly inconsistent potential propagation coefficient update in cellular
automation

  * •

Over-Simplified factors others than dual species competition between AGH and
Bees.

  * •

Sub-optimal interpretability at intermediate layer’s feature processing

  * •

Subjective Image Indicator for binary classification.

##  8 Memorandum

To: The Washington State Department of Agriculture

From: Team #2116473 of 2021 MCM

Date: Feb 9, 2021

Subject: Recommendations and Solutions to Potential Bioinvasion of Asian Giant
Hornet

Honorable Washington State Department of Agriculture,

Recently, Asian giant hornet was discovered around Washington State. According
to research, this species of hornet could be extremely dangerous and a serious
threat that will probably lead to a biological invasion in a short time. To
predict the pest’s spread and control it in the future, we make a serial of
analysis upon current data.

The report has analyzed the Asian giant hornet using statistics, mathematics,
and deep learning tools. First, we visualize the geographical distribution of
insects, investigate insect populations of different dates and examine the
their expansion period. Second, we calculate Asian giant hornet numbers at
different conditions and simulate its geographical spread. Third, we build an
advanced classification model and discuss the reason for widespread mistaken
reporting. Besides, we demonstrate the robustness of our models under various
input variations. Finally, we recommend measures for successfully addressing
Asian giant hornet.

There are insightful findings from our analysis.

  * •

Most Asian giant hornet is reported at west part of Washington. We use a deep
learning model with 0.975 testing accuracy for verification. However, we only
found 14 Asian Giant Hornet among 4400+ reports.

  * •

Local insects have a monthly active cycle of about 7.71 days. Asian giant
hornet breed at a rapid rate every Spring.

  * •

Asian giant hornet’s explosion rely on high potential propagation coeffcient,
according to our geographical spread analysis. However, AGH are inactive all
year along and can hardly reach that condition.

  * •

Asian giant hornet’s fragile at early explosion. Timely human interference
(about 30 days earlier than spreading) could efficiently control the AGH
population and protect other insects.

  * •

Light-weight deep learning models is far more accurate (¿0.97) for classifying
Asian giant hornet than human observations (¡0.01).

  * •

It is most efficient to update our models (with pretrained weight saved) after
every winter.

  * •

Reported image with high intensity change is a good metric for subjectively
classifying AGH, whereas images with high luminary is misleading.

We finally provide recommendations based on our analysis result.

  * •

Clear Asian giant hornet hives at early spring. You may use fire fighting
drones to accelerate that process.

  * •

Timely Classify and Report Asian giant hornet using light-weight convolutional
neural network. Once you detect the signs of an outbreak, you must take action
within 30 days.

  * •

Ignore reports that has low image intensity change or has noisy backgrounds.
They are almost certainly fake Asian giant hornet images. Be careful of
misguiding image with high contrast.

  * •

Update classification and explosion models at every winter. At that time, you
should have collected data till the end of a year.

  * •

There is no need for massive panic about Asian giant hornet. Even among large
numbers eyewitnesses reports over a long period of time, there are few AGH.

  * •

The public should be informed of AGH to reduce unnecessary review time.
Meanwhile, you could prioritize images with high TOPSIS score.

We suggest exploiting updating our model in the following ways:

First, you can use the non-reference image quality metrics and Topsis score to
determine the order of image audit. Then you can use SqueezeNet to identify if
a image contains an Asian giant hornet. For images that are difficult to
discriminate, you can use manual screening.

After screening the images, you should collect the data on the adaptation of
Asian Giant Hornets in Washington to update the differential equations for
simulation. You should also pay attention to the number of Asian giant hornet,
at early Spring. Finally, the government should take prompt action when there
is a evidence of outbreak.

We hope our analysis and advice will help to control spreading of Asian giant
hornets and generating better plans for future issues. You are welcome to
contact us at any time.

Yours Sincerely,

Team #2116473 of 2021 MCM

## References

  * [1] Jacqueline R Beggs, Eckehard G Brockerhoff, Juan C Corley, Marc Kenis, Maité Masciocchi, Franck Muller, Quentin Rome, and Claire Villemant.  Ecological effects and management of invasive alien vespidae.  BioControl, 56(4):505–526, 2011.
  * [2] Alberto J Alaniz, Mario A Carvajal, and Pablo M Vergara.  Giants are coming? predicting the potential spread and impacts of the giant asian hornet (vespa mandarinia, hymenoptera: Vespidae) in the usa.  Pest Management Science, 77(1):104–112, 2021.
  * [3] Gengping Zhu, Javier Gutierrez Illan, Chris Looney, and David W Crowder.  Assessing the ecological niche and invasion potential of the asian giant hornet.  Proceedings of the National Academy of Sciences, 117(40):24646–24648, 2020.
  * [4] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  Advances in neural information processing systems, 25:1097–1105, 2012.
  * [5] Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  arXiv preprint arXiv:1409.1556, 2014.
  * [6] Andrew G Howard, Menglong Zhu, Bo Chen, Dmitry Kalenichenko, Weijun Wang, Tobias Weyand, Marco Andreetto, and Hartwig Adam.  Mobilenets: Efficient convolutional neural networks for mobile vision applications.  arXiv preprint arXiv:1704.04861, 2017.
  * [7] Forrest N Iandola, Song Han, Matthew W Moskewicz, Khalid Ashraf, William J Dally, and Kurt Keutzer.  Squeezenet: Alexnet-level accuracy with 50x fewer parameters and¡ 0.5 mb model size.  arXiv preprint arXiv:1602.07360, 2016.
  * [8] Leo Breiman.  Random forests.  Machine learning, 45(1):5–32, 2001.
  * [9] Christopher Torrence and Gilbert P Compo.  A practical guide to wavelet analysis.  Bulletin of the American Meteorological society, 79(1):61–78, 1998\.
  * [10] Albert.  no.  https://github.com/Albert-l-w/no, 2021.
  * [11] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In 2009 IEEE conference on computer vision and pattern recognition, pages 248–255. Ieee, 2009.
  * [12] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: a simple way to prevent neural networks from overfitting.  The journal of machine learning research, 15(1):1929–1958, 2014\.
  * [13] Vinod Nair and Geoffrey E Hinton.  Rectified linear units improve restricted boltzmann machines.  In Icml, 2010.
  * [14] Anish Mittal, Rajiv Soundararajan, and Alan C Bovik.  Making a “completely blind” image quality analyzer.  IEEE Signal processing letters, 20(3):209–212, 2012.