# Modeling Cloud Reflectance Fields using Conditional Generative Adversarial
Networks

Victor Schmidt , Mustafa Alghali*, Kris Sankaran & Yoshua Bengio
Mila, Université de Montréal
{schmidtv<EMAIL_ADDRESS><EMAIL_ADDRESS>
&Tianle Yuan
NASA Goddard Space Flight Center
University of Maryland Baltimore County
<EMAIL_ADDRESS>
Equal contribution

###### Abstract

We introduce a conditional Generative Adversarial Network (cGAN) approach to
generate cloud reflectance fields (CRFs) conditioned on large scale
meteorological variables such as sea surface temperature and relative
humidity. We show that our trained model can generate realistic CRFs from the
corresponding meteorological observations, which represents a step towards a
data-driven framework for stochastic cloud parameterization.

##  1 Introduction

Global Climate Models (GCMs) are one of the most important tools available to
understand and anticipate the consequences of climate change, including
changes in precipitation, increases in temperatures, and acceleration in
glacial melting [22]. One of the key physical principles these models rely on
is the Earth’s energy balance [14]: in short, the difference between how much
energy the earth receives and how much it emits. In this context, it is
paramount to model clouds accurately as they both reflect energy coming to the
Earth and the infrared radiations it radiates [18]. However, as physical
processes at play in cloud composition and evolution typically range from
$10^{-6}$ to $10^{6}$m, direct simulation of their behavior can consume up to
20% of a GCM’s computations - depending on their time and spatial scales [1,
4, 21]. Various efforts have tried to address this challenge. This includes
traditional approaches that incorporate domain knowledge to build and validate
model hypotheses using observations as well as sub-grid cloud modeling (known
as super-parameterization). Alternatively, recent machine learning approaches
use meteorological variables to model sub-grid clouds, thereby reducing the
computational cost of super-parameterization [3, 19, 16, 25].

In this paper, we extend [25], a data-driven approach to contribute to cloud
modeling, focusing on one of the main features used in energy balance
calculations: reflectance fields. We use Conditional Generative Adversarial
Networks [11] to generate these reflectance fields conditioned on
meteorological variables111The code is available on Github:
https://github.com/krisrs1128/clouds_dist. We suggest using these generated
images to extract important cloud parameters such as optical depth. We believe
our approach is a step towards building a data-driven framework that can
reduce the computational complexity in traditional cloud modeling techniques.

Our goal is to model reflectance fields, which in turn could be used as a
proxy for cloud optical depth, a major component of GCMs’ energy balance
computations [9, 5]. To do so, we leverage 3100 aligned sample pairs
$\mathcal{X}=\\{r_{i},m_{i}\\}$, where meteorological data $m_{i}$ are
collocated with reflectances $r_{i}$. Each $m_{i}$ is a $44\times 256\times
256$ matrix, representing 42 measurements from MERRA-2 [6] (see Table 1) along
with longitude and latitude to account for the Earth’s movement relative to
the satellite222As the earth rotates, the actual geographical locations on
Earth change pixel position in the data.. On the other hand $r_{i}$, is a
$3\times 256\times 256$ matrix representing each location’s reflectance at RGB
wavelengths (680, 550 and 450 $nm$) as measured by the Aqua dataset [17]. One
could consider working in a Supervised Learning setting to learn a
deterministic mapping $f:m_{i}\mapsto r_{i}$ ; however given the chaotic
nature of climate, we need not a point estimate of the potential cloud
distribution on earth, but rather an ensemble of likely scenarios given
initial conditions. This motivates a generative approach using conditional
GANs.

##  2 Modeling reflectance fields

###  2.1 Network

Architecture: motivated by Ronneberger et al. in [20], we use a U-Net as
conditional generator. The U-Net architecture helps our generator capture
global context, and skip connections allow localization. All of the
convolution modules in our U-Net implementation consist of the same building
blocks: a 3x3 convolutional layer followed by padding - which eliminates the
need for cropping - followed by batch normalization, leaky ReLU, and a dropout
layer with 0.25 probability.

Source of stochasticity: we introduce stochasticity in the generator only
through the dropout layers at both training and test times, i.e we do not use
noise input vectors. As observed by Isola et al. [10] dropout introduces
diversity in the output of conditional GANs

Checkerboard artifacts: a direct implementation of this generator results in
checkerboard artifacts, a result of the use of transposed convolutions to
upsample the feature maps in the U-Net expansion path [15]. We solve this
problem by replacing transposed convolution with a resize operation of the
feature maps using 2d nearest neighbor interpolation, followed by a
convolution as proposed in [15].

Discriminator: we use a multi-scale discriminator as proposed by Wang et al.
in [24]. This introduces 3 discriminators with identical network structure
operating on different input scales: one discriminator operates on the raw
input image, while the other two operate on the raw image downsampled by
factors of 2 and 4, using average pooling with a stride of 2. The motivation
behind using discriminators at different scales is to provide the generator
with better guidance both in the scale of global context and finer details in
the image.

###  2.2 Training

#### Training objectives

To train our generator, we use a weighted objective function composed of two
losses: a non-saturating GAN loss and a matching loss:

  1. 1.

_Non-saturating adversarial loss_. We experimented with two types of
adversarial loss: the hinge loss of [12] and the least square loss (LSGAN) of
[13]. We observe better performance with LSGAN. In figure 6, we also see that
least squares loss is more stable during training.

  2. 2.

_$L^{1}$ matching loss_. We use $L^{1}$ loss between generated and true
reflectance images, which encourages the generator to produce outputs close to
the observed images from a regression perspective. The $L^{1}$ loss has been
found to produce less blurry outputs than $L^{2}$ loss [10].

#### Optimizer

. As we explored various optimization strategies and regularization methods,
we observed significant improvement both in terms of convergence and in the
quality of the generated output by using the Extra-Adam333See code at
https://github.com/GauthierGidel/Variational-Inequality-GAN method proposed
[7] compared to Adam and SGD, see Figure 7.

##  3 Results and Discussion

###  3.1 Visual analysis

Our U-Net generator, trained against a multi-scale discriminator and optimized
by Extra-Adam, is able to generate visually appealing CRFs that are difficult
to distinguish from true samples. On a validation set of ground truth images,
we obtain an $\ell^{2}$ loss $\sim 0.027$. Figure 1 shows 4 different pairs
$(G(m_{i}),r_{i})$: we can see that the model is able to pick up large-scale
cloud structures as well as the continents and oceans beneath them. Although
not as precise as the ground-truth $r_{i}$, the generated samples exhibit
similar global composition as well as local structures.

In Figure 2, we generate 3 reflectance fields from the same conditioning
measurements. We notice a consistent global pattern in the three samples, with
variations visible in finer details. In order to quantitatively measure
diversity across generations, we fix the validation set to 5 samples that are
selected manually to capture different regions of the rotating earth and
generate 15 samples in total: 3 for each validation sample. For each set, we
compute 3 metrics: pixel-wise mean, standard deviation and inter-quartile
range across samples. Figure 2 shows that the model can obtain high image
quality and proximity to the original distribution, but only the cost of low
diversity.

Our model still has limitations, such as blurriness and small size
checkerboard artifacts. We believe the reasons for this are:

  1. 1.

_More training samples_ are needed to represent such a high dimensional
distribution, i.e 3100 samples are not enough to train a deep U-Net generator
($\sim 1.4$ million parameters) and discriminator ($\sim 8.3$ million
parameters).

  2. 2.

_More hyperparameter tuning_ , including architectural choices of the
generator and discriminator to ensure the right capacity balance that lead to
a long lasting game and avoid prematurely saturated learning.

  3. 3.

_Further training_ – we can see that the discriminator loss still slightly
oscillates after saturation points and eventually decreases with number of
steps as shown in 10(c).

(a)

(b)

(c)

(d)

Figure 1: 4 inferences obtained from our trained model. Generated images are
on the left and the corresponding true images on the right. We can see how
composition is preserved, most large cloud fields have similar shapes but
differ in the details. Figure 2: Reflectance fields generated by conditioning
on the same input (noise comes from dropout, which is kept at test time).
Figure 3: Plot of the inter-quartile distance, mean, and standard deviation of
15 generated samples at different steps during the model’s training.

###  3.2 Spectral analysis

Although visual inspection techniques can give insight into GAN performance,
it is an expensive, cumbersome, and subjective measure [2]. We address this
issue by comparing the frequency spectrum of true and generated samples using
2D Discrete Fourier transform (DFT). This allows us to compare the images’
geometric structures by examining the contribution of frequency components
[8]. We compare the the magnitudes of the 2D DFT calculated from the grayscale
versions of the true and generated images, and compare the histograms of the
calculated magnitudes, their means, variances and the logarithmic average
$L^{2}$ distances. In figure 4 we observe that our generated images have
consistent and similar DFT distributions to those of their corresponding true
images, with a very small average $L^{2}$ distance.

Figure 4: Comparison of DFTs (from left to right: image, frequency magnitudes,
histogram of magnitudes), with real data in the top row and a generated
reflectance field (from the real data’s associated measurements)in the bottom
one

##  4 Conclusion and future work

We show that using conditional GANs to model CRFs can be an effective approach
towards building a data-driven framework. We think our approach could
significantly help improve the computation time of clouds modeling in global
climate models.
Future work includes increasing the size of the dataset and exploitation of
the temporal structure in our data in two ways: by adding date and time as
extra labels to the input variable, and by using temporal cross validation
[23] to validate our generator’s ability to predict possible changes in cloud
distribution over time. We also plan to increase the diversity in the
generated ensembles by incorporating input noise channels as an extra source
of stochasticity. To address what we suspect to be mode collapsing in our
network (the matching loss discourages the exploration of other potential
modes in the data) we suggest using staged training where an adaptive weight
for the matching loss encourages the generator to regress onto true images
during early stages of the training, eventually decreasing to zero as training
progresses.

## References

  * Arakawa [2004] Akio Arakawa.  The cumulus parameterization problem: Past, present, and future.  _Journal of Climate_ , 17(13):2493–2525, 2004\.
  * Borji [2019] Ali Borji.  Pros and cons of gan evaluation measures.  _Computer Vision and Image Understanding_ , 179:41–65, 2019\.
  * Brenowitz & Bretherton [2018] Noah D Brenowitz and Christopher S Bretherton.  Prognostic validation of a neural network unified physics parameterization.  _Geophysical Research Letters_ , 45(12):6289–6298, 2018.
  * Bretherton [2015] Christopher S Bretherton.  Insights into low-latitude cloud feedbacks from high-resolution models.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 373(2054):20140415, 2015.
  * Corti et al. [2005] T Corti, BP Luo, Th Peter, H Vömel, and Q Fu.  Mean radiative energy balance and vertical mass fluxes in the equatorial upper troposphere and lower stratosphere.  _Geophysical research letters_ , 32(6), 2005.
  * Gelaro et al. [2017] Ronald Gelaro, Will McCarty, Max J Suárez, Ricardo Todling, Andrea Molod, Lawrence Takacs, Cynthia A Randles, Anton Darmenov, Michael G Bosilovich, Rolf Reichle, et al.  The modern-era retrospective analysis for research and applications, version 2 (merra-2).  _Journal of Climate_ , 30(14):5419–5454, 2017\.
  * Gidel et al. [2018] Gauthier Gidel, Hugo Berard, Gaëtan Vignoud, Pascal Vincent, and Simon Lacoste-Julien.  A variational inequality perspective on generative adversarial networks.  _arXiv preprint arXiv:1802.10551_ , 2018.
  * Gonzalez & Woods [1991] CR Gonzalez and ER Woods.  Digital image processing, addison-wesley publishing company.  1991\.
  * Hartmann et al. [1992] Dennis L Hartmann, Maureen E Ockert-Bell, and Marc L Michelsen.  The effect of cloud type on earth’s energy balance: Global analysis.  _Journal of Climate_ , 5(11):1281–1304, 1992\.
  * Isola et al. [2017a] Phillip Isola, Jun-Yan Zhu, Tinghui Zhou, and Alexei A Efros.  Image-to-image translation with conditional adversarial networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 1125–1134, 2017a.
  * Isola et al. [2017b] Phillip Isola, Jun-Yan Zhu, Tinghui Zhou, and Alexei A Efros.  Image-to-image translation with conditional adversarial networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 1125–1134, 2017b.
  * Lim & Ye [2017] Jae Hyun Lim and Jong Chul Ye.  Geometric gan.  _arXiv preprint arXiv:1705.02894_ , 2017.
  * Mao et al. [2017] Xudong Mao, Qing Li, Haoran Xie, Raymond YK Lau, Zhen Wang, and Stephen Paul Smolley.  Least squares generative adversarial networks.  In _Proceedings of the IEEE International Conference on Computer Vision_ , pp. 2794–2802, 2017.
  * North et al. [1981] Gerald R North, Robert F Cahalan, and James A Coakley Jr.  Energy balance climate models.  _Reviews of Geophysics_ , 19(1):91–121, 1981\.
  * Odena et al. [2016] Augustus Odena, Vincent Dumoulin, and Chris Olah.  Deconvolution and checkerboard artifacts.  _Distill_ , 1(10):e3, 2016.
  * O’Gorman & Dwyer [2018] Paul A O’Gorman and John G Dwyer.  Using machine learning to parameterize moist convection: Potential for modeling of climate, climate change, and extreme events.  _Journal of Advances in Modeling Earth Systems_ , 10(10):2548–2563, 2018.
  * Platnick et al. [2017] S Platnick, KG Meyer, MD King, G Wind, N Amarasinghe, B Marchant, GT Arnold, Z Zhang, PA Hubanks, RE Holz, et al.  The modis cloud optical and microphysical products: Collection 6 updates and examples from terra and aqua, ieee t. geosci. remote, 55, 502–525, 2017.
  * Ramanathan et al. [1989] VLRD Ramanathan, RD Cess, EF Harrison, P Minnis, BR Barkstrom, E Ahmad, and D Hartmann.  Cloud-radiative forcing and climate: Results from the earth radiation budget experiment.  _Science_ , 243(4887):57–63, 1989.
  * Rasp et al. [2018] Stephan Rasp, Michael S Pritchard, and Pierre Gentine.  Deep learning to represent subgrid processes in climate models.  _Proceedings of the National Academy of Sciences_ , 115(39):9684–9689, 2018.
  * Ronneberger et al. [2015] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In _International Conference on Medical image computing and computer-assisted intervention_ , pp. 234–241. Springer, 2015.
  * Schneider et al. [2017] Tapio Schneider, João Teixeira, Christopher S Bretherton, Florent Brient, Kyle G Pressel, Christoph Schär, and A Pier Siebesma.  Climate goals and computing the future of clouds.  _Nature Climate Change_ , 7(1):3–5, 2017.
  * Stocker et al. [2013] Thomas F Stocker, Dahe Qin, Gian-Kasper Plattner, Melinda Tignor, Simon K Allen, Judith Boschung, Alexander Nauels, Yu Xia, Vincent Bex, Pauline M Midgley, et al.  Climate change 2013: The physical science basis.  _Contribution of working group I to the fifth assessment report of the intergovernmental panel on climate change_ , 1535, 2013.
  * Varma & Simon [2006] Sudhir Varma and Richard Simon.  Bias in error estimation when using cross-validation for model selection.  _BMC bioinformatics_ , 7(1):91, 2006.
  * Wang et al. [2018] Ting-Chun Wang, Ming-Yu Liu, Jun-Yan Zhu, Andrew Tao, Jan Kautz, and Bryan Catanzaro.  High-resolution image synthesis and semantic manipulation with conditional gans.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 8798–8807, 2018.
  * Yuan et al. [2019] T. Yuan, H. Song, D. Hall, V. Schmidt, K. Sankaran, and Y. Bengio.  Artificial intelligence based cloud distributor (ai-cd): probing clouds with generative adversarial networks.  _AGU Fall Meeting 2019_ , 2019.

##  Appendix A Data processing

During data processing, we Winsorize reflectance data to remove artifacts that
are present in some training samples due to sensor noise, clipping values to
the 95th percentile for each channel. This is necessary as meteorological
variables have different scales, see Figure 5. We standardize channels to have
values in $\left[-1,1\right]$ and zero mean. In order to to avoid introducing
unnecessary bias from values outside the earth disk, we first crop the images
to cut off most of these values and then upsample them again to their original
size using 2D nearest neighbors, replacing the remaining values with $-3$
(mean - 3x standard deviation) to avoid introducing any unnecessary bias in
the data distribution.

We used running statistics to compute the summary statistics of the data due
to the huge size of input tensors, which do not fit in 16GB GPU memory at one
time. We also increased the number of data loader workers to 12; this
accelerates the data loading process by $6\times$.

Table 1: Description of input components Name | Description | Number of channels
---|---|---
U, V | Wind components in 10 atmospheric levels | 20
T | Temperature in 10 atmospheric levels | 10
RH | Relative-humidity in 10 atmospheric levels | 10
SA | Scattering angle | 1
TS | Surface Temperature | 1
Lat, Long | Latitude and Longitude | 2
Figure 5: Histograms of six input variables shows the variance in scales.

##  Appendix B Hyper-parameter comparisons

(a) L1 matching loss

(b) Total weighted generator loss.

Figure 6: Comparison between the hinge loss (green) and the least squares loss
(purple) on model training stability and convergence, we observe that the
latter performs better both in optimization of the $L^{1}$ loss and the total
weighted generator loss. We configured Adam and ExtraAdam to use $\beta_{1}$ =
0.5, and $\beta_{2}$ = 0.99 in all experiments.

(a) $L^{1}$ matching loss.

(b) Generator adversarial loss.

(c) Discriminator loss.

(d) Adam

(e) ExtraSGD

(f) ExtraAdam

(g) Real earth

Figure 7: The losses of 3 experiments with 3 different optimizers: Adam,
ExtraSGD, and ExtraAdam, along with generated outputs for each experiment
conditioned on a fixed meteorological input. ExtraAdam shows better
convergence, less oscillating losses, and more visually appealing output
relative to Adam and ExtraSGD.

#### Regression vs. hallucinated features

The $\lambda_{1}/\lambda_{2}$ ratio in the generator’s weighted objective
function plays an important role in the behavior of our generator; experiments
with large ratios in the range of $\left[1,10\right]$ behave like supervised
models where we regress the generated images with $L^{1}$ loss, while small
ratios of $\leq 0.5$ tend to give the generator more freedom to explore the
distribution of interest, without being penalized for not matching low
frequency details. This causes the generator to hallucinate features that do
not exist in the true images (Figure 8). This behavior matches our
expectations: 3100 samples is not sufficient to learn the conditional
distribution of such variable and high-dimensional data.

(a) $L^{1}$ matching loss

(b) $\lambda_{1}/\lambda_{2}=5$

(c) $\lambda_{1}/\lambda_{2}=0.2$

(d) real

Figure 8: Two experiments with two different $\lambda_{1}/\lambda_{2}$ ratios
shows stability of the discriminator loss for larger ratios – resulting in
more realistic images, see (b) –, while small ratios produce hallucinated and
unrealistic images, see (c).

#### Architecture Summary

The main components of our architecture are summarized in Figure 9.

Figure 9: The overall architecture includes two models, a U-Net generator and
multiscale discriminator, and the optimization objective combines a cGAN loss
with an $L^{1}$ matching loss.

#### Sharpness of generated images

Generating sharp images that can show complicated and detailed clouds
structures such as spinning clouds is both important and challenging. We
address this challenge by carefully choosing the discriminator learning rate
to avoid saddle point convergence and non-convergence, which result in bad
generations. We hypothesize that the generator might have not been trained
long enough to learn such micro-level details and thus generate blurry output.
In most of these cases, we observe the generator is the dominant player, and
the discriminator is fooled at saturation points, see the example shown in
Figure 10.

(a) $L^{1}$ matching loss

(b) Generator loss

(c) Discriminator loss

(d) fake

(e) real

Figure 10: An example experiment with a large discriminator learning rate
(0.001) shows early saturation of the discriminator loss after $\sim$ 10K
steps, resulting in an idle generator loss starting from the same point, the
generated sample thus looks blurry when compared to the real one.


# Locality Sensitive Hashing for Network Traffic Fingerprinting

Nowfel Mashnoor, Jay Thom, Abdur Rouf, Shamik Sengupta, Batyr Charyyev
Department of Computer Science and Engineering, University of Nevada, Reno,
USA  Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The advent of the Internet of Things (IoT) has brought forth additional
intricacies and difficulties to computer networks. These gadgets are
particularly susceptible to cyber-attacks because of their simplistic design.
Therefore, it is crucial to recognise these devices inside a network for the
purpose of network administration and to identify any harmful actions. Network
traffic fingerprinting is a crucial technique for identifying devices and
detecting anomalies. Currently, the predominant methods for this depend
heavily on machine learning (ML). Nevertheless, machine learning (ML) methods
need the selection of features, adjustment of hyperparameters, and retraining
of models to attain optimal outcomes and provide resilience to concept drifts
detected in a network. In this research, we suggest using locality-sensitive
hashing (LSH) for network traffic fingerprinting as a solution to these
difficulties. Our study focuses on examining several design options for the
Nilsimsa LSH function. We then use this function to create unique fingerprints
for network data, which may be used to identify devices. We also compared it
with ML-based traffic fingerprinting and observed that our method increases
the accuracy of state-of-the-art by 12% achieving around 94% accuracy in
identifying devices in a network.

###### Index Terms:

Network Traffic Fingerprinting, Locality-Sensitive Hashing, Device
Identification, Internet of Things.

††publicationid: pubid: 979-8-3503-4693-0/23/$31.00 ©2023 IEEE

##  I Introduction

Network traffic fingerprinting studies inferential methods that take the
network traces of a group of devices as input and produce information about
those devices, their users, and the environment [1]. Network traffic
fingerprinting is an important tool for network security and management as it
enables to characterize traffic flows, identify the type of application such
as streaming or web browsing, and match the traffic to the application or the
device that generated it. Furthermore, it enables detecting and identifying
anomalous flows and prevent malicious activities in a network. Traffic
fingerprinting has been studied over a couple of decades focusing on operating
system identification [2], website detection [3], and device fingerprinting
[4], etc. However, it regained the interest of researchers with the emergence
of the Internet of Things (IoT) devices and the realization of smart systems.

Previous studies [5, 6, 7] use network traffic fingerprinting for device
identification. However, these traffic fingerprinting approaches rely on
machine learning (ML) and require careful selection of features from network
traffic data to prevent over-fitting/under-fitting the data and reduce
computational overhead. There exist studies on feature selection for ML-based
traffic fingerprinting that use optimization approaches such as Genetic
Algorithms [6], Grey Wolf Optimization [8], and Artificial Bee Colony [9]. ML-
based approaches also require hyperparameter tuning and retraining the models
if there exist newly available data, and concept drift in a network due to
firmware and configuration update. There exist studies on tuning and comparing
ML models for network traffic fingerprinting [10, 7, 11]. All these are
computationally costly and add complexity to ML-based networking tools.

In our previous studies, we proposed using locality-sensitive hash (LSH) based
network traffic fingerprinting for device identification [12, 13], anomaly
detection [14], and user interaction classification [15]. Locality-sensitive
hash is different from cryptographic hashes and produces similar hash values
for similar inputs. Thus we can use the hash similarity to identify devices by
comparing the hash value of the network traffic data to hash values from
different devices stored in a database. Since traffic fingerprinting with LSH
uses the digest (i.e., hash) of the network traffic data we do not need to
select and extract features from data. Also fingerprinting depends on hash
similarity thus traffic fingerprinting with LSH does not have complex
parameters that need to be tuned. In [12, 13] we proposed LSIF (Locality-
Sensitive IoT Fingerprinting) to identify IoT devices with LSH-based network
traffic fingerprinting. LSIF uses Nilsimsa [16] as the underlying LSH function
and achieves comparable results with state-of-the-art machine learning-based
methods. In this paper, we extend LSIF and propose LSIF-R (i.e., LSIF
Reloaded), by exploring the alternative design of Nilsimsa. Specifically, we
alter the implementation of Nilsimsa by trying different parameter
configurations such as window size, accumulator size, etc., making the hash
function more tunable. This will enable LSIF-R to be more robust to concept
drift observed in a network, in addition to increasing the accuracy of device
identification with traffic fingerprinting. We also compare LSIF-R with state-
of-the-art IoTSentinel [5] and make a preliminary analysis of overhead and
computational resource requirements.

Contributions of this paper are as follows:

  * •

We extensively explore design alternatives for LSH function Nilsimsa for
network traffic fingerprinting without feature extraction from the data,
tuning the parameters of the model, and retrain the model.

  * •

We propose a device identification system LSIF-R that uses LSH-based network
traffic fingerprinting. LSIF-R is lightweight and robust to concept drifts.

  * •

We compare LSIF-R with the state-of-the-art IoT device identification method
IoTSentinel.

##  II Related Works

In this section, we provide related works on network traffic fingerprinting,
IoT device identification, and studies that use locality-sensitive hashing for
various applications. We also briefly discuss how this work differs from our
prior studies and existing related works.

###  II-A Network Traffic Fingerprinting

Previous studies used network traffic fingerprinting for the operating system
and device identification [2, 5, 6], mobile application classification [17,
18], and monitoring/anomaly detection [19, 20]. Researchers have focused on
different aspects of traffic fingerprinting such as feature selection [6, 8,
9], model tuning [10, 7], user privacy [21, 22], fingerprinting in NATed
network, NAPT and VPN enabled networks [23, 24], fingerprinting with
background traffic [25], and fingerprinting on different states of the device
[6, 12].

###  II-B IoT Device Identification

Device identification is important for network management and monitoring [26].
Previous studies on IoT device identification with traffic fingerprinting use
machine learning models [10, 5, 11] and signature-based methods [27, 12, 13].
There exist studies [28, 29] that use natural language processing models to
identify devices with textual data extracted from mDNS, SSDP, DHCP, and HTTP
protocols, and from DNS queries. There also exist studies that use one/multi-
class classifiers and neural networks to identify devices with features
extracted from wireless communication [30, 7] and Bluetooth communication of
the devices [31, 32].

###  II-C Locality-Sensitive Hashing

Previous studies use locality-sensitive hashing (LSH) for spam detection [33],
malware classification [34], genome assembling [35], and content-based video
retrieval [36]. Brian et al. [37] show that LSH enables accurate and fast
performance with object classification, feature matching, and content-based
retrieval. Muhammad et al. propose a technique for spam email fingerprinting
at the packet level [33]. Locality-sensitive hashing was also used to recover
the biometric information of individuals in a corporate environment to detect
information leakage and similar events [38]. We also showed the effectiveness
of LSH function Nilsimsa in identifying IoT devices [12, 13], user interaction
classification [15], and monitoring and detecting malicious activities in a
network [14].

Different from existing studies on network traffic fingerprinting [2, 17, 20]
and existing studies on device identification [5, 6] we propose a new approach
with locality-sensitive hashing, which does not require feature
selection/extraction from data and does not have complex hyperparameters that
need to be tuned. Our study differs from existing studies on locality-
sensitive hashing [33, 34, 36] by focusing on different utilization (i.e.,
device identification) of the hash function. Also compared to our previous
studies [12, 15, 14] with Nilsimsa, our study differs by proposing an
alternative implementation of Nilsimsa making the device identification system
more flexible and at the same time robust to concept drifts observed in the
network.

Figure 1: Nilsimsa functionality

##  III Methodology

This section will provide a concise overview of locality-sensitive hashing,
with a special emphasis on Nilsimsa. We will provide comprehensive information
about the alternative configuration of Nilsimsa and its practical execution.
Next, we will provide comprehensive information on LSIF-R (Locality Sensitive
IoT Fingerprinting-Reloaded) as well as the specifics of the dataset used for
assessment purposes.

Nilsimsa is a locality-sensitive hash function that generates similar hash
values for similar inputs. There are alternative LSH functions such as ssdeep,
sdhash, and tlsh [39]. The functionality of Nilsimsa is presented in Figure 1.
Nilsimsa uses a sliding window of size 5 characters that slide one character
at a time. In each iteration it generates trigrams from data in the window and
passes those trigrams to tran54 hash function which generates value i between
0 and 255. Then i-th counter of the accumulator (i.e. array of integers of
size 256) is increased by 1. After iterating over the entire input file the
accumulator holds frequencies of the trigrams that have been found in the
input. For each bucket in the accumulator, the ratio of the value in that
bucket to the average values of the buckets will be assigned. Finally, if the
i-th ratio is greater than the median of the ratios, the i-th bit of the
Nilsimsa code will be set to 1 and it is set to 0, otherwise. The size of the
hash generated by Nilsimsa is 32 bytes and the similarity score of two hashes
ranges from -128 (i.e., completely different) to 128 (i.e., completely same).
The similarity score is calculated by subtracting 128 from the number of
similar bits in the hashes.

###  III-A LSIF-R

The Nilsimsa hashes generated from different traffic flows of similar devices
will have a high similarity score and a low score for hashes of traffic from
different devices. Using this similarity we can identify devices if we have
the signature (i.e., hash) of the traffic flow from devices in our database.
We altered the design of Nilsimsa and made it more tunable. Specifically, we
explored alternatives for parameters (i.e., window size, trigram, etc.). The
set of parameters that we explored is shown in Figure 1 with green color and
preceded with their parameter ID (PID), and parameter values are presented in
Table I. The default values of Nilsimsa are shown with blue color on Table I.
For PID 1 (window size) we explored the values of 6, 7, 8, 9, and 10
alternative to default value 5. For PID 2 we explored 3, 5, 7, 9 with the
restriction of $n-grams<windowsize$, as we can not generate n-gram where n is
less than the window size. In terms of the underlying hash function (PID 3)
and the threshold (PID 5) we used other classes of hashes such as SHA1-512,
MD4-5, etc., and statistical values such as first quartile (Q1), third
quartile (Q3), interquartile range (IQR), and standard deviation (std) as
alternative threshold metric. The default accumulator size in Nilsimsa is 256
and we used an accumulator size ranging from 16 to 1024 with an increment
value of 16.

Table I: Summary of the parameter set

Parameter ID | Values
---|---
window size (PID 1) |  5, 6, 7, 8, 9, 10
n-grams (PID 2) |  3, 5, 7, 9
hash (PID 3) |  TRAN53, SHA512, SHA384, SHA256, SHA224
SHA1, MD5, MD4, FNV164, FNV1A32, FNV1A64,
MMH3, CRC32, ADLER32, WHIRLPOOL
accumulator size (PID 4) | [16 : 1024], $\delta$ = 16 — 256
threshold (PID 5) | mean, median, mode, IQR, Q1, Q3, std

Figure 2: System Design of LSIF-R

In our previous studies [12, 13] we proposed a device identification system
that uses Nilsimsa with its default parameters. In this paper, we propose IoT
device identification system LSIF-R (Locality Sensitive IoT Fingerprinting
Reloaded), that uses the tunable Nilsimsa described above. The design of the
LSIF-R is presented in Figure 2. First, LSIF-R computes a set of hash values
of traffic flows generated by the devices and stores them as signatures in a
database, with the device names as labels of signatures. These signatures are
generated by tunable Nilsimsa. Here traffic flow is a sequence of network
packets sent and received by a particular IoT device in a certain time
interval. Note that creating the initial signature database is performed only
once, thus, it does not impose future overheads. When we try to identify the
device we generate a hash of its traffic flow and compare it to signatures
stored in a database. Then, the device with the highest average hash
similarity score is selected as the predicted label for the device. If the
highest similarity score is less than a threshold value (i.e., min threshold
for similarity score), then the device is labeled as unknown (i.e., the
signature of the device does not exist in the database). We can use average
traffic flow similarity as the min similarity threshold value or we can employ
a device-specific threshold based on the variance of the signatures of the
device.

Compared to existing machine learning-based device identification systems [6,
5], LSIF-R does not require feature selection and extraction from data, as it
relies on the hash of the traffic flow. This provides some advantages over
existing machine learning-based systems in terms of overhead as described in
the next sections of our paper. Also it is simple to update LSIF-R if there is
new available data. For instance, if we update the firmware of a particular
device then its network characteristics might change and it needs to be
reflected in the predictive model. If we want to integrate a new device into
our system, again with ML-based device identification we need to re-train the
model whereas with LSIF-R we only need to append signatures of the new device
to the end of the signature database. Thus, LSIF-R has clearly some benefits
over machine learning based methods in terms of updating the model.

###  III-B Dataset and Metrics

The dataset that we used for evaluation is from study [5]. It contains the
network traffic of 23 IoT devices collected during the setup phase of the
device. The dataset includes broad range of IoT devices ranging from sensors
to home appliances. For each device, there exist 20 measurements. Since we use
the hash of the traffic flow, we removed IP and MAC addresses associated with
the devices to remove any bias that may occur in classifying the devices. We
used 4-fold cross-validation for evaluations and performed 5 random
experiments and presented the average. The evaluations were in terms of
precision (i.e., $TP/(TP+FP)$), recall (i.e., $TP/(TP+FN)$), f1-score (i.e.,
$2/(1/precision+1/recall)$), and accuracy (i. e., $(TP+TN)/(TP+FP+TN+FN)$)
where TP, TN, FP, and FN stand for true positive, true negative, false
positive and false negative, respectively.

Figure 3: Device identification F1-score for default hash (Tran53) of Nilsimsa
and trying different values for rest of parameters. The red star is F1-score
for the default parameters of Nilsimsa. Figure 4: F1-score for different
classes of hash functions.

##  IV Evaluation and Results

In this section, we provide evaluation results of the proposed device
identification method LSIF-R in identifying 23 IoT devices described in
Section III-B. We first evaluate the LSIF-R by trying different parameters
from Table I. Then we focus on if there is any correlation between the
hyperparameter set and accuracy. Finally, we compare LSIF-R with the existing
state-of-the-art method IoTSentinel [5].

###  IV-A Design alternatives of Nilsimsa

First, we started the evaluation by trying different hyperparameters from
Table I while keeping the underlying hash (Tran53) function of Nilsimsa
unchanged. Figure 3 presents results in terms of the F1-score. When Tran53 is
stable, we can increase the average accuracy by 3% from the accuracy of 83%
achieved with default parameters of Nilsimsa (shown with a red star in Figure
3). The change in accuracy increases even more when we try different hash
values instead of Tran53. Figure 4 shows the results when we try different
hash functions from other hash classes such as SHA, MD, FNV, etc. We can
achieve an average F1-score above 90% with other classes of hash functions.
For instance SHA1 achieves F1-score of 93%, MD4 achieves 93%, and FNV164 and
CRC32 achieves 92% and 94% of F1-scores respectively. These values show that
we can increase the accuracy of LSIF-R by 10% by trying different
hyperparameters for Nilsimsa and changing its design and implementation. While
a majority of functions achieve similar results, not all hash functions
achieve good results. For instance in Figure 5 we can see that hash functions
MMH3 and ADLER32 perform very poorly in terms of device identification for all
parameter combinations. However, we can also see that hash functions FNV164
and FNVA164 achieve slightly better performance for parameters between 0-1000
compared to the rest of the hash functions. This suggests that we can use the
ensemble of FNV164 and FNVA164 with other hash functions to increase accuracy
for parameters between 0-1000. This might be useful if parameters between
0-1000 are more lightweight compared rest of the parameter combinations. For
instance assume the $900^{th}$ parameter has window size and n-gram equal to 5
and 3 respectively and assume these values for parameter $2000^{th}$ are 9 and
3. In this case, parameter $900^{th}$ will be more lightweight than
$2000^{th}$ because $Combination(5,3)<Combination(9,3)$, and
$Combination(5,3)$ indicates how many trigrams will be generated with window
size equal to 5 and n-gram equal to 3. LSIF-R achieves higher accuracy on the
parameter $2000^{th}$ with single hash functions. However, if it can achieve
comparable performance with a combination of multiple hash functions at
parameter $900^{th}$ then ensembling the hash functions needs to be considered
as using the parameter combination $900^{th}$ is more lightweight than
parameter $2000^{th}$. Overall we observed that if we tune the other
parameters of Nilsimsa and use the hash functions such as SHA, MD instead of
default TRAN53 the accuracy of Nilsimsa-based LSIF-R reaches up to 94%. From
alternative hash functions for Nilsimsa, we observed that CRC32 achieves the
best result thus in the rest of the paper our evaluations will be based on
CRC32.

Figure 5: F1-score of LSIF-R with other hash functions. Figure 6: Correlation
of individual parameter with accuracy.

###  IV-B Correlation with hyperparameters

In Figure 5, if we remove MMH3, ADLER32, FNV164, and FNVA164 hashes, we can
see that the rest of the hashes have very similar behavior. We can also see
that for parameters above $1500^{th}$ the accuracy of LSIF-R reaches around
90% with all hashes. Thus, it is important to analyze which parameters achieve
high accuracy and if there is a correlation between accuracy and any specific
parameter. Figure 6 shows the correlation between the accuracy of LSIF-R and
accumulator size one of the hyperparameters of Nilsimsa from Table I. We can
see that there is no correlation between accuracy and accumulator size. We
observed similar behavior for other parameters such as window size, n-grams,
and threshold type, which are not shown due to space constraints. While there
is no correlation between accuracy and individual parameters, there might
exist some correlation between accuracy and a combination of multiple
parameters. These analyzes are important as the overhead of Nilsimsa might
change depending on the values of hyperparameters. As we discussed in
subsection IV-A the parameter pair (window size=5, n-gram=3) will have low
overhead compared to the parameter pair (window size=9, n-gram=3). Overall, we
did not observe any correlation between individual parameters and accuracy
however exploring the overhead for different parameter combinations and the
trade-off between accuracy and imposed overhead from parameter combinations
need to be explored in the future.

###  IV-C Comparison with related works

In this subsection we compared LSIF-R with the state-of-the-art machine
learning-based IoT device identification system IoTSentinel [5]. IoTSentinel
works as follows: it extracts the 23 network traffic features from the first
12 packets of the device and generates a feature matrix and feeds it to the
machine learning model Random Forest. The features include the existence of
Link/Network/Transport layer protocols, data related to IP headers, packet
count, etc. A Random Forest model is trained for each device and it makes
binary classification (if traffic belongs to that device or not). If traffic
flow belongs to multiple or none of the devices as a result of random forest
classifications then edit distance discrimination is used to identify the
device. Different from the IoTSentinel, LSIF-R is simple and it is easy to
update the model if there is new available data. For instance, if we need to
update the device identification system due to a firmware update of that
device or if we want to integrate the new device into the system we will need
to update the whole underlying machine learning model of the system. Whereas
this is not the case with LSIF-R where we only need to update the signatures
of a particular device or append the signatures of a new device to the
signature database. Also, LSIF-R does not require feature selection as it only
takes the hash of the traffic flow. The feature selection requires expert
knowledge and it is computationally expensive. For instance, IoTSentinel uses
the first 23 packets to extract the features from data. The features include
the existence of a particular protocol in data or statistics about
source/destination port numbers and IP numbers. Since these features are
computed over the first 23 packets we need to store the information related to
23 packets and wait till those packets are generated, making feature selection
costly in terms of storage and time. This cost might be significant especially
if IoTSentinel is deployed on embedded devices with limited computational
resources. However, this is not the case with LSIF-R as it uses the hash of
the traffic flow and the hash can be generated as network traffic data is
generated. For instance in our previous study [13] we observed that the
processing time of Nilsimsa with default configuration is $3.3\times 10^{5}$
bytes/s which is significantly higher than the traffic rate of IoT devices
used in our evaluation (WeMoLink has the highest rate of the traffic 992
bytes/s). Since the processing speed of Nilsimsa is significantly higher than
the traffic rate of devices we can generate the hash of the traffic as it
comes with negligible overhead. These suggest that LSIF-R is more lightweight
than machine learning-based device identification systems in terms of feature
selection/extraction and updating the model.

Figure 7: Comparison of LSIF-R and IoTSentinel in identifying 23 devices. The
average and individual results are presented.

In terms of accuracy, IoTSentinel achieves average precision, recall, and
f1-score of 84%, 82%, and 81% respectively and these values for LSIF-R are
94%, 93%, and 93%. Figure 7 presents the average performance comparison of
IoTSentinel and LSIF-R and the performance comparison for individual devices.
We can see that LSIF-R achieves better accuracy in 17 devices out of 23, and
in terms of average accuracy LSIF-R has an average F1-score 12% higher than
the IoTSentinel. These findings show that LSIF-R is more lightweight and
achieves higher accuracy compared to the state-of-the-art machine learning-
based device identification system IoTSentinel.

##  V Conclusion

Network traffic fingerprinting is an important tool for network management and
monitoring. In this paper, we proposed network traffic fingerprinting with
locality-sensitive hashing. Specifically we explored design alternatives for
the hash function Nilsimsa and proposed LSIF-R a device identification method
that uses Nilsimsa-based network traffic fingerprinting. Since LSIF-R depends
on locality sensitive hashing it does not require feature selection and
extraction from traffic data. Also by design LSIF-R does not have complex
parameters that need to be tuned and it is relatively straightforward to
update it with new measurements. We compared the LSIF-R with IoTSentinel which
uses machine learning-based network traffic fingerprinting to identify
devices. We observed that LSIF-R has higher accuracy (by 12%) than IoTSentinel
and achieves 94% accuracy in identifying 23 IoT devices. In the future, we
want to explore the computational overhead of different parameter combinations
of Nilsimsa on LSIF-R and explore the trade-off between accuracy and imposed
overhead. We also want to explore self-tuning properties for LSIF-R to
mitigate the impact of concept drift observed in a network.

## References

  * [1] M. Conti, Q. Q. Li, A. Maragno, and R. Spolaor, “The dark side(-channel) of mobile devices: A survey on network traffic analysis,” _IEEE Commun. Surv. Tutor._ , vol. 20, no. 4, pp. 2658–2713, 2018.
  * [2] A. Aksoy, S. Louis, and M. H. Gunes, “Operating system fingerprinting via automated network traffic analysis,” in _2017 IEEE Congress on Evolutionary Computation_ , 2017, pp. 2502–2509.
  * [3] M. Shen, Y. Liu, S.Chen, L.Zhu, and Y.Zhang, “Webpage fingerprinting using only packet length information,” in _IEEE ICC 2019_ , pp. 1–6.
  * [4] S. V. Radhakrishnan, A. S. Uluagac, and R. Beyah, “Gtid: A technique for physical device and device type fingerprinting,” _IEEE Trans. Dependable Secure Comput._ , vol. 12, no. 5, pp. 519–532, 2015.
  * [5] M.Miettinen, S.Marchal, I.Hafeez, N.Asokan, A.Sadeghi, and S.Tarkoma, “Iot sentinel: Automated device-type identification for security enforcement in iot,” in _2017 IEEE 37th Int. Conf. Distrib. Comput. Syst._ , pp. 2177–2184.
  * [6] A. Aksoy and M. H. Gunes, “Automated iot device identification using network traffic,” in _2019 IEEE ICC_ , 2019, pp. 1–7.
  * [7] A. J. Pinheiro, J. de M. Bezerra, C. A. Burgardt, and D. R. Campelo, “Identifying iot devices and events based on packet length from encrypted traffic,” _Computer Communications_ , vol. 144, pp. 8 – 17, 2019.
  * [8] Al.Shorman, A.Faris, and H.Aljarah, “Unsupervised intelligent system based on one class support vector machine and grey wolf optimization for iot botnet detection,” _J Ambient Intell Humaniz Comput_ , vol. 11, p. 2809–2825, 2020.
  * [9] W. A. Ghanem and A. Jantan, “Training a neural network for cyberattack classification applications using hybridization of an artificial bee colony and monarch butterfly optimization,” _Neural Process. Lett._ , vol. 51, pp. 905–946, 2020.
  * [10] E. Anthi, L. Williams, M. Słowińska, G. Theodorakopoulos, and P. Burnap, “A supervised intrusion detection system for smart home iot devices,” _IEEE Internet Things J._ , vol. 6, no. 5, pp. 9042–9053, 2019.
  * [11] B. A. Desai, D. M. Divakaran, I. Nevat, G. W. Peter, and M. Gurusamy, “A feature-ranking framework for iot device classification,” in _COMSNETS 2019_ , 2019, pp. 64–71.
  * [12] B. Charyyev and M. H. Gunes, “Iot traffic flow identification using locality sensitive hashes,” in _IEEE ICC 2020_ , 2020, pp. 1–6.
  * [13] B. Charyyev and M. H. Gunes, “Locality-sensitive iot network traffic fingerprinting for device identification,” _IEEE Internet Things J._ , vol. 8, no. 3, pp. 1272–1281, 2021.
  * [14] B. Charyyev and M. H. Gunes, “Detecting anomalous iot traffic flow with locality sensitive hashes,” in _2020 IEEE Globecom_ , 2020.
  * [15] B. Charyyev and M. H. Gunes, “Iot event classification based on network traffic,” in _INFOCOM 2022 WKSHPS_ , 2020, pp. 854–859.
  * [16] E. Damiani, S. D. C. di Vimercati, S. Paraboschi, and P. Samarati, “An open digest-based technique for spam detection,” in _2004 Int. WKSHPS on Sec. in Paral. and Distrib. Systems_ , vol. 41, no. 8, 2004, pp. 74–83.
  * [17] V. Ede, Thijs, and et al., “Flowprint: Semi-supervised mobile-app fingerprinting on encrypted network traffic,” in _NDSS_ , vol. 27, 2020.
  * [18] H. D. Trinh and et.al., “Mobile traffic classification through physical control channel fingerprinting: A deep learning approach,” _IEEE Trans. Netw. Service Manag_ , vol. 18, no. 2, pp. 1946–1961, 2020.
  * [19] Y. Meidan and et.al., “N-baiot—network-based detection of iot botnet attacks using deep autoencoders,” _IEEE Pervasive Comput._ , vol. 17, no. 3, pp. 12–22, 2018.
  * [20] Y. Mirsky, T. Doitshman, Y. Elovici, and A. Shabtai, “Kitsune: An Ensemble of Autoencoders for Online Network Intrusion Detection,” _arXiv e-prints_ , p. arXiv:1802.09089, Feb. 2018.
  * [21] X. Zhang, J. Hamm, M. K. Reiter, and Y. Zhang, “Defeating traffic analysis via differential privacy: a case study on streaming traffic,” _Int. J. Inf. Secur._ , vol. 21, no. 3, pp. 689–706, 2022.
  * [22] K. Dietz, M. Mühlhauser, M. Seufert, N. Gray, T. Hoßfeld, and D. Herrmann, “Browser fingerprinting: How to protect machine learning models and data with differential privacy?” _Electron. Commun. EASST_ , vol. 80, 2021.
  * [23] S. Dong, Z. Li, D. Tang, J. Chen, M. Syn, and K. Zhang, “Your smart home can’t keep a secret: Towards automated fingerprinting of iot traffic with neural networks,” in _arXiv:1909.00104_.
  * [24] Y. Meidan, V. Sachidananda, Y. Elovici, and A. Shabtai, “Privacy-preserving detection of iot devices connected behind a nat in a smart home setup,” in _arXiv:1905.13430_.
  * [25] S. Marchal and et.al., “Audi: Toward autonomous iot device-type identification using periodic communication,” _IEEE J. Sel. Areas Commun._ , vol. 37, no. 6, pp. 1402–1412, 2019.
  * [26] J. Ortiz, C. Crawford, and F. Le, “Devicemien: Network device behavior modeling for identifying unknown iot devices,” in _Proceedings of the IoTDI_ , ser. IoTDI ’19, 2019, p. 106–117.
  * [27] L. Babun, H. Aksu, L. Ryan, K. Akkaya, E. S. Bentley, and A. S. Uluagac, “Z-iot: Passive device-class fingerprinting of zigbee and z-wave iot devices,” in _IEEE ICC 2020_ , 2020, pp. 1–7.
  * [28] N. Ammar, L. Noirie, and S. Tixeuil, “Network-protocol-based iot device identification,” in _2019 FMEC_ , pp. 204–209.
  * [29] F. Le, J. Ortiz, D. Verma, and D. Kandlur, _Policy-Based Identification of IoT Devices’ Vendor and Type by DNS Traffic Analysis_. Cham: Springer International Publishing, 2019, pp. 180–201.
  * [30] M. R. Shahid, G. Blanc, Z. Zhang, and H. Debar, “Iot devices recognition through network traffic analysis,” in _2018 IEEE Trans. Big Data_ , 2018, pp. 5187–5192.
  * [31] H. Aksu, A. S. Uluagac, and E. Bentley, “Identification of wearable devices with bluetooth,” _IEEE Trans. Sustain. Comput_ , pp. 1–1, 2018.
  * [32] T. Gu and et. al., “Bf-iot: Securing the iot networks via fingerprinting-based device authentication,” in _IEEE MASS 2018_ , pp. 254–262.
  * [33] M. N. Marsono, “Packet-level open-digest fingerprinting for spam detection on middleboxes,” _Int. J. Netw. Manag._ , vol. 22, no. 1, pp. 12–26, 2012.
  * [34] L. Friborg and S. C. Peiser, “Malware classification using locality sensitive hashing and neural networks,” 2019.
  * [35] K.Berlin, S.Koren, C.Chin, J.Drake, J. Landolin, and A. Phillippy, “Assembling large genomes with single-molecule sequencing and locality-sensitive hashing,” _Nat. Biotechnol_ , vol. 33, no. 6, pp. 623–630, 2015.
  * [36] S. Hu, “Efficient video retrieval by locality sensitive hashing,” in _2005 IEEE ICASSP_ , vol. 2. IEEE, 2005, pp. ii–449.
  * [37] B. Kulis and K. Grauman, “Kernelized locality-sensitive hashing for scalable image search,” in _2009 IEEE CV_ , pp. 2130–2137.
  * [38] A.Alruban, N.Clarke, F.Li, and S.Furnell, “Biometrically linking document leakage to the individuals responsible,” in _Trust, Priv. and Sec, in Dig. Business: 15th Int. Conf., TrustBus_. Springer, 2018, pp. 135–149.
  * [39] V. G. Martinez, F. H. Alvarez, and L. H. E. Fernando, “State of the art in similarity preserving hashing functions,” _digital.csic.es_ , 2014.


# Monte Carlo modelling of the linear Breit-Wheeler process within the GEANT4
framework

R. A. Watt  S. J. Rose  B. Kettle  S. P. D. Mangles  The John Adams Institute
for Accelerator Science, Imperial College London, London, SW7 2AZ, United
Kingdom

(February 9, 2023)

###### Abstract

A linear Breit-Wheeler module for the code Geant4 has been developed. This
allows signal-to-noise ratio calculations of linear Breit-Wheeler detection
experiments to be performed within a single framework. The interaction between
two photon sources is modelled by treating one as a static field, then photons
from the second source are sampled and tracked through the field. To increase
the efficiency of the module, we have used a Gaussian process regression,
which can lead to an increase in the calculation rate by a factor of up to
1000.

To demonstrate the capabilities of this module, we use it to perform a
parameter scan, modelling an experiment based on that recently reported by
Kettle et al. [1]. We show that colliding $50\,$fs duration $\gamma$-rays,
produced through bremsstrahlung emission of a $100\,$pC, $2\,$GeV laser
wakefield accelerator beam, with a $50\,$ps X-ray field, generated by a
germanium burn-through foil heated to temperatures $>\,150\,$eV, this
experiment is capable of producing $>1\,$ Breit-Wheeler pair per shot.

††preprint: APS/123-QED

##  I  Introduction

The linear Breit-Wheeler (BW) process is the annihilation of two photons to
produce an electron positron pair ($\gamma\gamma\rightarrow e^{+}e^{-}$) and
is the simplest mechanism by which matter can be generated from light [2]. The
process was first predicted in 1934, however, despite the long time that has
passed, the annihilation of real photons has never been directly observed in
the laboratory. This is due to the difficulty in generating the high energy-
density photon sources required to overcome the relatively small cross-
section. However, such photon sources are routinely found in astrophysical
environments, and the BW process is predicted to play an important role in a
range of phenomena, with examples including gamma ray bursts and the emission
from quasars [3, 4]. Also for photons propagating through the intergalactic
medium there is a high energy ($<100\,\mathrm{GeV}$) cut-off in the cosmic
gamma ray spectrum observed at Earth due to BW annihilation with the cosmic
microwave background [5, 6]. This mechanism was thought to be well understood.
However, recent observations have found a larger number of high energy photons
from quasar 3C 279 reaching Earth than expected [7]. This demonstrates the
need for more experimental work, to better understand these astrophysical
environments.

Recently, the STAR collaboration [8] have reported the observation of the BW
process from quasi-real photons in the peripheral collisions of high energy
ions. However, observing the annihilation of real photons remains elusive.
With the advance of high power laser systems, generating high energy-density
photon sources in the laboratory has become possible. This has led to the
proposal of several laser based BW detection experimental schemes utilising
current and future facilities. Pike et al. [9] have suggested a scheme
involving the generation of a thermal X-ray field with a temperature of $\sim
300\,$eV, using a laser heated hohlraum. A beam of GeV gamma rays, produced by
bremsstrahlung emission of a laser wakefield accelerated (LWFA) electron beam
in a high Z material target, then interacts with the X-ray photons. To produce
the thermal X-ray field requires a facility such as the National Ignition
Facility (NIF). Pike et al. predict up to $10^{5}$ pairs can be produced in a
single shot. A second scheme involving the interaction between two symmetrical
photon sources has been proposed by Ribeyre et al. [10]. With symmetrical
sources, the photon energy required to produces the rest mass of the
$e^{+}e^{-}$ pair is reduced to the MeV scale. Ribeyre et al. have compared
multiple sources and found synchrotron emission from a highly energetic
electron beam in an intense laser field and bremsstrahlung emission to be the
most efficient methods. A pair yield of $\sim 10^{4}$ per shot is predicted
with a relatively modest laser energy of $100\,\mathrm{J}$. Finally, building
on this symmetric setup, Drebot et al. [11] have investigated a scheme
involving the interaction of Compton sources produced by $\sim 260\,$MeV
electron beams and joule class lasers. The predicted yield of $\sim 10^{-4}$
pairs per shot is far lower than other experimental schemes, however, due to
the low laser energy required, a high repetition rate system can be used. The
BW process would then be detected over many shots.

In the Spring of 2018 a BW detection campaign took place at the Gemini laser
of the Central Laser Facility at the Rutherford Appleton Laboratory, in the
UK. The experiment was based on the scheme proposed by Pike et al. [9] and a
diagram of the experimental setup is shown in figure 1. A more in depth review
of this experiment can be found in Kettle et al. [1]. An electron beam
generated by a LWFA interacted with a $1\,$mm bismuth target, producing a
$\sim 50\,$fs pulse of gamma rays with a spectrum extending up to 100s of MeV.
This was directed through a $\sim 50\,$ps duration X-ray field, generated by
the thermal emission of a laser heated germanium foil which emits strongly
between 1 keV and 2 keV [12]. Both beams of the Gemini laser are capable of
providing 15 J, which is low compared to a high energy facility such as NIF at
$2\,$MJ. However, the Gemini laser has a much higher repetition rate of 0.05
Hz compared to 1 shot per day. With a large number of shots a better
characterisation of the background can be performed as well as a statistical
analysis of the interaction.

Figure 1: Schematic of Gemini linear BW detection experiment. Starting from
the left: a LWFA generates an electron beam which is converted into a gamma
ray beam through bremsstrahlung emission in a thin bismuth converter foil. A
tungsten collimator and block are placed in the beam path to remove highly
divergent gamma rays and those directed towards the X-ray foil respectively. A
large number of BH pairs are produced in the converter foil, collimator and
block. These are removed with an on-axis magnet before the interaction zone.
The gamma rays interact with an X-ray field, generated by a laser heated
germanium foil, producing BW pairs. The residual gamma rays continue on axis
to a spectrometer. The BW pairs pass through a magnetic chicane to single
particle detectors, situated behind lead shielding (not shown). Diagram
provided by E. Gerstmayr.

All linear BW detection experiments require the generation of high energy
($>0.511\,$MeV) photons. As a result, the environment is inherently noisy.
There are two main sources of noise: the direct detection of photons by the
particle detectors, and photons interacting with the experimental setup,
producing background $e^{+}e^{-}$ pairs through the Bethe-Heitler (BH) process
[13] (the annihilation of a photon with a virtual photon in the nuclear field
of an atom). For a successful experiment, the ratio of signal BW pairs to the
background noise should be maximised. An estimation of the background noise
can be obtained by performing full scale simulations of the passage of high
energy particles through the experimental setup. Several publicly available
Monte Carlo particle tracking codes exist for this purpose, including Geant4,
Fluka and MCNP [14, 15, 16]. For the 2018 Gemini campaign Geant4 was used to
analyse the background noise [1]. However, the linear BW process is not
included within the standard Geant4 physics package. It is advantageous to
perform both signal and noise calculations within a single framework, and the
aim of this work has been to develop a linear BW module for Geant4 to enable
this.

Geant4 models the passage of individual particles through matter. To include
the linear BW process within Geant4, it must be modelled within this same
framework. This has been achieved by treating one photon source as a static
photon field. Individual photons from the second source are then sampled and
tracked through this field. In this work, we will refer to a photon from the
field as static and a photon tracked through the field as dynamic. Using this
method, the temporal evolution of the static photon source cannot be accounted
for. This module is therefore suited for modelling experiments with asymmetric
photon sources, in which one is constant in time, to a good approximation,
over the full interaction. Figure 1 shows an example of such an experiment.

##  II Module Overview

As discussed in section I, this module involves tracking dynamic photons
through a static photon field. The static photon field is represented by a
physical volume within the Geant4 computational domain. It is fully defined by
$n(\epsilon,\theta,\phi)$, the number of photons per unit volume with energy
between $\epsilon$ and $\epsilon+\mathrm{d}\epsilon$, and travelling at an
angle between $\theta+\mathrm{d}\theta$ and $\phi+\mathrm{d}\phi$. For a
dynamic photon which enters a static photon field the calculation of the
interaction involves two steps:

  1. 1.

Calculating the probability of the interaction occurring.

  2. 2.

Calculating the dynamics of the interaction.

A dynamic photon enters the static photon field at some arbitrary angle in the
simulation frame. To perform the calculation of step 1, we transform
$n(\epsilon,\theta,\phi)$ into a frame where the z-axis is the direction of
propagation of the dynamic photon. We will refer to this frame as the dynamic
photon frame. For the linear BW process to occur, the total energy in the
centre-of-mass (c.m.) frame must be larger than twice the rest mass energy of
an electron. This sets the following interaction threshold

$s=2\epsilon E(1-\textrm{cos}\,\theta)>4m^{2}c^{4}$ (1)

where $E$ is the energy of the dynamic photon. If $n(\epsilon,\theta,\phi)=0$
for all $\epsilon$, $\theta$ and $\phi$ satisfying equation 1, the threshold
is not met and the probability of the dynamic photon interacting is 0. In this
case the dynamic photon will propagate through the static field unaffected.
However, if $n(\epsilon,\theta,\phi)\not=0$ there is a finite probability of
the dynamic photon interacting. To obtain this probability requires the mean
free path of the interaction, $\lambda$, which is obtained through the
following [17]

$\begin{split}\frac{1}{\lambda}=\int_{0}^{2\pi}\mathrm{d}\phi\int_{0}^{\pi}\mathrm{d}\theta\int_{0}^{\infty}\mathrm{d}\epsilon\,\,\sigma_{\textrm{BW}}(s)\,n(\phi,\theta,\epsilon)\,(1-\cos{\theta})\end{split}$
(2)

where the factor of $(1-\textrm{cos}(\theta))$ accounts for the relative
velocity between photons and $\sigma_{\mathrm{BW}}$ is the polarisation
averaged, total BW cross-section [18]:

$\sigma_{\textrm{BW}}=\frac{\pi}{2}r_{e}^{2}(1-\beta^{2})\bigg{[}-2\beta(2-\beta^{2})+(3-\beta^{4})\ln\frac{1+\beta}{1-\beta}\bigg{]}$
(3)

where $\beta=\sqrt{(1-4s^{-1})}$ and $r_{e}$ is the classical electron radius.
The probability that a dynamic photon propagates a distance $l$ is

$P(l)=\lambda^{-1}e^{-l/\lambda}.$ (4)

For each dynamic photon, a value of $l$ is sampled from this distribution. If
$l$ is longer than the static photon field, the dynamic photon will propagate
through unaffected. If $l$ is shorter, the dynamic photon will propagate $l$
and annihilate.

When a dynamic photon annihilates it is removed from the simulation. We then
move onto step 2, determining the dynamics of the interaction. First the
properties of the annihilating static photon ($\epsilon,\theta,\phi$) must be
calculated. These are obtained by sampling from the integrand in equation 2.
In the c.m. frame, both the electron and positron receive an energy,
$E_{+,-}$, of $\sqrt{s}/2$. The c.m. frame positron polar scattering angle,
$\theta_{+}$, is obtained by sampling from the BW differential cross-section
[18]

$\frac{d\sigma_{\textrm{BW}}}{d\Omega}=\frac{r_{e}^{2}\beta}{s}\bigg{[}\frac{1+2\beta\mathrm{sin}^{2}\theta-\beta^{4}-\beta^{4}\mathrm{sin}^{4}\theta}{(1-\beta^{2}\mathrm{cos}^{2}\theta)^{2}}\bigg{]}$
(5)

which is plotted in figure 2. The positron azimuthal scattering angle,
$\phi_{+}$, is sampled from the uniform distribution $\mathcal{U}(0,2\pi)$.
Using conservation of momentum, the electron polar and azimuthal scattering
angles are given by $\theta_{-}=\pi-\theta_{+}$ and $\phi_{-}=\phi_{+}-\pi$.
From $E_{+,-}$, $\theta_{+,-}$ and $\phi_{+,-}$ the momentum in the c.m.
frame, $\textbf{P}_{+,-}$, is calculated. Finally a Lorentz boost is performed
into the laboratory frame and an electron and positron are added with the
calculated momentum.

Figure 2: Polarisation averaged BW differential cross-section.

##  III Increasing Efficiency with Gaussian Process Regression

Due to the low photon density, length of field and value of
$\sigma_{\mathrm{BW}}$, the probability that an individual photon will undergo
the linear BW process is small. To analyse the interaction, many dynamic
photons must be simulated. For each dynamic photon $\lambda$ is calculated,
involving a triple integral over the static photon variables (equation 2).
This makes the method outlined in section II computationally expensive.

A common solution to this problem is to replace the direct, expensive
calculation of $\lambda$ with a lookup table. $\lambda$ is a function of three
variables of the dynamic photon, its energy, $E$, and direction, $\theta_{d}$
and $\phi_{d}$, through the rotation to the dynamic photon frame. A lookup
table would involve discretizing $E$, $\theta_{d}$ and $\phi_{d}$ to form a
three dimensional grid. At each grid point $\lambda$ is calculated and values
are obtained at runtime by interpolating from the grid. However, a lookup
table is not appropriate in this case for the following reasons

  1. 1.

$\lambda$ varies rapidly with $E$, $\theta_{d}$ and $\phi_{d}$. Therefore, to
avoid numerical errors, a fine grid resolution is required. The generation of
this lookup table is then computationally expensive.

  2. 2.

To avoid numerical errors through extrapolation, the limits of the table must
be known prior to the simulation runtime. This can be challenging if the
dynamic photons are generated through a multi-step process.

  3. 3.

The generation of a uniform dense lookup table does not represent the
distribution of the dynamic photon parameters, making it inefficient.

We have used an alternative approach where a Gaussian process regression (GPR)
model is trained dynamically and used to quickly calculate $\lambda$ when
required. GPR is a Bayesian method for approximating predictions from an
expensive physical model. Here, we will only discuss our implementation and
not the concepts of GPR. For a full review of GPR refer to Rasmussen and
Williams (2006) [19].

We can consider equation 2 as an expensive function mapping the dynamic photon
parameters $\textbf{x}\equiv(E,\theta_{d},\phi_{d})$ to $\lambda$. The goal is
to replace equation 2 with a cheaper function, $f$. However, $\lambda$ varies
rapidly with x, therefore, it is more numerically accurate to find a function
which returns the natural logarithm of the mean free path estimate

$\lambda_{a}=e^{f(\textbf{x})}$ (6)

where $f(\textbf{x})$ is provided by a GPR. By using a GPR we have assumed
that $f(\textbf{x})$ is Gaussian distributed. Therefore, $\lambda_{a}$ is log-
normal distributed with mean and variance given by

$\begin{split}\mu_{\lambda_{a}}&=\mathrm{exp}\Bigg{(}\mu_{f}-\frac{\sigma^{2}_{f}}{2}\Bigg{)}\\\
\sigma^{2}_{\lambda_{a}}&=\Big{[}\mathrm{exp}(\sigma^{2}_{f})-1\Big{]}\,\mathrm{exp}(2\mu_{f}+\sigma_{f}^{2})\end{split}$
(7)

where $\mu_{f}$ and $\sigma^{2}_{f}$ are the mean and variance of
$f(\textbf{x})$.

To train the GPR a data set,
$\mathcal{D}=\\{\mathbf{x}_{i},\mathrm{log}(\lambda_{i})\\}$, is required.
This is generated by calculating equation 2 for a limited number of events. In
the units system of Geant4, the dimensions of $\mathbf{x}$ vary over vastly
different scales. Therefore, if the regression model is trained directly on
$\mathcal{D}$ it would perform poorly. To solve this problem $\mathbf{x}$ is
normalised using min-max feature scaling.

In this project the open source library libgp [20], was used for the GPR and
the implementation can be broken down into the following three stages:

  1. 1.

Data accumulation stage: Equation 2 is solved in full for $n_{d}$ dynamic
photons and the result saved, generating $\mathcal{D}$.

  2. 2.

Training stage: After simulating $n_{d}$ dynamic photons the GPR is trained on
$\mathcal{D}$ by optimising the model hyperparameters.

  3. 3.

Acceleration stage: For subsequent events, ${\sigma_{\lambda_{a}}}$ is
calculated using equation 7 if ${\sigma_{\lambda_{a}}}<\sigma_{\mathrm{max}}$
where $\sigma_{\mathrm{max}}$ is a user defined limit, $\lambda$ is given by
$\mu_{\lambda_{a}}$. If ${\sigma_{\lambda_{a}}}>\sigma_{\mathrm{max}}$
equation 2 is solved and the result is appended to $\mathcal{D}$. After
another $n_{d}$ points are added to $\mathcal{D}$ the hyperparameters are
again optimised.

As the simulation progresses $\mathcal{D}$ grows. This causes the confidence
of the GPR to increase, so it is used more often than the full calculation of
equation 2. The rate at which dynamic photons are simulated then increases.
The values of $n_{d}$ and $\sigma_{\mathrm{max}}$ can be tuned to optimise the
efficiency for an acceptable level of error.

This GPR scheme overcomes all the shortfalls of a lookup table listed above.
As the GPR is trained dynamically, no expensive calculations prior to runtime
are required. Also, no limits on the dynamic photon parameters are set prior
to runtime. This avoids extrapolation as $\lambda$ is calculated in full when
the dynamic photon parameters fall well outside the data set. Finally, as the
GPR is trained on a sample from the dynamic photon population, it is fully
representative of the dynamic photon parameters distribution.

To demonstrate the increase in computational efficiency a test simulation has
been carried out for both the full calculation of equation 2 and the GPR
scheme. This test simulation involves a gamma ray beam with a divergence of
$10\,\mathrm{mrad}$ and a Gaussian energy spectrum with mean $1\,\mathrm{GeV}$
and standard deviation $0.25\,\mathrm{GeV}$, interacting with a $1\,$cm long,
$300\,$eV isotropic black body radiation field

$n(\epsilon,\theta,\phi)=\frac{\epsilon^{2}}{\pi^{2}}\frac{1}{e^{\epsilon/T}-1}\frac{\sin\theta}{4\pi}.$
(8)

For the GPR scheme the parameters $n_{d}$ and $\sigma_{\mathrm{max}}$ are set
to 100 and $10^{-2}$ respectively. Figure 3 shows the number of dynamic
photons simulated against runtime. The red curve corresponds to the full
calculation of equation 2 whereas the blue curve corresponds to the GPR
scheme. At the start of the simulation figure 3 shows both schemes have the
same rate of dynamic photon calculation. This is during the data accumulation
stage. Just after $500\,$s, the GPR enters the training stage. Beyond this
point the GPR model dynamic photon computation rate increases rapidly. After
$2000\,$s the GPR scheme has calculated 100 times more dynamic photons. When
performing experimental parameter optimisation, the same static photon field
is often used in multiple simulations. In this case the GPR model can be saved
after the end of a simulation and reused in a later one. This avoids the slow
data accumulation stage and provides an even greater increase in efficiency as
shown by the black curve in figure 3.

Figure 3: Number of simulated dynamic photons against runtime for the full
calculation of equation 2 (red), the accelerated GPR scheme (blue) and the
pre-trained GPR scheme (black).

To test if the GPR scheme is accurately emulating the full calculation of
equation 2, simulations can be carried out for both algorithms and the number
of pairs generated, $N_{\mathrm{BW}}$, compared. However, both algorithms
include uncertainty in $N_{\mathrm{BW}}$ due to the finite number of events
simulated. Therefore, an ensemble of simulations should be carried out and the
distributions over $N_{BW}$ compared. The result of this is shown in figure 4
for two different values of $\sigma_{\mathrm{max}}$. When
$\sigma_{\mathrm{max}}=0.01$ there is little difference between the
distributions of the two algorithms, suggesting that the GPR scheme is
accurately emulating 2. However, figure 4 also shows that if
$\sigma_{\mathrm{max}}$ is set too large, this is not the case and the GPR
algorithm introduces errors to the estimate.

Figure 4: Histograms of number of positrons generated per simulation. The red
curve shown the full calculation of equation 2 and the blue and black curves
show calculations performed with the GPR scheme with a $\sigma_{\mathrm{max}}$
of $0.01$ and $1.0$ respectively.

##  IV Module Discussion

Fully modelling the interaction between two photon sources where both are
treated dynamically is an N-body problem with complexity of
$\mathcal{O}(\mathrm{N}^{2})$. This is too computationally expensive to
directly solve for the large number of photons required to statistically
analyse the interaction. Through the use of a tree code, as demonstrated by
Jansen et al. [21], the problem complexity can be reduced to
$\mathcal{O}(\mathrm{N\,log\,N})$. Here, by treating one source as a static
field, the complexity is reduced further to $\mathcal{O}(\mathrm{N})$.
However, the type of photon-photon interactions which can be modelled are also
constrained. By defining the static photon field with
$n(\epsilon,\theta,\phi)$ neither spatial or temporal gradients are accounted
for. It is trivial to redefine the static field also as a function of
position, $n(\mathbf{r},\epsilon,\theta,\phi)$, and account for spatial
gradients by constructing the full field through multiple sub-fields with
different $n(\epsilon,\theta,\phi)$. But it is not possible to relax the
temporal gradient constraint due to the static nature of the Geant4
computational domain. Therefore, the module can only be used when one photon
source is effectively stationary over the full duration of the interaction.

Although it constrains the type of experiment which can be modelled,
developing the module within Geant4 allows us to make use of its extensive
toolkit. It has the capability to model particle tracking, the geometry of
experiments and the response of detectors to energetic particles. Combining
these capabilities with the module discussed here, start-to-end simulations
can be performed. This includes: the generation of the gamma ray beam in the
bremsstrahlung target; the background noise; the tracking of pairs through the
analyser magnets and the particle detector interaction. This can be used for
experimental parameter sensitivity analysis and optimisation. The ability to
fully visualise the setup with Geant4 (see figure 5) also makes experimental
design much quicker [1].

Figure 5: Geant4 visualisation of basic Breit-Wheeler experiment. An electron
(red) produces a high energy gamma ray (green) in a bismuth target. The
particles propagate through an on-axis magnet, diverting seed electron. The
gamma ray propagates into an X-ray field annihilating and forming an
$e^{+}e^{-}$ pair.

##  V Module Demonstration

Here we demonstrate how the module can be used to analyse a BW detection
experiment similar to the Kettle et al. experiment discussed in section I but
with a simplified geometry and representation of the X-ray source. The setup
is shown in figure 5. An electron beam is incident onto a $1\,\mathrm{mm}$
bismuth converter foil. This produces gamma rays through the bremsstrahlung
process. The particles then pass through an on-axis magnet, removing the seed
electrons and BH pairs generated in the converter foil. The gamma rays
propagate into a $1\,$mm long X-ray field, generated by a laser heated
germanium foil, and annihilate.

The X-ray field has spatial gradients, is anisotropic and has a non-blackbody
energy spectrum, making this an ideal test for the module. To model the
spatial gradients, the field is split up into an array of $10\times 10\times
10$ photon volumes, each with a different $n(\epsilon,\theta,\phi)$. We assume
that $n(\epsilon,\theta,\phi)$ is not correlated in energy and angle. This
allows us to separate the energy and angular dependant parts,
$n(\epsilon,\theta,\phi)=f(\epsilon)\Phi(\theta,\phi)$. $f(\epsilon)$ is then
modelled using the atomic physics code Flychk [22] and an example is shown in
figure 6 with a foil temperature of $350\,\mathrm{eV}$. $\Phi(\theta,\phi)$ is
modelled using a Monte-Carlo X-ray photon tracking code. Here, photons are
launched from an emitting plane and tracked through the photon volume array.
By launching many photons, $\Phi(\theta,\phi)$ can be estimated.

Figure 6: Radiation spectrum emitted from a static, solid density germanium
plasma with an electron temperature of 350 eV. This was generated using the
atomic physics code Flychk [22].

Figure 7 shows the angle and energy probability distribution for BW positrons
produced with a 1 GeV electron beam and a burn-through foil electron
temperature of 350 eV. From figure 7 we can see the positron beam is well
collimated with a divergence of $\sim 10\,\mathrm{mrad}$. The divergence of
the positrons decreases with increasing energy. The particle beaming is a
result of the gamma ray carrying almost all the momentum of the interaction.
This is a key advantage of the Pike et al. detection scheme. If the positron
beam had a large divergence, which is the case for two symmetric photon
sources, to prevent a low detection efficiency the interaction region would
have to be surrounded by detectors. A large number of high divergence
background BH pairs would also be detected. However, with a highly beamed
source, an analyser chicane can be used to select positrons within a specific
divergence and energy band and transport them to an area with low noise to be
detected.

Figure 7: Positron angle and energy distribution. Figure 8: Breit-Wheeler pair
yield against X-ray foil temperature and electron beam energy.

With the increased computational efficiency, as a result of implementing a
GPR, the module can be used to perform scans over experimental parameters.
This is demonstrated in figure 8 where we can see how the yield of BW pairs is
affected by the burn-through foil temperature and the LWFA electron beam
energy. Each simulation consisted of $10^{9}$ electrons (160 pC) and the
electron beam energy was varied from 500 MeV to 3000 MeV and foil temperature
from $50\,$eV to $500\,$eV. The Gemini laser facility is capable of producing
$2\,$GeV electron beams with a charge of $<100\,$pC. Therefore, these results
suggest that with a foil temperature of $>150\,$eV, it is possible to generate
$>1$ BW pair per shot. Utilising the relatively high repetition rate of the
Gemini laser facility, over multiple shots a detectable number of BW pairs
above background could be produced.

##  VI Conclusion

As the power of laser facilities has increased in recent years, it is now
theoretically possible to produce a detectable number of BW pairs in the
laboratory. However, the expected signal-to-noise ratio of such an experiment
is low, making detailed numerical modelling vital. Here, we have presented the
development of a new linear BW module for Geant4. This will allow signal-to-
noise ratio calculations to be performed within a single framework. We have
shown how a Gaussian process regression can be used to greatly increase the
rate at which photons are simulated without a loss in the accuracy of the
module. We have used this increase in efficiency to perform a parameter scan
which suggests that $>1$ BW pair per shot can be produced using the Gemini
laser facility.

The module presented here can be readily extended to include other particle-
photon interactions such as Compton scattering ($e^{-}\gamma\rightarrow
e^{-}\gamma$) and photon-photon scattering
($\gamma\gamma\rightarrow\gamma\gamma$). Photon-photon scattering from two
real photons is another QED process with astrophysical relevance that has
never been directly observed in the laboratory. The difficulties of observing
the process arise due to the relatively small cross-section ($\sim\alpha^{2}$
times smaller than $\sigma_{\mathrm{BW}}$) and the incoming and outgoing
particles all being photons, making background subtraction challenging. The
scattering of quasi-real photons in the collision of heavy ion beams has
recently been observed by the ATLAS detector at the Large Hadron Collider
[23]. With improvements to the experimental scheme from section I, it may also
be possible to study photon-photon scattering of two real photons, which could
be modelled using this module.

## References

  * Kettle  _et al._ [2021] B. Kettle, D. Hollatz, E. Gerstmayr, G. Samarin, A. Alejo, S. Astbury, C. Baird, S. Bohlen, M. Campbell, C. Colgan, _et al._ , New Journal of Physics 23, 115006 (2021).
  * Breit and Wheeler [1934] G. Breit and J. A. Wheeler, Phys. Rev. 46, 1087 (1934).
  * Piran [2005] T. Piran, Rev. Mod. Phys. 76, 1143 (2005).
  * Bonometto and Rees [1971] S. Bonometto and M. Rees, Monthly Notices of the Royal Astronomical Society 152, 21 (1971).
  * Nikishov [1961] A. Nikishov, Zhur. Eksptl’. i Teoret. Fiz. 46, (1961).
  * Gould and Schréder [1967a] R. J. Gould and G. P. Schréder, Phys. Rev. 155, 1408 (1967a).
  * Aleksić  _et al._ [2011] J. Aleksić, L. A. Antonelli, P. Antoranz, M. Backes, J. A. Barrio, D. Bastieri, J. Becerra González, W. Bednarek, A. Berdyugin, K. Berger, _et al._ , The Astrophysical Journal Letters 530, A4 (2011), arXiv:1101.2522 [astro-ph.CO] .
  * Adam  _et al._ [2021] J. Adam, L. Adamczyk, J. Adams, J. Adkins, G. Agakishiev, M. Aggarwal, Z. Ahammed, I. Alekseev, D. Anderson, A. Aparin, _et al._ , Physical review letters 127, 052302 (2021).
  * Pike  _et al._ [2014] O. J. Pike, F. Mackenroth, E. G. Hill, and S. J. Rose, Nature Photonics 8, 434 (2014).
  * Ribeyre  _et al._ [2016] X. Ribeyre, E. D’Humières, O. Jansen, S. Jequier, V. T. Tikhonchuk, and M. Lobet, Physical Review E 93, 013201 (2016).
  * Drebot  _et al._ [2017] I. Drebot, D. Micieli, E. Milotti, V. Petrillo, E. Tassi, and L. Serafini, Physical Review Accelerators and Beams 20, 043402 (2017).
  * Hoarty  _et al._ [2010] D. Hoarty, S. James, C. Brown, B. Williams, H. Chung, J. Harris, L. Upcraft, B. Crowley, C. Smith, and R. Lee, High Energy Density Physics 6, 105 (2010).
  * Bethe and Heitler [1934] H. Bethe and W. Heitler, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 146, 83 (1934).
  * Agostinelli  _et al._ [2003] S. Agostinelli, J. Allison, K. a. Amako, J. Apostolakis, H. Araujo, P. Arce, M. Asai, D. Axen, S. Banerjee, G. Barrand, _et al._ , Nuclear instruments and methods in physics research section A: Accelerators, Spectrometers, Detectors and Associated Equipment 506, 250 (2003).
  * Battistoni  _et al._ [2015] G. Battistoni, T. Boehlen, F. Cerutti, P. W. Chin, L. S. Esposito, A. Fassò, A. Ferrari, A. Lechner, A. Empl, A. Mairani, A. Mereghetti, P. G. Ortega, J. Ranft, S. Roesler, P. R. Sala, V. Vlachoudis, and G. Smirnov, Annals of Nuclear Energy 82, 10 (2015).
  * Briesmeister  _et al._ [2000] J. F. Briesmeister _et al._ , Version 4C, LA-13709-M, Los Alamos National Laboratory 2 (2000).
  * Gould and Schréder [1967b] R. J. Gould and G. P. Schréder, Phys. Rev. 155, 1404 (1967b).
  * Jauch  _et al._ [2009] J. M. Jauch, F. Rohrlich, and J. C. Polkinghorne, Physics Today 9, 34 (2009).
  * Rasmussen and Williams [2006] C. Rasmussen and C. Williams,  _Gaussian Processes for Machine Learning_ , Adaptive Computation and Machine Learning (MIT Press, Cambridge, MA, USA, 2006) p. 248.
  * Blum [2019] M. Blum, libgp, https://github.com/mblum/libgp (2019).
  * Jansen  _et al._ [2018] O. Jansen, E. D’Humières, X. Ribeyre, S. Jequier, and V. T. Tikhonchuk, Journal of Computational Physics 355, 582 (2018).
  * Chung  _et al._ [2005] H.-K. Chung, M. Chen, W. Morgan, Y. Ralchenko, and R. Lee, High Energy Density Physics 1, 3 (2005).
  * ATLAS collaboration [2017] ATLAS collaboration, Nature Physics 13, 852 (2017).