11institutetext: Institute of Frontier and Interdisciplinary Science, Shandong
University, Qingdao, Shandong, 266237, People’s Republic of China
11email<EMAIL_ADDRESS>22institutetext: Institute of Space Sciences,
Shandong University, Shandong, 264209, People’s Republic of China
33institutetext: Department of Earth and Space Sciences, Southern University
of Science and Technology, Shenzhen, Guangdong, People’s Republic of China
44institutetext: Institute of Solar-Terrestrial Physics, Irkutsk, 664033,
Russia

# Excitation of extraordinary modes inside the source of Saturn’s kilometric
radiation

Hao Ning 1122 Yao Chen 1122 Chuanyang Li 1122 Shengyi Ye 33 Alexey Kuznetsov
44 Siyuan Wu 33

(Received ; accepted )

The electron cyclotron maser instability (ECMI) of extraordinary mode waves
was investigated with the parameters observed in Saturn’s kilometric radiation
(SKR) sources. Previous studies employed simplified dispersion relations, and
did not consider the excitation of the relativistic (R) mode. This mode is
introduced by considering the relativistic effect in plasmas consisting of
both cold and hot electrons. Using particle-in-cell simulations, we
investigated the excitation of R and X modes based on the measured data. Using
the reported value of the density ratio of energetic to total electrons
$n_{e}/n_{0}=24\%$, the most unstable mode is the R mode. The escaping X-mode
emissions are amplified only if the energetic electrons are dominant with
$n_{e}/n_{0}\geq 90\%$. For these cases, only the X mode is excited and the R
mode disappears due to its strong coupling. The results are well in line with
the linear kinetic theory of ECMI. The properties of both the R and X modes
are consistent with the observed SKR emissions. This raises questions about
the nature of the measured electric field fluctuations within “presumed” SKR
sources. The study provides new insights into the ECMI process relevant to SKR
emission mechanisms.

###### Key Words.:

Planets and satellites: gaseous planets – Radio continuum: planetary systems –
Masers – Waves – Plasmas – Methods: numerical

##  1 Introduction

Saturn’s kilometric radiations (SKRs) represent intense radio emissions
originating from the kronian auroral region, considered to be the counterpart
of the auroral kilometric radiations (AKRs) on Earth (see, e.g., Kaiser et
al., 1980, 1981; Gurnett et al., 2004; Kurth et al., 2005; Lamy et al., 2009).
The radiations are observed in the kilohertz to megahertz frequency range, and
the radiated powers of $10^{8}$–$10^{10}$ watts. Furthermore, SKRs are closely
related to Saturn’s magnetosphere dynamics (Jackman, 2011; Kimura et al.,
2013), and have been extensively studied.

Saturn’s kilometric radiations are often strongly polarized with narrow
beaming angles and high intensity (Lamy et al., 2008, 2018), indicating
coherent emissions generated by electron cyclotron maser instability (ECMI)
(Wu & Lee, 1979; Wu, 1985; Treumann, 2006; Bingham et al., 2013). The
instability excites escaping modes at gyro-frequencies ($\sim\Omega_{ce}$) via
a gyro-resonant wave-particle interaction, driven by electrons with
anisotropic velocity distribution functions (VDFs; see Ning et al., 2021a, b,
for recent studies). The ECMI is now widely applied to planetary and stellar
radio emissions.

The Cassini spacecraft measures the radio waves with the Radio and Plasma Wave
Science (RPWS, Gurnett et al., 2004) instrument and the electrons with the
Cassini Plasma Spectrometer Electron Spectrometer (CAPS/ELS, Linder et al.,
1998; Young et al., 2004) in Saturn’s magnetosphere. Since 2010, researchers
have reported several crossings of SKR sources by Cassini, identified with the
cutoff frequency of the spectrum close to or slightly lower than the local
$\Omega_{ce}$ (e.g., Kurth et al., 2011; Mutel et al., 2010; Lamy et al.,
2010; Menietti et al., 2011; Lamy et al., 2018).

As has previously been reported, SKR sources are characterized by a low ratio
of plasma oscillation frequency to electron gyro-frequency
($\omega_{pe}/\Omega_{ce}<0.1$), in accordance with the ECMI conditions
($\omega_{pe}/\Omega_{ce}<1$, Wu & Lee, 1979). The emissions mainly propagate
quasi-perpendicularly to the background magnetic field with linear
polarization near the sources, consistent with the X-mode properties. The
energy conversion rate from electrons to radiations is $\sim 1\%$ (Lamy et
al., 2011). The abundances of energetic and cold electrons are comparable with
each other (Schippers et al., 2011).

The measurement of electron VDFs is crucial to the study of the ECMI process.
In the first encounter of the SKR source by Cassini in 2008, Mutel et al.
(2010) found that the cold and hot electrons can be fitted with a Kappa and a
Dory-Guest-Harri (DGH, Dory et al., 1965) distribution, respectively. They
further calculated the corresponding ECMI growth rates of X-mode emissions
with the linear kinetic theory, and suggested that the ECMI process can
efficiently amplify this escaping X mode below $\Omega_{ce}$, along directions
perpendicular to the magnetic field. Menietti et al. (2011) investigated the
ECMI process of another SKR source using a similar method and obtained the
same results. They also suggested that O-mode and Z-mode waves can be
generated.

However, these studies on the ECMI of SKRs followed Mutel et al. (2007) and
employed the dispersion relations of pure relativistic energetic electrons.
They estimated the correction to the dispersion relation with a scale factor
when calculating the growth rates of the X and Z modes, rather than solving
the dispersion matrix to calculate the real frequency and the growth rate in a
self-consistent manner. As reported by Mutel et al. (2010) and Menietti et al.
(2011), the density ratio of the energetic to the total electrons was fitted
as $n_{e}/n_{0}=20\%$–24%. With this range of $n_{e}/n_{0}$ (and other
comparable values), one obtains very different dispersion relations. For waves
perpendicular to the magnetic field, the gyro-resonance condition is
simplified as $\omega=n\Omega_{ce}/\gamma<n\Omega_{ce}$. When adopting the
cold-plasma dispersion relation deduced from the magneto-ionic theory, the
ECMI process cannot amplify the fundamental X mode with the cutoff frequency
($\omega_{\text{X}}$) being larger than $\Omega_{ce}$. For pure energetic
electrons, the X-mode cutoff is less than $\Omega_{ce}$ and thus can grow via
ECMI (e.g., Winglee, 1983; Pritchett, 1984a; Winglee, 1985; Robinson, 1986).

For plasmas consisting of both the cold and energetic electrons, another
important mode may emerge. This mode, first deduced by Pritchett (1984b),
appears in the range of $\Omega_{ce}/\gamma<\omega<\Omega_{ce}$. It is
relevant to the relativistic effect, and is therefore called the relativistic
(R) mode (Strangeway, 1985, 1986). Our present study shows that the R mode
becomes the most unstable mode, and is probably relevant to SKR. However, it
is trapped and cannot escape. This mode has not been considered in studies on
SKR sources, to the best of our knowledge.

According to the earlier reports on the crossings of SKR sources by Cassini,
both cold and hot components of electrons exist with comparable densities
(Schippers et al., 2011). This raises questions concerning the applicability
of the oversimplified dispersion relations adopted by Mutel et al. (2010) and
Menietti et al. (2011), which did not consider the existence of the R mode.
For this study, we conducted fully kinetic electromagnetic particle-in-cell
(PIC) simulations with the electron VDFs and plasma parameters reported by
Menietti et al. (2011) to investigate the extraordinary modes (including R, X,
and Z) and their coupling within SKR sources. We conducted a parametric study
on $n_{e}/n_{0}$, and applied the kinetic linear theory to verify and explain
the PIC simulations.

##  2 Numerical method and parameter setting

As reported by Menietti et al. (2011), the frequency ratio was
$\omega_{pe}/\Omega_{ce}=0.046$, and the VDFs of both the cold and hot
electrons were written as

$f_{0}=\frac{1}{(2\pi)^{3/2}v_{0}^{3}}\exp(-\frac{u^{2}}{2v_{0}^{2}})$ (1)

and

$f_{e}=\frac{1}{j!(2\pi)^{3}/2v_{\perp}^{2}v_{\parallel}}\left(\frac{u_{\perp}^{2}}{2v_{\perp}^{2}}\right)^{j}\exp\left(-\frac{u_{\perp}^{2}}{2v_{\perp}^{2}}-\frac{u_{\parallel}^{2}}{2v_{\parallel}^{2}}\right),$
(2)

respectively, where $j=11$, $u$ represents the momentum per mass, and
$v_{\parallel}\leavevmode\nobreak\ (=0.11c)$ and
$v_{\perp}\leavevmode\nobreak\ (=0.038c)$ represent the thermal velocity of
the DGH distribution along parallel and perpendicular directions. The density
ratio of the energetic to the total electrons was fitted to be
$n_{e}/n_{0}=24\%$. Considering the uncertainties of the source measurements,
we conducted a parametric study of $n_{e}/n_{0}$. The models and parameter
settings are presented below.

###  2.1 PIC simulations

We employed the Vector-PIC (vPIC) code developed by the Los Alamos National
Laboratory (LANL). The code employs a second-order, explicit, leapfrog
algorithm to update the positions and the velocities of charged particles,
along with a full Maxwell description of the electric and magnetic fields via
a second-order finite-difference time-domain solver (Bowers et al., 2008,
2009).

We performed 2D3V (two spatial dimensions with three velocity components)
simulations in the $xOz$ plane with periodical boundary conditions and a
uniform background magnetic field ($\@vec{B}_{0}=B_{0}\hat{e}_{z}$). The
simulation domain was set to $L_{x}=L_{z}=2048\Delta$, the grid length to
$\Delta=0.019\leavevmode\nobreak\ \lambda_{DE}$, and the time step to $\Delta
t=0.0097\leavevmode\nobreak\ \omega_{pe}^{-1}$.

Two components of electrons were included. Their distribution functions are
described by Eqs. (1) and (2), and the center of the VDF of energetic
electrons is at $v_{\perp}\sim 0.18\leavevmode\nobreak\ c$ (Fig. 1(a)). The
thermal velocity of the low-energy electrons was set to $v_{0}=0.01c$ ($\sim
50$ eV). The charge neutrality was maintained by including thermal protons
with a realistic proton-to-electron mass ratio of 1836 and the same
temperature (with the thermal velocity $v_{p}=v_{0}/\sqrt{1836}$). We included
1000 macroparticles per cell for each component of the electrons and 500 for
the protons.

Figure 1: Evolution of electron VDFs and energy profiles of the simulation
with $n_{e}/n_{0}=24\%$. Panels (a–c): Snapshots of the VDFs at the start (a,
$t=0\leavevmode\nobreak\ \omega_{pe}^{-1}$), middle (b, 100
$\omega_{pe}^{-1}$), and end (c, 200 $\omega_{pe}^{-1}$) of the simulation.
Panel (d): Temporal variations of energies of the six fluctuated field
components. The black line represents the decline of the electron kinetic
energy ($-\Delta E_{k}$). The energies are normalized to the initial kinetic
energy of the total electrons. The gray dashed line represents the exponential
fitting of the energy profiles. The vertical dotted line in (d) indicates the
separate time (120 $\omega_{pe}^{-1}$) of the linear stage and the saturation
stage. Figure 2: Wave dispersion diagrams of electric and magnetic fields with
$n_{e}/n_{0}=24\%$ and $\theta_{kB}=90^{\circ}$, during the interval of [100,
200] $\omega_{pe}^{-1}$. The six field components are provided in panels (a–f)
and the separated LH and RH components in panels (g–h).

###  2.2 Linear analysis of ECMI

Following the analysis of Pritchett (1984b), we employed the $\delta$ function
below to describe the two components of electrons,

$f=(1-\frac{n_{e}}{n_{0}})\delta(u_{\perp})\delta(u_{\parallel})+\frac{n_{e}}{n_{0}}\delta(u_{\perp}-u_{r})\delta(u_{\parallel}),$
(3)

where $u_{r}$ refers to the momentum of the energetic electrons. In accordance
with the PIC simulation, we assumed the waves propagate in the $xOz$ plane
with $\@vec{B}_{0}=B_{z}\hat{e}_{z}$. For simplicity, only the perpendicular
propagating waves with $k_{\parallel}=0$ were studied.

With the VDF presented above, we derived the elements of the plasma dispersion
tensor (Eq. (5)) and obtained the frequencies and the growth rates of the ECMI
modes by numerically solving their determinant. We analyzed the dispersion
relations, the ECMI growth rates, and the polarization of the X and R modes.
The results are presented in the appendix.

##  3 PIC simulations

We present the PIC simulations with the reported SKR source parameters as the
reference case with ($n_{e}/n_{0}=24\%$) to show the excitation and properties
of extraordinary modes (R, X, Z), followed by the parametric study on
$n_{e}/n_{0}$ to explore its effect.

###  3.1 Results with measured parameters ($n_{e}/n_{0}=0.24$)

According to the energy profiles presented in Fig. 1(d), the wave excitation
and further evolution can be divided into two stages, the linear (0–120
$\omega_{pe}^{-1}$) stage and the saturation (120–200 $\omega_{pe}^{-1}$)
stage. In the linear stage, the energies of $\Delta B_{z}\leavevmode\nobreak\
(=B_{z}-B_{0})$ and $E_{y}$ rise rapidly, reaching the level of
$10^{-2}\leavevmode\nobreak\ E_{k0}$ ($E_{k0}$ is the initial kinetic energy
of total electrons). The electron kinetic energy ($E_{k}$) declines in
accordance with the growth of the ECMI modes. The energies of other field
components are much weaker, being $\sim 10^{-5}\leavevmode\nobreak\ E_{k0}$.
In the saturation stage, the energies of all the fields rise gradually. In
total, $\sim$3.4% of $E_{k0}$ was converted to wave energies.

Figure 1(a–c) presents the evolution of the electron VDFs. The two components
are initially separated. Within 100–200 $\omega_{pe}^{-1}$, the shell-like hot
electrons diffuse in the phase space, later developing a butterfly-like
morphology. Meanwhile, the cool Maxwellian electrons are significantly
energized along the perpendicular direction.

We performed a fast Fourier transform (FFT) analysis on each field component.
In Fig. 2(a–f), we show the $\omega-k$ dispersion diagrams for perpendicular-
propagating waves with $\theta_{kB}=90^{\circ}$ (the angle between $\@vec{k}$
and $\@vec{B}_{0}$), along which the ECMI waves are the most amplified. The
strongest wave appears in the $E_{y}$ and $B_{z}$ components, with frequencies
of $\sim 0.98\leavevmode\nobreak\ \Omega_{ce}$, in accordance with the
dispersion relation of the R mode as shown in the Appendix. This mode can be
identified as the transverse electromagnetic wave with significant $E_{y}$ and
$B_{z}$ components and an insignificant $E_{x}$ component. The property of the
simulated R mode is in line with the linear kinetic theory (see Section A.2).

To separate the left-hand (LH) and right-hand (RH) polarized components of the
R mode, we processed ($E_{x}$, $E_{y}$) with the following equation:

$E_{LH,RH}=\frac{1}{\sqrt{2}}(E_{x}\pm iE_{y}).$ (4)

The results are presented in Fig. 2(g–h). The intensities of the LH and RH
components are comparable to each other, indicating a strong linear
polarization.

According to the energy profiles in Fig. 1(d), the total energy of the R mode
is $\sim$0.034 $E_{k0}$, and its growth rate can be fitted as
$\Gamma=5.57\times 10^{-3}\leavevmode\nobreak\ \Omega_{ce}$. The fundamental O
mode is slightly amplified in $E_{z}$ and $B_{y}$. Its total energy is weaker
than the R mode by approximately two orders of magnitude.

In this case, only the R mode can be amplified, with no significant growth of
X. The results can be clearly understood with the linear kinetic theory. As
shown in the Appendix, the R mode is the most unstable at $k\sim\Omega_{ce}/c$
with a maximal growth rate of $0.006\leavevmode\nobreak\ \Omega_{ce}$, and
linearly polarized in $E_{y}$ with $E_{x}/iE_{y}\sim 0$ (Fig. 7(a–c)). The
growth rates of the X mode remain at 0 with its cutoff being larger than
$\Omega_{ce}$.

###  3.2 Effect of $n_{e}/n_{0}$

Figure 3: Energy profiles and growth rates for different simulation cases.
Panels (a–e): Energy profiles of the six fluctuated field components for the
cases with $n_{e}/n_{0}=$ 0.1, 0.5, 0.7, 0.9, and 1.0. The energies are
normalized to the respective initial kinetic energy of total electrons of each
case. Panel (f) presents the growth rates with different density ratios
obtained by the PIC simulations (blue) and the linear kinetic theory (orange).
Figure 4: Snapshots of VDFs of simulations with $n_{e}/n_{0}=$ 10% (a), 50%
(b), 70% (c), 90% (d), and 100% (e) at moments representative of the initial
(left), the transition (middle), and the saturation (right) stages.

Considering the uncertainties of the electron measurement of Cassini CAPS/ELS,
we varied the values of $n_{e}/n_{0}$ from 10% to 100% to investigate the
conditions for X-mode excitations. For cases with a $n_{e}/n_{0}>10\%$, the
simulations lasted for 200 $\omega_{pe}^{-1}$, and for the $n_{e}/n_{0}=10\%$
case, the simulation lasted for 300 $\omega_{pe}^{-1}$.

In Fig. 3(a–e), we plotted the energy profiles of each case. As introduced in
Section 3.1, we divided the whole ECMI process into two stages, the linear
stage and the saturation stage. For each case we observed a sharp rise of
$E_{y}$ in the linear stage. The energy conversion rates from the electron
kinetic energy to the waves is $\sim 10^{-2}$. With a higher $n_{e}/n_{0}$,
the wave excitation is more impulsive, with a short linear stage. For cases
with $n_{e}/n_{0}\geq 50\%$, we observed a gradual decline of the field
energies and $-\Delta E_{k}\leavevmode\nobreak\ (=E_{k0}-E_{k})$, suggesting
further damping of the amplified waves. Figure 4 displays the evolution of the
VDFs in different cases. With a larger $n_{e}/n_{0}$, the VDF of the hot
electrons diffuses faster, and the cold electrons are more efficiently
accelerated. The VDFs do not change much after 50–100 $\omega_{pe}^{-1}$.

Figure 5: $\omega-k$ dispersion diagrams of $E_{y}$ at
$\theta_{kB}=90^{\circ}$ during the saturation stage for different cases.
Figure 6: Wave dispersion diagrams of the six field components (a–f) and the
separated LH and RH components (g–h) with $n_{e}/n_{0}=100\%$ and
$\theta_{kB}=90^{\circ}$ during [100, 200] $\omega_{pe}^{-1}$.

For each case, we performed a FFT analysis of the fluctuated fields to show
the amplified wave modes (Fig. 5). The waves in $E_{x}$ and $E_{z}$ are not
shown here since they are very weak. For cases with $n_{e}/n_{0}\leq 70\%$,
the R mode is the most intensive, while no significant X mode appears. With
increasing $n_{e}/n_{0}$, the R-mode dispersion curve becomes flatter and
closer to that of the X mode. With $n_{e}/n_{0}\geq 90\%$, the two dispersion
curves merge, indicating their strong coupling; this leads to the significant
excitation of the X mode with a lower cutoff frequency of
$\omega_{\text{X}}\sim 0.97\leavevmode\nobreak\ \Omega_{ce}$.

We plotted the growth rates ($\Gamma$) of the simulated R mode (or the X mode
for the cases with $n_{e}/n_{0}=90$ and 100%) by fitting the energy curves
(Fig. 3(f)). The values of $\Gamma$ rise from 0.0032 to 0.018 $\Omega_{ce}$
with $n_{e}/n_{0}$ from 10% to 100%. With the linear kinetic theory (see Fig.
8 of the Appendix), we obtained comparable growth rates for the R (for the
cases with $0<n_{e}/n_{0}<100\%$) and the X mode (for the case with
$n_{e}/n_{0}=100\%$), which rise from $4.59\times 10^{-3}$ to $1.14\times
10^{-2}\leavevmode\nobreak\ \Omega_{ce}$.

According to our PIC simulations, the X mode is excited only if the energetic
component is dominant ($n_{e}/n_{0}\geq 90\%$). The characteristics of the X
mode in the case of $n_{e}/n_{0}=100\%$, including the growth rate,
polarization, and wave fields (see Figs. 3(e), 4(e), and 6), are consistent
with the linear kinetic theory (Fig. 7(d–f)). In this case, the X mode becomes
saturated within the first 50 $\omega_{pe}^{-1}$, mainly in $E_{y}$ and
$B_{z}$, and much less so in $E_{x}$. The energy conversion rate from
energetic electrons to the X mode is $\sim$3%, and the growth rate is
$\sim$0.018 $\Omega_{ce}$, similar to that obtained by the linear kinetic
theory ($\sim 0.014\leavevmode\nobreak\ \Omega_{ce}$; see Fig. 7(e)).
According to Fig. 6(g–h), the amplified X-mode waves are linearly polarized
with comparable LH and RH intensities.

With the linear kinetic theory (see Section A.3), for the case with
$n_{e}/n_{0}=90\%$, the growth rate of the X mode is 0. This is inconsistent
with the corresponding PIC simulation presented above. This can be explained
by the evolution of the VDFs (Fig. 4(d)). Almost all of the cool electrons are
energized in the simulation within the first 50 $\omega_{pe}^{-1}$, and their
VDFs overlap with the hot component. This makes it difficult to distinguish
the two components in the phase space, leading to an equivalent case with
$n_{e}/n_{0}=100\%$. This explains why the simulation with $n_{e}/n_{0}=90\%$
is in line with that of $n_{e}/n_{0}=100\%$. In other cases with less
energetic electrons, the two components of the VDF are still separated until
the end of the simulation, and the X mode cannot grow. Further discussion is
provided in the following section.

##  4 Conclusion and discussion

The present study investigates the excitation of the R and X modes based on
the in situ observation of SKR sources reported by Menietti et al. (2011).
Both PIC simulations and linear theory analysis were performed, and the
results are in line with each other. With the given set of parameters, the
trapped R mode is amplified as the dominant mode, while the escaping X mode
can not grow. The growth rate of the R mode increases with a $n_{e}/n_{0}$ of
10%–70%. The X-mode emissions can only be amplified if the electrons are
dominated by hot components ($n_{e}/n_{0}=90\%$ and 100%) with the cutoff
frequency less than $\Omega_{ce}$ and the growth rate $\sim
0.018\leavevmode\nobreak\ \Omega_{ce}$. According to our PIC simulations, both
the X and R modes can be characterized with similar frequencies
($\sim\Omega_{ce}$), polarization (linearly polarized in $E_{y}$), growth
rates ($10^{-2}\leavevmode\nobreak\ \Omega_{ce}$), and energy conversion rates
(1%–3%). The properties can be well explained using the linear kinetic theory.

As the kinetic theory suggests, the X mode can be amplified only with
$n_{e}/n_{0}=100\%$. In other cases with $0<n_{e}/n_{0}\leq 90\%$, the R mode
is the most unstable and the X mode cannot grow. In our simulation with
$n_{e}/n_{0}=90\%$, we observed significant declines of wave energies after
the linear stage, with the cold electrons significantly accelerated at 50
$\omega_{pe}^{-1}$. We propose a multi-step process to explain the X-mode
excitation in this case (and its inconsistency with the linear theory): the
R-mode waves are first amplified via ECMI, then interact strongly with the
cool electrons and become damped soon via the gyro-resonant absorption. This
results in the bulk acceleration of the minor cool electrons, which then
become the “hot” component. This leads to an effective “$n_{e}/n_{0}=100\%$”
case, and thus the strong coupling of the R and X modes and the strong
excitation of the X mode. This means that the non-escaping R mode plays a
crucial role in the dynamic evolution of electron VDFs and the excitation of
the X mode.

However, with the actual measured parameters, we obtained significant
excitation of the trapped R mode rather than the escaping X-mode radiation.
This does not explain SKR observed from a distance. We propose two
possibilities to explain this problem:

1\. The measured electron VDFs are not the direct driver of ECMI due to the
poor resolution. The CAPS/ELS instrument measures the electrons of different
pitch angles using eight detectors, which rotate for approximately 3 min to
cover $2\pi$ sr of the sky. The temporal resolution is much larger than the
timescales of the wave growth ($\sim 100$$\omega_{pe}^{-1}$). It is hard to
measure electrons right before the excitation of ECMI since the VDFs diffuse
rapidly at a timescale of 50–100 $\omega_{pe}^{-1}$. Therefore, the ECMI of
SKRs may not be induced by the measured electron VDFs.

2\. The in situ observation cannot distinguish the R and X modes in the source
region. According to our simulations and the linear theory analysis, both the
R and X modes are identified as transverse electromagnetic waves with
$\omega\sim\Omega_{ce}$, linearly polarized in $E_{y}$. These properties are
consistent with the observed SKR emissions. With the Cassini in situ
measurements of the electric field only, it is difficult to distinguish the
two wave modes. In other words, it is possible that the enhancement of the
fluctuating electric fields observed by Cassini near or within the sources is
actually the trapped R mode, rather than the generally believed escaping X
mode.

Saturn’s kilometric radiation has long been considered to be the analogy of
the AKR from the geo-magnetosphere. Similar problems were present when
applying the ECMI mechanism to AKRs. It is believed that AKRs are generated
from low-density cavities with no cold electron population. According to the
PIC simulations of the finite source cavity (Pritchett, 1986; Pritchett &
Winglee, 1989; Pritchett et al., 2002), the generated X-mode radiations could
escape from the cavities via mode conversion at the boundaries. Although the
presence of such cavity structures has not been conclusively identified in SKR
sources, this could be a possible explanation for SKR.

In the present simulations, we employed periodical boundary conditions to
model the excitation of the wave modes, assuming the electrons are uniformly
distributed in a large spatial scale. No replenishment of the electron free
energy was included in this model. We will conduct simulations with various
setups of the source configurations to examine their effect on the generation
process of SKR .

The ECMI growth rates of the X and Z modes in SKR sources have been studied by
Mutel et al. (2010) and Menietti et al. (2011). As we pointed out in Section
1, previous studies employed oversimplified dispersion relations with
inappropriate assumptions and did not consider the R mode. In this paper, we
performed PIC simulations to study the excitation of the three extraordinary
modes (R, X, and Z), focusing mainly on the R mode and its effect on the
excitation of the X mode. We also performed linear analysis of the plasma
kinetic theory to obtain the dispersion relations and wave growths. The
results are in agreement with the PIC simulations.

The R mode is dominant in most cases with both the cold and hot electrons,
reaching energies of $\sim 0.03\leavevmode\nobreak\ E_{k0}$. This mode is
trapped between plasma layers delimited by the two electron gyro-frequencies
($\Omega_{ce}/\gamma,\leavevmode\nobreak\ \Omega_{ce}$), yet it may be
converted to escaping emissions when propagating in inhomogeneous plasmas. In
addition, the R mode can be absorbed by electrons, leading to the acceleration
and heating of plasmas, which, in turn, can modify the dispersion relations.
Further simulations should investigate the heating and acceleration process
due to the R mode and its coupling with the escaping X mode.

###### Acknowledgements.

This study is supported by NSFC grants (11790303 (11790300), 12203031, and
12103029), the China Postdoctoral Science Foundation (2022TQ0189), and the
Natural Science Foundation of Shandong Province (ZR2021QA033). We thank the
National Supercomputer Centers in Tianjin and the Beijing Super Cloud
Computing Center (BSCC, URL: http://www.blsc.cn/) for providing HPC resources,
and LANL for the open-source VPIC code.

## References

  * Bingham et al. (2013) Bingham, R., Speirs, D. C., Kellett, B. J., et al. 2013, Space Sci. Rev., 178, 695
  * Bowers et al. (2008) Bowers, K. J., Albright, B. J., Yin, L., Bergen, B., & Kwan, T. J. T. 2008, Physics of Plasmas, 15, 055703
  * Bowers et al. (2009) Bowers, K. J., Albright, B. J., Yin, L., et al. 2009, in Journal of Physics Conference Series, Vol. 180, Journal of Physics Conference Series, 012055
  * Dory et al. (1965) Dory, R. A., Guest, G. E., & Harris, E. G. 1965, Phys. Rev. Lett., 14, 131
  * Gurnett et al. (2004) Gurnett, D. A., Kurth, W. S., Kirchner, D. L., et al. 2004, Space Sci. Rev., 114, 395
  * Jackman (2011) Jackman, C. M. 2011, in Planetary, Solar and Heliospheric Radio Emissions (PRE VII), ed. H. O. Rucker, W. S. Kurth, P. Louarn, & G. Fischer, 1–12
  * Kaiser et al. (1981) Kaiser, M. L., Desch, M. D., & Lecacheux, A. 1981, Nature, 292, 731
  * Kaiser et al. (1980) Kaiser, M. L., Desch, M. D., Warwick, J. W., & Pearce, J. B. 1980, Science, 209, 1238
  * Kimura et al. (2013) Kimura, T., Lamy, L., Tao, C., et al. 2013, Journal of Geophysical Research (Space Physics), 118, 7019
  * Kurth et al. (2005) Kurth, W. S., Gurnett, D. A., Clarke, J. T., et al. 2005, Nature, 433, 722
  * Kurth et al. (2011) Kurth, W. S., Gurnett, D. A., Menietti, J. D., et al. 2011, in Planetary, Solar and Heliospheric Radio Emissions (PRE VII), ed. H. O. Rucker, W. S. Kurth, P. Louarn, & G. Fischer, 75–85
  * Lamy et al. (2009) Lamy, L., Cecconi, B., Prangé, R., et al. 2009, Journal of Geophysical Research (Space Physics), 114, A10212
  * Lamy et al. (2011) Lamy, L., Cecconi, B., Zarka, P., et al. 2011, Journal of Geophysical Research (Space Physics), 116, A04212
  * Lamy et al. (2010) Lamy, L., Schippers, P., Zarka, P., et al. 2010, Geochim. Res. Lett., 37, L12104
  * Lamy et al. (2008) Lamy, L., Zarka, P., Cecconi, B., et al. 2008, Journal of Geophysical Research (Space Physics), 113, A07201
  * Lamy et al. (2018) Lamy, L., Zarka, P., Cecconi, B., et al. 2018, Science, 362, aat2027
  * Linder et al. (1998) Linder, D. R., Coates, A. J., Woodliffe, R. D., et al. 1998, Geophysical Monograph Series, 102, 257
  * Menietti et al. (2011) Menietti, J. D., Mutel, R. L., Schippers, P., et al. 2011, Journal of Geophysical Research (Space Physics), 116, A12222
  * Mutel et al. (2010) Mutel, R. L., Menietti, J. D., Gurnett, D. A., et al. 2010, Geochim. Res. Lett., 37, L19105
  * Mutel et al. (2007) Mutel, R. L., Peterson, W. M., Jaeger, T. R., & Scudder, J. D. 2007, Journal of Geophysical Research (Space Physics), 112, A07211
  * Ning et al. (2021a) Ning, H., Chen, Y., Ni, S., et al. 2021a, A&A, 651, A118
  * Ning et al. (2021b) Ning, H., Chen, Y., Ni, S., et al. 2021b, ApJ, 920, L40
  * Pritchett (1984a) Pritchett, P. L. 1984a, Geochim. Res. Lett., 11, 143
  * Pritchett (1984b) Pritchett, P. L. 1984b, J. Geophys. Res., 89, 8957
  * Pritchett (1986) Pritchett, P. L. 1986, J. Geophys. Res., 91, 13569
  * Pritchett et al. (2002) Pritchett, P. L., Strangeway, R. J., Ergun, R. E., & Carlson, C. W. 2002, Journal of Geophysical Research (Space Physics), 107, 1437
  * Pritchett & Winglee (1989) Pritchett, P. L. & Winglee, R. M. 1989, J. Geophys. Res., 94, 129
  * Robinson (1986) Robinson, P. A. 1986, Journal of Plasma Physics, 35, 187
  * Schippers et al. (2011) Schippers, P., Arridge, C. S., Menietti, J. D., et al. 2011, Journal of Geophysical Research (Space Physics), 116, A05203
  * Strangeway (1985) Strangeway, R. J. 1985, J. Geophys. Res., 90, 9675
  * Strangeway (1986) Strangeway, R. J. 1986, J. Geophys. Res., 91, 3152
  * Treumann (2006) Treumann, R. A. 2006, A&A Rev., 13, 229
  * Winglee (1983) Winglee, R. M. 1983, Plasma Physics, 25, 217
  * Winglee (1985) Winglee, R. M. 1985, ApJ, 291, 160
  * Wu (1985) Wu, C. S. 1985, Space Sci. Rev., 41, 215
  * Wu & Lee (1979) Wu, C. S. & Lee, L. C. 1979, ApJ, 230, 621
  * Young et al. (2004) Young, D. T., Berthelier, J. J., Blanc, M., et al. 2004, Space Sci. Rev., 114, 1

##  Appendix A Linear analysis of plasma kinetic theory

###  A.1 Dispersion tensor

We followed Pritchett (1984b) to derive the dielectric tensors of plasmas with
the VDF given by the $\delta$ function (Eq. (3)). For perpendicular
propagating waves, the elements of plasma dispersion tensors were simplified
as

$\displaystyle\varsigma_{xx}=1-r^{\prime}\frac{\omega_{pe}^{2}}{\omega^{2}-1}-r\sum_{n=-\infty}^{\infty}\frac{n^{2}\omega_{pe}^{2}}{\xi^{2}\omega
b}\left(2\xi J_{n}J_{n}^{\prime}-G_{r}\frac{\omega J_{n}^{2}}{b}\right)$ (5)
$\displaystyle\varsigma_{xy}=ir^{\prime}\frac{\omega_{pe}^{2}}{(\omega^{2}-1)\omega}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
+ir\sum_{n=-\infty}^{\infty}\frac{n\omega_{pe}^{2}}{\omega\xi
b}\left[J_{n}J_{n}^{\prime}+\xi(J_{n}J_{n}^{\prime\prime}+J_{n}^{\prime
2})-G_{r}\frac{\omega J_{n}J_{n}^{\prime}}{b}\right]=-\varsigma_{yx}$
$\displaystyle\varsigma_{yy}=1-\frac{k^{2}c^{2}}{\omega^{2}}-r^{\prime}\frac{\omega_{pe}^{2}}{\omega^{2}-1}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
-r\sum_{n=-\infty}^{\infty}\frac{\omega_{pe}^{2}}{\omega
b}\left[2J_{n}^{\prime 2}+2\xi
J_{n}^{\prime}J_{n}^{\prime\prime}-G_{r}\frac{\omega J_{n}^{\prime
2}}{b}\right]$
$\displaystyle\varsigma_{xz}=\varsigma_{zx}=\varsigma_{yz}=\varsigma_{zy}=0$
$\displaystyle\varsigma_{zz}=1-\frac{k^{2}c^{2}}{\omega^{2}}-r^{\prime}\frac{\omega_{pe}^{2}}{\omega^{2}}-r\frac{\omega_{pe}^{2}}{\gamma_{r}\omega^{2}}$

,

where $r=n_{e}/n_{0}$, $\xi=ku_{r}/\Omega_{ce}$, $r^{\prime}=1-r$,
$G_{r}=\gamma_{r}-1/\gamma_{r}$, $\gamma_{r}=(1+u_{r}^{2}/c^{2})^{1/2}$,
$b=\gamma_{r}\omega-n\Omega_{ce}$, and $J_{n}(\xi)$ refers to the usual Bessel
function of order $n$.

The following equation was solved to obtain both the frequency and growth rate
of each extraordinary mode (X, Z, or R):

$\varsigma_{xx}\varsigma_{yy}-\varsigma_{xy}\varsigma_{yx}=0.$ (6)

###  A.2 ECMI of X and R with $n_{e}/n_{0}=24\%$

To understand the PIC simulations, we first conducted linear theory analysis
with the parameters of SKR sources reported by Menietti et al. (2011). We set
$u_{r}=0.18\leavevmode\nobreak\ c$, $\gamma_{r}=1.016$, $n_{e}/n_{0}=24\%$,
and $\omega_{pe}/\Omega_{ce}=0.046$.

With a given $k$, we obtained four sets of complex roots
$\omega=\omega_{r}+i\Gamma$ of Eq. (6), representing different branches of the
wave modes (see Fig. 7(a–b)). The dispersion curves of the X mode lie in the
range of $\omega_{r}>\Omega_{ce}$ with a zero growth rate. The slow branch of
the extraordinary mode (Z) is below $\Omega_{ce}/\gamma_{r}$. Within the range
of $\Omega_{ce}/\gamma_{r}<\omega_{r}<\Omega_{ce}$ and
$k\lesssim\Omega_{ce}/c$, there exist two conjugate roots of the R mode (R1
and R2). The maximum absolute value of their growth rates is 0.006
$\Omega_{ce}$. With increasing $k$ ($>\Omega_{ce}/c$), R1 and R2 are decoupled
and approach their respective resonant frequency, that is, $\sim\Omega_{ce}$
and $\Omega_{ce}/\gamma_{r}$.

Figure 7: Dispersion relation (a, d) and growth rates (b, e) of the
extraordinary wave modes according to the kinetic linear theory, with
$n_{e}/n_{0}=24\%$ (upper panels) and 100% (lower panels). The blue and red
lines refer to the X and Z modes; the purple and yellow lines refer to the two
branches of the relativistic modes (R1 and R2). Panels (c) and (f) present the
corresponding intensity ratio of $E_{x}/iE_{y}$ for each mode.

The perpendicular-propagating extraordinary waves are mainly carried by
$E_{x}$ and $E_{y}$. The polarization was determined by

$\frac{E_{x}}{E_{y}}=-\frac{\varsigma_{xy}}{\varsigma_{xx}}.$ (7)

Figure 7(c) displays the variations in $E_{x}/iE_{y}$ with $k$ for each
extraordinary mode. With increasing $k$, $E_{x}/iE_{y}$ rises from $-1$ to 0
for the X mode and from 0 to infinity for the Z mode. The $E_{x}/iE_{y}$ of R
coincides with that of the X mode, with $k<\Omega_{ce}/c$, and with that of
the Z mode, with $k>\Omega_{ce}/c$.

With $n_{e}/n_{0}=24\%$, R is the most unstable mode, with the highest growth
rate at $k\sim\Omega_{ce}/c$; where the three modes are basically linearly
polarized with $E_{x}/iE_{y}\sim 0$, the X mode cannot be amplified via the
ECMI process.

###  A.3 Parameter studies on $n_{e}/n_{0}$, $\omega_{pe}/\Omega_{ce}$, and
$u_{r}$

We varied $n_{e}/n_{0}$ from 0 to 100%, keeping the other parameters fixed
(Fig. 8). With $n_{e}/n_{0}=0$, we obtained the magnetoionic cold-plasma
dispersion relation with two branches (X and Z), with no instability.

With $n_{e}/n_{0}>0$, we obtained conjugate pairs of solutions for the R1, R2,
and Z modes, and only those with positive growth rates are presented. The
unstable R1 and R2 modes appear within $\Omega_{ce}/\gamma<\omega<\Omega_{ce}$
and their growth rates rise with $n_{e}/n_{0}$. With a higher abundance of
energetic electrons, the X mode has a lower cutoff ($\omega_{\text{X}}$) and
its dispersion curve becomes closer to that of the R mode. The X mode is
always stable with $\omega_{\text{X}}>\Omega_{ce}$ if $n_{e}/n_{0}<100\%$.

With $n_{e}/n_{0}=100\%,$ the dispersion curve of the R and X modes merge
together, with the cutoff $\omega_{\text{X}}<\Omega_{ce}$. In this case, the X
mode can be amplified via ECMI with the highest growth rate, $\sim
0.011\leavevmode\nobreak\ \Omega_{ce}$. As Figure 7(d–f) shows, the X mode is
linearly polarized in $E_{y}$ ($E_{x}/iE_{y}\sim 0$) with
$k\gtrsim\Omega_{ce}/c$, and has the highest growth rate at
$k\sim\Omega_{ce}/c$. In summary, the linear theory suggests that the R mode
is unstable if both cold and hot components are present, while the X mode can
only be amplified with $n_{e}/n_{0}=100\%$.

We further expanded the parameter study by varying $\omega_{pe}/\Omega_{ce}$
and $u_{r}$ while fixing $n_{e}/n_{0}$ (=24%). We find that the R mode is
always dominant. With increasing $\omega_{pe}/\Omega_{ce}$ (from 0.01 to 0.25)
or $u_{r}$ (from 0.05 to 0.3 $c$), the dispersion curves of the R mode appear
in a broader range of frequencies with increasing growth rates.

Figure 8: Dispersion relation (upper panels) and growth rates (lower panels)
of the extraordinary wave modes (X, R1, R2, and Z) with varying $n_{e}/n_{0}$.


$s_{1}$ from phase 1 and with the guess of $g$, the attacker obtains
$x_{i}+s_{2}\mod(N-1)=(\log_{2}g)^{-1}\log_{2}(F_{i}\oplus s_{1})\mod(N-1)$
for all $L$ values at the attacker’s disposal (the attacker does not need
pairs in this phase, just individual flow label values).

Now the attacker can use the fact that $0\leq x_{i}<M$, and that $M$ is
significantly lower than $N-1$. This means that for the correct $g$, the set
$\\{(\log_{2}g)^{-1}\log_{2}(F_{i}\oplus
s_{1})\mod(N-1)\\}=\\{x_{i}+s_{2}\mod(N-1)\\}$ will have a “hole” of size of
at least $N-1-M$ in its coverage of $[0,N-1)$. A-priori, the attacker does not
know where the hole is (because the attacker does not know $s_{2}$), but the
existence of a hole whose size is at least $E=N-1-M=244332$ can be detected
even without knowing $s_{2}$, by splitting the range of $N-1$ numbers into
segments whose size is no larger than $\frac{E}{2}$. It is then guaranteed
that one of the segments will be empty. In our case, the maximal segment size
is $122166$, and the attacker needs 5 segments. For simplicity, the attacker
can use equal-sized segments with $\lceil\frac{N-1}{5}\rceil=104854$ elements
each. For each guess of $g$, the attacker goes over the the list of
$(\log_{2}g)^{-1}\log_{2}(F_{i}\oplus s_{1})\mod(N-1)$ values, and each value
is divided by the segment size to find the corresponding segment. The attacker
stops if there are hits in all 5 segments. In such case, the attacker discards
the guess of $g$ and moves on to the next one. For the correct guess of $g$,
on the other hand, there will be at least one segment that is not a hit by all
$L$ values.

Interestingly, due to the algebraic properties of the series $(x_{i})$, there
may be incorrect $g^{\prime}$ values that cause clustering in the values

$(\log_{2}g^{\prime})^{-1}\log_{2}(F_{i}\oplus s_{1})\mod(N-1)=$
$(\log_{2}g^{\prime})^{-1}(\log_{2}g)(x_{i}+s_{2})\mod(N-1)$

This can result in “holes” large enough to become false positives. When the
$(x_{i})$ series is truly random, this is statistically impossible, given
large enough $L$ values (e.g. $L=100$), yet when $(x_{i})$ is generated via a
linear congruential generator, this does happen.

The following paragraph contains an in-depth explanation of this phenomenon,
and the reader can skip it if these details are not necessary.

To demonstrate such false positive, define:

$\gamma=(\log_{2}g^{\prime})^{-1}\log_{2}g\mod(N-1)$

Then the values that cluster are $(\gamma\cdot x_{i}\mod(N-1))$ (we can drop
$s_{2}$ as it is merely a shift, and so does not affect clustering). Now
consider a situation wherein $\gamma\cdot b$ and $\gamma\cdot M$ are “small”
modulo $N-1$ (by that we mean they are close to either 0 or to $N-1$).
Additionally, assume that

$\frac{M}{\gcd(a-1,M)}\big{(}\gamma\cdot\gcd(a-1,M)\mod(N-1)\big{)}$
$\ll(N-1)$

For simplicity, let us assume that $n_{i}=1$ and that the series is extracted
sequentially, i.e. $x_{i+1}=ax_{i}+b\mod M$. Now, $x_{i+1}-x_{i}\mod
M=(a-1)x_{i}+b\mod M$ and $x_{i}(a-1)\mod
M=x_{i}\frac{a-1}{\gcd(a-1,M)}\gcd(a-1,M)\mod M$. There are only
$\frac{M}{\gcd(a-1,M)}$ possible values for this expression:
$\\{0,\gcd(a-1,M),2\cdot\gcd(a-1,M),\cdots,(\frac{M}{\gcd(a-1,M)}-1)\gcd(a-1,M)\\}$.
Define $\Delta_{i}=x_{i+1}-x_{i}$ (over the integers). Then $\Delta_{i}$ is
either $t_{i}\cdot\gcd(a-1,M)+b$ or $t_{i}\cdot\gcd(a-1,M)+b-M$ or
$t_{i}\cdot\gcd(a-1,M)+b-2M$, for some $0\leq t_{i}<\frac{M}{\gcd(a-1,M)}$.
Note that these expressions are over the integers (not modulo $M$). When we
multiply this by $\gamma$ (modulo $N-1$), we get
$\gamma\cdot\Delta_{i}\mod(N-1)=t_{i}\cdot(\gamma\cdot\gcd(a-1,M)\mod(N-1))+(\gamma\cdot
b\mod(N-1))-u_{i}\cdot(\gamma\cdot M\mod(N-1))$ where $0\leq
t_{i}<\frac{M}{\gcd(a-1,M)}$ and $0\leq u_{i}\leq 2$. Each $\Delta_{i}$ is
therefore “small” (compared to $N-1$), so
$x_{i}=x_{0}+\sum_{j=0}^{i-1}\Delta_{j}\mod(N-1)$ does not get too far from
$x_{0}$. This explains why there may be (few) false positives.

At the end of phase 2, the attacker has the correct $g$, perhaps along with
some false positives. The attacker also knows which segment is inside the hole
(the empty segment) for each $g$ candidate. There may be more than one such
segment, and in such case the attacker can choose one segment arbitrarily as
the empty segment.

Phase 3 – Eliminating the $g$ False Positives, and Finding $s_{2}$
(Partially): For each $g$ candidate, the attacker finds the minimum $S_{min}$
and maximum $S_{max}$ of the set $\\{(\log_{2}g)^{-1}\log_{2}(F_{i}\oplus
s_{1})\mod(N-1)\\}=\\{x_{i}+s_{2}\mod(N-1)\\}$ such that the range
$[S_{min},S_{max}]$ form a continuous segment over the integers. If the set is
entirely below the empty segment, then the minimum and maximum are taken from
the set as usual. Likewise if the set is entirely above the empty segment. If
there are elements in the set both above and below the empty segment, then the
set starts above the empty segment, wraps around (courtesy of the modulo
$(N-1)$) and ends below the empty segment. Therefore, the minimum is taken
from the entries above the empty segment, and the maximum of the set is taken
from the entries below the empty segment, to which we add $N-1$ in order to
“undo” the wrap-around.

The attacker can easily eliminate almost all false positives for $g$ by
discarding ranges larger than $M$, i.e. discard $g$ if $S_{max}-S_{min}+1>M$.
Additionally, the attacker can apply a heuristic that eliminates additional,
very rare false positives wherein the cluster is too “tight”, i.e. form a
range which is “too small”. For example, the attacker can eliminate $g$ values
that suggest $S_{max}-S_{min}+1<\frac{M}{2}$. Empirically, this leaves the
attacker with a single (correct) $g$.

Now clearly, $s_{2}\leq S_{min}$ and $S_{max}\leq s_{2}+M-1$. Put differently,
$S_{max}-M+1\leq s_{2}\leq S_{min}$ (note that above, we discarded $g$
candidates for which $S_{max}-S_{min}+1>M$ i.e. $S_{max}-M+1>S_{min}$,
therefore we know that $S_{max}-M+1\leq S_{min}$; also, we allow negative
values for $s_{2}$ here). We see that there are $M+S_{min}-S_{max}$ possible
values for $s_{2}$. Since $S_{max}$ is distributed like a maximum of $L$
random variables uniformly distributed on $[0,M)$,
$E(S_{max})=\frac{L}{L+1}M$, and similarly $E(S_{min})=\frac{1}{L+1}M$.
Therefore, $E(M+S_{min}-S_{max})=\frac{2}{L+1}M$.

At the end of phase 3, the attacker has a correct $g$ and an average of
$\frac{2}{L+1}M$ candidates for $s_{2}$, which is the range $S_{max}-M+1\leq
s_{2}\leq S_{min}$.

Phase 4 – Finding $a$ and Possible $b$ and $s_{2}$ Values: First, we make
another observation on Algorithm 2. Concerning the intermediate values
generated in line 13, it is easy to see that the output sequence induced by
$(x,a,b)$ is at offset $+k$ (modulo $M$) from the output sequence induced by
$((x-k)\mod M,a,(b+(a-1)k)\mod M)$ for any $k$. Let us assume that $s_{2}$ is
correct, and so are $(x_{0},a,b)$. Let us now look at an incorrect guess
$s^{\prime}_{2}=s_{2}+k\mod(N-1)$. This guess will lead to an (incorrect)
internal state $x^{\prime}_{0}=(x_{0}-k)\mod(N-1)$. If $\forall_{0\leq
i<L}(0\leq x_{i}-k<M)$, then $x^{\prime}_{i}=(x_{i}-k)\mod M$, and we can pair
$x^{\prime}_{0}=x_{0}-k$ with $a^{\prime}=a,b^{\prime}=(b+(a-1)k)\mod M$ to
form a sequence such that $\forall_{0\leq
i<L}(x_{i}+s_{2}\mod(N-1)=x^{\prime}_{i}+s^{\prime}_{2}\mod(N-1))$. We see,
therefore, that for every $k$ for which $\forall_{0\leq i<L}(0\leq
x_{i}-k<M)$, it is impossible to distinguish the correct sequence (with
$s_{2}$) from a sequence generated with $s_{2}+k$. This gives rise, however,
to an attack optimization. Instead of going over all $S_{max}-M\leq s_{2}\leq
S_{min}$, we can choose one arbitrarily, because every value in this range
represents a possible $s^{\prime}_{2}$ which keeps all $x^{\prime}_{i}$ in
$[0,M)$.

There is an important caveat here. While an arbitrary $S_{max}-M+1\leq
s^{\prime}_{2}\leq S_{min}$ is indistinguishable from the correct $s_{2}$ for
all the flow label values already observed, it may happen that a future $x$
will fall outside the observed cyclic range $[S_{min}-s_{2},S_{max}-s_{2}]$.
In such case, an arbitrarily chosen $s^{\prime}_{2}$ may fail to predict the
correct PRNG output for $x$. One can think of this “out of bounds” $x$ value
as redefining – reducing – the range of possible $s^{\prime}_{2}$. At the
extreme, if the attacker has outputs for $x_{i}=0,x_{j}=M-1$, then (and only
then) the attacker will have a singleton “range” for $s^{\prime}_{2}$, and
thus $s^{\prime}_{2}=s_{2}$.

The larger $L$ is, the less likely a random $x$ is to fall outside the range.
In fact, the number of possible “out of range” values is identical to the size
of the $s_{2}$ range, i.e. in expectancy of $\frac{2}{L+1}M$, therefore the
probability for a random $x$ to fall “out of range” is $\approx\frac{2}{L+1}$.
In fact, the attacker can choose $s_{2}=\frac{S_{max}-M+1+S_{min}}{2}$ to
further reduce this probability. Note that since negative values of $s_{2}$
are possible if $S_{max}-M+1<0$, we may need to apply modulo $(N-1)$ to the
chosen $s_{2}$.

After choosing $s_{2}$, the attacker uses pairs again, to find $a$ (the
correct one) and $b$ (which matches the chosen $s_{2}$). Given $s_{2},g,s_{1}$
and pairs $(F_{i},F^{\prime}_{i})$ of consecutive flow label values, the
attacker can calculate $x_{i}=(\log_{2}g)^{-1}\log_{2}(F_{i}\oplus
s_{1})-s_{2}\mod(N-1)$ for every $F_{i}$, and likewise $x^{\prime}_{i}$ for
$F^{\prime}_{i}$. The attacker also knows $n_{i}$ for every such pair. The
attacker starts with $a\mod 12=1$ and $b\mod 12=(\log_{2}g)^{-1}I\mod 12$. At
each step, the attacker extends $a$ and $b$ by $\log_{2}3$ bits each. The
attacker does this by defining $\tilde{a}=12m_{a}+(a\mod
12),\tilde{b}=12m_{b}+(b\mod 12)$ where $m_{a},m_{b}\in\\{0,1,2\\}$ and going
over all 9 combinations of $m_{a},m_{b}$. For each combination of
$m_{a},m_{b}$, i.e. for each candidate of $a\mod 36=\tilde{a},b\mod
36=\tilde{b}$ the attacker goes over all $P$ pairs with $n_{i}\neq 3\land
n_{i}\neq 6$ (it is easy to see that for pairs with $n_{i}=3$ or $n_{i}=6$,
all such $\tilde{a}$ values yield correct results, hence they add no
information), and counts how many pairs yield a consistent result, i.e.

$\tilde{C}=$ $|\\{i|x^{\prime}_{i}\mod
36=\tilde{a}^{n_{i}}x_{i}+(\tilde{a}^{n_{i}-1}+\dots+1)\tilde{b}\mod 36\\}|$

Since the condition holds for all true pairs, if all pairs are true, we expect
$\tilde{C}=\tilde{P}$ where $\tilde{P}=|\\{i|n_{i}\neq 3\land n_{i}\neq
6\\}|$. If there are false pairs, we must use a statistical approach instead.
For a false pair (and also, when the guess of $m_{a},m_{b}$ is wrong), the
probability for the condition to hold, for a single pair, is $\frac{1}{3}$.
Hence, for a wrong guess,
$\tilde{C}\sim\mathit{N}(\frac{\tilde{P}}{3},\frac{\sqrt{2}}{3}\sqrt{\tilde{P}})$.
We can use a threshold of $5\sigma$ again, and determine that a guess of
$m_{a},m_{b}$ is correct if

$\tilde{C}\geq\frac{\tilde{P}}{3}+5\frac{\sqrt{2}}{3}\sqrt{\tilde{P}}$

The attacker applies this technique repeatedly. First, the attacker calculates
$a\mod 2^{2}3^{3}$ and $b\mod 2^{2}3^{3}$, then $a\mod 2^{2}3^{4}$ and $b\mod
2^{2}3^{4}$, etc., until the attacker has $a\mod 2^{2}3^{7},b\mod 2^{2}3^{7}$.
At this point, the attacker moves to applying the same technique with powers
of two (this time, going over all pairs with $n_{i}\neq 4\land n_{i}\neq 8$),
to get $a\mod 2^{3}3^{7},b\mod 2^{3}3^{7}$, etc., until the attacker gets
$a\mod 2^{7}3^{7}=a\mod M=a$, and $b\mod 2^{7}3^{7}=b\mod M=b$.

At the end of this phase, the attacker has $s_{1},g,s_{2},a,b$. With these
values, and a recent sample of the NetBSD flow label, $F$, the attacker can
predict the next value $F^{\prime}$ of the flow label – actually, the attacker
can provide 7 possible values, one of them is guaranteed to be correct. For
simplicity, we ignore the MSB prediction, which is trivial, and assume MSB=0.
Then the attacker can calculate $x$ which was used for generating $F$:
$x=(\log_{2}g)^{-1}\log_{2}(F\oplus s_{1})-s_{2}\mod(N-1)$. The attacker can
then enumerate over the 7 possible values of $n$, for each, the attacker
calculates $x^{\prime}=a^{n}x+(a^{n-1}+\dots+1)b\mod M$, and from that,
$F^{\prime}=s_{1}\oplus(g^{x^{\prime}+s_{2}}\mod N)$. In general, given
$g,s_{1},s_{2},a,b$ and a sample $F$, the attacker can roll forward the PRNG
indefinitely (up to reseeding, of course).

Covert Channel: This cryptanalytic result can be employed in order to
implement a covert channel between a sender and a receiver. The receiver first
establishes around 100 consecutive TCP connections with the target host,
extracts the IPv6 flow labels, and obtains the next PRNG internal outputs
using the cryptanalytic techniques. The sender then either establishes four
TCP connections with the target host, to signal the bit “1”, or does nothing,
to signal the bit “0”. Finally, the receiver establishes another TCP
connection with the target host, and extracts the final IPv6 flow label. The
receiver then matches this value to the series of eight possible future values
(corresponding to the PRNG internal steps). If a match is found in positions
2-8 of the list, then the receiver infers that the sender did not establish
four TCP connections with the target host (because each TCP connection incurs
at least two internal steps, and the receiver’s final connection also consumes
at least two internal steps – altogether 10 internal steps at the minimum).
This is interpreted by the receiver as bit “0”. Otherwise, the receiver
assumes the bit sent was “1”.

For host alias resolution use case, the attacker can improve the basic
technique a bit, since the attacker can observe the flow label received from
both endpoints, and match them to the same series (for endpoints that are
mapped to the same host) or not (for endpoints that are mapped to different
hosts).

####  A.2.3 A Covert Channel via NetBSD TCP/IPv6 Flow Label Most Significant
Bit

We now describe an attack against the NetBSD TCP/IPv6 flow label generation,
which does not make any use of the specific algorithm that generates the least
significant 19 bits of the flow label field. This is a more generic attack,
compared to Appendix A.2.2, but at the price of being much less efficient, in
terms of number of packets (TCP connections) required.

The attack is based on the observation that a “passive” side channel exists in
the flow label generator, in the form of the timing of the MSB flip. When two
TCP/IPv6 services are served from a single NetBSD host, the MSB in their flow
label fields will flip simultaneously. But when the services are served from
different NetBSD hosts, this is unlikely to happen.

We now calculate how many TCP connections are needed to force a reseed. In
general, let $S$ be the internal steps the flow label PRNG executed while
generating $n$ flow labels, then $S$ is a sum of $2n$ IIDs whose distribution
is uniform on $\\{1,2,3,4\\}$. Therefore,
$S\sim\mathit{N}(5n,\frac{\sqrt{10n}}{2})$, and $Prob(S\geq
200000)=\Phi(\frac{5n-200000}{\nicefrac{{\sqrt{10n}}}{{2}}})$. We can thus
force an immediate algorithm reseed and MSB flip with probability $\approx 1$
by establishing $n=41000$ TCP connections with the target host.

This can be turned into an “active” covert channel. The sender and receiver
agree on time $t_{0}$. Shortly before $t_{0}$, the sender establishes 41000
TCP connections with the target host and forces an algorithm reseed. Then
every $\Delta t$ seconds the sender either establishes 41000 TCP connections
with the target host to send the bit “1”, or does nothing to send the bit “0”.
Every $\Delta t$ seconds the receiver establishes a single TCP connection with
the target host, and observes whether the MSB bit flipped (signalling the bit
“1”) or not (“0”). $\Delta t$ should be chosen such that the sender can
establish 41000 TCP connections with the target host during this time period.
Special care must be taken if 180 seconds elapse since the last transmission
of the bit “1”, as this will cause the bit to flip due to a 180s timeout, and
may thus result in a false positive (reading bit “1” by the receiver). To
counteract this, the sender and receiver can decide to skip the time slice in
which 180s elapsed since the previous bit “1” was transmitted.

###  A.3 Linux NetFilter/nftables

NetFilter is a network filtering framework which is built into the Linux
kernel. NetFilter is controlled from user space by the nftables utility (the
nft command line). The NetFilter object hierarchy has tables that contain
chains of rules.

Assuming a table named override-ipid is created, e.g. by running:

nft add table override-ipid

And assuming an output filter chain is created at priority “raw” (i.e. at the
IPv4 level, but before fragmentation) is created, e.g. by running:

nft 'add chain override-ipid output {type filter hook output priority raw;}'

Then a safe configuration that only sets IPv4 ID=0 when DF=1 can be
implemented by running:

nft 'add rule override-ipid output ip frag-off & 0x4000 != 0 ip id set 0x0000'

A more aggressive configuration randomizes IP ID for ICMP packets, regardless
of DF. This can be implemented by running the following (in addition to the
above):

nft 'add rule override-ipid output ip protocol==1 ip id set numgen random mod
65536'

Finally, here is a very aggressive configuration, that unconditionally
randomizes IPv4 ID:

nft 'add rule override-ipid output ip id set numgen random mod 65536'

Note that until recently, numgen random was based on a weak Linux kernel
random number generator (this implementation was attacked in [30]). A more
secure version is available in kernel versions 5.19-rc4 and above, and in LTS
branches starting with versions 5.15.51 and 5.10.127. Using numgen random in
earlier versions is not recommended.


11institutetext: Federal University of Santa Catarina (UFSC), Joinville, Santa
Catarina, Brazil
11email<EMAIL_ADDRESS>
22institutetext: Institute for Computational Genomics, Joint Research Center
for Computational Biomedicine, RWTH Aachen University Medical School, 52074
Aachen, Germany
22email<EMAIL_ADDRESS>
# Clustering Validation with The Area Under Precision-Recall Curves

Pablo Andretta Jaskowiak 11 Ivan Gesteira Costa 22

###### Abstract

Confusion matrices and derived metrics provide a comprehensive framework for
the evaluation of model performance in machine learning. These are well-known
and extensively employed in the supervised learning domain, particularly
classification. Surprisingly, such a framework has not been fully explored in
the context of clustering validation. Indeed, just recently such a gap has
been bridged with the introduction of the Area Under the ROC Curve for
Clustering (AUCC), an internal/relative Clustering Validation Index (CVI) that
allows for clustering validation in real application scenarios. In this work
we explore the Area Under Precision-Recall Curve (and related metrics) in the
context of clustering validation. We show that these are not only appropriate
as CVIs, but should also be preferred in the presence of cluster imbalance. We
perform a comprehensive evaluation of proposed and state-of-art CVIs on real
and simulated data sets. Our observations corroborate towards an unified
validation framework for supervised and unsupervised learning, given that they
are consistent with existing guidelines established for the evaluation of
supervised learning models.

###### Keywords:

Clustering validation Clustering validation index CVI Area Under the
Precision-Recall Curve AUPR.

##  1 Introduction

Clustering is an unsupervised Machine Learning (ML) task, with a broad and
diverse spectrum of applications. Given a data set
$\mathbf{X}=\\{\mathbf{x}_{1},\dots,\mathbf{x}_{n}\\}$, with $n$ objects
embedded in a $m$-dimensional space, a typical clustering procedure consists
in partitioning $\mathbf{X}$ into a number of clusters $k$, usually unknown
beforehand. One expects that, at the end of the clustering procedure, objects
within each cluster are more similar to each other than objects that belong to
different clusters. Henceforth we consider hard clustering solutions, i.e.,
partitions, with $2\leq k\leq n-1$, expressed in the form
$\mathcal{C}=\\{C_{1},\dots,C_{k}\\}$, and subject to (i)
${C_{1}}\cup\dots\cup{C_{k}}=\mathbf{X}$; (ii) ${C_{i}}\neq\emptyset,\forall
i$; and (iii) ${C_{i}}\cap{C_{j}}=\emptyset,\forall i,j$ with $i\neq j$.

Given the unsupervised nature of clustering, the vast majority of its methods
will produce a result for any given set of input parameters. In such a
scenario, the validation of clustering results acquires critical importance.
By resorting to clustering validation procedures in the context of exploratory
data analysis, one is more likely to avoid meaningless solutions, selecting
only a narrow subset of relevant ones. Ideally, such solutions would be
further considered and inspected by a specialist of the subject under
investigation, which then might endorse and provide context to the findings.
Conversely, if clustering is employed as an automated pre-processing step
within a broader ML framework, e.g., feature selection, quantitative
validation procedures can provide the necessary degree of automation to select
the most significant clustering results that will be further considered
downstream in the pipeline, without the need of human intervention.

Clustering Validity Indices (CVIs) are the “tools” usually applied for the
quantitative validation of clustering results [24]. These mathematical indices
can be broadly divided into three categories: external, internal, and relative
[22]. External CVIs compare pairs of partitions, providing a score based on
their degree of correspondence [27]. Such indices are typically employed to
the development/evaluation of clustering algorithms, in which candidate
clustering solutions are compared against a reference partition (ground truth)
known beforehand. Internal validity criteria, take into account only the data
under analysis and the clustering outcome (partition). Such measures are
usually heavily influenced by the number of clusters from the partition under
evaluation and cannot be directly applied to distinguish among results with
different $k$. Finally, relative validity criteria comprehend a subset of
internal CVIs that can be employed to evaluate clustering results without
direct influence of their respective number of clusters (at least ideally).
Given that relative CVIs allow for the selection of “the best” clustering
outcome (and its corresponding number of clusters), they are ultimately
adopted in real clustering applications scenarios.

There is a significant number of relative CVIs available in the literature
[39], such as Silhouette, C-Index, Gamma, Dunn, and Davies-Boudin, just to
mention a few. This considerable number of indices suggests that no single
relative CVI can capture the notion of clustering, neither prevail in all
scenarios, as particular data characteristics might affect its performance in
unpredictable ways [12]. Indeed, such an universal view of clustering is, at
best, ill-advised, and has already been heavily criticized, with a call for
the identification of classes of clustering problems that might be treated
similarly [25]. Such an identification would not only avoid an universal and
simplified view of the clustering problem, but more importantly also prevent
the search for an (nonexistent) universal solution.

Just recently, it has been demonstrated that the Area Under the (ROC) Curve
(AUC), previously employed mostly to the validation of binary classification
models [14, 15, 16] also serve as an effective and computationally efficient
relative CVI [23]. In this work we further shorten the gap between the
evaluation of supervised and unsupervised learning models by introducing
relative CVIs based on Area Under the Precision-Recall Curve (AUPR). Our
motivation arises from the fact that AUPR has been shown to be more reliable
than AUC when dealing with heavily imbalanced classification problems [9, 36,
38]. We, therefore, propose the adoption of AUPR based relative CVIs, targeted
at the evaluation of clustering solutions with (heavy) imbalance. Such
problems are far from rare [7] and pose challenges in scenarios such as single
cell data analysis [41].

The remainder of the paper is organized as follows. In Section 2, we provide
the necessary background and revisit the Area Under the (ROC) Curve for
Clustering (AUUC), which will provide the basis for the introduction of AUPR
based CVIs, in Section 3. In Section 4 we detail the experimental design of
our empirical evaluation procedure, for which results are given in Section 5.
Final remarks and possible future work directions are discussed in Section 6.

##  2 Background

A relative CVI (Cluster Validity Index) allows for the quantitative evaluation
of a clustering solution. It provides a numerical score that relates to the
quality of the partition under evaluation, according, of course, to that
specific CVI. In order to select the most suited clustering solution among a
set of candidate ones (generated, for instance, with different clustering
algorithms and/or parameters), a user can apply a relative CVI and select the
partition that yields the best score.

Relative CVIs usually consider two basic characteristics of a partition to
quantify its overall quality, namely, intra-cluster compactness and inter
cluster separation [40]. A number of indices have been proposed based on
different definitions/formulations of these two concepts. A handful of reviews
and empirical evaluations of relative CVIs can be found in the literature [28,
29, 26, 39, 2]. Due to their practical appeal, implementations for a number of
such indices are readily available, e.g., in Python and R Packages [4, 6, 11,
32].

Just recently the Area Under the (ROC) Curve (AUC) has been borrowed from the
supervised learning domain and considered as a relative CVI for the validation
of clustering results [23]. The Area Under the (ROC) Curve of a clustering
solution which, in order to differentiate from AUC in the supervised setting,
is referred to as Area Under the (ROC) Curve for Clustering (AUCC), can be
computed based on: (i) a pairwise clustering solution and; (ii) the pairwise
similarities from the objects that compose such solution. Given a clustering
solution, i.e., a partition $\mathcal{C}$, a pairwise clustering solution
$\mathcal{C}^{p}$, with a total of $n(n-1)/2$ distinct pairs, can be readily
obtained with the following rule:

$\mathcal{C}^{p}(\mathbf{x}_{i},\mathbf{x}_{j})=\left\\{\begin{array}[]{l
l}1&\quad\mbox{if }\exists l:\mathbf{x}_{i},\mathbf{x}_{j}\in C_{l},\\\ \vskip
6.0pt plus 2.0pt minus 2.0pt\cr 0&\quad\mbox{otherwise}.\\\
\end{array}\right.$

Pairwise similarities can be readily obtained by the transformation of
distance metrics, i.e., Euclidean distance, into similarities.

Once the two pairwise relationships previously discussed are available, they
can be provided as input to a standard ROC Analysis procedure [14, 15, 16],
which, in the context of clustering validation, can be briefly summarized as
follows. Similarity values between all object pairs are sorted in a decreasing
fashion. Each numerically distinct similarity value is considered as a
decision threshold $\phi_{s}$, which establishes that: object pairs with a
similarity higher than $\phi_{s}$ are deemed “positive” (should be in a same
cluster), otherwise, they are deemed “negative” (should be in different
clusters). Binary predictions, as derived from a given similarity threshold
$\phi_{s}$, can be compared against the ones established by the pairwise
clustering solution $\mathcal{C}^{p}$. With the aid of a confusion matrix
(contingency table), as depicted by Figure 1, each pair of objects can, on the
basis of their pairwise similarity/threshold and clustering membership, be
regarded as: True Positive (TP), False Positive (FP), False Negative (FN), or
True Negative (TN). Finally, a ROC Plot can be obtained by relating all the
values of two particular metrics derived from the confusion matrices, namely
the True Positive rate ($TPR$) and false negatie rate ($FPR$). From such a ROC
Curve, the corresponding Area Under the Curve for Clustering (AUCC) can be
obtained, thus providing the quality score for the partition at hand.

Prediction ($\phi_{S}$) |  | Expected ($\mathcal{C}^{p}$) |
---|---|---|---
| Positive | Negative |
Positive |  $TP$ |  $FP$ |
Negative |  $FN$ |  $TN$ |

AUC/AUCC Metrics

$TPR=TP/(TP+FP)$

$FPR=FP/(FP+TN)$

AUPRC/AUPRCC/Related Metrics

$\text{Precision}=TP/(TP+FP)$

$\text{Recall}=TP/(TP+FN)$

$\text{Invese Precision}=TN/(TN+FN)$

$\text{Invese Recall}=TN/(TN+FP)$

Figure 1: A confusion matrix/contingency table (left) and some metrics (right)
that can be readily extracted from it for computing AUC (Area Under the ROC
Curve), AUPR (Area Under Precision-Recall curve) and related metrics. Note
that Inverse Precision and Inverse Recall are also known as Negative
Predictive Value and True Negative Rate, respectively. We opt for the
“Inverse” terminology from Powers [33], given that it emphasize the fact that
these are the values of Precision and Recall w.r.t. the negative class (Powers
refers to this as the “Inverse Problem”), in which one assumes that the
negative class is that of prevalent interest, contrary to the positive one.

The rationale behind this evaluation is that in a good clustering solution,
object pairs in a same cluster should have higher similarities than those in
different clusters. The relationship to ROC Analysis in classification is
straightforward: similarity values correspond to “classification thresholds”
whereas pairwise clustering memberships correspond to the “true classes”. It
has been shown that AUCC holds an expected value of 0.5 for random clustering
solutions and, moreover, that it is actually a linear transformation of the
Gamma relative validity index from Baker and Hubert [3], but with a much lower
computational cost [23].

##  3 The Area Under Precision-Recall Curve for Clustering

Even though results for AUCC are promising for relative clustering validation
in a general context, it is important to keep in mind that the ROC Curve and
its corresponding AUC are not the only quality measures that can be derived
from confusion matrices. In fact, in the context of supervised learning, AUC
is already known to perform poorly as a metric of quality for heavily
imbalanced problems [9, 36, 38]. These are our main motivations to explore
AUPR based metrics in the context of relative clustering validation. Bearing
that in mind, a counterpart of the Area Under Precision-Recall curve (AUPR)
for the relative validation of clustering results can be derived similarly as
shown before for AUCC. The major change arises from the fact that instead of
relating values of $TPR$ and $FPR$ for each one of the confusion matrices we
obtain during the process (one for each different similarity value), we now
relate values of Precision and Recall, which can easily computed as shown in
Figure 1.

It is worth noticing that in the context of binary classification problems
(supervised learning), the positive class usually has a somewhat greater
importance than the negative one. This is due to the fact that the positive
class is generally associated with the “case” of interest, which one aims to
correctly identify. Take, for instance, the problem of predicting whether a
patient is sick (positive) or healthy (negative). It can be argued that the
misclassification of a sick patient as healthy is more severe than the
misclassification of a healthy patient as sick. Therefore, in such problems,
most of the interest lies in: (i) the success rate of the classifier w.r.t.
its positive predictions, that is, from its positive predictions, how many are
actually positive (precision) and; (ii) the success rate of the classifier in
recovering/identifying actual positive cases (recall).

In some problems, however, such asymmetry might simply not exist. In other
problems, it may present itself reversed. In the latter case, the negative
class has more “importance” than the positive one, as discussed by Powers
[33], which refers to this as “Inverse Problem”. By taking these scenarios
into account, one may define counterpart metrics of Precision and Recall that
provide evaluation scores on the basis of the negative class, namely, Inverse
Precision (Negative Predicted Value) and Inverse Recall (True Negative Rate);
please, refer to Figure 1. Indeed, such an observation has already led to the
combination of measures/metrics and their inverse counterparts in the context
of external validation of clustering results, see for instance [1].

In the context of relative clustering validation, the problem we address in
the current paper, such an asymmetry between positive and negative classes is
neither observed nor desired. Notice, from the definition of AUCC, that one
deals with clustering memberships and similarities in terms of object pairs.
Therefore, all confusion matrices and statistics derived from them are also
w.r.t. object pairs. In this context, assigning a pair of objects from a same
cluster (positive pair) to different clusters (negative pair) is as bad as
assigning a pair of objects from different clusters (negative pair) to a same
cluster (positive pair).

Therefore, we consider for evaluation three relative CVIs based on the
concepts of the Area Under Precision-Recall curve (AUPR) and its Inverse
counterparts, namely: (i) the Area Under Precision-Recall curve for Clustering
(AUPRC), which relates the values of Precision and Recall; (ii) the Area Under
Inverse Precision-Recall curve for Clustering (AUIPRC), which relates the
values of Inverse Precision and Inverse Recall; and (iii) the average between
AUPRC and AUIPRC, which we refer to as Symmetric Area Under Precision-Recall
curves Clustering (SAUPRC). Note that the actual procedure adopted in order to
obtain such metrics is analogous to the one previously discussed in the
context of AUCC, except for the metrics which are derived from the confusion
matrices.

##  4 Experimental Design

###  4.1 Evaluation Procedure

In order to empirically assess the performance of the relative CVIs considered
in this study, we resort to a well established evaluation procedure from the
literature. This procedure was first discussed by [28] and more recently
formalized and employed by [39] to the empirical evaluation of a number of
relative CVIs. In brief, it is based on the reasoning that a good relative CVI
should provide evaluations that agree with that of an external CVI. Given a
data set with known cluster structure, i.e., a ground truth, the evaluation
procedure is as follows:

  1. 1.

Generate partitions of varied qualities for the data set. This may be achieved
by considering different clustering algorithms and configurations of their
parameters, ideally, leading to partitions with different numbers of clusters
($k$);

  2. 2.

For each one of the partitions from the previous step, compute its quality
according to: (i) the relative CVI under evaluation and; (ii) an external CVI
that quantifies its agreement w.r.t. the ground truth partition;

  3. 3.

Correlate the evaluation scores of each partition, as obtained from the
relative and external evaluations. This indicates the quality of the relative
CVI.

Four definitions are needed within this evaluation framework. The first one is
that of the clustering algorithms employed in order to generate a varied set
of partitions. In our study we considered k-means and agglomerative
hierarchical clustering algorithms with the following linkage functions:
Single, Average, Complete, and Ward’s [40]. The second one corresponds to the
range of the number of clusters ($k$) allowed for the partitions, for which we
consider $2\leq k\leq\lceil\sqrt{n}\rceil$, where $n$ is the number of objects
from the data set. The third and fourth definitions are related to the
external index and the correlation between external and relative evaluations.
For those we consider two well-known measures, i.e., the Adjusted Rand Index
(ARI) [20] and the Pearson correlation coefficient [17].

We evaluate within the same experimental setting several well established
relative CVIs from the literature [30, 42]. These are the same111In the
evaluation of AUCC [23], in addition to their original formulations, 3
variants of SWC and 17 variants of Dunn were considered. Given that such
variants provided similar results to their “original” counterparts, we believe
that their inclusion would bring more distractions than insights and,
therefore, do not consider them here.  measures employed in the previous study
that introduced AUCC [23]. In brief, besides the three relative CVIs we
introduce in this work, i.e., AUPRC, AUIPRC, and SAUPRC, we also take into
account AUCC, Silhouette Width Criterion (SWC) [35], Davies-Bouldin (DB) [8],
C Index [21], Dunn [13], PBM [31], Calinski-Harabasz (VRC) [5], Point Biserial
(PB) [28] and, Ratkowsky Lance (C/Sqrt(k) [34]. For a thorough review of such
measures we refer the reader to Vendramin et. al [39].

###  4.2 Evaluation Scenarios

Three data sources are considered in the evaluation, as we detail in the
following.

####  4.2.1 Julia Handl

This set comprises a benchmark of 80 data sets from [18], which account for
multivariate normally distributed clusters. The number of objects within each
data set ranges from around 800 to 4.000 (almost all datasets have different
numbers of objects). Given that this is a synthetic benchmark set, the
following parameters were employed by the aforementioned authors to generate
them. Regarding dimensions, data sets have 2 or 20 dimensions. The numbers of
clusters ($k$) considered are 4, 10, 20, and 40. Two different cluster sizes
are considered by authors, depending on the number of clusters present in the
data set. For data sets with 4 and 10 clusters, their sizes are uniformly
drawn within [50,500]; whereas for data sets with 20 and 40 clusters, their
sizes are uniformly drawn within [10,100]. Data sets were generated
iteratively, preventing overlap between clusters, in a trial and error
procedure. For more details regarding the data generation procedure, the
reader is referred to [18]. A detailed description of the generator and the
actual data sets are provided in one of the author’s webpage:
https://personalpages.manchester.ac.uk/staff/julia.handl/.

####  4.2.2 Tabula Muris

This set of data is comprised of 64 data sets from the Tabula Muris Compendium
[37], which is based on single cell experiments from different organs and
tissues of the _Mus musculus_ organism. The following 20 organs/tissues are
considered: aorta, bladder, brain (myeloid), brain (non-myeloid), diaphragm,
fat, heart, kidney, large intestine, limb muscle, liver, lung, mammary gland,
marrow, pancreas, skin, spleen, thymus, tongue and, trachea. Data from the two
available single cell isolation methods (FACS and droplet) were considered.
Previously normalized data were used [37]. In order to obtain the ground truth
for each data set, allowing thus for the use of an external index, we
considered the two available options: cluster IDs and cell ontology class.
Data sets were obtained with the Seurat Package [19].

####  4.2.3 Synthetic Data sets

This set of 2D Synthetic data was generated for this study, in an iterative
fashion. Cluster centers were constrained within the [0,500] interval in both
$x$ and $y$ axis. Cluster variances for both axis were identical, randomly
drawn from the the [2,50] interval. For a given data set, each cluster center
was randomly drawn in a sequential manner. If the newly generated cluster did
not overlap with any of the existing clusters (considering, their respective
variances), it was thus added to the data set. The process then continued
until all the desired clusters could be added to the data set. If a new
cluster could no be added to the data set after 10,000 attempts (center
draws), the maximum variance was decreased in 5% and the process would
continue, with an attempt to add a new cluster. Subsequent 5% decreases could
occur (each after another round of 10,000 failed attempts). It is wort
noticing, however, that this reduction was only necessary in a few cases
during the generation of data sets with large numbers of clusters. In total,
2,100 data sets were generated, from the following configurations: (i) number
of objects fixed in 500; (ii) 21 different numbers of clusters, in the
$[2,22]$ interval; (iii) 10 different cluster imbalances within the
$[0,10,20,\dots,80,90]$ interval, where each value corresponds to a percentage
of the number of objects assigned to one of the clusters, with the remaining
ones having the same size — 0(%) accounts for no imbalance. Combining number
of clusters and imbalances, we had a total of 210 configurations, for which we
generated 10 different data sets (replicates), thus totalling 2,100.

####  4.2.4 Data Characteristics

In order to provide some notion of the cluster imbalance present in these data
sets we depict, in Figure 2, boxplots regarding the cluster size imbalance
ratio (i.e., largest cluster size divided by smallest cluster size) for each
benchmark set (left), as well as a stratified analysis per imbalance for the
Synthetic data (right). Julia Handl data sets are, in general, well balanced.
As for Tabula Muris, a higher imbalance can be observed, with some extreme
cases, in which the largest cluster is approximately $150$ times larger than
the smallest one. For the Synthetic data, we have a wide range of cluster size
imbalance ratios, as defined by their design. These diverse data sets will
allow us to evaluate CVIs under distinct conditions, from completely balanced
to heavy imbalance.

Figure 2: Characterization of cluster size imbalance in the distinct
evaluation scenarios.

##  5 Results

###  5.1 Overview

Due to the distinct characteristics of each set of benchmark data, we
investigate them individually (Figure 3). We observe that Symmetric AUPR
(SAUPRC) obtains the highest correlation w.r.t. external evaluations from ARI
for both Julia Handl and Synthetic data. In the case of the Tabula Muris
benchmark set, VRC (Calinski Harabasz) provided the overall best results, in
terms of median correlation with ARI, followed by AUIPRC. Nonetheless, SAUPRC
appears as the fourth best CVI for this benchmark set. In general terms,
SAUPRC provides more consistent/robust evaluations than its base counterparts,
which are also evaluated here, namely, AUPRC and AUIPRC. Except for Julia
Handl data sets, for which AUPRC and AUIPRC appear next to each other, results
for these measure are quite contrasting for Tabula Muris and Synthetic data,
which can be observed by the change in their relative order w.r.t. the
boxplots from Figure 3. Taking such behaviours into account, by combining AUPR
and AUIPR, SAUPR provides robust evaluations in different settings.

Figure 3: Boxplots depict correlation between the evaluations of each relative
CVI and those of the Adjusted Rand Index (ARI), i.e., the external CVI.
Results are sorted according to the average correlation value, which is
indicated by a circle in each boxplot.

It is worth noticing that in all benchmark sets SAUPRC provided better results
than SWC (Silhouette), which is one of the most commonly employed CVIs from
the clustering literature. Moreover, results for SAUPRC are superior than
those of the recently introduced AUCC in all benchmark sets. Altogether, these
results indicate the power of CVIs based on metrics from confusion matrices,
given that AUCC also shows competitive results to SWC in all settings.

We applied statistical tests taking into account the results obtained within
each benchmark set. Given that we are dealing with the comparison of multiple
methods (CVIs) in multiple data sets, we resorted to the Friedman statistical
test followed by the Nemenyi posthoc test [10]. Results are shown in Figure 4.
Apart from the statistical differences themselves, these plots provide a nice
summary w.r.t. the mean ranks of each CVI. In terms of mean ranks, SAUPRC is
the best CVI for both Julia Handl and Synthetic data. As for Tabula Muris, it
ranks 4.47. In this latter set, VRC has a slightly higher mean rank than
AUIPRC (the lower the rank, the better the method). In terms of statistical
differences, the overall scenario considering each benchmark set is the
following: (i) for Julia Handl, SAUPRC is statistically superior than all the
remaining CVIs, except for C Index; (ii) for Tabula Muris, VRC and AUIPRC
provide the top results, but are not statistically superior than C/Sqrt(k),
SAUPRC, and C Index; and (iii) for the Syntetic set, SAUPRC is superiror to
all the remaining CVIs.

Figure 4: Results regarding Friedman-Nemenyi statistical tests for the three
benchmark sets: Julia Handl (left), Tabula Muris (center) and Synthetic data
(right). Mean rank values for each method are shown in the vertical axis
(lower mean ranks indicate better methods). The horizontal line for each
method accounts for its mean rank (central dot), as well as the critical
difference established by the test (highlighted in gray). A statistically
significant difference is observed between CVIs (in favor of the lower ranked
one) if the projections of their horizontal lines do not overlap in the $x$
axis.

###  5.2 Cluster Size Imbalance Effect

One advantage of simulated data is that we can check how different cluster
imbalances affect the performance of the CVIs. Here we present results for the
Synthetic data in a stratified manner, considering each imbalance separately.
As previously discussed, 10 different cluster imbalances within the
$[0,10,20,\dots,80,90]$ interval were considered. Each imbalance value
corresponds to the percentage of objects assigned to a single cluster, with
the remaining objects uniformly distributed among remaining clusters. The
value 0(%) accounts for no imbalance.

Results for such analysis are depicted by Fig. 5. Broadly speaking, three
distinct performance behaviours can be observed from the plot w.r.t.
imbalance.

The first one is that of CVIs with deteriorating performance as cluster
imbalance increases. This trend is noticeable for AUIPRC, C Index, PBM, PB,
and VRC. A quite extreme case is that of PB (Point Biserial), for which the
median correlation transitions from values close to $0.5$ in the case of
balanced data to $-1.0$ in the highest imbalance considered in our evaluation
design ($90\%$).

Figure 5: Correlation of relative CVIs with external index (ARI) considering
the different cluster imbalances from the Synthetic data sets. Balanced data
is indicated by “Bal.”.

The second distinctive behaviour is that of measures with a (mostly) stable
performance w.r.t. class imbalance. Within this category we can place AUCC,
DB, and SWC (arguably). In the case of SWC, performance is actually slightly
better with increasing imbalance. It is worth noticing, however, that this
particular CVI provides a high number of outlying evaluations for higher
imbalances.

Finally, the third observable behavior is that of CVIs that provide better
evaluations as cluster imbalance increases. In this category we can place
AUPRC, SAUPRC, Dunn, and C/Sqrt(k). The case of C/Sqrt(k) is quite peculiar,
as the evaluations of this particular CVI provide a negative agreement w.r.t.
those of the external index with balanced data. Its performance improves as
imbalance rises, but median correlation values stay below $0.5$, indicating
only a moderate agreement w.r.t. external evaluations. In the case of Dunn, a
similar analysis is possible. Its median correlations with ARI range, however,
from $0.0$ (no imbalance) to $0.75$ ($90\%$ imbalance). Here, SAUPRC and AUPRC
stand out as the superior CVIs within this category. Notice that even though
their performance is enhanced as imbalance increases, reaching median
correlation values close to $1.0$, for balanced data, their median correlation
values are already above $0.75$. Finally, it is interesting to observe that
SAUPRC combined the best behaviours of AUIPRC, which has its best performance
for balanced data, and AUPRC, for which best performance is observed with high
imbalance data.

##  6 Conclusions

The positive results obtained in the context of relative clustering validation
with the Area Under the ROC Curve for Clustering (AUCC) [23] opened new venues
of investigation. In the context of supervised learning, AUC is just one among
several metrics that can be derived from confusion matrices in order to
quantify the performance of classification models. In the realm of relative
clustering validation, however, most of this measures remain yet unexplored.
In this work, we proposed and investigated empirically three CVIs based on the
concept of the Area Under Precision-Recall (AUPR), taking a step further
towards a unified framework for the evaluation of supervised and unsupervised
models.

Our evaluation was conducted on synthetic data sets well established in the
clustering literature, i.e., the Julia Handl benchmark set, and data sets with
different cluster size imbalances, introduced by us, in order to isolate and
evaluate cluster size imbalance effect. We also considered data sets of Single
Cell experiments from the Tabula Muris project, which consists of a real
application scenario with different levels of imbalances in cluster sizes.

In general terms, the Symetric Area Under Precision-Recall for Clustering
(SAUPRC), which combines the evaluations obtained from AUPRC and AUIPRC
provided the best results for cases of high cluster size imbalance. This
particular measure also provided superior results for well balanced data, with
better overall results than established CVIs from the literature, such as
Silhouette (SWC) or AUCC. For such a reason, we believe that SAUPRC should be
favored, particularly in the presence of heavy cluster size imbalance.

As future work, we aim to investigate the effect of sub sampling in the
performance of metrics based on a given area under the curve (either ROC or
Precision-Recall Curves). Given that to compute such measures for a data set
with $n$ objects we need to derive a pairwise clustering solution of size
$n(n-1)$, it is possible that a sub sampling strategy may reduce the
computational time needed in order to compute such measure, without
deteriorating their performance. The actual degree of possible reduction,
however, needs to be investigated.

## Acknowledgements

PAJ thanks the German Academic Exchange Service (DAAD) for financial support
under the Research Stays for University Academics and Scientists 2022 Program,
allowing him to conduct this research at RWTH Aachen University.

## References

  * [1] Amigó, E., Gonzalo, J., Artiles, J., Verdejo, F.: A comparison of extrinsic clustering evaluation metrics based on formal constraints. Information Retrieval 12(5), 613 (2009). https://doi.org/10.1007/s10791-009-9106-z
  * [2] Arbelaitz, O., Gurrutxaga, I., Muguerza, J., Pérez, J.M., Perona, I.: An extensive comparative study of cluster validity indices. Pattern Recognition 46(1), 243 – 256 (2013). https://doi.org/https://doi.org/10.1016/j.patcog.2012.07.021, http://www.sciencedirect.com/science/article/pii/S003132031200338X
  * [3] Baker, F.B., Hubert, L.J.: Measuring the power of hierarchical cluster analysis. Journal of the American Statistical Association 70(349), 31–38 (1975)
  * [4] Brock, G., Pihur, V., Datta, S., Datta, S.: clValid: An R package for cluster validation. Journal of Statistical Software 25(4), 1–22 (2008), http://www.jstatsoft.org/v25/i04/
  * [5] Calinski, R., Harabasz, J.: A dentrite method for cluster analysis. Commun Stat 3, 1–27 (1974)
  * [6] Charrad, M., Ghazzali, N., Boiteau, V., Niknafs, A.: NbClust: An R package for determining the relevant number of clusters in a data set. J. of Statistical Software 61(6), 1–36 (2014), http://www.jstatsoft.org/v61/i06/
  * [7] Chawla, N.V., Japkowicz, N., Kotcz, A.: Editorial: Special issue on learning from imbalanced data sets. SIGKDD Explor. Newsl. 6(1), 1–6 (jun 2004). https://doi.org/10.1145/1007730.1007733
  * [8] Davies, D.L., Bouldin, D.W.: A cluster separation measure. IEEE Trans. Pattern Analysis and Machine Intelligence 1, 224–227 (1979)
  * [9] Davis, J., Goadrich, M.: The relationship between precision-recall and roc curves. In: Proceedings of the 23rd International Conference on Machine Learning. pp. 233–240. ICML ’06, ACM, New York, NY, USA (2006)
  * [10] Demšar, J.: Statistical Comparisons of Classifiers over Multiple Data Sets. Journal of Machine Learning Research 7, 1–30 (2006)
  * [11] Desgraupes, B.: clusterCrit: Clustering Indices (2016), https://CRAN.R-project.org/package=clusterCrit, r package version 1.2.7
  * [12] Dubes, R.C.: Cluster Analysis And Related Issues, pp. 3–32. Handbook of Pattern Recognition and Computer Vision, World Scientific (Aug 1993), 0
  * [13] Dunn, J.C.: Well separated clusters and optimal fuzzy-partitions. Journal of Cybernetics 4, 95–104 (1974)
  * [14] Fawcett, T.: An introduction to ROC analysis. Pattern Recognition Letters 27(8), 861–874 (2006)
  * [15] Flach, P., Hernández-Orallo, J., Ferri, C.: A coherent interpretation of AUC as a measure of aggregated classification performance. In: International Conference on Machine Learning (ICML) (2011)
  * [16] Flach, P.A.: Encyclopedia of Machine Learning, chap. ROC Analysis, pp. 869–875. Springer US, Boston, MA (2010). https://doi.org/10.1007/978-0-387-30164-8_733
  * [17] Freedman, D., Pisani, R., Purves, R.: Statistics. W.W. Norton (1998)
  * [18] Handl, J., Knowles, J.: An evolutionary approach to multiobjective clustering. IEEE Transactions on Evolutionary Computation 11(1), 56–76 (2007). https://doi.org/10.1109/TEVC.2006.877146
  * [19] Hao, Y., et al.: Integrated analysis of multimodal single-cell data. Cell (2021). https://doi.org/10.1016/j.cell.2021.04.048
  * [20] Hubert, L., Arabie, P.: Comparing Partitions. Journal of Classification 2, 193–218 (1985)
  * [21] Hubert, L.J., Levin, J.R., Hubert L, J.: A General Statistical Framework for Assessing Categorical Clustering in Free Recall. Psychological Bulletin 83(6), 1072–1080 (1976). https://doi.org/10.1037/0033-2909.83.6.1072
  * [22] Jain, A.K., Dubes, R.C.: Algorithms for Clustering Data. Prentice Hall (1988)
  * [23] Jaskowiak, P.A., Costa, I.G., Campello, R.J.G.B.: The area under the roc curve as a measure of clustering quality. Data Mining and Knowledge Discovery 36(3), 1219–1245 (May 2022). https://doi.org/10.1007/s10618-022-00829-0, https://doi.org/10.1007/s10618-022-00829-0
  * [24] Luna-Romera, J.M., Martínez-Ballesteros, M., García-Gutiérrez, J., Riquelme, J.C.: External clustering validity index based on chi-squared statistical test. Information Sciences 487, 1–17 (2019). https://doi.org/https://doi.org/10.1016/j.ins.2019.02.046
  * [25] von Luxburg, U., Williamson, R.C., Guyon, I.: Clustering: Science or art? In: Guyon, I., Dror, G., Lemaire, V., Taylor, G., Silver, D. (eds.) Proceedings of ICML Workshop on Unsupervised and Transfer Learning. Proceedings of Machine Learning Research, vol. 27, pp. 65–79. PMLR, Bellevue, Washington, USA (02 Jul 2012)
  * [26] Maulik, U., Bandyopadhyay, S.: Performance evaluation of some clustering algorithms and validity indices. IEEE Trans. Pattern Analysis and Machine Intelligence 24(12), 1650–1654 (2002). https://doi.org/10.1109/TPAMI.2002.1114856
  * [27] Meilǎ, M.: Comparing clusterings: an axiomatic view. In: Proceedings of the 22nd international conference on Machine learning - ICML ’05. pp. 577–584. ACM Press, New York, New York, USA (Aug 2005), http://dl.acm.org/citation.cfm?id=1102351.1102424
  * [28] Milligan, G.W.: A monte carlo study of thirty internal criterion measures for cluster analysis. Psychometrika 46(2), 187–199 (1981)
  * [29] Milligan, G.W., Cooper, M.C.: An examination of procedures for determining the number of clusters in a data set. Psychometrika 50(2), 159–179 (1985)
  * [30] Nguyen, T., Viehman, J., Yeboah, D., Olbricht, G.R., Obafemi-Ajayi, T.: Statistical comparative analysis and evaluation of validation indices for clustering optimization. In: 2020 IEEE Symposium Series on Computational Intelligence (SSCI). pp. 3081–3090 (2020). https://doi.org/10.1109/SSCI47803.2020.9308412
  * [31] Pakhira, M.K., Bandyopadhyay, S., Maulik, U.: Validity index for crisp and fuzzy clusters. Pattern Recognition 37, 487–501 (2004)
  * [32] Pedregosa, F., et al.: Scikit-learn: Machine learning in Python. Journal of Machine Learning Research 12, 2825–2830 (2011)
  * [33] Powers, D.M.W.: Evaluation: From precision, recall and f-measure to roc., informedness, markedness & correlation. Journal of Machine Learning Technologies 2(1), 37–63 (2011)
  * [34] Ratkowsky, D.A., Lance, G.N.: A criterion for determining the number of groups in a classification. Australian Computer Journal 10, 115–117 (1978)
  * [35] Rousseeuw, P.J.: Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. Journal of computational and applied mathematics 20, 53–65 (1987)
  * [36] Saito, T., Rehmsmeier, M.: The precision-recall plot is more informative than the roc plot when evaluating binary classifiers on imbalanced datasets. PloS one 10(3), e0118432–e0118432 (Mar 2015)
  * [37] Schaum, N., et al.: Single-cell transcriptomics of 20 mouse organs creates a tabula muris. Nature 562(7727), 367–372 (Oct 2018). https://doi.org/10.1038/s41586-018-0590-4, https://doi.org/10.1038/s41586-018-0590-4
  * [38] Sofaer, H.R., Hoeting, J.A., Jarnevich, C.S.: The area under the precision-recall curve as a performance metric for rare binary events. Methods in Ecology and Evolution 10(4), 565–577 (2019)
  * [39] Vendramin, L., Campello, R.J.G.B., Hruschka, E.R.: Relative clustering validity criteria: A comparative overview. Statistical Analysis and Data Mining 3(4), 209–235 (2010). https://doi.org/10.1002/sam.10080
  * [40] Xu, R., Wunsch II, D.: Clustering. IEEE Press (2009)
  * [41] Yu, L., Cao, Y., Yang, J.Y.H., Yang, P.: Benchmarking clustering algorithms on estimating the number of cell types from single-cell rna-sequencing data. Genome Biology 23(1), 49 (Feb 2022). https://doi.org/10.1186/s13059-022-02622-0, https://doi.org/10.1186/s13059-022-02622-0
  * [42] Zhou, S., Liu, F., Song, W.: Estimating the optimal number of clusters via internal validity index. Neural Processing Letters 53(2), 1013–1034 (Apr 2021). https://doi.org/10.1007/s11063-021-10427-8