# Generalized anisotropic $\kappa$-cookbook: 2D fitting of Ulysses electron
data

K. Scherer,1,2 and E. Husidic1,3, M. Lazar,1,3, H. Fichtner,1,2
1Institut für Theoretische Physik, Lehrstuhl IV: Plasma-Astroteilchenphysik,
Ruhr-Universität Bochum, D-44780 Bochum, Germany
2Research Department, Plasmas with Complex Interactions, Ruhr-Universität
Bochum, 44780 Bochum, Germany
3Centre for Mathematical Plasma Astrophysics, 3001 Leuven Belgium
<EMAIL_ADDRESS>
###### Abstract

Observations in space plasmas reveal particle velocity distributions out of
thermal equilibrium, with anisotropies (e.g., parallel drifts or/and different
temperatures, $T_{\parallel}$ \- parallel and $T_{\perp}$ \- perpendicular,
with respect to the background magnetic field), and multiple quasithermal and
suprathermal populations with different properties. The recently introduced
(isotropic) $\kappa$-cookbook is generalized in the present paper to cover all
these cases of anisotropic and multi-component distributions reported by the
observations. We derive general analytical expressions for the velocity
moments and show that the common (bi-)Maxwellian and
(bi-)$\kappa-$distributions are obtained as limiting cases of the generalized
anisotropic $\kappa$-cookbook (or recipes). Based on this generalization, a
new 2D fitting procedure is introduced, with an improved level of confidence
compared to the 1D fitting methods widely used to quantify the main properties
of the observed distributions. The nonlinear least-squares fit is applied to
electron data sets measured by the Ulysses spacecraft confirming the existence
of three different populations, a quasithermal core and two suprathermal (halo
and strahl) components. In general, the best overall fit is given by the sum
of a Maxwellian and two generalized $\kappa$-distributions.

###### keywords:

plasmas, Sun: heliosphere, solar wind, methods: data analysis

††pagerange: Generalized anisotropic $\kappa$-cookbook: 2D fitting of Ulysses
electron data–B.2

##  1 Introduction

Space plasmas such as the solar wind or planetary magnetospheres are dilute
and (nearly) collisionless (Marsch, 2006; Heikkila, 2011), and particle
velocity distributions measured in-situ reveal non-equilibrium states, both
for electrons (Štverák et al., 2008; Wilson et al., 2019a, 2020) and ions
(Gloeckler, 2003; Kasper et al., 2006; Lazar et al., 2012b). The distributions
exhibit typical non-Maxwellian features such as suprathermal tails due to an
excess of high-energy particles (Maksimovic et al., 1997; Štverák et al.,
2008), or deviations from isotropy, e.g., temperature anisotropies with
respect to the magnetic field direction (Marsch, 2006; Štverák et al., 2008),
or field-aligned beams (e.g., electron strahls) along the magnetic field lines
with speeds higher than the bulk plasma flow (Pierrard et al., 2001; Wilson et
al., 2019b).

Up to a few keV the electron velocity distributions in space plasmas reveal
three distinct components (Pierrard et al., 2001; Wilson et al., 2019b): at
low energies a quasi-thermal and highly dense core ($\approx$ 80% - 90% of the
total particle number density) approximated by a Maxwellian; and two
suprathermal populations, a hot, tenuous halo ($\approx$ 5% - 10% of the total
density) enhancing suprathermal tails, which are well described by power-law
functions such as the $\kappa$-distribution function (see below), and likewise
hot and tenuous beams ($\approx$ 5% - 10% of the total density), which can
also be fitted by $\kappa$-distributions, see e.g., Štverák et al. (2008);
Wilson et al. (2019a, 2020). Electrons in the solar wind are described by
different combinations of multi-component models, i.e., dual models to fit
only the core and halo when the strahl is not prominent enough, triple models
to include also the asymmetric strahl, or even a quadruple model to include
two counterbeaming strahls (double strahl) observed in closed magnetic field
topologies, associated to coronal loops or interplaneteraty shocks (Pilipp et
al., 1987b; Lazar et al., 2012a, 2014; Wilson et al., 2019b; Macneil et al.,
2020).

Kappa (or $\kappa$-) power laws have become a widely employed tool to describe
the observed distributions in non-equilibrium distributions, either as a
global fit to incorporate both the core and halo (Olbert, 1968; Vasyliunas,
1968), or as a partial fit to reproduce only the halo or strahl components,
see the review by Pierrard & Lazar (2010) and references therein. In the
present work we propose a generalized fitting model capable to reproduce any
of the aforementioned anisotropies or combinations of beams (or drifting
populations) with intrinsic anisotropic temperatures, as revealed by the
realistic distributions in space plasmas (Maksimovic et al., 2005; Štverák et
al., 2008; Pilipp et al., 1987a, b, c). In kinetic studies, e.g., when
exploring the stability properties of these distributions, it is important to
capture their non-Maxwellian features, and thus quantify the free energy
sources triggering instabilities and wave fluctuations (e.g. Astfalk & Jenko,
2016). In the absence of collisions space plasmas are governed by the wave-
particle interactions.

The anisotropic model mostly invoked is the bi-$\kappa$-distribution function
(BK) (Summers & Thorne, 1992; Lazar & Poedts, 2009)

$\displaystyle f_{\mathrm{BK}}(\vec{r},\vec{v},t)=$
$\displaystyle\frac{n(\vec{r},t)}{\mathchoice{{\hbox{$\displaystyle\sqrt{\pi^{3}\,\kappa\,}$}\lower
0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\textstyle\sqrt{\pi^{3}\,\kappa\,}$}\lower
0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\scriptstyle\sqrt{\pi^{3}\,\kappa\,}$}\lower
0.4pt\hbox{\vrule
height=4.30276pt,depth=-3.44223pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\pi^{3}\,\kappa\,}$}\lower
0.4pt\hbox{\vrule
height=3.44165pt,depth=-2.75334pt}}}\,\theta_{\parallel}\,\theta_{\perp}^{2}}\frac{\Gamma(\kappa)}{\Gamma\left(\kappa-\frac{1}{2}\right)}$
$\displaystyle\times\left(1+\frac{(v_{\parallel}-u_{\parallel
0})^{2}}{\kappa\,\theta_{\parallel}^{2}}+\frac{v_{\perp}^{2}}{\kappa\,\theta_{\perp}^{2}}\right)^{-(\kappa+1)},$
(1)

where $n(\vec{r},t)$ is the number density (depending on location $\vec{r}$
and time $t$), $\Gamma(\mu)$ is the (complete) Gamma function, $v$ the
particle speed, which is normalized to a nominal thermal speed $\theta$,
$u_{\parallel 0}$ is the drift (or bulk) speed in parallel direction, and
$\kappa$ is a free parameter determining the width/slope of the suprathermal
tails. The symbols $\parallel$ and $\perp$ denote directions parallel and
perpendicular to the background magnetic field, respectively. The BK is
defined only for values of $\kappa>3/2$, and as $\kappa\to\infty$, it
approaches the standard bi-Maxwellian distribution function (BM):

$\displaystyle
f_{\mathrm{BM}}(\vec{r},\vec{v},t)=\frac{n(\vec{r},t)}{\mathchoice{{\hbox{$\displaystyle\sqrt{\pi^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\textstyle\sqrt{\pi^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\scriptstyle\sqrt{\pi^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=4.30276pt,depth=-3.44223pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\pi^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=3.44165pt,depth=-2.75334pt}}}\,\theta_{\parallel}\,\theta_{\perp}^{2}}\exp\left(-\frac{(v_{\parallel}-u_{\parallel
0})^{2}}{\theta_{\parallel}^{2}}-\frac{v_{\perp}^{2}}{\theta_{\perp}^{2}}\right)\,.$
(2)

If this BM approximates the core of the BK, it is described by the same
thermal speeds $\theta_{\parallel,\perp}$, see Lazar et al. (2015, 2016) for
an extended discussion motivating a definition of the kinetic temperature for
$\kappa$-distributed plasmas, i.e., $\theta_{\parallel,\perp}$ independent of
$\kappa$.

Despite their successful applications, standard $\kappa$-distribution models
are undermined by certain unphysical implications, such as diverging higher-
order velocity moments (which restricts values of the $\kappa$-parameter to
always exceed some minimum limits, although the observations may suggest even
lower values) and unrealistic contributions from superluminal particles
(Scherer et al., 2019a; Fahr & Heyl, 2020), and thus prevent a macroscopic
modeling of $\kappa$-distributed plasmas. These issues have been fixed by the
regularized bi-$\kappa$-distribution function (RBK) (Scherer et al., 2017,
2019a; Lazar et al., 2020)

$\displaystyle f_{\mathrm{RBK}}(\vec{r},$
$\displaystyle\vec{v},t)=\frac{n(\vec{r},t)}{\mathchoice{{\hbox{$\displaystyle\sqrt{\pi^{3}\,\kappa\,}$}\lower
0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\textstyle\sqrt{\pi^{3}\,\kappa\,}$}\lower
0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\scriptstyle\sqrt{\pi^{3}\,\kappa\,}$}\lower
0.4pt\hbox{\vrule
height=4.30276pt,depth=-3.44223pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\pi^{3}\,\kappa\,}$}\lower
0.4pt\hbox{\vrule
height=3.44165pt,depth=-2.75334pt}}}\,\theta_{\parallel}\,\theta_{\perp}^{2}}\frac{\Gamma(\kappa)}{\Gamma\left(\kappa-\frac{1}{2}\right)}\mathcal{W}_{[0,0]}\times$
$\displaystyle\left(1+\frac{(v_{\parallel}-u_{\parallel
0})^{2}}{\kappa\,\theta_{\parallel}^{2}}+\frac{v_{\perp}^{2}}{\kappa\,\theta_{\perp}^{2}}\right)^{-(\kappa+1)}e^{-\xi_{\parallel}\frac{(v_{\parallel}-u_{\parallel
0})^{2}}{\theta_{\parallel}^{2}}-\xi_{\perp}\frac{v_{\perp}^{2}}{\theta_{\perp}^{2}}}\,,$
(3)

where $\xi_{\parallel,\perp}$ are the cut-off parameters and
$\mathcal{W}_{[0,0]}$ is given by

$\mathcal{W}_{[0,0]}=\left(\int\limits_{0}^{1}\mathrm{d}t\,U\left(\frac{3}{2},\frac{3-2\,\kappa}{2},\kappa[\xi_{\perp}^{2}+(\xi_{\parallel}-\xi_{\perp})t^{2}]\right)\right)^{-1}$
(4)

with $U(a,c,x)$ being the Kummer-U or Tricomi function (e.g., Oldham et al.
(2010)). For $\xi_{\parallel}=\xi_{\perp}=0$ the RBK turns into the BK, and if
in addition $\kappa\to\infty$, the BM is recovered.

A universally accepted version would be useful to identify and differentiate
physical conditions when each of these versions of $\kappa$-distributions
applies (see Scherer et al. (2020) and references therein). Scherer et al.
(2020) proposed an isotropic generalized $\kappa$-distribution, named the
(isotropic) $\kappa$-cookbook, and exploring the implications of various
parameter setups. Here, in Section 2 we introduce an extended version able to
cover cases of anisotropic and multi-component distributions reported by the
observations, by introducing the generalized anisotropic $\kappa$-distribution
function (GAK), calling it the anisotropic $\kappa$-cookbook. All the commonly
used anisotropic model distributions, i.e., the BM, the BK and the RBK, can be
thus derived as particular cases of the GAK. We provide the reader with
general analytical expressions for the velocity moments (see Appendix A), when
one is not only concerned about the quality of the fits to the observations,
but also about the implications for macroscopic quantities when using a
certain distribution function. Most fitting procedures invoke 1D model
functions such as the BM or the BK, which are used to fit 1D cuts along
parallel, perpendicular or any other oblique direction of the observed
distributions or individual components (Pilipp et al., 1987a, b, c; Štverák et
al., 2008; Wilson et al., 2019b). While this is a common procedure, the
accuracy in reproducing the overall shape of the distribution and unveil
details (e.g., counterbeams mimicking an excess of parallel temperature) is
questionable. In Section 3 we present a new 2D fitting method using the GAK
with a Levenberg-Marquardt algorithm (Levenberg, 1944; Marquardt, 1963) to
minimize the nonlinear least-squares, which here is applied to selected
electron data sets, measured by the SWOOPS instrument of the Ulysses
spacecraft (Bame et al., 1992). The method allows us to depict electron
populations and sum-up their model distributions to find the best overall fit
and parameters controlling the quality of the fits. In Section 4 we present
the conclusions, pointing to the importance and the potential of a generalized
model combined with a 2D fitting to provide realistic descriptions of the
particle velocity distributions observed in space plasmas.

##  2 The anisotropic $\kappa$-cookbook

Similarly to the isotropic generalized $\kappa$-distribution function (GKD) in
Scherer et al. (2020), here we introduce the generalized anisotropic
$\kappa$-distribution function (GAK) as

$\displaystyle f_{GAK}$
$\displaystyle(\eta_{\parallel},\eta_{\perp},\zeta,\xi_{\parallel},\xi_{\perp})=n_{0}\,N_{{GAK}}$
(5) $\displaystyle\times\left(1+\frac{(v_{\parallel}-u_{\parallel
0})^{2}}{\eta_{\parallel}\Theta_{\parallel}^{2}}+\frac{v_{\perp}^{2}}{\eta_{\perp}\Theta_{\perp}^{2}}\right)^{-\zeta}\mathrm{e}^{-\xi_{\parallel}\frac{(v_{\parallel}-u_{\parallel
0})^{2}}{\Theta_{\parallel}^{2}}-\xi_{\perp}\frac{v_{\perp}^{2}}{\Theta_{\perp}^{2}}}\,,$

which we will call accordingly the anisotropic $\kappa$-cookbook. The 5-tuples
$(\eta_{\parallel},\eta_{\perp},\zeta,\xi_{\parallel},\xi_{\perp})$, which
define the cookbook, we will call ’recipes’. The integrals for the velocity
moments are quite similar to those used in Scherer et al. (2019b, 2020) and
are given in detail in Appendix A. The normalization constant
$N_{\mathrm{GAK}}$ in Eq. (5) is such that
$\int\mathrm{d}^{3}v\,f_{\mathrm{GAK}}$ gives the number density $n_{0}$, and
is obtained by

$\displaystyle N_{\mathrm{GAK}}^{-1}=$
$\displaystyle\mathchoice{{\hbox{$\displaystyle\sqrt{\pi^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\textstyle\sqrt{\pi^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\scriptstyle\sqrt{\pi^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=4.30276pt,depth=-3.44223pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\pi^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=3.44165pt,depth=-2.75334pt}}}\Theta_{\parallel}\Theta_{\perp}^{2}\eta_{\parallel}^{\frac{1}{2}}\eta_{\perp}$
(6)
$\displaystyle\times\int\limits_{0}^{1}U\left(\frac{3}{2},\frac{5}{2}-\zeta,(a_{2}+(a_{1}-a_{2})t^{2}\right)\mathrm{d}t\,.$

With the definition of the GAK by a tuple
$(\eta_{1},\eta_{2},\zeta,\xi_{1},\xi_{2})$, we can state the recipes for the
bi-Maxwellian as $(1,1,0,1,1)$, for the standard bi-$\kappa$-distribution as
$(\kappa,\kappa,\kappa+1,0,0)$ and for the regularized
bi-$\kappa$-distribution as $(\kappa,\kappa,\kappa+1,\xi_{1},\xi_{2})$, see
for further discussion of these distributions Scherer et al. (2019b). Here we
will focus on the fitting of the high resolution electron data from the SWOOPS
instrument of the Ulysses mission111http://ufa.esac.esa.int/ufa/#data. We used
also hourly average electron data (Bame et al., 1992, same webpage) and the
magnetic field data (Balogh et al., 1992, same webpage) for comparison and to
determine the Alfvén speed.

##  3 Data analysis

###  3.1 The model distributions

In order to differentiate between the anisotropic models used in the fitting
procedure, we use the following notations $f_{\alpha}$, with
$\alpha=\\{a,b,c,d,e\\}$ presented in Table 1, where RBK∗ is a modified
version of the RBK, allowing for an anisotropic $\kappa$-parameterization,
i.e., $\kappa_{\perp}\neq\kappa_{\parallel}$.

Index | Recipe | Distribution
---|---|---
a | (1,1,0,1,1) | Bi-Maxwellian
b | $(\eta_{{}_{\parallel}},\eta_{{}_{\perp}},\zeta,\xi_{{}_{\parallel}},\xi_{{}_{\perp}})$ | GAK
c | $(\kappa,\kappa,\kappa+1,\xi_{{}_{\parallel}},\xi_{{}_{\perp}})$ | RBK
d | $(\kappa_{{}_{\parallel}},\kappa{{}_{\perp}},\kappa_{{}_{\parallel}}+\kappa_{{}_{\perp}},\xi_{{}_{\parallel}},\xi_{{}_{\perp}})$ | RBK∗
e | $(\kappa,\kappa,\kappa+1,0,0)$ | BK
Table 1: Anisotropic distribution functions used in the fitting procedure in
Sec. 3.2

In general the electron data up to a few keV reveal three populations, a
thermal core and two suprathermal populations called halo and strahl, and
model distributions can be combined as

$\displaystyle f_{\alpha_{1}\alpha_{2}\alpha_{3}}\equiv
f_{\alpha_{1}}+f_{\alpha_{2}}+f_{\alpha_{3}}\,,$ (7)

where $\alpha_{i}\in\alpha$ with $i=\\{1,2,3\\}$ and the ordering is according
to the decreasing number density from the core (most dense) to the halo and
strahl, though the strahl may sometimes be more dense than the halo, e.g., in
the outer corona (Halekas et al., 2020) or the CME (single or double) strahls
(Skoug et al., 2000; Anderson et al., 2012).

$\displaystyle n_{\alpha_{1}}>n_{\alpha_{2}}\geqslant n_{\alpha_{3}}$ (8)

For example, the fit for the observed distribution function $f_{abc}$ can be
the sum of the core fit with the recipe (1,1,0,1,1), and two suprathermal
components, one fitted by the recipe
$(\eta_{{}_{\parallel}},\eta_{{}_{\perp}},\zeta,\xi_{{}_{\parallel}},\xi_{{}_{\perp}})$,
and the other one by
$(\kappa,\kappa,\kappa+1,\xi_{{}_{\parallel}},\xi_{{}_{\perp}})$ (with the
lowest number density $n_{\alpha_{3}}$). The bi-$\kappa$-distribution function
is not used explicitly in the fitting procedure, but is included in the
application of the GAK and the RBK, e.g., if in the RBK
$\xi_{\parallel}=\xi_{\perp}=0$.

The velocity moments for a sum of distribution functions are already discussed
in Scherer et al. (2020):

$\displaystyle n_{\alpha_{1}\alpha_{2}\alpha_{3}}$
$\displaystyle=\sum\limits_{i=1}^{3}n_{\alpha_{i}}$ (9)
$\displaystyle\vec{u}_{cms}$
$\displaystyle\equiv\vec{u}_{\alpha_{1}\alpha_{2}\alpha_{3}}=\frac{\sum\limits_{i=1}^{3}n_{i}\vec{u}_{\alpha_{i}}}{\sum\limits_{i=1}^{3}n_{\alpha_{i}}}=(u_{{}_{\parallel},cms},\vec{u}_{{}_{\perp}})^{T}$
(10) $\displaystyle\overleftrightarrow{P}$
$\displaystyle=\sum\limits_{i=1}^{3}\int
f_{i}(\vec{v})\left[\vec{v}\otimes\vec{v}+(\vec{u}-\vec{u}_{i})\otimes(\vec{u}-\vec{u}_{i})\right]\mathrm{d}^{3}v$
(11) $\displaystyle P_{{}_{\parallel}}$
$\displaystyle\equiv\sum\limits_{i=1}^{3}P_{\alpha_{i},11}+\sum\limits_{i=1}^{3}n_{\alpha_{i}}(u_{{}_{\parallel},cms}-u_{{}_{\parallel},i})^{2}$
(12) $\displaystyle P_{{}_{\perp}}$
$\displaystyle\equiv\sum\limits_{i=1}^{3}P_{\alpha_{i},22}$ (13)
$\displaystyle\vec{H}$
$\displaystyle=\frac{1}{2}\sum\limits_{i=1}^{3}\left(2\overleftrightarrow{P}_{\alpha_{i}}+\mathrm{tr}(\overleftrightarrow{P})_{\alpha_{i}}-n_{\alpha_{i}}U_{\alpha_{i}}^{2}\right)\vec{U}_{\alpha{i}}$
(14) $\displaystyle{\color[rgb]{0,0,0}H_{{}_{\parallel}}}$
$\displaystyle{\color[rgb]{0,0,0}=\sum\limits_{i=1}^{3}\left(3P_{\alpha_{i},11}+2P_{\alpha_{i},22}-n_{\alpha_{i}}U_{{}_{\parallel},\alpha_{i}}^{2}\right)U_{{}_{\parallel},\alpha_{i}}}$
(15) $\displaystyle H_{{}_{\perp}}$ $\displaystyle=0$ (16)

with $\vec{U}_{\alpha{i}}=\vec{u}_{cms}-\vec{u}_{\alpha{i}}$ , where
$\vec{u}_{cms}$ is the center of mass velocity. $\overleftrightarrow{P}$
denotes the thermal pressure tensor and $\vec{H}$ is the heat flux vector. It
is always assumed that the perpendicular drift speed is zero (see discussion
below) and thus $\vec{u}_{{}_{\perp},\alpha{i}}\equiv 0$ and
$\vec{u}_{{}_{\perp},cms}=0$, and we are only left with the parallel
components in the center of mass velocity as well as in the heat flux
$\vec{H}$. The above moments are all normalized to the electron mass.

Figure 1: 3D representation of the data (2002, doy 14,22:09:58) by a
Maxwellian, a GAK and a Maxwellian: On the x-axis the parallel speed
$v_{{}_{\parallel}}$ and on the y-axis the perpendicular speed
$v_{{}_{\perp}}$ in units of [Mm/s] are given, while on the logarithmic z-axis
the distribution function in [cm-3(Mm/s)-3] is shown. The black crosses
represent the data, while the colored grids show the following fits: top left
panel a Maxwellian fit of the core, top right a GAK fit of the suprathermal
part, bottom left a shifted Maxwellian fit of the strahl, bottom right the sum
of all three model distributions. One can see that the data are best fitted
with the sum of three distribution functions. The fits are done
simultaneously, and the single distributions are shown in the plots above.

In Fig. 1 we show the data (black crosses) and the simultaneously fitted
distribution functions (colored grids). In the upper left panel only the core
is fitted with a Maxwellian (Eq. 2), in the upper right panel only the first
suprathermal part is fitted with a GAK (Eq. 5), in the lower left panel only
the second suprathermal part is fitted with a Maxwelian (Eq. 2), and finally
in the lower right panel the fit in form of the sum of all three model
distribution is presented.

From Fig. 1 it becomes evident that the data consists of three populations: a
core and two suprathermal populations, where the hotter population is called
strahl and the less hot population constitutes the halo. The best fit is then
obtained with the sum of three model distributions. We fitted all three
distribution functions at once, and the single distributions shown in Fig.1
are then the single distributions taken from the fit.

###  3.2 Data fitting

We use a nonlinear least-squares fit based on the Marquardt-Levenberg
(Levenberg, 1944; Marquardt, 1963) algorithm to find the best (combination of)
model distribution(s). During that procedure we fitted all possible
combinations of only one (global) fit $f_{\alpha_{i}}$, dual fits
$f_{\alpha_{i}\alpha_{j}}$, or triple fits combining three models
$f_{\alpha_{i}\alpha_{j}\alpha_{k}}$ as shown in Eq. (7) (for more details see
Appendix B). f we fit more than one distribution the fit to all distributions
is always done simultaneously. This prevents us from deciding which part of
the data belongs to a given distribution. For testing we fitted first a
Maxwellian (for $v<3$ Mm/s, then using that initial values two distributions,
and with the same procedure at least three distributions. There were only
negligible differences to the full fit, and thus we always do a full fit.

For the accuracy of fitting we defined the relative error

$\displaystyle E_{i}=\frac{|f_{fit,i}-f_{obs,i}|}{f_{obs,i}}$ (17)

for each data point $i$, and the mean relative error $<E>$, its standard
deviation $\sigma$ and the $R^{2}$ value as

$\displaystyle<E>$ $\displaystyle=\frac{1}{N}\sum\limits_{i=1}^{N}E_{i}\,,$
(18) $\displaystyle\sigma$
$\displaystyle=\mathchoice{{\hbox{$\displaystyle\sqrt{\frac{1}{N-1}\sum\limits_{i=1}^{N}(<E>-E_{i})^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=9.8611pt,depth=-7.88892pt}}}{{\hbox{$\textstyle\sqrt{\frac{1}{N-1}\sum\limits_{i=1}^{N}(<E>-E_{i})^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=9.30444pt,depth=-7.44359pt}}}{{\hbox{$\scriptstyle\sqrt{\frac{1}{N-1}\sum\limits_{i=1}^{N}(<E>-E_{i})^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=6.53888pt,depth=-5.23112pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\frac{1}{N-1}\sum\limits_{i=1}^{N}(<E>-E_{i})^{2}\,}$}\lower
0.4pt\hbox{\vrule height=5.03888pt,depth=-4.03113pt}}}\,,$ (19) $\displaystyle
R^{2}$
$\displaystyle=1-\frac{(\sum\limits_{i=1}^{N}E_{i}-<E>)^{2}}{\sum\limits_{i=1}^{N}E_{i}^{2}}\,.$
(20)

The values for the mean relative error $<E>$ and the standard deviation
$\sigma$ are required to be as small as possible for a satisfying fit. Here,
by direct visual inspection (comparison of the different fits) we define
values below 0.3 for $<E>$ and $\sigma$ as a “good” fit, and thus only such
fits are considered in the following. However, in order to evaluate the
goodness of the fits not only visually, but to have a quantitative measure, we
use the $R^{2}$ values. Once $<E>$ and $\sigma$ are below 0.3, the highest
$R^{2}$ value determines the best fit, while the limit $R^{2}=1$ is the
idealized best fit (i.e., $<E>\to 0$ and $\sigma\to 0$).

In the present analysis we compare two data sets arbitrarily chosen from the
year 2002: day 288, 00:08:17 (hh:mm:ss), and day 365, 00:13:45. Table 2
presents all tested model distributions for these data sets, the values
obtained for $<E>$ and $\sigma$ with the corresponding $R^{2}$ values. It is
to be observed that for day 365 (right column) not $f_{\mathrm{adb}}$ gives
the best fit (despite the lowest $<E>$ and $\sigma$ values), but
$f_{\mathrm{aba}}$ due to the highest $R^{2}$ value.

day 228 | day 365
---|---
$f_{\alpha_{i}\alpha_{j}\alpha_{k}}$ | $<E>$ | $\sigma$ | $R^{2}$ | $f_{\alpha_{i}\alpha_{j}\alpha_{k}}$ | $<E>$ | $\sigma$ | $R^{2}$
aaa | 0.28 | 0.24 | 0.57 | aaa | 0.21 | 0.18 | 0.58
aba | 0.23 | 0.20 | 0.56 | aba | 0.28 | 0.19 | 0.69
abb | 0.27 | 0.18 | 0.70 | - | - | - | -
- | - | - | - | abc | 0.16 | 0.15 | 0.55
aca | 0.29 | 0.27 | 0.54 | aca | 0.19 | 0.19 | 0.51
- | - | - | - | acb | 0.17 | 0.16 | 0.54
add | 0.29 | 0.24 | 0.61 | - | - | - | -
- | - | - | - | ada | 0.17 | 0.17 | 0.50
- | - | - | - | adb | 0.15 | 0.14 | 0.54

Table 2: The mean average error $<E>$, the standard deviation $\sigma$ and the
$R^{2}$. The model distribution which gives the best fit is marked in bold for
the particular day. The $-$ sign indicates that these values of $<E>$ and
$\sigma$ are above 0.3. All other combinations above 0.3 are left out.

In order to show differences of accuracy, in Fig. 2 we have compared the best
fits from Table 2 with another bad fit, arbitrarily chosen. In this figure
both electron data sets are used. The black contour lines show the model fits,
while the red contour lines show that of the data. For a good fit the black
and red contour lines should coincide. For a good fit (the top and bottom left
panels, respectively), these contours are close to each other, while in the
bad fits (the top and bottom right panels, respectively), they markedly
differ. Thus our choice first excludes all the $<E>$ and $\sigma$ values which
are larger than 0.3, and than chooses from the remaining models that for which
the $R^{2}$ value is the largest as a good choice. Nevertheless, one should
keep in mind that the nonlinear least-squares fit may find a local minimum,
which gives not the best fit. Therefore, we have tested different initial
values, and in rare cases, this leads to better fitting models. The presented
fits are robust against different initial values, and the best model is in all
tests the same with the same fitted values (up to rounding errors).

Figure 2: Contour plots of the fitted distribution functions. The upper panels
show two fitted models at the date 2002, day 288, 00:08:17: the upper left
panel shows the best model $f_{abb}$ with $<E>=0.27$, $\sigma=0.18$, and
$R^{2}=0.70$, while in the upper right panel a “bad” model $f_{acc}$ is shown
with $<E>=0.5$, $\sigma=0.37$, and $R^{2}=0.65$. In the lower panels two
fitted models at the date 2002, day 365, 00:13:45 are shown: the lower left
panel shows the best model $f_{aba}$ with $<E>=0.28$, $\sigma=0.19$, and
$R^{2}=0.69$, while on the lower right panel a “bad” model $f_{acc}$ is shown
with $<E>=0.81$, $\sigma=0.14$, and $R^{2}=0.97$. The bad fits are excluded
because of the large $<E>$ or $\sigma$ values (above our defined threshold of
0.3, see text for details). The black contour lines are those of the model
fit, while the red contours are those from the data. The blue contour lines
show the levels of the third distributions, which are in the above case too
small to contribute substantially to those of the total distribution function.

In Table 3 we show the best fit values for the two discussed cases. With help
of the formulas (10), (12), (14), (15), and (29) we can calculate the center
of mass speed, the pressure and the heat flow using the fitted values. The
results are given in Table 4 together with the integrated SWOOPS electron data
(see Bame et al., 1992), which give the isotropic moment for the core and halo
as well as the magnetic field and the resulting Alfvén speed. To make the
fitted anisotropic temperatures compatible to the isotropic temperatures, we
estimated the corresponding isotropic pressure as
$P_{iso}=P_{{}_{\parallel}}+2P_{{}_{\perp}}$, and from there the isotropic
temperatures. One should keep in mind that the ion data data are averaged over
1 hour, while the electron data are given in 2 min intervals (with a much
longer cadence time of 3-17 min.)

day | $f_{\alpha_{i}\alpha_{j}\alpha_{k}}$ |  $n_{i}$ [cm-3] |  $\Theta_{\parallel,i}$ [Mm/s] |  $\Theta_{\perp,i}$ [Mm/s] | $\eta_{\parallel,i}\kappa$ | $\eta_{\perp,i}$ | $\zeta$ | $\xi_{\parallel,i}$ | $\xi_{\perp,i}$ |  $u_{\parallel,i}$[Mm/s]
---|---|---|---|---|---|---|---|---|---|---
288 | $\alpha_{i}=a$ | $1.20\cdot 10^{-1}$ | $8.08\cdot 10^{-1}$ | $7.50\cdot 10^{-1}$ |  |  |  |  |  | $1.05\cdot 10^{-1}$
| $\alpha_{j}=b$ | $1.42\cdot 10^{-3}$ | $2.37\cdot 10^{-1}$ | $1.94\cdot 10^{-1}$ | $6.59\cdot 10^{-2}$ | $4.80\cdot 10^{-2}$ | $5.01\cdot 10^{-1}$ | $1.36\cdot 10^{-3}$ | $9.86\cdot 10^{-4}$ | $-5.25\cdot 10^{-1}$
| $\alpha_{k}=b$ | $9.68\cdot 10^{-3}$ | $1.17\cdot 10^{1}$ | $1.00\cdot 10^{1}$ | $3.98\cdot 10^{-1}$ | $3.08\cdot 10^{-1}$ | $1.21\cdot 10^{1}$ | $4.51\cdot 10^{-3}$ | $1.12\cdot 10^{-3}$ | $-1.72\cdot 10^{-1}$
365 | $\alpha_{i}=a$ | $1.34\cdot 10^{-1}$ | $1.02\cdot 10^{0}$ | $9.23\cdot 10^{-1}$ |  |  |  |  |  | $-4.81\cdot 10^{-2}$
| $\alpha_{j}=b$ | $5.63\cdot 10^{-3}$ | $2.34\cdot 10^{0}$ | $2.30\cdot 10^{0}$ | $3.31\cdot 10^{0}$ | $3.11\cdot 10^{0}$ | $4.58\cdot 10^{0}$ | $9.19\cdot 10^{-4}$ | $5.88\cdot 10^{-5}$ | $-6.97\cdot 10^{-1}$
| $\alpha_{k}=a$ | $1.47\cdot 10^{-4}$ | $1.45\cdot 10^{0}$ | $2.51\cdot 10^{0}$ |  |  |  |  |  | $2.61\cdot 10^{0}$

Table 3: Best fit parameters.

day | function |  $n_{tot}$ [cm-3] |  $n_{core}$ [cm-3] |  $n_{halo}$ [cm-3] |  $u_{cms}$ [Mm/s] |  $T_{tot}$ [K] |  $T_{core}$ [K] |  $T_{halo}$ [K] |  $H_{{}_{\parallel}}$ [gs-3] |  $|\vec{B}|$ [nT] |  $v_{A}$ [Mm/s]
---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  |  |  $u_{SW}$ [Mm/s] |  $T_{ion}$ [K] |  |  |  |  |
288 | abb | $1.7\cdot 10^{-1}$ | $2.7\cdot 10^{-2}$ | $1.2\cdot 10^{-6}$ | $7.8\cdot 10^{-2}$ | $1.7\cdot 10^{5}$ | $1.6\cdot 10^{5}$ | $3.1\cdot 10^{5}$ | $-7.6\cdot 10^{-5}$ |  |
| isotropic | $1.5\cdot 10^{-1}$ | $1.2\cdot 10^{-1}$ | $2.3\cdot 10^{-2}$ |  | $1.2\cdot 10^{5}$ | $4.4\cdot 10^{4}$ | $5.4\cdot 10^{6}$ |  | $7.7\cdot 10^{-1}$ | $4.2\cdot 10^{-2}$
| ion | $1.4\cdot 10^{-1}$ |  |  | $5.1\cdot 10^{-1}$ | $5.8\cdot 10^{5}$ |  |  |  |  |
365 | aba | $1.4\cdot 10^{-1}$ | $1.3\cdot 10^{-1}$ | $5.6\cdot 10^{-3}$ | $-7.1\cdot 10^{-1}$ | $7.3\cdot 10^{5}$ | $3.1\cdot 10^{5}$ | $4.6\cdot 10^{5}$ | $1.2\cdot 10^{-5}$ |  |
| isotropic | $3.8\cdot 10^{-1}$ | $3.7\cdot 10^{-1}$ | $1.2\cdot 10^{-2}$ |  | $6.3\cdot 10^{5}$ | $4.7\cdot 10^{5}$ | $5.2\cdot 10^{6}$ |  | $9.38\cdot 10^{-1}$ | $5.1\cdot 10^{-2}$
| ion | $4.1\cdot 10^{-1}$ |  |  | $4.2\cdot 10^{-2}$ | $6.5\cdot 10^{5}$ |  |  |  |  |

Table 4: The moments for the fitted anisotropic model distributions and those
from the isotropic data. The magnetic field and the Alvfén speed are also
displayed. Additionally, the rows indicated by “ion” give the proton number
densities, speeds and temperatures from the SWOOPS instrument. The temperature
is calculated from the pressure $P_{tot}=P_{{}_{\parallel}}+2P_{{}_{\perp}}$
using the ideal gas law for the core, halo and the total pressure (using Eq.
12). The position of the Ulysses spacecraft on day 288 is at $r=4.2$ au,
$\vartheta=29.88^{o}$, $\varphi=18.09^{o}$, and on day 365 at $r=4.5$ au,
$\vartheta=24.00^{o}$, $\varphi=8.94^{o}$.

From Table 3 one can also see that fits of the anisotropic distributions give
similar results as those for the isotropic data and the ion number density.
The match is not perfect, possibly for various reasons: 1. The data recording
does not occur at the same time, or 2\. the data are averaged over different
time intervals, or 3. there may be major differences between isotropic and
anisotropic distribution functions. However, the number densities and
temperatures are in the same range. The difference also affects the third
fitting model, which can be much hotter than expected (at day 288
$T_{{}_{\parallel},3}=1.6\cdot 10^{8}$ [K] and $T_{{}_{\perp},3}=1.2\cdot
10^{8}$ [K] and at day 365 $T_{{}_{\parallel},3}=5.4\cdot 10^{6}$ [K] and
$T_{{}_{\perp},3}=2.1\cdot 10^{6}$ [K] ).

Figure 3: Finally we show here the contour lines for distributions shown in
Fig 1. The best fit is given by the a Bi-Maxwellian, a GAK, and again a
shifted Bi-Maxwellian. The data are taken at 2002, doy 14, 22:09:52 see also
(Lazar et al., 2014) for the same data set. It can be seen that the red
contours for the observed data scatter a lot and our quality criteria fail. We
thus have chosen the fit with the lowest mean error (0.33). The black contours
show the fit with all three distributions, and the blue contours only the fit
for the strahl.

In Fig. 3 we show an event including a distinct strahl. Unfortunately, the
data scatter a lot, so that the fit has larger mean errors and standard
deviations. Therefore, we have chosen the lowest values for $<E>$ and $\sigma$
(neglecting $R^{2}$) for the best fit. The number density, temperature and
drift speed for the three fitted distributions (Maxwellian, GAK, Maxwellian)
are shown in Table 5. For convenience we have also added the temperature
anisotropy in Table 5. It can be seen that the core is a "cold" isotropic
distribution, the halo is a hot moderate anisotropic distribution, and the
strahl is an even hotter highly anisotropic distribution. A discussion of the
physical implications, like heat transfer, will go beyond the scope of the
recent discussion and will be object to further work.

$f_{a}$ |  $n$[cm${}^{-}3$] |  $T_{{}_{\parallel}}$[K] |  $T_{{}_{\perp}}$[K] |  $u_{{}_{\parallel}}$[Mm/s] | $A=T_{{}_{\perp}}/T_{{}_{\parallel}}$
---|---|---|---|---|---
a | $2.7\cdot 10^{-1}$ | 23194 | 22792 | 0.14 | 0.99
b | $1.7\cdot 10^{-2}$ | 270454 | 172650 | -0.36 | 0.64
a | $2.2\cdot 10^{-5}$ | 6584576 | 135541 | -8.80 | 0.02
Table 5: The moments for the distribution function shown in Fig. 3. For the
GAK the remaining parameter are $\eta_{{}_{\parallel}}=1.3$,
$\eta_{{}_{\perp}}=1.2$, $\zeta=3.1$, $\xi_{{}_{\parallel}}=2.6\cdot 10^{-2}$,
and $\xi_{{}_{\perp}}=9.5\cdot 10^{-5}$

To fit the distribution in Fig. 3 our above described quality parameter for
the goodness of the fit failed, because of the data scattering. In that case
one has either to weaken the conditions or to check the fits by hand.
Nevertheless, for quite smooth data our quality assessment is sufficient. In
the data set where the data points seem to scatter, an additional provided
error estimate for each data point would be helpful, because then one can use
a weighted fitting procedure.

##  4 Conclusions

In hot and dilute plasmas from space, studies of fundamental processes, e.g.,
particle heating, instabilities, etc., require laborious kinetic approaches
conditioned by the velocity distributions of plasma particles and their basic
features, e.g, anisotropies with respect to the magnetic field combined with
suprathermal high-energy tails described by the $\kappa$-distribution
functions. Recently, this empirical model has been regularized in order to
provide a well-defined macroscopic description of $\kappa$-distributed plasmas
as fluids (Scherer et al., 2017, 2019a; Lazar et al., 2020). Scherer et al.
(2020) have also introduced the $\kappa$-cookbook, a generalization such that
the multitude of forms of $\kappa$-distributions but also Maxwellian limits
invoked in the literature are recovered as particular cases or recipes. In the
present paper we have advanced this notion by including the anisotropic
distributions, often reported by the observations, in a new generalized
anisotropic $\kappa$-cookbook (GAK), and provided the reader with analytical
expressions of the integrals for the velocity moments.

In addition, this generalized model has been applied to describe the velocity
distributions measured in-situ in the solar wind, by introducing a new fitting
method with two-dimensional (2D) models. We demonstrated the quality of the
method for samples of high resolution electron data sets provided by Ulysses
missions. It turned out that the moments of our best fit distribution models
are in a good agreement with the isotropic electron and ion data. The best fit
may include various combinations of the GAK recipes, which can vary from a
data set to another. These can be directly related to different properties of
plasma species, or even their components, e.g., electron core, halo or strahl
populations, conditioned by the solar activity, heliocentric distance and
solar wind processes. Numerical reasons can also be invoked, because for
example, for large $\kappa$-values the distribution function becomes a
Maxwellian, which numerically can lead to a $\kappa-$model distribution
instead of a Maxwellian model distribution. From a fitting point of view these
differences are marginal, but from a theoretical point of view this can cause
unwanted complications (for example, with Maxwellian distributions analytic
solutions are often possible, while for non-Maxwellian they are not). These
results promise extended applications in the observations and modeling of the
observed distributions and their dependence with distance, latitude and
longitude.

##  5 Data availability statement

The high resolution electron data are from the SWOOPS instrument of the
Ulysses mission: http://ufa.esac.esa.int/ufa/#data. We used also hourly
average electron data (Bame et al., 1992, same webpage) and the magnetic field
data (Balogh et al., 1992, same webpage)

## References

  * Anderson et al. (2012) Anderson, B. R., Skoug, R. M., Steinberg, J. T., & McComas, D. J. 2012, Journal of Geophysical Research: Space Physics, 117, A04107
  * Astfalk & Jenko (2016) Astfalk, P. & Jenko, F. 2016, Journal of Geophysical Research (Space Physics), 121, 2842
  * Balogh et al. (1992) Balogh, A., Beek, T. J., Forsyth, R. J., et al. 1992, A&AS, 92, 221
  * Bame et al. (1992) Bame, S. J., McComas, D. J., Barraclough, B. L., et al. 1992, A&AS, 92, 237
  * Elzhov et al. (2016) Elzhov, T. V., Mullen, K. M., Spiess, A. N., & Bolker, B. 2016, Package ‘minpack.lm’, Website, online available under https://cran.r-project.org/web/packages/minpack.lm/index.html; accessed:
  * Fahr & Heyl (2020) Fahr, H. J. & Heyl, M. 2020, MNRAS, 491, 3967
  * Gloeckler (2003) Gloeckler, G. 2003, in American Institute of Physics Conference Series, Vol. 679, Solar Wind Ten, ed. M. Velli, R. Bruno, F. Malara, & B. Bucci, 583–588
  * Halekas et al. (2020) Halekas, J. S., Whittlesey, P., Larson, D. E., et al. 2020, ApJS, 246, 22
  * Heikkila (2011) Heikkila, W. J. 2011, Earth’s Magnetosphere, 1st edn. (Elsevier Publications)
  * Kasper et al. (2006) Kasper, J. C., Lazarus, A. J., Steinberg, J. T., Ogilvie, K. W., & Szabo, A. 2006, J. Geophys. Res., 111, A03105
  * Lazar et al. (2016) Lazar, M., Fichtner, H., & Yoon, P. H. 2016, A&A, 589, A39
  * Lazar et al. (2012a) Lazar, M., Pierrard, V., Poedts, S., & Schlickeiser, R. 2012a, Astrophysics and Space Science Proceedings, 33, 97
  * Lazar & Poedts (2009) Lazar, M. & Poedts, S. 2009, A&A, 494, 311
  * Lazar et al. (2015) Lazar, M., Poedts, S., & Fichtner, H. 2015, A&A, 582, A124
  * Lazar et al. (2014) Lazar, M., Pomoell, J., Poedts, S., Dumitrache, C., & Popescu, N. A. 2014, Sol. Phys., 289, 4239–4266
  * Lazar et al. (2020) Lazar, M., Scherer, K., Fichtner, H., & Pierrard, V. 2020, A&A, 634, A20
  * Lazar et al. (2012b) Lazar, M., Schlickeiser, R., & Poedts, S. 2012b, in Exploring the Solar Wind, ed. M. Lazar (Rijeka: IntechOpen)
  * Levenberg (1944) Levenberg, K. 1944, Quart. Appl. Math., 2, 164
  * Macneil et al. (2020) Macneil, A. R., Owen, M. J., Lockwood, M., Štverák, Š., & Owen, C. J. 2020, Sol. Phys., 295, 16
  * Maksimovic et al. (1997) Maksimovic, M., Pierrard, V., & Riley, P. 1997, Geophys. Res. Lett., 24, 1151
  * Maksimovic et al. (2005) Maksimovic, M., Zouganelis, I., Chaufray, J.-Y., et al. 2005, Journal of Geophysical Research (Space Physics), 110, A09104
  * Marquardt (1963) Marquardt, D. W. 1963, J. Soc. Ind. Appl. Math., 11, 431
  * Marsch (2006) Marsch, E. 2006, Living Rev. Sol. Phys., 3, 1
  * Olbert (1968) Olbert, S. 1968, in Astrophysics and Space Science Library, Vol. 10, Physics of the Magnetosphere, ed. R. D. L. Carovillano & J. F. McClay, 641
  * Oldham et al. (2010) Oldham, K., Myland, J., & Spanier, J. 2010, An Atlas of Functions: with Equator, the Atlas Function Calculator, An Atlas of Functions (Springer New York)
  * Pierrard & Lazar (2010) Pierrard, V. & Lazar, M. 2010, Sol. Phys., 267, 153
  * Pierrard et al. (2001) Pierrard, V., Maksimovic, M., & Lemaire, J. 2001, Astrophys. Space Sci., 277, 195–200
  * Pilipp et al. (1987a) Pilipp, W. G., Miggenrieder, H., Montgomery, M. D., et al. 1987a, J. Geophys. Res.: Space Phys., 92, 1075
  * Pilipp et al. (1987b) Pilipp, W. G., Miggenrieder, H., Montgomery, M. D., et al. 1987b, J. Geophys. Res.: Space Phys., 92, 1093
  * Pilipp et al. (1987c) Pilipp, W. G., Miggenrieder, H., Mühlhäuser, K. H., et al. 1987c, J. Geophys. Res.: Space Phys., 92, 1103
  * Scherer et al. (2019a) Scherer, K., Fichtner, H., Fahr, H. J., & Lazar, M. 2019a, ApJ, 881, 93
  * Scherer et al. (2017) Scherer, K., Fichtner, H., & Lazar, M. 2017, EPL (Europhysics Letters), 120, 50002
  * Scherer et al. (2020) Scherer, K., Husidic, E., Lazar, M., & Fichtner, H. 2020, MNRAS, 497, 1738–1756
  * Scherer et al. (2019b) Scherer, K., Lazar, M., Husidic, E., & Fichtner, H. 2019b, ApJ, 880, 118
  * Skoug et al. (2000) Skoug, R. M., Feldman, W. C., Gosling, J. T., McComas, D. J., & Smith, C. W. 2000, Journal of Geophysical Research: Space Physics, 105, 23069
  * Summers & Thorne (1992) Summers, D. & Thorne, R. M. 1992, J. Geophys. Res., 97, 16827
  * Štverák et al. (2008) Štverák, Š., Trávníček, P., Maksimovic, M., et al. 2008, Journal of Geophysical Research (Space Physics), 113, A03103
  * Vasyliunas (1968) Vasyliunas, V. M. 1968, in Astrophysics and Space Science Library, Vol. 10, Physics of the Magnetosphere, ed. R. D. L. Carovillano & J. F. McClay, 622
  * Wilson et al. (2019b) Wilson, L. B., Chen, L.-J., Wang, S., et al. 2019b, Astrophys. J. Suppl.S., 243, 8
  * Wilson et al. (2019a) Wilson, L. B., Chen, L.-J., Wang, S., et al. 2019a, Astrophys. J. Suppl.S., 245, 24
  * Wilson et al. (2020) Wilson, L. B., Chen, L.-J., Wang, S., et al. 2020, Astrophys. J., 893, 22

##  Appendix A The integrals

The integrals for the velocity moments $M^{\lambda,\mu}$ of order
$\lambda,\mu$ of the anistropic cookbook (Eq. 5) have the form (neglecting the
factors $n_{0}$ and $N_{\mathrm{GAK}}$ and having already taken into account
the drift speed)

$\displaystyle M^{\lambda,\mu}$
$\displaystyle=\int\limits_{0}^{2\pi}\int\limits_{-\infty}^{\infty}\int\limits_{0}^{\infty}f_{Gak}(\eta_{\parallel},\eta_{\perp},\zeta,\xi_{\parallel},\xi_{\perp})v_{\parallel}^{\lambda}v_{\perp}^{\mu+1}\mathrm{d}v_{\parallel}\mathrm{d}v_{\perp}\mathrm{d}\varphi$
(21)
$\displaystyle=\int\limits_{0}^{2\pi}\int\limits_{-\infty}^{\infty}\int\limits_{0}^{\infty}\left(1+\frac{v_{\parallel}^{2}}{\eta_{\parallel}\Theta_{\parallel}^{2}}+\frac{v_{\perp}^{2}}{\eta_{\perp}\Theta_{\perp}^{2}}\right)^{-\zeta}\mathrm{e}^{-\xi_{\parallel}\frac{v_{\parallel}^{2}}{\Theta_{\parallel}^{2}}-\xi_{\perp}\frac{v_{\perp}^{2}}{\Theta_{\perp}^{2}}}$
$\displaystyle\hskip 85.35826pt\times
v_{\parallel}^{\lambda}v_{\perp}^{\mu+1}\mathrm{d}v_{\parallel}\mathrm{d}v_{\perp}\mathrm{d}\varphi$
$\displaystyle=2\pi\Theta_{\parallel}^{1+\lambda}\Theta_{\perp}^{2+\mu}\eta_{\parallel}^{\frac{\lambda+1}{2}}\eta_{\perp}^{\frac{\mu+2}{2}}\int\limits_{-\infty}^{\infty}\int\limits_{0}^{\infty}\left(1+w_{\parallel}^{2}+w^{2}_{\perp}\right)^{-\zeta}$
$\displaystyle\hskip
56.9055pt\times\mathrm{e}^{-\xi_{\parallel}\eta_{\parallel}\Theta_{\parallel}^{2}w_{\parallel}^{2}-\xi_{\perp}\eta_{\perp}\Theta_{\perp}^{2}w_{\perp}^{2}}w_{\parallel}^{\lambda}w_{\perp}^{\mu+1}\mathrm{d}w_{\parallel}\mathrm{d}w_{\perp}$
$\displaystyle=2\pi\Theta_{\parallel}^{1+\lambda}\Theta_{\perp}^{2+\mu}\eta_{\parallel}^{\frac{\lambda+1}{2}}\eta_{\perp}^{\frac{\mu+2}{2}}\int\limits_{0}^{\pi}\int\limits_{0}^{\infty}r^{\lambda+\mu+2}\left(1+r^{2}\right)^{-\zeta}$
$\displaystyle\hskip
56.9055pt\times\mathrm{e}^{-r^{2}(a_{2}+(a_{1}-a_{2})\cos^{2}\vartheta)}\sin\vartheta^{\mu+1}\cos\vartheta^{\lambda}\mathrm{d}r\mathrm{d}\vartheta\,,$

where we have first replaced

$\displaystyle
v_{\parallel}=\mathchoice{{\hbox{$\displaystyle\sqrt{\eta_{\parallel}\,}$}\lower
0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\textstyle\sqrt{\eta_{\parallel}\,}$}\lower
0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\scriptstyle\sqrt{\eta_{\parallel}\,}$}\lower
0.4pt\hbox{\vrule
height=3.01389pt,depth=-2.41113pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\eta_{\parallel}\,}$}\lower
0.4pt\hbox{\vrule
height=2.15277pt,depth=-1.72223pt}}}\Theta_{\parallel}w_{\parallel}\,,\qquad
v_{\perp}=\mathchoice{{\hbox{$\displaystyle\sqrt{\eta_{\perp}\,}$}\lower
0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\textstyle\sqrt{\eta_{\perp}\,}$}\lower
0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\scriptstyle\sqrt{\eta_{\perp}\,}$}\lower
0.4pt\hbox{\vrule
height=3.01389pt,depth=-2.41113pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\eta_{\perp}\,}$}\lower
0.4pt\hbox{\vrule
height=2.15277pt,depth=-1.72223pt}}}\,,\Theta_{\perp}w_{\perp}$ (22)

and next

$\displaystyle w_{\parallel}=r\cos\vartheta\,,\qquad
w_{\perp}=r\sin\vartheta\,,\qquad
r=\mathchoice{{\hbox{$\displaystyle\sqrt{w_{\parallel}^{2}+w_{\perp}^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\textstyle\sqrt{w_{\parallel}^{2}+w_{\perp}^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\scriptstyle\sqrt{w_{\parallel}^{2}+w_{\perp}^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=4.30276pt,depth=-3.44223pt}}}{{\hbox{$\scriptscriptstyle\sqrt{w_{\parallel}^{2}+w_{\perp}^{2}\,}$}\lower
0.4pt\hbox{\vrule height=3.44165pt,depth=-2.75334pt}}}\,,$ (23)

and used the shorthand notation

$\displaystyle a_{1}=\xi_{\parallel}\eta_{\parallel}\,,\qquad
a_{2}=\xi_{\perp}\eta_{\perp}\,.$ (24)

Next we use $z=r^{2}$ and $t=\cos\vartheta$ to obtain

$\displaystyle M^{\lambda,\mu}=$ $\displaystyle
2\pi\Theta_{\parallel}^{1+\lambda}\Theta_{\perp}^{2+\mu}\eta_{\parallel}^{\frac{\lambda+1}{2}}\eta_{\perp}^{\frac{\mu+2}{2}}\int\limits_{0}^{1}\int\limits_{0}^{\infty}z^{\frac{\lambda+\mu+1}{2}}\left(1+z\right)^{-\zeta}$
(25)
$\displaystyle\times\mathrm{e}^{-z(a_{2}+(a_{1}-a_{2})t^{2})}(1-t^{2})^{\frac{\mu}{2}}t^{\lambda}\mathrm{d}z\mathrm{d}t$
$\displaystyle=2\pi\Theta_{\parallel}^{1+\lambda}\Theta_{\perp}^{2+\mu}\eta_{\parallel}^{\frac{\lambda+1}{2}}\eta_{\perp}^{\frac{\mu+2}{2}}\Gamma\left(\frac{3+\lambda+\mu}{2}\right)$
$\displaystyle\times\int\limits_{0}^{1}U\left(\frac{3+\lambda+\mu}{2},\frac{5+\lambda+\mu}{2}-\zeta,a_{2}+(a_{1}-a_{2})t^{2}\right)$
$\displaystyle\times(1-t^{2})^{\frac{\mu}{2}}t^{\lambda}\mathrm{d}t\,.$

For the normalization constant we thus find ($\lambda=0,\mu=0$):

$\displaystyle N_{\mathrm{GAK}}^{-1}=$
$\displaystyle\mathchoice{{\hbox{$\displaystyle\sqrt{\pi^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\textstyle\sqrt{\pi^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\scriptstyle\sqrt{\pi^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=4.30276pt,depth=-3.44223pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\pi^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=3.44165pt,depth=-2.75334pt}}}\Theta_{\parallel}\Theta_{\perp}^{2}\eta_{\parallel}^{\frac{1}{2}}\eta_{\perp}$
(26)
$\displaystyle\times\int\limits_{0}^{1}U\left(\frac{3}{2},\frac{5}{2}-\zeta,(a_{2}+(a_{1}-a_{2})t^{2}\right)\mathrm{d}t\,.$

By introducing the shorthand notations

$\ {}^{[\lambda,\mu]}\mathcal{W}^{[]}$
$\displaystyle(\eta_{\parallel}.\eta_{\perp},\zeta,\xi_{\parallel},\xi_{\perp})\equiv\mbox{$\
{}^{[\lambda,\mu]}\mathcal{W}^{[]}$}\equiv$ (27)
$\displaystyle\frac{\int\limits_{0}^{1}U\left(\frac{3+\lambda+\mu}{2},\frac{5+\lambda+\mu}{2}-\zeta,a_{2}+(a_{1}-a_{2})t^{2}\right)(1-t^{2})^{\frac{\mu}{2}}t^{\lambda}\mathrm{d}t}{\int\limits_{0}^{1}U\left(\frac{3}{2},\frac{5}{2}-\zeta,a_{2}+(a_{1}-a_{2})t^{2}\right)\mathrm{d}t}$

and

$\displaystyle\mbox{$\
{}^{[]}\mathcal{W}^{[]}$}\equiv\frac{1}{\int\limits_{0}^{1}U\left(\frac{3}{2},\frac{5}{2}-\zeta,a_{2}+(a_{1}-a_{2})t^{2}\right)\mathrm{d}t}\,,$
(28)

we can write the moments as

$\displaystyle M^{\lambda,\mu}$
$\displaystyle=\frac{2}{\mathchoice{{\hbox{$\displaystyle\sqrt{\pi\,}$}\lower
0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\textstyle\sqrt{\pi\,}$}\lower
0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\scriptstyle\sqrt{\pi\,}$}\lower
0.4pt\hbox{\vrule
height=3.01389pt,depth=-2.41113pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\pi\,}$}\lower
0.4pt\hbox{\vrule
height=2.15277pt,depth=-1.72223pt}}}}n_{0}\Theta_{\parallel}^{\lambda}\Theta_{\perp}^{\mu}\eta_{\parallel}^{\frac{\lambda}{2}}\eta_{\perp}^{\frac{\mu}{2}}\Gamma\left(\frac{3+\lambda+\mu}{2}\right)\mbox{$\
{}^{[\lambda,\mu]}\mathcal{W}^{[]}$}\,.$ (29)

The distribution function is then

$\displaystyle f_{\mathrm{GAK}}=$
$\displaystyle\frac{n_{0}}{\mathchoice{{\hbox{$\displaystyle\sqrt{\pi^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\textstyle\sqrt{\pi^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\scriptstyle\sqrt{\pi^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=4.30276pt,depth=-3.44223pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\pi^{3}\,}$}\lower
0.4pt\hbox{\vrule
height=3.44165pt,depth=-2.75334pt}}}\Theta_{\parallel}\Theta_{\perp}^{2}\eta_{\parallel}^{\frac{1}{2}}\eta_{\perp}}\mbox{$\
{}^{[]}\mathcal{W}^{[]}$}$ (30)
$\displaystyle\left(1+\frac{v_{\parallel}^{2}}{\eta_{\parallel}\Theta_{\parallel}^{2}}+\frac{v_{\perp}^{2}}{\eta_{\perp}\Theta_{\perp}^{2}}\right)^{-\zeta}\mathrm{e}^{-\xi_{\parallel}\frac{v_{\parallel}^{2}}{\Theta_{\parallel}^{2}}-\xi_{\perp}\frac{v_{\perp}^{2}}{\Theta_{\perp}^{2}}}\,.$

$\vec{c}(1:9,i)$ |  $n_{i}$ [cm-3] |  $\Theta_{\parallel,i}$ [Mm/s] |  $\Theta_{\perp,i}$ [Mm/s] | $\eta_{\parallel,i}$ | $\eta_{\perp,i}$ | $\zeta$ | $\xi_{\parallel,i}$ | $\xi_{\perp,i}$ |  $u_{\parallel,i}$ [Mm/s]
---|---|---|---|---|---|---|---|---|---
1 | $1\cdot 10^{-1}$ | 1.0 | 1.0 | 1.5 | 1.5 | 2.5 | $1\cdot 10^{-2}$ | $1\cdot 10^{-3}$ | $5\cdot 10^{-2}$
2 | $1\cdot 10^{-3}$ | 1.5 | 1.5 | 1.5 | 1.5 | 2.5 | $1\cdot 10^{-2}$ | $1\cdot 10^{-3}$ | $5\cdot 10^{-2}$
3 | $1\cdot 10^{-7}$ | 2.0 | 2.0 | 1.5 | 1.5 | 2.5 | $1\cdot 10^{-2}$ | $1\cdot 10^{-3}$ | $5\cdot 10^{-2}$
Table 6: The initial parameter set.

##  Appendix B Fitting procedure

###  B.1 Data and normalization

The electron high-resolution data from the Ulysses SWOOPS instrument are given
for 20 energy channels ($E_{kin}$): 1.69, 2.35, 3.25, 4.51, 6.26, 8.65, 12.1,
16.8, 23.2, 31.9, 43.9, 60.6, 84.0, 116, 163, 226, 312, 429, 591, 815 [eV],
and for 20 angle bins $\varphi$ in 9∘ steps in the range from 0∘ to 180∘,
which gives us a 20$\times$20 array.

The first angle bin is along the magnetic field line. For the speed $v$ we use
the classical energy-velocity relation

$\displaystyle\frac{1}{2}m_{e}v^{2}=E_{kin}\,,$ (31)

where $m_{\mathrm{e}}$ is the electron mass, and $v$ is the speed. The speed
can be decomposed into a parallel and a perpendicular component in the
following way:

$\displaystyle v=v_{\parallel}\cos\varphi+v_{\perp}\sin\varphi\,.$ (32)

Thus, with the above angles, $v_{{}_{\parallel}}$ has values in the interval
[$-15,\ldots,15$] Mm, and $v_{{}_{\perp}}$ in the interval [$0,\ldots,15$] Mm.

The distribution function is normalized to [s3/cm3/Mm3]. The nonlinear least-
squares fit is done by a Levenbergh-Marquardt algorithm from the minpack-
package (Elzhov et al. 2016).

We performed a couple of tests, where we defined an initial function from the
intial dataset, i.e., those data which are used as best guess. We fit always
the logarithm of the data and the model distributions, because otherwise only
the core is well fitted. The reason is that the least square fit minimizes the
total distances between the data points and the model. Thus in fitting the
original functions the distances of the smaller values of the distribution do
not contribute much to the total distance. For a logarithmic fit the distances
on all scales contribute in a similar manner to the total distances.

###  B.2 Input parameter

The array $\vec{c}(1:9,1:3)$ which is used as initial array for fitting, e.g.,
the assumed “best” initial values is given in Table 6.


# Semi-Supervised Noisy Student Pre-training on EfficientNet Architectures for
Plant Pathology Classification

Sedrick Scott Keh
The Hong Kong University of Science and Technology

<EMAIL_ADDRESS>
###### Abstract

In recent years, deep learning has vastly improved the identification and
diagnosis of various diseases in plants. In this report, we investigate the
problem of pathology classification using images of a single leaf. We explore
the use of standard benchmark models such as VGG16, ResNet101, and DenseNet
161 to achieve a 0.945 score on the task. Furthermore, we explore the use of
the newer EfficientNet model, improving the accuracy to 0.962. Finally, we
introduce the state-of-the-art idea of semi-supervised Noisy Student training
to the EfficientNet, resulting in significant improvements in both accuracy
and convergence rate. The final ensembled Noisy Student model performs very
well on this particular classification task, achieving a test score of 0.982.

##  1 Introduction

Plant pathology classification is a very challenging and important task. In
2019, more than 20$\%$ of all global crop productions were lost due to pests
and other pathogens. This results in massive economic losses and food deficits
in many parts of the world. Therefore, it is very important to be able to
develop systems that can detect such diseases. Such systems will have a
significant impact on the agricultural sector, especially in areas where
agriculture is a primary food source for millions of people.

Figure 1: Examples of fungal diseases in apples [22]

Currently, plant disease diagnosis is mostly done manually by human
examination. However, this is inefficient, as it is tedious, time-consuming,
and error-prone. Recently, deep learning and computer vision techniques have
shown great promise in this task. These systems have the ability to quickly
and reliably make predictions, and have greatly contributed to plant
classification in the past few years. However, there are still many
challenges. For instance, images of different leaves may look very similar
with each other. Meanwhile, two images from the same leaf species may have
very different images. Plant disease symptoms may also manifest in different
ways due to age of infected tissues, genetic variations, and light conditions.

In this report, we examine the use of pretrained deep learning models in plant
pathology classification. More specifically, we examine the state-of-the-art
EfficientNet [21] and a semi-supervised approach called Noisy Student training
[24]. We demonstrate that the Noisy Student EfficientNet performs very well on
our dataset, and it has significant improvements over previous attempts at the
task.

##  2 Related Work

###  2.1 Traditional Methods

Prior to the introduction of deep learning methods, plant pathology detection
was mostly done through traditional image processing and classification
techniques [12]. Until now, some of these traditional methods remain widely
used due to their simplicity.

####  2.1.1 Image Processing Methods

An image of a leaf can reveal many characteristics of a particular plant, such
as color, style, and shape. Therefore, many early attempts at this task were
focused on extracting features from these leaf characteristics. Using image
processing techniques such as edge detectors to analyze image gradients, one
can easily identify features such as discolorations, spots, or holes on the
leaf surfaces [13].

While this idea makes intuitive sense, it may not work that well in all
situations. This is because certain diseases may not be so easily visible from
the image due to variations in lighting, obstruction, and position.
Conversely, the discoloration seen in an image may not necessarily be from
diseases as well. It is possible that these changes in color are just due to
sunlight patterns or other external factors.

####  2.1.2 Clustering

There have also been various attempts to use clustering algorithms to group
similar leaves together. The biggest advantage of using a clustering algorithm
is that it is sensitive to patterns in local areas [11]. For example, some
diseases may be focused more on damaging the edges of a leaf, while others may
cause discolorations in the center of the leaf. The k-means clustering
algorithm will be able to capture these relationships since it will group
together leaves which have symptoms in the same general area. [10]

However, one disadvantage of clustering is that once again, it may fail to
pick up on certain patterns or be sensitive false patterns that are not
actually there. Because the distance metric is usually based on pixel vales or
other pixel features, it is prone to errors from variations in position,
occlusion, lighting, and many others.

####  2.1.3 SVM and Other Classification Models

To address the issues of the previous methodologies, many recent approaches
have explored the use of machine learning architectures [7]. For instance, a
report by Singh [18] explored the use of genetic algorithms in disease image
classification. In their report, they first performed image processing to find
certain features of the leaves, followed by a genetic optimization algorithm
to reach a classification accuracy of $95.71\%$. Another report by Rumpf [14]
uses SVMs and achieves a similar accuracy of $96\%$. These, accuracies,
however, are only binary predictions, and fail to identify specific diseases
or whether the plant has more than one disease. When predicting multiple
diseases, the SVM accuracy drops significantly to around $86\%$.

###  2.2 Deep Learning Based Methods

More recently, deep neural networks and CNN-based models have greatly aided in
plant disease classification [23], reaching accuracies as high as $97.4\%$ on
some datasets. Furthermore, these also offer the advantage that there is no
need to for manual feature extraction and engineering. Rather, these deep
learning models can take in the pure images as training data and make
predictions from them. This makes the process more scalable.

####  2.2.1 Transfer Learning on Plant Classification

There have been a few studies on using transfer learning for plant disease
classification. However, many of these studies are too highly specific and may
not generalize well. For instance, one study [16] focused specifically on just
rice crops, and achieved a score of around $91.37\%$ with a pre-trianed
AlexNet model.

Additionally, a report by Aravind $\&$ Raja [1] explores the use of transfer
learning models such as VGG16, VGG19, GoogLeNet, ResNet101 and DenseNet201,
achieving an accuracy of $97.3\%$ on four different species of plants.
However, the state-of-the-art has become much better since then, with new
models such as EfficientNet coming out. This means that the accuracy on the
task can likely be improved. Furthermore, the dataset used was also quite
limited. We aim to improve on this accuracy using more novel architectures and
other model training strategies.

##  3 Data

###  3.1 Dataset

The dataset used in this report is the same dataset from the Plant Pathology
FGVC7 workshop of CVPR 2020 [22]. The data was compiled by the Cornell
Initiative for Digital Agriculture (CIDA), and it contains image data of 1820
apple leaves, together with corresponding labels of their diseases. Each leaf
has 4 columns of binary values (“healthy”, “multiple$\\_$diseases”, “rust”,
“scab”). There are also 1820 testing images with no labels.

Four sample images are shown below

(a) Healthy Leaf

(b) ”Rust” Disease

(c) ”Scab” Disease

(d) Both ”Rust” & ”Scab”

Figure 2: Sample leaf images from the dataset

The target output prediction for each image is a vector of $4$ values,
representing the probabilities of the leaf being in each of the possibilities
(“healthy”, “multiple$\\_$diseases”, “rust”, “scab”).

###  3.2 Data Preprocessing

Before training, some augmentations and preprocessing techniques are first
performed on the images. All input images were first reshaped to have
dimensions $(545\times 545\times 3)$.

For the training data, the input image is randomly flipped horizontally with a
probability of $0.5$, as well as flipped vertically with a probability of
$0.5$. Next, there is a random scalar shift and rotation, with rotation limit
$25$ degrees and probability $0.7$. Furthermore, with a $0.5$ probability, we
randomly select one of the following to apply to the images:

  * •

Emboss the image

  * •

Sharpen the image

  * •

Apply blur on the image

Finally, there is a piecewise affine shift with probability $0.5$, followed by
channel-wise normalization. The purpose of these augmentations is to prevent
overfitting and to add variations to make the training task ”harder” so that
the model can learn a better representation.

Figure 3: Preprocessing, before (left) and after (right)

These steps are not done for the validation and testing data. The only
preprocessing done for validation and testing data is channel-wise
normalization.

##  4 Problem Statement

The main problem is as follows: Given a $(545\times 545\times 3)$ image of a
leaf, we want to identify whether or not the plant has a particular disease,
and which specific disease the plant has. Specifically, our output will be
four probabilities of the leaf being in each of our four classes.

###  4.1 Loss Function

During training, the loss function used is cross-entropy. This compares the
similarity of two probability distributions and is given by the formula

$\textbf{Loss}=\sum_{i=1}^{C}p_{i}log(y_{i})$

Here, $C$ is the number of classes ($C=4$ in our case). $p_{i}$ represents the
actual probabilities, and $y_{i}$ represents the predicted probabilities.

###  4.2 Evaluation Metrics

For evaluation, we follow the metric used by the official Kaggle competition,
which is the mean column-wise ROC AUC. In other words, the score is the
average of the individual AUCs of each predicted column. The ROC AUC is
dependent on certain values such as the true positive rate (TPR) and false
positive rate (FPR).

##  5 Methodology

###  5.1 Benchmark Transfer Learning Models

The idea of Transfer Learning is to use the weights of a large model pre-
trained on another dataset and apply them on a novel task. Some of these large
models include AlexNet [9], VGGNet [17], ResNet [3], and more recently,
DenseNet [5] and EfficientNet [21]. Since the early 2010s, one main dataset
used to train and evaluate these large models is the ImageNet dataset [2].
Since then, the models have achieved progressively better performances on
various tasks.

In this report, our main approach will be to begin with a baseline simple CNN
model, followed by using these benchmark models. Eventually, we will progress
to more state-of-the-art models such as EfficientNet and even semi-supervised
Noisy Student traning. We will then compare the performance of various models
and demonstrate that Noisy Student self-training model outperforms previous
state-of-the-art models on this particular dataset.

####  5.1.1 VGGNet

The first model we attempt will be the VGGNet. The VGGNet is basically a stack
of multiple CNN layers with dropout and pooling. Here, convolutional layers
are of size $2\times 2$ or $3\times 3$, which is smaller compared to the
previous AlexNet model. More specifically, we use the VGG16.

Figure 4: VGG16 Architecture

####  5.1.2 ResNet

After the VGGNet, the next model we try will be the ResNet [3]. The ResNet
introduces the new idea of ”skip connections”, as shown in the figure below.
These special connections allow the gradients to be able to propagate backward
more effectively, hence improving the performance.

Figure 5: ResNet Skip Connection

More specifically, we use the ResNet101 model.

####  5.1.3 DenseNet

After the ResNet, the next model we try on the dataset is the DenseNet [5].
Here, the basic building block is called a dense block. It follows a similar
idea to the ResNet, except that in the DenseNet, there are bottleneck
connections across multiple dense blocks. This helps weights and gradients
propagate even better.

Figure 6: Dense Block Architecture

More specifically, we use the DenseNet161 model.

###  5.2 EfficientNet

From the name, EfficientNet [21], released by Google last 2019, is not only
about accuracy, but also about efficiency. It scales the model to perform much
faster with less parameters than other previous SOTA models.

EfficientNets use the idea of model scaling (see Figure 7). In any model,
there are three types of scaling:

  1. 1.

Depth Scaling: This is the most common form of scaling. It basically means
making a model deeper. In most cases, this will increase the accuracy of the
model. However, as the model goes deeper, it may also lead to problems such as
vanishing gradients, so the accuracy gains start to diminish if we go too
deep.

  2. 2.

Width Scaling: This refers to how wide the network is (for example, how many
channels there are in a convolutional layer). The advantage of having a wider
model is that it’s smaller and able to capture fine-grained features. However,
the disadvantage is that the accuracy will saturate quickly.

  3. 3.

Resolution Scaling: This simply refers to the resolution of the input iamge.
Inuitively, higher images will be better for the model, but practically, in
high resolutions, there is no significant difference between, say, $500\times
500$ as opposed to $560\times 560$.

In conclusion, we see that depth, width, or resolution scaling may improve
performance, but the accuracy gains eventually diminish for larger models.

Figure 7: Combined Scaling for EfficientNet

####  5.2.1 Combined Scaling

EfficientNet uses the idea of combined scaling, which combines the three types
of scaling to achieve optimal performance. It achieves this balance by
considering the following constraints:

$\text{depth: }d=\alpha^{\phi}$ $\text{width: }w=\beta^{\phi}$
$\text{resolution: }r=\gamma^{\phi}$ $\text{s.t
}\quad\alpha\cdot\beta^{2}\cdot\gamma^{2}\approx 2$ $\alpha\geq 1,\beta\geq
1,\gamma\geq 1$

where $\phi$ is a user-specified parameter. These constraints are motivated by
the computational idea that the FLOPS of a convolution operation is usually
proportional to $d\cdot w^{2}\cdot r^{2}$, so implementing these constraints
effectively allows us to control the number of FLOPS.

####  5.2.2 Network Structure

For the building block of the EfficientNet, it actually uses MBConv layers,
which are Inverted Residual Blocks that were introduced in MobileNet V2 [15].
There are also some Squeeze and Excitation blocks [4]. To find the exact
architecture, the authors employed a Neural Architecture Search (NAS) to find
the optimal values of parameters $\alpha,\beta,\gamma,\phi$. Once these
parameters are found, the full EfficientNet model can then be constructed.

There are multiple configurations of EfficientNet. In this report, we use
Efficientnet B5.

###  5.3 Noisy Student Training

Released in November 2019, Noisy Student training [24] is the latest state-of-
the-art in ImageNet classification. It is a semi-supervised self-training
model that uses both labeled and unlabeled data. For the ImageNet task, the
unlabeled dataset used is the JFT-300M datset [20]. One can clearly observe
the massive scale of this unlabeled dataset (300M as compared to ImageNet’s
14M). To handle the unlabeled data, the algorithm trains two types of
networks: a teacher network and a student network.

####  5.3.1 Teacher-Student Architecture

The basic idea is to train a teacher model on labeled images to predict a
distribution of class labels, then applying this to unlabeled data to create
pseudo-labels. Then, a student model (usually larger) will be trained on both
originally labeled data and the new pseudo-labeled data. This is done
iteratively, meaning that the student model will then be used as a new teacher
model to re-predict the pseudo-labels, and a new student model will be
trained, and so on. Both student and teacher models are EfficientNet models,
and the new student models usually become progressively larger as they scale
up through the iterations.

Figure 8: Noisy Student Pipeline

####  5.3.2 Adding Noise to the Student Model

The authors of the Noisy Student paper also discovered that adding a noise to
the student networks was very significant to improving the performance. One
particular useful form of noise in this model is Stochastic Depth [6], which
means adding random shortcuts between layers. Other forms of noise include
data augmentation (rotation, horizontal and vertical shifts), and dropout [19]
(randomly discarding certain neurons). The main idea is to make the learning
more difficult for the student model, which will allow it to learn more
meaningful representations than the teacher, and therefore the model becomes
more robust.

##  6 Experiments

In this report, we demonstrate the different experiments done on the various
transfer learning models outlined in Section 5. Because these models were
quite large, training was done with a batch size of 4 with the Adam optimizer
[8] and a learning rate of 1e-3. Training was done for 30 epochs with a
learning rate decay of 1e-3, and the loss and accuracy were measured after
every epoch.

###  6.1 Benchmark Models

We first begin with a basic three-layer convolutional network. This is
followed by our three benchmark transfer learning models, namely VGGNet,
ResNet, and DenseNet. These models were pre-trained, and a linear classifier
was added at the end. The results of the experiments are shown below:

Model | Test ROC AUC |  |
---|---|---|---
3 Layer CNN | 0.615 |  |
VGG16 | 0.933 |  |
ResNet101 | 0.943 |  |
DenseNet161 | 0.945 |  |
Table 1: Test ROC AUC for the different models

From the table above, we see that the three-layer CNN performs very poorly on
the dataset, while the other three transfer learning models perform quite well
on this dataset, with VGG performing significantly better than the simple
three-layer CNN, then ResNet performing slightly better than VGGNet, and
finally DenseNet achieving an ROC AUC score of around $0.945$. However, while
these may seem quite high, it is actually not that high for this specific
task. The DenseNet score of $0.945$ ranks around 850 out of 1345 in the Kaggle
public leaderboard. We see later on that our EfficientNet model greatly
outperforms this.

In order to further examine the DenseNet model and its common failures, we
take the DenseNet and plot its confusion matrix. (See Figure 12 in next page.)

Figure 9: DenseNet Confusion Matrix

From this, confusion matrix, we see that most of the errors were in the second
row, i.e. when the leaf class was ”multiple diseases”. Only 6 images with
multiple diseases were classified correctly, while 12 were classified
incorrectly. These errors are consistent with our intuition, since we would
expect that plants with multiple diseases would be harder to classify.

###  6.2 EfficientNet Architectures

EfficientNet architectures perform significantly better than the previous
benchmark models. The baseline EfficientNet achieves a leaderboard test
accuracy of $0.962$. Meanwhile, using weights from the Noisy Student Training
achieves a leaderboard test accuracy of $0.967$. This is indeed a big jump
compared to DenseNet at $0.945$

Model | Test ROC AUC |  |
---|---|---|---
3 Layer CNN | 0.615 |  |
VGG16 | 0.933 |  |
ResNet101 | 0.943 |  |
DenseNet161 | 0.945 |  |
EfficientNet | 0.962 |  |
EfficientNet w/ Noisy Student | 0.967 |  |
Table 2: Summary of Final Results

These observations are consistent with what we expect based on their
performance on ImageNet and other similar tasks. The EfficientNet with Noisy
Student performance is actually very good on this task. With an ROC AUC
accuracy of $0.967$, this already outperforms most of the models in
literature, and it achieves this without extensive feature engineering or
model ensembling.

####  6.2.1 Loss Curve Analysis (Figure 11)

The improvement of EfficientNet models and Noisy Student Models can be seen
most evidently in the loss and accuracy curves plotted at Figure 11 (see next
page). In Figure 11 in the next page, we have two plots, one for training loss
across epochs, and another for validation ROC AUC accuracies.

In the curve of training losses, it is evident that EfficientNet and Noisy
Student have consistently lower losses than the other models. These
EfficientNet architectures also converge faster and are less prone to
fluctuations. Comparing the baseline EfficientNet with the Noisy Student
Training, we can also notice that using Noisy Student shows a small but clear
improvement.

A similar pattern is shown in the curves for validation accuracies across
epochs. We see that as they converge, the EfficientNet and Noisy Student
models have a much higher validation accuracy as compared to VGGNet, ResNet,
and DenseNet. Here, we also notice that there are many fluctuations in the
early epochs for EfficientNet (and ResNet). However, once we add the Noisy
Student weights, it becomes much more stable and converges faster, converging
at around epoch 18 as opposed to the baseline EfficientNet, which conveges at
around epoch 23.

####  6.2.2 Confusion Matrix

To further examine the performance of the EfficientNet with Noisy Student
training, we visualize the confusion matrix, just as we did with the DenseNet.

Figure 10: Noisy Student Training Confusion Matrix

Indeed, we can see the confusion matrix of the Noisy Student no longer suffers
from the same shortcomings of the DenseNet. The plants with multiple diseases
are classified accurately.

Figure 11: Train Loss (top) and Validation Accuracy (bottom)

###  6.3 Hyperparameter Tuning (For Noisy Student)

Since we are using pre-trained models, there are not that many hyperparameters
to tune. Furthermore, we are forced to use small batches due to GPU memory
restrictions. Here, we focus on two main ones: learning rate and optimizer.

####  6.3.1 Optimizer

For the optimizers, we compare Adam and SGD

Optimizer | Test ROC AUC |  |
---|---|---|---
SGD | 0.9669 |  |
Adam | 0.9672 |  |
Table 3: Test ROC AUC for the different models

From Table 3, we see that the final results of Adam and SGD were more or less
the same. However, we noticed during training that Adam converged faster, so
we select Adam as our optimizer.

####  6.3.2 Learning Rate

Below are the performances of different learning rates.

Learning Rate | Test ROC AUC |  |
---|---|---|---
1e-1 | 0.9437 |  |
1e-2 | 0.9653 |  |
1e-3 | 0.9672 |  |
1e-4 | 0.9667 |  |
Table 4: Performance of different learning rates

We see in Table 4 that the ROC AUC for various learning rates are more or less
the same (except for 1e-1, which is much lower). However, the performance for
1e-3 is slightly better than the rest, so that is the one we select.

###  6.4 Additional Techniques

####  6.4.1 Model Ensembling

Model Ensembling is the idea of combining the losses of multiple models into
one so that the models can compensate for the shortcomings of the other
models. Here, I combined my best EfficientNet model with my two best Noisy
Student models.

Figure 12: Model Ensembling

The final ensembled model had a test score of 0.982 and ranked 64 out of 1345
(top 5$\%$) on the public leaderboard.

##  7 Conclusion

In this report, we applied various CNN-based image classification transfer
learning models to the novel Plant Pathology Classification task. We
discovered that the EfficientNet model outperforms the other models
significantly. Additionally, we see that adding Noisy Student Training to the
EfficientNet model further improves the model performance and makes the model
converge faster and less prone to fluctuations.

##  8 Future Work

The field of plant pathology classification is a very significant one with
large global impacts. This paper showed the use of EfficientNet architectures
on apple leaves. In the future, more tests can be done on various other types
of plants and fruits. Since the accuracy for this task is already quite high,
one possible area of improvement would be the speed of detection and
classification (real-time detection, mobile devices, etc.)

## References

  * [1] K. R. Aravind and P. Raja.  Automated disease classification in (selected) agricultural crops using transfer learning.  Automatika, 61(2):260–272, 2020.
  * [2] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei.  ImageNet: A Large-Scale Hierarchical Image Database.  In CVPR09, 2009.
  * [3] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  CoRR, abs/1512.03385, 2015.
  * [4] J. Hu, L. Shen, and G. Sun.  Squeeze-and-excitation networks.  CoRR, abs/1709.01507, 2017.
  * [5] G. Huang, Z. Liu, and K. Q. Weinberger.  Densely connected convolutional networks.  CoRR, abs/1608.06993, 2016.
  * [6] G. Huang, Y. Sun, Z. Liu, D. Sedra, and K. Q. Weinberger.  Deep networks with stochastic depth.  CoRR, abs/1603.09382, 2016.
  * [7] S. D. Khirade and A. B. Patil.  Plant disease detection using image processing.  pages 768–771, 2015.
  * [8] D. Kingma and J. Ba.  Adam: A method for stochastic optimization.  International Conference on Learning Representations, 12 2014.
  * [9] A. Krizhevsky, I. Sutskever, and G. E. Hinton.  Imagenet classification with deep convolutional neural networks.  In Proceedings of the 25th International Conference on Neural Information Processing Systems - Volume 1, NIPS’12, page 1097–1105, Red Hook, NY, USA, 2012. Curran Associates Inc.
  * [10] C. U. Kumari, S. Jeevan Prasad, and G. Mounika.  Leaf disease detection: Feature extraction with k-means clustering and classification with ann.  In 2019 3rd International Conference on Computing Methodologies and Communication (ICCMC), pages 1095–1098, 2019.
  * [11] S. Maity, S. Sarkar, A. Tapadar, A. Dutta, S. Biswas, S. Nayek, and P. Saha.  Fault area detection in leaf diseases using k-means clustering.  10 2018.
  * [12] A. Patel and B. Joshi.  A survey on the plant leaf disease detection techniques.  IJARCCE, 6:229–231, 01 2017.
  * [13] S. Radha.  Leaf disease detection using image processing.  Journal of Chemical and Pharmaceutical Sciences, 03 2017.
  * [14] T. Rumpf, A.-K. Mahlein, U. Steiner, E.-C. Oerke, H.-W. Dehne, and L. Plümer.  Early detection and classification of plant diseases with support vector machines based on hyperspectral reflectance.  Computers and Electronics in Agriculture, 74:91–99, 10 2010.
  * [15] M. Sandler, A. G. Howard, M. Zhu, A. Zhmoginov, and L. Chen.  Inverted residuals and linear bottlenecks: Mobile networks for classification, detection and segmentation.  CoRR, abs/1801.04381, 2018.
  * [16] V. K. Shrivastava, M. K. Pradhan, S. Minz, and M. P. Thakur.  Rice Plant Disease Classification Using Transfer Learning of Deep Convolution Neural Network.  ISPRS - International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 423:631–635, July 2019.
  * [17] K. Simonyan and A. Zisserman.  Very deep convolutional networks for large-scale image recognition.  CoRR, abs/1409.1556, 2014.
  * [18] V. Singh and A. Misra.  Detection of plant leaf diseases using image segmentation and soft computing techniques.  Information Processing in Agriculture, 4(1):41 – 49, 2017.
  * [19] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov.  Dropout: A simple way to prevent neural networks from overfitting.  J. Mach. Learn. Res., 15(1):1929–1958, Jan. 2014.
  * [20] C. Sun, A. Shrivastava, S. Singh, and A. Gupta.  Revisiting unreasonable effectiveness of data in deep learning era.  CoRR, abs/1707.02968, 2017.
  * [21] M. Tan and Q. V. Le.  Efficientnet: Rethinking model scaling for convolutional neural networks, 2019.  cite arxiv:1905.11946Comment: Published in ICML 2019.
  * [22] R. Thapa, N. Snavely, S. Belongie, and A. Khan.  The plant pathology 2020 challenge dataset to classify foliar disease of apples, 2020.
  * [23] Y. Toda and F. Okura.  How convolutional neural networks diagnose plant disease.  Plant Phenomics, 2019.
  * [24] Q. Xie, E. Hovy, M.-T. Luong, and Q. Le.  Self-training with noisy student improves imagenet classification.  11 2019.