# A machine learning approach to the Berezinskii-Kosterlitz-Thouless
transition in classical and quantum models

M. Richter-Laskowska${}^{\textsf{\footnotesize{\mbox{?}}}}$, H.
Khan${}^{\textsf{\footnotesize{\mbox{?}}}}$, N.
Trivedi${}^{\textsf{\footnotesize{\mbox{?}}}}$, M.M.
Maśka${}^{\textsf{\footnotesize{\mbox{?}}}}$

(Received May 31, 2018, in final form August 16, 2018)

###### Abstract

The Berezinskii-Kosterlitz-Thouless transition is a very specific phase
transition where all thermodynamic quantities are smooth. Therefore, it is
difficult to determine the critical temperature in a precise way. In this
paper we demonstrate how neural networks can be used to perform this task. In
particular, we study how the accuracy of the transition identification depends
on the way the neural networks are trained. We apply our approach to three
different systems: (i) the classical XY model, (ii) the phase-fermion model,
where classical and quantum degrees of freedom are coupled and (iii) the
quantum XY model.

Key words: phase transitions, topological defects, XY model, artificial neural
networks, machine learning

PACS: 64.60.-i, 05.70.Fh, 07.05.Mh

###### Abstract

Ïåðåõä Áåðåçèíñüêîãî-Êîñòåðëöà-Òàóëåññà äóæå ñïåöèôíèì ôàçîâèì ïåðåõîäîì, ïðè
ÿêîìó âñ òåðìîäèíàìí âåëèèíè íåïåðåðâíèìè. Òîìó âàæêî òîíî âèçíàèòè êðèòèíó
òåìïåðàòóðó. Ó öé ñòàòò íàìè ïîêàçàíî, ÿê ìîæíà âèêîðèñòàòè íåéðîíí ìåðåæ äëÿ
ðîçâ’ÿçàííÿ öüîãî çàâäàííÿ. Çîêðåìà, äîñëäæåíî, äî ÿêî¿ ìðè òîíñòü
ðîçïçíàâàííÿ ïåðåõîäó çàëåæèòü âä ñïîñîáó íàâàííÿ íåéðîííèõ ìåðåæ. Ìè
çàñòîñîâóìî íàø ïäõä äî òðüîõ ðçíèõ ñèñòåì: (i) êëàñèíà XY ìîäåëü, (ii)
ôàçîâî-ôåðìîííà ìîäåëü ç âçàìîäþ ìæ êëàñèíèìè é êâàíòîâèìè ñòóïåíÿìè âëüíîñò
òà (iii) êâàíòîâà XY ìîäåëü.

Ключов слова: ôàçîâ ïåðåõîäè, òîïîëîãí äåôåêòè, XY ìîäåëü, øòóí íåéðîíí ìåðåæ,
ìàøèííå íàâàííÿ

##  1 Introduction

In many cases, thermodynamic phase transitions are clearly visible with well
defined and easily identifiable critical points. In the Landau picture, we
define a macroscopic order parameter that is non-zero in the ordered phase and
vanishes when we cross the critical temperature. Typically, the transition is
signaled by a discontinuity or divergence of some thermodynamic quantities
such as specific heat or magnetic susceptibility. There also exist
unconventional phase transitions which are much more difficult to find. One
example are topological phase transitions, connected to the formation of
topological defects, such as dislocations in two-dimensional crystals,
vortices in two-dimensional superconductors and so on. The proliferation of
these defects leads to the Berezinskii-Kosterlitz-Thouless (BKT) phase
transition [1, 2, 3, 4] where thermodynamic quantities behave smoothly.
Therefore, traditional detection methods based on, e.g., the divergence of the
specific heat or the spin susceptibility, cannot be used. In numerical
approaches, the main difficulty with the identification of the BKT transition
stems from the fact that the interaction between the topological charges
depends logarithmically on the spatial separation. Therefore, numerical
results converge very slowly with the size of the system, and a precise
determination of the critical temperature is a computationally challenging
task. The standard approach is based on the scaling properties of, e.g., the
spin stiffness or superfluid density. This is particularly difficult for
quantum models, where numerical methods are usually very involved and memory-
and time-consuming.

Besides the traditional methods which rely on the idea of the order parameter
or quantities accessible in experiments, such as the specific heat or magnetic
susceptibility, attempts to use artificial neutral networks to identify phases
in condensed matter and transition between them have been recently undertaken
[5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. These new methods
have proven to be accurate and reliable for classical Ising-type models [5, 6,
7]. Less spectacular progress has been made for quantum systems [15, 16] and
for systems with topological phase transitions [17, 18, 19, 20]. In the latter
case, the main difficulty comes from the fact that the topologically ordered
phase is not described by a local order parameter [1, 2, 3, 4]. Instead, its
formation is connected with suppression of non-local topological defects which
are difficult to identify.

In this paper, we demonstrate the application of machine learning approaches
to identify topological transitions in a few different types of two-
dimensional classical and quantum systems. In particular, we study the
classical XY (c-XY) model, the phase-fermion (PF) [21] where the interaction
is only between quantum and classical degrees of freedom, and the fully
quantum XY (q-XY) model.

The first of these models has already been thoroughly analyzed in [17]. The
authors have shown there that treating spin configurations as raw images in
the case of a feed-forward network does not lead to the correct value of the
critical temperature. Instead, they propose to preprocess the spin
configurations into vorticity and then use the results to train two kinds of
artificial neural networks (ANN): a one-layer feed-forward network and a deep
convolutional network. In both cases, the results are scaled with the system
size towards the correct value of the critical temperature, but the one-layer
network performed poorly for large systems. In the present approach, we do not
have convolutional layers, but we use a deep feed-forward network composed of
four fully-connected layers (we learned that the choice of particular meta-
parameters is not crucial for the network performance). What is important is
that, instead of using raw spin configurations where each spin is represented
by a number from 0 to $2\piup$ (which gives rather inaccurate results, as
demonstrated in [17]), we represent the configurations as vectors of sines and
cosines of the spin angles, which reflects the system’s symmetry.

##  2 Models

The first model describes classical spins of unit length in a square lattice
with a nearest neighbor interaction given by the following Hamiltonian

$H_{\text{c-XY}}=-J\sum_{\langle
i,j\rangle}\cos\left(\theta_{i}-\theta_{j}\right),$ (2.1)

where $J$ is the coupling between the spins, and $\theta_{i}$ describes the
direction of spin $i$. It is known that this model exhibits the BKT phase
transition, and precise finite-size scaling gives the critical temperature
$T_{\rm BKT}\approx 0.8935$ in units of $J$ [22].

In the next model, classical spins $\theta_{i}$ interact with fermions. The
Hamiltonian of the PF model is given by

$H_{\rm PF}=-t\sum_{\langle
i,j\rangle,\sigma}\hat{c}^{\dagger}_{i\sigma}\hat{c}_{j\sigma}+g\sum_{i}\left(\mathrm{e}^{\mathrm{i}\theta_{i}}\hat{c}_{i\uparrow}\hat{c}_{i\downarrow}+{\rm
h.c.}\right),$ (2.2)

where $\hat{c}^{\dagger}_{i\sigma}\ \left(\hat{c}_{i\sigma}\right)$ is an
operator that creates (annihilates) spin-$\sigma$ electron at the lattice site
$i$, and $g$ describes the strength of the interaction between the classical
and quantum degrees of freedom. We set the hopping integral $t$ as the energy
unit ($t=1$). In this model, the fermions mediate an effective interaction
between the classical spins $\theta_{i}$ which also leads to the BKT phase
transition. The critical temperature is a function of $g$, and for $g=2$ Monte
Carlo (MC) simulations give $T_{\rm BKT}\approx 0.12$ [21]. The PF model can
be treated as an approximation of the boson-fermion model [23], valid when
fluctuations of the number of bosons at a lattice site can be neglected.

The Hamiltonian of the last model, the quantum XY model, is given by

$H_{\text{q-XY}}=\frac{E_{\text{c}}}{2}\sum_{i}\hat{n}_{i}^{2}-J\sum_{\langle
i,j\rangle}\cos\left(\hat{\theta}_{i}-\hat{\theta}_{j}\right),$ (2.3)

where $\hat{n}_{i}$ is the number operator that is canonically conjugate to
the quantum phase operator $\hat{\theta}_{i}$, and $E_{\text{c}}$ is the
charging energy.

To train neural networks and to classify phases, one needs extensive sets of
spin configurations generated at different temperatures. They were produced
with the help of the MC simulations. In the case of the classical XY model, we
directly used the Metropolis algorithm. For the PF model, we also used the
Metropolis algorithm, but in each MC step (i.e., for each generated spin
configuration) we needed to diagonalize the fermionic Hamiltonian (2.2) [24].
For the quantum XY model, we used a quantum MC method in which the Hamiltonian
was mapped to a classical action of spins on an effective (2+1)D lattice. We
then used a Wolff cluster algorithm to simulate these spins. In all cases, the
simulations were performed on $16\times 16$ systems.

In all these models, the helicity modulus $\Upsilon$, at finite temperature
defined as the second derivative of the free energy with respect to an
externally imposed global twist across the sample [25], has a universal jump
at the BKT transition. In the thermodynamic limit, it drops at $T_{\rm BKT}$
from $\frac{2}{\piup}T_{\rm BKT}$ to zero. However, in finite systems it
evolves smoothly and converges only logarithmically to the thermodynamic
limit, as shown in figures 1 (a), 1 (c), and 1 (e).

Figure 1: (Colour online) Helicity modulus for the c-XY model (a), PF model
for $g=4t$ (c), and q-XY model for $E_{\text{c}}=0.1J$ (e). The corresponding
critical temperatures extrapolated to the thermodynamic limit are presented in
panels (b), (d), and (f). Since $T^{*}(L)-T_{\rm BKT}\propto\log(L)^{-2}$,
where $T^{*}(L)$ is the temperature at which $\Upsilon$ drops the most rapidly
in a $L\times L$ system [27, 28, 22], the temperature is presented in these
panels as a function of $\log(L)^{-2}$.

A finite size scaling of the temperature at which $\Upsilon(T)$ crosses
$\frac{2}{\piup}T$ can give a rough estimate of $T_{\rm BKT}$. This is
demonstrated in figures 1 (b), 1 (d), and 1 (f). It is also possible to
determine $T_{\rm BKT}$ in a more precise way. At the critical temperature,
the helicity modulus scales with the system size according to the Kosterlitz
renormalization group (RG) equations [26],

$\Upsilon_{L}=\Upsilon_{L\to\infty}\left[1+\frac{1}{2}\frac{1}{\ln(L)+C}\right],$
(2.4)

where $L$ is a linear size of the system and $C$ is a constant. Therefore, if
one fits $\Upsilon_{L}(T)$ calculated in MC simulations for different system
sizes $L$ to the RG predictions, the fitting errors drop almost to zero
exactly at $T_{\rm BKT}$ [29]. In figure 2 one can see the procedure.

Figure 2: (Colour online) The root-mean-square error $\delta$ for fitting MC
results for the c-XY model (a), the PF model (b), and the q-XY model (c) to
the RG predictions given by equation (2.4). The sharp minimum of the fitting
errors indicates $T_{\rm BKT}$. The model parameters are the same as in figure
1.

##  3 Artificial neural network

Here, however, we study how the critical temperature can be determined with
the help of ANN trained to identify the low- and high-temperature phases. When
a system approaches the critical temperature, thermal fluctuations increase
drastically so that the spin configurations can be very different from
configurations generated at very low temperatures. However, as long as $T$ is
below $T_{\rm BKT}$, the system is still in the same low-temperature phase.
The problem with systems where the BKT transition takes place, such as the
c-XY, PF or q-XY models, is that in the low-temperature phase, MC simulations
on finite clusters show finite magnetization, which according to the Mermin-
Wagner theorem cannot exist in the thermodynamic limit. Since the MC results
are used to train the ANN, it is possible that the network will learn finite-
size features of the spin configurations, such as magnetization. It was
proposed in [17] that this difficulty can be overcome by adding a
convolutional layer designed to identify topological defects in spin
configurations. Then, the network learns configurations of vortices and
antivortices instead of raw spin configurations.

Our aim here is to show how the features that can be used to identify phases
depend on the distance from the critical temperature $T_{\rm BKT}$. The
standard ANN-based approach to the problem of finding a phase transition is to
train the ANN to recognize features of the low- and high-temperature phases.
Within the scheme of supervised learning, one needs to feed the ANN with
labelled spin configurations generated at low and high temperatures. The
network is supposed to extract characteristic patterns and to learn how to
classify configurations which were not used at the training stage. Usually,
this is an easy task for configurations generated at very low or at very high
temperatures, since they are clearly distinctive. However, to precisely
determine the critical temperature, the ANN must distinguish configurations
generated slightly below $T_{\rm BKT}$ and slightly above $T_{\rm BKT}$.
Thermal fluctuations in this regime make those configurations very different
from the fully ordered low-temperature configurations and from completely
random high-temperature configurations. The question then is, whether an ANN
trained at extreme temperatures will be capable of classifying the phases
close to $T_{\rm BKT}$? In other words, do the distinctive features learned by
the ANN at the extremes balance each other out just right in the evaluation
process such that the correct $T_{\rm BKT}$ is predicted? To answer this
question we trained the ANN at different distances from $T_{\rm BKT}$ and
checked how the distance $|T-T_{\rm BKT}|$ affects the accuracy of finding
$T_{\rm BKT}$. To be precise, we used MC simulations to generate sets of
configurations $\\{C_{1}\\},\,\\{C_{2}\\},\,\ldots,\\{C_{N}\\}$ at
temperatures $T_{1},\,T_{2},\,\ldots,\,T_{N}$ ($T_{1}<T_{2}<\ldots<T_{N}$)
below and above $T_{\rm BKT}$. For the c-XY, PF, and q-XY models, the ranges
of temperatures were from $0.1J$ to $1.6J$, from $0.02t$ to $0.2t$, and from
$0.1J$ to $1.5J$, respectively. Then, we generated sets of configurations
representing different low- (${\cal L}_{m}$) and high-temperature (${\cal
H}_{m}$) ranges:

${\cal L}_{m}=\bigcup\limits_{i=1}^{m}C_{i}\hskip 14.22636pt\mbox{and}\hskip
14.22636pt{\cal H}_{m}=\bigcup\limits_{i=N-m+1}^{N}C_{i}\,,$ (3.1)

where $m<N/2$ is the number of temperatures in each range. We randomly removed
from ${\cal L}_{m}$ and ${\cal H}_{m}$ some number of configurations to keep
the cardinality of these sets fixed. In this way, our study, on how the
critical temperature predicted by the ANN depends on the range of temperatures
used to train it, was not affected by a different number of configurations for
different ranges. In order to connect $m$ with a temperature range used in the
trainings, we introduce $\tau$ as a measure of the relative temperature range:

$\tau=\frac{T_{\rm high}-T_{\rm low}}{T_{\rm BKT}-T_{\rm low}}\,,$ (3.2)

where $T_{\rm low}$ and $T_{\rm high}$ are the lowest and the highest
temperatures lower than $T_{\rm BKT}$ which were used to train the ANN.

Configurations from ${\cal L}_{m}$ and ${\cal H}_{m}$ were mixed together and
shuffled and then they were used to teach the ANN to identify the low- and
high-temperature phases.

The spin configurations generated in MC are stored as numbers $\theta_{i}$
from 0 to $2\piup$. However, in order to take into account the character of
classical two-dimensional spins, the configurations were rewritten as an array
composed of cosines and sines:
$\cos\theta_{1},\ldots,\cos\theta_{N},\sin\theta_{1},\ldots,\sin\theta_{N}$,
where $N=L\times L$ is the number of lattice sites. This is equivalent to a
representation by complex numbers and has the advantage that almost parallel
spins are represented by close numbers, which is not the case for the original
representation by the angles $\theta_{i}$.

Since the generation of extensive sets of spin configurations for different
system sizes and at different temperatures is a time-consuming task,
especially for the PF and q-XY models, we used a technique known from machine
learning-based image recognition to increase the number of configurations
without additional MC simulations. Namely, we transformed the original
configurations according to symmetries of the system. We used periodic
boundary conditions, which allowed us to apply translations to create new
configurations. Other transformations were reflections and rotations. Figure 3
schematically shows how the ANN is used to classify spin configurations.

Figure 3: (Colour online) A scheme of a binary classification process in a
feedforward ANN: a given spin configuration $\\{\theta_{i}\\},\ i=1,\ldots,N$,
is rewritten as a length-$2N$ vector
$[\cos\theta_{1},\ldots,\cos\theta_{N},\sin\theta_{1},\ldots,\sin\theta_{N}]$
and is presented to the input layer of the ANN (green circles). Then,
activations flow across fully connected hidden layers (red circles) to the
output layer composed of only one neuron (violet circle).

The ANN was implemented using the KERAS package with TENSORFLOW as the
computational backend. We used a deep feedforward network with four fully
connected hidden layers with 512, 192, 64, 16, and 16 neurons. Such a
structure was a balance between the number of training epochs required for
convergence and the time needed for one epoch. It turns out, however, that the
metaparameters are not crucial for the ANN performance. As the activation
function, a rectified linear unit (ReLU) was used in the hidden layers and
sigmoid function in the output layer. The network was trained to minimize the
distance between the MC data and the model predictions defined by the binary
cross entropy. The loss function is given by

$L=-\sum_{i}\left[y_{i}\log p_{i}+(1-y_{i})\log(1-p_{i})\right],$ (3.3)

where $y_{i}$ are labels and $p_{i}$ are the corresponding predictions.

We found that with the Tikhonov regularization (L2) [31], the network performs
better in the classification of the low- and high-temperature phases. Figure 4
shows an example of learned weights used to classify phases of a $16\times 16$
system. One can see that despite a rather large size of the network, most of
the neurons are activated.

Figure 4: (Colour online) Visualization of the activations in a deep ANN
composed of hidden layers of 512, 192, 64, 16, and 16 neurons. White colour
corresponds to zero weight, blue to negative weight and red to positive
weight. “Weights 1” connect the output of layer 1 to layer 2, so they are
represented by a rectangle with $512\times 192$ colour squares. Similarly,
“weights 2” are represented by a $192\times 64$ rectangle, etc. The network
was initiated with all weights close to zero (we assumed a finite value
$|w_{i}|<10^{-3}$ to avoid the vanishing gradient problem [30]), so blue and
red parts indicate neurons which were activated during training.

##  4 Results

Figure 5: (Colour online) Calculated by the ANN, the probability that a given
configuration belongs to the high-temperature phase of the c-XY model. The
network was trained on labelled data generated according to equation (3.1) for
$m=1$ (a), $m=10$ (b), and $m=16$ (c). Then, the network determined the
probabilities 100 times for configurations generated at different temperatures
each time starting from different weights and biases. The vertical error bars
show the standard deviation. The solid red line is the best fit
$P(T)=0.5\tanh\left[\alpha\left(T-T_{\rm BKT}\right)\right]+0.5$, where
$\alpha$ and $T_{\rm BKT}$ are fitting parameters. The dashed red line shows
$1-P(T)$ which is the probability of being classified as a low-temperature
phase. The black arrows indicate the temperatures (a) or the ranges of
temperatures (b), (c) used to train the network. The black vertical line
indicates $T_{\rm BKT}$ determined from the RG equations. The dashed green
line in panel (c) shows magnetization. Comparing panels (a)–(c) one can see
that the for the c-XY model, the average critical temperature does not change
significantly with an increasing range of training temperatures, but the
spread of the results shrinks. The temperatures at which the ANN was trained
and tested were (in units of $J$): from 0.10 to 0.70 and from 1.10 to 1.60
with stepsize 0.05 and from 0.750 to 1.050 with stepsize 0.025. Panel (d)
shows estimated $T_{\rm BKT}$ as a function of $\tau$. The horizontal red line
shows the critical temperature determined from fitting the MC results to the
RG equation (2.4). The right-hand vertical axis shows $(T_{\rm BKT}-T^{0}_{\rm
BKT})/T^{0}_{\rm BKT}\times 100\%$, where $T_{\rm BKT}$ is the average ANN
prediction and $T^{0}_{\rm BKT}$ is the actual critical temperature.

Figure 5 shows how the predictions for $T_{\rm BKT}$ of the c-XY model
calculated by the ANN depend on the way the network was trained. In each case,
we trained the neural network using a 10-fold cross validation technique [32]
and repeated this procedure 10 times. As a result, we obtained 100 possible
values of probability $P$ that a given configuration belongs to the high-
temperature phase (the probability that it belongs to the low-temperature
phase is $1-P$). Each time, the ANN was initialized with random weights and
biases, and a larger spread of the predictions indicates a greater difficulty
in an unambiguous classification of the phase. The same method of multiple
trainings starting from different random weights and biases was used in [13]
to determine the standard error of the predicted critical temperature. One can
see in figure 5 that even if the ANN was trained only at the extreme
temperatures ($T=0.1$ and $T=1.6$, $m=1$) corresponding to a fully ordered and
completely random configurations, the average predicted $T_{\rm BKT}$ is not
far from the actual value. This could be an accidental coincidence because
with an increasing $m$ the deviation slightly increases, but always remains
below 10%, which can be seen in figure 5 (d). For $m=16$, the deviation is
smaller than the line width. The problem, however, is that in this case, the
uncertainty is large. So, for a precise determination of $T_{\rm BKT}$, the
averaging over a large number of configurations generated at different
temperatures is necessary. For example, for $m=1$, the average prediction for
1000 statistically independent configurations is less than 1% off the actual
$T_{\rm BKT}$, but individual predictions are spread within the range of $\pm
18\%$ around this value. This means that if one saves the computational time
required to generate configurations for training, more configurations should
be generated for the evaluation stage. With an increasing width of the
temperature range used for training, the spread of the calculated
probabilities decreases significantly. Figure 5 (d) shows how the average
$T_{\rm BKT}$ depends on the number of different temperatures used in
training. On the right-hand vertical axis, showing the relative error, one can
see that the error is always below 10%, and with an increasing number $m$, it
converges to the actual value of $T_{\rm BKT}$ for the c-XY model.

Figure 6: (Colour online) The same as in figure 5, but for the PF model. In
this case, the critical temperature is affected by the width of the range
training temperatures and a relatively wide range is necessary to obtain
precise $T_{\rm BKT}$ [c.f. figure 6 (b)]. The temperatures at which the ANN
was trained and tested were (in units of $t$): from 0.02 to 0.20 with stepsize
0.01.

In the case of the PF model, the results are different. As can be seen in
figure 6, the spread of the calculated probabilities is less dependent on $m$,
but the average critical temperature strongly depends on $m$. This means that
for the PF model, an increase of the number of configurations used at the
stage of phase classification will not guarantee a more precise estimation of
$T_{\rm BKT}$. Instead, for this model, a sufficiently wide range of
temperatures at which configurations for training are generated is necessary.
From the physical point of view, this means that in the PF model, extremely
low-temperature configurations and extremely high-temperature configurations
are more different from the configurations close to the critical point than in
the case of the c-XY model. It can be seen in figure 6 (d) that for the PF
model, the relative error for $m=1$ extends to more than 20%.

Figure 7: (Colour online) The same as in figure 5, but for the q-XY model.
Similarly to the c-XY model, also here the critical temperature is rather
insensitive to the width of the training temperatures, but the spread of the
results decreases with an increasing $m$. The temperatures at which the ANN
was trained and tested were (in units of $J$): from 0.1 to 0.8 and from 1.2 to
1.5 with stepsize 0.1 and from 0.85 to 1.15 with stepsize 0.05.

The distribution of probabilities calculated for the q-XY model is similar to
that for its classical counterpart. It is presented in figure 7. Though for
$m=1$, the estimated critical temperature differs from its real value by 14%
[see figure 7 (d)], the difference decreases very quickly with an increasing
$m$ and already for $m\geqslant 3$, the relative error is around 2%.

The results show that in the case of the PF model, a much richer set of
configurations is required to properly train the ANN than for the XY models.
The reason can be connected to a different character of this model. In both
the classical and quantum XY models, the interaction range is limited to
nearest neighbors. On the other hand, fermions in the PF model mediate the
effective interactions between arbitrarily spaced lattice sites. This effect
is seen in figure 1: at low temperature the helicity modulus in the c-XY and
q-XY models converges even for very small systems. This is not the case for
the PF model, where even at very low temperatures (i.e., in an almost fully
ordered state) the energy per lattice site depends on the system size. This
results from the delocalization of fermions which in the ordered state are
similar to quantum particles in an infinite quantum well, with their energies
strongly dependent on the size of the well.

##  5 Summary

We have demonstrated how the accuracy of finding the BKT transition in three
different models depends on the range of temperatures at which the ANN was
trained. We used a simple feedforward network with densely connected hidden
layers. We did not perform any feature engineering of the spin configurations
generated in MC simulations and we did not use convolutional layers.
Therefore, the phase classification was based on raw spin configurations, not
on the explicitly extracted vortices as in [17]. Nevertheless, in figure 5 (c)
we compare the calculated probabilities and magnetization that results from
the finite size of the system. One can see there that the section of $P(T)$
that indicates the BKT transition is much steeper than the temperature
dependence of the magnetization, even if the network was trained on the
extreme temperatures [figure 5 (a)]. Therefore, we believe that the ANN learns
not only the magnetization (which would vanish in the thermodynamic limit),
but also some topological features connected with the BKT transition. One also
cannot exclude that the ANN is capable of learning the character of the spin-
spin correlations which change their behavior at the BKT transition. To
confirm this, however, at least a finite size analysis of the ML results would
be necessary, which has not been performed here. However, our aim was
different — we wanted to demonstrate how the critical temperature determined
by the ANN depends on the composition of the training set. As one can expect,
the larger is the variety of the configurations representing the low-
temperature and high-temperature phases, the better is the accuracy of the
critical temperature. We found, however, that for the c-XY and q-XY model, the
average $T_{\rm BKT}$ was close to the actual one even if the ANN was trained
relatively far from the critical point. Increasing the range of temperatures
at which the network was trained, only slightly improves the numerical
accuracy (i.e., the difference between the average $T_{\rm BKT}$ determined by
the ANN and the value found from the RG equations), but significantly reduces
the uncertainty. The situation is different for the PF model, where the
numerical accuracy is strongly dependent on the temperature range used at the
training stage. We attribute this behavior to the long-range effective
interaction present in the PF model which can lead to a longer range of the
spin-spin correlations and their different temperature dependence.

Despite the difference in the results for the XY models and the PF model, in
all cases it is important to train the ANN not only at very low and at very
high temperatures, but also as close as possible to the critical temperature.
The main problem in training at temperatures close to $T_{\rm BKT}$ is that
for supervised learning, the configurations must be labelled, so one should
know at least an approximate value of the critical temperature. One of the
ways to overcome this difficulty is to use the learning by confusion approach
[14] based on a combination of supervised and unsupervised techniques.

## Acknowledgements

M.M.M. acknowledges support by NCN (Poland) under grant 2016/23/B/ST3/00647.
H.K. and N.T. acknowledge funding from grant no. NSF DMR 1629382.

## References

  * [1] Berezinskii V.L., Sov. Phys. JETP, 1971, 32, 493.
  * [2] Berezinskii V.L., Sov. Phys. JETP, 1972, 34, 610.
  * [3] Kosterlitz J.M., Thouless D.J., J. Phys. C: Solid State Phys., 1972, 5, L124, doi:10.1088/0022-3719/5/11/002.
  * [4] Kosterlitz J.M., Thouless D.J., J. Phys. C: Solid State Phys., 1973, 6, 1181, doi:10.1088/0022-3719/6/7/010.
  * [5] Carrasquilla J., Melko R.G., Nat. Phys., 2017, 13, 431, doi:10.1038/nphys4035.
  * [6] Morningstar A., Melko R., J. Mach. Learn. Res., 2018, 18, 163.
  * [7] Ponte P., Melko R., Phys. Rev. B, 2017, 96, 205146, doi:10.1103/PhysRevB.96.205146.
  * [8] Zhang Y., Kim E.-A., Phys. Rev. Lett., 2017, 118, 216401, doi:10.1103/PhysRevLett.118.216401.
  * [9] Wang L., Phys. Rev. B, 2016, 94, 195105, doi:10.1103/PhysRevB.94.195105.
  * [10] Hu W., Singh R.R.P., Scalettar R.T., Phys. Rev. E, 2017, 95, 062122, doi:10.1103/PhysRevE.95.062122.
  * [11] Wetzel S.J., Phys. Rev. E, 2017, 96, 022140, doi:10.1103/PhysRevE.96.022140.
  * [12] Broecker P., Carrasquilla J., Melko R.G., Trebst S., Sci. Rep., 2017, 7, 8823,
doi:10.1038/s41598-017-09098-0.

  * [13] Ch’ng K., Carrasquilla J., Melko R.G., Khatami E., Phys. Rev. X, 2017, 7, 031038,
doi:10.1103/PhysRevX.7.031038.

  * [14] Van Nieuwenburg E.P.L., Liu Y.-H., Huber S.D., Nat. Phys., 2017, 13, 435, doi:10.1038/nphys4037.
  * [15] Torlai G., Mazzola G., Carrasquilla J., Troyer M., Melko R., Carleo G., Nat. Phys., 2018, 14, 447, doi:10.1038/s41567-018-0048-5.
  * [16] Carleo G., Troyer M., Science, 2017, 355, 602, doi:10.1126/science.aag2302.
  * [17] Beach M.J.S., Golubeva A., Melko R.G., Phys. Rev. B, 2018, 97, 045207, doi:10.1103/PhysRevB.97.045207.
  * [18] Deng D.-L., Li X., Sarma S.D., Phys. Rev. B, 2017, 96, 195145, doi:10.1103/PhysRevB.96.195145.
  * [19] Zhang W., Liu J., Wei T.-C., Preprint arXiv:1804.02709, 2018.
  * [20] Rodriguez-Nieva J.F., Scheurer M.S., Preprint arXiv:1805.05961, 2018.
  * [21] Maśka M.M., Trivedi N., Preprint arXiv:1706.04197, 2017.
  * [22] Hsieh Y.-D., Kao Y.-J., Sandvik A.W., J. Stat. Mech.: Theory Exp., 2013, 2013, P09001,
doi:10.1088/1742-5468/2013/09/P09001.

  * [23] Micnas R., Ranninger J., Robaszkiewicz S., Rev. Mod. Phys., 1990, 62, 113, doi:10.1103/RevModPhys.62.113.
  * [24] Maśka M.M., Czajka K., Phys. Rev. B, 2006, 74, 035109, doi:10.1103/PhysRevB.74.035109.
  * [25] Fisher M.E., Barber M.N., Jasnow D., Phys. Rev. A, 1973, 8, 1111, doi:10.1103/PhysRevA.8.1111.
  * [26] Kosterlitz J.M., J. Phys. C: Solid State Phys., 1974, 7, 1046, doi:10.1088/0022-3719/7/6/005.
  * [27] Schultka N., Manousakis E., Phys. Rev. B, 1994, 49, 12071, doi:10.1103/PhysRevB.49.12071.
  * [28] Tomita Y., Okabe Y., Phys. Rev. B, 2002, 65, 184405, doi:10.1103/PhysRevB.65.184405.
  * [29] Weber H., Minnhagen P., Phys. Rev. B, 1988, 37, 5986(R), doi:10.1103/PhysRevB.37.5986.
  * [30] Sussillo D., Abbott L.F., Preprint arXiv:1412.6558v3, 2015.
  * [31] Ng A.Y., In: Proceedings of the Twenty-First International Conference on Machine Learning (Banff, Canada, 2004), ACM, New York, 2004, 78, doi:10.1145/1015330.1015435.
  * [32] Gunasegaran T., Cheah Y.-N., In: Proceedings of the 8th International Conference on Information Technology (Amman, Jordan, 2017), IEEE, 2017, 89–95, doi:10.1109/ICITECH.2017.8079960.

Ukrainian 0

Çàñòîñóâàííÿ ìàøèííîãî íàâàííÿ äî ïåðåõîäó Áåðåçèíñüêîãî-Êîñòåðëöà-Òàóëåññà â
êëàñèíèõ êâàíòîâèõ ìîäåëÿõ Ì. Ðõòåð-
Ëÿñêîâñüêà${}^{\textsf{\footnotesize{\mbox{?}}}}$, Ã.
Õàí${}^{\textsf{\footnotesize{\mbox{?}}}}$, Í.
Òðâåä${}^{\textsf{\footnotesize{\mbox{?}}}}$, Ì.Ì.
Ìàñüêà${}^{\textsf{\footnotesize{\mbox{?}}}}$

${}^{\textsf{\footnotesize 1}}$ íñòèòóò ôçèêè, Ñëåçüêèé óíâåðñèòåò, âóë. 75-ãî
ïõîòíîãî ïîëêó, 1, 41-500 Õîæâ, Ïîëüùà ${}^{\textsf{\footnotesize 2}}$
Ôàêóëüòåò ôçèêè, óíâåðñèòåò øòàòó Îãàéî, ïðîñï. Â. Âóäðàôôà, 191, Êîëóìáóñ,
Îãàéî 43210, ÑØÀ


# Jet Quenching of the Heavy Quarks in the Quark-Gluon Plasma
and the Nonadditive Statistics

Trambak Bhattacharyya<EMAIL_ADDRESS>Bogoliubov Laboratory
of Theoretical Physics, Joint Institute for Nuclear Research, Dubna, 141980,
Moscow Region, Russia  Eugenio Megías<EMAIL_ADDRESS>Departamento de Física
Atómica, Molecular y Nuclear and Instituto Carlos I de Física Teórica y
Computacional, Universidad de Granada, Avenida de Fuente Nueva s/n, 18071
Granada, Spain  Airton Deppman<EMAIL_ADDRESS>Instituto de Física,
Universidade de São Paulo, Rua do Matão 1371, São Paulo 05508-090, Brazil

###### Abstract

Using the Plastino-Plastino (PP) equation, we calculate transport coefficients
of the heavy-quarks traversing inside the quark-gluon plasma, and generalize
their relationship with differential energy loss. The PP equation indicates
anomalous diffusion of the probe particles and yields a quasi-exponential
stationary distribution obtained also from the nonadditive statistics proposed
by C. Tsallis. We estimate energy loss in a nonadditive quark-gluon medium,
and calculate the jet-quenching parameter ($\hat{q}$) for the PP dynamics.
With the help of the estimate of $\hat{q}$, we calculate the nuclear
suppression factor ($R_{\text{AA}}$) of the heavy-quarks passing through a
nonadditive quark-gluon plasma using the model proposed by Dokshitzer and
Kharzeev. In many a case, the parameters in the analysis are fixed from the
experimental results to minimize arbitrariness. There is a good agreement
between the theoretical calculation and experimental $R_{\text{AA}}$ data,
indicating that fast heavy-quarks may be subjected to anomalous diffusion
inside the QGP.

##  I Introduction

Transport of high-energy particles inside the hot and dense quark-gluon plasma
(QGP) medium created in high-energy collisions (e.g. Pb-Pb) is an interesting,
and active field of research. One of the standard ways to study this
phenomenon is to utilize the Boltzmann transport equation (BTE), or its
approximations (e.g. the Fokker-Planck equation), that dictate the evolution
of a non-equilibrium distribution (say that of the heavy quarks) inside the
QGP medium of the light quarks (LQ), and gluons (g). Important inputs to these
transport equations are the transport coefficients like drag, and diffusion
coefficients that are also related to the energy loss inside the medium.
Transport equations are solved using these inputs like the transport
coefficients, and the solutions are used to establish a connection with the
experimental observables like the nuclear suppression factor.

Transport coefficients are also important inputs to establish the stationary
distributions of the transport equations rafelskiwaltonprl . Ref.
rafelskiwaltonprl  shows that using the perturbative QCD estimates of the
drag, and the diffusion coefficients of the heavy quarks (HQ), the stationary
solution of the Fokker-Planck equation (that the authors initially assumed to
belong to a broad class characterized by some parameters) is a quasi-
exponential distribution. On the other hand, it was also shown that, without
assuming a mathematical form of the stationary distribution, it is possible to
generalize the Boltzmann transport equation lavagnopla ; wilkosada ; biroepja
so that the same quasi-exponential function is obtained as the stationary
solution in a more straightforward way. This generalized Boltzmann transport
equation has been a subject matter in many other works, and the quasi-
exponential stationary distribution it yields also arises from the nonadditive
statistics (also known as the nonextensive statistics) proposed by C. Tsallis
Tsal88 .

The nonadditive statistics arises in systems with a fluctuating ambiance
Wilk00 ; Wilk09 . Imprints of nonadditivity are found in the experimental data
obtained in high-energy collision experiments, for example, where the particle
spectra can be described in terms of the following quasi-exponential function
jpg20 ; Cleymansplb ,

$\displaystyle\left(1+\frac{q-1}{T}E\right)^{-\frac{1}{q-1}},$ (1)

that owes its origin to the Tsallis statistics and converges to the
exponential Boltzmann-Gibbs distribution in the limit $q\rightarrow 1$. In the
above equation $E=\sqrt{P^{2}+m^{2}}$ ($P\equiv|\vec{P}|$) is the single
particle energy of a particle of the mass $m$ moving with the three-momentum
$\vec{P}$, $q$ is the entropic parameter, and $T$ is the temperature
parameter. It can be shown that the $q$ parameter is connected to the relative
variance in (inverse) temperature Wilk00 . Also, scaling properties of the
Yang-Mills theory shows that $q$ can be deduced from the field theory
parameters like $N_{c}$ (no. of colours) and $N_{f}$ (no. of flavours)
deppmanprdq . The $q$ parameter is also utilized to generalize the molecular
chaos hypothesis (stosszahlansatz) used in the conventional Boltzmann
transport equation. This generalized stosszahlansatz, that may also be a
requirement for systems with small number of particles deppmanqideal , appears
in the generalized Boltzmann transport equation (gBTE) for a non-equilibrium
distribution $f$ given below (see, e.g., wilkosada ),

$\displaystyle
P^{\mu}\partial_{\mu}f^{q}=\mathcal{C}_{q}[f]~{}~{}(\mu=0,1,2,3).$ (2)

It is noteworthy that the power index $q$ on the left hand side ensures
energy-momentum conservation wilkosada , and the generalized molecular chaos
hypothesis is convoluted inside the collision term $\mathcal{C}_{q}$ that
dictates the evolution of the distribution $f$. It is also possible to
establish a connection between the quantity $\partial_{\mu}f^{q}$ and the
fractal derivatives. In the fractal approach, the Boltzmann equation,

$\frac{\partial f}{\partial t}+{\bf v}\cdot\frac{\partial f}{\partial{\bf
r}}+{\bf F}\cdot\frac{\partial f}{\partial{\bf v}}=C[{\bf p}(t+dt),{\bf
p}(t)],$ (3)

is modified to incorporate the fractal momentum space by imposing
deppmanfractalfpe

$C[{\bf p}(t+dt),{\bf p}(t)]=D_{t}^{\alpha_{q}}f\,,$ (4)

where $D_{t}^{\alpha_{q}}f$ is the fractal derivative defined by Parvate and
Gangal parvatefractal . An appropriate continuous approximation of the fractal
derivative leads to the $q$-deformed derivative deppmanplb ; deppmanfractal2 ,

$D_{t}^{\alpha_{q}}f=f^{q-1}\frac{df}{dt}\,,$ (5)

where $\alpha_{q}-1=1-q$. The relaxation time approximation of the BTE yields

$\frac{\partial f}{\partial t}+{\bf v}\cdot\frac{\partial f}{\partial{\bf
r}}+{\bf F}\cdot\frac{\partial f}{\partial{\bf v}}=-\frac{f-f_{0}}{\tau}\,,$
(6)

where $\tau$ is the relaxation time and $f_{0}$ is the stationary state
distribution. The approximation is valid for states near the stationary one.
The term in the r.h.s. of the equation above is a finite difference
approximation of the derivative $df/dt$ near the stationary state. By
considering the time-momentum space is fractal, this derivative must be
substituted by the fractal derivative operator or, in the appropriate
continuous approximation, the substitution

$\frac{df_{q}}{dt}\rightarrow
f_{q}^{q-1}\frac{df_{q}}{dt}\sim\frac{df_{q}^{q}}{dt}\,,$ (7)

where $f_{q}$ is a distribution characterized by the $q$ parameter.
Considering that the l.h.s. of Eq. (6) is the time derivative of the
distribution $f$, we get

$\frac{df}{dt}=\frac{df_{q}^{q}}{dt}\,.$ (8)

A general solution for this equation can be search with the ansatz

$f({\bf p},t)=\sum_{q^{\prime}}f_{q^{\prime}}({\bf p},t)\,.$ (9)

However, for systems associated with Tsallis Statistics, it is usually found
that the system is characterized by a single and well defined value $q$, so
$f\rightarrow f_{q}$. Due to Eq. (6), this condition is satisfied only if
$f=f_{q}^{q}$. With this result, Eq. (6) needs to be rewritten as

$\frac{\partial f_{q}^{q}}{\partial t}+{\bf v}\cdot\frac{\partial
f_{q}^{q}}{\partial{\bf r}}+{\bf F}\cdot\frac{\partial f_{q}^{q}}{\partial{\bf
v}}=-\frac{f_{q}-f_{0}}{\tau}\,,$ (10)

which is the equation utilized in Ref. tbphysica .

The Boltzmann transport equation can be subjected to different approximations
like the relaxation time approximation considered above. Recently an iterative
scheme has been developed to obtain approximate analytical solutions of the
generalized Boltzmann transport equation tbphysica . Although the relaxation
time approximation of the collision term is a useful first approximation for
studying transport phenomena, its assumptions may be valid only for
restrictive scenarios ashcroftbook . So, attempts have been made to go beyond
the relaxation time approximation maciejprd . One may also utilize the Landau
kinetic approximation landaubook  that observes that the most of the quark-
gluon scattering is soft. In this approximation, it is possible to obtain the
Fokker-Planck (FP) equation, and the Plastino-Plastino (PP) equation from the
BTE, and the gBTE respectively tbphysica ; svetitsky ; deppmanplb .

In this work, we investigate the energy loss of the heavy quarks in the quark-
gluon plasma medium dictated by the Plastino-Plastino dynamics involving
anomalous diffusion of the heavy-quarks. Following are the highlights of our
paper:

  1. 1.

We evaluate the transport coefficients from the Plastino-Plastino equation
that arises from the generalized Boltzmann transport equation, and has a
stationary distribution also obtained from the nonadditive statistics.

  2. 2.

We generalize the Rafelski-Walton formula rafelskiwaltonprl  connecting
differential energy loss (-$dE/dr$), drag, and diffusion coefficients, and
estimate energy loss per unit length by the heavy quarks traversing through a
nonadditive QGP medium.

  3. 3.

Later, we compute the jet quenching transport coefficient that is connected to
the transverse diffusion coefficient, and provide an estimate of the nuclear
suppression factor using the Dokshitzer-Kharzeev (DK) formula DK  and compare
the results with experimental data. Using the DK formula, we also provide
analytical results connecting nuclear suppression factor ($R_{\text{AA}}$),
jet quenching parameter ($\hat{q}$), and differential energy loss (-$dE/dr$).

A comparative study of the Fokker-Planck, and the Plastino-Plastino dynamics
of the heavy quarks have been performed in a recent article megiasplb .
However, in this work transport coefficients are estimated following Ref.
svetitsky  for small incoming momentum of the heavy quarks, and the present
work generalizes those estimates. Some estimates of the nuclear suppression
factor with the nonextensive initial distributions have been the subject
matter in some recent works tbphysica ; tsallisraaepja  that consider the
relaxation time approximation. The present work generalizes that framework and
the estimates by considering the microscopic scenario of the parton-parton
interaction.

The paper is organized as follows: in the next section we provide the
formalism to estimate differential energy loss of the heavy quarks from the
Plastino-Plastino (PP) transport coefficients. In Section III, we compute the
jet quenching transport coefficient ($\hat{q}$), and provide an estimate of
the nuclear suppression factor in the Plastino-Plastino framework. Section IV
contains analytical results establishing a relationship among nuclear
suppression factor ($R_{\text{AA}}$), jet quenching parameter ($\hat{q}$), and
differential energy loss. Section V contains a summary of our work,
conclusions, and discussion.

##  II Differential energy loss, and the Plastino-Plastino transport
coefficients

###  II.1 PP drag and diffusion transport coefficients

In this section, we shall use an extension of the Rafelski-Walton relationship
connecting the transport coefficients, and energy loss. However, before that,
we review how the transport coefficients like drag, and diffusion can be
evaluated from the Plastino-Plastino equation. First, we expand the l.h.s. of
the gBTE given by Eq. (2) assuming absence of an external force, and obtain
the following,

$\displaystyle\frac{\partial f^{q}(\vec{P},\vec{r},t)}{\partial
t}+\vec{v}.\vec{\nabla}f^{q}(\vec{P},\vec{r},t)=E^{-1}\mathcal{C}_{q}[f(\vec{P},\vec{r},t)],$
(11)

where $f$ is a non-equilibrium distribution (of the heavy-quarks) dependent on
three-momentum $\vec{P}\equiv\\{\vec{P_{\text{T}}},P_{\text{L}}\\}$ (T:
transverse; L: longitudinal), position $\vec{r}$, and time $t$. For a
relativistic system, energy of the heavy-quark of the mass $m$ is given by,
$E=\sqrt{P^{2}+m^{2}}$, and its velocity is given by $\vec{v}=E^{-1}\vec{P}$.

Now, we assume that the scale of spatial variations along the longitudinal
direction is sufficiently large, and the initial motion is one-dimensional
along the $x$-direction. By using the Bjorken’s assumption of longitudinal
boost-invariance (at a time $t$, conditions at point $x$ is the same as those
at point $x=0$) at the central rapidity region baym ; japrd , $f$ obeys

$\displaystyle f(\vec{P}_{t},P_{x},x,t)=f(\vec{P}_{t},P_{x}^{\prime},t),$ (12)

where $P_{x}^{\prime}=(P_{x}-u_{x}E)\gamma,~{}u_{x}=x/t,$ and
$~{}\gamma=(1-u_{x}^{2})^{-1/2}$. Noting that at $x=0$, $\partial
P_{x}^{\prime}/\partial x=-E/t$, we find,

$\displaystyle v_{x}\frac{\partial f^{q}}{\partial
x}=-\frac{P_{x}}{t}\frac{\partial f^{q}}{\partial P_{x}}.$ (13)

Hence, it is possible to write the gBTE for the heavy-quark distributions in
the following form tbphysica ,

$\displaystyle\left(\frac{\partial}{\partial
t}-\frac{P_{x}}{t}\frac{\partial}{\partial
P_{x}}\right)f^{q}(\vec{P}_{t},P_{x},t)=E^{-1}\mathcal{C}_{q}[f(\vec{P}_{t},P_{x},t)].$
(14)

For the binary collisions of particles with four-momenta
$P^{\mu}\equiv(E,\vec{P})$ (having the momentum distribution $f$), and
$Q^{\mu}\equiv(\mathcal{E},\vec{Q})$ (having the momentum distribution $g$),
that change to $P^{\prime\mu}\equiv(E^{\prime},\vec{P^{\prime}})$, and
$Q^{\prime\mu}\equiv(\mathcal{E^{\prime}},\vec{Q^{\prime}})$ after scattering,
the collision term of the gBTE can be written as wilkosada ; biroepja ,

$\displaystyle\mathcal{C}_{q}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int\frac{d^{3}\vec{Q}}{(2\pi)^{3}2\mathcal{E}}\frac{d^{3}\vec{Q^{\prime}}}{(2\pi)^{3}2\mathcal{E^{\prime}}}\frac{d^{3}\vec{P^{\prime}}}{(2\pi)^{3}2E^{\prime}}|\overline{\mathcal{M}}|^{2}(2\pi)^{4}\delta^{4}(P+Q-P^{{}^{\prime}}-Q^{{}^{\prime}})\left[h_{q}\\{f(\vec{P^{\prime}}),g(\vec{Q^{\prime}})\\}-h_{q}\\{f(\vec{P}),g(\vec{Q})\\}\right].$

$|\overline{\mathcal{M}}|^{2}$ is the square-averaged amplitude of the heavy
quark collisional processes (with the light quarks, and gluons) in a quark-
gluon plasma medium given by Ref. combridge . The function $h_{q}$ provides an
ansatz for the generalization of the molecular chaos hypothesis given by,

$\displaystyle h_{q}\\{f,g\\}=\exp_{q}\left[\log_{q}(f)+\log_{q}(g)\right],$
(16)

where

$\displaystyle\log_{q}f=\frac{1-f^{1-q}}{q-1}.$ (17)

Here $\log_{q}$ is the $q$-logarithm function and becomes the conventional
logarithm in the limit $q\rightarrow 1$, which also implies,

$\displaystyle\lim_{q\rightarrow 1}h_{q}\\{f,g\\}$ $\displaystyle=$
$\displaystyle\lim_{q\rightarrow
1}\exp_{q}\left[\log_{q}(f)+\log_{q}(g)\right]$ (18) $\displaystyle=$
$\displaystyle\exp\left[\log(f)+\log(g)\right]$ $\displaystyle=$
$\displaystyle f\times g,$

i.e., the distribution functions are independent, which is a consequence of
the molecular chaos hypothesis. In this limit, one gets back the collision
term of the conventional BTE svetitsky ,

$\displaystyle\mathcal{C}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int\frac{d^{3}\vec{Q}}{(2\pi)^{3}2\mathcal{E}}\frac{d^{3}\vec{Q^{\prime}}}{(2\pi)^{3}2\mathcal{E^{\prime}}}\frac{d^{3}\vec{P^{\prime}}}{(2\pi)^{3}2E^{\prime}}|\overline{\mathcal{M}}|^{2}(2\pi)^{4}\delta^{4}(P+Q-P^{{}^{\prime}}-Q^{{}^{\prime}})\left[f(\vec{P^{\prime}})g(\vec{Q^{\prime}})-f(\vec{P})g(\vec{Q})\right].$

Using the Landau kinetic approximation, it is possible to obtain tbphysica
the Plastino-Plastino equation from Eq. (14):

$\displaystyle\left.\frac{\partial f}{\partial
t}\right|_{P_{x}t}=-\frac{\partial}{\partial
P_{i}}\left(A_{i}^{\text{PP}}f\right)+\frac{\partial}{\partial P_{i}\partial
P_{j}}\left(B_{ij}^{\text{PP}}f^{2-q}\right)~{}~{}~{}~{}~{}~{}~{}~{}(i,j=1,2,3).$
(20)

$B_{ij}^{\text{PP}}$, and $A_{i}^{\text{PP}}$, the Plastino-Plastino diffusion
and drag coefficients, are given by

$\displaystyle B_{ij}^{\text{PP}}$ $\displaystyle=$
$\displaystyle\frac{f^{q-1}}{2}\left[\frac{1}{2Eq}\int\frac{d^{3}\vec{Q}}{(2\pi)^{3}2\mathcal{E}}\frac{d^{3}\vec{Q^{\prime}}}{(2\pi)^{3}2\mathcal{E^{\prime}}}\frac{d^{3}\vec{P^{\prime}}}{(2\pi)^{3}2E^{\prime}}~{}|\overline{\mathcal{M}}|^{2}(2\pi)^{4}\delta^{4}(P+Q-P^{{}^{\prime}}-Q^{{}^{\prime}})f^{1-q}S_{\vec{P},\vec{Q}}K_{i}K_{j}\right],$
(21a) $\displaystyle\equiv\frac{f^{q-1}}{2}\langle\langle
K_{i}K_{j}\rangle\rangle^{PP}$ $\displaystyle A_{i}^{\text{PP}}$
$\displaystyle=$
$\displaystyle\frac{1}{2Eq}\int\frac{d^{3}\vec{Q}}{(2\pi)^{3}2\mathcal{E}}\frac{d^{3}\vec{Q^{\prime}}}{(2\pi)^{3}2\mathcal{E^{\prime}}}\frac{d^{3}\vec{P^{\prime}}}{(2\pi)^{3}2E^{\prime}}~{}|\overline{\mathcal{M}}|^{2}(2\pi)^{4}\delta^{4}(P+Q-P^{{}^{\prime}}-Q^{{}^{\prime}})f^{1-q}S_{\vec{P},\vec{Q}}K_{i}\equiv\left<\left<K_{i}\right>\right>^{PP},$
(21b)

where the 3-momentum transfer is $K_{i}\equiv(P-P^{\prime})_{i}$ and the
factor $S_{\vec{P},\vec{Q}}$, given by

$\displaystyle S_{\vec{P},\vec{Q}}=\left[1+\frac{\left(1+\delta
q_{b}\frac{|\vec{Q}|}{T_{b}}\right)^{\frac{\delta q}{\delta
q_{b}}}}{\left(1+\delta q_{h}\frac{E}{T}\right)^{\frac{\delta q}{\delta
q_{h}}}}-\frac{1}{\left(1+\delta q_{h}\frac{E}{T}\right)^{\frac{\delta
q}{\delta q_{h}}}}\right]^{-\frac{1}{\delta q}},$ (22)

represents the interplay (characterized by $\delta q\equiv q-1$) between the
incoming heavy quark distribution (temperature $T$, entropic parameter $\delta
q_{h}\equiv q_{h}-1$), and the medium distribution (temperature $T_{b}$,
entropic parameter $\delta q_{b}\equiv q_{b}-1$). Both the distributions are
parameterized by a function similar to Eq. (1) and for all the calculations in
the paper we use $q_{h}=q_{b}=1.10,~{}q=1.002,~{}T_{h}=0.4$ GeV, and
$T_{b}=0.2$ GeV.

The factor $S_{\vec{P},\vec{Q}}$ arises due to generalized molecular chaos in
the collision term, and converges to the exponential Maxwell-Boltzmann medium
distribution (hence giving rise to molecular chaos) when $q,~{}q_{h}$, and
$q_{b}$ approach 1.

###  II.2 Differential energy loss

Utilizing the notation used in Eq. (21a), differential energy loss of the
heavy quarks in an environment dictated by the non-additive statistics can be
written as below following Ref. thomaprd :

$\displaystyle-\frac{dE}{dr}$ $\displaystyle=$
$\displaystyle\left\langle\left\langle\frac{E-E^{\prime}}{v}\right\rangle\right\rangle^{PP}=\left\langle\left\langle\frac{E^{2}-EE^{\prime}}{P}\right\rangle\right\rangle^{PP}.$
(23)

By defining $2B_{00}^{\text{PP}}\equiv
f^{q-1}\langle\langle(E-E^{\prime})^{2}\rangle\rangle^{PP}$, and from Eqs.
(21a), and (21b) it is possible to establish a relationship between $-dE/dr$,
and the heavy quark transport coefficients as shown below:

$\displaystyle-\frac{dE}{dr}$ $\displaystyle=$
$\displaystyle\left\langle\left\langle\frac{E^{2}-EE^{\prime}}{P}\right\rangle\right\rangle^{PP}$
(24) $\displaystyle=$
$\displaystyle\left\langle\left\langle\frac{P^{2}+m^{2}-EE^{\prime}-\vec{P}\cdot\vec{P^{\prime}}+\vec{P}\cdot\vec{P^{\prime}}}{P}\right\rangle\right\rangle^{PP}$
$\displaystyle=$
$\displaystyle\left\langle\left\langle\frac{1}{P}\left[P^{2}-\vec{P}.\vec{P^{\prime}}+\frac{1}{2}\left\\{E^{2}-P^{2}+E^{\prime
2}-P^{\prime
2}-2EE^{\prime}+2\vec{P}\cdot\vec{P^{\prime}}\right\\}\right]\right\rangle\right\rangle^{PP}$
$\displaystyle=$
$\displaystyle\left\langle\left\langle\left[\frac{P_{i}K_{i}}{P}+\frac{1}{2P}\left\\{(E-E^{\prime})^{2}-(\vec{P}-\vec{P^{\prime}})^{2}\right\\}\right]\right\rangle\right\rangle^{PP}$
$\displaystyle=$ $\displaystyle\frac{P_{i}}{P}\left\langle\left\langle
K_{i}\right\rangle\right\rangle^{PP}+\frac{1}{2P}\left\langle\left\langle(E-E^{\prime})^{2}\right\rangle\right\rangle^{PP}-\frac{1}{2P}\left\langle\left\langle(\vec{P}-\vec{P^{\prime}})^{2}\right\rangle\right\rangle^{PP}$
$\displaystyle=$
$\displaystyle\frac{P_{i}A_{i}^{\text{PP}}+f^{1-q}B_{00}^{\text{PP}}-f^{1-q}B_{ii}^{\text{PP}}}{P}.$

In the above equation (also in Eqs. 28, and 42), repeated indices are summed
over. It is noteworthy that in the limit $q\rightarrow 1$, Eq. (24) converges
to a similar equation obtained in Ref. rafelskiwaltonprl  for the Boltzmann-
Gibbs statistics. In the high-energy regime ($E,E^{\prime}>>m$) where small
angle scattering is more dominant, the above equation can be approximated as,

$\displaystyle-\frac{dE}{dr}$ $\displaystyle\approx$
$\displaystyle\gamma^{\text{PP}}P,$ (25)

where $A_{i}^{\text{PP}}=\gamma^{\text{PP}}P_{i}$, and $P_{i}P_{i}=P^{2}$. We
evaluate the transport coefficients following the method elaborated in Ref.
svetitsky . Fig. 1 compares the differential energy loss in the Plastino-
Plastino, and the Fokker-Planck formalism. It shows that the Plastino-Plastino
energy loss starts surpassing the Fokker-Planck counterpart at around 2 GeV
momentum of the incoming charm quark.

Figure 1: Comparison of the energy loss of the charm quarks ($m$=1.5 GeV)
following the FP and PP dynamics in the QGP of $T_{b}$= 0.2 GeV.

Experimental evidences of energy loss in the QGP is manifest through measuring
the nuclear suppression factor defined by,

$R_{\mathrm{AA}}=\frac{\left(d^{2}N/dp_{\mathrm{T}}dy\right)^{\mathrm{A+A}}}{N_{\mathrm{coll}}\times\left(d^{2}N/dp_{\mathrm{T}}dy\right)^{\mathrm{p+p}}}.$
(26)

$N_{\mathrm{coll}}$ is the number of nucleon-nucleon binary collisions (p+p)
while a nucleus ‘A’ collides with another nucleus. $d^{2}N/dp_{\mathrm{T}}dy$
denotes the differential yield within the transverse momentum and rapidity
range $p_{\rm{T}}$ to $p_{\mathrm{T}}+dp_{\mathrm{T}}$ and $y$ to $y+dy$.
Theoretical estimates of $R_{\mathrm{AA}}$ are obtained by solving the
transport equation. However, in this work, we aim to estimate $R_{\text{AA}}$
by computing another transport coefficient $\hat{q}$, the jet-quenching
parameter, that yields the transverse momentum broadening per unit length.
Dokshitzer and Kharzeev DK  have derived an analytical result connecting
$\hat{q}$, and $R_{\mathrm{AA}}$, and the rest of the article will be devoted
to computing $\hat{q}$, and providing numerical estimates of
$R_{\mathrm{AA}}$. We also will discuss analytical results connecting the
transport coefficients and nuclear suppression factor.

##  III Calculating the jet-quenching parameter $\hat{q}$

In the previous section, we have derived the Plastino-Plastino equation from
the gBTE, and found the integral expressions for the drag and diffusion
coefficients, $A^{\text{PP}}_{i}$, and $B^{\text{PP}}_{ij}$. The diffusion
tensor $B^{\text{PP}}_{ij}$ can be written in terms of a longitudinal and a
transverse diffusion coefficient in the following way svetitsky ,

$\displaystyle B^{\text{PP}}_{ij}$ $\displaystyle=$
$\displaystyle\left(\delta_{ij}-\frac{{P_{i}}{P_{j}}}{{{P}^{2}}}\right)B^{\text{PP}}_{\text{T}}({P^{2}})+\frac{{P_{i}}{P_{j}}}{{P^{2}}}B^{\text{PP}}_{\text{L}}(P^{2})~{}~{}~{}(P^{2}\equiv|\vec{P}|^{2}).$
(27)

Hence, the trace of the diffusion tensor (that appears in Eq. 24) is
$B_{ii}^{\text{PP}}=2B^{\text{PP}}_{\text{T}}+B^{\text{PP}}_{\text{L}}$. From
Eq. (27) we obtain,

$\displaystyle B^{\text{PP}}_{\text{T}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\delta_{ij}-\frac{{P_{i}}{P_{j}}}{{P^{2}}}\right)B^{\text{PP}}_{ij}$
(28) $\displaystyle=$
$\displaystyle\frac{f^{q-1}}{4}\left(\delta_{ij}-\frac{{P_{i}}{P_{j}}}{{P^{2}}}\right)\langle\langle
K_{i}K_{j}\rangle\rangle^{PP}$ $\displaystyle=$
$\displaystyle\frac{f^{q-1}}{4}\left\langle\left\langle
K^{2}-\frac{(\vec{P}\cdot\vec{K})^{2}}{P^{2}}\right\rangle\right\rangle^{PP}$
$\displaystyle=$ $\displaystyle\frac{f^{q-1}}{4}\left\langle\left\langle
K_{\text{T}}^{2}\right\rangle\right\rangle^{PP},$

where
$\vec{K}\equiv(\vec{K}_{\text{T}},K_{\text{L}})\equiv(\vec{K}_{\text{T}},\vec{K}\cdot\vec{P}/P)$.
Further, $\langle\langle K_{\text{T}}^{2}\rangle\rangle$, that gives
transverse momentum broadening per unit time, can be connected to the jet-
quenching coefficients $\hat{q}$ that estimates average transverse momentum
broadening per unit length in the following way SolovevaPRD :

$\displaystyle\hat{q}=\frac{\langle\langle
K_{\text{T}}^{2}\rangle\rangle}{v_{\text{L}}}=\frac{E}{P_{\text{L}}}\langle\langle
K_{\text{T}}^{2}\rangle\rangle~{}~{}(v_{\text{L}}:~{}\text{longitudinal~{}velocity}).$
(29)

Hence, from Eq. (28)

$\displaystyle\hat{q}^{\text{PP}}=4f^{1-q}\frac{E}{P_{\text{L}}}B^{\text{PP}}_{\text{T}}.$
(30)

It is noteworthy that, in the limit $q\rightarrow 1$, a similar relationship
between $\hat{q}$, and $B_{\text{T}}$ has also been derived in Ref. TBPrd  for
the kinematic region $P_{\text{L}}>>m,P_{\text{T}}$ that implies
$\hat{q}\approx 4B_{\text{T}}$.

Now, the jet-quenching transport coefficient $\hat{q}$ is connected to the
nuclear suppression factor $R_{\text{AA}}$ that is an experimental observable.
So, with our estimates of $\hat{q}$, we are also able to provide an estimate
of $R_{\text{AA}}$ in a nonadditive medium. In particular, it will be
interesting to see what effect the entropic $q$ parameter (that parameterizes
the generalized molecular chaos) has on heavy-quark suppression. An analytical
formula relating $\hat{q}$ and the nuclear suppression factor
$R^{\text{h}}_{\text{AA}}$ for the heavy quarks (denoted by ‘h’) has been used
in Ref. laceyPRL , and derived in Ref. DK . This formula is given by (for a
path length $\ell$, strong coupling $\alpha_{s}$, and colour factor $C_{F}$),

$\displaystyle
R^{\text{h}}_{\text{AA}}(P_{\text{T}},\ell)\simeq\exp\left[-\frac{2\alpha_{s}C_{F}}{\sqrt{\pi}}\ell\sqrt{\frac{\hat{q}\mathcal{L}_{h}^{\text{abs}}}{P_{\text{T}}}}+\frac{16\alpha_{s}C_{F}}{9\sqrt{3}}\ell\left(\frac{\hat{q}m^{2}}{m^{2}+P_{\text{T}}^{2}}\right)^{1/3}\right].$
(31)

However, this formula for $R^{\text{h}}_{\text{AA}}$ is derived considering
radiative processes. Hence, we need to calculate the radiative transport
coefficients to be discussed in the following section, along with the
evaluation of the quantity $\mathcal{L}_{h}^{\text{abs}}$.

###  III.1 Estimating the jet quenching parameter $\hat{q}$: collisional and
radiative

Figure 2: Comparison of total (collisional+radiative), and collisional jet-
quenching parameters ($\hat{q}$) of the charm quarks ($m$=1.5 GeV;
$P_{\text{L}}$=10 GeV) following the FP and PP ($q_{b}$=1.10) dynamics in the
QGP of $T_{b}$= 0.2 GeV.

According to Eq. (30), $\hat{q}=4(E/P_{\text{L}})f^{1-q}B_{\text{T}}$, and is
given by the following integration in the PP case for the collisional (e.g.
HQ+LQ/g $\rightarrow$ HQ+LQ/g) processes.

$\displaystyle\hat{q}^{\text{PP}}_{\text{coll}}=\frac{1}{2qP_{\text{L}}}\int\frac{d^{3}\vec{Q}}{(2\pi)^{3}2\mathcal{E}}\frac{d^{3}\vec{Q^{\prime}}}{(2\pi)^{3}2\mathcal{E^{\prime}}}\frac{d^{3}\vec{P^{\prime}}}{(2\pi)^{3}2E^{\prime}}|\overline{\mathcal{M}}|_{2\rightarrow
2}^{2}(2\pi)^{4}\delta^{4}(P+Q-P^{\prime}-Q^{\prime})S_{\vec{P},\vec{Q}}\left(P^{\prime
2}-\frac{(\vec{P}.\vec{P^{\prime}})^{2}}{P^{2}}\right),$ (32)

For the radiative processes (e.g. HQ+LQ/g $\rightarrow$ HQ+LQ/g+g), we follow
Refs. TBPrd ; SKDPrd  to calculate the radiative transport coefficients from
the collisional ones. To be consistent with Eq. (31), we consider only
collinear, and soft gluon (4-momentum
$k_{5}^{\mu}\equiv\\{E_{5},\vec{k}_{5}\\}$) radiation without the Bose
enhancement factor in the present calculation. These considerations lead to
the following estimate of the radiative jet transport coefficient,

$\displaystyle\hat{q}^{\text{PP}}_{\text{rad}}$ $\displaystyle=$
$\displaystyle\frac{1}{2qP_{\text{L}}}\int\frac{d^{3}\vec{Q}}{(2\pi)^{3}2\mathcal{E}}\frac{d^{3}\vec{Q^{\prime}}}{(2\pi)^{3}2\mathcal{E^{\prime}}}\frac{d^{3}\vec{P^{\prime}}}{(2\pi)^{3}2E^{\prime}}\frac{d^{3}\vec{k_{5}}}{(2\pi)^{3}2E_{5}}(2\pi)^{4}\delta^{4}(P+Q-P^{\prime}-Q^{\prime}-k_{5})S_{\vec{P},\vec{Q}}\left(P^{\prime
2}-\frac{(\vec{P}.\vec{P^{\prime}})^{2}}{P^{2}}\right)$ (33)
$\displaystyle\times~{}|\overline{\mathcal{M}}|_{2\rightarrow
3}^{2}\theta_{1}\left(E-E_{5}\right)\theta_{2}\left(\tau-\tau_{\text{F}}\right),$

where the first $\theta$-function prohibits gluon radiation of energy greater
than that of the incoming heavy quark, and the second $\theta$-function
accounts for the additive region of the gluon spectra due to multiple
scattering in which the interaction time ($\tau$) in greater than the
formation time ($\tau_{\text{F}}$) of the gluons PlumerPRD . These two
$\theta$-functions put an upper-bound and a lower-bound respectively in the
gluon 3-momentum ($\vec{k}_{5}\equiv\\{\vec{k}_{\text{T}},k_{\text{L}}\\}$)
integration. The radiative matrix element for the soft, collinear gluon
emission at an angle $\theta_{g}$ can be written in the following way DK  (for
$g_{s}^{2}=4\pi\alpha_{s}$; $\alpha_{s}\equiv$ strong coupling),

$\displaystyle|\overline{\mathcal{M}}|_{2\rightarrow
3}^{2}\approx|\overline{\mathcal{M}}|_{2\rightarrow
2}^{2}\times\frac{12g_{s}^{2}}{k_{\text{T}}^{2}}\left(1+\frac{m^{2}}{E^{2}\theta_{g}^{2}}\right)^{-2}.$
(34)

Also, parameterizing the soft gluon 4-momentum in terms of its rapidity
($\eta$), we get,

$\displaystyle\hat{q}^{\text{PP}}_{\text{rad}}$ $\displaystyle\approx$
$\displaystyle\int\text{collisional~{}part}\times\frac{6\alpha_{s}}{\pi}\int
d\eta\int\frac{dk_{\text{T}}}{k_{\text{T}}}~{}\left(1+\frac{m^{2}}{E^{2}\theta_{g}^{2}}\right)^{-2}\theta_{1}\left(E-E_{5}\right)\theta_{2}\left(\tau-\tau_{\text{F}}\right).$
(35)

Now, using the parameterization $E_{5}=k_{\text{T}}\cosh\eta$, putting
$\tau_{\text{F}}=\cosh\eta/k_{\text{T}}$, and using
$\theta_{g}=2\tan^{-1}\left[\exp(-\eta)\right]$,

$\displaystyle\hat{q}^{\text{PP}}_{\text{rad}}$ $\displaystyle\approx$
$\displaystyle\int\text{collisional~{}part}\times\frac{6\alpha_{s}}{\pi}\int
d\eta\int_{\Lambda\cosh\eta}^{\frac{E}{\cosh\eta}}\frac{dk_{\text{T}}}{k_{\text{T}}}~{}\left[1+\frac{m^{2}}{E^{2}\left(2\tan^{-1}\left[\exp(-\eta)\right]\right)^{2}}\right]^{-2}.$
(36)

In the above integration, the interaction rate $\Lambda$ is evaluated from the
following integral,

$\displaystyle\Lambda=\frac{1}{2qE}\int\frac{d^{3}\vec{Q}}{(2\pi)^{3}\mathcal{E}}\frac{d^{3}\vec{Q^{\prime}}}{(2\pi)^{3}\mathcal{E^{\prime}}}\frac{d^{3}\vec{P^{\prime}}}{(2\pi)^{3}E^{\prime}}|\overline{\mathcal{M}}|_{2\rightarrow
2}^{2}(2\pi)^{4}\delta^{4}(P+Q-P^{\prime}-Q^{\prime})f^{1-q}S_{\vec{P},\vec{Q}}.$
(37)

In Fig. 2 we show the variation of the (Plastino-Plastino and Fokker-Planck)
jet transport coefficient
$\hat{q}^{\text{PP/FP}}=\hat{q}^{\text{PP/FP}}_{\text{coll}}+\hat{q}^{\text{PP/FP}}_{\text{rad}}$
with the incoming HQ transverse momentum $p_{\text{T}}$.

###  III.2 Estimating $\mathcal{L}_{h}^{\text{abs}}$, and $R_{\text{AA}}$ and
comparing with experimental data

Figure 3: Comparison of the nuclear suppression factor ($R_{\text{AA}}$) of
the charm quarks ($m$=1.5 GeV; $P_{\text{L}}$=10 GeV) following the FP and PP
dynamics in the QGP of $T_{b}$= 0.2 GeV calculated using Eq. (31). We put
$b=1.7,~{}n=11,~{}\ell=3$ fm.

For evaluating $\mathcal{L}_{h}^{\text{abs}}$, we shall follow Ref. DK  that
defines

$\displaystyle\mathcal{L}_{h}=\frac{d}{d\ln
p_{\text{T}}}\ln\left[\frac{d\sigma^{\text{vac}}}{dp_{\text{T}}^{2}}\left(p_{\text{T}}\right)\right],$
(38)

and uses the following vacuum differential cross-section for the D-mesons
dmeson

$\displaystyle\frac{d\sigma^{\text{vac}}}{dp_{\text{T}}^{2}}=C\left(\frac{1}{bm^{2}+p_{\text{T}}^{2}}\right)^{\frac{n}{2}},$
(39)

for $b=1.4\pm 0.3$, $n=10.0\pm 1.2$, and $m=1.5$ GeV. Hence, the quantity
$\mathcal{L}_{h}$ can be evaluated to be,

$\displaystyle\mathcal{L}_{h}=-\frac{np_{\text{T}}^{2}}{bm^{2}+p_{\text{T}}^{2}}\Rightarrow\mathcal{L}_{h}^{\text{abs}}=\left|\mathcal{L}_{h}\right|.$
(40)

We plot the variation of the charm quark nuclear suppression factor with its
transverse momentum in Fig. 3.

We also compare the theoretical model calculations with the experimental data
ALICED0  of the $D^{0}$ mesons in Fig. 4. To this end, we integrate Eq. (31)
with the following function peterson  fragmenting heavy-quarks (‘h’) to the
hadrons (‘H’) with the momentum fraction $z$ of the partons:

$\displaystyle
D^{\text{H}}_{\text{h}}(z)\propto\frac{1}{z\left(1-\frac{1}{z}-\frac{\epsilon_{\text{h}}}{(1-z)}\right)^{2}}.$
(41)

Figure 4: Comparison of nuclear suppression factor of the $D^{0}$ mesons
calculated using Eqs. (31), and (41) for the PP and FP dynamics with
experimental data obtained by the ALICE collaboration ALICED0 . We put
$T_{b}=0.2$ GeV, $b$=1.7, $n$=11, and $\epsilon_{\text{h}}$=0.056.

##  IV Relationship among jet quenching parameter, differential energy loss,
and nuclear suppression factor

From Eqs. (24), and (30), one can establish a relationship between $-dE/dr$,
and $\hat{q}$ in the following way:

$\displaystyle-\frac{dE}{dr}$ $\displaystyle=$
$\displaystyle\frac{P_{i}A_{i}^{\text{PP}}+f^{1-q}B_{00}^{\text{PP}}-f^{1-q}B_{ii}^{\text{PP}}}{P}$
(42) $\displaystyle=$
$\displaystyle\frac{P_{i}A_{i}^{\text{PP}}+f^{1-q}(B_{00}^{\text{PP}}-2B^{\text{PP}}_{\text{T}}-B^{\text{PP}}_{\text{L}})}{P}$
$\displaystyle=$ $\displaystyle
P\gamma^{\text{PP}}+f^{1-q}\left(\frac{B_{00}^{\text{PP}}-B^{\text{PP}}_{\text{L}}}{P}\right)-\frac{\hat{q}^{\text{PP}}P_{\text{L}}}{2PE}$

In the Boltzmann-Gibbs (BG) limit (in which case we simply remove the
superscript ‘PP’), the above equation becomes,

$\displaystyle-\frac{dE}{dr}=P\gamma+\frac{B_{00}-B_{\text{L}}}{P}-\frac{\hat{q}P_{\text{L}}}{2PE}.$
(43)

The transport coefficients $\gamma$, $B_{00}$, $\hat{q}$ and $B_{\text{L}}$
are the conventional Fokker-Planck transport coefficients. Eq. (42) also
implies,

$\displaystyle\hat{q}^{\text{PP}}$ $\displaystyle=$
$\displaystyle\frac{2PE}{P_{\text{L}}}\left[P\gamma^{\text{PP}}+f^{1-q}\left(\frac{B_{00}^{\text{PP}}-B^{\text{PP}}_{\text{L}}}{P}\right)+\frac{dE}{dr}\right]\text{(Plastino-
Plastino),}~{}~{}\text{and}$ $\displaystyle\hat{q}$ $\displaystyle=$
$\displaystyle\frac{2PE}{P_{\text{L}}}\left(P\gamma+\frac{B_{00}-B_{\text{L}}}{P}+\frac{dE}{dr}\right)~{}~{}~{}~{}\text{(Fokker-
Planck)}.$ (44)

From Eq. (44), it is possible to establish a connection between nuclear
suppression factor, and differential energy loss in the following way:

$\displaystyle
R^{\text{h}}_{\text{AA}}(P_{\text{T}},\ell)\simeq\exp\left[-\frac{2\alpha_{s}C_{F}}{\sqrt{\pi}}\ell\sqrt{\frac{\frac{2PE}{P_{\text{L}}}\left(P\gamma^{\text{PP}}+\frac{f^{1-q}(B_{00}^{\text{PP}}-B^{\text{PP}}_{\text{L}})}{P}+\frac{dE}{dr}\right)\mathcal{L}_{h}^{\text{abs}}}{P_{\text{T}}}}\right.$
$\displaystyle\left.+\frac{16\alpha_{s}C_{F}}{9\sqrt{3}}\ell\left(\frac{\frac{2PE}{P_{\text{L}}}\left(P\gamma^{\text{PP}}+\frac{f^{1-q}(B_{00}^{\text{PP}}-B^{\text{PP}}_{\text{L}})}{P}+\frac{dE}{dr}\right)m^{2}}{m^{2}+P_{\text{T}}^{2}}\right)^{1/3}\right],$
(45)

that, in the BG limit ($q\rightarrow 1$), becomes,

$\displaystyle
R^{\text{h}}_{\text{AA}}(P_{\text{T}},\ell)\simeq\exp\left[-\frac{2\alpha_{s}C_{F}}{\sqrt{\pi}}\ell\sqrt{\frac{\frac{2PE}{P_{\text{L}}}\left(P\gamma+\frac{B_{00}-B_{\text{L}}}{P}+\frac{dE}{dr}\right)\mathcal{L}_{h}^{\text{abs}}}{P_{\text{T}}}}\right.$
$\displaystyle\left.+\frac{16\alpha_{s}C_{F}}{9\sqrt{3}}\ell\left(\frac{\frac{2PE}{P_{\text{L}}}\left(P\gamma+\frac{B_{00}-B_{\text{L}}}{P}+\frac{dE}{dr}\right)m^{2}}{m^{2}+P_{\text{T}}^{2}}\right)^{1/3}\right].$
(46)

##  V Discussion and conclusion

In this work, we revisit the derivation of the Plastino-Plastino equation from
a generalized Boltzmann transport equation. The index $q$ also appears in Ref.
deppmanfractalfpe  (as a power of the test particle distribution) that derives
the Plastino-Plastino equation using the fractal space approach. We calculate
the transport coefficients like the diffusion and the jet-quenching parameter
of the charm quarks, and differential energy loss in a nonadditive QGP medium.
We have generalized some previously used relationships (Eqs. (24), and (30)),
and presented analytical formulae establishing a connection between
$R_{\text{AA}}$, an experimental observable, and the transport coefficients
using Ref. DK  (Eqs. (45), and (46)). From this paper it is also possible to
estimate the elliptic flow parameter by defining path-length-dependent
$R_{\text{AA}}$ for the in-plane and out-of-plane directions RAAv2 .

We observe that when there is an interplay between the incoming heavy-quarks
and the medium particle distribution, manifested by a generalized molecular
chaos, numerical values of transport coefficients, energy loss, and
suppression of the heavy-quarks increase compared to their Boltzmann-Gibbs
counterparts. Given that for finite systems generalized molecular chaos may be
a more realistic scenario deppmanqideal , our estimates in this paper will be
relevant for the QGP system. The $\hat{q}$ values obtained using the PP
dynamics is around 1 GeV2-fm-1 which is comparable to the value obtained in
Ref. laceyPRL  and less compared to the ‘experimentally favoured range’ for
the RHIC data armestoqhat . However, the DK formula considers soft, eikonal
and collinear gluon radiation approximation, and it will be interesting to
investigate how generalizing these approximation considered in Ref. DK ,
following Ref. tbnoneikonal , will impact the estimate of $\hat{q}$, and other
experimental observables. With these modifications, it can also be expected
that the contribution of radiation in $\hat{q}$ will increase, more so for the
PP case where $\hat{q}_{\text{rad}}$ is around 40-60% of
$\hat{q}_{\text{coll}}$. A few salient points that come out of the discussion
are listed below:

  * •

Calculation of $R_{\text{AA}}$ is based on calculating the jet transport
coefficient evaluated using the Plastino-Plastino equation that takes into
account anomalous diffusion of the heavy-quarks in the medium.

  * •

In many a case, the parameters in the analysis are fixed from the experimental
results to minimize arbitrariness.

  * •

For the given parameter values, there is a good agreement between the
theoretical calculation and $R_{\text{AA}}$ data, indicating that fast heavy-
quarks may be subjected to anomalous diffusion inside the QGP anodiff .

Our calculation of $R_{\text{AA}}$ using the PP dynamics captures the trend of
decreasing suppression (considering the central values) at higher transverse
momenta (Fig. 4), but the FP dynamics does not show this trend. It can be
intuitively understood that $R_{\text{AA}}$ depends on $\ell$ that may be
equalled to the average path length traversed by the heavy-quarks. This
quantity $\ell$ can introduce a centrality dependence in our calculations. We
observe that for a 0-20 % central collision between nuclei with the same
radius (e.g. Pb on Pb), the average path length of the heavy-quarks varies
between 55% to 75 % of the radius of the produced system. So, our analysis
with the PP dynamics indicates that the diameter of the QGP system produced
may be around 8-10 fm, but the FP dynamics indicates a larger diameter of
around 13-18 fm. Nevertheless, we think that a more detailed analysis
involving centrality dependent experimental data and interferometry may be
required. For the lower $p_{\text{T}}$ values, more involved investigation
(perhaps involving collectivity) may benefit the study. This also begs a more
rigorous approach to establish a connection between theory and experiment by
utilizng a numerical solution of the Plastino-Plastino equation with the
transport coefficients as inputs. In this extension, the collective behaviour
of medium (that is not considered in the present work) should be included for
a more reliable comparison with experimental data. We reserve these works for
the future.

## Acknowledgements

The work of EM is supported by the project PID2020-114767GB-I00 and by the
Ramón y Cajal Program under Grant RYC-2016-20678 funded by
MCIN/AEI/10.13039/501100011033 and by “FSE Investing in your future”, by Junta
de Andalucía under Grant FQM-225, and by the “Prórrogas de Contratos Ramón y
Cajal” Program of the University of Granada. This work is a part of the
project INCT-FNA Proc. No. 464898/2014-5. A.D. acknowledges the support from
the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq-
Brazil), grant 306093/2022-7.

## References

  * (1) J. Rafelski and D. B. Walton Phys. Rev. Lett. 84, 31 (1999).
  * (2) A. Lavagno, Phys. Lett. A 301, 13 (2002).
  * (3) T. Osada G. Wilk, Phys. Rev. C 77, 044903 (2009).
  * (4) T. S. Biro and E. Molnar, Eur. Phys. J. A 48, 172 (2012)
  * (5) C. Tsallis, J. Stat. Phys. 52, 479 (1988).
  * (6) G. Wilk, and Z. Włodarczyk, Phys. Rev. Lett. 84, 2770 (2000).
  * (7) G. Wilk, and Z. Włodarczyk, Phys. Rev. C 79, 054903 (2009).
  * (8) M. D. Azmi et al.,J. Phys. G 47, 045001 (2020).
  * (9) J. Cleymans et al., Phys. Lett. B 723, 351-354 (2013).
  * (10) A. Deppman, E. Megías and D. P. Menezes, Phys. Rev. D 101, 034019 (2019).
  * (11) J. S. Lima and A. Deppman, Phys. Rev. E 101, 040102 (R) (2020).
  * (12) E. Megías, A. K. Golmankhaneh and A. Deppman, Phys. Lett. B 848, 2024.
  * (13) A. Parvate and A. Gangal, Fractals 17 (01) 53 (2009).
  * (14) A. Deppman, E. Megías, R. Pasechnik, Entropy 25 (7), 1008 (2023).
  * (15) A. Deppman et al., Phys. Lett. B 839, 137752 (2023).
  * (16) T. Bhattacharyya, Physica A 624, 128910 (2023).
  * (17) N. W. Ashcroft, and N. D. Mermin, _Solid State Physics_ (USA: Harcourt College Publishers, 1975).
  * (18) M. Rybcziński, G. Wilk and Z. Włodarczyk, Phys. Rev. D 103, 114026 (2021).
  * (19) E. M. Lifshitz, and L. P. Pitaevskii, _Physical Kinetics_ (Pergamon, New York, 1981), Chap. 2.
  * (20) B. Svetitsky, Phys. Rev. D 37, 2484 (1988).
  * (21) Y. L. Dokshitzer, and D. E. Kharzeev, Phys. Lett. B 519, 199 (2001).
  * (22) E. Megías et al., Phys. Lett. B 845 138136 (2023).
  * (23) S. Tripathy et al., Eur. Phys. J. A 52, no. 9, 289 (2016).
  * (24) G. Baym, Phys. Lett. B 138, 18 (1984).
  * (25) P. Roy, A. K. Dutt-Mazumder, and J.-e Alam, Phys. Rev. D 73, 044911 (2006).
  * (26) B. L. Combridge, Nucl. Phys. B 151, 429 (1979).
  * (27) E. Braaten, and M. H. Thoma, Phys. Rev. D 44, R2625 (1991); 44, 1298 (1991).
  * (28) T. Song, I. Grishmanovskii, and O. Soloveva, Phys. Rev. D 107, 036009 (2023).
  * (29) S. Mazumder, T. Bhattacharyya, and J.-e Alam, Phys. Rev. D 89, 014002 (2014).
  * (30) R. A. Lacey et al., Phys. Rev. Lett. 103, 142302 (2009).
  * (31) A. Shaikh et al., Phys. Rev. D 104, 034017 (2021).
  * (32) X. N. Wang, M. Gyulassy, and M. Plumer, Phys. Rev. D 51, 3436 (1995).
  * (33) C. A. Alves et al., E769 Collaboration, Phys. Rev. Lett. 77 2392 (1996).
  * (34) B. Abelev et al., The ALICE Collaboration, Jour. of High En. Phys. 09, 112 (2012).
  * (35) C. Peterson et al., Phys. Rev. D 27, 105 (1983).
  * (36) D. Zigic et al. J. Phys. G: Nucl. Part. Phys. 46, 085101 (2019).
  * (37) N. Armesto et al., Phys. Lett. B 637 362 (2006).
  * (38) T. Bhattacharyya, S. Mazumder, and R. Abir, Adv. High Energy Phys. 2016, 1298986 (2016).
  * (39) Paul Caucal and Yacine Mehtar-Tani, Phys. Rev. D 106, L051501 (2022).