# Measuring densities of cold atomic clouds smaller than the resolution limit.

A. Litvinov1, P. Bataille1, E. Maréchal1,2, P. Pedri1, O. Gorceix1, M. Robert-
de-Saint-Vincent1, B. Laburthe-Tolra1 1Laboratoire de Physique des Lasers,
Université Sorbonne Paris Nord F-93430 Villetaneuse, France and LPL CNRS, UMR
7538, F-93430 Villetaneuse, France
2Laboratoire de Physique et d’Etude des Matériaux (LPEM), CNRS, Université
Paris Sciences et Lettres

###### Abstract

We propose and demonstrate an experimental method to measure by absorption
imaging the size and local column density of a cloud of atoms, even when its
smallest dimension is smaller than the resolution of the imaging system. To do
this, we take advantage of the fact that, for a given total number of atoms, a
smaller and denser cloud scatters less photons when the gas is optically
thick. The method relies on making an ansatz on the cloud shape along the
unresolved dimension(s), and on providing an additional information such as
the total number of atoms. We demonstrate the method on in-situ absorption
images of elongated 87Sr Fermi gases. We find significant non-linear
corrections to the estimated size and local density of the cloud compared to a
standard analysis. This allows us to recover an un-distorted longitudinal
density profile, and to measure transverse sizes as small as one fourth of our
imaging resolution. The ultimate limit of our method is the wavelength that is
used for imaging.

The cold atoms community typically uses absorption imaging to explore the
physics of dense clouds of atoms. In many cases, especially in the case of in
situ observations to investigate $e.g.$ the equation of state of quantum
fluids Shin2008 ; ho2009 ; Nascimbene2010 ; Yefsah2011  or 1D gases
cazalilla2011 , these objects are extremely small and sometimes even below the
resolution limit of the imaging system. More generally, quantum gases reveal
interesting phenomena in small, local features: vortices vortex  whose size is
set by condensates’ healing length and typically cannot be resolved in situ,
density fluctuations fluctuations , Wannier functions in optical lattices,
etc. Resolving such structures is a difficult but rewarding problem that
prompted important technical developments, using for example high-resolution
objectives in the quantum gas microscope approach gross2017  (which are
nevertheless still limited to the diffraction limit), the newly demonstrated
quantum gas magnifier Asteria2021 , super-resolution imaging Subhankar2019 ;
Donald2019 , scanning probes using electrons wurtz2009  or ions Veit2021 .

Here we focus on the in-situ imaging of small and dense objects using standard
absorption imaging. In practice, trapped gases can be both extremely absorbing
and smaller than the imaging resolution, which constitutes a severe difficulty
for image analysis due to $e.g.$ total light absorption Reinaudi2007  or
diffraction of the optical fields used for imaging. We nevertheless show that
it is still possible to accurately measure the size and local density of an
object, even when this object is smaller than the imaging resolution.

The decrease in light intensity induced by propagation through an atomic
sample can be linked to the column density thanks to the Beer-Lambert law.
However, for objects smaller than the imaging resolution and when light
absorption is strong, the information on the column density is partially lost.
Indeed, the Beer-Lambert law is non linear and cannot be averaged over the
imaging resolution. Crucially, the average number of absorbed photons per atom
depends on the cloud size since atoms in a smaller and therefore optically
thicker medium are exposed to a reduced average light intensity. Our main idea
is to take advantage of this to reconstruct information on size and local
density at scales below the imaging resolution.

Figure 1: Recovering density profiles from distorted absorption images of
tightly confined gases. a) Geometry. b) Absorption image of a gas tightly
confined and expanding in one dimension. c) Time of flight absorption image
when the gas is released in all three directions. d) Cross section along the
short axis of the column density shown in b), at peak value. The raw data
(column density deduced from the pixelated optical depth) shows diffraction
fringes locally leading to non-physical negative optical depth (black
squares). Our method does not use the optical depth, but rather estimates
density by counting missing photons and assuming a Gaussian shape to the cloud
in its shortest dimension. The result of our method is shown by the red solid
lines, delimiting a confidence interval (red hashes) resulting from
experimental fluctuations. The raw data underestimates the atom number by a
factor 1.5, and the peak column density by a factor of roughly 4 in this case,
while our method allows to infer an object size of only 0.25 px (pixel size).
e) Linear density along the long axis (Oy), $n_{1D}(y)$ deduced either by
integrating the pixelated optical depth along the short axis (Ox) (black
dashes), or by our method (red solid lines, delimiting the confidence
interval).

In the first section of this paper, we explain our image processing, which
estimates the size of atomic clouds by using the non linearity of the Beer
Lambert law and an additional information. This latter, we here take to be an
independent measurement of the total number of atoms. In the second section,
we use this method to analyze images of ultra-cold strontium clouds at
temperature ranging from $0.14\,T_{F}$ to $1.6\,T_{F}$ (where $T_{F}$ is the
Fermi temperature). We focus on cases where the ultra-cold gases have expanded
while being strongly confined in two dimensions, which results in a very
elongated cloud with a very small transverse dimension due to confinement, and
a longitudinal profile that reproduces the initial momentum distribution. This
very elongated cloud, and in particular the measurement of its transverse
size, is our subject of study. In practice, using the averaged Beer Lambert
law results in a loss of information on this transverse size (which is smaller
than the resolution limit of our imaging system), and a strongly distorted
shape along the longitudinal axis. We show that our method enables us to
deduce the transverse size, which in practice was down to 0.25 times the
imaging resolution at the lowest temperature. This transverse size is in good
agreement with theoretical predictions; the longitudinal density profile
deduced from our method is also in agreement with the expected Fermi
distribution.

##  I Beer-Lambert Law for objects below imaging resolution

We first describe our approach to analyze absorption images of atomic clouds
smaller than the imaging system resolution. We will focus on a situation where
the profile of the column density is smaller than the resolution limit in one
direction only. The approach can be generalized to structures that are smaller
than the imaging resolution in the two dimensions perpendicular to the
propagation axis of the imaging laser beam.

Under standard conditions of imaging, the propagation along $Oz$ of resonant
light at intensities $I(x,y,z)$ below the saturation intensity, through an
atomic cloud of density $n(x,y,z)$, follows the Beer-Lambert law

$\frac{dI(x,y,z)}{I(x,y,z)}=-n(x,y,z)\sigma_{0}dz$ (1)

where $\sigma_{0}=3\lambda^{2}/2\pi$ is the resonant light absorption cross-
section for a single atom. $\lambda$ is the photon wavelength. We will
consider the situation where the longitudinal size of the object $l$ along the
direction of the imaging laser beam is smaller than the Rayleigh length
associated with its transverse size $\sigma$, $i.e.$
$l<\pi\sigma^{2}/\lambda$. In this regime, diffraction of the imaging beam as
it propagates through the cloud can be neglected over the distance $l$, such
that the intensity variations along $x,y,z$ are independent. Then Eq.(1) can
be integrated to provide a measurement of the optical depth (OD):

$OD(x,y)=-\log\frac{I(x,y)}{I_{0}(x,y)}=\sigma_{0}\int_{-\infty}^{\infty}n(x,y,z)dz$
(2)

where $I_{0}(x,y)$ is the transmitted light intensity without atoms and
$I(x,y)$ the transmitted light intensity with atoms. $OD(x,y)/\sigma_{0}$ thus
provides the column density $\int_{\mathbb{R}}n(x,y,z)dz\equiv\tilde{n}(x,y)$.

Experimentally, information on light intensity is limited by the imaging
resolution, set by $e.g.$ the detector pixel size, aberrations, or the
diffraction limit. We thus consider the collected power over a size $a$, at
position $\left\\{x=ia,y=ja\right\\}$, characterized by the indices $(i,j)$

$P(i,j)=\iint_{D_{i,j}}I(x,y)dxdy$ (3)

where
$D_{i,j}=\left[\left\\{ia,(i+1)a\right\\}\times\left\\{ja,(j+1)a\right\\}\right]$
is the domain of integration. It is possible to extract density information
from $\frac{P(i,j)}{P_{0}(i,j)}$ (where $P_{0}(i,j)$ is the collected power
without atoms) and Eq. (2) when $I(x,y)$ varies slowly over the distance $a$,
or when absorption is negligible. However, in general, and in the present
situation,

$\log\left(\frac{P(i,j)}{P_{0}(i,j)}\right)\neq\frac{1}{a^{2}}\iint_{D_{i,j}}\log\left(\frac{I(x,y)}{I_{0}(x,y)}\right)dxdy$
(4)

such that only measuring $P(i,j)$ and $P_{0}(i,j)$ is insufficient to provide
information on either the peak or the averaged column density.

Nevertheless, for a given total atom number, the total absorbed light power
strongly depends on the extension of the imaged object in the imaged plane
when absorption is strong, irrespective of the imaging resolution. This is due
to a shadowing effect in which the first atoms met by the imaging light reduce
the light intensity for the subsequent atoms. This is strongest for small size
samples because of increased density. While this effect is well captured by
the Beer-Lambert law, Eq. (4) indicates deviations when information is derived
from images that are pixelated or blurred at the scale of the sample size.
Then, the number of absorbed photons still depends on the size of the cloud,
and, as we shall now see, it is sill possible to deduce this size and the
local density if one provides independent information, such as the total atom
number.

In our experiment, the column density varies rapidly along the short axis $Ox$
and slowly along $Oy$, see Fig. 1. We will now relate the total number of
absorbed photons integrated along a pixel line parallel to $Ox$ to the local
column density. We define:

$R_{ph}(j)=\sum_{i}\frac{P_{0}(i,j)-P(i,j)}{P_{0}(i,j)}$ (5)

For brevity, we introduce $\delta I(x,y)=I_{0}(x,y)-I(x,y)$. We then relate
the pixelated information to local intensities $I(x,y)$ using Eq. (3), and
then use Eq. (2).

$\displaystyle R_{ph}(j)$ $\displaystyle=$
$\displaystyle\sum_{i}\frac{1}{P_{0}(i,j)}\iint_{D_{i,j}}dxdy\,\delta I(x,y)$
$\displaystyle=$
$\displaystyle\sum_{i}\iint_{D_{i,j}}dxdy\frac{I_{0}(x,y)}{P_{0}(i,j)}$
$\displaystyle\times\Big{(}1-\exp(-\sigma_{0}\int_{\mathbb{R}}n(x,y,z)dz)\Big{)}$
(6)

The discrete indices $(i,j)$ and the size $a$ can be matched to physical
pixels of the camera, or to effective pixels after binning. In any case, Eq.
(6) is valid provided the incident intensity in the atomic plane, imaged to
form $I_{0}(x,y)$ on the camera, is homogeneous at the spatial scale of the
resolution limit of the imaging system.

To simplify Eq. (6), we note that
$\frac{I_{0}(x,y)}{P_{0}(i,j)}=\frac{1}{a^{2}}$, and that
$\sum_{i}\int_{ia}^{(i+1)a}dx[..]=\int_{\mathbb{R}}dx[..]$ as the integrand
$[..]$ is zero far away from the cloud. Furthermore, in our case, the atomic
density along $Oy$ varies slowly over a pixel length, such that we can also
replace $\int_{ja}^{(j+1)a}dyf(y)$ (where $f$ is the integrand over $y$ in Eq.
(6)) by $af(ja)$.

Our idea is to introduce an ansatz on the local column density profile along
$Ox$, which for simplicity is taken as Gaussian:

$\tilde{n}(x,y)=e^{-\frac{x^{2}}{2\sigma_{x}^{2}}}\tilde{n}(0,y),$ (7)

Then, in the integration along $x$ of Eq.(6) we apply a change of variable
$u=x/\sigma_{x}$. We obtain:

$\displaystyle R_{ph}(j)$ $\displaystyle=$
$\displaystyle\frac{\sigma_{x}}{a}\int_{\mathbb{R}}\Big{(}1-\exp\Big{(}-\sigma_{0}\tilde{n}(0,aj)e^{-\frac{u^{2}}{2})}\Big{)}du$
$\displaystyle\equiv$
$\displaystyle\frac{\sigma_{x}}{a}F\big{(}\sigma_{0}\tilde{n}(0,aj)\big{)},$
(8)

where we defined the transfer function $F$.

The total number of atoms $N_{at}$ can be computed according to the definition
of the column density and to Eq.(7), and related to $R_{ph}$ using Eq.(8):

$N_{at}=\frac{\sqrt{2\pi}a\sigma_{x}}{\sigma_{0}}\sum_{j}F^{-1}\big{(}\frac{a}{\sigma_{x}}R_{ph}(j)\big{)}$
(9)

To obtain these equations, four assumptions have been made. First, the density
along the short axis $Ox$ is assumed to be a Gaussian of width $\sigma_{x}$.
Second, the column density $\tilde{n}(x=0,y)$ along the long axis $Oy$ is
quasi-uniform along the width of a pixel. Third, the illumination is
homogeneous at the scale of the resolution limit. Finally we also assumed that
diffraction of the probe light by the atomic sample can be neglected over the
depth of the sample along $Oz$, so that the Beer-Lambert law can be applied
locally.

Both Eq.(8) and Eq.(9) give the possibility to reconstruct the actual density
profiles even when the Beer-Lambert law cannot be directly used on the raw
data. Indeed, $F$ is non linear so that a given density profile
$\tilde{n}(0,aj)$ results in distorted profiles $R_{ph}(j)$ that
parametrically depend on $\sigma_{x}$. Therefore, we can in principle retrieve
the transverse size $\sigma_{x}$ from Eq.(8) when the longitudinal density
profile is known a priori. This is the case for example in our experiment at
high temperature, where the longitudinal cloud shape along $Oy$ should obey
Boltzmann statistics and be Gaussian. Alternatively Eq.(9) simply relates the
free unknown parameter $\sigma_{x}$ to the total number of atoms, which is
useful provided this number can be measured independently.

In what follows, we obtain $N_{at}$ by an independent measurement on expanded
clouds for which the gas is larger than the resolution in 3D, so that the
usual analysis holds and
$N^{3D}_{at}=\frac{a^{2}}{\sigma_{0}}\sum_{i,j}OD^{3D}(i,j)$ ($OD^{3D}$ is the
optical depth for the images with 3D expansion). We point out that the
estimate of $\sigma_{x}$ that we thus perform is independent from
uncertainties in $\sigma_{0}$ (uncertainties that can arise for example for
atoms with an hyperfine structure due to optical pumping effects during the
imaging pulse). Indeed, setting $N_{at}=N^{3D}_{at}$, from the explicit
expressions of $N_{at}$ and $N^{3D}_{at}$ if follows immediately that
$\sigma_{x}$ is independent of $\sigma_{0}$.

##  II EXPERIMENTAL RESULTS

In our experiment, whose setup is detailed in bataille2019 , the 87Sr cold
atomic sample is obtained after evaporation to reach regimes from T$\simeq$ 2
$T_{F}$ down to T$\simeq 0.15\,T_{F}$. The optical dipole trap is made of an
horizontal anisotropic laser beam crossed by a second laser beam at 30 degrees
from vertical, see Fig.1.a. To produce a tightly confined gas in one
dimension, we switch the horizontal beam off. This allows for an expansion of
the gas channeled by the second laser beam, for times ranging from 0.1 ms to
20 ms. We then take an absorption image, see Fig.1.b, using a pulse of
circularly polarized resonant light whose intensity is about 50 times below
the saturation intensity, and a non-magnifying imaging system that uses a
telescope configuration to approximately conjugate the atomic plane to a CCD
chip with 6.5 $\mu$m-wide pixels. The telescope is made of two achromatic
lenses of local length $f=150$ mm, and diameter 50 mm. The opening diameter is
set to 30 mm to reduce geometrical aberrations, and the diffraction-limited
resolution corresponds to 2.8 $\mu$m. However, in practice, the resolution of
our imaging system is limited by the pixel size of the CCD chip. Furthermore,
as a function of time, the cloud is positioned at different depths $z$ along
the imaging axis, due to the orientation of the guiding beam (see Fig 1.a).
This results in images that can be slightly out-of-focus. To independently
calibrate the number of atoms, we release the gas by switching off both
trapping lasers and take absorption images after a three dimensional free
flight, see Fig.1.c.

As can be seen in Fig.1.b, images of tightly confined gases are affected by
pixelation and diffraction fringes perpendicular to the elongated direction.
Fig.1.d shows a cut through those fringes, which are visible because the
object is slightly out of focus. The only diffraction signature is along the
transverse axis, and arises because of fast variations of the atomic density
along this axis. Thus the light is redistributed along this direction only.
Therefore, we can measure $R_{ph}(j)$ on camera pixel lines without binning,
even when the pixel size $a$ is smaller than the actual optical resolution on
out-of-focus clouds. Using Eq. (8) we deduce the transverse size and the local
density profile of the gas, also shown in Fig.1.d.

Our approach thus recovers the local peak density from which we also infer the
actual density profile along the elongated axis. As shown in Fig.1.e the
linear density derived by integrating the pixelated optical depth has a
distorted shape compared to the linear density, as recovered by our method. As
far as the low density part is concerned the reconstructed profile matches
that deduced from integrating the pixelated optical depth, as expected for low
absorption. In contrast, the high density region, corresponding to larger
absorption, is underestimated when using the pixelated optical depth.

###  II.1 Analysis of the sub-resolution transverse size

Figure 2: Determination of $\sigma_{x}$ based on the independently measured
total atom number $N_{at}$, for two images (at 0.14 $T_{F}$ and 0.56 $T_{F}$)
whose peak optical depths significantly differ. The black line is the inverse
of Eq.(9) for a specific $R_{ph}$, i.e. a single image. The atom number
confidence intervals (delimited by the dotted lines) are due to atom number
fluctuations from shot to shot, and are the dominant cause for the
$\sigma_{x}$ confidence intervals (delimited by the dashed lines).
Qualitatively, a low optical depth results in a large slope, and thus large
$\sigma_{x}$ uncertainty. Figure 3: Time dependence of $\sigma_{x}$ from 2 ms
to 20 ms channeled expansion, at different degeneracies. Each point aggregates
typically 10 images, with error bars that show the standard deviations between
the individual measurements. Those standard deviations match the individual
picture confidence interval (see Fig. 2), and thus can be attributed to mostly
measurement uncertainty. The degeneracy is deduced from the reference images.

We now focus on our measurements of $\sigma_{x}$. Fig. 2 illustrates how the
atom number from reference images is used to recover $\sigma_{x}$. For a given
picture of the elongated cloud, we compute $R_{ph}(j)$, we plot the number of
atoms $N_{at}$ deduced from Eq.(9) as a function of an assumed value of
$\sigma_{x}$. From the reference measurement of $N_{at}^{3D}$ and its
uncertainty associated with shot to shot fluctuations, we deduce the
corresponding value and confidence interval of $\sigma_{x}$ that we will
present on the following figures.

On Fig. 3 we observe that $\sigma_{x}$ does not vary much with the guided time
of flight (its variations are similar to the approximately $10\%$ standard
deviations of the data). Such a small variation is expected because the gas
evolves while being confined in the guiding beam with a transverse confinement
frequency which we estimate to vary from 210 Hz at the initial position to 145
Hz after 20 ms of guided fall. Indeed, if we assume that the transverse
degrees of freedom follow adiabatically such an evolution, the decrease in
trapping frequency corresponds to an increase of only $20\%$ in the transverse
size. Our ability to measure $\sigma_{x}$ for a large range of expansion times
illustrates the robustness of the method as, from short to long time of
flight, the peak optical density varies by a factor 14 because of the
longitudinal expansion. Furthermore, the method provides accurate estimates of
$\sigma_{x}$ even for short expansion times (2 ms), which indicates that the
longitudinal profile varies slowly enough for the approximation made to derive
Eq.(8) to remain valid.

As shown in Fig. 4, $\sigma_{x}$ depends on the temperature and reaches 0.25
times the pixel size at lowest temperatures, $T\simeq 30\,$nK at
$0.14\,T_{F}$. This value, $i.e.$ $1.6\,\mu$m, is smaller than the typical
distortions introduced by diffraction of the imaging beam and out-of-focus
measurements. It is also smaller than the ultimate resolution of our imaging
system (set by the pixel size).

To compare our measurements to theoretical expectations, we use the reference
images of 3D expansions, and we measure the mean kinetic energy along the
longitudinal axis $E_{K}$. Based on the equi-partition of energy, and
neglecting inter-atomic interactions, we expect
$\frac{1}{2}m\omega^{2}\sigma_{x_{th}}^{2}=E_{K}$. In Fig. 4, we thus compare
the measured dependence of $\sigma_{x}$ as a function of the energy per atom
to the expected values of $\sigma_{x_{th}}$, using the independently measured
initial trap confinement frequency. We find an overall good agreement with
this rough theoretical model without free parameter, with discrepancies
smaller than 15 $\%$ throughout the curve. We point out that, based on the
resolution of the Boltzmann equations Menotti2002 , we expect interaction
effects on the measurement of $E_{K}$ during the expansion to be of the order
of 8 percent at the lowest temperature given our experimental parameters
Sonderhouse2020 , while the interaction effects on $\sigma_{x}$ are of the
order of 1 percent.

###  II.2 Longitudinal profile

Figure 4: Dependence of $\sigma_{x}$ as a function of kinetic energy along one
direction. The kinetic energy is deduced from the reference images. The black
line shows the expected size according to equipartition of energy. The error
bars show the standard deviations, aggregating data for all channeled
expansion times.

In order to further verify that our results match the ones that are expected
from the theoretical point of view, we now turn to the longitudinal profiles
that we obtain, and compare them to theory. In Fig. 5, we present both the
longitudinal density profile deduced from the pixelated optical depth, and the
1D density profile $n_{1D}$ deduced from our method.

Fig. 5.a shows a non degenerate gas at $1.6\,T_{F}$ after a channeled
expansion of 10 ms. For such a long expansion, the longitudinal profile should
reflect the momentum distribution along this axis before expansion. Both
profiles are fitted by Gaussian distributions since at this temperature the
momentum distribution is almost perfectly described by Boltzmann statistics.
The residuals to the fits show that the 1D profile deduced from the pixelated
optical depth (upper residuals) systematically distorts the shape of the
cloud, while our reconstruction method (lower residuals) naturally recovers
the expected Gaussian profile.

Figure 5: Density profiles along the long axis Oy. The black hashed surface is
the confidence interval on the measured $n_{1D}(y)$. The dashed line
corresponds to the density deduced from the pixelated optical depth. a) Non
degenerate gas at 1.6 $T_{F}$ after a 10 ms guided expansion. The red lines
are Gaussian fits for the linear density deduced from the pixelated optical
depth (upper normalized residual) and for $n_{1D}(y)$ (lower normalized
residual). The residuals show that the raw data differs significantly from the
expected Gaussian shape of a thermal gas. On the contrary, our method recovers
this expected Gaussian shape. b) Degenerate gas at 0.21 $T_{F}$ after a 18 ms
guided expansion. In blue, the expected density profile (with no free
parameter) matches the profile recovered by our method, as demonstrated by the
residuals shown below.

As shown in Fig. 1.a, this expansion is performed along an axis that makes an
angle with the imaging plane, which defines a parallax for our measurements.
In addition, during the expansion, the atoms feel a position-dependent
potential due to the divergence of the Gaussian laser beam. In practice this
results in an anti-confinement. The parallax and anti-confinement can both be
measured by comparing the size of the expanding cloud to its temperature, and
measuring the trajectory of the center of mass. For the analysis shown in Fig.
5, the parallax and anti-confinement simply result in a single multiplicative
correction parameter C to the cloud size for a given expansion time.

Fig. 5.b shows the linear density of a degenerate gas at $0.21\,T_{F}$. We
compare it to the expected density profile of a degenerate Fermi gas expanding
in 1D. There are no free parameters since this density profile only depends on
the number of atoms, the temperature and degeneracy, which are all measured
from the reference images, and $C$, measured from Fig. 5.a. Obviously, the
linear density profile from pixelated data is in severe disagreement with the
theoretical model, unlike the reconstructed profile. The agreement with this
latter profile is good throughout the curve as can be seen from the residuals.
However, for the lowest momenta, we also observe a strong noise that affects
the reconstruction. These irregularities are due to the strong non-linearity
of the transfer function $F$ at high density. This non-linearity may also bias
an initially symmetric noise. This constitutes a limitation to our
reconstruction scheme when the source of the noise is purely from technical
origin; however, it also signals that our method could unveil density
fluctuations that would otherwise be washed out by the limited optical
resolution.

###  II.3 Discussion and conclusion

We introduced a novel method to analyze absorption images and to measure
density profiles at scales below the imaging resolution limit that may occur
due to $e.g.$ pixelation, distortion of the images caused by aberrations, and
the diffraction by the numerical aperture. These limitations are well known in
the community and make the analysis of very small objects particularly
difficult. Here we go beyond the approach performed in bouchoule2006 ;
armijo2010  that corrects the absorption by a simple numerical factor. Since
the deviation to the Beer-Lambert law is non-linear, we find that it is
important to follow such a procedure as presented here, in order to recover
the exact local density and an un-distorted density profile, even along the
direction in which the cloud is larger than the imaging resolution.

In our experiment, we characterize objects as small as 0.25 times our pixel
size, despite diffraction fringes due to out of focus imaging as large as two
pixels - these objects are thus significantly below our resolution limit.
Comparisons with theoretical models of the cloud shape with no free parameter
confirm the validity of these results. We have further checked our method by
testing its sensitivity to the ansatz that is chosen for the transverse
profile. We found that using a Fermi distribution in the transverse profile
barely modifies our results. This good agreement relies on a certain
similarity between the actual density profile and the chosen ansatz. However,
our method would not be able to reveal unexpected features in the shape along
the non-resolved direction, since different density shapes could lead to the
same number of scattered photons.

We point out that although our method allows retrieving features below the
imaging resolution, its fundamental limit is $\simeq\lambda$, as $\lambda$
sets the scattering cross section of light by the atoms. As a consequence, an
object whose size would be below $\lambda$ would affect light propagation
identically irrespective of its actual size. In addition to this limit on the
transverse size $\sigma$, there is also a limitation associated with the depth
$l$ of the object along the axis of the imaging beam. For our equations to
hold, it is necessary that the diffraction of the light field can be neglected
over the distance $l$. Therefore is is important that
$l<\pi\sigma^{2}/\lambda$.

Finally, we point out that our method could easily be generalized to study
features smaller than the imaging resolution within clouds otherwise larger
than the imaging resolution. For example, using the appropriate ansatz, we
have in mind to measure the size of vortex cores within a 2D or a 3D
superfluid, that are typically below the imaging resolution, or other
hydrodynamic structures such as solitons Denschlag2000 .

###### Acknowledgements.

Acknowledgements: This research was funded by the Agence Nationale de la
Recherche (projects ANR-18-CE47-0004 and ANR-16-TERC-0015-01), the Conseil
Regional d’Ile-de-France, Institut FRancilien des Atomes Froids, DIM Nano’K
(projects METROSPIN and ACOST), DIM Sirteq (project SureSpin), and Labex
FIRST-TF (project CUSAS).

## References

  * (1) Y. Shin et al., Nature 451, 689 (2008), Y. Shin, Phys. Rev. A 77, 041603(R) (2008).
  * (2) T. Ho, Q. Zhou, Nature Physics 6, 131 (2009).
  * (3) S. Nascimbene et al., Nature 463, 1057 (2010)
  * (4) T. Yefsah et al., Phys. Rev. Lett. 107, 130401 (2011)
  * (5) M. A. Cazalilla et al., Rev. Mod. Phys. 83, 1405 (2011).
  * (6) M. R. Matthews et al., Phys. Rev. Lett. 83, 2498 (1999), B. P. Anderson et al., Phys. Rev. Lett. 85, 2857 (2000), K. W. Madison et al., Phys. Rev. Lett. 84, 806 (2000), J. R. Abo-Shaeer et al., Science 292, 476 (2001), E. Hodby et al., Phys. Rev. Lett. 88, 010405 (2001), C. Raman et al., Phys. Rev. Lett. 87, 210402 (2001), Zwierlein et al., Nature 435, 1047 (2005)
  * (7) R. N. Bisset, C. Ticknor, and P. B. Blakie, Phys. Rev. A 88, 063624 (2013)
  * (8) C. Gross and I. Bloch, Science, 357, 995 (2017)
  * (9) L. Asteria et al., arXiv:2104.10089 (2021)
  * (10) S. Subhankar et al., Phys. Rev. X 9, 021002 (2019)
  * (11) M. McDonald et al. Phys. Rev. X 9, 021001 (2019).
  * (12) P. Würtz et al., Phys. Rev. Lett. 103, 080404 (2009)
  * (13) C. Veit et al., Phys. Rev. X 11, 011036 (2021)
  * (14) G. Reinaudi et al., Opt. Lett. 32, 3143 (2007).
  * (15) P. Bataille et al., Phys. Rev. A 102, 013317 (2020)
  * (16) C. Menotti, P. Pedri, and S. Stringari, Phys. Rev. Lett. 89, 250402 (2002)
  * (17) L. Sonderhouse et al., Nature Physics 16, 1216 (2020)
  * (18) J. Esteve et al., Phys. Rev. Lett. 96, 130403 (2006)
  * (19) J. Armijo et al., Phys. Rev. Lett. 105, 230402 (2010)
  * (20) J. Denschlag et al., Science, 287, 97 (2000)


11institutetext: University of Delaware, Newark DE 19711, USA
11email<EMAIL_ADDRESS>
# Learning the Propagation of Worms in Wireless Sensor Networks

Yifan Wang  0000-0002-6312-4663 Siqi Wang  0000-0002-0089-1005 Guangmo Tong
0000-0003-3247-4019

###### Abstract

Wireless sensor networks (WSNs) are composed of spatially distributed sensors
and are considered vulnerable to attacks by worms and their variants. Due to
the distinct strategies of worms propagation, the dynamic behavior varies
depending on the different features of the sensors. Modeling the spread of
worms can help us understand the worm attack behaviors and analyze the
propagation procedure. In this paper, we design a communication model under
various worms. We aim to learn our proposed model to analytically derive the
dynamics of competitive worms propagation. We develop a new searching space
combined with complex neural network models. Furthermore, the experiment
results verified our analysis and demonstrated the performance of our proposed
learning algorithms.

###### Keywords:

Wireless sensor network Multi-worm propagation Complex neural network.

##  1 Introduction

Due to the recent advances in wireless and micro-mechanical sensor
technologies, wireless sensor networks (WSNs) have been widely used in various
real-time applications, for example, target tracking, biological detection,
health care monitoring, environmental sensing and etc. [15, 27, 24, 21, 16,
9]. WSNs are composed of the spatially distributed sensors that are low-cost,
small-sized, low-energy, and multi-functional [15, 33, 10]. Each sensor is
designed to sense and collect environmental information. It has the ability to
carry out a local computation to obtain the required data and will transfer
the required data to a higher-level host [3, 33]. Besides, the communication
paradigm of WSNs is based on the broadcasting such that each sensor can
communicate with its neighbors directly [30, 3].

The computer worms are the programs that can be disseminated across the
network by self-replicating to exploit the security or policy flaws [34, 26].
Although many advantages of the sensors are proposed, WSNs are vulnerable to
attacks by worms and variants. One reason is that sensor network protocols
primarily concentrate on power conservation, and each sensor is usually
assigned to a limited resource. Therefore it generally has a weak defense
mechanism [18, 28]. Besides, since the spread of worms usually relies on the
information involved in the infected sensors and WSNs often embed an efficient
propagation mechanism, the worms can achieve a rapid propagation speed [4, 15,
2]. Furthermore, some worms will upgrade and be more complicated in order to
propagate successfully, which will result in a widespread destruction [33].
Therefore, worms and variants become one of the main threats to network
security.

In order to control the spreading of worms, people seek to understand the
dynamic behaviors of worms. In practice, worms utilize different propagation
strategies. For example, scan-based techniques will identify the vulnerable
sensors based on the IP address [32, 37]. Topology-based techniques allow
worms to spread through the topology neighbors [33]. Many studies have modeled
these propagation procedures mathematically. Most researchers proposed
propagation models based on the classic epidemic models, for example, SIR and
SEIR models [14, 12]. In recent works, extensions to the epidemic models have
been proposed: Some studies considered the hybrid dynamic of worms and
antivirus propagation [1, 6, 11]. Some works investigated the virus behavior
containing the stability analysis of WSNs [36, 13, 25].

Considering a WSN composed of a number of sensors, the transmission rates for
different worms vary based on the distinct mechanism embedded in the sensors.
That means that for different worms and sensors, the dynamic behavior will
vary [7, 35]. Meanwhile, both scan-based and topology-based propagation
strategies select the targets based on the features related to the sensors
[33, 23]. Since there is a basic similarity between worm propagation through
wireless devices and the information diffusion process over the regular
network, from this standpoint, we propose a discrete communication model based
on the competitive diffusion models under multiple worms. We seek to learn our
proposed model to derive the dynamics of competitive worms propagation
analytically. Our contributions can be summarized as follows.

  * •

Propagation Model design. We design a communication model that simulates the
propagation of the competitive worm. Our model includes the special case with
the hybrid dynamics with virus and patch, which provides a new aspect to the
antivirus protection design.

  * •

Searching Space Analysis. We develop a new searching space combined with the
complex-valued neural network. To our best knowledge, this provides a starting
point for analyzing the worms propagation model. Besides, a back-propagation-
based learning algorithm for the complex-valued neural network has been
proposed.

  * •

Empirical Studies. We provide experimental validation of our model compared
with some baseline methods.

The paper is organized as follows. In the next section, we introduce some
related works of the different propagation models. Section 3 illustrates the
structure of the communication model and formulates the learning problem.
Section 4 describes the learning framework, including the searching space
design and the learning algorithm. Experiments are given in section 5.

##  2 Related Work

A large number of models have been proposed to analyze the dynamics of worm
propagation in WSNs based on epidemic models. Kephart et al. first proposed
the idea of formulating the worm’s propagation based on simple epidemic
models. They assumed each sensor or host stays in one of two states:
susceptible or infectious. And the infected sensor will not change its
infectious state [14]. In recent works, Abdel et al. proposed an extension SIR
model by considering antivirus dynamics that add the antidotal state to each
sensor [1]. And the authors present a unified study of the epidemic-antivirus
nonlinear dynamical system. The work in [2] designed a fractional-order
compartmental epidemic model by accounting for the global stability analysis.
Several studies added the quarantined constraints for each sensor to the
epidemic models [15, 19, 20].

For multiple worms propagation, authors in [30] proposed a model that supports
the sleep and work interleaving schedule policy for the sensor nodes. They
assigned different infect probabilities to each worm when the infected nodes
attempted to infect other nodes. However, they did not specify the
relationship among these worms, and it is hard to track the dynamic of
different worms. [22] developed a SEjIjR-V Epidemic Model with differential
exposure, infectivity, and communication Range. Besides, they generate the
epidemic thresholds for each sensor corresponding to different worms. In our
work, we formulate a worms propagation model based on the works in [31]. And
inspired by them, we develop an explicitly complex-valued neural network to
exactly simulate the worms’ propagation process.

##  3 Preliminaries

In this section, we present the communication model and problem settings.

###  3.1 Communication Model

We consider the wireless network as a directed graph $G=(V,E)$, where $V$ and
$E$ denote the sensors and the communication links, respectively. Let $|V|=N$
and $|E|=M$. For each node $v$, its in-neighbors is denoted by $N_{v}\subseteq
V$. It is assumed that $K$ different worms are spread over the network, and
each worm $k\in[K]$ propagates from a set of sensors $S^{k}\in V$. For each
node $v\in V$, there is a set of threshold
$\Theta_{v}=\\{\theta_{v}^{1},...,\theta_{v}^{K}\\}$ that represents its
vulnerabilities to different worms, where $\theta\in[0,\infty)$. For each edge
$(u,v)\in E$, there is a set of weights
$W_{(u,v)}=\\{w_{(u,v)}^{1},...,w_{(u,v)}^{K}\\}$ indicates its levels of
transmission with different worms, where $w\in[0,\infty)$. We assume that a
node can only be infected by one worm and suppose that the spread of worms
continuous step by step based on the following rules.

  * •

Initialization. For nodes $v\in S^{k}$, they are $k$-infected; Otherwise, the
nodes stay innocent.

  * •

Step $t$. There are two stages during each step. Let $S^{k}_{t-1}\subseteq V$
be the $k$-infected nodes at the end of step $t-1$.

    * –

Stage 1. For each innocent node $v$, it can be possibly infected by the worms
in set $P_{v}^{t}\subseteq[K]$, where

$P_{v}^{t}=\Big{\\{}k\Big{|}\sum_{u\in N_{v}\cap
S^{k}_{t-1}}w_{(u,v)}^{k}\geq\theta_{v}^{k}\Big{\\}}.$ (1)

    * –

Stage 2. Node $v$ will finally be infected by worm $k^{*}$ such that

$k^{*}=\operatorname*{arg\,max}_{k\in P_{v}^{t}}\sum_{u\in N_{v}\cap
S^{k}_{t-1}}w_{(u,v)}^{k}$ (2)

Once a node is infected, its status will not change during the propagation.
For tie breaking, each node is more vulnerable with the worm with a large
index.

  * •

Termination. The propagation process ends when no nodes can be infected. We
can observe that $t\leq N$.

###  3.2 Problem Settings

We begin by presenting some notations that are used in the rest of the paper.

Notations. For a complex number $a+ib$, $a$ is the real part and $b$ is the
imaginary part, where $a,b\in\mathbb{R}$. For convenience, we use $(a,b)$ to
simplify the notation. During the propagation process, for a step $t\in[N]$,
we denoted the status of each node $v\in V$ by a complex-valued vector
$I_{t}(v)\in\mathbb{C}^{K}$: If $v$ is $k$-infection, then the $k_{th}$
element in $I_{t}(v)$ is $1+i0$ and the rest elements are $0+i0$; If $v$ is a
innocent node, then all the elements in $I_{t}(v)$ are $0+i0$. The status over
all nodes, named all-infection status, is then formulated as a matrix
$\mathcal{I}_{t}\in\mathbb{C}^{K\times N}$ by concatenating $I_{t}$ for each
node, where $\mathcal{I}_{t,v}=I_{t}(v)$. Furthermore, to simplify, we assume
that $K=2^{P}$ where $p\in\mathbb{Z}^{+}$.

We are interested in learning propagation functions given the collected
propagation data.

######  Definition 1.

Given a network graph and a communication model, the propagation function
$F:\mathbb{C}^{K\times N}\to\mathbb{C}^{K\times N}$ maps from initial all-
infection status $\mathcal{I}_{0}$ to the final status $\mathcal{I}_{N}$.

######  Definition 2.

The propagation data set $D$ is a set of initial-final status pairs with size
$Q\in\mathbb{Z}^{+}$, where

$D=\\{(\mathcal{I}_{0},\mathcal{I}_{N})^{i}\\}_{i=1}^{Q}.$ (3)

Given a new initial-final status pair $(\mathcal{I}_{0},\mathcal{I}_{N})$, let
$\hat{\mathcal{I}}_{N}$ be the predicted final all-infection status, we use
$l$ to measuring the loss.

$l(\hat{\mathcal{I}}_{N},\mathcal{I}_{N})=\frac{1}{N}|\\{v:\hat{\mathcal{I}}_{N,v}\neq\mathcal{I}_{N,v},v\in
V\\}|$ (4)

The problem can be formally defined as follows.

######  Problem 1.

Given the graph structure $G$ and a set of propagation data $D$ draw from some
distribution $\mathcal{D}$, we wish to find a propagation function
$F:\mathbb{C}^{K\times N}\to\mathbb{C}^{K\times N}$, that minimizes the
generalization error $L$, where

$L(F)=\mathbb{E}_{(\mathcal{I}_{0},\mathcal{I}_{N})\sim\mathcal{D}}l(F(\mathcal{I}_{0}),\mathcal{I}_{N})$
(5)

##  4 A Learning Framework

In this section, for the communication model, we first design a searching
space that composed by the complex neural network. Then the learning algorithm
is proposed given such searching space.

###  4.1 Hypothesis Space Design

In this section, we construct the hypothesis space for problem 1 by simulating
the propagation process as a finite complex-valued Neural Network (CVNNs) from
local to global. Assume that for the complex-valued neuron $n$ at layer $l+1$,
it is linked with $m$ numbers of neurons in the previous layer, then the
output of neuron $n$ is described as follows:

$H^{l+1}_{n}=\sigma^{l+1}_{n}(\sum_{i=1}^{m}W^{l}_{in}H^{l}_{n}+T_{n}^{l})$
(6)

Here $W^{l}_{in}$ denotes the weight connecting the neuron $n$ and the neuron
$i$ from the previous layer, $\sigma^{l+1}_{n}$ denotes the activation
function on neuron $n$, $T_{n}^{l}$ is the bias term.

####  4.1.1 Local Propagation

For the local propagation, it simulates the diffusion process during every
step $t$. It starts by receiving the all-infection status
$\mathcal{I}_{t}\in\mathbb{C}^{K\times N}$ at the end of step $t$. Figure 1
illustrates an example of local propagation process.

Figure 1: Local propagation structure. Suppose that two worms are disseminated
in the network. The top half of the figure shows the whole propagation process
across a small WSN. The bottom half shows the structure of our proposed
complex neural network during step 1. Each stages are simulated by several
layers in the network. And three additional layers are added to control the
format of the output.

Stage 1. We aim to find out all the possible worms that can infect the each
innocent nodes. During this stage, each worm diffuses separately and follows a
linear threshold process. Therefore, we wish the links between the layers
represents the topology in the WSNs and the activation function can be simply
defined as the threshold functions. In addition, the communication model
indicates that once a node is affected, it will never change the status.
Therefore, for each node, we add one link towards it self at the next layer
with a relatively large value. The propagation can be simulated by

$\displaystyle\mathbf{H}^{2}=\sigma^{1}(\mathbf{W}^{1}\mathbf{H}^{1}+\mathbf{T}^{1}),$
(7)

where for $i,j\in[N]$ and $k,l\in[K]$, we have

$\displaystyle\mathbf{W}^{1}_{i\cdot K+k,j\cdot
K+l}\operatorname{\coloneqq}\begin{cases}(w_{(i,j)}^{k},0)&\text{if $(i,j)\in
E$ and $k=l$}\\\ 1000&\text{if $i=j$ and $k=l$}\\\
(0,0)&\text{otherwise}\end{cases}$

$\sigma^{1}_{j\cdot K+l}(x)\operatorname{\coloneqq}\begin{cases}x&\text{if
$\operatorname*{Re}(x)\geq\theta_{j}^{l}$}\\\
(0,0)&\text{Otherwise}\end{cases}.$ (8)

$T^{1}_{j\cdot K+l}(x)\operatorname{\coloneqq}(0,l)$ (9)

Stage 2. With the input $\mathbf{H}^{2}$ from phase 1, this stage determines
the final infect status of each node. That is for a node $i\in V$, we want to
find out the largest value among the sub-array

$\Big{[}\operatorname*{Re}(\mathbf{H}^{2}_{i+1}),\operatorname*{Re}(\mathbf{H}^{2}_{i+2}),...,\operatorname*{Re}(\mathbf{H}^{2}_{i+K})\Big{]}.$
(10)

Intuitively, finding the largest value among equation 10 can be implemented
through a hierarchy structure of pairwise comparisons. Starting from the sub-
array, each hierarchy compares every two neighboring elements and passing the
larger value to the next level, which leads to a $P$ hierarchy levels
structure.

Figure 2: Hierarchy structure of pairwise comparisons. Suppose that there are
four different worms: $a,b,c$ and $d$. There are two comparison process with
two layers. The nodes in the bottom layer represent the sub-array 10. All the
weights on the black solid links are $1+i0$ and the weights on the dash links
are $-1+i0$. During each comparison process in each two nodes, the first layer
implements the activation function in equation 11 such that: if the prior one
has the larger value of real part, the node returns the difference both in
real and imaginary part; Otherwise, it will return $(0,0)$. Therefore, by
combining the second layer, the output is exactly the larger value between two
nodes.

From bottom to top, for node $v\in[N]$, in each hierarchy level $p\in[P]$, the
size of input elements $I_{p}=\frac{K}{2^{p-1}}$ and the size of output
elements $O_{p}=\frac{K}{2^{p}}$. Each level can be decoded by two layers:
$l^{c_{1}}$ and $l^{c_{2}}$.

In $l^{c_{1}}$, for $i,j\in[I_{p}]$,

$\displaystyle\mathbf{W}^{c_{1}}_{i,j}\operatorname{\coloneqq}\begin{cases}(1,0)&\text{if
$i=j$}\\\ (0,-1)&\text{if $j=i-1$ for $j=2n+1$ where
$n\in\\{0,1,2,...,I_{p}-1\\}$}\\\ (0,0)&\text{otherwise}\end{cases}$

$\sigma^{l_{c1}}$ is defined as follows:

$\sigma^{l_{c_{1}}}(x)\operatorname{\coloneqq}\begin{cases}(\operatorname*{Im}(x),\operatorname*{Re}(x))&\text{if
$\operatorname*{Im}(x)\geq 0$}\\\ (0,0)&\text{otherwise}\end{cases}$ (11)

In $l^{c_{2}}$, for $i\in[I_{p}]$ and $j\in[O_{p}]$,

$\displaystyle\mathbf{W}^{l_{c_{2}}}_{i,j}\operatorname{\coloneqq}\begin{cases}(1,0)&\text{if
$i=2\cdot j$ or $i=2\cdot j-1$}\\\ (0,0)&\text{otherwise}\end{cases}$

$\sigma^{l_{c_{2}}}$ is defined as follows.

$\displaystyle\sigma^{l_{c_{2}}}(x)\operatorname{\coloneqq}x$

By repeatedly using these two layers $P$ times, the output of stage 2 is
$\mathbf{H}^{2P+2}\in\mathbb{C}^{N}$. For $i\in[N]$ and $k\in[K]$, each node
in $\mathbf{H}^{2P+2}$ holds the value as follows,

$\mathbf{H}^{2P+2}_{i}=\Big{(}\max_{k\in
K}\\{\operatorname*{Re}(\mathbf{H}^{2}_{i+k})\\},\operatorname*{arg\,max}_{k\in
K}\\{\operatorname*{Re}(\mathbf{H}^{2}_{i+k})\\}\Big{)},$ (12)

where the real part for each node $i$ indicates the largest incoming weights
among all the worms, and the imaginary part is the corresponding worm index.

Format Adjustment. In order to simulate the global worm propagation process,
we need to make the format of the output at the end of local propagation be
same as the input, that is
$\mathbf{H}^{2Z+3}=\mathcal{I}_{t+1}\in\mathbb{C}^{K\times N}$. Two layers are
added. Figure 3 shows an example of the structure of the format adjustment.
The two layers can be formulated as follows.

Figure 3: Structure of the format adjustment. For node $v\in V$, suppose that
there are four different worms: $a,b,c$ and $d$. All the weights on the black
solid links are $1+i0$ and the weights on the dash links are $-1+i0$. The
value on the nodes with shadow is $1+i0$; Otherwise, it is $0+i0$. In
addition, the threshold on each node $v^{k}$ where $k\in\\{a,b,c,d\\}$ is $k$.

$\displaystyle\mathbf{H}^{2P+2+i}=\sigma^{2Z+2+i}(\mathbf{W}^{2Z+2+(i-1)}\mathbf{H}^{2Z+2+(i-1)}),$
(13)

where $i=\\{1,2\\}$. The first layer can be seen as a transition layer. The
number of the output nodes will be recovered as the input, and the real number
of each node are set to $0$. For the imaginary part, its value will be equal
to $1$ if the node’s corresponding worm’s index is no large than the real part
in equation 12. Therefore, it is formulated as follows. For each $i,j\in[N]$,
$k\in[K]$, we have

$\displaystyle\mathbf{W}^{2P+3}_{i\cdot K,j\cdot
K+k}\operatorname{\coloneqq}\begin{cases}(1,0)&\text{if $i=j$}\\\
(0,0)&\text{otherwise}\end{cases}$

and

$\sigma_{j\cdot
K+k}^{2P+3}(x)\operatorname{\coloneqq}\begin{cases}(1,0)&\text{if
$\operatorname*{Im}(x)$ and $\operatorname*{Im}(x)\geq k$}\\\
(0,0)&\text{Otherwise}\end{cases}.$ (14)

The last layer are designed to extract the all-infection status at the end of
local propagation, and that is $\mathbf{H}^{2P+4}=\mathcal{I}_{t}$

where for each $i,j\in[N]$, $k,l\in[K]$, we have

$\displaystyle\mathbf{W}^{2P+4}_{i\cdot K+k,j\cdot
K+l}\operatorname{\coloneqq}\begin{cases}(1,0)&\text{if $i=j$ and $l\geq
k$}\\\ (-1,0)&\text{if $i=j$ and $l<k$}\\\ (0,0)&\text{otherwise}\end{cases}$

and

$\sigma_{j\cdot
K+l}^{2P+4}(x)\operatorname{\coloneqq}\begin{cases}(1,0)&\text{if
$\operatorname*{Re}(x)>l$}\\\ (0,0)&\text{Otherwise}\end{cases}.$ (15)

We can then obtain the following lemma.

######  Lemma 1.

The local propagation under transmission model can be simulated by a complex
neural network with $2P+4$ layers.

####  4.1.2 From Local to Global

The output of local propagation complex neural network can be viewed as the
input for the next time step. We can repeatedly using the local network to
exactly simulated the whole propagation process. We denoted such global
complex neural network by $M_{C}$. Therefore, we can obtain our hypothesis
space of problem 1.

######  Definition 3.

Taking the weights $W_{(u,v)}$ for $(u,v)\in E$ and thresholds $\Theta_{v}$
for $v\in V$ as the parameters, the hypothesis space $\mathcal{F}$ is the set
of all possible $M_{C}$.

###  4.2 Learning Algorithm

Given the hypothesis space $\mathcal{F}$, the propagation function $F$ based
on the propagation model is now implemented as complex neural network models.
Instead of directly learning the $F$, we seek to learn a complex neural
network model from $\mathcal{F}$ that can minimize the generalization error in
equation 5. Same as the neural network, learning in complex neural networks
refers to the process of modifying the parameters involved in the network such
that the learning objectives can be achieved. We choose the gradient-based
approach with the complex domain back-propagation [5]. During the training, We
utilize the complexPytorch package, which is a high-level API of PyTorch [17,
29].

##  5 Experiments

In this section, we evaluated our proposed method by experiments. We compare
our method with several baseline algorithm on synthetic and real worlds
networks.

Table 1: Main results. Each simulation is implemented with 600 training
samples, 400 testing samples and 4 worms.

|  | F1 score | Precision | Recall | Accuracy
---|---|---|---|---|---
ER graph | Proposed method | 0.469(6.E-03) | 0.419(6.E-03) | 0.629(6.E-03) | 0.727(3.E-03)
RegularNN | 0.408(4.E-03) | 0.329(4.E-03) | 0.590(6.E-03) | 0.695(3.E-03)
SVM | 0.184(0) | 0.192(0) | 0.201(0) | 0.498(1.E-03)
Random | 0.182(2.E-03) | 0.200(2.E-03) | 0.200(2.E-03) | 0.200(2.E-03)
Sensors network | Proposed method | 0.408(1.E-02) | 0.372(1.E-02) | 0.575(1.E-02) | 0.684(6.E-03)
RegularNN | 0.208(1.E-02) | 0.302(6.E-02) | 0.450(3.E-02) | 0.597(3.E-03)
SVM | 0.146(8.E-04) | 0.117(3.E-03) | 0.205(3.E-04) | 0.537(3.E-03)
Random | 0.176(7.E-03) | 0.200(4.E-03) | 0.200(8.E-03) | 0.199(7.E-03)

###  5.1 Settings

WSN data. We adopt two types of graph structures to evaluate our proposed
method: one synthetic graph and a real-world WSN. One synthetic graph is a
randomly generated Erdős–Rényi (ER) graph with 200 nodes and 0.2 edge
probability [8]. The real-world WSN is collected from 54 sensors deployed in
the Intel Berkeley Research lab 111
http://db.csail.mit.edu/labdata/labdata.html..

Communication model and sample generating The number of worms are select from
${2,4,8}$. For each node $v$, the threshold value of worm $k$,
$\theta_{v}^{k}$, is generated uniformly random from $[0,1]$. The weight on
each edge $(u,v)$ of worm $k$ is set as $w_{(u,v)}^{k}=\frac{1}{|N_{v}|+k}$.
For each graph, we generate a sample pool
$D=\\{(\mathcal{I}_{0},\mathcal{I}_{N})^{i}\\}_{i}^{1000}$ with $N/2$ seeds
and each seed is selected uniformly at random. Each selected seed node is
assigned to worm uniformly at random. Given $\mathcal{I}_{0}$,
$\mathcal{I}_{N})$ is generated follow the communication model.

Baseline method. We consider two baseline algorithms: RegularNN and Support
Vector Machine (SVM). RegularNN is a real-valued deep neural network model
with $4\cdot N$ layers. The activation functions are set as a combination of
linear, ReLU, and sigmoid. We aim to compare ComplexNN, which explicitly
simulates the diffusion process, with the regular designed neural network.
Since the proposed problem is a binary classification task, we implement the
SVM algorithm. Besides, the random method with a random prediction is also
considered.

Implementation The size of the training set is selected from
$\\{200,400,600\\}$ and the testing size is $400$. During each simulation, the
train and test samples are randomly selected from the sample pool. For each
method, we use F1 score, precision, recall, and accuracy to evaluate the
performance. Besides, the number of simulations for each setting is set to 5,
and we report the average performance together with the standard deviation.

Table 2: Impact of the different number of worms.  Each simulation is
implemented with 600 training samples, 400 testing samples, and 100 seeds.

Number of worms | F1 | Precision | Recall | Accuracy
---|---|---|---|---
2 | 0.717(3.E-02) | 0.657(1.E-02) | 0.948(3.E-02) | 0.803(2.E-02)
4 | 0.469(6.E-03) | 0.419(6.E-03) | 0.629(6.E-03) | 0.727(3.E-03)
8 | 0.277(4.E-03) | 0.206(5.E-03) | 0.313(5.E-03) | 0.656(3.E-03)

Table 3: Impace of the different number of seeds. Each simulation is
implemented with 600 training samples, 400 testing samples, and 4 worms.

Number of seeds | F1 score | Precision | Recall | Accuracy
---|---|---|---|---
50 | 0.332(4.E-03) | 0.277(4.E-03) | 0.579(7.E-03) | 0.676(4.E-03)
100 | 0.469(6.E-03) | 0.419(6.E-03) | 0.629(6.E-03) | 0.727(3.E-03)
150 | 0.549(5.E-03) | 0.524(4.E-03) | 0.626(4.E-03) | 0.733(3.E-03)

###  5.2 Observations

Table 1 shows the main results. Each cell includes the average value and its
corresponding standard deviation. We can observe that our method outperforms
other methods. Compared to the RegularNN, the proposed method are able to
improve the prediction errors. And this shows the advantage of using an
explicit model. Besides, it shows that the performance of each algorithm is
increased given a larger graph.

Impact of number of worms Table 2 gives the simulation results for our
proposed method under the different numbers of worms. We can observe that with
a lower number of worms, the performance is improved. One plausible reason is
that the number of worms determines the size of the number of input and output
features for our model. And it meets the nature of neural network training
properties.

Imbalanced class distribution The simulation results for our proposed method
under the different number of seeds are given in Table 2. We see that the
performance is increased when more seeds are provided. Since our model is a
binary classifier, a small number of seeds leads to an imbalance class
distribution. It is an inherent issue in predicting the node status in that
the majority of the nodes tend to be inactive.

##  6 Conclusion

In this paper, we formulate a propagation model for multiple worms and propose
a complex-valued neural network model. In particular, our designed complex-
valued neural network model simulates the propagation process exactly. Overall
our method outperforms classic supervised learning methods. However, it is not
sufficiently stable for different data within the different numbers of worms
and seeds. And it is possible that there exists an advanced representation
strategy for the neural network design that can achieve a better performance,
which is one direction for our future work. In addition, encouraged by the
possibility of using a complex-valued neural network to simulate the worms’
propagation models, we believe it would be interesting to apply our method to
the traditional epidemic model.

## References

  * [1] Abdel-Gawad, H.I., Baleanu, D., Abdel-Gawad, A.H.: Unification of the different fractional time derivatives: An application to the epidemic-antivirus dynamical system in computer networks. Chaos, Solitons & Fractals 142, 110416 (2021)
  * [2] Achar, S.J., Baishya, C., Kaabar, M.K.: Dynamics of the worm transmission in wireless sensor network in the framework of fractional derivatives. Mathematical Methods in the Applied Sciences 45(8), 4278–4294 (2022)
  * [3] Akyildiz, I.F., Su, W., Sankarasubramaniam, Y., Cayirci, E.: Wireless sensor networks: a survey. Computer networks 38(4), 393–422 (2002)
  * [4] Awasthi, S., Kumar, N., Srivastava, P.K.: A study of epidemic approach for worm propagation in wireless sensor network. In: Intelligent computing in engineering, pp. 315–326. Springer (2020)
  * [5] Bassey, J., Qian, L., Li, X.: A survey of complex-valued neural networks. arXiv preprint arXiv:2101.12249 (2021)
  * [6] Behal, K.S., Gakkhar, S., Srivastava, T.: Dynamics of virus-patch model with latent effect. International Journal of Computer Mathematics pp. 1–16 (2022)
  * [7] Chakrabarti, D., Leskovec, J., Faloutsos, C., Madden, S., Guestrin, C., Faloutsos, M.: Information survival threshold in sensor and p2p networks. In: IEEE INFOCOM 2007-26th IEEE International Conference on Computer Communications. pp. 1316–1324. IEEE (2007)
  * [8] Erdös, P.: Graph theory and probability. Canadian Journal of Mathematics 11, 34–38 (1959)
  * [9] Galluccio, L., Morabito, G.: Impact of worm propagation on vehicular sensor networks exploiting v2v communications. In: 2019 International Conference on Wireless and Mobile Computing, Networking and Communications (WiMob). pp. 1–6. IEEE (2019)
  * [10] Gulati, K., Boddu, R.S.K., Kapila, D., Bangare, S.L., Chandnani, N., Saravanan, G.: A review paper on wireless sensor network techniques in internet of things (iot). Materials Today: Proceedings (2021)
  * [11] Haghighi, M.S., Wen, S., Xiang, Y., Quinn, B., Zhou, W.: On the race of worms and patches: Modeling the spread of information in wireless sensor networks. IEEE Transactions on Information Forensics and Security 11(12), 2854–2865 (2016)
  * [12] Han, X., Tan, Q.: Dynamical behavior of computer virus on internet. Applied mathematics and computation 217(6), 2520–2526 (2010)
  * [13] Hu, Z., Wang, H., Liao, F., Ma, W.: Stability analysis of a computer virus model in latent period. Chaos, Solitons & Fractals 75, 20–28 (2015)
  * [14] Kephart, J.O., White, S.R., Chess, D.M.: Computers and epidemiology. IEEE spectrum 30(5), 20–26 (1993)
  * [15] Khanh, N.H.: Dynamics of a worm propagation model with quarantine in wireless sensor networks. Appl. Math. Inf. Sci 10(5), 1739–1746 (2016)
  * [16] Krishnamachari, B.: Networking wireless sensors. Cambridge University Press (2005)
  * [17] Matthès, M.W., Bromberg, Y., de Rosny, J., Popoff, S.M.: Learning and avoiding disorder in multimode fibers. Physical Review X 11(2), 021060 (2021)
  * [18] Mishra, B.K., Keshri, N.: Mathematical model on the transmission of worms in wireless sensor network. Applied Mathematical Modelling 37(6), 4103–4111 (2013)
  * [19] Mishra, B.K., Tyagi, I.: Defending against malicious threats in wireless sensor network: A mathematical model. International Journal of Information Technology and Computer Science 6(3), 12–19 (2014)
  * [20] Mishra, B.K., Srivastava, S.K., Mishra, B.K.: A quarantine model on the spreading behavior of worms in wireless sensor network. Transaction on IoT and Cloud Computing 2(1), 1–12 (2014)
  * [21] Nain, M., Goyal, N.: Localization techniques in underwater wireless sensor network. In: 2021 International Conference on Advance Computing and Innovative Technologies in Engineering (ICACITE). pp. 747–751. IEEE (2021)
  * [22] Nwokoye, C., Umeugoji, C., Umeh, I.: Evaluating degrees of differential infections on sensor networks’ features using the sejijr-v epidemic model. Egyptian Computer Science Journal 44(3) (2020)
  * [23] Nwokoye, C.H., Madhusudanan, V.: Epidemic models of malicious-code propagation and control in wireless sensor networks: An indepth review. Wireless Personal Communications pp. 1–30 (2022)
  * [24] Ojha, R.P., Srivastava, P.K., Sanyal, G., Gupta, N.: Improved model for the stability analysis of wireless sensor network against malware attacks. Wireless Personal Communications 116(3), 2525–2548 (2021)
  * [25] Qin, P.: Analysis of a model for computer virus transmission. Mathematical Problems in Engineering 2015 (2015)
  * [26] Rajesh, B., Reddy, Y.J., Reddy, B.D.K.: A survey paper on malicious computer worms. International Journal of Advanced Research in Computer Science and Technology 3(2), 161–167 (2015)
  * [27] Regin, R., Rajest, S.S., Singh, B.: Fault detection in wireless sensor network based on deep learning algorithms. EAI Transactions on Scalable Information Systems p. e8 (2021)
  * [28] Srivastava, A.P., Awasthi, S., Ojha, R.P., Srivastava, P.K., Katiyar, S.: Stability analysis of sidr model for worm propagation in wireless sensor network. Indian Journal of Science and Technology 9(31), 1–5 (2016)
  * [29] Trabelsi, C., Bilaniuk, O., Zhang, Y., Serdyuk, D., Subramanian, S., Santos, J.F., Mehri, S., Rostamzadeh, N., Bengio, Y., Pal, C.J.: Deep complex networks. In: International Conference on Learning Representations (2018), https://openreview.net/forum?id=H1T2hmZAb
  * [30] Wang, X., Li, Q., Li, Y.: Eisirs: a formal model to analyze the dynamics of worm propagation in wireless sensor networks. Journal of Combinatorial Optimization 20(1), 47–62 (2010)
  * [31] Wang, Y., Tong, G.: Learnability of competitive threshold models. arXiv preprint arXiv:2205.03750 (2022)
  * [32] Wang, Y., Wen, S., Cesare, S., Zhou, W., Xiang, Y.: The microcosmic model of worm propagation. The Computer Journal 54(10), 1700–1720 (2011)
  * [33] Wang, Y., Wen, S., Xiang, Y., Zhou, W.: Modeling the propagation of worms in networks: A survey. IEEE Communications Surveys & Tutorials 16(2), 942–960 (2013)
  * [34] Weaver, N., Paxson, V., Staniford, S., Cunningham, R.: A taxonomy of computer worms. In: Proceedings of the 2003 ACM workshop on Rapid Malcode. pp. 11–18 (2003)
  * [35] Yao, Y., Sheng, C., Fu, Q., Liu, H., Wang, D.: A propagation model with defensive measures for plc-pc worms in industrial networks. Applied Mathematical Modelling 69, 696–713 (2019)
  * [36] Zarin, R., Khaliq, H., Khan, A., Khan, D., Akgül, A., Humphries, U.W.: Deterministic and fractional modeling of a computer virus propagation. Results in Physics 33, 105130 (2022)
  * [37] Zou, C.C., Towsley, D., Gong, W., Cai, S.: Routing worm: A fast, selective attack worm based on ip address information. In: Workshop on Principles of Advanced and Distributed Simulation (PADS’05). pp. 199–206. IEEE (2005)