# Dynamic Graph Neural Network for
Traffic Forecasting in Wide Area Networks

Tanwi Mallick  Argonne National Laboratory
Lemont, Illinois
<EMAIL_ADDRESS>Mariam Kiran  Lawrence Berkeley
National Laboratory
Berkeley, California
<EMAIL_ADDRESS>Bashir Mohammed  Lawrence Berkeley
National Laboratory
Berkeley, California
<EMAIL_ADDRESS>Prasanna Balaprakash  Argonne National Laboratory
Lemont, Illinois
<EMAIL_ADDRESS>
###### Abstract

Wide area networking infrastructures (WANs), particularly science and research
WANs, are the backbone for moving large volumes of scientific data between
experimental facilities and data centers. With demands growing at exponential
rates, these networks are struggling to cope with large data volumes, real-
time responses, and overall network performance. Network operators are
increasingly looking for innovative ways to manage the limited underlying
network resources. Forecasting network traffic is a critical capability for
proactive resource management, congestion mitigation, and dedicated transfer
provisioning. To this end, we propose a nonautoregressive graph-based neural
network for multistep network traffic forecasting. Specifically, we develop a
dynamic variant of diffusion convolutional recurrent neural networks to
forecast traffic in research WANs. We evaluate the efficacy of our approach on
real traffic from ESnet, the U.S. Department of Energy’s dedicated science
network. Our results show that compared to classical forecasting methods, our
approach explicitly learns the dynamic nature of spatiotemporal traffic
patterns, showing significant improvements in forecasting accuracy. Our
technique can surpass existing statistical and deep learning approaches by
achieving $\approx$20% mean absolute percentage error for multiple hours of
forecasts despite dynamic network traffic settings.

##  I Introduction

Large scientific experimental facilities, with high-speed data production
rates, are bringing enormous data movement and transfer challenges to the
underlying network infrastructure that supports distributed science workflows.
With this relentless growth, there is a need to develop proactive planning,
resource allocation, and provisioning methods to manage the current backbone
bandwidth of research wide area networks (R-WANs), while keeping costs low
[1]. One of these proactive strategies could rely on network operators’
ability to forecast future traffic on the network based on the current state
of the network. For example, if we can forecast the traffic for several hours
in the future, we can minimize the effect of congestion by diverting flows
from congested links to free or less congested links; similarly, dedicated
lines can be assigned to large transfers so that they do not interfere with
other data transfers. The ability to forecast congestion patterns a few hours
ahead can enable new traffic management strategies with traditional network
protocol design, essentially allowing one to utilize underused bandwidth more
efficiently and minimize overall congestion in the network infrastructure.

Network traffic forecasting in R-WANs is a formidable task, owing to a lack of
regular patterns in how users access and perform data transfers over the
network [2]. Compared to Internet WANs, which have periodic patterns [3],
R-WANs have random traffic spikes that are difficult to understand and
anticipate. The traffic on R-WANs depends on which science experiments and
devices are running and which groups are involved, and it is characterized by
high-variability data transfers lasting minutes or even hours [4].

Network monitoring tools such as Simple Network Management Protocol (SNMP),
sflow, and netflow [5] allow collecting traffic information on network nodes
and flow transfers as time-stamped data recording gigabytes of log files (GB).
Most monitoring tools collect data at 30-sec intervals, giving a fine-grained
view that includes other key features such as protocols used (e.g., TCP, UDP),
interfaces, source and destination IP addresses, and even flow speeds. These
data sets can be leveraged for developing data-driven forecasting methods.
Classical statistical time-series forecasting methods such as ARIMA and Holt-
Winters have been investigated for network traffic forecasting [6], but have
been less effective for R-WANs, as regular and seasonal patterns do not exist
[1]. Forecasting methods based on classical ML methods such as random forest
and SVM have been developed to provide per-site forecasts [7], but they fail
to consider the whole network and its spatial patterns and connectivity,
making them less robust with respect to the dynamic nature of the R-WAN
traffic [8]. Recent improvements in data collection present new opportunities
for developing deep learning (DL) methods for forecasting R-WAN traffic [9].
For example, variants of convolutional neural networks (CNNs) and long-short-
term memory methods (LSTMs) have been shown to provide better accuracy than
classical statistical techniques for data sets with seasonality [10].
Nevertheless, existing DL methods for network forecasting are out-of-the-box
methods and are not particularly customized for the R-WAN. They do not take
into account the spatial correlation of the entire network and the dynamic
temporal correlation of R-WANs.

We model the R-WAN networking infrastructure as a graph, where nodes and edges
correspond to sites and their connectivity, modeling network traffic as time
series. We propose a nonautoregressive graph neural network approach to
forecast traffic for multiple time steps ahead. Specifically, we develop a
dynamic diffusion convolutional recurrent neural network (DDCRNN) that
considers data movement as a diffusion process from one node to another based
on connectivity. The spatial and temporal correlations are learned through
graph diffusion convolution and recurrent units, respectively. Consequently, a
single model is trained and used to forecast traffic on the entire R-WAN
network. The DDCRNN that we propose is based on the diffusion convolutional
recurrent neural network (DCRNN), which was originally developed for highway
transportation forecasting [11, 12]. The key difference between the highway
traffic DCRNN and our DDCRNN is the way in which the connectivity between the
nodes is considered. In the DCRNN, the connectivity is static and computed on
the basis of (driving) distance; in DDCRNN, the connectivity is dynamic and
computed on the basis of the current state of the traffic in the network. This
approach is designed to explicitly model the dynamic nature of the R-WAN
traffic. Our main contributions are as follows:

  * •

We develop a dynamic diffusion graph recurrent neural network architecture to
forecast traffic for multiple time steps on an R-WAN that is characterized by
a lack of regular patterns and seasonality.

  * •

We demonstrate the effectiveness of the proposed method on real data from the
Energy Sciences Network (ESnet), a high-performance R-WAN built by the U.S.
Department of Energy (DOE) to support U.S. scientific research. We show that
our approach explicitly learns the dynamic nature of spatiotemporal traffic
patterns in this R-WAN and outperforms existing statistical and DL methods
used for traffic forecasting by achieving $\approx$20% mean absolute
percentage error for multiple hours of forecasts.

Using real traffic data from ESnet, our model captures key patterns among
sites and links, highlighting several interesting behaviors and relationships
that were not previously known, such as some sites having patterns that are
more amenable for forecasting than others.

##  II Related Work

Wolski et al. [13] developed a TCP performance network forecasting method to
help schedule computations over distributed networks. These authors utilized
statistical approaches such as ARIMA, Holt-Winters and Hidden Markov, and were
successful in predicting a few time-steps by modeling traffic as stochastic
processes [14]. Nevertheless, these methods are dependent on finding
seasonality and using it to improve predictions. Network traffic spikes
randomly with sudden data transfers and lacks seasonal patterns [15],
affecting the forecasting accuracy of these methods [6]. Moreover, these
approaches also fail to learn long-range dependency [16].

Research in software-defined networking (SDN) promises to provide flexible
solutions for building agile networks and leveraging active monitoring,
prediction, and informed decision-making [17]. Google [18] used SDN to
optimize link usage by doing ”what-if” scenarios to schedule transfers. Using
Multi-Protocol Label Switching, advanced forwarding schemes can control and
optimize flows for packet forwarding [19]. However, the optimization
techniques discussed do not exploit forecasting algorithms to make decisions.

Most R-WAN network links are underutilized, especially as routing protocols
use path-finding algorithms rather than accounting for current traffic
utilization. In prior work, dynamic congestion-based routing algorithms have
been developed to adapt as traffic changes [20]. However, these algorithms are
susceptible to oscillatory behaviors and performance degradation, as shown in
[21].

The European R&E network, GEANT, used LSTM models to forecast traffic on
European links, but only demonstrated 15-min forecasts [21]. Similarly stacked
autoencoders [9] were used for 15-, 30- and 60-minute forecasts. However,
these authors did not use highly dynamic R-WANs such as ESnet and exploited
relatively simpler and shorter forecasts. Recently, DCRNN implementations,
without any modifications, were successfully applied to forecast congestion
events in a WAN [22]. Network traffic forecasting has explored statistical and
simpler ML models as seen in [6][7] [8], but prediction accuracy is seen to
suffer because of lack of seasonality in R-WAN traffic traces.

In this paper, we show that DCRNN alone is not effective, as it does not
consider the dynamic spatial and temporal traffic patterns in R-WAN networks.

##  III U.S. Research Network: ESnet

The DOE science network ESnet provides services to more than 50 research sites
and universities, including supercomputing facilities and major scientific
instruments. ESnet also connects to 140 research and commercial networks,
permitting geography-free collaboration around the world.

Figure 1: Full topology of ESnet. Figure 2: Traffic patterns in January 2018
on all transatlantic links in ESnet, including London, New York, Amsterdam,
Boston, AOFA, CERN, and Washington, D.C.

We use real traffic traces of two-way traffic data, collected from all of
2018, across the 48 sites in the U.S. and Europe. Figure 1 shows the full
network topology of the network, divided among five time zones. All data are
collected in GMT. Sample network traffic patterns on transatlantic links in
ESnet appear in Figure 2, which shows the lack of temporal patterns in the
network traffic.

###  III-A SNMP Data Collection

Network operators monitor link capacity and data movement across routers using
SNMP [23]. Table I shows a snapshot of these data during a two-way transfer
between Sunnyvale and Sacramento in California.

The traffic is collected in moving GBs across router interfaces at 30-sec
intervals. These data are then aggregated to 1-hour intervals and modeled as
discrete observations. Aggregating traffic compounds the burst patterns [20]
and reduces the model complexity, allowing faster training and predictions.

Timestamp | SACR_SUNN_in (GB) | SACR_SUNN_out (GB)
---|---|---
1514822400 | 14110930202 | 1025131246
1514826000 | 13453619303 | 9191557943
1514829600 | 12168879944 | 7793842045
1514833200 | 11231198033 | 7097237528
1514836800 | 10780847622 | 8048293939
TABLE I: Sample timestamped traffic trace collected from one router collecting
traffic in both directions between Sacramento and Sunnyvale.

The data that we used for our study spans from 1 January 2018 to 31 December
2018, covering 48 sites with 96 traffic traces of two-way traffic.

While building hourly traffic summaries, we found that router interfaces
sometimes miss recording traffic movements at specific intervals. To address
this issue, we calculated the average data values for the missing points from
the surrounding values to fill in the gaps. These missing values were only
seen in the Washington-Chicago link during a 1-week interval in November.

##  IV Dynamic Diffusion Convolution Recurrent Neural Network

We model the R-WAN as a graph $G=(V,E,A)$, where $V$ is a set of $N$ nodes
that represents the site; $E$ is a set of directed edges representing the
connection between nodes, and $A\in R^{N\times N}$ is the weighted adjacency
matrix representing the strength of connectivity between nodes. Given the
historical traffic observations at each node of the graph, the goal is to
learn a function f(.) that takes traffic observations for $T^{\prime}$ time
steps as input to forecast the traffic for the next $T$ time steps:

$X(t-T^{\prime}+1),...,X(t);G\xrightarrow{\text{f(.)}}X(t+1),...,X(t+T)$

Figure 3: DDCRNN model architecture. It takes an adjacency matrix computed
from the current state of the network traffic among the sites of the WAN
network topology and the traffic as a time-series at each node of the graph.
The encoder-decoder deep neural network is used to forecast the network
traffic for mutiple time steps.

DDCRNN utilizes a graph-based encoder-decoder architecture that models spatial
and temporal patterns through diffusion convolution operation on a graph and
gated recurrent units (GRUs), respectively. Specifically, the matrix
multiplication operations of the GRU cell are replaced by a diffusion
convolution operation to perform graph convolution on the time-series data. A
unit in the DDCRNN architecture is given by:

$\begin{array}[]{lcl}r^{t}&=&\sigma(W_{r\bigstar G}[X_{t},h_{t-1}]+b_{r})\\\
u^{t}&=&\sigma(W_{u\bigstar G}[X_{t},h_{t-1}]+b_{u})\\\
c^{t}&=&\tanh(W_{c\bigstar G}[X_{t}(r_{t}\odot h_{t-1}]+b_{c})\\\
h_{t}&=&u_{t}\odot h_{t-1}+(1-u_{t})\odot c_{t},\end{array}$ (1)

where $X_{t}$ and $h_{t}$ denote the input and final state at time $t$,
respectively; $r_{t}$, $u_{t}$, and $c_{t}$ are the reset gate, update gate,
and cell state at time $t$, respectively; $\bigstar G$ denotes the diffusion
convolution operation; and $W_{r},W_{u},$ and $W_{c}$ are parameters of the
GRU cell. The diffusion convolution operation [11] on the graph $G$ and the
input data $X$ is defined as:

$W_{\bigstar
G}X=\sum_{d=0}^{K-1}(W_{O}(D_{O}^{-1}A)^{d}+W_{I}(D_{I}^{-1}A)^{d})X,$ (2)

where $K$ is a maximum diffusion step; $D_{O}^{-1}A$ and $D_{I}^{-1}A$ are
transition matrices of the diffusion process and the reverse one,
respectively; $D_{O}$ and $D_{I}$ are the in-degree and out-degree diagonal
matrices, respectively, and $W_{O}$ and $W_{I}$ are the learnable filters for
the bidirectional diffusion process.

To model the dynamic nature of the R-WAN traffic, DDCRNN uses a dynamic
weighted adjacency matrix estimated from the time-series data. Specifically,
instead of keeping $A$ as static, for each input sequence of $T^{\prime}$ time
steps on $N$ nodes (for both training and inference), DDCRNN computes the
linear correlation coefficient values between sites that have data transfer
between them and uses it in Eq. 2. Consequently, the elements in the adjacency
matrix $\tilde{A}$ at time $T$ are given by:

$\tilde{A}_{ij}=\rho_{X_{i},X_{j}}=\frac{cov(X_{i},X_{j})}{\sigma_{X_{i}},\sigma_{X_{j}}},$
(3)

where, $A_{ij}$ represents the edge weight between the nodes $N_{i}$ and
$N_{j}$; $X_{i}$ and $X_{j}$ are the time-series data of nodes $N_{i}$ and
$N_{j}$; $\rho$ is the Pearson correlation coefficient between $X_{i}$ and
$X_{j}$; $cov$ is the covariance; $\sigma_{X_{i}}$ is the standard deviation
of $X_{i}$; and $\sigma_{X_{j}}$ is the standard deviation of $X_{j}$.

The overall architecture of DDCRNN is shown in Figure 3. There are two inputs
to the DDCRNN model: 1) an adjacency matrix representing the traffic
correlation among the sites of the WAN network topology and 2) the time-series
data or the traffic statistic at each node of the graph. The encoder of the
DDCRNN network encodes the input into a fixed-length vector and passes it to
the decoder. The decoder forecasts future traffic conditions. During the
training phase, the time series of $T^{\prime}$ time steps on $N$ nodes is fed
in as input; the correlation matrix is computed on the time-series data and
given as the weighted adjacency matrix for diffusion convolution. The encoder
architecture takes the data and the decoder is used to forecast the output of
the next $T$ time steps. The learnable weights of DDCRNN are trained using a
minibatch stochastic gradient, using a mean absolute error (MAE) loss
function.

The DDCRNN for R-WAN traffic forecasting is based on a DCRNN that was
originally developed for forecasting traffic on highway networks [11]. The
network traffic forecasting is similar to highway traffic forecasting, where
the traffic diffuses from one node to the other node based on the network
graph connectivity. However, the key difference stems from the spatial
connectivity and temporal regularity. In a highway setting, the traffic
exhibits regular patterns across weekdays, peak hours, regions (i.e.,
distribution of vehicles going from one region to another region), and
seasons. However, such patterns cannot be seen in a R-WAN network because
users can initiate large data transfers at any time. Therefore, DDCRNN
utilizes a dynamic adjacency matrix computed from the current state of the
network traffic. During training, DDCRNN learns to diffuse traffic on the
graph under different weighted adjacency matrix settings. Consequently, during
inference, the trained model can perform forecasts based on the time series of
$T^{\prime}$ steps and the weighted adjacency matrix computed on it.

##  V Experiments

We used ESnet traffic traces for 1 year with 1-hour time resolution. We
grouped data as follows: from 2018-01-01 to 2018-09-13 for training (70%);
from 2018-09-13 to 2018-10-20 (10%) for validation; and from 2018-10-20 to
2018-12-31 (20%) for testing. We measured the network traffic as bandwidth in
GB/s. Each site has two bandwidth values: incoming and outgoing. Therefore, we
mapped the 48 physical sites to a 96-node graph by considering two nodes
(incoming and outgoing) for each physical site, as they have separate fiber
optic links. Since data transfer varies significantly from site to site, we
computed min-max scalar transformation on the training data and use it
normalize bandwidth values of training, validation, and testing data.

Our DDCRNN implementation is based on the open-source DCRNN implementation
[11]. We used the following hyperparameter values for DDCRNN: batch size: 64,
number of epochs: 30, maximum diffusion steps: 2, number of RNN layers: 2,
number of RNN units per layers: 16, max_grad_norm: 5, initial learning rate:
0.01, and learning rate decay: 0.1. These hyperparameter values were set as
default for the open-source DCRNN implementation.

The DDCRNN training was performed on a single node of Cooley GPU cluster at
Argonne Leadership Computing Facility. The node consists of two 2.4 GHz Intel
Haswell E5-2620 v3 processors (6 cores per CPU, 12 cores total), one NVIDIA
Tesla K80 (two GPUs per node), 384 GB of RAM per node, and 24 GB GPU RAM per
node (12 GB per GPU). The software stack comprises Python 3.6.0, TensorFlow
1.3.1, NumPy 1.16.3, Pandas 0.19.2, and HDF5 1.8.17. The analysis was run on
the Haswell compute nodes of Cori supercomputer at NERSC where each node has
two sockets, and each socket is populated with a 2.3 GHz 16-core processor
(Intel Xeon Processor E5-2698 v3) and 128GB DDR4 2133MHz memory.

To compare the forecasting accuracy of different models, we utilize mean
absolute percentage error (MAPE) and coefficient of determination ($R^{2}$)
computed on the original bandwidth scale (after applying inverse min-max
scalar transformation).

###  V-A Exploratory analysis

Figure 4: Pearson correlation coefficient matrix for the 12 nodes. These nodes
exhibit high correlations and are selected for the exploratory analysis.

We conducted an exploratory analysis on DDCRNN to study the impact of input
horizon duration ($T^{\prime}$), forecasting horizon ($T$), and the type of
autoregression on the forecasting accuracy. Given a pair of sites, we computed
the Pearson correlation coefficient between their time series on the training
data. We selected the 12 out of 96 sites based on the strong correlation
values. The correlation matrix of the 12 sites are shown in Figure 4, where
_in and _out refer to incoming and outgoing nodes, respectively. The reason
for selecting the subset of nodes is to avoid bias in the experimental
comparison. In particular, if we use all the nodes and select the best
options, then it might introduce bias in favor of the DDCRNN while comparing
it with other methods. For the same reason, the forecasting accuracy for the
exploratory analysis was computed on the validation data and not on the test
data, thus avoiding test data leakage.

####  V-A1 Impact of input horizon duration

Input horizon duration | 6hrs | 12hrs | 18hrs | 24hrs | 30hrs | 36hrs | 42hrs | 48hrs
---|---|---|---|---|---|---|---|---
$R^{2}$ ($\mu$) | 0.58 | 0.60 | 0.70 | 0.72 | 0.76 | 0.69 | 0.71 | 0.77
$R^{2}$ ($\sigma$) | 0.02 | 0.09 | 0.10 | 0.12 | 0.07 | 0.11 | 0.05 | 0.08
MAPE ($\mu$) | 22.71 | 22.60 | 20.90 | 19.50 | 19.01 | 22.38 | 20.47 | 20.20
MAPE ($\sigma$) | 6.16 | 4.25 | 6.55 | 1.74 | 4.13 | 8.23 | 3.37 | 3.77
TABLE II: Mean ($\mu$) and standard deviation ($\sigma$) of $R^{2}$ and MAPE
values for varying input horizon duration.

We studied the impact of the length of the input horizon ($T^{\prime}$) on the
forecasting accuracy of the DDCRNN. For training DDCRNN, we varied the length
as 6, 12, 18, 24, 30, 36, 42, and 48 hours to forecast the traffic for the
next 24 hours.

The $R^{2}$ and MAPE values obtained on the validation data are shown in Table
II. For each node, we considered all the observed and their corresponding
predicted values from the DDCRNN on the validation set and computed $R^{2}$
and MAPE. The mean and standard deviation values are computed over 12 nodes.
We can observe that, up to 30 hours, the forecasting accuracy improves as the
sequence length increases. However, after 30 hours, the improvements are not
significant. Therefore, we use an input duration of 30 hours for the rest of
the experiments. Note that a larger input horizons will eventually increase
the training data size and time required for training.

####  V-A2 Impact of forecasting horizon

We analyzed the impact of forecasting horizon on DDCRNN performance. Given 30
hours of input horizon, we trained the DDCRNN model to forecast 6, 12, 18, 24
hours horizon. Note that by training the DDCRNN model for 24 hours, one can
get the forecasting results for all other horizons. However, this strategy may
not be optimal for 6, 12, or 18 hours forecast. Therefore, we trained DDCRNN
for each forecasting horizon.

Table V shows the $R^{2}$ and MAPE values for different forecasting horizons
on the validation data. From the results, as expected, we observe a trend in
which the increase in the forecasting horizon decreases the accuracy. We
select 24 hours as a forecasting horizon for further study, where the DDCRNN
model will be trained to forecast 24 hours. However, we will analyze the
forecasting accuracy for the intermediate time intervals as well.

Forecasting Time Horizon | $R^{2}$ | MAPE
---|---|---
|  $\mu$ ($\sigma$) |  $\mu$ ($\sigma$)
6hrs | 0.92 (0.07) | 11.31 (8.07)
12hrs | 0.82 (0.10) | 16.74 (6.35)
18hrs | 0.72 (0.11) | 20.61 (6.55)
24hrs | 0.76 (0.07) | 19.01 (4.13)
TABLE III: Mean ($\mu$) and standard deviation ($\sigma$) of $R^{2}$ and MAPE
values for varying forecasting horizon.

####  V-A3 Comparison between autoregressive and nonautoregressive DDCRNN

We compared the forecasting accuracy of autoregressive and nonautoregressive
DDCRNN variants. By default, DDCRNN adopts a nonautoregressive approach to
forecasting. An alternative approach is autoregressive forecasting, where
DDCRNN can be trained to forecast only a one-time step, which is then given to
the model recursively to obtain forecasting for 24 hours. This study is
motivated by a previous work [24], where autoregressive forecasting was
adopted within LSTM models for traffic forecasting.

Figure 5: Distribution of R2 values obtained by autoregressive and
nonautoregressive DDCRNN variants for the 24th hour forecast. Figure 6:
Distribution of MAPE values obtained by autoregressive and nonautoregressive
DDCRNN variants for different forecasting intervals.

(a) Autoregressive forecasting strategy

(b) Nonautoregressive forecasting strategy

Figure 7: Comparison between autoregressive and nonautoregressive DDCRNN
variants on ELPA_SUNN_in site.

Figure 5 shows the distribution of $R^{2}$ values obtained by the two
strategies at the 24th hour forecast. Figure 6 shows the distribution of MAPE
values for different forecasting horizon intervals. The result shows that the
default nonautoregressive strategy performs significantly better than the
autoregressive strategy. To take a close look, we plot the forecasting
accuracy ($R^{2}$) of a particular node ELPA_SUNN_in Figure 7. The results
show that the forecasting accuracy of autoregressive DDCRNN decreases rapidly
with an increase in forecasting time steps. This can be attributed to the
error introduced at each time step, which drastically reduces accuracy for the
next time step. This issue is less severe in the nonautoregressive strategy
because it was trained on the entire forecasting horizon.

###  V-B Comparison between DDCRNN and DCRNN

Here, we compare the forecasting accuracy of DDCRNN and DCRNN to show that our
dynamic approach that we proposed is critical for forecasting accuracy. These
two methods differ only with respect to how the weighted adjacency matrix is
used in the diffusion convolution. While in the former it is dynamically
computed based on the current traffic state, in the latter it is static and
based on the connectivity and distance. We also note that DCRNN has been
applied to network traffic forecasting previously in [22]. Both DDCRNN and
DCRNN used the same hyperparameters and the same training data with 96 nodes.
We used 30 hours of input horizon to forecast 24 hours of output horizon. The
forecasting accuracy values are computed on the test data.

Figure 8 shows the $R^{2}$ distribution obtained by DDCRNN and DCRNN for the
24th hour forecast. Figure 9 shows the MAPE distribution for different
forecasting horizon intervals. We can observe that DDCRNN achieves forecasting
accuracy, which is significantly better than DCRNN for all forecasting
horizons. The superior performance of DDCRNN can be attributed to its ability
to model the dynamic spatiotemporal traffic patterns; DCRNN did not have this
capability because it uses the static adjacency matrix, which leads to poor
forecasting accuracy. While a previous work [22] reported superior performance
of the direct application of DCRNN to network traffic, we hypothesize that it
was due to less dynamic traffic data with shorter forecasting horizon.

Figure 8: Distribution of R2 values obtained by DDCRNN and DCRNN for the 24th
hour forecast. Figure 9: Distribution of MAPE values obtained by DDCRNN and
DCRNN for different forecasting intervals.

###  V-C Comparison of DDCRNN with other methods

Here, we compare DDCRNN with other statistical, machine learning, and deep
learning methods including those previously studied for network traffic
forecasting, and show that DDCRNN outperforms all these methods.

We used ARIMA [25] (with order (2,1,2) implemented using statsmodel [26]
Python package), linear regression (LR) [7] (from scikit-learn), Random Forest
(RF) [27] (from scikit-learn with default hyperparameters), Gradient Boosting
(GB) [28] (from scikit-learn with default hyperparameters), Simple RNN (SRNN)
[29] with one hidden layers with 4 neurons (implemented in Keras), LSTM [30],
stacked LSTM (SLSTM): two hidden layers with 256 neurons per layer
(implemented in Keras), Gated Recurrent Unit [31]: one hidden layers with 4
neurons (implemented in Keras), and FireTS [32, 33], a Python package for
multivariate time series forecasting with linear regression as base estimator,
auto-regression order of 24, exogenous order of 24, and exogenous delay of 0.
For each node, we build a node-specific model using each method. Consequently,
for a given method, we have 96 models. For DDCRNN, we have a single model to
forecast the traffic on all 96 nodes. We used 30 hours of input horizon to
forecast 24 hours of output horizon. We computed the accuracy values on the
test data.

Figure 10 shows the distribution of $R^{2}$ obtained by different methods on
the test data for different forecasting intervals (1, 3, 6, 9, 12, and 24
hours). The results show that DDCRNN achieves forecasting accuracy values that
are significantly higher than all the other node-specific models. At 1st hour
forecast, node-specific models accuracy values are closer to DDCRNN, however,
they become poor with an increase in the forecasting horizon. Starting from
the 3rd hour forecast, DDCRNN results are significantly better than the node-
specific models. From 6th hour forecast, we can observe that only DDCRNN has
achieved reasonable accuracy. The superior performance of DDCRNN over node-
specific models can be attributed to the former’s ability to capture both
spatial and temporal patterns using gated recurrent units with diffusion
convolution defined on a dynamic graph.

(a) $1^{st}$ hour forecasting

(b) $3^{rd}$ hour forecasting

(c) $6^{th}$ hour forecasting

(d) $9^{th}$ hour forecasting

(e) $12^{th}$ hours forecasting

(f) $24^{th}$ hour forecasting

Figure 10: Distribution of $R^{2}$ values obtained on 96 nodes by DDCRNN and
other node-specific models for different forecasting horizon intervals.

(a) ACF for ANL-STAR-in

(b) PACF for ANL-STAR-in

(c) ACF for SLAC-SUNN-in

(d) PACF for SLAC-SUNN-in

(e) ACF for JGI-SACR-in

(f) PACF for JGI-SACR-in

(g) ACF for BNL-NEWY-in

(h) PACF for BNL-NEWY-in

(i) ACF for ALBQ-DENV-out

(j) PACF for ALBQ-DENV-out

(k) ACF for CERN513-WASH-in

(l) PACF for CERN513-WASH-in

Figure 11: Autocorrelation (ACF) and partial autocorrelation (PACF) plot
showing correlation of observations as a function of time lag for different
sites.

###  V-D Analysis and characterization

Here, we analyse the forecasting accuracy of DDCRNN and provide insights on
its strengths and limitations.

Table IV shows 10 nodes where DDCRNN obtained best and worst $R^{2}$ values,
respectively. For each node, we computed the autocorrelation and partial
autocorrealtion values with a lag of 10 to analyse to what extent the values
in the time series are correlated. Then, we computed the mean of the obtained
autocorrelation values. Table V shows the mean autocorrelation values for the
top three high and low accuracy nodes. For the same nodes, the autocorrelation
and partial autocorrelation plots are show in Figure 11.

S/n | High accuracy | Low accuracy
---|---|---
1 | ANL_STAR_in | BNL_NEWY_in
2 | ANL_STAR_out | ALBQ_DENV_out
3 | SLAC_SUNN_in | CERN-513_WASH_in
4 | JGI_SACR_in | NASH_WASH_in
5 | FNAL_STAR_out | CERN-513_WASH_out
6 | BNL_NEWY_out | ATLA_ORNL_out
7 | NERSC_SUNN_out | NASH_WASH_out
8 | NERSC_SUNN_in | ATLA_SRS_out
9 | SLAC_SUNN_out | ATLA_NASH_in
10 | LSVN_SUNN_out | ELPA_HOUS_out
TABLE IV: Nodes with high and low forecasting accuracy. Sites | Mean ACF Value
---|---
High accuracy
ANL_STAR_in | 0.783
SLAC_SUNN_in | 0.974
JGI_SACR_in | 0.842
Low accuracy
BNL_NEWY_in | 0.238
ALBQ_DENV_out | 0.526
CERN-513_WASH_in | 0.232
TABLE V: Mean autocorrelation values for the top three nodes with high and low
accuracy.

From the results, we observed that nodes, where DDCRNN obtained high accuracy,
exhibit large autocorrelation values. On the other hand, DDCRNN obtains poor
results when there is no strong temporal correlation, where the traffic is
random. Upon further analysis, we found that nodes with high forecasting
accuracy such as ANL_STAR_in, ANL_STAR_out, and SLAC_SUNN_in have large
experimental facilities. These include Argonne National Laboratory, Stanford
Linear Accelerator, Joint Genome Institute, and Linac Coherent Light Source.
In these user facilities, the data transfers can happen any time, however,
when a transfer happens, a large volume of experimental data will be moved to
data centers for further analysis. These transfers can last several hours and
thus exhibit spatial temporal patterns. DDCRNN can learn these patterns to
provide high forecasting accuracy. Specifically, our analysis reveal that
DDCRNN can model short-term transfers, particularly long connections leading
to predictable bulk data transfers especially driven by these large
facilities. With respect to low accuracy nodes, we see a multitude of factors
that prevented DDCRNN from reaching high accuracy over a long forecasting
horizon. For instance, the CERN facility was shut down most of the year for
maintenance which shows no transfers or very few patterns are picked up to and
from the facility. Additionally, sites such as DENV, BNL, and NASH run
performance tests which consist of very small transfers that take up less
bandwidth. These transfers are so small that the spatial-temporal patterns are
not highly correlated. Therefore, they are not amenable for modeling and
learning with DDCRNN.

##  VI Conclusion and Future work

We developed DDCRNN, a dynamic diffusion convolution recurrent neural network
for forecasting traffic bandwidth in a highly dynamic research wide area
network. These traffic traces show dynamic traffic and lack long term periodic
patterns. Our approach is built on a diffusion convolution recurrent neural
network that models spatial and temporal patterns using graph diffusion
convolution operations within recurrent units. In our approach, the dynamic
traffic is explicitly captured using a state-specific adjacency matrix
computed from the current traffic state in the network.

We evaluated our approach on real traffic traces from ESnet, a US research
network. Our approach outperformed several other alternative methods, first,
showing that the proposed DDCRNN achieves higher forecasting accuracy than
that of the static variant. And secondly, it shows that a single-model can be
used to achieve higher forecasting accuracy than the site-specific models
obtained from linear regression, ARIMA, random forest, gradient boosting, and
simple neural network variants that were previously used for traffic
forecasting.

Achieving higher prediction accuracy presents several advantages to performing
informed routing decisions for potentially new flows that could be impacted
with regular R-WAN behavior. Our DDCRNN approach exposes many potentials to
allow informed flow and routing allocations for network operations, such as
scheduling long-running flow on alternative routes to prevent congestion
points for other smaller flows. In the future, we will couple these decisions
with a controller to test how congestion-free routing will impact the average
utilization of a network in practice.

## Acknowledgements

This work was supported by the U.S. Department of Energy, Office of Science
Early Career Research Program for ‘Large-scale Deep Learning for Intelligent
Networks’ Contract no FP00006145. This material is based upon work supported
by the U.S. Department of Energy (DOE), Office of Science, Office of Advanced
Scientific Computing Research, under Contract DE-AC02-06CH11357. This research
used resources of the Argonne Leadership Computing Facility, which is a DOE
Office of Science User Facility.

## References

  * [1] Cisco, “Trending analysis,” _https://bit.ly/2G2zwwD_ , 2019.
  * [2] R. Boutaba, M. A. Salahuddin, N. Limam, S. Ayoubi, N. Shahriar, F. Estrada-Solano, and O. M. Caicedo, “A comprehensive survey on machine learning for networking: evolution, applications and research opportunities,” _Journal of Internet Services and Applications_ , vol. 9, no. 1, p. 16, 2018. [Online]. Available: https://doi.org/10.1186/s13174-018-0087-2
  * [3] k. claffy, “Internet traffic characterization,” Ph.D. dissertation, UC San Diego, Jun 1994.
  * [4] Q. Lu, L. Zhang, S. Sasidharan, W. Wu, P. DeMar, C. Guok, J. Macauley, I. Monga, S. Yu, J. H. Chen, J. Mambretti, J. Kim, S. Noh, X. Yang, T. Lehman, and G. Liu, “Bigdata express: Toward schedulable, predictable, and high-performance data transfer,” in _2018 IEEE/ACM Innovating the Network for Data-Intensive Science (INDIS)_ , 2018, pp. 75–84.
  * [5] R. Hofstede, P. Čeleda, B. Trammell, I. Drago, R. Sadre, A. Sperotto, and A. Pras, “Flow monitoring explained: From packet capture to data analysis with netflow and ipfix,” _IEEE Communications Surveys Tutorials_ , vol. 16, no. 4, pp. 2037–2064, 2014.
  * [6] W. Yoo and A. Sim, “Time-series forecast modeling on high-bandwidth network measurements,” _J. Grid Comput._ , vol. 14, no. 3, pp. 463–476, Sep. 2016\. [Online]. Available: http://dx.doi.org/10.1007/s10723-016-9368-9
  * [7] K. Lau and Q. Wu, “Local prediction of non-linear time series using support vector regression,” _Pattern recognition_ , vol. 41, no. 5, pp. 1539–1547, 2008.
  * [8] P. Zhao, J. Xia, Y. Dai, and J. He, “Wind speed prediction using support vector regression,” in _2010 5th IEEE Conference on Industrial Electronics and Applications_. IEEE, 2010, pp. 882–886.
  * [9] W. Wang, Y. Bai, C. Yu, Y. Gu, P. Feng, X. Wang, and R. Wang, “A network traffic flow prediction with deep learning approach for large-scale metropolitan area network,” _NOMS_ , 2018.
  * [10] B. Zhou, D. He, and Z. Sun, “Traffic modeling and prediction using arima/garch model,” in _Modeling and Simulation Tools for Emerging Telecommunication Networks_ , 2006, pp. 101–121.
  * [11] Y. Li, R. Yu, C. Shahabi, and Y. Liu, “Diffusion convolutional recurrent neural network: Data-driven traffic forecasting,” in _International Conference on Learning Representations (ICLR ’18)_ , 2018.
  * [12] T. Mallick, P. Balaprakash, E. Rask, and J. Macfarlane, “Graph-partitioning-based diffusion convolution recurrent neural network for large-scale traffic forecasting,” _arXiv preprint arXiv:1909.11197_ , 2019\.
  * [13] R. Wolski, N. T. Spring, and J. Hayes, “The network weather service: A distributed resource performance forecasting service for metacomputing,” _Future Gener. Comput. Syst._ , vol. 15, no. 5-6, pp. 757–768, Oct. 1999\. [Online]. Available: http://dx.doi.org/10.1016/S0167-739X(99)00025-4
  * [14] J. S. Armstrong, “Principles of forecasting: A handbook for researchers and practitioners,” _Kluwer Academic Publishers_ , 2001.
  * [15] H. Wang, H. Xie, L. Qiu, Y. R. Yang, Y. Zhang, and A. Greenberg, “Cope: Traffic engineering in dynamic networks,” _SIGCOMM Comput. Commun. Rev._ , vol. 36, no. 4, pp. 99–110, Aug. 2006. [Online]. Available: http://doi.acm.org/10.1145/1151659.1159926
  * [16] M. B. Ntlangu and A. Baghai-Wadji, “Modeling network traffic using time series analysis: A review,” in _Int. Conf. Big Data and IoT_ , 2017, pp. 209–215. [Online]. Available: http://doi.acm.org/10.1145/3175684.3175725
  * [17] N. McKeown, T. Anderson, H. Balakrishnan, G. Parulkar, L. Peterson, J. Rexford, S. Shenker, and J. Turner, “Openflow: Enabling innovation in campus networks,” _SIGCOMM Comput. Commun. Rev._ , vol. 38, no. 2, pp. 69–74, Mar. 2008. [Online]. Available: http://doi.acm.org/10.1145/1355734.1355746
  * [18] S. Jain, A. Kumar, S. Mandal, J. Ong, L. Poutievski, A. Singh, S. Venkata, J. Wanderer, J. Zhou, M. Zhu, J. Zolla, U. Hölzle, S. Stuart, and A. Vahdat, “B4: Experience with a globally-deployed software wan,” _SIGCOMM Comput. Commun. Rev._ , vol. 43, no. 4, pp. 3–14, 2013. [Online]. Available: http://doi.acm.org/10.1145/2534169.2486019
  * [19] A. Pathak, M. Zhang, Y. C. Hu, R. Mahajan, and D. Maltz, “Latency inflation with mpls-based traffic engineering,” in _Internet Measurement Conf_ , 2011, pp. 463–472.
  * [20] P. Cortez, M. Rio, M. Rocha, and P. Sousa, “Internet traffic forecasting using neural networks,” _IEEE Int. Joint Conf. Neural Network Proceedings_ , 2006\.
  * [21] A. Azzouni and G. Pujolle, “Neutm: A neural network-based framework for traffic matrix prediction in sdn,” _NOMS_ , 2018.
  * [22] D. Andreoletti, S. Troia, F. Musumeci, S. Giordano, G. Maier, and M. Tornatore, “Network traffic prediction based on diffusion convolutional recurrent neural networks,” in _IEEE INFOCOM 2019 - IEEE Conference on Computer Communications Workshops (INFOCOM WKSHPS)_ , 2019, pp. 246–251.
  * [23] S. M. Feit, _SNMP: A guide to network management_. McGraw-Hill Professional, 1993.
  * [24] N. Krishnaswamy, M. Kiran, K. Singh, and B. Mohammed, “Data-driven learning to predict wan network traffic,” in _Proceedings of the 3rd International Workshop on Systems and Network Telemetry and Analytics_ , ser. SNTA ’20. New York, NY, USA: Association for Computing Machinery, 2020, p. 11–18. [Online]. Available: https://doi.org/10.1145/3391812.3396268
  * [25] B. M. Williams and L. A. Hoel, “Modeling and forecasting vehicular traffic flow as a seasonal arima process: Theoretical basis and empirical results,” _Journal of Transportation Engineering_ , vol. 129, no. 6, pp. 664–672, 2003\.
  * [26] S. Seabold and J. Perktold, “statsmodels: Econometric and statistical modeling with python,” in _Python in Science Conference_ , 2010.
  * [27] G. Dudek, “Short-term load forecasting using random forests,” in _Intelligent Systems’ 2014_. Springer, 2015, pp. 821–828.
  * [28] S. B. Taieb and R. J. Hyndman, “A gradient boosting approach to the kaggle load forecasting competition,” _International journal of forecasting_ , vol. 30, no. 2, pp. 382–394, 2014.
  * [29] N. Laptev, J. Yosinski, L. E. Li, and S. Smyl, “Time-series extreme event forecasting with neural networks at uber,” in _International Conference on Machine Learning_ , vol. 34, 2017, pp. 1–5.
  * [30] S. Selvin, R. Vinayakumar, E. Gopalakrishnan, V. K. Menon, and K. Soman, “Stock price prediction using lstm, rnn and cnn-sliding window model,” in _International conference advances in computing, communications and informatics_ , 2017, pp. 1643–1647.
  * [31] R. Fu, Z. Zhang, and L. Li, “Using lstm and gru neural network methods for traffic flow prediction,” in _2016 31st Youth Academic Annual Conference of Chinese Association of Automation (YAC)_. IEEE, 2016, pp. 324–328.
  * [32] J. Xie and Q. Wang, “Benchmark machine learning approaches with classical time series approaches on the blood glucose level prediction challenge.” in _KHD@ IJCAI_ , 2018, pp. 97–102.
  * [33] J. Xie, “firets,” Aug. 2018. [Online]. Available: https://github.com/jxx123/fireTS

The submitted manuscript has been created by UChicago Argonne, LLC, Operator
of Argonne National Laboratory (“Argonne”). Argonne, a U.S. Department of
Energy Office of Science laboratory, is operated under Contract No. DE-
AC02-06CH11357. It is also supported by the Lawrence Berkeley National
Laboratory, a U.S. Department of Energy, Office of Science Early Career
Research Program for ‘Large-scale Deep Learning for Intelligent Networks’
Contract no. FP00006145. The U.S. Government retains for itself, and others
acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in
said article to reproduce, prepare derivative works, distribute copies to the
public, and perform publicly and display publicly, by or on behalf of the
Government. The Department of Energy will provide public access to these
results of federally sponsored research in accordance with the DOE Public
Access Plan. http://energy.gov/downloads/doe-public-access-plan


# Two-photon decay of para-positronium
within a composite approach††thanks: Based on the presentation of M.P. at the
$5^{\text{th}}$ Jagiellonian Symposium on ‘Advances in Particle Physics and
Medicine’, June 29 - July 7 2024, Cracow, Poland

M. Piotrowska1 F. Giacosa1,2
1 Institute of Physics, Jan Kochanowski University,
ul. Uniwersytecka 7, 25-406 Kielce, Poland.
2 Institute for Theoretical Physics, J. W. Goethe University,
Max-von-Laue-Str. 1, 60438 Frankfurt, Germany.

###### Abstract

The decay of the para-positronium into two photons is studied in the framework
of a composite Quantum Field Theoretical approach. This amounts to the
evaluation of the electron-positron dressing together with the Weinberg
compositeness condition for the positronium and the triangle-shaped diagram
with virtual electrons circulating in it, leading to the two-photon final
state. An important role is played by the positronium-electron-positron
vertex, which is linked to the wave function of the para-positronium. We show
how possible choices for the vertex function affect the $\gamma\gamma$ decay
rate. Outlooks to other decay channels and to other positronia are presented.

##  1 Introduction

The positronium is a bound state of an electron-positron pair
($e^{-}$-$e^{+}$), that arises in the framework of Quantum Electrodynamic
(QED). As such, it plays an important role in fundamental physics for being
the lightest ‘element’ [1, 2], as well as in medical physics, mainly due to
its applications in positron emission tomography technique (PET) [1, 3, 4, 5].
The ground state of the positronium, known under the name para-positronium
(p-Ps), corresponds to a spin-singlet state described by the spectroscopic
quantum numbers $n$ ${}^{2S+1}L_{J}=$ $1$ ${}^{1}S_{0}$, or equivalently by
using the relativistic notation as $J^{PC}=0^{-+}$, hence a pseudoscalar state
just as the pion in QCD.

The para-positronium is an unstable short-living state with a mean lifetime
$\sim 0.125$ ns. It decays into photons due to the annihilation of its
components. Charge-conjugation conservation requires that the decay of p-Ps
occurs only into an even number of photons, with the $\gamma\gamma$ mode
(p-Ps$\rightarrow\gamma\gamma$) being by far the dominant one. This decay
width can be expressed by (see e.g. Ref. [6]):

$\Gamma(Ps\rightarrow
n\gamma)=\frac{1}{2J+1}\left|\psi(0)\right|^{2}\lim_{v\to
0}\left[4v\sigma(e^{+}e^{-}\rightarrow n\gamma)\right]\text{ ,}$ (1)

where $v$ and $\sigma$ are relative velocity and annihilation cross-section of
$e^{-}$-$e^{+}$, and where
$\left|\psi(0)\right|^{2}=m_{e}^{2}\alpha^{3}/(8\pi)$ is the spatial wave
function at the origin (annihilation part of the amplitude). Namely, the
ground-state wave function is $\psi(\vec{x})=(\pi a^{3})^{-1/2}e^{-r/a}$,
where $r=|\vec{x}|$ and $a=2a_{0}$ (twice the Bohr radius of the hydrogen
atom). At lowest order it becomes:

$\displaystyle\Gamma\left({}^{1}S_{0}\rightarrow 2\gamma\right)$
$\displaystyle=$
$\displaystyle\frac{1}{2}\frac{e^{4}|\psi(\vec{x}=0)|^{2}}{\pi
m^{4}}\int\limits_{0}^{\infty}|\vec{k}_{1}|^{2}\delta(2m-2|\vec{k}_{1}|)d|\vec{k}_{1}|=$
(2) $\displaystyle=$ $\displaystyle\frac{e^{4}|\psi(\vec{x}=0)|^{2}}{4\pi
m^{2}}=\frac{4\pi\alpha^{2}}{m^{2}}|\psi(\vec{x}=0)|^{2}=\alpha^{5}\frac{m_{e}}{2}\text{
,}$

where the last result is the famous Wheeler-Pirenne formula [7, 8] involving
solely the fine-structure constant $\alpha$ and the electron mass $m_{e}$. In
the framework of QED higher orders can be systematically evaluated, see e.g.
Refs. [9, 10, 11] and Table 1 for a summary.

$\Gamma$(p-Ps $\rightarrow 2\gamma)$
---
Theory
formula | result [$\mu s^{-1}$] | comment | ref.
$\frac{\alpha^{5}m_{e}}{2}$ | $8032.5028(1)$ | lowest order (LO) | [7, 8]
$\Gamma_{0}\left\\{1+\frac{\alpha}{\pi}\left(\frac{\pi^{2}}{4}-5\right)\right\\}$ | $7985.249$ | LO+NLO | [9]
|  | corrections |
$\Gamma_{0}\left\\{1+\frac{\alpha}{\pi}\left(\frac{\pi^{2}}{4}-5\right)-\right.$ | 7989.6178(2) | LO+NLO+NNLO | [10, 11, 19]
$2\alpha^{2}ln\alpha+B_{2\gamma}\left(\frac{\alpha}{\pi}\right)^{2}-$ |  | corrections |
$\frac{3\alpha^{3}}{2\pi}ln^{2}\alpha+C\frac{\alpha^{3}}{\pi}ln\alpha+$ |  |  |
$\left.D\left(\frac{\alpha}{\pi}\right)^{3}\right\\}$ |  |  |
Experiment
result $[\mu s^{-1}]$ | ref.
$7990.9(1.7)$ | [20]
Table 1:  Summary of theoretical and experimental results for the decay rate
p-Ps$\rightarrow 2\gamma$. The constants $A$, $B$, $C$ and $D$ are reported in
[10].

Here, we intend to study the p-Ps within a composite model that makes use of
the compositeness condition, originally proposed to describe the deuteron [12,
13], in a way that resembles the treatment of QCD bound states [14, 15, 16].
To this end, we extend the scalar model described in Ref. [17]. The triangle
diagram of Fig. 1 leading to $\gamma\gamma$ decay is calculated within this
approach. Of course, the aim is not to go beyond the QED precision studies,
but to learn how to deal with the nonperturbative vertex linking the
positronium to its constituents [18]. In the following, a detailed discussion
on this will be provided.

##  2 The composite model

The Lagrangian of our model, describing the decay (and interaction) of para-
positronium (with mass $M_{P}$) into two massless photons, reads:

$\mathcal{L}_{int}=g_{P}P(x)\bar{\psi}(x)i\gamma^{5}\psi(x)-eA_{\mu}(x)\bar{\psi}(x)\gamma^{\mu}\psi(x)\text{
,}$ (3)

where $P(x)$ stands for the pseudoscalar positronium field, $\psi(x)$ is the
electron field, $A_{\mu}$ is the photon field, $e$ is the electric charge of
the proton, and finally $g_{P}$ is the positronium-constituent coupling
constant. Note, there is no direct coupling between p-Ps and photons: the
process p-Ps$\rightarrow\gamma\gamma$ is realized through triangle-shaped
diagram with virtual electrons circulating in it, see Fig. 1.

Figure 1: Triangle-shaped diagram for the process
p-Ps$\rightarrow\gamma\gamma$ .

At this point one needs to specify the kinematics of the two-body decay
illustrated in Fig. 1. For what concerns the external momenta one has:
$p^{\mu}=(M_{P},\vec{0})$, $k_{1}^{\mu}=(\omega,0,0,\omega)$,
$k_{2}^{\mu}=(\omega,0,0,-\omega)$, with $\omega=\frac{M_{P}}{2}$, while the
internal momenta are: $q_{1}=\frac{p}{2}+q$, $q_{2}=\frac{p}{2}-q$ and finally
$q_{3}=\frac{p}{2}+q-k_{1}$.

Let us now introduce one of the most important objects of our model, the
triangle amplitude $I$ related to diagram of the process depicted in Fig. 1,
that reads:

$I=i\int\frac{d^{4}q}{(2\pi)^{4}}\frac{\mathcal{F}(q,p)}{\left(q_{1}^{2}-m_{e}^{2}+i\varepsilon\right)\left(q_{2}^{2}-m_{e}^{2}+i\varepsilon\right)\left(q_{3}^{2}-m_{e}^{2}+i\varepsilon\right)}\text{
,}$ (4)

where $\mathcal{F}(q,p)$ is the ‘nonlocal’ vertex function proportional to the
Fourier transform of the wave function of the para-positronium. Note, this
function is formally not present in Eq. (3), since the latter is local. It
could be however easily included by rendering it nonlocal, see details in
Refs. [14, 15, 16, 21].

The evaluation of the integral of Eq. (4) has been done in two independent
ways. The first one by using a Wick rotation, the second by making use of
residue theorem. The basics of our formalism, including the discussion about
the convergence of integral of Eq. (4), has been presented in the QFT scalar
toy model of Ref. [17].

The triangle amplitude of Eq. (4) can be written as:

$I=i\int\frac{d^{4}q}{(2\pi)^{4}}\frac{\mathcal{F}(q,p)}{D_{1}D_{2}D_{3}}\text{
,}$ (5)

where:

$\displaystyle D_{1,2}$ $\displaystyle=\left(p/2\pm
q\right)^{2}-m_{e}^{2}+i\varepsilon=\left(M_{P}/2\pm
q^{0}\right)^{2}-\vec{q}^{2}-m_{e}^{2}+i\varepsilon,$ (6) $\displaystyle
D_{3}$
$\displaystyle=\left(M_{P}/2+q^{0}-k_{1}^{0}\right)^{2}-(\vec{q}-\vec{k_{1}})^{2}-m_{e}^{2}+i\varepsilon\text{
.}$

By setting $D_{1,2,3}=0$ one gets the corresponding poles:

$\textrm{Poles of
}D_{1}:\left\\{\begin{array}[]{l}L_{1}=-\frac{M_{P}}{2}-\sqrt{\rho^{2}+q_{z}^{2}+m_{e}^{2}}+i\delta\\\
R_{1}=-\frac{M_{P}}{2}+\sqrt{\rho^{2}+q_{z}^{2}+m_{e}^{2}}-i\delta\end{array}\right.\text{
,}$ (7)

$\textrm{Poles of
}D_{2}:\left\\{\begin{array}[]{l}L_{2}=\frac{M_{P}}{2}-\sqrt{\rho^{2}+q_{z}^{2}+m_{e}^{2}}+i\delta\\\
R_{2}=\frac{M_{P}}{2}+\sqrt{\rho^{2}+q_{z}^{2}+m_{e}^{2}}-i\delta\end{array}\right.\text{
,}$ (8)

$\textrm{Poles of
}D_{3}:\left\\{\begin{array}[]{l}L_{3}=-\sqrt{\rho^{2}+(q_{z}-k_{z})^{2}+m_{e}^{2}}+i\delta\\\
R_{3}=\sqrt{\rho^{2}+(q_{z}-k_{z})^{2}+m_{e}^{2}}-i\delta\end{array}\right.\text{
.}$ (9)

The resulting decay width into $\gamma\gamma$ within our approach reads:

$\Gamma_{P\rightarrow\gamma\gamma}=\frac{1}{2}\frac{|\vec{k}_{1}|}{8\pi
M_{P}^{2}}2\left|8m_{e}4\pi\alpha g_{P}I\frac{M_{P}^{2}}{4}\right|^{2}\text{
,}$ (10)

with $|\vec{k}_{1}|=\frac{M_{P}}{2}$.

The para-positronium is not an elementary particle, but an extended object
emerging as a bound state of one electron and one positron. Thus, the
positronium-constituent coupling constant $g_{P}$ entering the Lagrangian in
Eq. (3) is not a free parameter of our model and it can be obtained by using
the Weinberg’s compositeness condition [12, 13]:

$g_{P}=\sqrt{\frac{1}{\Sigma^{\prime}(s=M_{p}^{2})}}\text{ with
}\Sigma(s)=-i\int\frac{d^{4}q}{(2\pi)^{4}}\frac{\mathcal{F}^{2}(\vec{q})^{2}}{D_{1}D_{2}}(-p^{2}/4+q^{2}-m_{e}^{2})\text{
,}$ (11)

where $\Sigma(s)$ is the loop function depicted in Fig. 2 as a self-energy
loop diagram.

Figure 2: Loop diagram of the process p-Ps$\rightarrow e^{-}e^{+}\rightarrow$
p-Ps.

Within our approach the vertex function $\mathcal{F}(q,p)$ turns out to be
proportional to the wave function of the positronium in momentum coordinates,
$A(\vec{q})$:

$\mathcal{F}(q,p)=\mathcal{F}(\vec{q}^{2})\sim
A(\vec{q})=\left(1+\frac{4}{\alpha^{4}+m_{e}^{2}}\vec{q}^{2}\right)^{-2}\text{
.}$ (12)

However, by simply setting it as equal ($\mathcal{F}(q,p)=A(\vec{q})$) does
not work because, for this naive Ansatz, the value of the decay rate of the
process p-Ps$\rightarrow\gamma\gamma$ turns out to be too small by a factor 2
when compared to experimental data.

Moreover, at first sight it seems also that in Eq. (12) covariance is broken.
That is not necessarily the case, since we can interpret the vertex function
as resulting from the Lorentz-invariant object [22]:

$\frac{-(pq)^{2}+p^{2}q^{2}}{p^{2}}=\vec{q}^{2}\text{ ,}$

thus in the rest frame of the decaying particle, in which the four-momentum of
the positronium is $p=(\sqrt{s},\vec{0})$, reduces to

$\mathcal{F}(p,q)=\mathcal{F}\left(\frac{-(pq)^{2}+p^{2}q^{2}}{p^{2}}\right)=\mathcal{F}_{RF}(\vec{q}^{2})\text{
.}$ (13)

As a consequence of this setup, in the rest frame there is no $q^{0}$
dependence (and no additional pole) resulting from the vertex function.

Following Ref. [18] (see also the scalar model of Ref. [17]) let us now
consider another possibility for the vertex function, constructed as follow:

$\mathcal{F}(\vec{q})^{2}=\frac{1}{\left(1+\frac{\vec{q}^{2}}{\gamma^{2}}\right)^{2}}\left(\vec{q}^{2}+\gamma^{2}\right)$
(14)

with $\gamma^{2}=m^{2}-\frac{M_{P}^{2}}{4}$. First, in Fig. 3 we present the
loop function $\Sigma(s=M_{P}^{2})$ and the coupling constant for this choice.
It is visible that for increasing the variable $s$ the value of the loop
function $\Sigma(s)$ also increases and reaches a cusp at threshold. As a
consequence, the coupling constant $g_{P}$ as function of $s=M_{P}^{2}$
decreases for increasing $s$ and vanishes at threshold.

Figure 3:  Left panel: Dependence of loop function $\Sigma$ on the running
square mass of the positronium $s=M_{P}^{2}$. Right panel: dependence of
coupling constant $g_{P}$ on $s=M_{P}^{2}$. The red dots represent the
physical case.

In Table 2 we list the result of the decay width of the process
p-Ps$\rightarrow\gamma\gamma$ obtained in our model by making use of the
vertex function of Eq. (14). It is visible that the result involving the
contribution of all three poles are now closer to the experimental value of
7990.9(1.7) $\mu s^{-1}$ than lowest-order the Wheeler-Pirenne result. In this
respect, it can be shown that in the non-relativistic limit our model
correctly reduces to $\alpha^{5}m_{e}/2$.

Moreover, it is interesting to study the role of each pole, where in Eq. (10)
we split $I=I_{1}+I_{2}+I_{3}$, where $I_{k}$ arises from the pole $R_{k}$.
The contribution to the total decay rate from the first pole is by far the
dominant one. Namely, the ratio of the amplitude contribution of the second
pole w.r.t. the first one is $0.0017087$, while the contribution of the third
w.r.t. first is $-0.00472588$. Interestingly, the third pole gives a negative
contribution to the decay width. This goes in the right direction but the
contribution is even too strong, delivering a theoretical result smaller than
the experiment.

pole(s) contribution | result [$\mu s^{-1}$]
---|---
pole 1 | 8027.83
pole 1 + pole 2 | 8055.29
pole 1+pole 2+pole 3 | 7979.46
Table 2:  Pole contribution to the $\Gamma$ (p-Ps $\rightarrow\gamma\gamma)$.

Note, the result found in Ref. [18] for the same vertex function reads
$7.9527\mu s^{-1}$, the difference being due to the Weinberg compositeness
condition implemented in our approach. Still, we agree in the interpretation
that QFT approaches implicitly contain the resummation of a certain class (but
not all) of QED diagrams, see the discussion in Ref. [23].

##  3 Conclusions

After a brief recall of the theoretical and experimental results on the para-
positronium decay rate into two photons, we introduced a QFT composite model
for its description. Upon evaluating the triangle diagram with virtual
electrons and using the Weinberg compositeness condition, we calculated the
$\gamma\gamma$ decay width. We have shown that the role of the vertex function
is important. The choice of Ref. [18] delivers quite good results, but further
improvement is needed for a determination being in agreement with the very
precise experimental result.

In the future, one should test other choices for the vertex function, such as
the promising covariant Ansatz of Ref. [24]. In general the connection to
bound-state of quarks can also be useful [14, 16, 23, 25]. Also, studies of
the decay of radially excited states of para-positronium into photons as well
as the decay of ortho-positronium into three photons represent an outlook of
our model.

Acknowledgements: This work was supported by the Minister of Science (Poland)
under the ”Regional Excellence Initiative” program (project no.:
RID/SP/0015/2024/01).

## References

  * [1] G. S. Adkins, D. B. Cassidy and J. Pérez-Ríos, Phys. Rept. 975 (2022), 1-61 doi:10.1016/j.physrep.2022.05.002
  * [2] S. D. Bass, Acta Phys. Polon. B 50 (2019), 1319. doi:10.5506/APhysPolB.50.1319
  * [3] S. D. Bass, S. Mariazzi, P. Moskal and E. Stepien, Rev. Mod. Phys. 95 (2023) no.2, 021002. doi:10.1103/RevModPhys.95.021002
  * [4] Moskal, P., Jasinska, B., Stepien, E.Ł. et al., Nat Rev Phys 1, 527–529 (2019). https://doi.org/10.1038/s42254-019-0078-7
  * [5] M. D. Harpen, Med. Phys. 31 (2004), 57-61 doi:10.1118/1.1630494
  * [6] A. Sen and Z. K. Silagadze, Can. J. Phys. 97 (2019) no.7, 693-700 doi:10.1139/cjp-2018-0500
  * [7] J.A. Wheeler, Ann. N.Y. Acad. Sci. 48, 219 (1946).
  * [8] J. Pirenne, Arch. Sci. Phys. Nat. 29, 265 (1947)
  * [9] I. Harris and L.M. Brown, Phys. Rev. 105, 1656 (1957)
  * [10] G. S. Adkins, N. M. McGovern, R. N. Fell and J. Sapirstein, Phys. Rev. A 68 (2003), 032512
  * [11] A. Czarnecki, K. Melnikov and A. Yelkhovsky, Phys. Rev. A 61 (2000), 052502 doi:10.1103/PhysRevA.61.052502
  * [12] S. Weinberg, Phys. Rev. 130 (1963), 776-783. doi:10.1103/PhysRev.130.776.
  * [13] K. Hayashi et al. Fortsch. Phys. 15 (1967) no.10, 625-660. doi:10.1002/prop.19670151002
  * [14] A. Faessler, T. Gutsche, M. A. Ivanov, V. E. Lyubovitskij and P. Wang, Phys. Rev. D 68 (2003), 014011 doi:10.1103/PhysRevD.68.014011
  * [15] F. Giacosa, T. Gutsche and A. Faessler, Phys. Rev. C 71 (2005), 025202 doi:10.1103/PhysRevC.71.025202
  * [16] F. Giacosa, T. Gutsche and V. E. Lyubovitskij, Phys. Rev. D 77 (2008), 034007 doi:10.1103/PhysRevD.77.034007
  * [17] M. Piotrowska and F. Giacosa, Acta Phys. Polon. Supp. 17 (2024) no.1, 1-A7 doi:10.5506/APhysPolBSupp.17.1-A7
  * [18] J. Pestieau, C. Smith and S. Trine, Int. J. Mod. Phys. A 17 (2002), 1355-1398 doi:10.1142/S0217751X02009606
  * [19] B. A. Kniehl and A. A. Penin, Phys. Rev. Lett. 85 (2000), 1210. doi:10.1103/PhysRevLett.85.1210
  * [20] Al-Ramadhan, A. H., and D. Gidley (1994), Phys. Rev. Lett. 72, 1632.
  * [21] T. Wolkanowski, M. Sołtysiak and F. Giacosa, Nucl. Phys. B 909 (2016), 418-428 doi:10.1016/j.nuclphysb.2016.05.025
  * [22] M. Soltysiak and F. Giacosa, Acta Phys. Polon. Supp. 9 (2016), 467-472 doi:10.5506/APhysPolBSupp.9.467
  * [23] C. Smith, Int. J. Mod. Phys. A 19 (2004), 3905-3918 doi:10.1142/S0217751X04020166
  * [24] D. Gromes, Z. Phys. C 57 (1993), 631-638 doi:10.1007/BF01561482
  * [25] Z. P. Li, F. E. Close and T. Barnes, Phys. Rev. D 43 (1991), 2161-2170 doi:10.1103/PhysRevD.43.2161