# Application of Molecular Topology to the Prediction of Antioxidant Activity
in a Group of Phenolic Compounds

Jaime Barros Silva Filho
Fernando de Souza Bastos
Diogo da Silva Machado
Maria Luiza Ferreira Delfim

###### Abstract

The study of compounds with antioxidant capabilities is of great interest to
the scientific community, as it has implications in several areas, from
Agricultural Sciences to Biological Sciences, including Food Engineering,
Medicine and Pharmacy. In applications related to human health, it is known
that antioxidant activity can delay or inhibit oxidative damage to cells,
reducing damage caused by free radicals, helping in the treatment, or even
preventing or postponing the onset of various diseases. Among the compounds
that have antioxidant properties, there are several classes of Phenolic
Compounds, which include several compounds with different chemical structures.
In this work, based on the molecular branching of compounds and their
intramolecular charge distributions, and using Molecular Topology, we propose
a significant topological-mathematical model to evaluate the potential of
candidate compounds to have an antioxidant function.

##  1 Introduction

Studies and research on compounds with antioxidant activities have gained
prominence in recent years. This is due to the antioxidant activity delay or
inhibiting the oxidative damage of cells, decreasing the damage caused by free
radicals, preventing or postponing the onset of various diseases and also
helping in the treatment of diseases such as diabetes, inflammatory diseases,
cardiovascular disorders, associated with oxidative stress, some types of
cancer and aging (see eg [24, 19, 27, 28, 13, 20]).

Measuring the antioxidant activity/capacity of foods and biological compounds
is therefore essential not only to ensure food quality, but also to assess the
efficiency of dietary antioxidants in the prevention and treatment of diseases
related to oxidative stress.

Among the compounds that have antioxidant properties, there are several
classes of Phenolic Compounds, which are those that have one or more hydroxyl
groups directly linked to an aromatic ring. Such compounds represent a large
group of molecules with a variety of functions related to plant growth,
development and defense. They include signaling molecules, pigments and
flavors that can attract or repel, as well as compounds that can protect the
plant against insects, fungi, bacteria and viruses.

With the growing recognition of their bioactivity values, the research effort
on the antioxidant behavior of Phenolic Compounds, which include thousands of
compounds with different chemical structures, has increased significantly in
recent decades. Furthermore, the methods and instruments used to measure
antioxidant activity/capacity have made remarkable progress. According to
[21], there are several methods for assessing this ability, which generally
fall into three distinct categories, namely, spectrometry, electrochemical
assays and chromatography. However, all are based on chemical reactions that
demand time and financial resources.

Therefore, the implementation of mechanisms that can provide an early estimate
of the potential of a molecule to exert an antioxidant function is an
interesting research strategy. Molecular Topology based on Quantitative
Structure Activity Relations (QSAR) is a tool that has great potential in this
context. The QSAR technique, developed by [11], combines Topology and
Statistics concepts and its approach can be described as a statistical method
of data analysis to develop models that can correctly predict certain
biological activity or properties of compounds based on their structure
chemistry. QSAR techniques apply descriptors based on molecular structures and
use algorithms to correlate the obtained descriptors with the value of the
target property of interest, for more information we suggest [4].

The development of Molecular Topology is attributed to Randic [23] and Kier-
Hall [10, 15] and is based, above all, on the premise that, in many cases,
there is a close relationship between the structures of organic compounds and
many of its chemical and biological properties. From which, a set of suitable
numerical characterizations for the molecules of interest is obtained. Such
characterizations are designated by topological indices, which allow, after
statistical treatment (specifically, Linear Discriminant Analysis), the
classification into groups, with greater and lesser probabilities of having a
chemical function of interest. It should be noted that Molecular Topology has
been used satisfactorily in obtaining new compounds for the production of
drugs, cosmetics and agrochemicals of great interest to society, see more
details in [8, 18] and [1].

To obtain the topological indices, the molecules to be considered are
mathematically modeled as Topological Graphs, naturally taking their atoms as
vertices and atomic bonds as edges. Thus, from the graphs, numerical matrices
are obtained containing information about the molecules (such as the
structural arrangement, global charge transfer and between pairs of atoms,
atomic electronegativity, among others), which, through specific formulas,
lead to the survey of topological indices.

The main objective of this study was to develop a predictive QSAR model to
investigate the antioxidant activity of a group of Phenolic Compounds using
Molecular Topology and Linear Discriminant Analysis.

##  2 Materials and methods

###  2.1 Compounds Analyzed and Their Antioxidant Activities

In this work we consider a set of 23 Phenolic Compounds, whose structures can
be seen in Tables 1, 2 and 3.

[atom style=scale=0.3, rotate=80][:30]*6((-H3C)-(-CH3)=(*6(-O-(-[:-30]CH3)(-[:30](-[:-30](-[:30](-[:-30](-[:-90]CH3)(-[:30](-[:-30](-[:30](-[:-30](-[:-90]CH3)(-[:30](-[:-30](-[:30](-[:-30](-[:30]CH3)(-[:-90]CH3)))))))))))))—))-=(-H3C)-(-HO)=) |  [atom style=scale=0.3, rotate=80][:30]*6(-(-CH3)=(*6(-O-(-[:-30]CH3)(-[:30](-[:-30](-[:30](-[:-30](-[:-90]CH3)(-[:30](-[:-30](-[:30](-[:-30](-[:-90]CH3)(-[:30](-[:-30](-[:30](-[:-30](-[:30]CH3)(-[:-90]CH3)))))))))))))—))-=(-H3C)-(-HO)=) |  [atom style=scale=0.3, rotate=80][:30]*6((-H3C)-(-CH3)=(*6(-O-(-[:-30]CH3)(-[:30](-[:-30](-[:30](-[:-30](-[:-90]CH3)(-[:30](-[:-30](-[:30](-[:-30](-[:-90]CH3)(-[:30](-[:-30](-[:30](-[:-30](-[:30]CH3)(-[:-90]CH3)))))))))))))—))-=-(-HO)=) |  [atom style=scale=0.3, rotate=80][:30]*6(-(-CH3)=(*6(-O-(-[:-30]CH3)(-[:30](-[:-30](-[:30](-[:-30](-[:-90]CH3)(-[:30](-[:-30](-[:30](-[:-30](-[:-90]CH3)(-[:30](-[:-30](-[:30](-[:-30](-[:30]CH3)(-[:-90]CH3)))))))))))))—))-=-(-HO)=)
---|---|---|---
1 - $\alpha$-Tocoferol | 2 - $\beta$-Tocoferol | 3 - $\gamma$-Tocoferol | 4 - $\delta$-Tocoferol
[atom style=scale=0.5]*6((-O(-[:150]H3C))-(-OH)=(-O(-[:30]CH3))-=(-(-[:30]OH)(=[:150]O))-=-) |  [atom style=scale=0.4]*6((-HO)-=(-OH)-=(-(=[:30](-[:90]*6(-=-(-OH)=-=))))-=-) |  [atom style=scale=0.4]*6((-HO)-(-OH)=(-OH)-=(-(-[:30]OH)(=[:150]O))-=-) |  [atom style=scale=0.4]*6((-HO)-=-(-OH)=(-(-[:30]OH)(=[:150]O))-=-)
5 - Siringico | 6 - Reverastrol | 7 - Ácido Galico | 8 - Ácido Gentissico
[atom style=scale=0.4,
rotate=160]*6((-HO)=(-O(-[:-30]CH3))-=(-(=[:90](-(-[:-30](-[:30](=[:90]O)(-[:-30](=[:30](-[:-30](*6(=-=(-OH)-(-O(-[:-30]CH3))=-)))))))(=[:90]O))))-=-)
5 - Curcumina
Table 1: Structural formulas of Phenolic Compounds used in the present study with antioxidant capacity. [atom style=scale=0.5]*6(-(-OH)=-=(-(=[:150]-[:90](-[:30]OH)(=[:150]O)))-=-) |  [atom style=scale=0.5]*6(-(-OH)=(-O(-[:30]CH3))-=(-(-[:30]OH)(=[:150]O))-=-) |  [atom style=scale=0.5]*6(-(-OH)=(-O(-[:30]CH3))-=(-(=[:150]-[:90](-[:30]OH)(=[:150]O)))-=-) |  [atom style=scale=0.5]*6(-(-OH)=(-O(-[:30]CH3))-=(-(=[:30]O))-=-)
---|---|---|---
10 - Ácido $p$-cumarico | 11 - Ácido Vanilico | 12 - Ácido Ferulico | 13 - Vanilina
[atom style=scale=0.5]*6(-=-(-OH)=(-(=[:150]-[:90](-[:30]OH)(=[:150]O)))-=-) |  [atom style=scale=0.5]*6(-(-OH)=-=(-(-[:30]OH)(=[:150]O))-=-) |  [atom style=scale=0.5]*6(-(-OH)=(-OH)-=(-(-[:150](-[:90]OH)(=[:210]O)))-=-) |  [atom style=scale=0.5]*6(-=-(-OH)=(-(-[:30]OH)(=[:150]O))-=-)
14 - Ácido $o$-cumarico | 15 - Ácido $p$-hidroxibenzoico | 16 - Ácido Protocatequino | 17 - Ácido Salicilico
[atom style=scale=0.5,
rotate=80]*6(-(-OH)=(-O(-[:30]CH3))-=(-(-[:30]-[:90](-[:30]CH3)(=[:150]O)))-=-)
18 - Gengerona Zigerona
Table 2: Structural formulas of Phenolic Compounds used in the present study,
classified with non-active antioxidant capacity.

We used topological indices, described in the subsection 2.2, related to the
compounds in Tables 1 and 2 as vectors of observed characterists to estimate
the discriminant function. There are several indices in the literature, each
one more sensitive to a certain characteristic of the compounds. Table 1
contains all compounds that have proven antioxidant capacity. Table 2 has the
same amount of phenolic compounds as Table 1, however, such compounds either
have low or no antioxidant activity.

In Table 3 we present a group of five phenolic compounds with proven
antioxidant capacity, such compounds were used to apply the discriminant
function, that is, they are our test group.

[atom style=scale=0.5, rotate=80][:30]*6((-HO)-=(*6(-O-(-(*6(=-=(-OH)-(-OH)=-)))-(-OH)–))-=(-HO)-=) |  [atom style=scale=0.5, rotate=80]*6((-HO)=(-OH)-=(-(=[:-30](-[:30](=[:90]O)(-[:-30]O(-[:30]*6(-(-OH)-(-OH)-(-[:60]OH)(-[:0](=[:60]O)(-[:-60]OH))—))))))-=-=) |  [atom style=scale=0.5, rotate=80]*6((-HO)-=-(-(-[:-30](-O(-[:-30](-(=[:-30](-(*6(=-=(-OH)-(-OH)=-)))))(=[:-90]O)))(-[:-90](-[:-30]OH)(=[:-150]O))))=-(-HO)=) |  [atom style=scale=0.5, rotate=80]*6((-HO)=(-OH)-(*6(-O-(=O)-(*6(=-(-OH)=(-OH)-=-))—))=(*6(—O-(=O)-))-=-=)
---|---|---|---
19 - Catequina | 20 - Ácido Clorogenico | 21 - Ácido Rosmanico | 22 - Ácido Elagico
[atom style=scale=0.5,
rotate=-20]*6(-(-OH)=(*6(-(=O)-(-OH)-(-(*6(=-=(-OH)-(-OH)=-)))-O-))-=-(-HO)=)
23 - Quercetina
Table 3: Structural formulas of Phenolic Compounds used in the present study
as a test group, all with antioxidant capacity.

###  2.2 Calculation of Topological Indices

Topological indices are numerical functions of molecular graphs and are
considered important molecular descriptors. According to [26], in addition to
being graph invariant, topological indices do not consider information about
molecular geometry, such as bond lengths, bond angles or twist angles, but
instead encode information about adjacencies of atoms and branches within a
molecule. Also according to the same authors, since the computation of
topological indices uses less resources than the computation of those
molecular descriptors that also take into account molecular geometry,
topological indices have gained considerable popularity and many new
topological indices have been proposed and studied in the literature
specialized in recent years.

We assume that the following aspects of molecules are relevant to the
investigation of the occurrence of antioxidant activity: the type of molecular
branching and the distribution of intramolecular charge. Therefore, in this
work, we use two types of topological indices: the Randic index $\chi^{1}$ and
the topological load indices $G_{4},J_{2}$ and $J_{5}$.

All the topological indices used in this work were obtained with the aid of an
appropriate python algorithm that is available in the supplementary material
that accompanies this article.

####  2.2.1 Randic Index

According to [9], among the several hundred descriptors of molecular
structures based on graphs (see [25]), the Randic index is certainly the most
widely applied in chemistry and pharmacology. This index characterizes the
branching of the molecular graph and was introduced by Milan Randic [23]. If
$E(G)$ denotes the set of edges of the molecular graph $G$, then the Randic
index of $G$ is defined by

$\displaystyle\chi^{1}(G)=\displaystyle\sum_{e_{ij}\in E(G)}(deg_{i}\cdot
deg_{j})^{-1/2},$ (1)

where, $deg_{i}$ and $deg_{j}$ are the degrees of vertices $i$ and $j$,
respectively.

####  2.2.2 Topological Load Indexes

Topological Charge Indices were introduced in the literature by [7] and have
the ability to describe molecular charge distribution. In fact, given a
molecular graph $G$, let $A(e)$ be its adjacency matrix (modified), with the
relative electronegativity of each atom in the main diagonal entries and
$D^{\ast}$ the distance matrix inverse square of $G$, with entries on the main
diagonal taking on a value of zero. Consider $M=[m_{ij}]$ the square matrix of
order $N$ (where $N$ is the number of vertices of $G$) defined by

$\displaystyle M=A(e)\times D^{\ast}$ (2)

and take, for each $i,j$, with $1\leq i,j\leq N$ the charge term $CT_{ij}$
defined by

$\displaystyle CT_{ij}=m_{ij}-m_{ji}.$ (3)

The (valency) topological load indices, $G_{k}$, with $1\leq k\leq N-1$ is
defined by

$\displaystyle G_{k}=\sum^{i=N-1,j=N}_{i,j=i+1}\delta(k,d_{ij})\,\cdot\mid
CT_{ij}\mid,$ (4)

where $\delta$ is the Kronecker delta function ($\delta(p,q)=1$, if $p=q$ and
$0$ otherwise) and $d_{ij}$ denote the entries of the matrix of topological
distance. Remember that in the main diagonal entries of the matrix $A(e)$ the
relative electronegativity of an element $Q$ can be calculated by the formula

$\displaystyle Re^{-}(Q)=\lambda\cdot(e^{-}(Q)-e^{-}(C)),$ (5)

where $e^{-}(Q)$ and $e^{-}(C)$ denote, respectively, the Pauling
electronegativity of the element $Q$ and of the carbon atom $C$. The value
$\lambda$ constitutes the conversion factor. In this work, we consider two
values for $\lambda$: 2.2 (using hydrogen normalization = 0.77) and 3.28
(using chlorine normalization = 2). For this reason, to distinguish the two
respective situations, we will use the symbols $G_{k}^{2.2}$ and
$G_{k}^{3.28}$ to indicate the considered conversion factor.

On the other hand, the index $J_{k}$ is defined by,

$\displaystyle J_{k}=\frac{G_{k}}{N-1},\quad k=1,\cdots,N-1$ (6)

which measures the average value of charge transfer for each chemical bond in
the substance.

_Comp. no._ | _$\chi^{1}$_ | _$G_{4}^{3.28}$_ | _$J_{2}^{3.28}$_ | _$J_{5}^{3.28}$_ | _$G_{4}^{2.2}$_ | _$J_{2}^{2.2}$_ | _$J_{5}^{2.2}$_
---|---|---|---|---|---|---|---
| Antioxidant group training |
1 | 13,72456 | 3,59624 | 0,33965 | 0,05315 | 3,41601 | 0,33924 | 0,05059
2 | 13,19576 | 3,04046 | 0,3174 | 0,05595 | 2,92031 | 0,31262 | 0,05329
3 | 13,19576 | 3,07157 | 0,3174 | 0,05666 | 2,95142 | 0,31262 | 0,05401
4 | 12,77311 | 2,51578 | 0,26956 | 0,05174 | 2,45571 | 0,2603 | 0,049
5 | 5,35317 | 1,55516 | 0,41686 | 0,06137 | 1,14092 | 0,4076 | 0,04784
6 | 5,63077 | 0,59578 | 0,26229 | 0,01284 | 0,53571 | 0,3079 | 0,00993
7 | 4,37239 | 0,67669 | 0,46465 | 0,02008 | 0,59762 | 0,51897 | 0,01968
8 | 3,84159 | 0,84513 | 0,42525 | 0,0125 | 0,73126 | 0,43543 | 0,0125
9 | 9,83297 | 2,81071 | 0,5237 | 0,06254 | 2,5695 | 0,51444 | 0,04775
| Non-antioxidant group training |
10 | 4,0856 | 0,57801 | 0,54821 | 0,04045 | 0,51793 | 0,52525 | 0,02297
11 | 4,40225 | 1,18803 | 0,4169 | 0,03334 | 1,01364 | 0,40595 | 0,03293
12 | 5,06293 | 1,14313 | 0,50661 | 0,05282 | 1,102253 | 0,49734 | 0,03508
13 | 3,99562 | 1,17424 | 0,34116 | 0,0303 | 1,05364 | 0,35315 | 0,025
14 | 3,89175 | 0,76889 | 0,42424 | 0,03879 | 0,76889 | 0,45645 | 0,03237
15 | 3,42492 | 0,71179 | 0,44781 | 0,01066 | 0,53157 | 0,41975 | 0,01016
16 | 4,29869 | 1,01402 | 0,45592 | 0,0397 | 0,83379 | 0,47666 | 0,02371
17 | 3,4309 | 0,74291 | 0,45847 | 0,00694 | 0,56268 | 0,44444 | 0,00694
18 | 5,34935 | 1,43179 | 0,47242 | 0,06691 | 1,37126 | 0,46315 | 0,05508
| Test group |
19 | 8,1157 | 1,84068 | 0,31414 | 0,04914 | 1,84651 | 0,33093 | 0,03719
20 | 9,89595 | 2,43829 | 0,52694 | 0,03863 | 1,83233 | 0,43631 | 0,02781
21 | 9,40614 | 2,59738 | 0,50404 | 0,06857 | 2,11678 | 0,50306 | 0,04338
22 | 8,02602 | 3,50632 | 0,47763 | 0,05734 | 2,66436 | 0,5331 | 0,04339
23 | 7,76622 | 2,52803 | 0.33478 | 0,04599 | 246.000 | 0.40183 | 0,04106
Table 4: Topological indices related to Phenolic Compounds presented in Tables
1 and 2 broken down by group, as used in the Discriminant Analysis.

###  2.3 QSAR Algorithms: Linear Discriminant Analysis

With the objective of detecting the presence of the antioxidant function in
Phenolic Compounds, the QSAR prediction model by Linear Discriminant Analysis
(LDA) was established, a branch of multivariate statistics used in problems of
discrimination and classification of categories or objects, and whose idea
base appear in R. A. Fisher’s seminal work [6], and are currently a widespread
topic (see [12, 14]). In addition, several works use Discriminant Analysis in
establishing QSAR models, see, for example, [5, 3, 16, 2, 17].

Although Discriminant Analysis cannot provide concrete predicted values of the
antioxidant effect, it can determine the likelihood of classifying compounds
as active or inactive for antioxidant activity and thus aid in the discovery
and development of efficient antioxidants. Therefore, we obtained a
discriminant function (DF) that allows a classification between compounds with
active and inactive antioxidant activity. This function linearly depends on
the topological indices considered, which are assumed to be usable in
distinguishing compounds.

In the set of Phenolic Compounds used to obtain the discriminant function, we
separated each compound into two groups: the first formed by compounds with
proven antioxidant activity and the second by compounds without antioxidant
activity. The LDA was applied using the statistical software [22].

##  3 Results and Discussion

The Discriminant Analysis was applied to a training group formed by 18
compounds, of which 09 are classified as active for the antioxidant function
and 09 are classified as non-active for this function. For classification,
experimental data were taken from the literature and consisted of antioxidant
activity determined by DPPH free radical scavenging using lipid peroxidation
inhibitory effects expressed as the concentration of $50\%$ lipid peroxidation
inhibition ($IC_{50}/\ mumol\ ,L^{-}1$). In addition, we considered a test
group formed by 05 Phenolic Compounds, with active antioxidant function,
presented in Table 4. All phenolic compounds considered in this work are found
in foods used in the human diet.

As expected, the discriminant function (DF) found depends linearly on its
variables, which are given by each of the topological indices used in the
study:

$\displaystyle{\bf DF}$
$\displaystyle=-31.023674+4.4419931\,{\bf\chi^{1}}+41.2828013\,{\bf
G_{4}^{3.28}}-151.5908411\,{\bf J_{2}^{3.28}}-$
$\displaystyle-270.9518622\,{\bf J_{5}^{3.28}}-52.570601\,{\bf
G_{4}^{2.2}}+180.1409685\,{\bf J_{2}^{2.2}}+310.4418446\,{\bf J_{5}^{2,2}}$

The DF is used to characterize the presence of antioxidant function in
compounds: a certain compound will be selected as potential antioxidant if
$DF\geq 31.02367$, otherwise, if $DF<31.02367$, then it is classified as not
antioxidant.

The most important fact about the $DF$ function is determining its precision.
For this, we apply the discriminant function on the outer test group. The
result was: $80\%$ precision. Table 5 presents the results obtained for each
compound. In the test group compounds, only 01 compound (Catechin) was
incorrectly classified by the $DF$ function, as non-oxidant ($DF<31.02367$).
For the training set, sorting by the $DF$ function was also very significant.
In fact, $88.9\%$ correct prediction for the antioxidant group, 8 out of 9
classified correctly. In the non-antioxidant group, $100\%$ of correct
prediction, according to Table 5.

| Compostos | DF | Classificacao
---|---|---|---
1 | Catequina | 25.19 | NonAntioxidant
2 | Ácido Clorogenico | 45.18 | Antioxidant
3 | Ácido Rosmanico | 46.83 | Antioxidant
4 | Ácido Elagico | 61.90 | Antioxidant
5 | Quercitina | 31.46 | Antioxidant
6 | Alphatacoferol | 40.77 | Antioxidant
7 | Betatacoferol | 40.20 | Antioxidant
8 | Gammatacoferol | 39.88 | Antioxidant
9 | Deltatacoferol | 38.72 | Antioxidant
10 | Curcumina | 35.79 | Antioxidant
11 | Reverastrol | 36.75 | Antioxidant
12 | Ácido Galico | 39.66 | Antioxidant
13 | Ácido Gentisico | 27.98 | NonAntioxidant
14 | Ácido Sirigico | 36.46 | Antioxidant
15 | Ácido Pcumarico | 22.47 | NonAntioxidant
16 | Ácido Vanilico | 26.43 | NonAntioxidant
17 | Ácido Ferulico | 21.11 | NonAntioxidant
18 | Vanilina | 22.29 | NonAntioxidant
19 | Ácido Ocumarico | 26.06 | NonAntioxidant
20 | Ácido Phidroxibenzoico | 24.65 | NonAntioxidant
21 | Ácido Protocatequino | 30.48 | NonAntioxidant
22 | Ácido Salicilico | 27.17 | NonAntioxidant
23 | Gengerona Zigerona | 21.57 | NonAntioxidant
Table 5: Classification of compounds after application of the Discriminant
Function.

##  4 Conclusions

Discriminant Analysis was applied to establish the QSAR model incorporating
Topological Indices that proved satisfactory to predict the antioxidant
activity of 23 known Phenolic compounds. We made use of topological indices,
which are particularly sensitive to two types of molecular information,
namely, the molecular branching of the compound that is captured by $\chi^{1}$
and the intramolecular charge distributions of the compound through charge
indices $G_{4},J_{2},J_{5}$. We believe that our methodology can be applied to
assess the probability of existence of antioxidant activity in new molecular
compounds, before the application of chemical methods, saving researchers time
and financial resources.

##  5 acknowledgment

The third author thanks the University of California Riverside for the
hospitality.

## References

  * [1] José Manuel Amigó, Vicente Villar Amigó, Antonio Falcó Montesinos and J Gálvez  “Topología molecular”  In _Boletín de la Sociedad Española de Matemática Aplicada. N. 39 (2007)._ Sociedad Española de Matemática Aplicada (SeMA), 2007
  * [2] Lahmadi Ayoub et al.  “A specific QSAR model for proteasome inhibitors from Oleaeuropaea and Ficuscarica”  In _Bioinformation_ 14.7 Biomedical Informatics Publishing Group, 2018, pp. 384
  * [3] Joseph F Contrera, Philip MacLaughlin, Lowell H Hall and Lemont B Kier  “QSAR modeling of carcinogenic risk using discriminant analysis and topological molecular descriptors”  In _Current Drug Discovery Technologies_ 2.2 Bentham Science Publishers, 2005, pp. 55–67
  * [4] Richard D Cramer  “The inevitable QSAR renaissance”  In _Journal of computer-aided molecular design_ 26.1 Springer, 2012, pp. 35–38
  * [5] MTD Cronin and DA Basketter  “Multivariate QSAR analysis of a skin sensitization database”  In _SAR and QSAR in Environmental Research_ 2.3 Taylor & Francis, 1994, pp. 159–179
  * [6] Ronald A Fisher  “The use of multiple measurements in taxonomic problems”  In _Annals of eugenics_ 7.2 Wiley Online Library, 1936, pp. 179–188
  * [7] Jorge Galvez, R Garcia, MT Salabert and Rosa Soler  “Charge indexes. New topological descriptors”  In _Journal of Chemical Information and Computer Sciences_ 34.3 ACS Publications, 1994, pp. 520–525
  * [8] Jorge Galvez, Ramón Garcia-Domenech, V De Julian-Ortiz and Rosa Soler  “Topological approach to analgesia”  In _Journal of Chemical Information and Computer Sciences_ 34.5 ACS Publications, 1994, pp. 1198–1203
  * [9] Ivan Gutman, Boris Furtula and Vladimir Katanić  “Randić index and information”  In _AKCE International Journal of Graphs and Combinatorics_ 15.3, 2018, pp. 307–312  DOI: https://doi.org/10.1016/j.akcej.2017.09.006
  * [10] Lowell H Hall and Lemont B Kier  “Molecular connectivity and substructure analysis”  In _Journal of pharmaceutical sciences_ 67.12 Wiley Online Library, 1978, pp. 1743–1747
  * [11] Corwin Hansch and Toshio Fujita  “p-$\sigma$-$\pi$ Analysis. A Method for the Correlation of Biological Activity and Chemical Structure”  In _Journal of the American Chemical Society_ 86.8 ACS Publications, 1964, pp. 1616–1626
  * [12] Richard A Johnson and Dean W Wichern  “Applied multivariate statistical analysis”  In _New Jersey_ 405, 1992
  * [13] Thiruventhan Karunakaran et al.  “Nitric oxide inhibitory and anti-Bacillus activity of phenolic compounds and plant extracts from Mesua species”  In _Revista Brasileira de Farmacognosia_ 28 SciELO Brasil, 2018, pp. 231–234
  * [14] Ravindra Khattree and Dayanand N Naik  “Multivariate data reduction and discrimination”  In _SAS Institute, Cary, North Carolina_ , 2000
  * [15] Lemont Kier  “Molecular connectivity in chemistry and drug research”  Elsevier, 2012
  * [16] Dmitry A Konovalov, Lyndon E Llewellyn, Yvan Vander Heyden and Danny Coomans  “Robust cross-validation of linear regression QSAR models”  In _Journal of chemical information and modeling_ 48.10 ACS Publications, 2008, pp. 2081–2094
  * [17] Ang Lu et al.  “QSAR study of phenolic compounds and their anti-DPPH radical activity by discriminant analysis”  In _Scientific Reports_ 12.1 Nature Publishing Group, 2022, pp. 1–9
  * [18] Nassira Mahmoudi et al.  “Identification of new antimalarial drugs by linear discriminant analysis and topological virtual screening”  In _Journal of Antimicrobial Chemotherapy_ 57.3 Oxford University Press, 2006, pp. 489–497
  * [19] Natália Martins, Lillian Barros and Isabel CFR Ferreira  “In vivo antioxidant activity of phenolic compounds: Facts and gaps”  In _Trends in Food Science & Technology_ 48 Elsevier, 2016, pp. 1–12
  * [20] Ashley G Muller, Satyajit D Sarker, Imran Y Saleem and Gillian A Hutcheon  “Delivery of natural phenolic compounds for the potential treatment of lung cancer”  In _DARU Journal of Pharmaceutical Sciences_ 27.1 Springer, 2019, pp. 433–449
  * [21] Irina Georgiana Munteanu and Constantin Apetrei  “Analytical methods used in determining antioxidant activity: A review”  In _International Journal of Molecular Sciences_ 22.7 MDPI, 2021, pp. 3380
  * [22] R Core Team  “R: A Language and Environment for Statistical Computing”, 2022  R Foundation for Statistical Computing  URL: https://www.R-project.org/
  * [23] Milan Randic  “Characterization of molecular branching”  In _Journal of the American Chemical Society_ 97.23 ACS Publications, 1975, pp. 6609–6615
  * [24] Muhammed Asin Soobrattee et al.  “Phenolics as potential antioxidant therapeutic agents: mechanism and actions”  In _Mutation Research/Fundamental and Molecular mechanisms of mutagenesis_ 579.1-2 Elsevier, 2005, pp. 200–213
  * [25] Roberto Todeschini and Viviana Consonni  “Handbook of molecular descriptors”  John Wiley & Sons, 2008
  * [26] Alexander Vasilyev and Dragan Stevanovic  “MathChem: a Python package for calculating topological indices”  In _MATCH Commun. Math. Comput. Chem_ 71, 2014, pp. 657–680
  * [27] Muhammad Yasir, Bushra Sultana and Matthew Amicucci  “Biological activities of phenolic compounds extracted from Amaranthaceae plants and their LC/ESI-MS/MS profiling”  In _Journal of Functional Foods_ 26 Elsevier, 2016, pp. 645–656
  * [28] Sae-Rom Yoo et al.  “Simultaneous determination and anti-inflammatory effects of four phenolic compounds in Dendrobii Herba”  In _Natural Product Research_ 31.24 Taylor & Francis, 2017, pp. 2923–2926

##  6 Supplementary Material

    \noindent import networkx as nx\\
    \noindent import numpy as np\\
    \noindent import matplotlib.pyplot as plt\\
    \noindent import networkx.drawing\\
    \noindent import math
    \medskip
    \noindent def mod(x):
      if x>0 or x==0:
         return x
      elif x <0:
         return -x

    def find_shortest_path(graph,start,end, path=[]):
      path=path + [start]
      if start == end:
        return path
      shortest= None
      for node in graph[start]:
        if node not in path:
          newpath=find_shortest_path(graph, node, end, path)
          if newpath:
              if not shortest or len(newpath) < len(shortest):
                  shortest=newpath
      return shortest

    def Randic(MA,GRAU):
      r=0
      for i in range(0,n):
        for j in range(i,n):
          if MA[i][j]!=0 and GRAU[i]!=0 and GRAU[j]!=0:
            r = r + 1/(GRAU[i]*GRAU[j])**0.5
          else:
            r=r
      return r

    def galvez(k,MD,MC, c=0):
        for i in range(0,n):
          for j in range(i,n):
            if k == MD[i][j]:
               c = c + mod(MC[i][j])
        return c

    def Matriz_Adjacencia(n):
        MAD=np.zeros((n,n))
        i=0
        j=0
        resp=0
        print()
        print(’ATENÇÃO APENAS OS NUMEROS ACIMA DA DIAGONAL PRINCIPAL (MAij TEM J>I)
        E A ELETRONEGATIVIDADE NA DIAGONAL PRINCIPAL’)
        print()
        while True:
            i= int(input(’qual o índice da linha?’)) - 1
            while i>n-1 or i==’’:
                  print(’digite um valor menor que %d’ %(n+1))
                  i=int(input(’qual o índice da linha?’)) - 1
            j=int(input(’qual o índice da coluna?’)) - 1
            while j>n-1 or j==’’:
                print(’digite um valor menor que %d’ %(n+1))
                j=int(input(’qual o índice da coluna?’)) - 1
            MAD[i][j]=  float(input(’qual o valor?’))
            if MAD[i][j] == 0:
               resp = input(’voce deseja colocar mas algum valor?’)
               if resp == ’não’ or resp == ’nao’ or resp==’n’:
                   break

        for i in range(0,len(MAD)):
          for j in range (i,len(MAD)):
            MAD[j][i]=MAD[i][j]
        return MAD

    def print_matriz(M):
        print()
        for i in range(0,len(M)):
            for j in range(0,len(M)):
                print(’%2.4f’ %M[i][j], end=’ ’)
            print()

    def verificar(M):
        pup=0
        for i in range(0,len(M)):
          for j in range(0,len(M)):
            if M[i][j]==0:
                pup=pup+1
          if pup == len(M):
              print(’os valores da linha/coluna ’, i+1,’ estao errados. Digite novamente’)
              MAD=Matriz_Adjacencia(len(M))
          pup=0

    #MONTAR MATRIZ ADJACÊNCIA MODIFICADA:

    n=int(input(’numero de elementos do grafo’))
    MAM= Matriz_Adjacencia(n)

    for i in range(0,len(MAM)):
        for j in range (0,len(MAM)):
            if i == j and MAM[i][j]!=0:
                MAM[i][j]= 3.28*(MAM[i][j]-2.55)
            else:
                MAM[i][j]=MAM[i][j]
    print()
    print(’MATRIZ DE ADJACÊNCIA MODIFICADA 3,28:’)
    print()
    print_matriz(MAM)

    #MONTAR MATRIZ DE ADJACÊNCIA:

    MA=np.zeros((n,n))
    for i in range(0,n):
        for j in range(0,n):
            if i!=j:
                MA[i][j]=MAM[i][j]
            else:
                MA[i][j]=0

    #MONTAR DICIONARIO DE LISTA DE ADJACENCIA:

    M=np.zeros((n,n))
    for i in range(0,n):
     p=0
     for j in range(0,n):
      if MA[i][j]!=0:
          M[i][p]=j
          p=p+1

    print()

    dicio={}
    lista=[]
    for i in range(0,n):
        for j in range(0,n):
           if M[i][j]!=0 or j==0:
              lista.append(M[i][j])
           if j==n-1:
             dicio[i]=lista[:]
             lista=[]

    #MONTAR MATRIZ DE DISTÂNCIA

    MD=np.empty((n,n))

    for i in range(0,n):
      for j in range(0,n):
            MD[i][j]=len(find_shortest_path(dicio,i,j)) -1

    print()
    print(’MATRIZ DISTÂNCIA:’)
    print()
    print_matriz(MD)

    #MONTAR MATRIZ DISTÂNCIA MODIFICADA:

    MDM=np.zeros((n,n))

    for i in range(0,n):
      for j in range(0,n):
          if MD[i][j]!=0:
             MDM[i][j]= 1/((MD[i][j]))**2

    #CALCULAR O GRAU DOS VERTICES:

    GRAU=np.zeros(n)

    for i in range (0,n):
        for j in range(0,n):
          GRAU[i]= GRAU[i]+MA[i][j]

    #MUlTIPLICAR A MATRIZ CRIANDO MATRIZ INTERMEDIARIA:

    Mint=np.empty((n,n))

    for i in range(0,n):
        p=0
        for j in range(0,n):
            p=0
            for k in range(0,n):
                p = p + (MAM[i][k])*(MDM[k][j])
            Mint[i][j]=p

    print()
    print(’MATRIZ INTERMEDIÁRIA 3.28’)
    print()
    print_matriz(Mint)

    #MONTAR MATRIZ DE TRANSFERÊNCIA DE CARGA

    MC=np.zeros((n,n))

    for i in range(0,n):
        for j in range(0,n):
            MC[i][j]=Mint[i][j]-Mint[j][i]

    print()
    print(’MATRIZ TRANSFERENCIA DE CARGA: 3.28’)
    print()
    print_matriz(MC)

    #CALCULAR DO INDICE DE RANDIC:

    print()
    print(’1X = %5.5f’ %(Randic(MA,GRAU)))

    #CALCULAR DO ÍNDICE DE GALVEZ:

    print(’Gv,4 3,28 = %5.5f’ %(galvez(4, MD, MC)))
    print(’Jv,2 3,28 = %5.5f’ %((galvez(2,MD,MC))/(n - 1)))
    print(’Jv,5 3,28 = %5.5f’ %((galvez(5,MD,MC))/(n - 1)))

    #MONTAR MATRIZ ADJACÊNCIA MODIFICADA:

    for i in range(0,len(MAM)):
        for j in range (0,len(MAM)):
            if i == j and MAM[i][j]!=0:
                MAM[i][j]= 2.2*(MAM[i][j])/3.28
            else:
                MAM[i][j]=MAM[i][j]
    print()
    print(’MATRIZ DE ADJACÊNCIA MODIFICADA 2,2:’)
    print()
    print_matriz(MAM)

    #MONTAR MATRIZ DE ADJACÊNCIA:

    MA=np.zeros((n,n))
    for i in range(0,n):
        for j in range(0,n):
            if i!=j:
                MA[i][j]=MAM[i][j]
            else:
                MA[i][j]=0

    #MONTAR DICIONARIO DE LISTA DE ADJACENCIA:

    M=np.zeros((n,n))
    for i in range(0,n):
     p=0
     for j in range(0,n):
      if MA[i][j]!=0:
          M[i][p]=j
          p=p+1

    print()

    dicio={}
    lista=[]
    for i in range(0,n):
        for j in range(0,n):
           if M[i][j]!=0 or j==0:
              lista.append(M[i][j])
           if j==n-1:
             dicio[i]=lista[:]
             lista=[]

    #MONTAR MATRIZ DE DISTÂNCIA

    MD=np.empty((n,n))

    for i in range(0,n):
      for j in range(0,n):
            MD[i][j]=len(find_shortest_path(dicio,i,j)) -1

    print()
    print(’MATRIZ DISTÂNCIA:’)
    print()
    print_matriz(MD)

    #MONTAR MATRIZ DISTÂNCIA MODIFICADA:

    MDM=np.zeros((n,n))

    for i in range(0,n):
      for j in range(0,n):
          if MD[i][j]!=0:
             MDM[i][j]= 1/((MD[i][j]))**2

    #CALCULAR O GRAU DOS VERTICES:

    GRAU=np.zeros(n)

    for i in range (0,n):
        for j in range(0,n):
          GRAU[i]= GRAU[i]+MA[i][j]

    #MUlTIPLICAR A MATRIZ CRIANDO MATRIZ INTERMEDIARIA:

    Mint=np.empty((n,n))

    for i in range(0,n):
        p=0
        for j in range(0,n):
            p=0
            for k in range(0,n):
                p = p + (MAM[i][k])*(MDM[k][j])
            Mint[i][j]=p

    print()
    print(’MATRIZ INTERMEDIÁRIA 2.2’)
    print()
    print_matriz(Mint)

    #MONTAR MATRIZ DE TRANSFERÊNCIA DE CARGA

    MC=np.zeros((n,n))

    for i in range(0,n):
        for j in range(0,n):
            MC[i][j]=Mint[i][j]-Mint[j][i]

    print()
    print(’MATRIZ TRANSFERENCIA DE CARGA:2,2’)
    print()
    print_matriz(MC)

    #CALCULAR O INDICE DE RANDIC:

    print()
    print(’1X = %5.5f’ %(Randic(MA,GRAU)))

    #CALCULAR O ÍNDICE DE GALVEZ:

    print($’Gv,4 2,2 = \%5.5f’ \%(galvez(4,MD,MC))$)
    print($’Jv,2 2,2 = %5.5f’ %((galvez(2,MD,MC))/(n - 1))$$)
    print($’Jv,5 2,2 = \%5.5f’ \%((galvez(5,MD,MC))/(n - 1))$$)}


# Construction and Analysis of a Discrete Heat Equation Using Dynamic
Consistency: The Meso-scale Limit

Ronald E. Mickensa Department of Physics, Clark Atlanta University, Atlanta,
GA, 30314  and  Talitha Washington${}^{b},^{*}$ Department of Mathematical
Sciences, Clark Atlanta University, Atlanta, GA, 30314  Corresponding author

###### Abstract.

We present and analyze a new derivation of the meso-level behavior of a
discrete microscopic model of heat transfer. This construction is based on the
principle of dynamic consistency. Our work reproduces and corrects, when
needed, all the major previous expressions which provide modifications to the
standard heat PDE. However, unlike earlier efforts, we do not allow the
microscopic level parameters to have zero limiting values. We also give
insight into the difficulties of constructing physically valid heat equations
within the framework of the general mathematically inequivalent of difference
and differential equations.

###### Key words and phrases:

Heat equation; Dynamic consistency; Random walk; Asymptotic analysis;
Continuum limit

###### 2020 Mathematics Subject Classification:

35K05,; 39A14

##  1\. Introduction

The purpose of this paper is to examine the meso-scale limit of a discrete
micro-level mathematical model constructed to represent simple heat transfer.
The creation of the microscopic model is based on the application of the
principle of dynamic consistency (Mickens (2005, 2015, 2021)). Under the
appropriate mathematical assumptions, we are able to obtain the previous
results of Maxwell (1867), Cattaneo (1948), and Vernotte (1958) regarding the
replacement of the standard heat equation

(1) $u_{t}=Du_{xx},$

by the generalization

(2) $\tau u_{tt}+u_{t}=Du_{xx},$

where $D$ is the temperature diffusion constant, $\tau$ is a time-lag
parameter, and

$u_{t}=\frac{\partial u}{\partial t},\quad
u_{tt}=\frac{\partial^{2}u}{\partial t^{2}},\quad u_{x}=\frac{\partial
u}{\partial x},\quad u_{xx}=\frac{\partial^{2}u}{\partial x^{2}}.$

Note that $(x,t)$ are the one-dimensional space and time independent
variables, and $u=u(x,t)$ is defined over appropriate intervals of $(x,t)$
with suitable boundary conditions and initial values. The need for a
generalization of Eq. (1) comes from the fact that the solutions of Eq. (1)
transmit information at an infinite speed (Ali et al. (2005); Christov et al.
(2005); Dreher et al. (2009); Guyer et al. (1966); Joseph et al. (1989)), a
condition which violates the principle of causality (Ali et al. (2005);
Christov et al. (2005); Dreher et al. (2009); Guyer et al. (1966); Joseph et
al. (1989)).

For convenience, we now provide a resume of our major results:

  * (i)

Based on our formulation of the discrete mirco-level model, the drived first
modified form of the standard heat PDE has the mathematical structure

(3) $\epsilon_{1}u_{tt}+u_{t}=Du_{xx}+(\epsilon_{2}D)u_{xxxx}.$

where $(\epsilon_{1},\epsilon_{2})$ are non-negative parameters.

  * (ii)

Finite-order approximate models of discrete micro-level models can never be
“exact.” By this, we mean that the meso-level equations will always require
knowledge of certain properties of the solutions which in practice and/or
theoretically can not be experimentally measured.

This paper is organized as follows: Section 2 provides a brief historical
summary of work on the standard heat PDE, its major defect, and previous
attempts to resolve the associated issue. It also includes several of the
mathematical models put forth to solve the causality problem and discusses why
they all have failed in this quest. Section 3 gives several general comments
on the modeling process and introduces the concept of dynamic consistency. In
Section 4, we use dynamic consistency to construct perhaps the simplest
discrete mathematical model for heat conduction; it corresponds to a
generalized random walk model. We take this to be our micro-level model. In
Section 5, we show that by using a Taylor series expansion, a meso-level
hybrid difference-differential equation is obtained. Special cases of this
equation includes all previous cases of the modified standard heat equation.
Finally, in Section 6, we summarize our major conclusions and provide reasons
why the usual methods for generalizing the heat PDE may never produce valid
results.

A last comment. We have not tried to reference every relevant publication on
the subjects of this work. Such references are generally well known by active
researchers in these areas and are readily accessible. This view is consistent
with the fact that the current article is not a review paper.

##  2\. Preliminaries

Elementary heat conduction phenomena have been modeled by for approximately
150 years by the linear partial differential given in Eq. (1). However, this
PDE has associated with it two difficulties. First, Eq. (1) allows for the
transformation of information between two separate space points with infinite
speed (Joseph et al. (1989)). Second, this PDE does not provide a valid or
accurate model for many heat conduction problems arising in a broad range of
scientific and engineering situations (Dreher et al. (2009); Guyer et al.
(1966); Joseph et al. (1989)). As a result of the existence of these issues,
there has been a large research effort to resolve and/or explain these issues;
see for example (Ali et al. (2005); Christov et al. (2005); Dreher et al.
(2009); Guyer et al. (1966); Joseph et al. (1989)).

In general, the modifications to Eq. (1) have included adding a second-
derivative term in time. The simplest such model is the Maxwell-Cattaneo
representation (Maxwell (1867); Cattaneo (1948); Vernotte (1958)),

(4) $\tau u_{tt}+u_{t}=Du_{xx}$

where the constant, positive parameter $\tau$ is given the interpretation of a
delay or lag time. Note that this PDE has exactly the same mathematical
structure as the damped wave equation (Ali et al. (2005); Guyer et al. (1966);
Maxwell (1867)). Note that since the dependent variable, $u(x,t)$, represents
the temperature as measured with respect to the absolute Kelvin scale, then
only solutions satisfying the condition

(5) $u(x,t)\geq 0,\quad t\geq 0,$

are physically meaningful. Inspection of Eq. (4) shows that it is a hyperbolic
PDE and, as a consequence, can produce negative valued solutions. So this is a
difficulty that must be overcome. However, Eq. (4) does resolve the issue of
infinite speed for information transfer. We expect this speed to be
approximately

(6) $c=\sqrt{\frac{D}{\tau}}.$

Other modifications of Eq. (1) produce the linear PDE

$\tau u_{tt}+u_{t}=Du_{xx}+D_{1}u_{xxt},$

where $(\tau,D,D_{1})$ are non-negative, constant parameters; see, for
examples, the paper of Guyer et al. (1966). Here, the issues also include
whether the data can actually be measured to obtain the appropriate initial
and boundary conditions.

Our analysis of these modified heat equations leads us to the conclusion that
these efforts have no fundamental groundings in an overall consistent physical
theory. Further, they are related to other efforts involving the examination
of macroscopic limits of microscopic mathematical models (Abeyaratne (2014);
Ascher (2020)) and the validity of expansions in a (supposedly) small
parameter (Mazanov et al. (1974)). An accurate, but colorful way of
characterizing this situation is to look at it from the perspective of a
“wack-a-mole” situation where the things required to resolve one issue
provides the opportunity for the creation of another unresolved issue.

##  3\. Modeling and dynamic consistency

A mathematical model of a system is a symbolic representation incorporating
some aspects of its properties and features. These aspects concern properties
and important features of immediate interest to the modeler and the relevant
scientific community. However, it must not be forgotten that a mathematical
model, in general, can never capture all aspects and details of a system. This
is because at a given moment, no realistic system has all of its properties
known and as a consequence, the unknown ones can not be included in the model.
Further, even if a particular feature or aspect of a system is known, the
currently available mathematical structures may not be able to incorporate
them in an appropriate and viable model.

In the next Section we present arguments which form the bases for the
construction of a discrete microscopic model of heat conduction in one-space
dimension. Our dependent variable, $u$, represents the temperature measured on
the absolute Kelvin scale. This model will be derived by using the principle
of dynamic consistency (Mickens (2005, 2021)). We now provide a concise
definition of dynamic consistency and explain how it may be applied to the
construction of mathematical models.

######  Definition 3.1.

Consider two systems, $S_{1}$ and $S_{2}$. Let $S_{1}$ have some property,
$P$. If $S_{2}$ also has property $P$, then $S_{2}$ is said to be dynamically
consistent to $S_{1}$ with respect to property $P$.

Observe that the systems can be essentially anything. An explicit example is
where $S_{1}$ is a differential equation and $S_{2}$ is a discretization of
$S_{1}$. The case where $S_{2}$ is a finite-difference discretization has been
extensively studied within the context of the “nonstandard finite differences”
methodology (Mickens (2005, 2021)).

Another example is where $S_{1}$ is a physical system and $S_{2}$ is a
mathematical model of $S_{1}$. For this case, $P$ might include any or all of
the following features (Mickens (2005)):

  1. (i)

the dependent variables are non-negative

  2. (ii)

the dependent variables are bounded

  3. (iii)

conservation laws exist

  4. (iv)

the dependent variables have special asymptotic behaviors.

Note that the mathematical model may not have all of the properties $P$. This
is a general consequence of the modeling process. Again, considering the
situation where $S_{1}$ is a differential equation and $S_{2}$ its
discretization, then numerical instabilities may occur in the equations of
discretizations (Mickens (2005, 2021)).

Since the major goals of mathematical modeling are to provide an understanding
of the system, help to make valid predictions of the evolution of the physical
system, and otherwise provide insights into the fundamental nature of the
system, we must be thoughtful in our construction of such models. In
particular, we must always be aware that the results of mathematical modeling
depend strongly on the mathematical structures we use to create the models. An
insightful example is the modeling of single-population systems where discrete
models give solution behaviors entirely different from those provided by
differential equations (Mickens (2015)).

Another difficulty in creating mathematical models is to take as exact
particular functional forms which are based on experimental results valid only
at low concentrations, or a limited range of temperatures, etc. A prime
illustration of this is Fick’s law for heat conduction (Dreher et al. (2009);
Joseph et al. (1989)).

Finally, we give a concise summary of the dynamic consistency methodology for
constructing mathematical models of a system:

  1. (a)

Understand and analyze the system as fully possible;

  2. (b)

from (a), list the major properties and constraints that should be
incorporated into the model;

  3. (c)

using a particular mathematical structure (difference equations, differential
equations, integral equations, etc.), construct a mathematical model that
includes as many as possible of the features and/or items listed in (b).

Note that different mathematical models will arise from the effort to include
some, all, or a separate set of properties of the system and its associated
constraints, such as symmetries and conservation laws. In general, the
application of the dynamic consistency methodology does not yield a unique
mathematical model (Lax (1954); Mickens (2005)], and often this ambiguity can
be used to our advantage. It should be clear that the inclusion of an
increasing number of dynamic, consistent properties should yield better
mathematical models. However, the modeler may not, in general, know all the
important properties of a system. This implies that it is unlikely that the
modeling process will generate equations that are “exact,” either
mathematically or physically.

In the next section, we apply the dynamic consistency methodology to the issue
of constructing a discrete microscopic level model for simple heat conduction.

##  4\. A discrete microscopic level model

Before proceeding to construct the microscopic level model equation, we need
to provide some brief comments on several topics.

First, we assume that it makes sense to divide the physical universe into
three separate relevant domains in terms of space-time characterization. This
is depicted in Fig. (1) and Table (1). The smallest or lowest level is where
the atomic structure of matter dominates our understanding of physical
systems. The top level, sometimes called the continuum limit, is where macro-
scale phenomena is often observed, measured, and interpreted. Generally, at
this level, the graininess of the atomic aspects of matter does not play an
important role. Finally, the meso-level corresponds to considering an a priori
atomic system, but averaged over its atomic distinctness to produce variables
that have macroscopic meaning.

Second, we associate with each of the three levels characteristic time and
space physical scales; see Fig. (1). and Table (1). Thus, $(t_{a},x_{a})$ are
constant parameters associated with the evolution of the system dynamics at
the micro-level; $(\delta_{t},\delta_{x})$ are the constant parameters
associated with physical phenomena at the meso-level; and $(T^{*},L^{*})$ are
the time and length scales for what occurs at the macro-level. Relations
between their magnitudes are indicated in Fig. (1) and Table (1).

Figure 1. Ordering of time and space scales Table 1. Definition of time and space related variables $u_{s}^{r}$ | $u_{m}^{k}$ | $u(x,t)$
---|---|---
$t_{r}=t_{a}\cdot r$ | $t_{k}=\delta_{t}\cdot k$ | $(x,t)\textrm{ continuous}$
$x_{s}=x_{a}\cdot s$ | $x_{m}=\delta_{x}\cdot m$ |
$r=0,1,2,\dots$ | $k=0,1,2,\dots$ |
$s=\textrm{ integers}$ | $m=\textrm{ integers}$ |

Our main purpose is to construct a discrete space-time model for the micro-
level and then extend it to the meso-level. Next, we examine the meso-level
model and interpret it as a modified or corrected version of what is usually
considered to be a mathematical model for the continuum limit.

The space and time variables at the micro-level are discrete and are defined
in Fig. (1) and Table (1). We show that the derived meso-level evolution
equation has to be considered a hybrid differential-difference equation. This
we explain below. Note that $t_{r}$ and $x_{s}$ are discrete time and space
variables. The discrete meso-scale variables and their relationship to the
continuum variables $(t,x)$ will also be clarified later.

To proceed in the derivative of a discrete micro-level model for heat
conduction in one-dimension, we assume that our mathematical model satisfies
the following conditions:

  1. (A)

be a linear difference equation,

  2. (B)

satisfy the constraint that information can only travel with a finite speed in
the discrete space-time;

  3. (C)

be invariant under the parity transformation; i.e., the discrete equation
should not change its form under the interchange

(7) $(s+i)\ \longleftrightarrow\ (s-i),\\\ s=\textrm{ integer},\
i=1,2,3,\dots.$

  4. (D)

not be invariant under time inversion;

  5. (E)

satisfy a non-negativity condition for its solutions;

  6. (F)

possess the property of measurability.

The following is a brief discussion of the implications of these requirements:

  1. (A1)

Assuming that the micro-level equation is linear provides a great
simplification in that there are many more possible nonlinear equations and,
generally, we have no clue as to which one(s) to select.

  2. (B1)

The micro-level is characterized by the time scale to $t_{a}$ and the space
scale to $x_{a}$. Therefore, the speed of information flow, $c$, will be

(8) $c=\frac{x_{a}}{t_{a}}$

  3. (C1)

For a discrete equation the condition given by Eq. (7) ensures parity
conservation (Mickens (2005)).

  4. (D1)

The simplest way to have this condition hold is to have the difference
equation depend only on the discrete times $t_{r}$ and $t_{r+1}$. Physically,
this requirement implies physical dissipation.

  5. (E1)

Since our dependent variable corresponds to temperature on the absolution
Kelvin scale, it must be always non-negative (Dreher et al. (2009); Maxwell
(1867)).

  6. (F1)

By having the property of measurability, we mean that the model should not
contain any parameter and/or dependent variables that can not be readily
measured using standard experimental techniques.

Let $u_{s}^{r}$ be the temperature at discrete time $t_{r}$ and location
$x_{s}$. The simplest difference equation that satisfies all of the above
requirements is given by the expression

(9) $\begin{cases}u_{s}^{r+1}=&pu_{s+1}^{r}+(1-2p)u_{s}^{r}+pu_{s-1}^{r}\\\
&0<p\leq\frac{1}{2}.\end{cases}$

In the work to follow, we specialize to the case $p=1/3$ and Eq. (9) reduces
to

(10)
$u_{s}^{r+1}=\left(\frac{1}{3}\right)\left(u_{s+1}^{r}+u_{s}^{r}+u_{s-1}^{r}\right).$

This special value of $p$ does not change any of our conclusion.

It should be indicated that Eq. (10) is a partial difference equation, first-
order in the discrete-time and second-order in the discrete space variables.
Furthermore, Eq. (10) can be solved exactly by a number of standard techniques
(Mickens (2015)). We will not provide such solutions since their explicit
expressions play no role in the purpose of this article, namely, the
derivation of the meso-level equation and how it should be interpreted.
However, we do assume that there exists an analytic function, $U(x,t)$, such
that

(11) $\begin{cases}u_{s}^{r}=U(x_{s},t^{r}),\\\ s=\textrm{ integer };\quad
r=(0,2,2,\dots);\quad x_{a}\textrm{ and }t_{a}\textrm{ are fixed.}\end{cases}$

Note that Eq. (10) can be rewritten to the form

(12)
$u_{s}^{r+1}-u_{s}^{r}=\left(\frac{1}{3}\right)\left(u_{s+1}^{r}-2u_{s}^{r}+u_{s-1}^{r}\right),$

and from this we can obtain the expression

(13)
$\frac{u_{s}^{r+1}-u_{s}^{r}}{t_{a}}=D\left(\frac{u_{s+1}^{r}-2u_{s}^{r}+u_{s-1}^{r}}{x_{a}^{2}}\right),$

where

(14) $D=\left(\frac{1}{3}\right)\left(\frac{x_{a}^{2}}{t_{a}}\right).$

If $D$ is assumed constant and if the limits

(15) $t_{a}\rightarrow 0,\quad x_{a}\rightarrow 0,\quad D=\textrm{ constant},$

then Eq. (13) becomes the following partial differential equation

(16) $u_{t}=Du_{xx}.$

######  Comment 4.1.

In more detail, the following are the correct limits:

(17) $\begin{cases}t_{a}\rightarrow 0,x_{a}\rightarrow 0;\quad
r\to\infty,\quad s\to\infty;\\\
\left(\frac{1}{3}\right)\left(\frac{x_{a}^{2}}{t_{a}}\right)=D=\textrm{
constant};\\\ t_{a}r=t=\textrm{ fixed},\quad x_{a}s=x=\textrm{
fixed.}\end{cases}$

Observe that Eq. (16) is the standard heat equation (Joseph et al. (1989)).
Also, the procedure used to derive this result opens up a number of important
issues related to the macroscopic limits of microscopic models (Abeyaratne
(2014); Ascher (2020); Lax (1954)). As we shall see in the next section, an
alternative analysis and interpretation is possible.

##  5\. Macro-level approximation

Consider Eq. (10) and as stated above, assume that there exists a function,
$U(x,t)$, with appropriate analytical properties in $x$ and $t$, such that

(18) $u_{s}^{r}=U(x_{x},t_{r}).$

This implies that $U(x,t)$, for fixed $\overline{x}$ and $\overline{t}$,
satisfies the same recurrence/different equation as $u_{s}^{r}$, i.e.,

(19)
$U(x,t+\overline{t})=\left(\frac{1}{3}\right)\left[U(x+\overline{x},t)+U(x,t)+U(x-\overline{x},t)\right]$

Thus, at the meso-level where $(\delta_{x},\delta_{t})$ are the length and
time (physical) scales, we may identify $(\overline{x},\overline{t})$,
respectively, with $(\delta_{x},\delta_{t})$, and obtain

(20)
$U(x,t+\delta_{t})=\left(\frac{1}{3}\right)\left[U(x+\delta_{x},t)+U(x,t)+U(x-\delta_{x},t)\right]$

To proceed, carry out Taylor series expansions of the individual terms in Eq.
(20) and retain only terms up to $O(\delta_{t}^{3})+O(\delta_{x}^{6})$. Doing
this gives, after division by $\delta_{t}$, the expression

(21)
$U_{t}+\left(\frac{\delta_{t}}{2}\right)U_{tt}+O(\delta_{t}^{2})=DU_{xx}+\left[\left(\frac{D}{36}\right)\delta_{x}^{2}\right]U_{xxxx}+O(D\delta_{x}^{4}),$

where

(22) $D=\frac{\delta_{x}^{2}}{3\delta_{t}}.$

Note by assuming that $D$ is fixed in value, we have

(23) $\delta_{t}=O(\delta_{x}^{2}).$

Consequently, to terms of $O(1)$ and $O(\delta_{t})+O(\delta_{x}^{2})$, the
evolution equations for $U(x,t)$ are

(24) $\displaystyle O(1)$ $\displaystyle:\overline{U}_{t}=D\overline{U}_{xx},$
(25) $\displaystyle O(1)+\mathcal{O}(\delta_{t})+O(\delta_{x}^{2})$
$\displaystyle:\left(\frac{\delta_{t}}{2}\right)\overline{U}_{tt}+\overline{U}_{t}=D\overline{U}_{xx}+\left[\left(\frac{D}{36}\right)\delta_{x}^{2}\right]\overline{U}_{xxxx}.$

We have placed bars over the $U$’s to indicate that these equations only yield
approximations to the exact PDE for $U(x,t)$.

Note that the lowest order expansion, i.e., retaining terms to $O(1)$, the
evolution equation is the standard heat conduction PDE. Since

(26) $\partial_{t}=D\partial_{xx}+O(\delta_{t})+O(\delta_{x}^{2}),$

we have

(27) $\displaystyle\partial_{xxxx}=$
$\displaystyle\partial_{xx}\cdot\partial_{xx}$ (28) $\displaystyle=$
$\displaystyle\partial_{xx}\left[\left(\frac{1}{D}\right)\partial_{t}\right]+O(\delta_{t})+O(\delta_{x}^{2}).$

Substituting this latter result into the second term on the right side of Eq.
(25) gives

(29)
$\left(\frac{\delta_{t}}{2}\right)\overline{U}_{tt}+\overline{U}_{t}=D\overline{U}_{x}x+\left(\frac{\delta_{x}^{2}D}{36}\right)\overline{U}_{xxt}.$

This PDE is the same PDE that was derived by Guyer et al. (1966) under
different assumptions.

##  6\. Discussion

In this work, we began by assuming that we have the exact micro-level
evolution equation for simple heat conduction, given by Eq. (10). Assuming
that the solution of this discrete equation can be represented as the values
of an analytic function, $U(x,t)$, evaluated on the lattice of discrete space
and time points, we used Taylor series expansion on this expression at the
meso-level to obtain what we will call modified hybrid-PDEs, see for example
Eq. (29). Dropping the bar over $\overline{U}$, the following two lowest order
PDEs were found

(30) $U_{t}=DU_{xx},$ (31a) $\displaystyle\tau U_{tt}+U_{t}=$ $\displaystyle
DU_{xx}+D_{1}U_{xxt},$ (31b) $\displaystyle\tau U_{tt}+U_{t}=$ $\displaystyle
DU_{xx}+D_{2}U_{xxxx},$

where in this representation, we have expressed them, especially the Eqs.
(31a, 31b), in a general way suitable for use as phenomenological equations,
where $(\tau,D,D_{1})$ are to be determined from experiment. In principle, as
indicated in the previous section, these parameters are determined by the time
and length scales at the meso-level.

Observe that Eq. (30) is the standard heat equation. Within our
interpretation, this PDE should be considered the macroscopic limiting
equation of the microscopic Eq. (10). We do not designate it by the name
continuum limit (Abeyaratne (2014); Ascher (2020)), since no micro-level or
meso-level parameters are required to go to zero. Also, note that within the
context of this work, relationships should exist between the space and time
scales at the micro- and meso-levels, e.g.,

(32) $\delta_{x}=N_{1}x_{a},\quad\delta_{t}=N_{2}t_{a},$

where $N_{1}$ and $N_{2}$ are large integers.

Since neither the micro- or meso-levels parameters go to zero, we expect the
small parameters $(\epsilon_{1},\epsilon_{2})$, where

(33)
$\epsilon_{2}=\frac{\delta_{t}}{T^{*}},\quad\epsilon=\frac{\delta_{x}}{L^{*}},$

to play important roles in the expression of the first modification to the
macroscopic equation. To show this, carry out the scaling of independent
variables in Eq. (29) as follows:

(34) $t=T^{*}\overline{t},\quad x=L^{*}\overline{x},$

where $(L^{*},T^{*})$ are the macro-level scales, and
$(\overline{x},\overline{t})$ are the dimensionless independent variables.
Since the PDE is linear, we need to make no explicit scaling for $U$.
Substituting Eq. (32) into Eq. (29) and multiplying each term by $T^{*}$ gives

(35)
$\left[\left(\frac{1}{2}\right)\left(\frac{\delta_{t}}{T^{*}}\right)\right]\frac{\partial^{2}U}{\partial\bar{t}^{2}}+\frac{\partial
U}{\partial\bar{t}}=\left(\frac{DT^{*}}{(L^{*})^{2}}\right)\frac{\partial^{2}U}{\partial\bar{x}^{2}}+\left(\frac{\delta_{x}^{2}}{36}\right)\left[\frac{T^{*}D}{(L^{*})^{4}}\right]\
\frac{\partial^{4}U}{\partial\bar{x}^{4}}.$

Keeping in mind the definition of $D$ given in Eq. (22), we define
$\overline{D}$, $\epsilon_{1}$, and $\epsilon_{2}$ as follows

(36)
$\overline{D}=\frac{DT^{*}}{(L^{*})^{2}},\quad\epsilon_{1}=\frac{1}{2}\left(\frac{\delta_{t}}{T^{*}}\right),\quad\epsilon_{2}=\left(\frac{\delta_{x}}{L^{*}}\right)^{2},$

where it should be noted that these three parameters are dimensionless, i.e.,
pure numbers. Also, it is expected that $\epsilon_{1}$ and $\epsilon_{2}$ are
small, i.e.,

(37) $0<\epsilon_{1}<<1,\quad 0<\epsilon_{2}<<1.$

With this notation, Eq. (29) in dimensionless form is

(38)
$\epsilon_{1}U_{\bar{t}\,\bar{t}}+U_{\bar{t}}=\overline{D}U_{\bar{x}\bar{x}}+\epsilon_{2}\overline{D}U_{\bar{x}\bar{x}\bar{x}\bar{x}}.$

From the fact that $\delta_{t}=O(\delta_{x}^{2})$, we can expect an infinite
hierarchy of modified PDEs, starting from

(39) $\begin{cases}O(1):U_{t}=DU_{xx},\\\
O(\delta_{t})+O(\delta_{x}^{2}):\textrm{Eq.~{}(\ref{eq37})}\end{cases}$

Figure 2. Hiearchy of modified PDE.

Keep in mind the fact the $N$-th modified PDE contains terms from all the
previous equations. The $N=0$ and $N=1$ cases are the two expressions
presented in Eq. (39).

Observe that at the $N$-th level of the modified PDEs hierarchy, the time
derivatives up to order $(N+1)$ appear. This means also that spatial
derivatives to even order $2(N+1)$ must occur. Thus, for $N=1$, the time
derivatives $U_{t}$ and $U_{tt}$ appear only along with the spatial
derivatives, $U_{xx}$ and $U_{xxxx}$ or $U_{xxt})$. However, for most if not
all one-dimension space systems, the initial condition, $U(x,0)$, is known,
along with related boundary conditions, but $U_{t}(x,0)$ or higher time
derivative are not known, it follows that none of the modified PDEs is
measurable, i.e., their mathematical solutions require information that can
not be physically determined. An important implication is that this way or
methodology of investigating heat condition is invalid or at least not
productive. One way to view this situation is to remember that if we begin
with a discrete space and time model, then the mathematical structure
currently used corresponds to difference equations. However, it is well known
that difference equations do not model exactly differential equations; see the
books by (Mickens (2015, 2021)). In particular, the following relation holds

(40) $\textrm{ difference equation }=\textrm{ infinite-order differential
equation.}$

This can easily be seen by considering the following argument:

  * (1)

Define the (difference) shift and derivative operators by the following
relations

(41) $E_{a}f(x)\equiv f(x+a),\quad Df(x)\equiv\frac{df(x)}{dx}.$

  * (2)

The Taylor expansion of $f(x+a)$ is

(42) $\displaystyle f(x+a)$
$\displaystyle=\overset{\infty}{\underset{k=1}{\sum}}\left(\frac{a^{k}}{k!}\right)D^{k}f(x)$
(43) $\displaystyle=e^{aD}f(x).$

  * (3)

From the definition of the shift operator, it follows that the operators
satisfy

(44) $E_{a}=e^{aD}$

  * (4)

The difference operator is defined by the relation

(45) $\Delta_{a}=E_{a}-1.$

Comparison of Eq. (44) and Eq. (45) gives

(46) $\Delta_{a}=e^{aD}-1.$

  * (5)

Since $\Delta_{a}/a$ and $D$, are the discrete and continuous forms of the
derivative, Eq. (46) shows that they nonlinearly related and any finite order
difference equation corresponds to an infinite order differential equation.
Also, if we solve Eq. (46) for $D$, we obtain

(47) $\displaystyle D$
$\displaystyle\equiv\left(\frac{1}{a}\right)\ln(1+\Delta_{a})$ (48)
$\displaystyle=\left(\frac{1}{a}\right)\left[\Delta_{a}-\frac{\Delta_{a}^{2}}{2}+\frac{\Delta_{a}^{3}}{3}-\frac{\Delta_{a}^{4}}{4}+\cdots\right].$

Consequently, it can be concluded that in general a finite order differential
equation is equivalent to an infinite order difference equation.

In summary, we have used the dynamic consistency methodology (Mickens (2005,
2015)) to construct a discrete space and time microscopic mathematical model
for one-space dimension heat conduction. To obtain solutions to this
difference equation, only initial and/or boundary values are required to be
known. After some manipulations, we were able to derive an infinite number of
mesoscopic PDEs. These PDEs form a hierarchy of equations, each of which
modifies the equations on previous levels; see Fig. (2) However, the modified
equations at the first and higher levels generally require knowledge of time
derivatives for which experimental data is either difficult or impossible to
measure. Our major conclusion is that the procedures used to construct and
justify previous generalizations of the basic heat conduction equation,
$u_{t}=Du_{xx}$, are not valid and that other mathematical structures must be
discovered or created to resolve issues such as infinite speed of propagation
for the transfer of information (Ali et al. (2005); Cattaneo (1948); Christov
et al. (2005); Dreher et al. (2009); Guyer et al. (1966); Joseph et al.
(1989); Maxwell (1867); Vernotte (1958)). We are currently investigating this
problem.

##  7\. Bibliography

## References

  * Abeyaratne (2014) R. Abeyaratne, Macroscopic limits of microscopic models, International Journal of Mechanical Engineering Education, 42 (2014) 185-198.
  * Ali et al. (2005) Y.M. Ali, L.C. Zhang, Relativisitc heat conduction, International Journal of Heat and Mass Transfer, 48 (2005) 2397-2406.
  * Ascher (2020) U.M. Ascher, Discrete processes and their continuous limits, Journal of Dynamics and Games, 7 (2020) 123-140.
  * Cattaneo (1948) C. Cattaneo, Salla conduzione del calore, Atti del Seminario Matematico e Fisico dell’Universita di Moderna, 3 (1948) 3-21.
  * Christov et al. (2005) C.I. Christov, P.M. Jordan, Heat conduction paradox involving second-sound propagation in moving media, Physical Review Letters, 94 (2005) 154301.
  * Dreher et al. (2009) M. Dreher, R. Quintanilla, R. Racke, Ill-posed problems in thermomechanics, Applied Mathematics Letters, 22 (2009) 1374-1379.
  * Guyer et al. (1966) R.A. Guyer, J.A. Ktumhans, Thermal conductivity second-sound and phonon hydrodynamic phenomena in non-metallic crystals, Physical Review, 148 (1966) 778-788.
  * Joseph et al. (1989) D. Joseph, L. Preziosi, Heat waves, Reviews of Modern Physics, 61 (1989) 41-78.
  * Lax (1954) P.D. Lax, Weak solutions of nonlinear hyperbolic equations and their numerical computations, Communications on Pure and Applied Matheamtics, 7 (1954) 159-193.
  * Levy et al. (1961) H. Levy, F . Lessman, Finite Difference Equations, Macmillan, New York, 1961. See Chapter 8.
  * Maxwell (1867) J.C. Maxwell, On the dynamical theory of gases, Philosophical Transactions of the Royal Society of London, 157 (1867) 49-88.
  * Mazanov et al. (1974) A. Mazanov, K.P. Tognetti, Taylor series expansions of delay differential equations - a warning, Journal of Theoretical Biology, 46 (1974) 271-282.
  * Mickens (2005) R.E. Mickens, Dynamic consistency, Journal of Difference Equations and Applications, 11 (2005) 645-653.
  * Mickens (2015) R.E. Mickens, Difference Equations: Theory, Applications, and Advanced Topics, third edition, CRC Press, New York, 2015. See Sections 5.3 and 5.4.
  * Mickens (2021) R.E. Mickens, Nonstandard Finite Difference Schemes: Methodology and Applications, World Scientific, Singapore, 2021. See Sections 9.2 and 9.4.
  * Vernotte (1958) P. Vernotte, Les paradoxes de la theorie continue de l’equation de la chaleur, Comptes Rendus de l’Academie des Sciences, 246 (1958) 3154-3155.


# Gluon shadowing in nuclei and the role of the $c\bar{c}g$ state in the
coherent photoproduction of $J/\psi$ in nucleus-nucleus collisions

Agnieszka Łuszczak<EMAIL_ADDRESS>T. Kosciuszko Cracow
University of Technology, PL-30-084 Cracow, Poland  Wolfgang Schäfer
<EMAIL_ADDRESS>Institute of Nuclear Physics Polish Academy of
Sciences, ul. Radzikowskiego 152, PL-31-342 Cracow, Poland

###### Abstract

In this brief note, we confront our results on diffractive photoproduction of
$J/\psi$ mesons with the putative gluon shadowing ratio defined as
$R_{g}=\sqrt{\sigma(\gamma A\to J/\psi A)/\sigma_{\rm IA.}}$, where
$\sigma_{\rm IA}$ is the result in impulse approximation. Building on our
earlier description of the process in the color-dipole approach, where we took
into account the rescattering of $c\bar{c}$ states only, we demonstrate that
the inclusion of $c\bar{c}g$-Fock states improves the description at small $x$
commonly associated with gluon shadowing.

##  I Introduction

The recent measurements Abelev et al. (2013); Acharya et al. (2021);
Khachatryan et al. (2017); Bursche (2019); Kryshen (2017); Tumasyan et al.
(2023) of exclusive production of $J/\psi$ mesons in ultraperipheral heavy-ion
collisions at the LHC have attracted much attention.

At high energies, the diffractive photoproduction of the $J/\psi$ meson can be
described in the color dipole picture, see for example numerous approaches
Goncalves and Machado (2005); Ayala Filho et al. (2008); Ducati et al. (2013);
Cisek et al. (2012); Lappi and Mantysaari (2013); Sampaio dos Santos and
Machado (2015); Gonçalves et al. (2017); Xie and Chen (2016); Kopp and Machado
(2018); Łuszczak and Schäfer (2018, 2019); Henkels et al. (2020). In the
dipole approach the imaginary part of the diffractive production amplitude in
forward direction takes the form

$\displaystyle\Im m\mathcal{A}(\gamma A\rightarrow
VA;W,\mbox{\boldmath$q$}=0)$ $\displaystyle=$ $\displaystyle
2\int_{0}^{1}dz\int
d^{2}\mbox{\boldmath$r$}\,\Psi^{*}_{V}(z,\mbox{\boldmath$r$})\Psi_{\gamma}(z,\mbox{\boldmath$r$})\,\int
d^{2}\mbox{\boldmath$b$}\,\Gamma_{A}(x,\mbox{\boldmath$b$},\mbox{\boldmath$r$}).$
(1)

Here $\Psi_{\gamma},\Psi_{V}$ stand for the light-front wave functions (LFWFs)
of the the virtual photon and $J/\psi$ meson respectively. The crucial input
is the scattering amplitude
$\Gamma_{A}(x,\mbox{\boldmath$r$},\mbox{\boldmath$b$})$ of a color dipole of
size $r$ in impact parameter space Nikolaev and Zakharov (1991); Mueller
(1990). Here the color dipole are treated as a concrete set of diffraction
scattering eigenstates Good and Walker (1960); Miettinen and Pumplin (1978);
Nikolaev (1981); Nikolaev and Zakharov (1991). For heavy nuclei, the dipole
amplitude can be written in the exponential from familiar from Glauber theory
Glauber (1959).

$\displaystyle\Gamma_{A}(x,\mbox{\boldmath$b$},\mbox{\boldmath$r$})=1-S_{A}(x,\mbox{\boldmath$b$},\mbox{\boldmath$r$})\,,\,\,{\rm{with}}\,\,S_{A}(x,\mbox{\boldmath$b$},\mbox{\boldmath$r$})=\exp\Big{[}-{1\over
2}\sigma(x,\mbox{\boldmath$r$})T_{A}(\mbox{\boldmath$b$})\Big{]}\,.$ (2)

Above, $T_{A}(\mbox{\boldmath$b$})$ is the optical thickness of the nucleus of
mass number $A$ normalized as $\int
d^{2}\mbox{\boldmath$b$}T_{A}(\mbox{\boldmath$b$})=A$. The dipole amplitude
depends on the $\gamma A$ per-nucleon cm-energy $W$ via $x=M_{V}^{2}/W^{2}$.
Here $\sigma(x,\mbox{\boldmath$r$})$ is the dipole cross section for the free
nucleon target. An important feature of diffractive production on the free
nucleon is the “scanning radius” property of the overlap integral Eq. 1, which
entails that the process is dominated by a typical dipole size $r=r_{s}$, in
this case $r_{s}\sim 1/m_{c}$, where $m_{c}$ is the charm quark mass. At small
dipole sizes, we have

$\displaystyle\sigma(x,\mbox{\boldmath$r$})=\frac{\pi^{2}}{N_{c}}\alpha_{s}(\mu^{2})r^{2}\,xg(x,\mu^{2}),\quad\mu^{2}=\frac{A}{r^{2}}+\mu_{0}^{2}\,,$
(3)

which leads to the diffractive cross section proportional to the square of the
proton’s gluon distribution Ryskin (1993),

$\displaystyle\sigma(\gamma p\to J\psi
p)\propto\sigma^{2}(x,r_{s})\propto[xg(x,\mu^{2})]^{2}\,.$ (4)

Clearly this motivates an attempt to find information on the nuclear gluon
distribution from diffractive $J/\psi$ production on nuclear targets. The
strongly absorbing nuclear target however introduces another scale, the
saturation scale (for a review see Gelis et al. (2010))

$\displaystyle Q_{A}^{2}(x,\mbox{\boldmath$b$})={4\pi^{2}\over
N_{c}}\alpha_{s}(Q_{A}^{2})xg(x,Q_{A}^{2})T_{A}(\mbox{\boldmath$b$}),$ (5)

in terms of which the dipole-nucleus S-matrix Eq.2 takes the form Mueller
(1999)

$\displaystyle
S_{A}(x,\mbox{\boldmath$r$},\mbox{\boldmath$b$})=\exp\Big{[}-{1\over
8}Q_{A}^{2}(x,\mbox{\boldmath$b$})r^{2}\Big{]}\,,$ (6)

Therefore the scale

$\displaystyle r_{A}={2\sqrt{2}\over Q_{A}}\,,$ (7)

separates standard hard processes with $r_{s}\ll r_{A}$, for which the leading
twist pQCD factorization works, from processes in the strongly absorptive
saturation regime $r_{s}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}r_{A}$.

An estimate in Ref.Nikolaev et al. (2003) gives
$\langle{Q_{A}^{2}(\mbox{\boldmath$b$})}\rangle\sim 0.9\,\rm{GeV}^{2}$ for
$A^{1/3}=6$ and $x\sim 0.01$. For diffractive $J/\psi$ photoproduction on lead
we are therefore in a regime $r_{s}/r_{A}\sim 0.2$, where we may regard the
rescatterings included in the Glauber-Gribov form as a summation of higher-
twist corrections, but not yet as a fully developed saturated regime.

In Ref.Łuszczak and Schäfer (2019) based on Glauber-Gribov rescattering of a
color dipole were compared to the rapidity dependent cross sections for the
exclusive production of $J/\psi$ in ultraperipheral collisions of lead ions. A
rather good agreement was obtained for various input dipole cross sections
when compared to the results of LHCb, ALICE and CMS at forward rapidities
$|y|\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise
1.0pt\hbox{$>$}}2$.

##  II Beyond the $c\bar{c}$ rescattering

The Glauber form of the dipole amplitude however applies only in range of not
too small $x=x_{A}\sim 0.01$. Indeed, in the dipole picture of high-energy
scattering, higher Fock states containing multiple additional gluons are
responsible for the energy dependence of the dipole cross section. Multigluon
states strongly ordered in light cone momentum fraction $z_{n}\ll
z_{n-1}\ll\dots\ll z_{1}\ll 1$ give rise to the BFKL evolution of the dipole
cross section Mueller (1994); Nikolaev et al. (1994a); Nikolaev and Zakharov
(1994a). The dipole model applies in a limit of large coherence length of the
relevant Fock state. For example, in the limit of large photon energy $\omega$
in the rest frame of the nucleus, the coherence length
$l_{c}=2\omega/M_{c\bar{c}}^{2}$ for $c\bar{c}$ invariant mass
$M_{c\bar{c}}\sim M_{J/\psi}$ becomes much larger than the size of the nucleus
$l_{c}\gg R_{A}$ Kopeliovich and Zakharov (1991); Nikolaev (1992). This
justifies the application of the Glauber approach for color dipoles Mueller
(1990); Nikolaev and Zakharov (1991) for the $c\bar{c}$ component. With
increasing energy, the coherency condition $l_{c}\gg R_{A}$ will be satisfied
not only by the $c\bar{c}$-state, but also by higher
$c\bar{c}g,c\bar{c}gg,\dots$ states. In the language of Glauber–Gribov theory,
these correspond to inelastic shadowing corrections induced by high–mass
diffractive states. Below we restrict ourselves to $c\bar{c}g$ states to
represent high mass diffractive contributions. This is motivated by a
phenomenological success of color dipole approaches such as Bartels et al.
(1999); Golec-Biernat and Łuszczak (2009) which can describe inclusive
diffractive structure functions at HERA just including $q\bar{q}$ and
$q\bar{q}g$ states. We therefore follow previous considerations for inclusive
and diffractive DIS Nikolaev et al. (2007), where a similar approach had been
adopted. We also omit possible contributions of $c\bar{c}q\bar{q}$-Fock states
addressed in Nikolaev et al. (2004).

We very closely follow the approach outlined in Nikolaev et al. (2007) adopted
to the problem at hand in Łuszczak and Schäfer (2022), the details must not be
repeated here.

Assuming the Glauber form to hold at $x\sim x_{A}$, account for the
$c\bar{c}g$ states leads to a dipole amplitude

$\displaystyle\Gamma_{A}(x,\mbox{\boldmath$r$},\mbox{\boldmath$b$})=\Gamma_{A}(x_{A},\mbox{\boldmath$r$},\mbox{\boldmath$b$})+\log\Big{(}{x_{A}\over
x}\Big{)}\Delta\Gamma_{A}(x_{A},\mbox{\boldmath$r$},\mbox{\boldmath$b$})\,,$
(8)

where

$\displaystyle\Delta\Gamma_{A}(x_{A},\mbox{\boldmath$r$},\mbox{\boldmath$b$})$
$\displaystyle=$ $\displaystyle\int
d^{2}\mbox{\boldmath$\rho$}_{1}|\psi(\mbox{\boldmath$\rho$}_{1})-\psi(\mbox{\boldmath$\rho$}_{2})|^{2}\Big{\\{}\Gamma_{A,c\bar{c}g}(x_{A},\mbox{\boldmath$\rho$}_{1},\mbox{\boldmath$\rho$}_{2},\mbox{\boldmath$b$})-\Gamma_{A,c\bar{c}}(x_{A},\mbox{\boldmath$r$},\mbox{\boldmath$b$})\Big{\\}}$
(9)

and

$\displaystyle\Gamma_{A,c\bar{c}g}(x_{A},\mbox{\boldmath$\rho$}_{1},\mbox{\boldmath$\rho$}_{2},\mbox{\boldmath$b$})$
$\displaystyle=$
$\displaystyle\Gamma_{A}(x_{A},\mbox{\boldmath$\rho$}_{1},\mbox{\boldmath$b$}+{\mbox{\boldmath$\rho$}_{2}\over
2})+\Gamma_{A}(x_{A},\mbox{\boldmath$\rho$}_{2},\mbox{\boldmath$b$}+{\mbox{\boldmath$\rho$}_{1}\over
2})$ (10)
$\displaystyle-\Gamma_{A}(x_{A},\mbox{\boldmath$\rho$}_{1},\mbox{\boldmath$b$}+{\mbox{\boldmath$\rho$}_{2}\over
2})\Gamma_{A}(x_{A},\mbox{\boldmath$\rho$}_{2},\mbox{\boldmath$b$}+{\mbox{\boldmath$\rho$}_{1}\over
2})$

The transverse separation of the gluon from quark/antiquark, encoded in the
radial part of the $q\to qg$ light front WF:

$\displaystyle\psi(\mbox{\boldmath$\rho$})={\sqrt{C_{F}\alpha_{s}}\over\pi}{\mbox{\boldmath$\rho$}\over\rho^{2}}\,\mu_{G}\rho
K_{1}(\mu_{G}\rho)\,$ (11)

In Eq.(9) the integration extend over all dipole sizes, including the infrared
domain of large dipoles, where perturbation theory does not apply. Here, we
follow Nikolaev et al. (1994b, 2007) and introduce the parameter $\mu_{G}$
which ensures a the finite propagation radius $R_{c}=1/\mu_{G}\sim 0.2\div
0.3\,\mathrm{fm}$ of gluons. In this work, we employ $\mu_{G}=0.7\,\rm{GeV}$.
We use the running coupling

$\displaystyle\alpha_{s}(r)=\frac{4\pi}{9}\frac{1}{\log\Big{(}\frac{C}{r^{2}\Lambda^{2}}\Big{)}},$
(12)

with $C=1.5$, which we freeze for $r>1/\mu_{G}$. It enters the WF of Eq.(11)
as $\alpha_{s}(\min\\{r,\rho\\})$. In order to obtain the impulse
approximation result, i.e. the production on a quasifree nucleon in the
nucleus, we switch off the nonlinear piece in $\Gamma_{A,c\bar{c}g}$ and put

$\displaystyle\Gamma_{\rm
IA}(x_{A},\mbox{\boldmath$r$},\mbox{\boldmath$b$})={1\over
2}\sigma(x_{A},\mbox{\boldmath$r$})T_{A}(\mbox{\boldmath$b$})\,,$ (13)

We then follow Cisek et al. (2012) and evaluate

$\displaystyle R_{\rm coh}(x)={\displaystyle\int
d^{2}\mbox{\boldmath$b$}\,\Big{|}\langle{J/\psi}|\Gamma_{A}(x,\mbox{\boldmath$r$},\mbox{\boldmath$b$})|{\gamma}\rangle\Big{|}^{2}\over\displaystyle\int
d^{2}\mbox{\boldmath$b$}\,\Big{|}\langle{J/\psi}|\Gamma_{\rm
IA}(x,\mbox{\boldmath$r$},\mbox{\boldmath$b$})|{\gamma}\rangle\Big{|}^{2}}\,,$
(14)

and finally obtain the nuclear cross section as

$\displaystyle\sigma(\gamma A\to J/\psi A;W)=R_{\rm coh}(x)\,\sigma(\gamma
p\to J/\psi p;W)\,.$ (15)

The cross section $\sigma(\gamma p\to J/\psi p;W)$ is taken from our previous
work Łuszczak and Schäfer (2019).

##  III Comparison to experiment

In our numerical calculations we employed an easy to use dipole cross section
designed for the use with heavy flavours from Ref.Golec-Biernat and Sapeta
(2006). The light front wave function overlap is taken from Kowalski et al.
(2006), as in Łuszczak and Schäfer (2019) We refer the reader to these
references for further details.

In Fig.1 we show our results for the total diffractive photoproduction cross
section of $J/\psi$ on lead as a function of $\gamma A$ per-nucleon cm-energy.
The data were obtained by the CMS Tumasyan et al. (2023) and ALICE Acharya et
al. (2021) collaborations. We also show data points extracted by Contreras in
Ref. Contreras (2017) from data obtained in ultraperipheral heavy-ion
collisions. We show separately, the result of the impulse approximation, the
Glauber-Gribov rescattering of $c\bar{c}$ pairs, as well as the cross section
including rescattering of the $c\bar{c}g$ state. To be precise, the solid
curve is obtained by adding

$\displaystyle\sigma(\gamma A\to J/\psi A;W)=(1-f(W))\sigma_{c\bar{c}}(\gamma
A\to J/\psi A;W)+f(W)\sigma_{c\bar{c}+c\bar{c}g}(\gamma A\to J/\psi A;W)$

with

$\displaystyle f(W)={1\over
2}\Big{\\{}1+\tanh\Big{(}\frac{W-W_{0}}{W_{1}}\Big{)}\Big{\\}}\,,$ (17)

with $W_{0}=50\,\rm{GeV}$, $W_{1}=5\,\rm{GeV}$. Clearly the impulse
approximation fails dramatically, which illustrates the scale of nuclear
effects. A large part of the nuclear suppression can be explained by Glauber-
Gribov rescattering of the $c\bar{c}$ state alone. We also observe that the
calculations including the effect of the $c\bar{c}g$ state show an additional
suppression of the nuclear cross section, as required by experimental data.

Figure 1: The total diffractive photoproduction cross section of $J/\psi$ on
${}^{208}\rm Pb$. The data points are taken from Refs.Acharya et al. (2021);
Tumasyan et al. (2023); Contreras (2017) Figure 2: The “gluon shadowing ratio”
$R_{g}(x)=\sqrt{\frac{\sigma(\gamma A\to J/\psi A;W)}{\tilde{\sigma}_{\rm
IA}(\gamma A\to J/\psi A;W)}}\,$ for ${}^{208}\rm Pb$. Notice that the
reference cross section in the denominator is the parametrization taken from
Guzey et al. (2013).

Finally, in Fig.2 we show the so-called gluon shadowing ratio $R_{g}(x)$. The
latter is defined as

$\displaystyle R_{g}(x)=\sqrt{\frac{\sigma(\gamma A\to J/\psi
A;W)}{\tilde{\sigma}_{\rm IA}(\gamma A\to J/\psi A;W)}}\,.$ (18)

It is important to note, that here the impulse approximation cross section in
the denominator, $\tilde{\sigma}_{\rm IA}$ is taken from Ref. Guzey et al.
(2013), as this was also done in the analysis of the CMS collaboration
Tumasyan et al. (2023). This particular ratio is motivated by the leading-
twist limit of the diffractive cross section, where it is expected that

$\displaystyle R_{g}(x)\propto\frac{g_{A}(x,\mu^{2})}{Ag(x,\mu^{2})}\,,$ (19)

and therefore $R_{g}(x)$ is supposed to quantify the suppression (shadowing)
of the per-nucleon glue in the nucleus at small-$x$. Here the scale
$\mu^{2}\sim M^{2}_{J/\psi}/4\sim 2.4\,\rm{GeV}^{2}$.

As can be seen from Fig.2, again including the $c\bar{c}g$ state contributes
to an additional suppression improving the agreement with data. This is also
in qualitative agreement with the results for the nuclear structure function
$F_{2A}$ obtained in Ref. Nikolaev et al. (2007).

At larger $x$, but still well within the range of applicability of the dipole
approach we observe a large discrepancy with experimental data. Note that a
similar discrepancy is also visible in the rapidity dependence of cross
sections from ultraperipheral heavy-ion collisions (UPCs) shown in Łuszczak
and Schäfer (2019, 2022). There we found good agreement with LHCb data, which
have not been analyzed in the form of Fig.2, but obtained a substantial
underprediction of the ALICE forward rapidity data.

We note that this is the case for all dipole models used by us, and as can be
seen from Ref.Tumasyan et al. (2023) the same is true for slightly different
dipole model calculations by other groups. An analysis of exclusive $J/\psi$
production in UPCs based on the NLO collinear factorization approach Eskola et
al. (2023) also finds an ambiguity of gluon distributions fitting either set
of data.

From the point of view of the dipole approach some amount of overprediction of
shadowing in the frozen dipole size approximation is indeed expected Zakharov
(1998) in the range of $x\sim 0.01$. We have accounted for the finite
coherence length only in the most crude way, by including the nuclear form
factor suppression Łuszczak and Schäfer (2019) depending on longitudinal
momentum transfer.

##  IV Conclusions

In this letter we have addressed the diffractive photoproduction of $J/\psi$
mesons at the highest available energies. Data extracted from UPCs are well
described at high energies/small-$x$ after including additional shadowing that
stems from the $c\bar{c}g$ Fock state.

Admittedly, our calculation relies on a model for the behaviour of the $c\to
cg$ light front wave function in the nonperturbative region at large
transverse distances. This is an inevitable element of any color dipole
approach, though. We believe that our modeling of the essentially
nonperturbative physics is well motivated by a phenomenological success of
earlier works in the color dipole approach Nikolaev et al. (1994b, 2007). Our
restriction to the $c\bar{c}g$-system is to some extent justified by the fact,
that diffractive structure functions of the proton measured at HERA are well
described by the inclusion of $q\bar{q}$ and $q\bar{q}g$-states Bartels et al.
(1999); Golec-Biernat and Łuszczak (2009).

It is interesting to ask, whether our inelastic shadowing induced by
$q\bar{q}g$ states can be directly related to a shadowing of the leading
twist, DGLAP evolving gluon distribution of a nucleus. While some early papers
on the color dipole approach Nikolaev and Zakharov (1994a, b); Zakharov (1998)
indeed suggest so, such a conclusion would be ad odds with the treatment of
parton saturation effects along the lines of Balitsky (1996); Kovchegov
(1999). A closer look at the triple Pomeron coupling Nikolaev and Schäfer
(2006); Bartels and Kutak (2008) also casts some doubts on this
interpretation.

A precise relation to gluon shadowing would require a removal of higher-twist
corrections which are also present in our $c\bar{c}g$ contribution. In any
case, DGLAP evolution can only make sense at a hard scale $Q^{2}$ much larger
than the saturation scale. In the case of gluon shadowing, the relevant
saturation scale is the one for an gluon-gluon dipole, which is enhanced by a
factor $C_{A}/C_{F}=9/4$ over the one of the triplet-antitriplet dipole.
Roughly, one expects $\langle{Q^{2}_{8A}(\mbox{\boldmath$b$})}\rangle\sim
2\,\rm{GeV}^{2}$, in a similar ballpark as the hard scale for $J/\psi$
production $Q^{2}\sim M^{2}_{J/\psi}/4\sim 2.4\,\rm{GeV}^{2}$.

## References

  * Abelev et al. (2013) B. Abelev et al. (ALICE), Phys. Lett. B 718, 1273 (2013), eprint 1209.3715.
  * Acharya et al. (2021) S. Acharya et al. (ALICE), Eur. Phys. J. C 81, 712 (2021), eprint 2101.04577.
  * Khachatryan et al. (2017) V. Khachatryan et al. (CMS), Phys. Lett. B 772, 489 (2017), eprint 1605.06966.
  * Bursche (2019) A. Bursche (LHCb), Nucl. Phys. A 982, 247 (2019).
  * Kryshen (2017) E. L. Kryshen (ALICE), Nucl. Phys. A 967, 273 (2017), eprint 1705.06872.
  * Tumasyan et al. (2023) A. Tumasyan et al. (CMS), Phys. Rev. Lett. 131, 262301 (2023), eprint 2303.16984.
  * Goncalves and Machado (2005) V. P. Goncalves and M. V. T. Machado, Eur. Phys. J. C 40, 519 (2005), eprint hep-ph/0501099.
  * Ayala Filho et al. (2008) A. L. Ayala Filho, V. P. Goncalves, and M. T. Griep, Phys. Rev. C 78, 044904 (2008), eprint 0808.0366.
  * Ducati et al. (2013) M. B. G. Ducati, M. T. Griep, and M. V. T. Machado, Phys. Rev. C 88, 014910 (2013), eprint 1305.2407.
  * Cisek et al. (2012) A. Cisek, W. Schäfer, and A. Szczurek, Phys. Rev. C 86, 014905 (2012), eprint 1204.5381.
  * Lappi and Mantysaari (2013) T. Lappi and H. Mantysaari, Phys. Rev. C 87, 032201 (2013), eprint 1301.4095.
  * Sampaio dos Santos and Machado (2015) G. Sampaio dos Santos and M. V. T. Machado, J. Phys. G 42, 105001 (2015), eprint 1411.7918.
  * Gonçalves et al. (2017) V. P. Gonçalves, M. V. T. Machado, B. D. Moreira, F. S. Navarra, and G. S. dos Santos, Phys. Rev. D 96, 094027 (2017), eprint 1710.10070.
  * Xie and Chen (2016) Y.-p. Xie and X. Chen, Eur. Phys. J. C 76, 316 (2016), eprint 1602.00937.
  * Kopp and Machado (2018) F. Kopp and M. V. T. Machado, Phys. Rev. D 98, 014010 (2018), eprint 1806.06701.
  * Łuszczak and Schäfer (2018) A. Łuszczak and W. Schäfer, Phys. Rev. C 97, 024903 (2018), eprint 1712.04502.
  * Łuszczak and Schäfer (2019) A. Łuszczak and W. Schäfer, Phys. Rev. C 99, 044905 (2019), eprint 1901.07989.
  * Henkels et al. (2020) C. Henkels, E. G. de Oliveira, R. Pasechnik, and H. Trebien, Phys. Rev. D 102, 014024 (2020), eprint 2004.00607.
  * Nikolaev and Zakharov (1991) N. N. Nikolaev and B. G. Zakharov, Z. Phys. C 49, 607 (1991).
  * Mueller (1990) A. H. Mueller, Nucl. Phys. B 335, 115 (1990).
  * Good and Walker (1960) M. L. Good and W. D. Walker, Phys. Rev. 120, 1857 (1960).
  * Miettinen and Pumplin (1978) H. I. Miettinen and J. Pumplin, Phys. Rev. D 18, 1696 (1978).
  * Nikolaev (1981) N. N. Nikolaev, Zh.Eksp.Teor.Fiz. 81, 814 (1981).
  * Glauber (1959) R. J. Glauber, _High-Energy Collision Theory_ (Interscience, New York, NY, 1959), pp. 315–414.
  * Ryskin (1993) M. G. Ryskin, Z. Phys. C 57, 89 (1993).
  * Gelis et al. (2010) F. Gelis, E. Iancu, J. Jalilian-Marian, and R. Venugopalan, Ann. Rev. Nucl. Part. Sci. 60, 463 (2010), eprint 1002.0333.
  * Mueller (1999) A. H. Mueller, Nucl. Phys. B 558, 285 (1999), eprint hep-ph/9904404.
  * Nikolaev et al. (2003) N. N. Nikolaev, W. Schäfer, B. G. Zakharov, and V. R. Zoller, J. Exp. Theor. Phys. 97, 441 (2003), eprint hep-ph/0303024.
  * Mueller (1994) A. H. Mueller, Nucl. Phys. B 415, 373 (1994).
  * Nikolaev et al. (1994a) N. N. Nikolaev, B. G. Zakharov, and V. R. Zoller, JETP Lett. 59, 6 (1994a), eprint hep-ph/9312268.
  * Nikolaev and Zakharov (1994a) N. N. Nikolaev and B. G. Zakharov, Z. Phys. C 64, 631 (1994a), eprint hep-ph/9306230.
  * Kopeliovich and Zakharov (1991) B. Z. Kopeliovich and B. G. Zakharov, Phys. Rev. D 44, 3466 (1991).
  * Nikolaev (1992) N. N. Nikolaev, Comments Nucl. Part. Phys. 21, 41 (1992).
  * Bartels et al. (1999) J. Bartels, J. R. Ellis, H. Kowalski, and M. Wüsthoff, Eur. Phys. J. C 7, 443 (1999), eprint hep-ph/9803497.
  * Golec-Biernat and Łuszczak (2009) K. Golec-Biernat and A. Łuszczak, Phys. Rev. D 79, 114010 (2009), eprint 0812.3090.
  * Nikolaev et al. (2007) N. N. Nikolaev, W. Schäfer, B. G. Zakharov, and V. R. Zoller, JETP Lett. 84, 537 (2007), eprint hep-ph/0610319.
  * Nikolaev et al. (2004) N. N. Nikolaev, W. Schäfer, B. G. Zakharov, and V. R. Zoller, JETP Lett. 80, 371 (2004), eprint hep-ph/0408113.
  * Łuszczak and Schäfer (2022) A. Łuszczak and W. Schäfer, SciPost Phys. Proc. 8, 109 (2022), eprint 2108.06788.
  * Nikolaev et al. (1994b) N. N. Nikolaev, B. G. Zakharov, and V. R. Zoller, Phys. Lett. B 328, 486 (1994b), eprint hep-th/9401052.
  * Golec-Biernat and Sapeta (2006) K. J. Golec-Biernat and S. Sapeta, Phys. Rev. D 74, 054032 (2006), eprint hep-ph/0607276.
  * Kowalski et al. (2006) H. Kowalski, L. Motyka, and G. Watt, Phys. Rev. D 74, 074016 (2006), eprint hep-ph/0606272.
  * Contreras (2017) J. G. Contreras, Phys. Rev. C 96, 015203 (2017), eprint 1610.03350.
  * Guzey et al. (2013) V. Guzey, E. Kryshen, M. Strikman, and M. Zhalov, Phys. Lett. B 726, 290 (2013), eprint 1305.1724.
  * Eskola et al. (2023) K. J. Eskola, C. A. Flett, V. Guzey, T. Löytäinen, and H. Paukkunen, Phys. Rev. C 107, 044912 (2023), eprint 2210.16048.
  * Zakharov (1998) B. G. Zakharov, Phys. Atom. Nucl. 61, 838 (1998), eprint hep-ph/9807540.
  * Nikolaev and Zakharov (1994b) N. N. Nikolaev and B. G. Zakharov, J. Exp. Theor. Phys. 78, 598 (1994b).
  * Balitsky (1996) I. Balitsky, Nucl. Phys. B 463, 99 (1996), eprint hep-ph/9509348.
  * Kovchegov (1999) Y. V. Kovchegov, Phys. Rev. D 60, 034008 (1999), eprint hep-ph/9901281.
  * Nikolaev and Schäfer (2006) N. N. Nikolaev and W. Schäfer, Phys. Rev. D 74, 014023 (2006), eprint hep-ph/0604117.
  * Bartels and Kutak (2008) J. Bartels and K. Kutak, Eur. Phys. J. C 53, 533 (2008), eprint 0710.3060.