# NeuroLGP-SM: Scalable Surrogate-Assisted Neuroevolution for Deep Neural
Networks

Fergal Stapleton∗ ∗joint-first authorship Naturally Inspired Computation
Research Group,
Department of Computer Science, Hamilton Institute
Maynooth University, Ireland
<EMAIL_ADDRESS>Edgar Galván∗ Naturally Inspired Computation
Research Group,
Department of Computer Science, Hamilton Institute
Maynooth University, Lero, Ireland
<EMAIL_ADDRESS>
###### Abstract

Evolutionary Algorithms (EAs) play a crucial role in the architectural
configuration and training of Artificial Deep Neural Networks (DNNs), a
process known as neuroevolution. However, neuroevolution is hindered by its
inherent computational expense, requiring multiple generations, a large
population, and numerous epochs. The most computationally intensive aspect
lies in evaluating the fitness function of a single candidate solution. To
address this challenge, we employ Surrogate-assisted EAs (SAEAs). While a few
SAEAs approaches have been proposed in neuroevolution, none have been applied
to truly large DNNs due to issues like intractable information usage. In this
work, drawing inspiration from Genetic Programming semantics, we use
phenotypic distance vectors, outputted from DNNs, alongside Kriging Partial
Least Squares (KPLS), an approach that is effective in handling these large
vectors, making them suitable for search. Our proposed approach, named Neuro-
Linear Genetic Programming surrogate model (NeuroLGP-SM), efficiently and
accurately estimates DNN fitness without the need for complete evaluations.
NeuroLGP-SM demonstrates competitive or superior results compared to 12 other
methods, including NeuroLGP without SM, convolutional neural networks, support
vector machines, and autoencoders. Additionally, it is worth noting that
NeuroLGP-SM is 25% more energy-efficient than its NeuroLGP counterpart. This
efficiency advantage adds to the overall appeal of our proposed NeuroLGP-SM in
optimising the configuration of large DNNs.

###### Index Terms:

Neuroevolution, Linear Genetic Programming, Surrogate-assisted Evolutionary
Algorithms

##  I Introduction

Evolutionary Algorithms (EAs) [10] have proven to be effective in both the
crafting of architectures and hyperparameter optimisation of Deep Neural
Networks (DDNs) [25]. This application is commonly known as neuroevolution
[11], and have been applied to numerous problem domains, such as autonomous
vehicles [14, 39] and face recognition [9, 8]. The pursuit of optimal DNN
architectures has led to various methodologies, including EAs [11],
reinforcement learning [44], etc. However, a significant challenge persists
across these methods: the substantial computational resources required to
identify high-performing networks.

The rise of GPU-accelerated hardware has helped alleviate some of this
computational cost, however, a significant proportion of research in DNNs is
based on incremental improvements on DNN algorithms for benchmark problems,
where there is a significant correlation between network complexity for
incremental gains in terms of additional performance. In fact, when looking at
very large models of hundreds of billions of parameters, it can cost millions
of dollars for a single iteration [26]. This energy consumption is further
compounded when considering population-based neuroevolutionary techniques
which require many networks to be trained and evaluated in order to find
suitable architectures.

One way to address this significant issue is with the use of surrogate-
assisted evolutionary algorithms (SAEAs). SAEAs can be used to estimate the
fitness of DNNs without the need to fully train each network. In particular,
surrogate modelling strategies that employ Bayesian optimisation have shown
much promise [26]. However, a major challenge remains in how best to deal with
the surrogate representation. For instance, using genotype information to
build surrogates often requires complex encoding strategies [41], and in some
instances, constructing adequate distance metrics to compare different network
topologies is not feasible [40]. Using phenotypic information on the other
hand has shown some promise [19, 40], but a challenge remains in scaling to
deeper and more complex networks which inherently requires a high-dimensional
representation [38].

In this work, we analyse a novel population-based Neurovolutionary technique,
referred to as NeuroLGP, and its surrogate model variant NeuroLGP-SM [36].
Using a robust model management strategy, we use phenotypic distance vectors
to estimate the performance of partially trained DNNs. These vectors are
comparatively large for the optimisation problem at hand [38] and, as such, we
use an approach that is designed for handling high-dimensional data, known as
Kriging Partial Least Squares (KPLS). This approach allows for a novel and
scalable surrogate-assisted technique that is skilfully adept at handling
neuroevolution of DNNs and to the the best of our knowledge, this method of
surrogate-assisted neuroevolution has not been studied before.

The aim of this study is to apply Surrogate-assisted Evolutionary Algorithms
(SAEA) in neuroevolution, using Kriging Partial Least Squares (KPLS) on
phenotypic distance vectors inspired by Genetic Programming Semantics. Our key
contributions are outlined as follows, 1) Efficient Fitness Estimation: We
demonstrate the accurate estimation of DNN fitness values without full
evaluations through our proposed approach, NeuroLGP-SM, employing KPLS on
phenotypic distance vectors. 2) Performance Metrics: We employ three well-
defined metrics to assess the predictive capabilities of NeuroLGP-SM in terms
of model fitness. These results align with the competitive or superior
performance of our NeuroLGP-SM approach compared to 12 popular techniques,
including convolutional neural networks, autoencoders, and support vector
machines, across four challenging classification tasks. 3) Energy Consumption
Analysis: We provide a reliable formula to gauge the energy consumption of our
approaches, revealing that NeuroLGP-SM is 25% more efficient compared to its
counterpart that does not use surrogate models and 4) Encoding for Analysis:
Through clever encoding, we allow easy access to analyse the internal
structures of the architectures, enabling us to conduct an in-depth analysis
of the networks discovered by our proposed approaches.

##  II Related Work

###  II-A Neuroevolution using Surrogate Models

Santos et al. [33] proposed a novel approach that makes use of semantics in
neuroevolution. They do so by using Geometric Semantic Genetic Programming, in
conjunction with a neuroevolutionary approach called Semantic Learning
Machines (SLM) [16], a form of neuroevolution technique. Furthermore, Hagg et
al. [19], made the connection between semantic distances in GP and phenotypic
distances within the context of surrogate-assisted evolutionary algorithms for
neuroevolution. Similarly, Stork et al. [40] extended CGP to use a surrogate-
assisted neuroevolution approach that makes use of phenotyic distance vectors.
A limitation of Stork’s work is the scalability of using Kriging on high-
dimensional data. For instance, in the recent work by Stapleton and Galván,
highlighted that traditional approaches such as the Kriging approach suffer
with high-dimensionality [38] and may not be suitable for DNN architectures.

Freeze-Thaw Bayesian Optimisation (FBO), proposed by Swersky et al. [43], uses
Bayesian optimisation to determine whether a particular neural network that
has been partially trained should be fully evaluated. This approach is novel
in that the system stops training (or freezing) of less promising networks,
instead spending valuable resources on the most promising networks. Of note,
is the fact that the FBO approach relies on the phenotypic behaviour of DNNs.
It is important to note the majority of works rely on genotypic information
when building surrogate models for neuroevolution [41, 17]. For reference, a
major work in this regard is that of Sun et al. [41] referred to as the End-
to-End Performance Predictor. This approach uses an offline surrogate model
based on random forests to accelerate learning of CNN architectures. The CNN
architecture is encoded such that it maps to a numerical decision variable,
which is passed to the random forest-based surrogate-model. Not only does this
approach not require large amounts of training time but, also alleviates the
requirement of a smooth learning curve [43].

##  III Background

###  III-A Surrogate assisted models: Kriging and Kriging Partial Least
Squares

Surrogate modelling has many use cases but in the context of this work, the
aim is to effectively estimate the fitness values for candidate solutions
while simultaneously reducing the run time of the evolutionary process. To
this end, not only must the surrogate model be well-posed, but also, the
evolutionary process must interact with a robust surrogate model management
strategy [20]. The surrogate model differs from the parent model that instead
of training directly on the data the surrogate model is trained on the design
space, where the aim is to identify and further explore regions of this design
space that will produce preferable parameters. As such, interpolation-based
approaches may be used to simulate the unknown regions of the parameter space,
such as Gaussian processes, commonly referred to as Kriging. Another benefit
to the Kriging approach is that it allows for estimates of the uncertainty of
predictions.

Kriging is an interpolation-based technique that assumes spatial correlation
exists between known data points, based on the distance, and variation between
these points. We aim to use observations $\\{y(x_{1}),y(x_{2}),...y(x_{n})\\}$
to help estimate an unknown function value $\hat{y}(x^{*})$ for the unknown
data point $x^{*}$, where $n$ is the number of individuals in the training
data. A kernel function $K(\cdot)$ is used to express the spatial correlation
between two samples $x_{i}$ and $x^{{}^{\prime}}_{i}$ as shown in Eq. 1,

$k(x,x^{\prime})=\prod_{i=1}^{m}exp(-\theta(\mathcal{D}(x_{i},x^{{}^{\prime}}_{i})))$
(1)

where the $\theta$ parameter controls the rate at which the correlation decays
to zero between $x_{i}$ and $x^{{}^{\prime}}_{i}$ and
$\mathcal{D}(x_{i},x^{{}^{\prime}}_{i})=(x_{i}-x^{{}^{\prime}}_{i})^{2}$. The
$\theta$ parameter is determined using the Maximum Likelihood Estimator (MLE).
A significant drawback, however, is that for large dimensions $m$, the
computational cost increases dramatically since the MLE algorithm calculates
the inverse of the correlation matrix multiple times.

Kriging Partial Least Squares (KPLS) helps to alleviate this limitation by
reducing the number of parameters calculated [4]. It does so by using partial
least squares which project the high-dimensional data into a lower dimension
using principal components. Eq. 2 details the KPLS kernel,

$K(x,x^{\prime})=\prod_{k=1}^{h}\prod_{i=1}^{m}exp(-\theta_{k}(w^{(k)}_{i}x_{i}-w^{(k)}_{i}x^{\prime}_{i})^{2})$
(2)

where $w$ are rotated principal directions which maximise the covariance and
are a measure of how important each principal component is. Typically, the
number of principal components $h$ is much less than the number of dimensions
$m$ and, as such, allows for the substantial improvement in computational cost
associated with the KPLS approach. See [4] for details.

###  III-B Phenotypic Distance

The inspiration for using phenotypic distance is by works such as Stork et al.
[40] and Galván et al. [13, 15, 12, 37], whose previous work in GP centred on
semantics, in particular, the use of semantic distance metrics. In traditional
GP, semantics can be understood as the behaviour of program given a finite set
of inputs. A definition by Pawlak et al. [31] states

Def 1. _The semantics $s(p)$ of a program $p$ is the vector of values from the
output set $O$ obtained by computing $p$ on all inputs from the input set
$I$_. This is expressed in Eq. 3,

$s(p)=[p(in_{1}),p(in_{2}),...,p(in_{l})]$ (3)

where $l=|I|$ is the size of the input set and where the inputs $in$ from can
be some subset of the dataset. In terms of neuroevolution, we can define a
solution sample $x$ as having the semantic or phenotypic behaviour of the
$i^{th}$ program such that $x_{i}=s(p_{i})$, where the semantics $s(p)$ is the
vector of values from the output. From Eq. 1 we can then define our phenotypic
distance $\mathcal{D}$ as,

$\mathcal{D}(x_{i},x_{j})=\mathcal{D}(s(p_{i}),s(p_{j})))$ (4)

As such the distance metric is dependent on the outputs of each network, where
the $x_{i}$ is a flattened vector containing the outputs from the final layer
for all data instances. Thes phenotypic distance vector length, as such, is
the number of images of the test dataset times the number of classes. Notably,
this approach can be extended to any deep learning model architecture that
uses a vectorised output (i.e., transformers), however, the limitations and
scalability of other deep learning models with this method are yet to be
investigated.

##  IV Methodology

###  IV-A NeuroLGP

The original Linear Genetic Programming (LGP) encoding is based on the concept
of using registers, which are units computer memory storage for manipulating
data while executing instructions, in a low-level programming context. The
content of these registers are altered using instruction operations. There are
three main components to instruction; an operand which performs a specific
function on one or more registers which store the result in a destination
register. In the case of 2-register instruction encoding the operand operates
on a single instruction and for a 3-register instruction encoding operates on
two instructions. The left of Fig. 1 highlights an example of LGP written in C
code, where r[i] denotes the $i^{th}$ register. This example contains both
effective and non-effective lines of code, where the non-effective lines of
code are commented out and subsequently not compiled. Each line of code is
executed imperatively. The register r[0] is a specially designated register
for the final output of the program. With the NeuroLGP approach, instead of
using registers to store small amounts of data, we instead use the idea of
registers as pointers to much larger amounts of data. As such, the registers
instead control the flow of the initial and intermediary data from each
outputted layer of our evolvable DNN. The right of Fig. 1, demonstrates how
the expected representation would look, where the general form of this code is
similar to the figure on the left.

⬇

void gp(r)

double r[8];

{ …

r[0] = r[5] + 71;

// r[7] = r[0] - 59;

if (r[1] > 0)

if (r[5] > 2)

r[4] = r[2] * r[1];

// r[2] = r[5] + r[4];

r[6] = r[4] * 13;

r[0] = r[6];

}

⬇

def neuro_gp(…)

{ …

r[0] := Conv(r[1])

// r[4] := MaxPool(r[3])

// r[7] := Conv(r[4])

r[5] := MaxPool(r[2])

r[6] := Conv(r[5])

r[3] := AvgPool(r[6])

// r[5] := MaxPool(r[8])

r[0] := Dense(r[3])

}

Figure 1: left: LGP example in C language. Based on example from [5]. right:
NeuroLGP psuedocode for python.

Traditionally, in LGP, registers can be of different types (input registers,
calculation registers and constant registers). In this work, each register is
initialised with the inputs of various layers within the neural network and so
act as just input registers (though it may be useful to have special registers
specifically for handling skip connections and multiple branches). While this
may seem wasteful in terms of memory, multiple registers are only used in the
genotypic representation and when calculating the fitness, the number of
effective registers can be reduced as part of a repair process. Likewise, the
non-effective code can also be removed at the same time. As such, the
phenotypic representation of the network will not only be much simpler than
the genotypic representation, but also much smaller. Fig.2 demonstrates the
genotype-to-phenotype mapping.

Figure 2: Diagram of the NeuroLGP genotype-to-phenotype mapping. The
pseudocode for the set of instructions (left-hand side) can be represented as
the genotype with effective and non-effective code (top) and produces the
resulting phenotype (bottom right-hand side) as a specific neural network
architecture. Note that the non-effective coding is not present in the
phenotype.

###  IV-B NeuroLGP with Surrogate Model (NeuroLGP-SM)

To identify individuals requiring full evaluation, we use the Expected
Improvement (EI) criteria [21]. EI guides the selection of candidate solutions
for evaluation by estimating the improvement over the current best solution.
This enables us to prioritise solutions from areas in the search space
expected to exhibit the most significant improvement. The calculation of EI is
shown in Eq. 5,

EI $\displaystyle=\left\\{\begin{aligned}
&(\hat{f}(x)-f(x^{*}))\Phi(Z)+\sigma(x)\phi(Z)&&\quad\text{if}\ \sigma(x)>0\\\
&0&&\quad\text{if}\ \sigma(x)=0\end{aligned}\right.$ (5) $\displaystyle Z$
$\displaystyle=\frac{\hat{f}(x)-f(x^{*})}{\sigma(x)}$

where $\hat{f}(x)$ is the model’s predicted performance of the surrogate for
the phenotypic distance vector $x$, where $f(x^{*})$ is the best-known value
of the objective function so far (in this case maximum) and $\Phi$ and $\phi$
are the cumulative distribution function (CDF) and probability density
function (PDF) of the standard normal distribution, respectively. Fig. 3
summarises the surrogate model management strategy. The annotations are: (i)
Split: first, the population is split with a 40/60 split for individuals to be
fully evaluated vs. the partially trained individuals, respectively, (ii)
Estimate fitness: the fitness is estimated using the KPLS approach as informed
by the previous generation, (iii) Add data: the phenotypic distance vector for
the fully evaluated portion of the population is added to the surrogate
training data (note: the phenotypic distance vector is taken before the full
fitness evaluation), (iv) Extract data and train: the phenotypic distance
vectors are collectively used to train the KPLS approach as detailed in
Section III-A, and (v) Calculate EI: the EI criteria is calculated for the
incoming population of individuals.

Figure 3: Model management strategy in terms of fitness evaluation, model
management and training the surrogate model.

##  V Experimental setup

The Breast Cancer Histopathological Image Classification (BreakHis) [35] is a
binary classification dataset consisting of microscopic images containing
2,480 benign and 5,429 malignant tumours, where images were obtained as four
different magnifications ($\times$40, $\times$100, $\times$200 and
$\times$400) and are split into four individual datasets for this work. Each
image consists of 64x64 pixels (down-scaled resolution from 700x460 pixels)
and 3-channel RGB with 8-bit depth in each channel. The split for training,
validation, test and test2 is approximately (63.5%, 12.5%, 12.5%, 12.5%),
where the first three splits conform to a ratio of (70:15:15), conforming to
the general training/test split as seen in other approaches (Table I). The
first three splits are used purely for training and evaluating networks and
the last split is retained for an unbiased analysis of the evolutionary
process. The dataset is imbalanced and the synthetic minority oversampling
technique [6] is used to up-sample the minority class.

We employ two approaches to investigate the validity and quality of our
proposed surrogate-assisted neuroevolution model. These are: (i) The NeuroLGP
approach, which uses a novel encoding to evolve the structure of our networks,
training to the full number of epochs as outlined in Section IV-A, and (ii)
The NeuroLGP-SM approach where we integrate the surrogate modelling strategy,
as outlined in Section IV-B, into the NeuroLGP approach. Experiments were
conducted using Kay supercomputer provided by the Irish Centre for High-End
Computing (ICHEC). Experiments were run in parallel, with each run assigned to
a single Nvidia Tesla V100 GPU with 16GB Ram. Additionally, each run has
access to a 20-core 2.4 GHz Intel Xeon Gold 6148 (Skylake) CPU processor which
is used during training of the surrogate model portion of the NeuroLGP-SM
approach. Overall, the NeuroLGP approach took $\sim$28 GPU days and the
NeuroLGP-SM approach $\sim$21 GPU days, for 8 runs across the 4 datasets.

##  VI Results

###  VI-A Comparison of Models

Table I summarises approaches that have used the BreakHis data set from the
last few years, including our proposed approaches. Most of these works have
been taken from a comprehensive 2020 review paper from Benhammou et al. [3].
For fair comparison, we have selected works that also do not use transfer
learning as these represent manually crafted networks that have zero pre-
training. Furthermore, many of these approaches use specialised feature
extraction and pre-processing steps for histopathological data. As such, our
goal is not necessarily to outperform previous state-of-the-art works, but
rather to show that our neuroevolutionary technique can achieve similar
results to hand-crafted architectures, even without specialised knowledge of
the problem domain.

Looking at the last four columns from the right, we list the accuracies for
the various magnifications. The last four rows correspond to our proposed
approaches: NeuroLGP and its surrogate-assisted variant (NeuroLGP-SM). We can
see that our results are in good accordance with the other works, in some
cases outperforming other approaches. Notably, the results yielded by our two
approaches on the $\times$200 magnification are very competitive in terms of
accuracy.

If we turn our attention to specifically comparing the NeuroLGP-SM and
NeuroLGP models, we can see that each has very similar performances. For
instance, the mean values (third and fourth last lines from bottom), the
NeuroLGP-SM and NeuroLGP models are typically within 0.2$\%$ of each other and
for the best (last two lines) the NeuroLGP-SM and NeuroLGP models are within
1-2$\%$.

TABLE I: Accuracy results for approaches using BreakHis dataset ($\times 40$,
$\times 100$, $\times 200$ and $\times 400$). Results from this work are
highlighted in boldface and are presented in the last four rows. WSI: wholes
slide image, CNN: convolutional neural network, SVM: support vector machine,
AE: auto-encoder, DBN: deep-belief network.

Work | Preprocessing | Patch/slide | Model | Training/Test |  $\times$40 |  $\times$100 |  $\times$200 |  $\times$400
---|---|---|---|---|---|---|---|---
Gupta [18] | None | WSI | Ensemble | 70% / 30% | 88.9 | 88.9 | 88.9 | 88.9
Sharma [34] | Mixed | WSI | Ensemble | Not specified | 81.7$\pm$2.8 | 81.2$\pm$2.7 | 80.7$\pm$3.4 | 81.5$\pm$3.1
Nahid [27] | None | WSI | CNN | Not specified | 94.4 | 95.93 | 97.19 | 96
Nahid [28] | K-Mean Clustering | WSI | CNN | Not specified | 85 | 90 | 90 | 90
Karthiga [22] | K-Mean Clustering | WSI | SVM | Not specified | 93.3 | 93.3 | 93.3 | 93.3
Pratiher [32] | Gray Scale, Data Aug | WSI | AE | Not specified | 96.8 | 98.1 | 98.2 | 97.6
Badejo [2] | None | WSI | SVM | 70% / 30% | 91.1 | 90.7 | 86.2 | 84.3
Nahid [29] | Contrast Enhancement | WSI | DBN | 70% / 30% | 88.7 | 89.06 | 88.84 | 87.67
Das [7] | Resize (370x230) | Patches (224x224) | CNN | 80% / 20% | 89.52 | 85.3 | 88.6 | 88.4
Kumar [23] | Stain Normalisation | Patches (64x64, 32x32) | CNN | 70% / 30% | 82$\pm$2.8 | 86.2$\pm$4.6 | 84.6$\pm$3.0 | 84$\pm$4.0
Aswathy [1] | Mixed | WSI | SVM | 90% / 10% | 89.1 | 89.1 | 89.1 | 89.1
NeuroLGP (Mean) | Resize (64x64), Data Aug | WSI | CNN | See Section V | 89.7$\pm$2.1 | 0.872$\pm$3.0 | 92.6$\pm$2.1 | 90.8$\pm$1.0
NeuroLGP-SM (Mean) | Resize (64x64), Data Aug | WSI | CNN | See Section V | 89.8$\pm$1.4 | 0.873$\pm$1.8 | 92.8$\pm$2.4 | 91.0$\pm$0.7
NeuroLGP (Best) | Resize (64x64), Data Aug | WSI | CNN | See Section V | 93 | 91.6 | 96 | 92.5
NeuroLGP-SM (Best) | Resize (64x64), Data Aug | WSI | CNN | See Section V | 91.3 | 90.3 | 97 | 92.5

###  VI-B Analysis of the NeuroLGP-SM Model

Three metrics are used to determine how well the surrogate model performs in
predicting the fitness of our partially trained models: (i) the Mean Squared
Error (MSE) gives a measure of how accurate our model is in terms of
predicting fitness where values close to 0 are preferable, (ii) Kendall’s Tau
is used to measure the correlation between the predicted fitness and the
actual fitness [42], where values close to 0 indicates no correlation, -1 a
perfect negative correlation, +1 a perfect positive correlation, and (iii) the
R2 score is a measure of how close the predicted value is to its true value
and ranges from -$\infty$ to 1 and explains how much variability is in the
prediction model, where values closer to 1 are preferable.

Table II summarises the effectiveness of the NeuroLGP-SM model for each of the
datasets. The low MSE values, as seen in the first row, show a relatively low
error between the predicted and actual fitness across the four dataset. The
Kendall’s Tau values range from 0.5647 to 0.6791, indicating a strong positive
correlation between the actual and predicted fitness. We can see that for
$\times$40 we have R2 of 0.5026 indicating a moderate level of fit being
captured by the surrogate model while for $\times$100, $\times$200 and
$\times$400, we have R2 values of 0.6665, 0.7079 and 0.7786. There is also a
general trend showing the R2 score increasing with higher magnifications. This
may be a result of more noise/artefacts present in the lower magnification
images, making it more difficult for the surrogate model to predict based on
the network outputs but further studies would be required to confirm this.
Overall, the quality of fit and performance of the surrogate models based on
the three metrics is very strong for each dataset.

Next, we plot the average MSE per generation to get an idea of how the
predictive capability of the NeuroLGP-SM approach changes over time as seen in
Fig. 4. MSE values that either decrease or remain stable are preferable, as
increasing MSE values would indicate our surrogate model is losing its
predictive capability as new individuals are introduced. We can see that after
an initial drop in average MSE from the first couple of generations (x-axis),
while there are some fluctuations in the average MSE, in general values remain
relatively stable, demonstrating the robustness of our NeuroLGP-SM.

Figure 4: Avg. MSE between predicted vs. actual fitness over 15 generations for the NeuroLGP-SM approach. The relative stability shows the robustness of the NeuroLGP-SM approach. TABLE II: Average MSE, Kendall’s Tau and R2. Metric | Magnification Size
---|---
$\times$40 |  $\times$100 |  $\times$200 |  $\times$400
MSE | 0.0037 | 0.0017 | 0.0014 | 0.0009
Kendall’s Tau | 0.6019 | 0.6791 | 0.6225 | 0.5647
R2 | 0.5026 | 0.6665 | 0.7079 | 0.7786

Fig. 5 details the quality of fit in terms of how well the predicted
accuracies relate to the actual accuracies. In each figure, values closer to
the red diagonal line, represent individuals that have a better quality of
fit. We can see that the higher the accuracy the closer our predicted values
are to the red line. This is a by-product of the surrogate model management
strategy we employ. In essence, we are less concerned with the accuracies of
poorer-performing models, instead, we have taken a more greedy approach
favouring models that are better performing.

Figure 5(a) informs us why the $\times$40 magnification slightly under-
performed in terms of the metrics as shown in Table II. We can see that for
the predicted values, on the y-axis, there are quite several individuals
predicted to have 80% accuracy when the actual accuracies, for these
individuals range from 30-80% in terms of actual accuracy, as shown on the
x-axis. A deeper dive revealed that this was a result of one particularly bad-
performing run and was not present in the other 7 runs.

(a) $\times$40

(b) $\times$100

(c) $\times$200

(d) $\times$400

Figure 5:  Predicted vs. actual accuracies (black points), for $\times$40,
$\times$100, $\times$200, and $\times$400 datasets, shown in (a) – (d),
respectively, across 8 independent runs for the surrogate-assisted approach
(NeuroLGP-SM). The red line denotes where the accuracy for both predicted and
actual are the same, where points closer to this line are preferential.

###  VI-C Analysis of the Energy Consumption

We perform a comparison of the estimated energy consumption of the NeuroLGP-SM
and the NeuroLGP models using the Green-algorithm by Lannelongue et al. [24].
Eq. 6 is used to get this estimate,

$\displaystyle E=r_{t}*(P_{c}*U+P_{m})*PUE*PSF$ (6)

where $E$ is the energy consumed in KW/h, $r_{t}$ is the runtime, $P_{c}$ is
the power draw for the computing cores and depends on the CPU/GPU model, $U$
is the usage factor, $P_{m}$ is the power draw for the memory. $PUE$ is the
Power Usage Efficiency and is a measure of how much additional energy is
needed to operate the data centre. $PSF$ is the Pragmatic Scaling Factor and
can be used to estimate the performance when multiple runs are taken into
consideration. A conservative $PUE$ value of 1.67 was used, which represents
the worldwide average power usage $PUE$ across all data centres [24]. For
simplicity, $PSF$ was kept as 1 as we are interested in the average runtime
rather than the cumulative runtime.

The real usage factor was minimal for both the NeuroLGP-SM and the NeuroLGP
models, even during training for the NeuroLGP-SM model and additionally, the
CPU energy consumption is negligible compared to that of the GPU, so we have
reported just the energy saved by the GPU, overall 2.79 KW/h are saved using
the NeuroLGP-SM model. Considering the total energy saved across all 4
magnifications for 8 runs each, we save approximately 89.28 KW/h: 25% less
energy is consumed using the NeuroLGP-SM approach.

###  VI-D Analysis of Genotype

Figs. 6 (a) and (b) show the proportions of various genes for two of the four
magnifications $\times$100 and $\times$400, respectively. Due to page
constraints, we have shown only these two datasets, but our analysis for the
same for the other two, unless stated otherwise in the text. Initial
proportions (blue) represent the proportions of specific genes at
initialisation (first generation). The NeuroLGP-SM (green) and NeuroLGP
(orange) model are also represented and denote the proportions in the final
generation across all 8 runs. The comparison we make will aim to show how the
proportions the NeuroLGP-SM and NeuroLGP for specific layers change from
initialisation (i.e., we compare green and orange proportions against blue).

For the initial proportions, genes were proportioned based on their functional
grouping. The four groups are dropout, batch normalisation, pooling, and
convolutional layers. As such, each group makes up a quarter of the initial
population and are subsequently divided again by specific genes. For instance,
as there are 6 convolutional layer genes we divide $\frac{0.25}{6}$ to get the
proportion for each convolutional layer.

Turning our attention to the NeuroLGP-SM and NeuroLGP model, we can see
general trends in how the various gene groupings change by the final
generation shown in green and orange for the NeuroLGP-SM and NeuroLGP,
respectively. For instance, the dropout layer proportions sizes tend to
decrease significantly by the final generation in the NeuroLGP-SM and NeuroLGP
models. Similarly, the proportion of batch normalisation layers tends to
increase. The pooling layers are more specific to the dataset. For instance,
for $\times$100 there is an increase specifically for the max pooling layer
(also found in $\times$200 but not shown). On the other hand, for $\times$400
there is a general decrease, albeit the two models differ on the specific
pooling layer they decrease. Again, for the convolutional layers there is a
more specific trend for particular datasets. While there are some slight
increases and decreases in proportion size for $\times$40, $\times$100 and
$\times$200, for $\times$400 there are notable increases in proportion sizes
for $\times$400 for convolutional layers using a 3$\times$3 filter (Figure
6(b), CONV_32_3x3, CONV_64_3x3 and CONV_128_3x3).

Some of these trends are unsurprising, we would expect that the proportion of
dropout layers would be less, for instance, as CNN architectures are not as
prone to overfitting [30]. It would seem that the decrease in pooling layers
in the $\times$400 is being compensated by an increase in convolutional
layers, in other words, some of the dimension reduction is being handled by
the convolutional layers rather than the pooling. Further research could
investigate how these proportions change relative to different resolution
sizes. Interestingly, while in many cases the change in proportions are
consistent for both the NeuroLGP-SM and NeuroLGP, in some instances the
proportion sizes differ significantly, for instance, CONV_64_5x5 as seen in
$\times$400\. This is likely a result of solutions converging at different
local optima, however, further analysis of the fitness landscape would be
beneficial in gaining more insight.

###  VI-E Limitations of our analysis

While a single run is the norm in terms of neuroevolution, we endeavoured to
perform as many runs as possible, to add better confidence to our analysis.
While we ran individual runs in parallel, both the NeuroLGP-SM and NeuroLGP
approaches, for 8 runs each and for the 4 datasets, the total runtime amassed
$\sim$50 GPU days. To conduct a full statistical analysis we would require
another +50 GPU days of experimentation. Future work, could look into using
different problem domains and different Deep learning models.

(a) Breakhis $\times$100 dataset

(b) Breakhis $\times$400 dataset

Figure 6: Proportion of evolved network layers for initialisation (blue) and
final generations for NueroLGP-SM (green) and NeuroLGP (orange) for
$\times$100, and $\times$400 datasets

##  VII Conclusion

In this work, we demonstrate a novel surrogate-assisted neurovolutionary
approach, named Neuro-Linear Genetic Programming LGP surrogate model
(NeuroLGP-SM). This approach makes use of Kriging Partial Least Squares (KPLS)
to estimate the fitness of partially trained Deep Neural Networks (DNNs) using
phenotypic distance vectors in a high-dimensional context. This scalable
approach to surrogate-assisted neuroevolution is adept at finding robust and
efficient architectures. We demonstrate that our approach is competitive or
superior to other state-of-the-art handcrafted networks. Furthermore, we
demonstrate that the proposed NeuroLGP-SM approach is 25% more efficient in
energy consumption, compared to a NeuroLGP variant that does not use
surrogacy. Additionally, due to the unique encoding properties of our NeuroLGP
approach, we can easily analyse the internal structures of the architectures,
giving greater insight into the favorable components of the discovered
architectures.

## Acknowledgements

This publication has emanated from research conducted with the financial
support of Science Foundation Ireland under Grant number 18/CRT/6049. The
author has applied a CC BY public copyright licence to any Author Accepted
Manuscript version arising from this submission.

## References

  * [1] M. Aswathy and M. Jagannath.  An svm approach towards breast cancer classification from h&e-stained histopathology images based on integrated features.  Medical & biological engineering & computing, 59(9):1773–1783, 2021.
  * [2] J. A. Badejo, E. Adetiba, A. Akinrinmade, and M. B. Akanle.  Medical image classification with hand-designed or machine-designed texture descriptors: a performance evaluation.  In Bioinformatics and Biomedical Engineering: 6th International Work-Conference, IWBBIO 2018, Granada, Spain, April 25–27, 2018, Proceedings, Part II 6, pages 266–275. Springer, 2018.
  * [3] Y. Benhammou, B. Achchab, F. Herrera, and S. Tabik.  Breakhis based breast cancer automatic diagnosis using deep learning: Taxonomy, survey and insights.  Neurocomputing, 375:9–24, 2020.
  * [4] M. A. Bouhlel, N. Bartoli, A. Otsmane, and J. Morlier.  Improving kriging surrogates of high-dimensional design models by partial least squares dimension reduction.  Structural and Multidisciplinary Optimization, 53:935–952, 2016\.
  * [5] M. Brameier and W. Banzhaf.  Linear genetic programming, volume 1.  Springer, 2007.
  * [6] N. V. Chawla, K. W. Bowyer, L. O. Hall, and W. P. Kegelmeyer.  Smote: synthetic minority over-sampling technique.  Journal of artificial intelligence research, 16:321–357, 2002.
  * [7] K. Das, S. Conjeti, A. G. Roy, J. Chatterjee, and D. Sheet.  Multiple instance learning of deep convolutional neural networks for breast histopathology whole slide classification.  In 2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018), pages 578–581. IEEE, 2018.
  * [8] S. Deng, Z. Lv, E. Galván, and Y. Sun.  Evolutionary neural architecture search for facial expression recognition.  IEEE Transactions on Emerging Topics in Computational Intelligence, 7(5):1405–1419, 2023.
  * [9] S. Deng, Y. Sun, and E. Galvan.  Neural architecture search using genetic algorithm for facial expression recognition.  In Proceedings of the Genetic and Evolutionary Computation Conference Companion, GECCO ’22, page 423–426, New York, NY, USA, 2022. Association for Computing Machinery.
  * [10] A. E. Eiben and J. E. Smith.  Introduction to Evolutionary Computing.  Springer Verlag, 2003.
  * [11] E. Galván and P. Mooney.  Neuroevolution in deep neural networks: Current trends and future challenges.  IEEE Transactions on Artificial Intelligence, 2:476–493, 2021.
  * [12] E. Galván and M. Schoenauer.  Promoting semantic diversity in multi-objective genetic programming.  In A. Auger and T. Stützle, editors, Proceedings of the Genetic and Evolutionary Computation Conference, GECCO 2019, Prague, Czech Republic, July 13-17, 2019, pages 1021–1029. ACM, 2019.
  * [13] E. Galván and F. Stapleton.  Semantic-based distance approaches in multi-objective genetic programming.  In 2020 IEEE Symposium Series on Computational Intelligence (SSCI), pages 149–156. IEEE, 2020.
  * [14] E. Galván and F. Stapleton.  Evolutionary multi-objective optimisation in neurotrajectory prediction.  Applied Soft Computing, 146:110693, 2023.
  * [15] E. Galván, L. Trujillo, and F. Stapleton.  Semantics in multi-objective genetic programming.  Applied Soft Computing, 115:108143, 2022.
  * [16] I. Gonçalves, S. Silva, and C. M. Fonseca.  Semantic learning machine: a feedforward neural network construction algorithm inspired by geometric semantic genetic programming.  In Progress in Artificial Intelligence: 17th Portuguese Conference on Artificial Intelligence, EPIA 2015, Coimbra, Portugal, September 8-11, 2015. Proceedings 17, pages 280–285. Springer, 2015.
  * [17] B. Greenwood and T. McDonnell.  Surrogate-assisted neuroevolution.  In Proceedings of the Genetic and Evolutionary Computation Conference, pages 1048–1056, 2022.
  * [18] V. Gupta and A. Bhavsar.  An integrated multi-scale model for breast cancer histopathological image classification with joint colour-texture features.  In Computer Analysis of Images and Patterns: 17th International Conference, CAIP 2017, Ystad, Sweden, August 22-24, 2017, Proceedings, Part II 17, pages 354–366. Springer, 2017.
  * [19] A. Hagg, M. Zaefferer, J. Stork, and A. Gaier.  Prediction of neural network performance by phenotypic modeling.  In Proceedings of the Genetic and Evolutionary Computation Conference Companion, pages 1576–1582, 2019.
  * [20] Y. Jin.  Surrogate-assisted evolutionary computation: Recent advances and future challenges.  Swarm and Evolutionary Computation, 1(2):61–70, 2011.
  * [21] D. R. Jones, M. Schonlau, and W. J. Welch.  Efficient global optimization of expensive black-box functions.  Journal of Global optimization, 13:455–492, 1998.
  * [22] R. Karthiga and K. Narasimhan.  Automated diagnosis of breast cancer using wavelet based entropy features.  In 2018 Second international conference on electronics, communication and aerospace technology (ICECA), pages 274–279. IEEE, 2018.
  * [23] K. Kumar and A. C. S. Rao.  Breast cancer classification of image using convolutional neural network.  In 2018 4th International Conference on Recent Advances in Information Technology (RAIT), pages 1–6. IEEE, 2018.
  * [24] L. Lannelongue, J. Grealey, and M. Inouye.  Green algorithms: quantifying the carbon footprint of computation.  Advanced science, 8(12):2100707, 2021.
  * [25] Y. LeCun, Y. Bengio, and G. Hinton.  Deep learning.  Nature, 521(7553):436–444, 2015.
  * [26] G. Menghani.  Efficient deep learning: A survey on making deep learning models smaller, faster, and better.  ACM Computing Surveys, 55(12):1–37, 2023.
  * [27] A.-A. Nahid and Y. Kong.  Histopathological breast-image classification using local and frequency domains by convolutional neural network.  Information, 9(1):19, 2018.
  * [28] A.-A. Nahid, M. A. Mehrabi, Y. Kong, et al.  Histopathological breast cancer image classification by deep neural network techniques guided by local clustering.  BioMed research international, 2018, 2018.
  * [29] A.-A. Nahid, A. Mikaelian, and Y. Kong.  Histopathological breast-image classification with restricted boltzmann machine along with backpropagation.  Biomedical Research, 29(10):2068–2077, 2018.
  * [30] S. Park and N. Kwak.  Analysis on the dropout effect in convolutional neural networks.  In Computer Vision–ACCV 2016: 13th Asian Conference on Computer Vision, Taipei, Taiwan, November 20-24, 2016, Revised Selected Papers, Part II 13, pages 189–204. Springer, 2017.
  * [31] T. P. Pawlak, B. Wieloch, and K. Krawiec.  Semantic backpropagation for designing search operators in genetic programming.  IEEE Transactions on Evolutionary Computation, 19(3):326–340, June 2015.
  * [32] S. Pratiher and S. Chattoraj.  Diving deep onto discriminative ensemble of histological hashing & class-specific manifold learning for multi-class breast carcinoma taxonomy.  In ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 1025–1029. IEEE, 2019.
  * [33] F. J. Santos, I. Gonçalves, and M. Castelli.  Neuroevolution with box mutation: An adaptive and modular framework for evolving deep neural networks.  Applied Soft Computing, page 110767, 2023.
  * [34] M. Sharma, R. Singh, and M. Bhattacharya.  Classification of breast tumors as benign and malignant using textural feature descriptor.  In 2017 Ieee international conference on bioinformatics and biomedicine (bibm), pages 1110–1113. IEEE, 2017.
  * [35] F. A. Spanhol, L. S. Oliveira, C. Petitjean, and L. Heutte.  Breast cancer histopathological image classification using convolutional neural networks.  In 2016 international joint conference on neural networks (IJCNN), pages 2560–2567. IEEE, 2016.
  * [36] F. Stapleton, B. Cody-Kenny, and E. Galván.  Neurolgp-sm: A surrogate-assisted neuroevolution approach using linear genetic programming.  In International Conference on Optimization and Learning (OLA), 2024\.
  * [37] F. Stapleton and E. Galván.  Semantic neighborhood ordering in multi-objective genetic programming based on decomposition.  In 2021 IEEE Congress on Evolutionary Computation (CEC), pages 580–587. IEEE, 2021.
  * [38] F. Stapleton and E. Galvan.  Initial steps towards tackling high-dimensional surrogate modeling for neuroevolution using kriging partial least squares.  In Proceedings of the Companion Conference on Genetic and Evolutionary Computation, GECCO ’23 Companion, page 83–84, New York, NY, USA, 2023. Association for Computing Machinery.
  * [39] F. Stapleton, E. Galván, G. Sistu, and S. Yogamani.  Neuroevolutionary multi-objective approaches to trajectory prediction in autonomous vehicles.  In Proceedings of the Genetic and Evolutionary Computation Conference Companion, GECCO ’22, page 675–678, New York, NY, USA, 2022. Association for Computing Machinery.
  * [40] J. Stork, M. Zaefferer, and T. Bartz-Beielstein.  Improving neuroevolution efficiency by surrogate model-based optimization with phenotypic distance kernels.  In International Conference on the Applications of Evolutionary Computation (Part of EvoStar), pages 504–519. Springer, 2019.
  * [41] Y. Sun, H. Wang, B. Xue, Y. Jin, G. G. Yen, and M. Zhang.  Surrogate-assisted evolutionary deep learning using an end-to-end random forest-based performance predictor.  IEEE Transactions on Evolutionary Computation, 24(2):350–364, 2020\.
  * [42] Y. Sun, B. Xue, M. Zhang, and G. G. Yen.  Evolving deep convolutional neural networks for image classification.  IEEE Transactions on Evolutionary Computation, 24(2):394–407, 2019\.
  * [43] K. Swersky, J. Snoek, and R. P. Adams.  Freeze-thaw bayesian optimization.  arXiv preprint arXiv:1406.3896, 2014.
  * [44] B. Zoph and Q. V. Le.  Neural architecture search with reinforcement learning.  CoRR, abs/1611.01578, 2016.


# Thermal Pair Production from Photon-Photon Collision:
Breit-Wheeler Process at Finite Temperature

D. S. Cabral<EMAIL_ADDRESS>Instituto de Física, Universidade
Federal de Mato Grosso,
78060-900, Cuiabá, Mato Grosso, Brazil  A. F. Santos
<EMAIL_ADDRESS>Instituto de Física, Universidade Federal de
Mato Grosso,
78060-900, Cuiabá, Mato Grosso, Brazil  R. Bufalo<EMAIL_ADDRESS>Departamento de Física, Universidade Federal de Lavras,
Caixa Postal 3037, 37200-900 Lavras, MG, Brazil

###### Abstract

In this paper we examine the pair production through the Breit-Wheeler process
$\displaystyle\gamma\leavevmode\nobreak\ \gamma\to e^{+}e^{-}$ in a thermal
background. We compute the thermal contribution to the Breit-Wheeler
differential cross section within the thermofield dynamics formalism. We
evaluate in details the cross section for this process, which possess a
surprisingly simple expression valid for any temperature $\displaystyle\beta$,
from which we discuss some physically relevant aspects. We also consider the
high temperature regime of the cross section in order to have a better
understanding about its thermal behavior.

##  I Introduction

Although quantum electrodynamics (QED) can be regarded as one of the main
triumphs in physics with unprecedented significance, by describing the atomic
spectra and particle scattering, as well as being the cornerstone of the
standard model of particle physics, it is still an active area of research. We
can cite as one of the main ongoing interest the study of QED in intense
electromagnetic fields, or simply strong-field QED [1, 2, 3, 4], which had its
beginning with the study of Schwinger-Sauter pair-production from an external
field, which consists in the spontaneous creation of charged particle pairs by
electromagnetic fields [5, 6]. On the experimental side of studying QED at
extreme conditions, the continued development in laser technology and maximum
focused intensities led to the growing impact of laser-plasma interactions in
many facilities: SLAC [7], the European XFEL [8], DESY [9] and CERN [10].

In addition to the Schwinger’s effect, one can also produce charged particle
pairs from the collision of photons, for instance: Positron Production (from
linear and nonlinear trident process) [11] and Breit-Wheeler process [12]. The
Breit-Wheeler pair production is a physical process in which a positron-
electron pair is created from the collision of two photons,
$\displaystyle\gamma\leavevmode\nobreak\ \gamma\to e^{+}e^{-}$ (which is
basically the inverse process of electron-positron annihilation), which has
recently been experimentally observed [13] 111The detection of this Breit-
Wheeler process, with quasi-real photons, occurred in the peripheral
collisions of high Z nuclei at a conventional accelerator [13].. This process
can be extended to the multiphoton (nonlinear) Breit-Wheeler process [14, 15],
when the collision of two photons is replaced by a high-energetic probe photon
(a single $\displaystyle\gamma$ with several eV laser photons) decaying into
pairs propagating through an electromagnetic field (for example, a laser
pulse). Both linear and nonlinear Breit-Wheeler processes have been
extensively studied in many scenarios: extreme astrophysical environments,
like those surrounding pulsars and gamma-ray bursts, but also in heavy-ion
colliders, and high power laser facilities, among others [16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34].

One interesting aspect to be analyzed, on theoretical grounds, in the
production of charged pairs is the role played by thermal effects and possible
implications. Although the Schwinger effect at finite temperature has been
widely studied [36, 35, 37, 38, 39, 40, 41, 42, 43], as well as the Breit-
Wheeler process [41, 42] (see ref. [41] for a detailed list of works) 222These
should not be thought of as a competing thermal processes, actually, applying
the Optical Theorem, one can see that a unitarity cut of the dominant Feynman
diagram for thermal Schwinger pair production gives the Breit-Wheeler process
[42]., pair production in an electric field at finite temperature (with
thermal fermions and/or thermal photons) has a long history with many
conflicting papers [2]. Realizing that the thermal pair production is a
somewhat controversial topic, in which many papers (results) disagree with
each other, we approach the Breit-Wheeler process in terms of the thermofield
dynamics, with the expectation to shed some new light and clarify (some of)
the controversies.

The close relation of the technical aspects of the thermofield dynamics [44,
45, 46, 47] with the case of $\displaystyle T=0$ field theory, considering the
analysis and calculations of scattering amplitudes and decay processes, makes
this formalism very interesting to describe finite temperature field theory,
presenting some practical advantages in regard to Matsubara’s imaginary time
formalism and Schwinger-Keldysh closed time path method [48, 49, 50]. The main
feature of the thermofield dynamics method is that the Fock space is doubled,
in which a new set of operators is introduced, designed tilde operators,
acting on the second Fock space, tilde space. This can be understood as the
following: physically speaking, the second Fock space is interpreted as a heat
bath that ensures the dynamical system to stay in equilibrium. These two Fock
spaces are connected by means of a Bogoliubov transformation, which consists
of a rotation between the two spaces, original and tilde.

Since finite temperature effects have profound implications in the study of
high energy physics, we examine the pair production of Breit-Wheeler process
in the presence of a heat bath within the thermofield dynamics. We start Sec.
II by reviewing the main aspects of the Breit-Wheeler process and establishing
the framework of thermofield dynamics to compute the respective amplitude
related with the differential cross section. In Sec. III we evaluate
explicitly the squared transition amplitude of the Breit-Wheeler process in
the center of mass frame, where the kinematics of the process takes a simple
form. Moreover, since the angular dependence is such that we can obtain the
respective cross section valid for any value of $\displaystyle\beta$, from
which we discuss some physically relevant aspects, including comments about
the presence of off-shell contributions. In Sec. IV we summarize the results,
and present our final remarks.

##  II Breit-Wheeler Amplitude at finite temperature

The main purpose of this section is to calculate the differential cross
section for the Breit-Wheeler (BW) process at finite temperature within the
thermofield dynamics (TFD) [44, 45, 46, 47]. The BW process consists of
creating an electron-positron pair through the collision of two photons, i.e.
$\displaystyle\gamma+\gamma=e^{+}+e^{-}$. The Feynman diagrams that describe
this process are depicted in Figure 1.

Figure 1: The Feynman diagrams contributing to the pair creation in the BW
process.

Since we wish to compute the thermal cross section of the BW process, our
starting point is

$\frac{d\sigma}{d\Omega}=\frac{1}{64\pi^{2}(4\omega^{2})}\langle|\widehat{\mathcal{M}}_{fi}|^{2}\rangle$
(1)

in which the temperature effects are solely contained in the transition
amplitude $\displaystyle\mathcal{M}(\beta)$. This thermal amplitude is
formally defined as

$\mathcal{\mathcal{M}}\left(\beta\right)=\left\langle
f,\beta\left|\hat{S}\right|i,\beta\right\rangle$ (2)

where the $\displaystyle\hat{S}$-matrix is given

$\hat{S}=\sum_{n=0}^{\infty}\frac{\left(-i\right)^{n}}{n!}\int
dx_{1}dx_{2}\cdots
dx_{n}\tau\left[\hat{\mathcal{L}}_{i}\left(x_{1}\right)\hat{\mathcal{L}}_{\rm
int}\left(x_{2}\right)\cdots\hat{\mathcal{L}}_{\rm
int}\left(x_{n}\right)\right],$ (3)

with $\displaystyle\tau$ being the time ordering operator and

$\hat{\mathcal{L}}_{\rm int}\left(x\right)=\mathcal{L}_{\rm
int}\left(x_{n}\right)-\tilde{\mathcal{L}}_{\rm int}\left(x_{n}\right)$ (4)

describes the Lagrangian interaction part in the doubled notation of the TFD
formalism. In this context, the final and initial thermal states for the
process $\displaystyle\gamma(k)+\gamma(k^{\prime})=e^{+}(p^{\prime})+e^{-}(p)$
are given as

$\displaystyle\displaystyle\ket{i}$
$\displaystyle\displaystyle=\ket{\gamma,\gamma^{\prime}}=d^{\dagger}_{\beta}(k)d^{\dagger}_{\beta}(k^{\prime})\ket{0(\beta)},$
(5) $\displaystyle\displaystyle\ket{f}$
$\displaystyle\displaystyle=\ket{e^{+},e^{-}}=b^{\dagger}_{\beta}(p)a^{\dagger}_{\beta}(p^{\prime})\ket{0(\beta)},$
(6)

with $\displaystyle d^{\dagger}_{\beta}$, $\displaystyle b^{\dagger}_{\beta}$
and $\displaystyle a^{\dagger}_{\beta}$ being the thermal creation operators
for the photons, electrons and positrons, respectively.

In order to calculate the thermal transition amplitude for the BW process, the
second-order term of the scattering matrix (3) is considered,

$\widehat{S}^{(2)}=\int d^{4}xd^{4}y\bra{f}\tau\left[\hat{\mathcal{L}}_{\rm
int}(x)\hat{\mathcal{L}}_{\rm int}(y)\right]\ket{i},$ (7)

and the thermal transition amplitude for the BW process is defined as

$\displaystyle\displaystyle\widehat{\mathcal{M}}_{fi}=\bra{f,\beta}\widehat{S}^{(2)}\ket{i,\beta}=$
$\displaystyle\displaystyle-\frac{e^{2}}{2}\int
d^{4}xd^{4}y\left\\{\bra{f}\tau\left[\mathcal{L}_{\rm int}(x)\mathcal{L}_{\rm
int}(y)\right]\ket{i}+\bra{f}\tau\left[\widetilde{\mathcal{L}}_{\rm
int}(x)\widetilde{\mathcal{L}}_{\rm int}(y)\right]\ket{i}\right.$ (8)
$\displaystyle\displaystyle-$
$\displaystyle\displaystyle\left.\bra{f}\tau\left[\mathcal{L}_{\rm
int}(x)\widetilde{\mathcal{L}}_{\rm
int}(y)\right]\ket{i}-\bra{f}\tau\left[\widetilde{\mathcal{L}}_{\rm
int}(x)\mathcal{L}_{\rm int}(y)\right]\ket{i}\right\\}.$ (9)

In the quantum electrodynamics, the interaction part of the Lagrangian is
given as

$\displaystyle\displaystyle\mathcal{L}_{\rm int}(x)$
$\displaystyle\displaystyle=e\bar{\psi}(x)\gamma^{\mu}\psi(x)A_{\mu}(x),$ (10)
$\displaystyle\displaystyle\tilde{\mathcal{L}}_{\rm int}(x)$
$\displaystyle\displaystyle=e\tilde{\bar{\psi}}(x)\gamma^{\mu}\tilde{\psi}(x)\tilde{A}_{\mu}(x),$
(11)

where $\displaystyle e$ is the electron charge, $\displaystyle\gamma^{\mu}$
are the Dirac matrices, $\displaystyle\psi(x)$ is the fermion field and
$\displaystyle A_{\mu}(x)$ is the photon field, while the tilde fields are the
TFD dual (conjugated) fields.

In order to evaluate explicitly the matrix elements in (9), we should also
consider the solution for the fermionic field (and its tilde dual)

$\psi(x)=\int
d^{3}pN_{p}\sum_{s}\left[b(p,s)u(p,s)e^{-ipx}+a^{\dagger}(p,s)\bar{u}(p,s)e^{ipx}\right],$
(12)

where $\displaystyle N_{p}$ is the normalization constant and $\displaystyle
px=\xi t-\vec{p}\cdot\vec{x}$, with the dispersion relation
$\displaystyle\xi=\sqrt{\vec{p}^{2}+m^{2}}$ and metric signature
$\displaystyle\eta_{\mu\nu}={\rm diag}(+1,-1,-1,-1)$. In the TFD formalism it
is worth to mention the (tilde) conjugation operations for the fermions
$\displaystyle\tilde{u}(p,s)=u^{\dagger}(p,s)$. Moreover, for the photon field

$A_{\mu}(x)=\int
d^{3}kM_{k}\sum_{\lambda}\left[d(k,\lambda)\epsilon_{\mu}(k,\lambda)e^{-ikx}+d^{\dagger}(k,\lambda)\epsilon^{\dagger}_{\mu}(k,\lambda)e^{ikx}\right]$
(13)

with $\displaystyle\epsilon_{\mu}(k,\lambda)$ being the polarization vector of
the physical polarization states. With these elements we are now in a position
to compute the amplitude (9). Let us start with the first term

$\displaystyle\displaystyle\bra{f}\tau\left[\mathcal{L}_{\rm
int}(x)\mathcal{L}_{\rm int}(y)\right]\ket{i}=\leavevmode\nobreak\ $
$\displaystyle\displaystyle
2\bar{u}(p^{\prime},s^{\prime})\gamma^{\mu}S(x-y)\gamma^{\nu}v(p,s)\varepsilon_{\mu}(k,\lambda)\varepsilon_{\nu}(k^{\prime},\lambda^{\prime})$
$\displaystyle\displaystyle\times\left[e^{i(p-k)x+i(p^{\prime}-k^{\prime})y}+e^{i(p-k^{\prime})x+i(p^{\prime}-k)y}\right]$
$\displaystyle\displaystyle\times\bra{f}\tau\left[a^{\dagger}(p^{\prime},s^{\prime})b^{\dagger}(p,s)d(k,\lambda)d(k^{\prime},\lambda)\right]\ket{i}.$
(14)

One should pay attention to the element matrix present in (14) because it
involves both $\displaystyle T=0$ and thermal creation (annihilation)
operators (recall that the in and out states are thermal states, see eqs. (5)
and (6)). In order to evaluate this element matrix we should use the relation
among these two sets of operators, which formally is given in terms of the
following Bogoliubov transformation [47]

$a(p)=Ua_{\beta}(p)+V\widetilde{a}^{\dagger}_{\beta}(p);\quad\quad\quad
a^{\dagger}(p)=Ua^{\dagger}_{\beta}(p)+V\widetilde{a}_{\beta}(p),$ (15)

for the case of fermions (analogous relations hold for $\displaystyle d(p,s)$)
and

$d(p)=U^{\prime}d_{\beta}(p)+V^{\prime}\widetilde{d}^{\dagger}_{\beta}(p);\quad\quad\quad
d^{\dagger}(p)=U^{\prime}d^{\dagger}_{\beta}(p)+V^{\prime}\widetilde{d}_{\beta}(p)$
(16)

for bosons. The set of functions $\displaystyle(U,V)$ and
$\displaystyle(U^{\prime},V^{\prime})$ are related with the Fermi-Dirac
$\displaystyle n_{F}(E)$ and Bose-Einstein $\displaystyle n_{B}(E)$
distributions, respectively, defined as the following

$\displaystyle\displaystyle U$
$\displaystyle\displaystyle=\frac{1}{\sqrt{1+e^{-\beta E}}}=\sqrt{e^{\beta
E}n_{F}\left(E\right)},\quad V=\frac{1}{\sqrt{1+e^{\beta
E}}}=\sqrt{n_{F}\left(E\right)},$ (17) $\displaystyle\displaystyle U^{\prime}$
$\displaystyle\displaystyle=\frac{1}{\sqrt{1-e^{-\beta E}}}=\sqrt{e^{\beta
E}n_{B}\left(E\right)},\quad V^{\prime}=\frac{e^{-\beta
E/2}}{\sqrt{1-e^{-\beta E}}}=\sqrt{n_{B}\left(E\right)}.$ (18)

Thus, with help of these informations we cast the element matrix of Eq. (14)
as

$\displaystyle\displaystyle\bra{f}\tau\left[a^{\dagger}(p^{\prime},s^{\prime})b^{\dagger}(p,s)d(k,\lambda)d(k^{\prime},\lambda)\right]\ket{i}=U^{2}(U^{\prime})^{2}\bra{f}\tau\left[a_{\beta}^{\dagger}(p^{\prime},s^{\prime})b_{\beta}^{\dagger}(p,s)d_{\beta}(k,\lambda)d_{\beta}(k^{\prime},\lambda)\right]\ket{i},$
(19)

in which all operators are now thermal and this element can be evaluated.

From this analysis we can compute the remaining matrix elements of Eq. (9),
i.e. those proportional to
$\displaystyle\left\langle\widetilde{\mathcal{L}}\widetilde{\mathcal{L}}\right\rangle$,
$\displaystyle\left\langle\widetilde{\mathcal{L}}{\mathcal{L}}\right\rangle$
and
$\displaystyle\left\langle{\mathcal{L}}\widetilde{\mathcal{L}}\right\rangle$,
and it follows that they contribute with the following thermal functions
$\displaystyle V^{2}(V^{\prime})^{2}$, $\displaystyle-
VU(V^{\prime})(U^{\prime})$ and $\displaystyle-UV(U^{\prime})(V^{\prime})$,
respectively. Hence, taking all these results into consideration, the
transition amplitude (9) becomes in momentum space

$\displaystyle\displaystyle\widehat{\mathcal{M}}_{fi}=$
$\displaystyle\displaystyle-
ie^{2}F(\beta)\left(2\pi\right)^{4}\Big{[}\bar{u}(p^{\prime},s^{\prime})\gamma^{\mu}S_{\beta}\left(k-p\right)\gamma^{\nu}v(p,s)\varepsilon_{\mu}(k,\lambda)\varepsilon_{\nu}(k^{\prime},\lambda^{\prime})$
(20)
$\displaystyle\displaystyle+\bar{u}(p^{\prime},s^{\prime})\gamma^{\mu}S_{\beta}\left(k^{\prime}-p\right)\gamma^{\nu}v(p,s)\varepsilon_{\mu}(k,\lambda)\varepsilon_{\nu}(k^{\prime},\lambda^{\prime})\Big{]}$
(21)
$\displaystyle\displaystyle\times\bar{u}(p^{\prime},s^{\prime})\gamma^{\mu}S_{\beta}(q)\gamma^{\nu}v(p,s)\varepsilon_{\mu}(k,\lambda)\varepsilon_{\nu}(k^{\prime},\lambda^{\prime})$
(22)
$\displaystyle\displaystyle\times\bra{f}\tau\left[a_{\beta}^{\dagger}(p^{\prime},s^{\prime})b_{\beta}^{\dagger}(p,s)d_{\beta}(k,\lambda)d_{\beta}(k^{\prime},\lambda)\right]\ket{i},$
(23)

where thermal distributions are contained in the factor

$F(\beta)=(UU^{\prime}+VV^{\prime})^{2}$ (24)

which can be written with help of (17) and (18) in terms of the thermal
distributions

$F(\beta)=\frac{1}{2}\left[\tanh\left(\frac{\beta\omega}{2}\right)+2\sqrt{\frac{1}{e^{\beta\omega}-1}}\sqrt{\frac{1}{\cosh(\beta\omega)+1}}\sqrt{\coth\left(\frac{\beta\omega}{2}\right)+1}+\coth\left(\frac{\beta\omega}{2}\right)\right].$
(25)

Moreover, we have that $\displaystyle S_{\beta}(q)$ is the fermion propagator
at finite temperature defined as

$S_{\beta}(q)=S^{(0)}(q)+S^{(\beta)}(q),$ (26)

in which $\displaystyle S^{(0)}(q)=\frac{\not{q}+m}{q^{2}-m^{2}}$ is the usual
$\displaystyle T=0$ propagator, while $\displaystyle S^{(\beta)}(q)$ is the
finite temperature defined as [47]

$\displaystyle\displaystyle S^{(\beta)}(q)=\frac{\pi i}{(e^{\beta
q^{0}}+1)}\left[\frac{(\gamma^{0}\xi-\gamma^{i}q_{i}+m)}{2\xi}\Delta_{1}\delta(q^{0}-\xi)+\frac{(\gamma^{0}\xi+\gamma^{i}q_{i}-m)}{2\xi}\Delta_{2}\delta(q^{0}+\xi)\right],$
(27)

with the definitions

$\Delta_{1}=\begin{pmatrix}1&e^{\beta q_{0}/2}\\\ e^{\beta
q_{0}/2}&-1\end{pmatrix};\quad\quad\Delta_{2}=\begin{pmatrix}-1&e^{\beta
q_{0}/2}\\\ e^{\beta q_{0}/2}&1\end{pmatrix}.$ (28)

Finally, taking these results into account, we are able to compute the squared
modulus of the transition amplitude

$\displaystyle\displaystyle\langle|\widehat{\mathcal{M}}_{fi}|^{2}\rangle=F^{2}(\beta)\left\\{\langle|\widehat{\mathcal{M}}_{a}|^{2}\rangle_{\beta}+\langle|\widehat{\mathcal{M}}_{b}|^{2}\rangle_{\beta}+\left\langle
2\real\left[\widehat{\mathcal{M}}_{a}^{\dagger}\widehat{\mathcal{M}}_{b}\right]\right\rangle_{\beta}\right\\}$
(29)

with

$\displaystyle\displaystyle\langle|\widehat{\mathcal{M}}_{a}|^{2}\rangle_{\beta}=$
$\displaystyle\displaystyle\frac{e^{4}}{4}\eta_{\mu\rho}\eta_{\nu\alpha}\Tr{(\not{p}^{\prime}+m)\gamma^{\nu}S_{\beta}(k-p)\gamma^{\mu}(\not{p}-m)\gamma^{\rho}S_{\beta}^{\dagger}(k-p)\gamma^{\alpha}};$
(30)
$\displaystyle\displaystyle\langle|\widehat{\mathcal{M}}_{b}|^{2}\rangle_{\beta}=$
$\displaystyle\displaystyle\frac{e^{4}}{4}\eta_{\nu\rho}\eta_{\mu\alpha}\Tr{(\not{p}^{\prime}+m)\gamma^{\mu}S_{\beta}(k^{\prime}-p)\gamma^{\nu}(\not{p}-m)\gamma^{\rho}S_{\beta}^{\dagger}(k^{\prime}-p)\gamma^{\alpha}};$
(31) $\displaystyle\displaystyle\left\langle
2\real\left[\widehat{\mathcal{M}}_{a}^{\dagger}\widehat{\mathcal{M}}_{b}\right]\right\rangle_{\beta}=$
$\displaystyle\displaystyle\frac{e^{4}}{2}\eta_{\mu\rho}\eta_{\nu\alpha}\Re{\Tr\left[(\not{p}-m)\gamma^{\rho}S_{\beta}(k-p)\gamma^{\alpha}(\not{p}^{\prime}+m)\gamma^{\mu}S_{\beta}^{\dagger}(k^{\prime}-p)\gamma^{\nu}\right]}.$
(32)

where we have averaged over the polarization and spin states, with use of the
completeness relations,

$\sum_{s}u(p,s)\bar{u}(p,s)={\not{p}+m},\quad\quad\sum_{s}v(p,s)\bar{v}(p,s)={\not{p}-m}$
(33)

and

$\sum_{\lambda}\epsilon_{\mu}(k,\lambda)\epsilon^{\dagger}_{\rho}(k,\lambda)=\eta_{\mu\rho}.$
(34)

Our next step to compute the cross section for the Breit-Wheeler process is to
compute explicitly (29).

##  III Cross section at finite temperature

In order to evaluate explicitly the squared transition amplitude (29) we
consider the center of mass frame, so that the kinematical variables are
defined as

$\displaystyle\displaystyle p=$
$\displaystyle\displaystyle\,(\omega,A\sin{\theta},0,A\cos{\theta}),\qquad
p^{\prime}=(\omega,-A\sin{\theta},0,-A\cos{\theta})$ (35)
$\displaystyle\displaystyle k=$
$\displaystyle\displaystyle\,(\omega,0,0,\omega),\qquad
k^{\prime}=(\omega,0,0,-\omega),$ (36)

where $\displaystyle A=\sqrt{\omega^{2}-m^{2}}$, in this case the Mandelstam
variables become

$s=2m^{2}-t-u,\qquad t=m^{2}-2\omega(\omega+A\cos{\theta}),\qquad
u=m^{2}-2\omega(\omega-A\cos{\theta}).$ (37)

Within this framework we can now evaluate the trace in Eqs. (30), (31) and
(32), so that the transition amplitude (29) is cast as

$\langle|\widehat{\mathcal{M}}_{fi}|^{2}\rangle=e^{4}F^{2}(\beta)\left[\mathcal{M}^{(0)}+\frac{1}{(1+e^{2\beta\omega})^{2}}\mathcal{M}^{(\beta)}\right]$
(38)

in which the thermal factor $\displaystyle F(\beta)$ is given by (25), and the
$\displaystyle T=0$ part reads

$\mathcal{M}^{(0)}=\frac{-8m^{4}-4\cos^{4}(\theta)\left(m^{2}-\omega^{2}\right)^{2}+8m^{2}\cos^{2}(\theta)\left(m^{2}-\omega^{2}\right)+8m^{2}\omega^{2}+4\omega^{4}}{\left[\cos^{2}(\theta)\left(m^{2}-\omega^{2}\right)+\omega^{2}\right]^{2}}$
(39)

while the thermal amplitude has the following expression

$\displaystyle\displaystyle\mathcal{M}^{(\beta)}=$
$\displaystyle\displaystyle-4\pi\frac{\Gamma_{1}+3\xi_{t}^{2}m^{2}}{\xi_{t}^{2}}\left[\delta^{2}(2\omega+\xi_{t})+\delta^{2}(2\omega-\xi_{u})\right]$
$\displaystyle\displaystyle-4\pi\frac{\Gamma_{1}+3\xi_{u}^{2}m^{2}}{\xi_{u}^{2}}\left[\delta^{2}(2\omega+\xi_{u})+\delta^{2}(2\omega-\xi_{u})\right]$
$\displaystyle\displaystyle+8\pi^{2}\frac{\left(3\xi_{t}^{2}m^{2}-\Gamma_{1}\right)}{\xi_{t}^{2}}\delta(2\omega+\xi_{t})\delta(2\omega-\xi_{t})+8\pi^{2}\frac{\left(3\xi_{u}^{2}m^{2}-\Gamma_{1}\right)}{\xi_{u}^{2}}\delta(2\omega+\xi_{u})\delta(2\omega-\xi_{u})$
$\displaystyle\displaystyle+8\pi^{2}\frac{\Gamma_{2}+(m^{2}-4\omega^{2})\xi_{t}\xi_{u}}{\xi_{t}\xi_{u}}\left[\delta(2\omega+\xi_{t})\delta(2\omega+\xi_{u})+\delta(2\omega-\xi_{t})\delta(2\omega-\xi_{u})\right]$
$\displaystyle\displaystyle+8\pi^{2}\frac{\Gamma_{2}-(m^{2}-4\omega^{2})\xi_{t}\xi_{u}}{\xi_{t}\xi_{u}}\left[\delta(2\omega+\xi_{t})\delta(2\omega-\xi_{u})+\delta(2\omega-\xi_{t})\delta(2\omega+\xi_{u})\right],$
(40)

where we have defined, by simplicity of notation, the following quantities

$\displaystyle\displaystyle\Gamma_{1}$
$\displaystyle\displaystyle=18m^{4}-19\omega^{2}m^{2}-18\omega
m^{2}\sqrt{\omega^{2}-m^{2}}\cos(\theta)-\omega^{2}m^{2}\cos(2\theta)-\omega^{4}+\omega^{4}\cos(2\theta),$
(41) $\displaystyle\displaystyle\Gamma_{2}$
$\displaystyle\displaystyle=6m^{4}-8\omega^{2}m^{2}+8\omega^{4}+2\omega\sqrt{\omega^{2}-m^{2}}\left(4\omega^{2}-3m^{2}\right)\cos(\theta),$
(42) $\displaystyle\displaystyle\xi_{\text{t}}$
$\displaystyle\displaystyle=\sqrt{2\omega\cos(\theta)\sqrt{\omega^{2}-m^{2}}+2m^{2}-2\omega^{2}},$
(43) $\displaystyle\displaystyle\xi_{\text{u}}$
$\displaystyle\displaystyle=\sqrt{-2\omega\cos(\theta)\sqrt{\omega^{2}-m^{2}}+2m^{2}-2\omega^{2}}.$
(44)

Hence, in terms of the definition (1) the differential cross section for the
Breit-Wheeler process at finite temperature is given as

$\frac{d\sigma}{d\Omega}=\frac{1}{64\pi^{2}(4\omega^{2})}\langle|\widehat{\mathcal{M}}_{fi}|^{2}\rangle=\frac{e^{4}F^{2}(\beta)}{256\pi^{2}\omega^{2}}\left[\left.\frac{d\sigma}{d\Omega}\right|_{(0)}+\frac{1}{(1+e^{2\beta\omega})^{2}}\left.\frac{d\sigma}{d\Omega}\right|_{(\beta)}\right]$
(45)

where the contribution at $\displaystyle T=0$ is given by

$\left.\frac{d\sigma}{d\Omega}\right|_{(0)}=\mathcal{M}^{(0)},$ (46)

while the finite temperature contribution is

$\left.\frac{d\sigma}{d\Omega}\right|_{(\beta)}=\mathcal{M}^{(\beta)}$ (47)

where the zero and finite temperature amplitudes
$\displaystyle\mathcal{M}^{(0)}$ and $\displaystyle\mathcal{M}^{(\beta)}$ are
given by Eqs. (39) and (40), respectively.

Some comments are in order about the expression (45): it is interesting to
observe that at the limit when the temperature goes to zero
($\displaystyle\beta\rightarrow\infty$) it implies that
$\displaystyle\frac{1}{(1+e^{2\beta\omega})^{2}}\to 0$ and $\displaystyle
F(\beta)\to 1$, then, differential cross section (45) is purely the known
result at zero temperature [16, 41, 24].

Since all angular dependence of the differential cross section (45) is present
in the amplitudes squared $\displaystyle\mathcal{M}^{(0)}$ and
$\displaystyle\mathcal{M}^{(\beta)}$, one can evaluate the angular
integrations and obtain the total cross section for any value of
$\displaystyle\beta$

$\sigma_{\rm
BW}=\frac{e^{4}}{32\pi\omega^{2}}F^{2}\left(\beta\right)\left[\mathbb{I}_{1}+\frac{1}{\left(1+e^{2\beta\omega}\right)^{2}}\mathbb{I}_{2}\right]$
(48)

in which it was defined

$\displaystyle\displaystyle\mathbb{I}_{1}=-\frac{2}{\chi}-2+\frac{3\chi+\frac{1}{\chi}-2}{\sqrt{\chi\left(\chi-1\right)}}\ln\left(\frac{1+\sqrt{1-\frac{1}{\chi}}}{1-\sqrt{1-\frac{1}{\chi}}}\right)$
(49)

and

$\displaystyle\displaystyle\mathbb{I}_{2}=$
$\displaystyle\displaystyle\frac{\pi^{2}m}{\sqrt{\chi\left(\chi-1\right)}}\left(-10+21\chi-8\chi^{2}\right)\delta\left(\sqrt{\chi}\right)$
(50)
$\displaystyle\displaystyle+\frac{4\pi^{2}}{\sqrt{\chi\left(\chi-1\right)}}\left(9+20\chi^{2}-24\chi\right)\delta\left(3\chi-1\right)$
(51)

in which we have introduced the parameter
$\displaystyle\chi=\frac{\omega^{2}}{m^{2}}$ by means of simplicity of
notation. One can observe that the term $\displaystyle\mathbb{I}_{1}$
corresponds to the usual $\displaystyle T=0$ contribution, while
$\displaystyle\mathbb{I}_{2}$ corresponds to the purely thermal part.

It is important to observe that the delta functions present in
$\displaystyle\mathbb{I}_{2}$ in the expression (48) play the role of filters
for such values, this is a consequence of the approximation that has been
used, a tree-level approximation. One can naively find a singularity in
$\displaystyle\mathbb{I}_{2}$ present in the term proportional to
$\displaystyle\delta\left(\sqrt{\chi}\right)$ or even that
$\displaystyle\delta\left(3\chi-1\right)$ leads to a imaginary contribution to
the cross section. However, these observations are simply apparent, because
when we analyze the full expression for the cross-section (48) a natural
condition arises upon these delta functions arise. Once the cross-section
represents a measurable physical quantity it must be real. Therefore, applying
this condition on (48), we realize that the terms
$\displaystyle\mathbb{I}_{1}$ and $\displaystyle\mathbb{I}_{2}$ must be real.
This implies a strong constraint upon the parameter $\displaystyle\chi$, i.e.
$\displaystyle\chi>1$ (i.e. $\displaystyle\omega^{2}>m^{2}$). Imposing this
condition on $\displaystyle\mathbb{I}_{2}$, we realize that both delta
functions $\displaystyle\delta\left(\sqrt{\chi}\right)$ and
$\displaystyle\delta\left(3\chi-1\right)$ become zero and the cross-section
depends only on the term $\displaystyle\mathbb{I}_{1}$. Hence, the real and
finite cross-section at finite temperature becomes

$\sigma_{\rm
BW}=\frac{e^{4}}{32\pi\omega^{2}}F^{2}\left(\beta\right)\mathbb{I}_{1}.$ (52)

It is worth stress that the expression for the Breit-Wheeler process cross
section (52) is valid for any value of $\displaystyle\beta$. However, this
thermal dependence is rather complicated in the factor $\displaystyle
F^{2}\left(\beta\right)$. Thus, one can consider, by means of elucidation, the
high temperature limit, where one can get a better grasp of the thermal
behavior of the cross section (52). Hence, at the high temperature regime, we
obtain that

$\displaystyle\displaystyle\sigma_{\rm BW}^{\rm
high}\simeq\frac{e^{4}}{8\pi\omega^{2}}\mathbb{I}_{1}e^{-\frac{3\beta\omega}{4}}.$
(53)

This expression is accordance with our expectation, since it is well known
that at the high-temperature limit the Fermi-Dirac and Bose-Einstein
distributions reproduce the Maxwell-Boltzmann statistics:
$\displaystyle\frac{1}{e^{\beta\omega}\mp 1}\simeq e^{-\beta\omega}$.

##  IV Conclusion

In this paper we have studied thermal contributions to the Breit-Wheeler
process within the thermofield dynamics approach. Due to phenomenological
interest (QED in intense electromagnetic fields) the Breit-Wheeler process,
which consists of a pair production through the collision of two photons, has
been considered in many contexts. However, its description in a thermal bath
is somewhat controversial, since many papers (results) disagree with each
other. Hence, we approach this topic within the thermofield dynamics
formalism, which is very robust in order to shed some new light and clarity
(some of) the controversies.

We started our analysis by reviewing the main aspects of the thermofield
dynamics applied to the evaluation of the (thermal) differential cross
section. With this purpose, we discussed the second-order thermal transition
amplitude related with the process $\displaystyle\gamma+\gamma\to
e^{-}+e^{+}$, where the interacting part of the process is described by the
usual QED coupling $\displaystyle e\bar{\psi}(x)\gamma^{\mu}\psi(x)A_{\mu}(x)$
(and its dual). One important point in the calculation of the transition
amplitude is the use of the Bogoliubov transformation between the
$\displaystyle T=0$ and thermal creation/annihilation operators, this is
necessary because the in/out states are thermal states (operators), while the
operators in the interacting Lagrangian are not thermal. However, this is not
the only thermal contribution to the scattering matrix, but the fermionic
propagator also has a finite temperature part.

After establishing the main operations of the thermofield dynamics formalism,
we considered the center of mass frame to compute the squared transition
amplitude. With this step we finally obtained the desired expression for the
differential cross section Breit-Wheeler process at finite temperature. In the
zero temperature limit ($\displaystyle\beta\rightarrow\infty$) our results
reproduced the known literature [16, 41, 24]. Remarkably, the thermal
dependence of the differential cross section (45) had a simple form which
allowed to compute the angular integration in order to obtain the Breit-
Wheeler process cross section valid for any value of $\displaystyle\beta$.
Moreover, in order to have a better grasp about the thermal behavior of the
cross section, we considered the high temperature regime, which led to a
thermal dependence in terms of the Maxwell-Boltzmann distribution as we would
expect.

## Acknowledgments

The authors thank Professor Ademir E. Santana for his interest and comments in
the work. This work by A. F. S. is partially supported by National Council for
Scientific and Technological Development - CNPq project No. 313400/2020-2. R.
B. acknowledges partial support from Conselho Nacional de Desenvolvimento
Científico e Tecnológico (CNPq Project No. 306769/2022-0).

## References

  * [1] A. Gonoskov, T. G. Blackburn, M. Marklund and S. S. Bulanov, “Charged particle motion and radiation in strong electromagnetic fields,” Rev. Mod. Phys. 94 (2022) no.4, 045001
  * [2] A. Fedotov, A. Ilderton, F. Karbstein, B. King, D. Seipt, H. Taya and G. Torgrimsson, “Advances in QED with intense background fields,” Phys. Rept. 1010 (2023), 1-138
  * [3] A. Di Piazza, C. Muller, K. Z. Hatsagortsyan and C. H. Keitel, “Extremely high-intensity laser interactions with fundamental quantum systems,” Rev. Mod. Phys. 84 (2012), 1177
  * [4] C. G. W. D. Colgan, “Laser-plasma interactions as tools for studying processes in quantum electrodynamics,” PhD Thesis Imperial College (2022)
  * [5] J. S. Schwinger, “On gauge invariance and vacuum polarization,” Phys. Rev. 82 (1951) 664-679
  * [6] F. Sauter, “Uber das Verhalten eines Elektrons im homogenen elektrischen Feld nach der relativistischen Theorie Diracs,” Z. Phys. 69 (1931) 742-764
  * [7] I. Blumenfeld, C. E. Clayton, F. J. Decker, M. J. Hogan, C. Huang, R. Ischebeck, R. Iverson, C. Joshi, T. C. Katsouleas and N. Kirby, et al. “Energy doubling of 42 GeV electrons in a metre-scale plasma wakefield accelerator,” Nature 445 (2007), 741-744
  * [8] H. Abramowicz, U. Acosta, M. Altarelli, R. Aßmann, Z. Bai, T. Behnke, Y. Benhammou, T. Blackburn, S. Boogert and O. Borysov, et al. “Conceptual design report for the LUXE experiment,” Eur. Phys. J. ST 230 (2021) no.11, 2445-2560
  * [9] A. Aschikhin, C. Behrens, S. Bohlen, J. Dale, N. Delbos, L. di Lucchio, E. Elsen, J. H. Erbe, M. Felber and B. Foster, et al. “The FLASHForward Facility at DESY,” Nucl. Instrum. Meth. A 806 (2016), 175-183
  * [10] E. Adli et al. [AWAKE], “Acceleration of electrons in the plasma wakefield of a proton bunch,” Nature 561 (2018) no.7723, 363-367
  * [11] D. L. Burke, R. C. Field, G. Horton-Smith, T. Kotseroglou, J. E. Spencer, D. Walz, S. C. Berridge, W. M. Bugg, K. Shmakov and A. W. Weidemann, et al. “Positron production in multi - photon light by light scattering,” Phys. Rev. Lett. 79 (1997), 1626-1629
  * [12] G. Breit and J. A. Wheeler, “Collision of two light quanta,” Phys. Rev. 46 (1934) no.12, 1087-1091
  * [13] J. Adam et al. [STAR], “Measurement of $\displaystyle e^{+}e^{-}$ Momentum and Angular Distributions from Linearly Polarized Photon Collisions,” Phys. Rev. Lett. 127 (2021) no.5, 052302
  * [14] V. I. Ritus, “Quantum effects of the interaction of elementary particles with an intense electromagnetic field,” J. Russ. Laser Res. 6, 497-617 (1985)
  * [15] C. Bula et al. [E144], “Observation of nonlinear effects in Compton scattering,” Phys. Rev. Lett. 76 (1996), 3116-3119
  * [16] M. Bottcher and R. Schlickeiser, “The Pair production spectrum from photon-photon annihilation,” Astron. Astrophys. 325 (1997), 866
  * [17] G. Baur, K. Hencken and D. Trautmann, “Photon-photon physics in very peripheral collisions of relativistic heavy ions,” J. Phys. G 24 (1998), 1657-1692
  * [18] M. Kłusek-Gawenda, R. McNulty, R. Schicker and A. Szczurek, “Light-by-light scattering in ultraperipheral heavy-ion collisions at low diphoton masses,” Phys. Rev. D 99 (2019) no.9, 093013
  * [19] T. Heinzl, “Strong-Field QED and High Power Lasers,” Int. J. Mod. Phys. A 27 (2012), 1260010
  * [20] T. G. Blackburn, C. P. Ridgers, J. G. Kirk and A. R. Bell, “Quantum radiation reaction in laser-electron beam collisions,” Phys. Rev. Lett. 112 (2014), 015001
  * [21] A. Di Piazza, “Nonlinear Breit-Wheeler pair production in a tightly focused laser beam,” Phys. Rev. Lett. 117 (2016) no.21, 213201
  * [22] D. Seipt, T. Heinzl, M. Marklund and S. S. Bulanov, “Depletion of Intense Fields,” Phys. Rev. Lett. 118 (2017) no.15, 154803
  * [23] A. I. Titov, B. Kämpfer and H. Takabe, “Nonlinear Breit-Wheeler process in short laser double pulses,” Phys. Rev. D 98 (2018) no.3, 036022
  * [24] X. Ribeyre et al., “Effect of differential cross section in Breit-Wheeler pair production,” Plasma Phys. Control. Fusion 60, 104001 (2018).
  * [25] F. C. Salgado, K. Grafenstein, A. Golub, A. Döpp, A. Eckey, D. Hollatz, C. Müller, A. Seidel, D. Seipt and S. Karsch, et al. “Towards pair production in the non-perturbative regime,” New J. Phys. 23 (2021) no.10, 105002
  * [26] A. Golub, S. Villalba-Chávez, H. Ruhl and C. Müller, “Linear Breit-Wheeler pair production by high-energy bremsstrahlung photons colliding with an intense X-ray laser pulse,” Phys. Rev. D 103 (2021) no.1, 016009
  * [27] A. Golub, S. Villalba-Chávez and C. Müller, “Strong-field Breit-Wheeler pair production in QED2+1,” Phys. Rev. D 103 (2021) no.9, 096002
  * [28] T. Adamo, A. Ilderton and A. J. MacLeod, “Particle-beam scattering from strong-field QED,” Phys. Rev. D 104 (2021) no.11, 116013
  * [29] B. Kettle, D. Hollatz, E. Gerstmayr, G. M. Samarin, A. Alejo, S. Astbury, C. Baird, S. Bohlen, M. Campbell and C. Colgan, et al. “A laser-plasma platform for photon-photon physics: the two photon Breit-Wheeler process,” New J. Phys. 23 (2021) no.11, 115006
  * [30] Q. Zhao, L. Tang, F. Wan, B. C. Liu, R. Y. Liu, R. Z. Yang, J. Q. Yu, X. G. Ren, Z. F. Xu and Y. T. Zhao, et al. “Signatures of linear Breit-Wheeler pair production in polarized $\displaystyle\gamma$$\displaystyle\gamma$ collisions,” Phys. Rev. D 105 (2022) no.7, L071902
  * [31] T. G. Blackburn and B. King, “Higher fidelity simulations of nonlinear Breit-Wheeler pair creation in intense laser pulses,” Eur. Phys. J. C 82 (2022) no.1, 44
  * [32] J. D. Brandenburg, J. Seger, Z. Xu and W. Zha, “Report on Progress in Physics: Observation of the Breit-Wheeler Process and Vacuum Birefringence in Heavy-Ion Collisions,” arXiv:2208.14943 [hep-ph].
  * [33] M. Turner, S. S. Bulanov, C. Benedetti, A. J. Gonsalves, W. P. Leemans, K. Nakamura, J. van Tilborg, C. B. Schroeder, C. G. R. Geddes and E. Esarey, “Strong-field QED experiments using the BELLA PW laser dual beamlines,” Eur. Phys. J. D 76 (2022) no.11, 205
  * [34] A. J. MacLeod, P. Hadjisolomou, T. M. Jeong and S. V. Bulanov, “All-optical nonlinear Breit-Wheeler pair production with $\displaystyle\gamma$-flash photons,” Phys. Rev. A 107 (2023) no.1, 012215
  * [35] L. Medina and M. C. Ogilvie, “Schwinger pair production at finite temperature,” Phys. Rev. D 95 (2017) 056006
  * [36] B. King, H. Gies and A. Di Piazza, “Pair production in a plane wave by thermal background photons,” Phys. Rev. D 86 (2012), 125007 [erratum: Phys. Rev. D 87 (2013) no.6, 069905]
  * [37] O. Gould and A. Rajantie, “Thermal Schwinger pair production at arbitrary coupling,” Phys. Rev. D 96 (2017) 076002
  * [38] A. R. Brown, “Schwinger pair production at nonzero temperatures or in compact directions,” Phys. Rev. D 98 (2018) 036008
  * [39] O. Gould, A. Rajantie, and C. Xie, “Worldline sphaleron for thermal Schwinger pair production,” Phys. Rev. D 98 (2018) 056022
  * [40] M. Korwar and A. M. Thalapillil, “Finite temperature Schwinger pair production in coexistent electric and magnetic fields,”  Phys. Rev. D 98 (2018) 076016
  * [41] G. Torgrimsson, “Thermally versus dynamically assisted Schwinger pair production,” Phys. Rev. D 99 (2019) 096007
  * [42] O. Gould, S. Mangles, A. Rajantie, S. Rose, and C. Xie, “Observing thermal Schwinger pair production,” Phys. Rev. A 99 (5) (2019) 052120
  * [43] Y. L. Wang, H. B. Sang, and B. S. Xie, “Schwinger pair production correction in thermal system,” Phys. Rev. D 100 (11) (2019) 116016
  * [44] Y. Takahashi and H. Umezawa, “Thermo Field Dynamics,” Collective Phenomena 2 (1975), 55; Reprinted in Int. J. Mod. Phys.B 10 (1996), 1755.
  * [45] H. Umezawa, H. Matsumoto and M. Tachiki, “Thermofield Dynamics and Condensed States,” North-Holland, Amsterdam (Netherlands) (1982).
  * [46] H. Umezawa, “Advanced field theory: Micro, macro, and thermal physics,” American Institute of Physics (1995)
  * [47] F. C. Khanna, A. P. C. Malbouisson, J. M. C. Malbouisson and A. R. Santana, “Thermal quantum field theory - Algebraic aspects and applications,” World Scientific Publishing Company (2009)
  * [48] T. Matsubara, “A New Approach to Quantum-Statistical Mechanics,” Prog. Theor. Phys. 14 (1955), 351.
  * [49] J. Schwinger, “Brownian Motion of a Quantum Oscillator,”  J. Math. Phys. 2 (1961), 407.
  * [50] J. I. Kapusta and C. Gale, “Finite-Temperature Field Theory: Principles and Applications,” 2nd ed. Cambridge University Press (2011).