# Particle-In-Cell Simulations of the Cassini Spacecraft’s Interaction with
Saturn’s Ionosphere during the Grand Finale

Zeqi Zhang1, Ravindra T. Desai1, Yohei Miyake2, Hideyuki Usui2, Oleg
Shebanits1,3

###### Abstract

A surprising and unexpected phenomenon observed during Cassini’s Grand Finale
was the spacecraft charging to positive potentials in Saturn’s ionosphere.
Here, the ionospheric plasma was depleted of free electrons with negatively
charged ions and dust accumulating up to over 95 % of the negative charge
density. To further understand the spacecraft-plasma interaction, we perform a
three dimensional Particle-In-Cell study of a model Cassini spacecraft
immersed in plasma representative of Saturn’s ionosphere. The simulations
reveal complex interaction features such as electron wings and a highly
structured wake containing spacecraft-scale vortices. The results show how a
large negative ion concentration combined with a large negative to positive
ion mass ratio is able to drive the spacecraft to the observed positive
potentials. Despite the high electron depletions, the electron properties are
found as a significant controlling factor for the spacecraft potential
together with the magnetic field orientation which induces a potential
gradient directed across Cassini’s asymmetric body. This study reveals the
global spacecraft interaction experienced by Cassini during the Grand Finale
and how this is influenced by the unexpected negative ion and dust
populations.

1Blackett Laboratory, Imperial College London, London, UK

2Education Center on Computational Science and Engineering, Kobe University,
Kobe, Japan

3Swedish Institute of Space Physics, Uppsala, Sweden

Key words: space vehicles – plasmas – planets and satellites: individual:
Saturn – planets and satellites: atmospheres – planets and satellites:
composition

Accepted 2021 March 8. Received 2021 March 5; in original form 2021 January 22

##  1 Introduction

Spacecraft immersed within a plasma will become electrically charged due to
currents of incident electron and ion species [Whipple, 1981, Garrett, 1981].
When the net current induced by these species are not zero, net charge is
accumulated and there is a potential difference between the spacecraft and the
surrounding plasma. Since the plasma species have different charges, each
current is either decreased or increased by the change of the spacecraft
potential, which continues until an equilibrium potential is reached, where
the net current of all plasma species sum up to zero. Understanding how the
spacecraft potential is affected by its environment is important for
interpretting the surrounding plasma conditions and on-board plasma
measurements which can be significantly affected by the potential difference.

Near a given spacecraft surface a sheath boundary layer screens the potential
of the surface over distances of the order of the Debye length [Robertson,
2013] and, for sufficiently high relative velocities between the spacecraft
and surrounding plasma, a wake of depleted plasma is produced behind the
obstacle [Al’pert et al., 1966, Ludwig et al., 2012, Miloch, 2014]. As the
spacecraft absorbs incident plasma, a shock is unable to form upstream as
would for an obstacle able to withstand the incoming flow. If the relative
velocity exceeds the Alfvénic or acoustic velocities, the associated
discontinuities will instead trail downstream and form an approximately
conical wake region of disturbed flow commonly referred to as a Mach cone
[Willis et al., 2011]. Significant asymmetries are introduced into the
spacecraft-plasma interaction by the presence of a strong ambient magnetic
field which generates a convective electric field in the spacecraft frame
[Marklund et al., 1994, Pecseli, 2012] and modifies the structure of the wake
[Darian et al., 2017, Usui et al., 2019]. The high velocities of the
significantly lighter electrons compared with ions mean that a spacecraft
immersed within a typical space plasma will charge to negative potentials
[Spitzer, 1941] although processes such as photoelectron and secondary
electron emission, can shift the potential to positive values [Roussel and
Berthelier, 2004, Engwall et al., 2006, Miloch and Vladimirov, 2009,
Yaroshenko et al., 2011].

When passing through Saturn’s ionosphere, the Cassini spacecraft’s floating
potential was, surprisingly, observed as positive on all encounters below 3000
km altitude [Morooka et al., 2019]. The presence of negative ions and dust
grains in Saturn’s ionosphere is evident from the LP measurements [Morooka et
al., 2019], which showed significant concentrations up to over 95% of the
negative charge density, at altitudes of 3200 km down to the closest
measurement at $\approx$1600km. These appear an intrinsic part of the giant
planet’s ionosphere and distinct from electrons depletions associated with
Saturn’s main rings [e.g. Farrell et al., 2018] and transient negative ion
populations observed near Saturn’s icy satellites [Coates et al., 2010, Desai
et al., 2018, Nordheim et al., 2020]. Unfortunately, it was not possible to
obtain a mass distribution of these negative ions or dust grains due to
Cassini’s plasma spectrometers being offline and the Grand Finale plasma
datasets therefore lack crucial pieces of information.

A body immersed in a plasma with large quantities of negative ions can attain
a positive potential due to the reduced electron currents, as was shown for a
dust grain by Kim and Merlino [2006] using orbital motion limit theory, and
provides some indication for what might be occurring with Cassini. Similar
electron depletions of up to $\approx$96 % were, however, observed by Cassini
within Titan’s ionosphere where the spacecraft potential was consistently
negative [Wahlund et al., 2005, Crary et al., 2009, Shebanits et al., 2016,
Desai et al., 2017a]. Significant unknowns therefore remain regarding how
spacecraft interact with these outer solar system plasmas.

Plasmas with a significant negative ion content (ion-ion or dusty plasmas)
possess very different characteristics to typical electron-ion plasmas. The
reduced mobility of the heavier negative charge carriers alters the electric
field screening phenomenon and increases Debye length scales. Plasma
conductivities can reverse [Muralikrishna and Kulkarni, 2006, Shebanits et
al., 2020] and the plasma can host a variety of altered instabilities and wave
phenomenon [Shukla and Mamun, 2002, Desai et al., 2017b]. The interaction of
Cassini with Saturn’s ionosphere therefore represents a class of physics quite
different to the classic view of spacecraft charging within space plasmas.

In this article, we describe a self-consistent three dimensional Particle-In-
Cell (PIC) study where a model Cassini spacecraft is simulated immersed within
plasmas representative of Saturn’s ionosphere as observed during the Grand
Finale. Section 2 describes the simulation approach, how Cassini is modelled,
and the available measurements of Saturn’s ionospheric plasma. Section 3 then
applies these simulations to Saturn’s ionosphere with regards to specific and
general parameterisations and conducts a parametric survey to assess the
sensitivity of the results to the measured and inferred plasma properties.
Section 4 concludes with a summary of the key findings.

##  2 Simulation Development

###  2.1 EMSES

The three dimensional Particle-In-Cell simulation code Electro-Magnetic
Spacecraft Environment Simulation (EMSES) has been developed for the self-
consistent analysis of spacecraft-plasma interactions on either an
electromagnetic or electrostatic basis [Miyake and Usui, 2009]. The
electrostatic version of EMSES is utilised as the typical Alfvén velocities in
Saturn’s ionosphere are significantly greater than the spacecraft velocity and
in this regime Alfvénic perturbations are assumed to form only a small
contribution to the plasma currents [Rehman and Marchand, 2014].

Within the PIC approximation [Hockney and Eastwood, 1981, Birdsall and
Langdon, 1991], individual particles are represented by large quantities of
“super-particles” which are integrated through the same equations of motion as
a real particle [Boris, 1970]. The plasma is defined to consist of an
arbitrary number of plasma species in a drifting Maxwellian velocity
distribution to represent in-flowing plasma. Each species has mass and charge
normalized to the proton scale with a real ion-to-electron mass ratio. In this
study a negatively charged ion component is included, which is required to
study the target plasma environments. The positions of the super-particles are
tracked continuously in the simulation box, but the electric and the magnetic
fields are assigned and updated only on the grid points based on the Yee
[1966] algorithm, and interpolated onto the super-particles’ positions during
calculation. The charge density profile at each timestep is used to solve
Poisson’s equation for the electrostatic potential with Dirichlet boundary
conditions.

The simulations are run in the spacecraft frame and consists of a three
dimensional box with inflow and outflow boundary conditions along a specified
direction of plasma flow, and periodic boundary conditions orthogonal to
these. Within the domain the spacecraft is considered as a perfect conductor
with separate boundary treatments for both longitudinal and transverse
electric fields. The spacecraft surface can accumulate charge caused by
impinging super-particles, and redistributes this to ensure the spacecraft is
equipotential using the capacity matrix method [Hockney and Eastwood, 1981,
Miyake and Usui, 2009].

Figure 1: Simulation configuration for Cassini during Grand Finale Rev 292
ingress at 2500 km Saturn altitude. The X, Y and Z axes correspond to the
Zs/c, Xs/c, and Ys/c axes of the Cassini spacecraft attitude coordinate
system, respectively, and the precise simulation and plasma parameters are
provided in Table 1.

###  2.2 Application to Cassini

####  2.2.1 Main Body

The Cassini spacecraft is a three-axis stabilised spacecraft of approximately
six metres in length. Previous spacecraft charging simulations have considered
Cassini in two dimensions [Olson et al., 2010] and as a cylinder in three
dimensions during the 2004 Saturn orbit insertion [Yaroshenko et al., 2011].
In this study we model Cassini’s non-uniform shape encompassed within the
simulation domain of 12.8$\times$12.8$\times$12.8 m3, across
128$\times$128$\times$128 grid cells. Cassini is approximated using three
structures: a large thin cylinder representing the antenna dish with an
approximated diameter of 4 metres and width of 0.55 metres, and a longer
cylinder representing the main body with an approximated diameter of 2.2
metres and length of 3.8 metres. These are separated by a 0.55 metres gap
between them, as per the real Cassini spacecraft. Figure 1 shows a schematic
of the simulation geometry where the X, Y and Z axes correspond to the Zs/c,
Xs/c, and Ys/c axes of the Cassini spacecraft attitude coordinate system,
respectively: where X is the direction of main thrusters, Y is opposite to the
Langmuir probe direction and Z completes the right-hand set.

The dish and the main body of the spacecraft are considered as a single
perfect conductor, as Cassini was designed. Although this approximation does
not account for the curvature of the antenna dish, or other instruments
attached onto the body, the most important factors for spacecraft charging are
the surface area of structures larger than the Debye length scales as well as
the ram profile for a fast moving object. As far as the net surface current of
the spacecraft and its potential evolution are concerned, this approximation
is judged to reasonably approximate dynamics associated with Cassini’s
asymmetric shape interacting with Saturn’s ionosphere.

####  2.2.2 Langmuir Probe

The Langmuir Probe (LP) is represented by a small sphere at the side of
Cassini, the precise extent of which is defined subgrid. The LP is held at a
bias to the main spacecraft, the precise voltage difference specified
according to which point in the LP voltage sweep between –4 V and 4 V that is
simulated. The time required for the currents to equilibriate is small
compared to the sweep time scale of 0.5 s [Gurnett et al., 2004], and a fixed
bias within the sweep is chosen. The whole Cassini model is floating within
the plasma environment and can therefore become charged relative to it but the
voltage between Cassini and the probe remains fixed.

For a Langmuir Probe at a net negative potential within an electron-ion
plasma, the total external plasma currents can be considered to derive purely
from the positive ion current. This is because the lighter electrons are
repelled and therefore their current contribution can be neglected. Similarly,
at a net positive potential the external plasma currents is saturated by the
faster moving electrons and can be considered to purely consist of electrons.
In the presence of large negative ions, however, the situation is different.
When the LP is negatively biased, the large negative ions and dust grains are
still collected due to their high inertia (mass). At positive biases the
current will be a combination of electrons and negative ions but the
significantly smaller electrons still dominate the total current. The
simulations are therefore run with the Langmuir Probe biased at a negative
potential of –3 V to ensure no electrons are collected as the bias potential
is higher than the electron thermal energy. The self-consistent addition of
the negative ion component therefore allows us to constrain their contribution
to the Langmuir Probe currents and to the overall potential accumulated by the
spacecraft.

###  2.3 Saturn’s Ionospheric Plasma

Within Saturn’s ionosphere the total current onto the spacecraft is a function
of the positive ions, electrons and negative ions, charged dust, secondary
electron emission and photoelectron emission. The photoelectron current at
Saturn’s orbit are several orders of magnitude lower than the ion and electron
currents [Holmberg et al., 2017, Shebanits et al., 2017], and the secondary
electron emission currents are assumed to be of a similar magnitude to this
[Morooka et al., 2019]. They are therefore not included in this analysis. The
simulated currents can therefore be expressed as

$I_{total}=I_{electron}+I_{ion^{+}}+I_{ion^{-}},$ (1)

where I${}_{ion^{+}}$ represents the positive ions and dust and
I${}_{ion^{-}}$, represents the negative ions and dust. While multiple
positive and negative plasma components can be included the lack of
information on the mass distribution function leads us to use a mean mass to
represent these respective components.

####  2.3.1 Cassini Observations

The measured composition of Saturn’s ionosphere provides the inputs to the
simulations which produce the currents in Equation 1. Post–2012 CAPS was,
unfortunately, turned off and the mass distribution on the ambient plasma in
Saturn’s ionosphere was therefore unknown. Cassini’s Ion and Neutral Mass
Spectrometer was able to provide some information on positive ion composition
but this was limited up to low masses of several amu due to the high
spacecraft velocity [Waite et al., 2018]. The Langmuir Probe was, however,
able to differentiate between the bulk ion, electron and negative ion currents
[Wahlund et al., 2018, Hadid et al., 2019, Shebanits et al., 2020] and
therefore provide estimates of their densities [Morooka et al., 2019]. The
electron density was determined from the positive bias side of the Langmuir
Probe sweep [Morooka et al., 2019]. The positive ion density could, however,
only be constrained as a lower limit due to the unknown mass distribution
[Shebanits et al., 2013, Morooka et al., 2019]. A lower limit on the negative
ion density was therefore also available by assuming quasi-neutrality.

Table 1: Environmental and System Simulations Parameters Environmental
Parameters
---
Plasma ion density, $n_{0}$ | 854 cm-3
Negative ion concentration, $n_{ni}$ | 40.2 %
Ion mass, $m_{i}$ | 1.4 amu
Negative ion mass, $m^{-}_{i}$ | 3.0 amu
Electron temperature, $T_{e}$ | 0.093 eV
Ion temperature, $T_{i}$ | 0.093 eV
Negative ion temperature, $T^{-}_{i}$ | 0.093 eV
Magnetic field, $\vec{B}$ |  $\mu$T
Flow velocity, $\vec{v}_{flow}$ | km s-1
Alfvén speed, $v_{A}$ | 11,149 km s-1
Ion acoustic speed, $v_{S}$ | 2.56 km s-1
Debye length, $\lambda_{D}$ | 10.05 cm
Electron gyroperiod, $\tau_{ge}$ | 2.38 $\mu$s
Electron plasma period, $\tau_{pe}$ | 3.81 $\mu$s
Ion gyroperiod, $\tau_{gi}$ | 6.12 ms
Ion plasma period, $\tau_{pi}$ | 0.193 ms
Negative ion gyroperiod $\tau^{-}_{gi}$ | 13.1 ms
Negative ion plasma period, $\tau^{-}_{pi}$ | 0.283 ms
System Parameters
Grid width, $\Delta$r | 10 cm
Time step, $\Delta$t | 0.033 $\mu$s
Simulation time, $t$ | 0.67 ms
Particles per cell | 20
Domain size | 12.83 m3
Probe Bias, $\phi_{LP}$ | –3 V

The electron density was observed to increase with decreasing altitude as
Cassini sampled denser regions of Saturns ionosphere [Persoon et al., 2019,
Morooka et al., 2019]. The estimated negative ion density was also observed to
increase and at a greater rate than the electrons. As electrons are lost to
the negative ions/dust grains, the electron depletion also increased with
decreasing altitude, and reached an estimated lower bound of $94$ % [Morooka
et al., 2019]. The electron depletion is traditionally quoted as a ratio of
negative to positive ion densities such that a 94 % depletion corresponds to 6
% of the total negatively charged species being electrons, with the rest being
negative ions. The Langmuir Probe determined a mean electron temperature of
$\approx$0.1 eV and, based upon current balance, an estimate of the minimum
mean positive ion mass of $\approx$5 amu at the deepest point sampled [Morooka
et al., 2019]. The electron temperature was used as a conservative upper limit
on the ion temperatures and reducing the ion temperatures below this was later
found not to have a significant influence on the simulation results. Within
Saturn’s ionosphere the spacecraft potential varied between –1.25 and +0.75 V
with a positive floating potential on every encounter below 3000 km.

####  2.3.2 Simulation Inputs

Rev 292 is selected as a representative flyby with sufficient data from which
input parameters are derived. In the first instance, Cassini is simulated at
an altitude of 2500 km during ingress where the LP observations revealed an
estimated plasma density of 854 cm-3, temperature of 0.093 $eV$, and the lower
bound on the estimated electron depletion was 40 %. A higher electron density
decreases the Debye length which is the primary constraint on the simulation
grid size. This altitude therefore provided an optimal balance between
computational load and studying regions of interest with large negative ion
concentrations. The magnetic field was measured to be 14,925 nT which results
in the electrons being highly magnetised with gyroradii of the order of
centimetres and therefore significantly less than the spacecraft size. The
positive and negative ions are only weakly magnetised, with gyroradii of tens
to hundreds of metres, and therefore much larger than the spacecraft and the
simulation domain. Collision frequencies in Saturn’s ionosphere are much less
than a herz [Shebanits et al., 2020] and the associated period is
significantly greater than the total simulation time of less than a
millisecond. During Rev 292 the spacecraft potential was –0.12 V at this
altitude, and varied between –0.75 and +0.65 V. The precise simulations input
parameters are provided in Table 1.

The relative orientation of the spacecraft to the plasma and magnetic field
varies very slowly throughout a flyby and is different for different flybys.
The simulations described herein do not attempt a produce a precise
reproduction of the Cassini measurements, instead attempting to understand the
key features of the plasma interaction and associated charge accumulated.

##  3 Results

###  3.1 Global Interaction

The simulated interaction between Cassini and Saturn’s ionosphere is shown at
Figure 2. The left-hand (a–b), centre (c–d) and right-hand (e–f) panels show
electron, ion and negative ion densities respectively, with the upper (a, c,
e) and lower panels (b, d, f) respectively showing x–y and y–z slices through
the simulation. The results are displayed at the end of the simulation run.
The two cylinders representing Cassini’s antenna dish and main body are
apparent as distinct regions absent of plasma. The plasma velocity is arriving
predominantly along the y-direction at an angle of $\approx 168^{\circ}$
degrees to the near-oppositely directed magnetic field. The Langmuir Probe is
pointed in the spacecraft ram direction and therefore directly samples the
incoming plasma flow. The probe is most visible in Figure 2(a) at x=4 and y=4
as a local depletion in the electron density due to the negative bias of –3 V.

Figure 2: Two dimensional slices of the electron (a–b), ion (c–d) and negative
ion density (e–f) for simulation parameters outlined in Table 1 and described
in Section 3.1. The magnetic field is oriented along 1.48$\hat{x}$,
–14.8$\hat{y}$, 1.24$\hat{z}$ and a schematic of the simulation geometry is
displayed in Figure 1.

.

The high relative velocity of the plasma interaction produces an extended wake
of depleted plasma. The length of the wake is greater for positive and
negative ions due to their larger inertia resulting in slower refilling. When
looking at a slice through Cassini’s main body, panels (b, d, f), the wake
exhibits the characteristic structure associated with supersonic flow around a
cylinder. The slice in the x-y plane, however, in panels (a, c, e) reveals the
wake as non-uniform and highly structured behind the antenna dish and the gap
between the antenna dish and the main body

The electron density reveals wing-like structures of enhanced and then
depleted density attached to the spacecraft in the y-z plane in panel (b).
Similar wing-like structures have been identified at moving bodies produced by
Alfvén and whistler waves [Drell et al., 1965, Neuebauer, 1980, Stenzel and
Urrutia, 1989], which propagate at characteristic velocities associated with
their relative speed to the spacecraft such that they advect downstream.
Electron-wing structures, similar to these reported herein, have been
identified as consisting of propagating Langmuir waves [Miyake et al., 2020]
produced by electrons reflected from a negatively charged spacecraft which are
then guided by the magnetic field lines. Calculation of the Langmuir Probe
group velocity, cL=$\sqrt{3k_{B}T_{e}m_{e}}\cong$ 1,200 km s-1, and the angle
to the magnetic field, $\theta=\arctan({v_{flow}/c_{L}})\approx 2^{\circ}$,
confirms this mode propagates at small angles to the magnetic field in the
spacecraft frame, as can be seen in Figure 2. Figure 2a shows these wings
structures striking the inflow boundary condition and modifying the electron
density. Figure 2b however reveals that the effects of this are shifted
towards positive z-regions and mostly miss the spacecraft interaction and
therefore have a negligible numerical impact. The electron-wings propagate
upstream of Cassini and, for flybys where the magnetic field is more closely
aligned with the spacecraft velocity, they may therefore have influenced the
properties of the assumed pristine plasma ahead of Cassini.

Different density distributions are also visible for the different species and
the electrons particularly show enhanced spatial variations around the
spacecraft. To either side behind the antenna dish in Figure 2a two distinct
regions of depleted flows are visible which appears to produce a vortex-type
structure in the x-z plane. This is analysed further in Section 3.5.

Figure 3: Currents onto the spacecraft within Saturn simulations. (a) shows
the spacecraft potential over time, (b) shows the current onto the spacecraft
over time, (c) shows the current onto the dish over time, and (d) shows the
current decomposed into the positive and negative parts.

###  3.2 Spacecraft Potential

In this Section we compare the simulations results to Cassini Langmuir Probe
measurements of the spacecraft potential and plasma currents during ingress on
Rev 292. Figure 3 shows the evolution of the spacecraft potential and plasma
currents through this simulation run for input parameters derived for when
Cassini was at a Saturn altitude of 2500 km in the northern hemisphere.

The net spacecraft potential is shown in panel (a) and incident plasma
currents in panel (b), both of which start near zero. As the spacecraft
charges, the main spacecraft body and the probe remain at a fixed bias
relative to one another of 3 V. The net current becomes strongly negative
before returning to zero when an equilibrium is reached after $\approx$1
$\mu$s. The spacecraft has at this point accumulated a negative floating
potential of just under –0.49 V with the probe biased at –3.49 V. The
simulation is, however, run for significantly longer to ensure steady state.
The simulated potential is –0.42 V, and 0.3 V more negative than the observed
potential of –0.12 V. The sensitivity of the spacecraft potential to
variations of this magnitude are analysed within following Sections 3.4.

Figure 4: Currents onto the simulated Langmuir Probe decomposed into the
positive and negative parts and compared to the Cassini Langmuir Probe
readings during Rev 292 at 2500 km altitude ingress.

The current decomposition for each of the electrons, positive and negative
ions onto each part of the spacecraft is shown in panels (c–d). The spike in
the negative current in panel (b) corresponds to electron current which also
peaks at the beginning of the simulation. The spacecraft surfaces accumulate
the associated negative charge and increasingly repel electrons away, thereby
reducing the electron current. These initial electron dynamics reflects how
the smaller electrons are much more sensitive to the change in the potential
of the spacecraft compared to the negative ions which, due to their larger
momentum, are not as easily deflected. At the equilibrium potential, all three
species of the plasma contribute comparable level of currents onto the
spacecraft relative to each other. This current balance is significantly
different to within typical electron-ion plasmas where the electrons
constitute the majority of the incident current [Miyake et al., 2020, and
references therein]. Comparing between the current received from the dish and
from the main body, the current composition and relative magnitudes are also
different. This is due to the difference in their surface areas but also the
deflection of the particle’s trajectories by the ambient magnetic field
resulting in increased quantities of ions and negative ions striking the sides
of the dish despite the incident flow being approximately parallel to it.

###  3.3 Langmuir Probe Currents

Figure 4 shows the currents onto the simulated Langmuir Probe. The positive
ion current stabilises after $\approx$ 10 $\mu$s at $\approx$ 125 nA, slightly
lower than the observed total current of $\approx$ 150 nA. However, when
summing the net current over all plasma species, the current is around 30 %
less. The initial analysis of the LP current was conducted assuming the total
current derived from positive ions and that the negative current was
negligible. As a result of the simulation it is now clear that the negative
ion current is not negligible and can be approximated as an almost constant
current, regardless of the potential difference, due to their inertia. These
results highlight that while the electron density can be accurately
determined, the total plasma density cannot due to the ion and negative ion
components. Therefore, by scaling the plasma concentration the net currents
onto the probe will increase. There could also be other effects at play such
as from the impact and break-up of large dust particles onto the probe and
spacecraft which are not included in these simulations.

###  3.4 Parametric Survey

To further understand the sensitivity of the results to the unknown or
inferred parameters, we carry out a parametric study and systematically varied
the following parameters; the electron depletion, the ion and negative ion
masses, the electron temperature and also the ion and negative ion
temperatures which were previously assumed to be in thermal equilibrium with
the electrons. The total plasma density, however, remains unchanged and the
densities are scaled within this. The positive ion mass is now set at 5 amu to
represent the measured lower bound in the deep ionosphere [Morooka et al.,
2019].

Figure 5: Parametric scan of electron density depletion for a range of
negative ion masses (a) and positive ion masses (b).

Figure 5(a) shows the spacecraft potential as a function of electron
depletion. The electron density is varied from 1 % up to 50 % of the total ion
density, for negative ion masses of 4, 6, 50 and 100 amu. For all negative ion
masses simulated, as the electron density tends to zero, the spacecraft
potential becomes less negative and eventually reaches positive values. For
larger negative ion masses the spacecraft become positive at slightly higher
electron densities although this change becomes increasingly smaller with the
difference between 50 and 100 amu being significantly smaller than that
between 4 and 50 amu.

The effect of a varying positive ion mass on the spacecraft potential is shown
in Figure 5b where the negative ion mass is held constant at 16 amu and the
spacecraft similarly tends to positive potentials for increased electron
depletions. Smaller positive ion masses produced positive potentials at lower
masses than larger ones and this change also becomes smaller for the larger
masses. The spacecraft potential also appears more sensitive to the positive
ion mass than to the negative ion mass. These trends are explained by the
variation in the relative mobilities of the positive and negative charge
carriers. The production of a positive potential is due to the positive ions
actually becoming the most mobile charge carriers which result in the
spacecraft accumulating a net positive current. Kim and Merlino [2006] indeed
predict that a body immersed in an electron-depleted negative ion plasma can
gain a positive potential when the electron density reaches just a few percent
of the total plasma density. Figure 5a demonstrates this is possible for the
Cassini spacecraft in Saturn’s ionosphere.

Figure 6: Parametric scan of ion (a) and electron (b) temperatures.

To study the sensitivity to the electron and ion temperatures, the negative
ion mass is held at 16 amu and positive ion mass at 5 amu but the depletion
rate is now set at 95 % to represent the parameter regime close to where
Cassini gains a positive potential. The ion temperature variation, shown in
Figure 6a, reveals that the spacecraft potential increases as the temperature
increases from the measured mean temperature of $\approx$0.1 eV, reaching
–0.21 V when Ti=1 eV. For temperatures below the measured 0.1 eV, the
spacecraft potential initially follows this trend but then, surprisingly,
increases again when T${}_{i}\leq$0.5 eV. This slight effect is attributed to
the ions becoming increasingly magnetised as the mean ion gyroradius decreases
until it is of a comparable size to the spacecraft. The electron temperature
variation, shown in Figure 6b, has a much larger impact on the spacecraft
potential due the electron’s high relative mobility. The potential
consequently varies from -0.05 V to -1.5 V, despite the electrons constituting
just a few percent of the total plasma density.

###  3.5 Potential gradient

The relative mobility of the positive and negative ions explains the positive
potential. However, in Figure 5a, the spacecraft still charges to a positive
potential when using negative ions lighter than positive ions (4 and 5 amu
respectively), so that the negatively charged particles remain the more mobile
charge carriers. This therefore appear to contradicts this explanation. Beyond
the variation in the plasma parameters examined in Section 3.4, the effect of
the ambient magnetic field remains to be evaluated.

Within the next simulations the magnetic field direction is oriented
perpendicularly to the plasma flow direction, instead along the +z and –z
axis. The influence of the magnetic field is felt in the generation of the
convective electric field, Ec, in the spacecraft frame, namely,

$\boldsymbol{E_{c}}=\boldsymbol{-v}_{flow}\times\boldsymbol{B_{0}},$ (2)

which produces a potential gradient along the main axis of the spacecraft
[Pecseli, 2012]. Figure 7 shows the electron density and potential within the
two simulations with the two different field orientation. The electrons can be
seen to be accumulating on one sides of the spacecraft and the potential
gradient across the spacecraft interaction is clearly visible. The potential
distribution controls where the electrons impact the spacecraft and alters the
effective spacecraft cross section due to Cassni’s asymmetric shape. The
spacecraft potential is consequently +0.38 V in the initial magnetic field
configuration but –0.42 V with the field reversed. Further tests run without
the presence of a magnetic field (not shown) also show that the spacecraft is
unable to attain a positive potential when the positive ions are less massive
than the negative ions.

Figure 7: Electron density (a & c) and potential (b & d) of Cassini immersed
in a plasma with a negative ion mass of 4 amu and positive ion mass of 5 amu.
The left-hand panels (a & b) show the case where the magnetic field is
oriented along Bz which produces a spacecraft potential of –0.42 V whereas the
right-hand panels (c & d) show the case with the magnetic field reversed where
the spacecraft potential is +0.38 V.

The perpendicular magnetic field enhances the magnetisation of the plasma
interaction. The differential electron flows around Cassini can be seen to
produce highly structured vortices at a scale that is comparable to Cassini
itself. It is interesting to note that electrostatic spacecraft-scale modes
have been identified in Saturn’s deep ionosphere which are thought to be
modulated by Cassini’s presence [Sulaiman et al., 2017].

##  4 Conclusions

This study simulated the spacecraft-plasma interaction experienced by the
Cassini within Saturn’s ionosphere during the Grand Finale using a self-
consistent three-dimensional PIC code. The asymmetric spacecraft shape was
taken into account and a Langmuir Probe was also modelled which enabled
preliminary comparisons between the simulated and observed currents.

The global interaction of the Cassini spacecraft was found to be strongly
influenced by its asymmetric shape. Differential flows around Cassini’s
antenna produced a highly structured wake and spacecraft-scale electron
vortices behind the spacecraft. Electron-wings associated with propagating
Langmuir waves were also identified around the spacecraft which produced a
sharp enhancements to the in-flowing electrons. Similar electron-wings were
previously identified in simulations of spacecraft moving through the polar
regions of the Earth’s ionosphere [Miyake et al., 2020] and these simulations
indicate this phenomena occurred at Cassini too.

An electron density of 40 % of the total negative charge density was initially
used in this simulations, which was determined to be a lower bound by
Cassini’s Langmuir Probe at 2500 km Saturn altitude during Rev 292 [Morooka et
al., 2019]. The spacecraft potential was found to be negative in the same
range as that observed at –0.42 V compared with the observed –0.12 V. The
simulations indicated that in Saturn’s ionosphere, the ion currents onto the
spacecraft and Langmuir Probe were comparable to those from the electrons, a
situation unique to plasmas depleted of electrons. The simulated currents onto
the Langmuir probe also revealed that the total currents were in the same
range, although slightly lower, than those measured. This is consistent with
the 40 % electron depletion being near the lower bound and indicates that the
actual positive and negative ion density may have been higher.

Following these initial simulations a parametric study was carried out to
examine how the floating potential is affected by the ambient conditions. To
this end, we studied how the spacecraft potential changes with positive and
negative ion mass and temperature as well as electron depletion rate. As the
electron depletion rate increased the spacecraft potential was observed to
become less negative and even attained positive potentials of up to 0.55 V
when just a few percent of electrons were contained within the plasma. Varying
the plasma masses and temperatures had a similar effect in terms of charge
mobility with the electron temperature found to have the largest influence.

The magnetic field orientation was also varied to examine how the induced
electric field along Cassini’s main axis affected the charge accumulated. The
direction of the magnetic field, therefore electric field, produces a strong
potential gradient in the potential of the plasma and results in oppositely
charge currents preferentially collecting at different parts of the
spacecraft. Given the large different in surface area across Cassini due to
its large antenna dish, reversing the magnetic field actually shifted the
spacecraft through zero from an initial 0.38 V to negative values of a
comparable magnitude of –0.42 V.

As well as providing an insight into the spacecraft plasma interactions
experienced by Cassini during the Grand Finale, these simulations have
provided an explanation for the positive potentials attained through
considering classical charging theory and reversed charge mobility. Despite
this, further effects such as secondary electron emission and the break-up of
large dust or ice particles on the spacecraft, and an even more detailed
Cassini model incorporating effects such as the curvature of the antenna dish,
remain to be constrained. The study does indicate, however, that this approach
presents a valuable method for understanding spacecraft charging in the outer
solar system and can assist in interpretting in-situ measurements of these
exotic environments.

## Acknowledgements

ZZ acknowledges an Undergraduate Summer Research Bursary from the Royal
Astronomical Society. RTD acknowledges funding from NERC grant NE/P017347/1.
YM and HU acknowledge grant 20K04041 from the Japan Society for the Promotion
of Science: JSPS, and support from the innovative High-Performance Computing
Infrastructure (HPCI: hp200032) in Japan. OS acknowledges RS grant RP EA180014
and SNSA grant Dnr:195/20. This work has benefited from discussions with
International Space Science Institute (ISSI) International Team 437. This work
used the Imperial College High Performance Computing Service (doi:
10.14469/hpc/2232).

## Data Availability

All simulation data presented in this study can be retrieved from the Zenodo
open-access repository at https://doi.org/10.5281/zenodo.4592954. Cassini
Langmuir Probe data is available from the NASA Planetary Data System archive.

## References

  * Al’pert et al. [1966] Y. L. Al’pert, A. V. Gurevich, and L. P. Pitaevskii.  Space physics with artificial satellites.  _American Journal of Physics_ , 34(6):544–544, 1966.  doi: 10.1119/1.1973101.  URL https://doi.org/10.1119/1.1973101.
  * Birdsall and Langdon [1991] C. K. Birdsall and A. B. Langdon.  _Plasma Physics via Computer Simulation_.  1991\.
  * Boris [1970] J. P. Boris.  A physically motivated solution of the alfven problem.  _Tech. Rep. 2167_ , 1970.
  * Coates et al. [2010] A. J. Coates, G. H. Jones, G. R. Lewis, A. Wellbrock, D. T. Young, F. J. Crary, R. E. Johnson, T. A. Cassidy, and T. W. Hill.  Negative ions in the Enceladus plume.  _Icarus_ , 206(2):618–622, Apr. 2010.  doi: 10.1016/j.icarus.2009.07.013.
  * Crary et al. [2009] F. J. Crary, B. A. Magee, K. Mandt, J. H. Waite, J. Westlake, and D. T. Young.  Heavy ions, temperatures and winds in Titan’s ionosphere: Combined Cassini CAPS and INMS observations.  _Planetary and Space Science_ , 57:1847–1856, Dec. 2009\.  doi: 10.1016/j.pss.2009.09.006.
  * Darian et al. [2017] D. Darian, S. Marholm, J. J. P. Paulsson, Y. Miyake, H. Usui, M. Mortensen, and W. J. Miloch.  Numerical simulations of a sounding rocket in ionospheric plasma: Effects of magnetic field on the wake formation and rocket potential.  _Journal of Geophysical Research (Space Physics)_ , 122(9):9603–9621, Sept. 2017.  doi: 10.1002/2017JA024284.
  * Desai et al. [2017a] R. T. Desai, A. J. Coates, A. Wellbrock, V. Vuitton, F. J. Crary, D. González-Caniulef, O. Shebanits, G. H. Jones, G. R. Lewis, J. H. Waite, M. Cordiner, S. A. Taylor, D. O. Kataria, J.-E. Wahlund, N. J. T. Edberg, and E. C. Sittler.  Carbon chain anions and the growth of complex organic molecules in titan’s ionosphere.  _The Astrophysical Journal Letters_ , 844(2):L18, 2017a.  URL http://stacks.iop.org/2041-8205/844/i=2/a=L18.
  * Desai et al. [2017b] R. T. Desai, M. M. Cowee, H. Wei, X. Fu, S. P. Gary, M. Volwerk, and A. J. Coates.  Hybrid simulations of positively and negatively charged pickup ions and cyclotron wave generation at europa.  _Journal of Geophysical Research: Space Physics_ , 122(10):10,408–10,420, 2017b.  ISSN 2169-9402.  doi: 10.1002/2017JA024479.  URL http://dx.doi.org/10.1002/2017JA024479.  2017JA024479.
  * Desai et al. [2018] R. T. Desai, S. A. Taylor, L. H. Regoli, A. J. Coates, T. A. Nordheim, M. A. Cordiner, B. D. Teolis, M. F. Thomsen, R. E. Johnson, G. H. Jones, M. M. Cowee, and J. H. Waite.  Cassini caps identification of pickup ion compositions at rhea.  _Geophysical Research Letters_ , pages n/a–n/a, 2018.  ISSN 1944-8007.  doi: 10.1002/2017GL076588.  URL http://dx.doi.org/10.1002/2017GL076588.  2017GL076588.
  * Drell et al. [1965] S. D. Drell, H. M. Foley, and M. A. Ruderman.  Drag and Propulsion of Large Satellites in the Ionosphere: An Alfvén Propulsion Engine in Space.  _Journ. Geophys. Res._ , 70(13):3131–3145, July 1965.  doi: 10.1029/JZ070i013p03131.
  * Engwall et al. [2006] E. Engwall, A. I. Eriksson, and J. Forest.  Wake formation behind positively charged spacecraft in flowing tenuous plasmas.  _Physics of Plasmas_ , 13(6):062904, 2006.  doi: 10.1063/1.2199207.  URL https://doi.org/10.1063/1.2199207.
  * Farrell et al. [2018] W. M. Farrell, L. Z. Hadid, M. W. Morooka, W. S. Kurth, J. E. Wahlund, R. J. MacDowall, A. H. Sulaiman, A. M. Persoon, and D. A. Gurnett.  Saturn’s Plasma Density Depletions Along Magnetic Field Lines Connected to the Main Rings.  _Geophys. Res. Lett._ , 45(16):8104–8110, Aug. 2018.  doi: 10.1029/2018GL078137.
  * Garrett [1981] H. B. Garrett.  The charging of spacecraft surfaces.  _Reviews of Geophysics_ , 19(4):577–616, 1981\.  doi: 10.1029/RG019i004p00577.  URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/RG019i004p00577.
  * Gurnett et al. [2004] D. A. Gurnett, W. S. Kurth, D. L. Kirchner, G. B. Hospodarsky, T. F. Averkamp, P. Zarka, A. Lecacheux, R. Manning, A. Roux, P. Canu, N. Cornilleau-Wehrlin, P. Galopeau, A. Meyer, R. Boström, G. Gustafsson, J. E. Wahlund, L. Åhlen, H. O. Rucker, H. P. Ladreiter, W. Macher, L. J. C. Woolliscroft, H. Alleyne, M. L. Kaiser, M. D. Desch, W. M. Farrell, C. C. Harvey, P. Louarn, P. J. Kellogg, K. Goetz, and A. Pedersen.  The Cassini Radio and Plasma Wave Investigation.  _Space Science Reviews_ , 114(1-4):395–463, Sept. 2004.  doi: 10.1007/s11214-004-1434-0.
  * Hadid et al. [2019] L. Z. Hadid, M. W. Morooka, J. E. Wahlund, A. M. Persoon, D. J. Andrews, O. Shebanits, W. S. Kurth, E. Vigren, N. J. T. Edberg, A. F. Nagy, and A. I. Eriksson.  Saturn’s Ionosphere: Electron Density Altitude Profiles and D-Ring Interaction From The Cassini Grand Finale.  _Geophys. Res. Lett._ , 46(16):9362–9369, Aug. 2019.  doi: 10.1029/2018GL078004.
  * Hockney and Eastwood [1981] R. W. Hockney and J. W. Eastwood.  _Computer Simulation Using Particles_.  1981\.
  * Holmberg et al. [2017] M. K. G. Holmberg, O. Shebanits, J. E. Wahlund, M. W. Morooka, E. Vigren, N. André, P. Garnier, A. M. Persoon, V. Génot, and L. K. Gilbert.  Density Structures, Dynamics, and Seasonal and Solar Cycle Modulations of Saturn’s Inner Plasma Disk.  _Journal of Geophysical Research (Space Physics)_ , 122(12):12,258–12,273, Dec. 2017.  doi: 10.1002/2017JA024311.
  * Kim and Merlino [2006] S.-H. Kim and R. L. Merlino.  Charging of dust grains in a plasma with negative ions.  _Physics of Plasmas_ , 13(5):052118, May 2006\.  doi: 10.1063/1.2204830.
  * Ludwig et al. [2012] P. Ludwig, W. J. Miloch, H. Kählert, and M. Bonitz.  On the wake structure in streaming complex plasmas.  _New Journal of Physics_ , 14(5):053016, May 2012\.  doi: 10.1088/1367-2630/14/5/053016.
  * Marklund et al. [1994] G. T. Marklund, L. G. Blomberg, P. A. Lindqvist, C. G. Fälthammar, G. Haerendel, F. S. Mozer, A. Pedersen, and P. Tanskanen.  The double probe electric field experiment on Freja: Experiment description and first results.  _Space Science Reviews_ , 70(3-4):483–508, Nov. 1994.  doi: 10.1007/BF00756883.
  * Miloch [2014] W. J. Miloch.  Numerical simulations of dust charging and wakefield effects.  _Journal of Plasma Physics_ , 80(6):795–801, Dec. 2014.  doi: 10.1017/S0022377814000300.
  * Miloch and Vladimirov [2009] W. J. Miloch and S. V. Vladimirov.  Charging of spinning insulating objects by plasma and photoemission.  _Geophysical Research Letters_ , 36(18), 2009.  doi: 10.1029/2009GL040337.  URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2009GL040337.
  * Miyake and Usui [2009] Y. Miyake and H. Usui.  New electromagnetic particle simulation code for the analysis of spacecraft-plasma interactions.  _Physics of Plasmas_ , 16(6):062904, 2009.  doi: 10.1063/1.3147922.  URL https://doi.org/10.1063/1.3147922.
  * Miyake et al. [2020] Y. Miyake, W. J. Miloch, S. H. Kjus, and H. L. Pécseli.  Electron Wing-Like Structures Formed at a Negatively Charged Spacecraft Moving in a Magnetized Plasma.  _Journal of Geophysical Research (Space Physics)_ , 125(2):e27379, Feb. 2020.  doi: 10.1029/2019JA027379.
  * Morooka et al. [2019] M. W. Morooka, J. E. Wahlund, L. Z. Hadid, A. I. Eriksson, N. J. T. Edberg, E. Vigren, D. J. Andrews, A. M. Persoon, W. S. Kurth, D. A. Gurnett, W. M. Farrell, J. H. Waite, R. S. Perryman, and M. Perry.  Saturn’s Dusty Ionosphere.  _Journal of Geophysical Research (Space Physics)_ , 124(3):1679–1697, Mar. 2019.  doi: 10.1029/2018JA026154.
  * Muralikrishna and Kulkarni [2006] P. Muralikrishna and V. H. Kulkarni.  On the height variation of the E-region cowling conductivity - effect of charged dust particles.  _Annales Geophysicae_ , 24(11):2949–2957, Nov. 2006.  doi: 10.5194/angeo-24-2949-2006.
  * Neuebauer [1980] F. M. Neuebauer.  Nonlinear standing Alfvén wave current system at Io: Theory.  _Journ. Geophys. Res._ , 85(A3):1171–1178, Mar. 1980.  doi: 10.1029/JA085iA03p01171.
  * Nordheim et al. [2020] T. A. Nordheim, A. Wellbrock, G. H. Jones, R. T. Desai, A. J. Coates, B. D. Teolis, and J. M. Jasinski.  Detection of Negative Pickup Ions at Saturn’s Moon Dione.  _Geophys. Res. Lett._ , 47(7):e87543, Apr. 2020\.  doi: 10.1029/2020GL087543.
  * Olson et al. [2010] J. Olson, W. J. Miloch, S. Ratynskaia, and V. Yaroshenko.  Potential structure around the Cassini spacecraft near the orbit of Enceladus.  _Physics of Plasmas_ , 17(10):102904–102904, Oct. 2010.  doi: 10.1063/1.3486523.
  * Pecseli [2012] H. L. Pecseli.  _Waves and Oscillations in Plasmas_.  Series in Plasma Physics. CRC Press, Hoboken, NJ, 2012.  URL https://cds.cern.ch/record/1488873.
  * Persoon et al. [2019] A. M. Persoon, W. S. Kurth, D. A. Gurnett, J. B. Groene, A. H. Sulaiman, J. E. Wahlund, M. W. Morooka, L. Z. Hadid, A. F. Nagy, J. H. Waite, and T. E. Cravens.  Electron Density Distributions in Saturn’s Ionosphere.  _Geophys. Res. Lett._ , 46(6):3061–3068, Mar. 2019.  doi: 10.1029/2018GL078020.
  * Rehman and Marchand [2014] S. Rehman and R. Marchand.  Plasma-satellite interaction driven magnetic field perturbations.  _Physics of Plasmas_ , 21(9):090701, 2014.  doi: 10.1063/1.4894678.  URL https://doi.org/10.1063/1.4894678.
  * Robertson [2013] S. Robertson.  Sheaths in laboratory and space plasmas.  _Plasma Physics and Controlled Fusion_ , 55(9):093001, Sept. 2013.  doi: 10.1088/0741-3335/55/9/093001.
  * Roussel and Berthelier [2004] J.-F. Roussel and J.-J. Berthelier.  A study of the electrical charging of the rosetta orbiter: 1. numerical model.  _Journal of Geophysical Research: Space Physics_ , 109(A1), 2004.  doi: 10.1029/2003JA009836.  URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2003JA009836.
  * Shebanits et al. [2013] O. Shebanits, J.-E. Wahlund, K. Mandt, K. Ågren, N. J. T. Edberg, and J. H. Waite.  Negative ion densities in the ionosphere of Titan-Cassini RPWS/LP results.  _Planetary and Space Science_ , 84:153–162, Aug. 2013.  doi: 10.1016/j.pss.2013.05.021.
  * Shebanits et al. [2016] O. Shebanits, J.-E. Wahlund, N. J. T. Edberg, F. J. Crary, A. Wellbrock, D. J. Andrews, E. Vigren, R. T. Desai, A. J. Coates, K. E. Mandt, and J. H. Waite.  Ion and aerosol precursor densities in Titan’s ionosphere: A multi-instrument case study.  _Journal of Geophysical Research (Space Physics)_ , 121(A10):10, Oct. 2016.  doi: 10.1002/2016JA022980.
  * Shebanits et al. [2017] O. Shebanits, E. Vigren, J. E. Wahlund, M. K. G. Holmberg, M. Morooka, N. J. T. Edberg, K. E. Mand t, and J. H. Waite.  Titan’s ionosphere: A survey of solar EUV influences.  _Journal of Geophysical Research (Space Physics)_ , 122(7):7491–7503, July 2017.  doi: 10.1002/2017JA023987.
  * Shebanits et al. [2020] O. Shebanits, L. Z. Hadid, H. Cao, M. W. Morooka, G. J. Hunt, M. K. Dougherty, J. E. Wahlund, J. H. Waite, and I. Müller-Wodarg.  Saturn’s near-equatorial ionospheric conductivities from in situ measurements.  _Scientific Reports_ , 10:7932, May 2020.  doi: 10.1038/s41598-020-64787-7.
  * Shukla and Mamun [2002] P. K. Shukla and A. A. Mamun.  _Introduction to dusty plasma physics_.  2002\.
  * Spitzer [1941] J. Spitzer, Lyman.  The Dynamics of the Interstellar Medium. I. Local Equilibrium.  _Astrophys. Journ. Lett._ , 93:369, May 1941.  doi: 10.1086/144273.
  * Stenzel and Urrutia [1989] R. L. Stenzel and J. M. Urrutia.  Whistler wings from moving electrodes in a magnetized laboratory plasma.  _Geophys. Res. Lett._ , 16(5):361–364, May 1989\.  doi: 10.1029/GL016i005p00361.
  * Sulaiman et al. [2017] A. H. Sulaiman, W. S. Kurth, A. M. Persoon, J. D. Menietti, W. M. Farrell, S. Y. Ye, G. B. Hospodarsky, D. A. Gurnett, and L. Z. Hadid.  Intense Harmonic Emissions Observed in Saturn’s Ionosphere.  _Geophys. Res. Lett._ , 44(24):12,049–12,056, Dec. 2017.  doi: 10.1002/2017GL076184.
  * Usui et al. [2019] H. Usui, Y. Miyake, W. J. Miloch, and K. Ito.  Numerical study of plasma depletion region in a satellite wake.  _IEEE Transactions on Plasma Science_ , 47(8):3717–3723, 2019.  doi: 10.1109/TPS.2019.2918789.
  * Wahlund et al. [2005] J. E. Wahlund, R. Boström, G. Gustafsson, D. A. Gurnett, W. S. Kurth, A. Pedersen, T. F. Averkamp, G. B. Hospodarsky, A. M. Persoon, P. Canu, F. M. Neubauer, M. K. Dougherty, A. I. Eriksson, M. W. Morooka, R. Gill, M. André, L. Eliasson, and I. Müller-Wodarg.  Cassini Measurements of Cold Plasma in the Ionosphere of Titan.  _Science_ , 308(5724):986–989, May 2005.  doi: 10.1126/science.1109807.
  * Wahlund et al. [2018] J. E. Wahlund, M. W. Morooka, L. Z. Hadid, A. M. Persoon, W. M. Farrell, D. A. Gurnett, G. Hospodarsky, W. S. Kurth, S. Y. Ye, D. J. Andrews, N. J. T. Edberg, A. I. Eriksson, and E. Vigren.  In situ measurements of Saturn’s ionosphere show that it is dynamic and interacts with the rings.  _Science_ , 359(6371):66–68, Jan. 2018.  doi: 10.1126/science.aao4134.
  * Waite et al. [2018] J. H. Waite, R. S. Perryman, M. E. Perry, K. E. Miller, J. Bell, T. E. Cravens, C. R. Glein, J. Grimes, M. Hedman, J. Cuzzi, T. Brockwell, B. Teolis, L. Moore, D. G. Mitchell, A. Persoon, W. S. Kurth, J. E. Wahlund, M. Morooka, L. Z. Hadid, S. Chocron, J. Walker, A. Nagy, R. Yelle, S. Ledvina, R. Johnson, W. Tseng, O. J. Tucker, and W. H. Ip.  Chemical interactions between Saturn’s atmosphere and its rings.  _Science_ , 362(6410):aat2382, Oct. 2018.  doi: 10.1126/science.aat2382.
  * Whipple [1981] E. C. Whipple.  Potentials of surfaces in space.  _Reports on Progress in Physics_ , 44:1197–1250, Nov. 1981\.  doi: 10.1088/0034-4885/44/11/002.
  * Willis et al. [2011] C. T. N. Willis, J. E. Allen, M. Coppins, and M. Bacharis.  Wakes formed by dust grains in supersonically flowing plasmas.  _Physical Review E_ , 84(4):046410, Oct. 2011\.  doi: 10.1103/PhysRevE.84.046410.
  * Yaroshenko et al. [2011] V. V. Yaroshenko, W. J. Miloch, S. Vladimirov, H. M. Thomas, and G. E. Morfill.  Modeling of cassini’s charging at saturn orbit insertion flyby.  _Journal of Geophysical Research: Space Physics_ , 116(A12), 2011.  doi: 10.1029/2011JA016775.  URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2011JA016775.
  * Yee [1966] K. S. Yee.  Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media.  _IEEE Trans. Antennas and Propagation_ , pages 302–307, 1966.


11institutetext: INAF - Osservatorio Astronomico di Capodimonte, via
Moiariello 16, 80131 Napoli, Italy 22institutetext: INFN – Section of Naples,
via Cinthia 9, 80126 Napoli, Italy 33institutetext: Department of Physics,
University of Napoli “Federico II”, via Cinthia 9, 80126 Napoli, Italy
44institutetext: Millennium Institute of Astrophysics (MAS), Nuncio Monseñor
Sotero Sanz 100, Providencia, Santiago, Chile 55institutetext: AIMI, ARTORG
Center, University of Bern, Murtenstrasse 50, Bern, CH-3008, Switzerland
66institutetext: Department of Physics and Astronomy ‘Augusto Righi’,
University of Bologna, via Piero Gobetti 93/2, I-40129, Bologna, Italy;
77institutetext: INAF — Osservatorio di Astrofisica e Scienza dello Spazio di
Bologna, via Piero Gobetti 101, I-40129, Bologna, Italy;

# Identification of problematic epochs in astronomical time series through
transfer learning

Stefano Cavuoti , The first two authors share the first authorship and are in
alphabetical order.
e-mail addresses:
<EMAIL_ADDRESS>and<EMAIL_ADDRESS>0000-0002-3787-4196
Demetra De Cicco,⋆ 331144 0000-0001-7208-5101 Lars Doorenbos 55
0000-0002-0231-9950
Massimo Brescia 331122 0000-0001-9506-5680 Olena Torbaniuk 6677
0000-0003-4465-2564 Giuseppe Longo 33 0000-0002-9182-8414 Maurizio Paolillo
331122 0000-0003-4210-7693

(Received March 28, 2024; accepted May 06, 2024)

###### Abstract

Aims. We present a novel method for detecting outliers in astronomical time
series based on the combination of a deep neural network and a k-nearest
neighbor algorithm with the aim of identifying and removing problematic epochs
in the light curves of astronomical objects.

Methods. We used an EfficientNet network pretrained on ImageNet as a feature
extractor and performed a k-nearest neighbor search in the resulting feature
space to measure the distance from the first neighbor for each image. If the
distance was above the one obtained for a stacked image, we flagged the image
as a potential outlier.

Results. We applied our method to a time series obtained from the VLT Survey
Telescope monitoring campaign of the Deep Drilling Fields of the Vera C. Rubin
Legacy Survey of Space and Time ††thanks: Observations were provided by the
ESO programs 088.D-4013, 092.D-0370, and 094.D-0417 (PI G. Pignata).. We show
that our method can effectively identify and remove artifacts from the VST
time series and improve the quality and reliability of the data. This approach
may prove very useful in light of the amount of data that will be provided by
the LSST, which will prevent the inspection of individual light curves. We
also discuss the advantages and limitations of our method and suggest possible
directions for future work.
The code is available here.

###### Key Words.:

Machine learning – Galaxies: active – Techniques: image processing – Methods:
data analysis – Artificial Intelligence

##  1 Introduction

Astronomical time series are sequences of data that represent the variation of
a physical quantity over time in the observation of celestial objects or
phenomena. They are essential for studying and discovering the properties and
dynamics of various astrophysical sources, such as stars, planets, supernovae,
black holes, and galaxies (Scargle, 1997; Aigrain & Foreman-Mackey, 2023).

However, an astronomical time series analysis poses several challenges, such
as the presence of noise, gaps, outliers, and artifacts in the data. Artifacts
are unwanted signals that do not reflect the true nature of the observed
object and are caused by external factors, such as instrumental errors,
atmospheric conditions, peculiar events (e.g., cosmic rays, satellite tracks),
and contaminant sources (Li et al., 2018; Malz et al., 2019). Artifacts can
severely affect the quality and reliability of data and may lead to false or
misleading conclusions. Therefore, it is crucial to detect and remove
artifacts from astronomical time series before performing any further
analysis.

One specific case of astronomical time series are those obtained from the VLT
Survey Telescope (VST), a 2.6-meter optical telescope located at the Paranal
Observatory in Chile. The VST is designed to perform wide-field imaging
surveys of the southern sky, covering a field of view (FoV) of one square
degree with a pixel scale of 0.214″ (Capaccioli & Schipani, 2011). The VST
time series are useful for studying various topics, such as variable stars,
transient events, Solar System objects, and cosmology (Cappellaro et al.,
2015; Falocco et al., 2015; De Cicco et al., 2015, 2019, 2021; Botticella et
al., 2017; Fu et al., 2018; Liu et al., 2020; Poulain et al., 2020).

The VST time series, as all astronomical data, are affected by artifacts, such
as bad pixels, saturation, cosmic rays, ghost images, and background
fluctuations. These artifacts can hamper the detection and characterization of
faint or fast-varying sources and decrease the scientific value of the data.
Data reduction and calibration methods (see, e.g., Grado et al. 2012) are
expected to correct and mask such artifacts or, when this proves difficult,
flag them to allow easy filtering in the data analysis phase. However, it is
not unusual to find residual defects in the reduced data, especially in time
series, where we cannot use time averaging to remove many of the
aforementioned problems. Therefore, it is important to develop improved
methods for identifying and flagging these artifacts in order to limit their
impact on scientific results.

One possible way to tackle this problem is to use transfer learning, a
technique in machine learning that allows one to leverage the knowledge
learned from a source domain to improve the performance on a target domain.
Transfer learning can be useful when the target domain has limited or scarce
data or when the source domain has some similarities or relevance to the
target domain, as happens in astronomy (at least in terms of labels; Awang
Iskandar et al. 2020; Martinazzo et al. 2021), malware classification (Prima &
Bouhorma, 2020), Earth science (Zou & Zhong, 2018), and medicine (Ding et al.,
2019; Esteva et al., 2017; Kim et al., 2021; Menegola et al., 2017).

In this work, we propose the use of an EfficientNet network trained on
ImageNet as an anomaly detector for astronomical time series. An EfficientNet
network (Tan & Le, 2019) is a convolutional neural network (CNN) architecture
that scales up the network depth, width, and resolution in a balanced and
efficient way using a compound coefficient. ImageNet is a large-scale dataset
often used for such a purpose (Deng et al., 2009), and it contains around 1.3
million images, where the original task was to classify each image into one of
1,000 classes.

The idea is to use the EfficientNet network as a feature extractor, as
proposed in Doorenbos et al. (2022), in the detection of active galactic
nuclei (AGN) and to then use an algorithm for the detection of sources
deviating from the standard behavior.

This approach can improve and automatize the identification and removal of
problematic epochs in the light curves without the need for visual inspection.
This process can therefore improve the reliability of the light curves and the
results obtained with their analysis. We expect this to be highly relevant in
light of the Vera C. Rubin Legacy Survey of Space and Time (LSST; Ivezić et
al. 2019).

The paper is organized as follows: In Sec. 2, we describe the data used for
this work, while in Sec. 3, we present the method. In Sec. 4, we discuss the
results, and finally in Sec. 5 we draw our conclusions.

##  2 Data

As a benchmark, we used the same dataset as De Cicco et al. (2019, 2021,
2022). The dataset consists of $r$-band observations of the COSMOS field
performed with the VST spanning three observing seasons (hereafter, seasons)
from December 2011 to March 2015 and are part of a long-term effort to monitor
the LSST Deep Drilling Fields before the start of the Vera Rubin Telescope
operations. The VST is a 2.6-meter optical telescope that covers a FoV of
$1^{\circ}\times 1^{\circ}$ with a single pointing (the pixel scale is
$0.214$″/pixel). The three seasons include 54 visits in total and with two
gaps; more details about the dataset are given in Table 1 (adapted from Table
1 of De Cicco et al. (2019)). As explained in De Cicco et al. (2019), they
excluded 11 of the 65 visits that constitute the full dataset.

Table 1: VST-COSMOS dataset used in this work. The four columns report the
visit number, OB identification number, observing date, and seeing FWHM,
respectively, for the 54 visits used for the present analysis; visits are
listed in chronological order. All the visits were obtained by combining five
exposures, for a total exposure time of 1,800 s. However, visit 53 was
obtained by combining ten exposures, for a total exposure time of 3,600 s.
Adapted from Table 1 of De Cicco et al. (2019).

visit | OB-ID | obs. date | seeing (FWHM)
---|---|---|---
|  |  | (arcsec)
1 | 611279 | 2011-Dec-18 | $0.64$
2 | 611283 | 2011-Dec-22 | $0.94$
3 | 611287 | 2011-Dec-27 | $1.04$
4 | 611291 | 2011-Dec-31 | $1.15$
5 | 611295 | 2012-Jan-02 | $0.67$
6 | 611299 | 2012-Jan-06 | $0.58$
7 | 611311 | 2012-Jan-18 | $0.62$
8 | 611315 | 2012-Jan-20 | $0.88$
9 | 611319 | 2012-Jan-22 | $0.81$
10 | 611323 | 2012-Jan-24 | $0.67$
11 | 611327 | 2012-Jan-27 | $0.98$
12 | 611331 | 2012-Jan-29 | $0.86$
13 | 611335 | 2012-Feb-02 | $0.86$
14 | 611351 | 2012-Feb-16 | $0.50$
15 | 611355 | 2012-Feb-19 | $0.99$
16 | 611359 | 2012-Feb-21 | $0.79$
17 | 611363 | 2012-Feb-23 | $0.73$
18 | 611367 | 2012-Feb-26 | $0.83$
19 | 611371 | 2012-Feb-29 | $0.90$
20 | 611375 | 2012-Mar-03 | $0.97$
21 | 611387 | 2012-Mar-13 | $0.70$
22 | 611391 | 2012-Mar-15 | $1.08$
23 | 611395 | 2012-Mar-17 | $0.91$
24 | 768813 | 2012-May-08 | $0.74$
25 | 768817 | 2012-May-11 | $0.85$
26 | 768820 | 2012-May-17 | $0.77$
27 | 986611 | 2013-Dec-27 | $0.72$

visit | OB-ID | obs. date | seeing (FWHM)
---|---|---|---
|  |  | (arcsec)
28 | 986614 | 2013-Dec-30 | $1.00$
29 | 986617 | 2014-Jan-03 | $0.86$
30 | 986620 | 2014-Jan-05 | $0.81$
31 | 986626 | 2014-Jan-12 | $0.73$
32 | 986630 | 2014-Jan-21 | $1.18$
33 | 986633 | 2014-Jan-24 | $0.80$
34 | 986648 | 2014-Feb-09 | $1.28$
35 | 986652 | 2014-Feb-19 | $0.89$
36 | 986655 | 2014-Feb-21 | $0.93$
37 | 986658 | 2014-Feb-23 | $0.81$
38 | 986661 | 2014-Feb-26 | $0.81$
39 | 986664 | 2014-Feb-28 | $0.77$
40 | 986670 | 2014-Mar-08 | $0.91$
41 | 986674 | 2014-Mar-21 | $0.96$
42 | 986677 | 2014-Mar-23 | $0.92$
43 | 986680 | 2014-Mar-25 | $0.66$
44 | 1095777 | 2014-Mar-29 | $0.89$
45 | 1095783 | 2014-Apr-04 | $0.58$
46 | 986683 | 2014-Apr-07 | $0.61$
47 | 1136410 | 2014-Dec-03 | $1.00$
48 | 1136457 | 2015-Jan-10 | $0.71$
49 | 1136481 | 2015-Jan-28 | $0.90$
50 | 1136490 | 2015-Jan-31 | $0.73$
51 | 1136503 | 2015-Feb-15 | $0.70$
52 | 1136531 | 2015-Mar-10 | $0.80$
53 | 1136540 | 2015-Mar-14 | $0.84$
54 | 1136543 | 2015-Mar-19 | $1.00$
stacked | - | - | $0.67$

The observations in the $r$-band were designed to have a three-day observing
cadence, although the exact final cadence depended on the observational
constraints. The single-visit depth is $r\lesssim 24.6$ mag for point sources,
at a ${\sim}5\sigma$ confidence level. This makes our dataset particularly
interesting for studies aimed at forecasting the performance of the LSST, as
its single-visit depth is expected to be similar to our stacked COSMOS images.
We refer the reader to De Cicco et al. (2015) for information on the reduction
and combination of the exposures performed using the VST-Tube pipeline (Grado
et al., 2012) as well as the source extraction and sample assembly. The VST-
Tube magnitudes are in the AB system.

The sample contains 22,927 sources detected in at least 50% of the visits in
the dataset (i.e., they have at least 27 points in their light curves) with an
average magnitude of $r\leq 23.5$ mag within a $1\arcsec$-radius aperture. We
started with the same sample for our analysis, but while De Cicco et al.
(2019) focused on the sub-sample of sources selected by either variability or
multiwavelength properties, we used the entire sample, as done in De Cicco et
al. (2021). In fact, our aim is to efficiently identify the unreliable
measurements in the full dataset before engaging in any classification effort
that could be biased by photometric or aesthetic artifacts. These problematic
sources are mainly satellite tracks, aesthetic defects, cosmic rays,
saturation, and bad seeing together with objects blended with a neighbor in at
least some of the visits (those with poor seeing), making it substantially
harder to measure fluxes correctly. We decided to analyze and compare this
sample, as it comes from the pipeline and dataset already presented in De
Cicco et al. (2021).

As the reference for each object, we used the stacked image, that is, the
median of all the exposures in the first season – which covers a five-month
baseline – with a seeing full width at half maximum (FWHM) of less than
$0.80\arcsec$. This choice ensured that the reference image is of higher
quality than each individual epoch and that it does not include all the
defects present in the whole dataset. A visual representation of the resulting
image from single epochs is shown in Fig. 1.

${54}$${42}$${37}$${61}$${36}$${71}$${91}$${96}$${31}$${32}$${120}$${65}$${89}$${110}$${45}$${83}$${93}$${105}$${105}$${36}$${62}$${47}$${38}$${63}$${55}$Image
1${14}$${50}$${27}$${73}$${86}$${71}$${81}$${253}$${33}$${45}$${73}$${75}$${250}$${103}$${78}$${37}$${94}$${101}$${96}$${55}$${26}$${73}$${45}$${96}$${68}$Image
2${35}$${26}$${8}$${53}$${66}$${55}$${71}$${102}$${21}$${29}$${65}$${58}$${88}$${85}$${92}$${16}$${90}$${68}$${106}$${74}$${9}$${51}$${50}$${75}$${81}$Image
3$\cdots$${34}$${32}$${37}$${63}$${36}$${53}$${91}$${116}$${18}$${16}$${118}$${95}$${87}$${113}$${45}$${93}$${99}$${110}$${103}$${26}$${72}$${37}$${76}$${76}$${61}$Image
N-1${31}$${46}$${56}$${77}$${15}$${50}$${108}$${108}$${15}$${28}$${111}$${89}$${79}$${118}$${61}$${103}$${103}$${92}$${83}$${10}$${75}$${38}$${65}$${96}$${69}$Image
N$\Longrightarrow$${34}$${42}$${37}$${63}$${36}$${53}$${91}$${96}$${31}$${29}$${111}$${75}$${88}$${118}$${61}$${83}$${94}$${101}$${103}$${36}$${62}$${47}$${50}$${76}$${68}$Stacked
Image Figure 1: Example derivation of our stacked image starting from
individual images. The individual images are represented by simple 5x5
matrices. The median value for each cell is selected in order to derive the
stacked image. This process is suitable to avoiding wrong values in single
epochs. For example, the values above 250 in Image 2 have been removed
completely, as when using the mean value, the presence of such values,
although mitigated by the number of independent images, would affect the final
result much more.

##  3 Method

We propose a novel method for detecting outliers in astronomical time series
based on combining a deep neural network and a k-nearest neighbor (k-NN)
algorithm. The preliminary phase is similar to what has been done in Doorenbos
et al. (2022). The main idea is to use an EfficientNet, a CNN (Tan & Le,
2019), pretrained on ImageNet as a feature extractor and to perform a k-NN
search in the resulting feature space to measure the distance to the first
neighbor for each image. If the distance is above a certain threshold, we flag
the image as a potential outlier.

Convolutional neural networks Lecun et al. 1990; LeCun et al. 2015 are a class
of deep neural networks designed for processing structured grid data, such as
images. They have been proven to be highly effective in computer vision tasks,
including image classification, object detection, and segmentation. A typical
CNN consists of multiple layers, including:

  * •

Convolutional layers: These layers apply convolution operations to their
input, using learnable filters to extract features. Convolutional layers help
capture hierarchical patterns in the input.

  * •

Pooling layers: Pooling layers reduce the spatial dimensions of the input by
downsampling. Common pooling operations include max pooling, which retains the
maximum value in a region, and average pooling, which computes the average.

  * •

Fully connected layers: These layers connect every neuron to every neuron in
the previous and subsequent layers, and they help make predictions based on
the extracted features.

  * •

Activation functions: Non-linear activation functions, such as a rectified
linear unit (ReLU), introduce non-linearity, allowing the network to learn
complex relationships.

An example of the typical architecture of a CNN is shown in Fig. 2.

input image
---
layer
convolutional layer
---
pooling layer
---
convolutional layer
---
pooling layer
---
fully connected layer
---
fully connected
---
output layer
Figure 2: Architecture of a traditional CNN in which there is a sequence of
convolutional and pooling layers, and after a final subsampling layer, there
is a fully connected one that performs the classification or regression task.
In this specific work, we removed the fully connected part, and we used only
the feature extraction, which is performed up to the last pooling layer.

Scaling in CNNs involves adjusting various architectural dimensions to control
the model complexity and computational efficiency. In particular, classical
CNNs perform one of the following kinds of scaling:

  * •

Depth scaling: Increasing the depth of a CNN involves adding more layers.
Deeper networks can capture more complex features but require more
computational resources.

  * •

Width scaling: Wider networks increase the number of filters in each layer.
This helps in capturing more diverse features but also leads to higher
computational requirements.

  * •

Resolution scaling: Resolution scaling involves adjusting the input image
size. Lower resolutions reduce computational demands but may result in loss of
information.

The EfficientNet network is a CNN architecture that scales up the network
depth, width, and resolution in a balanced and efficient way. This method is
called compound scaling, and it combines depth, width, and resolution scaling
in a balanced manner, maintaining a constant ratio. It achieved state-of-the-
art accuracy on ImageNet, a large-scale image dataset comprising more than 1.3
million images belonging to more than 1,000 categories. By using a pretrained
EfficientNet network, we can leverage the general image processing skills
learned from the ImageNet dataset and adapt them to the specific
characteristics of astronomical time series. For this work, we made use only
of the first part of the network that extracts the features until the last
pooling layer, while we completely discarded the final fully connected layer
that in EfficientNet works as classifier (see Fig. 2).

The k-NN (Cover & Hart, 1967) algorithm is a non-parametric supervised
learning method that classifies an object based on the label of its k closest
neighbors through a majority vote. In our case, we represented images by their
features extracted from the EfficientNet network and measured the Euclidean
distance of an image to its closest neighbor in this feature space (k=1). This
distance reflects the similarity to the other images, and a large distance
indicates a possible anomaly.

Since we needed a threshold to distinguish a real anomaly from something that
is simply less similar, we used a stacked image (see Sec. 2) as a reference
for the distance calculation. The stacked image is one that minimizes the
noise and the presence of defects or spurious features with respect to the
individual epochs, and it results in a sharper and cleaner reference template.
We could thus safely assume that the stacked image represents the normal state
of the scene and that any image that deviates significantly from it is an
outlier.

In practice, as depicted in Fig. 3, our method consists of the following
steps:

  1. a.

Input data preparation: Time series data from astronomical observations
consisting of multiple images of the same celestial object and its stacked
image are collected.

  2. b.

Feature extraction: A pre-trained EfficientNet network on ImageNet is used.
EfficientNet is a CNN architecture that efficiently scales depth, width, and
resolution. The features are collected from input images using the
EfficientNet as a feature extractor.

  3. c.

Identification of first neighbor: The k-NN algorithm (k=1) is applied to
measure the Euclidean distance among features extracted by the EfficientNet
network.

  4. d.

Thresholding: The threshold for anomaly detection is identified. We use a
stacked image as a reference for distance calculation.

  5. e.

Anomaly detection: If the distance from the first neighbor is above the
predefined threshold, the image is flagged as a potential outlier. We identify
images deviating significantly from the normal state represented by the
stacked image.

  6. f.

Output: A list of potential outlier images identified during the anomaly
detection process is produced.

Input Data Feature Extraction (EfficientNet) Identification of 1st Neighbor
distance (k-NN Algorithm) Thresholding with respect to stacked image Anomaly
Detection Output Figure 3: Schematic description of the algorithm.

##  4 Experiments and results

We decided to perform two testing campaigns, one on a larger dataset in order
to have a wider case study and one on a smaller dataset, as we could then
manually verify the results in a more detailed way. We decided to use as the
larger dataset the Unlabel Set and as the smaller dataset the Label Set, both
derived from De Cicco et al. (2021).

###  4.1 Large dataset

We conducted a thorough analysis of the results obtained from applying our
algorithm to the Unlabel Set dataset by De Cicco et al. (2021). The dataset
comprises 17,995 sources observed multiple times (Table 1), yielding a total
of 989,725 individual observations.

As a first step, we used the same procedure of De Cicco et al. (2019) to
expunge bad epochs from the light curves, applying a sigma-clipping algorithm
with a threshold of five times the value of the normalized median absolute
deviation (NMAD) from the median magnitude of each source. This conservative
threshold enabled us to exclude extreme outliers, preserving the most real
variable events of astrophysical interest, such as AGN, supernovae, and other
transients, provided the transient was not an extremely bright event in an
otherwise quiescent galaxy (of which we have found no evidence in our data so
far). This process refined the dataset to 914,156 epochs. Despite this initial
cleaning step, the final light curves can still be affected by many
photometric problems and artifacts, as discussed in Sec. 1, that may require
additional ad hoc refinements to obtain fully reliable light curves. In our
case, we directly applied the algorithm described in Sec. 3, finding 336
additional anomalous epochs.

Upon initial analysis, we observed that a specific epoch, namely epoch 32,
according to the nomenclature provided in Table 1, exhibited 44 problematic
images. A swift inspection revealed them to be images with a low S/N, as this
is one of the epochs with the worst seeing values (1.18″, see Table 1, Fig. 4,
14).

Figure 4: Examples of problematic objects from epoch 32. The full set of
sources flagged as problematic in epoch 32 is available in Fig. 14.

On four occasions, epoch 34 produced outliers linked to a low S/N ratio (see
Fig. 5). Epoch 34 actually has a worse seeing (i.e., 1.28″) than Epoch 32.

Figure 5: Sources flagged as outliers in epoch 34 due to low S/N ratio with
regard to the reference epoch.

Another noteworthy case concerns source 23560, where nine out of 54 epochs
were identified as anomalous, likely due to the contamination from a bright
neighbor (see Fig. 6). In fact, source 23560 has a bright companion that can
be seen in the bottom panel of Fig. 6, so we are in the presence of a
configuration that, due also to the fluctuation in the point spread function,
affects the companion outskirts that contaminate this source and is spotted by
our algorithm.

Figure 6: Thumbnails of the flagged epochs for source 23560 (including epoch
32, which was excluded from Fig. 14). The last panel is a broader view of the
stacked image showing the nearby contaminant; the size of the previous
thumbnails is shown by the red dashed square.

Twelve additional epochs showed contamination because of defects from
electronics or saturated neighbors. Even when the central object flux does not
seem directly affected, the algorithm background estimation might be distorted
(Fig. 7).

Figure 7: Objects with a bright source nearby. Each row corresponds to a
different object with some flagged epochs, with the exception of the last two
lines, as they correspond to the same object.

Additionally, we identified 185 instances of astronomical tracks (asteroids,
satellites, etc.) that, while not directly influencing the central source
flux, could distort background estimation and thus the magnitude determination
(Figs. 8, 15).

Figure 8: Some examples of objects with tracks. (See Fig. 15 for the whole
sample.)

In three cases, an epoch was flagged due to the appearance of a secondary
object in a single epoch (Fig. 9). While this is an intriguing collateral
study, it is unlikely to impact variability measurements; nonetheless, the
algorithm correctly identified them as outliers.

Figure 9: Sources with an additional object that appears nearby. For each
object, the epoch flagged by the method is shown in the central column. The
epoch before and the one after are also shown in the first and last column,
respectively. For the first object, the new object detection corresponds to
the faint object located approximately at (10,25).

The most significant finding pertains to 52 events where source measurements
were undoubtedly influenced to some extent by tracks or other issues (Figs.
10, 19, 20). These epochs should be excluded from any further study involving
these objects. In particular, we would like to emphasize that, with some
exceptions, those epochs could not be spotted with a different threshold of
the sigma clip since they lie in the core of the distribution of the light
curve despite being spoiled. There are some cases, such as the last line of
Fig. 10, where although there is an additional flux coming from a track, the
measured flux from the pipeline is lower. In this specific case, we observed
that due to the track, the pipeline provided a wrong estimate of the centroid
of the object, and therefore the aperture magnitude does not contain all the
flux of the actual source. Since our method does not rely on centroid
measurements but makes use of the original images, it is able to identify
these outliers anyway.

Figure 10: Some of the problematic objects identified by our method. In most
cases, if a track appears, it is closer to the object with respect to the
aperture photometry, so it slightly affects the flux measurement. In the left
column, one can see the image as identified by our algorithm, while in the
right we report the light curve with the flagged epoch marked with red. The
horizontal lines represent the threshold of the sigma clip. (See Fig. 19 and
20 for the whole sample.)

Finally there remain 27 epochs for which, despite being flagged as outliers, a
clear interpretation could not be found. They represent approximately 8% of
the flagged sources (that corresponds to 0.003% of the whole set of images
presented to the algorithm) and could be considered as false positives
identified by our algorithm (Fig. 11, 21, and 22). However, we observed that a
significant portion (approximately 22% of these cases) belonged to one of the
two epochs with the best seeing conditions (epochs 14 and 6, with seeing of
0.5″and 0.58″ respectively; see Table.1). This suggests that given the
exceptionally good seeing conditions, these objects exhibit even greater
dissimilarities than others when compared to the stacked image. However,
considering the minimal information loss (0.003% relative to the original
dataset), we find this result to be more than acceptable.

Figure 11: Three of the objects flagged as anomalies for which we do not have
a proper motivation. Some of them can still be low S/N objects. For each
source, the epoch flagged is the central one, while on the left and on the
right, we have the previous and the next epoch, respectively. (See Figs. 21
and 22 for the whole sample.)

Finally, we performed a test on all 989,725 individual observations, including
epochs rejected by the sigma-clipping criterion, which resulted in 1,715
flagged images. A large part of these are already included in the $\sim
85,000$ objects discarded by the sigma clipping, while the majority of the
remaining objects are the same ones selected by our approach after applying
the sigma clipping. It is clear, however, that our algorithm is not capable of
reproducing the sigma-clipping results, and therefore it is better to use it
after the preliminary cleaning procedure as a further refinement in order to
flag and eliminate most of the remaining outliers.

###  4.2 Small dataset

We then transitioned to utilizing the Label Set dataset. This smaller dataset
enabled us to visually explore data filtered by sigma clipping. It was derived
by De Cicco et al. (2021) and consists of all the objects in the field for
which a classification of the sources, in terms of star, galaxy, and AGN, is
available. The size of this dataset made it an ideal case for testing for any
possible issue with our new algorithm. Our method did not identify all the
anomalies in this dataset, comprising 2,675 sources and a total of 147,125
individual observations. Post sigma clipping, this number decreased to
140,920, with 6,205 epochs eliminated by the sigma clip.

By executing our algorithm on the complete Label Set dataset, we identified
289 anomalies, of which 202 were already flagged by the sigma clip. For the
anomalies not identified by the sigma clip, a significant number was found to
be tracks, in some cases overlapping the source.

Our focus then shifted to the 6,003 sources that were clipped but not
identified by our algorithm. Notably, these sources did not exhibit images
entirely distorted by artifacts that we know are present based on the
unfiltered dataset, suggesting both methods effectively detect and filter such
image-related issues.

However, there are some cases where something was overlooked, although we
expected our algorithm to easily spot the epoch as anomalous. We explore two
examples of the few cases where our algorithm did not identify the problem.

For source 0035, by manually inspecting the light curve, we identified the
following issue: In epochs 22 and 23, there is a large contamination coming
from a defect of the image on the left, while in epoch 42 there is a track. In
addition, there are epochs with a lower S/N (for instance epoch 32). The
algorithm spotted epoch 23, where the contamination is bigger, but completely
missed the two other occurrences. By looking at the histogram of the distances
(see Fig. 12) as computed by our algorithm, we noticed that all the
problematic epochs were in fact detected as outliers of the distribution but
are simply not different enough to meet our exclusion criteria.

Figure 12: Source 0035. First row: Histogram of the distances from the stacked
image. Second row: Stacked image. Last row: Epochs 22, 23, 32, and 42.

Source 0076 exhibits a track in epoch 51. In Fig. 13, we show two more epochs
with epochs 51 and 32 for comparison. From the histogram, it appears that the
track, although being anomalous, is considered more similar to the stacked
image; we spotted a similar behavior in other sources. In all of these cases,
the source with the track is outside the distribution of standard epochs, but
the stacked image is still above it. This effect could be mitigated by
considering a ”weirdness” threshold. Regardless of the position of the stacked
image, the image was flagged as to be rejected, and clearly this, on one hand,
reduces the number of wrong epochs not flagged, but, on the other hand, it
increases the number of good epochs rejected. Although we added this
possibility to the code, we consider that performing the sigma clip before
implementing our algorithm provides the best results.

Figure 13: Source 0076. First row: Histogram of the distances. Second row:
Stacked image. Last row: Epochs 1, 32, 47, and 51.

##  5 Discussion and conclusion

In this study, we propose a novel method for detecting outliers in
astronomical time series, using VST monitoring campaigns of the COSMOS field
as a test bed. Astronomical time series, crucial for understanding celestial
phenomena, present challenges such as noise, gaps, contaminants, and
artifacts. Our approach, which falls under so-called transfer learning
methods, combines a deep neural network, specifically an EfficientNet network
trained on ImageNet, with a k-NN algorithm. The EfficientNet network serves as
a feature extractor, and the k-NN algorithm measures the distance in the
feature space to identify potential outliers. Anomalies are flagged based on
the distance from the first neighbor and comparison with the same distance
obtained using the stacked image, representing the normal state of the scene
that has, by construction, a higher S/N and hence represents an outlier by
definition.

Our experiments were separated into two cases extracted from the same dataset,
one involving a larger dataset (Unlabel Set) and one involving a smaller
dataset (Label Set), in order to verify the effectiveness of our method in
detecting various anomalies, including a low S/N ratio, contamination from
neighboring sources, data reduction errors, bright objects affecting
background estimation, and tracks. The performance of the algorithm was
evaluated in comparison to a classic sigma-clipping procedure commonly used in
astronomical data analysis.

In the case of the larger dataset (Unlabel Set), our method identified
anomalies that were missed by the sigma-clipping algorithm, including
instances of a low S/N ratio, contamination, and tracks. Additionally, it
successfully flagged epochs affected by bright saturated objects, data
reduction errors, and off-center tracks (regarding the latter aspects, the
method seems to be very efficient). Although there were some epochs that were
flagged by the method for which we do not have a proper explanation, even when
considering all of them as false positives, they are negligible in number with
respect to the whole dataset.

In the case of the smaller dataset (Label Set), we did the opposite by looking
at the anomalies identified by the sigma clip that were not flagged by our
algorithm. In most cases, the distribution of the distances showed that the
epochs that should clearly be removed visually, although deviating from the
distribution, are below the threshold identified by the stacked image. This
suggests that both of the analyses need to be performed, that is, the sigma
clipping first and then a refinement with the method proposed in this work.

One key aspect of our method is that unlike traditional approaches, such as
sigma clipping, that focus on identifying outliers in a limited parameter
space (e.g., extreme magnitude variations in the light curve), our method
examines the overall structure of the images across the time series,
considering the relationships and patterns in the feature space. This is
particularly evident in the light curves presented in the appendix (see Fig.
20), where anomalies flagged by our method go beyond the typical peak-
detection scenarios. Peaks alone do not capture the diverse range of anomalies
present in astronomical time series, and our method provides a more
comprehensive approach to outlier detection. In substance, a simple adjustment
of the sigma-clipping parameters would not suffice to capture the nuanced
anomalies identified by our method. This highlights the complementary nature
of our approach, offering a distinct and valuable perspective in the realm of
artifact detection in astronomical time series.

The presence of defected epochs in astronomical time series can significantly
impact statistical indicators commonly derived from light curves. One notable
example is the pair slope, a metric used for instance in De Cicco et al.
(2021). The pair slope measures the trend of variability between consecutive
data points in a light curve, providing insights into the underlying
astrophysical processes. However, the accuracy of such indicators is
compromised when defected epochs are present, and these epochs may hamper the
ability to detect and characterize a variety of important astrophysical
events, such as tidal disruption events, supernovae, AGN, and blazars. Our
method, by effectively identifying and flagging these anomalies, contributes
to the preservation of the integrity of statistical measures derived from
light curves, ensuring that the results are more reliable and reflective of
the true astrophysical phenomena.

As a further development in the future, we plan to explore different
pretrained networks, such as vision transformers Dosovitskiy et al. 2020.
Additionally, a possible improvement could be found by exploring adaptive
thresholding methods based on the distribution of the distance meant to work
with the threshold determined by the stacked image in order to make the
algorithm more robust to variations in data characteristics.

In addition to the effectiveness of the proposed method in detecting outliers
in astronomical time series, it is noteworthy to discuss the computational
efficiency and parallelization potential. The nature of the algorithm lends
itself to an “embarrassingly” parallel paradigm. Each time series is entirely
independent from the others, allowing for parallel execution on separate
computational nodes. This parallelization capability makes the method well
suited for deployment on high-performance computing clusters or distributed
computing environments.

During our testing phase, the algorithm demonstrated impressive efficiency
even on a machine without a GPU – we utilized an Intel® Core™ i9-10980XE CPU @
3.00GHz.111https://ark.intel.com/content/www/us/en/ark/products/198017/intel-
core-i9-10980xe-extreme-edition-processor-24-75m-cache-3-00-ghz.html On
average, processing each time series consisting of a maximum of 54 epochs took
less than 1 second. It is important to highlight that this timing was achieved
on a CPU, and the use of a GPU is expected to further accelerate the
processing time.

The scalability and adaptability of the algorithm to parallel processing
environments make it a promising solution to efficiently handle large-scale
astronomical datasets, such as the ones that will be obtained from the LSST.
Future implementations may leverage GPU resources to achieve even faster
processing times, which is particularly beneficial when dealing with extensive
datasets generated by modern astronomical surveys.

## Acknowledgements

DD acknowledges PON R&I 2021, CUP E65F21002880003. DD and MP also acknowledge
the financial contribution from PRIN-MIUR 2022 and from the Timedomes grant
within the “INAF 2023 Finanziamento della Ricerca Fondamentale”. MB and SC
acknoweldge the ASI-INAF TI agreement, 2018-23-HH.0 ”Attività scientifica per
la missione Euclid - fase D”

Topcat (Taylor, 2005) and STILTS (Taylor, 2006) have been used for this work.
Some of the resources from Stutz (2022) has been used for this work. The k-NN
used in this work is part of the Scikit package (Pedregosa et al., 2011).

## References

  * Aigrain & Foreman-Mackey (2023) Aigrain, S. & Foreman-Mackey, D. 2023, ARA&A, 61, 329
  * Awang Iskandar et al. (2020) Awang Iskandar, D. N. F., Zijlstra, A. A., McDonald, I., et al. 2020, Galaxies, 8, 88
  * Botticella et al. (2017) Botticella, M. T., Cappellaro, E., Greggio, L., et al. 2017, A&A, 598, A50
  * Capaccioli & Schipani (2011) Capaccioli, M. & Schipani, P. 2011, The Messenger, 146, 2
  * Cappellaro et al. (2015) Cappellaro, E., Botticella, M. T., Pignata, G., et al. 2015, A&A, 584, A62
  * Cover & Hart (1967) Cover, T. & Hart, P. 1967, IEEE Transactions on Information Theory, 13, 21
  * De Cicco et al. (2021) De Cicco, D., Bauer, F. E., Paolillo, M., et al. 2021, A&A, 645, A103
  * De Cicco et al. (2022) De Cicco, D., Bauer, F. E., Paolillo, M., et al. 2022, A&A, 664, A117
  * De Cicco et al. (2015) De Cicco, D., Paolillo, M., Covone, G., et al. 2015, A&A, 574, A112
  * De Cicco et al. (2019) De Cicco, D., Paolillo, M., Falocco, S., et al. 2019, A&A, 627, A33
  * Deng et al. (2009) Deng, J., Dong, W., Socher, R., et al. 2009, in 2009 IEEE conference on computer vision and pattern recognition, Ieee, 248–255
  * Ding et al. (2019) Ding, Y., Sohn, J. H., Kawczynski, M. G., et al. 2019, Radiology, 290, 456
  * Doorenbos et al. (2022) Doorenbos, L., Torbaniuk, O., Cavuoti, S., et al. 2022, A&A, 666, A171
  * Dosovitskiy et al. (2020) Dosovitskiy, A., Beyer, L., Kolesnikov, A., et al. 2020, arXiv e-prints, arXiv:2010.11929
  * Esteva et al. (2017) Esteva, A., Kuprel, B., Novoa, R. A., et al. 2017, nature, 542, 115
  * Falocco et al. (2015) Falocco, S., Paolillo, M., Covone, G., et al. 2015, A&A, 579, A115
  * Fu et al. (2018) Fu, L., Liu, D., Radovich, M., et al. 2018, MNRAS, 479, 3858
  * Grado et al. (2012) Grado, A., Capaccioli, M., Limatola, L., & Getman, F. 2012, Memorie della Societa Astronomica Italiana Supplementi, 19, 362
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111
  * Kim et al. (2021) Kim, Y. J., Bae, J. P., Chung, J.-W., et al. 2021, Scientific Reports, 11, 3605
  * LeCun et al. (2015) LeCun, Y., Bengio, Y., & Hinton, G. 2015, nature, 521, 436
  * Lecun et al. (1990) Lecun, Y., Boser, B., Denker, J. S., et al. 1990, in Advances in Neural Information Processing Systems 2, ed. D. S. Touretzky (Morgan Kaufmann), 396–404
  * Li et al. (2018) Li, B., Yu, C., Hu, X., et al. 2018, in Algorithms and Architectures for Parallel Processing, ed. J. Vaidya & J. Li (Cham: Springer International Publishing), 284–298
  * Liu et al. (2020) Liu, D., Deng, W., Fan, Z., et al. 2020, MNRAS, 493, 3825
  * Malz et al. (2019) Malz, A. I., Hložek, R., Allam, T., et al. 2019, The Astronomical Journal, 158, 171
  * Martinazzo et al. (2021) Martinazzo, A., Espadoto, M., & Hirata, N. S. 2021, in 2020 25th International Conference on Pattern Recognition (ICPR), IEEE, 4176–4182
  * Menegola et al. (2017) Menegola, A., Fornaciali, M., Pires, R., et al. 2017, in 2017 IEEE 14th international symposium on biomedical imaging (ISBI 2017), IEEE, 297–300
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Poulain et al. (2020) Poulain, M., Paolillo, M., De Cicco, D., et al. 2020, A&A, 634, A50
  * Prima & Bouhorma (2020) Prima, B. & Bouhorma, M. 2020, ISPRS - International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 4443, 343
  * Scargle (1997) Scargle, J. D. 1997, in Astronomical Time Series, ed. D. Maoz, A. Sternberg, & E. M. Leibowitz (Dordrecht: Springer Netherlands), 1–12
  * Stutz (2022) Stutz, D. 2022, Collection of LaTeX resources and examples, https://github.com/davidstutz/latex-resources, accessed on 01.01.2024
  * Tan & Le (2019) Tan, M. & Le, Q. V. 2019, arXiv e-prints, arXiv:1905.11946
  * Taylor (2005) Taylor, M. B. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert, 29
  * Taylor (2006) Taylor, M. B. 2006, in Astronomical Society of the Pacific Conference Series, Vol. 351, Astronomical Data Analysis Software and Systems XV, ed. C. Gabriel, C. Arviset, D. Ponz, & S. Enrique, 666
  * Zou & Zhong (2018) Zou, M. & Zhong, Y. 2018, Sensing and Imaging, 19, 6

##  Appendix A Complete figures

Figure 14: Full set of sources flagged as problematic in epoch 32.

Figure 15: Objects with track.

Figure 16: continued

Figure 17: continued

Figure 18: continued

Figure 19: Problematic objects. In most of the cases, if a track appears, it
is closer to the object with respect to the aperture photometry, so it
slightly affects the flux measurement.

Figure 20: Light curves for the 52 problematic objects shown in Fig. 19. The
flagged epoch is marked in red. The horizontal lines represent the threshold
of the sigma clip.

Figure 21: Objects flagged as anomalies for which we do not have a proper
motivation. Some of the anomalies can still be low S/N objects. For each epoch
flagged, we added the one before and the one after for comparison. In two
cases (second line, left column and fourth line, right column), epoch 32 is
selected as the “after” epoch.

Figure 22: Objects flagged as anomalies for which we do not have a proper
motivation. Some of the anomalies can still be low S/N objects. For each epoch
flagged, we added the one before and the one after for comparison. In the last
two lines, six epochs appear since the method selected more than one epoch
(epochs 31 and 35; 34 was removed by the sigma clipping), and this causes a
series (from epoch 30 to 36). In the last line, there are only two epochs
since the selected one is the last.